source: trunk/GSASIIimage.py @ 1694

Last change on this file since 1694 was 1694, checked in by vondreele, 7 years ago

some missing stuff in PWDR[0] dictionary
add new fxn to G2img - meanAZm - calc effective azimuth for a azm range for correct polarization calc.
fix some FlexGridSizers? in G2phsGUI
fix '+'/'=' key response for Anaconda vs Enthought versions; one sees '=' the other sees '+'!
fix L/R shifting in waterfall plots.

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