source: trunk/GSASIIimage.py @ 400

Last change on this file since 400 was 397, checked in by vondreele, 14 years ago

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