source: trunk/GSASIIimage.py @ 4569

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

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