source: trunk/GSASIIimage.py @ 2671

Last change on this file since 2671 was 2671, checked in by toby, 5 years ago

Set calib d-min as option for xfer angles; fix trig error

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