source: trunk/GSASIIimage.py @ 98

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

new fortran polymask used

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