source: trunk/GSASIIimage.py @ 606

Last change on this file since 606 was 584, checked in by vondreele, 13 years ago

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