Changeset 4685


Ignore:
Timestamp:
Dec 29, 2020 12:37:31 PM (9 months ago)
Author:
vondreele
Message:

Add gain map calculator

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIdataGUI.py

    r4682 r4685  
    38713871                        self.ErrorDialog('Image sum error','No nonzero image multipliers found')
    38723872                        return
    3873                        
    3874                        
     3873                                               
    38753874                    newImage = np.array(newImage,dtype=np.int32)                       
    38763875                    outname = 'IMG '+result[-1]
  • trunk/GSASIIimage.py

    r4671 r4685  
    12621262    return useMask
    12631263
     1264def MakeGainMap(image,Ix,Iy,data,masks,blkSize=128):
     1265    import scipy.ndimage.filters as sdif
     1266    Iy *= npcosd(Ix[:-1])**3       #undo parallax 3 or 4?
     1267    Iy *= (1000./data['distance'])**2    #undo r^2 effect
     1268    Iy /= np.array(G2pwd.Polarization(data['PolaVal'][0],Ix[:-1],0.)[0])    #undo polarization
     1269    if data['Oblique'][1]:
     1270        Iy *= G2pwd.Oblique(data['Oblique'][0],Ix[:-1])     #undo penetration
     1271    IyInt = scint.interp1d(Ix[:-1],Iy[0],bounds_error=False)
     1272    GainMap = np.zeros_like(image,dtype=float)
     1273    #do interpolation on all points - fills in the empty bins; leaves others the same
     1274    Nx,Ny = data['size']
     1275    nXBlks = (Nx-1)//blkSize+1
     1276    nYBlks = (Ny-1)//blkSize+1
     1277    for iBlk in range(nYBlks):
     1278        iBeg = iBlk*blkSize
     1279        iFin = min(iBeg+blkSize,Ny)
     1280        for jBlk in range(nXBlks):
     1281            jBeg = jBlk*blkSize
     1282            jFin = min(jBeg+blkSize,Nx)
     1283            TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
     1284            Ipix = IyInt(TA[0])
     1285            GainMap[iBeg:iFin,jBeg:jFin] = image[iBeg:iFin,jBeg:jFin]/(Ipix*TA[3])
     1286    GainMap = np.nan_to_num(GainMap)
     1287    GainMap = sdif.gaussian_filter(GainMap,3.,order=0)
     1288    return 1./GainMap
     1289
    12641290def ImageIntegrate(image,data,masks,blkSize=128,returnN=False,useTA=None,useMask=None):
    12651291    'Integrate an image; called from OnIntegrateAll and OnIntegrate in G2imgGUI'    #for q, log(q) bins need data['binType']
     
    13731399        H0 *= np.array(G2pwd.Polarization(data['PolaVal'][0],H2[:-1],0.)[0])
    13741400    H0 /= np.abs(npcosd(H2[:-1]-np.abs(data['det2theta'])))           #parallax correction
     1401    H0 *= (data['distance']/1000.)**2    #remove r^2 effect
    13751402    if 'SASD' in data['type']:
    13761403        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
  • trunk/GSASIIimgGUI.py

    r4672 r4685  
    122122            GMid = G2gd.GetGPXtreeItemId(G2frame, G2frame.root, gainMap)
    123123            if GMid:
    124                 Npix,gainfile,imagetag = G2IO.GetCheckImageFile(G2frame,GMid)
    125                 GMimage = G2IO.GetImageData(G2frame,gainfile,True,ImageTag=imagetag,FormatName=formatName)
     124                Npix,gainfile,imagetag = G2IO.GetCheckImageFile(G2frame,GMid)               
     125                Gdata = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,GMid,'Image Controls'))
     126                gformat = Gdata['formatName']
     127                GMimage = G2IO.GetImageData(G2frame,gainfile,True,ImageTag=imagetag,FormatName=gformat)
    126128                sumImg = sumImg*GMimage/1000
    127129    sumImg -= int(data.get('Flat Bkg',0))
     
    154156        dlg.Destroy()
    155157       
     158    def OnMakeGainMap(event):
     159        sumImg = GetImageZ(G2frame,data)
     160        masks = copy.deepcopy(G2frame.GPXtree.GetItemPyData(
     161            G2gd.GetGPXtreeItemId(G2frame,G2frame.Image,'Masks')))
     162        Data = copy.deepcopy(data)
     163        #force defaults for GainMap calc
     164        Data['IOtth'] = [0.1,60.0]
     165        Data['outAzimuths'] = 1
     166        Data['LRazimuth'] = [0.,360.]
     167        Data['outChannels'] = 2500
     168        Data['SampleAbs'] = [0.0,False]
     169        Data['binType'] = '2-theta'
     170        G2frame.Integrate = G2img.ImageIntegrate(sumImg,Data,masks,blkSize)           
     171        Iy,azms,Ix = G2frame.Integrate[:3]
     172        GainMap = G2img.MakeGainMap(sumImg,Ix,Iy,Data,masks,blkSize)*1000.
     173        Npix,imagefile,imagetag = G2IO.GetCheckImageFile(G2frame,G2frame.Image)
     174        pth = os.path.split(os.path.abspath(imagefile))[0]
     175        outname = 'GainMap'
     176        dlg = wx.FileDialog(G2frame, 'Choose gain map filename', pth,outname,
     177            'G2img files (*.G2img)|*.G2img',wx.FD_SAVE|wx.FD_OVERWRITE_PROMPT)
     178        if dlg.ShowModal() == wx.ID_OK:
     179            newimagefile = dlg.GetPath()
     180            newimagefile = G2IO.FileDlgFixExt(dlg,newimagefile)
     181            Data['formatName'] = 'GSAS-II image'
     182            Data['range'] = [(800,1200),[800,1200]]
     183            G2IO.PutG2Image(newimagefile,[],data,Npix,GainMap)
     184            Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text='IMG '+os.path.split(newimagefile)[1])
     185            G2frame.GPXtree.SetItemPyData(Id,[Npix,newimagefile])
     186            G2frame.GPXtree.SetItemPyData(G2frame.GPXtree.AppendItem(Id,text='Comments'),[])
     187            G2frame.GPXtree.SetItemPyData(G2frame.GPXtree.AppendItem(Id,text='Image Controls'),Data)
     188            G2frame.GPXtree.SetItemPyData(G2frame.GPXtree.AppendItem(Id,text='Masks'),masks)
     189            G2frame.GPXtree.Expand(Id)
     190            G2frame.GPXtree.SelectItem(Id)      #to show the gain map & put it in the list
     191
    156192    G2frame.dataWindow.ClearData()
    157193    G2frame.ImageZ = GetImageZ(G2frame,data)
     
    194230    tthSizer.Add(wx.StaticText(G2frame.dataWindow,label=' Sample changer position %.2f mm '%data['samplechangerpos']),0,WACV)
    195231    mainSizer.Add(tthSizer,0)
     232    if not data['Gain map']:
     233        makeGain = wx.Button(G2frame.dataWindow,label='Make gain map from this image? NB: use only an amorphous pattern for this')
     234        makeGain.Bind(wx.EVT_BUTTON,OnMakeGainMap)
     235        mainSizer.Add(makeGain)
    196236    G2frame.dataWindow.SetDataSize()
    197237
     
    473513        '''Integrate image in response to a menu event or from the AutoIntegrate
    474514        dialog. In the latter case, event=None.
    475         TODO: think of not repeating x,y-->2th,azm calc for unchanging calibrations?
    476515        '''
    477516        CleanupMasks(masks)
     
    13111350            ResetThresholds()
    13121351            wx.CallAfter(G2plt.PlotExposedImage,G2frame,event=event)
    1313        
     1352           
    13141353        backSizer = wx.FlexGridSizer(0,6,5,5)
    13151354        oldFlat = data.get('Flat Bkg',0.)
     
    13441383        backSizer.Add((5,5),0)
    13451384        backSizer.Add(wx.StaticText(G2frame.dataWindow,-1,' Gain map'),0,WACV)
    1346         gainMap = wx.ComboBox(parent=G2frame.dataWindow,value=data['Gain map'],choices=Choices,
     1385        gainMap = wx.ComboBox(G2frame.dataWindow,value=data['Gain map'],choices=Choices,
    13471386            style=wx.CB_READONLY|wx.CB_DROPDOWN)
    13481387        gainMap.Bind(wx.EVT_COMBOBOX,OnGainMap)
Note: See TracChangeset for help on using the changeset viewer.