source: trunk/GSASIIimage.py @ 661

Last change on this file since 661 was 661, checked in by toby, 10 years ago

complete revision number tracking

  • Property svn:keywords set to Date Author Revision URL Id
File size: 25.9 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2012-06-29 21:24:59 +0000 (Fri, 29 Jun 2012) $
5# $Author: toby $
6# $Revision: 661 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 661 2012-06-29 21:24:59Z toby $
9########### SVN repository information ###################
10import math
11import wx
12import time
13import numpy as np
14import numpy.linalg as nl
15import GSASIIpath
16GSASIIpath.SetVersionNumber("$Revision: 661 $")
17import GSASIIplot as G2plt
18import GSASIIlattice as G2lat
19import fellipse as fel
20
21# trig functions in degrees
22sind = lambda x: math.sin(x*math.pi/180.)
23asind = lambda x: 180.*math.asin(x)/math.pi
24tand = lambda x: math.tan(x*math.pi/180.)
25atand = lambda x: 180.*math.atan(x)/math.pi
26atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
27cosd = lambda x: math.cos(x*math.pi/180.)
28acosd = lambda x: 180.*math.acos(x)/math.pi
29rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
30#numpy versions
31npsind = lambda x: np.sin(x*np.pi/180.)
32npasind = lambda x: 180.*np.arcsin(x)/np.pi
33npcosd = lambda x: np.cos(x*np.pi/180.)
34nptand = lambda x: np.tan(x*np.pi/180.)
35npatand = lambda x: 180.*np.arctan(x)/np.pi
36npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
37   
38def pointInPolygon(pXY,xy):
39    #pXY - assumed closed 1st & last points are duplicates
40    Inside = False
41    N = len(pXY)
42    p1x,p1y = pXY[0]
43    for i in range(N+1):
44        p2x,p2y = pXY[i%N]
45        if (max(p1y,p2y) >= xy[1] > min(p1y,p2y)) and (xy[0] <= max(p1x,p2x)):
46            if p1y != p2y:
47                xinters = (xy[1]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
48            if p1x == p2x or xy[0] <= xinters:
49                Inside = not Inside
50        p1x,p1y = p2x,p2y
51    return Inside
52       
53def makeMat(Angle,Axis):
54    #Make rotation matrix from Angle in degrees,Axis =0 for rotation about x, =1 for about y, etc.
55    cs = cosd(Angle)
56    ss = sind(Angle)
57    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
58    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
59                   
60def FitRing(ring,delta):
61    parms = []
62    if delta:
63        err,parms = FitEllipse(ring)
64        errc,parmsc = FitCircle(ring)
65        errc = errc[0]/(len(ring)*parmsc[2][0]**2)
66        if not parms or errc < .1:
67            parms = parmsc
68    else:
69        err,parms = FitCircle(ring)
70    return parms
71       
72def FitCircle(ring):
73   
74    def makeParmsCircle(B):
75        cent = [-B[0]/2,-B[1]/2]
76        phi = 0.
77        sr1 = sr2 = math.sqrt(cent[0]**2+cent[1]**2-B[2])
78        return cent,phi,[sr1,sr2]
79       
80    ring = np.array(ring)
81    x = np.asarray(ring.T[0])
82    y = np.asarray(ring.T[1])
83   
84    M = np.array((x,y,np.ones_like(x)))
85    B = np.array(-(x**2+y**2))
86    result = nl.lstsq(M.T,B)
87    return result[1],makeParmsCircle(result[0])
88       
89def FitEllipse(ring):
90           
91    def makeParmsEllipse(B):
92        det = 4.*(1.-B[0]**2)-B[1]**2
93        if det < 0.:
94            print 'hyperbola!'
95            return 0
96        elif det == 0.:
97            print 'parabola!'
98            return 0
99        cent = [(B[1]*B[3]-2.*(1.-B[0])*B[2])/det, \
100            (B[1]*B[2]-2.*(1.+B[0])*B[3])/det]
101        phi = 0.5*atand(0.5*B[1]/B[0])
102       
103        a = (1.+B[0])/cosd(2*phi)
104        b = 2.-a
105        f = (1.+B[0])*cent[0]**2+(1.-B[0])*cent[1]**2+B[1]*cent[0]*cent[1]-B[4]
106        if f/a < 0 or f/b < 0:
107            return 0
108        sr1 = math.sqrt(f/a)
109        sr2 = math.sqrt(f/b)
110        if sr1 > sr2:
111            sr1,sr2 = sr2,sr1
112            phi -= 90.
113            if phi < -90.:
114                phi += 180.
115        return cent,phi,[sr1,sr2]
116               
117    ring = np.array(ring)
118    x = np.asarray(ring.T[0])
119    y = np.asarray(ring.T[1])
120    M = np.array((x**2-y**2,x*y,x,y,np.ones_like(x)))
121    B = np.array(-(x**2+y**2))
122    bb,err = fel.fellipse(len(x),x,y,1.E-7)
123#    print nl.lstsq(M.T,B)[0]
124#    print bb
125    return err,makeParmsEllipse(bb)
126   
127def FitDetector(rings,p0,wave):
128    from scipy.optimize import leastsq
129       
130    def CalibPrint(ValSig):
131        print 'Image Parameters:'
132        ptlbls = 'names :'
133        ptstr =  'values:'
134        sigstr = 'esds  :'
135        for name,fmt,value,sig in ValSig:
136            ptlbls += "%s" % (name.rjust(12))
137            if name == 'rotate':
138                ptstr += fmt % (value-90.)      #kluge to get rotation from vertical - see GSASIIimgGUI
139            else:
140                ptstr += fmt % (value)
141            if sig:
142                sigstr += fmt % (sig)
143            else:
144                sigstr += 12*' '
145        print ptlbls
146        print ptstr
147        print sigstr       
148       
149    def ellipseCalcD(B,xyd,wave):
150        x = xyd[0]
151        y = xyd[1]
152        dsp = xyd[2]
153        dist,x0,y0,phi,tilt = B
154        tth = 2.0*npasind(wave/(2.*dsp))
155        ttth = nptand(tth)
156        radius = dist*ttth
157        stth = npsind(tth)
158        cosb = npcosd(tilt)
159        R1 = dist*stth*npcosd(tth)*cosb/(cosb**2-stth**2)
160        R0 = np.sqrt(R1*radius*cosb)
161        zdis = R1*ttth*nptand(tilt)
162        X = x-x0+zdis*npsind(phi)
163        Y = y-y0-zdis*npcosd(phi)
164        XR = X*npcosd(phi)-Y*npsind(phi)
165        YR = X*npsind(phi)+Y*npcosd(phi)
166        return (XR/R0)**2+(YR/R1)**2-1
167       
168    def ellipseCalcW(C,xyd):
169        dist,x0,y0,phi,tilt,wave = C
170        B = dist,x0,y0,phi,tilt
171        return ellipseCalcD(B,xyd,wave)
172       
173    names = ['distance','det-X0','det-Y0','rotate','tilt','wavelength']
174    fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.5f']   
175    result = leastsq(ellipseCalcD,p0,args=(rings.T,wave),full_output=True)
176    result[0][3] = np.mod(result[0][3],360.0)               #remove random full circles
177    vals = list(result[0])
178    vals.append(wave)
179    chi = np.sqrt(np.sum(ellipseCalcD(result[0],rings.T,wave)**2))
180    sig = list(chi*np.sqrt(np.diag(result[1])))
181    sig.append(0.)
182    ValSig = zip(names,fmt,vals,sig)
183    CalibPrint(ValSig)
184    try:
185        print 'Trial refinement of wavelength - not used for calibration'
186        p0 = result[0]
187        p0 = np.append(p0,wave)
188        resultW = leastsq(ellipseCalcW,p0,args=(rings.T,),full_output=True)
189        resultW[0][3] = np.mod(result[0][3],360.0)          #remove random full circles
190        sig = np.sqrt(np.sum(ellipseCalcW(resultW[0],rings.T)**2))
191        ValSig = zip(names,fmt,resultW[0],sig*np.sqrt(np.diag(resultW[1])))
192        CalibPrint(ValSig)
193        return result[0],resultW[0][-1]       
194    except ValueError:
195        print 'Bad wavelength refinement - no result'
196        return result[0],wave
197           
198def ImageLocalMax(image,w,Xpix,Ypix):
199    w2 = w*2
200    sizey,sizex = image.shape
201    xpix = int(Xpix)            #get reference corner of pixel chosen
202    ypix = int(Ypix)
203    if (w < xpix < sizex-w) and (w < ypix < sizey-w) and image[ypix,xpix]:
204        Z = image[ypix-w:ypix+w,xpix-w:xpix+w]
205        Zmax = np.argmax(Z)
206        Zmin = np.argmin(Z)
207        xpix += Zmax%w2-w
208        ypix += Zmax/w2-w
209        return xpix,ypix,np.ravel(Z)[Zmax],max(0.0001,np.ravel(Z)[Zmin])   #avoid neg/zero minimum
210    else:
211        return 0,0,0,0     
212   
213def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
214    cent,phi,radii = ellipse
215    cphi = cosd(phi)
216    sphi = sind(phi)
217    ring = []
218    amin = 0
219    amax = 360
220    for a in range(0,360,1):
221        x = radii[0]*cosd(a)
222        y = radii[1]*sind(a)
223        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
224        Y = (sphi*x+cphi*y+cent[1])*scaley
225        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
226        if I and J and I/J > reject:
227            X += .5                             #set to center of pixel
228            Y += .5
229            X /= scalex                         #convert to mm
230            Y /= scaley
231            amin = min(amin,a)
232            amax = max(amax,a)
233            ring.append([X,Y,dsp])
234    delt = amax-amin
235    if len(ring) < 20:             #want more than 20 deg
236        return [],delt > 90
237    return ring,delt > 90
238   
239def makeIdealRing(ellipse,azm=None):
240    cent,phi,radii = ellipse
241    cphi = cosd(phi)
242    sphi = sind(phi)
243    if azm:
244        aR = azm[0]-90,azm[1]-90,azm[1]-azm[0]
245        if azm[1]-azm[0] > 180:
246            aR[2] /= 2
247    else:
248        aR = 0,362,181
249       
250    a = np.linspace(aR[0],aR[1],aR[2])
251    x = radii[0]*npcosd(a-phi)
252    y = radii[1]*npsind(a-phi)
253    X = (cphi*x-sphi*y+cent[0])
254    Y = (sphi*x+cphi*y+cent[1])
255    return zip(X,Y)
256               
257def calcDist(radii,tth):
258    stth = sind(tth)
259    ctth = cosd(tth)
260    ttth = tand(tth)
261    return math.sqrt(radii[0]**4/(ttth**2*((radii[0]*ctth)**2+(radii[1]*stth)**2)))
262   
263def calcZdisCosB(radius,tth,radii):
264    cosB = sinb = radii[0]**2/(radius*radii[1])
265    if cosB > 1.:
266        return 0.,1.
267    else:
268        cosb = math.sqrt(1.-sinb**2)
269        ttth = tand(tth)
270        zdis = radii[1]*ttth*cosb/sinb
271        return zdis,cosB
272   
273def GetEllipse(dsp,data):
274    dist = data['distance']
275    cent = data['center']
276    tilt = data['tilt']
277    phi = data['rotation']
278    radii = [0,0]
279    tth = 2.0*asind(data['wavelength']/(2.*dsp))
280    ttth = tand(tth)
281    stth = sind(tth)
282    ctth = cosd(tth)
283    cosb = cosd(tilt)
284    radius = dist*ttth
285    radii[1] = dist*stth*ctth*cosb/(cosb**2-stth**2)
286    if radii[1] > 0:
287        radii[0] = math.sqrt(radii[1]*radius*cosb)
288        zdis = radii[1]*ttth*tand(tilt)
289        elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)]
290        return elcent,phi,radii
291    else:
292        return False
293       
294def GetDetectorXY(dsp,azm,data):
295    from scipy.optimize import fsolve
296    def func(xy,*args):
297       azm,phi,R0,R1,A,B = args
298       cp = cosd(phi)
299       sp = sind(phi)
300       x,y = xy
301       out = []
302       out.append(y-x*tand(azm))
303       out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)
304       return out
305    elcent,phi,radii = GetEllipse(dsp,data)
306    cent = data['center']
307    tilt = data['tilt']
308    phi = data['rotation']
309    wave = data['wavelength']
310    dist = data['distance']
311    tth = 2.0*asind(wave/(2.*dsp))
312    ttth = tand(tth)
313    radius = dist*ttth
314    stth = sind(tth)
315    cosb = cosd(tilt)
316    R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2)
317    R0 = math.sqrt(R1*radius*cosb)
318    zdis = R1*ttth*tand(tilt)
319    A = zdis*sind(phi)
320    B = -zdis*cosd(phi)
321    xy0 = [radius*cosd(azm),radius*sind(azm)]
322    xy = fsolve(func,xy0,args=(azm,phi,R0,R1,A,B))+cent
323    return xy
324                   
325def GetTthAzmDsp(x,y,data):
326    wave = data['wavelength']
327    dist = data['distance']
328    cent = data['center']
329    tilt = data['tilt']
330    phi = data['rotation']
331    LRazim = data['LRazimuth']
332    azmthoff = data['azmthOff']
333    dx = np.array(x-cent[0],dtype=np.float32)
334    dy = np.array(y-cent[1],dtype=np.float32)
335    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
336    X = np.dot(X,makeMat(phi,2))
337    Z = np.dot(X,makeMat(tilt,0)).T[2]
338    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
339    dsp = wave/(2.*npsind(tth/2.))
340    azm = (npatan2d(dx,-dy)+azmthoff+720.)%360.
341    return tth,azm,dsp
342   
343def GetTth(x,y,data):
344    return GetTthAzmDsp(x,y,data)[0]
345   
346def GetTthAzm(x,y,data):
347    return GetTthAzmDsp(x,y,data)[0:2]
348   
349def GetDsp(x,y,data):
350    return GetTthAzmDsp(x,y,data)[2]
351       
352def GetAzm(x,y,data):
353    return GetTthAzmDsp(x,y,data)[1]
354       
355def ImageCompress(image,scale):
356    if scale == 1:
357        return image
358    else:
359        return image[::scale,::scale]
360       
361def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
362    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
363    curr = np.array([dist,x,y])
364    return abs(avg-curr)/avg < .02
365
366def EdgeFinder(image,data):          #this makes list of all x,y where I>edgeMin suitable for an ellipse search?
367    import numpy.ma as ma
368    Nx,Ny = data['size']
369    pixelSize = data['pixelSize']
370    edgemin = data['edgemin']
371    scalex = pixelSize[0]/1000.
372    scaley = pixelSize[1]/1000.   
373    tay,tax = np.mgrid[0:Nx,0:Ny]
374    tax = np.asfarray(tax*scalex,dtype=np.float32)
375    tay = np.asfarray(tay*scaley,dtype=np.float32)
376    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
377    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
378    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
379    return zip(tax,tay)
380   
381def ImageRecalibrate(self,data):
382    import copy
383    import ImageCalibrants as calFile
384    print 'Image recalibration:'
385    time0 = time.time()
386    pixelSize = data['pixelSize']
387    scalex = 1000./pixelSize[0]
388    scaley = 1000./pixelSize[1]
389    pixLimit = data['pixLimit']
390    cutoff = data['cutoff']
391    data['rings'] = []
392    data['ellipses'] = []
393    if not data['calibrant']:
394        print 'no calibration material selected'
395        return True
396   
397    skip = data['calibskip']
398    dmin = data['calibdmin']
399    Bravais,Cells = calFile.Calibrants[data['calibrant']][:2]
400    HKL = []
401    for bravais,cell in zip(Bravais,Cells):
402        A = G2lat.cell2A(cell)
403        hkl = G2lat.GenHBravais(dmin,bravais,A)[skip:]
404        HKL += hkl
405    HKL = G2lat.sortHKLd(HKL,True,False)
406    wave = data['wavelength']
407    cent = data['center']   
408    dist = data['distance']
409    cent = data['center']
410    tilt = data['tilt']
411    phi = data['rotation']
412    for H in HKL: 
413        dsp = H[3]
414        ellipse = GetEllipse(dsp,data)
415        Ring,delt = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
416        if Ring:
417            numZ = len(Ring)
418            data['rings'].append(np.array(Ring))
419            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
420        else:
421            continue
422    rings = np.concatenate((data['rings']),axis=0)
423    p0 = [dist,cent[0],cent[1],phi,tilt]
424    result,newWave = FitDetector(rings,p0,wave)
425    data['distance'] = result[0]
426    data['center'] = result[1:3]
427    data['rotation'] = np.mod(result[3],360.0)
428    data['tilt'] = result[4]
429    N = len(data['ellipses'])
430    data['ellipses'] = []           #clear away individual ellipse fits
431    for H in HKL[:N]:
432        ellipse = GetEllipse(H[3],data)
433        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
434   
435    print 'calibration time = ',time.time()-time0
436    G2plt.PlotImage(self,newImage=True)       
437    return True
438   
439           
440def ImageCalibrate(self,data):
441    import copy
442    import ImageCalibrants as calFile
443    print 'Image calibration:'
444    time0 = time.time()
445    ring = data['ring']
446    pixelSize = data['pixelSize']
447    scalex = 1000./pixelSize[0]
448    scaley = 1000./pixelSize[1]
449    pixLimit = data['pixLimit']
450    cutoff = data['cutoff']
451    if len(ring) < 5:
452        print 'not enough inner ring points for ellipse'
453        return False
454       
455    #fit start points on inner ring
456    data['ellipses'] = []
457    outE = FitRing(ring,True)
458    if outE:
459#        print 'start ellipse:',outE
460        ellipse = outE
461    else:
462        return False
463       
464    #setup 360 points on that ring for "good" fit
465    Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
466    if Ring:
467        ellipse = FitRing(Ring,delt)
468        Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
469        ellipse = FitRing(Ring,delt)
470    else:
471        print '1st ring not sufficiently complete to proceed'
472        return False
473#    print 'inner ring:',ellipse
474    data['center'] = copy.copy(ellipse[0])           #not right!! (but useful for now)
475    data['ellipses'].append(ellipse[:]+('r',))
476    G2plt.PlotImage(self,newImage=True)
477   
478    #setup for calibration
479    data['rings'] = []
480    data['ellipses'] = []
481    if not data['calibrant']:
482        print 'no calibration material selected'
483        return True
484   
485    skip = data['calibskip']
486    dmin = data['calibdmin']
487    Bravais,Cells = calFile.Calibrants[data['calibrant']][:2]
488    HKL = []
489    for bravais,cell in zip(Bravais,Cells):
490        A = G2lat.cell2A(cell)
491        hkl = G2lat.GenHBravais(dmin,bravais,A)[skip:]
492        HKL += hkl
493    HKL = G2lat.sortHKLd(HKL,True,False)
494    wave = data['wavelength']
495    cent = data['center']
496    elcent,phi,radii = ellipse
497    dsp = HKL[0][3]
498    tth = 2.0*asind(wave/(2.*dsp))
499    ttth = tand(tth)
500    data['distance'] = dist = calcDist(radii,tth)
501    radius = dist*tand(tth)
502    zdis,cosB = calcZdisCosB(radius,tth,radii)
503    cent1 = []
504    cent2 = []
505    xSum = 0
506    ySum = 0
507    zxSum = 0
508    zySum = 0
509    phiSum = 0
510    tiltSum = 0
511    distSum = 0
512    Zsum = 0
513    for i,H in enumerate(HKL):
514        dsp = H[3]
515        tth = 2.0*asind(0.5*wave/dsp)
516        stth = sind(tth)
517        ctth = cosd(tth)
518        ttth = tand(tth)
519        radius = dist*ttth
520        elcent,phi,radii = ellipse
521        if i:
522            radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2)
523            radii[0] = math.sqrt(radii[1]*radius*cosB)
524            zdis,cosB = calcZdisCosB(radius,tth,radii)
525            zsinp = zdis*sind(phi)
526            zcosp = zdis*cosd(phi)
527            cent = data['center']
528            elcent = [cent[0]-zsinp,cent[1]+zcosp]
529            ellipse = (elcent,phi,radii)
530        ratio = radii[1]/radii[0]
531        Ring,delt = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
532        if Ring:
533            numZ = len(Ring)
534            data['rings'].append(np.array(Ring))
535            newellipse = FitRing(Ring,delt)
536            elcent,phi,radii = newellipse               
537            if abs(phi) > 45. and phi < 0.:
538                phi += 180.
539            dist = calcDist(radii,tth)
540            distR = 1.-dist/data['distance']
541            if abs(distR) > 0.1:
542#                print dsp,dist,data['distance'],distR,len(Ring),delt
543                break
544            if distR > 0.001:
545                print 'Wavelength too large?'
546            elif distR < -0.001:
547                print 'Wavelength too small?'
548            else:
549                ellipse = newellipse
550            zdis,cosB = calcZdisCosB(radius,tth,radii)
551            Tilt = acosd(cosB)          # 0 <= tilt <= 90
552            zsinp = zdis*sind(ellipse[1])
553            zcosp = zdis*cosd(ellipse[1])
554            cent1.append(np.array([elcent[0]+zsinp,elcent[1]-zcosp]))
555            cent2.append(np.array([elcent[0]-zsinp,elcent[1]+zcosp]))
556            if i:
557                d1 = cent1[-1]-cent1[-2]        #get shift of 2 possible center solutions
558                d2 = cent2[-1]-cent2[-2]
559                if np.dot(d2,d2) > np.dot(d1,d1):  #right solution is the larger shift
560                    data['center'] = cent1[-1]
561                else:
562                    data['center'] = cent2[-1]
563                Zsum += numZ
564                phiSum += numZ*phi
565                distSum += numZ*dist
566                xSum += numZ*data['center'][0]
567                ySum += numZ*data['center'][1]
568                tiltSum += numZ*abs(Tilt)
569                if not np.all(checkEllipse(Zsum,distSum,xSum,ySum,dist,data['center'][0],data['center'][1])):
570                    print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i) 
571            cent = data['center']
572#            print ('for ring # %2i @ d-space %.4f: dist %.3f rotate %6.2f tilt %6.2f Xcent %.3f Ycent %.3f Npts %d'
573#                %(i,dsp,dist,phi-90.,Tilt,cent[0],cent[1],numZ))
574            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
575        else:
576            break
577    G2plt.PlotImage(self,newImage=True)
578    fullSize = len(self.ImageZ)/scalex
579    if 2*radii[1] < .9*fullSize:
580        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
581    if not Zsum:
582        print 'Only one ring fitted. Check your wavelength.'
583        return False
584    cent = data['center'] = [xSum/Zsum,ySum/Zsum]
585    wave = data['wavelength']
586    dist = data['distance'] = distSum/Zsum
587   
588    #possible error if no. of rings < 3! Might need trap here
589    d1 = cent1[-1]-cent1[1]             #compare last ring to 2nd ring
590    d2 = cent2[-1]-cent2[1]
591    Zsign = 1
592    len1 = math.sqrt(np.dot(d1,d1))
593    len2 = math.sqrt(np.dot(d2,d2))
594    t1 = d1/len1
595    t2 = d2/len2
596    if len2 > len1:
597        if -135. < atan2d(t2[1],t2[0]) < 45.:
598            Zsign = -1
599    else:
600        if -135. < atan2d(t1[1],t1[0]) < 45.:
601            Zsign = -1
602   
603    tilt = data['tilt'] = Zsign*tiltSum/Zsum
604    phi = data['rotation'] = phiSum/Zsum
605    rings = np.concatenate((data['rings']),axis=0)
606    p0 = [dist,cent[0],cent[1],phi,tilt]
607    result,newWave = FitDetector(rings,p0,wave)
608    data['distance'] = result[0]
609    data['center'] = result[1:3]
610    data['rotation'] = np.mod(result[3],360.0)
611    data['tilt'] = result[4]
612    N = len(data['ellipses'])
613    data['ellipses'] = []           #clear away individual ellipse fits
614    for H in HKL[:N]:
615        ellipse = GetEllipse(H[3],data)
616        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
617    print 'calibration time = ',time.time()-time0
618    G2plt.PlotImage(self,newImage=True)       
619    return True
620   
621def Make2ThetaAzimuthMap(data,masks,iLim,jLim):
622    import numpy.ma as ma
623    import polymask as pm
624    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
625    pixelSize = data['pixelSize']
626    scalex = pixelSize[0]/1000.
627    scaley = pixelSize[1]/1000.
628   
629    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
630    tax = np.asfarray(tax*scalex,dtype=np.float32)
631    tay = np.asfarray(tay*scaley,dtype=np.float32)
632    nI = iLim[1]-iLim[0]
633    nJ = jLim[1]-jLim[0]
634    #make position masks here
635    spots = masks['Points']
636    polygons = masks['Polygons']
637    tam = ma.make_mask_none((nI,nJ))
638    for polygon in polygons:
639        if polygon:
640            tamp = ma.make_mask_none((nI*nJ))
641            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
642                tay.flatten(),len(polygon),polygon,tamp))
643            tam = ma.mask_or(tam.flatten(),tamp)
644    if tam.shape: tam = np.reshape(tam,(nI,nJ))
645    for X,Y,diam in spots:
646        tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2))
647        tam = ma.mask_or(tam,tamp)
648    TA = np.array(GetTthAzm(tax,tay,data))
649    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
650    return np.array(TA),tam           #2-theta & azimuth arrays & position mask
651
652def Fill2ThetaAzimuthMap(masks,TA,tam,image):
653    import numpy.ma as ma
654    Zlim = masks['Thresholds'][1]
655    rings = masks['Rings']
656    arcs = masks['Arcs']
657    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0])))    #azimuth, 2-theta
658    tax,tay = np.dsplit(TA,2)    #azimuth, 2-theta
659    for tth,thick in rings:
660        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
661    for tth,azm,thick in arcs:
662        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
663        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
664        tam = ma.mask_or(tam.flatten(),tamt*tama)
665    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
666    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
667    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
668    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
669    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))
670    del(tam)
671    return tax,tay,taz
672   
673def ImageIntegrate(image,data,masks):
674    import histogram2d as h2d
675    print 'Begin image integration'
676    LUtth = data['IOtth']
677    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
678    numAzms = data['outAzimuths']
679    numChans = data['outChannels']
680    Dtth = (LUtth[1]-LUtth[0])/numChans
681    Dazm = (LRazm[1]-LRazm[0])/numAzms
682    if data['centerAzm']:
683        LRazm += Dazm/2.
684    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
685    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
686    imageN = len(image)
687    Nx,Ny = data['size']
688    nXBlks = (Nx-1)/1024+1
689    nYBlks = (Ny-1)/1024+1
690    Nup = nXBlks*nYBlks*3+3
691    dlg = wx.ProgressDialog("Elapsed time","2D image integration",Nup,
692        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
693    try:
694        t0 = time.time()
695        Nup = 0
696        dlg.Update(Nup)
697        for iBlk in range(nYBlks):
698            iBeg = iBlk*1024
699            iFin = min(iBeg+1024,Ny)
700            for jBlk in range(nXBlks):
701                jBeg = jBlk*1024
702                jFin = min(jBeg+1024,Nx)               
703                print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin
704                TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
705               
706                Nup += 1
707                dlg.Update(Nup)
708                Block = image[iBeg:iFin,jBeg:jFin]
709                tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
710                del TA,tam
711                Nup += 1
712                dlg.Update(Nup)
713                tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
714                tax = np.where(tax < LRazm[0],tax+360.,tax)
715                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
716                print 'block done'
717                del tax,tay,taz
718                Nup += 1
719                dlg.Update(Nup)
720        NST = np.array(NST)
721        H0 = np.divide(H0,NST)
722        H0 = np.nan_to_num(H0)
723        del NST
724        if Dtth:
725            H2 = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
726        else:
727            H2 = np.array(LUtth)
728        if Dazm:       
729            H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
730        else:
731            H1 = LRazm
732        H0[0] /= npcosd(H2[:-1])
733        Nup += 1
734        dlg.Update(Nup)
735        t1 = time.time()
736    finally:
737        dlg.Destroy()
738    print 'Integration complete'
739    print "Elapsed time:","%8.3f"%(t1-t0), "s"
740    return H0,H1,H2
Note: See TracBrowser for help on using the repository browser.