source: trunk/GSASIIimage.py @ 2694

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

after a repeat of a sequential fit, clear away old plot. Removes possibly bad parameter references
make note in G2grid of possible problem of pseudo vars storage in Controls
set new patch value for penetration corr to 0.5; above that convert to new value
fix peak picking problem when pattern has excluded regions
remove RaisePageNoRefresh? from Seq result plotting - irrelevant if plot is deleted

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 46.0 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2017-02-07 20:50:13 +0000 (Tue, 07 Feb 2017) $
5# $Author: vondreele $
6# $Revision: 2694 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 2694 2017-02-07 20:50:13Z vondreele $
9########### SVN repository information ###################
10'''
11*GSASIIimage: Image calc module*
12================================
13
14Ellipse fitting & image integration
15
16'''
17
18import math
19import time
20import numpy as np
21import numpy.linalg as nl
22import numpy.ma as ma
23import polymask as pm
24from scipy.optimize import leastsq
25import copy
26import GSASIIpath
27GSASIIpath.SetVersionNumber("$Revision: 2694 $")
28import GSASIIplot as G2plt
29import GSASIIlattice as G2lat
30import GSASIIpwd as G2pwd
31import GSASIIspc as G2spc
32
33# trig functions in degrees
34sind = lambda x: math.sin(x*math.pi/180.)
35asind = lambda x: 180.*math.asin(x)/math.pi
36tand = lambda x: math.tan(x*math.pi/180.)
37atand = lambda x: 180.*math.atan(x)/math.pi
38atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
39cosd = lambda x: math.cos(x*math.pi/180.)
40acosd = lambda x: 180.*math.acos(x)/math.pi
41rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
42#numpy versions
43npsind = lambda x: np.sin(x*np.pi/180.)
44npasind = lambda x: 180.*np.arcsin(x)/np.pi
45npcosd = lambda x: np.cos(x*np.pi/180.)
46npacosd = lambda x: 180.*np.arccos(x)/np.pi
47nptand = lambda x: np.tan(x*np.pi/180.)
48npatand = lambda x: 180.*np.arctan(x)/np.pi
49npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
50debug = False
51   
52def pointInPolygon(pXY,xy):
53    'Needs a doc string'
54    #pXY - assumed closed 1st & last points are duplicates
55    Inside = False
56    N = len(pXY)
57    p1x,p1y = pXY[0]
58    for i in range(N+1):
59        p2x,p2y = pXY[i%N]
60        if (max(p1y,p2y) >= xy[1] > min(p1y,p2y)) and (xy[0] <= max(p1x,p2x)):
61            if p1y != p2y:
62                xinters = (xy[1]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
63            if p1x == p2x or xy[0] <= xinters:
64                Inside = not Inside
65        p1x,p1y = p2x,p2y
66    return Inside
67   
68def peneCorr(tth,dep,dist,tilt=0.,azm=0.):
69    'Needs a doc string'
70#    return dep*(1.-npcosd(abs(tilt*npsind(azm))-tth*npcosd(azm)))  #something wrong here
71    return dep*(1.-npcosd(tth))*dist**2/1000.         #best one
72#    return dep*npsind(tth)             #not as good as 1-cos2Q
73       
74def makeMat(Angle,Axis):
75    '''Make rotation matrix from Angle and Axis
76
77    :param float Angle: in degrees
78    :param int Axis: 0 for rotation about x, 1 for about y, etc.
79    '''
80    cs = npcosd(Angle)
81    ss = npsind(Angle)
82    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
83    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
84                   
85def FitEllipse(xy):
86   
87    def ellipse_center(p):
88        ''' gives ellipse center coordinates
89        '''
90        b,c,d,f,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[0]
91        num = b*b-a*c
92        x0=(c*d-b*f)/num
93        y0=(a*f-b*d)/num
94        return np.array([x0,y0])
95   
96    def ellipse_angle_of_rotation( p ):
97        ''' gives rotation of ellipse major axis from x-axis
98        range will be -90 to 90 deg
99        '''
100        b,c,a = p[1]/2, p[2], p[0]
101        return 0.5*npatand(2*b/(a-c))
102   
103    def ellipse_axis_length( p ):
104        ''' gives ellipse radii in [minor,major] order
105        '''
106        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
107        up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
108        down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
109        down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
110        res1=np.sqrt(up/down1)
111        res2=np.sqrt(up/down2)
112        return np.array([ res2,res1])
113   
114    xy = np.array(xy)
115    x = np.asarray(xy.T[0])[:,np.newaxis]
116    y = np.asarray(xy.T[1])[:,np.newaxis]
117    D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
118    S = np.dot(D.T,D)
119    C = np.zeros([6,6])
120    C[0,2] = C[2,0] = 2; C[1,1] = -1
121    E, V =  nl.eig(np.dot(nl.inv(S), C))
122    n = np.argmax(np.abs(E))
123    a = V[:,n]
124    cent = ellipse_center(a)
125    phi = ellipse_angle_of_rotation(a)
126    radii = ellipse_axis_length(a)
127    phi += 90.
128    if radii[0] > radii[1]:
129        radii = [radii[1],radii[0]]
130        phi -= 90.
131    return cent,phi,radii
132
133def FitDetector(rings,varyList,parmDict,Print=True):
134    'Needs a doc string'
135       
136    def CalibPrint(ValSig,chisq,Npts):
137        print 'Image Parameters: chi**2: %12.3g, Np: %d'%(chisq,Npts)
138        ptlbls = 'names :'
139        ptstr =  'values:'
140        sigstr = 'esds  :'
141        for name,value,sig in ValSig:
142            ptlbls += "%s" % (name.rjust(12))
143            if name == 'phi':
144                ptstr += Fmt[name] % (value%360.)
145            else:
146                ptstr += Fmt[name] % (value)
147            if sig:
148                sigstr += Fmt[name] % (sig)
149            else:
150                sigstr += 12*' '
151        print ptlbls
152        print ptstr
153        print sigstr       
154       
155    def ellipseCalcD(B,xyd,varyList,parmDict):
156       
157        x,y,dsp = xyd
158        varyDict = dict(zip(varyList,B))
159        parms = {}
160        for parm in parmDict:
161            if parm in varyList:
162                parms[parm] = varyDict[parm]
163            else:
164                parms[parm] = parmDict[parm]
165        phi = parms['phi']-90.               #get rotation of major axis from tilt axis
166        tth = 2.0*npasind(parms['wave']/(2.*dsp))
167        phi0 = npatan2d(y-parms['det-Y'],x-parms['det-X'])
168        dxy = peneCorr(tth,parms['dep'],parms['dist'],parms['tilt'],phi0)
169        stth = npsind(tth)
170        cosb = npcosd(parms['tilt'])
171        tanb = nptand(parms['tilt'])       
172        tbm = nptand((tth-parms['tilt'])/2.)
173        tbp = nptand((tth+parms['tilt'])/2.)
174        d = parms['dist']+dxy
175        fplus = d*tanb*stth/(cosb+stth)
176        fminus = d*tanb*stth/(cosb-stth)
177        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
178        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
179        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
180        R1 = (vplus+vminus)/2.                                    #major axis
181        zdis = (fplus-fminus)/2.
182        Robs = np.sqrt((x-parms['det-X'])**2+(y-parms['det-Y'])**2)
183        rsqplus = R0**2+R1**2
184        rsqminus = R0**2-R1**2
185        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
186        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
187        P = 2.*R0**2*zdis*npcosd(phi0-phi)
188        Rcalc = (P+Q)/R
189        M = (Robs-Rcalc)*10.        #why 10? does make "chi**2" more reasonable
190        return M
191       
192    names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
193    fmt = ['%12.3f','%12.3f','%12.3f','%12.3f','%12.3f','%12.4f','%12.6f']
194    Fmt = dict(zip(names,fmt))
195    p0 = [parmDict[key] for key in varyList]
196    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
197    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0))   #reduced chi^2 = M/(Nobs-Nvar)
198    parmDict.update(zip(varyList,result[0]))
199    vals = list(result[0])
200    sig = list(np.sqrt(chisq*np.diag(result[1])))
201    sigList = np.zeros(7)
202    for i,name in enumerate(varyList):
203        sigList[i] = sig[varyList.index(name)]
204    ValSig = zip(varyList,vals,sig)
205    if Print:
206        CalibPrint(ValSig,chisq,rings.shape[0])
207    return [chisq,vals,sigList]
208
209def ImageLocalMax(image,w,Xpix,Ypix):
210    'Needs a doc string'
211    w2 = w*2
212    sizey,sizex = image.shape
213    xpix = int(Xpix)            #get reference corner of pixel chosen
214    ypix = int(Ypix)
215    if not w:
216        ZMax = np.sum(image[ypix-1:ypix+1,xpix-1:xpix+1])
217        return xpix,ypix,ZMax,0.0001
218    if (w2 < xpix < sizex-w2) and (w2 < ypix < sizey-w2) and image[ypix,xpix]:
219        ZMax = image[ypix-w:ypix+w,xpix-w:xpix+w]
220        Zmax = np.argmax(ZMax)
221        ZMin = image[ypix-w2:ypix+w2,xpix-w2:xpix+w2]
222        Zmin = np.argmin(ZMin)
223        xpix += Zmax%w2-w
224        ypix += Zmax/w2-w
225        return xpix,ypix,np.ravel(ZMax)[Zmax],max(0.0001,np.ravel(ZMin)[Zmin])   #avoid neg/zero minimum
226    else:
227        return 0,0,0,0     
228   
229def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
230    'Needs a doc string'
231    def ellipseC():
232        'compute estimate of ellipse circumference'
233        if radii[0] < 0:        #hyperbola
234            theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
235            print theta
236            return 0
237        apb = radii[1]+radii[0]
238        amb = radii[1]-radii[0]
239        return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
240       
241    cent,phi,radii = ellipse
242    cphi = cosd(phi-90.)        #convert to major axis rotation
243    sphi = sind(phi-90.)
244    ring = []
245    C = int(ellipseC())         #ring circumference
246    azm = []
247    for i in range(0,C,1):      #step around ring in 1mm increments
248        a = 360.*i/C
249        x = radii[1]*cosd(a)        #major axis
250        y = radii[0]*sind(a)
251        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
252        Y = (sphi*x+cphi*y+cent[1])*scaley
253        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
254        if I and J and float(I)/J > reject:
255            X += .5                             #set to center of pixel
256            Y += .5
257            X /= scalex                         #convert back to mm
258            Y /= scaley
259            if [X,Y,dsp] not in ring:           #no duplicates!
260                ring.append([X,Y,dsp])
261                azm.append(a)
262    if len(ring) < 10:
263        ring = []
264        azm = []
265    return ring,azm
266   
267def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
268    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
269    on output
270    radii[0] (b-minor axis) set < 0. for hyperbola
271   
272    '''
273    radii = [0,0]
274    stth = sind(tth)
275    cosb = cosd(tilt)
276    tanb = tand(tilt)
277    tbm = tand((tth-tilt)/2.)
278    tbp = tand((tth+tilt)/2.)
279    sinb = sind(tilt)
280    d = dist+dxy
281    if tth+abs(tilt) < 90.:      #ellipse
282        fplus = d*tanb*stth/(cosb+stth)
283        fminus = d*tanb*stth/(cosb-stth)
284        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
285        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
286        radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
287        radii[1] = (vplus+vminus)/2.                                    #major axis
288        zdis = (fplus-fminus)/2.
289    else:   #hyperbola!
290        f = d*abs(tanb)*stth/(cosb+stth)
291        v = d*(abs(tanb)+tand(tth-abs(tilt)))
292        delt = d*stth*(1.+stth*cosb)/(abs(sinb)*cosb*(stth+cosb))
293        eps = (v-f)/(delt-v)
294        radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.)                     #-minor axis
295        radii[1] = eps*(delt-f)/(eps**2-1.)                             #major axis
296        if tilt > 0:
297            zdis = f+radii[1]*eps
298        else:
299            zdis = -f
300#NB: zdis is || to major axis & phi is rotation of minor axis
301#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
302    elcent = [cent[0]+zdis*sind(phi),cent[1]-zdis*cosd(phi)]
303    return elcent,phi,radii
304   
305def GetEllipse(dsp,data):
306    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
307    as given in image controls dictionary (data) and a d-spacing (dsp)
308    '''
309    cent = data['center']
310    tilt = data['tilt']
311    phi = data['rotation']
312    dep = data.get('DetDepth',0.0)
313    tth = 2.0*asind(data['wavelength']/(2.*dsp))
314    dist = data['distance']
315    dxy = peneCorr(tth,dep,dist,tilt)
316    return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
317       
318def GetDetectorXY(dsp,azm,data):
319    'Needs a doc string'
320   
321    elcent,phi,radii = GetEllipse(dsp,data)
322    phi = data['rotation']-90.          #to give rotation of major axis
323    tilt = data['tilt']
324    dist = data['distance']
325    cent = data['center']
326    tth = 2.0*asind(data['wavelength']/(2.*dsp))
327    stth = sind(tth)
328    cosb = cosd(tilt)
329    if radii[0] > 0.:
330        sinb = sind(tilt)
331        tanb = tand(tilt)
332        fplus = dist*tanb*stth/(cosb+stth)
333        fminus = dist*tanb*stth/(cosb-stth)
334        zdis = (fplus-fminus)/2.
335        rsqplus = radii[0]**2+radii[1]**2
336        rsqminus = radii[0]**2-radii[1]**2
337        R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus
338        Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm-phi)**2)
339        P = 2.*radii[0]**2*zdis*cosd(azm-phi)
340        radius = (P+Q)/R
341        xy = np.array([radius*cosd(azm),radius*sind(azm)])
342        xy += cent
343    else:   #hyperbola - both branches (one is way off screen!)
344        sinb = abs(sind(tilt))
345        tanb = abs(tand(tilt))
346        f = dist*tanb*stth/(cosb+stth)
347        v = dist*(tanb+tand(tth-abs(tilt)))
348        delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
349        ecc = (v-f)/(delt-v)
350        R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm))
351        if tilt > 0.:
352            offset = 2.*radii[1]*ecc+f      #select other branch
353            xy = [-R*cosd(azm)-offset,R*sind(azm)]
354        else:
355            offset = -f
356            xy = [-R*cosd(azm)-offset,-R*sind(azm)]
357        xy = -np.array([xy[0]*cosd(phi)+xy[1]*sind(phi),xy[0]*sind(phi)-xy[1]*cosd(phi)])
358        xy += cent
359    return xy
360   
361def GetDetXYfromThAzm(Th,Azm,data):
362    '''Computes a detector position from a 2theta angle and an azimultal
363    angle (both in degrees)
364    '''
365    dsp = data['wavelength']/(2.0*npsind(Th))   
366    return GetDetectorXY(dsp,Azm,data)
367                   
368def GetTthAzmDsp(x,y,data): #expensive
369    '''Computes a 2theta, etc. from a detector position and calibration constants
370    - checked OK for ellipses & hyperbola
371    :returns: np.array(tth,azm,G,dsp) where tth is 2theta, azm is the azimutal angle,
372       G is ? and dsp is the d-space
373    '''
374    wave = data['wavelength']
375    cent = data['center']
376    tilt = data['tilt']
377    dist = data['distance']/cosd(tilt)
378    x0 = dist*tand(tilt)
379    phi = data['rotation']
380    dep = data['DetDepth']
381    azmthoff = data['azmthOff']
382    dx = np.array(x-cent[0],dtype=np.float32)
383    dy = np.array(y-cent[1],dtype=np.float32)
384    D = ((dx-x0)**2+dy**2+dist**2)      #sample to pixel distance
385    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
386    X = np.dot(X,makeMat(phi,2))
387    Z = np.dot(X,makeMat(tilt,0)).T[2]
388    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
389    dxy = peneCorr(tth,dep,dist,tilt,npatan2d(dy,dx))
390    DX = dist-Z+dxy
391    DY = np.sqrt(dx**2+dy**2-Z**2)
392    tth = npatan2d(DY,DX) 
393    dsp = wave/(2.*npsind(tth/2.))
394    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
395    G = D/dist**2       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
396    return np.array([tth,azm,G,dsp])
397   
398def GetTth(x,y,data):
399    'Give 2-theta value for detector x,y position; calibration info in data'
400    return GetTthAzmDsp(x,y,data)[0]
401   
402def GetTthAzm(x,y,data):
403    'Give 2-theta, azimuth values for detector x,y position; calibration info in data'
404    return GetTthAzmDsp(x,y,data)[0:2]
405   
406def GetTthAzmG(x,y,data):
407    '''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
408     calibration info in data - only used in integration
409    '''
410    'Needs a doc string - checked OK for ellipses & hyperbola'
411    tilt = data['tilt']
412    dist = data['distance']/npcosd(tilt)
413    x0 = data['distance']*nptand(tilt)
414    MN = -np.inner(makeMat(data['rotation'],2),makeMat(tilt,0))
415    distsq = data['distance']**2
416    dx = x-data['center'][0]
417    dy = y-data['center'][1]
418    G = ((dx-x0)**2+dy**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
419    X = np.dstack([dx.T,dy.T,np.zeros_like(dx.T)])
420    Z = np.dot(X,MN).T[2]
421    xyZ = dx**2+dy**2-Z**2   
422    tth = npatand(np.sqrt(xyZ)/(dist-Z))
423    dxy = peneCorr(tth,data['DetDepth'],dist,tilt,npatan2d(dy,dx))
424    tth = npatan2d(np.sqrt(xyZ),dist-Z+dxy) 
425    azm = (npatan2d(dy,dx)+data['azmthOff']+720.)%360.
426    return tth,azm,G
427
428def GetDsp(x,y,data):
429    'Give d-spacing value for detector x,y position; calibration info in data'
430    return GetTthAzmDsp(x,y,data)[3]
431       
432def GetAzm(x,y,data):
433    'Give azimuth value for detector x,y position; calibration info in data'
434    return GetTthAzmDsp(x,y,data)[1]
435   
436def meanAzm(a,b):
437    AZM = lambda a,b: npacosd(0.5*(npsind(2.*b)-npsind(2.*a))/(np.pi*(b-a)/180.))/2.
438    azm = AZM(a,b)
439#    quad = int((a+b)/180.)
440#    if quad == 1:
441#        azm = 180.-azm
442#    elif quad == 2:
443#        azm += 180.
444#    elif quad == 3:
445#        azm = 360-azm
446    return azm     
447       
448def ImageCompress(image,scale):
449    ''' Reduces size of image by selecting every n'th point
450    param: image array: original image
451    param: scale int: intervsl between selected points
452    returns: array: reduced size image
453    '''
454    if scale == 1:
455        return image
456    else:
457        return image[::scale,::scale]
458       
459def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
460    'Needs a doc string'
461    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
462    curr = np.array([dist,x,y])
463    return abs(avg-curr)/avg < .02
464
465def EdgeFinder(image,data):
466    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
467    Not currently used but might be useful in future?
468    '''
469    import numpy.ma as ma
470    Nx,Ny = data['size']
471    pixelSize = data['pixelSize']
472    edgemin = data['edgemin']
473    scalex = pixelSize[0]/1000.
474    scaley = pixelSize[1]/1000.   
475    tay,tax = np.mgrid[0:Nx,0:Ny]
476    tax = np.asfarray(tax*scalex,dtype=np.float32)
477    tay = np.asfarray(tay*scaley,dtype=np.float32)
478    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
479    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
480    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
481    return zip(tax,tay)
482   
483def MakeFrameMask(data,frame):
484    pixelSize = data['pixelSize']
485    scalex = pixelSize[0]/1000.
486    scaley = pixelSize[1]/1000.
487    blkSize = 512
488    Nx,Ny = data['size']
489    nXBlks = (Nx-1)/blkSize+1
490    nYBlks = (Ny-1)/blkSize+1
491    tam = ma.make_mask_none(data['size'])
492    for iBlk in range(nXBlks):
493        iBeg = iBlk*blkSize
494        iFin = min(iBeg+blkSize,Nx)
495        for jBlk in range(nYBlks):
496            jBeg = jBlk*blkSize
497            jFin = min(jBeg+blkSize,Ny)               
498            nI = iFin-iBeg
499            nJ = jFin-jBeg
500            tax,tay = np.mgrid[iBeg+0.5:iFin+.5,jBeg+.5:jFin+.5]         #bin centers not corners
501            tax = np.asfarray(tax*scalex,dtype=np.float32)
502            tay = np.asfarray(tay*scaley,dtype=np.float32)
503            tamp = ma.make_mask_none((1024*1024))
504            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
505                tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
506            if tamp.shape:
507                tamp = np.reshape(tamp[:nI*nJ],(nI,nJ))
508                tam[iBeg:iFin,jBeg:jFin] = ma.mask_or(tamp[0:nI,0:nJ],tam[iBeg:iFin,jBeg:jFin])
509            else:
510                tam[iBeg:iFin,jBeg:jFin] = True
511    return tam.T
512   
513def ImageRecalibrate(G2frame,data,masks):
514    '''Called to repeat the calibration on an image, usually called after
515    calibration is done initially to improve the fit.
516    '''
517    import ImageCalibrants as calFile
518    print 'Image recalibration:'
519    time0 = time.time()
520    pixelSize = data['pixelSize']
521    scalex = 1000./pixelSize[0]
522    scaley = 1000./pixelSize[1]
523    pixLimit = data['pixLimit']
524    cutoff = data['cutoff']
525    data['rings'] = []
526    data['ellipses'] = []
527    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
528        data['DetDepth'] /= data['distance']
529    if not data['calibrant']:
530        print 'no calibration material selected'
531        return []   
532    skip = data['calibskip']
533    dmin = data['calibdmin']
534    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
535    HKL = []
536    for bravais,sg,cell in zip(Bravais,SGs,Cells):
537        A = G2lat.cell2A(cell)
538        if sg:
539            SGData = G2spc.SpcGroup(sg)[1]
540            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
541            HKL += hkl
542        else:
543            hkl = G2lat.GenHBravais(dmin,bravais,A)
544            HKL += hkl
545    HKL = G2lat.sortHKLd(HKL,True,False)
546    varyList = [item for item in data['varyList'] if data['varyList'][item]]
547    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
548        'setdist':data.get('setdist',data['distance']),
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 []   
578       
579    rings = np.concatenate((data['rings']),axis=0)
580    [chisq,vals,sigList] = 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 [vals,varyList,sigList,parmDict]
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    ctth = npcosd(tth)
683#1st estimate of tilt; assume ellipse - don't know sign though
684    if varyDict['tilt']:
685        tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
686        if not tilt:
687            print 'WARNING - selected ring was fitted as a circle'
688            print ' - if detector was tilted we suggest you skip this ring - WARNING'
689    else:
690        tilt = data['tilt']
691#1st estimate of dist: sample to detector normal to plane
692    if varyDict['dist']:
693        data['distance'] = dist = radii[0]**2/(ttth*radii[1])
694    else:
695        dist = data['distance']
696    if varyDict['tilt']:
697#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
698        zdisp = radii[1]*ttth*tand(tilt)
699        zdism = radii[1]*ttth*tand(-tilt)
700#cone axis position; 2 choices. Which is right?     
701#NB: zdisp is || to major axis & phi is rotation of minor axis
702#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
703        centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
704        centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
705#check get same ellipse parms either way
706#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
707        fail = True
708        i2 = 1
709        while fail:
710            dsp = HKL[i2][3]
711            print '2nd ring: try %.4f'%(dsp)
712            tth = 2.0*asind(wave/(2.*dsp))
713            ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
714            print fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1])
715            Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
716            parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
717                'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}       
718            varyList = [item for item in varyDict if varyDict[item]]
719            if len(Ringp) > 10:
720                chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)[0]
721                tiltp = parmDict['tilt']
722                phip = parmDict['phi']
723                centp = [parmDict['det-X'],parmDict['det-Y']]
724                fail = False
725            else:
726                chip = 1e6
727            ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
728            print fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1])
729            Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
730            if len(Ringm) > 10:
731                parmDict['tilt'] *= -1
732                chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)[0]
733                tiltm = parmDict['tilt']
734                phim = parmDict['phi']
735                centm = [parmDict['det-X'],parmDict['det-Y']]
736                fail = False
737            else:
738                chim = 1e6
739            if fail:
740                i2 += 1
741        if chip < chim:
742            data['tilt'] = tiltp
743            data['center'] = centp
744            data['rotation'] = phip
745        else:
746            data['tilt'] = tiltm
747            data['center'] = centm
748            data['rotation'] = phim
749        data['ellipses'].append(ellipsep[:]+('b',))
750        data['rings'].append(np.array(Ringp))
751        data['ellipses'].append(ellipsem[:]+('r',))
752        data['rings'].append(np.array(Ringm))
753        G2plt.PlotImage(G2frame,newImage=True)
754    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
755        data['DetDepth'] /= data['distance']
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)[0]
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)[0]
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,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    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 data['DetDepth'] > 0.5:          #patch - redefine DetDepth
892        data['DetDepth'] /= data['distance']
893    if 'SASD' in data['type']:
894        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
895    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
896    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
897    Nx,Ny = data['size']
898    nXBlks = (Nx-1)/blkSize+1
899    nYBlks = (Ny-1)/blkSize+1
900    tbeg = time.time()
901    times = [0,0,0,0,0]
902    for iBlk in range(nYBlks):
903        iBeg = iBlk*blkSize
904        iFin = min(iBeg+blkSize,Ny)
905        for jBlk in range(nXBlks):
906            jBeg = jBlk*blkSize
907            jFin = min(jBeg+blkSize,Nx)
908            # next is most expensive step!
909            TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin),times)           #2-theta & azimuth arrays & create position mask
910            Block = image[iBeg:iFin,jBeg:jFin]
911            t0 = time.time()
912            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
913            del TA; del tam
914            times[2] += time.time()-t0
915            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
916            tax = np.where(tax < LRazm[0],tax+360.,tax)
917            if data.get('SampleAbs',[0.0,''])[1]:
918                if 'Cylind' in data['SampleShape']:
919                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))      #adjust for additional thickness off sample normal
920                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
921                elif 'Fixed' in data['SampleShape']:    #assumes flat plate sample normal to beam
922                    tabs = G2pwd.Absorb('Fixed',muT,tay)
923            if 'log(q)' in data.get('binType',''):
924                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
925            elif 'Q' == data.get('binType',''):
926                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
927            t0 = time.time()
928            taz = np.array((taz*tad/tabs),dtype='float32')
929            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
930                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,
931                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
932            times[3] += time.time()-t0
933            del tax; del tay; del taz; del tad; del tabs
934    t0 = time.time()
935    NST = np.array(NST,dtype=np.float)
936    H0 = np.divide(H0,NST)
937    H0 = np.nan_to_num(H0)
938    H2 = np.array([tth for tth in np.linspace(lutth[0],lutth[1],numChans+1)])
939    if 'log(q)' in data.get('binType',''):
940        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
941    elif 'Q' == data.get('binType',''):
942        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
943    if Dazm:       
944        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
945    else:
946        H1 = LRazm
947    H0 /= npcosd(H2[:-1])           #**2? I don't think so, **1 is right for powders
948    if 'SASD' in data['type']:
949        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
950    if data['Oblique'][1]:
951        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
952    if 'SASD' in data['type'] and data['PolaVal'][1]:
953        #NB: in G2pwd.Polarization azm is defined from plane of polarization, not image x axis!
954        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm-90.)[0] for azm in (H1[:-1]+np.diff(H1)/2.)])
955    times[4] += time.time()-t0
956    print 'Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
957        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4])
958    print "Elapsed time:","%8.3fs"%(time.time()-tbeg)
959    print 'Integration complete'
960    if returnN:     #As requested by Steven Weigand
961        return H0,H1,H2,NST,CancelPressed
962    else:
963        return H0,H1,H2,CancelPressed
964   
965def MakeStrStaRing(ring,Image,Controls):
966    ellipse = GetEllipse(ring['Dset'],Controls)
967    pixSize = Controls['pixelSize']
968    scalex = 1000./pixSize[0]
969    scaley = 1000./pixSize[1]
970    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
971    if len(Ring):
972        ring['ImxyObs'] = copy.copy(Ring[:2])
973        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
974        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
975        ring['ImtaObs'] = TA
976        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
977        Ring[0] = TA[0]
978        Ring[1] = TA[1]
979        return Ring,ring
980    else:
981        ring['ImxyObs'] = [[],[]]
982        ring['ImtaObs'] = [[],[]]
983        ring['ImtaCalc'] = [[],[]]
984        return [],[]    #bad ring; no points found
985   
986def FitStrSta(Image,StrSta,Controls):
987    'Needs a doc string'
988   
989    StaControls = copy.deepcopy(Controls)
990    phi = StrSta['Sample phi']
991    wave = Controls['wavelength']
992    pixelSize = Controls['pixelSize']
993    scalex = 1000./pixelSize[0]
994    scaley = 1000./pixelSize[1]
995    StaType = StrSta['Type']
996    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
997
998    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
999        dset = ring['Dset']
1000        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1001        if len(Ring):
1002            ring.update(R)
1003            p0 = ring['Emat']
1004            val,esd = FitStrain(Ring,p0,dset,wave,phi,StaType)
1005            ring['Emat'] = val
1006            ring['Esig'] = esd
1007            ellipse = FitEllipse(R['ImxyObs'].T)
1008            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1009            ringint = np.array([float(Image[int(y*scaley),int(x*scalex)]) for x,y in np.array(ringxy)[:,:2]])
1010            ringint /= np.mean(ringint)
1011            ring['Ivar'] = np.var(ringint)
1012            print 'Variance in normalized ring intensity: %.3f'%(ring['Ivar'])
1013    CalcStrSta(StrSta,Controls)
1014   
1015def IntStrSta(Image,StrSta,Controls):
1016    StaControls = copy.deepcopy(Controls)
1017    pixelSize = Controls['pixelSize']
1018    scalex = 1000./pixelSize[0]
1019    scaley = 1000./pixelSize[1]
1020    phi = StrSta['Sample phi']
1021    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1022    RingsAI = []
1023    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1024        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1025        if len(Ring):
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            print ' %s %.3f %s %.3f'%('d-spacing',ring['Dcalc'],'var(MRD):',np.var(ringint))
1031            RingsAI.append(np.array(zip(ringazm,ringint)).T)
1032#            GSASIIpath.IPyBreak()
1033    return RingsAI
1034   
1035def CalcStrSta(StrSta,Controls):
1036
1037    wave = Controls['wavelength']
1038    phi = StrSta['Sample phi']
1039    StaType = StrSta['Type']
1040    for ring in StrSta['d-zero']:
1041        Eij = ring['Emat']
1042        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
1043        th,azm = ring['ImtaObs']
1044        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
1045        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
1046        if StaType == 'True':
1047            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
1048        else:
1049            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
1050        dmin = np.min(ring['ImtaCalc'][0])
1051        dmax = np.max(ring['ImtaCalc'][0])
1052        if abs(Eij[0]) < abs(Eij[2]):         #tension
1053            ring['Dcalc'] = dmin+(dmax-dmin)/4.
1054        else:                       #compression
1055            ring['Dcalc'] = dmin+3.*(dmax-dmin)/4.
1056#        ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
1057
1058def calcFij(omg,phi,azm,th):
1059    '''    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
1060
1061    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
1062        90 when perp. to sample surface
1063    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
1064    :param azm: his chi = azimuth around incident beam
1065    :param th:  his theta = theta
1066    '''
1067    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
1068    b = -npcosd(azm)*npcosd(th)
1069    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
1070    d = a*npsind(phi)+b*npcosd(phi)
1071    e = a*npcosd(phi)-b*npsind(phi)
1072    Fij = np.array([
1073        [d**2,d*e,c*d],
1074        [d*e,e**2,c*e],
1075        [c*d,c*e,c**2]])
1076    return -Fij
1077
1078def FitStrain(rings,p0,dset,wave,phi,StaType):
1079    'Needs a doc string'
1080    def StrainPrint(ValSig,dset):
1081        print 'Strain tensor for Dset: %.6f'%(dset)
1082        ptlbls = 'names :'
1083        ptstr =  'values:'
1084        sigstr = 'esds  :'
1085        for name,fmt,value,sig in ValSig:
1086            ptlbls += "%s" % (name.rjust(12))
1087            ptstr += fmt % (value)
1088            if sig:
1089                sigstr += fmt % (sig)
1090            else:
1091                sigstr += 12*' '
1092        print ptlbls
1093        print ptstr
1094        print sigstr
1095       
1096    def strainCalc(p,xyd,dset,wave,phi,StaType):
1097        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1098        dspo,azm,dsp = xyd
1099        th = npasind(wave/(2.0*dspo))
1100        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1101        if StaType == 'True':
1102            dspc = dset*np.exp(V)
1103        else:
1104            dspc = dset*(V+1.)
1105        return dspo-dspc
1106       
1107    names = ['e11','e12','e22']
1108    fmt = ['%12.2f','%12.2f','%12.2f']
1109    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1110    vals = list(result[0])
1111    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1112    sig = list(np.sqrt(chisq*np.diag(result[1])))
1113    ValSig = zip(names,fmt,vals,sig)
1114    StrainPrint(ValSig,dset)
1115    return vals,sig
1116   
1117def AutoSpotMasks(Image,Masks,Controls):
1118    rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
1119    print 'auto spot search'
1120    pixelSize = Controls['pixelSize']
1121    spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
1122    indices = (-1,0,1)
1123    rolls = np.array([[ix,iy] for ix in indices for iy in indices])
1124    for roll in rolls:
1125        if np.any(roll):        #avoid [0,0]
1126            spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.))
1127    mags = spotMask[spotMask.nonzero()]
1128    indx = np.transpose(spotMask.nonzero())
1129    nx,ny = Image.shape
1130    jndx = []
1131    for [ind,mag] in zip(indx,mags):
1132        if (1 < ind[0] < nx-2) and (1 < ind[1] < ny-2):
1133            cent = np.zeros((3,3))
1134            cent[1,1] = mag
1135            msk = np.array(Image[ind[0]-2:ind[0]+3,ind[1]-2:ind[1]+3])
1136            msk = ma.array(msk,mask=(msk==mag))
1137            if mag > 1.5*ma.mean(msk):
1138                jndx.append([ind[1]+.5,ind[0]+.5])
1139    print 'Spots found: ',len(jndx)
1140    jndx = np.array(jndx)
1141    peaks = jndx*pixelSize/1000.
1142    tth = GetTth(peaks.T[0],peaks.T[1],Controls)
1143    histth,bins = np.histogram(tth,2500)
1144    for shft in [-1,1]:
1145        histmsk = ma.array(histth,mask=-(histth-np.roll(histth,shft)<2.))
1146    histmsk = ma.array(histmsk,mask=(histmsk>5))
1147    binmsk = np.array(ma.nonzero(histmsk))
1148    digts = np.digitize(tth,bins)-1
1149    PeaksList = []
1150    for digt,peak in zip(digts,peaks):
1151        if digt in binmsk:
1152            PeaksList.append(peak)
1153    #should be able to filter out spotty Bragg rings here
1154    PeaksList = np.array(PeaksList)
1155    Points = np.ones((PeaksList.shape[0],3))
1156    Points[:,:2] = PeaksList
1157    Masks['Points'] = list(Points)
1158    return None
1159                   
1160   
1161       
Note: See TracBrowser for help on using the repository browser.