source: trunk/GSASIIimage.py @ 706

Last change on this file since 706 was 706, checked in by vondreele, 9 years ago

fix sumimage - missing code!
more work on strain fitting

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