source: trunk/GSASIIimage.py @ 2421

Last change on this file since 2421 was 2421, checked in by vondreele, 5 years ago

name sequential results by type of sequence (peak fit, SASD, etc.) so multiple tables are kept in project
add reverse & copyforward to seq. strain fitting
trap cancel from various sequential fitting

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