source: trunk/GSASIIimage.py @ 4415

Last change on this file since 4415 was 4415, checked in by vondreele, 18 months ago

change export_PDB to handle fullrmc pdb files
implement Wenqian's mod for the geometric image correction
TransformPhase? & FillUnitCell? now have option to not force atoms to new unit cell (default=True)
refactor FillAtomLookup?
change FindAllNeighbors? to optionally give short output
remove result of double click on Uiso column heading
many additions & changes for fulrmc GUI, etc.
comment out obsolete MPLsubplots from G2plot - no longer needed (pre mpl 2.2)
modify PDB importer to handle fullrmc pdb files.

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