Granetz, R.S.

Wang, L.

DSpace@MIT

PFC/JA-92-22A Series Expansion Method inThree Dimensional TomographyLing Wang and Robert S. GranetzPlasma Fusion CenterMassachusetts Institute of TechnologyCambridge, MA 02139August, 1992Submitted to Journal of Optical Society of America A.This work was supported by the U. S. Department of Energy Contract No. DE-AC02-78ET51013. Reproduction, translation, publication, use and disposal, in whole or in partby or for the United States government is permitted.A Series Expansion Method in Three DimensionalTomographyLing Wang and Robert S. GranetzPlasma Fusion CenterMassachusetts Institute of TechnologyCambridge, MA02139United States of AmericaAbbreviated title:A Series Expansion Method in 3-D TomographyKey words. X-ray transform, three dimensional tomography, spherical harmonic series ex-pansion, Euler angles.AMS(MOS) subject classifications. 44A15,42C15ABSTRACTA series expansion method in three dimensional tomography is presented. While the functionto be reconstructed is expanded in terms of spherical harmonics and orthogonal polynomials,the projection of the function is expanded in terms of Euler angles in the angular part andorthogonal polynomials in the radial part. The computer simulation of the reconstructiondemonstrates the applicability of the method.11. IntroductionIn two dimensional tomography since the pioneering work of Cormack [1], the method ofseries expansion in terms of Fourier harmonics is well known to have great advantages whenthe number of views is limited. Deans has given an extensive review of the method and itsapplications in many branches of science and engineering [2]. It is especial useful in somephysical problems, plasma fusion for example, where identifying dominant low order harmoniccomponents is essential [3]. It is conceivable that this method will play an even more importantrole in three dimensions not only because stacking of two dimensional slices becomes difficult ifnot impossible due to the constraints on viewing geometry but also because series expansion byspherical harmonics is a useful tool in many physical problems. There have been some importantdevelopments in the series expansion method in three dimensions. As in two dimensions, theanalytical expressions in both real and projection spaces have been found for certain classesof radial orthogonal polynomial basis functions [4, 5, 6] and the formulas for forward andinverse transforms of general radial functions have also been derived [7]. The choice of Eulerangles not only directly represents viewing geometry but also makes the angular basis functionsexplicitly analytical [7]. However, the above achievements have been mostly focused on thetheoretical side and their numerical application has just started [8, 9]. Many practical andimportant aspects have yet to be explored and this is the motivation of the present paper. Themajor difference between the proposed method of this work and that of references [8] and [9]is that while choosing the same basis functions in the radial part, we use in the projectionspace the Euler angle representation in the angular part. This choice enables us to perform theleast squares fit directly in the projection space and therefore no additional repeated Fouriertransforms are required. This is especially useful not only in terms of saving computing time2but even more importantly, it is suitable to the cases where one has only a limited number ofviews, for which discrete Fourier transforms are not appropriate. We will demonstrate for thefirst time the applicability of this method by numerical reconstruction of a test function. Thedetailed algorithms and procedures are given in Section 2, the reconstruction is performed inSection 3 and Section 4 draws the conclusion.2 Euler angle representation and series expansionThree dimensional tomography starts from the X-ray transform of a function f(f) definedby [10]g = f(ff + s)ds (1)where i is a unit vector along the integration path and fl= F- (F i)i. Hence fT L i, where F isan arbitrary point on the integration path. In tomographic applications, g(, i) is a measuredprojection and f(F) is the three dimensional image of interest to be reconstructed.As shown in Figure 1, for a given integration path along a unit vector i, we can alwaysrotate the coordinate system in such a way that the new coordinate axis can be aligned asi' 11 i, V' 11 P. (2)It is well known that this rotation can be parameterized by the Euler angles, (a,#, ) as shownin Figure 1. We adopt in this paper the "y convention" of Reference [11]. Using the Eulerangles, a vector E in XYZ system is related to the same vector S' in X'Y'Z' system byV '= Azf, :F= A-iz' = A TE' (3)3where the matrix A is given byCos a cos # cos - - sin a sin y sin a cos 0 cos y + cos a sin - - sin 0 cos yA = - cos ssiny - sin a cos y - sin a cos # sin y + cos a cos y sin 0 sin 7 (4)cos a sin sin a sin3 cos )and0 < a<27, 0 S/ r, 0 y <27r. (5)We can therefore always represent i and #in terms of Euler angles (a, 0, y) and p Ip1. UsingEquation (3) their Euler angles are found to becos 0 = iz (6)cos a - (7)sin 3sin a = (8)sinsin y = . (9)sin,#cos- -psin/# + psin a cos3 (10)when / # 0 and p $ 0. For / = 0 we can set y = 0 and a is determined bycos a = p, sin a = py. (11)while for p = 0, -y is not needed and (a, 0) are from Equations (6) to (8).We may now rewrite Equation (1), the x-ray transform of f(r, 0, #), asg(, t) = g(p, a, /, 7) = f(r, 0, 4)ds (12)where L(p, a,3, 7) represents the integration path as shown in Figure 1.4For any given function f(r, 0, 4) we can always make the following expansion in terms ofspherical harmonics:f(r, ,4) = > (fl.(r) Sfm(, 0)+ flm(r) ( (13)1=0 m=0where ff,(r) are the radial expansion components and the real spherical harmonics are simplySIm(0, 4) Re(Yim(0, )) = NimP m (cos 0) cos m4 (14)Sf.(0, Im(Ym(O,4)) = NimPm (cos 0) sin m (15)where Yim (9, 0) are the well known spherical harmonics defined byYim(0,q$) NmPim(cos 0)e im (16)where(21 + 1)(l - m)!Nim3 (1)" 47r(l + m)!andP(x) = (1 - X2)M/2 d 2I+ _ 1)' (18)is the associated Legendre function which can be rewritten in a more useful form for numericalevaluation as [12])PIM"(cos 0) = (l + m)! sin 0 M P(Mm)(Cos ) (19)where P('") (x) is the Jacobi polynomial of degree n and type (v, p) [13].The radial components ff;"(r) also have to be expanded in terms of a complete set of basisfunctions {h3 (r)}'o0 (to be specified later):ff;s(r) = ZA"'h (r), (20)j=05where Ajf are the expansion coefficients to be determined. Equation (13) then becomesf(r,,4) = E (Aj.Sf.(0, #) + AjimSIm(9, #))hj3 (r)j=0 1=0 m=0We define Wj', (a,3, y, p) as follows and it has been proven that [7]Wi',(a, 3, p) ] dsS,"(O, #)h 3 (r) = V,', (a, 0, -y)gjim,(p)J(PA#OP,-) MI=O(21)(22)whereVinmm,(a, 0, -) = m, [cos(ma + m'(y ++(-1)m'(cos(ma - m'(y + 7r/2))dd,_m,(#)]VImm,(a,,y) = bm' sin(ma + m'(-y + 7r/2))d(')(3)+(-1)m'(sin(ma - m'(y + 7r/2))d$,- ,(/)](23)(24)where1/26M =m = 0m 540(25)andM I-M M I+M It -MM'+M)Cos M, (cos,3)2) (26)(for m' > m)and(for m > m'), (27)and gjIm(p) are the radial basis functions in the projection space (a, 3, y, p), obtained fromhjl(r) by [7]gjIm(p) = Nim(1 + (-1)1+m)J rdh i " (F I ( 1 (28)6d~(l ) 0)= I!sin )dmm,(#3) = (-1)"-'~"d,()Taking the x-ray transform of f(r, 0, 4) in Equation (21) and employing Equation (22) we cantherefore expand the projection g(a,#3, y,p) defined in Equation (12) aso0 00g(a,3, -y, p) = (A imW-m(Q,/,y, p)+ A-imWim(a, ,y,p)). (29)j=O 1=0 m=OSo far we have not specified h3 (r) yet. In principle once it is given, by a certain physical modelfor example, one can always use Equation (28) to get gjim(p) at least by numerical integration.In this paper, we choose the orthogonal polynomialshji(r) = rl(1 - r2)p,(+1/2,1)(1 - 2r 2), (30)because its gjl,(p) can be carried out analytically [4, 6], namely{ Gjimp m (1 _ P2)3/2p('' 3 /2 ) (1 - 2p 2) I + m = evengl (p) =1m/+ (31)0 l + m = odd,whereGjlm (-1)'(j + 1)((l2- + j)! V(21 + 1)(1 - m)!(I + m)! (32)21+'F('- + j+ 2 + ) -mf)!('+-)!It should be noted that although the above formulas look complicated, the only special functionwhich appears in Equations (19), (30), (31) and (26) is the Jacobi polynomial PV'M"')(x). Itsevaluation can be performed by using the stable recursion formula suggested by Abramowitzand Stegun (Eq. 22.7.1 of [13]) that2(n + 1)(n + v + it + 1)(2n + v + p)P("')(x)= [(2n + v + p + 1)(v 2 _ 12) + (2n + v + p)(2n + v + p + 1)(2n+ nup +2)]P )()-2(n + v)(n + p)(2n + v + y + 2)P,(')(x) (33)with0')(x) = 1, and P( ')(x) = [v - p + (v + p + 2)x]/2. (34)73 ReconstructionWe take the following function with an ellipsoidal iso-valued surface as our test functionf'(x, y, z) = exp - + + . ) (35)where a = 0.2 b = 0.3 c = 0.4and its projection gt(a, /, y, p) is from the numerical integration, Equation (12).In our simulation, 21 x 21 element array detectors are used. We have two detection setups.There are four arrays in the first setup. They are located symmetrically at the four corners ofan equal-sided tetrahedron which is centered at the origin of a coordinate system. The straightline through the center and each corner of the tetrahedron passes through the center of eacharray and is normal to it as well. Each array receives signals through a pin hole aperture. Themaximum chord distance, p, is about one. In the second case, the setup is the same except fourmore arrays are added. The eight arrays are now located symmetrically at the eight corners ofa cube.The number of series expansion terms in Equations (21) and (29) is infinite, however inpractical applications the number of measured projections is finite, and therefore the series hasto be truncated. After a proper truncation we haveJ L If(r, 0, 4) = (A. Si.(0, 0) + AjimS8.(0, 4))h 1(r) (36)j=O 1=0 m=OandJ L Ig(ak, Ok,,p) = Z (A ImWjim(ak, 3 , Yk, pk) + AjlmW mk, (ak, ,k, Pk)) (37)j=0 1=0 m=0k = 1, 2, -,N.where N is the total number of measurements. Equations (37) are just linear equations for theunknowns Ajf and this system can be solved by the method of least squares. If the samples are8equally or approximately equally spaced then the Nyquist principle can be used to estimate themaximum number of expansion harmonics which may be resolved [3]. However, this is not thecase in our sample space, at least in the angular subspace (a,,3, -y). How to directly determinethe truncation in this case is still an open question and therefore numerical reconstruction testswere conducted in lieu of this. We have found by trial and error thatL = 3 J = 9 (for 4 arrays) and L = 6 J = 9 (for 8 arrays) (38)gives optimal reconstruction results. We'd like to point out that the maximum radial number,J = 9, from numerical tests does agree with what the Nyquist principle requires; namely thatthe number of samples should be at least twice the frequency of the highest spatial harmonicsin the radial direction. In our detector geometry (21 samples in roughly radial direction), thechord radial spacing is uniform enough to apply the Nyquist theorem. (This points out anotherone of the reasons for deriving this analytical expansion technique - the radial and angularparts of the transform problem are separable.)After solving Equation (37) for the expansion coefficients, A', the reconstructed functionis then given by Equation (36). In Figure 2 (4 array case) and Figure 3 (8 array case) we showa few three dimensional iso-valued surfaces of the reconstructed function in comparison withthat of the test function. As one can see while the reconstruction is reasonably good for thecase with four arrays it is greatly improved when eight arrays are used. The slight deformationof the surfaces in Figure 2 is due to the fact that, although the test function Equation (35)looks simple, it does contain an infinite number of spherical harmonics. This can be seen byrewriting ft(x, y, z) in spherical coordinates asr 2 (cos #sin0) 2  (sin 4 sin 0)2 (cos 0)2 (39)ft(rq) ,= exp -r a2 + + (39)9and making the expansion e- = 1- z+z 2/2+.... It is immediately clear that higher harmonicsare there, especially at large radii. The slight deformation is an "aliasing" effect, which showsup as expected as a result of not having a sufficient number of samples.4 ConclusionWe have presented a series expansion method in which the Euler angles and orthogonalpolynomials were employed. We have also show the detail algorithms and procedures of themethod. The applicability is demonstrated through a successful reconstruction of a test func-tion. We hope to see the use of this method to a real problem in the field of plasma physics,for example laser fusion [14, 15], in the near future.AcknowledgementsWe wish to thank Professor Ian Hutchinson for his valuable suggestions.10References[1] A. M. Cormack, Representation of a function by its line integrals, with some radiologicalapplications, J. Appl. Phys., 10 (1964), pp. 2908-2913.[2] S. R. Deans, A unified Radon inversion formula, J. Math. Phys., 19 (1978), pp. 2346-2349.[3] R. S. Granetz, P. Smeulders, X-ray tomography on Jet, Nuclear Fusion, 28 (1988), pp.457-476.[4] P. Maass, The X-ray transform: singular value decomposition and resolution, Inverse Prob-lems 3 (1987), pp. 727-741.[5] S. H. Izen, A series inversion for the X-ray transform in n dimensions, Inverse Problems,4 (1988), pp. 725-748.[6] S. H. Izen, Inversion of the k-plane transform by orthogonal function series expansions,Inverse Problems, 4 (1988), pp. 725-748.[7] L. Wang The X-ray transform and its inversion for the series expansion basis functions inthree dimensional tomography, to appear on SIAM, J. Math. Anal., 52 (1992). (Oct. Issue5)[8] S. H. Izen, An application for a limited solid angle x-ray transform, Contemporary Math-ematics, 113 (1990), pp. 151-170.[9] S. H. Izen, Inversion of the x-ray transform from data in a limited angular range, in SignalProcessing, Part 2, IMA Volumes in Mathematics and its Applications #23, SpringerVerlag, 1990.11[10] D. C. Solmon, The X-ray transform, J. Math. Anal. Appl., 56 (1976), pp. 61-83.[11] H. Goldstein, Classical mechanics, Addison-Wesley, Mass, 1980.[12] L. C. Biedenharn and J. D. Louck, Encyclopedia of Mathematics and its applications. Vol.8:Angular Momentum in Quantum Physics, Adison-Wesley, Mass., 1981.[13] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover Publica-tions, New York, 1970.[14] M. Nakai, et al, Development of x-ray emission computed tomography for ICF research,Rev. Sci. Instrum. 61 (1990) pp. 2783-2785.[15] Y. W. Chen, et al, Three-dimensional imaging of laser imploded targets, J. Appl. Phys. 68(1990) pp. 1483-1488.12zz-- *- - ---. 4.. - - .*... - . e-. --. * e. .- .. . -.--. - .... *#, 1Jy)xThe Euler angle representation of the integrationpath L(p,a,3,y).13xFIG. 1.L (, PReconstructed function Test functionA. f(r,6, 4) = 0.4B. f(r, 6, #) = 0.2surfacesurfaceFIG. 2. 3-D iso-valued surfaces of the reconstructed function(left) and test function (right) when 4 arrays are used.14(DReconstructed functionA. f(r, 0, 4) = 0.4B. f(r, 0, 4) = 0.2surfacesurfaceFIG. 3. 3-D iso-valued surfaces of the reconstructed function(left) and test function (right) when 8 arrays are used.15Test function

A series expansion method in three dimensional tomography

http://dspace.mit.edu/bitstream/1721.1/95153/1/92ja022_full.pdf

A series expansion method in three dimensional tomography

Abstract

Similar works

Full text

Available Versions

DSpace@MIT