source: trunk/GSASIIimage.py @ 4585

Last change on this file since 4585 was 4585, checked in by vondreele, 3 years ago

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