source: trunk/GSASIIimage.py @ 2125

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

a bit of simplification to G2sfact for HKLF5 data
Add sample shape (Cylinder & Fixed flat plate) to ImageGUI & integration

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