source: trunk/GSASIIimage.py @ 2468

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

add dzero calc from compression 2D images
plot magnetic structures; rework a number of draw atom routines to accomodate magnetic structures
fix EXP phase importer for nonmagnetic atoms in magnetic structure

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