Changeset 3186


Ignore:
Timestamp:
Dec 9, 2017 8:54:59 AM (5 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.