source: trunk/GSASIIimage.py @ 3104

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

achieve some speed (~10%) improvement in integration of images

  • 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 17:25:46 +0000 (Fri, 29 Sep 2017) $
5# $Author: vondreele $
6# $Revision: 3104 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 3104 2017-09-29 17:25:46Z 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: 3104 $")
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    X = np.dstack([dx.T,dy.T,np.zeros_like(dx.T)])
426    Z = np.dot(X,MN).T[2]
427    xyZ = dx**2+dy**2-Z**2   
428    tth = npatand(np.sqrt(xyZ)/(dist-Z))
429    dxy = peneCorr(tth,data['DetDepth'],dist,tilt,npatan2d(dy,dx))
430    tth = npatan2d(np.sqrt(xyZ),dist-Z+dxy) 
431    azm = (npatan2d(dy,dx)+data['azmthOff']+720.)%360.
432    return tth,azm,G
433
434def GetDsp(x,y,data):
435    'Give d-spacing value for detector x,y position; calibration info in data'
436    return GetTthAzmDsp(x,y,data)[3]
437       
438def GetAzm(x,y,data):
439    'Give azimuth value for detector x,y position; calibration info in data'
440    return GetTthAzmDsp(x,y,data)[1]
441   
442def meanAzm(a,b):
443    AZM = lambda a,b: npacosd(0.5*(npsind(2.*b)-npsind(2.*a))/(np.pi*(b-a)/180.))/2.
444    azm = AZM(a,b)
445#    quad = int((a+b)/180.)
446#    if quad == 1:
447#        azm = 180.-azm
448#    elif quad == 2:
449#        azm += 180.
450#    elif quad == 3:
451#        azm = 360-azm
452    return azm     
453       
454def ImageCompress(image,scale):
455    ''' Reduces size of image by selecting every n'th point
456    param: image array: original image
457    param: scale int: intervsl between selected points
458    returns: array: reduced size image
459    '''
460    if scale == 1:
461        return image
462    else:
463        return image[::scale,::scale]
464       
465def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
466    'Needs a doc string'
467    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
468    curr = np.array([dist,x,y])
469    return abs(avg-curr)/avg < .02
470
471def EdgeFinder(image,data):
472    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
473    Not currently used but might be useful in future?
474    '''
475    import numpy.ma as ma
476    Nx,Ny = data['size']
477    pixelSize = data['pixelSize']
478    edgemin = data['edgemin']
479    scalex = pixelSize[0]/1000.
480    scaley = pixelSize[1]/1000.   
481    tay,tax = np.mgrid[0:Nx,0:Ny]
482    tax = np.asfarray(tax*scalex,dtype=np.float32)
483    tay = np.asfarray(tay*scaley,dtype=np.float32)
484    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
485    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
486    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
487    return zip(tax,tay)
488   
489def MakeFrameMask(data,frame):
490    pixelSize = data['pixelSize']
491    scalex = pixelSize[0]/1000.
492    scaley = pixelSize[1]/1000.
493    blkSize = 512
494    Nx,Ny = data['size']
495    nXBlks = (Nx-1)/blkSize+1
496    nYBlks = (Ny-1)/blkSize+1
497    tam = ma.make_mask_none(data['size'])
498    for iBlk in range(nXBlks):
499        iBeg = iBlk*blkSize
500        iFin = min(iBeg+blkSize,Nx)
501        for jBlk in range(nYBlks):
502            jBeg = jBlk*blkSize
503            jFin = min(jBeg+blkSize,Ny)               
504            nI = iFin-iBeg
505            nJ = jFin-jBeg
506            tax,tay = np.mgrid[iBeg+0.5:iFin+.5,jBeg+.5:jFin+.5]         #bin centers not corners
507            tax = np.asfarray(tax*scalex,dtype=np.float32)
508            tay = np.asfarray(tay*scaley,dtype=np.float32)
509            tamp = ma.make_mask_none((1024*1024))
510            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
511                tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
512            if tamp.shape:
513                tamp = np.reshape(tamp[:nI*nJ],(nI,nJ))
514                tam[iBeg:iFin,jBeg:jFin] = ma.mask_or(tamp[0:nI,0:nJ],tam[iBeg:iFin,jBeg:jFin])
515            else:
516                tam[iBeg:iFin,jBeg:jFin] = True
517    return tam.T
518   
519def ImageRecalibrate(G2frame,data,masks):
520    '''Called to repeat the calibration on an image, usually called after
521    calibration is done initially to improve the fit.
522    '''
523    import ImageCalibrants as calFile
524    print 'Image recalibration:'
525    time0 = time.time()
526    pixelSize = data['pixelSize']
527    scalex = 1000./pixelSize[0]
528    scaley = 1000./pixelSize[1]
529    pixLimit = data['pixLimit']
530    cutoff = data['cutoff']
531    data['rings'] = []
532    data['ellipses'] = []
533    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
534        data['DetDepth'] /= data['distance']
535    if not data['calibrant']:
536        print 'no calibration material selected'
537        return []   
538    skip = data['calibskip']
539    dmin = data['calibdmin']
540    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
541    HKL = []
542    for bravais,sg,cell in zip(Bravais,SGs,Cells):
543        A = G2lat.cell2A(cell)
544        if sg:
545            SGData = G2spc.SpcGroup(sg)[1]
546            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
547            HKL += hkl
548        else:
549            hkl = G2lat.GenHBravais(dmin,bravais,A)
550            HKL += hkl
551    HKL = G2lat.sortHKLd(HKL,True,False)
552    varyList = [item for item in data['varyList'] if data['varyList'][item]]
553    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
554        'setdist':data.get('setdist',data['distance']),
555        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
556    Found = False
557    wave = data['wavelength']
558    frame = masks['Frames']
559    tam = ma.make_mask_none(G2frame.ImageZ.shape)
560    if frame:
561        tam = ma.mask_or(tam,MakeFrameMask(data,frame))
562    for iH,H in enumerate(HKL):
563        if debug:   print H
564        dsp = H[3]
565        tth = 2.0*asind(wave/(2.*dsp))
566        if tth+abs(data['tilt']) > 90.:
567            print 'next line is a hyperbola - search stopped'
568            break
569        ellipse = GetEllipse(dsp,data)
570        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(G2frame.ImageZ,mask=tam))[0]
571        if Ring:
572            if iH >= skip:
573                data['rings'].append(np.array(Ring))
574            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
575            Found = True
576        elif not Found:         #skipping inner rings, keep looking until ring found
577            continue
578        else:                   #no more rings beyond edge of detector
579            data['ellipses'].append([])
580            continue
581    if not data['rings']:
582        print 'no rings found; try lower Min ring I/Ib'
583        return []   
584       
585    rings = np.concatenate((data['rings']),axis=0)
586    [chisq,vals,sigList] = FitDetector(rings,varyList,parmDict)
587    data['wavelength'] = parmDict['wave']
588    data['distance'] = parmDict['dist']
589    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
590    data['rotation'] = np.mod(parmDict['phi'],360.0)
591    data['tilt'] = parmDict['tilt']
592    data['DetDepth'] = parmDict['dep']
593    data['chisq'] = chisq
594    N = len(data['ellipses'])
595    data['ellipses'] = []           #clear away individual ellipse fits
596    for H in HKL[:N]:
597        ellipse = GetEllipse(H[3],data)
598        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))   
599    print 'calibration time = %.3f'%(time.time()-time0)
600    G2plt.PlotImage(G2frame,newImage=True)       
601    return [vals,varyList,sigList,parmDict]
602           
603def ImageCalibrate(G2frame,data):
604    '''Called to perform an initial image calibration after points have been
605    selected for the inner ring.
606    '''
607    import copy
608    import ImageCalibrants as calFile
609    print 'Image calibration:'
610    time0 = time.time()
611    ring = data['ring']
612    pixelSize = data['pixelSize']
613    scalex = 1000./pixelSize[0]
614    scaley = 1000./pixelSize[1]
615    pixLimit = data['pixLimit']
616    cutoff = data['cutoff']
617    varyDict = data['varyList']
618    if varyDict['dist'] and varyDict['wave']:
619        print 'ERROR - you can not simultaneously calibrate distance and wavelength'
620        return False
621    if len(ring) < 5:
622        print 'ERROR - not enough inner ring points for ellipse'
623        return False
624       
625    #fit start points on inner ring
626    data['ellipses'] = []
627    data['rings'] = []
628    outE = FitEllipse(ring)
629    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
630    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi**2: %.3f, Np: %d'
631    if outE:
632        print fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1])
633        ellipse = outE
634    else:
635        return False
636       
637    #setup 360 points on that ring for "good" fit
638    data['ellipses'].append(ellipse[:]+('g',))
639    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
640    if Ring:
641        ellipse = FitEllipse(Ring)
642        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]    #do again
643        ellipse = FitEllipse(Ring)
644    else:
645        print '1st ring not sufficiently complete to proceed'
646        return False
647    if debug:
648        print fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],
649            ellipse[2][0],ellipse[2][1],0.,len(Ring))     #cent,phi,radii
650    data['ellipses'].append(ellipse[:]+('r',))
651    data['rings'].append(np.array(Ring))
652    G2plt.PlotImage(G2frame,newImage=True)
653   
654#setup for calibration
655    data['rings'] = []
656    if not data['calibrant']:
657        print 'no calibration material selected'
658        return True
659   
660    skip = data['calibskip']
661    dmin = data['calibdmin']
662#generate reflection set
663    Bravais,SGs,Cells = calFile.Calibrants[data['calibrant']][:3]
664    HKL = []
665    for bravais,sg,cell in zip(Bravais,SGs,Cells):
666        A = G2lat.cell2A(cell)
667        if sg:
668            SGData = G2spc.SpcGroup(sg)[1]
669            hkl = G2pwd.getHKLpeak(dmin,SGData,A)
670            HKL += hkl
671        else:
672            hkl = G2lat.GenHBravais(dmin,bravais,A)
673            HKL += hkl
674    HKL = G2lat.sortHKLd(HKL,True,False)[skip:]
675#set up 1st ring
676    elcent,phi,radii = ellipse              #from fit of 1st ring
677    dsp = HKL[0][3]
678    print '1st ring: try %.4f'%(dsp)
679    if varyDict['dist']:
680        wave = data['wavelength']
681        tth = 2.0*asind(wave/(2.*dsp))
682    else:   #varyDict['wave']!
683        dist = data['distance']
684        tth = npatan2d(radii[0],dist)
685        data['wavelength'] = wave =  2.0*dsp*sind(tth/2.0)
686    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
687    ttth = nptand(tth)
688    ctth = npcosd(tth)
689#1st estimate of tilt; assume ellipse - don't know sign though
690    if varyDict['tilt']:
691        tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
692        if not tilt:
693            print 'WARNING - selected ring was fitted as a circle'
694            print ' - if detector was tilted we suggest you skip this ring - WARNING'
695    else:
696        tilt = data['tilt']
697#1st estimate of dist: sample to detector normal to plane
698    if varyDict['dist']:
699        data['distance'] = dist = radii[0]**2/(ttth*radii[1])
700    else:
701        dist = data['distance']
702    if varyDict['tilt']:
703#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
704        zdisp = radii[1]*ttth*tand(tilt)
705        zdism = radii[1]*ttth*tand(-tilt)
706#cone axis position; 2 choices. Which is right?     
707#NB: zdisp is || to major axis & phi is rotation of minor axis
708#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
709        centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
710        centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
711#check get same ellipse parms either way
712#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
713        fail = True
714        i2 = 1
715        while fail:
716            dsp = HKL[i2][3]
717            print '2nd ring: try %.4f'%(dsp)
718            tth = 2.0*asind(wave/(2.*dsp))
719            ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
720            print fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1])
721            Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
722            parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
723                'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}       
724            varyList = [item for item in varyDict if varyDict[item]]
725            if len(Ringp) > 10:
726                chip = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)[0]
727                tiltp = parmDict['tilt']
728                phip = parmDict['phi']
729                centp = [parmDict['det-X'],parmDict['det-Y']]
730                fail = False
731            else:
732                chip = 1e6
733            ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
734            print fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1])
735            Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,G2frame.ImageZ)[0]
736            if len(Ringm) > 10:
737                parmDict['tilt'] *= -1
738                chim = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)[0]
739                tiltm = parmDict['tilt']
740                phim = parmDict['phi']
741                centm = [parmDict['det-X'],parmDict['det-Y']]
742                fail = False
743            else:
744                chim = 1e6
745            if fail:
746                i2 += 1
747        if chip < chim:
748            data['tilt'] = tiltp
749            data['center'] = centp
750            data['rotation'] = phip
751        else:
752            data['tilt'] = tiltm
753            data['center'] = centm
754            data['rotation'] = phim
755        data['ellipses'].append(ellipsep[:]+('b',))
756        data['rings'].append(np.array(Ringp))
757        data['ellipses'].append(ellipsem[:]+('r',))
758        data['rings'].append(np.array(Ringm))
759        G2plt.PlotImage(G2frame,newImage=True)
760    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
761        data['DetDepth'] /= data['distance']
762    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
763        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
764    varyList = [item for item in varyDict if varyDict[item]]
765    data['rings'] = []
766    data['ellipses'] = []
767    for i,H in enumerate(HKL):
768        dsp = H[3]
769        tth = 2.0*asind(wave/(2.*dsp))
770        if tth+abs(data['tilt']) > 90.:
771            print 'next line is a hyperbola - search stopped'
772            break
773        if debug:   print 'HKLD:',H[:4],'2-theta: %.4f'%(tth)
774        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
775        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
776        if debug:   print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1])
777        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,G2frame.ImageZ)[0]
778        if Ring:
779            data['rings'].append(np.array(Ring))
780            rings = np.concatenate((data['rings']),axis=0)
781            if i:
782                chisq = FitDetector(rings,varyList,parmDict,False)[0]
783                data['distance'] = parmDict['dist']
784                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
785                data['rotation'] = parmDict['phi']
786                data['tilt'] = parmDict['tilt']
787                data['DetDepth'] = parmDict['dep']
788                data['chisq'] = chisq
789                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
790                if debug:   print fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings))
791            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
792#            G2plt.PlotImage(G2frame,newImage=True)
793        else:
794            if debug:   print 'insufficient number of points in this ellipse to fit'
795#            break
796    G2plt.PlotImage(G2frame,newImage=True)
797    fullSize = len(G2frame.ImageZ)/scalex
798    if 2*radii[1] < .9*fullSize:
799        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
800    N = len(data['ellipses'])
801    if N > 2:
802        FitDetector(rings,varyList,parmDict)[0]
803        data['wavelength'] = parmDict['wave']
804        data['distance'] = parmDict['dist']
805        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
806        data['rotation'] = parmDict['phi']
807        data['tilt'] = parmDict['tilt']
808        data['DetDepth'] = parmDict['dep']
809    for H in HKL[:N]:
810        ellipse = GetEllipse(H[3],data)
811        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
812    print 'calibration time = %.3f'%(time.time()-time0)
813    G2plt.PlotImage(G2frame,newImage=True)       
814    return True
815   
816def Make2ThetaAzimuthMap(data,masks,iLim,jLim,tamp,times): #most expensive part of integration!
817    'Needs a doc string'
818    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
819    pixelSize = data['pixelSize']
820    scalex = pixelSize[0]/1000.
821    scaley = pixelSize[1]/1000.
822   
823    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
824    tax = np.asfarray(tax*scalex,dtype=np.float32).flatten()
825    tay = np.asfarray(tay*scaley,dtype=np.float32).flatten()
826    nI = iLim[1]-iLim[0]
827    nJ = jLim[1]-jLim[0]
828    t0 = time.time()
829    #make position masks here
830    frame = masks['Frames']
831    tam = ma.make_mask_none((nI*nJ))
832    if frame:
833        tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
834            tay,len(frame),frame,tamp)[:nI*nJ])-True)
835    polygons = masks['Polygons']
836    for polygon in polygons:
837        if polygon:
838            tam = ma.mask_or(tam,ma.make_mask(pm.polymask(nI*nJ,tax,
839                tay,len(polygon),polygon,tamp)[:nI*nJ]))
840    for X,Y,rsq in masks['Points'].T:
841        tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,rsq)))
842    if tam.shape: tam = np.reshape(tam,(nI,nJ))
843    times[0] += time.time()-t0
844    t0 = time.time()
845    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
846    times[1] += time.time()-t0
847    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
848    return np.array(TA),tam           #2-theta, azimuth & geom. corr. arrays & position mask
849
850def Fill2ThetaAzimuthMap(masks,TA,tam,image):
851    'Needs a doc string'
852    Zlim = masks['Thresholds'][1]
853    rings = masks['Rings']
854    arcs = masks['Arcs']
855    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2])))    #azimuth, 2-theta, dist
856    tax,tay,tad = np.dsplit(TA,3)    #azimuth, 2-theta, dist**2/d0**2
857    for tth,thick in rings:
858        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
859    for tth,azm,thick in arcs:
860        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
861        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
862        tam = ma.mask_or(tam.flatten(),tamt*tama)
863    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
864    tabs = np.ones_like(taz)
865    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
866    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))   #azimuth
867    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))   #2-theta
868    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))   #intensity
869    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))   #dist**2/d0**2
870    tabs = ma.compressed(ma.array(tabs.flatten(),mask=tam)) #ones - later used for absorption corr.
871    return tax,tay,taz,tad,tabs
872
873def ImageIntegrate(image,data,masks,blkSize=128,returnN=False):
874    'Integrate an image; called from OnIntegrateAll and OnIntegrate in G2imgGUI'    #for q, log(q) bins need data['binType']
875    import histogram2d as h2d
876    print 'Begin image integration'
877    CancelPressed = False
878    LUtth = np.array(data['IOtth'])
879    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
880    numAzms = data['outAzimuths']
881    numChans = data['outChannels']
882    Dazm = (LRazm[1]-LRazm[0])/numAzms
883    if '2-theta' in data.get('binType','2-theta'):
884        lutth = LUtth               
885    elif 'log(q)' in data['binType']:
886        lutth = np.log(4.*np.pi*npsind(LUtth/2.)/data['wavelength'])
887    elif 'q' == data['binType'].lower():
888        lutth = 4.*np.pi*npsind(LUtth/2.)/data['wavelength']
889    dtth = (lutth[1]-lutth[0])/numChans
890    muT = data.get('SampleAbs',[0.0,''])[0]
891    if data['DetDepth'] > 0.5:          #patch - redefine DetDepth
892        data['DetDepth'] /= data['distance']
893    if 'SASD' in data['type']:
894        muT = -np.log(muT)/2.       #Transmission to 1/2 thickness muT
895    Masks = copy.deepcopy(masks)
896    Masks['Points'] = np.array(Masks['Points']).T           #get spots as X,Y,R arrays
897    if np.any(masks['Points']):
898        Masks['Points'][2] = np.square(Masks['Points'][2]/2.)
899    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
900    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
901    Nx,Ny = data['size']
902    nXBlks = (Nx-1)/blkSize+1
903    nYBlks = (Ny-1)/blkSize+1
904    tbeg = time.time()
905    times = [0,0,0,0,0]
906    tamp = ma.make_mask_none((1024*1024))       #NB: this array size used in the fortran histogram2d
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),tamp,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.