Commit fcd42a0c authored by Alessandro Mirone's avatar Alessandro Mirone
Browse files

codato raffinamento per fattori per rois

parent 06c9d0e6
......@@ -1391,18 +1391,20 @@ class Functor_Run_Yaml:
print " RETURN "
return
yamltext = self.get_yaml() ##
methodname = self.get_method_name(yamltext)
if self.isrunning():
self.identity_widget.shadok_message_signal.emit("This Job is already running")
return
if methodname != "view_Volume_myavi" :
self.identity_widget.shadok_message_signal.emit("This Job is already running %s"%methodname)
return
if hasattr(self,"sysout"):
oldsysout = sys.stdout
sys.stdout = self.sysout
yamltext = self.get_yaml() ##
methodname = self.get_method_name(yamltext)
where = os.path.join("wizardRUNS", methodname)
if not os.path.exists(where):
......
......@@ -36,8 +36,6 @@ def swissknife_runner( yamltext, where, ret_dico ):
ret_dico["input"] =inputname
ret_dico["stdout"]=stdname
ret_dico["stderr"]=errname
open(inputname,"w").write(yamltext)
stdout = open(stdname,"w")
......@@ -362,6 +360,15 @@ def getMethod(options):
defaults2 = ("1", simplecopy),nmin=1, nmax=40
)
scal["optional_solution"] = Hdf5FilePath(doc=("OPTIONAL : The hdf5 file where scalar volumes has been written\n"
"If given, used for balancing analyzer factors\n"),
defaults2 = ( "", simplecopy )
, fileme=0.5, gnameme=0.5)
scal["target_address"] = Hdf5FilePath(doc=("The hdf5 file where scalar products will be written\n"
"Must be a free name\n"),
......@@ -423,6 +430,30 @@ def getMethod(options):
gv["target_address"] .isresult = 1
gv["swissknife_runner"] = swissknife_runner
# ==============================================================================
vis = collections.OrderedDict()
metodo["VIS_VOLUME"] = vis
vis["getyaml"] = dic2yaml_generic
vis["method_name"] = "view_Volume_myavi"
vis["volume_address"] = Hdf5FilePath(doc=" The hdf5 file and datagroup where volume is \n",
defaults2 = ( gv["target_address"] ,simplecopy),
fileme=1, gnameme=1)
vis["isolevel"] = aNumber( doc="The isolevel : a number between 0 and 1 ( 0 represents minimum, 1, maximum)", float_also=True,
defaults2 = ("0.1", simplecopy),nmin=0, nmax=1.0
)
vis["opacity"] = aNumber( doc="The opacity : a number between 0 and 1 ", float_also=True,
defaults2 = ("0.6", simplecopy),nmin=0, nmax=1.0)
vis["HELP"] = "Visualising teh volume \n"
vis["swissknife_runner"] = swissknife_runner
return metodo
......
......@@ -44,7 +44,7 @@ and even include graphics.
"""
# from mayavi import mlab
from mayavi import mlab
import string
......@@ -561,23 +561,28 @@ def volume_from_2Dimages(mydata):
def view_Volume_myavi(mydata):
"""
volumeaddress : "myfile.hdf5:/path/to/hdf5/group" # the target destination for volume.
volume_address : "myfile.hdf5:/path/to/hdf5/group" # the target destination for volume.
"""
filename, groupname = split_hdf5_address( mydata['volumeaddress'] )
filename, groupname = split_hdf5_address( mydata['volume_address'] )
h5=h5py.File(filename,'r')
Volume = h5[groupname] [:]
h5.close()
h5=None
view_Volume_myavi_(Volume)
h5=None
isolevel = mydata['isolevel']
opacity = mydata['opacity']
view_Volume_myavi_(Volume, isolevel, opacity)
def view_Volume_myavi_(V) :
from mayavi import mlab
def view_Volume_myavi_(V, isolevel, opacity) :
src = mlab.pipeline.scalar_field(V)
mlab.pipeline.iso_surface(src, contours=[V.min()+0.1*V.ptp(), ], opacity=0.6)
mlab.pipeline.iso_surface(src, contours=[V.min()+isolevel*V.ptp(), ], opacity=opacity)
mlab.show()
src = mlab.pipeline.scalar_field(V)
mlab.pipeline.volume(src,vmin=1000.0, vmax=2000.0)
......@@ -1219,14 +1224,10 @@ def superR_getVolume(mydata):
scalprods_address = mydata["scalprods_address"]
scalprods_filename, scalprods_groupname = split_hdf5_address(scalprods_address)
target_filename = mydata["target_filename"]
if os.path.exists(target_filename):
sys.stdout.write("Error : file %s exists already. Remove it yourself\n"%target_filename)
sys.stderr.write("Error : file %s exists already. Remove it yourself\n"%target_filename)
sys.exit(1)
target_address = mydata["target_address"]
target_filename, target_groupname = split_hdf5_address(target_address)
niter = mydata["niter"]
beta = mydata["beta"]
eps = mydata["eps"]
......@@ -1257,8 +1258,15 @@ def superR_getVolume(mydata):
Volume = superr.superr( scalDD, scalDS, scalSS, niter=niter, beta=beta)
h5 = h5py.File(target_filename,"w")
h5["Volume"] = Volume
if os.path.exists(target_filename):
h5 = h5py.File(target_filename,"a")
else:
h5 = h5py.File(target_filename,"w")
if target_groupname in h5.keys():
del h5[target_groupname]
h5[target_groupname] = Volume
h5.close()
......@@ -1318,6 +1326,12 @@ def superR_scal_deltaXimages(mydata):
# defaults to None. If given the integrated image and the average line will be written
# to check the superposability between the foil scans and the sample scans
###
## optional
optional_solution : /mntdirect/_scisoft/users/mirone/WORKS/Christoph/XRSTools/volumes_gasket.h5:/ref408_bis_423_548/Volume0p03
## a solution with dimensions [ZDIM,YDIM,XDIM]
## If given, will be used to balance analyzer factors
If nbin is given the dimensios of the superresolution axis, will be reduced or increased,
by binning together the foil PSFs.
What the program will produce, under *target_address* datagroup, is
......@@ -1338,6 +1352,8 @@ def superR_scal_deltaXimages(mydata):
for a given Z,Y position of the sample, with the reponse function for a given X position
of the foil. The sum over the ROIs is implicitely done.
"""
delta_address = mydata["delta_address"]
......@@ -1357,7 +1373,18 @@ def superR_scal_deltaXimages(mydata):
else:
width = mydata['nbin']
if not mydata.has_key('optional_solution'):
solution = None
else:
solution_address = mydata["orig_delta_address"]
solution_filename, solution_groupname = split_hdf5_address(solution_address)
solution = h5py.File(solution_filename,"r")
solution = solution[solution_groupname][:]
if not mydata.has_key('roi_keys'):
roi_keys = list(h5.keys())
else:
......@@ -1430,25 +1457,76 @@ def superR_scal_deltaXimages(mydata):
m = h5[zscan_keys[0]][roi_keys[0]]["matrix"][:]
YDIM = m.shape[0]
fattori = {}
for i,rk in enumerate(roi_keys):
fattori[rk] = 1.0
if solution is not None:
for rk in roi_keys:
scal_dd=np.array([0.0],"d")
scal_ds = np.array([0.0],"d")
scal_ss = np.array([0.0],"d")
probes = sonde [rk]
SS = np.tensordot( probes, probes, axes = [ [1,2], [1,2] ] )
for iz in range(myZDIM):
zkey = zscan_keys[myrange[iz] ]
m = np.array(h5[ zkey ][ rk ]["matrix"][:],"d")
msum = m.sum(axis=0)
probes = sonde [rk]
assert( probes.shape[1:] == m.shape[1:])
assert( XDIM == probes.shape[0] )
assert( YDIM == m.shape[0] )
plane_contrib = np.tensordot( m, probes, axes = [ [1,2], [1,2] ] )
scal_dd += (m*m).sum()
scal_ds[:] = scal_ds + np.tensordot( plane_contrib , solution[zkey], axes = [ [0,1], [0,1] ] )
scal_ss[:] = scal_ss + np.tensordot( np.tensordot(SS,solution[zkey],axes=[[1],[1]]) , solution[zkey],axes=[[0,1],[1,0]]) )
if nprocs>1:
comm.Allreduce([np.array(scal_ss), MPI.DOUBLE],
[scal_ss, MPI.DOUBLE],
op=MPI.SUM)
comm.Allreduce([np.array(scal_dd), MPI.DOUBLE],
[scal_dd, MPI.DOUBLE],
op=MPI.SUM)
comm.Allreduce([np.array(scal_ds), MPI.DOUBLE],
[scal_ds, MPI.DOUBLE],
op=MPI.SUM)
comm.Barrier()
fattori[rk] = scal_ds/scal_ss
sum = 0.0
for rk in roi_keys:
sum = sum + fattori[rk]*fattori[rk]
for rk in roi_keys:
fattori[rk] = fattori[rk]/np.sqrt( sum/len(toi_keys) )
fattori = fattori/( np.sqrt( np.dot(fattori,fattori) / len(fattori)) )
scalDS = np.zeros( [myZDIM,YDIM,XDIM] ,"d" )
scalDD = 0.0
# Volume = np.zeros( [ZDIM,YDIM,XDIM] ,"d" )
scalSS = np.zeros( [XDIM,XDIM] ,"d" )
for i,rk in enumerate(roi_keys):
if i%nprocs == myrank:
probes = sonde [rk]
scalSS[:] += np.tensordot( probes, probes, axes = [ [1,2], [1,2] ] )
scalSS[:] = scalSS[:] + np.tensordot( probes, probes, axes = [ [1,2], [1,2] ] ) *fattori[rk]*fattori[rk]
if nprocs>1:
comm.Reduce([np.array(scalSS), MPI.DOUBLE],
[scalSS, MPI.DOUBLE],
op=MPI.SUM, root=0)
for iz in range(myZDIM):
comm.Barrier()
for iz in range(myZDIM):
zkey = zscan_keys[myrange[iz] ]
print " process %d analyzing scan : "%myrank , zkey
for rk in roi_keys:
......@@ -1459,28 +1537,40 @@ def superR_scal_deltaXimages(mydata):
msg = " ERROR : the yscan elements have different shapes.\n selects homogeneous scans."
print msg
raise Exception, msg
integrated_images[rk][1] = integrated_images[rk][1]+msum
cornerpos = np.array(h5 [ zkey ][ rk ]["cornerpos"][:])
integrated_images [rk][3] = cornerpos
probes = sonde [rk]
assert( probes.shape[1:] == m.shape[1:])
assert( XDIM == probes.shape[0] )
assert( YDIM == m.shape[0] )
plane_contrib = np.tensordot( m, probes, axes = [ [1,2], [1,2] ] )
scalDS[iz] += plane_contrib
scalDS[iz] = scalDS[iz]+ plane_contrib*fattori[rk]
scalDD += (m*m).sum()
if nprocs>1:
comm.Reduce([np.array(scalSS), MPI.DOUBLE],
[scalSS, MPI.DOUBLE],
op=MPI.SUM, root=0)
comm.Barrier()
h5f.close()
##
## ######################
if nprocs>1:
for n in list(integrated_images.keys()):
if myrank:
......
......@@ -20,6 +20,10 @@ def superr( scalDD, scalDS, scalSS, niter=15, beta=1.0e-8):
scalSS = scalSS.astype("f")
scalDS = scalDS.astype("f")
print " SHAPES ", scalDS.shape, scalSS.shape
Fista ( scalDD, scalDS, scalSS, Volume,niter=niter, beta=beta)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment