source: trunk/GSASIIimage.py @ 1165

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

image calibration much better now - works on transposed data
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 34.9 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASII image calculations: ellipse fitting & image integration       
3########### SVN repository information ###################
4# $Date: 2013-12-14 14:54:05 +0000 (Sat, 14 Dec 2013) $
5# $Author: vondreele $
6# $Revision: 1165 $
7# $URL: trunk/GSASIIimage.py $
8# $Id: GSASIIimage.py 1165 2013-12-14 14:54:05Z vondreele $
9########### SVN repository information ###################
10'''
11*GSASIIimage: Image calc module*
12================================
13
14Ellipse fitting & image integration
15
16'''
17
18import math
19import time
20import numpy as np
21import numpy.linalg as nl
22import numpy.ma as ma
23import polymask as pm
24from scipy.optimize import leastsq
25import copy
26import GSASIIpath
27GSASIIpath.SetVersionNumber("$Revision: 1165 $")
28import GSASIIplot as G2plt
29import GSASIIlattice as G2lat
30import GSASIIpwd as G2pwd
31import fellipse as fel
32
33# trig functions in degrees
34sind = lambda x: math.sin(x*math.pi/180.)
35asind = lambda x: 180.*math.asin(x)/math.pi
36tand = lambda x: math.tan(x*math.pi/180.)
37atand = lambda x: 180.*math.atan(x)/math.pi
38atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi
39cosd = lambda x: math.cos(x*math.pi/180.)
40acosd = lambda x: 180.*math.acos(x)/math.pi
41rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p)
42#numpy versions
43npsind = lambda x: np.sin(x*np.pi/180.)
44npasind = lambda x: 180.*np.arcsin(x)/np.pi
45npcosd = lambda x: np.cos(x*np.pi/180.)
46npacosd = lambda x: 180.*np.arccos(x)/np.pi
47nptand = lambda x: np.tan(x*np.pi/180.)
48npatand = lambda x: 180.*np.arctan(x)/np.pi
49npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
50   
51def pointInPolygon(pXY,xy):
52    'Needs a doc string'
53    #pXY - assumed closed 1st & last points are duplicates
54    Inside = False
55    N = len(pXY)
56    p1x,p1y = pXY[0]
57    for i in range(N+1):
58        p2x,p2y = pXY[i%N]
59        if (max(p1y,p2y) >= xy[1] > min(p1y,p2y)) and (xy[0] <= max(p1x,p2x)):
60            if p1y != p2y:
61                xinters = (xy[1]-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
62            if p1x == p2x or xy[0] <= xinters:
63                Inside = not Inside
64        p1x,p1y = p2x,p2y
65    return Inside
66   
67def peneCorr(tth,dep):
68    'Needs a doc string'
69    return dep*(1.-npcosd(tth))         #best one
70#    return dep*npsind(tth)             #not as good as 1-cos2Q
71       
72def makeMat(Angle,Axis):
73    '''Make rotation matrix from Angle and Axis
74
75    :param float Angle: in degrees
76    :param int Axis: 0 for rotation about x, 1 for about y, etc.
77    '''
78    cs = cosd(Angle)
79    ss = sind(Angle)
80    M = np.array(([1.,0.,0.],[0.,cs,-ss],[0.,ss,cs]),dtype=np.float32)
81    return np.roll(np.roll(M,Axis,axis=0),Axis,axis=1)
82                   
83def FitEllipse(xy):
84   
85    def ellipse_center(p):
86        ''' gives ellipse center coordinates
87        '''
88        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
89        num = b*b-a*c
90        x0=(c*d-b*f)/num
91        y0=(a*f-b*d)/num
92        return np.array([x0,y0])
93   
94    def ellipse_angle_of_rotation( p ):
95        ''' gives rotation of ellipse major axis from x-axis
96        range will be -90 to 90 deg
97        '''
98        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
99        return 0.5*npatand(2*b/(a-c))
100   
101    def ellipse_axis_length( p ):
102        ''' gives ellipse radii in [minor,major] order
103        '''
104        b,c,d,f,g,a = p[1]/2, p[2], p[3]/2, p[4]/2, p[5], p[0]
105        up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
106        down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
107        down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
108        res1=np.sqrt(up/down1)
109        res2=np.sqrt(up/down2)
110        return np.array([ res2,res1])
111   
112    xy = np.array(xy)
113    x = np.asarray(xy.T[0])[:,np.newaxis]
114    y = np.asarray(xy.T[1])[:,np.newaxis]
115    D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
116    S = np.dot(D.T,D)
117    C = np.zeros([6,6])
118    C[0,2] = C[2,0] = 2; C[1,1] = -1
119    E, V =  nl.eig(np.dot(nl.inv(S), C))
120    n = np.argmax(np.abs(E))
121    a = V[:,n]
122    cent = ellipse_center(a)
123    phi = ellipse_angle_of_rotation(a)
124    radii = ellipse_axis_length(a)
125    phi += 90.
126    if radii[0] > radii[1]:
127        radii = [radii[1],radii[0]]
128        phi -= 90.
129    return cent,phi,radii
130
131def FitDetector(rings,varyList,parmDict,Print=True):
132    'Needs a doc string'
133       
134    def CalibPrint(ValSig):
135        print 'Image Parameters:'
136        ptlbls = 'names :'
137        ptstr =  'values:'
138        sigstr = 'esds  :'
139        for name,fmt,value,sig in ValSig:
140            ptlbls += "%s" % (name.rjust(12))
141            if name == 'phi':
142                ptstr += fmt % (value%360.)
143            else:
144                ptstr += fmt % (value)
145            if sig:
146                sigstr += fmt % (sig)
147            else:
148                sigstr += 12*' '
149        print ptlbls
150        print ptstr
151        print sigstr       
152       
153    def ellipseCalcD(B,xyd,varyList,parmDict):
154        x,y,dsp = xyd
155        wave = parmDict['wave']
156        if 'dep' in varyList:
157            dist,x0,y0,tilt,chi,dep = B[:6]
158        else:
159            dist,x0,y0,tilt,chi = B[:5]
160            dep = parmDict['dep']
161        if 'wave' in varyList:
162            wave = B[-1]
163        phi = chi-90.               #get rotation of major axis from tilt axis
164        tth = 2.0*npasind(wave/(2.*dsp))
165        dxy = peneCorr(tth,dep)
166        ttth = nptand(tth)
167        stth = npsind(tth)
168        cosb = npcosd(tilt)
169        tanb = nptand(tilt)       
170        tbm = nptand((tth-tilt)/2.)
171        tbp = nptand((tth+tilt)/2.)
172        sinb = npsind(tilt)
173        d = dist+dxy
174        fplus = d*tanb*stth/(cosb+stth)
175        fminus = d*tanb*stth/(cosb-stth)
176        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
177        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
178        R0 = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
179        R1 = (vplus+vminus)/2.                                    #major axis
180        zdis = (fplus-fminus)/2.
181        Robs = np.sqrt((x-x0)**2+(y-y0)**2)
182        phi0 = npatan2d(y-y0,x-x0)
183        rsqplus = R0**2+R1**2
184        rsqminus = R0**2-R1**2
185        R = rsqminus*npcosd(2.*phi0-2.*phi)+rsqplus
186        Q = np.sqrt(2.)*R0*R1*np.sqrt(R-2.*zdis**2*npsind(phi0-phi)**2)
187        P = 2.*R0**2*zdis*npcosd(phi0-phi)
188        Rcalc = (P+Q)/R
189        M = Robs-Rcalc
190        return M
191       
192    names = ['dist','det-X','det-Y','tilt','phi','dep','wave']
193    fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.3f','%12.5f']
194    p0 = [parmDict[key] for key in varyList]
195    result = leastsq(ellipseCalcD,p0,args=(rings.T,varyList,parmDict),full_output=True,ftol=1.e-8)
196    chisq = np.sum(result[2]['fvec']**2)
197    parmDict.update(zip(varyList,result[0]))
198    vals = list(result[0])
199    chi = np.sqrt(np.sum(ellipseCalcD(result[0],rings.T,varyList,parmDict)**2))
200    sig = list(chi*np.sqrt(np.diag(result[1])))
201    sigList = np.zeros(7)
202    for i,name in enumerate(names):
203        if name in varyList:
204            sigList[i] = sig[varyList.index(name)]
205    ValSig = zip(names,fmt,vals,sig)
206    if Print:
207        CalibPrint(ValSig)
208    return chi,chisq
209                   
210def ImageLocalMax(image,w,Xpix,Ypix):
211    'Needs a doc string'
212    w2 = w*2
213    sizey,sizex = image.shape
214    xpix = int(Xpix)            #get reference corner of pixel chosen
215    ypix = int(Ypix)
216    if (w < xpix < sizex-w) and (w < ypix < sizey-w) and image[ypix,xpix]:
217        Z = image[ypix-w:ypix+w,xpix-w:xpix+w]
218        Zmax = np.argmax(Z)
219        Zmin = np.argmin(Z)
220        xpix += Zmax%w2-w
221        ypix += Zmax/w2-w
222        return xpix,ypix,np.ravel(Z)[Zmax],max(0.0001,np.ravel(Z)[Zmin])   #avoid neg/zero minimum
223    else:
224        return 0,0,0,0     
225   
226def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
227    'Needs a doc string'
228    def ellipseC():
229        'compute estimate of ellipse circumference'
230        if radii[0] < 0:        #hyperbola
231            theta = npacosd(1./np.sqrt(1.+(radii[0]/radii[1])**2))
232            print theta
233            return 0
234        apb = radii[1]+radii[0]
235        amb = radii[1]-radii[0]
236        return np.pi*apb*(1+3*(amb/apb)**2/(10+np.sqrt(4-3*(amb/apb)**2)))
237    cent,phi,radii = ellipse
238    cphi = cosd(phi)
239    sphi = sind(phi)
240    ring = []
241    sumI = 0
242    amin = 0
243    amax = 360
244    C = int(ellipseC())
245    for i in range(0,C,1):
246        a = 360.*i/C
247        x = radii[0]*cosd(a)
248        y = radii[1]*sind(a)
249        X = (cphi*x+sphi*y+cent[0])*scalex      #convert mm to pixels
250        Y = (-sphi*x+cphi*y+cent[1])*scaley
251        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
252        if I and J and I/J > reject:
253            sumI += I/J
254            X += .5                             #set to center of pixel
255            Y += .5
256            X /= scalex                         #convert to mm
257            Y /= scaley
258            amin = min(amin,a)
259            amax = max(amax,a)
260            if [X,Y,dsp] not in ring:
261                ring.append([X,Y,dsp])
262    if len(ring) < 10:             #want more than 10 deg
263        ring = []
264    return ring
265   
266def GetEllipse2(tth,dxy,dist,cent,tilt,phi):
267    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
268    on output
269    radii[0] (b-minor axis) set < 0. for hyperbola
270   
271    '''
272    radii = [0,0]
273    ttth = tand(tth)
274    stth = sind(tth)
275    ctth = cosd(tth)
276    cosb = cosd(tilt)
277    tanb = tand(tilt)
278    tbm = tand((tth-tilt)/2.)
279    tbp = tand((tth+tilt)/2.)
280    sinb = sind(tilt)
281    d = dist+dxy
282    if tth+abs(tilt) < 90.:      #ellipse
283        fplus = d*tanb*stth/(cosb+stth)
284        fminus = d*tanb*stth/(cosb-stth)
285        vplus = d*(tanb+(1+tbm)/(1-tbm))*stth/(cosb+stth)
286        vminus = d*(tanb+(1-tbp)/(1+tbp))*stth/(cosb-stth)
287        radii[0] = np.sqrt((vplus+vminus)**2-(fplus+fminus)**2)/2.      #+minor axis
288        radii[1] = (vplus+vminus)/2.                                    #major axis
289        zdis = (fplus-fminus)/2.
290    else:   #hyperbola!
291        f = d*tanb*stth/(cosb+stth)
292        v = d*(tanb+tand(tth-tilt))
293        delt = d*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
294        eps = (v-f)/(delt-v)
295        radii[0] = -eps*(delt-f)/np.sqrt(eps**2-1.)                     #-minor axis
296        radii[1] = eps*(delt-f)/(eps**2-1.)                             #major axis
297        zdis = f+radii[1]*eps
298#NB: zdis is || to major axis & phi is rotation of minor axis
299#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
300    elcent = [cent[0]+zdis*sind(phi),cent[1]-zdis*cosd(phi)]
301    return elcent,phi,radii
302   
303def GetEllipse(dsp,data):
304    '''uses Dandelin spheres to find ellipse or hyperbola parameters from detector geometry
305    as given in image controls dictionary (data) and a d-spacing (dsp)
306    '''
307    cent = data['center']
308    tilt = data['tilt']
309    phi = data['rotation']
310    dep = data['DetDepth']
311    tth = 2.0*asind(data['wavelength']/(2.*dsp))
312    dxy = peneCorr(tth,dep)
313    dist = data['distance']
314    return GetEllipse2(tth,dxy,dist,cent,tilt,phi)
315       
316def GetDetectorXY(dsp,azm,data):
317    'Needs a doc string'
318   
319    elcent,phi,radii = GetEllipse(dsp,data)
320    phi = data['rotation']-90.          #to give rotation of major axis
321    tilt = data['tilt']
322    dist = data['distance']
323    cent = data['center']
324    tth = 2.0*asind(data['wavelength']/(2.*dsp))
325    ttth = tand(tth)
326    stth = sind(tth)
327    ctth = cosd(tth)
328    sinb = sind(tilt)
329    cosb = cosd(tilt)
330    tanb = tand(tilt)
331    if radii[0] > 0.:
332        fplus = dist*tanb*stth/(cosb+stth)
333        fminus = dist*tanb*stth/(cosb-stth)
334        zdis = (fplus-fminus)/2.
335        rsqplus = radii[0]**2+radii[1]**2
336        rsqminus = radii[0]**2-radii[1]**2
337        R = rsqminus*cosd(2.*azm-2.*phi)+rsqplus
338        Q = np.sqrt(2.)*radii[0]*radii[1]*np.sqrt(R-2.*zdis**2*sind(azm-phi)**2)
339        P = 2.*radii[0]**2*zdis*cosd(azm-phi)
340        radius = (P+Q)/R
341        xy = np.array([radius*cosd(azm),radius*sind(azm)])
342        xy += cent
343    else:   #hyperbola - both branches (one is way off screen!)
344        f = dist*tanb*stth/(cosb+stth)
345        v = dist*(tanb+tand(tth-abs(tilt)))
346        delt = dist*stth*(1+stth*cosb)/(sinb*cosb*(stth+cosb))
347        ecc = (v-f)/(delt-v)
348        R = radii[1]*(ecc**2-1)/(1-ecc*cosd(azm))
349        if tilt > 0.:
350            offset = 2.*radii[1]*ecc+f      #select other branch
351            xy = [-R*cosd(azm)-offset,R*sind(azm)]
352        else:
353            offset = -f
354            xy = [-R*cosd(azm)-offset,-R*sind(azm)]
355        xy = -np.array([xy[0]*cosd(phi)-xy[1]*sind(phi),xy[0]*sind(phi)+xy[1]*cosd(phi)])
356        xy += cent
357    return xy
358   
359def GetDetXYfromThAzm(Th,Azm,data):
360    'Needs a doc string'
361    dsp = data['wavelength']/(2.0*npsind(Th))   
362    return GetDetectorXY(dsp,azm,data)
363                   
364def GetTthAzmDsp(x,y,data):
365    'Needs a doc string - checked OK for ellipses dont know about hyperbola'
366    wave = data['wavelength']
367    cent = data['center']
368    tilt = data['tilt']
369    dist = data['distance']/cosd(tilt)
370    phi = data['rotation']
371    dep = data['DetDepth']
372    LRazim = data['LRazimuth']
373    azmthoff = data['azmthOff']
374    dx = np.array(x-cent[0],dtype=np.float32)
375    dy = np.array(y-cent[1],dtype=np.float32)
376    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
377    X = np.dot(X,makeMat(phi,2))
378    Z = np.dot(X,makeMat(tilt,0)).T[2]
379    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
380    dxy = peneCorr(tth,dep)
381    DX = dist-Z+dxy
382    DY = np.sqrt(dx**2+dy**2-Z**2)
383    D = (DX**2+DY**2)/dist**2       #for geometric correction = 1/cos(2theta)^2 if tilt=0.
384    tth = npatan2d(DY,DX) #depth corr not correct for tilted detector
385    dsp = wave/(2.*npsind(tth/2.))
386    azm = (npatan2d(dy,dx)+azmthoff+720.)%360.
387    return tth,azm,D,dsp
388   
389def GetTth(x,y,data):
390    'Give 2-theta value for detector x,y position; calibration info in data'
391    return GetTthAzmDsp(x,y,data)[0]
392   
393def GetTthAzm(x,y,data):
394    'Give 2-theta, azimuth values for detector x,y position; calibration info in data'
395    return GetTthAzmDsp(x,y,data)[0:2]
396   
397def GetTthAzmD(x,y,data):
398    '''Give 2-theta, azimuth & geometric correction values for detector x,y position;
399     calibration info in data
400    '''
401    return GetTthAzmDsp(x,y,data)[0:3]
402
403def GetDsp(x,y,data):
404    'Give d-spacing value for detector x,y position; calibration info in data'
405    return GetTthAzmDsp(x,y,data)[3]
406       
407def GetAzm(x,y,data):
408    'Give azimuth value for detector x,y position; calibration info in data'
409    return GetTthAzmDsp(x,y,data)[1]
410       
411def ImageCompress(image,scale):
412    'Needs a doc string'
413    if scale == 1:
414        return image
415    else:
416        return image[::scale,::scale]
417       
418def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
419    'Needs a doc string'
420    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
421    curr = np.array([dist,x,y])
422    return abs(avg-curr)/avg < .02
423
424def EdgeFinder(image,data):
425    '''this makes list of all x,y where I>edgeMin suitable for an ellipse search?
426    Not currently used but might be useful in future?
427    '''
428    import numpy.ma as ma
429    Nx,Ny = data['size']
430    pixelSize = data['pixelSize']
431    edgemin = data['edgemin']
432    scalex = pixelSize[0]/1000.
433    scaley = pixelSize[1]/1000.   
434    tay,tax = np.mgrid[0:Nx,0:Ny]
435    tax = np.asfarray(tax*scalex,dtype=np.float32)
436    tay = np.asfarray(tay*scaley,dtype=np.float32)
437    tam = ma.getmask(ma.masked_less(image.flatten(),edgemin))
438    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
439    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
440    return zip(tax,tay)
441   
442def MakeFrameMask(data,frame):
443    pixelSize = data['pixelSize']
444    scalex = pixelSize[0]/1000.
445    scaley = pixelSize[1]/1000.
446    blkSize = 512
447    Nx,Ny = data['size']
448    nXBlks = (Nx-1)/blkSize+1
449    nYBlks = (Ny-1)/blkSize+1
450    tam = ma.make_mask_none(data['size'])
451    for iBlk in range(nXBlks):
452        iBeg = iBlk*blkSize
453        iFin = min(iBeg+blkSize,Nx)
454        for jBlk in range(nYBlks):
455            jBeg = jBlk*blkSize
456            jFin = min(jBeg+blkSize,Ny)               
457            nI = iFin-iBeg
458            nJ = jFin-jBeg
459            tax,tay = np.mgrid[iBeg+0.5:iFin+.5,jBeg+.5:jFin+.5]         #bin centers not corners
460            tax = np.asfarray(tax*scalex,dtype=np.float32)
461            tay = np.asfarray(tay*scaley,dtype=np.float32)
462            tamp = ma.make_mask_none((1024*1024))
463            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
464                tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
465            if tamp.shape:
466                tamp = np.reshape(tamp[:nI*nJ],(nI,nJ))
467                tam[iBeg:iFin,jBeg:jFin] = ma.mask_or(tamp[0:nI,0:nJ],tam[iBeg:iFin,jBeg:jFin])
468            else:
469                tam[iBeg:iFin,jBeg:jFin] = True
470    return tam.T
471   
472def ImageRecalibrate(self,data,masks):
473    'Needs a doc string'
474    import ImageCalibrants as calFile
475    print 'Image recalibration:'
476    time0 = time.time()
477    pixelSize = data['pixelSize']
478    scalex = 1000./pixelSize[0]
479    scaley = 1000./pixelSize[1]
480    pixLimit = data['pixLimit']
481    cutoff = data['cutoff']
482    data['rings'] = []
483    data['ellipses'] = []
484    if not data['calibrant']:
485        print 'no calibration material selected'
486        return True   
487    skip = data['calibskip']
488    dmin = data['calibdmin']
489    Bravais,Cells = calFile.Calibrants[data['calibrant']][:2]
490    HKL = []
491    for bravais,cell in zip(Bravais,Cells):
492        A = G2lat.cell2A(cell)
493        hkl = G2lat.GenHBravais(dmin,bravais,A)
494        HKL += hkl
495    HKL = G2lat.sortHKLd(HKL,True,False)
496    varyList = ['dist','det-X','det-Y','tilt','phi']
497    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
498        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
499    Found = False
500    wave = data['wavelength']
501    frame = masks['Frames']
502    tam = ma.make_mask_none(self.ImageZ.shape)
503    if frame:
504        tam = ma.mask_or(tam,MakeFrameMask(data,frame))
505    for iH,H in enumerate(HKL): 
506        dsp = H[3]
507        tth = 2.0*asind(wave/(2.*dsp))
508        if tth+abs(data['tilt']) > 90.:
509            print 'next line is a hyperbola - search stopped'
510            break
511        ellipse = GetEllipse(dsp,data)
512        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,ma.array(self.ImageZ,mask=tam))
513        if Ring:
514            if iH >= skip:
515                data['rings'].append(np.array(Ring))
516            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
517            Found = True
518        elif not Found:         #skipping inner rings, keep looking until ring found
519            continue
520        else:                   #no more rings beyond edge of detector
521            break
522    rings = np.concatenate((data['rings']),axis=0)
523    if data['DetDepthRef']:
524        varyList.append('dep')
525    FitDetector(rings,varyList,parmDict)
526    data['distance'] = parmDict['dist']
527    data['center'] = [parmDict['det-X'],parmDict['det-Y']]
528    data['rotation'] = np.mod(parmDict['phi'],360.0)
529    data['tilt'] = parmDict['tilt']
530    data['DetDepth'] = parmDict['dep']
531    N = len(data['ellipses'])
532    data['ellipses'] = []           #clear away individual ellipse fits
533    for H in HKL[:N]:
534        ellipse = GetEllipse(H[3],data)
535        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
536   
537    print 'calibration time = ',time.time()-time0
538    G2plt.PlotImage(self,newImage=True)       
539    return True
540           
541def ImageCalibrate(self,data):
542    'Needs a doc string'
543    import copy
544    import ImageCalibrants as calFile
545    print 'Image calibration:'
546    time0 = time.time()
547    ring = data['ring']
548    pixelSize = data['pixelSize']
549    scalex = 1000./pixelSize[0]
550    scaley = 1000./pixelSize[1]
551    pixLimit = data['pixLimit']
552    cutoff = data['cutoff']
553    if len(ring) < 5:
554        print 'not enough inner ring points for ellipse'
555        return False
556       
557    #fit start points on inner ring
558    data['ellipses'] = []
559    data['rings'] = []
560    outE = FitEllipse(ring)
561    fmt  = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f'
562    fmt2 = '%s X: %.3f, Y: %.3f, phi: %.3f, R1: %.3f, R2: %.3f, chi^2: %.3f, N: %d'
563    if outE:
564        print fmt%('start ellipse: ',outE[0][0],outE[0][1],outE[1],outE[2][0],outE[2][1])
565        ellipse = outE
566    else:
567        return False
568       
569    #setup 360 points on that ring for "good" fit
570    data['ellipses'].append(ellipse[:]+('g',))
571    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
572    if Ring:
573        ellipse = FitEllipse(ring)
574        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
575        ellipse = FitEllipse(ring)
576    else:
577        print '1st ring not sufficiently complete to proceed'
578        return False
579    print fmt2%('inner ring:    ',ellipse[0][0],ellipse[0][1],ellipse[1],ellipse[2][0],ellipse[2][1],0.,len(Ring))     #cent,phi,radii
580    data['ellipses'].append(ellipse[:]+('r',))
581    data['rings'].append(np.array(Ring))
582    G2plt.PlotImage(self,newImage=True)
583   
584#setup for calibration
585    data['rings'] = []
586    if not data['calibrant']:
587        print 'no calibration material selected'
588        return True
589   
590    skip = data['calibskip']
591    dmin = data['calibdmin']
592#generate reflection set
593    Bravais,Cells = calFile.Calibrants[data['calibrant']][:2]
594    HKL = []
595    for bravais,cell in zip(Bravais,Cells):
596        A = G2lat.cell2A(cell)
597        hkl = G2lat.GenHBravais(dmin,bravais,A)[skip:]
598        HKL += hkl
599    HKL = G2lat.sortHKLd(HKL,True,False)
600    wave = data['wavelength']
601#set up 1st ring
602    elcent,phi,radii = ellipse              #from fit of 1st ring
603    dsp = HKL[0][3]
604    tth = 2.0*asind(wave/(2.*dsp))
605    Ring0 = makeRing(dsp,ellipse,3,cutoff,scalex,scaley,self.ImageZ)
606    ttth = nptand(tth)
607    stth = npsind(tth)
608    ctth = npcosd(tth)
609#1st estimate of tilt; assume ellipse - don't know sign though
610    tilt = npasind(np.sqrt(max(0.,1.-(radii[0]/radii[1])**2))*ctth)
611    if not tilt:
612        print 'WARNING - selected ring was fitted as a circle'
613        print ' - if detector was tilted we suggest you skip this ring - WARNING'
614#1st estimate of dist: sample to detector normal to plane
615    data['distance'] = dist = radii[0]**2/(ttth*radii[1])
616#ellipse to cone axis (x-ray beam); 2 choices depending on sign of tilt
617    zdisp = radii[1]*ttth*tand(tilt)
618    zdism = radii[1]*ttth*tand(-tilt)
619#cone axis position; 2 choices. Which is right?     
620#NB: zdisp is || to major axis & phi is rotation of minor axis
621#thus shift from beam to ellipse center is [Z*sin(phi),-Z*cos(phi)]
622    centp = [elcent[0]+zdisp*sind(phi),elcent[1]-zdisp*cosd(phi)]
623    centm = [elcent[0]+zdism*sind(phi),elcent[1]-zdism*cosd(phi)]
624#check get same ellipse parms either way
625#now do next ring; estimate either way & do a FitDetector each way; best fit is correct one
626    dsp = HKL[1][3]
627    tth = 2.0*asind(wave/(2.*dsp))
628    ellipsep = GetEllipse2(tth,0.,dist,centp,tilt,phi)
629    print fmt%('plus ellipse :',ellipsep[0][0],ellipsep[0][1],ellipsep[1],ellipsep[2][0],ellipsep[2][1])
630    Ringp = makeRing(dsp,ellipsep,3,cutoff,scalex,scaley,self.ImageZ)
631    parmDict = {'dist':dist,'det-X':centp[0],'det-Y':centp[1],
632        'tilt':tilt,'phi':phi,'wave':wave,'dep':0.0}
633    varyList = ['dist','det-X','det-Y','tilt','phi']
634    if len(Ringp) > 10:
635        chip,chisqp = FitDetector(np.array(Ring0+Ringp),varyList,parmDict,True)
636        tiltp = parmDict['tilt']
637        phip = parmDict['phi']
638        centp = [parmDict['det-X'],parmDict['det-Y']]
639    else:
640        chip = 1e6
641    ellipsem = GetEllipse2(tth,0.,dist,centm,-tilt,phi)
642    print fmt%('minus ellipse:',ellipsem[0][0],ellipsem[0][1],ellipsem[1],ellipsem[2][0],ellipsem[2][1])
643    Ringm = makeRing(dsp,ellipsem,3,cutoff,scalex,scaley,self.ImageZ)
644    if len(Ringm) > 10:
645        parmDict['tilt'] *= -1
646        chim,chisqm = FitDetector(np.array(Ring0+Ringm),varyList,parmDict,True)
647        tiltm = parmDict['tilt']
648        phim = parmDict['phi']
649        centm = [parmDict['det-X'],parmDict['det-Y']]
650    else:
651        chim = 1e6
652    if chip < chim:
653        data['tilt'] = tiltp
654        data['center'] = centp
655        data['rotation'] = phip
656    else:
657        data['tilt'] = tiltm
658        data['center'] = centm
659        data['rotation'] = phim
660    data['ellipses'].append(ellipsep[:]+('b',))
661    data['rings'].append(np.array(Ringp))
662    data['ellipses'].append(ellipsem[:]+('r',))
663    data['rings'].append(np.array(Ringm))
664    G2plt.PlotImage(self,newImage=True)
665    parmDict = {'dist':data['distance'],'det-X':data['center'][0],'det-Y':data['center'][1],
666        'tilt':data['tilt'],'phi':data['rotation'],'wave':data['wavelength'],'dep':data['DetDepth']}
667    varyList = ['dist','det-X','det-Y','tilt','phi']
668    if data['DetDepthRef']:
669        varyList.append('dep')
670    data['rings'] = []
671    data['ellipses'] = []
672    for i,H in enumerate(HKL):
673        dsp = H[3]
674        tth = 2.0*asind(wave/(2.*dsp))
675        if tth+abs(data['tilt']) > 90.:
676            print 'next line is a hyperbola - search stopped'
677            break
678        print 'HKLD:',H[:4],'2-theta: %.4f'%(tth)
679        elcent,phi,radii = ellipse = GetEllipse(dsp,data)
680        data['ellipses'].append(copy.deepcopy(ellipse+('g',)))
681        print fmt%('predicted ellipse:',elcent[0],elcent[1],phi,radii[0],radii[1])
682        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
683        if Ring:
684            data['rings'].append(np.array(Ring))
685            rings = np.concatenate((data['rings']),axis=0)
686            if i:
687                chi,chisq = FitDetector(rings,varyList,parmDict,False)
688                data['distance'] = parmDict['dist']
689                data['center'] = [parmDict['det-X'],parmDict['det-Y']]
690                data['rotation'] = parmDict['phi']
691                data['tilt'] = parmDict['tilt']
692                data['DetDepth'] = parmDict['dep']
693                elcent,phi,radii = ellipse = GetEllipse(dsp,data)
694                print fmt2%('fitted ellipse:   ',elcent[0],elcent[1],phi,radii[0],radii[1],chisq,len(rings))
695            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
696#            G2plt.PlotImage(self,newImage=True)
697        else:
698            print 'insufficient number of points in this ellipse to fit'
699            break
700    G2plt.PlotImage(self,newImage=True)
701    fullSize = len(self.ImageZ)/scalex
702    if 2*radii[1] < .9*fullSize:
703        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
704    N = len(data['ellipses'])
705    if N > 2:
706        FitDetector(rings,varyList,parmDict)
707        data['distance'] = parmDict['dist']
708        data['center'] = [parmDict['det-X'],parmDict['det-Y']]
709        data['rotation'] = parmDict['phi']
710        data['tilt'] = parmDict['tilt']
711        data['DetDepth'] = parmDict['dep']
712    for H in HKL[:N]:
713        ellipse = GetEllipse(H[3],data)
714        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
715    print 'calibration time = ',time.time()-time0
716    G2plt.PlotImage(self,newImage=True)       
717    return True
718   
719def Make2ThetaAzimuthMap(data,masks,iLim,jLim):
720    'Needs a doc string'
721    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
722    pixelSize = data['pixelSize']
723    scalex = pixelSize[0]/1000.
724    scaley = pixelSize[1]/1000.
725   
726    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
727    tax = np.asfarray(tax*scalex,dtype=np.float32)
728    tay = np.asfarray(tay*scaley,dtype=np.float32)
729    nI = iLim[1]-iLim[0]
730    nJ = jLim[1]-jLim[0]
731    #make position masks here
732    frame = masks['Frames']
733    tam = ma.make_mask_none((nI,nJ))
734    if frame:
735        tamp = ma.make_mask_none((1024*1024))
736        tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
737            tay.flatten(),len(frame),frame,tamp)[:nI*nJ])-True  #switch to exclude around frame
738        tam = ma.mask_or(tam.flatten(),tamp)
739    polygons = masks['Polygons']
740    for polygon in polygons:
741        if polygon:
742            tamp = ma.make_mask_none((1024*1024))
743            tamp = ma.make_mask(pm.polymask(nI*nJ,tax.flatten(),
744                tay.flatten(),len(polygon),polygon,tamp)[:nI*nJ])
745            tam = ma.mask_or(tam.flatten(),tamp)
746    if tam.shape: tam = np.reshape(tam,(nI,nJ))
747    spots = masks['Points']
748    for X,Y,diam in spots:
749        tamp = ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2))
750        tam = ma.mask_or(tam,tamp)
751    TA = np.array(GetTthAzmD(tax,tay,data))     #includes geom. corr.
752    TA[1] = np.where(TA[1]<0,TA[1]+360,TA[1])
753    return np.array(TA),tam           #2-theta, azimuth & geom. corr. arrays & position mask
754
755def Fill2ThetaAzimuthMap(masks,TA,tam,image):
756    'Needs a doc string'
757    Zlim = masks['Thresholds'][1]
758    rings = masks['Rings']
759    arcs = masks['Arcs']
760    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0]),ma.getdata(TA[2])))    #azimuth, 2-theta, dist
761    tax,tay,tad = np.dsplit(TA,3)    #azimuth, 2-theta, dist
762    for tth,thick in rings:
763        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
764    for tth,azm,thick in arcs:
765        tamt = ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))
766        tama = ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1]))
767        tam = ma.mask_or(tam.flatten(),tamt*tama)
768    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
769    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
770    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
771    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
772    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))
773    tad = ma.compressed(ma.array(tad.flatten(),mask=tam))
774    return tax,tay,taz      #*tad**2 wrong - why?
775   
776def ImageIntegrate(image,data,masks,blkSize=128,dlg=None):
777    'Needs a doc string'
778    import histogram2d as h2d
779    print 'Begin image integration'
780    LUtth = data['IOtth']
781    LRazm = np.array(data['LRazimuth'],dtype=np.float64)
782    numAzms = data['outAzimuths']
783    numChans = data['outChannels']
784    Dtth = (LUtth[1]-LUtth[0])/numChans
785    Dazm = (LRazm[1]-LRazm[0])/numAzms
786    if data['centerAzm']:
787        LRazm += Dazm/2.
788    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
789    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
790    imageN = len(image)
791    Nx,Ny = data['size']
792    nXBlks = (Nx-1)/blkSize+1
793    nYBlks = (Ny-1)/blkSize+1
794    Nup = nXBlks*nYBlks*3+3
795    t0 = time.time()
796    Nup = 0
797    if dlg:
798        dlg.Update(Nup)
799    for iBlk in range(nYBlks):
800        iBeg = iBlk*blkSize
801        iFin = min(iBeg+blkSize,Ny)
802        for jBlk in range(nXBlks):
803            jBeg = jBlk*blkSize
804            jFin = min(jBeg+blkSize,Nx)               
805            TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
806           
807            Nup += 1
808            if dlg:
809                dlg.Update(Nup)
810            Block = image[iBeg:iFin,jBeg:jFin]
811            tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
812            Nup += 1
813            if dlg:
814                dlg.Update(Nup)
815            tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
816            tax = np.where(tax < LRazm[0],tax+360.,tax)
817            if any([tax.shape[0],tay.shape[0],taz.shape[0]]):
818                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
819            Nup += 1
820            if dlg:
821                dlg.Update(Nup)
822    NST = np.array(NST)
823    H0 = np.divide(H0,NST)
824    H0 = np.nan_to_num(H0)
825    del NST
826    if Dtth:
827        H2 = np.array([tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)])
828    else:
829        H2 = np.array(LUtth)
830    if Dazm:       
831        H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
832    else:
833        H1 = LRazm
834    H0 /= npcosd(H2[:-1])**2
835    if data['Oblique'][1]:
836        H0 /= G2pwd.Oblique(data['Oblique'][0],H2[:-1])
837    if 'SASD' in data['type'] and data['PolaVal'][1]:
838        H0 /= np.array([G2pwd.Polarization(data['PolaVal'][0],H2[:-1],Azm=azm)[0] for azm in H1[:-1]+np.diff(H1)/2.])
839    Nup += 1
840    if dlg:
841        dlg.Update(Nup)
842    t1 = time.time()
843    print 'Integration complete'
844    print "Elapsed time:","%8.3f"%(t1-t0), "s"
845    return H0,H1,H2
846   
847def FitStrSta(Image,StrSta,Controls,Masks):
848    'Needs a doc string'
849   
850#    print 'Masks:',Masks
851    StaControls = copy.deepcopy(Controls)
852    phi = StrSta['Sample phi']
853    wave = Controls['wavelength']
854    StaControls['distance'] += StrSta['Sample z']*cosd(phi)
855    pixSize = StaControls['pixelSize']
856    scalex = 1000./pixSize[0]
857    scaley = 1000./pixSize[1]
858    rings = []
859
860    for ring in StrSta['d-zero']:       #get observed x,y,d points for the d-zeros
861        ellipse = GetEllipse(ring['Dset'],StaControls)
862        Ring = makeRing(ring['Dset'],ellipse,ring['pixLimit'],ring['cutoff'],scalex,scaley,Image)
863        Ring = np.array(Ring).T
864        ring['ImxyObs'] = np.array(Ring[:2])      #need to apply masks to this to eliminate bad points
865        Ring[:2] = GetTthAzm(Ring[0],Ring[1],StaControls)       #convert x,y to tth,azm
866        Ring[0] /= 2.                                           #convert to theta
867        if len(rings):
868            rings = np.concatenate((rings,Ring),axis=1)
869        else:
870            rings = np.array(Ring)
871    E = StrSta['strain']
872    p0 = [E[0][0],E[1][1],E[2][2],E[0][1],E[0][2],E[1][2]]
873    E = FitStrain(rings,p0,wave,phi)
874    StrSta['strain'] = E
875
876def calcFij(omg,phi,azm,th):
877    '''Does something...
878
879    Uses parameters as defined by Bob He & Kingsley Smith, Adv. in X-Ray Anal. 41, 501 (1997)
880
881    :param omg: his omega = sample omega rotation; 0 when incident beam || sample surface, 90 when perp. to sample surface
882    :param phi: his phi = sample phi rotation; usually = 0, axis rotates with omg.
883    :param azm: his chi = azimuth around incident beam
884    :param th:  his theta = theta
885    '''
886    a = npsind(th)*npcosd(omg)+npsind(azm)*npcosd(th)*npsind(omg)
887    b = -npcosd(azm)*npcosd(th)
888    c = npsind(th)*npsind(omg)-npsind(azm)*npcosd(th)*npcosd(omg)
889    d = a*npcosd(phi)-b*npsind(phi)
890    e = a*npsind(phi)+b*npcosd(phi)
891    Fij = np.array([
892        [d**2,d*e,c*d],
893        [d*e,e**2,c*e],
894        [c*d,c*e,c**2]])
895    return -Fij*nptand(th)
896
897def FitStrain(rings,p0,wave,phi):
898    'Needs a doc string'
899    def StrainPrint(ValSig):
900        print 'Strain tensor:'
901        ptlbls = 'names :'
902        ptstr =  'values:'
903        sigstr = 'esds  :'
904        for name,fmt,value,sig in ValSig:
905            ptlbls += "%s" % (name.rjust(12))
906            ptstr += fmt % (value)
907            if sig:
908                sigstr += fmt % (sig)
909            else:
910                sigstr += 12*' '
911        print ptlbls
912        print ptstr
913        print sigstr
914       
915    def strainCalc(p,xyd,wave,phi):
916#        E = np.array([[p[0],p[3],p[4]],[p[3],p[1],p[5]],[p[4],p[5],p[2]]])
917        E = np.array([[p[0],0,0],[0,p[1],0],[0,0,0]])
918        th,azm,dsp = xyd
919        th0 = npasind(wave/(2.*dsp))
920        dth = 180.*np.sum(E*calcFij(phi,0.,azm,th).T)/(np.pi*1.e6) #in degrees & microstrain units
921        th0 += dth
922        return (th-th0)**2
923       
924    names = ['e11','e22','e33','e12','e13','e23']
925    fmt = ['%12.2f','%12.2f','%12.2f','%12.2f','%12.2f','%12.5f']
926    p1 = [p0[0],p0[1]]   
927    result = leastsq(strainCalc,p1,args=(rings,wave,phi),full_output=True)
928    vals = list(result[0])
929    chi = np.sqrt(np.sum(strainCalc(result[0],rings,wave,phi)**2))
930    sig = list(chi*np.sqrt(np.diag(result[1])))
931    ValSig = zip(names,fmt,vals,sig)
932    StrainPrint(ValSig)
933#    return np.array([[vals[0],vals[3],vals[4]],[vals[3],vals[1],vals[5]],[vals[4],vals[5],vals[2]]])
934    return np.array([[vals[0],0,0],[0,vals[1],0],[0,0,0]])
935   
936       
Note: See TracBrowser for help on using the repository browser.