source: trunk/GSASIIimage.py @ 3105

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

another slight integration improvement

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 48.8 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2017-09-29 18:10:11 +0000 (Fri, 29 Sep 2017) $
5# $Author: vondreele $
6# $Revision: 3105 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 3105 2017-09-29 18:10:11Z 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: 3105 $")
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
52nxs = np.newaxis
53debug = False
54   
55def pointInPolygon(pXY,xy):
56    'Needs a doc string'
57    #pXY - assumed closed 1st & last points are duplicates
58    Inside = False
59    N = len(pXY)
60    p1x,p1y = pXY[0]
61    for i in range(N+1):
62        p2x,p2y = pXY[i%N]
63        if (max(p1y,p2y) >= xy[1] > min(p1y,p2y)) and (xy[0] <= max(p1x,p2x)):
64            if p1y != p2y:
65                xinters = (xy[1]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
66            if p1x == p2x or xy[0] <= xinters:
67                Inside = not Inside
68        p1x,p1y = p2x,p2y
69    return Inside
70   
71def peneCorr(tth,dep,dist,tilt=0.,azm=0.):
72    'Needs a doc string'
73#    return dep*(1.-npcosd(abs(tilt*npsind(azm))-tth*npcosd(azm)))  #something wrong here
74    return dep*(1.-npcosd(tth))*dist**2/1000.         #best one
75#    return dep*npsind(tth)             #not as good as 1-cos2Q
76       
77def makeMat(Angle,Axis):
78    '''Make rotation matrix from Angle and Axis
79
80    :param float Angle: in degrees
81    :param int Axis: 0 for rotation about x, 1 for about y, etc.
82    '''
83    cs = npcosd(Angle)
84    ss = npsind(Angle)
85    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
86    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
87                   
88def FitEllipse(xy):
89   
90    def ellipse_center(p):
91        ''' gives ellipse center coordinates
92        '''
93        b,c,d,f,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[0]
94        num = b*b-a*c
95        x0=(c*d-b*f)/num
96        y0=(a*f-b*d)/num
97        return np.array([x0,y0])
98   
99    def ellipse_angle_of_rotation( p ):
100        ''' gives rotation of ellipse major axis from x-axis
101        range will be -90 to 90 deg
102        '''
103        b,c,a = p[1]/2, p[2], p[0]
104        return 0.5*npatand(2*b/(a-c))
105   
106    def ellipse_axis_length( p ):
107        ''' gives ellipse radii in [minor,major] order
108        '''
109        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
110        up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
111        down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
112        down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
113        res1=np.sqrt(up/down1)
114        res2=np.sqrt(up/down2)
115        return np.array([ res2,res1])
116   
117    xy = np.array(xy)
118    x = np.asarray(xy.T[0])[:,np.newaxis]
119    y = np.asarray(xy.T[1])[:,np.newaxis]
120    D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
121    S = np.dot(D.T,D)
122    C = np.zeros([6,6])
123    C[0,2] = C[2,0] = 2; C[1,1] = -1
124    E, V =  nl.eig(np.dot(nl.inv(S), C))
125    n = np.argmax(np.abs(E))
126    a = V[:,n]
127    cent = ellipse_center(a)
128    phi = ellipse_angle_of_rotation(a)
129    radii = ellipse_axis_length(a)
130    phi += 90.
131    if radii[0] > radii[1]:
132        radii = [radii[1],radii[0]]
133        phi -= 90.
134    return cent,phi,radii
135
136def FitDetector(rings,varyList,parmDict,Print=True):
137    'Needs a doc string'
138       
139    def CalibPrint(ValSig,chisq,Npts):
140        print 'Image Parameters: chi**2: %12.3g, Np: %d'%(chisq,Npts)
141        ptlbls = 'names :'
142        ptstr =  'values:'
143        sigstr = 'esds  :'
144        for name,value,sig in ValSig:
145            ptlbls += "%s" % (name.rjust(12))
146            if name == 'phi':
147                ptstr += Fmt[name] % (value%360.)
148            else:
149                ptstr += Fmt[name] % (value)
150            if sig:
151                sigstr += Fmt[name] % (sig)
152            else:
153                sigstr += 12*' '
154        print ptlbls
155        print ptstr
156        print sigstr       
157       
158    def ellipseCalcD(B,xyd,varyList,parmDict):
159       
160        x,y,dsp = xyd
161        varyDict = dict(zip(varyList,B))
162        parms = {}
163        for parm in parmDict:
164            if parm in varyList:
165                parms[parm] = varyDict[parm]
166            else:
167                parms[parm] = parmDict[parm]
168        phi = parms['phi']-90.               #get rotation of major axis from tilt axis
169        tth = 2.0*npasind(parms['wave']/(2.*dsp))
170        phi0 = npatan2d(y-parms['det-Y'],x-parms['det-X'])
171        dxy = peneCorr(tth,parms['dep'],parms['dist'],parms['tilt'],phi0)
172        stth = npsind(tth)
173        cosb = npcosd(parms['tilt'])
174        tanb = nptand(parms['tilt'])       
175        tbm = nptand((tth-parms['tilt'])/2.)
176        tbp = nptand((tth+parms['tilt'])/2.)
177        d = parms['dist']+dxy
178        fplus = d*tanb*stth/(cosb+stth)
179        fminus = d*tanb*stth/(cosb-stth)
180        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
181        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
182        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
183        R1 = (vplus+vminus)/2.                                    #major axis
184        zdis = (fplus-fminus)/2.
185        Robs = np.sqrt((x-parms['det-X'])**2+(y-parms['det-Y'])**2)
186        rsqplus = R0**2+R1**2
187        rsqminus = R0**2-R1**2
188        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
189        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
190        P = 2.*R0**2*zdis*npcosd(phi0-phi)
191        Rcalc = (P+Q)/R
192        M = (Robs-Rcalc)*10.        #why 10? does make "chi**2" more reasonable
193        return M
194       
195    names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
196    fmt = ['%12.3f','%12.3f','%12.3f','%12.3f','%12.3f','%12.4f','%12.6f']
197    Fmt = dict(zip(names,fmt))
198    p0 = [parmDict[key] for key in varyList]
199    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
200    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[0]-len(p0))   #reduced chi^2 = M/(Nobs-Nvar)
201    parmDict.update(zip(varyList,result[0]))
202    vals = list(result[0])
203    sig = list(np.sqrt(chisq*np.diag(result[1])))
204    sigList = np.zeros(7)
205    for i,name in enumerate(varyList):
206        sigList[i] = sig[varyList.index(name)]
207    ValSig = zip(varyList,vals,sig)
208    if Print:
209        CalibPrint(ValSig,chisq,rings.shape[0])
210    return [chisq,vals,sigList]
211
212def ImageLocalMax(image,w,Xpix,Ypix):
213    'Needs a doc string'
214    w2 = w*2
215    sizey,sizex = image.shape
216    xpix = int(Xpix)            #get reference corner of pixel chosen
217    ypix = int(Ypix)
218    if not w:
219        ZMax = np.sum(image[ypix-2:ypix+2,xpix-2:xpix+2])
220        return xpix,ypix,ZMax,0.0001
221    if (w2 < xpix < sizex-w2) and (w2 < ypix < sizey-w2) and image[ypix,xpix]:
222        ZMax = image[ypix-w:ypix+w,xpix-w:xpix+w]
223        Zmax = np.argmax(ZMax)
224        ZMin = image[ypix-w2:ypix+w2,xpix-w2:xpix+w2]
225        Zmin = np.argmin(ZMin)
226        xpix += Zmax%w2-w
227        ypix += Zmax/w2-w
228        return xpix,ypix,np.ravel(ZMax)[Zmax],max(0.0001,np.ravel(ZMin)[Zmin])   #avoid neg/zero minimum
229    else:
230        return 0,0,0,0     
231   
232def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
233    'Needs a doc string'
234    def ellipseC():
235        'compute estimate of ellipse circumference'
236        if radii[0] < 0:        #hyperbola
237            theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
238            print theta
239            return 0
240        apb = radii[1]+radii[0]
241        amb = radii[1]-radii[0]
242        return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
243       
244    cent,phi,radii = ellipse
245    cphi = cosd(phi-90.)        #convert to major axis rotation
246    sphi = sind(phi-90.)
247    ring = []
248    C = int(ellipseC())         #ring circumference
249    azm = []
250    for i in range(0,C,1):      #step around ring in 1mm increments
251        a = 360.*i/C
252        x = radii[1]*cosd(a-phi+90.)        #major axis
253        y = radii[0]*sind(a-phi+90.)
254        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
255        Y = (sphi*x+cphi*y+cent[1])*scaley
256        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
257        if I and J and float(I)/J > reject:
258            X += .5                             #set to center of pixel
259            Y += .5
260            X /= scalex                         #convert back to mm
261            Y /= scaley
262            if [X,Y,dsp] not in ring:           #no duplicates!
263                ring.append([X,Y,dsp])
264                azm.append(a)
265    if len(ring) < 10:
266        ring = []
267        azm = []
268    return ring,azm
269   
270def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
271    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
272    on output
273    radii[0] (b-minor axis) set < 0. for hyperbola
274   
275    '''
276    radii = [0,0]
277    stth = sind(tth)
278    cosb = cosd(tilt)
279    tanb = tand(tilt)
280    tbm = tand((tth-tilt)/2.)
281    tbp = tand((tth+tilt)/2.)
282    sinb = sind(tilt)
283    d = dist+dxy
284    if tth+abs(tilt) < 90.:      #ellipse
285        fplus = d*tanb*stth/(cosb+stth)
286        fminus = d*tanb*stth/(cosb-stth)
287        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
288        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
289        radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
290        radii[1] = (vplus+vminus)/2.                                    #major axis
291        zdis = (fplus-fminus)/2.
292    else:   #hyperbola!
293        f = d*abs(tanb)*stth/(cosb+stth)
294        v = d*(abs(tanb)+tand(tth-abs(tilt)))
295        delt = d*stth*(1.+stth*cosb)/(abs(sinb)*cosb*(stth+cosb))
296        eps = (v-f)/(delt-v)
297        radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.)                     #-minor axis
298        radii[1] = eps*(delt-f)/(eps**2-1.)                             #major axis
299        if tilt > 0:
300            zdis = f+radii[1]*eps
301        else:
302            zdis = -f
303#NB: zdis is || to major axis & phi is rotation of minor axis
304#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
305    elcent = [cent[0]+zdis*sind(phi),cent[1]-zdis*cosd(phi)]
306    return elcent,phi,radii
307   
308def GetEllipse(dsp,data):
309    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
310    as given in image controls dictionary (data) and a d-spacing (dsp)
311    '''
312    cent = data['center']
313    tilt = data['tilt']
314    phi = data['rotation']
315    dep = data.get('DetDepth',0.0)
316    tth = 2.0*asind(data['wavelength']/(2.*dsp))
317    dist = data['distance']
318    dxy = peneCorr(tth,dep,dist,tilt)
319    return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
320       
321def GetDetectorXY(dsp,azm,data):
322    '''Get detector x,y position from d-spacing (dsp), azimuth (azm,deg)
323    & image controls dictionary (data)
324    it seems to be only used in plotting
325    '''   
326    elcent,phi,radii = GetEllipse(dsp,data)
327    phi = data['rotation']-90.          #to give rotation of major axis
328    tilt = data['tilt']
329    dist = data['distance']
330    cent = data['center']
331    tth = 2.0*asind(data['wavelength']/(2.*dsp))
332    stth = sind(tth)
333    cosb = cosd(tilt)
334    if radii[0] > 0.:
335        sinb = sind(tilt)
336        tanb = tand(tilt)
337        fplus = dist*tanb*stth/(cosb+stth)
338        fminus = dist*tanb*stth/(cosb-stth)
339        zdis = (fplus-fminus)/2.
340        rsqplus = radii[0]**2+radii[1]**2
341        rsqminus = radii[0]**2-radii[1]**2
342        R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus
343        Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm-phi)**2)
344        P = 2.*radii[0]**2*zdis*cosd(azm-phi)
345        radius = (P+Q)/R
346        xy = np.array([radius*cosd(azm),radius*sind(azm)])
347        xy += cent
348    else:   #hyperbola - both branches (one is way off screen!)
349        sinb = abs(sind(tilt))
350        tanb = abs(tand(tilt))
351        f = dist*tanb*stth/(cosb+stth)
352        v = dist*(tanb+tand(tth-abs(tilt)))
353        delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
354        ecc = (v-f)/(delt-v)
355        R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm))
356        if tilt > 0.:
357            offset = 2.*radii[1]*ecc+f      #select other branch
358            xy = [-R*cosd(azm)-offset,-R*sind(azm)]
359        else:
360            offset = -f
361            xy = [-R*cosd(azm)-offset,R*sind(azm)]
362        xy = -np.array([xy[0]*cosd(phi)+xy[1]*sind(phi),xy[0]*sind(phi)-xy[1]*cosd(phi)])
363        xy += cent
364    return xy
365   
366def GetDetXYfromThAzm(Th,Azm,data):
367    '''Computes a detector position from a 2theta angle and an azimultal
368    angle (both in degrees) - apparently not used!
369    '''
370    dsp = data['wavelength']/(2.0*npsind(Th))   
371    return GetDetectorXY(dsp,Azm,data)
372                   
373def GetTthAzmDsp(x,y,data): #expensive
374    '''Computes a 2theta, etc. from a detector position and calibration constants - checked
375    OK for ellipses & hyperbola.
376
377    :returns: np.array(tth,azm,G,dsp) where tth is 2theta, azm is the azimutal angle,
378       G is ? and dsp is the d-space
379    '''
380    wave = data['wavelength']
381    cent = data['center']
382    tilt = data['tilt']
383    dist = data['distance']/cosd(tilt)
384    x0 = dist*tand(tilt)
385    phi = data['rotation']
386    dep = data.get('DetDepth',0.)
387    azmthoff = data['azmthOff']
388    dx = np.array(x-cent[0],dtype=np.float32)
389    dy = np.array(y-cent[1],dtype=np.float32)
390    D = ((dx-x0)**2+dy**2+dist**2)      #sample to pixel distance
391    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
392    X = np.dot(X,makeMat(phi,2))
393    Z = np.dot(X,makeMat(tilt,0)).T[2]
394    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
395    dxy = peneCorr(tth,dep,dist,tilt,npatan2d(dy,dx))
396    DX = dist-Z+dxy
397    DY = np.sqrt(dx**2+dy**2-Z**2)
398    tth = npatan2d(DY,DX) 
399    dsp = wave/(2.*npsind(tth/2.))
400    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
401    G = D/dist**2       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
402    return np.array([tth,azm,G,dsp])
403   
404def GetTth(x,y,data):
405    'Give 2-theta value for detector x,y position; calibration info in data'
406    return GetTthAzmDsp(x,y,data)[0]
407   
408def GetTthAzm(x,y,data):
409    'Give 2-theta, azimuth values for detector x,y position; calibration info in data'
410    return GetTthAzmDsp(x,y,data)[0:2]
411   
412def GetTthAzmG(x,y,data):
413    '''Give 2-theta, azimuth & geometric corr. values for detector x,y position;
414     calibration info in data - only used in integration
415    '''
416    'Needs a doc string - checked OK for ellipses & hyperbola'
417    tilt = data['tilt']
418    dist = data['distance']/npcosd(tilt)
419    x0 = data['distance']*nptand(tilt)
420    MN = -np.inner(makeMat(data['rotation'],2),makeMat(tilt,0))
421    distsq = data['distance']**2
422    dx = x-data['center'][0]
423    dy = y-data['center'][1]
424    G = ((dx-x0)**2+dy**2+distsq)/distsq       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
425    Z = np.dot(np.dstack([dx.T,dy.T,np.zeros_like(dx.T)]),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,tamp,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).flatten()
824    tay = np.asfarray(tay*scaley,dtype=np.float32).flatten()
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        tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
833            tay,len(frame),frame,tamp)[:nI*nJ])-True)
834    polygons = masks['Polygons']
835    for polygon in polygons:
836        if polygon:
837            tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
838                tay,len(polygon),polygon,tamp)[:nI*nJ]))
839    for X,Y,rsq in masks['Points'].T:
840        tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,rsq)))
841    if tam.shape: tam = np.reshape(tam,(nI,nJ))
842    times[0] += time.time()-t0
843    t0 = time.time()
844    TA = np.array(GetTthAzmG(np.reshape(tax,(nI,nJ)),np.reshape(tay,(nI,nJ)),data))     #includes geom. corr. as dist**2/d0**2 - most expensive step
845    times[1] += time.time()-t0
846    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
847    return np.array(TA),tam           #2-theta, azimuth & geom. corr. arrays & position mask
848
849def Fill2ThetaAzimuthMap(masks,TA,tam,image):
850    'Needs a doc string'
851    Zlim = masks['Thresholds'][1]
852    rings = masks['Rings']
853    arcs = masks['Arcs']
854    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2])))    #azimuth, 2-theta, dist
855    tax,tay,tad = np.dsplit(TA,3)    #azimuth, 2-theta, dist**2/d0**2
856    for tth,thick in rings:
857        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
858    for tth,azm,thick in arcs:
859        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
860        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
861        tam = ma.mask_or(tam.flatten(),tamt*tama)
862    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
863    tabs = np.ones_like(taz)
864    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
865    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
866    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
867    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
868    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
869    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
870    return tax,tay,taz,tad,tabs
871
872def ImageIntegrate(image,data,masks,blkSize=128,returnN=False):
873    'Integrate an image; called from OnIntegrateAll and OnIntegrate in G2imgGUI'    #for q, log(q) bins need data['binType']
874    import histogram2d as h2d
875    print 'Begin image integration'
876    CancelPressed = False
877    LUtth = np.array(data['IOtth'])
878    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
879    numAzms = data['outAzimuths']
880    numChans = data['outChannels']
881    Dazm = (LRazm[1]-LRazm[0])/numAzms
882    if '2-theta' in data.get('binType','2-theta'):
883        lutth = LUtth               
884    elif 'log(q)' in data['binType']:
885        lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
886    elif 'q' == data['binType'].lower():
887        lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
888    dtth = (lutth[1]-lutth[0])/numChans
889    muT = data.get('SampleAbs',[0.0,''])[0]
890    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
891        data['DetDepth'] /= data['distance']
892    if 'SASD' in data['type']:
893        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
894    Masks = copy.deepcopy(masks)
895    Masks['Points'] = np.array(Masks['Points']).T           #get spots as X,Y,R arrays
896    if np.any(masks['Points']):
897        Masks['Points'][2] = np.square(Masks['Points'][2]/2.)
898    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
899    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
900    Nx,Ny = data['size']
901    nXBlks = (Nx-1)/blkSize+1
902    nYBlks = (Ny-1)/blkSize+1
903    tbeg = time.time()
904    times = [0,0,0,0,0]
905    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
906    for iBlk in range(nYBlks):
907        iBeg = iBlk*blkSize
908        iFin = min(iBeg+blkSize,Ny)
909        for jBlk in range(nXBlks):
910            jBeg = jBlk*blkSize
911            jFin = min(jBeg+blkSize,Nx)
912            # next is most expensive step!
913            TA,tam = Make2ThetaAzimuthMap(data,Masks,(iBeg,iFin),(jBeg,jFin),tamp,times)           #2-theta & azimuth arrays & create position mask
914            Block = image[iBeg:iFin,jBeg:jFin]
915            t0 = time.time()
916            tax,tay,taz,tad,tabs = Fill2ThetaAzimuthMap(Masks,TA,tam,Block)    #and apply masks
917            del TA; del tam
918            times[2] += time.time()-t0
919            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
920            tax = np.where(tax < LRazm[0],tax+360.,tax)
921            if data.get('SampleAbs',[0.0,''])[1]:
922                if 'Cylind' in data['SampleShape']:
923                    muR = muT*(1.+npsind(tax)**2/2.)/(npcosd(tay))      #adjust for additional thickness off sample normal
924                    tabs = G2pwd.Absorb(data['SampleShape'],muR,tay)
925                elif 'Fixed' in data['SampleShape']:    #assumes flat plate sample normal to beam
926                    tabs = G2pwd.Absorb('Fixed',muT,tay)
927            if 'log(q)' in data.get('binType',''):
928                tay = np.log(4.*np.pi*npsind(tay/2.)/data['wavelength'])
929            elif 'q' == data.get('binType','').lower():
930                tay = 4.*np.pi*npsind(tay/2.)/data['wavelength']
931            t0 = time.time()
932            taz = np.array((taz*tad/tabs),dtype='float32')
933            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
934                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,
935                    numAzms,numChans,LRazm,lutth,Dazm,dtth,NST,H0)
936            times[3] += time.time()-t0
937            del tax; del tay; del taz; del tad; del tabs
938    t0 = time.time()
939    NST = np.array(NST,dtype=np.float)
940    H0 = np.divide(H0,NST)
941    H0 = np.nan_to_num(H0)
942    H2 = np.array([tth for tth in np.linspace(lutth[0],lutth[1],numChans+1)])
943    if 'log(q)' in data.get('binType',''):
944        H2 = 2.*npasind(np.exp(H2)*data['wavelength']/(4.*np.pi))
945    elif 'q' == data.get('binType','').lower():
946        H2 = 2.*npasind(H2*data['wavelength']/(4.*np.pi))
947    if Dazm:       
948        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
949    else:
950        H1 = LRazm
951    H0 /= npcosd(H2[:-1])           #**2? I don't think so, **1 is right for powders
952    if 'SASD' in data['type']:
953        H0 /= npcosd(H2[:-1])           #one more for small angle scattering data?
954    if data['Oblique'][1]:
955        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
956    if 'SASD' in data['type'] and data['PolaVal'][1]:
957        #NB: in G2pwd.Polarization azm is defined from plane of polarization, not image x axis!
958        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm-90.)[0] for azm in (H1[:-1]+np.diff(H1)/2.)])
959    times[4] += time.time()-t0
960    print 'Step times: \n apply masks  %8.3fs xy->th,azm   %8.3fs fill map     %8.3fs \
961        \n binning      %8.3fs cleanup      %8.3fs'%(times[0],times[1],times[2],times[3],times[4])
962    print "Elapsed time:","%8.3fs"%(time.time()-tbeg)
963    print 'Integration complete'
964    if returnN:     #As requested by Steven Weigand
965        return H0,H1,H2,NST,CancelPressed
966    else:
967        return H0,H1,H2,CancelPressed
968   
969def MakeStrStaRing(ring,Image,Controls):
970    ellipse = GetEllipse(ring['Dset'],Controls)
971    pixSize = Controls['pixelSize']
972    scalex = 1000./pixSize[0]
973    scaley = 1000./pixSize[1]
974    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
975    if len(Ring):
976        ring['ImxyObs'] = copy.copy(Ring[:2])
977        TA = GetTthAzm(Ring[0],Ring[1],Controls)       #convert x,y to tth,azm
978        TA[0] = Controls['wavelength']/(2.*npsind(TA[0]/2.))      #convert 2th to d
979        ring['ImtaObs'] = TA
980        ring['ImtaCalc'] = np.zeros_like(ring['ImtaObs'])
981        Ring[0] = TA[0]
982        Ring[1] = TA[1]
983        return Ring,ring
984    else:
985        ring['ImxyObs'] = [[],[]]
986        ring['ImtaObs'] = [[],[]]
987        ring['ImtaCalc'] = [[],[]]
988        return [],[]    #bad ring; no points found
989   
990def FitStrSta(Image,StrSta,Controls):
991    'Needs a doc string'
992   
993    StaControls = copy.deepcopy(Controls)
994    phi = StrSta['Sample phi']
995    wave = Controls['wavelength']
996    pixelSize = Controls['pixelSize']
997    scalex = 1000./pixelSize[0]
998    scaley = 1000./pixelSize[1]
999    StaType = StrSta['Type']
1000    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1001
1002    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1003        dset = ring['Dset']
1004        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1005        if len(Ring):
1006            ring.update(R)
1007            p0 = ring['Emat']
1008            val,esd,covMat = FitStrain(Ring,p0,dset,wave,phi,StaType)
1009            ring['Emat'] = val
1010            ring['Esig'] = esd
1011            ellipse = FitEllipse(R['ImxyObs'].T)
1012            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1013            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1014            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1015            ringint /= np.mean(ringint)
1016            ring['Ivar'] = np.var(ringint)
1017            ring['covMat'] = covMat
1018            print 'Variance in normalized ring intensity: %.3f'%(ring['Ivar'])
1019    CalcStrSta(StrSta,Controls)
1020   
1021def IntStrSta(Image,StrSta,Controls):
1022    StaControls = copy.deepcopy(Controls)
1023    pixelSize = Controls['pixelSize']
1024    scalex = 1000./pixelSize[0]
1025    scaley = 1000./pixelSize[1]
1026    phi = StrSta['Sample phi']
1027    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
1028    RingsAI = []
1029    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
1030        Ring,R = MakeStrStaRing(ring,Image,StaControls)
1031        if len(Ring):
1032            ellipse = FitEllipse(R['ImxyObs'].T)
1033            ringxy,ringazm = makeRing(ring['Dcalc'],ellipse,0,0.,scalex,scaley,Image)
1034            ring['ImxyCalc'] = np.array(ringxy).T[:2]
1035            ringint = np.array([float(Image[int(x*scalex),int(y*scaley)]) for y,x in np.array(ringxy)[:,:2]])
1036            ringint /= np.mean(ringint)
1037            print ' %s %.3f %s %.3f'%('d-spacing',ring['Dcalc'],'sig(MRD):',np.sqrt(np.var(ringint)))
1038            RingsAI.append(np.array(zip(ringazm,ringint)).T)
1039    return RingsAI
1040   
1041def CalcStrSta(StrSta,Controls):
1042
1043    wave = Controls['wavelength']
1044    phi = StrSta['Sample phi']
1045    StaType = StrSta['Type']
1046    for ring in StrSta['d-zero']:
1047        Eij = ring['Emat']
1048        E = [[Eij[0],Eij[1],0],[Eij[1],Eij[2],0],[0,0,0]]
1049        th,azm = ring['ImtaObs']
1050        th0 = np.ones_like(azm)*npasind(wave/(2.*ring['Dset']))
1051        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th0).T/1.e6,axis=2),axis=1)
1052        if StaType == 'True':
1053            ring['ImtaCalc'] = np.array([np.exp(V)*ring['Dset'],azm])
1054        else:
1055            ring['ImtaCalc'] = np.array([(V+1.)*ring['Dset'],azm])
1056        dmin = np.min(ring['ImtaCalc'][0])
1057        dmax = np.max(ring['ImtaCalc'][0])
1058        if ring.get('fixDset',True):
1059            if abs(Eij[0]) < abs(Eij[2]):         #tension
1060                ring['Dcalc'] = dmin+(dmax-dmin)/4.
1061            else:                       #compression
1062                ring['Dcalc'] = dmin+3.*(dmax-dmin)/4.
1063        else:
1064            ring['Dcalc'] = np.mean(ring['ImtaCalc'][0])
1065
1066def calcFij(omg,phi,azm,th):
1067    '''    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
1068
1069    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface,
1070        90 when perp. to sample surface
1071    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
1072    :param azm: his chi = azimuth around incident beam
1073    :param th:  his theta = theta
1074    '''
1075    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
1076    b = -npcosd(azm)*npcosd(th)
1077    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
1078    d = a*npsind(phi)+b*npcosd(phi)
1079    e = a*npcosd(phi)-b*npsind(phi)
1080    Fij = np.array([
1081        [d**2,d*e,c*d],
1082        [d*e,e**2,c*e],
1083        [c*d,c*e,c**2]])
1084    return -Fij
1085
1086def FitStrain(rings,p0,dset,wave,phi,StaType):
1087    'Needs a doc string'
1088    def StrainPrint(ValSig,dset):
1089        print 'Strain tensor for Dset: %.6f'%(dset)
1090        ptlbls = 'names :'
1091        ptstr =  'values:'
1092        sigstr = 'esds  :'
1093        for name,fmt,value,sig in ValSig:
1094            ptlbls += "%s" % (name.rjust(12))
1095            ptstr += fmt % (value)
1096            if sig:
1097                sigstr += fmt % (sig)
1098            else:
1099                sigstr += 12*' '
1100        print ptlbls
1101        print ptstr
1102        print sigstr
1103       
1104    def strainCalc(p,xyd,dset,wave,phi,StaType):
1105        E = np.array([[p[0],p[1],0],[p[1],p[2],0],[0,0,0]])
1106        dspo,azm,dsp = xyd
1107        th = npasind(wave/(2.0*dspo))
1108        V = -np.sum(np.sum(E*calcFij(90.,phi,azm,th).T/1.e6,axis=2),axis=1)
1109        if StaType == 'True':
1110            dspc = dset*np.exp(V)
1111        else:
1112            dspc = dset*(V+1.)
1113        return dspo-dspc
1114       
1115    names = ['e11','e12','e22']
1116    fmt = ['%12.2f','%12.2f','%12.2f']
1117    result = leastsq(strainCalc,p0,args=(rings,dset,wave,phi,StaType),full_output=True)
1118    vals = list(result[0])
1119    chisq = np.sum(result[2]['fvec']**2)/(rings.shape[1]-3)     #reduced chi^2 = M/(Nobs-Nvar)
1120    covM = result[1]
1121    covMatrix = covM*chisq
1122    sig = list(np.sqrt(chisq*np.diag(result[1])))
1123    ValSig = zip(names,fmt,vals,sig)
1124    StrainPrint(ValSig,dset)
1125    return vals,sig,covMatrix
1126   
1127def AutoSpotMasks(Image,Masks,Controls):
1128    print 'auto spot search'
1129    maxPks = 500
1130    dmax = 1.5         #just beyond diamond 111 @ 2.0596
1131    wave = Controls['wavelength']
1132    tthMax = 2.0*asind(wave/(2.*dmax))
1133    pixelSize = Controls['pixelSize']
1134    lap = np.array([[0,1,0],[1,-4,1],[0,1,0]],np.float32)
1135    ck = scsg.cspline2d(Image.T,24.0)
1136    spotMask1=-scsg.convolve2d(ck,lap,mode='same',boundary='symm')
1137    for mul in range(10):
1138        spotMask = ma.array(Image,mask=(spotMask1<(mul+1)*1.e4))
1139        indx = np.transpose(spotMask.nonzero())
1140        if not len(indx):       #smooth image - no sharp points
1141            break
1142        pts,distort = scvq.kmeans(np.array(indx,np.float32),maxPks,thresh=1.0)
1143        peaks = pts*pixelSize/1000.
1144        tths = GetTth(peaks.T[0],peaks.T[1],Controls)
1145        masktths = ma.getmask(ma.array(tths,mask=tths>tthMax))
1146        pts = pts[masktths,:]
1147        A,B = np.mgrid[0:len(pts),0:len(pts)]   
1148        M = np.sqrt(np.sum((pts[A]-pts[B])**2,axis=-1))
1149        MX = ma.array(M,mask=(M<10.))
1150        Indx = ma.nonzero(ma.getmask(MX))
1151        print 'Spots found: ',mul,len(pts),distort,len(Indx[0])
1152        if distort < .3:
1153            break
1154    if not len(Indx[0]):
1155        print 'no auto search spots found'
1156        return None
1157    if distort > 2.:
1158        print 'no big spots found'
1159        return None
1160    #use clustered points to get position & spread
1161   
1162    Indx = np.sort(Indx,axis=0).T
1163    Nmax = np.max(Indx)
1164    nums = [list(set([ind[1] for ind in Indx if ind[1]>=jnd and ind[0] == jnd])) for jnd in range(Nmax+1)]
1165    delList = []
1166    for Num in nums:
1167        delList += list(np.sort(Num))[1:]
1168    delList = list(set(delList))
1169    Nums = [num for inum,num in enumerate(nums) if inum not in delList]
1170    points = []
1171    esds = []
1172    for num in Nums:
1173        points.append(np.mean(pts[num],axis=0))
1174        esds.append(np.mean(np.var(pts[num],axis=0)))
1175    peaks = (np.array(points)+np.array([.5,.5]))*pixelSize/1000.
1176    Points = np.ones((peaks.shape[0],3))*2.
1177    Points[:,2] += np.array(esds)*pixelSize[0]/1000.
1178    Points[:,:2] = peaks
1179    Masks['Points'] = list(Points)
1180    return None
1181
1182#    rollImage = lambda rho,roll: np.roll(np.roll(rho,roll[0],axis=0),roll[1],axis=1)
1183#    print 'auto spot search'
1184#    pixelSize = Controls['pixelSize']
1185#    spotMask = ma.array(Image,mask=(Image<np.mean(Image)))
1186#    indices = (-1,0,1)
1187#    rolls = np.array([[ix,iy] for ix in indices for iy in indices])
1188#    for roll in rolls:
1189#        if np.any(roll):        #avoid [0,0]
1190#            spotMask = ma.array(spotMask,mask=(spotMask-rollImage(Image,roll)<0.))
1191#    mags = spotMask[spotMask.nonzero()]
1192#    indx = np.transpose(spotMask.nonzero())
1193#    nx,ny = Image.shape
1194#    jndx = []
1195#    for [ind,mag] in zip(indx,mags):
1196#        if (1 < ind[0] < nx-2) and (1 < ind[1] < ny-2):
1197#            cent = np.zeros((3,3))
1198#            cent[1,1] = mag
1199#            msk = np.array(Image[ind[0]-2:ind[0]+3,ind[1]-2:ind[1]+3])
1200#            msk = ma.array(msk,mask=(msk==mag))
1201#            if mag > 1.5*ma.mean(msk):
1202#                jndx.append([ind[1]+.5,ind[0]+.5])
1203#    print 'Spots found: ',len(jndx)
1204#    jndx = np.array(jndx)
1205#    peaks = jndx*pixelSize/1000.
1206#    tth = GetTth(peaks.T[0],peaks.T[1],Controls)
1207#    histth,bins = np.histogram(tth,2500)
1208#    for shft in [-1,1]:
1209#        histmsk = ma.array(histth,mask=-(histth-np.roll(histth,shft)<2.))
1210#    histmsk = ma.array(histmsk,mask=(histmsk>5))
1211#    binmsk = np.array(ma.nonzero(histmsk))
1212#    digts = np.digitize(tth,bins)-1
1213#    PeaksList = []
1214#    for digt,peak in zip(digts,peaks):
1215#        if digt in binmsk:
1216#            PeaksList.append(peak)
1217#    #should be able to filter out spotty Bragg rings here
1218#    PeaksList = np.array(PeaksList)
1219#    Points = np.ones((PeaksList.shape[0],3))
1220#    Points[:,:2] = PeaksList
1221#    Masks['Points'] = list(Points)
1222#    return None
1223       
1224                   
1225   
1226       
Note: See TracBrowser for help on using the repository browser.