source:trunk/GSASIIimage.py@248

Last change on this file since 248 was 248, checked in by vondreele, 12 years ago

histogram2d.for - cleaned up
add azmthOff - as for a rotation of a 2D detector
fix image integration problems

• 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-17 19:33:02 +0000 (Thu, 17 Feb 2011) \$
4# \$Author: vondreele \$
5# \$Revision: 248 \$
6# \$URL: trunk/GSASIIimage.py \$
7# \$Id: GSASIIimage.py 248 2011-02-17 19:33:02Z 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    azmthoff = data['azmthOff']
284    dx = np.array(x-cent[0],dtype=np.float32)
285    dy = np.array(y-cent[1],dtype=np.float32)
286    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
287    X = np.dot(X,makeMat(phi,2))
288    Z = np.dot(X,makeMat(tilt,0)).T[2]
289    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
290    dsp = wave/(2.*npsind(tth/2.))
291    azm = npatan2d(dy,dx)+azmthoff
292    return tth,azm,dsp
293
294def GetTth(x,y,data):
295    return GetTthAzmDsp(x,y,data)[0]
296
297def GetTthAzm(x,y,data):
298    return GetTthAzmDsp(x,y,data)[0:2]
299
300def GetDsp(x,y,data):
301    return GetTthAzmDsp(x,y,data)[2]
302
303def GetAzm(x,y,data):
304    return GetTthAzmDsp(x,y,data)[1]
305
306def ImageCompress(image,scale):
307    if scale == 1:
308        return image
309    else:
310        return image[::scale,::scale]
311
312def ImageCalibrate(self,data):
313    import copy
314    import ImageCalibrants as calFile
315    print 'image calibrate'
316    ring = data['ring']
317    pixelSize = data['pixelSize']
318    scalex = 1000./pixelSize[0]
319    scaley = 1000./pixelSize[1]
320    pixLimit = data['pixLimit']
321    cutoff = data['cutoff']
322    if len(ring) < 5:
323        print 'not enough inner ring points for ellipse'
324        return False
325
326    #fit start points on inner ring
327    data['ellipses'] = []
328    outE = FitRing(ring)
329    if outE:
330        print 'start ellipse:',outE
331        ellipse = outE
332    else:
333        return False
334
335    #setup 180 points on that ring for "good" fit
336    Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
337    if Ring:
338        ellipse = FitRing(Ring)
339        Ring = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
340        ellipse = FitRing(Ring)
341    else:
342        print '1st ring not sufficiently complete to proceed'
343        return False
344    print 'inner ring:',ellipse
345    data['center'] = copy.copy(ellipse[0])           #not right!! (but useful for now)
346    data['ellipses'].append(ellipse[:]+('r',))
347    G2plt.PlotImage(self)
348
349    #setup for calibration
350    data['rings'] = []
351    data['ellipses'] = []
352    if not data['calibrant']:
353        print 'no calibration material selected'
354        return True
355
356    skip = data['calibskip']
357    dmin = data['calibdmin']
358    Bravais,cell = calFile.Calibrants[data['calibrant']][:2]
359    A = G2lat.cell2A(cell)
360    wave = data['wavelength']
361    cent = data['center']
362    elcent,phi,radii = ellipse
363    HKL = G2lat.GenHBravais(dmin,Bravais,A)[skip:]
364    dsp = HKL[0][3]
365    tth = 2.0*asind(wave/(2.*dsp))
366    ttth = tand(tth)
367    data['distance'] = dist = calcDist(radii,tth)
368    radius = dist*tand(tth)
369    zdis,cosB = calcZdisCosB(radius,tth,radii)
370    cent1 = []
371    cent2 = []
372    xSum = 0
373    ySum = 0
374    zxSum = 0
375    zySum = 0
376    phiSum = 0
377    tiltSum = 0
378    distSum = 0
379    Zsum = 0
380    for i,H in enumerate(HKL):
381        dsp = H[3]
382        tth = 2.0*asind(0.5*wave/dsp)
383        stth = sind(tth)
384        ctth = cosd(tth)
385        ttth = tand(tth)
386        radius = dist*ttth
387        elcent,phi,radii = ellipse
388        radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2)
389        radii[0] = math.sqrt(radii[1]*radius*cosB)
390        zdis,cosB = calcZdisCosB(radius,tth,radii)
391        zsinp = zdis*sind(phi)
392        zcosp = zdis*cosd(phi)
393        cent = data['center']
394        elcent = [cent[0]+zsinp,cent[1]-zcosp]
395        ratio = radii[1]/radii[0]
396        Ring = makeRing(dsp,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
397        if Ring:
398            numZ = len(Ring)
399            data['rings'].append(np.array(Ring))
400            newellipse = FitRing(Ring)
401            elcent,phi,radii = newellipse
402            if abs(phi) > 45. and phi < 0.:
403                phi += 180.
404            dist = calcDist(radii,tth)
405            distR = 1.-dist/data['distance']
406            if abs(distR) > 0.1:
407                print distR,dist,data['distance']
408                del data['rings'][-1]
409                continue
410            if distR > 0.001:
411                print 'Wavelength too large?'
412            elif distR < -0.001:
413                print 'Wavelength too small?'
414            else:
415                ellipse = newellipse
416                if abs((radii[1]/radii[0]-ratio)/ratio) > 0.01:
417                    print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i)
418                    return False
419            zdis,cosB = calcZdisCosB(radius,tth,radii)
420            Tilt = acosd(cosB)          # 0 <= tilt <= 90
421            zsinp = zdis*sind(ellipse[1])
422            zcosp = zdis*cosd(ellipse[1])
423            cent1.append(np.array([elcent[0]+zsinp,elcent[1]-zcosp]))
424            cent2.append(np.array([elcent[0]-zsinp,elcent[1]+zcosp]))
425            if i:
426                d1 = cent1[-1]-cent1[-2]        #get shift of 2 possible center solutions
427                d2 = cent2[-1]-cent2[-2]
428                if np.dot(d2,d2) > np.dot(d1,d1):  #right solution is the larger shift
429                    data['center'] = cent1[-1]
430                else:
431                    data['center'] = cent2[-1]
432                Zsum += numZ
433                phiSum += numZ*phi
434                distSum += numZ*dist
435                xSum += numZ*data['center'][0]
436                ySum += numZ*data['center'][1]
437                tiltSum += numZ*abs(Tilt)
438            cent = data['center']
439            print ('for ring # %2i dist %.3f rotate %6.2f tilt %6.2f Xcent %.3f Ycent %.3f Npts %d'
440                %(i,dist,phi,Tilt,cent[0],cent[1],numZ))
441            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
442            G2plt.PlotImage(self)
443        else:
444            break
445    fullSize = len(self.ImageZ)/scalex
446    if 2*radii[1] < .9*fullSize:
447        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
448    if not Zsum:
449        print 'Only one ring fitted. Check your wavelength.'
450        return False
451    cent = data['center'] = [xSum/Zsum,ySum/Zsum]
452    wave = data['wavelength']
453    dist = data['distance'] = distSum/Zsum
454
455    #possible error if no. of rings < 3! Might need trap here
456    d1 = cent1[-1]-cent1[1]             #compare last ring to 2nd ring
457    d2 = cent2[-1]-cent2[1]
458    Zsign = 1
459    len1 = math.sqrt(np.dot(d1,d1))
460    len2 = math.sqrt(np.dot(d2,d2))
461    t1 = d1/len1
462    t2 = d2/len2
463    if len2 > len1:
464        if -135. < atan2d(t2[1],t2[0]) < 45.:
465            Zsign = -1
466    else:
467        if -135. < atan2d(t1[1],t1[0]) < 45.:
468            Zsign = -1
469
470    tilt = data['tilt'] = Zsign*tiltSum/Zsum
471    phi = data['rotation'] = phiSum/Zsum
472    rings = np.concatenate((data['rings']),axis=0)
473    p0 = [dist,cent[0],cent[1],phi,tilt]
474    result,newWave = FitDetector(rings,p0,wave)
475    print 'Suggested new wavelength = ',('%.5f')%(newWave),' (not reliable if distance > 2m)'
476    data['distance'] = result[0]
477    data['center'] = result[1:3]
478    data['rotation'] = np.mod(result[3],360.0)
479    data['tilt'] = result[4]
480    N = len(data['ellipses'])
481    data['ellipses'] = []           #clear away individual ellipse fits
482    for H in HKL[:N]:
483        ellipse = GetEllipse(H[3],data)
484        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
485    G2plt.PlotImage(self)
486    return True
487
488def Make2ThetaAzimuthMap(data,masks,iLim,jLim):
489    import numpy.ma as ma
490    import polymask as pm
491    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
492    pixelSize = data['pixelSize']
493    scalex = pixelSize[0]/1000.
494    scaley = pixelSize[1]/1000.
495
496    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
497    tax = np.asfarray(tax*scalex,dtype=np.float32)
498    tay = np.asfarray(tay*scaley,dtype=np.float32)
499    nI = iLim[1]-iLim[0]
500    nJ = jLim[1]-jLim[0]
501    #make position masks here
502    spots = masks['Points']
503    polygons = masks['Polygons']
504    tam = ma.make_mask_none((iLim[1]-iLim[0],jLim[1]-jLim[0]))
505    for polygon in polygons:
506        tamp = ma.make_mask_none((nI*nJ))
507        tam = ma.mask_or(tam.flatten(),ma.make_mask(pm.polymask(nI*nJ,
508            tax.flatten(),tay.flatten(),len(polygon),polygon,tamp)))
509    if tam.shape: tam = np.reshape(tam,(nI,nJ))
510    for X,Y,diam in spots:
511        tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2)))
512    return np.array(GetTthAzm(tax,tay,data)),tam           #2-theta & azimuth arrays & position mask
513
514def Fill2ThetaAzimuthMap(masks,TA,tam,image):
515    import numpy.ma as ma
516    Zlim = masks['Thresholds'][1]
517    rings = masks['Rings']
518    arcs = masks['Arcs']
519    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0])))    #azimuth, 2-theta
520    tax,tay = np.dsplit(TA,2)    #azimuth, 2-theta
521    for tth,thick in rings:
522        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
523    for tth,azm,thick in arcs:
524        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))* \
525            ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1])))
526    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
527    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
528    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
529    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
530    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))
531    del(tam)
532    return tax,tay,taz
533
534def ImageIntegrate(image,data,masks):
535    import histogram2d as h2d
536    print 'Begin image integration'
537    LUtth = data['IOtth']
538    if data['fullIntegrate']:
539        LRazm = [-180,180]
540    else:
541        LRazm = data['LRazimuth']
542    numAzms = data['outAzimuths']
543    numChans = data['outChannels']
544    Dtth = (LUtth[1]-LUtth[0])/numChans
545    Dazm = (LRazm[1]-LRazm[0])/numAzms
546    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
547    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
548    imageN = len(image)
549    Nx,Ny = data['size']
550    nXBlks = (Nx-1)/1024+1
551    nYBlks = (Ny-1)/1024+1
552    Nup = nXBlks*nYBlks*3+3
553    dlg = wx.ProgressDialog("Elapsed time","2D image integration",Nup,
554        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
555    try:
556        t0 = time.time()
557        Nup = 0
558        dlg.Update(Nup)
559        for iBlk in range(nYBlks):
560            iBeg = iBlk*1024
561            iFin = min(iBeg+1024,Ny)
562            for jBlk in range(nXBlks):
563                jBeg = jBlk*1024
564                jFin = min(jBeg+1024,Nx)
565                print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin
566                TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
567                Nup += 1
568                dlg.Update(Nup)
569                Block = image[iBeg:iFin,jBeg:jFin]
570                tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
571                del TA,tam
572                Nup += 1
573                dlg.Update(Nup)
574                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
575                del tax,tay,taz
576                Nup += 1
577                dlg.Update(Nup)
578        NST = np.array(NST)
579        H0 = np.divide(H0,NST)
580        H0 = np.nan_to_num(H0)
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.