source: trunk/GSASIIimage.py @ 94

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

move all fortran to fsource
add gray for readonly textctrl
do polygon insert, show in imageGUI and make mask
fix to ReadPDBphase to use correct crystal to cartesian transformation matrix

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