source: trunk/GSASIIimage.py @ 3345

Last change on this file since 3345 was 3345, checked in by toby, 4 years ago

rework ValidatedTextCtrl? so that typeHint overrides the initial variable type; fix display of invalid numbers; use nDig to drive typeHint to float; postpone import of pyspg & polymask when not needed

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