source: trunk/GSASIIimage.py @ 2788

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

make I/Ib floats in ImageCalibrants?.py
fix issue for SASD seq. refinement
fix some issues with SASD image integrations

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