source: trunk/GSASIIimage.py @ 1164

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

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