source: trunk/GSASIIimage.py @ 2183

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

fix General - atomic wt. now updated on isotope change
add new Stacking fault option to make sequence of DIFFaX runs varying one variable - under development
fix problem in G2image & G2imageGUI where change in FlatBkg?, Background or Dark didn't change thresholds or image plot limits - these are now set to new image limits.
Some of these used tricky changes to sizer values

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 42.1 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2016-03-24 14:58:40 +0000 (Thu, 24 Mar 2016) $
5# $Author: vondreele $
6# $Revision: 2183 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 2183 2016-03-24 14:58:40Z 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: 2183 $")
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    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
417    dxy = peneCorr(tth,data['DetDepth'],tilt,npatan2d(dy,dx))
418    tth = npatan2d(np.sqrt(dx**2+dy**2-Z**2),dist-Z+dxy) 
419    azm = (npatan2d(dy,dx)+data['azmthOff']+720.)%360.
420    return tth,azm,G
421
422def GetDsp(x,y,data):
423    'Give d-spacing value for detector x,y position; calibration info in data'
424    return GetTthAzmDsp(x,y,data)[3]
425       
426def GetAzm(x,y,data):
427    'Give azimuth value for detector x,y position; calibration info in data'
428    return GetTthAzmDsp(x,y,data)[1]
429   
430def meanAzm(a,b):
431    AZM = lambda a,b: npacosd(0.5*(npsind(2.*b)-npsind(2.*a))/(np.pi*(b-a)/180.))/2.
432    azm = AZM(a,b)
433#    quad = int((a+b)/180.)
434#    if quad == 1:
435#        azm = 180.-azm
436#    elif quad == 2:
437#        azm += 180.
438#    elif quad == 3:
439#        azm = 360-azm
440    return azm     
441       
442def ImageCompress(image,scale):
443    ''' Reduces size of image by selecting every n'th point
444    param: image array: original image
445    param: scale int: intervsl between selected points
446    returns: array: reduced size image
447    '''
448    if scale == 1:
449        return image
450    else:
451        return image[::scale,::scale]
452       
453def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
454    'Needs a doc string'
455    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
456    curr = np.array([dist,x,y])
457    return abs(avg-curr)/avg < .02
458
459def EdgeFinder(image,data):
460    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
461    Not currently used but might be useful in future?
462    '''
463    import numpy.ma as ma
464    Nx,Ny = data['size']
465    pixelSize = data['pixelSize']
466    edgemin = data['edgemin']
467    scalex = pixelSize[0]/1000.
468    scaley = pixelSize[1]/1000.   
469    tay,tax = np.mgrid[0:Nx,0:Ny]
470    tax = np.asfarray(tax*scalex,dtype=np.float32)
471    tay = np.asfarray(tay*scaley,dtype=np.float32)
472    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
473    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
474    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
475    return zip(tax,tay)
476   
477def MakeFrameMask(data,frame):
478    pixelSize = data['pixelSize']
479    scalex = pixelSize[0]/1000.
480    scaley = pixelSize[1]/1000.
481    blkSize = 512
482    Nx,Ny = data['size']
483    nXBlks = (Nx-1)/blkSize+1
484    nYBlks = (Ny-1)/blkSize+1
485    tam = ma.make_mask_none(data['size'])
486    for iBlk in range(nXBlks):
487        iBeg = iBlk*blkSize
488        iFin = min(iBeg+blkSize,Nx)
489        for jBlk in range(nYBlks):
490            jBeg = jBlk*blkSize
491            jFin = min(jBeg+blkSize,Ny)               
492            nI = iFin-iBeg
493            nJ = jFin-jBeg
494            tax,tay = np.mgrid[iBeg+0.5:iFin+.5,jBeg+.5:jFin+.5]         #bin centers not corners
495            tax = np.asfarray(tax*scalex,dtype=np.float32)
496            tay = np.asfarray(tay*scaley,dtype=np.float32)
497            tamp = ma.make_mask_none((1024*1024))
498            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
499                tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
500            if tamp.shape:
501                tamp = np.reshape(tamp[:nI*nJ],(nI,nJ))
502                tam[iBeg:iFin,jBeg:jFin] = ma.mask_or(tamp[0:nI,0:nJ],tam[iBeg:iFin,jBeg:jFin])
503            else:
504                tam[iBeg:iFin,jBeg:jFin] = True
505    return tam.T
506   
507def ImageRecalibrate(G2frame,data,masks):
508    '''Called to repeat the calibration on an image, usually called after
509    calibration is done initially to improve the fit.
510    '''
511    import ImageCalibrants as calFile
512    print 'Image recalibration:'
513    time0 = time.time()
514    pixelSize = data['pixelSize']
515    scalex = 1000./pixelSize[0]
516    scaley = 1000./pixelSize[1]
517    pixLimit = data['pixLimit']
518    cutoff = data['cutoff']
519    data['rings'] = []
520    data['ellipses'] = []
521    if not data['calibrant']:
522        print 'no calibration material selected'
523        return True   
524    skip = data['calibskip']
525    dmin = data['calibdmin']
526    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
527    HKL = []
528    for bravais,sg,cell in zip(Bravais,SGs,Cells):
529        A = G2lat.cell2A(cell)
530        if sg:
531            SGData = G2spc.SpcGroup(sg)[1]
532            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
533            HKL += hkl
534        else:
535            hkl = G2lat.GenHBravais(dmin,bravais,A)
536            HKL += hkl
537    HKL = G2lat.sortHKLd(HKL,True,False)
538    varyList = [item for item in data['varyList'] if data['varyList'][item]]
539    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
540        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
541    Found = False
542    wave = data['wavelength']
543    frame = masks['Frames']
544    tam = ma.make_mask_none(G2frame.ImageZ.shape)
545    if frame:
546        tam = ma.mask_or(tam,MakeFrameMask(data,frame))
547    for iH,H in enumerate(HKL):
548        if debug:   print H
549        dsp = H[3]
550        tth = 2.0*asind(wave/(2.*dsp))
551        if tth+abs(data['tilt']) > 90.:
552            print 'next line is a hyperbola - search stopped'
553            break
554        ellipse = GetEllipse(dsp,data)
555        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(G2frame.ImageZ,mask=tam))
556        if Ring:
557            if iH >= skip:
558                data['rings'].append(np.array(Ring))
559            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
560            Found = True
561        elif not Found:         #skipping inner rings, keep looking until ring found
562            continue
563        else:                   #no more rings beyond edge of detector
564            data['ellipses'].append([])
565            continue
566    if not data['rings']:
567        print 'no rings found; try lower Min ring I/Ib'
568        return True   
569       
570    rings = np.concatenate((data['rings']),axis=0)
571    chisq = FitDetector(rings,varyList,parmDict)
572    data['wavelength'] = parmDict['wave']
573    data['distance'] = parmDict['dist']
574    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
575    data['rotation'] = np.mod(parmDict['phi'],360.0)
576    data['tilt'] = parmDict['tilt']
577    data['DetDepth'] = parmDict['dep']
578    data['chisq'] = chisq
579    N = len(data['ellipses'])
580    data['ellipses'] = []           #clear away individual ellipse fits
581    for H in HKL[:N]:
582        ellipse = GetEllipse(H[3],data)
583        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))   
584    print 'calibration time = ',time.time()-time0
585    G2plt.PlotImage(G2frame,newImage=True)       
586    return True
587           
588def ImageCalibrate(G2frame,data):
589    '''Called to perform an initial image calibration after points have been
590    selected for the inner ring.
591    '''
592    import copy
593    import ImageCalibrants as calFile
594    print 'Image calibration:'
595    time0 = time.time()
596    ring = data['ring']
597    pixelSize = data['pixelSize']
598    scalex = 1000./pixelSize[0]
599    scaley = 1000./pixelSize[1]
600    pixLimit = data['pixLimit']
601    cutoff = data['cutoff']
602    varyDict = data['varyList']
603    if varyDict['dist'] and varyDict['wave']:
604        print 'ERROR - you can not simultaneously calibrate distance and wavelength'
605        return False
606    if len(ring) < 5:
607        print 'ERROR - not enough inner ring points for ellipse'
608        return False
609       
610    #fit start points on inner ring
611    data['ellipses'] = []
612    data['rings'] = []
613    outE = FitEllipse(ring)
614    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
615    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi**2: %.3f, Np: %d'
616    if outE:
617        print fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1])
618        ellipse = outE
619    else:
620        return False
621       
622    #setup 360 points on that ring for "good" fit
623    data['ellipses'].append(ellipse[:]+('g',))
624    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)
625    if Ring:
626        ellipse = FitEllipse(Ring)
627        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)    #do again
628        ellipse = FitEllipse(Ring)
629    else:
630        print '1st ring not sufficiently complete to proceed'
631        return False
632    if debug:
633        print fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],
634            ellipse[2][0],ellipse[2][1],0.,len(Ring))     #cent,phi,radii
635    data['ellipses'].append(ellipse[:]+('r',))
636    data['rings'].append(np.array(Ring))
637    G2plt.PlotImage(G2frame,newImage=True)
638   
639#setup for calibration
640    data['rings'] = []
641    if not data['calibrant']:
642        print 'no calibration material selected'
643        return True
644   
645    skip = data['calibskip']
646    dmin = data['calibdmin']
647#generate reflection set
648    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
649    HKL = []
650    for bravais,sg,cell in zip(Bravais,SGs,Cells):
651        A = G2lat.cell2A(cell)
652        if sg:
653            SGData = G2spc.SpcGroup(sg)[1]
654            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
655            HKL += hkl
656        else:
657            hkl = G2lat.GenHBravais(dmin,bravais,A)
658            HKL += hkl
659    HKL = G2lat.sortHKLd(HKL,True,False)[skip:]
660#set up 1st ring
661    elcent,phi,radii = ellipse              #from fit of 1st ring
662    dsp = HKL[0][3]
663    print '1st ring: try %.4f'%(dsp)
664    if varyDict['dist']:
665        wave = data['wavelength']
666        tth = 2.0*asind(wave/(2.*dsp))
667    else:   #varyDict['wave']!
668        dist = data['distance']
669        tth = npatan2d(radii[0],dist)
670        data['wavelength'] = wave =  2.0*dsp*sind(tth/2.0)
671    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,G2frame.ImageZ)
672    ttth = nptand(tth)
673    stth = npsind(tth)
674    ctth = npcosd(tth)
675#1st estimate of tilt; assume ellipse - don't know sign though
676    if varyDict['tilt']:
677        tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
678        if not tilt:
679            print 'WARNING - selected ring was fitted as a circle'
680            print ' - if detector was tilted we suggest you skip this ring - WARNING'
681    else:
682        tilt = data['tilt']
683#1st estimate of dist: sample to detector normal to plane
684    if varyDict['dist']:
685        data['distance'] = dist = radii[0]**2/(ttth*radii[1])
686    else:
687        dist = data['distance']
688    if varyDict['tilt']:
689#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
690        zdisp = radii[1]*ttth*tand(tilt)
691        zdism = radii[1]*ttth*tand(-tilt)
692#cone axis position; 2 choices. Which is right?     
693#NB: zdisp is || to major axis & phi is rotation of minor axis
694#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
695        centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
696        centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
697#check get same ellipse parms either way
698#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
699        fail = True
700        i2 = 1
701        while fail:
702            dsp = HKL[i2][3]
703            print '2nd ring: try %.4f'%(dsp)
704            tth = 2.0*asind(wave/(2.*dsp))
705            ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
706            print fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1])
707            Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,G2frame.ImageZ)
708            parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
709                'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}       
710            varyList = [item for item in varyDict if varyDict[item]]
711            if len(Ringp) > 10:
712                chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)
713                tiltp = parmDict['tilt']
714                phip = parmDict['phi']
715                centp = [parmDict['det-X'],parmDict['det-Y']]
716                fail = False
717            else:
718                chip = 1e6
719            ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
720            print fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1])
721            Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,G2frame.ImageZ)
722            if len(Ringm) > 10:
723                parmDict['tilt'] *= -1
724                chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)
725                tiltm = parmDict['tilt']
726                phim = parmDict['phi']
727                centm = [parmDict['det-X'],parmDict['det-Y']]
728                fail = False
729            else:
730                chim = 1e6
731            if fail:
732                i2 += 1
733        if chip < chim:
734            data['tilt'] = tiltp
735            data['center'] = centp
736            data['rotation'] = phip
737        else:
738            data['tilt'] = tiltm
739            data['center'] = centm
740            data['rotation'] = phim
741        data['ellipses'].append(ellipsep[:]+('b',))
742        data['rings'].append(np.array(Ringp))
743        data['ellipses'].append(ellipsem[:]+('r',))
744        data['rings'].append(np.array(Ringm))
745        G2plt.PlotImage(G2frame,newImage=True)
746    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
747        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
748    varyList = [item for item in varyDict if varyDict[item]]
749    data['rings'] = []
750    data['ellipses'] = []
751    for i,H in enumerate(HKL):
752        dsp = H[3]
753        tth = 2.0*asind(wave/(2.*dsp))
754        if tth+abs(data['tilt']) > 90.:
755            print 'next line is a hyperbola - search stopped'
756            break
757        if debug:   print 'HKLD:',H[:4],'2-theta: %.4f'%(tth)
758        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
759        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
760        if debug:   print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1])
761        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)
762        if Ring:
763            data['rings'].append(np.array(Ring))
764            rings = np.concatenate((data['rings']),axis=0)
765            if i:
766                chisq = FitDetector(rings,varyList,parmDict,False)
767                data['distance'] = parmDict['dist']
768                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
769                data['rotation'] = parmDict['phi']
770                data['tilt'] = parmDict['tilt']
771                data['DetDepth'] = parmDict['dep']
772                data['chisq'] = chisq
773                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
774                if debug:   print fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings))
775            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
776#            G2plt.PlotImage(G2frame,newImage=True)
777        else:
778            if debug:   print 'insufficient number of points in this ellipse to fit'
779#            break
780    G2plt.PlotImage(G2frame,newImage=True)
781    fullSize = len(G2frame.ImageZ)/scalex
782    if 2*radii[1] < .9*fullSize:
783        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
784    N = len(data['ellipses'])
785    if N > 2:
786        FitDetector(rings,varyList,parmDict)
787        data['wavelength'] = parmDict['wave']
788        data['distance'] = parmDict['dist']
789        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
790        data['rotation'] = parmDict['phi']
791        data['tilt'] = parmDict['tilt']
792        data['DetDepth'] = parmDict['dep']
793    for H in HKL[:N]:
794        ellipse = GetEllipse(H[3],data)
795        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
796    print 'calibration time = ',time.time()-time0
797    G2plt.PlotImage(G2frame,newImage=True)       
798    return True
799   
800def Make2ThetaAzimuthMap(data,masks,iLim,jLim,times): #most expensive part of integration!
801    'Needs a doc string'
802    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
803    pixelSize = data['pixelSize']
804    scalex = pixelSize[0]/1000.
805    scaley = pixelSize[1]/1000.
806   
807    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
808    tax = np.asfarray(tax*scalex,dtype=np.float32)
809    tay = np.asfarray(tay*scaley,dtype=np.float32)
810    nI = iLim[1]-iLim[0]
811    nJ = jLim[1]-jLim[0]
812    t0 = time.time()
813    #make position masks here
814    frame = masks['Frames']
815    tam = ma.make_mask_none((nI,nJ))
816    if frame:
817        tamp = ma.make_mask_none((1024*1024))
818        tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
819            tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
820        tam = ma.mask_or(tam.flatten(),tamp)
821    polygons = masks['Polygons']
822    for polygon in polygons:
823        if polygon:
824            tamp = ma.make_mask_none((1024*1024))
825            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
826                tay.flatten(),len(polygon),polygon,tamp)[:nI*nJ])
827            tam = ma.mask_or(tam.flatten(),tamp)
828    if tam.shape: tam = np.reshape(tam,(nI,nJ))
829    spots = masks['Points']
830    for X,Y,diam in spots:
831        tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2))
832        tam = ma.mask_or(tam,tamp)
833    times[0] += time.time()-t0
834    t0 = time.time()
835    TA = np.array(GetTthAzmG(tax,tay,data))     #includes geom. corr. as dist**2/d0**2 - most expensive step
836    times[1] += time.time()-t0
837    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
838    return np.array(TA),tam           #2-theta, azimuth & geom. corr. arrays & position mask
839
840def Fill2ThetaAzimuthMap(masks,TA,tam,image):
841    'Needs a doc string'
842    Zlim = masks['Thresholds'][1]
843    rings = masks['Rings']
844    arcs = masks['Arcs']
845    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2])))    #azimuth, 2-theta, dist
846    tax,tay,tad = np.dsplit(TA,3)    #azimuth, 2-theta, dist**2/d0**2
847    for tth,thick in rings:
848        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
849    for tth,azm,thick in arcs:
850        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
851        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
852        tam = ma.mask_or(tam.flatten(),tamt*tama)
853    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
854    tabs = np.ones_like(taz)
855    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
856    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
857    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
858    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
859    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
860    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
861    return tax,tay,taz,tad,tabs
862   
863def ImageIntegrate(image,data,masks,blkSize=128,dlg=None,returnN=False):
864    'Needs a doc string'    #for q, log(q) bins need data['binType']
865    import histogram2d as h2d
866    print 'Begin image integration'
867    LUtth = np.array(data['IOtth'])
868    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
869    numAzms = data['outAzimuths']
870    numChans = data['outChannels']
871    azmOff = data['azmthOff']
872    Dazm = (LRazm[1]-LRazm[0])/numAzms
873    if '2-theta' in data.get('binType','2-theta'):
874        lutth = LUtth               
875    elif 'Q' == data['binType']:
876        lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
877    elif 'log(q)' in data['binType']:
878        lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
879    dtth = (lutth[1]-lutth[0])/numChans
880    muT = data.get('SampleAbs',[0.0,''])[0]
881    if 'SASD' in data['type']:
882        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
883    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
884    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
885    imageN = len(image)
886    Nx,Ny = data['size']
887    nXBlks = (Nx-1)/blkSize+1
888    nYBlks = (Ny-1)/blkSize+1
889    Nup = nXBlks*nYBlks*3+3
890    tbeg = time.time()
891    Nup = 0
892    if dlg:
893        dlg.Update(Nup)
894    times = [0,0,0,0,0]
895    for iBlk in range(nYBlks):
896        iBeg = iBlk*blkSize
897        iFin = min(iBeg+blkSize,Ny)
898        for jBlk in range(nXBlks):
899            jBeg = jBlk*blkSize
900            jFin = min(jBeg+blkSize,Nx)
901            # next is most expensive step!
902            TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin),times)           #2-theta & azimuth arrays & create position mask
903            Nup += 1
904            if dlg:
905                dlg.Update(Nup)
906            Block = image[iBeg:iFin,jBeg:jFin]
907            t0 = time.time()
908            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
909            times[2] += time.time()-t0
910            Nup += 1
911            if dlg:
912                dlg.Update(Nup)
913            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
914            tax = np.where(tax < LRazm[0],tax+360.,tax)
915            if data.get('SampleAbs',[0.0,''])[1]:
916                if 'Cylind' in data['SampleShape']:
917                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))      #adjust for additional thickness off sample normal
918                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
919                elif 'Fixed' in data['SampleShape']:    #assumes flat plate sample normal to beam
920                    tabs = G2pwd.Absorb('Fixed',muT,tay)
921            if 'log(q)' in data.get('binType',''):
922                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
923            elif 'Q' == data.get('binType',''):
924                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
925            t0 = time.time()
926            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
927                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz*tad/tabs,
928                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
929            times[3] += time.time()-t0
930            Nup += 1
931            if dlg:
932                dlg.Update(Nup)
933    t0 = time.time()
934    NST = np.array(NST,dtype=np.float)
935    H0 = np.divide(H0,NST)
936    H0 = np.nan_to_num(H0)
937    H2 = np.array([tth for tth in np.linspace(lutth[0],lutth[1],numChans+1)])
938    if 'log(q)' in data.get('binType',''):
939        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
940    elif 'Q' == data.get('binType',''):
941        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
942    if Dazm:       
943        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
944    else:
945        H1 = LRazm
946    H0 /= npcosd(H2[:-1])           #**2? I don't think so, **1 is right for powders
947    if 'SASD' in data['type']:
948        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
949    if data['Oblique'][1]:
950        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
951    if 'SASD' in data['type'] and data['PolaVal'][1]:
952        #NB: in G2pwd.Polarization azm is defined from plane of polarization, not image x axis!
953        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm-90.)[0] for azm in (H1[:-1]+np.diff(H1)/2.)])
954    Nup += 1
955    if dlg:
956        dlg.Update(Nup)
957    times[4] += time.time()-t0
958    print 'Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
959        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4])
960    print "Elapsed time:","%8.3fs"%(time.time()-tbeg)
961    print 'Integration complete'
962    if returnN:     #As requested by Steven Weigand
963        return H0,H1,H2,NST
964    else:
965        return H0,H1,H2
966   
967def MakeStrStaRing(ring,Image,Controls):
968    ellipse = GetEllipse(ring['Dset'],Controls)
969    pixSize = Controls['pixelSize']
970    scalex = 1000./pixSize[0]
971    scaley = 1000./pixSize[1]
972    Ring = np.array(makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)).T   #returns x,y,dsp for each point in ring
973    if len(Ring):
974        ring['ImxyObs'] = copy.copy(Ring[:2])
975        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
976        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
977        ring['ImtaObs'] = TA
978        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
979        Ring[0] = TA[0]
980        Ring[1] = TA[1]
981        return Ring,ring
982    else:
983        ring['ImxyObs'] = [[],[]]
984        ring['ImtaObs'] = [[],[]]
985        ring['ImtaCalc'] = [[],[]]
986        return [],[]    #bad ring; no points found
987   
988def FitStrSta(Image,StrSta,Controls):
989    'Needs a doc string'
990   
991    StaControls = copy.deepcopy(Controls)
992    phi = StrSta['Sample phi']
993    wave = Controls['wavelength']
994    StaType = StrSta['Type']
995    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
996
997    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
998        dset = ring['Dset']
999        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1000        if len(Ring):
1001            ring.update(R)
1002            p0 = ring['Emat']
1003            val,esd = FitStrain(Ring,p0,dset,wave,phi,StaType)
1004            ring['Emat'] = val
1005            ring['Esig'] = esd
1006    CalcStrSta(StrSta,Controls)
1007   
1008def CalcStrSta(StrSta,Controls):
1009
1010    wave = Controls['wavelength']
1011    phi = StrSta['Sample phi']
1012    StaType = StrSta['Type']
1013    for ring in StrSta['d-zero']:
1014        Eij = ring['Emat']
1015        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
1016        th,azm = ring['ImtaObs']
1017        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
1018        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
1019        if StaType == 'True':
1020            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
1021        else:
1022            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
1023        ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
1024
1025def calcFij(omg,phi,azm,th):
1026    '''Does something...
1027
1028    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
1029
1030    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
1031        90 when perp. to sample surface
1032    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
1033    :param azm: his chi = azimuth around incident beam
1034    :param th:  his theta = theta
1035    '''
1036    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
1037    b = -npcosd(azm)*npcosd(th)
1038    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
1039    d = a*npsind(phi)+b*npcosd(phi)
1040    e = a*npcosd(phi)-b*npsind(phi)
1041    Fij = np.array([
1042        [d**2,d*e,c*d],
1043        [d*e,e**2,c*e],
1044        [c*d,c*e,c**2]])
1045    return -Fij
1046
1047def FitStrain(rings,p0,dset,wave,phi,StaType):
1048    'Needs a doc string'
1049    def StrainPrint(ValSig,dset):
1050        print 'Strain tensor for Dset: %.6f'%(dset)
1051        ptlbls = 'names :'
1052        ptstr =  'values:'
1053        sigstr = 'esds  :'
1054        for name,fmt,value,sig in ValSig:
1055            ptlbls += "%s" % (name.rjust(12))
1056            ptstr += fmt % (value)
1057            if sig:
1058                sigstr += fmt % (sig)
1059            else:
1060                sigstr += 12*' '
1061        print ptlbls
1062        print ptstr
1063        print sigstr
1064       
1065    def strainCalc(p,xyd,dset,wave,phi,StaType):
1066        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1067        dspo,azm,dsp = xyd
1068        th = npasind(wave/(2.0*dspo))
1069        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1070        if StaType == 'True':
1071            dspc = dset*np.exp(V)
1072        else:
1073            dspc = dset*(V+1.)
1074        return dspo-dspc
1075       
1076    names = ['e11','e12','e22']
1077    fmt = ['%12.2f','%12.2f','%12.2f']
1078    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1079    vals = list(result[0])
1080    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1081    sig = list(np.sqrt(chisq*np.diag(result[1])))
1082    ValSig = zip(names,fmt,vals,sig)
1083    StrainPrint(ValSig,dset)
1084    return vals,sig
1085   
1086       
Note: See TracBrowser for help on using the repository browser.