Changeset 32
- Timestamp:
- Feb 25, 2010 3:29:55 PM (14 years ago)
- Location:
- trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASII.py
r31 r32 415 415 Data['tilt'] = 0.0 416 416 Data['rotation'] = 0.0 417 Data['refine'] = [True,False,True,True,True]418 417 Data['showLines'] = False 419 418 Data['ring'] = [] 420 419 Data['rings'] = [] 420 Data['cutoff'] = 10 421 Data['pixLimit'] = 20 421 422 Data['ellipses'] = [] 422 423 Data['masks'] = [] … … 426 427 Data['outChannels'] = 2500 427 428 Data['fullIntegrate'] = False 429 Data['setRings'] = False 428 430 Data['setDefault'] = False 429 431 Data['range'] = [(Imin,Imax),[Imin,Imax]] … … 1155 1157 return 1156 1158 Ypos = int(event.ydata)*self.imScale 1157 if event.key == 'c' and Data['refine'][0]:1159 if event.key == 'c': 1158 1160 cent = Data['center'] = [Xpos*pixelSize[0]/1000.,Ypos*pixelSize[1]/1000.] #convert to mm 1159 self.centText.SetValue(("%8.3f,%8.3f" % (cent[0],cent[1])))1160 1161 elif event.key == 'm': 1161 1162 xpos = Xpos*pixelSize[0]/1000. … … 1314 1315 yring *= scaley 1315 1316 ax.text(xring,yring,'+',color='r',ha='center',va='center',picker=3) 1316 # for ring in Data['rings']: 1317 # for xring,yring in ring: 1318 # xring *= scalex 1319 # yring *= scaley 1320 # ax.text(xring,yring,'+',ha='center',va='center') 1317 if Data['setRings']: 1318 for ring in Data['rings']: 1319 for xring,yring in ring: 1320 xring *= scalex 1321 yring *= scaley 1322 ax.text(xring,yring,'+',ha='center',va='center') 1321 1323 for ellipse in Data['ellipses']: 1322 1324 cent,phi,[width,height] = ellipse 1323 ax.add_artist(Ellipse([cent[0]*scalex,cent[1]*scaley],2*width*scalex,2*height*scalex,phi, fc=None))1325 ax.add_artist(Ellipse([cent[0]*scalex,cent[1]*scaley],2*width*scalex,2*height*scalex,phi,ec='r',fc=None)) 1324 1326 self.Img.axes.set_xlim(xlim) 1325 1327 self.Img.axes.set_ylim(ylim) -
trunk/GSASIIcomp.py
r31 r32 1208 1208 cent = [(ell[1]*ell[4]-2.0*ell[2]*ell[3])/det, \ 1209 1209 (ell[1]*ell[3]-2.0*ell[0]*ell[4])/det] 1210 phi = 0.5*atand(ell[1]/(ell[0]-ell[2] ))1210 phi = 0.5*atand(ell[1]/(ell[0]-ell[2]+1e-32)) 1211 1211 1212 1212 a3 = 0.5*(ell[0]+ell[2]+(ell[0]-ell[2])/cosd(2.0*phi)) … … 1215 1215 if f3/a3 < 0 or f3/b3 < 0: 1216 1216 return 0 1217 if phi > 0.:1218 sr1 = math.sqrt(f3/a3)1219 sr2 = math.sqrt(f3/b3)1220 else:1221 sr1 = math.sqrt(f3/a3)1222 sr2 = math.sqrt(f3/b3)1223 phi *= -1.1217 sr1 = math.sqrt(f3/a3) 1218 sr2 = math.sqrt(f3/b3) 1219 if sr1 > sr2: 1220 sr1,sr2 = SwapXY(sr1,sr2) 1221 phi -= 90. 1222 if phi < -90.: 1223 phi += 180. 1224 1224 return cent,phi,[sr1,sr2] 1225 1225 … … 1245 1245 w2 = w*2 1246 1246 size = len(image) 1247 if (w < Xpos < size-w) and (w < Ypos < size-w) :1247 if (w < Xpos < size-w) and (w < Ypos < size-w) and image[Ypos,Xpos]: 1248 1248 Z = image[Ypos-w:Ypos+w,Xpos-w:Xpos+w] 1249 1249 Zmax = np.argmax(Z) … … 1255 1255 return 0,0,0,0 1256 1256 1257 def makeRing(ellipse, scalex,scaley,imScale,image):1258 cent,phi, semiradii = ellipse1257 def makeRing(ellipse,pix,reject,scalex,scaley,imScale,image): 1258 cent,phi,radii = ellipse 1259 1259 cphi = cosd(phi) 1260 1260 sphi = sind(phi) 1261 1261 ring = [] 1262 for a in range(0,360 ):1263 x = semiradii[0]*cosd(a)1264 y = semiradii[1]*sind(a)1262 for a in range(0,360,2): 1263 x = radii[0]*cosd(a) 1264 y = radii[1]*sind(a) 1265 1265 X = (cphi*x-sphi*y+cent[0])*scalex*imScale 1266 1266 Y = (sphi*x+cphi*y+cent[1])*scaley*imScale 1267 X,Y,I,J = ImageLocalMax(image, 20,X,Y)1268 if I and J :1267 X,Y,I,J = ImageLocalMax(image,pix,X,Y) 1268 if I and J and I/J > reject: 1269 1269 X /= scalex*imScale 1270 1270 Y /= scaley*imScale 1271 1271 ring.append([X,Y]) 1272 if len(ring) < 75: #want more than 1/3of a circle1272 if len(ring) < 45: #want more than 1/4 of a circle 1273 1273 return [] 1274 1274 return ring 1275 1276 def calcDist(radii,tth): 1277 stth = sind(tth) 1278 ctth = cosd(tth) 1279 ttth = tand(tth) 1280 return math.sqrt(radii[0]**4/(ttth**2*((radii[0]*ctth)**2+(radii[1]*stth)**2))) 1281 1282 def calcZdisCosB(radius,tth,radii): 1283 sinb = cosB = min(1.0,radii[0]**2/(radius*radii[1])) 1284 cosb = math.sqrt(1.-sinb**2) 1285 ttth = tand(tth) 1286 zdis = radii[1]*ttth*cosb/sinb 1287 return zdis,cosB 1275 1288 1276 1289 def ImageCalibrate(self,data): … … 1281 1294 scalex = data['scalex'] # = 1000./(pixelSize[0]*self.imScale) 1282 1295 scaley = data['scaley'] 1296 cutoff = data['cutoff'] 1283 1297 if len(ring) < 5: 1284 1298 print 'not enough inner ring points for ellipse' … … 1294 1308 return False 1295 1309 1296 #setup 360 points on that ring for "good" fit1297 Ring = makeRing(ellipse, scalex,scaley,self.imScale,self.ImageZ)1310 #setup 180 points on that ring for "good" fit 1311 Ring = makeRing(ellipse,20,cutoff,scalex,scaley,self.imScale,self.ImageZ) 1298 1312 if Ring: 1299 1313 ellipse = FitEllipse(Ring) … … 1302 1316 return False 1303 1317 print 'inner ring:',ellipse 1304 data['center'] = ellipse[0]#not right!! (but useful for now)1318 data['center'] = copy.copy(ellipse[0]) #not right!! (but useful for now) 1305 1319 data['ellipses'].append(ellipse[:]) 1306 1320 self.PlotImage() … … 1315 1329 A = cell2A(cell) 1316 1330 wave = data['wavelength'] 1317 dist = data['distance'] 1318 refine = data['refine'] #flags for center, wavelength,distance, tilt angle, tilt rotation 1331 cent = data['center'] 1332 pixLimit = data['pixLimit'] 1333 elcent,phi,radii = ellipse 1319 1334 HKL = GenHBravais(0.5,Bravais,A) 1320 dList = [] 1335 dsp = HKL[0][3] 1336 tth = 2.0*asind(wave/(2.*dsp)) 1337 ttth = tand(tth) 1338 data['distance'] = dist = calcDist(radii,tth) 1339 radius = dist*tand(tth) 1340 zdis,cosB = calcZdisCosB(radius,tth,radii) 1341 sinp = sind(ellipse[1]) 1342 cosp = cosd(ellipse[1]) 1343 cent1 = [] 1344 cent2 = [] 1345 Zsign = 1. 1346 xSum = 0 1347 ySum = 0 1348 phiSum = 0 1349 tiltSum = 0 1350 distSum = 0 1351 Zsum = 0 1321 1352 for i,H in enumerate(HKL): 1322 cent,phi,semiradii = ellipse 1323 meanRadius = semiradii[0]**2/semiradii[1] 1324 dsp = 1.0/math.sqrt(calc_rDsq(H,A)) 1325 tth = 2.0*asind(wave/(2.*dsp)) 1326 radius = dist*tand(tth) 1327 dList.append([dsp,radius]) 1328 semiradii[0] *= radius/meanRadius 1329 semiradii[1] *= radius/meanRadius 1330 Ring = makeRing(ellipse,scalex,scaley,self.imScale,self.ImageZ) 1353 dsp = H[3] 1354 tth = 2.0*asind(0.5*wave/dsp) 1355 stth = sind(tth) 1356 ctth = cosd(tth) 1357 ttth = tand(tth) 1358 radius = dist*ttth 1359 elcent,phi,radii = ellipse 1360 radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2) 1361 radii[0] = math.sqrt(radii[1]*radius*cosB) 1362 zdis,cosB = calcZdisCosB(radius,tth,radii) 1363 zdis *= Zsign 1364 sinp = sind(phi) 1365 cosp = cosd(phi) 1366 cent = data['center'] 1367 elcent = [cent[0]+zdis*sinp,cent[1]-zdis*cosp] 1368 ratio = radii[1]/radii[0] 1369 Ring = makeRing(ellipse,pixLimit,cutoff,scalex,scaley,self.imScale,self.ImageZ) 1331 1370 if Ring: 1371 numZ = len(ring) 1332 1372 data['rings'].append(Ring) 1333 ellipse = FitEllipse(Ring) 1334 cosb = (1.-ellipse[2][0]**2/ellipse[2][1]**2)*cosd(tth)**2 1335 print 'for ring #',i,ellipse[1],cosb,(ellipse[2][0]**2/ellipse[2][1])/tand(tth) 1373 elcent,phi,radii = ellipse = FitEllipse(Ring) 1374 if abs(phi) > 45. and phi < 0.: 1375 phi += 180. 1376 dist = calcDist(radii,tth) 1377 distR = 1.-dist/data['distance'] 1378 if distR > 0.001: 1379 print 'Wavelength too large?' 1380 elif distR < -0.001: 1381 print 'Wavelength too small?' 1382 else: 1383 if abs((radii[1]/radii[0]-ratio)/ratio) > 0.01: 1384 print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i) 1385 return False 1386 zdis,cosB = calcZdisCosB(radius,tth,radii) 1387 Tilt = acosd(cosB) 1388 sinp = sind(ellipse[1]) 1389 cosp = cosd(ellipse[1]) 1390 cent1.append(np.array([elcent[0]+zdis*sinp,elcent[1]-zdis*cosp])) 1391 cent2.append(np.array([elcent[0]-zdis*sinp,elcent[1]+zdis*cosp])) 1392 dist1 = dist2 = 0 1393 if i: 1394 d1 = cent1[-1]-cent1[-2] 1395 dist1 = np.dot(d1,d1) 1396 d2 = cent2[-1]-cent2[-2] 1397 dist2 = np.dot(d2,d2) 1398 if dist2 > dist1: 1399 data['center'] = cent1[-1] 1400 Zsign *= -1. 1401 else: 1402 data['center'] = cent2[-1] 1403 Zsign = 1. 1404 Zsum += numZ 1405 phiSum += numZ*phi 1406 distSum += numZ*dist 1407 xSum += numZ*data['center'][0] 1408 ySum += numZ*data['center'][1] 1409 tiltSum += numZ*Tilt 1410 cent = data['center'] 1411 print 'for ring # %i dist %.3f rotate %.2f tilt %.2f Xcent %.3f Ycent %.3f' \ 1412 %(i,dist,phi,Tilt,cent[0],cent[1]) 1336 1413 data['ellipses'].append(copy.deepcopy(ellipse)) 1337 1414 self.PlotImage() 1338 1415 else: 1339 1416 break 1417 fullSize = len(self.ImageZ)/(self.imScale*scalex) 1418 if 2*radii[1] < .9*fullSize: 1419 print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib' 1420 if Zsum: 1421 data['center'] = [xSum/Zsum,ySum/Zsum] 1422 data['tilt'] = tiltSum/Zsum 1423 data['distance'] = distSum/Zsum 1424 data['rotation'] = phiSum/Zsum 1425 print data['center'],data['tilt'],data['distance'],data['rotation'] 1426 else: 1427 print 'Only one ring fitted. Check your wavelength.' 1428 return False 1340 1429 self.PlotImage() 1341 1342 1343 1344 1345 1346 1347 1348 1430 1349 1431 return True -
trunk/GSASIIgrid.py
r31 r32 1109 1109 data['calibrant'] = calSel.GetValue() 1110 1110 1111 def OnPixLimit(event): 1112 data['pixLimit'] = int(pixLimit.GetValue()) 1113 1111 1114 def OnMaxSlider(event): 1112 1115 imax = max(data['range'][1][0],int(maxSel.GetValue())) … … 1145 1148 distSel.SetValue("%8.3f"%(data['distance'])) #reset in case of error 1146 1149 1147 def On ImRefine(event):1148 ImageCalibRef[0] = centRef.GetValue()1149 ImageCalibRef[1] = waveRef.GetValue()1150 ImageCalibRef[2] = distRef.GetValue()1151 ImageCalibRef[3] = tiltRef.GetValue()1152 ImageCalibRef[4] = rotRef.GetValue()1153 SetStatusLine()1154 1150 def OnCutOff(event): 1151 try: 1152 cutoff = float(cutOff.GetValue()) 1153 data['cutoff'] = cutoff 1154 except ValueError: 1155 pass 1156 cutOff.SetValue("%.1f"%(data['cutoff'])) #reset in case of error 1157 1155 1158 def OnShowLines(event): 1156 1159 if data['showLines']: … … 1180 1183 data['setDefault'] = True 1181 1184 1185 def OnSetRings(event): 1186 if data['setRings']: 1187 data['setRings'] = False 1188 else: 1189 data['setRings'] = True 1190 setRings.SetValue(data['setRings']) 1191 self.PlotImage() 1192 1182 1193 def OnClearCalib(event): 1183 1194 data['ring'] = [] … … 1187 1198 1188 1199 def OnCalibrate(event): 1200 data['setRings'] = False 1201 setRings.SetValue(data['setRings']) 1189 1202 msg = \ 1190 '''Select > 5points on inner ring of image pattern.1203 '''Select > 4 points on inner ring of image pattern. 1191 1204 Click right mouse button to select point. 1192 1205 Use left mouse button to delete point. … … 1200 1213 Status.SetStatusText('Calibration successful') 1201 1214 cent = data['center'] 1202 self.centText.SetValue(("%8.3f,%8.3f" % (cent[0],cent[1]))) 1215 centText.SetValue(("%8.3f,%8.3f" % (cent[0],cent[1]))) 1216 distSel.SetValue("%8.3f"%(data['distance'])) 1217 tiltSel.SetValue("%9.3f"%(data['tilt'])) 1218 rotSel.SetValue("%9.3f"%(data['rotation'])) 1203 1219 else: 1204 1220 Status.SetStatusText('Calibration failed') … … 1208 1224 1209 1225 def SetStatusLine(): 1210 if data['refine'][0]: 1211 Status.SetStatusText("On Image: key 'c' to mark center") 1212 else: 1213 Status.SetStatusText("") 1226 Status.SetStatusText("On Image: key 'c' to mark center") 1214 1227 1215 1228 colorList = [m for m in mpl.cm.datad.keys() if not m.endswith("_r")] … … 1243 1256 mainSizer.Add(maxSizer,1,wx.EXPAND|wx.RIGHT) 1244 1257 1245 comboSizer = wx. FlexGridSizer(2,4,5,5)1258 comboSizer = wx.BoxSizer(wx.HORIZONTAL) 1246 1259 comboSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Color bar '),0, 1247 1260 wx.ALIGN_CENTER_VERTICAL) … … 1257 1270 calSel.Bind(wx.EVT_COMBOBOX, OnNewCalibrant) 1258 1271 comboSizer.Add(calSel,0,wx.ALIGN_CENTER_VERTICAL) 1272 comboSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Pixel search range '),0, 1273 wx.ALIGN_CENTER_VERTICAL) 1274 pixLimit = wx.ComboBox(parent=self.dataDisplay,value=str(data['pixLimit']),choices=['5','10','15','20'], 1275 style=wx.CB_READONLY|wx.CB_DROPDOWN) 1276 pixLimit.Bind(wx.EVT_COMBOBOX, OnPixLimit) 1277 comboSizer.Add(pixLimit,0,wx.ALIGN_CENTER_VERTICAL) 1278 comboSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Min ring I/Ib '),0, 1279 wx.ALIGN_CENTER_VERTICAL) 1280 cutOff = wx.TextCtrl(parent=self.dataDisplay,value=("%.1f" % (data['cutoff'])), 1281 style=wx.TE_PROCESS_ENTER) 1282 cutOff.Bind(wx.EVT_TEXT_ENTER,OnCutOff) 1283 comboSizer.Add(cutOff,0,wx.ALIGN_CENTER_VERTICAL) 1284 1259 1285 mainSizer.Add(comboSizer,0,wx.ALIGN_CENTER_HORIZONTAL) 1260 1286 mainSizer.Add((5,5),0) 1261 1287 1262 dataSizer = wx.FlexGridSizer(6,5,5,5) 1263 dataSizer.Add((5,0),0) 1288 dataSizer = wx.FlexGridSizer(6,4,5,5) 1264 1289 dataSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Calibration coefficients'),0, 1265 1290 wx.ALIGN_CENTER_VERTICAL) … … 1269 1294 dataSizer.Add((5,0),0) 1270 1295 1271 ImageCalibRef = data['refine']1272 1296 dataSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Beam center X,Y'),0, 1273 1297 wx.ALIGN_CENTER_VERTICAL) 1274 1298 cent = data['center'] 1275 self.centText = wx.TextCtrl(parent=self.dataDisplay,value=("%8.3f,%8.3f" % (cent[0],cent[1])),style=wx.TE_READONLY) 1276 dataSizer.Add(self.centText,0,wx.ALIGN_CENTER_VERTICAL) 1277 centRef = wx.CheckBox(parent=self.dataDisplay,label='refine?') 1278 centRef.Bind(wx.EVT_CHECKBOX, OnImRefine) 1279 centRef.SetValue(ImageCalibRef[0]) 1280 dataSizer.Add(centRef,0,wx.ALIGN_CENTER_VERTICAL) 1299 centText = wx.TextCtrl(parent=self.dataDisplay,value=("%8.3f,%8.3f" % (cent[0],cent[1])),style=wx.TE_READONLY) 1300 dataSizer.Add(centText,0,wx.ALIGN_CENTER_VERTICAL) 1281 1301 1282 1302 dataSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Inner/Outer radii'),0, … … 1293 1313 waveSel.Bind(wx.EVT_TEXT_ENTER,OnWavelength) 1294 1314 dataSizer.Add(waveSel,0,wx.ALIGN_CENTER_VERTICAL) 1295 waveRef = wx.CheckBox(parent=self.dataDisplay,label='refine?')1296 waveRef.Bind(wx.EVT_CHECKBOX, OnImRefine)1297 waveRef.SetValue(ImageCalibRef[1])1298 dataSizer.Add(waveRef,0,wx.ALIGN_CENTER_VERTICAL)1299 1315 1300 1316 dataSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Start/End azimuth'),0, … … 1310 1326 distSel.Bind(wx.EVT_TEXT_ENTER,OnDistance) 1311 1327 dataSizer.Add(distSel,0,wx.ALIGN_CENTER_VERTICAL) 1312 distRef = wx.CheckBox(parent=self.dataDisplay,label='refine?')1313 distRef.Bind(wx.EVT_CHECKBOX, OnImRefine)1314 distRef.SetValue(ImageCalibRef[2])1315 dataSizer.Add(distRef,0,wx.ALIGN_CENTER_VERTICAL)1316 1328 1317 1329 dataSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' No. bins'),0, … … 1323 1335 dataSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Tilt angle'),0, 1324 1336 wx.ALIGN_CENTER_VERTICAL) 1325 self.tiltSel = wx.TextCtrl(parent=self.dataDisplay,value=("%9.3f"%(data['tilt'])),style=wx.TE_READONLY) 1326 dataSizer.Add(self.tiltSel,0,wx.ALIGN_CENTER_VERTICAL) 1327 tiltRef = wx.CheckBox(parent=self.dataDisplay,label='refine?') 1328 tiltRef.Bind(wx.EVT_CHECKBOX, OnImRefine) 1329 tiltRef.SetValue(ImageCalibRef[3]) 1330 dataSizer.Add(tiltRef,0,wx.ALIGN_CENTER_VERTICAL) 1337 tiltSel = wx.TextCtrl(parent=self.dataDisplay,value=("%9.3f"%(data['tilt'])),style=wx.TE_READONLY) 1338 dataSizer.Add(tiltSel,0,wx.ALIGN_CENTER_VERTICAL) 1331 1339 showLines = wx.CheckBox(parent=self.dataDisplay,label='Show integration limits?') 1332 1340 dataSizer.Add(showLines,0) … … 1340 1348 dataSizer.Add(wx.StaticText(parent=self.dataDisplay,label=' Tilt rotation'),0, 1341 1349 wx.ALIGN_CENTER_VERTICAL) 1342 self.rotSel = wx.TextCtrl(parent=self.dataDisplay,value=("%9.3f"%(data['rotation'])),style=wx.TE_READONLY) 1343 dataSizer.Add(self.rotSel,0,wx.ALIGN_CENTER_VERTICAL) 1344 rotRef = wx.CheckBox(parent=self.dataDisplay,label='refine?') 1345 rotRef.Bind(wx.EVT_CHECKBOX, OnImRefine) 1346 rotRef.SetValue(ImageCalibRef[4]) 1347 dataSizer.Add(rotRef,0,) 1350 rotSel = wx.TextCtrl(parent=self.dataDisplay,value=("%9.3f"%(data['rotation'])),style=wx.TE_READONLY) 1351 dataSizer.Add(rotSel,0,wx.ALIGN_CENTER_VERTICAL) 1348 1352 setDefault = wx.CheckBox(parent=self.dataDisplay,label='Use as default for all images?') 1349 1353 dataSizer.Add(setDefault,0) 1350 1354 setDefault.Bind(wx.EVT_CHECKBOX, OnSetDefault) 1351 1355 setDefault.SetValue(data['setDefault']) 1352 dataSizer.Add((10,5),0) 1356 setRings = wx.CheckBox(parent=self.dataDisplay,label='Show ring picks?') 1357 dataSizer.Add(setRings,0) 1358 setRings.Bind(wx.EVT_CHECKBOX, OnSetRings) 1359 setRings.SetValue(data['setRings']) 1353 1360 1354 1361 mainSizer.Add(dataSizer,0) … … 1741 1748 else: 1742 1749 self.dataFrame = DataFrame(parent=self.mainPanel) 1743 1750 1751 self.dataFrame.Raise() 1744 1752 self.PickId = 0 1745 1753 self.PatternId = 0
Note: See TracChangeset
for help on using the changeset viewer.