source:branch/MPbranch/GSASIIimage.py@5028

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

Add another variation on background - individual peaks (not finished)
fix image calibration problem

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