source: trunk/GSASIIimage.py @ 4316

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

slight streamline of AutoSpotMask2 - a bit faster

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