Changeset 4104


Ignore:
Timestamp:
Aug 20, 2019 9:06:16 PM (2 years ago)
Author:
toby
Message:

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

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/GSASIIdataGUI.py

    r4103 r4104  
    55095509        self.ImageEdit.Append(G2G.wxID_IMRECALIBRATE,'Recalibrate','Recalibrate detector by fitting to calibrant lines')
    55105510        self.ImageEdit.Append(G2G.wxID_IMRECALIBALL,'Recalibrate all','Recalibrate all images by fitting to calibrant lines')
     5511        G2G.Define_wxId('wxID_IMDISTRECALIB')
     5512        self.ImageEdit.Append(G2G.wxID_IMDISTRECALIB,'Multi-distance Recalibrate','Recalibrate all images varying delta-distance and fitting wavelength')
    55115513        self.ImageEdit.Append(G2G.wxID_IMCLEARCALIB,'Clear calibration','Clear calibration data points and rings')
    55125514       
  • trunk/GSASIIimage.py

    r4102 r4104  
    138138
    139139def FitDetector(rings,varyList,parmDict,Print=True,covar=False):
    140     '''Fit detector calibration paremeters
     140    '''Fit detector calibration parameters
    141141
    142142    :param np.array rings: vector of ring positions
     
    225225        return [chisq,vals,sigList]
    226226
     227def FitMultiDist(rings,varyList,parmDict,Print=True,covar=False):
     228    '''Fit detector calibration parameters with multi-distance data
     229
     230    :param np.array rings: vector of ring positions (x,y,dist,d-space)
     231    :param list varyList: calibration parameters to be refined
     232    :param dict parmDict: calibration parameters
     233    :param bool Print: set to True (default) to print the results
     234    :param bool covar: set to True to return the covariance matrix (default is False)
     235    :returns: [chisq,vals,sigDict] unless covar is True, then
     236        [chisq,vals,sigDict,coVarMatrix] is returned
     237    '''
     238       
     239    def CalibPrint(parmDict,sigDict,chisq,Npts):
     240        print ('Image Parameters: chi**2: %12.3g, Np: %d'%(chisq,Npts))
     241        ptlbls = 'names :'
     242        ptstr =  'values:'
     243        sigstr = 'esds  :'
     244        names = ['wavelength', 'dep', 'phi', 'tilt']
     245        if 'deltaDist' in parmDict:
     246            names += ['deltaDist']
     247        for name in names:
     248            if name == 'wavelength':
     249                fmt = '%12.6f'
     250            elif name == 'dep':
     251                fmt = '%12.4f'
     252            else:
     253                fmt = '%12.3f'
     254               
     255            ptlbls += "%s" % (name.rjust(12))
     256            if name == 'phi':
     257                ptstr += fmt % (parmDict[name]%360.)
     258            else:
     259                ptstr += fmt % (parmDict[name])
     260            if name in sigDict:
     261                sigstr += fmt % (sigDict[name])
     262            else:
     263                sigstr += 12*' '
     264        print (ptlbls)
     265        print (ptstr)
     266        print (sigstr)
     267        ptlbls = 'names :'
     268        ptstr =  'values:'
     269        sigstr = 'esds  :'
     270        for d in sorted(set([i[5:] for i in parmDict.keys() if 'det-X' in i]),key=lambda x:int(x)):
     271            fmt = '%12.3f'
     272            for key in 'det-X','det-Y','delta':
     273                name = key+d
     274                if name not in parmDict: continue
     275                ptlbls += "%12s" % name
     276                ptstr += fmt % (parmDict[name])
     277                if name in sigDict:
     278                    sigstr += fmt % (sigDict[name])
     279                else:
     280                    sigstr += 12*' '
     281                if len(ptlbls) > 68:
     282                    print()
     283                    print (ptlbls)
     284                    print (ptstr)
     285                    print (sigstr)
     286                    ptlbls = 'names :'
     287                    ptstr =  'values:'
     288                    sigstr = 'esds  :'
     289        if len(ptlbls) > 0:
     290            print()
     291            print (ptlbls)
     292            print (ptstr)
     293            print (sigstr)
     294
     295    def ellipseCalcD(B,xyd,varyList,parmDict):
     296        x,y,dist,dsp = xyd
     297        varyDict = dict(zip(varyList,B))
     298        parms = {}
     299        for parm in parmDict:
     300            if parm in varyList:
     301                parms[parm] = varyDict[parm]
     302            else:
     303                parms[parm] = parmDict[parm]
     304        # create arrays with detector center values
     305        detX = np.array([parms['det-X'+str(int(d))] for d in dist])
     306        detY = np.array([parms['det-Y'+str(int(d))] for d in dist])
     307        if 'deltaDist' in parms:
     308            deltaDist = parms['deltaDist']
     309        else:
     310            deltaDist = np.array([parms['delta'+str(int(d))] for d in dist])
     311               
     312        phi = parms['phi']-90.               #get rotation of major axis from tilt axis
     313        tth = 2.0*npasind(parms['wavelength']/(2.*dsp))
     314        phi0 = npatan2d(y-detY,x-detX)
     315        dxy = peneCorr(tth,parms['dep'],dist-deltaDist,parms['tilt'],phi0)
     316        stth = npsind(tth)
     317        cosb = npcosd(parms['tilt'])
     318        tanb = nptand(parms['tilt'])       
     319        tbm = nptand((tth-parms['tilt'])/2.)
     320        tbp = nptand((tth+parms['tilt'])/2.)
     321        d = (dist-deltaDist)+dxy
     322        fplus = d*tanb*stth/(cosb+stth)
     323        fminus = d*tanb*stth/(cosb-stth)
     324        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
     325        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
     326        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
     327        R1 = (vplus+vminus)/2.                                    #major axis
     328        zdis = (fplus-fminus)/2.
     329        Robs = np.sqrt((x-detX)**2+(y-detY)**2)
     330        rsqplus = R0**2+R1**2
     331        rsqminus = R0**2-R1**2
     332        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
     333        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
     334        P = 2.*R0**2*zdis*npcosd(phi0-phi)
     335        Rcalc = (P+Q)/R
     336        return (Robs-Rcalc)*10.        #why 10? does make "chi**2" more reasonable
     337       
     338    p0 = [parmDict[key] for key in varyList]
     339    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
     340    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0))   #reduced chi^2 = M/(Nobs-Nvar)
     341    parmDict.update(zip(varyList,result[0]))
     342    vals = list(result[0])
     343    if chisq > 1:
     344        sig = list(np.sqrt(chisq*np.diag(result[1])))
     345    else:
     346        sig = list(np.sqrt(np.diag(result[1])))
     347    sigDict = {name:s for name,s in zip(varyList,sig)}
     348    if Print:
     349        CalibPrint(parmDict,sigDict,chisq,rings.shape[0])
     350    if covar:
     351        return [chisq,vals,sigDict,result[1]]
     352    else:
     353        return [chisq,vals,sigDict]
     354   
    227355def ImageLocalMax(image,w,Xpix,Ypix):
    228356    'Needs a doc string'
     
    556684    return tam.T
    557685   
    558 def ImageRecalibrate(G2frame,ImageZ,data,masks):
     686def ImageRecalibrate(G2frame,ImageZ,data,masks,getRingsOnly=False):
    559687    '''Called to repeat the calibration on an image, usually called after
    560688    calibration is done initially to improve the fit.
     
    564692    :param dict data: the Controls dict for the image
    565693    :param dict masks: a dict with masks
    566     :returns: a list containing vals,varyList,sigList,parmDict,covar
     694    :returns: a list containing vals,varyList,sigList,parmDict,covar or rings
     695      (with an array of x, y, and d-space values) if getRingsOnly is True
     696      or an empty list, in case of an error
    567697    '''
    568698    import ImageCalibrants as calFile
    569     G2fil.G2Print ('Image recalibration:')
     699    if not getRingsOnly:
     700        G2fil.G2Print ('Image recalibration:')
    570701    time0 = time.time()
    571702    pixelSize = data['pixelSize']
     
    640771       
    641772    rings = np.concatenate((data['rings']),axis=0)
     773    if getRingsOnly:
     774        return rings
    642775    [chisq,vals,sigList,covar] = FitDetector(rings,varyList,parmDict,True,True)
    643776    data['wavelength'] = parmDict['wave']
  • trunk/GSASIIimgGUI.py

    r4102 r4104  
    193193               
    194194    def OnRecalibrate(event):
     195        '''Use existing calibration values as starting point for a calibration
     196        fit
     197        '''
    195198        G2img.ImageRecalibrate(G2frame,G2frame.ImageZ,data,masks)
    196199        wx.CallAfter(UpdateImageControls,G2frame,data,masks)
    197200       
    198201    def OnRecalibAll(event):
     202        '''Use existing calibration values as starting point for a calibration
     203        fit for a selected series of images
     204        '''
    199205        Names = G2gd.GetGPXtreeDataNames(G2frame,['IMG ',])
    200206        dlg = G2G.G2MultiChoiceDialog(G2frame,'Image calibration controls','Select images to recalibrate:',Names)
     
    247253        G2plt.PlotExposedImage(G2frame,event=None)
    248254        G2frame.GPXtree.SelectItem(Id)
     255
     256    def OnDistRecalib(event):
     257        '''Assemble rings & calibration input for a series of images with
     258        differing distances
     259        '''
     260        obsArr = np.array([]).reshape(0,4)
     261        parmDict = {}
     262        varList = []
     263        Names = G2gd.GetGPXtreeDataNames(G2frame,['IMG ',])
     264        dlg = G2G.G2MultiChoiceDialog(G2frame,'Image calibration controls','Select images to recalibrate:',Names)
     265        try:
     266            wx.BeginBusyCursor()
     267            if dlg.ShowModal() == wx.ID_OK:
     268                items = dlg.GetSelections()
     269                print('Scanning for ring picks...')
     270#                G2frame.EnablePlot = False
     271                for item in items:
     272                    name = Names[item]
     273                    print ('getting rings for',name)
     274                    G2frame.Image = G2gd.GetGPXtreeItemId(G2frame,G2frame.root,name)
     275                    Data = G2frame.GPXtree.GetItemPyData(G2gd.GetGPXtreeItemId(G2frame,G2frame.Image,'Image Controls'))
     276                    G2frame.ImageZ = GetImageZ(G2frame,Data)
     277                    Data['setRings'] = True
     278                    Mid = G2gd.GetGPXtreeItemId(G2frame,G2frame.Image,'Masks')
     279                    Masks = G2frame.GPXtree.GetItemPyData(Mid)
     280                    result = G2img.ImageRecalibrate(G2frame,G2frame.ImageZ,Data,Masks,getRingsOnly=True)
     281                    if not len(result):
     282                        print('calibrant missing from local image calibrants files')
     283                        return
     284                    # add detector set dist into data array, create a single really large array
     285                    distarr = np.zeros_like(result[:,2:3])
     286                    if 'setdist' not in Data:
     287                        print('Distance (setdist) not in image metadata')
     288                        return
     289                    distarr += Data['setdist']
     290                    obsArr = np.concatenate((
     291                        obsArr,
     292                        np.concatenate((result[:,0:2],distarr,result[:,2:3]),axis=1)),axis=0)
     293                    # create a parameter dict for combined fit
     294                    if 'wavelength' not in parmDict:
     295                        parmDict['wavelength'] = Data['wavelength']
     296                        if Data['varyList']['wave']:
     297                            varList += ['wavelength']
     298                        parmDict['dep'] = Data['DetDepth']
     299                        if Data['varyList']['dep']:
     300                            varList += ['dep']
     301                        # distance flag determines if individual values are refined
     302                        if not Data['varyList']['dist']:
     303                            # starts as zero, single variable, always refined
     304                            parmDict['deltaDist'] = 0.
     305                            varList += ['deltaDist']
     306                        parmDict['phi'] = Data['rotation']
     307                        if Data['varyList']['phi']:
     308                            varList += ['phi']
     309                        parmDict['tilt'] = Data['tilt']
     310                        if Data['varyList']['tilt']:
     311                            varList += ['tilt']
     312                    key = str(int(Data['setdist']))
     313                    if 'deltaDist' not in parmDict:
     314                        # starts as zero, variable refined for each image
     315                         parmDict['delta'+key] = 0
     316                         varList += ['delta'+key]
     317                    for i,z in enumerate(['X','Y']):
     318                        v = 'det-'+z
     319                        if v+key in parmDict:
     320                            print('Error: two images with setdist ~=',key)
     321                            return
     322                        parmDict[v+key] = Data['center'][i]
     323                        if Data['varyList'][v]:
     324                            varList += [v+key]
     325                #GSASIIpath.IPyBreak()
     326                print('\nFitting',len(obsArr.shape[0]),'ring picks...')
     327                result = G2img.FitMultiDist(obsArr,varList,parmDict)
     328                # create a sequential table?
     329#                Id =  G2gd.GetGPXtreeItemId(G2frame,G2frame.root,'Sequential image calibration results')
     330#                if Id:
     331#                    SeqResult = G2frame.GPXtree.GetItemPyData(Id)
     332#                else:
     333#                    Id = G2frame.GPXtree.AppendItem(parent=G2frame.root,text='Sequential image calibration results')
     334                #SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
     335#                    vals,varyList,sigList,parmDict,covar = result
     336#                    sigList = list(sigList)
     337#                    if 'dist' not in varyList:
     338#                        vals.append(parmDict['dist'])
     339#                        varyList.append('dist')
     340#                        sigList.append(None)
     341#                    vals.append(Data.get('setdist',Data['distance']))
     342#                    # add setdist to varylist etc. so that it is displayed in Seq Res table
     343#                    varyList.append('setdist')
     344#                    sigList.append(None)
     345#                    covar = np.lib.pad(covar, (0,1), 'constant')
     346#                    vals.append(Data.get('samplechangerpos',Data['samplechangerpos']))
     347#                    varyList.append('chgrpos')
     348#                    sigList.append(None)
     349                   
     350#                    SeqResult[name] = {'variables':vals,'varyList':varyList,'sig':sigList,'Rvals':[],
     351#                        'covMatrix':covar,'title':name,'parmDict':parmDict}
     352#                SeqResult['histNames'] = Names               
     353#                G2frame.GPXtree.SetItemPyData(Id,SeqResult)
     354        finally:
     355            dlg.Destroy()
     356            wx.EndBusyCursor()
     357
     358#        print ('All selected images recalibrated - results in Sequential image calibration results')
     359#        G2frame.G2plotNB.Delete('Sequential refinement')    #clear away probably invalid plot
     360#        G2plt.PlotExposedImage(G2frame,event=None)
     361#        G2frame.GPXtree.SelectItem(Id)
     362
    249363       
    250364    def OnClearCalib(event):
     
    12871401    G2frame.Bind(wx.EVT_MENU, OnRecalibrate, id=G2G.wxID_IMRECALIBRATE)
    12881402    G2frame.Bind(wx.EVT_MENU, OnRecalibAll, id=G2G.wxID_IMRECALIBALL)
     1403    G2frame.Bind(wx.EVT_MENU, OnDistRecalib, id=G2G.wxID_IMDISTRECALIB)
    12891404    G2frame.Bind(wx.EVT_MENU, OnClearCalib, id=G2G.wxID_IMCLEARCALIB)
    12901405#    if data.get('calibrant'):
Note: See TracChangeset for help on using the changeset viewer.