Changeset 4104

Timestamp:

Aug 20, 2019 9:06:16 PM (4 years ago)

Author:

toby

Message:

new mode for combined fit of wavelength to a set of images

Location:

trunk

Files:

• GSASIIdataGUI.py (modified) (1 diff)
• GSASIIimage.py (modified) (5 diffs)
• GSASIIimgGUI.py (modified) (3 diffs)

trunk/GSASIIdataGUI.py

@@ -5509 +5509 @@
         self.ImageEdit.Append(G2G.wxID_IMRECALIBRATE,'Recalibrate','Recalibrate detector by fitting to calibrant lines')
         self.ImageEdit.Append(G2G.wxID_IMRECALIBALL,'Recalibrate all','Recalibrate all images by fitting to calibrant lines')
+        G2G.Define_wxId('wxID_IMDISTRECALIB')
+        self.ImageEdit.Append(G2G.wxID_IMDISTRECALIB,'Multi-distance Recalibrate','Recalibrate all images varying delta-distance and fitting wavelength')
         self.ImageEdit.Append(G2G.wxID_IMCLEARCALIB,'Clear calibration','Clear calibration data points and rings')
         

trunk/GSASIIimage.py

@@ -138 +138 @@
 
 def FitDetector(rings,varyList,parmDict,Print=True,covar=False):
-    '''Fit detector calibration paremeters
+    '''Fit detector calibration parameters
 
     :param np.array rings: vector of ring positions
@@ -225 +225 @@
         return [chisq,vals,sigList]
 
+def FitMultiDist(rings,varyList,parmDict,Print=True,covar=False):
+    '''Fit detector calibration parameters with multi-distance data
+
+    :param np.array rings: vector of ring positions (x,y,dist,d-space)
+    :param list varyList: calibration parameters to be refined
+    :param dict parmDict: calibration parameters
+    :param bool Print: set to True (default) to print the results
+    :param bool covar: set to True to return the covariance matrix (default is False)
+    :returns: [chisq,vals,sigDict] unless covar is True, then 
+        [chisq,vals,sigDict,coVarMatrix] is returned
+    '''
+        
+    def CalibPrint(parmDict,sigDict,chisq,Npts):
+        print ('Image Parameters: chi**2: %12.3g, Np: %d'%(chisq,Npts))
+        ptlbls = 'names :'
+        ptstr =  'values:'
+        sigstr = 'esds  :'
+        names = ['wavelength', 'dep', 'phi', 'tilt']
+        if 'deltaDist' in parmDict:
+            names += ['deltaDist']
+        for name in names:
+            if name == 'wavelength':
+                fmt = '%12.6f'
+            elif name == 'dep':
+                fmt = '%12.4f'
+            else:
+                fmt = '%12.3f'
+                
+            ptlbls += "%s" % (name.rjust(12))
+            if name == 'phi':
+                ptstr += fmt % (parmDict[name]%360.)
+            else:
+                ptstr += fmt % (parmDict[name])
+            if name in sigDict:
+                sigstr += fmt % (sigDict[name])
+            else:
+                sigstr += 12*' '
+        print (ptlbls)
+        print (ptstr)
+        print (sigstr)
+        ptlbls = 'names :'
+        ptstr =  'values:'
+        sigstr = 'esds  :'
+        for d in sorted(set([i[5:] for i in parmDict.keys() if 'det-X' in i]),key=lambda x:int(x)):
+            fmt = '%12.3f'
+            for key in 'det-X','det-Y','delta':
+                name = key+d
+                if name not in parmDict: continue
+                ptlbls += "%12s" % name
+                ptstr += fmt % (parmDict[name])
+                if name in sigDict:
+                    sigstr += fmt % (sigDict[name])
+                else:
+                    sigstr += 12*' '
+                if len(ptlbls) > 68:
+                    print()
+                    print (ptlbls)
+                    print (ptstr)
+                    print (sigstr)
+                    ptlbls = 'names :'
+                    ptstr =  'values:'
+                    sigstr = 'esds  :'
+        if len(ptlbls) > 0:
+            print()
+            print (ptlbls)
+            print (ptstr)
+            print (sigstr)
+
+    def ellipseCalcD(B,xyd,varyList,parmDict):
+        x,y,dist,dsp = xyd
+        varyDict = dict(zip(varyList,B))
+        parms = {}
+        for parm in parmDict:
+            if parm in varyList:
+                parms[parm] = varyDict[parm]
+            else:
+                parms[parm] = parmDict[parm]
+        # create arrays with detector center values
+        detX = np.array([parms['det-X'+str(int(d))] for d in dist])
+        detY = np.array([parms['det-Y'+str(int(d))] for d in dist])
+        if 'deltaDist' in parms:
+            deltaDist = parms['deltaDist']
+        else:
+            deltaDist = np.array([parms['delta'+str(int(d))] for d in dist])
+                
+        phi = parms['phi']-90.               #get rotation of major axis from tilt axis
+        tth = 2.0*npasind(parms['wavelength']/(2.*dsp))
+        phi0 = npatan2d(y-detY,x-detX)
+        dxy = peneCorr(tth,parms['dep'],dist-deltaDist,parms['tilt'],phi0)
+        stth = npsind(tth)
+        cosb = npcosd(parms['tilt'])
+        tanb = nptand(parms['tilt'])        
+        tbm = nptand((tth-parms['tilt'])/2.)
+        tbp = nptand((tth+parms['tilt'])/2.)
+        d = (dist-deltaDist)+dxy
+        fplus = d*tanb*stth/(cosb+stth)
+        fminus = d*tanb*stth/(cosb-stth)
+        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
+        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
+        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
+        R1 = (vplus+vminus)/2.                                    #major axis
+        zdis = (fplus-fminus)/2.
+        Robs = np.sqrt((x-detX)**2+(y-detY)**2)
+        rsqplus = R0**2+R1**2
+        rsqminus = R0**2-R1**2
+        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
+        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
+        P = 2.*R0**2*zdis*npcosd(phi0-phi)
+        Rcalc = (P+Q)/R
+        return (Robs-Rcalc)*10.        #why 10? does make "chi**2" more reasonable
+        
+    p0 = [parmDict[key] for key in varyList]
+    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
+    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0))   #reduced chi^2 = M/(Nobs-Nvar)
+    parmDict.update(zip(varyList,result[0]))
+    vals = list(result[0])
+    if chisq > 1:
+        sig = list(np.sqrt(chisq*np.diag(result[1])))
+    else:
+        sig = list(np.sqrt(np.diag(result[1])))
+    sigDict = {name:s for name,s in zip(varyList,sig)}
+    if Print:
+        CalibPrint(parmDict,sigDict,chisq,rings.shape[0])
+    if covar:
+        return [chisq,vals,sigDict,result[1]]
+    else:
+        return [chisq,vals,sigDict]
+    
 def ImageLocalMax(image,w,Xpix,Ypix):
     'Needs a doc string'
