source: trunk/GSASIIimage.py @ 2836

Last change on this file since 2836 was 2836, checked in by vondreele, 4 years ago

fixes to scripts & fxye importer

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