Changeset 3104


Ignore:
Timestamp:
Sep 29, 2017 12:25:46 PM (4 years ago)
Author:
vondreele
Message:

achieve some speed (~10%) improvement in integration of images

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIimage.py

    r2836 r3104  
    5050npatand = lambda x: 180.*np.arctan(x)/np.pi
    5151npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
     52nxs = np.newaxis
    5253debug = False
    5354   
     
    813814    return True
    814815   
    815 def Make2ThetaAzimuthMap(data,masks,iLim,jLim,times): #most expensive part of integration!
     816def Make2ThetaAzimuthMap(data,masks,iLim,jLim,tamp,times): #most expensive part of integration!
    816817    'Needs a doc string'
    817818    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
     
    821822   
    822823    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
    823     tax = np.asfarray(tax*scalex,dtype=np.float32)
    824     tay = np.asfarray(tay*scaley,dtype=np.float32)
     824    tax = np.asfarray(tax*scalex,dtype=np.float32).flatten()
     825    tay = np.asfarray(tay*scaley,dtype=np.float32).flatten()
    825826    nI = iLim[1]-iLim[0]
    826827    nJ = jLim[1]-jLim[0]
     
    828829    #make position masks here
    829830    frame = masks['Frames']
    830     tam = ma.make_mask_none((nI,nJ))
     831    tam = ma.make_mask_none((nI*nJ))
    831832    if frame:
    832         tamp = ma.make_mask_none((1024*1024))
    833         tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
    834             tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
    835         tam = ma.mask_or(tam.flatten(),tamp)
     833        tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
     834            tay,len(frame),frame,tamp)[:nI*nJ])-True)
    836835    polygons = masks['Polygons']
    837836    for polygon in polygons:
    838837        if polygon:
    839             tamp = ma.make_mask_none((1024*1024))
    840             tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
    841                 tay.flatten(),len(polygon),polygon,tamp)[:nI*nJ])
    842             tam = ma.mask_or(tam.flatten(),tamp)
     838            tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
     839                tay,len(polygon),polygon,tamp)[:nI*nJ]))
     840    for X,Y,rsq in masks['Points'].T:
     841        tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,rsq)))
    843842    if tam.shape: tam = np.reshape(tam,(nI,nJ))
    844     spots = masks['Points']
    845     for X,Y,diam in spots:
    846         tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2))
    847         tam = ma.mask_or(tam,tamp)
    848843    times[0] += time.time()-t0
    849844    t0 = time.time()
    850     TA = np.array(GetTthAzmG(tax,tay,data))     #includes geom. corr. as dist**2/d0**2 - most expensive step
     845    TA = np.array(GetTthAzmG(np.reshape(tax,(nI,nJ)),np.reshape(tay,(nI,nJ)),data))     #includes geom. corr. as dist**2/d0**2 - most expensive step
    851846    times[1] += time.time()-t0
    852847    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
     
    875870    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
    876871    return tax,tay,taz,tad,tabs
    877    
     872
    878873def ImageIntegrate(image,data,masks,blkSize=128,returnN=False):
    879874    'Integrate an image; called from OnIntegrateAll and OnIntegrate in G2imgGUI'    #for q, log(q) bins need data['binType']
     
    898893    if 'SASD' in data['type']:
    899894        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
     895    Masks = copy.deepcopy(masks)
     896    Masks['Points'] = np.array(Masks['Points']).T           #get spots as X,Y,R arrays
     897    if np.any(masks['Points']):
     898        Masks['Points'][2] = np.square(Masks['Points'][2]/2.)
    900899    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
    901900    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
     
    905904    tbeg = time.time()
    906905    times = [0,0,0,0,0]
     906    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
    907907    for iBlk in range(nYBlks):
    908908        iBeg = iBlk*blkSize
     
    912912            jFin = min(jBeg+blkSize,Nx)
    913913            # next is most expensive step!
    914             TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin),times)           #2-theta & azimuth arrays & create position mask
     914            TA,tam = Make2ThetaAzimuthMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp,times)           #2-theta & azimuth arrays & create position mask
    915915            Block = image[iBeg:iFin,jBeg:jFin]
    916916            t0 = time.time()
    917             tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
     917            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(Masks,TA,tam,Block)    #and apply masks
    918918            del TA; del tam
    919919            times[2] += time.time()-t0
Note: See TracChangeset for help on using the changeset viewer.