source: trunk/GSASIIimage.py @ 4325

Last change on this file since 4325 was 4325, checked in by vondreele, 20 months ago

small RMCProfile fixes
New auto spot mask version - works

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