@@ -556 +684 @@
     return tam.T
     
-def ImageRecalibrate(G2frame,ImageZ,data,masks):
+def ImageRecalibrate(G2frame,ImageZ,data,masks,getRingsOnly=False):
     '''Called to repeat the calibration on an image, usually called after
     calibration is done initially to improve the fit.
@@ -564 +692 @@
     :param dict data: the Controls dict for the image
     :param dict masks: a dict with masks
-    :returns: a list containing vals,varyList,sigList,parmDict,covar
+    :returns: a list containing vals,varyList,sigList,parmDict,covar or rings
+      (with an array of x, y, and d-space values) if getRingsOnly is True
+      or an empty list, in case of an error
     '''
     import ImageCalibrants as calFile
-    G2fil.G2Print ('Image recalibration:')
+    if not getRingsOnly:
+        G2fil.G2Print ('Image recalibration:')
     time0 = time.time()
     pixelSize = data['pixelSize']
@@ -640 +771 @@
         
     rings = np.concatenate((data['rings']),axis=0)
+    if getRingsOnly:
+        return rings
     [chisq,vals,sigList,covar] = FitDetector(rings,varyList,parmDict,True,True)
     data['wavelength'] = parmDict['wave']

trunk/GSASIIimgGUI.py

@@ -193 +193 @@
                 
     def OnRecalibrate(event):
+        '''Use existing calibration values as starting point for a calibration
+        fit
+        '''
         G2img.ImageRecalibrate(G2frame,G2frame.ImageZ,data,masks)
         wx.CallAfter(UpdateImageControls,G2frame,data,masks)
         
     def OnRecalibAll(event):
+        '''Use existing calibration values as starting point for a calibration
+        fit for a selected series of images
+        '''
         Names = G2gd.GetGPXtreeDataNames(G2frame,['IMG ',])
         dlg = G2G.G2MultiChoiceDialog(G2frame,'Image calibration controls','Select images to recalibrate:',Names)
@@ -247 +253 @@
         G2plt.PlotExposedImage(G2frame,event=None)
         G2frame.GPXtree.SelectItem(Id)
