source: trunk/GSASIIimage.py @ 251

Last change on this file since 251 was 251, checked in by vondreele, 11 years ago

Add NaCl? even hkl to ImageCalibrants?.py
Fix bug in polygon mask in GSASIIimage.py
Fix copy bugs in GSASIIimgGUI.py
Make float image for float tiff file in GSASIIIO.py

  • 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-03-01 19:16:27 +0000 (Tue, 01 Mar 2011) $
4# $Author: vondreele $
5# $Revision: 251 $
6# $URL: trunk/GSASIIimage.py $
7# $Id: GSASIIimage.py 251 2011-03-01 19:16:27Z 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        if polygon:
507            tamp = ma.make_mask_none((nI*nJ))
508            tam = ma.mask_or(tam.flatten(),ma.make_mask(pm.polymask(nI*nJ,
509                tax.flatten(),tay.flatten(),len(polygon),polygon,tamp)))
510    if tam.shape: tam = np.reshape(tam,(nI,nJ))
511    for X,Y,diam in spots:
512        tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2)))
513    return np.array(GetTthAzm(tax,tay,data)),tam           #2-theta & azimuth arrays & position mask
514
515def Fill2ThetaAzimuthMap(masks,TA,tam,image):
516    import numpy.ma as ma
517    Zlim = masks['Thresholds'][1]
518    rings = masks['Rings']
519    arcs = masks['Arcs']
520    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0])))    #azimuth, 2-theta
521    tax,tay = np.dsplit(TA,2)    #azimuth, 2-theta
522    for tth,thick in rings:
523        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
524    for tth,azm,thick in arcs:
525        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))* \
526            ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1])))
527    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
528    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
529    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
530    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
531    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))
532    del(tam)
533    return tax,tay,taz
534   
535def ImageIntegrate(image,data,masks):
536    import histogram2d as h2d
537    print 'Begin image integration'
538    LUtth = data['IOtth']
539    if data['fullIntegrate']:
540        LRazm = [-180,180]
541    else:
542        LRazm = data['LRazimuth']
543    numAzms = data['outAzimuths']
544    numChans = data['outChannels']
545    Dtth = (LUtth[1]-LUtth[0])/numChans
546    Dazm = (LRazm[1]-LRazm[0])/numAzms
547    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
548    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
549    imageN = len(image)
550    Nx,Ny = data['size']
551    nXBlks = (Nx-1)/1024+1
552    nYBlks = (Ny-1)/1024+1
553    Nup = nXBlks*nYBlks*3+3
554    dlg = wx.ProgressDialog("Elapsed time","2D image integration",Nup,
555        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
556    try:
557        t0 = time.time()
558        Nup = 0
559        dlg.Update(Nup)
560        for iBlk in range(nYBlks):
561            iBeg = iBlk*1024
562            iFin = min(iBeg+1024,Ny)
563            for jBlk in range(nXBlks):
564                jBeg = jBlk*1024
565                jFin = min(jBeg+1024,Nx)               
566                print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin
567                TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
568                Nup += 1
569                dlg.Update(Nup)
570                Block = image[iBeg:iFin,jBeg:jFin]
571                tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
572                del TA,tam
573                Nup += 1
574                dlg.Update(Nup)
575                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
576                del tax,tay,taz
577                Nup += 1
578                dlg.Update(Nup)
579        NST = np.array(NST)
580        H0 = np.divide(H0,NST)
581        H0 = np.nan_to_num(H0)
582        del NST
583        if Dtth:
584            H2 = [tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)]
585        else:
586            H2 = LUtth
587        if Dazm:       
588            H1 = [azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)]
589        else:
590            H1 = LRazm
591        Nup += 1
592        dlg.Update(Nup)
593        t1 = time.time()
594    finally:
595        dlg.Destroy()
596    print 'Integration complete'
597    print "Elapsed time:","%8.3f"%(t1-t0), "s"
598    return H0,H1,H2
Note: See TracBrowser for help on using the repository browser.