source: trunk/GSASIIimage.py @ 2556

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

reorder contents of PDF Controls data page to put sample stuff first, then files & finally PDF calc. controls

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