source: trunk/GSASIIimage.py

Last change on this file was 5256, checked in by vondreele, 6 months ago

comments in G2image
a fix to REFD plotting
new help entries for SASD & REFD

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