Changeset 3154

Timestamp:

Nov 18, 2017 6:14:55 PM (6 years ago)

Author:

vondreele

Message:

new & imporved spot mask routines:
manual spot mask now fits position & size for selected spot.
changed mouse controls LB drag moves spot, shift-LB changes size & RB deletes spot
autospotmask improved - new algorithm

Location:

trunk

Files:

• GSASIIdataGUI.py (modified) (1 diff)
• GSASIIimage.py (modified) (3 diffs)
• GSASIIimgGUI.py (modified) (4 diffs)
• GSASIIpath.py (modified) (1 diff)
• GSASIIplot.py (modified) (9 diffs)

trunk/GSASIIdataGUI.py

@@ -3012 +3012 @@
         '''Allows user to supply scale factor(s) when summing data - 
         TODO: CAN WE PREVIEW RESULT HERE?'''
+        #TODO: also add filter box for selection of items
         def __init__(self,parent,title,text,dataType,data,dataList):
             wx.Dialog.__init__(self,parent,-1,title,size=(400,250),

trunk/GSASIIimage.py

@@ -33 +33 @@
 import GSASIIpwd as G2pwd
 import GSASIIspc as G2spc
+import GSASIImath as G2mth
 
 # trig functions in degrees
@@ -838 +839 @@
             tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
                 tay,len(polygon),polygon,tamp)[:nI*nJ]))
-    try:
-        import spotmask as sm
-        spots = masks['Points'].T
-        if len(spots):
-            tam = ma.mask_or(tam,ma.getmask(sm.spotmask(nI*nJ,tax,
-                    tay,len(spots),spots,tamp)[:nI*nJ]))
-    except:
+    if True:
+#    try:
+#        import spotmask as sm
+#        spots = masks['Points'].T
+#        if len(spots):
+#            tam = ma.mask_or(tam,ma.getmask(sm.spotmask(nI*nJ,tax,
+#                    tay,len(spots),spots,tamp)[:nI*nJ]))
+#    except:
         for X,Y,rsq in masks['Points'].T:
             tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,rsq)))
@@ -1138 +1140 @@
     StrainPrint(ValSig,dset)
     return vals,sig,covMatrix
