Changeset 1149
- Timestamp:
- Nov 24, 2013 2:50:12 PM (9 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
-
GSASIIimage.py (modified) (8 diffs)
-
GSASIIimgGUI.py (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/GSASIIimage.py
r1148 r1149 14 14 Ellipse fitting & image integration 15 15 16 needs some minor refactoring to remove wx code17 actually not easy in this case wx.ProgressDialog needs #of blocks to process when started18 not known in G2imageGUI.19 16 ''' 20 17 21 18 import math 22 import wx23 19 import time 24 20 import numpy as np … … 510 506 'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']} 511 507 Found = False 508 wave = data['wavelength'] 512 509 for H in HKL: 513 510 dsp = H[3] 511 tth = 2.0*asind(wave/(2.*dsp)) 512 if tth+abs(data['tilt']) > 90.: 513 print 'next line is a hyperbola - search stopped' 514 break 514 515 ellipse = GetEllipse(dsp,data) 515 516 Ring,delt,sumI = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ) … … 585 586 data['rings'] = [] 586 587 data['ellipses'] = [] 587 if not data['calibrant'] or 'None' in data['calibrant']:588 if not data['calibrant']: 588 589 print 'no calibration material selected' 589 590 return True … … 609 610 #1st estimate of tilt; assume ellipse - don't know sign though 610 611 tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth) 612 if not tilt: 613 print 'WARNING - selected ring was fitted as a circle' 614 print ' - if detector was tilted we suggest you skip this ring - WARNING' 611 615 #1st estimate of dist: sample to detector normal to plane 612 616 data['distance'] = dist = radii[0]**2/(ttth*radii[1]) … … 623 627 ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi) 624 628 Ringp,delt,sumIp = makeRing(dsp,ellipsep,2,cutoff,scalex,scaley,self.ImageZ) 625 outEp,errp = FitRing(Ringp,True) 626 print '+',ellipsep,errp 629 if len(Ringp) > 10: 630 outEp,errp = FitRing(Ringp,True) 631 else: 632 errp = 1e6 633 # print '+',ellipsep,errp 627 634 ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi) 628 635 Ringm,delt,sumIm = makeRing(dsp,ellipsem,2,cutoff,scalex,scaley,self.ImageZ) 629 outEm,errm = FitRing(Ringm,True) 630 print '-',ellipsem,errm 636 if len(Ringm) > 10: 637 outEm,errm = FitRing(Ringm,True) 638 else: 639 errm = 1e6 640 # print '-',ellipsem,errm 631 641 if errp < errm: 632 642 data['tilt'] = tilt … … 646 656 dsp = H[3] 647 657 tth = 2.0*asind(wave/(2.*dsp)) 658 if tth+abs(data['tilt']) > 90.: 659 print 'next line is a hyperbola - search stopped' 660 break 648 661 print 'HKLD:',H[:4],'2-theta: %.4f'%(tth) 649 662 elcent,phi,radii = ellipse = GetEllipse(dsp,data) … … 746 759 return tax,tay,taz 747 760 748 def ImageIntegrate(image,data,masks ):761 def ImageIntegrate(image,data,masks,dlg=None): 749 762 'Needs a doc string' 750 763 import histogram2d as h2d … … 766 779 nYBlks = (Ny-1)/blkSize+1 767 780 Nup = nXBlks*nYBlks*3+3 768 dlg = wx.ProgressDialog("Elapsed time","2D image integration",Nup, 769 style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE) 770 try: 771 t0 = time.time() 772 Nup = 0 781 t0 = time.time() 782 Nup = 0 783 if dlg: 773 784 dlg.Update(Nup) 774 for iBlk in range(nYBlks): 775 iBeg = iBlk*blkSize 776 iFin = min(iBeg+blkSize,Ny) 777 for jBlk in range(nXBlks): 778 jBeg = jBlk*blkSize 779 jFin = min(jBeg+blkSize,Nx) 780 # print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin 781 TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin)) #2-theta & azimuth arrays & create position mask 782 783 Nup += 1 785 for iBlk in range(nYBlks): 786 iBeg = iBlk*blkSize 787 iFin = min(iBeg+blkSize,Ny) 788 for jBlk in range(nXBlks): 789 jBeg = jBlk*blkSize 790 jFin = min(jBeg+blkSize,Nx) 791 # print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin 792 TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin)) #2-theta & azimuth arrays & create position mask 793 794 Nup += 1 795 if dlg: 784 796 dlg.Update(Nup) 785 Block = image[iBeg:iFin,jBeg:jFin] 786 tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block) #and apply masks 787 del TA,tam 788 Nup += 1 797 Block = image[iBeg:iFin,jBeg:jFin] 798 tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block) #and apply masks 799 del TA,tam 800 Nup += 1 801 if dlg: 789 802 dlg.Update(Nup) 790 tax = np.where(tax > LRazm[1],tax-360.,tax) #put azm inside limits if possible 791 tax = np.where(tax < LRazm[0],tax+360.,tax) 792 if any([tax.shape[0],tay.shape[0],taz.shape[0]]): 793 NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0) 794 # print 'block done' 795 del tax,tay,taz 796 Nup += 1 803 tax = np.where(tax > LRazm[1],tax-360.,tax) #put azm inside limits if possible 804 tax = np.where(tax < LRazm[0],tax+360.,tax) 805 if any([tax.shape[0],tay.shape[0],taz.shape[0]]): 806 NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0) 807 # print 'block done' 808 del tax,tay,taz 809 Nup += 1 810 if dlg: 797 811 dlg.Update(Nup) 798 NST = np.array(NST) 799 H0 = np.divide(H0,NST) 800 H0 = np.nan_to_num(H0) 801 del NST 802 if Dtth: 803 H2 = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)]) 804 else: 805 H2 = np.array(LUtth) 806 if Dazm: 807 H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)]) 808 else: 809 H1 = LRazm 810 H0 /= npcosd(H2[:-1])**2 811 if data['Oblique'][1]: 812 H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1]) 813 if 'SASD' in data['type'] and data['PolaVal'][1]: 814 H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm)[0] for azm in H1[:-1]+np.diff(H1)/2.]) 815 Nup += 1 812 NST = np.array(NST) 813 H0 = np.divide(H0,NST) 814 H0 = np.nan_to_num(H0) 815 del NST 816 if Dtth: 817 H2 = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)]) 818 else: 819 H2 = np.array(LUtth) 820 if Dazm: 821 H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)]) 822 else: 823 H1 = LRazm 824 H0 /= npcosd(H2[:-1])**2 825 if data['Oblique'][1]: 826 H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1]) 827 if 'SASD' in data['type'] and data['PolaVal'][1]: 828 H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm)[0] for azm in H1[:-1]+np.diff(H1)/2.]) 829 Nup += 1 830 if dlg: 816 831 dlg.Update(Nup) 817 t1 = time.time() 818 finally: 819 dlg.Destroy() 832 t1 = time.time() 820 833 print 'Integration complete' 821 834 print "Elapsed time:","%8.3f"%(t1-t0), "s" -
trunk/GSASIIimgGUI.py
r1148 r1149 73 73 74 74 def OnIntegrate(event): 75 if data['background image'][0]: 76 maskCopy = copy.deepcopy(masks) 77 backImg = data['background image'][0] 78 backScale = data['background image'][1] 79 id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, backImg) 80 Npix,imagefile = G2frame.PatternTree.GetItemPyData(id) 81 backImage = G2IO.GetImageData(G2frame,imagefile,True)*backScale 82 sumImage = G2frame.ImageZ+backImage 83 sumMin = np.min(sumImage) 84 sumMax = np.max(sumImage) 85 maskCopy['Thresholds'] = [(sumMin,sumMax),[sumMin,sumMax]] 86 G2frame.Integrate = G2img.ImageIntegrate(sumImage,data,maskCopy) 87 else: 88 G2frame.Integrate = G2img.ImageIntegrate(G2frame.ImageZ,data,masks) 89 # G2plt.PlotIntegration(G2frame,newPlot=True) 90 G2IO.SaveIntegration(G2frame,G2frame.PickId,data) 75 dlg = wx.ProgressDialog("Elapsed time","2D image integration",Nup, 76 style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE) 77 try: 78 if data['background image'][0]: 79 maskCopy = copy.deepcopy(masks) 80 backImg = data['background image'][0] 81 backScale = data['background image'][1] 82 id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, backImg) 83 Npix,imagefile = G2frame.PatternTree.GetItemPyData(id) 84 backImage = G2IO.GetImageData(G2frame,imagefile,True)*backScale 85 sumImage = G2frame.ImageZ+backImage 86 sumMin = np.min(sumImage) 87 sumMax = np.max(sumImage) 88 maskCopy['Thresholds'] = [(sumMin,sumMax),[sumMin,sumMax]] 89 G2frame.Integrate = G2img.ImageIntegrate(sumImage,data,maskCopy,dlg) 90 else: 91 G2frame.Integrate = G2img.ImageIntegrate(G2frame.ImageZ,data,masks,dlg) 92 # G2plt.PlotIntegration(G2frame,newPlot=True) 93 G2IO.SaveIntegration(G2frame,G2frame.PickId,data) 94 finally: 95 dlg.Destroy() 91 96 for item in G2frame.MakePDF: item.Enable(True) 92 97 … … 116 121 ifintegrate,name,id = item 117 122 if ifintegrate: 118 id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, name) 119 Npix,imagefile = G2frame.PatternTree.GetItemPyData(id) 120 image = G2IO.GetImageData(G2frame,imagefile,True) 121 Id = G2gd.GetPatternTreeItemId(G2frame,id, 'Image Controls') 122 Data = G2frame.PatternTree.GetItemPyData(Id) 123 backImage = [] 124 if Data['background image'][0]: 125 backImg = Data['background image'][0] 126 backScale = Data['background image'][1] 127 id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, backImg) 123 dlgp = wx.ProgressDialog("Elapsed time","2D image integration",Nup, 124 style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE) 125 try: 126 id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, name) 128 127 Npix,imagefile = G2frame.PatternTree.GetItemPyData(id) 129 backImage = G2IO.GetImageData(G2frame,imagefile,True)*backScale 130 try: 131 Masks = G2frame.PatternTree.GetItemPyData( 132 G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks')) 133 except TypeError: #missing Masks 134 Imin,Imax = Data['Range'] 135 Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Frames':[],'Thresholds':[(Imin,Imax),[Imin,Imax]]} 136 G2frame.PatternTree.SetItemPyData( 137 G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'),Masks) 138 if len(backImage): 139 G2frame.Integrate = G2img.ImageIntegrate(image+backImage,Data,Masks) 140 else: 141 G2frame.Integrate = G2img.ImageIntegrate(image,Data,Masks) 142 # G2plt.PlotIntegration(G2frame,newPlot=True,event=event) 143 G2IO.SaveIntegration(G2frame,Id,Data) 128 image = G2IO.GetImageData(G2frame,imagefile,True) 129 Id = G2gd.GetPatternTreeItemId(G2frame,id, 'Image Controls') 130 Data = G2frame.PatternTree.GetItemPyData(Id) 131 backImage = [] 132 if Data['background image'][0]: 133 backImg = Data['background image'][0] 134 backScale = Data['background image'][1] 135 id = G2gd.GetPatternTreeItemId(G2frame, G2frame.root, backImg) 136 Npix,imagefile = G2frame.PatternTree.GetItemPyData(id) 137 backImage = G2IO.GetImageData(G2frame,imagefile,True)*backScale 138 try: 139 Masks = G2frame.PatternTree.GetItemPyData( 140 G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks')) 141 except TypeError: #missing Masks 142 Imin,Imax = Data['Range'] 143 Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Frames':[],'Thresholds':[(Imin,Imax),[Imin,Imax]]} 144 G2frame.PatternTree.SetItemPyData( 145 G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'),Masks) 146 if len(backImage): 147 G2frame.Integrate = G2img.ImageIntegrate(image+backImage,Data,Masks,dlgp) 148 else: 149 G2frame.Integrate = G2img.ImageIntegrate(image,Data,Masks,dlgp) 150 # G2plt.PlotIntegration(G2frame,newPlot=True,event=event) 151 G2IO.SaveIntegration(G2frame,Id,Data) 152 finally: 153 dlgp.Destroy() 144 154 finally: 145 155 dlg.Destroy()
Note: See TracChangeset
for help on using the changeset viewer.
