source: trunk/GSASIIimage.py @ 652

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

Add delete map peaks menu item
correct integration result for flat detector
remove print of map peak list
more work on findOffset

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