source: trunk/GSASIIimage.py @ 81

Last change on this file since 81 was 81, checked in by vondreel, 12 years ago

split out image stuff into GSASIIimage.py
File size: 16.8 KB
Line 
1#GSASII image calculations: ellipse fitting & image integration       
2import math
3import wx
4import time
5import numpy as np
6import numpy.linalg as nl
7import GSASIIplot as G2plt
8import GSASIIlattice as G2lat
9
10# trig functions in degrees
11sind = lambda x: math.sin(x*math.pi/180.)
12asind = lambda x: 180.*math.asin(x)/math.pi
13tand = lambda x: math.tan(x*math.pi/180.)
14atand = lambda x: 180.*math.atan(x)/math.pi
15atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
16cosd = lambda x: math.cos(x*math.pi/180.)
17acosd = lambda x: 180.*math.acos(x)/math.pi
18rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
19#numpy versions
20npsind = lambda x: np.sin(x*np.pi/180.)
21npasind = lambda x: 180.*np.arcsin(x)/np.pi
22npcosd = lambda x: np.cos(x*np.pi/180.)
23nptand = lambda x: np.tan(x*np.pi/180.)
24npatand = lambda x: 180.*np.arctan(x)/np.pi
25npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
26       
27def makeMat(Angle,Axis):
28    #Make rotation matrix from Angle in degrees,Axis =0 for rotation about x, =1 for about y, etc.
29    cs = cosd(Angle)
30    ss = sind(Angle)
31    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
32    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
33                   
34def FitRing(ring):
35    Err,parms = FitCircle(ring)
36    Err /= len(ring)
37#    print 'circle error:','%8f'%(Err)
38    if Err > 20000.:
39        eparms = FitEllipse(ring)
40        if eparms:
41            parms = eparms
42    return parms
43       
44def FitCircle(ring):
45   
46    def makeParmsCircle(B):
47        cent = [-B[0]/2,-B[1]/2]
48        phi = 0.
49        sr1 = sr2 = math.sqrt(cent[0]**2+cent[1]**2-B[2])
50        return cent,phi,[sr1,sr2]
51       
52    ring = np.array(ring)
53    x = np.asarray(ring.T[0])
54    y = np.asarray(ring.T[1])
55   
56    M = np.array((x,y,np.ones_like(x)))
57    B = np.array(-(x**2+y**2))
58    result = nl.lstsq(M.T,B)
59    return result[1],makeParmsCircle(result[0])
60       
61def FitEllipse(ring):
62           
63    def makeParmsEllipse(B):
64        det = 4.*(1.-B[0]**2)-B[1]**2
65        if det < 0.:
66            print 'hyperbola!'
67            return 0
68        elif det == 0.:
69            print 'parabola!'
70            return 0
71        cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \
72            (B[1]*B[2]-2.*(1.+B[0])*B[3])/det]
73        phi = 0.5*atand(0.5*B[1]/B[0])
74       
75        a = (1.+B[0])/cosd(2*phi)
76        b = 2.-a
77        f = (1.+B[0])*cent[0]**2+(1.-B[0])*cent[1]**2+B[1]*cent[0]*cent[1]-B[4]
78        if f/a < 0 or f/b < 0:
79            return 0
80        sr1 = math.sqrt(f/a)
81        sr2 = math.sqrt(f/b)
82        if sr1 > sr2:
83            sr1,sr2 = sr2,sr1
84            phi -= 90.
85            if phi < -90.:
86                phi += 180.
87        return cent,phi,[sr1,sr2]
88               
89    ring = np.array(ring)
90    x = np.asarray(ring.T[0])
91    y = np.asarray(ring.T[1])
92    M = np.array((x**2-y**2,x*y,x,y,np.ones_like(x)))
93    B = np.array(-(x**2+y**2))
94    result = nl.lstsq(M.T,B)
95    return makeParmsEllipse(result[0])
96   
97def FitDetector(rings,p0,wave):
98    from scipy.optimize import leastsq
99    def ellipseCalc(B,xyd,wave):
100        x = xyd[0]
101        y = xyd[1]
102        dsp = xyd[2]
103        dist,x0,y0,phi,tilt = B
104        tth = 2.0*npasind(wave/(2.*dsp))
105        ttth = nptand(tth)
106        radius = dist*ttth
107        stth = npsind(tth)
108        cosb = npcosd(tilt)
109        R1 = dist*stth*npcosd(tth)*cosb/(cosb**2-stth**2)
110        R0 = np.sqrt(R1*radius*cosb)
111        zdis = R1*ttth*nptand(tilt)
112        X = x-x0+zdis*npsind(phi)
113        Y = y-y0-zdis*npcosd(phi)
114        XR = X*npcosd(phi)-Y*npsind(phi)
115        YR = X*npsind(phi)+Y*npcosd(phi)
116        return (XR/R0)**2+(YR/R1)**2-1
117    result = leastsq(ellipseCalc,p0,args=(rings.T,wave))
118    return result[0]
119           
120def ImageLocalMax(image,w,Xpix,Ypix):
121    w2 = w*2
122    size = len(image)
123    xpix = int(Xpix)            #get reference corner of pixel chosen
124    ypix = int(Ypix)
125    if (w < xpix < size-w) and (w < ypix < size-w) and image[ypix,xpix]:
126        Z = image[ypix-w:ypix+w,xpix-w:xpix+w]
127        Zmax = np.argmax(Z)
128        Zmin = np.argmin(Z)
129        xpix += Zmax%w2-w
130        ypix += Zmax/w2-w
131        return xpix,ypix,np.ravel(Z)[Zmax],np.ravel(Z)[Zmin]
132    else:
133        return 0,0,0,0
134   
135def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
136    cent,phi,radii = ellipse
137    cphi = cosd(phi)
138    sphi = sind(phi)
139    ring = []
140    for a in range(-180,180,2):
141        x = radii[0]*cosd(a)
142        y = radii[1]*sind(a)
143        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
144        Y = (sphi*x+cphi*y+cent[1])*scaley
145        X,Y,I,J = ImageLocalMax(image,pix,X,Y)     
146        if I and J and I/J > reject:
147            X += .5                             #set to center of pixel
148            Y += .5
149            X /= scalex                         #convert to mm
150            Y /= scaley
151            ring.append([X,Y,dsp])
152    if len(ring) < 45:             #want more than 1/4 of a circle
153        return []
154    return ring
155   
156def makeIdealRing(ellipse):
157    cent,phi,radii = ellipse
158    cphi = cosd(phi)
159    sphi = sind(phi)
160    ring = []
161    for a in range(0,360,2):
162        x = radii[0]*cosd(a)
163        y = radii[1]*sind(a)
164        X = (cphi*x-sphi*y+cent[0])
165        Y = (sphi*x+cphi*y+cent[1])
166        ring.append([X,Y])
167    return ring
168   
169def calcDist(radii,tth):
170    stth = sind(tth)
171    ctth = cosd(tth)
172    ttth = tand(tth)
173    return math.sqrt(radii[0]**4/(ttth**2*((radii[0]*ctth)**2+(radii[1]*stth)**2)))
174   
175def calcZdisCosB(radius,tth,radii):
176    cosB = sinb = radii[0]**2/(radius*radii[1])
177    if cosB > 1.:
178        return 0.,1.
179    else:
180        cosb = math.sqrt(1.-sinb**2)
181        ttth = tand(tth)
182        zdis = radii[1]*ttth*cosb/sinb
183        return zdis,cosB
184   
185def GetEllipse(dsp,data):
186    dist = data['distance']
187    cent = data['center']
188    tilt = data['tilt']
189    phi = data['rotation']
190    radii = [0,0]
191    tth = 2.0*asind(data['wavelength']/(2.*dsp))
192    ttth = tand(tth)
193    stth = sind(tth)
194    ctth = cosd(tth)
195    cosb = cosd(tilt)
196    radius = dist*ttth
197    radii[1] = dist*stth*ctth*cosb/(cosb**2-stth**2)
198    if radii[1] > 0:
199        radii[0] = math.sqrt(radii[1]*radius*cosb)
200        zdis = radii[1]*ttth*tand(tilt)
201        elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)]
202        return elcent,phi,radii
203    else:
204        return False
205       
206def GetDetectorXY(dsp,azm,data):
207    from scipy.optimize import fsolve
208    def func(xy,*args):
209       azm,phi,R0,R1,A,B = args
210       cp = cosd(phi)
211       sp = sind(phi)
212       x,y = xy
213       out = []
214       out.append(y-x*tand(azm))
215       out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)
216       return out
217    elcent,phi,radii = GetEllipse(dsp,data)
218    cent = data['center']
219    tilt = data['tilt']
220    phi = data['rotation']
221    wave = data['wavelength']
222    dist = data['distance']
223    tth = 2.0*asind(wave/(2.*dsp))
224    ttth = tand(tth)
225    radius = dist*ttth
226    stth = sind(tth)
227    cosb = cosd(tilt)
228    R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2)
229    R0 = math.sqrt(R1*radius*cosb)
230    zdis = R1*ttth*tand(tilt)
231    A = zdis*sind(phi)
232    B = -zdis*cosd(phi)
233    xy0 = [radius*cosd(azm),radius*sind(azm)]
234    xy = fsolve(func,xy0,args=(azm,phi,R0,R1,A,B))+cent
235    return xy
236                   
237def GetTthAzmDsp(x,y,data):
238    wave = data['wavelength']
239    dist = data['distance']
240    cent = data['center']
241    tilt = data['tilt']
242    phi = data['rotation']
243    dx = np.array(x-cent[0],dtype=np.float32)
244    dy = np.array(y-cent[1],dtype=np.float32)
245    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
246    X = np.dot(X,makeMat(phi,2))
247    Z = np.dot(X,makeMat(tilt,0)).T[2]
248    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
249    dsp = wave/(2.*npsind(tth/2.))
250    azm = npatan2d(dy,dx)
251    return tth,azm,dsp
252   
253def GetTth(x,y,data):
254    return GetTthAzmDsp(x,y,data)[0]
255   
256def GetTthAzm(x,y,data):
257    return GetTthAzmDsp(x,y,data)[0:2]
258   
259def GetDsp(x,y,data):
260    return GetTthAzmDsp(x,y,data)[2]
261       
262def ImageCompress(image,scale):
263    if scale == 1:
264        return image
265    else:
266        return image[::scale,::scale]
267       
268def ImageCalibrate(self,data):
269    import copy
270    import ImageCalibrants as calFile
271    print 'image calibrate'
272    ring = data['ring']
273    pixelSize = data['pixelSize']
274    scalex = 1000./pixelSize[0]
275    scaley = 1000./pixelSize[1]
276    cutoff = data['cutoff']
277    if len(ring) < 5:
278        print 'not enough inner ring points for ellipse'
279        return False
280       
281    #fit start points on inner ring
282    data['ellipses'] = []
283    outE = FitRing(ring)
284    if outE:
285        print 'start ellipse:',outE
286        ellipse = outE
287    else:
288        return False
289       
290    #setup 180 points on that ring for "good" fit
291    Ring = makeRing(1.0,ellipse,20,cutoff,scalex,scaley,self.ImageZ)
292    if Ring:
293        ellipse = FitRing(Ring)
294        Ring = makeRing(1.0,ellipse,20,cutoff,scalex,scaley,self.ImageZ)    #do again
295        ellipse = FitRing(Ring)
296    else:
297        print '1st ring not sufficiently complete to proceed'
298        return False
299    print 'inner ring:',ellipse
300    data['center'] = copy.copy(ellipse[0])           #not right!! (but useful for now)
301    data['ellipses'].append(ellipse[:]+('r',))
302    G2plt.PlotImage(self)
303   
304    #setup for calibration
305    data['rings'] = []
306    data['ellipses'] = []
307    if not data['calibrant']:
308        print 'no calibration material selected'
309        return True
310       
311    Bravais,cell = calFile.Calibrants[data['calibrant']]
312    A = G2lat.cell2A(cell)
313    wave = data['wavelength']
314    cent = data['center']
315    pixLimit = data['pixLimit']
316    elcent,phi,radii = ellipse
317    HKL = G2lat.GenHBravais(0.5,Bravais,A)
318    dsp = HKL[0][3]
319    tth = 2.0*asind(wave/(2.*dsp))
320    ttth = tand(tth)
321    data['distance'] = dist = calcDist(radii,tth)
322    radius = dist*tand(tth)
323    zdis,cosB = calcZdisCosB(radius,tth,radii)
324    cent1 = []
325    cent2 = []
326    xSum = 0
327    ySum = 0
328    zxSum = 0
329    zySum = 0
330    phiSum = 0
331    tiltSum = 0
332    distSum = 0
333    Zsum = 0
334    for i,H in enumerate(HKL):
335        dsp = H[3]
336        tth = 2.0*asind(0.5*wave/dsp)
337        stth = sind(tth)
338        ctth = cosd(tth)
339        ttth = tand(tth)
340        radius = dist*ttth
341        elcent,phi,radii = ellipse
342        radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2)
343        radii[0] = math.sqrt(radii[1]*radius*cosB)
344        zdis,cosB = calcZdisCosB(radius,tth,radii)
345        zsinp = zdis*sind(phi)
346        zcosp = zdis*cosd(phi)
347        cent = data['center']
348        elcent = [cent[0]+zsinp,cent[1]-zcosp]
349        ratio = radii[1]/radii[0]
350        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
351        if Ring:
352            numZ = len(Ring)
353            data['rings'].append(np.array(Ring))
354            ellipse = FitRing(Ring)
355            elcent,phi,radii = ellipse               
356            if abs(phi) > 45. and phi < 0.:
357                phi += 180.
358            dist = calcDist(radii,tth)
359            distR = 1.-dist/data['distance']
360            if distR > 0.001:
361                print 'Wavelength too large?'
362            elif distR < -0.001:
363                print 'Wavelength too small?'
364            else:
365                if abs((radii[1]/radii[0]-ratio)/ratio) > 0.01:
366                    print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i)
367                    return False
368            zdis,cosB = calcZdisCosB(radius,tth,radii)
369            Tilt = acosd(cosB)          # 0 <= tilt <= 90
370            zsinp = zdis*sind(ellipse[1])
371            zcosp = zdis*cosd(ellipse[1])
372            cent1.append(np.array([elcent[0]+zsinp,elcent[1]-zcosp]))
373            cent2.append(np.array([elcent[0]-zsinp,elcent[1]+zcosp]))
374            if i:
375                d1 = cent1[-1]-cent1[-2]        #get shift of 2 possible center solutions
376                d2 = cent2[-1]-cent2[-2]
377                if np.dot(d2,d2) > np.dot(d1,d1):  #right solution is the larger shift
378                    data['center'] = cent1[-1]
379                else:
380                    data['center'] = cent2[-1]
381                Zsum += numZ
382                phiSum += numZ*phi
383                distSum += numZ*dist
384                xSum += numZ*data['center'][0]
385                ySum += numZ*data['center'][1]
386                tiltSum += numZ*abs(Tilt)
387            cent = data['center']
388            print ('for ring # %2i dist %.3f rotate %6.2f tilt %6.2f Xcent %.3f Ycent %.3f Npts %d' 
389                %(i,dist,phi,Tilt,cent[0],cent[1],numZ))
390            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
391            G2plt.PlotImage(self)
392        else:
393            break
394    fullSize = len(self.ImageZ)/scalex
395    if 2*radii[1] < .9*fullSize:
396        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
397    if not Zsum:
398        print 'Only one ring fitted. Check your wavelength.'
399        return False
400    cent = data['center'] = [xSum/Zsum,ySum/Zsum]
401    wave = data['wavelength']
402    dist = data['distance'] = distSum/Zsum
403   
404    #possible error if no. of rings < 3! Might need trap here
405    d1 = cent1[-1]-cent1[1]             #compare last ring to 2nd ring
406    d2 = cent2[-1]-cent2[1]
407    Zsign = 1
408    len1 = math.sqrt(np.dot(d1,d1))
409    len2 = math.sqrt(np.dot(d2,d2))
410    t1 = d1/len1
411    t2 = d2/len2
412    if len2 > len1:
413        if -135. < atan2d(t2[1],t2[0]) < 45.:
414            Zsign = -1
415    else:
416        if -135. < atan2d(t1[1],t1[0]) < 45.:
417            Zsign = -1
418   
419    tilt = data['tilt'] = Zsign*tiltSum/Zsum
420    phi = data['rotation'] = phiSum/Zsum
421    rings = np.concatenate((data['rings']),axis=0)
422    p0 = [dist,cent[0],cent[1],phi,tilt]
423    result = FitDetector(rings,p0,wave)
424    data['distance'] = result[0]
425    data['center'] = result[1:3]
426    data['rotation'] = np.mod(result[3],360.0)
427    data['tilt'] = result[4]
428    N = len(data['ellipses'])
429    data['ellipses'] = []           #clear away individual ellipse fits
430    for H in HKL[:N]:
431        ellipse = GetEllipse(H[3],data)
432        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
433    G2plt.PlotImage(self)       
434    return True
435   
436def Make2ThetaAzimuthMap(data,imageN):
437    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
438    pixelSize = data['pixelSize']
439    scalex = pixelSize[0]/1000.
440    scaley = pixelSize[1]/1000.
441    tax,tay = np.mgrid[0.5:imageN+.5,0.5:imageN+.5]         #bin centers not corners
442    tax = np.asfarray(tax*scalex,dtype=np.float32)
443    tay = np.asfarray(tay*scaley,dtype=np.float32)
444    return GetTthAzm(tay,tax,data)           #2-theta & azimuth arrays
445
446def Fill2ThetaAzimuthMap(data,TA,image):
447    import numpy.ma as ma
448    LUtth = data['IOtth']
449    if data['fullIntegrate']:
450        LRazm = [-180,180]
451    else:
452        LRazm = data['LRazimuth']
453    imageN = len(image)
454    TA = np.reshape(TA,(2,imageN,imageN))
455    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0])))    #azimuth, 2-theta
456    tax,tay = np.dsplit(TA,2)    #azimuth, 2-theta
457    tax = ma.masked_outside(tax.flatten(),LRazm[0],LRazm[1])
458    tay = ma.masked_outside(tay.flatten(),LUtth[0],LUtth[1])
459    tam = ma.getmask(tax)+ma.getmask(tay)
460    taz = ma.masked_where(tam,image.flatten())
461    return tax,tay,taz,tam
462   
463def Bin2ThetaAzimuthMap(data,tax,tay,taz):
464    import numpy.ma as ma
465    LUtth = data['IOtth']
466    if data['fullIntegrate']:
467        LRazm = [-180,180]
468    else:
469        LRazm = data['LRazimuth']
470    numAzms = data['outAzimuths']
471    numChans = data['outChannels']
472    NST = np.histogram2d(tax,tay,normed=False,bins=(numAzms,numChans),range=[LRazm,LUtth])
473    HST = np.histogram2d(tax,tay,normed=False,bins=(numAzms,numChans),range=[LRazm,LUtth],weights=taz)
474    return NST,HST
475
476def ImageIntegrate(self,data):
477    dlg = wx.ProgressDialog("Elapsed time","2D image integration",5,
478        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
479    try:
480        print 'Begin image integration'
481        print 'Create 2-theta,azimuth map'
482        t0 = time.time()
483        dlg.Update(0)
484        imageN = len(self.ImageZ)
485        TA = Make2ThetaAzimuthMap(data,imageN)           #2-theta & azimuth arrays
486        dlg.Update(1)
487        print 'Fill map with 2-theta/azimuth values'
488        tax,tay,taz,tam = Fill2ThetaAzimuthMap(data,TA,self.ImageZ)
489        del TA
490        dlg.Update(2)
491        print 'Bin image by 2-theta/azimuth intervals'
492        NST,HST = Bin2ThetaAzimuthMap(data,tax,tay,taz)
493        del tax,tay,taz
494        dlg.Update(3)
495        print 'Form normalized 1D pattern(s)'
496        self.Integrate = [HST[0]/NST[0],HST[1],HST[2]]
497        del NST,HST
498        dlg.Update(4)
499        t1 = time.time()
500        print 'Integration complete'
501        print "Elapsed time:","%8.3f"%(t1-t0), "s"
502    finally:
503        dlg.Destroy()
Note: See TracBrowser for help on using the repository browser.