Changeset 1153


Ignore:
Timestamp:
Nov 26, 2013 2:32:25 PM (9 years ago)
Author:
vondreele
Message:

fix to multiimage Integrate
fix to calibration
new geometric correction - takes account for tilted detectors

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/GSASIIimage.py

    r1152 r1153  
    258258    amin = 0
    259259    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
    261263        x = radii[0]*cosd(a)
    262264        y = radii[1]*sind(a)
    … …  
    331333   
    332334    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               
    342353    elcent,phi,radii = GetEllipse(dsp,data)
    343354    cent = data['center']
    … …  
    353364    stth = sind(tth)
    354365    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
    362382    return xy
    363383   
    … …  
    384404    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
    385405    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
    387410    dsp = wave/(2.*npsind(tth/2.))
    388411    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
    389     return tth,azm,dsp
     412    return tth,azm,D,dsp
    390413   
    391414def GetTth(x,y,data):
    392     'Needs a doc string'
     415    'Give 2-theta value for detector x,y position; calibration info in data'
    393416    return GetTthAzmDsp(x,y,data)[0]
    394417   
    395418def GetTthAzm(x,y,data):
    396     'Needs a doc string'
     419    'Give 2-theta, azimuth values for detector x,y position; calibration info in data'
    397420    return GetTthAzmDsp(x,y,data)[0:2]
    398421   
     422def 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
    399428def 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]
    402431       
    403432def GetAzm(x,y,data):
    404     'Needs a doc string'
     433    'Give azimuth value for detector x,y position; calibration info in data'
    405434    return GetTthAzmDsp(x,y,data)[1]
    406435       
    … …  
    692721        tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2))
    693722        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.
    695724    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
    696     return np.array(TA),tam           #2-theta & azimuth arrays & position mask
     725    return np.array(TA),tam           #2-theta, azimuth & geom. corr. arrays & position mask
    697726
    698727def Fill2ThetaAzimuthMap(masks,TA,tam,image):
    … …  
    702731    rings = masks['Rings']
    703732    arcs = masks['Arcs']
    704     TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0])))    #azimuth, 2-theta
    705     tax,tay = np.dsplit(TA,2)    #azimuth, 2-theta
     733    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
    706735    for tth,thick in rings:
    707736        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
    … …  
    715744    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
    716745    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
    719748   
    720749def ImageIntegrate(image,data,masks,blkSize=128,dlg=None):
    … …  
    747776            jBeg = jBlk*blkSize
    748777            jFin = min(jBeg+blkSize,Nx)               
    749 #            print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin
    750778            TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
    751779           
    … …  
    755783            Block = image[iBeg:iFin,jBeg:jFin]
    756784            tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
    757             del TA,tam
    758785            Nup += 1
    759786            if dlg:
    … …  
    763790            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
    764791                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,taz
    767792            Nup += 1
    768793            if dlg:
    … …  
    780805    else:
    781806        H1 = LRazm
    782     H0 /= npcosd(H2[:-1])**2
     807#    H0 /= npcosd(H2[:-1])**2
    783808    if data['Oblique'][1]:
    784809        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])

  • trunk/GSASIIimgGUI.py

    r1152 r1153  
    126126                        ifintegrate,name,id = item
    127127                        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
    128135                            dlgp = wx.ProgressDialog("Elapsed time","2D image integration",Nup,
    129136                                style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
    … …  
    132139                                Npix,imagefile = G2frame.PatternTree.GetItemPyData(id)
    133140                                image = G2IO.GetImageData(G2frame,imagefile,True)
    134                                 Id = G2gd.GetPatternTreeItemId(G2frame,id, 'Image Controls')
    135                                 Data = G2frame.PatternTree.GetItemPyData(Id)
    136141                                backImage = []
    137142                                if Data['background image'][0]:
    … …  
    150155                                        G2gd.GetPatternTreeItemId(G2frame,G2frame.Image, 'Masks'),Masks)
    151156                                if len(backImage):                               
    152                                     G2frame.Integrate = G2img.ImageIntegrate(image+backImage,Data,Masks,dlgp)
     157                                    G2frame.Integrate = G2img.ImageIntegrate(image+backImage,Data,Masks,blkSize,dlgp)
    153158                                else:
    154                                     G2frame.Integrate = G2img.ImageIntegrate(image,Data,Masks,dlgp)
     159                                    G2frame.Integrate = G2img.ImageIntegrate(image,Data,Masks,blkSize,dlgp)
    155160#                               G2plt.PlotIntegration(G2frame,newPlot=True,event=event)
    156161                                G2IO.SaveIntegration(G2frame,Id,Data)

  • trunk/GSASIIplot.py

    r1151 r1153  
    21602160                if (0 <= xpix <= sizexy[0]) and (0 <= ypix <= sizexy[1]):
    21612161                    Int = G2frame.ImageZ[ypix][xpix]
    2162                 tth,azm,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
     2162                tth,azm,D,dsp = G2img.GetTthAzmDsp(xpos,ypos,Data)
    21632163                Q = 2.*math.pi/dsp
    21642164                if G2frame.MaskKey in ['p','f']:
    … …  
    23142314                        rings.remove(ring)
    23152315            else:
    2316                 tth,azm,dsp = G2img.GetTthAzmDsp(Xpos,Ypos,Data)
     2316                tth,azm,dsp = G2img.GetTthAzmDsp(Xpos,Ypos,Data)[:3]
    23172317                itemPicked = str(G2frame.itemPicked)
    23182318                if 'Line2D' in itemPicked and PickName == 'Image Controls':
Note: See TracChangeset for help on using the changeset viewer.