source: trunk/GSASIIimage.py @ 4587

Last change on this file since 4587 was 4587, checked in by vondreele, 13 months ago

New version of GetTthAzmDsp? to use det2theta

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 65.8 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2020-10-08 14:59:00 +0000 (Thu, 08 Oct 2020) $
5# $Author: vondreele $
6# $Revision: 4587 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 4587 2020-10-08 14:59:00Z 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: 4587 $")
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)
475    it seems to be only used in plotting
476    '''   
477    elcent,phi,radii = GetEllipse(dsp,data)
478    phi = data['rotation']-90.          #to give rotation of major axis
479    tilt = data['tilt']
480    dist = data['distance']
481    cent = data['center']
482    tth = 2.0*asind(data['wavelength']/(2.*dsp))
483    stth = sind(tth)
484    cosb = cosd(tilt)
485    if radii[0] > 0.:
486        sinb = sind(tilt)
487        tanb = tand(tilt)
488        fplus = dist*tanb*stth/(cosb+stth)
489        fminus = dist*tanb*stth/(cosb-stth)
490        zdis = (fplus-fminus)/2.
491        rsqplus = radii[0]**2+radii[1]**2
492        rsqminus = radii[0]**2-radii[1]**2
493        R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus
494        Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm-phi)**2)
495        P = 2.*radii[0]**2*zdis*cosd(azm-phi)
496        radius = (P+Q)/R
497        xy = np.array([radius*cosd(azm),radius*sind(azm)])
498        xy += cent
499    else:   #hyperbola - both branches (one is way off screen!)
500        sinb = abs(sind(tilt))
501        tanb = abs(tand(tilt))
502        f = dist*tanb*stth/(cosb+stth)
503        v = dist*(tanb+tand(tth-abs(tilt)))
504        delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
505        ecc = (v-f)/(delt-v)
506        R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm))
507        if tilt > 0.:
508            offset = 2.*radii[1]*ecc+f      #select other branch
509            xy = [-R*cosd(azm)-offset,-R*sind(azm)]
510        else:
511            offset = -f
512            xy = [-R*cosd(azm)-offset,R*sind(azm)]
513        xy = -np.array([xy[0]*cosd(phi)+xy[1]*sind(phi),xy[0]*sind(phi)-xy[1]*cosd(phi)])
514        xy += cent
515    return xy
516   
517def GetDetXYfromThAzm(Th,Azm,data):
518    '''Computes a detector position from a 2theta angle and an azimultal
519    angle (both in degrees) - apparently not used!
520    '''
521    dsp = data['wavelength']/(2.0*npsind(Th))   
522    return GetDetectorXY(dsp,Azm,data)
523                   
524def GetTthAzmDsp2(x,y,data): #expensive
525    '''Computes a 2theta, etc. from a detector position and calibration constants - checked
526    OK for ellipses & hyperbola.
527
528    :returns: np.array(tth,azm,G,dsp) where tth is 2theta, azm is the azimutal angle,
529       G is ? and dsp is the d-space
530    '''
531    wave = data['wavelength']
532    cent = data['center']
533    tilt = data['tilt']
534    dist = data['distance']/cosd(tilt)
535    x0 = dist*tand(tilt)
536    phi = data['rotation']
537    dep = data.get('DetDepth',0.)
538    azmthoff = data['azmthOff']
539    dx = np.array(x-cent[0],dtype=np.float32)
540    dy = np.array(y-cent[1],dtype=np.float32)
541    D = ((dx-x0)**2+dy**2+dist**2)      #sample to pixel distance
542    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
543    X = np.dot(X,makeMat(phi,2))
544    Z = np.dot(X,makeMat(tilt,0)).T[2]
545    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
546    dxy = peneCorr(tth,dep,dist)
547    DX = dist-Z+dxy
548    DY = np.sqrt(dx**2+dy**2-Z**2)
549    tth = npatan2d(DY,DX) 
550    dsp = wave/(2.*npsind(tth/2.))
551    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
552    G = D/dist**2       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
553    return np.array([tth,azm,G,dsp])
554   
555def GetTthAzmDsp(x,y,data): #expensive
556    '''Computes a 2theta, etc. from a detector position and calibration constants - checked
557    OK for ellipses & hyperbola.
558
559    :returns: np.array(tth,azm,G,dsp) where tth is 2theta, azm is the azimutal angle,
560       G is ? and dsp is the d-space
561    '''
562   
563    def costth(xyz):
564        u = xyz/nl.norm(xyz,axis=-1)[:,:,nxs]
565        return np.dot(u,np.array([0.,0.,1.]))
566       
567#zero detector 2-theta: tested with tilted images - perfect integrations
568    wave = data['wavelength']
569    dx = x-data['center'][0]
570    dy = y-data['center'][1]
571    tilt = data['tilt']
572    dist = data['distance']/npcosd(tilt)    #sample-beam intersection point
573    T = makeMat(tilt,0)
574    R = makeMat(data['rotation'],2)
575    MN = np.inner(R,np.inner(R,T))
576    dxyz0 = np.inner(np.dstack([dx,dy,np.zeros_like(dx)]),MN)    #correct for 45 deg tilt
577    dxyz0 += np.array([0.,0.,dist])
578    if data['DetDepth']:
579        ctth0 = costth(dxyz0)
580        tth0 = npacosd(ctth0)
581        dzp = peneCorr(tth0,data['DetDepth'],dist)
582        dxyz0[:,:,2] += dzp
583#non zero detector 2-theta:
584    if data['det2theta']:       
585        tthMat = makeMat(data['det2theta'],1)
586        dxyz = np.inner(dxyz0,tthMat.T)
587    else:
588        dxyz = dxyz0
589    ctth = costth(dxyz)
590    tth = npacosd(ctth)
591    dsp = wave/(2.*npsind(tth/2.))
592    azm = (npatan2d(dxyz[:,:,1],dxyz[:,:,0])+data['azmthOff']+720.)%360.       
593# G-calculation       
594    x0 = data['distance']*nptand(tilt)
595    x0x = x0*npcosd(data['rotation'])
596    x0y = x0*npsind(data['rotation'])
597    distsq = data['distance']**2
598    G = ((dx-x0x)**2+(dy-x0y)**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
599    return np.array([tth,azm,G,dsp])
600   
601def GetTth(x,y,data):
602    'Give 2-theta value for detector x,y position; calibration info in data'
603    return GetTthAzmDsp(x,y,data)[0]
604   
605def GetTthAzm(x,y,data):
606    'Give 2-theta, azimuth values for detector x,y position; calibration info in data'
607    return GetTthAzmDsp(x,y,data)[0:2]
608   
609def GetTthAzmG2(x,y,data):
610    '''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
611     calibration info in data - only used in integration - old version
612    '''
613    'Needs a doc string - checked OK for ellipses & hyperbola'
614    tilt = data['tilt']
615    dist = data['distance']/npcosd(tilt)
616    MN = -np.inner(makeMat(data['rotation'],2),makeMat(tilt,0))
617    dx = x-data['center'][0]
618    dy = y-data['center'][1]
619    dz = np.dot(np.dstack([dx.T,dy.T,np.zeros_like(dx.T)]),MN).T[2]   
620    xyZ = dx**2+dy**2-dz**2   
621    tth0 = npatand(np.sqrt(xyZ)/(dist-dz))
622    dzp = peneCorr(tth0,data['DetDepth'],dist)
623    tth = npatan2d(np.sqrt(xyZ),dist-dz+dzp) 
624    azm = (npatan2d(dy,dx)+data['azmthOff']+720.)%360.
625   
626    distsq = data['distance']**2
627    x0 = data['distance']*nptand(tilt)
628    x0x = x0*npcosd(data['rotation'])
629    x0y = x0*npsind(data['rotation'])
630    G = ((dx-x0x)**2+(dy-x0y)**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
631    return tth,azm,G
632
633def GetTthAzmG(x,y,data):
634    '''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
635     calibration info in data - only used in integration
636     checked OK for ellipses & hyperbola
637     This is the slow step in image integration
638     '''     
639    def costth(xyz):
640        u = xyz/nl.norm(xyz,axis=-1)[:,:,nxs]
641        return np.dot(u,np.array([0.,0.,1.]))
642       
643#zero detector 2-theta: tested with tilted images - perfect integrations
644    dx = x-data['center'][0]
645    dy = y-data['center'][1]
646    tilt = data['tilt']
647    dist = data['distance']/npcosd(tilt)    #sample-beam intersection point
648    T = makeMat(tilt,0)
649    R = makeMat(data['rotation'],2)
650    MN = np.inner(R,np.inner(R,T))
651    dxyz0 = np.inner(np.dstack([dx,dy,np.zeros_like(dx)]),MN)    #correct for 45 deg tilt
652    dxyz0 += np.array([0.,0.,dist])
653    if data['DetDepth']:
654        ctth0 = costth(dxyz0)
655        tth0 = npacosd(ctth0)
656        dzp = peneCorr(tth0,data['DetDepth'],dist)
657        dxyz0[:,:,2] += dzp
658#non zero detector 2-theta:
659    if data['det2theta']:       
660        tthMat = makeMat(data['det2theta'],1)
661        dxyz = np.inner(dxyz0,tthMat.T)
662    else:
663        dxyz = dxyz0
664    ctth = costth(dxyz)
665    tth = npacosd(ctth)
666    azm = (npatan2d(dxyz[:,:,1],dxyz[:,:,0])+data['azmthOff']+720.)%360.       
667# G-calculation       
668    x0 = data['distance']*nptand(tilt)
669    x0x = x0*npcosd(data['rotation'])
670    x0y = x0*npsind(data['rotation'])
671    distsq = data['distance']**2
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 GetDsp(x,y,data):
676    'Give d-spacing value for detector x,y position; calibration info in data'
677    return GetTthAzmDsp(x,y,data)[3]
678       
679def GetAzm(x,y,data):
680    'Give azimuth value for detector x,y position; calibration info in data'
681    return GetTthAzmDsp(x,y,data)[1]
682   
683def meanAzm(a,b):
684    AZM = lambda a,b: npacosd(0.5*(npsind(2.*b)-npsind(2.*a))/(np.pi*(b-a)/180.))/2.
685    azm = AZM(a,b)
686#    quad = int((a+b)/180.)
687#    if quad == 1:
688#        azm = 180.-azm
689#    elif quad == 2:
690#        azm += 180.
691#    elif quad == 3:
692#        azm = 360-azm
693    return azm     
694       
695def ImageCompress(image,scale):
696    ''' Reduces size of image by selecting every n'th point
697    param: image array: original image
698    param: scale int: intervsl between selected points
699    returns: array: reduced size image
700    '''
701    if scale == 1:
702        return image
703    else:
704        return image[::scale,::scale]
705       
706def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
707    'Needs a doc string'
708    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
709    curr = np.array([dist,x,y])
710    return abs(avg-curr)/avg < .02
711
712def GetLineScan(image,data):
713    Nx,Ny = data['size']
714    pixelSize = data['pixelSize']
715    scalex = 1000./pixelSize[0]         #microns --> 1/mm
716    scaley = 1000./pixelSize[1]
717    wave = data['wavelength']
718    numChans = data['outChannels']
719    LUtth = np.array(data['IOtth'],dtype=np.float)
720    azm = data['linescan'][1]-data['azmthOff']
721    Tx = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
722    Ty = np.zeros_like(Tx)
723    dsp = wave/(2.0*npsind(Tx/2.0))
724    xy = np.array([GetDetectorXY(d,azm,data) for d in dsp]).T
725    xy[1] *= scalex
726    xy[0] *= scaley
727    xy = np.array(xy,dtype=int)
728    Xpix = ma.masked_outside(xy[1],0,Ny-1)
729    Ypix = ma.masked_outside(xy[0],0,Nx-1)
730    xpix = Xpix[~(Xpix.mask+Ypix.mask)].compressed()
731    ypix = Ypix[~(Xpix.mask+Ypix.mask)].compressed()
732    Ty = image[xpix,ypix]
733    Tx = ma.array(Tx,mask=Xpix.mask+Ypix.mask).compressed()
734    return [Tx,Ty]
735
736def EdgeFinder(image,data):
737    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
738    Not currently used but might be useful in future?
739    '''
740    import numpy.ma as ma
741    Nx,Ny = data['size']
742    pixelSize = data['pixelSize']
743    edgemin = data['edgemin']
744    scalex = pixelSize[0]/1000.
745    scaley = pixelSize[1]/1000.   
746    tay,tax = np.mgrid[0:Nx,0:Ny]
747    tax = np.asfarray(tax*scalex,dtype=np.float32)
748    tay = np.asfarray(tay*scaley,dtype=np.float32)
749    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
750    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
751    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
752    return zip(tax,tay)
753   
754def MakeFrameMask(data,frame):
755    import polymask as pm
756    pixelSize = data['pixelSize']
757    scalex = pixelSize[0]/1000.
758    scaley = pixelSize[1]/1000.
759    blkSize = 512
760    Nx,Ny = data['size']
761    nXBlks = (Nx-1)//blkSize+1
762    nYBlks = (Ny-1)//blkSize+1
763    tam = ma.make_mask_none(data['size'])
764    for iBlk in range(nXBlks):
765        iBeg = iBlk*blkSize
766        iFin = min(iBeg+blkSize,Nx)
767        for jBlk in range(nYBlks):
768            jBeg = jBlk*blkSize
769            jFin = min(jBeg+blkSize,Ny)               
770            nI = iFin-iBeg
771            nJ = jFin-jBeg
772            tax,tay = np.mgrid[iBeg+0.5:iFin+.5,jBeg+.5:jFin+.5]         #bin centers not corners
773            tax = np.asfarray(tax*scalex,dtype=np.float32)
774            tay = np.asfarray(tay*scaley,dtype=np.float32)
775            tamp = ma.make_mask_none((1024*1024))
776            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
777                tay.flatten(),len(frame),frame,tamp)[:nI*nJ])^True  #switch to exclude around frame
778            if tamp.shape:
779                tamp = np.reshape(tamp[:nI*nJ],(nI,nJ))
780                tam[iBeg:iFin,jBeg:jFin] = ma.mask_or(tamp[0:nI,0:nJ],tam[iBeg:iFin,jBeg:jFin])
781            else:
782                tam[iBeg:iFin,jBeg:jFin] = True
783    return tam.T
784   
785def ImageRecalibrate(G2frame,ImageZ,data,masks,getRingsOnly=False):
786    '''Called to repeat the calibration on an image, usually called after
787    calibration is done initially to improve the fit.
788
789    :param G2frame: The top-level GSAS-II frame or None, to skip plotting
790    :param np.Array ImageZ: the image to calibrate
791    :param dict data: the Controls dict for the image
792    :param dict masks: a dict with masks
793    :returns: a list containing vals,varyList,sigList,parmDict,covar or rings
794      (with an array of x, y, and d-space values) if getRingsOnly is True
795      or an empty list, in case of an error
796    '''
797    import ImageCalibrants as calFile
798    if not getRingsOnly:
799        G2fil.G2Print ('Image recalibration:')
800    time0 = time.time()
801    pixelSize = data['pixelSize']
802    scalex = 1000./pixelSize[0]
803    scaley = 1000./pixelSize[1]
804    pixLimit = data['pixLimit']
805    cutoff = data['cutoff']
806    data['rings'] = []
807    data['ellipses'] = []
808    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
809        data['DetDepth'] /= data['distance']
810    if not data['calibrant']:
811        G2fil.G2Print ('warning: no calibration material selected')
812        return []   
813    skip = data['calibskip']
814    dmin = data['calibdmin']
815    if data['calibrant'] not in calFile.Calibrants:
816        G2fil.G2Print('Warning: %s not in local copy of image calibrants file'%data['calibrant'])
817        return []
818    calibrant = calFile.Calibrants[data['calibrant']]
819    Bravais,SGs,Cells = calibrant[:3]
820    HKL = []
821    for bravais,sg,cell in zip(Bravais,SGs,Cells):
822        A = G2lat.cell2A(cell)
823        if sg:
824            SGData = G2spc.SpcGroup(sg)[1]
825            hkl = G2pwd.getHKLpeak(dmin,SGData,A,Inst=None,nodup=True)
826            HKL += list(hkl)
827        else:
828            hkl = G2lat.GenHBravais(dmin,bravais,A)
829            HKL += list(hkl)
830    if len(calibrant) > 5:
831        absent = calibrant[5]
832    else:
833        absent = ()
834    HKL = G2lat.sortHKLd(HKL,True,False)
835    varyList = [item for item in data['varyList'] if data['varyList'][item]]
836    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
837        'setdist':data.get('setdist',data['distance']),
838        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
839    Found = False
840    wave = data['wavelength']
841    frame = masks['Frames']
842    tam = ma.make_mask_none(ImageZ.shape)
843    if frame:
844        tam = ma.mask_or(tam,MakeFrameMask(data,frame))
845    for iH,H in enumerate(HKL):
846        if debug:   print (H) 
847        dsp = H[3]
848        tth = 2.0*asind(wave/(2.*dsp))
849        if tth+abs(data['tilt']) > 90.:
850            G2fil.G2Print ('next line is a hyperbola - search stopped')
851            break
852        ellipse = GetEllipse(dsp,data)
853        if iH not in absent and iH >= skip:
854            Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
855        else:
856            Ring = makeRing(dsp,ellipse,pixLimit,1000.0,scalex,scaley,ma.array(ImageZ,mask=tam))[0]
857        if Ring:
858            if iH not in absent and iH >= skip:
859                data['rings'].append(np.array(Ring))
860            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
861            Found = True
862        elif not Found:         #skipping inner rings, keep looking until ring found
863            continue
864        else:                   #no more rings beyond edge of detector
865            data['ellipses'].append([])
866            continue
867    if not data['rings']:
868        G2fil.G2Print ('no rings found; try lower Min ring I/Ib',mode='warn')
869        return []   
870       
871    rings = np.concatenate((data['rings']),axis=0)
872    if getRingsOnly:
873        return rings,HKL
874    [chisq,vals,sigList,covar] = FitDetector(rings,varyList,parmDict,True,True)
875    data['wavelength'] = parmDict['wave']
876    data['distance'] = parmDict['dist']
877    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
878    data['rotation'] = np.mod(parmDict['phi'],360.0)
879    data['tilt'] = parmDict['tilt']
880    data['DetDepth'] = parmDict['dep']
881    data['chisq'] = chisq
882    N = len(data['ellipses'])
883    data['ellipses'] = []           #clear away individual ellipse fits
884    for H in HKL[:N]:
885        ellipse = GetEllipse(H[3],data)
886        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))   
887    G2fil.G2Print ('calibration time = %.3f'%(time.time()-time0))
888    if G2frame:
889        G2plt.PlotImage(G2frame,newImage=True)       
890    return [vals,varyList,sigList,parmDict,covar]
891
892def ImageCalibrate(G2frame,data):
893    '''Called to perform an initial image calibration after points have been
894    selected for the inner ring.
895    '''
896    import ImageCalibrants as calFile
897    G2fil.G2Print ('Image calibration:')
898    time0 = time.time()
899    ring = data['ring']
900    pixelSize = data['pixelSize']
901    scalex = 1000./pixelSize[0]
902    scaley = 1000./pixelSize[1]
903    pixLimit = data['pixLimit']
904    cutoff = data['cutoff']
905    varyDict = data['varyList']
906    if varyDict['dist'] and varyDict['wave']:
907        G2fil.G2Print ('ERROR - you can not simultaneously calibrate distance and wavelength')
908        return False
909    if len(ring) < 5:
910        G2fil.G2Print ('ERROR - not enough inner ring points for ellipse')
911        return False
912       
913    #fit start points on inner ring
914    data['ellipses'] = []
915    data['rings'] = []
916    outE = FitEllipse(ring)
917    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
918    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi**2: %.3f, Np: %d'
919    if outE:
920        G2fil.G2Print (fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1]))
921        ellipse = outE
922    else:
923        return False
924       
925    #setup 360 points on that ring for "good" fit
926    data['ellipses'].append(ellipse[:]+('g',))
927    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
928    if Ring:
929        ellipse = FitEllipse(Ring)
930        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]    #do again
931        ellipse = FitEllipse(Ring)
932    else:
933        G2fil.G2Print ('1st ring not sufficiently complete to proceed',mode='warn')
934        return False
935    if debug:
936        G2fil.G2Print (fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],
937            ellipse[2][0],ellipse[2][1],0.,len(Ring)))     #cent,phi,radii
938    data['ellipses'].append(ellipse[:]+('r',))
939    data['rings'].append(np.array(Ring))
940    G2plt.PlotImage(G2frame,newImage=True)
941   
942#setup for calibration
943    data['rings'] = []
944    if not data['calibrant']:
945        G2fil.G2Print ('Warning: no calibration material selected')
946        return True
947   
948    skip = data['calibskip']
949    dmin = data['calibdmin']
950#generate reflection set
951    calibrant = calFile.Calibrants[data['calibrant']]
952    Bravais,SGs,Cells = calibrant[:3]
953    HKL = []
954    for bravais,sg,cell in zip(Bravais,SGs,Cells):
955        A = G2lat.cell2A(cell)
956        if sg:
957            SGData = G2spc.SpcGroup(sg)[1]
958            hkl = G2pwd.getHKLpeak(dmin,SGData,A,Inst=None,nodup=True)
959            HKL += list(hkl)
960        else:
961            hkl = G2lat.GenHBravais(dmin,bravais,A)
962            HKL += list(hkl)
963    HKL = G2lat.sortHKLd(HKL,True,False)[skip:]
964#set up 1st ring
965    elcent,phi,radii = ellipse              #from fit of 1st ring
966    dsp = HKL[0][3]
967    G2fil.G2Print ('1st ring: try %.4f'%(dsp))
968    if varyDict['dist']:
969        wave = data['wavelength']
970        tth = 2.0*asind(wave/(2.*dsp))
971    else:   #varyDict['wave']!
972        dist = data['distance']
973        tth = npatan2d(radii[0],dist)
974        data['wavelength'] = wave =  2.0*dsp*sind(tth/2.0)
975    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
976    ttth = nptand(tth)
977    ctth = npcosd(tth)
978#1st estimate of tilt; assume ellipse - don't know sign though
979    if varyDict['tilt']:
980        tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
981        if not tilt:
982            G2fil.G2Print ('WARNING - selected ring was fitted as a circle')
983            G2fil.G2Print (' - if detector was tilted we suggest you skip this ring - WARNING')
984    else:
985        tilt = data['tilt']
986#1st estimate of dist: sample to detector normal to plane
987    if varyDict['dist']:
988        data['distance'] = dist = radii[0]**2/(ttth*radii[1])
989    else:
990        dist = data['distance']
991    if varyDict['tilt']:
992#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
993        zdisp = radii[1]*ttth*tand(tilt)
994        zdism = radii[1]*ttth*tand(-tilt)
995#cone axis position; 2 choices. Which is right?     
996#NB: zdisp is || to major axis & phi is rotation of minor axis
997#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
998        centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
999        centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
1000#check get same ellipse parms either way
1001#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
1002        fail = True
1003        i2 = 1
1004        while fail:
1005            dsp = HKL[i2][3]
1006            G2fil.G2Print ('2nd ring: try %.4f'%(dsp))
1007            tth = 2.0*asind(wave/(2.*dsp))
1008            ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
1009            G2fil.G2Print (fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1]))
1010            Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
1011            parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
1012                'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}       
1013            varyList = [item for item in varyDict if varyDict[item]]
1014            if len(Ringp) > 10:
1015                chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)[0]
1016                tiltp = parmDict['tilt']
1017                phip = parmDict['phi']
1018                centp = [parmDict['det-X'],parmDict['det-Y']]
1019                fail = False
1020            else:
1021                chip = 1e6
1022            ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
1023            G2fil.G2Print (fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1]))
1024            Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
1025            if len(Ringm) > 10:
1026                parmDict['tilt'] *= -1
1027                chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)[0]
1028                tiltm = parmDict['tilt']
1029                phim = parmDict['phi']
1030                centm = [parmDict['det-X'],parmDict['det-Y']]
1031                fail = False
1032            else:
1033                chim = 1e6
1034            if fail:
1035                i2 += 1
1036        if chip < chim:
1037            data['tilt'] = tiltp
1038            data['center'] = centp
1039            data['rotation'] = phip
1040        else:
1041            data['tilt'] = tiltm
1042            data['center'] = centm
1043            data['rotation'] = phim
1044        data['ellipses'].append(ellipsep[:]+('b',))
1045        data['rings'].append(np.array(Ringp))
1046        data['ellipses'].append(ellipsem[:]+('r',))
1047        data['rings'].append(np.array(Ringm))
1048        G2plt.PlotImage(G2frame,newImage=True)
1049    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
1050        data['DetDepth'] /= data['distance']
1051    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
1052        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
1053    varyList = [item for item in varyDict if varyDict[item]]
1054    data['rings'] = []
1055    data['ellipses'] = []
1056    for i,H in enumerate(HKL):
1057        dsp = H[3]
1058        tth = 2.0*asind(wave/(2.*dsp))
1059        if tth+abs(data['tilt']) > 90.:
1060            G2fil.G2Print ('next line is a hyperbola - search stopped')
1061            break
1062        if debug:   print ('HKLD:'+str(H[:4])+'2-theta: %.4f'%(tth))
1063        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
1064        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
1065        if debug:   print (fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1]))
1066        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
1067        if Ring:
1068            data['rings'].append(np.array(Ring))
1069            rings = np.concatenate((data['rings']),axis=0)
1070            if i:
1071                chisq = FitDetector(rings,varyList,parmDict,False)[0]
1072                data['distance'] = parmDict['dist']
1073                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
1074                data['rotation'] = parmDict['phi']
1075                data['tilt'] = parmDict['tilt']
1076                data['DetDepth'] = parmDict['dep']
1077                data['chisq'] = chisq
1078                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
1079                if debug:   print (fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings)))
1080            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
1081#            G2plt.PlotImage(G2frame,newImage=True)
1082        else:
1083            if debug:   print ('insufficient number of points in this ellipse to fit')
1084#            break
1085    G2plt.PlotImage(G2frame,newImage=True)
1086    fullSize = len(G2frame.ImageZ)/scalex
1087    if 2*radii[1] < .9*fullSize:
1088        G2fil.G2Print ('Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib')
1089    N = len(data['ellipses'])
1090    if N > 2:
1091        FitDetector(rings,varyList,parmDict)[0]
1092        data['wavelength'] = parmDict['wave']
1093        data['distance'] = parmDict['dist']
1094        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
1095        data['rotation'] = parmDict['phi']
1096        data['tilt'] = parmDict['tilt']
1097        data['DetDepth'] = parmDict['dep']
1098    for H in HKL[:N]:
1099        ellipse = GetEllipse(H[3],data)
1100        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
1101    G2fil.G2Print ('calibration time = %.3f'%(time.time()-time0))
1102    G2plt.PlotImage(G2frame,newImage=True)       
1103    return True
1104   
1105def Make2ThetaAzimuthMap(data,iLim,jLim): #most expensive part of integration!
1106    'Needs a doc string'
1107    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
1108    pixelSize = data['pixelSize']
1109    scalex = pixelSize[0]/1000.
1110    scaley = pixelSize[1]/1000.
1111    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
1112    tax = np.asfarray(tax*scalex,dtype=np.float32).flatten()
1113    tay = np.asfarray(tay*scaley,dtype=np.float32).flatten()
1114    nI = iLim[1]-iLim[0]
1115    nJ = jLim[1]-jLim[0]
1116    TA = np.empty((4,nI,nJ))
1117    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
1118    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
1119    TA[3] = G2pwd.Polarization(data['PolaVal'][0],TA[0],TA[1]-90.)[0]
1120    return TA           #2-theta, azimuth & geom. corr. arrays
1121
1122def MakeMaskMap(data,masks,iLim,jLim,tamp):
1123    import polymask as pm
1124    pixelSize = data['pixelSize']
1125    scalex = pixelSize[0]/1000.
1126    scaley = pixelSize[1]/1000.
1127   
1128    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
1129    tax = np.asfarray(tax*scalex,dtype=np.float32).flatten()
1130    tay = np.asfarray(tay*scaley,dtype=np.float32).flatten()
1131    nI = iLim[1]-iLim[0]
1132    nJ = jLim[1]-jLim[0]
1133    #make position masks here
1134    frame = masks['Frames']
1135    tam = ma.make_mask_none((nI*nJ))
1136    if frame:
1137        tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
1138            tay,len(frame),frame,tamp)[:nI*nJ])^True)
1139    polygons = masks['Polygons']
1140    for polygon in polygons:
1141        if polygon:
1142            tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
1143                tay,len(polygon),polygon,tamp)[:nI*nJ]))
1144    for X,Y,rsq in masks['Points'].T:
1145        tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,rsq)))
1146    if tam.shape: 
1147        tam = np.reshape(tam,(nI,nJ))
1148    else:
1149        tam = ma.make_mask_none((nI,nJ))
1150    for xline in masks.get('Xlines',[]):    #a y pixel position
1151        if iLim[0] <= xline <= iLim[1]:
1152            tam[xline-iLim[0],:] = True
1153    for yline in masks.get('Ylines',[]):    #a x pixel position
1154        if jLim[0] <= yline <= jLim[1]:
1155            tam[:,yline-jLim[0]] = True           
1156    return tam           #position mask
1157
1158def Fill2ThetaAzimuthMap(masks,TA,tam,image):
1159    'Needs a doc string'
1160    Zlim = masks['Thresholds'][1]
1161    rings = masks['Rings']
1162    arcs = masks['Arcs']
1163    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2]),ma.getdata(TA[3])))    #azimuth, 2-theta, dist, pol
1164    tax,tay,tad,pol = np.dsplit(TA,4)    #azimuth, 2-theta, dist**2/d0**2, pol
1165    for tth,thick in rings:
1166        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
1167    for tth,azm,thick in arcs:
1168        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
1169        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
1170        tam = ma.mask_or(tam.flatten(),tamt*tama)
1171    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
1172    tabs = np.ones_like(taz)
1173    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
1174    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
1175    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
1176    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
1177    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
1178    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
1179    pol = ma.compressed(ma.array(pol.flatten(),mask=tam))   #polarization
1180    return tax,tay,taz/pol,tad,tabs
1181
1182def MakeUseTA(data,blkSize=128):
1183    Nx,Ny = data['size']
1184    nXBlks = (Nx-1)//blkSize+1
1185    nYBlks = (Ny-1)//blkSize+1
1186    useTA = []
1187    for iBlk in range(nYBlks):
1188        iBeg = iBlk*blkSize
1189        iFin = min(iBeg+blkSize,Ny)
1190        useTAj = []
1191        for jBlk in range(nXBlks):
1192            jBeg = jBlk*blkSize
1193            jFin = min(jBeg+blkSize,Nx)
1194            TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))          #2-theta & azimuth arrays & create position mask
1195            useTAj.append(TA)
1196        useTA.append(useTAj)
1197    return useTA
1198
1199def MakeUseMask(data,masks,blkSize=128):
1200    Masks = copy.deepcopy(masks)
1201    Masks['Points'] = np.array(Masks['Points']).T           #get spots as X,Y,R arrays
1202    if np.any(masks['Points']):
1203        Masks['Points'][2] = np.square(Masks['Points'][2]/2.)
1204    Nx,Ny = data['size']
1205    nXBlks = (Nx-1)//blkSize+1
1206    nYBlks = (Ny-1)//blkSize+1
1207    useMask = []
1208    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
1209    for iBlk in range(nYBlks):
1210        iBeg = iBlk*blkSize
1211        iFin = min(iBeg+blkSize,Ny)
1212        useMaskj = []
1213        for jBlk in range(nXBlks):
1214            jBeg = jBlk*blkSize
1215            jFin = min(jBeg+blkSize,Nx)
1216            mask = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp)          #2-theta & azimuth arrays & create position mask
1217            useMaskj.append(mask)
1218        useMask.append(useMaskj)
1219    return useMask
1220
1221def ImageIntegrate(image,data,masks,blkSize=128,returnN=False,useTA=None,useMask=None):
1222    'Integrate an image; called from OnIntegrateAll and OnIntegrate in G2imgGUI'    #for q, log(q) bins need data['binType']
1223    import histogram2d as h2d
1224    G2fil.G2Print ('Begin image integration; image range: %d %d'%(np.min(image),np.max(image)))
1225    CancelPressed = False
1226    LUtth = np.array(data['IOtth'])
1227    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
1228    numAzms = data['outAzimuths']
1229    numChans = (data['outChannels']//4)*4
1230    Dazm = (LRazm[1]-LRazm[0])/numAzms
1231    if '2-theta' in data.get('binType','2-theta'):
1232        lutth = LUtth               
1233    elif 'log(q)' in data['binType']:
1234        lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
1235    elif 'q' == data['binType'].lower():
1236        lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
1237    dtth = (lutth[1]-lutth[0])/numChans
1238    muT = data.get('SampleAbs',[0.0,''])[0]
1239    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
1240        data['DetDepth'] /= data['distance']
1241    if 'SASD' in data['type']:
1242        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
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    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
1248    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
1249    H2 = np.linspace(lutth[0],lutth[1],numChans+1)
1250    Nx,Ny = data['size']
1251    nXBlks = (Nx-1)//blkSize+1
1252    nYBlks = (Ny-1)//blkSize+1
1253    tbeg = time.time()
1254    times = [0,0,0,0,0]
1255    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
1256    for iBlk in range(nYBlks):
1257        iBeg = iBlk*blkSize
1258        iFin = min(iBeg+blkSize,Ny)
1259        for jBlk in range(nXBlks):
1260            jBeg = jBlk*blkSize
1261            jFin = min(jBeg+blkSize,Nx)
1262            # next is most expensive step!
1263
1264            t0 = time.time()
1265            if useTA:
1266                TA = useTA[iBlk][jBlk]
1267            else:
1268                TA = Make2ThetaAzimuthMap(data,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
1269            times[1] += time.time()-t0      #xy->th,azm
1270
1271            t0 = time.time()
1272            if useMask:
1273                tam = useMask[iBlk][jBlk]
1274            else:
1275                tam = MakeMaskMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp)
1276            times[0] += time.time()-t0      #apply masks
1277
1278            t0 = time.time()
1279            Block = image[iBeg:iFin,jBeg:jFin]          #apply image spotmask here
1280            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(Masks,TA,tam,Block)    #and apply masks
1281            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
1282            tax = np.where(tax < LRazm[0],tax+360.,tax)
1283            if data.get('SampleAbs',[0.0,''])[1]:
1284                if 'Cylind' in data['SampleShape']:
1285                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))      #adjust for additional thickness off sample normal
1286                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
1287                elif 'Fixed' in data['SampleShape']:    #assumes flat plate sample normal to beam
1288                    tabs = G2pwd.Absorb('Fixed',muT,tay)
1289            if 'log(q)' in data.get('binType',''):
1290                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
1291            elif 'q' == data.get('binType','').lower():
1292                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
1293            times[2] += time.time()-t0          #fill map
1294
1295            t0 = time.time()
1296            taz = np.array((taz*tad/tabs),dtype='float32')
1297            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
1298                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,
1299                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
1300            times[3] += time.time()-t0          #binning
1301    G2fil.G2Print('End integration loops')
1302
1303    t0 = time.time()
1304    #prepare masked arrays of bins with pixels for interpolation setup
1305    H2msk = [ma.array(H2[:-1],mask=np.logical_not(nst)) for nst in NST]
1306#    H0msk = [ma.array(np.divide(h0,nst),mask=np.logical_not(nst)) for nst,h0 in zip(NST,H0)]
1307    H0msk = [ma.array(h0/nst,mask=np.logical_not(nst)) for nst,h0 in zip(NST,H0)]
1308    #make linear interpolators; outside limits give NaN
1309    H0 = []
1310    problemEntries = []
1311    for i,(h0msk,h2msk) in enumerate(zip(H0msk,H2msk)):
1312        try:
1313            h0int = scint.interp1d(h2msk.compressed(),h0msk.compressed(),bounds_error=False)
1314            #do interpolation on all points - fills in the empty bins; leaves others the same
1315            h0 = h0int(H2[:-1])
1316            H0.append(h0)
1317        except ValueError:
1318            problemEntries.append(i)
1319            H0.append(np.zeros(numChans,order='F',dtype=np.float32))
1320    H0 = np.nan_to_num(np.array(H0))       
1321    if 'log(q)' in data.get('binType',''):
1322        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
1323    elif 'q' == data.get('binType','').lower():
1324        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
1325    if Dazm:       
1326        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
1327    else:
1328        H1 = LRazm
1329    if 'SASD' not in data['type']:
1330        H0 *= np.array(G2pwd.Polarization(data['PolaVal'][0],H2[:-1],0.)[0])
1331    H0 /= np.abs(npcosd(H2[:-1]-np.abs(data['det2theta'])))           #parallax correction
1332    if 'SASD' in data['type']:
1333        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
1334    if data['Oblique'][1]:
1335        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
1336    times[4] += time.time()-t0          #cleanup
1337
1338    G2fil.G2Print ('Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
1339        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4]))
1340    G2fil.G2Print ("Elapsed time:","%8.3fs"%(time.time()-tbeg))
1341    G2fil.G2Print ('Integration complete')
1342    if problemEntries:
1343        msg = ""
1344        for i in problemEntries:
1345            if msg: msg += ', '
1346            msg += "{:.2f}".format((H1[i+1]+H1[i])/2.)
1347        print('\nWarning, integrations have no pixels at these azimuth(s)\n\t',msg,'\n')
1348    if returnN:     #As requested by Steven Weigand
1349        return H0,H1,H2,NST,CancelPressed
1350    else:
1351        return H0,H1,H2,CancelPressed
1352   
1353def MakeStrStaRing(ring,Image,Controls):
1354    ellipse = GetEllipse(ring['Dset'],Controls)
1355    pixSize = Controls['pixelSize']
1356    scalex = 1000./pixSize[0]
1357    scaley = 1000./pixSize[1]
1358    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
1359    if len(Ring):
1360        ring['ImxyObs'] = copy.copy(Ring[:2])
1361        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
1362        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
1363        ring['ImtaObs'] = TA
1364        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
1365        Ring[0] = TA[0]
1366        Ring[1] = TA[1]
1367        return Ring,ring
1368    else:
1369        ring['ImxyObs'] = [[],[]]
1370        ring['ImtaObs'] = [[],[]]
1371        ring['ImtaCalc'] = [[],[]]
1372        return [],[]    #bad ring; no points found
1373   
1374def FitStrSta(Image,StrSta,Controls):
1375    'Needs a doc string'
1376   
1377    StaControls = copy.deepcopy(Controls)
1378    phi = StrSta['Sample phi']
1379    wave = Controls['wavelength']
1380    pixelSize = Controls['pixelSize']
1381    scalex = 1000./pixelSize[0]
1382    scaley = 1000./pixelSize[1]
1383    StaType = StrSta['Type']
1384    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1385
1386    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1387        dset = ring['Dset']
1388        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1389        if len(Ring):
1390            ring.update(R)
1391            p0 = ring['Emat']
1392            val,esd,covMat = FitStrain(Ring,p0,dset,wave,phi,StaType)
1393            ring['Emat'] = val
1394            ring['Esig'] = esd
1395            ellipse = FitEllipse(R['ImxyObs'].T)
1396            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1397            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1398            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1399            ringint /= np.mean(ringint)
1400            ring['Ivar'] = np.var(ringint)
1401            ring['covMat'] = covMat
1402            G2fil.G2Print ('Variance in normalized ring intensity: %.3f'%(ring['Ivar']))
1403    CalcStrSta(StrSta,Controls)
1404   
1405def IntStrSta(Image,StrSta,Controls):
1406    StaControls = copy.deepcopy(Controls)
1407    pixelSize = Controls['pixelSize']
1408    scalex = 1000./pixelSize[0]
1409    scaley = 1000./pixelSize[1]
1410    phi = StrSta['Sample phi']
1411    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1412    RingsAI = []
1413    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1414        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1415        if len(Ring):
1416            ellipse = FitEllipse(R['ImxyObs'].T)
1417            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image,5)
1418            XY = np.array(ringxy).T
1419            Th,Azm = GetTthAzm(XY[0],XY[1],Controls)
1420            pola = G2pwd.Polarization(Controls['PolaVal'][0],Th,Azm-90.)[0]
1421            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1422            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1423            ringint /= np.mean(ringint)
1424            ringint /= pola
1425            G2fil.G2Print (' %s %.3f %s %.3f %s %d'%('d-spacing',ring['Dcalc'],'sig(MRD):',np.sqrt(np.var(ringint)),'# points:',len(ringint)))
1426            RingsAI.append(np.array(list(zip(ringazm,ringint))).T)
1427    return RingsAI
1428   
1429def CalcStrSta(StrSta,Controls):
1430
1431    wave = Controls['wavelength']
1432    phi = StrSta['Sample phi']
1433    StaType = StrSta['Type']
1434    for ring in StrSta['d-zero']:
1435        Eij = ring['Emat']
1436        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
1437        th,azm = ring['ImtaObs']
1438        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
1439        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
1440        if StaType == 'True':
1441            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
1442        else:
1443            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
1444        dmin = np.min(ring['ImtaCalc'][0])
1445        dmax = np.max(ring['ImtaCalc'][0])
1446        if ring.get('fixDset',True):
1447            if abs(Eij[0]) < abs(Eij[2]):         #tension
1448                ring['Dcalc'] = dmin+(dmax-dmin)/4.
1449            else:                       #compression
1450                ring['Dcalc'] = dmin+3.*(dmax-dmin)/4.
1451        else:
1452            ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
1453
1454def calcFij(omg,phi,azm,th):
1455    '''    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
1456
1457    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
1458        90 when perp. to sample surface
1459    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
1460    :param azm: his chi = azimuth around incident beam
1461    :param th:  his theta = theta
1462    '''
1463    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
1464    b = -npcosd(azm)*npcosd(th)
1465    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
1466    d = a*npsind(phi)+b*npcosd(phi)
1467    e = a*npcosd(phi)-b*npsind(phi)
1468    Fij = np.array([
1469        [d**2,d*e,c*d],
1470        [d*e,e**2,c*e],
1471        [c*d,c*e,c**2]])
1472    return -Fij
1473
1474def FitStrain(rings,p0,dset,wave,phi,StaType):
1475    'Needs a doc string'
1476    def StrainPrint(ValSig,dset):
1477        print ('Strain tensor for Dset: %.6f'%(dset))
1478        ptlbls = 'names :'
1479        ptstr =  'values:'
1480        sigstr = 'esds  :'
1481        for name,fmt,value,sig in ValSig:
1482            ptlbls += "%s" % (name.rjust(12))
1483            ptstr += fmt % (value)
1484            if sig:
1485                sigstr += fmt % (sig)
1486            else:
1487                sigstr += 12*' '
1488        print (ptlbls)
1489        print (ptstr)
1490        print (sigstr)
1491       
1492    def strainCalc(p,xyd,dset,wave,phi,StaType):
1493        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1494        dspo,azm,dsp = xyd
1495        th = npasind(wave/(2.0*dspo))
1496        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1497        if StaType == 'True':
1498            dspc = dset*np.exp(V)
1499        else:
1500            dspc = dset*(V+1.)
1501        return dspo-dspc
1502       
1503    names = ['e11','e12','e22']
1504    fmt = ['%12.2f','%12.2f','%12.2f']
1505    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1506    vals = list(result[0])
1507    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1508    covM = result[1]
1509    covMatrix = covM*chisq
1510    sig = list(np.sqrt(chisq*np.diag(result[1])))
1511    ValSig = zip(names,fmt,vals,sig)
1512    StrainPrint(ValSig,dset)
1513    return vals,sig,covMatrix
1514
1515def FitImageSpots(Image,ImMax,ind,pixSize,nxy,spotSize=1.0):
1516   
1517    def calcMean(nxy,pixSize,img):
1518        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1519        gdx = ma.array((gdx-nxy//2)*pixSize[0]/1000.,mask=~ma.getmaskarray(ImBox))
1520        gdy = ma.array((gdy-nxy//2)*pixSize[1]/1000.,mask=~ma.getmaskarray(ImBox))
1521        posx = ma.sum(gdx)/ma.count(gdx)
1522        posy = ma.sum(gdy)/ma.count(gdy)
1523        return posx,posy
1524   
1525    def calcPeak(values,nxy,pixSize,img):
1526        back,mag,px,py,sig = values
1527        peak = np.zeros([nxy,nxy])+back
1528        nor = 1./(2.*np.pi*sig**2)
1529        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1530        gdx = (gdx-nxy//2)*pixSize[0]/1000.
1531        gdy = (gdy-nxy//2)*pixSize[1]/1000.
1532        arg = (gdx-px)**2+(gdy-py)**2       
1533        peak += mag*nor*np.exp(-arg/(2.*sig**2))
1534        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))
1535   
1536    def calc2Peak(values,nxy,pixSize,img):
1537        back,mag,px,py,sigx,sigy,rho = values
1538        peak = np.zeros([nxy,nxy])+back
1539        nor = 1./(2.*np.pi*sigx*sigy*np.sqrt(1.-rho**2))
1540        gdx,gdy = np.mgrid[0:nxy,0:nxy]
1541        gdx = (gdx-nxy//2)*pixSize[0]/1000.
1542        gdy = (gdy-nxy//2)*pixSize[1]/1000.
1543        argnor = -1./(2.*(1.-rho**2))
1544        arg = (gdx-px)**2/sigx**2+(gdy-py)**2/sigy**2-2.*rho*(gdx-px)*(gdy-py)/(sigx*sigy)       
1545        peak += mag*nor*np.exp(argnor*arg)
1546        return ma.compressed(img-peak)/np.sqrt(ma.compressed(img))       
1547   
1548    nxy2 = nxy//2
1549    ImBox = Image[ind[1]-nxy2:ind[1]+nxy2+1,ind[0]-nxy2:ind[0]+nxy2+1]
1550    back = np.min(ImBox)
1551    mag = np.sum(ImBox-back)
1552    vals = [back,mag,0.,0.,0.2,0.2,0.]
1553    ImBox = ma.array(ImBox,dtype=float,mask=ImBox>0.75*ImMax)
1554    px = (ind[0]+.5)*pixSize[0]/1000.
1555    py = (ind[1]+.5)*pixSize[1]/1000.
1556    if ma.any(ma.getmaskarray(ImBox)):
1557        vals = calcMean(nxy,pixSize,ImBox)
1558        posx,posy = [px+vals[0],py+vals[1]]
1559        return [posx,posy,spotSize]
1560    else:
1561        result = leastsq(calc2Peak,vals,args=(nxy,pixSize,ImBox),full_output=True)
1562        vals = result[0]
1563        ratio = vals[4]/vals[5]
1564        if 0.5 < ratio < 2.0 and vals[2] < 2. and vals[3] < 2.:
1565            posx,posy = [px+vals[2],py+vals[3]]
1566            return [posx,posy,min(6.*vals[4],spotSize)]
1567        else:
1568            return None
1569   
1570def AutoSpotMask(Image,Masks,Controls,numChans,dlg=None):
1571   
1572    frame = Masks['Frames']
1573    tam = ma.make_mask_none(Image.shape)
1574    if frame:
1575        tam = ma.mask_or(tam,MakeFrameMask(Controls,frame))
1576    LUtth = np.array(Controls['IOtth'])
1577    dtth = (LUtth[1]-LUtth[0])/numChans
1578    esdMul = Masks['SpotMask']['esdMul']
1579    prob = 100.*sc.erf(esdMul/np.sqrt(2.))
1580    print(' Spots greater than %.2f of band intensity are masked'%prob)
1581    mask = ma.make_mask_none(Image.shape)
1582    band = ma.array(Image,mask=ma.nomask)
1583    TA = Make2ThetaAzimuthMap(Controls,(0,Image.shape[0]),(0,Image.shape[1]))[0]    #2-theta array
1584    TThs = np.linspace(LUtth[0],LUtth[1],numChans,False)
1585    for it,TTh in enumerate(TThs):
1586        band.mask = ma.masked_outside(TA,TTh,TTh+dtth).mask+tam
1587        pcmax = np.percentile(band.compressed(),[prob,50.])
1588        mband = ma.masked_greater(band,pcmax[0])
1589        std = ma.std(mband)
1590        anom = ma.masked_greater((band-pcmax[1])/std,esdMul)
1591        mask ^= (anom.mask^band.mask)
1592        if not dlg is None:
1593            GoOn = dlg.Update(it,newmsg='Processed 2-theta rings = %d'%(it))
1594            if not GoOn[0]:
1595                break
1596    return mask
1597
1598def DoPolaCalib(ImageZ,imageData,arcTth):
1599    ''' Determine image polarization by successive integrations with & without preset arc mask.
1600    After initial search, does a set of five with offset azimuth to get mean(std) result.
1601    '''
1602    from scipy.optimize import minimize_scalar
1603    print(' Polarization fit for mask at %.1f deg 2-theta'%arcTth)
1604    data = copy.deepcopy(imageData)
1605    data['IOtth'] = [arcTth-2.,arcTth+2.]
1606    data['fullIntegrate'] = True
1607    data['LRazimuth'] = [0.,360.]
1608    data['outChannels'] = 200
1609    data['outAzimuths'] = 1
1610    Arc = [arcTth,[85.,95.],2.0]
1611    print(' Integration 2-theta test range %.1f - %.1f in 200 steps'%(data['IOtth'][0],data['IOtth'][1]))
1612    print(' Mask azimuth range: %.1f - %.1f'%(Arc[1][0],Arc[1][1]))
1613    Masks = {'Points':[],'Rings':[],'Arcs':[],'Polygons':[],'Frames':[],
1614        'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
1615    AMasks = {'Points':[],'Rings':[],'Arcs':[Arc,],'Polygons':[],'Frames':[],
1616        'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
1617   
1618    def func(p):
1619        p = min(1.0,max(p,0.0))
1620        data['PolaVal'][0] = p
1621        H0 = ImageIntegrate(ImageZ,data,Masks)[0]
1622        H0A = ImageIntegrate(ImageZ,data,AMasks)[0]
1623        M = np.sum(H0)-np.sum(H0A)
1624        print(' Polarization %.4f, fxn: %.1f'%(p,M))
1625        return M**2
1626    time0 = time.time()
1627    res = minimize_scalar(func,bracket=[1.,.999],tol=.0001)
1628    print(res)
1629    pola = min(1.0,max(res.x,.0))
1630   
1631    Pola = []
1632    for arc in [75,80,85,90,95]:
1633        Arc = [arcTth,[arc,arc+10.],2.0]
1634        AMasks = {'Points':[],'Rings':[],'Arcs':[Arc,],'Polygons':[],'Frames':[],
1635            'Thresholds':imageData['range'],'SpotMask':{'esdMul':3.,'spotMask':None}}
1636        res = minimize_scalar(func,bracket=[pola-.001,pola],tol=.0001)
1637        print(' Mask azimuth range: %.1f - %.1f'%(Arc[1][0],Arc[1][1]))
1638        print(' pola: %.5f'%res.x)
1639        Pola.append(res.x)
1640    Pola = np.array(Pola)
1641    mean = np.mean(Pola)
1642    std = int(10000*np.std(Pola))
1643    print(' Polarization: %.4f(%d)'%(mean,std))
1644    print(' time: %.2fs'%(time.time()-time0))
1645    imageData['PolaVal'][0] = mean
1646   
Note: See TracBrowser for help on using the repository browser.