Changeset 48
- Timestamp:
- Apr 12, 2010 1:56:43 PM (13 years ago)
- Location:
- trunk
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASII.py
r45 r48 409 409 Data['calibrant'] = '' 410 410 Data['IOtth'] = [2.0,5.0] 411 Data['LRazimuth'] = [- 45,45]411 Data['LRazimuth'] = [-135,-45] 412 412 Data['outChannels'] = 2500 413 413 Data['fullIntegrate'] = False -
trunk/GSASIIcomp.py
r47 r48 1206 1206 return False,0,0 1207 1207 1208 def FitRing(ring): 1209 Err,parms = FitCircle(ring) 1210 Err /= len(ring) 1211 print 'circle error:','%8f'%(Err) 1212 if Err > 20000.: 1213 eparms = FitEllipse(ring) 1214 if eparms: 1215 parms = eparms 1216 return parms 1217 1218 def FitCircle(ring): 1219 import numpy.linalg as nl 1220 1221 def makeParmsCircle(B): 1222 cent = [-B[0]/2,-B[1]/2] 1223 phi = 0. 1224 sr1 = sr2 = math.sqrt(cent[0]**2+cent[1]**2-B[2]) 1225 return cent,phi,[sr1,sr2] 1226 1227 ring = np.array(ring) 1228 x = np.asarray(ring.T[0]) 1229 y = np.asarray(ring.T[1]) 1230 1231 M = np.array((x,y,np.ones_like(x))) 1232 B = np.array(-(x**2+y**2)) 1233 result = nl.lstsq(M.T,B) 1234 return result[1],makeParmsCircle(result[0]) 1235 1208 1236 def FitEllipse(ring): 1237 import numpy.linalg as nl 1238 1239 def makeParmsEllipse(B): 1240 det = 4.*(1.-B[0]**2)-B[1]**2 1241 if det < 0.: 1242 print 'hyperbola!' 1243 return 0 1244 elif det == 0.: 1245 print 'parabola!' 1246 return 0 1247 cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \ 1248 (B[1]*B[2]-2.*(1.+B[0])*B[3])/det] 1249 phi = 0.5*atand(0.5*B[1]/B[0]) 1250 1251 a = (1.+B[0])/cosd(2*phi) 1252 b = 2.-a 1253 f = (1.+B[0])*cent[0]**2+(1.-B[0])*cent[1]**2+B[1]*cent[0]*cent[1]-B[4] 1254 if f/a < 0 or f/b < 0: 1255 return 0 1256 sr1 = math.sqrt(f/a) 1257 sr2 = math.sqrt(f/b) 1258 if sr1 > sr2: 1259 sr1,sr2 = SwapXY(sr1,sr2) 1260 phi -= 90. 1261 if phi < -90.: 1262 phi += 180. 1263 return cent,phi,[sr1,sr2] 1264 1265 ring = np.array(ring) 1266 x = np.asarray(ring.T[0]) 1267 y = np.asarray(ring.T[1]) 1268 M = np.array((x**2-y**2,x*y,x,y,np.ones_like(x))) 1269 B = np.array(-(x**2+y**2)) 1270 result = nl.lstsq(M.T,B) 1271 return makeParmsEllipse(result[0]) 1272 1273 def FitDetector(rings,p0,wave): 1209 1274 from scipy.optimize import leastsq 1210 def ellipseCalc(A,xy): 1211 x = xy[0] 1212 y = xy[1] 1213 return (1.+A[0])*x**2+(1.-A[0])*y**2+A[1]*x*y+A[2]*x+A[3]*y+A[4] 1214 ring = np.array(ring) 1215 p0 = np.zeros(shape=5) 1216 result = leastsq(ellipseCalc,p0,args=(ring.T,)) 1217 B = result[0] 1218 1219 ell = [1.+B[0],B[1],1.-B[0],B[2],B[3],B[4]] 1220 # ell is [A,B,C,D,E,F] for Ax^2+Bxy+Cy^2+Dx+Ey+F = 0 1221 det = 4.*ell[0]*ell[2]-ell[1]**2 1222 if det < 0.: 1223 print 'hyperbola!' 1224 return 0 1225 elif det == 0.: 1226 print 'parabola!' 1227 return 0 1228 cent = [(ell[1]*ell[4]-2.0*ell[2]*ell[3])/det, \ 1229 (ell[1]*ell[3]-2.0*ell[0]*ell[4])/det] 1230 phi = 0.5*atand(ell[1]/(ell[0]-ell[2]+1e-32)) 1231 1232 a3 = 0.5*(ell[0]+ell[2]+(ell[0]-ell[2])/cosd(2.0*phi)) 1233 b3 = ell[0]+ell[2]-a3 1234 f3 = ell[0]*cent[0]**2+ell[2]*cent[1]**2+ell[1]*cent[0]*cent[1]-ell[5] 1235 if f3/a3 < 0 or f3/b3 < 0: 1236 return 0 1237 sr1 = math.sqrt(f3/a3) 1238 sr2 = math.sqrt(f3/b3) 1239 if sr1 > sr2: 1240 sr1,sr2 = SwapXY(sr1,sr2) 1241 phi -= 90. 1242 if phi < -90.: 1243 phi += 180. 1244 return cent,phi,[sr1,sr2] 1275 def ellipseCalc(B,xyd,wave): 1276 sind = lambda x: np.sin(x*math.pi/180.) 1277 asind = lambda x: 180.*np.arcsin(x)/math.pi 1278 cosd = lambda x: np.cos(x*math.pi/180.) 1279 tand = lambda x: np.tan(x*math.pi/180.) 1280 x = xyd[0] 1281 y = xyd[1] 1282 dsp = xyd[2] 1283 dist,x0,y0,phi,tilt = B 1284 tth = 2.0*asind(wave/(2.*dsp)) 1285 ttth = tand(tth) 1286 radius = dist*ttth 1287 stth = sind(tth) 1288 cosb = cosd(tilt) 1289 R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2) 1290 R0 = np.sqrt(R1*radius*cosb) 1291 zdis = R1*ttth*tand(tilt) 1292 X = x-x0+zdis*sind(phi) 1293 Y = y-y0-zdis*cosd(phi) 1294 XR = X*cosd(phi)-Y*sind(phi) 1295 YR = X*sind(phi)+Y*cosd(phi) 1296 return (XR/R0)**2+(YR/R1)**2-1 1297 result = leastsq(ellipseCalc,p0,args=(rings.T,wave)) 1298 return result[0] 1299 1245 1300 1246 1301 def ImageLocalMax(image,w,Xpix,Ypix): … … 1257 1312 return 0,0,0,0 1258 1313 1259 def makeRing( ellipse,pix,reject,scalex,scaley,image):1314 def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image): 1260 1315 cent,phi,radii = ellipse 1261 1316 cphi = cosd(phi) … … 1271 1326 X /= scalex #convert to mm 1272 1327 Y /= scaley 1273 ring.append([X,Y ])1328 ring.append([X,Y,dsp]) 1274 1329 if len(ring) < 45: #want more than 1/4 of a circle 1275 1330 return [] … … 1302 1357 stth = sind(tth) 1303 1358 ctth = cosd(tth) 1304 cosb = sind(tilt) 1305 sinb = math.sqrt(1.-cosb**2) 1306 sinp = sind(phi) 1307 cosp = cosd(phi) 1359 cosb = cosd(tilt) 1308 1360 radius = dist*ttth 1309 radii[1] = dist*stth*ctth* sinb/(sinb**2-stth**2)1361 radii[1] = dist*stth*ctth*cosb/(cosb**2-stth**2) 1310 1362 if radii[1] > 0: 1311 radii[0] = math.sqrt(radii[1]*radius* sinb)1312 zdis = radii[1]*ttth* cosb/sinb1313 elcent = [cent[0]-zdis*sin p,cent[1]+zdis*cosp]1363 radii[0] = math.sqrt(radii[1]*radius*cosb) 1364 zdis = radii[1]*ttth*tand(tilt) 1365 elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)] 1314 1366 return elcent,phi,radii 1315 1367 else: 1316 1368 return False 1317 1369 1370 def GetDetectorXY(dsp,azm,data): 1371 from scipy.optimize import fsolve 1372 def func(xy,*args): 1373 azm,phi,R0,R1,A,B = args 1374 cp = cosd(phi) 1375 sp = sind(phi) 1376 x,y = xy 1377 out = [] 1378 out.append(y-x*tand(azm)) 1379 out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2) 1380 return out 1381 elcent,phi,radii = GetEllipse(dsp,data) 1382 cent = data['center'] 1383 tilt = data['tilt'] 1384 phi = data['rotation'] 1385 wave = data['wavelength'] 1386 dist = data['distance'] 1387 tth = 2.0*asind(wave/(2.*dsp)) 1388 ttth = tand(tth) 1389 radius = dist*ttth 1390 stth = sind(tth) 1391 cosb = cosd(tilt) 1392 R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2) 1393 R0 = math.sqrt(R1*radius*cosb) 1394 zdis = R1*ttth*tand(tilt) 1395 A = zdis*sind(phi) 1396 B = -zdis*cosd(phi) 1397 xy0 = [radius*cosd(azm),radius*sind(azm)] 1398 xy = fsolve(func,xy0,args=(azm,phi,R0,R1,A,B))+cent 1399 return xy 1400 1318 1401 def GetTthDspAzm(x,y,data): 1402 #make numpy array compliant 1403 atand = lambda x: 180.*np.arctan(x)/np.pi 1404 sind = lambda x: np.sin(x*np.pi/180.) 1405 atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi 1319 1406 wave = data['wavelength'] 1320 1407 dist = data['distance'] … … 1327 1414 X = np.sum(X*makeMat(-phi,2),axis=1) 1328 1415 X = np.sum(X*makeMat(-tilt,0),axis=1) 1329 tth = atand( math.sqrt(dx**2+dy**2-X[2]**2)/(dist-X[2]))1416 tth = atand(np.sqrt(dx**2+dy**2-X[2]**2)/(dist-X[2])) 1330 1417 dsp = wave/(2.*sind(tth/2.)) 1331 1418 azm = atan2d(dy,dx) … … 1359 1446 #fit start points on inner ring 1360 1447 data['ellipses'] = [] 1361 out B = FitEllipse(ring)1362 if out B:1363 ellipse = outB1364 print 'start:',ellipse1448 outE = FitRing(ring) 1449 if outE: 1450 print 'start ellipse:',outE 1451 ellipse = outE 1365 1452 else: 1366 1453 return False 1367 1454 1368 1455 #setup 180 points on that ring for "good" fit 1369 Ring = makeRing( ellipse,20,cutoff,scalex,scaley,self.ImageZ)1456 Ring = makeRing(1.0,ellipse,20,cutoff,scalex,scaley,self.ImageZ) 1370 1457 if Ring: 1371 ellipse = Fit Ellipse(Ring)1372 Ring = makeRing( ellipse,20,cutoff,scalex,scaley,self.ImageZ) #do again1373 ellipse = Fit Ellipse(Ring)1458 ellipse = FitRing(Ring) 1459 Ring = makeRing(1.0,ellipse,20,cutoff,scalex,scaley,self.ImageZ) #do again 1460 ellipse = FitRing(Ring) 1374 1461 else: 1375 1462 print '1st ring not sufficiently complete to proceed' … … 1400 1487 radius = dist*tand(tth) 1401 1488 zdis,cosB = calcZdisCosB(radius,tth,radii) 1402 sinp = sind(ellipse[1])1403 cosp = cosd(ellipse[1])1404 1489 cent1 = [] 1405 1490 cent2 = [] … … 1428 1513 elcent = [cent[0]+zsinp,cent[1]-zcosp] 1429 1514 ratio = radii[1]/radii[0] 1430 Ring = makeRing( ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)1515 Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) 1431 1516 if Ring: 1432 1517 numZ = len(Ring) 1433 data['rings'].append(np.column_stack((np.array(Ring),dsp*np.ones(numZ)))) 1434 elcent,phi,radii = ellipse = FitEllipse(Ring) 1518 data['rings'].append(np.array(Ring)) 1519 ellipse = FitRing(Ring) 1520 elcent,phi,radii = ellipse 1435 1521 if abs(phi) > 45. and phi < 0.: 1436 1522 phi += 180. … … 1499 1585 phi = data['rotation'] = phiSum/Zsum 1500 1586 rings = np.concatenate((data['rings']),axis=0) 1501 print wave,dist,cent,phi,tilt 1502 1503 1587 p0 = [dist,cent[0],cent[1],phi,tilt] 1588 result = FitDetector(rings,p0,wave) 1589 data['distance'] = result[0] 1590 data['center'] = result[1:3] 1591 data['rotation'] = np.mod(result[3],360.0) 1592 data['tilt'] = result[4] 1504 1593 N = len(data['ellipses']) 1594 data['ellipses'] = [] #clear away individual ellipse fits 1505 1595 for H in HKL[:N]: 1506 1596 ellipse = GetEllipse(H[3],data) … … 1509 1599 1510 1600 return True 1511 1512 1601 1513 1602 def test(): 1514 1603 cell = [5,5,5,90,90,90] -
trunk/GSASIIgrid.py
r43 r48 1164 1164 if data['fullIntegrate']: 1165 1165 data['fullIntegrate'] = False 1166 data['LRazimuth'] = [-45,45] 1167 self.Lazim.SetValue("%6d" % (-45)) 1168 self.Razim.SetValue("%6d" % (45)) 1166 self.Lazim.SetEditable(True) 1167 self.Razim.SetEditable(True) 1169 1168 else: 1170 1169 data['fullIntegrate'] = True 1171 data['LRazimuth'] = [0,360] 1172 self.Lazim.SetValue("%6d" % (0)) 1173 self.Razim.SetValue("%6d" % (360)) 1170 self.Lazim.SetEditable(False) 1171 self.Razim.SetEditable(False) 1174 1172 G2plt.PlotImage(self) 1175 1173 -
trunk/GSASIIplot.py
r47 r48 186 186 ypos = event.ydata-ycent 187 187 if 'line3' in str(item) or 'line4' in str(item) and not Data['fullIntegrate']: 188 ang = int(atan2d( -ypos,xpos))188 ang = int(atan2d(xpos,ypos)) 189 189 imgPage.canvas.SetToolTipString('%6d deg'%(ang)) 190 elif ' inner' in str(item.get_gid()) or 'outer' in str(item.get_gid()):190 elif 'line1' in str(item) or 'line2' in str(item): 191 191 tth = G2cmp.GetTth(event.xdata,event.ydata,Data) 192 192 imgPage.canvas.SetToolTipString('%8.3fdeg'%(tth)) … … 335 335 IOtth = Data['IOtth'] 336 336 wave = Data['wavelength'] 337 dsp = wave/(2.0*sind(IOtth[0]/2.0))338 ell 0 = G2cmp.GetEllipse(dsp,Data) #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola)339 dsp = wave/(2.0*sind(IOtth[1]/2.0))340 ell 1 = G2cmp.GetEllipse(dsp,Data) #Ditto & more likely for outer ellipse337 dspI = wave/(2.0*sind(IOtth[0]/2.0)) 338 ellI = G2cmp.GetEllipse(dspI,Data) #=False if dsp didn't yield an ellipse (ugh! a parabola or a hyperbola) 339 dspO = wave/(2.0*sind(IOtth[1]/2.0)) 340 ellO = G2cmp.GetEllipse(dspO,Data) #Ditto & more likely for outer ellipse 341 341 if Data['fullIntegrate']: 342 arcxI = arcyI = np.array(range(0,361)) 343 arcxO = arcyO = np.array(range(0,361)) 342 Azm = np.array(range(0,361)) 344 343 else: 345 arcxI = arcyI= np.array(range(LRAzim[0],LRAzim[1]+1))346 arcxO = arcyO = np.array(range(LRAzim[0],LRAzim[1]+1))347 if ell0:348 cent0,phi0,radii0 = ell0349 radii0 = np.sum(G2cmp.makeMat(phi0,2).T*np.array(radii0+[0.,]),axis=1)350 arcxI = np.cos(arcxI*math.pi/180.)*radii0[0]+cent0[0]351 arc yI = np.sin(arcyI*math.pi/180.)*radii0[1]+cent0[1]344 Azm = np.array(range(LRAzim[0],LRAzim[1]+1)) 345 if ellI: 346 xyI = [] 347 for azm in Azm: 348 xyI.append(G2cmp.GetDetectorXY(dspI,azm,Data)) 349 xyI = np.array(xyI) 350 arcxI,arcyI = xyI.T 352 351 imgPlot.plot(arcxI,arcyI,picker=3) 353 if ell1: 354 cent1,phi1,radii1 = ell1 355 radii1 = np.sum(G2cmp.makeMat(phi1,2).T*np.array(radii1+[0.,]),axis=1) 356 arcxO = np.cos(arcxO*math.pi/180.)*radii1[0]+cent1[0] 357 arcyO = np.sin(arcyO*math.pi/180.)*radii1[1]+cent1[1] 358 imgPlot.plot(arcxO,arcyO,picker=3) 359 if ell0 and ell1 and not Data['fullIntegrate']: 352 if ellO: 353 xyO = [] 354 for azm in Azm: 355 xyO.append(G2cmp.GetDetectorXY(dspO,azm,Data)) 356 xyO = np.array(xyO) 357 arcxO,arcyO = xyO.T 358 imgPlot.plot(arcxO,arcyO,picker=3) 359 if ellO and ellI and not Data['fullIntegrate']: 360 360 imgPlot.plot([arcxI[0],arcxO[0]],[arcyI[0],arcyO[0]],picker=3) 361 361 imgPlot.plot([arcxI[-1],arcxO[-1]],[arcyI[-1],arcyO[-1]],picker=3) … … 363 363 imgPlot.text(xring,yring,'+',color='b',ha='center',va='center',picker=3) 364 364 if Data['setRings']: 365 rings = Data['rings']366 for xring,yring, tthin rings:365 rings = np.concatenate((Data['rings']),axis=0) 366 for xring,yring,dsp in rings: 367 367 imgPlot.text(xring,yring,'+',ha='center',va='center') 368 368 for ellipse in Data['ellipses']:
Note: See TracChangeset
for help on using the changeset viewer.