Changeset 3154 for trunk/GSASIIimage.py


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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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        
Note: See TracChangeset for help on using the changeset viewer.