Changeset 1148
- Timestamp:
- Nov 23, 2013 5:29:28 PM (9 years ago)
- Location:
- trunk
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIIO.py
r1147 r1148 410 410 File = open(filename,'rb') 411 411 dataType = 5 412 center = [None,None] 413 wavelength = None 414 distance = None 412 415 try: 413 416 Meta = open(filename+'.metadata','Ur') … … 464 467 pixy = (marFrame.pixelsizeX/1000.0,marFrame.pixelsizeY/1000.0) 465 468 head = marFrame.outputHead() 469 # extract resonable wavelength from header 470 wavelength = marFrame.sourceWavelength*1e-5 471 wavelength = (marFrame.opticsWavelength > 0) and marFrame.opticsWavelength*1e-5 or wavelength 472 wavelength = (wavelength <= 0) and None or wavelength 473 # extract resonable distance from header 474 distance = (marFrame.startXtalToDetector+marFrame.endXtalToDetector)*5e-4 475 distance = (distance <= 0) and marFrame.xtalToDetector*1e-3 or distance 476 distance = (distance <= 0) and None or distance 477 # extract resonable center from header 478 center = [marFrame.beamX*marFrame.pixelsizeX*1e-9,marFrame.beamY*marFrame.pixelsizeY*1e-9] 479 center = (center[0] != 0 and center[1] != 0) and center or [None,None] 480 #print head,tifType,pixy 466 481 elif 272 in IFD: 467 482 ifd = IFD[272] … … 561 576 562 577 image = np.reshape(image,(sizexy[1],sizexy[0])) 563 center = [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000] 564 data = {'pixelSize':pixy,'wavelength':0.10,'distance':100.0,'center':center,'size':sizexy} 578 center = (not center[0]) and [pixy[0]*sizexy[0]/2000,pixy[1]*sizexy[1]/2000] or center 579 wavelength = (not wavelength) and 0.10 or wavelength 580 distance = (not distance) and 100.0 or distance 581 data = {'pixelSize':pixy,'wavelength':wavelength,'distance':distance,'center':center,'size':sizexy} 565 582 File.close() 566 583 if imageOnly: -
trunk/GSASIIimage.py
r1133 r1148 15 15 16 16 needs some minor refactoring to remove wx code 17 actually not easy in this case wx.ProgressDialog needs #of blocks to process when started 18 not known in G2imageGUI. 17 19 ''' 18 20 … … 44 46 npasind = lambda x: 180.*np.arcsin(x)/np.pi 45 47 npcosd = lambda x: np.cos(x*np.pi/180.) 48 npacosd = lambda x: 180.*np.arccos(x)/np.pi 46 49 nptand = lambda x: np.tan(x*np.pi/180.) 47 50 npatand = lambda x: 180.*np.arctan(x)/np.pi … … 91 94 else: 92 95 err,parms = FitCircle(ring) 93 return parms 96 return parms,err 94 97 95 98 def FitCircle(ring): … … 118 121 if det < 0.: 119 122 print 'hyperbola!' 123 print B 120 124 return 0 121 125 elif det == 0.: 122 126 print 'parabola!' 127 print B 123 128 return 0 124 129 cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \ … … 150 155 return err,makeParmsEllipse(bb) 151 156 152 def FitDetector(rings,varyList,parmDict ):157 def FitDetector(rings,varyList,parmDict,Print=True): 153 158 'Needs a doc string' 154 159 … … 173 178 174 179 def ellipseCalcD(B,xyd,varyList,parmDict): 175 x = xyd[0] 176 y = xyd[1] 177 dsp = xyd[2] 180 x,y,dsp = xyd 178 181 wave = parmDict['wave'] 179 182 if 'dep' in varyList: … … 187 190 dxy = peneCorr(tth,dep) 188 191 ttth = nptand(tth) 189 radius = (dist+dxy)*ttth190 192 stth = npsind(tth) 191 193 cosb = npcosd(tilt) 192 R1 = (dist+dxy)*stth*npcosd(tth)*cosb/(cosb**2-stth**2) 193 R0 = np.sqrt(R1*radius*cosb) 194 zdis = R1*ttth*nptand(tilt) 195 X = x-x0+zdis*npsind(phi) 194 tanb = nptand(tilt) 195 tbm = nptand((tth-tilt)/2.) 196 tbp = nptand((tth+tilt)/2.) 197 sinb = npsind(tilt) 198 d = dist+dxy 199 fplus = d*tanb*stth/(cosb+stth) 200 fminus = d*tanb*stth/(cosb-stth) 201 vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth) 202 vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth) 203 R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2. #+minor axis 204 R1 = (vplus+vminus)/2. #major axis 205 zdis = (fplus-fminus)/2. 206 X = x-x0-zdis*npsind(phi) 196 207 Y = y-y0-zdis*npcosd(phi) 197 XR = X*npcosd(phi) -Y*npsind(phi)208 XR = X*npcosd(phi)+Y*npsind(phi) 198 209 YR = X*npsind(phi)+Y*npcosd(phi) 199 210 return (XR/R0)**2+(YR/R1)**2-1 … … 202 213 fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.3f','%12.5f'] 203 214 p0 = [parmDict[key] for key in varyList] 204 result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True) 215 result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8) 216 chisq = np.sum(result[2]['fvec']**2) 205 217 parmDict.update(zip(varyList,result[0])) 206 218 vals = list(result[0]) … … 212 224 sigList[i] = sig[varyList.index(name)] 213 225 ValSig = zip(names,fmt,vals,sig) 214 CalibPrint(ValSig) 215 # try: 216 # print 'Trial refinement of wavelength - not used for calibration' 217 # p0 = result[0] 218 # print p0 219 # print parms 220 # p0 = np.append(p0,parms[0]) 221 # resultW = leastsq(ellipseCalcW,p0,args=(rings.T,parms[1:]),full_output=True) 222 # resultW[0][3] = np.mod(result[0][3],360.0) #remove random full circles 223 # sig = np.sqrt(np.sum(ellipseCalcW(resultW[0],rings.T)**2)) 224 # ValSig = zip(names,fmt,resultW[0],sig*np.sqrt(np.diag(resultW[1]))) 225 # CalibPrint(ValSig) 226 # return result[0],resultW[0][-1] 227 # except ValueError: 228 # print 'Bad refinement - no result' 229 # return result[0],wave 226 if Print: 227 CalibPrint(ValSig) 228 return chi,chisq 230 229 231 230 def ImageLocalMax(image,w,Xpix,Ypix): … … 247 246 def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image): 248 247 'Needs a doc string' 248 def ellipseC(): 249 'compute estimate of ellipse circumference' 250 if radii[0] < 0: #hyperbola 251 theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2)) 252 print theta 253 return 0 254 apb = radii[1]+radii[0] 255 amb = radii[1]-radii[0] 256 return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2))) 249 257 cent,phi,radii = ellipse 250 258 cphi = cosd(phi) 251 259 sphi = sind(phi) 252 260 ring = [] 261 sumI = 0 253 262 amin = 0 254 263 amax = 360 255 for a in range(0, 360,1):264 for a in range(0,int(ellipseC()),1): 256 265 x = radii[0]*cosd(a) 257 266 y = radii[1]*sind(a) … … 260 269 X,Y,I,J = ImageLocalMax(image,pix,X,Y) 261 270 if I and J and I/J > reject: 271 sumI += I/J 262 272 X += .5 #set to center of pixel 263 273 Y += .5 … … 269 279 ring.append([X,Y,dsp]) 270 280 delt = amax-amin 271 if len(ring) < 10: #want more than 20 deg272 return [], delt > 90273 return ring, delt > 90281 if len(ring) < 10: #want more than 10 deg 282 return [],(delt > 90),0 283 return ring,(delt > 90),sumI/len(ring) 274 284 275 285 def makeIdealRing(ellipse,azm=None): … … 286 296 287 297 a = np.linspace(aR[0],aR[1],aR[2]) 288 x = radii[0]*npcosd(a-phi) 289 y = radii[1]*npsind(a-phi) 290 X = (cphi*x-sphi*y+cent[0]) 291 Y = (sphi*x+cphi*y+cent[1]) 292 return zip(X,Y) 293 298 slr = radii[0]**2/radii[1] 299 rat2 = (radii[0]/radii[1])**2 300 if radii[0] > 0.: #ellipse 301 ecc = np.sqrt(max(1.-rat2,0.)) 302 else: #hyperbola 303 ecc = np.sqrt(1.+rat2) 304 rad = slr/(1.+ecc*npcosd(a)) 305 xy = np.array([rad*npsind(a)+cent[0],rad*npcosd(a)+cent[1]]) 306 return xy 307 294 308 def calcDist(radii,tth): 295 309 'Needs a doc string' … … 309 323 zdis = radii[1]*ttth*cosb/sinb 310 324 return zdis,cosB 325 326 def GetEllipse2(tth,dxy,dist,cent,tilt,phi): 327 '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry 328 on output 329 radii[0] (b-minor axis) set < 0. for hyperbola 330 331 ''' 332 radii = [0,0] 333 ttth = tand(tth) 334 stth = sind(tth) 335 ctth = cosd(tth) 336 cosb = cosd(tilt) 337 tanb = tand(tilt) 338 tbm = tand((tth-tilt)/2.) 339 tbp = tand((tth+tilt)/2.) 340 sinb = sind(tilt) 341 d = dist+dxy 342 if tth+tilt < 90.: #ellipse 343 fplus = d*tanb*stth/(cosb+stth) 344 fminus = d*tanb*stth/(cosb-stth) 345 vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth) 346 vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth) 347 radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2. #+minor axis 348 radii[1] = (vplus+vminus)/2. #major axis 349 zdis = (fplus-fminus)/2. 350 else: #hyperbola! 351 f = d*tanb*stth/(cosb+stth) 352 v = d*(tanb+tand(tth-tilt)) 353 delt = d*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb)) 354 eps = (v-f)/(delt-v) 355 radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.) #-minor axis 356 radii[1] = eps*(delt-f)/(eps**2-1.) #major axis 357 zdis = f+radii[1]*eps 358 elcent = [cent[0]+zdis*sind(phi),cent[1]+zdis*cosd(phi)] 359 return elcent,phi,radii 311 360 312 361 def GetEllipse(dsp,data): 313 'Needs a doc string' 314 dist = data['distance'] 362 '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry 363 as given in image controls dictionary (data) and a d-spacing (dsp) 364 ''' 315 365 cent = data['center'] 316 366 tilt = data['tilt'] 317 367 phi = data['rotation'] 318 368 dep = data['DetDepth'] 319 radii = [0,0]320 369 tth = 2.0*asind(data['wavelength']/(2.*dsp)) 321 ttth = tand(tth)322 stth = sind(tth)323 ctth = cosd(tth)324 cosb = cosd(tilt)325 370 dxy = peneCorr(tth,dep) 326 radius = ttth*(dist+dxy) 327 radii[1] = (dist+dxy)*stth*ctth*cosb/(cosb**2-stth**2) 328 if radii[1] > 0: 329 radii[0] = math.sqrt(radii[1]*radius*cosb) 330 zdis = radii[1]*ttth*tand(tilt) 331 elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)] 332 return elcent,phi,radii 333 else: 334 print 'bad ellipse - radii:',radii 335 return False 371 dist = data['distance'] 372 return GetEllipse2(tth,dxy,dist,cent,tilt,phi) 336 373 337 374 def GetDetectorXY(dsp,azm,data): 338 375 'Needs a doc string' 376 339 377 from scipy.optimize import fsolve 340 378 def func(xy,*args): … … 349 387 elcent,phi,radii = GetEllipse(dsp,data) 350 388 cent = data['center'] 351 tilt = data['tilt']352 389 phi = data['rotation'] 353 390 wave = data['wavelength'] 354 dist = data['distance'] 391 tilt = data['tilt'] 392 dist = data['distance']/cosd(tilt) 355 393 dep = data['DetDepth'] 356 394 tth = 2.0*asind(wave/(2.*dsp)) … … 371 409 def GetDetXYfromThAzm(Th,Azm,data): 372 410 'Needs a doc string' 373 dsp = data['wavelength']/(2.0*npsind(Th)) 374 411 dsp = data['wavelength']/(2.0*npsind(Th)) 375 412 return GetDetectorXY(dsp,azm,data) 376 413 377 414 def GetTthAzmDsp(x,y,data): 378 'Needs a doc string '415 'Needs a doc string - checked OK for ellipses dont know about hyperbola' 379 416 wave = data['wavelength'] 380 dist = data['distance']381 417 cent = data['center'] 382 418 tilt = data['tilt'] 419 dist = data['distance']/cosd(tilt) 383 420 phi = data['rotation'] 384 421 dep = data['DetDepth'] … … 392 429 tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z)) 393 430 dxy = peneCorr(tth,dep) 394 tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z+dxy)) 431 tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z+dxy)) #depth corr not correct for tilted detector 395 432 dsp = wave/(2.*npsind(tth/2.)) 396 azm = (npatan2d(d x,-dy)+azmthoff+720.)%360.433 azm = (npatan2d(dy,dx)+azmthoff+720.)%360. 397 434 return tth,azm,dsp 398 435 … … 428 465 def EdgeFinder(image,data): 429 466 '''this makes list of all x,y where I>edgeMin suitable for an ellipse search? 467 Not currently used but might be useful in future? 430 468 ''' 431 469 import numpy.ma as ma … … 465 503 for bravais,cell in zip(Bravais,Cells): 466 504 A = G2lat.cell2A(cell) 467 hkl = G2lat.GenHBravais(dmin,bravais,A) [skip:]505 hkl = G2lat.GenHBravais(dmin,bravais,A) 468 506 HKL += hkl 469 507 HKL = G2lat.sortHKLd(HKL,True,False) … … 471 509 parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1], 472 510 'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']} 511 Found = False 473 512 for H in HKL: 474 513 dsp = H[3] 475 514 ellipse = GetEllipse(dsp,data) 476 Ring,delt = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)515 Ring,delt,sumI = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) 477 516 if Ring: 478 # numZ = len(Ring)479 517 data['rings'].append(np.array(Ring)) 480 518 data['ellipses'].append(copy.deepcopy(ellipse+('r',))) 481 else: 519 Found = True 520 elif not Found: #skipping inner rings, keep looking until ring found 482 521 continue 522 else: #no more rings beyond edge of detector 523 break 483 524 rings = np.concatenate((data['rings']),axis=0) 484 525 if data['DetDepthRef']: … … 518 559 #fit start points on inner ring 519 560 data['ellipses'] = [] 520 outE = FitRing(ring,True) 561 outE,err = FitRing(ring,True) 562 fmt = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f' 563 fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi^2: %.3f, N: %d' 521 564 if outE: 522 print 'start ellipse:',outE565 print fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1]) 523 566 ellipse = outE 524 567 else: … … 526 569 527 570 #setup 360 points on that ring for "good" fit 528 Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)571 Ring,delt,sumI = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) 529 572 if Ring: 530 ellipse = FitRing(Ring,delt)531 Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) #do again532 ellipse = FitRing(Ring,delt)573 ellipse,err = FitRing(Ring,delt) 574 Ring,delt,sumI = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) #do again 575 ellipse,err = FitRing(Ring,delt) 533 576 else: 534 577 print '1st ring not sufficiently complete to proceed' 535 578 return False 536 print 'inner ring:',ellipse #cent,phi,radii 537 data['center'] = copy.copy(ellipse[0]) #not right!! (but useful for now) 579 print fmt2%('inner ring: ',ellipse[0][0],ellipse[0][1],ellipse[1],ellipse[2][0],ellipse[2][1],err,len(Ring)) #cent,phi,radii 538 580 data['ellipses'].append(ellipse[:]+('r',)) 581 data['rings'].append(np.array(Ring)) 539 582 G2plt.PlotImage(self,newImage=True) 540 583 … … 542 585 data['rings'] = [] 543 586 data['ellipses'] = [] 544 if not data['calibrant'] :587 if not data['calibrant'] or 'None' in data['calibrant']: 545 588 print 'no calibration material selected' 546 589 return True … … 548 591 skip = data['calibskip'] 549 592 dmin = data['calibdmin'] 593 #generate reflection set 550 594 Bravais,Cells = calFile.Calibrants[data['calibrant']][:2] 551 595 HKL = [] … … 556 600 HKL = G2lat.sortHKLd(HKL,True,False) 557 601 wave = data['wavelength'] 558 cent = data['center'] 559 elcent,phi,radii = ellipse 602 #set up 1st ring 603 elcent,phi,radii = ellipse #from fit of 1st ring 560 604 dsp = HKL[0][3] 561 605 tth = 2.0*asind(wave/(2.*dsp)) 562 ttth = tand(tth) 563 data['distance'] = dist = calcDist(radii,tth) 564 radius = dist*tand(tth) 565 zdis,cosB = calcZdisCosB(radius,tth,radii) 566 cent1 = [] 567 cent2 = [] 568 xSum = 0 569 ySum = 0 570 zxSum = 0 571 zySum = 0 572 phiSum = 0 573 tiltSum = 0 574 distSum = 0 575 Zsum = 0 606 ttth = nptand(tth) 607 stth = npsind(tth) 608 ctth = npcosd(tth) 609 #1st estimate of tilt; assume ellipse - don't know sign though 610 tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth) 611 #1st estimate of dist: sample to detector normal to plane 612 data['distance'] = dist = radii[0]**2/(ttth*radii[1]) 613 #ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt 614 zdisp = radii[1]*ttth*tand(tilt) 615 zdism = radii[1]*ttth*tand(-tilt) 616 #cone axis position; 2 choices. Which is right? 617 centp = [elcent[0]+zdisp*sind(phi),elcent[1]+zdisp*cosd(phi)] 618 centm = [elcent[0]+zdism*sind(phi),elcent[1]+zdism*cosd(phi)] 619 #check get same ellipse parms either way 620 #now do next ring; estimate either way & check sum of Imax/Imin in 3x3 block around each point 621 dsp = HKL[1][3] 622 tth = 2.0*asind(wave/(2.*dsp)) 623 ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi) 624 Ringp,delt,sumIp = makeRing(dsp,ellipsep,2,cutoff,scalex,scaley,self.ImageZ) 625 outEp,errp = FitRing(Ringp,True) 626 print '+',ellipsep,errp 627 ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi) 628 Ringm,delt,sumIm = makeRing(dsp,ellipsem,2,cutoff,scalex,scaley,self.ImageZ) 629 outEm,errm = FitRing(Ringm,True) 630 print '-',ellipsem,errm 631 if errp < errm: 632 data['tilt'] = tilt 633 data['center'] = centp 634 negTilt = 1 635 else: 636 data['tilt'] = -tilt 637 data['center'] = centm 638 negTilt = -1 639 data['rotation'] = phi 640 parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1], 641 'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']} 642 varyList = ['dist','det-X','det-Y','tilt','phi'] 643 if data['DetDepthRef']: 644 varyList.append('dep') 576 645 for i,H in enumerate(HKL): 577 646 dsp = H[3] 578 tth = 2.0*asind(0.5*wave/dsp) 579 stth = sind(tth) 580 ctth = cosd(tth) 581 ttth = tand(tth) 582 radius = dist*ttth 583 elcent,phi,radii = ellipse 584 if i: 585 radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2) 586 radii[0] = math.sqrt(radii[1]*radius*cosB) 587 zdis,cosB = calcZdisCosB(radius,tth,radii) 588 zsinp = zdis*sind(phi) 589 zcosp = zdis*cosd(phi) 590 cent = data['center'] 591 elcent = [cent[0]-zsinp,cent[1]+zcosp] 592 ellipse = (elcent,phi,radii) 593 ratio = radii[1]/radii[0] 594 Ring,delt = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) 647 tth = 2.0*asind(wave/(2.*dsp)) 648 print 'HKLD:',H[:4],'2-theta: %.4f'%(tth) 649 elcent,phi,radii = ellipse = GetEllipse(dsp,data) 650 data['ellipses'].append(copy.deepcopy(ellipse+('g',))) 651 print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1]) 652 Ring,delt,sumI = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) 595 653 if Ring: 596 numZ = len(Ring)597 654 data['rings'].append(np.array(Ring)) 598 newellipse = FitRing(Ring,delt) 599 elcent,phi,radii = newellipse 600 if abs(phi) > 45. and phi < 0.: 601 phi += 180. 602 dist = calcDist(radii,tth) 603 distR = 1.-dist/data['distance'] 604 if abs(distR) > 0.1: 605 # print dsp,dist,data['distance'],distR,len(Ring),delt 606 break 607 if distR > 0.001: 608 print 'Wavelength too large?' 609 elif distR < -0.001: 610 print 'Wavelength too small?' 611 else: 612 ellipse = newellipse 613 zdis,cosB = calcZdisCosB(radius,tth,radii) 614 Tilt = acosd(cosB) # 0 <= tilt <= 90 615 zsinp = zdis*sind(ellipse[1]) 616 zcosp = zdis*cosd(ellipse[1]) 617 cent1.append(np.array([elcent[0]+zsinp,elcent[1]-zcosp])) 618 cent2.append(np.array([elcent[0]-zsinp,elcent[1]+zcosp])) 655 rings = np.concatenate((data['rings']),axis=0) 619 656 if i: 620 d1 = cent1[-1]-cent1[-2] #get shift of 2 possible center solutions 621 d2 = cent2[-1]-cent2[-2] 622 if np.dot(d2,d2) > np.dot(d1,d1): #right solution is the larger shift 623 data['center'] = cent1[-1] 624 else: 625 data['center'] = cent2[-1] 626 Zsum += numZ 627 phiSum += numZ*phi 628 distSum += numZ*dist 629 xSum += numZ*data['center'][0] 630 ySum += numZ*data['center'][1] 631 tiltSum += numZ*abs(Tilt) 632 if not np.all(checkEllipse(Zsum,distSum,xSum,ySum,dist,data['center'][0],data['center'][1])): 633 print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i) 634 cent = data['center'] 635 # print ('for ring # %2i @ d-space %.4f: dist %.3f rotate %6.2f tilt %6.2f Xcent %.3f Ycent %.3f Npts %d' 636 # %(i,dsp,dist,phi-90.,Tilt,cent[0],cent[1],numZ)) 657 chi,chisq = FitDetector(rings,varyList,parmDict,False) 658 data['distance'] = parmDict['dist'] 659 data['center'] = [parmDict['det-X'],parmDict['det-Y']] 660 data['rotation'] = np.mod(parmDict['phi'],360.0) 661 data['tilt'] = parmDict['tilt'] 662 data['DetDepth'] = parmDict['dep'] 663 elcent,phi,radii = ellipse = GetEllipse(dsp,data) 664 print fmt2%('fitted ellipse: ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings)) 637 665 data['ellipses'].append(copy.deepcopy(ellipse+('r',))) 638 666 else: 667 print 'insufficient number of points in this ellipse to fit' 639 668 break 640 669 G2plt.PlotImage(self,newImage=True) … … 642 671 if 2*radii[1] < .9*fullSize: 643 672 print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib' 644 if not Zsum: 645 print 'Only one ring fitted. Check your wavelength.' 646 return False 647 data['center'] = [xSum/Zsum,ySum/Zsum] 648 data['distance'] = distSum/Zsum 649 650 #possible error if no. of rings < 3! Might need trap here 651 d1 = cent1[-1]-cent1[1] #compare last ring to 2nd ring 652 d2 = cent2[-1]-cent2[1] 653 Zsign = 1 654 len1 = math.sqrt(np.dot(d1,d1)) 655 len2 = math.sqrt(np.dot(d2,d2)) 656 t1 = d1/len1 657 t2 = d2/len2 658 if len2 > len1: 659 if -135. < atan2d(t2[1],t2[0]) < 45.: 660 Zsign = -1 661 else: 662 if -135. < atan2d(t1[1],t1[0]) < 45.: 663 Zsign = -1 664 665 data['tilt'] = Zsign*tiltSum/Zsum 666 data['rotation'] = phiSum/Zsum 667 varyList = ['dist','det-X','det-Y','tilt','phi'] 668 parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1], 669 'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']} 670 rings = np.concatenate((data['rings']),axis=0) 671 if data['DetDepthRef']: 672 varyList.append('dep') 673 N = len(data['ellipses']) 674 if N > 2: 673 675 FitDetector(rings,varyList,parmDict) 674 data['distance'] = parmDict['dist'] 675 data['center'] = [parmDict['det-X'],parmDict['det-Y']] 676 data['rotation'] = np.mod(parmDict['phi'],360.0) 677 data['tilt'] = parmDict['tilt'] 678 data['DetDepth'] = parmDict['dep'] 679 N = len(data['ellipses']) 680 data['ellipses'] = [] #clear away individual ellipse fits 676 data['distance'] = parmDict['dist'] 677 data['center'] = [parmDict['det-X'],parmDict['det-Y']] 678 data['rotation'] = np.mod(parmDict['phi'],360.0) 679 data['tilt'] = parmDict['tilt'] 680 data['DetDepth'] = parmDict['dep'] 681 681 for H in HKL[:N]: 682 682 ellipse = GetEllipse(H[3],data) … … 837 837 for ring in StrSta['d-zero']: #get observed x,y,d points for the d-zeros 838 838 ellipse = GetEllipse(ring['Dset'],StaControls) 839 Ring,delt = makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)839 Ring,delt,sumI = makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image) 840 840 Ring = np.array(Ring).T 841 841 ring['ImxyObs'] = np.array(Ring[:2]) #need to apply masks to this to eliminate bad points -
trunk/GSASIIimgGUI.py
r1135 r1148 87 87 else: 88 88 G2frame.Integrate = G2img.ImageIntegrate(G2frame.ImageZ,data,masks) 89 G2plt.PlotIntegration(G2frame,newPlot=True)89 # G2plt.PlotIntegration(G2frame,newPlot=True) 90 90 G2IO.SaveIntegration(G2frame,G2frame.PickId,data) 91 91 for item in G2frame.MakePDF: item.Enable(True) … … 203 203 for key in keys: 204 204 if key in ['rotation']: 205 File.write(key+':'+str(data[key] -90.)+'\n')205 File.write(key+':'+str(data[key])+'\n') 206 206 else: 207 207 File.write(key+':'+str(data[key])+'\n') … … 229 229 save[key] = val 230 230 elif key in ['rotation']: 231 save[key] = float(val) +90.231 save[key] = float(val) 232 232 elif key in ['center',]: 233 233 if ',' in val: … … 416 416 calibSizer.Add(wx.StaticText(parent=G2frame.dataDisplay,label=' Tilt rotation'),0, 417 417 wx.ALIGN_CENTER_VERTICAL) 418 rotSel = wx.TextCtrl(parent=G2frame.dataDisplay,value=("%9.3f"%(data['rotation'] -90.)),style=wx.TE_READONLY) #kluge to get rotation from vertical - see GSASIIimage418 rotSel = wx.TextCtrl(parent=G2frame.dataDisplay,value=("%9.3f"%(data['rotation'])),style=wx.TE_READONLY) 419 419 rotSel.SetBackgroundColour(VERY_LIGHT_GREY) 420 420 calibSizer.Add(rotSel,0,wx.ALIGN_CENTER_VERTICAL) -
trunk/GSASIIplot.py
r1134 r1148 2494 2494 xyI = [] 2495 2495 for azm in Azm: 2496 xyI.append(G2img.GetDetectorXY(dspI,azm -90.,Data))2496 xyI.append(G2img.GetDetectorXY(dspI,azm,Data)) 2497 2497 xyI = np.array(xyI) 2498 2498 arcxI,arcyI = xyI.T … … 2501 2501 xyO = [] 2502 2502 for azm in Azm: 2503 xyO.append(G2img.GetDetectorXY(dspO,azm -90.,Data))2503 xyO.append(G2img.GetDetectorXY(dspO,azm,Data)) 2504 2504 xyO = np.array(xyO) 2505 2505 arcxO,arcyO = xyO.T … … 2526 2526 N += 1 2527 2527 for ellipse in Data['ellipses']: 2528 # print 'plot:',ellipse 2528 2529 cent,phi,[width,height],col = ellipse 2529 2530 Plot.add_artist(Ellipse([cent[0],cent[1]],2*width,2*height,phi,ec=col,fc='none')) … … 2650 2651 extent=[ysc[0],ysc[-1],xsc[-1],xsc[0]],aspect='auto') 2651 2652 colorBar = Page.figure.colorbar(Img) 2652 if Data['ellipses']: 2653 for ellipse in Data['ellipses']: 2654 ring = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color 2655 x,y = np.hsplit(ring,2) 2656 tth,azm = G2img.GetTthAzm(x,y,Data) 2657 # azm = np.where(azm < 0.,azm+360,azm) 2658 Plot.plot(tth,azm,'b,') 2653 # if Data['ellipses']: 2654 # for ellipse in Data['ellipses']: 2655 # x,y = np.array(G2img.makeIdealRing(ellipse[:3])) #skip color 2656 # tth,azm = G2img.GetTthAzm(x,y,Data) 2657 ## azm = np.where(azm < 0.,azm+360,azm) 2658 # Plot.plot(tth,azm,'b,') 2659 2659 if not newPlot: 2660 2660 Page.toolbar.push_current() -
trunk/fsource/SConstruct
r771 r1148 62 62 elif GFORTpath != "": 63 63 FCompiler='gfortran' 64 # LDFLAGS = '-static-libgfortran -static-libgcc' 64 65 else: 65 66 print 'No Fortran compiler in path'
Note: See TracChangeset
for help on using the changeset viewer.