source: trunk/GSASIIimage.py @ 893

Last change on this file since 893 was 893, checked in by vondreele, 10 years ago

final check on 2D penetration correction

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