Changeset 1153
- Timestamp:
- Nov 26, 2013 2:32:25 PM (9 years ago)
- Location:
- trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIimage.py
r1152 r1153 258 258 amin = 0 259 259 amax = 360 260 for a in range(0,int(ellipseC()),1): 260 C = int(ellipseC()) 261 for i in range(0,C,1): 262 a = 360.*i/C 261 263 x = radii[0]*cosd(a) 262 264 y = radii[1]*sind(a) … … 331 333 332 334 from scipy.optimize import fsolve 333 def func(xy,*args): 334 azm,phi,R0,R1,A,B = args 335 cp = cosd(phi) 336 sp = sind(phi) 337 x,y = xy 338 out = [] 339 out.append(y-x*tand(azm)) 340 out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2) 341 return out 335 def ellip(xy,*args): 336 azm,phi,R0,R1,A,B = args 337 cp = cosd(phi) 338 sp = sind(phi) 339 x,y = xy 340 out = [] 341 out.append(y-x*tand(azm)) 342 out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2) 343 return out 344 def hyperb(xy,*args): 345 azm,phi,R0,R1,A,B = args 346 cp = cosd(phi) 347 sp = sind(phi) 348 x,y = xy 349 out = [] 350 out.append(y+x*tand(azm)) 351 out.append(R0**2*((x+A)*sp-(y+B)*cp)**2-R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2) 352 return out 342 353 elcent,phi,radii = GetEllipse(dsp,data) 343 354 cent = data['center'] … … 353 364 stth = sind(tth) 354 365 cosb = cosd(tilt) 355 R1 = (dist+dxy)*stth*cosd(tth)*cosb/(cosb**2-stth**2) 356 R0 = math.sqrt(R1*radius*cosb) 357 zdis = R1*ttth*tand(tilt) 358 A = zdis*sind(phi) 359 B = -zdis*cosd(phi) 360 xy0 = [radius*cosd(azm),radius*sind(azm)] 361 xy = fsolve(func,xy0,args=(azm,phi,R0,R1,A,B))+cent 366 if tth+abs(tilt) < 90.: 367 R1 = (dist+dxy)*stth*cosd(tth)*cosb/(cosb**2-stth**2) 368 R0 = math.sqrt(R1*radius*cosb) 369 zdis = R1*ttth*tand(tilt) 370 A = zdis*sind(phi) 371 B = -zdis*cosd(phi) 372 xy0 = [radius*cosd(azm),radius*sind(azm)] 373 xy = fsolve(ellip,xy0,args=(azm,phi,R0,R1,A,B))+cent 374 else: 375 R1 = (dist+dxy)*stth*cosd(tth)*cosb/(cosb**2+stth**2) 376 R0 = math.sqrt(R1*radius*cosb) 377 zdis = R1*ttth*tand(tilt) 378 A = zdis*sind(phi) 379 B = -zdis*cosd(phi) 380 xy0 = [radius*cosd(azm),radius*sind(azm)] 381 xy = fsolve(hyperb,xy0,args=(azm,phi,R0,R1,A,B))+cent 362 382 return xy 363 383 … … 384 404 tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z)) 385 405 dxy = peneCorr(tth,dep) 386 tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z+dxy)) #depth corr not correct for tilted detector 406 DX = dist-Z+dxy 407 DY = np.sqrt(dx**2+dy**2-Z**2) 408 D = (DX**2+DY**2)/dist**2 #for geometric correction = 1/cos(2theta)^2 if tilt=0. 409 tth = npatand(DY/DX) #depth corr not correct for tilted detector 387 410 dsp = wave/(2.*npsind(tth/2.)) 388 411 azm = (npatan2d(dy,dx)+azmthoff+720.)%360. 389 return tth,azm, dsp412 return tth,azm,D,dsp 390 413 391 414 def GetTth(x,y,data): 392 ' Needs a doc string'415 'Give 2-theta value for detector x,y position; calibration info in data' 393 416 return GetTthAzmDsp(x,y,data)[0] 394 417 395 418 def GetTthAzm(x,y,data): 396 ' Needs a doc string'419 'Give 2-theta, azimuth values for detector x,y position; calibration info in data' 397 420 return GetTthAzmDsp(x,y,data)[0:2] 398 421 422 def GetTthAzmD(x,y,data): 423 '''Give 2-theta, azimuth & geometric correction values for detector x,y position; 424 calibration info in data 425 ''' 426 return GetTthAzmDsp(x,y,data)[0:3] 427 399 428 def GetDsp(x,y,data): 400 ' Needs a doc string'401 return GetTthAzmDsp(x,y,data)[ 2]429 'Give d-spacing value for detector x,y position; calibration info in data' 430 return GetTthAzmDsp(x,y,data)[3] 402 431 403 432 def GetAzm(x,y,data): 404 ' Needs a doc string'433 'Give azimuth value for detector x,y position; calibration info in data' 405 434 return GetTthAzmDsp(x,y,data)[1] 406 435 … … 692 721 tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2)) 693 722 tam = ma.mask_or(tam,tamp) 694 TA = np.array(GetTthAzm (tax,tay,data))723 TA = np.array(GetTthAzmD(tax,tay,data)) #includes geom. corr. 695 724 TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1]) 696 return np.array(TA),tam #2-theta & azimutharrays & position mask725 return np.array(TA),tam #2-theta, azimuth & geom. corr. arrays & position mask 697 726 698 727 def Fill2ThetaAzimuthMap(masks,TA,tam,image): … … 702 731 rings = masks['Rings'] 703 732 arcs = masks['Arcs'] 704 TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]) )) #azimuth, 2-theta705 tax,tay = np.dsplit(TA,2) #azimuth, 2-theta733 TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2]))) #azimuth, 2-theta, dist 734 tax,tay,tad = np.dsplit(TA,3) #azimuth, 2-theta, dist 706 735 for tth,thick in rings: 707 736 tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))) … … 715 744 tay = ma.compressed(ma.array(tay.flatten(),mask=tam)) 716 745 taz = ma.compressed(ma.array(taz.flatten(),mask=tam)) 717 del(tam)718 return tax,tay,taz 746 tad = ma.compressed(ma.array(tad.flatten(),mask=tam)) 747 return tax,tay,taz*tad**2 719 748 720 749 def ImageIntegrate(image,data,masks,blkSize=128,dlg=None): … … 747 776 jBeg = jBlk*blkSize 748 777 jFin = min(jBeg+blkSize,Nx) 749 # print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin750 778 TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin)) #2-theta & azimuth arrays & create position mask 751 779 … … 755 783 Block = image[iBeg:iFin,jBeg:jFin] 756 784 tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block) #and apply masks 757 del TA,tam758 785 Nup += 1 759 786 if dlg: … … 763 790 if any([tax.shape[0],tay.shape[0],taz.shape[0]]): 764 791 NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0) 765 # print 'block done'766 del tax,tay,taz767 792 Nup += 1 768 793 if dlg: … … 780 805 else: 781 806 H1 = LRazm 782 H0 /= npcosd(H2[:-1])**2807 # H0 /= npcosd(H2[:-1])**2 783 808 if data['Oblique'][1]: 784 809 H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1]) -
trunk/GSASIIimgGUI.py
r1152 r1153 126 126 ifintegrate,name,id = item 127 127 if ifintegrate: 128 Id = G2gd.GetPatternTreeItemId(G2frame,id, 'Image Controls') 129 Data = G2frame.PatternTree.GetItemPyData(Id) 130 blkSize = 128 #this seems to be optimal; will break in polymask if >1024 131 Nx,Ny = Data['size'] 132 nXBlks = (Nx-1)/blkSize+1 133 nYBlks = (Ny-1)/blkSize+1 134 Nup = nXBlks*nYBlks*3+3 128 135 dlgp = wx.ProgressDialog("Elapsed time","2D image integration",Nup, 129 136 style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE) … … 132 139 Npix,imagefile = G2frame.PatternTree.GetItemPyData(id) 133 140 image = G2IO.GetImageData(G2frame,imagefile,True) 134 Id = G2gd.GetPatternTreeItemId(G2frame,id, 'Image Controls')135 Data = G2frame.PatternTree.GetItemPyData(Id)136 141 backImage = [] 137 142 if Data['background image'][0]: … … 150 155 G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'),Masks) 151 156 if len(backImage): 152 G2frame.Integrate = G2img.ImageIntegrate(image+backImage,Data,Masks, dlgp)157 G2frame.Integrate = G2img.ImageIntegrate(image+backImage,Data,Masks,blkSize,dlgp) 153 158 else: 154 G2frame.Integrate = G2img.ImageIntegrate(image,Data,Masks, dlgp)159 G2frame.Integrate = G2img.ImageIntegrate(image,Data,Masks,blkSize,dlgp) 155 160 # G2plt.PlotIntegration(G2frame,newPlot=True,event=event) 156 161 G2IO.SaveIntegration(G2frame,Id,Data) -
trunk/GSASIIplot.py
r1151 r1153 2160 2160 if (0 <= xpix <= sizexy[0]) and (0 <= ypix <= sizexy[1]): 2161 2161 Int = G2frame.ImageZ[ypix][xpix] 2162 tth,azm, dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)2162 tth,azm,D,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data) 2163 2163 Q = 2.*math.pi/dsp 2164 2164 if G2frame.MaskKey in ['p','f']: … … 2314 2314 rings.remove(ring) 2315 2315 else: 2316 tth,azm,dsp = G2img.GetTthAzmDsp(Xpos,Ypos,Data) 2316 tth,azm,dsp = G2img.GetTthAzmDsp(Xpos,Ypos,Data)[:3] 2317 2317 itemPicked = str(G2frame.itemPicked) 2318 2318 if 'Line2D' in itemPicked and PickName == 'Image Controls':
Note: See TracChangeset
for help on using the changeset viewer.