source: trunk/GSASIIimage.py @ 1151

Last change on this file since 1151 was 1151, checked in by vondreele, 9 years ago

fix ring & arc masks
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 31.6 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2013-11-25 20:01:08 +0000 (Mon, 25 Nov 2013) $
5# $Author: vondreele $
6# $Revision: 1151 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 1151 2013-11-25 20:01:08Z 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
22from scipy.optimize import leastsq
23import copy
24import GSASIIpath
25GSASIIpath.SetVersionNumber("$Revision: 1151 $")
26import GSASIIplot as G2plt
27import GSASIIlattice as G2lat
28import GSASIIpwd as G2pwd
29import fellipse as fel
30
31# trig functions in degrees
32sind = lambda x: math.sin(x*math.pi/180.)
33asind = lambda x: 180.*math.asin(x)/math.pi
34tand = lambda x: math.tan(x*math.pi/180.)
35atand = lambda x: 180.*math.atan(x)/math.pi
36atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
37cosd = lambda x: math.cos(x*math.pi/180.)
38acosd = lambda x: 180.*math.acos(x)/math.pi
39rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
40#numpy versions
41npsind = lambda x: np.sin(x*np.pi/180.)
42npasind = lambda x: 180.*np.arcsin(x)/np.pi
43npcosd = lambda x: np.cos(x*np.pi/180.)
44npacosd = lambda x: 180.*np.arccos(x)/np.pi
45nptand = lambda x: np.tan(x*np.pi/180.)
46npatand = lambda x: 180.*np.arctan(x)/np.pi
47npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
48   
49def pointInPolygon(pXY,xy):
50    'Needs a doc string'
51    #pXY - assumed closed 1st & last points are duplicates
52    Inside = False
53    N = len(pXY)
54    p1x,p1y = pXY[0]
55    for i in range(N+1):
56        p2x,p2y = pXY[i%N]
57        if (max(p1y,p2y) >= xy[1] > min(p1y,p2y)) and (xy[0] <= max(p1x,p2x)):
58            if p1y != p2y:
59                xinters = (xy[1]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
60            if p1x == p2x or xy[0] <= xinters:
61                Inside = not Inside
62        p1x,p1y = p2x,p2y
63    return Inside
64   
65def peneCorr(tth,dep):
66    'Needs a doc string'
67    return dep*(1.-npcosd(tth))         #best one
68#    return dep*npsind(tth)             #not as good as 1-cos2Q
69       
70def makeMat(Angle,Axis):
71    '''Make rotation matrix from Angle and Axis
72
73    :param float Angle: in degrees
74    :param int Axis: 0 for rotation about x, 1 for about y, etc.
75    '''
76    cs = cosd(Angle)
77    ss = sind(Angle)
78    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
79    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
80                   
81def FitRing(ring,delta):
82    'Needs a doc string'
83    parms = []
84    if delta:
85        err,parms = FitEllipse(ring)
86        errc,parmsc = FitCircle(ring)
87        errc = errc[0]/(len(ring)*parmsc[2][0]**2)
88        if not parms or errc < .1:
89            parms = parmsc
90    else:
91        err,parms = FitCircle(ring)
92    return parms,err
93       
94def FitCircle(ring):
95    'Needs a doc string'
96   
97    def makeParmsCircle(B):
98        cent = [-B[0]/2,-B[1]/2]
99        phi = 0.
100        sr1 = sr2 = math.sqrt(cent[0]**2+cent[1]**2-B[2])
101        return cent,phi,[sr1,sr2]
102       
103    ring = np.array(ring)
104    x = np.asarray(ring.T[0])
105    y = np.asarray(ring.T[1])
106   
107    M = np.array((x,y,np.ones_like(x)))
108    B = np.array(-(x**2+y**2))
109    result = nl.lstsq(M.T,B)
110    return result[1],makeParmsCircle(result[0])
111       
112def FitEllipse(ring):
113    'Needs a doc string'
114           
115    def makeParmsEllipse(B):
116        det = 4.*(1.-B[0]**2)-B[1]**2
117        if det < 0.:
118            print 'hyperbola!'
119            print B
120            return 0
121        elif det == 0.:
122            print 'parabola!'
123            print B
124            return 0
125        cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \
126            (B[1]*B[2]-2.*(1.+B[0])*B[3])/det]
127        phi = 0.5*atand(0.5*B[1]/B[0])
128       
129        a = (1.+B[0])/cosd(2*phi)
130        b = 2.-a
131        f = (1.+B[0])*cent[0]**2+(1.-B[0])*cent[1]**2+B[1]*cent[0]*cent[1]-B[4]
132        if f/a < 0 or f/b < 0:
133            return 0
134        sr1 = math.sqrt(f/a)
135        sr2 = math.sqrt(f/b)
136        if sr1 > sr2:
137            sr1,sr2 = sr2,sr1
138            phi -= 90.
139            if phi < -90.:
140                phi += 180.
141        return cent,phi,[sr1,sr2]
142               
143    ring = np.array(ring)
144    x = np.asarray(ring.T[0])
145    y = np.asarray(ring.T[1])
146    M = np.array((x**2-y**2,x*y,x,y,np.ones_like(x)))
147    B = np.array(-(x**2+y**2))
148    bb,err = fel.fellipse(len(x),x,y,1.E-7)
149#    print nl.lstsq(M.T,B)[0]
150#    print bb
151    return err,makeParmsEllipse(bb)
152   
153def FitDetector(rings,varyList,parmDict,Print=True):
154    'Needs a doc string'
155       
156    def CalibPrint(ValSig):
157        print 'Image Parameters:'
158        ptlbls = 'names :'
159        ptstr =  'values:'
160        sigstr = 'esds  :'
161        for name,fmt,value,sig in ValSig:
162            ptlbls += "%s" % (name.rjust(12))
163            if name == 'rotate':
164                ptstr += fmt % (value-90.)      #kluge to get rotation from vertical - see GSASIIimgGUI
165            else:
166                ptstr += fmt % (value)
167            if sig:
168                sigstr += fmt % (sig)
169            else:
170                sigstr += 12*' '
171        print ptlbls
172        print ptstr
173        print sigstr       
174       
175    def ellipseCalcD(B,xyd,varyList,parmDict):
176        x,y,dsp = xyd
177        wave = parmDict['wave']
178        if 'dep' in varyList:
179            dist,x0,y0,tilt,phi,dep = B[:6]
180        else:
181            dist,x0,y0,tilt,phi = B[:5]
182            dep = parmDict['dep']
183        if 'wave' in varyList:
184            wave = B[-1]
185        tth = 2.0*npasind(wave/(2.*dsp))
186        dxy = peneCorr(tth,dep)
187        ttth = nptand(tth)
188        stth = npsind(tth)
189        cosb = npcosd(tilt)
190        tanb = nptand(tilt)       
191        tbm = nptand((tth-tilt)/2.)
192        tbp = nptand((tth+tilt)/2.)
193        sinb = npsind(tilt)
194        d = dist+dxy
195        fplus = d*tanb*stth/(cosb+stth)
196        fminus = d*tanb*stth/(cosb-stth)
197        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
198        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
199        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
200        R1 = (vplus+vminus)/2.                                    #major axis
201        zdis = (fplus-fminus)/2.
202        X = x-x0-zdis*npsind(phi)
203        Y = y-y0-zdis*npcosd(phi)
204        XR = X*npcosd(phi)+Y*npsind(phi)
205        YR = X*npsind(phi)+Y*npcosd(phi)
206        return (XR/R0)**2+(YR/R1)**2-1
207       
208    names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
209    fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.3f','%12.5f']
210    p0 = [parmDict[key] for key in varyList]
211    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
212    chisq = np.sum(result[2]['fvec']**2)
213    parmDict.update(zip(varyList,result[0]))
214    vals = list(result[0])
215    chi = np.sqrt(np.sum(ellipseCalcD(result[0],rings.T,varyList,parmDict)**2))
216    sig = list(chi*np.sqrt(np.diag(result[1])))
217    sigList = np.zeros(7)
218    for i,name in enumerate(names):
219        if name in varyList:
220            sigList[i] = sig[varyList.index(name)]
221    ValSig = zip(names,fmt,vals,sig)
222    if Print:
223        CalibPrint(ValSig)
224    return chi,chisq
225                   
226def ImageLocalMax(image,w,Xpix,Ypix):
227    'Needs a doc string'
228    w2 = w*2
229    sizey,sizex = image.shape
230    xpix = int(Xpix)            #get reference corner of pixel chosen
231    ypix = int(Ypix)
232    if (w < xpix < sizex-w) and (w < ypix < sizey-w) and image[ypix,xpix]:
233        Z = image[ypix-w:ypix+w,xpix-w:xpix+w]
234        Zmax = np.argmax(Z)
235        Zmin = np.argmin(Z)
236        xpix += Zmax%w2-w
237        ypix += Zmax/w2-w
238        return xpix,ypix,np.ravel(Z)[Zmax],max(0.0001,np.ravel(Z)[Zmin])   #avoid neg/zero minimum
239    else:
240        return 0,0,0,0     
241   
242def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
243    'Needs a doc string'
244    def ellipseC():
245        'compute estimate of ellipse circumference'
246        if radii[0] < 0:        #hyperbola
247            theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
248            print theta
249            return 0
250        apb = radii[1]+radii[0]
251        amb = radii[1]-radii[0]
252        return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
253    cent,phi,radii = ellipse
254    cphi = cosd(phi)
255    sphi = sind(phi)
256    ring = []
257    sumI = 0
258    amin = 0
259    amax = 360
260    for a in range(0,int(ellipseC()),1):
261        x = radii[0]*cosd(a)
262        y = radii[1]*sind(a)
263        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
264        Y = (sphi*x+cphi*y+cent[1])*scaley
265        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
266        if I and J and I/J > reject:
267            sumI += I/J
268            X += .5                             #set to center of pixel
269            Y += .5
270            X /= scalex                         #convert to mm
271            Y /= scaley
272            amin = min(amin,a)
273            amax = max(amax,a)
274            if [X,Y,dsp] not in ring:
275                ring.append([X,Y,dsp])
276    delt = amax-amin
277    if len(ring) < 10:             #want more than 10 deg
278        return [],(delt > 90),0
279    return ring,(delt > 90),sumI/len(ring)
280   
281def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
282    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
283    on output
284    radii[0] (b-minor axis) set < 0. for hyperbola
285   
286    '''
287    radii = [0,0]
288    ttth = tand(tth)
289    stth = sind(tth)
290    ctth = cosd(tth)
291    cosb = cosd(tilt)
292    tanb = tand(tilt)
293    tbm = tand((tth-tilt)/2.)
294    tbp = tand((tth+tilt)/2.)
295    sinb = sind(tilt)
296    d = dist+dxy
297    if tth+tilt < 90.:      #ellipse
298        fplus = d*tanb*stth/(cosb+stth)
299        fminus = d*tanb*stth/(cosb-stth)
300        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
301        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
302        radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
303        radii[1] = (vplus+vminus)/2.                                    #major axis
304        zdis = (fplus-fminus)/2.
305    else:   #hyperbola!
306        f = d*tanb*stth/(cosb+stth)
307        v = d*(tanb+tand(tth-tilt))
308        delt = d*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
309        eps = (v-f)/(delt-v)
310        radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.)                     #-minor axis
311        radii[1] = eps*(delt-f)/(eps**2-1.)                             #major axis
312        zdis = f+radii[1]*eps
313    elcent = [cent[0]+zdis*sind(phi),cent[1]+zdis*cosd(phi)]
314    return elcent,phi,radii
315   
316def GetEllipse(dsp,data):
317    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
318    as given in image controls dictionary (data) and a d-spacing (dsp)
319    '''
320    cent = data['center']
321    tilt = data['tilt']
322    phi = data['rotation']
323    dep = data['DetDepth']
324    tth = 2.0*asind(data['wavelength']/(2.*dsp))
325    dxy = peneCorr(tth,dep)
326    dist = data['distance']
327    return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
328       
329def GetDetectorXY(dsp,azm,data):
330    'Needs a doc string'
331   
332    from scipy.optimize import fsolve
333    def func(xy,*args):
334       azm,phi,R0,R1,A,B = args
335       cp = cosd(phi)
336       sp = sind(phi)
337       x,y = xy
338       out = []
339       out.append(y-x*tand(azm))
340       out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)
341       return out
342    elcent,phi,radii = GetEllipse(dsp,data)
343    cent = data['center']
344    phi = data['rotation']
345    wave = data['wavelength']
346    tilt = data['tilt']
347    dist = data['distance']/cosd(tilt)
348    dep = data['DetDepth']
349    tth = 2.0*asind(wave/(2.*dsp))
350    dxy = peneCorr(tth,dep)
351    ttth = tand(tth)
352    radius = (dist+dxy)*ttth
353    stth = sind(tth)
354    cosb = cosd(tilt)
355    R1 = (dist+dxy)*stth*cosd(tth)*cosb/(cosb**2-stth**2)
356    R0 = math.sqrt(R1*radius*cosb)
357    zdis = R1*ttth*tand(tilt)
358    A = zdis*sind(phi)
359    B = -zdis*cosd(phi)
360    xy0 = [radius*cosd(azm),radius*sind(azm)]
361    xy = fsolve(func,xy0,args=(azm,phi,R0,R1,A,B))+cent
362    return xy
363   
364def GetDetXYfromThAzm(Th,Azm,data):
365    'Needs a doc string'
366    dsp = data['wavelength']/(2.0*npsind(Th))   
367    return GetDetectorXY(dsp,azm,data)
368                   
369def GetTthAzmDsp(x,y,data):
370    'Needs a doc string - checked OK for ellipses dont know about hyperbola'
371    wave = data['wavelength']
372    cent = data['center']
373    tilt = data['tilt']
374    dist = data['distance']/cosd(tilt)
375    phi = data['rotation']
376    dep = data['DetDepth']
377    LRazim = data['LRazimuth']
378    azmthoff = data['azmthOff']
379    dx = np.array(x-cent[0],dtype=np.float32)
380    dy = np.array(y-cent[1],dtype=np.float32)
381    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
382    X = np.dot(X,makeMat(phi,2))
383    Z = np.dot(X,makeMat(tilt,0)).T[2]
384    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
385    dxy = peneCorr(tth,dep)
386    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z+dxy)) #depth corr not correct for tilted detector
387    dsp = wave/(2.*npsind(tth/2.))
388    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
389    return tth,azm,dsp
390   
391def GetTth(x,y,data):
392    'Needs a doc string'
393    return GetTthAzmDsp(x,y,data)[0]
394   
395def GetTthAzm(x,y,data):
396    'Needs a doc string'
397    return GetTthAzmDsp(x,y,data)[0:2]
398   
399def GetDsp(x,y,data):
400    'Needs a doc string'
401    return GetTthAzmDsp(x,y,data)[2]
402       
403def GetAzm(x,y,data):
404    'Needs a doc string'
405    return GetTthAzmDsp(x,y,data)[1]
406       
407def ImageCompress(image,scale):
408    'Needs a doc string'
409    if scale == 1:
410        return image
411    else:
412        return image[::scale,::scale]
413       
414def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
415    'Needs a doc string'
416    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
417    curr = np.array([dist,x,y])
418    return abs(avg-curr)/avg < .02
419
420def EdgeFinder(image,data):
421    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
422    Not currently used but might be useful in future?
423    '''
424    import numpy.ma as ma
425    Nx,Ny = data['size']
426    pixelSize = data['pixelSize']
427    edgemin = data['edgemin']
428    scalex = pixelSize[0]/1000.
429    scaley = pixelSize[1]/1000.   
430    tay,tax = np.mgrid[0:Nx,0:Ny]
431    tax = np.asfarray(tax*scalex,dtype=np.float32)
432    tay = np.asfarray(tay*scaley,dtype=np.float32)
433    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
434    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
435    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
436    return zip(tax,tay)
437   
438def ImageRecalibrate(self,data):
439    'Needs a doc string'
440    import ImageCalibrants as calFile
441    print 'Image recalibration:'
442    time0 = time.time()
443    pixelSize = data['pixelSize']
444    scalex = 1000./pixelSize[0]
445    scaley = 1000./pixelSize[1]
446    pixLimit = data['pixLimit']
447    cutoff = data['cutoff']
448    data['rings'] = []
449    data['ellipses'] = []
450    if not data['calibrant']:
451        print 'no calibration material selected'
452        return True
453   
454    skip = data['calibskip']
455    dmin = data['calibdmin']
456    Bravais,Cells = calFile.Calibrants[data['calibrant']][:2]
457    HKL = []
458    for bravais,cell in zip(Bravais,Cells):
459        A = G2lat.cell2A(cell)
460        hkl = G2lat.GenHBravais(dmin,bravais,A)
461        HKL += hkl
462    HKL = G2lat.sortHKLd(HKL,True,False)
463    varyList = ['dist','det-X','det-Y','tilt','phi']
464    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
465        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
466    Found = False
467    wave = data['wavelength']
468    for H in HKL: 
469        dsp = H[3]
470        tth = 2.0*asind(wave/(2.*dsp))
471        if tth+abs(data['tilt']) > 90.:
472            print 'next line is a hyperbola - search stopped'
473            break
474        ellipse = GetEllipse(dsp,data)
475        Ring,delt,sumI = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
476        if Ring:
477            data['rings'].append(np.array(Ring))
478            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
479            Found = True
480        elif not Found:         #skipping inner rings, keep looking until ring found
481            continue
482        else:                   #no more rings beyond edge of detector
483            break
484    rings = np.concatenate((data['rings']),axis=0)
485    if data['DetDepthRef']:
486        varyList.append('dep')
487    FitDetector(rings,varyList,parmDict)
488    data['distance'] = parmDict['dist']
489    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
490    data['rotation'] = np.mod(parmDict['phi'],360.0)
491    data['tilt'] = parmDict['tilt']
492    data['DetDepth'] = parmDict['dep']
493    N = len(data['ellipses'])
494    data['ellipses'] = []           #clear away individual ellipse fits
495    for H in HKL[:N]:
496        ellipse = GetEllipse(H[3],data)
497        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
498   
499    print 'calibration time = ',time.time()-time0
500    G2plt.PlotImage(self,newImage=True)       
501    return True
502           
503def ImageCalibrate(self,data):
504    'Needs a doc string'
505    import copy
506    import ImageCalibrants as calFile
507    print 'Image calibration:'
508    time0 = time.time()
509    ring = data['ring']
510    pixelSize = data['pixelSize']
511    scalex = 1000./pixelSize[0]
512    scaley = 1000./pixelSize[1]
513    pixLimit = data['pixLimit']
514    cutoff = data['cutoff']
515    if len(ring) < 5:
516        print 'not enough inner ring points for ellipse'
517        return False
518       
519    #fit start points on inner ring
520    data['ellipses'] = []
521    outE,err = FitRing(ring,True)
522    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
523    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi^2: %.3f, N: %d'
524    if outE:
525        print fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1])
526        ellipse = outE
527    else:
528        return False
529       
530    #setup 360 points on that ring for "good" fit
531    Ring,delt,sumI = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
532    if Ring:
533        ellipse,err = FitRing(Ring,delt)
534        Ring,delt,sumI = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
535        ellipse,err = FitRing(Ring,delt)
536    else:
537        print '1st ring not sufficiently complete to proceed'
538        return False
539    print fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],ellipse[2][0],ellipse[2][1],err,len(Ring))     #cent,phi,radii
540    data['ellipses'].append(ellipse[:]+('r',))
541    data['rings'].append(np.array(Ring))
542    G2plt.PlotImage(self,newImage=True)
543   
544    #setup for calibration
545    data['rings'] = []
546    data['ellipses'] = []
547    if not data['calibrant']:
548        print 'no calibration material selected'
549        return True
550   
551    skip = data['calibskip']
552    dmin = data['calibdmin']
553#generate reflection set
554    Bravais,Cells = calFile.Calibrants[data['calibrant']][:2]
555    HKL = []
556    for bravais,cell in zip(Bravais,Cells):
557        A = G2lat.cell2A(cell)
558        hkl = G2lat.GenHBravais(dmin,bravais,A)[skip:]
559        HKL += hkl
560    HKL = G2lat.sortHKLd(HKL,True,False)
561    wave = data['wavelength']
562#set up 1st ring
563    elcent,phi,radii = ellipse              #from fit of 1st ring
564    dsp = HKL[0][3]
565    tth = 2.0*asind(wave/(2.*dsp))
566    ttth = nptand(tth)
567    stth = npsind(tth)
568    ctth = npcosd(tth)
569#1st estimate of tilt; assume ellipse - don't know sign though
570    tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
571    if not tilt:
572        print 'WARNING - selected ring was fitted as a circle'
573        print ' - if detector was tilted we suggest you skip this ring - WARNING'
574#1st estimate of dist: sample to detector normal to plane
575    data['distance'] = dist = radii[0]**2/(ttth*radii[1])
576#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
577    zdisp = radii[1]*ttth*tand(tilt)
578    zdism = radii[1]*ttth*tand(-tilt)
579#cone axis position; 2 choices. Which is right?
580    centp = [elcent[0]+zdisp*sind(phi),elcent[1]+zdisp*cosd(phi)]
581    centm = [elcent[0]+zdism*sind(phi),elcent[1]+zdism*cosd(phi)]
582#check get same ellipse parms either way
583#now do next ring; estimate either way & check sum of Imax/Imin in 3x3 block around each point
584    dsp = HKL[1][3]
585    tth = 2.0*asind(wave/(2.*dsp))
586    ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
587    Ringp,delt,sumIp = makeRing(dsp,ellipsep,2,cutoff,scalex,scaley,self.ImageZ)
588    if len(Ringp) > 10:
589        outEp,errp = FitRing(Ringp,True)
590    else:
591        errp = 1e6
592#    print '+',ellipsep,errp
593    ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
594    Ringm,delt,sumIm = makeRing(dsp,ellipsem,2,cutoff,scalex,scaley,self.ImageZ)
595    if len(Ringm) > 10:
596        outEm,errm = FitRing(Ringm,True)
597    else:
598        errm = 1e6
599#    print '-',ellipsem,errm
600    if errp < errm:
601        data['tilt'] = tilt
602        data['center'] = centp
603        negTilt = 1
604    else:
605        data['tilt'] = -tilt
606        data['center'] = centm
607        negTilt = -1
608    data['rotation'] = phi
609    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
610        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
611    varyList = ['dist','det-X','det-Y','tilt','phi']
612    if data['DetDepthRef']:
613        varyList.append('dep')
614    for i,H in enumerate(HKL):
615        dsp = H[3]
616        tth = 2.0*asind(wave/(2.*dsp))
617        if tth+abs(data['tilt']) > 90.:
618            print 'next line is a hyperbola - search stopped'
619            break
620        print 'HKLD:',H[:4],'2-theta: %.4f'%(tth)
621        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
622        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
623        print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1])
624        Ring,delt,sumI = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
625        if Ring:
626            data['rings'].append(np.array(Ring))
627            rings = np.concatenate((data['rings']),axis=0)
628            if i:
629                chi,chisq = FitDetector(rings,varyList,parmDict,False)
630                data['distance'] = parmDict['dist']
631                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
632                data['rotation'] = np.mod(parmDict['phi'],360.0)
633                data['tilt'] = parmDict['tilt']
634                data['DetDepth'] = parmDict['dep']
635                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
636                print fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings))
637            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
638        else:
639            print 'insufficient number of points in this ellipse to fit'
640            break
641    G2plt.PlotImage(self,newImage=True)
642    fullSize = len(self.ImageZ)/scalex
643    if 2*radii[1] < .9*fullSize:
644        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
645    N = len(data['ellipses'])
646    if N > 2:
647        FitDetector(rings,varyList,parmDict)
648        data['distance'] = parmDict['dist']
649        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
650        data['rotation'] = np.mod(parmDict['phi'],360.0)
651        data['tilt'] = parmDict['tilt']
652        data['DetDepth'] = parmDict['dep']
653    for H in HKL[:N]:
654        ellipse = GetEllipse(H[3],data)
655        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
656    print 'calibration time = ',time.time()-time0
657    G2plt.PlotImage(self,newImage=True)       
658    return True
659   
660def Make2ThetaAzimuthMap(data,masks,iLim,jLim):
661    'Needs a doc string'
662    import numpy.ma as ma
663    import polymask as pm
664    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
665    pixelSize = data['pixelSize']
666    scalex = pixelSize[0]/1000.
667    scaley = pixelSize[1]/1000.
668   
669    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
670    tax = np.asfarray(tax*scalex,dtype=np.float32)
671    tay = np.asfarray(tay*scaley,dtype=np.float32)
672    nI = iLim[1]-iLim[0]
673    nJ = jLim[1]-jLim[0]
674    #make position masks here
675    frame = masks['Frames']
676    tam = ma.make_mask_none((nI,nJ))
677    if frame:
678        tamp = ma.make_mask_none((1024*1024))
679        tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
680            tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
681        tam = ma.mask_or(tam.flatten(),tamp)
682    polygons = masks['Polygons']
683    for polygon in polygons:
684        if polygon:
685            tamp = ma.make_mask_none((1024*1024))
686            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
687                tay.flatten(),len(polygon),polygon,tamp)[:nI*nJ])
688            tam = ma.mask_or(tam.flatten(),tamp)
689    if tam.shape: tam = np.reshape(tam,(nI,nJ))
690    spots = masks['Points']
691    for X,Y,diam in spots:
692        tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2))
693        tam = ma.mask_or(tam,tamp)
694    TA = np.array(GetTthAzm(tax,tay,data))
695    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
696    return np.array(TA),tam           #2-theta & azimuth arrays & position mask
697
698def Fill2ThetaAzimuthMap(masks,TA,tam,image):
699    'Needs a doc string'
700    import numpy.ma as ma
701    Zlim = masks['Thresholds'][1]
702    rings = masks['Rings']
703    arcs = masks['Arcs']
704    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0])))    #azimuth, 2-theta
705    tax,tay = np.dsplit(TA,2)    #azimuth, 2-theta
706    for tth,thick in rings:
707        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
708    for tth,azm,thick in arcs:
709        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
710        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
711        tam = ma.mask_or(tam.flatten(),tamt*tama)
712    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
713    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
714    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
715    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
716    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))
717    del(tam)
718    return tax,tay,taz
719   
720def ImageIntegrate(image,data,masks,dlg=None):
721    'Needs a doc string'
722    import histogram2d as h2d
723    print 'Begin image integration'
724    blkSize = 128   #this seems to be optimal; will break in polymask if >1024
725    LUtth = data['IOtth']
726    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
727    numAzms = data['outAzimuths']
728    numChans = data['outChannels']
729    Dtth = (LUtth[1]-LUtth[0])/numChans
730    Dazm = (LRazm[1]-LRazm[0])/numAzms
731    if data['centerAzm']:
732        LRazm += Dazm/2.
733    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
734    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
735    imageN = len(image)
736    Nx,Ny = data['size']
737    nXBlks = (Nx-1)/blkSize+1
738    nYBlks = (Ny-1)/blkSize+1
739    Nup = nXBlks*nYBlks*3+3
740    t0 = time.time()
741    Nup = 0
742    if dlg:
743        dlg.Update(Nup)
744    for iBlk in range(nYBlks):
745        iBeg = iBlk*blkSize
746        iFin = min(iBeg+blkSize,Ny)
747        for jBlk in range(nXBlks):
748            jBeg = jBlk*blkSize
749            jFin = min(jBeg+blkSize,Nx)               
750#            print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin
751            TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
752           
753            Nup += 1
754            if dlg:
755                dlg.Update(Nup)
756            Block = image[iBeg:iFin,jBeg:jFin]
757            tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
758            del TA,tam
759            Nup += 1
760            if dlg:
761                dlg.Update(Nup)
762            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
763            tax = np.where(tax < LRazm[0],tax+360.,tax)
764            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
765                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
766#            print 'block done'
767            del tax,tay,taz
768            Nup += 1
769            if dlg:
770                dlg.Update(Nup)
771    NST = np.array(NST)
772    H0 = np.divide(H0,NST)
773    H0 = np.nan_to_num(H0)
774    del NST
775    if Dtth:
776        H2 = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
777    else:
778        H2 = np.array(LUtth)
779    if Dazm:       
780        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
781    else:
782        H1 = LRazm
783    H0 /= npcosd(H2[:-1])**2
784    if data['Oblique'][1]:
785        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
786    if 'SASD' in data['type'] and data['PolaVal'][1]:
787        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm)[0] for azm in H1[:-1]+np.diff(H1)/2.])
788    Nup += 1
789    if dlg:
790        dlg.Update(Nup)
791    t1 = time.time()
792    print 'Integration complete'
793    print "Elapsed time:","%8.3f"%(t1-t0), "s"
794    return H0,H1,H2
795   
796def FitStrSta(Image,StrSta,Controls,Masks):
797    'Needs a doc string'
798   
799#    print 'Masks:',Masks
800    StaControls = copy.deepcopy(Controls)
801    phi = StrSta['Sample phi']
802    wave = Controls['wavelength']
803    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
804    pixSize = StaControls['pixelSize']
805    scalex = 1000./pixSize[0]
806    scaley = 1000./pixSize[1]
807    rings = []
808
809    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
810        ellipse = GetEllipse(ring['Dset'],StaControls)
811        Ring,delt,sumI = makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)
812        Ring = np.array(Ring).T
813        ring['ImxyObs'] = np.array(Ring[:2])      #need to apply masks to this to eliminate bad points
814        Ring[:2] = GetTthAzm(Ring[0],Ring[1],StaControls)       #convert x,y to tth,azm
815        Ring[0] /= 2.                                           #convert to theta
816        if len(rings):
817            rings = np.concatenate((rings,Ring),axis=1)
818        else:
819            rings = np.array(Ring)
820    E = StrSta['strain']
821    p0 = [E[0][0],E[1][1],E[2][2],E[0][1],E[0][2],E[1][2]]
822    E = FitStrain(rings,p0,wave,phi)
823    StrSta['strain'] = E
824
825def calcFij(omg,phi,azm,th):
826    '''Does something...
827
828    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
829
830    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface, 90 when perp. to sample surface
831    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
832    :param azm: his chi = azimuth around incident beam
833    :param th:  his theta = theta
834    '''
835    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
836    b = -npcosd(azm)*npcosd(th)
837    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
838    d = a*npcosd(phi)-b*npsind(phi)
839    e = a*npsind(phi)+b*npcosd(phi)
840    Fij = np.array([
841        [d**2,d*e,c*d],
842        [d*e,e**2,c*e],
843        [c*d,c*e,c**2]])
844    return -Fij*nptand(th)
845
846def FitStrain(rings,p0,wave,phi):
847    'Needs a doc string'
848    def StrainPrint(ValSig):
849        print 'Strain tensor:'
850        ptlbls = 'names :'
851        ptstr =  'values:'
852        sigstr = 'esds  :'
853        for name,fmt,value,sig in ValSig:
854            ptlbls += "%s" % (name.rjust(12))
855            ptstr += fmt % (value)
856            if sig:
857                sigstr += fmt % (sig)
858            else:
859                sigstr += 12*' '
860        print ptlbls
861        print ptstr
862        print sigstr
863       
864    def strainCalc(p,xyd,wave,phi):
865#        E = np.array([[p[0],p[3],p[4]],[p[3],p[1],p[5]],[p[4],p[5],p[2]]])
866        E = np.array([[p[0],0,0],[0,p[1],0],[0,0,0]])
867        th,azm,dsp = xyd
868        th0 = npasind(wave/(2.*dsp))
869        dth = 180.*np.sum(E*calcFij(phi,0.,azm,th).T)/(np.pi*1.e6) #in degrees & microstrain units
870        th0 += dth
871        return (th-th0)**2
872       
873    names = ['e11','e22','e33','e12','e13','e23']
874    fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.5f']
875    p1 = [p0[0],p0[1]]   
876    result = leastsq(strainCalc,p1,args=(rings,wave,phi),full_output=True)
877    vals = list(result[0])
878    chi = np.sqrt(np.sum(strainCalc(result[0],rings,wave,phi)**2))
879    sig = list(chi*np.sqrt(np.diag(result[1])))
880    ValSig = zip(names,fmt,vals,sig)
881    StrainPrint(ValSig)
882#    return np.array([[vals[0],vals[3],vals[4]],[vals[3],vals[1],vals[5]],[vals[4],vals[5],vals[2]]])
883    return np.array([[vals[0],0,0],[0,vals[1],0],[0,0,0]])
884   
885       
Note: See TracBrowser for help on using the repository browser.