source: trunk/GSASIIimage.py @ 271

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

make g77 versions of fortran for binwin2.7 - faster than gfortran!
fix to histogram2d.for
GSASII.py - allow 2D offset in multiplots, fix for large file lists
GSASIIimgGUI.py - more stuff in save/load controls
GSASIIplot.py - 2D offset in multiplots, fix bad behavior in structure drawings
GSASIIpwd.py - more mods to PDF calcs
  • Property svn:keywords set to Date Author Revision URL Id
File size: 21.4 KB
Line 
1#GSASII image calculations: ellipse fitting & image integration       
2########### SVN repository information ###################
3# $Date: 2011-04-28 18:11:32 +0000 (Thu, 28 Apr 2011) $
4# $Author: vondreele $
5# $Revision: 271 $
6# $URL: trunk/GSASIIimage.py $
7# $Id: GSASIIimage.py 271 2011-04-28 18:11:32Z 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 ellipseCalcD(B,xyd,wave):
129        x = xyd[0]
130        y = xyd[1]
131        dsp = xyd[2]
132        dist,x0,y0,phi,tilt = B
133        tth = 2.0*npasind(wave/(2.*dsp))
134        ttth = nptand(tth)
135        radius = dist*ttth
136        stth = npsind(tth)
137        cosb = npcosd(tilt)
138        R1 = dist*stth*npcosd(tth)*cosb/(cosb**2-stth**2)
139        R0 = np.sqrt(R1*radius*cosb)
140        zdis = R1*ttth*nptand(tilt)
141        X = x-x0+zdis*npsind(phi)
142        Y = y-y0-zdis*npcosd(phi)
143        XR = X*npcosd(phi)-Y*npsind(phi)
144        YR = X*npsind(phi)+Y*npcosd(phi)
145        return (XR/R0)**2+(YR/R1)**2-1
146    def ellipseCalcW(C,xyd):
147        dist,x0,y0,phi,tilt,wave = C
148        B = dist,x0,y0,phi,tilt
149        return ellipseCalcD(B,xyd,wave)
150    result = leastsq(ellipseCalcD,p0,args=(rings.T,wave))
151    if len(rings) > 1:
152        p0 = result[0]
153        p0 = np.append(p0,wave)
154        resultW = leastsq(ellipseCalcW,p0,args=(rings.T,))
155        return result[0],resultW[0][-1]
156    else:
157        return result[0],wave
158           
159def ImageLocalMax(image,w,Xpix,Ypix):
160    w2 = w*2
161    sizey,sizex = image.shape
162    xpix = int(Xpix)            #get reference corner of pixel chosen
163    ypix = int(Ypix)
164    if (w < xpix < sizex-w) and (w < ypix < sizey-w) and image[ypix,xpix]:
165        Z = image[ypix-w:ypix+w,xpix-w:xpix+w]
166        Zmax = np.argmax(Z)
167        Zmin = np.argmin(Z)
168        xpix += Zmax%w2-w
169        ypix += Zmax/w2-w
170        return xpix,ypix,np.ravel(Z)[Zmax],np.ravel(Z)[Zmin]
171    else:
172        return 0,0,0,0
173   
174def makeRing(dsp,ellipse,pix,reject,scalex,scaley,image):
175    cent,phi,radii = ellipse
176    cphi = cosd(phi)
177    sphi = sind(phi)
178    ring = []
179    amin = 180
180    amax = -180
181    for a in range(-180,180,1):
182        x = radii[0]*cosd(a)
183        y = radii[1]*sind(a)
184        X = (cphi*x-sphi*y+cent[0])*scalex      #convert mm to pixels
185        Y = (sphi*x+cphi*y+cent[1])*scaley
186        X,Y,I,J = ImageLocalMax(image,pix,X,Y)
187        if I and J and I/J > reject:
188            X += .5                             #set to center of pixel
189            Y += .5
190            X /= scalex                         #convert to mm
191            Y /= scaley
192            amin = min(amin,a)
193            amax = max(amax,a)
194            ring.append([X,Y,dsp])
195    delt = amax-amin
196    if len(ring) < 20:             #want more than 20 deg
197        return [],delt > 90
198    return ring,delt > 90
199   
200def makeIdealRing(ellipse,azm=None):
201    cent,phi,radii = ellipse
202    cphi = cosd(phi)
203    sphi = sind(phi)
204    ring = []
205    if azm:
206        aR = azm[0]-90,azm[1]-90,1
207        if azm[1]-azm[0] > 180:
208            aR[2] = 2
209    else:
210        aR = 0,362,2
211    for a in range(aR[0],aR[1],aR[2]):
212        x = radii[0]*cosd(a-phi)
213        y = radii[1]*sind(a-phi)
214        X = (cphi*x-sphi*y+cent[0])
215        Y = (sphi*x+cphi*y+cent[1])
216        ring.append([X,Y])
217    return ring
218   
219def calcDist(radii,tth):
220    stth = sind(tth)
221    ctth = cosd(tth)
222    ttth = tand(tth)
223    return math.sqrt(radii[0]**4/(ttth**2*((radii[0]*ctth)**2+(radii[1]*stth)**2)))
224   
225def calcZdisCosB(radius,tth,radii):
226    cosB = sinb = radii[0]**2/(radius*radii[1])
227    if cosB > 1.:
228        return 0.,1.
229    else:
230        cosb = math.sqrt(1.-sinb**2)
231        ttth = tand(tth)
232        zdis = radii[1]*ttth*cosb/sinb
233        return zdis,cosB
234   
235def GetEllipse(dsp,data):
236    dist = data['distance']
237    cent = data['center']
238    tilt = data['tilt']
239    phi = data['rotation']
240    radii = [0,0]
241    tth = 2.0*asind(data['wavelength']/(2.*dsp))
242    ttth = tand(tth)
243    stth = sind(tth)
244    ctth = cosd(tth)
245    cosb = cosd(tilt)
246    radius = dist*ttth
247    radii[1] = dist*stth*ctth*cosb/(cosb**2-stth**2)
248    if radii[1] > 0:
249        radii[0] = math.sqrt(radii[1]*radius*cosb)
250        zdis = radii[1]*ttth*tand(tilt)
251        elcent = [cent[0]-zdis*sind(phi),cent[1]+zdis*cosd(phi)]
252        return elcent,phi,radii
253    else:
254        return False
255       
256def GetDetectorXY(dsp,azm,data):
257    from scipy.optimize import fsolve
258    def func(xy,*args):
259       azm,phi,R0,R1,A,B = args
260       cp = cosd(phi)
261       sp = sind(phi)
262       x,y = xy
263       out = []
264       out.append(y-x*tand(azm))
265       out.append(R0**2*((x+A)*sp-(y+B)*cp)**2+R1**2*((x+A)*cp+(y+B)*sp)**2-(R0*R1)**2)
266       return out
267    elcent,phi,radii = GetEllipse(dsp,data)
268    cent = data['center']
269    tilt = data['tilt']
270    phi = data['rotation']
271    wave = data['wavelength']
272    dist = data['distance']
273    tth = 2.0*asind(wave/(2.*dsp))
274    ttth = tand(tth)
275    radius = dist*ttth
276    stth = sind(tth)
277    cosb = cosd(tilt)
278    R1 = dist*stth*cosd(tth)*cosb/(cosb**2-stth**2)
279    R0 = math.sqrt(R1*radius*cosb)
280    zdis = R1*ttth*tand(tilt)
281    A = zdis*sind(phi)
282    B = -zdis*cosd(phi)
283    xy0 = [radius*cosd(azm),radius*sind(azm)]
284    xy = fsolve(func,xy0,args=(azm,phi,R0,R1,A,B))+cent
285    return xy
286                   
287def GetTthAzmDsp(x,y,data):
288    wave = data['wavelength']
289    dist = data['distance']
290    cent = data['center']
291    tilt = data['tilt']
292    phi = data['rotation']
293    LRazim = data['LRazimuth']
294    azmthoff = data['azmthOff']
295    dx = np.array(x-cent[0],dtype=np.float32)
296    dy = np.array(y-cent[1],dtype=np.float32)
297    X = np.array(([dx,dy,np.zeros_like(dx)]),dtype=np.float32).T
298    X = np.dot(X,makeMat(phi,2))
299    Z = np.dot(X,makeMat(tilt,0)).T[2]
300    tth = npatand(np.sqrt(dx**2+dy**2-Z**2)/(dist-Z))
301    dsp = wave/(2.*npsind(tth/2.))
302    azm = (npatan2d(dx,-dy)+azmthoff+720.)%360.
303    azm = np.where(azm>180,azm-360.,azm)
304    return tth,azm,dsp
305   
306def GetTth(x,y,data):
307    return GetTthAzmDsp(x,y,data)[0]
308   
309def GetTthAzm(x,y,data):
310    return GetTthAzmDsp(x,y,data)[0:2]
311   
312def GetDsp(x,y,data):
313    return GetTthAzmDsp(x,y,data)[2]
314       
315def GetAzm(x,y,data):
316    return GetTthAzmDsp(x,y,data)[1]
317       
318def ImageCompress(image,scale):
319    if scale == 1:
320        return image
321    else:
322        return image[::scale,::scale]
323       
324def checkEllipse(Zsum,distSum,xSum,ySum,dist,x,y):
325    avg = np.array([distSum/Zsum,xSum/Zsum,ySum/Zsum])
326    curr = np.array([dist,x,y])
327    return abs(avg-curr)/avg < .02
328
329def ImageCalibrate(self,data):
330    import copy
331    import ImageCalibrants as calFile
332    print 'image calibrate'
333    time0 = time.time()
334    ring = data['ring']
335    pixelSize = data['pixelSize']
336    scalex = 1000./pixelSize[0]
337    scaley = 1000./pixelSize[1]
338    pixLimit = data['pixLimit']
339    cutoff = data['cutoff']
340    if len(ring) < 5:
341        print 'not enough inner ring points for ellipse'
342        return False
343       
344    #fit start points on inner ring
345    data['ellipses'] = []
346    outE = FitRing(ring,True)
347    if outE:
348        print 'start ellipse:',outE
349        ellipse = outE
350    else:
351        return False
352       
353    #setup 360 points on that ring for "good" fit
354    Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)
355    if Ring:
356        ellipse = FitRing(Ring,delt)
357        Ring,delt = makeRing(1.0,ellipse,pixLimit,cutoff,scalex,scaley,self.ImageZ)    #do again
358        ellipse = FitRing(Ring,delt)
359    else:
360        print '1st ring not sufficiently complete to proceed'
361        return False
362    print 'inner ring:',ellipse
363    data['center'] = copy.copy(ellipse[0])           #not right!! (but useful for now)
364    data['ellipses'].append(ellipse[:]+('r',))
365    G2plt.PlotImage(self,newImage=True)
366   
367    #setup for calibration
368    data['rings'] = []
369    data['ellipses'] = []
370    if not data['calibrant']:
371        print 'no calibration material selected'
372        return True
373   
374    skip = data['calibskip']
375    dmin = data['calibdmin']
376    Bravais,cell = calFile.Calibrants[data['calibrant']][:2]
377    A = G2lat.cell2A(cell)
378    wave = data['wavelength']
379    cent = data['center']
380    elcent,phi,radii = ellipse
381    HKL = G2lat.GenHBravais(dmin,Bravais,A)[skip:]
382    dsp = HKL[0][3]
383    tth = 2.0*asind(wave/(2.*dsp))
384    ttth = tand(tth)
385    data['distance'] = dist = calcDist(radii,tth)
386    radius = dist*tand(tth)
387    zdis,cosB = calcZdisCosB(radius,tth,radii)
388    cent1 = []
389    cent2 = []
390    xSum = 0
391    ySum = 0
392    zxSum = 0
393    zySum = 0
394    phiSum = 0
395    tiltSum = 0
396    distSum = 0
397    Zsum = 0
398    for i,H in enumerate(HKL):
399        dsp = H[3]
400        tth = 2.0*asind(0.5*wave/dsp)
401        stth = sind(tth)
402        ctth = cosd(tth)
403        ttth = tand(tth)
404        radius = dist*ttth
405        elcent,phi,radii = ellipse
406        radii[1] = dist*stth*ctth*cosB/(cosB**2-stth**2)
407        radii[0] = math.sqrt(radii[1]*radius*cosB)
408        zdis,cosB = calcZdisCosB(radius,tth,radii)
409        zsinp = zdis*sind(phi)
410        zcosp = zdis*cosd(phi)
411        cent = data['center']
412        elcent = [cent[0]+zsinp,cent[1]-zcosp]
413        ratio = radii[1]/radii[0]
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            newellipse = FitRing(Ring,delt)
419            elcent,phi,radii = newellipse               
420            if abs(phi) > 45. and phi < 0.:
421                phi += 180.
422            dist = calcDist(radii,tth)
423            distR = 1.-dist/data['distance']
424            if abs(distR) > 0.1:
425                print dsp,dist,data['distance'],distR,len(Ring),delt
426                break
427            if distR > 0.001:
428                print 'Wavelength too large?'
429            elif distR < -0.001:
430                print 'Wavelength too small?'
431            else:
432                ellipse = newellipse
433            zdis,cosB = calcZdisCosB(radius,tth,radii)
434            Tilt = acosd(cosB)          # 0 <= tilt <= 90
435            zsinp = zdis*sind(ellipse[1])
436            zcosp = zdis*cosd(ellipse[1])
437            cent1.append(np.array([elcent[0]+zsinp,elcent[1]-zcosp]))
438            cent2.append(np.array([elcent[0]-zsinp,elcent[1]+zcosp]))
439            if i:
440                d1 = cent1[-1]-cent1[-2]        #get shift of 2 possible center solutions
441                d2 = cent2[-1]-cent2[-2]
442                if np.dot(d2,d2) > np.dot(d1,d1):  #right solution is the larger shift
443                    data['center'] = cent1[-1]
444                else:
445                    data['center'] = cent2[-1]
446                Zsum += numZ
447                phiSum += numZ*phi
448                distSum += numZ*dist
449                xSum += numZ*data['center'][0]
450                ySum += numZ*data['center'][1]
451                tiltSum += numZ*abs(Tilt)
452                if not np.all(checkEllipse(Zsum,distSum,xSum,ySum,dist,data['center'][0],data['center'][1])):
453                    print 'Bad fit for ring # %i. Try reducing Pixel search range'%(i) 
454            cent = data['center']
455            print ('for ring # %2i @ d-space %.4f: dist %.3f rotate %6.2f tilt %6.2f Xcent %.3f Ycent %.3f Npts %d' 
456                %(i,dsp,dist,phi-90.,Tilt,cent[0],cent[1],numZ))
457            data['ellipses'].append(copy.deepcopy(ellipse+('r',)))
458        else:
459            break
460    G2plt.PlotImage(self,newImage=True)
461    fullSize = len(self.ImageZ)/scalex
462    if 2*radii[1] < .9*fullSize:
463        print 'Are all usable rings (>25% visible) used? Try reducing Min ring I/Ib'
464    if not Zsum:
465        print 'Only one ring fitted. Check your wavelength.'
466        return False
467    cent = data['center'] = [xSum/Zsum,ySum/Zsum]
468    wave = data['wavelength']
469    dist = data['distance'] = distSum/Zsum
470   
471    #possible error if no. of rings < 3! Might need trap here
472    d1 = cent1[-1]-cent1[1]             #compare last ring to 2nd ring
473    d2 = cent2[-1]-cent2[1]
474    Zsign = 1
475    len1 = math.sqrt(np.dot(d1,d1))
476    len2 = math.sqrt(np.dot(d2,d2))
477    t1 = d1/len1
478    t2 = d2/len2
479    if len2 > len1:
480        if -135. < atan2d(t2[1],t2[0]) < 45.:
481            Zsign = -1
482    else:
483        if -135. < atan2d(t1[1],t1[0]) < 45.:
484            Zsign = -1
485   
486    tilt = data['tilt'] = Zsign*tiltSum/Zsum
487    phi = data['rotation'] = phiSum/Zsum
488    rings = np.concatenate((data['rings']),axis=0)
489    p0 = [dist,cent[0],cent[1],phi,tilt]
490    result,newWave = FitDetector(rings,p0,wave)
491    print 'Suggested new wavelength = ',('%.5f')%(newWave),' (not reliable if distance > 2m)'
492    data['distance'] = result[0]
493    data['center'] = result[1:3]
494    data['rotation'] = np.mod(result[3],360.0)
495    data['tilt'] = result[4]
496    N = len(data['ellipses'])
497    data['ellipses'] = []           #clear away individual ellipse fits
498    for H in HKL[:N]:
499        ellipse = GetEllipse(H[3],data)
500        data['ellipses'].append(copy.deepcopy(ellipse+('b',)))
501    print 'calibration time = ',time.time()-time0
502    G2plt.PlotImage(self,newImage=True)       
503    return True
504   
505def Make2ThetaAzimuthMap(data,masks,iLim,jLim):
506    import numpy.ma as ma
507    import polymask as pm
508    #transforms 2D image from x,y space to 2-theta,azimuth space based on detector orientation
509    pixelSize = data['pixelSize']
510    scalex = pixelSize[0]/1000.
511    scaley = pixelSize[1]/1000.
512   
513    tay,tax = np.mgrid[iLim[0]+0.5:iLim[1]+.5,jLim[0]+.5:jLim[1]+.5]         #bin centers not corners
514    tax = np.asfarray(tax*scalex,dtype=np.float32)
515    tay = np.asfarray(tay*scaley,dtype=np.float32)
516    nI = iLim[1]-iLim[0]
517    nJ = jLim[1]-jLim[0]
518    #make position masks here
519    spots = masks['Points']
520    polygons = masks['Polygons']
521    tam = ma.make_mask_none((iLim[1]-iLim[0],jLim[1]-jLim[0]))
522    for polygon in polygons:
523        if polygon:
524            tamp = ma.make_mask_none((nI*nJ))
525            tam = ma.mask_or(tam.flatten(),ma.make_mask(pm.polymask(nI*nJ,
526                tax.flatten(),tay.flatten(),len(polygon),polygon,tamp)))
527    if tam.shape: tam = np.reshape(tam,(nI,nJ))
528    for X,Y,diam in spots:
529        tam = ma.mask_or(tam,ma.getmask(ma.masked_less((tax-X)**2+(tay-Y)**2,(diam/2.)**2)))
530    return np.array(GetTthAzm(tax,tay,data)),tam           #2-theta & azimuth arrays & position mask
531
532def Fill2ThetaAzimuthMap(masks,TA,tam,image):
533    import numpy.ma as ma
534    Zlim = masks['Thresholds'][1]
535    rings = masks['Rings']
536    arcs = masks['Arcs']
537    TA = np.dstack((ma.getdata(TA[1]),ma.getdata(TA[0])))    #azimuth, 2-theta
538    tax,tay = np.dsplit(TA,2)    #azimuth, 2-theta
539    for tth,thick in rings:
540        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.)))
541    for tth,azm,thick in arcs:
542        tam = ma.mask_or(tam.flatten(),ma.getmask(ma.masked_inside(tay.flatten(),max(0.01,tth-thick/2.),tth+thick/2.))* \
543            ma.getmask(ma.masked_inside(tax.flatten(),azm[0],azm[1])))
544    taz = ma.masked_outside(image.flatten(),int(Zlim[0]),Zlim[1])
545    tam = ma.mask_or(tam.flatten(),ma.getmask(taz))
546    tax = ma.compressed(ma.array(tax.flatten(),mask=tam))
547    tay = ma.compressed(ma.array(tay.flatten(),mask=tam))
548    taz = ma.compressed(ma.array(taz.flatten(),mask=tam))
549    del(tam)
550    return tax,tay,taz
551   
552def ImageIntegrate(image,data,masks):
553    import histogram2d as h2d
554    print 'Begin image integration'
555    LUtth = data['IOtth']
556    if data['fullIntegrate']:
557        LRazm = [0,360]
558    else:
559        LRazm = data['LRazimuth']
560    numAzms = data['outAzimuths']
561    numChans = data['outChannels']
562    Dtth = (LUtth[1]-LUtth[0])/numChans
563    Dazm = (LRazm[1]-LRazm[0])/numAzms
564    NST = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
565    H0 = np.zeros(shape=(numAzms,numChans),order='F',dtype=np.float32)
566    imageN = len(image)
567    Nx,Ny = data['size']
568    nXBlks = (Nx-1)/1024+1
569    nYBlks = (Ny-1)/1024+1
570    Nup = nXBlks*nYBlks*3+3
571    dlg = wx.ProgressDialog("Elapsed time","2D image integration",Nup,
572        style = wx.PD_ELAPSED_TIME|wx.PD_AUTO_HIDE)
573    try:
574        t0 = time.time()
575        Nup = 0
576        dlg.Update(Nup)
577        for iBlk in range(nYBlks):
578            iBeg = iBlk*1024
579            iFin = min(iBeg+1024,Ny)
580            for jBlk in range(nXBlks):
581                jBeg = jBlk*1024
582                jFin = min(jBeg+1024,Nx)               
583                print 'Process map block:',iBlk,jBlk,' limits:',iBeg,iFin,jBeg,jFin
584                TA,tam = Make2ThetaAzimuthMap(data,masks,(iBeg,iFin),(jBeg,jFin))           #2-theta & azimuth arrays & create position mask
585                Nup += 1
586                dlg.Update(Nup)
587                Block = image[iBeg:iFin,jBeg:jFin]
588                tax,tay,taz = Fill2ThetaAzimuthMap(masks,TA,tam,Block)    #and apply masks
589                del TA,tam
590                Nup += 1
591                dlg.Update(Nup)
592                tax = np.where(tax > LRazm[1],tax-360.,tax)                 #put azm inside limits if possible
593                tax = np.where(tax < LRazm[0],tax+360.,tax)
594                NST,H0 = h2d.histogram2d(len(tax),tax,tay,taz,numAzms,numChans,LRazm,LUtth,Dazm,Dtth,NST,H0)
595                print 'block done'
596                del tax,tay,taz
597                Nup += 1
598                dlg.Update(Nup)
599        NST = np.array(NST)
600        H0 = np.divide(H0,NST)
601        H0 = np.nan_to_num(H0)
602        del NST
603        if Dtth:
604            H2 = [tth for tth in np.linspace(LUtth[0],LUtth[1],numChans+1)]
605        else:
606            H2 = LUtth
607        if Dazm:       
608            H1 = np.array([azm for azm in np.linspace(LRazm[0],LRazm[1],numAzms+1)])
609        else:
610            H1 = LRazm
611        Nup += 1
612        dlg.Update(Nup)
613        t1 = time.time()
614    finally:
615        dlg.Destroy()
616    print 'Integration complete'
617    print "Elapsed time:","%8.3f"%(t1-t0), "s"
618    return H0,H1,H2
Note: See TracBrowser for help on using the repository browser.