source: trunk/GSASIIimage.py @ 675

Last change on this file since 675 was 675, checked in by vondreele, 10 years ago

refactor image controls GUI
add oblique detector absorption correction to integration

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