#!/usr/bin/env python
# -*- coding: utf-8 -*-
#/*##########################################################################
# Copyright (C) 2011-2017 European Synchrotron Radiation Facility
#
# TDS2EL
# Author :
# Alessandro Mirone, Bjorn Wehinger
#
# European Synchrotron Radiation Facility, Grenoble,France
#
# TDS2EL is developed at
# the ESRF by the Scientific Software staff.
# Principal author for TDS2EL: Alessandro Mirone, Bjorn Wehinger
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation; either version 2 of the License, or (at your option)
# any later version.
#
# TDS2EL is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TDS2EL; if not, write to the Free Software Foundation, Inc., 59 Temple Place,
# Suite 330, Boston, MA 02111-1307, USA.
#
# TDS2EL follows the dual licensing model of Trolltech's Qt and Riverbank's PyQt
# and cannot be used as a free plugin for a non-free program.
#
# Please contact the ESRF industrial unit (industry@esrf.fr) if this license
# is a problem for you.
#############################################################################*/
r"""
* tds2el_extract_maxima is always runned with two arguments ::
tds2el_extract_maxima input_file
* The maxima search operation profits from parallelism ::
mpirun -n 12 tds2el_extract_maxima input_file
* Three different kind of task can be activated by setting to True the corresponding switch parameter
* By default the 3D peaks are searched into a stack of images
* drawmask = True. This option sums up a small numer of images of the stack. It displays them and a mask can be drawn.
The mask is then written to a file and can be used during the peak search to set to zero anomalous pixels(filter_file option).
* review = True. Visual navigation in the found maxima. They can be tagged as good or bad. In the future
the possibility of introducing a mask-per-blob could be considered. But it is maybe better to do it in the reciprocal space.
* example de fichier d' input ::
images_prefix = "G19K45/set1_0001p_"
# Let the following images_prefix_alt to None and dont set it ( default None)
# It is used for reviewing with comparison of same position locations
# versus an alternative dataset
# images_prefix_alt = "G19K45/set1_0001p_"
images_postfix = ".cbf"
## peaks are searched, having a radius peak_size
# and with maximum > threshold_abs and > image_max*threshold_rel
threshold_abs=40000
threshold_rel=0.00
peak_size = 2
spot_view_size = 10
# Filtering, OPTIONAL or set it to None for no mask filtering of bad pixels
# filter_file = None
filter_file = "filter.edf"
## medianize can be 1 or 0 ( True or False)
## The peak is always searched into a stack
## of spot_view_size + 1 + spot_view_size slices
## When the option is activated each slice will be beforehand
## subtracted with the median of the slices going from the spot_view_size^th preceeding
## to the spot_view_size^th following. Spot_view_size should be greaer than peak_size
medianize = 1
## Let this advanced option as it is for now, or dive into the code.
## The rationale here would be that in the first step on works
## locally on spot_view_size + 1 + spot_view_size slices,
## in the second step one can work globally on the found peaks
globalize = (0,0)
## The file where found peaks are written
## or, when reviewing, read.
file_maxima = "maxima.p"
### if yes, a graphical window allows for the selection of the mask
### Closing the window will propose to write on newmask.edf
drawmask = False
## Allows to have a simultaneous and animated view of the blobs.
## Animated thumnails.
## They can be selected of deselected. When clicking on the
## thumbnails a larger 3D explorer is charged with
## a (review_size+1+review_size )^3 stack
review = False
review_size = 50
* The following image illustrates the use of the viewer. Peaks can be selected and deselected. When closing the windows
you are asked if yes or not you want to save the maxima on a new file. The file can be reused for the further steps: geometry alignement
and elastic constants fitting. The alignements steps create a new file where peaks are tagged by hkl. The image below is the run of the
viewer on such a file, hkl is displayed. Each icon is animated movie sweeping the small (2*spot_view_size+1)**3 cube around the peak in image space: each frame correspond
to an experimental image. If you let the mouse pointer idle on an icon a toolitip appear mentioning further metadata od the spot. If you click the 3D stack viewer from Silx is charged
with a larger (2*review_size+1)**3 data volume of the data. Whe you have loaded an image and working on it, its checkbox is red-framed to remember what you are working on. The alt-data
if any will be displayed as an additional companion icon on the right of the checkbox
* The thumbnails are generated in extraction phase using a Silx widget. At present the widget needs to have acces to a graphical display and interact with it,
this might be a performance issue if the extraction is runned remotely in some cases of slow network. With NX ( no-machine) the performances are
good even remotely.
.. image:: viewer.png
"""
import glob
# import fftw3f
import math
import sys
from silx import sx
import silx
import silx.gui
from silx.gui.plot import PlotWidget
from silx.gui import qt
from silx.gui.colors import Colormap
import h5py
import pylab
import random
# import Fit
import minimiser
import numpy
import os
import scipy
## from skimage.feature import peak_local_max
import pickle
if(sys.argv[0][-12:]!="sphinx-build"):
import spotpicker_cy
global observed
observed = None
try:
from mpi4py import MPI
myrank = MPI.COMM_WORLD.Get_rank()
nprocs = MPI.COMM_WORLD.Get_size()
procnm = MPI.Get_processor_name()
comm = MPI.COMM_WORLD
print( "MPI LOADED , nprocs = ", nprocs)
except:
nprocs=1
myrank = 0
print( "MPI NOT LOADED ")
# import sift
from math import pi, sin, cos
try :
import fabio
except :
print( 'Warning : fabio module could not be initialised')
import numpy as np
import sys, time
import scipy.ndimage as ndi
qapp = qt.QApplication.instance() or qt.QApplication([])
w=PlotWidget()
timer = qt.QTimer()
tmpw=PlotWidget()
__version__ = 0.7
def mylocalmaxs(image , min_dist=1 , tabs=None, trel=None , nobord=False ):
size = 2 * min_dist + 1
# print " filter "
# mask = mask.swapaxes(0, i)
image_max = np.array(image,"f")
tmp0 = np.array(image,"f")
spotpicker_cy.window_maximum( image_max, tmp0, 0, min_dist )
spotpicker_cy.window_maximum( tmp0, image_max, 1, min_dist )
spotpicker_cy.window_maximum( image_max, tmp0, 2, min_dist )
#print " filter2 "
#image_max = ndi.maximum_filter(image, size=size, mode='constant')
#print " O? "
mask = image == image_max
# if nobord:
# for i in range(mask.ndim):
# mask = mask.swapaxes(0, i)
# mask[:min_dist] = mask[-min_dist:] = False
# mask = mask.swapaxes(0, i)
if nobord:
mask[:nobord] = mask[-nobord:] = False
thresholds = []
if tabs is None:
tabs = image.min()
thresholds.append(tabs)
if trel is not None:
thresholds.append(trel * image.max())
if thresholds:
mask &= image > max(thresholds)
coord = np.nonzero(mask)
coord = np.column_stack(coord)
return coord
[docs]class PeakoScope(qt.QTabWidget):
def __init__(self,maxima_data, review_size, flist, flist_alt, parent = None):
super(PeakoScope, self).__init__(parent)
files_name = maxima_data.keys()
files_name.sort()
picchi = [ ]
for fn in files_name:
d=maxima_data[fn]
for pos, dico in d.items():
dico["fn"] = fn
picchi.append(dico)
nps = len(picchi)
ntabs = nps/100
if ntabs*100<nps:
ntabs+=1
self.tabs={}
self.maxima_data = maxima_data
self.review_size = review_size
self.flist = flist
self.flist_alt = flist_alt
for i in range(ntabs):
nmin = i*100
nmax = min((i+1)*100, nps)
mypeaks = picchi[nmin:nmax]
self.tabs[i] = qt.QWidget()
self.addTab( self.tabs[i], os.path.basename(mypeaks[0]["fn"]) +"--"+os.path.basename(mypeaks[-1]["fn"]) )
self.tabFill(self.tabs[i], mypeaks )
img = fabio.open(flist[0])
self.DIM1orig, self.DIM2orig = img.data.shape
def tabFill(self, tab, peaks):
grid = qt.QGridLayout()
for i in range(len(peaks)):
print( " aggiungo ", i)
r,c = i/10, i%10
grid.addWidget(PmapDisplayer(peaks[i], self.flist, self.flist_alt, self.review_size),r,c)
tab.setLayout(grid)
[docs] def closeEvent(self, event):
print("event")
reply = qt.QMessageBox.question(self, 'Message',
"Now exiting, also writing on newpeaks.p?\n ", qt.QMessageBox.Yes, qt.QMessageBox.No, qt.QMessageBox.Discard )
if reply == qt.QMessageBox.No:
sys.exit(0)
event.accept()
elif reply == qt.QMessageBox.Yes:
f = open("newpeaks.p","w")
pickle.dump( self.maxima_data ,f )
pickle.dump( [self.DIM1orig, self.DIM2orig] ,f )
f.close()
sys.exit(0)
else:
event.ignore()
[docs]class myPlot2D(sx.Plot2D):
[docs] def closeEvent(self, event):
print("event")
data = self.getMaskToolsDockWidget().getSelectionMask()
if len(data)==0:
sys.exit(0)
reply = qt.QMessageBox.question(self, 'Message',
"Now exiting, also writing on newmask.edf?\n ", qt.QMessageBox.Yes, qt.QMessageBox.No, qt.QMessageBox.Discard )
if reply == qt.QMessageBox.No:
sys.exit(0)
event.accept()
elif reply == qt.QMessageBox.Yes:
i=fabio.edfimage.edfimage()
i.data=(1-data).astype("f")
i.write("newmask.edf")
sys.exit(0)
else:
event.ignore()
def mainDrawmask( flist):
wscope = myPlot2D()
d=0
for f in flist[:20]:
d=d+fabio.open(f).data
wscope.addImage(d)
wscope.show()
qapp.exec_()
[docs]class MyStack(sx.StackView):
[docs] def closeEvent(self, event):
self.hide()
event.ignore()
def mainReview(params, flist, flist_alt):
maxima_data = pickle.load( open(params.file_maxima, 'r'))
# files_name = maxima_data.keys()
# files_name.sort()
global wscope
wscope = MyStack()
wscope.setStack(numpy.random.random(64**3).reshape(64, 64, 64))
wscope.show()
ps=PeakoScope(maxima_data, params.review_size, flist, flist_alt, parent = None)
ps.show()
qapp.exec_()
# def __init__(self, *args):
# super(CentralWidget, self).__init__(*args)
# self.setStyleSheet("background-color: rgb(255,0,0); margin:5px; border:1px solid rgb(0, 255, 0); ")
[docs]class ScaledLabel(qt.QLabel):
def __init__(self, *args, **kwargs):
qt.QLabel.__init__(self)
self._pixmap = qt.QPixmap(self.pixmap())
self.setSizePolicy(qt.QSizePolicy.Expanding, qt.QSizePolicy.Expanding)
[docs] def resizeEvent(self, event):
self.setPixmap(self._pixmap.scaled(
self.width(), self.height(),
qt.Qt.KeepAspectRatio))
def convertArrayToQImage(image):
image = np.array(image, copy=False, order='C', dtype=np.uint8)
qimage = qt.QImage(
image.data,
image.shape[1],
image.shape[0],
image.strides[0], # bytesPerLine
qt.QImage.Format_RGB888).copy()
return qimage
[docs]class PmapDisplayer(qt.QFrame):
def __init__(self, peak, flist, flist_alt, rsize):
self.peak=peak
self.rsize = rsize
self.flist = flist
self.flist_alt = flist_alt
try:
pmaps, pmaps_alt = peak["pmap"]
except:
pmaps = peak["pmap"]
pmaps_alt = None
peak["accepted"]
qt.QFrame.__init__(self)
# self.setAttribute(Qt.Qt.WA_NoSystemBackground)
self.setAutoFillBackground(True)
self.pmaps = []
for p in pmaps:
ima = convertArrayToQImage(p)
self.pmaps.append(qt.QPixmap(ima))
self.resize(self.pmaps[0].width(), self.pmaps[0].height())
print( " size " , self.pmaps[0].width(), self.pmaps[0].height())
self.pmaps_alt = []
if pmaps_alt is not None:
for p in pmaps_alt:
ima = convertArrayToQImage(p)
self.pmaps_alt.append(qt.QPixmap(ima))
self.grid = qt.QGridLayout(self)
# if "hkl" in peak:
# print "%-- s y,x = %s"%(peak["fn"], peak["pos"] )
# self.setToolTip("%s y,x = %s hkl=%s"%(peak["fn"], peak["pos"], peak["hkl"] ) )
# else:
# print "%s y,x = %s"%(peak["fn"], peak["pos"] )
# self.setToolTip("%s y,x = %s"%(peak["fn"], peak["pos"] ) )
self.label = qt.QLabel( )
self.label.setSizePolicy(qt.QSizePolicy.Expanding, qt.QSizePolicy.Expanding)
if "hkl" in peak:
print( "%-- s y,x = %s"%(peak["fn"], peak["pos"] ))
self.label.setToolTip("%s y,x = %s hkl=%s"%(peak["fn"], peak["pos"], peak["hkl"] ) )
else:
print( "%s y,x = %s"%(peak["fn"], peak["pos"] ))
self.label.setToolTip("%s y,x = %s"%(peak["fn"], peak["pos"] ) )
self.checkbox = qt.QCheckBox()
self.checkbox.setChecked( self.peak["accepted"] )
self.checkbox.stateChanged.connect(self.changecheck)
if "hkl_rounded" in peak:
self.grid.addWidget(self.label,0, 0, 2,1)
# self.tag_w = qt.QWidget()
# self.tag_w.setSizePolicy(qt.QSizePolicy.Expanding, qt.QSizePolicy.Expanding)
# self.grid.addWidget(self.tag_w,0,1)
# self.grid2 = qt.QGridLayout(self.tag_w)
self.tag = qt.QLabel("%s"% peak["hkl_rounded"] )
# self.tag.setSizePolicy(qt.QSizePolicy.Expanding, qt.QSizePolicy.Expanding)
# self.grid2.addWidget(self.checkbox,0,0)
#self.grid2.addWidget(self.tag,1,0)
self.grid.addWidget(self.checkbox,0,1)
self.grid.addWidget(self.tag ,1,1)
else:
self.grid.addWidget(self.label,0, 0, 1,1)
self.grid.addWidget(self.checkbox,0,1)
if len(self.pmaps_alt):
self.label_alt = qt.QLabel( )
self.label_alt.setSizePolicy(qt.QSizePolicy.Expanding, qt.QSizePolicy.Expanding)
self.grid.addWidget(self.label_alt,0,2)
fn = peak["fn"]
pp = self.flist.index(fn)
fn_alt = self.flist_alt[pp]
if "hkl" in peak:
print( "%-- s y,x = %s"%(peak["fn"], peak["pos"] ))
self.label_alt.setToolTip("%s y,x = %s hkl=%s"%(fn_alt, peak["pos"], peak["hkl"] ) )
else:
print( "%s y,x = %s"%(peak["fn"], peak["pos"] ))
self.label_alt.setToolTip("%s y,x = %s"%(fn_alt, peak["pos"] ) )
self.imgnumber=0
self.changeFrame()
self.label.resize(self.pmaps[0].width(), self.pmaps[0].height())
self.label.setPixmap(self.pmaps[0].scaled( self.label.width(), self.label.height(), qt.Qt.KeepAspectRatio))
timer.timeout.connect(self.changeFrame)
self.fixborder()
[docs] def mousePressEvent(self, QMouseEvent):
print( "PRESS " , QMouseEvent.pos())
print( "checkbos pos ", self.checkbox.pos())
m_pos = QMouseEvent.pos().x()
c_pos = self.checkbox.pos().x()
pos = self.peak["pos"]
Gi = self.peak["Gi"]
slice0 = slice( max(0,pos[0]-self.rsize), pos[0]+self.rsize+1 )
slice1 = slice( max(0,pos[1]-self.rsize), pos[1]+self.rsize+1 )
print( slice0)
print( slice1)
if m_pos< c_pos:
print (" PRENDO UNO STACK ", len(self.flist), Gi, max( 0, Gi-self.rsize),min(Gi+self.rsize+1,len(self.flist) ) )
stack=np.array([ fabio.open(fn).data[slice0, slice1 ] for fn
in self.flist[ max( 0, Gi-self.rsize):min(Gi+self.rsize+1,len(self.flist) ) ]])
elif self.flist_alt is not None:
print (" PRENDO UNO STACK from ALT", len(self.flist_alt), Gi, max( 0, Gi-self.rsize),min(Gi+self.rsize+1,len(self.flist_alt) ) )
stack=np.array([ fabio.open(fn).data[slice0, slice1 ] for fn
in self.flist_alt[ max( 0, Gi-self.rsize):min(Gi+self.rsize+1,len(self.flist) ) ]])
print( " RISULTATO " , stack.shape)
global wscope
wscope.setStack(stack)
wscope.show()
global observed
if observed is not None:
observed.fixborder()
observed = self
self.checkbox.setStyleSheet("margin:1px; border:1px solid rgb(255, 0, 0); ")
def changecheck(self, status):
print( status )
self.peak["accepted"] = (status!=0)
self.fixborder()
def fixborder(self):
if self.peak["accepted"] :
self.checkbox.setStyleSheet("margin:1px; border:1px solid rgb(0, 255, 0);")
else:
self.checkbox.setStyleSheet("margin:1px; border:1px solid rgb(200, 200, 200); ")
def changeFrame(self):
if not self.peak["accepted"]:
self.imgnumber = len(self.pmaps)/2
self.label.setPixmap(self.pmaps[self.imgnumber%len(self.pmaps)].scaled( self.label.width(), self.label.height(), qt.Qt.KeepAspectRatio))
if len(self.pmaps_alt):
self.label_alt.setPixmap(self.pmaps_alt[self.imgnumber%len(self.pmaps)].scaled( self.label.width(), self.label.height(), qt.Qt.KeepAspectRatio))
## self.label.setPixmap(self.pmaps[self.imgnumber%len(self.pmaps)])
self.imgnumber+=1
# def paintEvent(self,event):
# self.setAttribute(Qt.Qt.WA_NoSystemBackground)
def mainRead(Filter_fname,flist,flist_alt, t_abs, t_rel, binning, peak_size ,spot_view_size, medianize, globalize, file_maxima):
img = fabio.open(flist[0])
DIM1orig, DIM2orig = img.data.shape
data = img.data[::binning, ::binning]
dim1, dim2 = data.shape
print( 'Image Dimension :', dim1, dim2)
total = len(flist)
if Filter_fname is not None:
print( 'Loading images filter ...')
Filter = fabio.open(Filter_fname ).data
assert( Filter.shape == (dim1,dim2) )
else:
Filter = np.ones( (dim1,dim2) ,"f")
massimo=0.0
i=0
Nring = 2*spot_view_size+1
middle = spot_view_size
ring = np.zeros([Nring, data.shape[0], data.shape[1]],"f")
ring_orig = np.zeros([Nring, data.shape[0], data.shape[1]],"f")
ring_median = np.zeros([Nring, data.shape[0], data.shape[1]],"f")
ring_images_name = [None]*Nring
ring_alt = np.zeros([Nring, data.shape[0], data.shape[1]],"f")
nblockfl = len(flist)/nprocs+1
res={}
count = 0
count_median = 0
comm.Barrier()
marge = Nring-middle
tutti = []
for iim, fname in enumerate(flist[:]):
if flist_alt is not None:
# print flist_alt
fname_alt = flist_alt[iim]
if ((iim/nblockfl)%nprocs != myrank) and ( ((iim-marge)/nblockfl ) %nprocs != myrank ) and ( ((iim+marge) /nblockfl ) %nprocs != myrank ) :
continue
newim = (fabio.open(fname).data.astype(np.float32)*Filter)[::binning, ::binning]
if medianize:
ring_median[count_median%Nring] = newim
count_median +=1
imorig = newim
newim = newim -np.median(ring_median[:count_median], axis=0)
print( 'Working on image %s' % fname )
ring[count%Nring] = newim
if medianize :
ring_orig[count%Nring] = imorig
if flist_alt is not None:
ring_alt[count%Nring] = (fabio.open(fname_alt).data.astype(np.float32)*Filter)[::binning, ::binning]
ring_images_name[count%Nring] = fname
if count==Nring-1:
iplan=middle
elif count>=Nring:
iplan+=1
else:
count+=1
continue
count+=1
centered_ring = np.roll( ring, middle-(iplan%Nring),axis=0)
if medianize :
centered_ring_orig = np.roll( ring_orig, middle-(iplan%Nring),axis=0)
if flist_alt is not None:
centered_ring_alt = np.roll( ring_alt, middle-(iplan%Nring),axis=0)
# picchi = peak_local_max( centered_ring[middle-1:middle+2] , min_distance=peak_size , exclude_border = True) # threshold_abs=t_abs,threshold_rel=t_rel)
# picchi = peak_local_max( centered_ring , min_distance=peak_size , exclude_border = True, threshold_abs=t_abs,threshold_rel=t_rel)
picchi = mylocalmaxs( centered_ring , min_dist = peak_size , nobord = spot_view_size, tabs=t_abs , trel = t_rel)
picchiI = []
picchiPmap = []
picchiPmap_alt = []
picchi_again=[]
for i,j,k in picchi:
newpos = np.array([iim,j,k])
if len(tutti) and np.max(np.abs( np.array(tutti)-newpos), axis=1).min()<10:
continue
tutti.append([iim,j,k])
picchiI.append(centered_ring[i,j,k ] )
picchi_again.append( (i,j,k) )
tok=[]
tok_alt=[]
for II in range( len(centered_ring)):
cm = Colormap( name='viridis', normalization='log' )
# tmpw=PlotWidget()
if not medianize:
rgbA = cm.applyToData( centered_ring[II,j-spot_view_size:j+spot_view_size,k-spot_view_size:k+spot_view_size ])
# tmpw.addImage(centered_ring[II,j-spot_view_size:j+spot_view_size,k-spot_view_size:k+spot_view_size ], colormap={'name':'viridis', 'autoscale':True, 'normalization':'log'} )
else:
rgbA = cm.applyToData( centered_ring_orig[II,j-spot_view_size:j+spot_view_size,k-spot_view_size:k+spot_view_size ] )
# tmpw.addImage(centered_ring_orig[II,j-spot_view_size:j+spot_view_size,k-spot_view_size:k+spot_view_size ], colormap={'name':'viridis', 'autoscale':True, 'normalization':'log'} )
adata = rgbA[:, :, :3]
tok.append(adata)
if flist_alt is not None:
# tmpw=PlotWidget()
rgbA = cm.applyToData( centered_ring_alt[II,j-spot_view_size:j+spot_view_size,k-spot_view_size:k+spot_view_size ])
adata = rgbA[:, :, :3]
tok_alt.append(adata)
picchiPmap.append(tok)
if flist_alt is not None:
picchiPmap_alt.append(tok_alt)
# ima = silx.gui._utils.convertArrayToQImage(adata)
# p=qt.QPixmap(ima)
myimage = ring_images_name[ (iplan%Nring) ]
Giplan = iim -((Nring-1) -middle)
if (myrank==1) :print ("==== ", Giplan, ring_images_name[iplan%Nring])
# if (Giplan/nblockfl)%nprocs==myrank:
# nameout=params.airscattering_debug_dir+"/"+os.path.basename(myimage)+".blur.tif"
# TiffIO(nameout,"w+").writeImage( blurred_ring[middle] )
picchi = np.array(picchi_again)
print( "IMAGE ", myimage,iim, Giplan, picchi, picchiI)
if len(picchi):
picchi=picchi[:,1:]
picchiL=picchi.T.tolist()
# if Filter is not None:
# good = np.where( Filter[picchiL]>0 )
# print " good " , good
# if len(good[0]):
# res[ myimage]=((picchi[good].T*binning).tolist() ,Giplan , picchiI , picchiPmap)
# else:
# res[ myimage]=((picchi.T*binning).tolist() ,Giplan , picchiI, picchiPmap )
tok = {}
if flist_alt is not None:
for pos , II, pmap, pmap_alt in zip((picchi*binning).tolist() , picchiI, picchiPmap, picchiPmap_alt ):
tok[tuple(pos)] = { "pos":pos, "Gi":Giplan, "II":II, "pmap":[pmap, pmap_alt] , "accepted":True, "mask":None}
else:
for pos , II, pmap in zip((picchi*binning).tolist() , picchiI, picchiPmap ):
tok[tuple(pos)] = { "pos":pos, "Gi":Giplan, "II":II, "pmap":[pmap, None] , "accepted":True, "mask":None}
res[ myimage]=tok
# for i in range(1,nprocs):
# if myrank==0:
# edict = comm.recv( source=i, tag=777+i)
# res.update(edict)
# else:
# if i==myrank:
# comm.send( res ,dest=0, tag=777+i)
comm.Barrier()
for i in range(1,nprocs):
if myrank==0:
edict = comm.recv( source=i, tag=777+i)
res.update(edict)
else:
if i==myrank:
comm.send( res ,dest=0, tag=777+i)
if myrank==0:
if globalize[0]>0 and globalize[1]>0:
bin = globalize[0]
dist = globalize[1]
globVol = np.zeros( [len(flist)/bin, dim1/bin , dim2/bin ],"f")
count=0
for key, (poss, nseq , intes ) in res.items():
if True or (nseq<29+10 and nseq>29-10):
globVol[nseq/bin][(np.array(poss,"i")/bin).tolist()] = intes
count+=len(intes)
# globVol[nseq/bin][poss.tolist()] = intes
print( " GLOBALIZING ", count , " PEAKS ")
# picchi = peak_local_max( globVol , min_distance = 2 , exclude_border = True, threshold_abs=1.0e-30)
picchi = mylocalmaxs( globVol , min_dist = dist , nobord = True, tabs=1.0e-30)
print( " GLOBALIZING shrinked to ", picchi.shape[0] , " PEAKS ")
picchi = picchi.T
zpos = set(picchi[0] )
newdic = {}
for z in zpos:
ys = picchi[1][ picchi[0]==z]
xs = picchi[2][ picchi[0]==z]
intes = globVol[z][ [ ys.tolist(),xs.tolist()] ]
newdic[ flist[z] ]= ( [ ys*bin,xs*bin], z*bin, intes )
res = newdic
maxima_file = open(file_maxima, 'w')
pickle.dump(res, maxima_file)
pickle.dump((DIM1orig ,DIM2orig), maxima_file)
maxima_file.close()
def main(params):
Filter_fname = params.filter_file
flist = glob.glob(params.images_prefix+"*"+params.images_postfix )
flist.sort()
if params.images_prefix_alt is not None:
flist_alt = glob.glob(params.images_prefix_alt+"*"+params.images_postfix )
flist_alt.sort()
else:
flist_alt = None
if params.drawmask == True:
mainDrawmask( flist )
elif params.review == False:
mainRead( params.filter_file,flist,flist_alt, params.threshold_abs, params.threshold_rel, params.bin ,
params.peak_size, params.spot_view_size, params.medianize, params.globalize, params.file_maxima )
print( " processo ", myrank, " aspetta ")
comm.Barrier()
else:
assert(nprocs==1)
timer.start(100)
mainReview( params , flist, flist_alt)
USAGE = """ USAGE :
extract_maxima inputfile
example input file containing
images_prefix = "Diffuse1/bronze7_0001p_"
images_postfix = ".cbf"
threshold_abs=10
threshold_rel=0.05
peak_size = 10
# OPTIONAL or set it to None for no mask filtering of bad pixels
filter_file = "newMask.edf"
medianize = 0
globalize = (0,0)
file_maxima = "maxima.p"
review = False
review_size = 10
"""
print( __name__ )
if(sys.argv[0][-12:]!="sphinx-build"):
if len(sys.argv) < 2:
print( USAGE)
else:
class Parameters:
filter_file = None
images_prefix = None
images_prefix_alt = None
volume_reconstruction = True
review = False
file_maxima="massimi.p"
review_size = 10
bin = 1
exec(open(sys.argv[1]).read())
pp = Parameters()
main(pp )
print( "exiting")
sys.exit(0)
# if __name__ == "__main__":