source: trunk/GSASIIimage.py @ 2045

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

trap no rings found on Recalibrate, was cause of crash.

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