Changeset 2760
- Timestamp:
- Mar 24, 2017 12:20:38 PM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIimage.py
r2731 r2760 23 23 import polymask as pm 24 24 from scipy.optimize import leastsq 25 import scipy.signal as scsg 26 import scipy.cluster.vq as scvq 25 27 import copy 26 28 import GSASIIpath … … 1116 1118 1117 1119 def AutoSpotMasks(Image,Masks,Controls): 1118 rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)1119 1120 print 'auto spot search' 1121 maxPks = 500 1122 dmax = 1.5 #just beyond diamond 111 @ 2.0596 1123 wave = Controls['wavelength'] 1124 tthMax = 2.0*asind(wave/(2.*dmax)) 1120 1125 pixelSize = Controls['pixelSize'] 1121 spotMask = ma.array(Image,mask=(Image<np.mean(Image))) 1122 indices = (-1,0,1) 1123 rolls = np.array([[ix,iy] for ix in indices for iy in indices]) 1124 for roll in rolls: 1125 if np.any(roll): #avoid [0,0] 1126 spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.)) 1127 mags = spotMask[spotMask.nonzero()] 1128 indx = np.transpose(spotMask.nonzero()) 1129 nx,ny = Image.shape 1130 jndx = [] 1131 for [ind,mag] in zip(indx,mags): 1132 if (1 < ind[0] < nx-2) and (1 < ind[1] < ny-2): 1133 cent = np.zeros((3,3)) 1134 cent[1,1] = mag 1135 msk = np.array(Image[ind[0]-2:ind[0]+3,ind[1]-2:ind[1]+3]) 1136 msk = ma.array(msk,mask=(msk==mag)) 1137 if mag > 1.5*ma.mean(msk): 1138 jndx.append([ind[1]+.5,ind[0]+.5]) 1139 print 'Spots found: ',len(jndx) 1140 jndx = np.array(jndx) 1141 peaks = jndx*pixelSize/1000. 1142 tth = GetTth(peaks.T[0],peaks.T[1],Controls) 1143 histth,bins = np.histogram(tth,2500) 1144 for shft in [-1,1]: 1145 histmsk = ma.array(histth,mask=-(histth-np.roll(histth,shft)<2.)) 1146 histmsk = ma.array(histmsk,mask=(histmsk>5)) 1147 binmsk = np.array(ma.nonzero(histmsk)) 1148 digts = np.digitize(tth,bins)-1 1149 PeaksList = [] 1150 for digt,peak in zip(digts,peaks): 1151 if digt in binmsk: 1152 PeaksList.append(peak) 1153 #should be able to filter out spotty Bragg rings here 1154 PeaksList = np.array(PeaksList) 1155 Points = np.ones((PeaksList.shape[0],3)) 1156 Points[:,:2] = PeaksList 1126 lap = np.array([[0,1,0],[1,-4,1],[0,1,0]],np.float32) 1127 ck = scsg.cspline2d(Image.T,24.0) 1128 spotMask1=-scsg.convolve2d(ck,lap,mode='same',boundary='symm') 1129 for mul in range(10): 1130 spotMask = ma.array(Image,mask=(spotMask1<(mul+1)*1.e4)) 1131 indx = np.transpose(spotMask.nonzero()) 1132 if not len(indx): #smooth image - no sharp points 1133 break 1134 pts,distort = scvq.kmeans(np.array(indx,np.float32),maxPks,thresh=1.0) 1135 peaks = pts*pixelSize/1000. 1136 tths = GetTth(peaks.T[0],peaks.T[1],Controls) 1137 masktths = ma.getmask(ma.array(tths,mask=tths>tthMax)) 1138 pts = pts[masktths,:] 1139 A,B = np.mgrid[0:len(pts),0:len(pts)] 1140 M = np.sqrt(np.sum((pts[A]-pts[B])**2,axis=-1)) 1141 MX = ma.array(M,mask=(M<10.)) 1142 Indx = ma.nonzero(ma.getmask(MX)) 1143 print 'Spots found: ',mul,len(pts),distort,len(Indx[0]) 1144 if distort < .3: 1145 break 1146 if not len(Indx[0]): 1147 print 'no auto search spots found' 1148 return None 1149 if distort > 2.: 1150 print 'no big spots found' 1151 return None 1152 #use clustered points to get position & spread 1153 1154 Indx = np.sort(Indx,axis=0).T 1155 Nmax = np.max(Indx) 1156 nums = [list(set([ind[1] for ind in Indx if ind[1]>=jnd and ind[0] == jnd])) for jnd in range(Nmax+1)] 1157 delList = [] 1158 for Num in nums: 1159 delList += list(np.sort(Num))[1:] 1160 delList = list(set(delList)) 1161 Nums = [num for inum,num in enumerate(nums) if inum not in delList] 1162 points = [] 1163 esds = [] 1164 for num in Nums: 1165 points.append(np.mean(pts[num],axis=0)) 1166 esds.append(np.mean(np.var(pts[num],axis=0))) 1167 peaks = (np.array(points)+np.array([.5,.5]))*pixelSize/1000. 1168 Points = np.ones((peaks.shape[0],3))*2. 1169 Points[:,2] += np.array(esds)*pixelSize[0]/1000. 1170 Points[:,:2] = peaks 1157 1171 Masks['Points'] = list(Points) 1158 1172 return None 1173 1174 # rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1) 1175 # print 'auto spot search' 1176 # pixelSize = Controls['pixelSize'] 1177 # spotMask = ma.array(Image,mask=(Image<np.mean(Image))) 1178 # indices = (-1,0,1) 1179 # rolls = np.array([[ix,iy] for ix in indices for iy in indices]) 1180 # for roll in rolls: 1181 # if np.any(roll): #avoid [0,0] 1182 # spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.)) 1183 # mags = spotMask[spotMask.nonzero()] 1184 # indx = np.transpose(spotMask.nonzero()) 1185 # nx,ny = Image.shape 1186 # jndx = [] 1187 # for [ind,mag] in zip(indx,mags): 1188 # if (1 < ind[0] < nx-2) and (1 < ind[1] < ny-2): 1189 # cent = np.zeros((3,3)) 1190 # cent[1,1] = mag 1191 # msk = np.array(Image[ind[0]-2:ind[0]+3,ind[1]-2:ind[1]+3]) 1192 # msk = ma.array(msk,mask=(msk==mag)) 1193 # if mag > 1.5*ma.mean(msk): 1194 # jndx.append([ind[1]+.5,ind[0]+.5]) 1195 # print 'Spots found: ',len(jndx) 1196 # jndx = np.array(jndx) 1197 # peaks = jndx*pixelSize/1000. 1198 # tth = GetTth(peaks.T[0],peaks.T[1],Controls) 1199 # histth,bins = np.histogram(tth,2500) 1200 # for shft in [-1,1]: 1201 # histmsk = ma.array(histth,mask=-(histth-np.roll(histth,shft)<2.)) 1202 # histmsk = ma.array(histmsk,mask=(histmsk>5)) 1203 # binmsk = np.array(ma.nonzero(histmsk)) 1204 # digts = np.digitize(tth,bins)-1 1205 # PeaksList = [] 1206 # for digt,peak in zip(digts,peaks): 1207 # if digt in binmsk: 1208 # PeaksList.append(peak) 1209 # #should be able to filter out spotty Bragg rings here 1210 # PeaksList = np.array(PeaksList) 1211 # Points = np.ones((PeaksList.shape[0],3)) 1212 # Points[:,:2] = PeaksList 1213 # Masks['Points'] = list(Points) 1214 # return None 1215 1159 1216 1160 1217
Note: See TracChangeset
for help on using the changeset viewer.