+
+    def OnDistRecalib(event):
+        '''Assemble rings & calibration input for a series of images with
+        differing distances
+        '''
+        obsArr = np.array([]).reshape(0,4)
+        parmDict = {}
+        varList = []
+        Names = G2gd.GetGPXtreeDataNames(G2frame,['IMG ',])
+        dlg = G2G.G2MultiChoiceDialog(G2frame,'Image calibration controls','Select images to recalibrate:',Names)
+        try:
+            wx.BeginBusyCursor()
+            if dlg.ShowModal() == wx.ID_OK:
+                items = dlg.GetSelections()
+                print('Scanning for ring picks...')
+#                G2frame.EnablePlot = False
+                for item in items:
+                    name = Names[item]
+                    print ('getting rings for',name)
+                    G2frame.Image = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,name)
+                    Data = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,G2frame.Image,'Image Controls'))
+                    G2frame.ImageZ = GetImageZ(G2frame,Data)
+                    Data['setRings'] = True
+                    Mid = G2gd.GetGPXtreeItemId(G2frame,G2frame.Image,'Masks')
+                    Masks = G2frame.GPXtree.GetItemPyData(Mid)
+                    result = G2img.ImageRecalibrate(G2frame,G2frame.ImageZ,Data,Masks,getRingsOnly=True)
+                    if not len(result):
+                        print('calibrant missing from local image calibrants files')
+                        return
+                    # add detector set dist into data array, create a single really large array
+                    distarr = np.zeros_like(result[:,2:3])
+                    if 'setdist' not in Data:
+                        print('Distance (setdist) not in image metadata')
+                        return
+                    distarr += Data['setdist']
+                    obsArr = np.concatenate((
+                        obsArr,
+                        np.concatenate((result[:,0:2],distarr,result[:,2:3]),axis=1)),axis=0)
+                    # create a parameter dict for combined fit
+                    if 'wavelength' not in parmDict:
+                        parmDict['wavelength'] = Data['wavelength']
+                        if Data['varyList']['wave']:
+                            varList += ['wavelength']
+                        parmDict['dep'] = Data['DetDepth']
+                        if Data['varyList']['dep']:
+                            varList += ['dep']
+                        # distance flag determines if individual values are refined
+                        if not Data['varyList']['dist']:
+                            # starts as zero, single variable, always refined
+                            parmDict['deltaDist'] = 0.
+                            varList += ['deltaDist']
+                        parmDict['phi'] = Data['rotation']
+                        if Data['varyList']['phi']:
+                            varList += ['phi']
+                        parmDict['tilt'] = Data['tilt']
+                        if Data['varyList']['tilt']:
+                            varList += ['tilt']
+                    key = str(int(Data['setdist']))
+                    if 'deltaDist' not in parmDict:
+                        # starts as zero, variable refined for each image
+                         parmDict['delta'+key] = 0
+                         varList += ['delta'+key]
+                    for i,z in enumerate(['X','Y']):
+                        v = 'det-'+z
+                        if v+key in parmDict:
+                            print('Error: two images with setdist ~=',key)
+                            return
+                        parmDict[v+key] = Data['center'][i]
+                        if Data['varyList'][v]:
+                            varList += [v+key]
+                #GSASIIpath.IPyBreak()
+                print('\nFitting',len(obsArr.shape[0]),'ring picks...')
+                result = G2img.FitMultiDist(obsArr,varList,parmDict)
+                # create a sequential table?
+#                Id =  G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Sequential image calibration results')
+#                if Id:
+#                    SeqResult = G2frame.GPXtree.GetItemPyData(Id)
+#                else:
+#                    Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text='Sequential image calibration results')
+                #SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
+#                    vals,varyList,sigList,parmDict,covar = result
+#                    sigList = list(sigList)
+#                    if 'dist' not in varyList:
+#                        vals.append(parmDict['dist'])
+#                        varyList.append('dist')
+#                        sigList.append(None)
+#                    vals.append(Data.get('setdist',Data['distance']))
+#                    # add setdist to varylist etc. so that it is displayed in Seq Res table
+#                    varyList.append('setdist')
+#                    sigList.append(None)
+#                    covar = np.lib.pad(covar, (0,1), 'constant')
+#                    vals.append(Data.get('samplechangerpos',Data['samplechangerpos']))
+#                    varyList.append('chgrpos')
+#                    sigList.append(None)
+                    
+#                    SeqResult[name] = {'variables':vals,'varyList':varyList,'sig':sigList,'Rvals':[],
+#                        'covMatrix':covar,'title':name,'parmDict':parmDict}
+#                SeqResult['histNames'] = Names                
+#                G2frame.GPXtree.SetItemPyData(Id,SeqResult)
+        finally:
+            dlg.Destroy()
+            wx.EndBusyCursor()
+
+#        print ('All selected images recalibrated - results in Sequential image calibration results')
+#        G2frame.G2plotNB.Delete('Sequential refinement')    #clear away probably invalid plot
+#        G2plt.PlotExposedImage(G2frame,event=None)
+#        G2frame.GPXtree.SelectItem(Id)
+
         
     def OnClearCalib(event):
@@ -1287 +1401 @@
     G2frame.Bind(wx.EVT_MENU, OnRecalibrate, id=G2G.wxID_IMRECALIBRATE)
     G2frame.Bind(wx.EVT_MENU, OnRecalibAll, id=G2G.wxID_IMRECALIBALL)
+    G2frame.Bind(wx.EVT_MENU, OnDistRecalib, id=G2G.wxID_IMDISTRECALIB)
     G2frame.Bind(wx.EVT_MENU, OnClearCalib, id=G2G.wxID_IMCLEARCALIB)
 #    if data.get('calibrant'):