source: trunk/GSASIIimage.py @ 245

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