Changeset 114


Ignore:
Timestamp:
Jul 16, 2010 10:13:33 AM (11 years ago)
Author:
vondreel
Message:

new integration by blocks

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIimage.py

    r108 r114  
    459459    return True
    460460   
    461 def Make2ThetaAzimuthMap(data,masks,imageN):
     461def Make2ThetaAzimuthMap(data,masks,iLim,jLim):
    462462    import numpy.ma as ma
    463463    import polymask as pm
     
    466466    scalex = pixelSize[0]/1000.
    467467    scaley = pixelSize[1]/1000.
    468     tay,tax = np.mgrid[0.5:imageN+.5,0.5:imageN+.5]         #bin centers not corners
     468   
     469    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
    469470    tax = np.asfarray(tax*scalex,dtype=np.float32)
    470471    tay = np.asfarray(tay*scaley,dtype=np.float32)
     472    nI = iLim[1]-iLim[0]
     473    nJ = jLim[1]-jLim[0]
    471474    #make position masks here
    472475    spots = masks['Points']
    473476    polygons = masks['Polygons']
    474     tam = ma.make_mask_none((imageN,imageN))
    475     if polygons:
    476         print 'Generate polygon mask'
     477    tam = ma.make_mask_none((iLim[1]-iLim[0],jLim[1]-jLim[0]))
    477478    for polygon in polygons:
    478         tamp = ma.make_mask_none((imageN*imageN))
    479         tam = ma.mask_or(tam.flatten(),ma.make_mask(pm.polymask(imageN*imageN,
     479        tamp = ma.make_mask_none((nI*nJ))
     480        tam = ma.mask_or(tam.flatten(),ma.make_mask(pm.polymask(nI*nJ,
    480481            tax.flatten(),tay.flatten(),len(polygon),polygon,tamp)))
    481     if spots:
    482         print 'Generate spot mask'
    483     tam = np.reshape(tam,(imageN,imageN))
     482    if tam.shape: tam = np.reshape(tam,(nI,nJ))
    484483    for X,Y,diam in spots:
    485484        tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2)))
     
    491490    rings = masks['Rings']
    492491    arcs = masks['Arcs']
    493     imageN = len(image)
    494     TA = np.reshape(TA,(2,imageN,imageN))
     492    nJ = len(image)
     493    nI = len(image[0])
     494    TA = np.reshape(TA,(2,nI,nJ))
    495495    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0])))    #azimuth, 2-theta
    496496    tax,tay = np.dsplit(TA,2)    #azimuth, 2-theta
    497     if rings:
    498         print 'Generate ring mask'
    499497    for tth,thick in rings:
    500498        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
    501     if arcs:
    502         print 'Generate arc mask'
    503499    for tth,azm,thick in arcs:
    504500        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))* \
     
    512508    return tax,tay,taz
    513509   
    514 def Bin2ThetaAzimuthMap(data,tax,tay,taz):
    515     import numpy.ma as ma
     510def ImageIntegrate(self,data,masks):
    516511    import histogram2d as h2d
     512    print 'Begin image integration'
    517513    LUtth = data['IOtth']
    518514    if data['fullIntegrate']:
     
    522518    numAzms = data['outAzimuths']
    523519    numChans = data['outChannels']
    524     t0 = time.time()
     520    Dtth = (LUtth[1]-LUtth[0])/numChans
     521    Dazm = (LRazm[1]-LRazm[0])/numAzms
    525522    NST = np.zeros(shape=(numAzms,numChans),dtype=np.int,order='F')
    526523    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
    527     H1 = np.zeros(shape=(numAzms+1,))
    528     H2 = np.zeros(shape=(numChans+1,))   
    529     HST = [H0,H1,H2]
    530     NST,H0,H1,H2 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,NST,H0,H1,H2)
    531     return NST,HST
    532 
    533 def ImageIntegrate(self,data,masks):
    534     dlg = wx.ProgressDialog("Elapsed time","2D image integration",5,
     524    imageN = len(self.ImageZ)
     525    nBlks = (imageN-1)/1024+1
     526    dlg = wx.ProgressDialog("Elapsed time","2D image integration",nBlks*nBlks*3+3,
    535527        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
    536528    try:
    537         print 'Begin image integration'
    538         print 'Create 2-theta,azimuth map'
    539529        t0 = time.time()
    540         dlg.Update(0)
    541         imageN = len(self.ImageZ)
    542         TA,tam = Make2ThetaAzimuthMap(data,masks,imageN)           #2-theta & azimuth arrays & create position mask
    543         dlg.Update(1)
    544         print 'Fill map with 2-theta/azimuth values'
    545         tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,self.ImageZ)    #and apply masks
    546         del TA
    547         dlg.Update(2)
    548         print 'Bin image by 2-theta/azimuth intervals'
    549         NST,HST = Bin2ThetaAzimuthMap(data,tax,tay,taz)
     530        Nup = 0
     531        dlg.Update(Nup)
     532        for iBlk in range(nBlks):
     533            iBeg = iBlk*1024
     534            iFin = min(iBeg+1024,imageN)
     535            for jBlk in range(nBlks):
     536                jBeg = jBlk*1024
     537                jFin = min(jBeg+1024,imageN)               
     538                print 'Process map block',iBlk,jBlk,iBeg,iFin,jBeg,jFin
     539                TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
     540                Nup += 1
     541                dlg.Update(Nup)
     542                Block = self.ImageZ[iBeg:iFin,jBeg:jFin]
     543                tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
     544                del TA
     545                Nup += 1
     546                dlg.Update(Nup)
     547                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
     548                Nup += 1
     549                dlg.Update(Nup)
     550        H0 = np.nan_to_num(np.divide(H0,NST))
     551        del NST
     552        if Dtth:
     553            H2 = [tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)]
     554        else:
     555            H2 = LUtth
     556        if Dazm:       
     557            H1 = [azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)]
     558        else:
     559            H1 = LRazm
     560        self.Integrate = [H0,H1,H2]
    550561        print 'Binning complete'
    551         del tax,tay,taz
    552         dlg.Update(3)
    553         print 'Form normalized 1D pattern(s)'
    554         HST[0] = np.divide(HST[0],NST)
    555         print 'Normalize done'
    556         self.Integrate = HST
    557         del NST
    558         dlg.Update(4)
     562        Nup += 1
     563        dlg.Update(Nup)
    559564        t1 = time.time()
    560565        print "Elapsed time:","%8.3f"%(t1-t0), "s"
Note: See TracChangeset for help on using the changeset viewer.