# source:trunk/GSASIIimage.py@114

Last change on this file since 114 was 114, checked in by vondreel, 13 years ago

new integration by blocks

File size: 19.5 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)
123        stth = npsind(tth)
124        cosb = npcosd(tilt)
125        R1 = dist*stth*npcosd(tth)*cosb/(cosb**2-stth**2)
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):
153    cphi = cosd(phi)
154    sphi = sind(phi)
155    ring = []
156    for a in range(-180,180,2):
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):
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]):
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
192    stth = sind(tth)
193    ctth = cosd(tth)
194    ttth = tand(tth)
196
199    if cosB > 1.:
200        return 0.,1.
201    else:
202        cosb = math.sqrt(1.-sinb**2)
203        ttth = tand(tth)
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']
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)
223        elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)]
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
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)
248    stth = sind(tth)
249    cosb = cosd(tilt)
250    R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2)
252    zdis = R1*ttth*tand(tilt)
253    A = zdis*sind(phi)
254    B = -zdis*cosd(phi)
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']
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)
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)
370        zsinp = zdis*sind(phi)
371        zcosp = zdis*cosd(phi)
372        cent = data['center']
373        elcent = [cent[0]+zsinp,cent[1]-zcosp]
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)
381            if abs(phi) > 45. and phi < 0.:
382                phi += 180.
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:
391                    print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i)
392                    return False
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
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
462    import numpy.ma as ma
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
469    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
470    tax = np.asfarray(tax*scalex,dtype=np.float32)
471    tay = np.asfarray(tay*scaley,dtype=np.float32)
472    nI = iLim[1]-iLim[0]
473    nJ = jLim[1]-jLim[0]
478    for polygon in polygons:
481            tax.flatten(),tay.flatten(),len(polygon),polygon,tamp)))
482    if tam.shape: tam = np.reshape(tam,(nI,nJ))
483    for X,Y,diam in spots:
485    return GetTthAzm(tax,tay,data),tam           #2-theta & azimuth arrays & position mask
486
488    import numpy.ma as ma
492    nJ = len(image)
493    nI = len(image[0])
494    TA = np.reshape(TA,(2,nI,nJ))
495    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0])))    #azimuth, 2-theta
496    tax,tay = np.dsplit(TA,2)    #azimuth, 2-theta
497    for tth,thick in rings:
499    for tth,azm,thick in arcs:
507    del(tam)
508    return tax,tay,taz
509
511    import histogram2d as h2d
512    print 'Begin image integration'
513    LUtth = data['IOtth']
514    if data['fullIntegrate']:
515        LRazm = [-180,180]
516    else:
517        LRazm = data['LRazimuth']
518    numAzms = data['outAzimuths']
519    numChans = data['outChannels']
520    Dtth = (LUtth[1]-LUtth[0])/numChans
521    Dazm = (LRazm[1]-LRazm[0])/numAzms
522    NST = np.zeros(shape=(numAzms,numChans),dtype=np.int,order='F')
523    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
524    imageN = len(self.ImageZ)
525    nBlks = (imageN-1)/1024+1
526    dlg = wx.ProgressDialog("Elapsed time","2D image integration",nBlks*nBlks*3+3,
527        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
528    try:
529        t0 = time.time()
530        Nup = 0
531        dlg.Update(Nup)
532        for iBlk in range(nBlks):
533            iBeg = iBlk*1024
534            iFin = min(iBeg+1024,imageN)
535            for jBlk in range(nBlks):
536                jBeg = jBlk*1024
537                jFin = min(jBeg+1024,imageN)
538                print 'Process map block',iBlk,jBlk,iBeg,iFin,jBeg,jFin
539                TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
540                Nup += 1
541                dlg.Update(Nup)
542                Block = self.ImageZ[iBeg:iFin,jBeg:jFin]
544                del TA
545                Nup += 1
546                dlg.Update(Nup)
547                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
548                Nup += 1
549                dlg.Update(Nup)
550        H0 = np.nan_to_num(np.divide(H0,NST))
551        del NST
552        if Dtth:
553            H2 = [tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)]
554        else:
555            H2 = LUtth
556        if Dazm:
557            H1 = [azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)]
558        else:
559            H1 = LRazm
560        self.Integrate = [H0,H1,H2]
561        print 'Binning complete'
562        Nup += 1
563        dlg.Update(Nup)
564        t1 = time.time()
565        print "Elapsed time:","%8.3f"%(t1-t0), "s"
566        print 'Integration complete'
567    finally:
568        dlg.Destroy()
Note: See TracBrowser for help on using the repository browser.