source: trunk/GSASIIimage.py @ 3147

Last change on this file since 3147 was 3147, checked in by vondreele, 4 years ago

for image integration use linear interpolation to fill in bins with no pixels based on nearest filled bins.

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