# source:trunk/GSASIIimage.py@265

Last change on this file since 265 was 265, checked in by vondreele, 12 years ago

further progress on implementing pdf calculations
optionally put legends on the pdf plots
attempt implementation of a rotation of the azimuth ranges for multiazimuth integrations -not fully successful

• Property svn:keywords set to `Date Author Revision URL Id`
File size: 21.4 KB
Line
1#GSASII image calculations: ellipse fitting & image integration
2########### SVN repository information ###################
3# \$Date: 2011-04-19 20:12:37 +0000 (Tue, 19 Apr 2011) \$
4# \$Author: vondreele \$
5# \$Revision: 265 \$
6# \$URL: trunk/GSASIIimage.py \$
7# \$Id: GSASIIimage.py 265 2011-04-19 20:12:37Z vondreele \$
8########### SVN repository information ###################
9import math
10import wx
11import time
12import numpy as np
13import numpy.linalg as nl
14import GSASIIpath
15import GSASIIplot as G2plt
16import GSASIIlattice as G2lat
17import fellipse as fel
18
19# trig functions in degrees
20sind = lambda x: math.sin(x*math.pi/180.)
21asind = lambda x: 180.*math.asin(x)/math.pi
22tand = lambda x: math.tan(x*math.pi/180.)
23atand = lambda x: 180.*math.atan(x)/math.pi
24atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
25cosd = lambda x: math.cos(x*math.pi/180.)
26acosd = lambda x: 180.*math.acos(x)/math.pi
27rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
28#numpy versions
29npsind = lambda x: np.sin(x*np.pi/180.)
30npasind = lambda x: 180.*np.arcsin(x)/np.pi
31npcosd = lambda x: np.cos(x*np.pi/180.)
32nptand = lambda x: np.tan(x*np.pi/180.)
33npatand = lambda x: 180.*np.arctan(x)/np.pi
34npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
35
36def pointInPolygon(pXY,xy):
37    #pXY - assumed closed 1st & last points are duplicates
38    Inside = False
39    N = len(pXY)
40    p1x,p1y = pXY[0]
41    for i in range(N+1):
42        p2x,p2y = pXY[i%N]
43        if (max(p1y,p2y) >= xy[1] > min(p1y,p2y)) and (xy[0] <= max(p1x,p2x)):
44            if p1y != p2y:
45                xinters = (xy[1]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
46            if p1x == p2x or xy[0] <= xinters:
47                Inside = not Inside
48        p1x,p1y = p2x,p2y
49    return Inside
50
51def makeMat(Angle,Axis):
52    #Make rotation matrix from Angle in degrees,Axis =0 for rotation about x, =1 for about y, etc.
53    cs = cosd(Angle)
54    ss = sind(Angle)
55    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
56    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
57
58def FitRing(ring,delta):
59    parms = []
60    if delta:
61        err,parms = FitEllipse(ring)
62        errc,parmsc = FitCircle(ring)
63        errc = errc[0]/(len(ring)*parmsc[2][0]**2)
64        if not parms or errc < .1:
65            parms = parmsc
66    else:
67        err,parms = FitCircle(ring)
68    return parms
69
70def FitCircle(ring):
71
72    def makeParmsCircle(B):
73        cent = [-B[0]/2,-B[1]/2]
74        phi = 0.
75        sr1 = sr2 = math.sqrt(cent[0]**2+cent[1]**2-B[2])
76        return cent,phi,[sr1,sr2]
77
78    ring = np.array(ring)
79    x = np.asarray(ring.T[0])
80    y = np.asarray(ring.T[1])
81
82    M = np.array((x,y,np.ones_like(x)))
83    B = np.array(-(x**2+y**2))
84    result = nl.lstsq(M.T,B)
85    return result[1],makeParmsCircle(result[0])
86
87def FitEllipse(ring):
88
89    def makeParmsEllipse(B):
90        det = 4.*(1.-B[0]**2)-B[1]**2
91        if det < 0.:
92            print 'hyperbola!'
93            return 0
94        elif det == 0.:
95            print 'parabola!'
96            return 0
97        cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \
98            (B[1]*B[2]-2.*(1.+B[0])*B[3])/det]
99        phi = 0.5*atand(0.5*B[1]/B[0])
100
101        a = (1.+B[0])/cosd(2*phi)
102        b = 2.-a
103        f = (1.+B[0])*cent[0]**2+(1.-B[0])*cent[1]**2+B[1]*cent[0]*cent[1]-B[4]
104        if f/a < 0 or f/b < 0:
105            return 0
106        sr1 = math.sqrt(f/a)
107        sr2 = math.sqrt(f/b)
108        if sr1 > sr2:
109            sr1,sr2 = sr2,sr1
110            phi -= 90.
111            if phi < -90.:
112                phi += 180.
113        return cent,phi,[sr1,sr2]
114
115    ring = np.array(ring)
116    x = np.asarray(ring.T[0])
117    y = np.asarray(ring.T[1])
118    M = np.array((x**2-y**2,x*y,x,y,np.ones_like(x)))
119    B = np.array(-(x**2+y**2))
120    bb,err = fel.fellipse(len(x),x,y,1.E-7)
121#    print nl.lstsq(M.T,B)[0]
122#    print bb
123    return err,makeParmsEllipse(bb)
124
125def FitDetector(rings,p0,wave):
126    from scipy.optimize import leastsq
127
128    def ellipseCalcD(B,xyd,wave):
129        x = xyd[0]
130        y = xyd[1]
131        dsp = xyd[2]
132        dist,x0,y0,phi,tilt = B
133        tth = 2.0*npasind(wave/(2.*dsp))
134        ttth = nptand(tth)
136        stth = npsind(tth)
137        cosb = npcosd(tilt)
138        R1 = dist*stth*npcosd(tth)*cosb/(cosb**2-stth**2)
140        zdis = R1*ttth*nptand(tilt)
141        X = x-x0+zdis*npsind(phi)
142        Y = y-y0-zdis*npcosd(phi)
143        XR = X*npcosd(phi)-Y*npsind(phi)
144        YR = X*npsind(phi)+Y*npcosd(phi)
145        return (XR/R0)**2+(YR/R1)**2-1
146    def ellipseCalcW(C,xyd):
147        dist,x0,y0,phi,tilt,wave = C
148        B = dist,x0,y0,phi,tilt
149        return ellipseCalcD(B,xyd,wave)
150    result = leastsq(ellipseCalcD,p0,args=(rings.T,wave))
151    if len(rings) > 1:
152        p0 = result[0]
153        p0 = np.append(p0,wave)
154        resultW = leastsq(ellipseCalcW,p0,args=(rings.T,))
155        return result[0],resultW[0][-1]
156    else:
157        return result[0],wave
158
159def ImageLocalMax(image,w,Xpix,Ypix):
160    w2 = w*2
161    sizey,sizex = image.shape
162    xpix = int(Xpix)            #get reference corner of pixel chosen
163    ypix = int(Ypix)
164    if (w < xpix < sizex-w) and (w < ypix < sizey-w) and image[ypix,xpix]:
165        Z = image[ypix-w:ypix+w,xpix-w:xpix+w]
166        Zmax = np.argmax(Z)
167        Zmin = np.argmin(Z)
168        xpix += Zmax%w2-w
169        ypix += Zmax/w2-w
170        return xpix,ypix,np.ravel(Z)[Zmax],np.ravel(Z)[Zmin]
171    else:
172        return 0,0,0,0
173
174def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
176    cphi = cosd(phi)
177    sphi = sind(phi)
178    ring = []
179    amin = 180
180    amax = -180
181    for a in range(-180,180,1):
184        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
185        Y = (sphi*x+cphi*y+cent[1])*scaley
186        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
187        if I and J and I/J > reject:
188            X += .5                             #set to center of pixel
189            Y += .5
190            X /= scalex                         #convert to mm
191            Y /= scaley
192            amin = min(amin,a)
193            amax = max(amax,a)
194            ring.append([X,Y,dsp])
195    delt = amax-amin
196    if len(ring) < 20:             #want more than 20 deg
197        return [],delt > 90
198    return ring,delt > 90
199
200def makeIdealRing(ellipse,azm=None):
202    cphi = cosd(phi)
203    sphi = sind(phi)
204    ring = []
205    if azm:
206        aR = azm[0]-90,azm[1]-90,1
207        if azm[1]-azm[0] > 180:
208            aR[2] = 2
209    else:
210        aR = 0,362,2
211    for a in range(aR[0],aR[1],aR[2]):
214        X = (cphi*x-sphi*y+cent[0])
215        Y = (sphi*x+cphi*y+cent[1])
216        ring.append([X,Y])
217    return ring
218
220    stth = sind(tth)
221    ctth = cosd(tth)
222    ttth = tand(tth)
224
227    if cosB > 1.:
228        return 0.,1.
229    else:
230        cosb = math.sqrt(1.-sinb**2)
231        ttth = tand(tth)
233        return zdis,cosB
234
235def GetEllipse(dsp,data):
236    dist = data['distance']
237    cent = data['center']
238    tilt = data['tilt']
239    phi = data['rotation']
241    tth = 2.0*asind(data['wavelength']/(2.*dsp))
242    ttth = tand(tth)
243    stth = sind(tth)
244    ctth = cosd(tth)
245    cosb = cosd(tilt)
251        elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)]
253    else:
254        return False
255
256def GetDetectorXY(dsp,azm,data):
257    from scipy.optimize import fsolve
258    def func(xy,*args):
259       azm,phi,R0,R1,A,B = args
260       cp = cosd(phi)
261       sp = sind(phi)
262       x,y = xy
263       out = []
264       out.append(y-x*tand(azm))
265       out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)
266       return out
268    cent = data['center']
269    tilt = data['tilt']
270    phi = data['rotation']
271    wave = data['wavelength']
272    dist = data['distance']
273    tth = 2.0*asind(wave/(2.*dsp))
274    ttth = tand(tth)
276    stth = sind(tth)
277    cosb = cosd(tilt)
278    R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2)
280    zdis = R1*ttth*tand(tilt)
281    A = zdis*sind(phi)
282    B = -zdis*cosd(phi)
284    xy = fsolve(func,xy0,args=(azm,phi,R0,R1,A,B))+cent
285    return xy
286
287def GetTthAzmDsp(x,y,data):
288    wave = data['wavelength']
289    dist = data['distance']
290    cent = data['center']
291    tilt = data['tilt']
292    phi = data['rotation']
293    LRazim = data['LRazimuth']
294    azmthoff = data['azmthOff']
295    azmRot = data['azmthRotate']
296    Full = data['fullIntegrate']
297    dx = np.array(x-cent[0],dtype=np.float32)
298    dy = np.array(y-cent[1],dtype=np.float32)
299    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
300    X = np.dot(X,makeMat(phi,2))
301    Z = np.dot(X,makeMat(tilt,0)).T[2]
302    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
303    dsp = wave/(2.*npsind(tth/2.))
304    azm = (npatan2d(dx,-dy)+azmthoff+720.)%360.
305    if Full:
306        azm = (azm+azmRot+720.)%360.
307    else:
308        azm = np.where(azm<LRazim[0],azm+360.,azm)
309    return tth,azm,dsp
310
311def GetTth(x,y,data):
312    return GetTthAzmDsp(x,y,data)[0]
313
314def GetTthAzm(x,y,data):
315    return GetTthAzmDsp(x,y,data)[0:2]
316
317def GetDsp(x,y,data):
318    return GetTthAzmDsp(x,y,data)[2]
319
320def GetAzm(x,y,data):
321    return GetTthAzmDsp(x,y,data)[1]
322
323def ImageCompress(image,scale):
324    if scale == 1:
325        return image
326    else:
327        return image[::scale,::scale]
328
329def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
330    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
331    curr = np.array([dist,x,y])
332    return abs(avg-curr)/avg < .02
333
334def ImageCalibrate(self,data):
335    import copy
336    import ImageCalibrants as calFile
337    print 'image calibrate'
338    time0 = time.time()
339    ring = data['ring']
340    pixelSize = data['pixelSize']
341    scalex = 1000./pixelSize[0]
342    scaley = 1000./pixelSize[1]
343    pixLimit = data['pixLimit']
344    cutoff = data['cutoff']
345    if len(ring) < 5:
346        print 'not enough inner ring points for ellipse'
347        return False
348
349    #fit start points on inner ring
350    data['ellipses'] = []
351    outE = FitRing(ring,True)
352    if outE:
353        print 'start ellipse:',outE
354        ellipse = outE
355    else:
356        return False
357
358    #setup 360 points on that ring for "good" fit
359    Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
360    if Ring:
361        ellipse = FitRing(Ring,delt)
362        Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
363        ellipse = FitRing(Ring,delt)
364    else:
365        print '1st ring not sufficiently complete to proceed'
366        return False
367    print 'inner ring:',ellipse
368    data['center'] = copy.copy(ellipse[0])           #not right!! (but useful for now)
369    data['ellipses'].append(ellipse[:]+('r',))
370    G2plt.PlotImage(self,newImage=True)
371
372    #setup for calibration
373    data['rings'] = []
374    data['ellipses'] = []
375    if not data['calibrant']:
376        print 'no calibration material selected'
377        return True
378
379    skip = data['calibskip']
380    dmin = data['calibdmin']
381    Bravais,cell = calFile.Calibrants[data['calibrant']][:2]
382    A = G2lat.cell2A(cell)
383    wave = data['wavelength']
384    cent = data['center']
386    HKL = G2lat.GenHBravais(dmin,Bravais,A)[skip:]
387    dsp = HKL[0][3]
388    tth = 2.0*asind(wave/(2.*dsp))
389    ttth = tand(tth)
390    data['distance'] = dist = calcDist(radii,tth)
393    cent1 = []
394    cent2 = []
395    xSum = 0
396    ySum = 0
397    zxSum = 0
398    zySum = 0
399    phiSum = 0
400    tiltSum = 0
401    distSum = 0
402    Zsum = 0
403    for i,H in enumerate(HKL):
404        dsp = H[3]
405        tth = 2.0*asind(0.5*wave/dsp)
406        stth = sind(tth)
407        ctth = cosd(tth)
408        ttth = tand(tth)
414        zsinp = zdis*sind(phi)
415        zcosp = zdis*cosd(phi)
416        cent = data['center']
417        elcent = [cent[0]+zsinp,cent[1]-zcosp]
419        Ring,delt = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
420        if Ring:
421            numZ = len(Ring)
422            data['rings'].append(np.array(Ring))
423            newellipse = FitRing(Ring,delt)
425            if abs(phi) > 45. and phi < 0.:
426                phi += 180.
428            distR = 1.-dist/data['distance']
429            if abs(distR) > 0.1:
430                print dsp,dist,data['distance'],distR,len(Ring),delt
431                break
432            if distR > 0.001:
433                print 'Wavelength too large?'
434            elif distR < -0.001:
435                print 'Wavelength too small?'
436            else:
437                ellipse = newellipse
439            Tilt = acosd(cosB)          # 0 <= tilt <= 90
440            zsinp = zdis*sind(ellipse[1])
441            zcosp = zdis*cosd(ellipse[1])
442            cent1.append(np.array([elcent[0]+zsinp,elcent[1]-zcosp]))
443            cent2.append(np.array([elcent[0]-zsinp,elcent[1]+zcosp]))
444            if i:
445                d1 = cent1[-1]-cent1[-2]        #get shift of 2 possible center solutions
446                d2 = cent2[-1]-cent2[-2]
447                if np.dot(d2,d2) > np.dot(d1,d1):  #right solution is the larger shift
448                    data['center'] = cent1[-1]
449                else:
450                    data['center'] = cent2[-1]
451                Zsum += numZ
452                phiSum += numZ*phi
453                distSum += numZ*dist
454                xSum += numZ*data['center'][0]
455                ySum += numZ*data['center'][1]
456                tiltSum += numZ*abs(Tilt)
457                if not np.all(checkEllipse(Zsum,distSum,xSum,ySum,dist,data['center'][0],data['center'][1])):
458                    print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i)
459            cent = data['center']
460            print ('for ring # %2i @ d-space %.4f: dist %.3f rotate %6.2f tilt %6.2f Xcent %.3f Ycent %.3f Npts %d'
461                %(i,dsp,dist,phi-90.,Tilt,cent[0],cent[1],numZ))
462            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
463        else:
464            break
465    G2plt.PlotImage(self,newImage=True)
466    fullSize = len(self.ImageZ)/scalex
468        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
469    if not Zsum:
470        print 'Only one ring fitted. Check your wavelength.'
471        return False
472    cent = data['center'] = [xSum/Zsum,ySum/Zsum]
473    wave = data['wavelength']
474    dist = data['distance'] = distSum/Zsum
475
476    #possible error if no. of rings < 3! Might need trap here
477    d1 = cent1[-1]-cent1[1]             #compare last ring to 2nd ring
478    d2 = cent2[-1]-cent2[1]
479    Zsign = 1
480    len1 = math.sqrt(np.dot(d1,d1))
481    len2 = math.sqrt(np.dot(d2,d2))
482    t1 = d1/len1
483    t2 = d2/len2
484    if len2 > len1:
485        if -135. < atan2d(t2[1],t2[0]) < 45.:
486            Zsign = -1
487    else:
488        if -135. < atan2d(t1[1],t1[0]) < 45.:
489            Zsign = -1
490
491    tilt = data['tilt'] = Zsign*tiltSum/Zsum
492    phi = data['rotation'] = phiSum/Zsum
493    rings = np.concatenate((data['rings']),axis=0)
494    p0 = [dist,cent[0],cent[1],phi,tilt]
495    result,newWave = FitDetector(rings,p0,wave)
496    print 'Suggested new wavelength = ',('%.5f')%(newWave),' (not reliable if distance > 2m)'
497    data['distance'] = result[0]
498    data['center'] = result[1:3]
499    data['rotation'] = np.mod(result[3],360.0)
500    data['tilt'] = result[4]
501    N = len(data['ellipses'])
502    data['ellipses'] = []           #clear away individual ellipse fits
503    for H in HKL[:N]:
504        ellipse = GetEllipse(H[3],data)
505        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
506    print 'calibration time = ',time.time()-time0
507    G2plt.PlotImage(self,newImage=True)
508    return True
509
511    import numpy.ma as ma
513    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
514    pixelSize = data['pixelSize']
515    scalex = pixelSize[0]/1000.
516    scaley = pixelSize[1]/1000.
517
518    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
519    tax = np.asfarray(tax*scalex,dtype=np.float32)
520    tay = np.asfarray(tay*scaley,dtype=np.float32)
521    nI = iLim[1]-iLim[0]
522    nJ = jLim[1]-jLim[0]
527    for polygon in polygons:
528        if polygon:
531                tax.flatten(),tay.flatten(),len(polygon),polygon,tamp)))
532    if tam.shape: tam = np.reshape(tam,(nI,nJ))
533    for X,Y,diam in spots:
535    return np.array(GetTthAzm(tax,tay,data)),tam           #2-theta & azimuth arrays & position mask
536
538    import numpy.ma as ma
542    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0])))    #azimuth, 2-theta
543    tax,tay = np.dsplit(TA,2)    #azimuth, 2-theta
544    for tth,thick in rings:
546    for tth,azm,thick in arcs:
554    del(tam)
555    return tax,tay,taz
556
558    import histogram2d as h2d
559    print 'Begin image integration'
560    LUtth = data['IOtth']
561    if data['fullIntegrate']:
562        LRazm = [0,360]
563    else:
564        LRazm = data['LRazimuth']
565    numAzms = data['outAzimuths']
566    numChans = data['outChannels']
567    azmRot = data['azmthRotate']
568    Full = data['fullIntegrate']
569    Dtth = (LUtth[1]-LUtth[0])/numChans
570    Dazm = (LRazm[1]-LRazm[0])/numAzms
571    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
572    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
573    imageN = len(image)
574    Nx,Ny = data['size']
575    nXBlks = (Nx-1)/1024+1
576    nYBlks = (Ny-1)/1024+1
577    Nup = nXBlks*nYBlks*3+3
578    dlg = wx.ProgressDialog("Elapsed time","2D image integration",Nup,
579        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
580    try:
581        t0 = time.time()
582        Nup = 0
583        dlg.Update(Nup)
584        for iBlk in range(nYBlks):
585            iBeg = iBlk*1024
586            iFin = min(iBeg+1024,Ny)
587            for jBlk in range(nXBlks):
588                jBeg = jBlk*1024
589                jFin = min(jBeg+1024,Nx)
590                print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin
591                TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
592                Nup += 1
593                dlg.Update(Nup)
594                Block = image[iBeg:iFin,jBeg:jFin]
596                del TA,tam
597                Nup += 1
598                dlg.Update(Nup)
599                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
600                del tax,tay,taz
601                Nup += 1
602                dlg.Update(Nup)
603        NST = np.array(NST)
604        H0 = np.divide(H0,NST)
605        H0 = np.nan_to_num(H0)
606        del NST
607        if Dtth:
608            H2 = [tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)]
609        else:
610            H2 = LUtth
611        if Dazm:
612            H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
613            if Full:
614                H1 = H1+azmRot
615        else:
616            H1 = LRazm
617        Nup += 1
618        dlg.Update(Nup)
619        t1 = time.time()
620    finally:
621        dlg.Destroy()
622    print 'Integration complete'
623    print "Elapsed time:","%8.3f"%(t1-t0), "s"
624    return H0,H1,H2
Note: See TracBrowser for help on using the repository browser.