Changeset 1165 for trunk/GSASIIimage.py
 Timestamp:
 Dec 14, 2013 8:54:05 AM (9 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/GSASIIimage.py
r1164 r1165 81 81 return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1) 82 82 83 def FitRing(ring,delta): 84 'Needs a doc string' 85 parms = [] 86 if delta: 87 err,parms = FitEllipse(ring) 88 errc,parmsc = FitCircle(ring) 89 errc = errc[0]/(len(ring)*parmsc[2][0]**2) 90 if not parms or errc < .1: 91 parms = parmsc 92 else: 93 err,parms = FitCircle(ring) 94 return parms,err 95 96 def FitCircle(ring): 97 'Needs a doc string' 98 99 def makeParmsCircle(B): 100 cent = [B[0]/2,B[1]/2] 101 phi = 0. 102 sr1 = sr2 = math.sqrt(cent[0]**2+cent[1]**2B[2]) 103 return cent,phi,[sr1,sr2] 104 105 ring = np.array(ring) 106 x = np.asarray(ring.T[0]) 107 y = np.asarray(ring.T[1]) 108 109 M = np.array((x,y,np.ones_like(x))) 110 B = np.array((x**2+y**2)) 111 result = nl.lstsq(M.T,B) 112 return result[1],makeParmsCircle(result[0]) 113 114 def FitEllipse(ring): 115 'Needs a doc string' 116 117 def makeParmsEllipse(B): 118 det = 4.*(1.B[0]**2)B[1]**2 119 if det < 0.: 120 print 'hyperbola!' 121 print B 122 return 0 123 elif det == 0.: 124 print 'parabola!' 125 print B 126 return 0 127 cent = [(B[1]*B[3]2.*(1.B[0])*B[2])/det, \ 128 (B[1]*B[2]2.*(1.+B[0])*B[3])/det] 129 phi = 0.5*atand(0.5*B[1]/B[0]) 130 131 a = (1.+B[0])/cosd(2*phi) 132 b = 2.a 133 f = (1.+B[0])*cent[0]**2+(1.B[0])*cent[1]**2+B[1]*cent[0]*cent[1]B[4] 134 if f/a < 0 or f/b < 0: 135 return 0 136 sr1 = math.sqrt(f/a) 137 sr2 = math.sqrt(f/b) 138 if sr1 > sr2: 139 sr1,sr2 = sr2,sr1 140 phi = 90. 141 if phi < 90.: 142 phi += 180. 143 return cent,phi,[sr1,sr2] 144 145 ring = np.array(ring) 146 x = np.asarray(ring.T[0]) 147 y = np.asarray(ring.T[1]) 148 M = np.array((x**2y**2,x*y,x,y,np.ones_like(x))) 149 B = np.array((x**2+y**2)) 150 bb,err = fel.fellipse(len(x),x,y,1.E7) 151 # print nl.lstsq(M.T,B)[0] 152 # print bb 153 return err,makeParmsEllipse(bb) 154 83 def FitEllipse(xy): 84 85 def ellipse_center(p): 86 ''' gives ellipse center coordinates 87 ''' 88 b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0] 89 num = b*ba*c 90 x0=(c*db*f)/num 91 y0=(a*fb*d)/num 92 return np.array([x0,y0]) 93 94 def ellipse_angle_of_rotation( p ): 95 ''' gives rotation of ellipse major axis from xaxis 96 range will be 90 to 90 deg 97 ''' 98 b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0] 99 return 0.5*npatand(2*b/(ac)) 100 101 def ellipse_axis_length( p ): 102 ''' gives ellipse radii in [minor,major] order 103 ''' 104 b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0] 105 up = 2*(a*f*f+c*d*d+g*b*b2*b*d*fa*c*g) 106 down1=(b*ba*c)*( (ca)*np.sqrt(1+4*b*b/((ac)*(ac)))(c+a)) 107 down2=(b*ba*c)*( (ac)*np.sqrt(1+4*b*b/((ac)*(ac)))(c+a)) 108 res1=np.sqrt(up/down1) 109 res2=np.sqrt(up/down2) 110 return np.array([ res2,res1]) 111 112 xy = np.array(xy) 113 x = np.asarray(xy.T[0])[:,np.newaxis] 114 y = np.asarray(xy.T[1])[:,np.newaxis] 115 D = np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x))) 116 S = np.dot(D.T,D) 117 C = np.zeros([6,6]) 118 C[0,2] = C[2,0] = 2; C[1,1] = 1 119 E, V = nl.eig(np.dot(nl.inv(S), C)) 120 n = np.argmax(np.abs(E)) 121 a = V[:,n] 122 cent = ellipse_center(a) 123 phi = ellipse_angle_of_rotation(a) 124 radii = ellipse_axis_length(a) 125 phi += 90. 126 if radii[0] > radii[1]: 127 radii = [radii[1],radii[0]] 128 phi = 90. 129 return cent,phi,radii 130 155 131 def FitDetector(rings,varyList,parmDict,Print=True): 156 132 'Needs a doc string' … … 163 139 for name,fmt,value,sig in ValSig: 164 140 ptlbls += "%s" % (name.rjust(12)) 165 if name == ' rotate':166 ptstr += fmt % (value 90.) #kluge to get rotation from vertical  see GSASIIimgGUI141 if name == 'phi': 142 ptstr += fmt % (value%360.) 167 143 else: 168 144 ptstr += fmt % (value) … … 179 155 wave = parmDict['wave'] 180 156 if 'dep' in varyList: 181 dist,x0,y0,tilt, phi,dep = B[:6]157 dist,x0,y0,tilt,chi,dep = B[:6] 182 158 else: 183 dist,x0,y0,tilt, phi = B[:5]159 dist,x0,y0,tilt,chi = B[:5] 184 160 dep = parmDict['dep'] 185 161 if 'wave' in varyList: 186 162 wave = B[1] 163 phi = chi90. #get rotation of major axis from tilt axis 187 164 tth = 2.0*npasind(wave/(2.*dsp)) 188 165 dxy = peneCorr(tth,dep) … … 202 179 R1 = (vplus+vminus)/2. #major axis 203 180 zdis = (fplusfminus)/2. 204 X = xx0zdis*npsind(phi) #shift obsd. ellipses to ellipse center= 0,0 205 Y = yy0zdis*npcosd(phi) 206 XR = X*npcosd(phi)+Y*npsind(phi) #rotate by phi to put major axis along x 207 YR = X*npsind(phi)+Y*npcosd(phi) 208 M = (XR/R0)**2+(YR/R1)**21 #=0 for all points exactly on ellipse 181 Robs = np.sqrt((xx0)**2+(yy0)**2) 182 phi0 = npatan2d(yy0,xx0) 183 rsqplus = R0**2+R1**2 184 rsqminus = R0**2R1**2 185 R = rsqminus*npcosd(2.*phi02.*phi)+rsqplus 186 Q = np.sqrt(2.)*R0*R1*np.sqrt(R2.*zdis**2*npsind(phi0phi)**2) 187 P = 2.*R0**2*zdis*npcosd(phi0phi) 188 Rcalc = (P+Q)/R 189 M = RobsRcalc 209 190 return M 210 191 … … 266 247 x = radii[0]*cosd(a) 267 248 y = radii[1]*sind(a) 268 X = (cphi*x sphi*y+cent[0])*scalex #convert mm to pixels269 Y = ( sphi*x+cphi*y+cent[1])*scaley249 X = (cphi*x+sphi*y+cent[0])*scalex #convert mm to pixels 250 Y = (sphi*x+cphi*y+cent[1])*scaley 270 251 X,Y,I,J = ImageLocalMax(image,pix,X,Y) 271 252 if I and J and I/J > reject: … … 279 260 if [X,Y,dsp] not in ring: 280 261 ring.append([X,Y,dsp]) 281 delt = amaxamin282 262 if len(ring) < 10: #want more than 10 deg 283 r eturn [],(delt > 90)284 return ring ,(delt > 90)263 ring = [] 264 return ring 285 265 286 266 def GetEllipse2(tth,dxy,dist,cent,tilt,phi): … … 316 296 radii[1] = eps*(deltf)/(eps**21.) #major axis 317 297 zdis = f+radii[1]*eps 318 elcent = [cent[0]+zdis*sind(phi),cent[1]+zdis*cosd(phi)] 298 #NB: zdis is  to major axis & phi is rotation of minor axis 299 #thus shift from beam to ellipse center is [Z*sin(phi),Z*cos(phi)] 300 elcent = [cent[0]+zdis*sind(phi),cent[1]zdis*cosd(phi)] 319 301 return elcent,phi,radii 320 302 … … 350 332 fplus = dist*tanb*stth/(cosb+stth) 351 333 fminus = dist*tanb*stth/(cosbstth) 352 zdis = 0. 353 if tilt: 354 zdis = tilt*(fplusfminus)/(2.*abs(tilt)) 334 zdis = (fplusfminus)/2. 355 335 rsqplus = radii[0]**2+radii[1]**2 356 336 rsqminus = radii[0]**2radii[1]**2 357 337 R = rsqminus*cosd(2.*azm2.*phi)+rsqplus 358 Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R2.*zdis**2*sind(azm )**2)359 P = zdis*(rsqminus*cosd(azm2.*phi)+rsqplus*cosd(azm))338 Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R2.*zdis**2*sind(azmphi)**2) 339 P = 2.*radii[0]**2*zdis*cosd(azmphi) 360 340 radius = (P+Q)/R 361 341 xy = np.array([radius*cosd(azm),radius*sind(azm)]) … … 363 343 else: #hyperbola  both branches (one is way off screen!) 364 344 f = dist*tanb*stth/(cosb+stth) 365 v = dist*(tanb+tand(tth tilt))345 v = dist*(tanb+tand(tthabs(tilt))) 366 346 delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb)) 367 347 ecc = (vf)/(deltv) 368 348 R = radii[1]*(ecc**21)/(1ecc*cosd(azm)) 369 xy = [R*cosd(azm)f,R*sind(azm)] 349 if tilt > 0.: 350 offset = 2.*radii[1]*ecc+f #select other branch 351 xy = [R*cosd(azm)offset,R*sind(azm)] 352 else: 353 offset = f 354 xy = [R*cosd(azm)offset,R*sind(azm)] 370 355 xy = np.array([xy[0]*cosd(phi)xy[1]*sind(phi),xy[0]*sind(phi)+xy[1]*cosd(phi)]) 371 356 xy += cent … … 525 510 break 526 511 ellipse = GetEllipse(dsp,data) 527 Ring ,delt= makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(self.ImageZ,mask=tam))512 Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(self.ImageZ,mask=tam)) 528 513 if Ring: 529 514 if iH >= skip: … … 572 557 #fit start points on inner ring 573 558 data['ellipses'] = [] 574 outE,err = FitRing(ring,True) 559 data['rings'] = [] 560 outE = FitEllipse(ring) 575 561 fmt = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f' 576 562 fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi^2: %.3f, N: %d' … … 582 568 583 569 #setup 360 points on that ring for "good" fit 584 Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) 570 data['ellipses'].append(ellipse[:]+('g',)) 571 Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) 585 572 if Ring: 586 ellipse ,err = FitRing(Ring,delt)587 Ring ,delt= makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) #do again588 ellipse ,err = FitRing(Ring,delt)573 ellipse = FitEllipse(ring) 574 Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) #do again 575 ellipse = FitEllipse(ring) 589 576 else: 590 577 print '1st ring not sufficiently complete to proceed' 591 578 return False 592 print fmt2%('inner ring: ',ellipse[0][0],ellipse[0][1],ellipse[1],ellipse[2][0],ellipse[2][1], err,len(Ring)) #cent,phi,radii579 print fmt2%('inner ring: ',ellipse[0][0],ellipse[0][1],ellipse[1],ellipse[2][0],ellipse[2][1],0.,len(Ring)) #cent,phi,radii 593 580 data['ellipses'].append(ellipse[:]+('r',)) 594 581 data['rings'].append(np.array(Ring)) 595 582 G2plt.PlotImage(self,newImage=True) 596 583 597 584 #setup for calibration 598 585 data['rings'] = [] 599 data['ellipses'] = []600 586 if not data['calibrant']: 601 587 print 'no calibration material selected' … … 617 603 dsp = HKL[0][3] 618 604 tth = 2.0*asind(wave/(2.*dsp)) 619 Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,self.ImageZ) [0]605 Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,self.ImageZ) 620 606 ttth = nptand(tth) 621 607 stth = npsind(tth) … … 626 612 print 'WARNING  selected ring was fitted as a circle' 627 613 print '  if detector was tilted we suggest you skip this ring  WARNING' 628 tilt = 0.01629 614 #1st estimate of dist: sample to detector normal to plane 630 615 data['distance'] = dist = radii[0]**2/(ttth*radii[1]) … … 632 617 zdisp = radii[1]*ttth*tand(tilt) 633 618 zdism = radii[1]*ttth*tand(tilt) 634 #cone axis position; 2 choices. Which is right? 635 centp = [elcent[0]+zdisp*sind(phi),elcent[1]+zdisp*cosd(phi)] 636 centm = [elcent[0]+zdism*sind(phi),elcent[1]+zdism*cosd(phi)] 619 #cone axis position; 2 choices. Which is right? 620 #NB: zdisp is  to major axis & phi is rotation of minor axis 621 #thus shift from beam to ellipse center is [Z*sin(phi),Z*cos(phi)] 622 centp = [elcent[0]+zdisp*sind(phi),elcent[1]zdisp*cosd(phi)] 623 centm = [elcent[0]+zdism*sind(phi),elcent[1]zdism*cosd(phi)] 637 624 #check get same ellipse parms either way 638 625 #now do next ring; estimate either way & do a FitDetector each way; best fit is correct one … … 640 627 tth = 2.0*asind(wave/(2.*dsp)) 641 628 ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi) 642 ellipsem = GetEllipse2(tth,0.,dist,centm,tilt,phi) 643 Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,self.ImageZ)[0] 644 Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,self.ImageZ)[0] 629 print fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1]) 630 Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,self.ImageZ) 645 631 parmDict = {'dist':dist,'detX':centp[0],'detY':centp[1], 646 632 'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0} 647 633 varyList = ['dist','detX','detY','tilt','phi'] 648 634 if len(Ringp) > 10: 649 chip,chisqp = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,False) 635 chip,chisqp = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True) 636 tiltp = parmDict['tilt'] 637 phip = parmDict['phi'] 638 centp = [parmDict['detX'],parmDict['detY']] 650 639 else: 651 640 chip = 1e6 641 ellipsem = GetEllipse2(tth,0.,dist,centm,tilt,phi) 642 print fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1]) 643 Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,self.ImageZ) 652 644 if len(Ringm) > 10: 653 645 parmDict['tilt'] *= 1 654 chim,chisqm = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,False) 646 chim,chisqm = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True) 647 tiltm = parmDict['tilt'] 648 phim = parmDict['phi'] 649 centm = [parmDict['detX'],parmDict['detY']] 655 650 else: 656 651 chim = 1e6 657 652 if chip < chim: 658 data['tilt'] = tilt 653 data['tilt'] = tiltp 659 654 data['center'] = centp 655 data['rotation'] = phip 660 656 else: 661 data['tilt'] = tilt657 data['tilt'] = tiltm 662 658 data['center'] = centm 663 data['rotation'] = phi 659 data['rotation'] = phim 660 data['ellipses'].append(ellipsep[:]+('b',)) 661 data['rings'].append(np.array(Ringp)) 662 data['ellipses'].append(ellipsem[:]+('r',)) 663 data['rings'].append(np.array(Ringm)) 664 G2plt.PlotImage(self,newImage=True) 664 665 parmDict = {'dist':data['distance'],'detX':data['center'][0],'detY':data['center'][1], 665 666 'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']} … … 668 669 varyList.append('dep') 669 670 data['rings'] = [] 671 data['ellipses'] = [] 670 672 for i,H in enumerate(HKL): 671 673 dsp = H[3] … … 678 680 data['ellipses'].append(copy.deepcopy(ellipse+('g',))) 679 681 print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1]) 680 Ring ,delt= makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)682 Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) 681 683 if Ring: 682 684 data['rings'].append(np.array(Ring)) … … 686 688 data['distance'] = parmDict['dist'] 687 689 data['center'] = [parmDict['detX'],parmDict['detY']] 688 data['rotation'] = np.mod(parmDict['phi']+720.,360.0)690 data['rotation'] = parmDict['phi'] 689 691 data['tilt'] = parmDict['tilt'] 690 if data['rotation'] >180.:691 data['rotation'] = 180.692 data['tilt'] *= 1.693 692 data['DetDepth'] = parmDict['dep'] 694 693 elcent,phi,radii = ellipse = GetEllipse(dsp,data) 695 694 print fmt2%('fitted ellipse: ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings)) 696 695 data['ellipses'].append(copy.deepcopy(ellipse+('r',))) 696 # G2plt.PlotImage(self,newImage=True) 697 697 else: 698 698 print 'insufficient number of points in this ellipse to fit' … … 707 707 data['distance'] = parmDict['dist'] 708 708 data['center'] = [parmDict['detX'],parmDict['detY']] 709 data['rotation'] = np.mod(parmDict['phi']+720.,360.0)709 data['rotation'] = parmDict['phi'] 710 710 data['tilt'] = parmDict['tilt'] 711 if data['rotation'] >180.:712 data['rotation'] = 180.713 data['tilt'] *= 1.714 711 data['DetDepth'] = parmDict['dep'] 715 712 for H in HKL[:N]: … … 863 860 for ring in StrSta['dzero']: #get observed x,y,d points for the dzeros 864 861 ellipse = GetEllipse(ring['Dset'],StaControls) 865 Ring ,delt= makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)862 Ring = makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image) 866 863 Ring = np.array(Ring).T 867 864 ring['ImxyObs'] = np.array(Ring[:2]) #need to apply masks to this to eliminate bad points
Note: See TracChangeset
for help on using the changeset viewer.