Changeset 3186


Ignore:
Timestamp:
Dec 9, 2017 8:54:59 AM (4 years ago)
Author:
vondreele
Message:

revisions to integrate images so that masks are reused if not changed in integrate all & auto integrate

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIimage.py

    r3173 r3186  
    2323import polymask as pm
    2424from scipy.optimize import leastsq
    25 import scipy.signal as scsg
    26 import scipy.cluster.vq as scvq
    2725import scipy.interpolate as scint
    2826import copy
     
    237235        'compute estimate of ellipse circumference'
    238236        if radii[0] < 0:        #hyperbola
    239             theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
     237#            theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
    240238#            print (theta)
    241239            return 0
     
    815813    return True
    816814   
    817 def Make2ThetaAzimuthMap(data,iLim,jLim,times): #most expensive part of integration!
     815def Make2ThetaAzimuthMap(data,iLim,jLim): #most expensive part of integration!
    818816    'Needs a doc string'
    819817    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
     
    827825    nI = iLim[1]-iLim[0]
    828826    nJ = jLim[1]-jLim[0]
    829     t0 = time.time()
    830827    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
    831     times[1] += time.time()-t0
    832828    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
    833829    return TA           #2-theta, azimuth & geom. corr. arrays
    834830
    835 def MakeMaskMap(data,masks,iLim,jLim,tamp,times):
     831def MakeMaskMap(data,masks,iLim,jLim,tamp):
    836832    pixelSize = data['pixelSize']
    837833    scalex = pixelSize[0]/1000.
     
    843839    nI = iLim[1]-iLim[0]
    844840    nJ = jLim[1]-jLim[0]
    845     t0 = time.time()
    846841    #make position masks here
    847842    frame = masks['Frames']
     
    859854            tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,rsq)))
    860855    if tam.shape: tam = np.reshape(tam,(nI,nJ))
    861     times[0] += time.time()-t0
    862856    return tam           #position mask
    863857
     
    886880
    887881def MakeUseTA(data,blkSize=128):
    888     times = [0,0,0,0,0]
    889882    Nx,Ny = data['size']
    890883    nXBlks = (Nx-1)//blkSize+1
     
    898891            jBeg = jBlk*blkSize
    899892            jFin = min(jBeg+blkSize,Nx)
    900             TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin),times)          #2-theta & azimuth arrays & create position mask
     893            TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))          #2-theta & azimuth arrays & create position mask
    901894            useTAj.append(TA)
    902895        useTA.append(useTAj)
    903896    return useTA
    904897
    905 def ImageIntegrate(image,data,masks,blkSize=128,returnN=False,useTA=None):
     898def MakeUseMask(data,masks,blkSize=128):
     899    Masks = copy.deepcopy(masks)
     900    Masks['Points'] = np.array(Masks['Points']).T           #get spots as X,Y,R arrays
     901    if np.any(masks['Points']):
     902        Masks['Points'][2] = np.square(Masks['Points'][2]/2.)
     903    Nx,Ny = data['size']
     904    nXBlks = (Nx-1)//blkSize+1
     905    nYBlks = (Ny-1)//blkSize+1
     906    useMask = []
     907    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
     908    for iBlk in range(nYBlks):
     909        iBeg = iBlk*blkSize
     910        iFin = min(iBeg+blkSize,Ny)
     911        useMaskj = []
     912        for jBlk in range(nXBlks):
     913            jBeg = jBlk*blkSize
     914            jFin = min(jBeg+blkSize,Nx)
     915            mask = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp)          #2-theta & azimuth arrays & create position mask
     916            useMaskj.append(mask)
     917        useMask.append(useMaskj)
     918    return useMask
     919
     920def ImageIntegrate(image,data,masks,blkSize=128,returnN=False,useTA=None,useMask=None):
    906921    'Integrate an image; called from OnIntegrateAll and OnIntegrate in G2imgGUI'    #for q, log(q) bins need data['binType']
    907922    import histogram2d as h2d
     
    944959            jFin = min(jBeg+blkSize,Nx)
    945960            # next is most expensive step!
     961            t0 = time.time()
    946962            if useTA:
    947963                TA = useTA[iBlk][jBlk]
    948964            else:
    949                 TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin),times)           #2-theta & azimuth arrays & create position mask
    950             tam = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp,times)
     965                TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
     966            times[1] += time.time()-t0
     967            t0 = time.time()
     968            if useMask:
     969                tam = useMask[iBlk][jBlk]
     970            else:
     971                tam = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp)
    951972            Block = image[iBeg:iFin,jBeg:jFin]
     973            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(Masks,TA,tam,Block)    #and apply masks
     974            times[0] += time.time()-t0
    952975            t0 = time.time()
    953             tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(Masks,TA,tam,Block)    #and apply masks
    954             times[2] += time.time()-t0
    955976            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
    956977            tax = np.where(tax < LRazm[0],tax+360.,tax)
     
    965986            elif 'q' == data.get('binType','').lower():
    966987                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
     988            times[2] += time.time()-t0
    967989            t0 = time.time()
    968990            taz = np.array((taz*tad/tabs),dtype='float32')
  • trunk/GSASIIimgGUI.py

    r3184 r3186  
    1919import copy
    2020import glob
     21import time
    2122import re
    2223import math
     
    144145################################################################################                   
    145146blkSize = 1024   #this seems to be optimal; will break in polymask if >1024
    146 def UpdateImageControls(G2frame,data,masks,useTA=None,IntegrateOnly=False):
     147def UpdateImageControls(G2frame,data,masks,useTA=None,useMask=None,IntegrateOnly=False):
    147148    '''Shows and handles the controls on the "Image Controls"
    148149    data tree entry
     
    238239        G2frame.slideSizer.GetChildren()[4].Window.SetValue(Imin)   #tricky
    239240         
    240     def OnIntegrate(event,useTA=None):
     241    def OnIntegrate(event,useTA=None,useMask=None):
    241242        '''Integrate image in response to a menu event or from the AutoIntegrate
    242243        dialog. In the latter case, event=None.
     
    247248        wx.BeginBusyCursor()
    248249        try:
    249             G2frame.Integrate = G2img.ImageIntegrate(sumImg,data,masks,blkSize,useTA=useTA)           
     250            G2frame.Integrate = G2img.ImageIntegrate(sumImg,data,masks,blkSize,useTA=useTA,useMask=useMask)           
    250251        finally:
    251252            wx.EndBusyCursor()   
     
    270271                    pId = 0
    271272                    oldData = {'tilt':0.,'distance':0.,'rotation':0.,'center':[0.,0.],'DetDepth':0.,'azmthOff':0.}
     273                    oldMhash = 0
    272274                    for icnt,item in enumerate(items):
    273275                        GoOn = dlgp.Update(icnt)
     
    283285                                same = False
    284286                        if not same:
    285                             print('Use new image controls')
     287                            t0 = time.time()
    286288                            useTA = G2img.MakeUseTA(Data,blkSize)
     289                            print(' Use new image controls; new xy -> th,azm time %.3f'%(time.time()-t0))
    287290                        Masks = G2frame.GPXtree.GetItemPyData(
    288291                            G2gd.GetGPXtreeItemId(G2frame,G2frame.Image,'Masks'))
     292                        Mhash = hash(str(Masks))
     293                        if  Mhash != oldMhash:
     294                            t0 = time.time()
     295                            useMask = G2img.MakeUseMask(Data,Masks,blkSize)
     296                            print(' Use new mask; make mask time: %.3f'%(time.time()-t0))
     297                            oldMhash = Mhash
    289298                        image = GetImageZ(G2frame,Data)
    290                         G2frame.Integrate = G2img.ImageIntegrate(image,Data,Masks,blkSize,useTA=useTA)
     299                        G2frame.Integrate = G2img.ImageIntegrate(image,Data,Masks,blkSize,useTA=useTA,useMask=useMask)
    291300                        del image   #force cleanup
    292301                        pId = G2IO.SaveIntegration(G2frame,CId,Data)
     
    11931202   
    11941203    if IntegrateOnly:
    1195         OnIntegrate(None,useTA=useTA)
     1204        OnIntegrate(None,useTA=useTA,useMask=useMask)
    11961205        return
    11971206   
     
    27712780        self.Pause = True
    27722781           
    2773     def IntegrateImage(self,img,useTA=None):
     2782    def IntegrateImage(self,img,useTA=None,useMask=None):
    27742783        '''Integrates a single image. Ids for created PWDR entries (more than one is possible)
    27752784        are placed in G2frame.IntgOutList
     
    27972806        # simulate a Image Controls press, since that is where the
    27982807        # integration is hidden
    2799         UpdateImageControls(G2frame,data,masks,useTA=useTA,IntegrateOnly=True)
     2808        UpdateImageControls(G2frame,data,masks,useTA=useTA,useMask=useMask,IntegrateOnly=True)
    28002809        G2frame.IntegratedList.append(img) # note this as integrated
    28012810        # split name and control number
     
    29742983        This is called only after the "Start" button is pressed (then its label reads "Pause").
    29752984        '''
    2976         def AutoIntegrateImage(imgId,useTA=None):
     2985        def AutoIntegrateImage(imgId,useTA=None,useMask=None):
    29772986            '''Integrates an image that has been read into the data tree and updates the
    29782987            AutoInt window.
     
    29923001            self.EnableButtons(False)
    29933002            try:
    2994                 self.IntegrateImage(img,useTA=useTA)
     3003                self.IntegrateImage(img,useTA=useTA,useMask=useMask)
    29953004            finally:
    29963005                self.EnableButtons(True)
     
    30793088        # have not yet been processed           
    30803089        oldData = {'tilt':0.,'distance':0.,'rotation':0.,'center':[0.,0.],'DetDepth':0.,'azmthOff':0.}
    3081         self.useTA = None
     3090        oldMhash = 0
     3091        if 'useTA' not in dir(self):    #initial definition; reuse if after Resume
     3092            self.useTA = None
     3093            self.useMask = None
    30823094        for img in G2gd.GetGPXtreeDataNames(G2frame,['IMG ']):
    30833095            imgId = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,img)
     
    30893101            Data = G2frame.GPXtree.GetItemPyData(
    30903102                G2gd.GetGPXtreeItemId(G2frame,imgId, 'Image Controls'))
    3091             same = True
     3103            sameTA = True
    30923104            for item in ['tilt','distance','rotation','center','DetDepth','azmthOff']:
    30933105                if Data[item] != oldData[item]:
    3094                     same = False
    3095             if not same:
    3096                 print('Use new image controls')
     3106                    sameTA = False
     3107            if not sameTA:
     3108                t0 = time.time()
    30973109                self.useTA = G2img.MakeUseTA(Data,blkSize)
    3098             AutoIntegrateImage(imgId,self.useTA)
     3110                print(' Use new image controls; xy->th,azm mtime: %.3f'%(time.time()-t0))
     3111            Mask = G2frame.GPXtree.GetItemPyData(
     3112                G2gd.GetGPXtreeItemId(G2frame,imgId, 'Masks'))
     3113            Mhash = hash(str(Mask))
     3114            if  Mhash != oldMhash:
     3115                t0 = time.time()
     3116                self.useMask = G2img.MakeUseMask(Data,Mask,blkSize)
     3117                print(' Use new mask; make mask time: %.3f'%(time.time()-t0))
     3118                oldMhash = Mhash
     3119            AutoIntegrateImage(imgId,self.useTA,self.useMask)
    30993120            oldData = Data
    31003121            if self.pdfControls: AutoComputePDF(imgId)
     
    31073128
    31083129        # loop over image files matching glob, reading in any new ones
     3130        if self.useTA is None or self.useMask is None:
     3131            print('Integration will not be fast; there is no beginning image controls')     #TODO: work this out??
    31093132        for newImage in self.currImageList:
    31103133            if newImage in imageFileList or self.Pause: continue # already read?
    31113134            for imgId in G2IO.ReadImages(G2frame,newImage):
    3112                 AutoIntegrateImage(imgId,self.useTA)
     3135                AutoIntegrateImage(imgId,self.useTA,self.useMask)
    31133136                if self.pdfControls: AutoComputePDF(imgId)
    31143137                self.Pause |= G2frame.PauseIntegration
  • trunk/GSASIIspc.py

    r3176 r3186  
    28492849                return spc + rhomb
    28502850    # how about the post-2002 orthorhombic names?
    2851     for i,spc in sgequiv_2002_orthorhombic:
    2852         if rspc == i.replace(' ','').upper():
    2853             return spc
     2851    if rspc in sgequiv_2002_orthorhombic:
     2852        return sgequiv_2002_orthorhombic[rspc]
     2853    else:
    28542854    # not found
    2855     return ''
     2855        return ''
    28562856
    28572857spgbyNum = []
     
    29452945    'C2/m':('C 2','C m','C c','C n',
    29462946        'C 2/m','C 2/c','C 2/n',),
    2947     'Pmmm':('P 2 2 2',
     2947    'A2/m':('A 2','A m','A a','A n',
     2948        'A 2/m','A 2/a','A 2/n',),
     2949   'Pmmm':('P 2 2 2',
    29482950        'P 2 2 21','P 21 2 2','P 2 21 2',
    29492951        'P 21 21 2','P 2 21 21','P 21 2 21',
     
    30193021        'F -4 3 c','F m -3 m','F m 3 m','F m -3 c','F d -3 m','F d -3 c',),
    30203022}
    3021 
     3023sgequiv_2002_orthorhombic = {}
     3024''' A dictionary of orthorhombic space groups that were renamed in the 2002 Volume A,
     3025 along with the pre-2002 name. The e designates a double glide-plane
     3026'''
     3027sgequiv_2002_orthorhombic = {
     3028        'AEM2':'A b m 2','B2EM':'B 2 c m','CM2E':'C m 2 a',
     3029        'AE2M':'A c 2 m','BME2':'B m a 2','C2ME':'C 2 m b',
     3030        'AEA2':'A b a 2','B2EB':'B 2 c b','CC2E':'C c 2 a',
     3031        'AE2A':'A c 2 a','BBE2':'B b a 2','C2CE':'C 2 c b',
     3032        'CMCE':'C m c a','AEMA':'A b m a','BBEM':'B b c m',
     3033        'BMEB':'B m a b','CCME':'C c m b','AEAM':'A c a m',
     3034        'CMME':'C m m a','AEMM':'A b m m','BMEM':'B m c m',
     3035        'CCCE':'C c c a','AEAA':'A b a a','BBEB':'B b c b'}
    30223036ptssdict = {}
    30233037'''A dictionary of superspace group symbols allowed for each point group
     
    31663180                      'R 3 c r','R -3 c r','R -3 m r',),
    31673181
    3168 #A list of orthorhombic space groups that were renamed in the 2002 Volume A,
    3169 # along with the pre-2002 name. The e designates a double glide-plane'''
    3170 sgequiv_2002_orthorhombic= (('A e m 2', 'A b m 2',),
    3171                             ('A e a 2', 'A b a 2',),
    3172                             ('C m c e', 'C m c a',),
    3173                             ('C m m e', 'C m m a',),
    3174                             ('C c c e', 'C c c a'),)
    31753182#Use the space groups types in this order to list the symbols in the
    31763183#order they are listed in the International Tables, vol. A'''
  • trunk/imports/G2phase_CIF.py

    r3172 r3186  
    152152            SpGrp = SpGrp.replace('_','')
    153153            # try normalizing the space group, to see if we can pick the space group out of a table
    154             SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
    155             if SpGrpNorm:
    156                 E,SGData = G2spc.SpcGroup(SpGrpNorm)
     154            E,SGData = G2spc.SpcGroup(SpGrp)
     155            if E and SpGrp:
     156                SpGrpNorm = G2spc.StandardizeSpcName(SpGrp)
     157                if SpGrpNorm:
     158                    E,SGData = G2spc.SpcGroup(SpGrpNorm)
    157159            # nope, try the space group "out of the Box"
    158             if E and SpGrp:
    159                 E,SGData = G2spc.SpcGroup(SpGrp)
    160160            if E:
    161161                if not SpGrp:
     
    196196                waveloop = blk.GetLoop('_cell_wave_vector_seq_id')
    197197                waveDict = dict(waveloop.items())
    198                 SuperVec = [[float(waveDict['_cell_wave_vector_x'][0].replace('?','0')),
    199                     float(waveDict['_cell_wave_vector_y'][0].replace('?','0')),
    200                     float(waveDict['_cell_wave_vector_z'][0].replace('?','0'))],False,4]
     198                SuperVec = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0],
     199                    cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0],
     200                    cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[0]],False,4]
    201201            # read in atoms
    202202            self.errors = 'Error during reading of atoms'
     
    263263            ranIdlookup = {}
    264264            for aitem in atomloop:
     265                mc = 0
    265266                if magnetic:
    266                     atomlist = ['','','',0.,0.,0.,1.0,0.,0.,0.,'',0.,'I',0.01,0.,0.,0.,0.,0.,0.]
    267                 else:
    268                     atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01,0.,0.,0.,0.,0.,0.]
     267                    atomlist = ['','','',0.,0.,0.,1.0, 0.,0.,0.,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
     268                    mc = 3
     269                else:
     270                    atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,]
    269271                for val,key in zip(aitem,atomkeys):
    270272                    col = G2AtomDict.get(key,-1)
    271273                    if col >= 3:
    272274                        atomlist[col] = cif.get_number_with_esd(val)[0]
    273                         if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag
    274                     elif col is not None:
     275                        if col >= 11: atomlist[9+mc] = 'A' # if any Aniso term is defined, set flag
     276                    elif col is not None and col != -1:
    275277                        atomlist[col] = val
    276278                    elif key in ('_atom_site_thermal_displace_type',
    277279                               '_atom_site_adp_type'):   #Iso or Aniso?
    278280                        if val.lower() == 'uani':
    279                             if magnetic:
    280                                 atomlist[12] = 'A'
    281                             else:
    282                                 atomlist[9] = 'A'
     281                            atomlist[9+mc] = 'A'
    283282                    elif key == '_atom_site_u_iso_or_equiv':
    284283                        uisoval = cif.get_number_with_esd(val)[0]
    285284                        if uisoval is not None:
    286                             if magnetic:
    287                                 atomlist[13] = uisoval                               
    288                             else:   
    289                                 atomlist[10] = uisoval
     285                            atomlist[10+mc] = uisoval
    290286                if not atomlist[1] and atomlist[0]:
    291287                    typ = atomlist[0].rstrip('0123456789-+')
     
    296292                        self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n'
    297293                if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop?
    298                     if magnetic:
    299                         atomlist[12] = 'A'
    300                     else:
    301                         atomlist[9] = 'A'
    302                     for val,key in zip( # load the values
    303                             anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),
    304                             anisokeys):
     294                    atomlist[9+mc] = 'A'
     295                    for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys):
    305296                        col = G2AtomDict.get(key)
    306                         if magnetic:
    307                             col += 3
    308297                        if col:
    309                             atomlist[col] = cif.get_number_with_esd(val)[0]
     298                            atomlist[col+mc] = cif.get_number_with_esd(val)[0]
    310299                if magnetic:
    311300                    for mitem in magatomloop:
     
    317306                                    atomlist[mcol] = cif.get_number_with_esd(mval)[0]
    318307                            break                           
    319                     atomlist[10],atomlist[11] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
    320                 else:
    321                     atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
     308                atomlist[7+mc],atomlist[8+mc] = G2spc.SytSym(atomlist[3:6],SGData)[:2]
    322309                atomlist[1] = G2elem.FixValence(atomlist[1])
    323310                atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id
Note: See TracChangeset for help on using the changeset viewer.