source: trunk/GSASIIimage.py @ 4671

Last change on this file since 4671 was 4671, checked in by vondreele, 11 months ago

add refinable multiplier for a fixed background; multiplier is now normally > 0
usable for peak fitting and Rietved refinement
add the fixed background entry to all background defaults
fix bug in GetDetectorXY fo when cursor outside image - returns [0,0] not None; changes elsewhere to use this
GetTthAzmDsp? now returns explicit list not assumed tuple
put the abs in the nl.qr test for singularities in HessianLSQ
Add 'BF mult' to name list in G2obj
put a try - except TypeError? around setting plot style stuff in PlotPatterns?
Remove the setting of Pattern[0]BackFile? - this was redundant for PWDR
remove picker/pickradius from linescan plot - failed
remove the alternate fixed background definition ('_fixedVary', etc.)
clear the PhaseReaderClass?Drawing? dictionary upon import of phase from a gpx file

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 67.5 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2020-12-12 19:30:31 +0000 (Sat, 12 Dec 2020) $
5# $Author: vondreele $
6# $Revision: 4671 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 4671 2020-12-12 19:30:31Z vondreele $
9########### SVN repository information ###################
10'''
11*GSASIIimage: Image calc module*
12================================
13
14Ellipse fitting & image integration
15
16'''
17from __future__ import division, print_function
18import math
19import time
20import numpy as np
21import numpy.linalg as nl
22import numpy.ma as ma
23from scipy.optimize import leastsq
24import scipy.interpolate as scint
25import scipy.special as sc
26import copy
27import GSASIIpath
28GSASIIpath.SetVersionNumber("$Revision: 4671 $")
29try:
30    import GSASIIplot as G2plt
31except ImportError: # expected in scriptable w/o matplotlib and/or wx
32    pass
33import GSASIIlattice as G2lat
34import GSASIIpwd as G2pwd
35import GSASIIspc as G2spc
36#import GSASIImath as G2mth
37import GSASIIfiles as G2fil
38
39# trig functions in degrees
40sind = lambda x: math.sin(x*math.pi/180.)
41asind = lambda x: 180.*math.asin(x)/math.pi
42tand = lambda x: math.tan(x*math.pi/180.)
43atand = lambda x: 180.*math.atan(x)/math.pi
44atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
45cosd = lambda x: math.cos(x*math.pi/180.)
46acosd = lambda x: 180.*math.acos(x)/math.pi
47rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
48#numpy versions
49npsind = lambda x: np.sin(x*np.pi/180.)
50npasind = lambda x: 180.*np.arcsin(x)/np.pi
51npcosd = lambda x: np.cos(x*np.pi/180.)
52npacosd = lambda x: 180.*np.arccos(x)/np.pi
53nptand = lambda x: np.tan(x*np.pi/180.)
54npatand = lambda x: 180.*np.arctan(x)/np.pi
55npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
56nxs = np.newaxis
57debug = False
58   
59def pointInPolygon(pXY,xy):
60    'Needs a doc string'
61    #pXY - assumed closed 1st & last points are duplicates
62    Inside = False
63    N = len(pXY)
64    p1x,p1y = pXY[0]
65    for i in range(N+1):
66        p2x,p2y = pXY[i%N]
67        if (max(p1y,p2y) >= xy[1] > min(p1y,p2y)) and (xy[0] <= max(p1x,p2x)):
68            if p1y != p2y:
69                xinters = (xy[1]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
70            if p1x == p2x or xy[0] <= xinters:
71                Inside = not Inside
72        p1x,p1y = p2x,p2y
73    return Inside
74   
75def peneCorr(tth,dep,dist):
76    'Needs a doc string'
77    return dep*(1.-npcosd(tth))*dist**2/1000.         #best one
78       
79def makeMat(Angle,Axis):
80    '''Make rotation matrix from Angle and Axis
81
82    :param float Angle: in degrees
83    :param int Axis: 0 for rotation about x, 1 for about y, etc.
84    '''
85    cs = npcosd(Angle)
86    ss = npsind(Angle)
87    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
88    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
89                   
90def FitEllipse(xy):
91   
92    def ellipse_center(p):
93        ''' gives ellipse center coordinates
94        '''
95        b,c,d,f,a = p[1]/2., p[2], p[3]/2., p[4]/2., p[0]
96        num = b*b-a*c
97        x0=(c*d-b*f)/num
98        y0=(a*f-b*d)/num
99        return np.array([x0,y0])
100   
101    def ellipse_angle_of_rotation( p ):
102        ''' gives rotation of ellipse major axis from x-axis
103        range will be -90 to 90 deg
104        '''
105        b,c,a = p[1]/2., p[2], p[0]
106        return 0.5*npatand(2*b/(a-c))
107   
108    def ellipse_axis_length( p ):
109        ''' gives ellipse radii in [minor,major] order
110        '''
111        b,c,d,f,g,a = p[1]/2., p[2], p[3]/2., p[4]/2, p[5], p[0]
112        up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
113        down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
114        down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
115        res1=np.sqrt(up/down1)
116        res2=np.sqrt(up/down2)
117        return np.array([ res2,res1])
118   
119    xy = np.array(xy)
120    x = np.asarray(xy.T[0])[:,np.newaxis]
121    y = np.asarray(xy.T[1])[:,np.newaxis]
122    D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
123    S = np.dot(D.T,D)
124    C = np.zeros([6,6])
125    C[0,2] = C[2,0] = 2; C[1,1] = -1
126    E, V =  nl.eig(np.dot(nl.inv(S), C))
127    n = np.argmax(np.abs(E))
128    a = V[:,n]
129    cent = ellipse_center(a)
130    phi = ellipse_angle_of_rotation(a)
131    radii = ellipse_axis_length(a)
132    phi += 90.
133    if radii[0] > radii[1]:
134        radii = [radii[1],radii[0]]
135        phi -= 90.
136    return cent,phi,radii
137
138def FitDetector(rings,varyList,parmDict,Print=True,covar=False):
139    '''Fit detector calibration parameters
140
141    :param np.array rings: vector of ring positions
142    :param list varyList: calibration parameters to be refined
143    :param dict parmDict: all calibration parameters
144    :param bool Print: set to True (default) to print the results
145    :param bool covar: set to True to return the covariance matrix (default is False)
146    :returns: [chisq,vals,sigList] unless covar is True, then
147        [chisq,vals,sigList,coVarMatrix] is returned
148    '''
149       
150    def CalibPrint(ValSig,chisq,Npts):
151        print ('Image Parameters: chi**2: %12.3g, Np: %d'%(chisq,Npts))
152        ptlbls = 'names :'
153        ptstr =  'values:'
154        sigstr = 'esds  :'
155        for name,value,sig in ValSig:
156            ptlbls += "%s" % (name.rjust(12))
157            if name == 'phi':
158                ptstr += Fmt[name] % (value%360.)
159            else:
160                ptstr += Fmt[name] % (value)
161            if sig:
162                sigstr += Fmt[name] % (sig)
163            else:
164                sigstr += 12*' '
165        print (ptlbls)
166        print (ptstr)
167        print (sigstr)       
168       
169    def ellipseCalcD(B,xyd,varyList,parmDict):
170       
171        x,y,dsp = xyd
172        varyDict = dict(zip(varyList,B))
173        parms = {}
174        for parm in parmDict:
175            if parm in varyList:
176                parms[parm] = varyDict[parm]
177            else:
178                parms[parm] = parmDict[parm]
179        phi = parms['phi']-90.               #get rotation of major axis from tilt axis
180        tth = 2.0*npasind(parms['wave']/(2.*dsp))
181        phi0 = npatan2d(y-parms['det-Y'],x-parms['det-X'])
182        dxy = peneCorr(tth,parms['dep'],parms['dist'])
183        stth = npsind(tth)
184        cosb = npcosd(parms['tilt'])
185        tanb = nptand(parms['tilt'])       
186        tbm = nptand((tth-parms['tilt'])/2.)
187        tbp = nptand((tth+parms['tilt'])/2.)
188        d = parms['dist']+dxy
189        fplus = d*tanb*stth/(cosb+stth)
190        fminus = d*tanb*stth/(cosb-stth)
191        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
192        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
193        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
194        R1 = (vplus+vminus)/2.                                    #major axis
195        zdis = (fplus-fminus)/2.
196        Robs = np.sqrt((x-parms['det-X'])**2+(y-parms['det-Y'])**2)
197        rsqplus = R0**2+R1**2
198        rsqminus = R0**2-R1**2
199        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
200        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
201        P = 2.*R0**2*zdis*npcosd(phi0-phi)
202        Rcalc = (P+Q)/R
203        M = (Robs-Rcalc)*25.        #why 25? does make "chi**2" more reasonable
204        return M
205       
206    names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
207    fmt = ['%12.3f','%12.3f','%12.3f','%12.3f','%12.3f','%12.4f','%12.6f']
208    Fmt = dict(zip(names,fmt))
209    p0 = [parmDict[key] for key in varyList]
210    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
211    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0))   #reduced chi^2 = M/(Nobs-Nvar)
212    parmDict.update(zip(varyList,result[0]))
213    vals = list(result[0])
214    if not len(vals):
215        sig = []
216        ValSig = []
217        sigList = []
218    else:   
219        sig = list(np.sqrt(chisq*np.diag(result[1])))
220        sigList = np.zeros(7)
221        for i,name in enumerate(varyList):
222            sigList[i] = sig[varyList.index(name)]
223        ValSig = zip(varyList,vals,sig)
224    if Print:
225        if len(sig):
226            CalibPrint(ValSig,chisq,rings.shape[0])
227        else:
228            print(' Nothing refined')
229    if covar:
230        return [chisq,vals,sigList,result[1]]
231    else:
232        return [chisq,vals,sigList]
233
234def FitMultiDist(rings,varyList,parmDict,Print=True,covar=False):
235    '''Fit detector calibration parameters with multi-distance data
236
237    :param np.array rings: vector of ring positions (x,y,dist,d-space)
238    :param list varyList: calibration parameters to be refined
239    :param dict parmDict: calibration parameters
240    :param bool Print: set to True (default) to print the results
241    :param bool covar: set to True to return the covariance matrix (default is False)
242    :returns: [chisq,vals,sigDict] unless covar is True, then
243        [chisq,vals,sigDict,coVarMatrix] is returned
244    '''
245       
246    def CalibPrint(parmDict,sigDict,chisq,Npts):
247        ptlbls = 'names :'
248        ptstr =  'values:'
249        sigstr = 'esds  :'
250        for d in sorted(set([i[5:] for i in parmDict.keys() if 'det-X' in i]),key=lambda x:int(x)):
251            fmt = '%12.3f'
252            for key in 'det-X','det-Y','delta':
253                name = key+d
254                if name not in parmDict: continue
255                ptlbls += "%12s" % name
256                ptstr += fmt % (parmDict[name])
257                if name in sigDict:
258                    sigstr += fmt % (sigDict[name])
259                else:
260                    sigstr += 12*' '
261                if len(ptlbls) > 68:
262                    print()
263                    print (ptlbls)
264                    print (ptstr)
265                    print (sigstr)
266                    ptlbls = 'names :'
267                    ptstr =  'values:'
268                    sigstr = 'esds  :'
269        if len(ptlbls) > 8:
270            print()
271            print (ptlbls)
272            print (ptstr)
273            print (sigstr)
274        print ('\nImage Parameters: chi**2: %12.3g, Np: %d'%(chisq,Npts))
275        ptlbls = 'names :'
276        ptstr =  'values:'
277        sigstr = 'esds  :'
278        names = ['wavelength', 'dep', 'phi', 'tilt']
279        if 'deltaDist' in parmDict:
280            names += ['deltaDist']
281        for name in names:
282            if name == 'wavelength':
283                fmt = '%12.6f'
284            elif name == 'dep':
285                fmt = '%12.4f'
286            else:
287                fmt = '%12.3f'
288               
289            ptlbls += "%s" % (name.rjust(12))
290            if name == 'phi':
291                ptstr += fmt % (parmDict[name]%360.)
292            else:
293                ptstr += fmt % (parmDict[name])
294            if name in sigDict:
295                sigstr += fmt % (sigDict[name])
296            else:
297                sigstr += 12*' '
298        print (ptlbls)
299        print (ptstr)
300        print (sigstr)
301        print()
302
303    def ellipseCalcD(B,xyd,varyList,parmDict):
304        x,y,dist,dsp = xyd
305        varyDict = dict(zip(varyList,B))
306        parms = {}
307        for parm in parmDict:
308            if parm in varyList:
309                parms[parm] = varyDict[parm]
310            else:
311                parms[parm] = parmDict[parm]
312        # create arrays with detector center values
313        detX = np.array([parms['det-X'+str(int(d))] for d in dist])
314        detY = np.array([parms['det-Y'+str(int(d))] for d in dist])
315        if 'deltaDist' in parms:
316            deltaDist = parms['deltaDist']
317        else:
318            deltaDist = np.array([parms['delta'+str(int(d))] for d in dist])
319               
320        phi = parms['phi']-90.               #get rotation of major axis from tilt axis
321        tth = 2.0*npasind(parms['wavelength']/(2.*dsp))
322        phi0 = npatan2d(y-detY,x-detX)
323        dxy = peneCorr(tth,parms['dep'],dist-deltaDist)
324        stth = npsind(tth)
325        cosb = npcosd(parms['tilt'])
326        tanb = nptand(parms['tilt'])       
327        tbm = nptand((tth-parms['tilt'])/2.)
328        tbp = nptand((tth+parms['tilt'])/2.)
329        d = (dist-deltaDist)+dxy
330        fplus = d*tanb*stth/(cosb+stth)
331        fminus = d*tanb*stth/(cosb-stth)
332        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
333        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
334        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
335        R1 = (vplus+vminus)/2.                                    #major axis
336        zdis = (fplus-fminus)/2.
337        Robs = np.sqrt((x-detX)**2+(y-detY)**2)
338        rsqplus = R0**2+R1**2
339        rsqminus = R0**2-R1**2
340        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
341        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
342        P = 2.*R0**2*zdis*npcosd(phi0-phi)
343        Rcalc = (P+Q)/R
344        return (Robs-Rcalc)*25.        #why 25? does make "chi**2" more reasonable
345       
346    p0 = [parmDict[key] for key in varyList]
347    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
348    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0))   #reduced chi^2 = M/(Nobs-Nvar)
349    parmDict.update(zip(varyList,result[0]))
350    vals = list(result[0])
351    if chisq > 1:
352        sig = list(np.sqrt(chisq*np.diag(result[1])))
353    else:
354        sig = list(np.sqrt(np.diag(result[1])))
355    sigDict = {name:s for name,s in zip(varyList,sig)}
356    if Print:
357        CalibPrint(parmDict,sigDict,chisq,rings.shape[0])
358    if covar:
359        return [chisq,vals,sigDict,result[1]]
360    else:
361        return [chisq,vals,sigDict]
362   
363def ImageLocalMax(image,w,Xpix,Ypix):
364    'Needs a doc string'
365    w2 = w*2
366    sizey,sizex = image.shape
367    xpix = int(Xpix)            #get reference corner of pixel chosen
368    ypix = int(Ypix)
369    if not w:
370        ZMax = np.sum(image[ypix-2:ypix+2,xpix-2:xpix+2])
371        return xpix,ypix,ZMax,0.0001
372    if (w2 < xpix < sizex-w2) and (w2 < ypix < sizey-w2) and image[ypix,xpix]:
373        ZMax = image[ypix-w:ypix+w,xpix-w:xpix+w]
374        Zmax = np.argmax(ZMax)
375        ZMin = image[ypix-w2:ypix+w2,xpix-w2:xpix+w2]
376        Zmin = np.argmin(ZMin)
377        xpix += Zmax%w2-w
378        ypix += Zmax//w2-w
379        return xpix,ypix,np.ravel(ZMax)[Zmax],max(0.0001,np.ravel(ZMin)[Zmin])   #avoid neg/zero minimum
380    else:
381        return 0,0,0,0     
382   
383def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image,mul=1):
384    'Needs a doc string'
385    def ellipseC():
386        'compute estimate of ellipse circumference'
387        if radii[0] < 0:        #hyperbola
388#            theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
389#            print (theta)
390            return 0
391        apb = radii[1]+radii[0]
392        amb = radii[1]-radii[0]
393        return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
394       
395    cent,phi,radii = ellipse
396    cphi = cosd(phi-90.)        #convert to major axis rotation
397    sphi = sind(phi-90.)
398    ring = []
399    C = int(ellipseC())*mul         #ring circumference in mm
400    azm = []
401    for i in range(0,C,1):      #step around ring in 1mm increments
402        a = 360.*i/C
403        x = radii[1]*cosd(a-phi+90.)        #major axis
404        y = radii[0]*sind(a-phi+90.)
405        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
406        Y = (sphi*x+cphi*y+cent[1])*scaley
407        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
408        if I and J and float(I)/J > reject:
409            X += .5                             #set to center of pixel
410            Y += .5
411            X /= scalex                         #convert back to mm
412            Y /= scaley
413            if [X,Y,dsp] not in ring:           #no duplicates!
414                ring.append([X,Y,dsp])
415                azm.append(a)
416    if len(ring) < 10:
417        ring = []
418        azm = []
419    return ring,azm
420   
421def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
422    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
423    on output
424    radii[0] (b-minor axis) set < 0. for hyperbola
425   
426    '''
427    radii = [0,0]
428    stth = sind(tth)
429    cosb = cosd(tilt)
430    tanb = tand(tilt)
431    tbm = tand((tth-tilt)/2.)
432    tbp = tand((tth+tilt)/2.)
433    sinb = sind(tilt)
434    d = dist+dxy
435    if tth+abs(tilt) < 90.:      #ellipse
436        fplus = d*tanb*stth/(cosb+stth)
437        fminus = d*tanb*stth/(cosb-stth)
438        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
439        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
440        radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
441        radii[1] = (vplus+vminus)/2.                                    #major axis
442        zdis = (fplus-fminus)/2.
443    else:   #hyperbola!
444        f = d*abs(tanb)*stth/(cosb+stth)
445        v = d*(abs(tanb)+tand(tth-abs(tilt)))
446        delt = d*stth*(1.+stth*cosb)/(abs(sinb)*cosb*(stth+cosb))
447        eps = (v-f)/(delt-v)
448        radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.)                     #-minor axis
449        radii[1] = eps*(delt-f)/(eps**2-1.)                             #major axis
450        if tilt > 0:
451            zdis = f+radii[1]*eps
452        else:
453            zdis = -f
454#NB: zdis is || to major axis & phi is rotation of minor axis
455#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
456    elcent = [cent[0]+zdis*sind(phi),cent[1]-zdis*cosd(phi)]
457    return elcent,phi,radii
458   
459def GetEllipse(dsp,data):
460    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
461    as given in image controls dictionary (data) and a d-spacing (dsp)
462    '''
463    cent = data['center']
464    tilt = data['tilt']
465    phi = data['rotation']
466    dep = data.get('DetDepth',0.0)
467    tth = 2.0*asind(data['wavelength']/(2.*dsp))
468    dist = data['distance']
469    dxy = peneCorr(tth,dep,dist)
470    return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
471
472def GetDetectorXY(dsp,azm,data):
473    '''Get detector x,y position from d-spacing (dsp), azimuth (azm,deg)
474    & image controls dictionary (data) - new version
475    it seems to be only used in plotting
476    '''
477    def LinePlaneCollision(planeNormal, planePoint, rayDirection, rayPoint, epsilon=1e-6):
478     
479        ndotu = planeNormal.dot(rayDirection)
480        if ndotu < epsilon:   
481                return None
482     
483        w = rayPoint - planePoint
484        si = -planeNormal.dot(w) / ndotu
485        Psi = w + si * rayDirection + planePoint
486        return Psi
487   
488   
489    dist = data['distance']
490    cent = data['center']
491    T = makeMat(data['tilt'],0)
492    R = makeMat(data['rotation'],2)
493    MN = np.inner(R,np.inner(R,T))
494    iMN= nl.inv(MN)
495    tth = 2.0*npasind(data['wavelength']/(2.*dsp))
496    vect = np.array([npsind(tth)*npcosd(azm),npsind(tth)*npsind(azm),npcosd(tth)])
497    dxyz0 = np.inner(np.array([0.,0.,1.0]),MN)    #tilt detector normal
498    dxyz0 += np.array([0.,0.,dist])                 #translate to distance
499    dxyz0 = np.inner(dxyz0,makeMat(data['det2theta'],1).T)   #rotate on 2-theta
500    dxyz1 = np.inner(np.array([cent[0],cent[1],0.]),MN)    #tilt detector cent
501    dxyz1 += np.array([0.,0.,dist])                 #translate to distance
502    dxyz1 = np.inner(dxyz1,makeMat(data['det2theta'],1).T)   #rotate on 2-theta
503    xyz = LinePlaneCollision(dxyz0,dxyz1,vect,2.*dist*vect)
504    if xyz is None:
505        return np.zeros(2)
506#        return None
507    xyz = np.inner(xyz,makeMat(-data['det2theta'],1).T)
508    xyz -= np.array([0.,0.,dist])                 #translate back
509    xyz = np.inner(xyz,iMN)
510    return np.squeeze(xyz)[:2]+cent
511       
512def GetDetectorXY2(dsp,azm,data):
513    '''Get detector x,y position from d-spacing (dsp), azimuth (azm,deg)
514    & image controls dictionary (data)
515    it seems to be only used in plotting
516    '''   
517    elcent,phi,radii = GetEllipse(dsp,data)
518    phi = data['rotation']-90.          #to give rotation of major axis
519    tilt = data['tilt']
520    dist = data['distance']
521    cent = data['center']
522    tth = 2.0*asind(data['wavelength']/(2.*dsp))
523    stth = sind(tth)
524    cosb = cosd(tilt)
525    if radii[0] > 0.:
526        sinb = sind(tilt)
527        tanb = tand(tilt)
528        fplus = dist*tanb*stth/(cosb+stth)
529        fminus = dist*tanb*stth/(cosb-stth)
530        zdis = (fplus-fminus)/2.
531        rsqplus = radii[0]**2+radii[1]**2
532        rsqminus = radii[0]**2-radii[1]**2
533        R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus
534        Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm-phi)**2)
535        P = 2.*radii[0]**2*zdis*cosd(azm-phi)
536        radius = (P+Q)/R
537        xy = np.array([radius*cosd(azm),radius*sind(azm)])
538        xy += cent
539    else:   #hyperbola - both branches (one is way off screen!)
540        sinb = abs(sind(tilt))
541        tanb = abs(tand(tilt))
542        f = dist*tanb*stth/(cosb+stth)
543        v = dist*(tanb+tand(tth-abs(tilt)))
544        delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
545        ecc = (v-f)/(delt-v)
546        R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm))
547        if tilt > 0.:
548            offset = 2.*radii[1]*ecc+f      #select other branch
549            xy = [-R*cosd(azm)-offset,-R*sind(azm)]
550        else:
551            offset = -f
552            xy = [-R*cosd(azm)-offset,R*sind(azm)]
553        xy = -np.array([xy[0]*cosd(phi)+xy[1]*sind(phi),xy[0]*sind(phi)-xy[1]*cosd(phi)])
554        xy += cent
555    if data['det2theta']:
556        xy[0] += dist*nptand(data['det2theta']+data['tilt']*npsind(data['rotation']))
557    return xy
558   
559def GetDetXYfromThAzm(Th,Azm,data):
560    '''Computes a detector position from a 2theta angle and an azimultal
561    angle (both in degrees) - apparently not used!
562    '''
563    dsp = data['wavelength']/(2.0*npsind(Th))   
564    return GetDetectorXY(dsp,Azm,data)
565                   
566def GetTthAzmDsp2(x,y,data): #expensive
567    '''Computes a 2theta, etc. from a detector position and calibration constants - checked
568    OK for ellipses & hyperbola.
569
570    :returns: np.array(tth,azm,G,dsp) where tth is 2theta, azm is the azimutal angle,
571       G is ? and dsp is the d-space
572    '''
573    wave = data['wavelength']
574    cent = data['center']
575    tilt = data['tilt']
576    dist = data['distance']/cosd(tilt)
577    x0 = dist*tand(tilt)
578    phi = data['rotation']
579    dep = data.get('DetDepth',0.)
580    azmthoff = data['azmthOff']
581    dx = np.array(x-cent[0],dtype=np.float32)
582    dy = np.array(y-cent[1],dtype=np.float32)
583    D = ((dx-x0)**2+dy**2+dist**2)      #sample to pixel distance
584    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
585    X = np.dot(X,makeMat(phi,2))
586    Z = np.dot(X,makeMat(tilt,0)).T[2]
587    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
588    dxy = peneCorr(tth,dep,dist)
589    DX = dist-Z+dxy
590    DY = np.sqrt(dx**2+dy**2-Z**2)
591    tth = npatan2d(DY,DX) 
592    dsp = wave/(2.*npsind(tth/2.))
593    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
594    G = D/dist**2       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
595    return np.array([tth,azm,G,dsp])
596   
597def GetTthAzmDsp(x,y,data): #expensive
598    '''Computes a 2theta, etc. from a detector position and calibration constants - checked
599    OK for ellipses & hyperbola.
600
601    :returns: np.array(tth,azm,G,dsp) where tth is 2theta, azm is the azimutal angle,
602       G is ? and dsp is the d-space
603    '''
604   
605    def costth(xyz):
606        u = xyz/nl.norm(xyz,axis=-1)[:,:,nxs]
607        return np.dot(u,np.array([0.,0.,1.]))
608       
609#zero detector 2-theta: tested with tilted images - perfect integrations
610    wave = data['wavelength']
611    dx = x-data['center'][0]
612    dy = y-data['center'][1]
613    tilt = data['tilt']
614    dist = data['distance']/npcosd(tilt)    #sample-beam intersection point
615    T = makeMat(tilt,0)
616    R = makeMat(data['rotation'],2)
617    MN = np.inner(R,np.inner(R,T))
618    dxyz0 = np.inner(np.dstack([dx,dy,np.zeros_like(dx)]),MN)    #correct for 45 deg tilt
619    dxyz0 += np.array([0.,0.,dist])
620    if data['DetDepth']:
621        ctth0 = costth(dxyz0)
622        tth0 = npacosd(ctth0)
623        dzp = peneCorr(tth0,data['DetDepth'],dist)
624        dxyz0[:,:,2] += dzp
625#non zero detector 2-theta:
626    if data['det2theta']:       
627        tthMat = makeMat(data['det2theta'],1)
628        dxyz = np.inner(dxyz0,tthMat.T)
629    else:
630        dxyz = dxyz0
631    ctth = costth(dxyz)
632    tth = npacosd(ctth)
633    dsp = wave/(2.*npsind(tth/2.))
634    azm = (npatan2d(dxyz[:,:,1],dxyz[:,:,0])+data['azmthOff']+720.)%360.       
635# G-calculation       
636    x0 = data['distance']*nptand(tilt)
637    x0x = x0*npcosd(data['rotation'])
638    x0y = x0*npsind(data['rotation'])
639    distsq = data['distance']**2
640    G = ((dx-x0x)**2+(dy-x0y)**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
641    return [tth,azm,G,dsp]
642   
643def GetTth(x,y,data):
644    'Give 2-theta value for detector x,y position; calibration info in data'
645    return GetTthAzmDsp(x,y,data)[0]
646   
647def GetTthAzm(x,y,data):
648    'Give 2-theta, azimuth values for detector x,y position; calibration info in data'
649    return GetTthAzmDsp(x,y,data)[0:2]
650   
651def GetTthAzmG2(x,y,data):
652    '''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
653     calibration info in data - only used in integration - old version
654    '''
655    'Needs a doc string - checked OK for ellipses & hyperbola'
656    tilt = data['tilt']
657    dist = data['distance']/npcosd(tilt)
658    MN = -np.inner(makeMat(data['rotation'],2),makeMat(tilt,0))
659    dx = x-data['center'][0]
660    dy = y-data['center'][1]
661    dz = np.dot(np.dstack([dx.T,dy.T,np.zeros_like(dx.T)]),MN).T[2]   
662    xyZ = dx**2+dy**2-dz**2   
663    tth0 = npatand(np.sqrt(xyZ)/(dist-dz))
664    dzp = peneCorr(tth0,data['DetDepth'],dist)
665    tth = npatan2d(np.sqrt(xyZ),dist-dz+dzp) 
666    azm = (npatan2d(dy,dx)+data['azmthOff']+720.)%360.
667   
668    distsq = data['distance']**2
669    x0 = data['distance']*nptand(tilt)
670    x0x = x0*npcosd(data['rotation'])
671    x0y = x0*npsind(data['rotation'])
672    G = ((dx-x0x)**2+(dy-x0y)**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
673    return tth,azm,G
674
675def GetTthAzmG(x,y,data):
676    '''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
677     calibration info in data - only used in integration
678     checked OK for ellipses & hyperbola
679     This is the slow step in image integration
680     '''     
681    def costth(xyz):
682        u = xyz/nl.norm(xyz,axis=-1)[:,:,nxs]
683        return np.dot(u,np.array([0.,0.,1.]))
684       
685#zero detector 2-theta: tested with tilted images - perfect integrations
686    dx = x-data['center'][0]
687    dy = y-data['center'][1]
688    tilt = data['tilt']
689    dist = data['distance']/npcosd(tilt)    #sample-beam intersection point
690    T = makeMat(tilt,0)
691    R = makeMat(data['rotation'],2)
692    MN = np.inner(R,np.inner(R,T))
693    dxyz0 = np.inner(np.dstack([dx,dy,np.zeros_like(dx)]),MN)    #correct for 45 deg tilt
694    dxyz0 += np.array([0.,0.,dist])
695    if data['DetDepth']:
696        ctth0 = costth(dxyz0)
697        tth0 = npacosd(ctth0)
698        dzp = peneCorr(tth0,data['DetDepth'],dist)
699        dxyz0[:,:,2] += dzp
700#non zero detector 2-theta:
701    if data['det2theta']:       
702        tthMat = makeMat(data['det2theta'],1)
703        dxyz = np.inner(dxyz0,tthMat.T)
704    else:
705        dxyz = dxyz0
706    ctth = costth(dxyz)
707    tth = npacosd(ctth)
708    azm = (npatan2d(dxyz[:,:,1],dxyz[:,:,0])+data['azmthOff']+720.)%360.       
709# G-calculation       
710    x0 = data['distance']*nptand(tilt)
711    x0x = x0*npcosd(data['rotation'])
712    x0y = x0*npsind(data['rotation'])
713    distsq = data['distance']**2
714    G = ((dx-x0x)**2+(dy-x0y)**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
715    return tth,azm,G
716
717def GetDsp(x,y,data):
718    'Give d-spacing value for detector x,y position; calibration info in data'
719    return GetTthAzmDsp(x,y,data)[3]
720       
721def GetAzm(x,y,data):
722    'Give azimuth value for detector x,y position; calibration info in data'
723    return GetTthAzmDsp(x,y,data)[1]
724   
725def meanAzm(a,b):
726    AZM = lambda a,b: npacosd(0.5*(npsind(2.*b)-npsind(2.*a))/(np.pi*(b-a)/180.))/2.
727    azm = AZM(a,b)
728#    quad = int((a+b)/180.)
729#    if quad == 1:
730#        azm = 180.-azm
731#    elif quad == 2:
732#        azm += 180.
733#    elif quad == 3:
734#        azm = 360-azm
735    return azm     
736       
737def ImageCompress(image,scale):
738    ''' Reduces size of image by selecting every n'th point
739    param: image array: original image
740    param: scale int: intervsl between selected points
741    returns: array: reduced size image
742    '''
743    if scale == 1:
744        return image
745    else:
746        return image[::scale,::scale]
747       
748def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
749    'Needs a doc string'
750    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
751    curr = np.array([dist,x,y])
752    return abs(avg-curr)/avg < .02
753
754def GetLineScan(image,data):
755    Nx,Ny = data['size']
756    pixelSize = data['pixelSize']
757    scalex = 1000./pixelSize[0]         #microns --> 1/mm
758    scaley = 1000./pixelSize[1]
759    wave = data['wavelength']
760    numChans = data['outChannels']
761    LUtth = np.array(data['IOtth'],dtype=np.float)
762    azm = data['linescan'][1]-data['azmthOff']
763    Tx = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
764    Ty = np.zeros_like(Tx)
765    dsp = wave/(2.0*npsind(Tx/2.0))
766    xy = [GetDetectorXY(d,azm,data) for d in dsp]
767    xy = np.array(xy).T
768    xy[1] *= scalex
769    xy[0] *= scaley
770    xy = np.array(xy,dtype=int)
771    Xpix = ma.masked_outside(xy[1],0,Ny-1)
772    Ypix = ma.masked_outside(xy[0],0,Nx-1)
773    xpix = Xpix[~(Xpix.mask+Ypix.mask)].compressed()
774    ypix = Ypix[~(Xpix.mask+Ypix.mask)].compressed()
775    Ty = image[xpix,ypix]
776    Tx = ma.array(Tx,mask=Xpix.mask+Ypix.mask).compressed()
777    return [Tx,Ty]
778
779def EdgeFinder(image,data):
780    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
781    Not currently used but might be useful in future?
782    '''
783    import numpy.ma as ma
784    Nx,Ny = data['size']
785    pixelSize = data['pixelSize']
786    edgemin = data['edgemin']
787    scalex = pixelSize[0]/1000.
788    scaley = pixelSize[1]/1000.   
789    tay,tax = np.mgrid[0:Nx,0:Ny]
790    tax = np.asfarray(tax*scalex,dtype=np.float32)
791    tay = np.asfarray(tay*scaley,dtype=np.float32)
792    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
793    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
794    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
795    return zip(tax,tay)
796   
797def MakeFrameMask(data,frame):
798    import polymask as pm
799    pixelSize = data['pixelSize']
800    scalex = pixelSize[0]/1000.
801    scaley = pixelSize[1]/1000.
802    blkSize = 512
803    Nx,Ny = data['size']
804    nXBlks = (Nx-1)//blkSize+1
805    nYBlks = (Ny-1)//blkSize+1
806    tam = ma.make_mask_none(data['size'])
807    for iBlk in range(nXBlks):
808        iBeg = iBlk*blkSize
809        iFin = min(iBeg+blkSize,Nx)
810        for jBlk in range(nYBlks):
811            jBeg = jBlk*blkSize
812            jFin = min(jBeg+blkSize,Ny)               
813            nI = iFin-iBeg
814            nJ = jFin-jBeg
815            tax,tay = np.mgrid[iBeg+0.5:iFin+.5,jBeg+.5:jFin+.5]         #bin centers not corners
816            tax = np.asfarray(tax*scalex,dtype=np.float32)
817            tay = np.asfarray(tay*scaley,dtype=np.float32)
818            tamp = ma.make_mask_none((1024*1024))
819            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
820                tay.flatten(),len(frame),frame,tamp)[:nI*nJ])^True  #switch to exclude around frame
821            if tamp.shape:
822                tamp = np.reshape(tamp[:nI*nJ],(nI,nJ))
823                tam[iBeg:iFin,jBeg:jFin] = ma.mask_or(tamp[0:nI,0:nJ],tam[iBeg:iFin,jBeg:jFin])
824            else:
825                tam[iBeg:iFin,jBeg:jFin] = True
826    return tam.T
827   
828def ImageRecalibrate(G2frame,ImageZ,data,masks,getRingsOnly=False):
829    '''Called to repeat the calibration on an image, usually called after
830    calibration is done initially to improve the fit.
831
832    :param G2frame: The top-level GSAS-II frame or None, to skip plotting
833    :param np.Array ImageZ: the image to calibrate
834    :param dict data: the Controls dict for the image
835    :param dict masks: a dict with masks
836    :returns: a list containing vals,varyList,sigList,parmDict,covar or rings
837      (with an array of x, y, and d-space values) if getRingsOnly is True
838      or an empty list, in case of an error
839    '''
840    import ImageCalibrants as calFile
841    if not getRingsOnly:
842        G2fil.G2Print ('Image recalibration:')
843    time0 = time.time()
844    pixelSize = data['pixelSize']
845    scalex = 1000./pixelSize[0]
846    scaley = 1000./pixelSize[1]
847    pixLimit = data['pixLimit']
848    cutoff = data['cutoff']
849    data['rings'] = []
850    data['ellipses'] = []
851    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
852        data['DetDepth'] /= data['distance']
853    if not data['calibrant']:
854        G2fil.G2Print ('warning: no calibration material selected')
855        return []   
856    skip = data['calibskip']
857    dmin = data['calibdmin']
858    if data['calibrant'] not in calFile.Calibrants:
859        G2fil.G2Print('Warning: %s not in local copy of image calibrants file'%data['calibrant'])
860        return []
861    calibrant = calFile.Calibrants[data['calibrant']]
862    Bravais,SGs,Cells = calibrant[:3]
863    HKL = []
864    for bravais,sg,cell in zip(Bravais,SGs,Cells):
865        A = G2lat.cell2A(cell)
866        if sg:
867            SGData = G2spc.SpcGroup(sg)[1]
868            hkl = G2pwd.getHKLpeak(dmin,SGData,A,Inst=None,nodup=True)
869            HKL += list(hkl)
870        else:
871            hkl = G2lat.GenHBravais(dmin,bravais,A)
872            HKL += list(hkl)
873    if len(calibrant) > 5:
874        absent = calibrant[5]
875    else:
876        absent = ()
877    HKL = G2lat.sortHKLd(HKL,True,False)
878    varyList = [item for item in data['varyList'] if data['varyList'][item]]
879    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
880        'setdist':data.get('setdist',data['distance']),
881        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
882    Found = False
883    wave = data['wavelength']
884    frame = masks['Frames']
885    tam = ma.make_mask_none(ImageZ.shape)
886    if frame:
887        tam = ma.mask_or(tam,MakeFrameMask(data,frame))
888    for iH,H in enumerate(HKL):
889        if debug:   print (H) 
890        dsp = H[3]
891        tth = 2.0*asind(wave/(2.*dsp))
892        if tth+abs(data['tilt']) > 90.:
893            G2fil.G2Print ('next line is a hyperbola - search stopped')
894            break
895        ellipse = GetEllipse(dsp,data)
896        if iH not in absent and iH >= skip:
897            Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
898        else:
899            Ring = makeRing(dsp,ellipse,pixLimit,1000.0,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
900        if Ring:
901            if iH not in absent and iH >= skip:
902                data['rings'].append(np.array(Ring))
903            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
904            Found = True
905        elif not Found:         #skipping inner rings, keep looking until ring found
906            continue
907        else:                   #no more rings beyond edge of detector
908            data['ellipses'].append([])
909            continue
910    if not data['rings']:
911        G2fil.G2Print ('no rings found; try lower Min ring I/Ib',mode='warn')
912        return []   
913       
914    rings = np.concatenate((data['rings']),axis=0)
915    if getRingsOnly:
916        return rings,HKL
917    [chisq,vals,sigList,covar] = FitDetector(rings,varyList,parmDict,True,True)
918    data['wavelength'] = parmDict['wave']
919    data['distance'] = parmDict['dist']
920    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
921    data['rotation'] = np.mod(parmDict['phi'],360.0)
922    data['tilt'] = parmDict['tilt']
923    data['DetDepth'] = parmDict['dep']
924    data['chisq'] = chisq
925    N = len(data['ellipses'])
926    data['ellipses'] = []           #clear away individual ellipse fits
927    for H in HKL[:N]:
928        ellipse = GetEllipse(H[3],data)
929        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))   
930    G2fil.G2Print ('calibration time = %.3f'%(time.time()-time0))
931    if G2frame:
932        G2plt.PlotImage(G2frame,newImage=True)       
933    return [vals,varyList,sigList,parmDict,covar]
934
935def ImageCalibrate(G2frame,data):
936    '''Called to perform an initial image calibration after points have been
937    selected for the inner ring.
938    '''
939    import ImageCalibrants as calFile
940    G2fil.G2Print ('Image calibration:')
941    time0 = time.time()
942    ring = data['ring']
943    pixelSize = data['pixelSize']
944    scalex = 1000./pixelSize[0]
945    scaley = 1000./pixelSize[1]
946    pixLimit = data['pixLimit']
947    cutoff = data['cutoff']
948    varyDict = data['varyList']
949    if varyDict['dist'] and varyDict['wave']:
950        G2fil.G2Print ('ERROR - you can not simultaneously calibrate distance and wavelength')
951        return False
952    if len(ring) < 5:
953        G2fil.G2Print ('ERROR - not enough inner ring points for ellipse')
954        return False
955       
956    #fit start points on inner ring
957    data['ellipses'] = []
958    data['rings'] = []
959    outE = FitEllipse(ring)
960    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
961    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi**2: %.3f, Np: %d'
962    if outE:
963        G2fil.G2Print (fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1]))
964        ellipse = outE
965    else:
966        return False
967       
968    #setup 360 points on that ring for "good" fit
969    data['ellipses'].append(ellipse[:]+('g',))
970    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
971    if Ring:
972        ellipse = FitEllipse(Ring)
973        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]    #do again
974        ellipse = FitEllipse(Ring)
975    else:
976        G2fil.G2Print ('1st ring not sufficiently complete to proceed',mode='warn')
977        return False
978    if debug:
979        G2fil.G2Print (fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],
980            ellipse[2][0],ellipse[2][1],0.,len(Ring)))     #cent,phi,radii
981    data['ellipses'].append(ellipse[:]+('r',))
982    data['rings'].append(np.array(Ring))
983    G2plt.PlotImage(G2frame,newImage=True)
984   
985#setup for calibration
986    data['rings'] = []
987    if not data['calibrant']:
988        G2fil.G2Print ('Warning: no calibration material selected')
989        return True
990   
991    skip = data['calibskip']
992    dmin = data['calibdmin']
993#generate reflection set
994    calibrant = calFile.Calibrants[data['calibrant']]
995    Bravais,SGs,Cells = calibrant[:3]
996    HKL = []
997    for bravais,sg,cell in zip(Bravais,SGs,Cells):
998        A = G2lat.cell2A(cell)
999        if sg:
1000            SGData = G2spc.SpcGroup(sg)[1]
1001            hkl = G2pwd.getHKLpeak(dmin,SGData,A,Inst=None,nodup=True)
1002            HKL += list(hkl)
1003        else:
1004            hkl = G2lat.GenHBravais(dmin,bravais,A)
1005            HKL += list(hkl)
1006    HKL = G2lat.sortHKLd(HKL,True,False)[skip:]
1007#set up 1st ring
1008    elcent,phi,radii = ellipse              #from fit of 1st ring
1009    dsp = HKL[0][3]
1010    G2fil.G2Print ('1st ring: try %.4f'%(dsp))
1011    if varyDict['dist']:
1012        wave = data['wavelength']
1013        tth = 2.0*asind(wave/(2.*dsp))
1014    else:   #varyDict['wave']!
1015        dist = data['distance']
1016        tth = npatan2d(radii[0],dist)
1017        data['wavelength'] = wave =  2.0*dsp*sind(tth/2.0)
1018    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
1019    ttth = nptand(tth)
1020    ctth = npcosd(tth)
1021#1st estimate of tilt; assume ellipse - don't know sign though
1022    if varyDict['tilt']:
1023        tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
1024        if not tilt:
1025            G2fil.G2Print ('WARNING - selected ring was fitted as a circle')
1026            G2fil.G2Print (' - if detector was tilted we suggest you skip this ring - WARNING')
1027    else:
1028        tilt = data['tilt']
1029#1st estimate of dist: sample to detector normal to plane
1030    if varyDict['dist']:
1031        data['distance'] = dist = radii[0]**2/(ttth*radii[1])
1032    else:
1033        dist = data['distance']
1034    if varyDict['tilt']:
1035#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
1036        zdisp = radii[1]*ttth*tand(tilt)
1037        zdism = radii[1]*ttth*tand(-tilt)
1038#cone axis position; 2 choices. Which is right?     
1039#NB: zdisp is || to major axis & phi is rotation of minor axis
1040#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
1041        centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
1042        centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
1043#check get same ellipse parms either way
1044#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
1045        fail = True
1046        i2 = 1
1047        while fail:
1048            dsp = HKL[i2][3]
1049            G2fil.G2Print ('2nd ring: try %.4f'%(dsp))
1050            tth = 2.0*asind(wave/(2.*dsp))
1051            ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
1052            G2fil.G2Print (fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1]))
1053            Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
1054            parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
1055                'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}       
1056            varyList = [item for item in varyDict if varyDict[item]]
1057            if len(Ringp) > 10:
1058                chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)[0]
1059                tiltp = parmDict['tilt']
1060                phip = parmDict['phi']
1061                centp = [parmDict['det-X'],parmDict['det-Y']]
1062                fail = False
1063            else:
1064                chip = 1e6
1065            ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
1066            G2fil.G2Print (fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1]))
1067            Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
1068            if len(Ringm) > 10:
1069                parmDict['tilt'] *= -1
1070                chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)[0]
1071                tiltm = parmDict['tilt']
1072                phim = parmDict['phi']
1073                centm = [parmDict['det-X'],parmDict['det-Y']]
1074                fail = False
1075            else:
1076                chim = 1e6
1077            if fail:
1078                i2 += 1
1079        if chip < chim:
1080            data['tilt'] = tiltp
1081            data['center'] = centp
1082            data['rotation'] = phip
1083        else:
1084            data['tilt'] = tiltm
1085            data['center'] = centm
1086            data['rotation'] = phim
1087        data['ellipses'].append(ellipsep[:]+('b',))
1088        data['rings'].append(np.array(Ringp))
1089        data['ellipses'].append(ellipsem[:]+('r',))
1090        data['rings'].append(np.array(Ringm))
1091        G2plt.PlotImage(G2frame,newImage=True)
1092    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
1093        data['DetDepth'] /= data['distance']
1094    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
1095        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
1096    varyList = [item for item in varyDict if varyDict[item]]
1097    data['rings'] = []
1098    data['ellipses'] = []
1099    for i,H in enumerate(HKL):
1100        dsp = H[3]
1101        tth = 2.0*asind(wave/(2.*dsp))
1102        if tth+abs(data['tilt']) > 90.:
1103            G2fil.G2Print ('next line is a hyperbola - search stopped')
1104            break
1105        if debug:   print ('HKLD:'+str(H[:4])+'2-theta: %.4f'%(tth))
1106        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
1107        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
1108        if debug:   print (fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1]))
1109        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
1110        if Ring:
1111            data['rings'].append(np.array(Ring))
1112            rings = np.concatenate((data['rings']),axis=0)
1113            if i:
1114                chisq = FitDetector(rings,varyList,parmDict,False)[0]
1115                data['distance'] = parmDict['dist']
1116                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
1117                data['rotation'] = parmDict['phi']
1118                data['tilt'] = parmDict['tilt']
1119                data['DetDepth'] = parmDict['dep']
1120                data['chisq'] = chisq
1121                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
1122                if debug:   print (fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings)))
1123            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
1124#            G2plt.PlotImage(G2frame,newImage=True)
1125        else:
1126            if debug:   print ('insufficient number of points in this ellipse to fit')
1127#            break
1128    G2plt.PlotImage(G2frame,newImage=True)
1129    fullSize = len(G2frame.ImageZ)/scalex
1130    if 2*radii[1] < .9*fullSize:
1131        G2fil.G2Print ('Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib')
1132    N = len(data['ellipses'])
1133    if N > 2:
1134        FitDetector(rings,varyList,parmDict)[0]
1135        data['wavelength'] = parmDict['wave']
1136        data['distance'] = parmDict['dist']
1137        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
1138        data['rotation'] = parmDict['phi']
1139        data['tilt'] = parmDict['tilt']
1140        data['DetDepth'] = parmDict['dep']
1141    for H in HKL[:N]:
1142        ellipse = GetEllipse(H[3],data)
1143        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
1144    G2fil.G2Print ('calibration time = %.3f'%(time.time()-time0))
1145    G2plt.PlotImage(G2frame,newImage=True)       
1146    return True
1147   
1148def Make2ThetaAzimuthMap(data,iLim,jLim): #most expensive part of integration!
1149    'Needs a doc string'
1150    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
1151    pixelSize = data['pixelSize']
1152    scalex = pixelSize[0]/1000.
1153    scaley = pixelSize[1]/1000.
1154    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
1155    tax = np.asfarray(tax*scalex,dtype=np.float32).flatten()
1156    tay = np.asfarray(tay*scaley,dtype=np.float32).flatten()
1157    nI = iLim[1]-iLim[0]
1158    nJ = jLim[1]-jLim[0]
1159    TA = np.empty((4,nI,nJ))
1160    TA[:3] = 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
1161    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
1162    TA[3] = G2pwd.Polarization(data['PolaVal'][0],TA[0],TA[1]-90.)[0]
1163    return TA           #2-theta, azimuth & geom. corr. arrays
1164
1165def MakeMaskMap(data,masks,iLim,jLim,tamp):
1166    import polymask as pm
1167    pixelSize = data['pixelSize']
1168    scalex = pixelSize[0]/1000.
1169    scaley = pixelSize[1]/1000.
1170   
1171    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
1172    tax = np.asfarray(tax*scalex,dtype=np.float32).flatten()
1173    tay = np.asfarray(tay*scaley,dtype=np.float32).flatten()
1174    nI = iLim[1]-iLim[0]
1175    nJ = jLim[1]-jLim[0]
1176    #make position masks here
1177    frame = masks['Frames']
1178    tam = ma.make_mask_none((nI*nJ))
1179    if frame:
1180        tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
1181            tay,len(frame),frame,tamp)[:nI*nJ])^True)
1182    polygons = masks['Polygons']
1183    for polygon in polygons:
1184        if polygon:
1185            tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
1186                tay,len(polygon),polygon,tamp)[:nI*nJ]))
1187    for X,Y,rsq in masks['Points'].T:
1188        tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,rsq)))
1189    if tam.shape: 
1190        tam = np.reshape(tam,(nI,nJ))
1191    else:
1192        tam = ma.make_mask_none((nI,nJ))
1193    for xline in masks.get('Xlines',[]):    #a y pixel position
1194        if iLim[0] <= xline <= iLim[1]:
1195            tam[xline-iLim[0],:] = True
1196    for yline in masks.get('Ylines',[]):    #a x pixel position
1197        if jLim[0] <= yline <= jLim[1]:
1198            tam[:,yline-jLim[0]] = True           
1199    return tam           #position mask
1200
1201def Fill2ThetaAzimuthMap(masks,TA,tam,image):
1202    'Needs a doc string'
1203    Zlim = masks['Thresholds'][1]
1204    rings = masks['Rings']
1205    arcs = masks['Arcs']
1206    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2]),ma.getdata(TA[3])))    #azimuth, 2-theta, dist, pol
1207    tax,tay,tad,pol = np.dsplit(TA,4)    #azimuth, 2-theta, dist**2/d0**2, pol
1208    for tth,thick in rings:
1209        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
1210    for tth,azm,thick in arcs:
1211        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
1212        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
1213        tam = ma.mask_or(tam.flatten(),tamt*tama)
1214    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
1215    tabs = np.ones_like(taz)
1216    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
1217    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
1218    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
1219    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
1220    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
1221    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
1222    pol = ma.compressed(ma.array(pol.flatten(),mask=tam))   #polarization
1223    return tax,tay,taz/pol,tad,tabs
1224
1225def MakeUseTA(data,blkSize=128):
1226    Nx,Ny = data['size']
1227    nXBlks = (Nx-1)//blkSize+1
1228    nYBlks = (Ny-1)//blkSize+1
1229    useTA = []
1230    for iBlk in range(nYBlks):
1231        iBeg = iBlk*blkSize
1232        iFin = min(iBeg+blkSize,Ny)
1233        useTAj = []
1234        for jBlk in range(nXBlks):
1235            jBeg = jBlk*blkSize
1236            jFin = min(jBeg+blkSize,Nx)
1237            TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))          #2-theta & azimuth arrays & create position mask
1238            useTAj.append(TA)
1239        useTA.append(useTAj)
1240    return useTA
1241
1242def MakeUseMask(data,masks,blkSize=128):
1243    Masks = copy.deepcopy(masks)
1244    Masks['Points'] = np.array(Masks['Points']).T           #get spots as X,Y,R arrays
1245    if np.any(masks['Points']):
1246        Masks['Points'][2] = np.square(Masks['Points'][2]/2.)
1247    Nx,Ny = data['size']
1248    nXBlks = (Nx-1)//blkSize+1
1249    nYBlks = (Ny-1)//blkSize+1
1250    useMask = []
1251    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
1252    for iBlk in range(nYBlks):
1253        iBeg = iBlk*blkSize
1254        iFin = min(iBeg+blkSize,Ny)
1255        useMaskj = []
1256        for jBlk in range(nXBlks):
1257            jBeg = jBlk*blkSize
1258            jFin = min(jBeg+blkSize,Nx)
1259            mask = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp)          #2-theta & azimuth arrays & create position mask
1260            useMaskj.append(mask)
1261        useMask.append(useMaskj)
1262    return useMask
1263
1264def ImageIntegrate(image,data,masks,blkSize=128,returnN=False,useTA=None,useMask=None):
1265    'Integrate an image; called from OnIntegrateAll and OnIntegrate in G2imgGUI'    #for q, log(q) bins need data['binType']
1266    import histogram2d as h2d
1267    G2fil.G2Print ('Begin image integration; image range: %d %d'%(np.min(image),np.max(image)))
1268    CancelPressed = False
1269    LUtth = np.array(data['IOtth'])
1270    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
1271    numAzms = data['outAzimuths']
1272    numChans = (data['outChannels']//4)*4
1273    Dazm = (LRazm[1]-LRazm[0])/numAzms
1274    if '2-theta' in data.get('binType','2-theta'):
1275        lutth = LUtth               
1276    elif 'log(q)' in data['binType']:
1277        lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
1278    elif 'q' == data['binType'].lower():
1279        lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
1280    dtth = (lutth[1]-lutth[0])/numChans
1281    muT = data.get('SampleAbs',[0.0,''])[0]
1282    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
1283        data['DetDepth'] /= data['distance']
1284    if 'SASD' in data['type']:
1285        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
1286    Masks = copy.deepcopy(masks)
1287    Masks['Points'] = np.array(Masks['Points']).T           #get spots as X,Y,R arrays
1288    if np.any(masks['Points']):
1289        Masks['Points'][2] = np.square(Masks['Points'][2]/2.)
1290    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
1291    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
1292    H2 = np.linspace(lutth[0],lutth[1],numChans+1)
1293    Nx,Ny = data['size']
1294    nXBlks = (Nx-1)//blkSize+1
1295    nYBlks = (Ny-1)//blkSize+1
1296    tbeg = time.time()
1297    times = [0,0,0,0,0]
1298    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
1299    for iBlk in range(nYBlks):
1300        iBeg = iBlk*blkSize
1301        iFin = min(iBeg+blkSize,Ny)
1302        for jBlk in range(nXBlks):
1303            jBeg = jBlk*blkSize
1304            jFin = min(jBeg+blkSize,Nx)
1305            # next is most expensive step!
1306
1307            t0 = time.time()
1308            if useTA:
1309                TA = useTA[iBlk][jBlk]
1310            else:
1311                TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
1312            times[1] += time.time()-t0      #xy->th,azm
1313
1314            t0 = time.time()
1315            if useMask:
1316                tam = useMask[iBlk][jBlk]
1317            else:
1318                tam = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp)
1319            times[0] += time.time()-t0      #apply masks
1320
1321            t0 = time.time()
1322            Block = image[iBeg:iFin,jBeg:jFin]          #apply image spotmask here
1323            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(Masks,TA,tam,Block)    #and apply masks
1324            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
1325            tax = np.where(tax < LRazm[0],tax+360.,tax)
1326            if data.get('SampleAbs',[0.0,''])[1]:
1327                if 'Cylind' in data['SampleShape']:
1328                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))      #adjust for additional thickness off sample normal
1329                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
1330                elif 'Fixed' in data['SampleShape']:    #assumes flat plate sample normal to beam
1331                    tabs = G2pwd.Absorb('Fixed',muT,tay)
1332            if 'log(q)' in data.get('binType',''):
1333                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
1334            elif 'q' == data.get('binType','').lower():
1335                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
1336            times[2] += time.time()-t0          #fill map
1337
1338            t0 = time.time()
1339            taz = np.array((taz*tad/tabs),dtype='float32')
1340            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
1341                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,
1342                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
1343            times[3] += time.time()-t0          #binning
1344    G2fil.G2Print('End integration loops')
1345
1346    t0 = time.time()
1347    #prepare masked arrays of bins with pixels for interpolation setup
1348    H2msk = [ma.array(H2[:-1],mask=np.logical_not(nst)) for nst in NST]
1349#    H0msk = [ma.array(np.divide(h0,nst),mask=np.logical_not(nst)) for nst,h0 in zip(NST,H0)]
1350    H0msk = [ma.array(h0/nst,mask=np.logical_not(nst)) for nst,h0 in zip(NST,H0)]
1351    #make linear interpolators; outside limits give NaN
1352    H0 = []
1353    problemEntries = []
1354    for i,(h0msk,h2msk) in enumerate(zip(H0msk,H2msk)):
1355        try:
1356            h0int = scint.interp1d(h2msk.compressed(),h0msk.compressed(),bounds_error=False)
1357            #do interpolation on all points - fills in the empty bins; leaves others the same
1358            h0 = h0int(H2[:-1])
1359            H0.append(h0)
1360        except ValueError:
1361            problemEntries.append(i)
1362            H0.append(np.zeros(numChans,order='F',dtype=np.float32))
1363    H0 = np.nan_to_num(np.array(H0))       
1364    if 'log(q)' in data.get('binType',''):
1365        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
1366    elif 'q' == data.get('binType','').lower():
1367        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
1368    if Dazm:       
1369        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
1370    else:
1371        H1 = LRazm
1372    if 'SASD' not in data['type']:
1373        H0 *= np.array(G2pwd.Polarization(data['PolaVal'][0],H2[:-1],0.)[0])
1374    H0 /= np.abs(npcosd(H2[:-1]-np.abs(data['det2theta'])))           #parallax correction
1375    if 'SASD' in data['type']:
1376        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
1377    if data['Oblique'][1]:
1378        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
1379    times[4] += time.time()-t0          #cleanup
1380
1381    G2fil.G2Print ('Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
1382        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4]))
1383    G2fil.G2Print ("Elapsed time:","%8.3fs"%(time.time()-tbeg))
1384    G2fil.G2Print ('Integration complete')
1385    if problemEntries:
1386        msg = ""
1387        for i in problemEntries:
1388            if msg: msg += ', '
1389            msg += "{:.2f}".format((H1[i+1]+H1[i])/2.)
1390        print('\nWarning, integrations have no pixels at these azimuth(s)\n\t',msg,'\n')
1391    if returnN:     #As requested by Steven Weigand
1392        return H0,H1,H2,NST,CancelPressed
1393    else:
1394        return H0,H1,H2,CancelPressed
1395   
1396def MakeStrStaRing(ring,Image,Controls):
1397    ellipse = GetEllipse(ring['Dset'],Controls)
1398    pixSize = Controls['pixelSize']
1399    scalex = 1000./pixSize[0]
1400    scaley = 1000./pixSize[1]
1401    Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)[0]).T   #returns x,y,dsp for each point in ring
1402    if len(Ring):
1403        ring['ImxyObs'] = copy.copy(Ring[:2])
1404        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
1405        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
1406        ring['ImtaObs'] = TA
1407        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
1408        Ring[0] = TA[0]
1409        Ring[1] = TA[1]
1410        return Ring,ring
1411    else:
1412        ring['ImxyObs'] = [[],[]]
1413        ring['ImtaObs'] = [[],[]]
1414        ring['ImtaCalc'] = [[],[]]
1415        return [],[]    #bad ring; no points found
1416   
1417def FitStrSta(Image,StrSta,Controls):
1418    'Needs a doc string'
1419   
1420    StaControls = copy.deepcopy(Controls)
1421    phi = StrSta['Sample phi']
1422    wave = Controls['wavelength']
1423    pixelSize = Controls['pixelSize']
1424    scalex = 1000./pixelSize[0]
1425    scaley = 1000./pixelSize[1]
1426    StaType = StrSta['Type']
1427    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1428
1429    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1430        dset = ring['Dset']
1431        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1432        if len(Ring):
1433            ring.update(R)
1434            p0 = ring['Emat']
1435            val,esd,covMat = FitStrain(Ring,p0,dset,wave,phi,StaType)
1436            ring['Emat'] = val
1437            ring['Esig'] = esd
1438            ellipse = FitEllipse(R['ImxyObs'].T)
1439            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1440            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1441            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1442            ringint /= np.mean(ringint)
1443            ring['Ivar'] = np.var(ringint)
1444            ring['covMat'] = covMat
1445            G2fil.G2Print ('Variance in normalized ring intensity: %.3f'%(ring['Ivar']))
1446    CalcStrSta(StrSta,Controls)
1447   
1448def IntStrSta(Image,StrSta,Controls):
1449    StaControls = copy.deepcopy(Controls)
1450    pixelSize = Controls['pixelSize']
1451    scalex = 1000./pixelSize[0]
1452    scaley = 1000./pixelSize[1]
1453    phi = StrSta['Sample phi']
1454    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1455    RingsAI = []
1456    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1457        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1458        if len(Ring):
1459            ellipse = FitEllipse(R['ImxyObs'].T)
1460            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image,5)
1461            XY = np.array(ringxy).T
1462            Th,Azm = GetTthAzm(XY[0],XY[1],Controls)
1463            pola = G2pwd.Polarization(Controls['PolaVal'][0],Th,Azm-90.)[0]     #get pola not dpola
1464            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1465            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1466            ringint /= np.mean(ringint)
1467            ringint /= pola[0]      #just 1st column
1468            G2fil.G2Print (' %s %.3f %s %.3f %s %d'%('d-spacing',ring['Dcalc'],'sig(MRD):',np.sqrt(np.var(ringint)),'# points:',len(ringint)))
1469            RingsAI.append(np.array(list(zip(ringazm,ringint))).T)
1470    return RingsAI
1471   
1472def CalcStrSta(StrSta,Controls):
1473
1474    wave = Controls['wavelength']
1475    phi = StrSta['Sample phi']
1476    StaType = StrSta['Type']
1477    for ring in StrSta['d-zero']:
1478        Eij = ring['Emat']
1479        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
1480        th,azm = ring['ImtaObs']
1481        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
1482        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
1483        if StaType == 'True':
1484            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
1485        else:
1486            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
1487        dmin = np.min(ring['ImtaCalc'][0])
1488        dmax = np.max(ring['ImtaCalc'][0])
1489        if ring.get('fixDset',True):
1490            if abs(Eij[0]) < abs(Eij[2]):         #tension
1491                ring['Dcalc'] = dmin+(dmax-dmin)/4.
1492            else:                       #compression
1493                ring['Dcalc'] = dmin+3.*(dmax-dmin)/4.
1494        else:
1495            ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
1496
1497def calcFij(omg,phi,azm,th):
1498    '''    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
1499
1500    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
1501        90 when perp. to sample surface
1502    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
1503    :param azm: his chi = azimuth around incident beam
1504    :param th:  his theta = theta
1505    '''
1506    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
1507    b = -npcosd(azm)*npcosd(th)
1508    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
1509    d = a*npsind(phi)+b*npcosd(phi)
1510    e = a*npcosd(phi)-b*npsind(phi)
1511    Fij = np.array([
1512        [d**2,d*e,c*d],
1513        [d*e,e**2,c*e],
1514        [c*d,c*e,c**2]])
1515    return -Fij
1516
1517def FitStrain(rings,p0,dset,wave,phi,StaType):
1518    'Needs a doc string'
1519    def StrainPrint(ValSig,dset):
1520        print ('Strain tensor for Dset: %.6f'%(dset))
1521        ptlbls = 'names :'
1522        ptstr =  'values:'
1523        sigstr = 'esds  :'
1524        for name,fmt,value,sig in ValSig:
1525            ptlbls += "%s" % (name.rjust(12))
1526            ptstr += fmt % (value)
1527            if sig:
1528                sigstr += fmt % (sig)
1529            else:
1530                sigstr += 12*' '
1531        print (ptlbls)
1532        print (ptstr)
1533        print (sigstr)
1534       
1535    def strainCalc(p,xyd,dset,wave,phi,StaType):
1536        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1537        dspo,azm,dsp = xyd
1538        th = npasind(wave/(2.0*dspo))
1539        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1540        if StaType == 'True':
1541            dspc = dset*np.exp(V)
1542        else:
1543            dspc = dset*(V+1.)
1544        return dspo-dspc
1545       
1546    names = ['e11','e12','e22']
1547    fmt = ['%12.2f','%12.2f','%12.2f']
1548    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1549    vals = list(result[0])
1550    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1551    covM = result[1]
1552    covMatrix = covM*chisq
1553    sig = list(np.sqrt(chisq*np.diag(result[1])))
1554    ValSig = zip(names,fmt,vals,sig)
1555    StrainPrint(ValSig,dset)
1556    return vals,sig,covMatrix
1557
1558def FitImageSpots(Image,ImMax,ind,pixSize,nxy,spotSize=1.0):
1559   
1560    def calcMean(nxy,pixSize,img):
1561        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1562        gdx = ma.array((gdx-nxy//2)*pixSize[0]/1000.,mask=~ma.getmaskarray(ImBox))
1563        gdy = ma.array((gdy-nxy//2)*pixSize[1]/1000.,mask=~ma.getmaskarray(ImBox))
1564        posx = ma.sum(gdx)/ma.count(gdx)
1565        posy = ma.sum(gdy)/ma.count(gdy)
1566        return posx,posy
1567   
1568    def calcPeak(values,nxy,pixSize,img):
1569        back,mag,px,py,sig = values
1570        peak = np.zeros([nxy,nxy])+back
1571        nor = 1./(2.*np.pi*sig**2)
1572        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1573        gdx = (gdx-nxy//2)*pixSize[0]/1000.
1574        gdy = (gdy-nxy//2)*pixSize[1]/1000.
1575        arg = (gdx-px)**2+(gdy-py)**2       
1576        peak += mag*nor*np.exp(-arg/(2.*sig**2))
1577        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))
1578   
1579    def calc2Peak(values,nxy,pixSize,img):
1580        back,mag,px,py,sigx,sigy,rho = values
1581        peak = np.zeros([nxy,nxy])+back
1582        nor = 1./(2.*np.pi*sigx*sigy*np.sqrt(1.-rho**2))
1583        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1584        gdx = (gdx-nxy//2)*pixSize[0]/1000.
1585        gdy = (gdy-nxy//2)*pixSize[1]/1000.
1586        argnor = -1./(2.*(1.-rho**2))
1587        arg = (gdx-px)**2/sigx**2+(gdy-py)**2/sigy**2-2.*rho*(gdx-px)*(gdy-py)/(sigx*sigy)       
1588        peak += mag*nor*np.exp(argnor*arg)
1589        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))       
1590   
1591    nxy2 = nxy//2
1592    ImBox = Image[ind[1]-nxy2:ind[1]+nxy2+1,ind[0]-nxy2:ind[0]+nxy2+1]
1593    back = np.min(ImBox)
1594    mag = np.sum(ImBox-back)
1595    vals = [back,mag,0.,0.,0.2,0.2,0.]
1596    ImBox = ma.array(ImBox,dtype=float,mask=ImBox>0.75*ImMax)
1597    px = (ind[0]+.5)*pixSize[0]/1000.
1598    py = (ind[1]+.5)*pixSize[1]/1000.
1599    if ma.any(ma.getmaskarray(ImBox)):
1600        vals = calcMean(nxy,pixSize,ImBox)
1601        posx,posy = [px+vals[0],py+vals[1]]
1602        return [posx,posy,spotSize]
1603    else:
1604        result = leastsq(calc2Peak,vals,args=(nxy,pixSize,ImBox),full_output=True)
1605        vals = result[0]
1606        ratio = vals[4]/vals[5]
1607        if 0.5 < ratio < 2.0 and vals[2] < 2. and vals[3] < 2.:
1608            posx,posy = [px+vals[2],py+vals[3]]
1609            return [posx,posy,min(6.*vals[4],spotSize)]
1610        else:
1611            return None
1612   
1613def AutoSpotMask(Image,Masks,Controls,numChans,dlg=None):
1614   
1615    frame = Masks['Frames']
1616    tam = ma.make_mask_none(Image.shape)
1617    if frame:
1618        tam = ma.mask_or(tam,MakeFrameMask(Controls,frame))
1619    LUtth = np.array(Controls['IOtth'])
1620    dtth = (LUtth[1]-LUtth[0])/numChans
1621    esdMul = Masks['SpotMask']['esdMul']
1622    prob = 100.*sc.erf(esdMul/np.sqrt(2.))
1623    print(' Spots greater than %.2f of band intensity are masked'%prob)
1624    mask = ma.make_mask_none(Image.shape)
1625    band = ma.array(Image,mask=ma.nomask)
1626    TA = Make2ThetaAzimuthMap(Controls,(0,Image.shape[0]),(0,Image.shape[1]))[0]    #2-theta array
1627    TThs = np.linspace(LUtth[0],LUtth[1],numChans,False)
1628    for it,TTh in enumerate(TThs):
1629        band.mask = ma.masked_outside(TA,TTh,TTh+dtth).mask+tam
1630        pcmax = np.percentile(band.compressed(),[prob,50.])
1631        mband = ma.masked_greater(band,pcmax[0])
1632        std = ma.std(mband)
1633        anom = ma.masked_greater((band-pcmax[1])/std,esdMul)
1634        mask ^= (anom.mask^band.mask)
1635        if not dlg is None:
1636            GoOn = dlg.Update(it,newmsg='Processed 2-theta rings = %d'%(it))
1637            if not GoOn[0]:
1638                break
1639    return mask
1640
1641def DoPolaCalib(ImageZ,imageData,arcTth):
1642    ''' Determine image polarization by successive integrations with & without preset arc mask.
1643    After initial search, does a set of five with offset azimuth to get mean(std) result.
1644    '''
1645    from scipy.optimize import minimize_scalar
1646    print(' Polarization fit for mask at %.1f deg 2-theta'%arcTth)
1647    data = copy.deepcopy(imageData)
1648    data['IOtth'] = [arcTth-2.,arcTth+2.]
1649    data['fullIntegrate'] = True
1650    data['LRazimuth'] = [0.,360.]
1651    data['outChannels'] = 200
1652    data['outAzimuths'] = 1
1653    Arc = [arcTth,[85.,95.],2.0]
1654    print(' Integration 2-theta test range %.1f - %.1f in 200 steps'%(data['IOtth'][0],data['IOtth'][1]))
1655    print(' Mask azimuth range: %.1f - %.1f'%(Arc[1][0],Arc[1][1]))
1656    Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Frames':[],
1657        'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
1658    AMasks = {'Points':[],'Rings':[],'Arcs':[Arc,],'Polygons':[],'Frames':[],
1659        'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
1660   
1661    def func(p):
1662        p = min(1.0,max(p,0.0))
1663        data['PolaVal'][0] = p
1664        H0 = ImageIntegrate(ImageZ,data,Masks)[0]
1665        H0A = ImageIntegrate(ImageZ,data,AMasks)[0]
1666        M = np.sum(H0)-np.sum(H0A)
1667        print(' Polarization %.4f, fxn: %.1f'%(p,M))
1668        return M**2
1669    time0 = time.time()
1670    res = minimize_scalar(func,bracket=[1.,.999],tol=.0001)
1671    print(res)
1672    pola = min(1.0,max(res.x,.0))
1673   
1674    Pola = []
1675    for arc in [75,80,85,90,95]:
1676        Arc = [arcTth,[arc,arc+10.],2.0]
1677        AMasks = {'Points':[],'Rings':[],'Arcs':[Arc,],'Polygons':[],'Frames':[],
1678            'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
1679        res = minimize_scalar(func,bracket=[pola-.001,pola],tol=.0001)
1680        print(' Mask azimuth range: %.1f - %.1f'%(Arc[1][0],Arc[1][1]))
1681        print(' pola: %.5f'%res.x)
1682        Pola.append(res.x)
1683    Pola = np.array(Pola)
1684    mean = np.mean(Pola)
1685    std = int(10000*np.std(Pola))
1686    print(' Polarization: %.4f(%d)'%(mean,std))
1687    print(' time: %.2fs'%(time.time()-time0))
1688    imageData['PolaVal'][0] = mean
1689   
Note: See TracBrowser for help on using the repository browser.