+
+def FitImageSpots(Image,ImMax,ind,pixSize,nxy):
+    
+    def calcMean(nxy,pixSize,img):
+        gdx,gdy = np.mgrid[0:nxy,0:nxy]
+        gdx = ma.array((gdx-nxy//2)*pixSize[0]/1000.,mask=~ma.getmaskarray(ImBox))
+        gdy = ma.array((gdy-nxy//2)*pixSize[1]/1000.,mask=~ma.getmaskarray(ImBox))
+        posx = ma.sum(gdx)/ma.count(gdx)
+        posy = ma.sum(gdy)/ma.count(gdy)
+        return posx,posy
+    
+    def calcPeak(values,nxy,pixSize,img):
+        back,mag,px,py,sig = values
+        peak = np.zeros([nxy,nxy])+back
+        nor = 1./(2.*np.pi*sig**2)
+        gdx,gdy = np.mgrid[0:nxy,0:nxy]
+        gdx = (gdx-nxy//2)*pixSize[0]/1000.
+        gdy = (gdy-nxy//2)*pixSize[1]/1000.
+        arg = (gdx-px)**2+(gdy-py)**2        
+        peak += mag*nor*np.exp(-arg/(2.*sig**2))
+        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))
+    
+    def calc2Peak(values,nxy,pixSize,img):
+        back,mag,px,py,sigx,sigy,rho = values
+        peak = np.zeros([nxy,nxy])+back
+        nor = 1./(2.*np.pi*sigx*sigy*np.sqrt(1.-rho**2))
+        gdx,gdy = np.mgrid[0:nxy,0:nxy]
+        gdx = (gdx-nxy//2)*pixSize[0]/1000.
+        gdy = (gdy-nxy//2)*pixSize[1]/1000.
+        argnor = -1./(2.*(1.-rho**2))
+        arg = (gdx-px)**2/sigx**2+(gdy-py)**2/sigy**2-2.*rho*(gdx-px)*(gdy-py)/(sigx*sigy)        
+        peak += mag*nor*np.exp(argnor*arg)
+        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))        
+    
+    nxy2 = nxy//2
+    ImBox = Image[ind[1]-nxy2:ind[1]+nxy2+1,ind[0]-nxy2:ind[0]+nxy2+1]
+    back = np.min(ImBox)
+    mag = np.sum(ImBox-back)
+    vals = [back,mag,0.,0.,0.2,0.2,0.]
+    ImBox = ma.array(ImBox,dtype=float,mask=ImBox>0.75*ImMax)
+    px = (ind[0]+.5)*pixSize[0]/1000.
+    py = (ind[1]+.5)*pixSize[1]/1000.
+    if ma.any(ma.getmaskarray(ImBox)):
+        vals = calcMean(nxy,pixSize,ImBox)
+        posx,posy = [px+vals[0],py+vals[1]]
+        return [posx,posy,6.]
+    else:
+        result = leastsq(calc2Peak,vals,args=(nxy,pixSize,ImBox),full_output=True)
+        vals = result[0]
+        ratio = vals[4]/vals[5]
+        if 0.5 < ratio < 2.0 and vals[2] < 2. and vals[3] < 2.:
+            posx,posy = [px+vals[2],py+vals[3]]
+            return [posx,posy,min(6.*vals[4],6.)]
+        else:
+            return None
     
 def AutoSpotMasks(Image,Masks,Controls):
+    
     print ('auto spot search')
-    maxPks = 500
-    dmax = 1.5         #just beyond diamond 111 @ 2.0596
-    wave = Controls['wavelength']
-    tthMax = 2.0*asind(wave/(2.*dmax))
+    nxy = 15
+    rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
     pixelSize = Controls['pixelSize']
-    lap = np.array([[0,1,0],[1,-4,1],[0,1,0]],np.float32)
-    ck = scsg.cspline2d(Image.T,24.0)
-    spotMask1=-scsg.convolve2d(ck,lap,mode='same',boundary='symm')
-    for mul in range(10):
-        spotMask = ma.array(Image,mask=(spotMask1<(mul+1)*1.e4))
-        indx = np.transpose(spotMask.nonzero())
-        if not len(indx):       #smooth image - no sharp points
-            break
-        pts,distort = scvq.kmeans(np.array(indx,np.float32),maxPks,thresh=1.0)
-        peaks = pts*pixelSize/1000.
-        tths = GetTth(peaks.T[0],peaks.T[1],Controls)
-        masktths = ma.getmask(ma.array(tths,mask=tths>tthMax))
-        pts = pts[masktths,:]
-        A,B = np.mgrid[0:len(pts),0:len(pts)]   
-        M = np.sqrt(np.sum((pts[A]-pts[B])**2,axis=-1))
-        MX = ma.array(M,mask=(M<10.))
-        Indx = ma.nonzero(ma.getmask(MX))
-        print ('Spots found: %d %d %.2f %d'%(mul,len(pts),distort,len(Indx[0])))
-        if distort < .3:
-            break
-    if not len(Indx[0]):
-        print ('no auto search spots found')
-        return None
-    if distort > 2.:
-        print ('no big spots found')
-        return None
-    #use clustered points to get position & spread
-    
-    Indx = np.sort(Indx,axis=0).T
-    Nmax = np.max(Indx)
-    nums = [list(set([ind[1] for ind in Indx if ind[1]>=jnd and ind[0] == jnd])) for jnd in range(Nmax+1)]
-    delList = []
-    for Num in nums:
-        delList += list(np.sort(Num))[1:]
-    delList = list(set(delList))
-    Nums = [num for inum,num in enumerate(nums) if inum not in delList]
-    points = []
-    esds = []
-    for num in Nums:
-        points.append(np.mean(pts[num],axis=0))
-        esds.append(np.mean(np.var(pts[num],axis=0)))
-    peaks = (np.array(points)+np.array([.5,.5]))*pixelSize/1000.
-    Points = np.ones((peaks.shape[0],3))*2.
-    Points[:,2] += np.array(esds)*pixelSize[0]/1000.
-    Points[:,:2] = peaks
-    Masks['Points'] = list(Points)
+    spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
+    indices = (-1,0,1)
+    rolls = np.array([[ix,iy] for ix in indices for iy in indices])
+    time0 = time.time()
+    for roll in rolls:
+        if np.any(roll):        #avoid [0,0]
+            spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.),dtype=float)
+    mags = spotMask[spotMask.nonzero()]
+    indx = np.transpose(spotMask.nonzero())
+    size1 = mags.shape[0]
+    magind = [[indx[0][0],indx[0][1],mags[0]],]
+    for ind,mag in list(zip(indx,mags))[1:]:        #remove duplicates
+#            ind[0],ind[1],I,J = ImageLocalMax(Image,nxy,ind[0],ind[1])
+        if (magind[-1][0]-ind[0])**2+(magind[-1][1]-ind[1])**2 > 16:
+            magind.append([ind[0],ind[1],Image[ind[0],ind[1]]])  
+    magind = np.array(magind).T
+    indx = np.array(magind[0:2],dtype=np.int32)
+    mags = magind[2]
+    size2 = mags.shape[0]
+    print ('Initial search done: %d -->%d %.2fs'%(size1,size2,time.time()-time0))
+    nx,ny = Image.shape
+    ImMax = np.max(Image)
+    peaks = []
+    nxy2 = nxy//2
+    mult = 0.001
+    num = 1e6
+    while num>500:
+        mult += .0001            
+        minM = mult*np.max(mags)
+        num = ma.count(ma.array(mags,mask=mags<=minM))
+        print('try',mult,minM,num)
+    minM = mult*np.max(mags)
+    print ('Find biggest spots:',mult,num,minM)
+    for i,mag in enumerate(mags):
+        if mag > minM:
+            if (nxy2 < indx[0][i] < nx-nxy2-1) and (nxy2 < indx[1][i] < ny-nxy2-1):
+#                    print ('try:%d %d %d %.2f'%(i,indx[0][i],indx[1][i],mags[i]))
+                peak = FitImageSpots(Image,ImMax,[indx[1][i],indx[0][i]],pixelSize,nxy)
+                if peak and not any(np.isnan(np.array(peak))):
+                    peaks.append(peak)
+#                    print (' Spot found: %s'%str(peak))
+    peaks = G2mth.sortArray(G2mth.sortArray(peaks,1),0)
+    Peaks = [peaks[0],]
+    for peak in peaks[1:]:
+        if GetDsp(peak[0],peak[1],Controls) >= 1.:      #toss non-diamond low angle spots
+            continue
+        if (peak[0]-Peaks[-1][0])**2+(peak[1]-Peaks[-1][1])**2 > peak[2]*Peaks[-1][2] :
+            Peaks.append(peak)
+#            print (' Spot found: %s'%str(peak))
+    print ('Spots found: %d time %.2fs'%(len(Peaks),time.time()-time0))
+    Masks['Points'] = Peaks
     return None
-
-#    rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
-#    print 'auto spot search'
-#    pixelSize = Controls['pixelSize']
-#    spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
-#    indices = (-1,0,1)
-#    rolls = np.array([[ix,iy] for ix in indices for iy in indices])
-#    for roll in rolls:
-#        if np.any(roll):        #avoid [0,0]
-#            spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.))
-#    mags = spotMask[spotMask.nonzero()]
-#    indx = np.transpose(spotMask.nonzero())
-#    nx,ny = Image.shape
-#    jndx = []
-#    for [ind,mag] in zip(indx,mags):
-#        if (1 < ind[0] < nx-2) and (1 < ind[1] < ny-2):
-#            cent = np.zeros((3,3))
-#            cent[1,1] = mag
-#            msk = np.array(Image[ind[0]-2:ind[0]+3,ind[1]-2:ind[1]+3])
-#            msk = ma.array(msk,mask=(msk==mag))
-#            if mag > 1.5*ma.mean(msk):
-#                jndx.append([ind[1]+.5,ind[0]+.5])
-#    print 'Spots found: ',len(jndx)
-#    jndx = np.array(jndx)
-#    peaks = jndx*pixelSize/1000.
-#    tth = GetTth(peaks.T[0],peaks.T[1],Controls)
-#    histth,bins = np.histogram(tth,2500)
-#    for shft in [-1,1]:
-#        histmsk = ma.array(histth,mask=-(histth-np.roll(histth,shft)<2.))
-#    histmsk = ma.array(histmsk,mask=(histmsk>5))
-#    binmsk = np.array(ma.nonzero(histmsk))
-#    digts = np.digitize(tth,bins)-1
-#    PeaksList = []
-#    for digt,peak in zip(digts,peaks):
-#        if digt in binmsk:
-#            PeaksList.append(peak)
-#    #should be able to filter out spotty Bragg rings here
-#    PeaksList = np.array(PeaksList)
-#    Points = np.ones((PeaksList.shape[0],3))
-#    Points[:,:2] = PeaksList
-#    Masks['Points'] = list(Points)
-#    return None
-        
-                    
-    
-        

trunk/GSASIIimgGUI.py

@@ -234 +234 @@
         '''
         CleanupMasks(masks)
-        blkSize = 128   #this seems to be optimal; will break in polymask if >1024
+        blkSize = 1024   #this seems to be optimal; will break in polymask if >1024
+#        blkSize = 128   #this seems to be optimal; will break in polymask if >1024
         sumImg = GetImageZ(G2frame,data)
         wx.BeginBusyCursor()
@@ -1514 +1515 @@
         #   (Imin0, Imax0) => Range[0] = data['range'][0] # lowest to highest pixel intensity
         #   [Imin, Imax] => Range[1] = data['range'][1] #   lowest to highest pixel intensity on cmap scale
+        scaleChoices = ("100%","99%","95%","90%","80%","?")
+        scaleSel = wx.Choice(G2frame.dataWindow,choices=scaleChoices,size=(-1,-1))
+        if (Range[1][0] == Range[0][0] and Range[1][1] == Range[0][1]):
+            scaleSel.SetSelection(0)
+        else:
+            scaleSel.SetSelection(len(scaleChoices)-1)
+        scaleSel.Bind(wx.EVT_CHOICE,OnAutoSet)
+        
         maxSizer = wx.GridBagSizer(0,0)
         r = c = 0
@@ -1533 +1542 @@
         maxSizer.Add(maxVal,(r,c))
         c += 1
-        scaleChoices = ("100%","99%","95%","90%","80%","?")
-        scaleSel = wx.Choice(G2frame.dataWindow,choices=scaleChoices,size=(-1,-1))
-        if (Range[1][0] == Range[0][0] and
-            Range[1][1] == Range[0][1]):
-            scaleSel.SetSelection(0)
-        else:
-            scaleSel.SetSelection(len(scaleChoices)-1)
-        scaleSel.Bind(wx.EVT_CHOICE,OnAutoSet)
         maxSizer.Add(scaleSel,(r,c),(2,1),flag=wx.ALIGN_CENTER)
         c = 0
@@ -1624 +1625 @@
     mainSizer.Add(littleSizer,0,)
     if len(Spots):
-        lbl = wx.StaticText(parent=G2frame.dataWindow,label=' Spot masks')
+        lbl = wx.StaticText(parent=G2frame.dataWindow,label=' Spot masks(on plot LB drag to move, shift-LB drag to resize, RB to delete)')
         lbl.SetBackgroundColour(wx.Colour(200,200,210))
         mainSizer.Add(lbl,0,wx.EXPAND|wx.ALIGN_CENTER,0)

trunk/GSASIIpath.py

@@ -629 +629 @@
     This routine is only used when debug=True is set in config.py    
     '''
-    #from IPython.core import ultratb
+    from IPython.core import ultratb
     if 'win' in sys.platform:
         ultratb.FormattedTB(call_pdb=False,color_scheme='NoColor')(*args)

trunk/GSASIIplot.py

@@ -4613 +4613 @@
     if G2frame.MaskKey == 's':
         Page.figure.suptitle('Multiple spot mode on (size={}), press s or right-click to end'
-                             .format(G2frame.spotSize),color='r',fontweight='bold')
+            .format(G2frame.spotSize),color='r',fontweight='bold')
     else:
         Page.figure.suptitle('New spot size={}'.format(G2frame.spotSize),
@@ -4666 +4666 @@
     from matplotlib.patches import Ellipse,Circle
     import numpy.ma as ma
+    G2frame.ShiftDown = False
     G2frame.cid = None
     #Dsp = lambda tth,wave: wave/(2.*npsind(tth/2.))
@@ -4729 +4730 @@
                      G2frame.G2plotNB.status.SetStatusText( \
                         'Radius=%.3fmm, 2-th=%.3fdeg, dsp=%.3fA, Q=%.5fA-1, azm=%.2fdeg, I=%6d'%(radius,tth,dsp,Q,azm,Int),1)
+
+    def OnImPlotKeyRelease(event):
+        if event.key == 'shift':
+            G2frame.ShiftDown = False
 
     def OnImPlotKeyPress(event):
@@ -4784 +4789 @@
                 G2frame.MskDelete = True
                 OnStartMask(G2frame)
+            elif event.key == 'shift':
+                G2frame.ShiftDown = True
                 
         elif treeItem == 'Stress/Strain':
@@ -4908 +4915 @@
                 itemNum = G2frame.itemPicked.itemNumber
                 if event.button == 1:
-                    x = Masks['Points'][itemNum][0]+Xpos-XposBeforeDrag
-                    y = Masks['Points'][itemNum][1]+Ypos-YposBeforeDrag
-                    pick.center=[x,y]
-                elif event.button == 3:
-                    r = math.sqrt((Xpos-Masks['Points'][itemNum][0])**2+
-                              (Ypos-Masks['Points'][itemNum][1])**2)
-                    pick.radius = r
+                    if G2frame.ShiftDown:
+                        r = math.sqrt((Xpos-Masks['Points'][itemNum][0])**2+
+                                  (Ypos-Masks['Points'][itemNum][1])**2)
+                        pick.radius = r
+                    else:
+                        x = Masks['Points'][itemNum][0]+Xpos-XposBeforeDrag
+                        y = Masks['Points'][itemNum][1]+Ypos-YposBeforeDrag
+                        pick.center=[x,y]
                 Page.figure.gca().draw_artist(pick)
             elif pickType.startswith('Ring'):
@@ -5132 +5140 @@
                     ToggleMultiSpotMask(G2frame)
                 else:
-                    #TODO - try a fit of spot to image?
-                    spot = [Xpos,Ypos,G2frame.spotSize]
+                    if G2frame.ShiftDown:   #force user selection
+                        sig = G2frame.spotSize
+                    else:                   #optimize spot pick
+                        pixLimit = 5
+                        Xpix,Ypix = Xpos*scalex,Ypos*scaley
+                        Xpix,Ypix,I,J = G2img.ImageLocalMax(G2frame.ImageZ,pixLimit,Xpix,Ypix)
+                        ind = [int(Xpix),int(Ypix)]
+                        nxy = 15
+                        ImMax = np.max(G2frame.ImageZ)
+                        result = G2img.FitImageSpots(G2frame.ImageZ,ImMax,ind,pixelSize,nxy)
+                        if result:
+                            Xpos,Ypos,sig = result
+                        else:
+                            print ('Not a spot')
+                            return
+                    spot = [Xpos,Ypos,sig]
                     Masks['Points'].append(spot)
                     artist = Circle((Xpos,Ypos),radius=spot[2]/2,fc='none',ec='r',picker=3)
@@ -5142 +5164 @@
                     Page.canvas.draw()
                 return 
-            elif G2frame.MaskKey == 's':
-                if event.button == 1:
-                    spot = [Xpos,Ypos,G2frame.spotSize]
-                    Masks['Points'].append(spot)
-                    G2imG.UpdateMasks(G2frame,Masks)
-                G2frame.MaskKey = ''
-                wx.CallAfter(PlotImage,G2frame,newImage=True)
-                return
             elif G2frame.MaskKey == 'r':
                 if event.button == 1:
@@ -5274 +5288 @@
                 G2frame.Razim.SetValue(Data['LRazimuth'][1])
             elif pickType == "Spot" and treeItem == 'Masks':
+                spotnum = G2frame.itemPicked.itemNumber
                 if event.button == 1:
-                    spotnum = G2frame.itemPicked.itemNumber
-                    Masks['Points'][spotnum][0:2] = G2frame.itemPicked.center
+                    if G2frame.ShiftDown:
+                        Masks['Points'][spotnum] = list(G2frame.itemPicked.center) + [2.*G2frame.itemPicked.radius,]
+                    else:
+                    # update the selected circle mask with the last drawn values
+                        Masks['Points'][spotnum][0:2] = G2frame.itemPicked.center
                 elif event.button == 3:
-                    # update the selected circle mask with the last drawn values
-                    spotnum = G2frame.itemPicked.itemNumber
-                    Masks['Points'][spotnum] = list(G2frame.itemPicked.center) + [
-                        2.*G2frame.itemPicked.radius]
+                    del Masks['Points'][spotnum]
                 G2imG.UpdateMasks(G2frame,Masks)
             elif pickType.startswith('Ring') and treeItem == 'Masks':
@@ -5312 +5327 @@
     else:
         Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
+        Page.canvas.mpl_connect('key_release_event', OnImPlotKeyRelease)
         Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
         Page.canvas.mpl_connect('pick_event', OnImPick)