source: trunk/GSASIIimage.py @ 4685

Last change on this file since 4685 was 4685, checked in by vondreele, 10 months ago

Add gain map calculator

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 68.8 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2020-12-29 18:37:31 +0000 (Tue, 29 Dec 2020) $
5# $Author: vondreele $
6# $Revision: 4685 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 4685 2020-12-29 18:37: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: 4685 $")
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 MakeGainMap(image,Ix,Iy,data,masks,blkSize=128):
1265    import scipy.ndimage.filters as sdif
1266    Iy *= npcosd(Ix[:-1])**3       #undo parallax 3 or 4?
1267    Iy *= (1000./data['distance'])**2    #undo r^2 effect
1268    Iy /= np.array(G2pwd.Polarization(data['PolaVal'][0],Ix[:-1],0.)[0])    #undo polarization
1269    if data['Oblique'][1]:
1270        Iy *= G2pwd.Oblique(data['Oblique'][0],Ix[:-1])     #undo penetration
1271    IyInt = scint.interp1d(Ix[:-1],Iy[0],bounds_error=False)
1272    GainMap = np.zeros_like(image,dtype=float)
1273    #do interpolation on all points - fills in the empty bins; leaves others the same
1274    Nx,Ny = data['size']
1275    nXBlks = (Nx-1)//blkSize+1
1276    nYBlks = (Ny-1)//blkSize+1
1277    for iBlk in range(nYBlks):
1278        iBeg = iBlk*blkSize
1279        iFin = min(iBeg+blkSize,Ny)
1280        for jBlk in range(nXBlks):
1281            jBeg = jBlk*blkSize
1282            jFin = min(jBeg+blkSize,Nx)
1283            TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
1284            Ipix = IyInt(TA[0])
1285            GainMap[iBeg:iFin,jBeg:jFin] = image[iBeg:iFin,jBeg:jFin]/(Ipix*TA[3])
1286    GainMap = np.nan_to_num(GainMap)
1287    GainMap = sdif.gaussian_filter(GainMap,3.,order=0)
1288    return 1./GainMap
1289
1290def ImageIntegrate(image,data,masks,blkSize=128,returnN=False,useTA=None,useMask=None):
1291    'Integrate an image; called from OnIntegrateAll and OnIntegrate in G2imgGUI'    #for q, log(q) bins need data['binType']
1292    import histogram2d as h2d
1293    G2fil.G2Print ('Begin image integration; image range: %d %d'%(np.min(image),np.max(image)))
1294    CancelPressed = False
1295    LUtth = np.array(data['IOtth'])
1296    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
1297    numAzms = data['outAzimuths']
1298    numChans = (data['outChannels']//4)*4
1299    Dazm = (LRazm[1]-LRazm[0])/numAzms
1300    if '2-theta' in data.get('binType','2-theta'):
1301        lutth = LUtth               
1302    elif 'log(q)' in data['binType']:
1303        lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
1304    elif 'q' == data['binType'].lower():
1305        lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
1306    dtth = (lutth[1]-lutth[0])/numChans
1307    muT = data.get('SampleAbs',[0.0,''])[0]
1308    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
1309        data['DetDepth'] /= data['distance']
1310    if 'SASD' in data['type']:
1311        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
1312    Masks = copy.deepcopy(masks)
1313    Masks['Points'] = np.array(Masks['Points']).T           #get spots as X,Y,R arrays
1314    if np.any(masks['Points']):
1315        Masks['Points'][2] = np.square(Masks['Points'][2]/2.)
1316    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
1317    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
1318    H2 = np.linspace(lutth[0],lutth[1],numChans+1)
1319    Nx,Ny = data['size']
1320    nXBlks = (Nx-1)//blkSize+1
1321    nYBlks = (Ny-1)//blkSize+1
1322    tbeg = time.time()
1323    times = [0,0,0,0,0]
1324    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
1325    for iBlk in range(nYBlks):
1326        iBeg = iBlk*blkSize
1327        iFin = min(iBeg+blkSize,Ny)
1328        for jBlk in range(nXBlks):
1329            jBeg = jBlk*blkSize
1330            jFin = min(jBeg+blkSize,Nx)
1331            # next is most expensive step!
1332
1333            t0 = time.time()
1334            if useTA:
1335                TA = useTA[iBlk][jBlk]
1336            else:
1337                TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
1338            times[1] += time.time()-t0      #xy->th,azm
1339
1340            t0 = time.time()
1341            if useMask:
1342                tam = useMask[iBlk][jBlk]
1343            else:
1344                tam = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp)
1345            times[0] += time.time()-t0      #apply masks
1346
1347            t0 = time.time()
1348            Block = image[iBeg:iFin,jBeg:jFin]          #apply image spotmask here
1349            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(Masks,TA,tam,Block)    #and apply masks
1350            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
1351            tax = np.where(tax < LRazm[0],tax+360.,tax)
1352            if data.get('SampleAbs',[0.0,''])[1]:
1353                if 'Cylind' in data['SampleShape']:
1354                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))      #adjust for additional thickness off sample normal
1355                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
1356                elif 'Fixed' in data['SampleShape']:    #assumes flat plate sample normal to beam
1357                    tabs = G2pwd.Absorb('Fixed',muT,tay)
1358            if 'log(q)' in data.get('binType',''):
1359                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
1360            elif 'q' == data.get('binType','').lower():
1361                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
1362            times[2] += time.time()-t0          #fill map
1363
1364            t0 = time.time()
1365            taz = np.array((taz*tad/tabs),dtype='float32')
1366            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
1367                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,
1368                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
1369            times[3] += time.time()-t0          #binning
1370    G2fil.G2Print('End integration loops')
1371
1372    t0 = time.time()
1373    #prepare masked arrays of bins with pixels for interpolation setup
1374    H2msk = [ma.array(H2[:-1],mask=np.logical_not(nst)) for nst in NST]
1375#    H0msk = [ma.array(np.divide(h0,nst),mask=np.logical_not(nst)) for nst,h0 in zip(NST,H0)]
1376    H0msk = [ma.array(h0/nst,mask=np.logical_not(nst)) for nst,h0 in zip(NST,H0)]
1377    #make linear interpolators; outside limits give NaN
1378    H0 = []
1379    problemEntries = []
1380    for i,(h0msk,h2msk) in enumerate(zip(H0msk,H2msk)):
1381        try:
1382            h0int = scint.interp1d(h2msk.compressed(),h0msk.compressed(),bounds_error=False)
1383            #do interpolation on all points - fills in the empty bins; leaves others the same
1384            h0 = h0int(H2[:-1])
1385            H0.append(h0)
1386        except ValueError:
1387            problemEntries.append(i)
1388            H0.append(np.zeros(numChans,order='F',dtype=np.float32))
1389    H0 = np.nan_to_num(np.array(H0))       
1390    if 'log(q)' in data.get('binType',''):
1391        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
1392    elif 'q' == data.get('binType','').lower():
1393        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
1394    if Dazm:       
1395        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
1396    else:
1397        H1 = LRazm
1398    if 'SASD' not in data['type']:
1399        H0 *= np.array(G2pwd.Polarization(data['PolaVal'][0],H2[:-1],0.)[0])
1400    H0 /= np.abs(npcosd(H2[:-1]-np.abs(data['det2theta'])))           #parallax correction
1401    H0 *= (data['distance']/1000.)**2    #remove r^2 effect
1402    if 'SASD' in data['type']:
1403        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
1404    if data['Oblique'][1]:
1405        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
1406    times[4] += time.time()-t0          #cleanup
1407
1408    G2fil.G2Print ('Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
1409        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4]))
1410    G2fil.G2Print ("Elapsed time:","%8.3fs"%(time.time()-tbeg))
1411    G2fil.G2Print ('Integration complete')
1412    if problemEntries:
1413        msg = ""
1414        for i in problemEntries:
1415            if msg: msg += ', '
1416            msg += "{:.2f}".format((H1[i+1]+H1[i])/2.)
1417        print('\nWarning, integrations have no pixels at these azimuth(s)\n\t',msg,'\n')
1418    if returnN:     #As requested by Steven Weigand
1419        return H0,H1,H2,NST,CancelPressed
1420    else:
1421        return H0,H1,H2,CancelPressed
1422   
1423def MakeStrStaRing(ring,Image,Controls):
1424    ellipse = GetEllipse(ring['Dset'],Controls)
1425    pixSize = Controls['pixelSize']
1426    scalex = 1000./pixSize[0]
1427    scaley = 1000./pixSize[1]
1428    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
1429    if len(Ring):
1430        ring['ImxyObs'] = copy.copy(Ring[:2])
1431        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
1432        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
1433        ring['ImtaObs'] = TA
1434        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
1435        Ring[0] = TA[0]
1436        Ring[1] = TA[1]
1437        return Ring,ring
1438    else:
1439        ring['ImxyObs'] = [[],[]]
1440        ring['ImtaObs'] = [[],[]]
1441        ring['ImtaCalc'] = [[],[]]
1442        return [],[]    #bad ring; no points found
1443   
1444def FitStrSta(Image,StrSta,Controls):
1445    'Needs a doc string'
1446   
1447    StaControls = copy.deepcopy(Controls)
1448    phi = StrSta['Sample phi']
1449    wave = Controls['wavelength']
1450    pixelSize = Controls['pixelSize']
1451    scalex = 1000./pixelSize[0]
1452    scaley = 1000./pixelSize[1]
1453    StaType = StrSta['Type']
1454    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1455
1456    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1457        dset = ring['Dset']
1458        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1459        if len(Ring):
1460            ring.update(R)
1461            p0 = ring['Emat']
1462            val,esd,covMat = FitStrain(Ring,p0,dset,wave,phi,StaType)
1463            ring['Emat'] = val
1464            ring['Esig'] = esd
1465            ellipse = FitEllipse(R['ImxyObs'].T)
1466            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1467            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1468            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1469            ringint /= np.mean(ringint)
1470            ring['Ivar'] = np.var(ringint)
1471            ring['covMat'] = covMat
1472            G2fil.G2Print ('Variance in normalized ring intensity: %.3f'%(ring['Ivar']))
1473    CalcStrSta(StrSta,Controls)
1474   
1475def IntStrSta(Image,StrSta,Controls):
1476    StaControls = copy.deepcopy(Controls)
1477    pixelSize = Controls['pixelSize']
1478    scalex = 1000./pixelSize[0]
1479    scaley = 1000./pixelSize[1]
1480    phi = StrSta['Sample phi']
1481    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1482    RingsAI = []
1483    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1484        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1485        if len(Ring):
1486            ellipse = FitEllipse(R['ImxyObs'].T)
1487            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image,5)
1488            XY = np.array(ringxy).T
1489            Th,Azm = GetTthAzm(XY[0],XY[1],Controls)
1490            pola = G2pwd.Polarization(Controls['PolaVal'][0],Th,Azm-90.)[0]     #get pola not dpola
1491            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1492            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1493            ringint /= np.mean(ringint)
1494            ringint /= pola[0]      #just 1st column
1495            G2fil.G2Print (' %s %.3f %s %.3f %s %d'%('d-spacing',ring['Dcalc'],'sig(MRD):',np.sqrt(np.var(ringint)),'# points:',len(ringint)))
1496            RingsAI.append(np.array(list(zip(ringazm,ringint))).T)
1497    return RingsAI
1498   
1499def CalcStrSta(StrSta,Controls):
1500
1501    wave = Controls['wavelength']
1502    phi = StrSta['Sample phi']
1503    StaType = StrSta['Type']
1504    for ring in StrSta['d-zero']:
1505        Eij = ring['Emat']
1506        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
1507        th,azm = ring['ImtaObs']
1508        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
1509        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
1510        if StaType == 'True':
1511            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
1512        else:
1513            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
1514        dmin = np.min(ring['ImtaCalc'][0])
1515        dmax = np.max(ring['ImtaCalc'][0])
1516        if ring.get('fixDset',True):
1517            if abs(Eij[0]) < abs(Eij[2]):         #tension
1518                ring['Dcalc'] = dmin+(dmax-dmin)/4.
1519            else:                       #compression
1520                ring['Dcalc'] = dmin+3.*(dmax-dmin)/4.
1521        else:
1522            ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
1523
1524def calcFij(omg,phi,azm,th):
1525    '''    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
1526
1527    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
1528        90 when perp. to sample surface
1529    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
1530    :param azm: his chi = azimuth around incident beam
1531    :param th:  his theta = theta
1532    '''
1533    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
1534    b = -npcosd(azm)*npcosd(th)
1535    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
1536    d = a*npsind(phi)+b*npcosd(phi)
1537    e = a*npcosd(phi)-b*npsind(phi)
1538    Fij = np.array([
1539        [d**2,d*e,c*d],
1540        [d*e,e**2,c*e],
1541        [c*d,c*e,c**2]])
1542    return -Fij
1543
1544def FitStrain(rings,p0,dset,wave,phi,StaType):
1545    'Needs a doc string'
1546    def StrainPrint(ValSig,dset):
1547        print ('Strain tensor for Dset: %.6f'%(dset))
1548        ptlbls = 'names :'
1549        ptstr =  'values:'
1550        sigstr = 'esds  :'
1551        for name,fmt,value,sig in ValSig:
1552            ptlbls += "%s" % (name.rjust(12))
1553            ptstr += fmt % (value)
1554            if sig:
1555                sigstr += fmt % (sig)
1556            else:
1557                sigstr += 12*' '
1558        print (ptlbls)
1559        print (ptstr)
1560        print (sigstr)
1561       
1562    def strainCalc(p,xyd,dset,wave,phi,StaType):
1563        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1564        dspo,azm,dsp = xyd
1565        th = npasind(wave/(2.0*dspo))
1566        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1567        if StaType == 'True':
1568            dspc = dset*np.exp(V)
1569        else:
1570            dspc = dset*(V+1.)
1571        return dspo-dspc
1572       
1573    names = ['e11','e12','e22']
1574    fmt = ['%12.2f','%12.2f','%12.2f']
1575    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1576    vals = list(result[0])
1577    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1578    covM = result[1]
1579    covMatrix = covM*chisq
1580    sig = list(np.sqrt(chisq*np.diag(result[1])))
1581    ValSig = zip(names,fmt,vals,sig)
1582    StrainPrint(ValSig,dset)
1583    return vals,sig,covMatrix
1584
1585def FitImageSpots(Image,ImMax,ind,pixSize,nxy,spotSize=1.0):
1586   
1587    def calcMean(nxy,pixSize,img):
1588        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1589        gdx = ma.array((gdx-nxy//2)*pixSize[0]/1000.,mask=~ma.getmaskarray(ImBox))
1590        gdy = ma.array((gdy-nxy//2)*pixSize[1]/1000.,mask=~ma.getmaskarray(ImBox))
1591        posx = ma.sum(gdx)/ma.count(gdx)
1592        posy = ma.sum(gdy)/ma.count(gdy)
1593        return posx,posy
1594   
1595    def calcPeak(values,nxy,pixSize,img):
1596        back,mag,px,py,sig = values
1597        peak = np.zeros([nxy,nxy])+back
1598        nor = 1./(2.*np.pi*sig**2)
1599        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1600        gdx = (gdx-nxy//2)*pixSize[0]/1000.
1601        gdy = (gdy-nxy//2)*pixSize[1]/1000.
1602        arg = (gdx-px)**2+(gdy-py)**2       
1603        peak += mag*nor*np.exp(-arg/(2.*sig**2))
1604        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))
1605   
1606    def calc2Peak(values,nxy,pixSize,img):
1607        back,mag,px,py,sigx,sigy,rho = values
1608        peak = np.zeros([nxy,nxy])+back
1609        nor = 1./(2.*np.pi*sigx*sigy*np.sqrt(1.-rho**2))
1610        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1611        gdx = (gdx-nxy//2)*pixSize[0]/1000.
1612        gdy = (gdy-nxy//2)*pixSize[1]/1000.
1613        argnor = -1./(2.*(1.-rho**2))
1614        arg = (gdx-px)**2/sigx**2+(gdy-py)**2/sigy**2-2.*rho*(gdx-px)*(gdy-py)/(sigx*sigy)       
1615        peak += mag*nor*np.exp(argnor*arg)
1616        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))       
1617   
1618    nxy2 = nxy//2
1619    ImBox = Image[ind[1]-nxy2:ind[1]+nxy2+1,ind[0]-nxy2:ind[0]+nxy2+1]
1620    back = np.min(ImBox)
1621    mag = np.sum(ImBox-back)
1622    vals = [back,mag,0.,0.,0.2,0.2,0.]
1623    ImBox = ma.array(ImBox,dtype=float,mask=ImBox>0.75*ImMax)
1624    px = (ind[0]+.5)*pixSize[0]/1000.
1625    py = (ind[1]+.5)*pixSize[1]/1000.
1626    if ma.any(ma.getmaskarray(ImBox)):
1627        vals = calcMean(nxy,pixSize,ImBox)
1628        posx,posy = [px+vals[0],py+vals[1]]
1629        return [posx,posy,spotSize]
1630    else:
1631        result = leastsq(calc2Peak,vals,args=(nxy,pixSize,ImBox),full_output=True)
1632        vals = result[0]
1633        ratio = vals[4]/vals[5]
1634        if 0.5 < ratio < 2.0 and vals[2] < 2. and vals[3] < 2.:
1635            posx,posy = [px+vals[2],py+vals[3]]
1636            return [posx,posy,min(6.*vals[4],spotSize)]
1637        else:
1638            return None
1639   
1640def AutoSpotMask(Image,Masks,Controls,numChans,dlg=None):
1641   
1642    frame = Masks['Frames']
1643    tam = ma.make_mask_none(Image.shape)
1644    if frame:
1645        tam = ma.mask_or(tam,MakeFrameMask(Controls,frame))
1646    LUtth = np.array(Controls['IOtth'])
1647    dtth = (LUtth[1]-LUtth[0])/numChans
1648    esdMul = Masks['SpotMask']['esdMul']
1649    prob = 100.*sc.erf(esdMul/np.sqrt(2.))
1650    print(' Spots greater than %.2f of band intensity are masked'%prob)
1651    mask = ma.make_mask_none(Image.shape)
1652    band = ma.array(Image,mask=ma.nomask)
1653    TA = Make2ThetaAzimuthMap(Controls,(0,Image.shape[0]),(0,Image.shape[1]))[0]    #2-theta array
1654    TThs = np.linspace(LUtth[0],LUtth[1],numChans,False)
1655    for it,TTh in enumerate(TThs):
1656        band.mask = ma.masked_outside(TA,TTh,TTh+dtth).mask+tam
1657        pcmax = np.percentile(band.compressed(),[prob,50.])
1658        mband = ma.masked_greater(band,pcmax[0])
1659        std = ma.std(mband)
1660        anom = ma.masked_greater((band-pcmax[1])/std,esdMul)
1661        mask ^= (anom.mask^band.mask)
1662        if not dlg is None:
1663            GoOn = dlg.Update(it,newmsg='Processed 2-theta rings = %d'%(it))
1664            if not GoOn[0]:
1665                break
1666    return mask
1667
1668def DoPolaCalib(ImageZ,imageData,arcTth):
1669    ''' Determine image polarization by successive integrations with & without preset arc mask.
1670    After initial search, does a set of five with offset azimuth to get mean(std) result.
1671    '''
1672    from scipy.optimize import minimize_scalar
1673    print(' Polarization fit for mask at %.1f deg 2-theta'%arcTth)
1674    data = copy.deepcopy(imageData)
1675    data['IOtth'] = [arcTth-2.,arcTth+2.]
1676    data['fullIntegrate'] = True
1677    data['LRazimuth'] = [0.,360.]
1678    data['outChannels'] = 200
1679    data['outAzimuths'] = 1
1680    Arc = [arcTth,[85.,95.],2.0]
1681    print(' Integration 2-theta test range %.1f - %.1f in 200 steps'%(data['IOtth'][0],data['IOtth'][1]))
1682    print(' Mask azimuth range: %.1f - %.1f'%(Arc[1][0],Arc[1][1]))
1683    Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Frames':[],
1684        'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
1685    AMasks = {'Points':[],'Rings':[],'Arcs':[Arc,],'Polygons':[],'Frames':[],
1686        'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
1687   
1688    def func(p):
1689        p = min(1.0,max(p,0.0))
1690        data['PolaVal'][0] = p
1691        H0 = ImageIntegrate(ImageZ,data,Masks)[0]
1692        H0A = ImageIntegrate(ImageZ,data,AMasks)[0]
1693        M = np.sum(H0)-np.sum(H0A)
1694        print(' Polarization %.4f, fxn: %.1f'%(p,M))
1695        return M**2
1696    time0 = time.time()
1697    res = minimize_scalar(func,bracket=[1.,.999],tol=.0001)
1698    print(res)
1699    pola = min(1.0,max(res.x,.0))
1700   
1701    Pola = []
1702    for arc in [75,80,85,90,95]:
1703        Arc = [arcTth,[arc,arc+10.],2.0]
1704        AMasks = {'Points':[],'Rings':[],'Arcs':[Arc,],'Polygons':[],'Frames':[],
1705            'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
1706        res = minimize_scalar(func,bracket=[pola-.001,pola],tol=.0001)
1707        print(' Mask azimuth range: %.1f - %.1f'%(Arc[1][0],Arc[1][1]))
1708        print(' pola: %.5f'%res.x)
1709        Pola.append(res.x)
1710    Pola = np.array(Pola)
1711    mean = np.mean(Pola)
1712    std = int(10000*np.std(Pola))
1713    print(' Polarization: %.4f(%d)'%(mean,std))
1714    print(' time: %.2fs'%(time.time()-time0))
1715    imageData['PolaVal'][0] = mean
1716   
Note: See TracBrowser for help on using the repository browser.