Changeset 3154


Ignore:
Timestamp:
Nov 18, 2017 6:14:55 PM (4 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:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIdataGUI.py

    r3151 r3154  
    30123012        '''Allows user to supply scale factor(s) when summing data -
    30133013        TODO: CAN WE PREVIEW RESULT HERE?'''
     3014        #TODO: also add filter box for selection of items
    30143015        def __init__(self,parent,title,text,dataType,data,dataList):
    30153016            wx.Dialog.__init__(self,parent,-1,title,size=(400,250),
  • trunk/GSASIIimage.py

    r3147 r3154  
    3333import GSASIIpwd as G2pwd
    3434import GSASIIspc as G2spc
     35import GSASIImath as G2mth
    3536
    3637# trig functions in degrees
     
    838839            tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
    839840                tay,len(polygon),polygon,tamp)[:nI*nJ]))
    840     try:
    841         import spotmask as sm
    842         spots = masks['Points'].T
    843         if len(spots):
    844             tam = ma.mask_or(tam,ma.getmask(sm.spotmask(nI*nJ,tax,
    845                     tay,len(spots),spots,tamp)[:nI*nJ]))
    846     except:
     841    if True:
     842#    try:
     843#        import spotmask as sm
     844#        spots = masks['Points'].T
     845#        if len(spots):
     846#            tam = ma.mask_or(tam,ma.getmask(sm.spotmask(nI*nJ,tax,
     847#                    tay,len(spots),spots,tamp)[:nI*nJ]))
     848#    except:
    847849        for X,Y,rsq in masks['Points'].T:
    848850            tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,rsq)))
     
    11381140    StrainPrint(ValSig,dset)
    11391141    return vals,sig,covMatrix
     1142
     1143def FitImageSpots(Image,ImMax,ind,pixSize,nxy):
     1144   
     1145    def calcMean(nxy,pixSize,img):
     1146        gdx,gdy = np.mgrid[0:nxy,0:nxy]
     1147        gdx = ma.array((gdx-nxy//2)*pixSize[0]/1000.,mask=~ma.getmaskarray(ImBox))
     1148        gdy = ma.array((gdy-nxy//2)*pixSize[1]/1000.,mask=~ma.getmaskarray(ImBox))
     1149        posx = ma.sum(gdx)/ma.count(gdx)
     1150        posy = ma.sum(gdy)/ma.count(gdy)
     1151        return posx,posy
     1152   
     1153    def calcPeak(values,nxy,pixSize,img):
     1154        back,mag,px,py,sig = values
     1155        peak = np.zeros([nxy,nxy])+back
     1156        nor = 1./(2.*np.pi*sig**2)
     1157        gdx,gdy = np.mgrid[0:nxy,0:nxy]
     1158        gdx = (gdx-nxy//2)*pixSize[0]/1000.
     1159        gdy = (gdy-nxy//2)*pixSize[1]/1000.
     1160        arg = (gdx-px)**2+(gdy-py)**2       
     1161        peak += mag*nor*np.exp(-arg/(2.*sig**2))
     1162        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))
     1163   
     1164    def calc2Peak(values,nxy,pixSize,img):
     1165        back,mag,px,py,sigx,sigy,rho = values
     1166        peak = np.zeros([nxy,nxy])+back
     1167        nor = 1./(2.*np.pi*sigx*sigy*np.sqrt(1.-rho**2))
     1168        gdx,gdy = np.mgrid[0:nxy,0:nxy]
     1169        gdx = (gdx-nxy//2)*pixSize[0]/1000.
     1170        gdy = (gdy-nxy//2)*pixSize[1]/1000.
     1171        argnor = -1./(2.*(1.-rho**2))
     1172        arg = (gdx-px)**2/sigx**2+(gdy-py)**2/sigy**2-2.*rho*(gdx-px)*(gdy-py)/(sigx*sigy)       
     1173        peak += mag*nor*np.exp(argnor*arg)
     1174        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))       
     1175   
     1176    nxy2 = nxy//2
     1177    ImBox = Image[ind[1]-nxy2:ind[1]+nxy2+1,ind[0]-nxy2:ind[0]+nxy2+1]
     1178    back = np.min(ImBox)
     1179    mag = np.sum(ImBox-back)
     1180    vals = [back,mag,0.,0.,0.2,0.2,0.]
     1181    ImBox = ma.array(ImBox,dtype=float,mask=ImBox>0.75*ImMax)
     1182    px = (ind[0]+.5)*pixSize[0]/1000.
     1183    py = (ind[1]+.5)*pixSize[1]/1000.
     1184    if ma.any(ma.getmaskarray(ImBox)):
     1185        vals = calcMean(nxy,pixSize,ImBox)
     1186        posx,posy = [px+vals[0],py+vals[1]]
     1187        return [posx,posy,6.]
     1188    else:
     1189        result = leastsq(calc2Peak,vals,args=(nxy,pixSize,ImBox),full_output=True)
     1190        vals = result[0]
     1191        ratio = vals[4]/vals[5]
     1192        if 0.5 < ratio < 2.0 and vals[2] < 2. and vals[3] < 2.:
     1193            posx,posy = [px+vals[2],py+vals[3]]
     1194            return [posx,posy,min(6.*vals[4],6.)]
     1195        else:
     1196            return None
    11401197   
    11411198def AutoSpotMasks(Image,Masks,Controls):
     1199   
    11421200    print ('auto spot search')
    1143     maxPks = 500
    1144     dmax = 1.5         #just beyond diamond 111 @ 2.0596
    1145     wave = Controls['wavelength']
    1146     tthMax = 2.0*asind(wave/(2.*dmax))
     1201    nxy = 15
     1202    rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
    11471203    pixelSize = Controls['pixelSize']
    1148     lap = np.array([[0,1,0],[1,-4,1],[0,1,0]],np.float32)
    1149     ck = scsg.cspline2d(Image.T,24.0)
    1150     spotMask1=-scsg.convolve2d(ck,lap,mode='same',boundary='symm')
    1151     for mul in range(10):
    1152         spotMask = ma.array(Image,mask=(spotMask1<(mul+1)*1.e4))
    1153         indx = np.transpose(spotMask.nonzero())
    1154         if not len(indx):       #smooth image - no sharp points
    1155             break
    1156         pts,distort = scvq.kmeans(np.array(indx,np.float32),maxPks,thresh=1.0)
    1157         peaks = pts*pixelSize/1000.
    1158         tths = GetTth(peaks.T[0],peaks.T[1],Controls)
    1159         masktths = ma.getmask(ma.array(tths,mask=tths>tthMax))
    1160         pts = pts[masktths,:]
    1161         A,B = np.mgrid[0:len(pts),0:len(pts)]   
    1162         M = np.sqrt(np.sum((pts[A]-pts[B])**2,axis=-1))
    1163         MX = ma.array(M,mask=(M<10.))
    1164         Indx = ma.nonzero(ma.getmask(MX))
    1165         print ('Spots found: %d %d %.2f %d'%(mul,len(pts),distort,len(Indx[0])))
    1166         if distort < .3:
    1167             break
    1168     if not len(Indx[0]):
    1169         print ('no auto search spots found')
    1170         return None
    1171     if distort > 2.:
    1172         print ('no big spots found')
    1173         return None
    1174     #use clustered points to get position & spread
    1175    
    1176     Indx = np.sort(Indx,axis=0).T
    1177     Nmax = np.max(Indx)
    1178     nums = [list(set([ind[1] for ind in Indx if ind[1]>=jnd and ind[0] == jnd])) for jnd in range(Nmax+1)]
    1179     delList = []
    1180     for Num in nums:
    1181         delList += list(np.sort(Num))[1:]
    1182     delList = list(set(delList))
    1183     Nums = [num for inum,num in enumerate(nums) if inum not in delList]
    1184     points = []
    1185     esds = []
    1186     for num in Nums:
    1187         points.append(np.mean(pts[num],axis=0))
    1188         esds.append(np.mean(np.var(pts[num],axis=0)))
    1189     peaks = (np.array(points)+np.array([.5,.5]))*pixelSize/1000.
    1190     Points = np.ones((peaks.shape[0],3))*2.
    1191     Points[:,2] += np.array(esds)*pixelSize[0]/1000.
    1192     Points[:,:2] = peaks
    1193     Masks['Points'] = list(Points)
     1204    spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
     1205    indices = (-1,0,1)
     1206    rolls = np.array([[ix,iy] for ix in indices for iy in indices])
     1207    time0 = time.time()
     1208    for roll in rolls:
     1209        if np.any(roll):        #avoid [0,0]
     1210            spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.),dtype=float)
     1211    mags = spotMask[spotMask.nonzero()]
     1212    indx = np.transpose(spotMask.nonzero())
     1213    size1 = mags.shape[0]
     1214    magind = [[indx[0][0],indx[0][1],mags[0]],]
     1215    for ind,mag in list(zip(indx,mags))[1:]:        #remove duplicates
     1216#            ind[0],ind[1],I,J = ImageLocalMax(Image,nxy,ind[0],ind[1])
     1217        if (magind[-1][0]-ind[0])**2+(magind[-1][1]-ind[1])**2 > 16:
     1218            magind.append([ind[0],ind[1],Image[ind[0],ind[1]]]) 
     1219    magind = np.array(magind).T
     1220    indx = np.array(magind[0:2],dtype=np.int32)
     1221    mags = magind[2]
     1222    size2 = mags.shape[0]
     1223    print ('Initial search done: %d -->%d %.2fs'%(size1,size2,time.time()-time0))
     1224    nx,ny = Image.shape
     1225    ImMax = np.max(Image)
     1226    peaks = []
     1227    nxy2 = nxy//2
     1228    mult = 0.001
     1229    num = 1e6
     1230    while num>500:
     1231        mult += .0001           
     1232        minM = mult*np.max(mags)
     1233        num = ma.count(ma.array(mags,mask=mags<=minM))
     1234        print('try',mult,minM,num)
     1235    minM = mult*np.max(mags)
     1236    print ('Find biggest spots:',mult,num,minM)
     1237    for i,mag in enumerate(mags):
     1238        if mag > minM:
     1239            if (nxy2 < indx[0][i] < nx-nxy2-1) and (nxy2 < indx[1][i] < ny-nxy2-1):
     1240#                    print ('try:%d %d %d %.2f'%(i,indx[0][i],indx[1][i],mags[i]))
     1241                peak = FitImageSpots(Image,ImMax,[indx[1][i],indx[0][i]],pixelSize,nxy)
     1242                if peak and not any(np.isnan(np.array(peak))):
     1243                    peaks.append(peak)
     1244#                    print (' Spot found: %s'%str(peak))
     1245    peaks = G2mth.sortArray(G2mth.sortArray(peaks,1),0)
     1246    Peaks = [peaks[0],]
     1247    for peak in peaks[1:]:
     1248        if GetDsp(peak[0],peak[1],Controls) >= 1.:      #toss non-diamond low angle spots
     1249            continue
     1250        if (peak[0]-Peaks[-1][0])**2+(peak[1]-Peaks[-1][1])**2 > peak[2]*Peaks[-1][2] :
     1251            Peaks.append(peak)
     1252#            print (' Spot found: %s'%str(peak))
     1253    print ('Spots found: %d time %.2fs'%(len(Peaks),time.time()-time0))
     1254    Masks['Points'] = Peaks
    11941255    return None
    1195 
    1196 #    rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
    1197 #    print 'auto spot search'
    1198 #    pixelSize = Controls['pixelSize']
    1199 #    spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
    1200 #    indices = (-1,0,1)
    1201 #    rolls = np.array([[ix,iy] for ix in indices for iy in indices])
    1202 #    for roll in rolls:
    1203 #        if np.any(roll):        #avoid [0,0]
    1204 #            spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.))
    1205 #    mags = spotMask[spotMask.nonzero()]
    1206 #    indx = np.transpose(spotMask.nonzero())
    1207 #    nx,ny = Image.shape
    1208 #    jndx = []
    1209 #    for [ind,mag] in zip(indx,mags):
    1210 #        if (1 < ind[0] < nx-2) and (1 < ind[1] < ny-2):
    1211 #            cent = np.zeros((3,3))
    1212 #            cent[1,1] = mag
    1213 #            msk = np.array(Image[ind[0]-2:ind[0]+3,ind[1]-2:ind[1]+3])
    1214 #            msk = ma.array(msk,mask=(msk==mag))
    1215 #            if mag > 1.5*ma.mean(msk):
    1216 #                jndx.append([ind[1]+.5,ind[0]+.5])
    1217 #    print 'Spots found: ',len(jndx)
    1218 #    jndx = np.array(jndx)
    1219 #    peaks = jndx*pixelSize/1000.
    1220 #    tth = GetTth(peaks.T[0],peaks.T[1],Controls)
    1221 #    histth,bins = np.histogram(tth,2500)
    1222 #    for shft in [-1,1]:
    1223 #        histmsk = ma.array(histth,mask=-(histth-np.roll(histth,shft)<2.))
    1224 #    histmsk = ma.array(histmsk,mask=(histmsk>5))
    1225 #    binmsk = np.array(ma.nonzero(histmsk))
    1226 #    digts = np.digitize(tth,bins)-1
    1227 #    PeaksList = []
    1228 #    for digt,peak in zip(digts,peaks):
    1229 #        if digt in binmsk:
    1230 #            PeaksList.append(peak)
    1231 #    #should be able to filter out spotty Bragg rings here
    1232 #    PeaksList = np.array(PeaksList)
    1233 #    Points = np.ones((PeaksList.shape[0],3))
    1234 #    Points[:,:2] = PeaksList
    1235 #    Masks['Points'] = list(Points)
    1236 #    return None
    1237        
    1238                    
    1239    
    1240        
  • trunk/GSASIIimgGUI.py

    r3142 r3154  
    234234        '''
    235235        CleanupMasks(masks)
    236         blkSize = 128   #this seems to be optimal; will break in polymask if >1024
     236        blkSize = 1024   #this seems to be optimal; will break in polymask if >1024
     237#        blkSize = 128   #this seems to be optimal; will break in polymask if >1024
    237238        sumImg = GetImageZ(G2frame,data)
    238239        wx.BeginBusyCursor()
     
    15141515        #   (Imin0, Imax0) => Range[0] = data['range'][0] # lowest to highest pixel intensity
    15151516        #   [Imin, Imax] => Range[1] = data['range'][1] #   lowest to highest pixel intensity on cmap scale
     1517        scaleChoices = ("100%","99%","95%","90%","80%","?")
     1518        scaleSel = wx.Choice(G2frame.dataWindow,choices=scaleChoices,size=(-1,-1))
     1519        if (Range[1][0] == Range[0][0] and Range[1][1] == Range[0][1]):
     1520            scaleSel.SetSelection(0)
     1521        else:
     1522            scaleSel.SetSelection(len(scaleChoices)-1)
     1523        scaleSel.Bind(wx.EVT_CHOICE,OnAutoSet)
     1524       
    15161525        maxSizer = wx.GridBagSizer(0,0)
    15171526        r = c = 0
     
    15331542        maxSizer.Add(maxVal,(r,c))
    15341543        c += 1
    1535         scaleChoices = ("100%","99%","95%","90%","80%","?")
    1536         scaleSel = wx.Choice(G2frame.dataWindow,choices=scaleChoices,size=(-1,-1))
    1537         if (Range[1][0] == Range[0][0] and
    1538             Range[1][1] == Range[0][1]):
    1539             scaleSel.SetSelection(0)
    1540         else:
    1541             scaleSel.SetSelection(len(scaleChoices)-1)
    1542         scaleSel.Bind(wx.EVT_CHOICE,OnAutoSet)
    15431544        maxSizer.Add(scaleSel,(r,c),(2,1),flag=wx.ALIGN_CENTER)
    15441545        c = 0
     
    16241625    mainSizer.Add(littleSizer,0,)
    16251626    if len(Spots):
    1626         lbl = wx.StaticText(parent=G2frame.dataWindow,label=' Spot masks')
     1627        lbl = wx.StaticText(parent=G2frame.dataWindow,label=' Spot masks(on plot LB drag to move, shift-LB drag to resize, RB to delete)')
    16271628        lbl.SetBackgroundColour(wx.Colour(200,200,210))
    16281629        mainSizer.Add(lbl,0,wx.EXPAND|wx.ALIGN_CENTER,0)
  • trunk/GSASIIpath.py

    r3136 r3154  
    629629    This routine is only used when debug=True is set in config.py   
    630630    '''
    631     #from IPython.core import ultratb
     631    from IPython.core import ultratb
    632632    if 'win' in sys.platform:
    633633        ultratb.FormattedTB(call_pdb=False,color_scheme='NoColor')(*args)
  • trunk/GSASIIplot.py

    r3153 r3154  
    46134613    if G2frame.MaskKey == 's':
    46144614        Page.figure.suptitle('Multiple spot mode on (size={}), press s or right-click to end'
    4615                              .format(G2frame.spotSize),color='r',fontweight='bold')
     4615            .format(G2frame.spotSize),color='r',fontweight='bold')
    46164616    else:
    46174617        Page.figure.suptitle('New spot size={}'.format(G2frame.spotSize),
     
    46664666    from matplotlib.patches import Ellipse,Circle
    46674667    import numpy.ma as ma
     4668    G2frame.ShiftDown = False
    46684669    G2frame.cid = None
    46694670    #Dsp = lambda tth,wave: wave/(2.*npsind(tth/2.))
     
    47294730                     G2frame.G2plotNB.status.SetStatusText( \
    47304731                        'Radius=%.3fmm, 2-th=%.3fdeg, dsp=%.3fA, Q=%.5fA-1, azm=%.2fdeg, I=%6d'%(radius,tth,dsp,Q,azm,Int),1)
     4732
     4733    def OnImPlotKeyRelease(event):
     4734        if event.key == 'shift':
     4735            G2frame.ShiftDown = False
    47314736
    47324737    def OnImPlotKeyPress(event):
     
    47844789                G2frame.MskDelete = True
    47854790                OnStartMask(G2frame)
     4791            elif event.key == 'shift':
     4792                G2frame.ShiftDown = True
    47864793               
    47874794        elif treeItem == 'Stress/Strain':
     
    49084915                itemNum = G2frame.itemPicked.itemNumber
    49094916                if event.button == 1:
    4910                     x = Masks['Points'][itemNum][0]+Xpos-XposBeforeDrag
    4911                     y = Masks['Points'][itemNum][1]+Ypos-YposBeforeDrag
    4912                     pick.center=[x,y]
    4913                 elif event.button == 3:
    4914                     r = math.sqrt((Xpos-Masks['Points'][itemNum][0])**2+
    4915                               (Ypos-Masks['Points'][itemNum][1])**2)
    4916                     pick.radius = r
     4917                    if G2frame.ShiftDown:
     4918                        r = math.sqrt((Xpos-Masks['Points'][itemNum][0])**2+
     4919                                  (Ypos-Masks['Points'][itemNum][1])**2)
     4920                        pick.radius = r
     4921                    else:
     4922                        x = Masks['Points'][itemNum][0]+Xpos-XposBeforeDrag
     4923                        y = Masks['Points'][itemNum][1]+Ypos-YposBeforeDrag
     4924                        pick.center=[x,y]
    49174925                Page.figure.gca().draw_artist(pick)
    49184926            elif pickType.startswith('Ring'):
     
    51325140                    ToggleMultiSpotMask(G2frame)
    51335141                else:
    5134                     #TODO - try a fit of spot to image?
    5135                     spot = [Xpos,Ypos,G2frame.spotSize]
     5142                    if G2frame.ShiftDown:   #force user selection
     5143                        sig = G2frame.spotSize
     5144                    else:                   #optimize spot pick
     5145                        pixLimit = 5
     5146                        Xpix,Ypix = Xpos*scalex,Ypos*scaley
     5147                        Xpix,Ypix,I,J = G2img.ImageLocalMax(G2frame.ImageZ,pixLimit,Xpix,Ypix)
     5148                        ind = [int(Xpix),int(Ypix)]
     5149                        nxy = 15
     5150                        ImMax = np.max(G2frame.ImageZ)
     5151                        result = G2img.FitImageSpots(G2frame.ImageZ,ImMax,ind,pixelSize,nxy)
     5152                        if result:
     5153                            Xpos,Ypos,sig = result
     5154                        else:
     5155                            print ('Not a spot')
     5156                            return
     5157                    spot = [Xpos,Ypos,sig]
    51365158                    Masks['Points'].append(spot)
    51375159                    artist = Circle((Xpos,Ypos),radius=spot[2]/2,fc='none',ec='r',picker=3)
     
    51425164                    Page.canvas.draw()
    51435165                return
    5144             elif G2frame.MaskKey == 's':
    5145                 if event.button == 1:
    5146                     spot = [Xpos,Ypos,G2frame.spotSize]
    5147                     Masks['Points'].append(spot)
    5148                     G2imG.UpdateMasks(G2frame,Masks)
    5149                 G2frame.MaskKey = ''
    5150                 wx.CallAfter(PlotImage,G2frame,newImage=True)
    5151                 return
    51525166            elif G2frame.MaskKey == 'r':
    51535167                if event.button == 1:
     
    52745288                G2frame.Razim.SetValue(Data['LRazimuth'][1])
    52755289            elif pickType == "Spot" and treeItem == 'Masks':
     5290                spotnum = G2frame.itemPicked.itemNumber
    52765291                if event.button == 1:
    5277                     spotnum = G2frame.itemPicked.itemNumber
    5278                     Masks['Points'][spotnum][0:2] = G2frame.itemPicked.center
     5292                    if G2frame.ShiftDown:
     5293                        Masks['Points'][spotnum] = list(G2frame.itemPicked.center) + [2.*G2frame.itemPicked.radius,]
     5294                    else:
     5295                    # update the selected circle mask with the last drawn values
     5296                        Masks['Points'][spotnum][0:2] = G2frame.itemPicked.center
    52795297                elif event.button == 3:
    5280                     # update the selected circle mask with the last drawn values
    5281                     spotnum = G2frame.itemPicked.itemNumber
    5282                     Masks['Points'][spotnum] = list(G2frame.itemPicked.center) + [
    5283                         2.*G2frame.itemPicked.radius]
     5298                    del Masks['Points'][spotnum]
    52845299                G2imG.UpdateMasks(G2frame,Masks)
    52855300            elif pickType.startswith('Ring') and treeItem == 'Masks':
     
    53125327    else:
    53135328        Page.canvas.mpl_connect('key_press_event', OnImPlotKeyPress)
     5329        Page.canvas.mpl_connect('key_release_event', OnImPlotKeyRelease)
    53145330        Page.canvas.mpl_connect('motion_notify_event', OnImMotion)
    53155331        Page.canvas.mpl_connect('pick_event', OnImPick)
Note: See TracChangeset for help on using the changeset viewer.