source: trunk/GSASIIimage.py @ 2305

Last change on this file since 2305 was 2305, checked in by vondreele, 6 years ago

in G2image.ImageIntegrate? force taz to be float-32 (it was float-64 & fortran histogram2d wanted real*4.

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