Vol. 72, No. 12/December 1982/J. Opt. Soc. Am. R. A. Niland 1677 Maximum reliable information obtainable by tomography R. A. Niland School of Physics, University of Sydney, New South Wales 2006, Australia Received March 1, 1982; revised manuscript received June 30, 1982 m equispaced views of a two-dimensional 2 object yield precisely M real numbers characterizing the object. These are generalized moments of the object and are free of aliasing contamination. A reconstruction that matches these moments is positively constrained and maximizes entropy can be generated. Certain types of a priori knowledge can be incorporated, and the reconstruction is not overly sensitive to noise on the data. INTRODUCTION This study was motivated by the experimental work of Meyers and Levine,1 who used 4-viewtomography in an attempt to obtain spatially resolved spectral measurements in a TORMAC plasma. A view is defined as the collection of all line integrals through the object along lines parallel to a given line. Investigation of the various algebraic reconstruction techniques led to the disturbing result that the reconstructed intensity distribution depended on the choice of algorithm, the number of iterations, and the intensity distribution of the object, in a generally unpredictable way. It became of interest only to those moments of the projections known to be free of aliasing contamination. Q TABLE REPRESENTATION OF AN OBJECT Suppose the object 4(r, 0) is contained within a circle of unit radius. We expand 4'in terms set of functions rheilo, of the complete nonorthogonal 1 =O. to answer the question: Given good quality projection data, 4 is fully described (1) k= 0,1, 2,..., 1,,2.... precisely what can be asserted about the object? The answer, which is shown here, is that m equispaced views determine M 2 independent real numbers, generalized moments of the object. A positively constrained function with the same moments and maximizing entropy can be generated and would seem to be so that a good candidate for the optimum reconstruction. This paper is related to other work as follows: It has been This table is a sampled form of the Fourier-Mellin transform known since 1917 (Radon, see Ref. 2) that the complete col- lection of line integrals through an object ip(r, 0) contained within the unit circle holds all the information required to reconstruct A' exactly, and an explicit reconstruction formula can be given that, however, is numerically ill conditioned. Cormack 3 independently solved the same problem by expanding 4'in orthogonal functions, which allowed a reasonable numerical solution. There are now many computer algorithms for reconstructing 4'from its projections. They fall into the general classes of algebraic and Fourier (or convolution) methods. If there are many projections available, as in medical tomography, the algorithms produce good results in general and differ only in speed, cost, or complexity. For few views, however, the results are variable, as noted above. Simple arguments 4 indicate a resolution cell size of diameter about 1/m, where m is the number of views. A more elaborate treatments yields about the same result and further demonstrates that the mechanism that limits resolution is angular aliasing. Their results are made sharp by the present paper. by a table of complex generalized moments 1 Qlk = fordr 2fr (2) fod~rkeilop. 4,which has a deep connection with tomography. We shall take 4'real so that only the half-table I > 0 need be considered. The Q table could be written equally well in of terms of functions orthogonal on the unit circle, and this will be done later for numerical calculations. However, the fun- damental results are clearer in the nonorthogonal representation. GENERATION OF THE Q TABLE FROM COMPLETE PROJECTION DATA Suppose we have a complete set of line integrals of 4': G(t, 0) 0 < t < 1, 0 < 0 < 2-r, where G(t, 4) is the integral alongy cos b - x sin / = t. By Radon's theorem 2 G establishes 4'uniquely. To connect G to the Q table we need two preliminary results whose proofs may be found in Refs. 3 and 4: 4 = fl(r)ei1 0 has a projection in the 0 direction of (Alf 1)(t) exp[il(o + 7r/2)], where AlfV(r)] = 2 f I fl (r) Ti(t/r)rdr (3) Minerbo 6 uses a maximum entropy method for recon- struction from few views. The maximum entropy solution is constrained to match the actual projections. By contrast, the present paper matches the maximum entropy solution 0030-3941/82/121677-06$01.00 and T1 (x) is a Chebyshev polynomial. We also require a result connecting moments of a function with moments of its projection: © 1982 Optical Society of America 1678 J. Opt. Soc. Am./Vol. 72, No. 12/December 1982 Sp1 Alf tkdt l R. A. Niland 1 = to retrieve functions of r from their high moments alone. This rPk o fl(r)rkrdr, reflects the well-known ill condition of the Radon inversion formula, which involves a derivative of the data. where r7(k + 1) = h Notice that r1k (5) 2(k +1+) (k -+2) Suppose that G (t, 4) is known only at m discrete equispaced values of 0 such that vanishes if k <1 and k + 1is even. Now expand the object as 1 4(r, 0) = I ao(r) + 2 c a [ai(r)cos 10 + 1=1 TV= w7r/m, bi(r)sin 101. (6) Its projection is G(t, 0) = -Aoao(r) + L Aal(r)cos I GENERATION OF THE Q TABLE FROM INCOMPLETE PROJECTION DATA v = 0, 1, ... , This is m-view tomography. we write *(t) = +2 L 12 2m 1. - Analogous to Eqs. (8) and (9) G(t, 0,) exp [ii 1( + 7r (10) using the discrete Fourier transform by necessity now, and + Albl(r)sin l +-2)] (7) Qlk*=-F1hk f dtthg,*(t). (11) The main result of this paper is that, for a certain set of 1, Qik = ffrdrdOrke"O k, Qlk* It can be seen as follows: = Qlk. 1 = 7r- rdrrh [al(r) + ibi(r)] gi*(t) = gj(t) + E gi,(t), [from Eq. (6)] (12) 1' = k Fik = 7r we k 4'° Ai(ai + ibi)tkdt where 1' ranges over the set of aliases of 1:1' = 2m ± 1, 4m + [from Eq. (4)] X thkdtg,(t ) 1.... [Proof: in Eq. (10) replace G by its inverse Fourier transform from Eq. (9).] [from Eq. (7)], (8) where gi (0 Inserting Eq. (12) into Eq. (11) yields Qlk* = Qlk + LA I' Pik fO =f 2 G+ =- l G(t, 0) exp[il(o + 7r/2)do. (9) = Qlk + Because of Eq. (5) not all the Qlk can be extracted from the complete projection data, as indicated in Fig. 1. However, according to Radon, those that can be obtained are sufficient to reconstruct 4 exactly. This follows from the Szasz-Muintz theorem, which says that Irk Ih = 0, k1 , k2,.... I is complete on [0, 1] provided that Z1 /kj = a. It is clear, however, that, as I increases, the reconstruction process becomes increasingly ill conditioned as one attempts 4 views I -Qo6 /-7 5 views I Qs8 dttkgl,(t) 6 views I Qolo = Qik E 1' PI'kQl'klrlk 1 =0,1,...I k = 1,1 +2,.... for m - 1, (13) This result is illustrated in Fig. 1 for several values of m, the number of views. In all, M2 real quantities can be obtained from the incomplete data. They lie in a triangular section of the Q table (Fig. 1) and are defined as Qi = ff rdrd0Oui, i = 1, 2 .... A For example, if m = 4, the 16 functions ui are '19 1, r2 , r4 , r6 , 2 M . (14) r cos 0, r sin 0, r3 cos 0, r 3 sin 0, r 5 cos 0, r 5 sin 0, r2 cos 20, r2 sin 20, H Cos 20, r4 sin 20, r 3 cos 30, r 3 sin 30. (15) If 1and k lie outside the range indicated in Eq. (13), then Q1k* 7- Qlk because not all the 1I'k = 0. For example (m = 4), inaccessible to tomography Q08* = Q08 + (P8 8/P08 )Q88, Fig. 1. Q table illustrating element retrievable by 4-, 5-, and 6-view tomography and those inaccessible to tomography. With 4 equi- and Qo8* differs from the true value by an unknown quantity that, since it reflects a component of ' with angular dependence ei80, cannot be measured by 4-view tomography. It is remarkable that a block of the Qlk remains totally free of this effect. spaced views, the elements in the triangle bounded by Q00, Qo6, can be retrieved. With 5 views the triangle extends to Qoo,Q08, etc. Elements on and to the left of the line Q20, Q64 Q33 Q44, cannot be ob- tained exactly by tomography with any finite number of views. R. A. Niland Vol. 72, No. 12/December 1982/J. Opt. Soc. Am. 1679 m2 and the procedure re- Until now essentially no assumptions have been made about Q'. If a priori knowledge is available, other Qlk* may be as- This is computed for i = 1, 2,... serted equal to their correct values Qlk and the number of parameters describing the object increased beyond M2 . An If Eq. (16) is discretized on a square grid with cell area equal to A, we have extreme example is if 4 is known to have cylindrical symmetry; then Qlk = 0 except for 1 = 0 and hence Qok* = Qok for all k, allowing an infinite number of parameters to be retrieved (in principle). This merely states that the Abel transform exists. In practice, noise would limit the highest usable k. Again, if 4'> 0, the constraint is reflected in the Q table: QOk > 0 for all k, and, in fact, these moments define a point peated to convergence. A uijfjA - Qi = 0, Jo where fj and uij are representative values of f and ui in cell j. Then Eq. (17) becomes Z uijfje8uiiA - Qi that lies in the convex moment body.7 More-complicated j. = (i = 1, 2,... M2), 0 (18) relations exist for I > 0. It is not clear how to incorporate this which is Bregman's method 9 for maximizing entropy of a information into the reconstruction process. Beyond these simple examples a priori knowledge is not treated further here. vector subject to linear equality constraints, and its converform, by Halley's method. Both are convergent. These methods and convergence results remain true for more sophisticated discretizations of Eq. (16). MAXIMUM ENTROPY RECONSTRUCTION Tomography with m equispaced views gives then M 2 real numbers pertaining to the object 4'. They are generalized moments of 4. For some applications these numbers could be used directly, e.g., to test a given hypothesis concerning 4'. In other cases an image of 4 may be desired. From now on suppose that 4 > 0. Because of its robust propertiess we will reconstruct a function f (r, 6), which is positively constrained, matches the observed Qi, and maximizes entropy: SS rdrdOflnf. f- is not reasonably diagonally dominant, as may be seen from standard treatments of linear Gauss Seidel. NEWTON METHOD written down explicitly and evaluated at moderate cost for small m: Jij = aFi1/Xj = ,ffrdrdOfuiuj. (19) m2 = exp E Xiui Then i=1 Jd = F with the m2 constraints of the form Fi(X) The convergence of Gauss Seidel will be slow if the ith equation is not a strong function of Xi-i.e., if the Jacobian The nonlinear system [Eq. (16) or (18)] was also solved by Newton's method by using the full Jacobian, which can be Standard variational techniques 6 show that f must be of the form f gence is guaranteed. Equation (18) may be solved for each i by Newton's method or, more efficiently, given the special rdrdtfui - Qi = 0, and i = 1, 2,... X(new)= m2 (16) determining the unknown parameters Xi. Note that f (r, 0) will generate a Q table matching that of 4'in a triangular block containing m 2 real values and making estimates of the other Qlk's based on the totality of reliable information. -rd, where Te[0, 1] permits a one-dimensional search for the min- imum of FTF in the Newton direction.1 0 In practice, with four views, 10 X 10 cells, and orthogonalized ui, the fixed choice r = 1 produced convergence in three to five iterations for a variety of artificial objects. REDUCTION OF DATA GAUSS SEIDEL ALGORITHM Two approaches to the solution of the nonlinear system [Eq. 2 (16)] will be described. If the number of equations M is small enough, Newton's method is probably superior. It is treated For practical computing, the ui form too skewed a basis, and it is desirable to orthogonalizethem over the disk. Equation (14) becomes in the next section. However, it is useful to know that the Qi= ffrrdO4ui nonlinear Gauss Seidel algorithm is convergent for this system and may be efficiently implemented. The ith equation of the nonlinear system [Eq. (16)] is solved for Xi. If we set Xi (new) = Xi + p, the equation becomes if rdrd~elluifui - Qi = 0, (17) which is a one-dimensional nonlinear equation in A whose root gives f(new) = ezuif. = dtvi(t)gi(t), (20) where if rdrduiz2j = 27rbij This ensures good initial condition for the Newton method 1680 J. Opt. Soc. Am./Vol. 72, No. 12/December 1982 R. A. Niland Table 1. Polynomials Relevant to 4-View Tomography ui 2i vi (s = 2t) l=0 11 2 3 4 5, 6 7, 8 9,10 = 2 11, 12 13, 14 1 =3 15, 16 r2 r4 r6 [2( 1)11/2 [6( 2r 2 - 1)11/2 [10( 6r 4 - 6r 2 + 1)]1/2 [14(20r 6 - 30r4 + 12r 2 - 1)11/2 [2( 1)11/2 [6( S2- 1)]1/2 [10( s4 - 3 2 + 1)11/2 [14(s 6 - 5s4 + 6s2 - 1)11/2 r r3 r5 2 ( r) 2[2( 3r3 -2r)]1/ 2 2[3(10r 5 - 12r3 + 3r)]1/2 2( 2[2( 2[3(s 5 r2 r 2 [6( r )J1/ [10(4r 4 - 3r 2 )]1/2 [6( [10( r3 2(2 r 3 )1/2 2(2 2 s) s 3 -2s)] 1 / 2 4s3 + 3s)]1/2 - 4 s S2 )]1/2 -32)11/2 S3)1/2 (if initial X = 0, J = 2rI) and moderately good condition for the extraction of moments from the data. For concreteness the various polynomials relevant to 4-view tomography are more than small changes in f and hence Q. So the inverse mapping Q - X is nearly singular for this f. In practice, the listed in Table 1. For I > 1, each entry ui, di represents two functions [as Eq. (15)], e.g., U5 = 2r cos 0, a6 = 2r sin 0. The di were formed by during the Newton process. It will manifest only if 4'has large zero or near-zero areas. In such critical cases a sure way to proceed is to add a constant offset b to q4(by altering Qok, all the Gram-Schmidt process from the ui but turn out also to be renormalized Zernike polynomials. 3 The formation of the vi is illustrated by if rdrd042Vd(3r 3 - d1 dt2Vr k), reconstruct, and study f(b) - b as b - 0, thus embedding the problem in a family parameterized by b. For large b the reconstruction may not be positive, and for b = 0 the search may fail. Such an embedding is computationally cheap because a good guess for X is available for each new value of b. Notice also that, if b , the problem can be solved explic- 2r)ei° = Q7+ iQ8 = r problem can be detected by monitoring the condition of J t3 _ r2t jgi(t) =f dt2v¶(8t3 - 4t)g1 (t) = dtv 7g1 (t) = r dtv7gi*(t). Note that V7= vg, etc., and it is convenient to define the vi as functions of s or t and not 0. The vi form a much better itly. NUMERICAL RESULTS --1 This section illustrates the performance of a straightforward implementation of the maximum entropy algorithm. For m views we take the set of alias-free functions orthogonalized over the unit disk as R 1 P cos 10, R1 sin 10, = 0,1, ... m-1 v = 0,1, .. . m conditioned basis than powers of t. (21) STABILITY OF THE RECONSTRUCTION where It is important to know that the reconstruction f will not change dramatically if 4 or Q is changed slightly. Overall we have a sequence of mappings 4 - Q - X - f from the object to the vector of its reliable moments to the maximum entropy solution. From numerical experiment and intuition, the map 4 - f is well behaved: f is a smoothed-out version of 4' matching it in a collection of moments, and there are no singularities in the map. Notice that it is projective (i.e., the map applied to f produces f) and the projective part resides in 4 - Q, which is also well behaved. Consider now the function Q(X), where Qi(X) = ff rdrdOuif. This is not invertible everywhere because, if at least one Xi is large in magnitude, it can be made larger without producing Rj"(r) = Aiv (-)s(l+ 2v-s)!rl+2-2s s=0 s!(l + v - s)!(v -s)! a Zernike polynomial, and the required normalizing factor Al is [2(2v + 1)]1/2for 1 = 0 and 2(1 + 2v + 1)1/2 otherwise. In passing note that the v polynomial corresponding to Riv(r) is A1" s=0 (-)5(1 + 2v-s)! s!(1 + 2v - 2s)! (2t)1+ 2p-2s. We take 8 equispaced views, so that there are 64 functions in Eq. (21). Integrals over the disk were done by a product trapezoidal Gaussian rule with 400 points: I JJ rdrd0p(r, 0) = - 1 20 20 F Z wip(ri, 0j), 20OJ=1i=1 Vol. 72, No. 12/December 1982/J. Opt. Soc. Am. R. A. Niland 1681 thogonality of the functions in Eq. (21). Given an artificial object 4(r, 0), the program generates the vector of moments Q, then taking as an initial guess X = 0 searches with Newton's method for a maximum entropy function possessingthe same moments. The Jacobian was used three times after each update. Although simple artificial objects were used, no a a priori knowledge was assumed. - - I- N ,1 ObjectA: - II I if r• 0.3 =0.2 if r>0.3. I X '= 1 The 8-view maximum entropy reconstruction is shown in Fig. 2(a), and details of convergence are shown in Table 2. The quantities given are (a) IQ - Q(X) 12, where X is the latest estimate of the maximum entropy parameters, and (b) the estimated condition of the Jacobian, as given by the subroutine DECOMP. 1 2 b Object B: 4= a = 0.5 = 1.0 0 x l Fig. 2. a, Slice of object A along any diameter and its 8-view maximum entropy reconstruction. b, Slice of object B along x axis and (= 0.2, 0.05, 0.005) if r • 0.8 if r > 0.8 2 2 0.25. if (x - O.5) - y Figure 2(b) shows the reconstruction for a = 0.2 along the slice y = 0. The progress of convergence is illustrated in Table 2. For low a it is clear that convergence is slower and also that the final condition of the Jacobian becomes large. The overshoots at the discontinuities also becomemore pronounced, its reconstruction. being about 50% for a = 0.005. where CONCLUSION -1)7r/10, j = 1, 2,... 20, ri is the ith node of a 20-point Gaussian rule on [0, 1] with weight r, and wi is the corresponding weight." Oj Coarser discretizations did not accurately reproduce the or- A finite number m of equispaced views of an object 4(r, 0) confined to the unit circle clearly does not contain all the information contained in the object. This paper explicitly identifies the features of the object that are not lost in the projections. In fact, precisely m2 generalized moments of ' can be extracted from the projections. Table 2. Convergence of Maximum Entropy Algorithm for 4 Objects Iteration 1 Error 0.2877 2.78 -2 Object A Condition 1.00 1.01 -2 4.64 -3 2 1.10 -4 1.05 -5 2.03 -11 1.63 -15 0.1728 1.44 -2 Object B a = 0.05 1.00 4.81 -2 2.08-3 3.90 1.17 -6 3 Object B a = 0.2 4.41 -5 4.15 -6 9.53 -12 4.00 -15 1.00 1.04-2 5.68-3 4.93 4.64 -7 5.55 0.2271 2.46 -2 Object B a = 0.005 4.86 -4 1.67 -4 1.52 -6 1.40 -7 6.99 3.17 -13 8.87 -4 4.03 -4 7.85 2.33-4 23.45 1.55-8 4 1.00 1.28-2 7.38-3 7.30 -5 7.49 0.2451 2.85 -2 2.75 -5 1.23 -5 40.47 7.07 -6 33.44 7.90 -7 172.89 3.00 -7 1.38-7 5 3.31 -9 210.81 2.90 -10 3.16 -11 6 1.09 -14 395.84 1682 J. Opt. Soc. Am./Vol. 72, No. 12/December 1982 In addition, for the case 6> 0, a convergent algorithm can be given to generate an object f that has the same m2 moments as if and maximizes entropy, i.e., is maximally noncommittal about those features of 4,that cannot be obtained from the projections. Some possible further lines of inquiry will now be indi- cated: (a) The algorithm that employs only m2 moments is conservative, and it is of interest to develop strategies to expand this number on the basis of either a priori knowledgeor judicious risk taking based on estimates of the unknown Ql'k in Eq. (13). (b) As mentioned in the introduction, the chief use for maximum entropy reconstruction would seem to be for small m. However, if m is large it may be possible to produce a conventional reconstruction and then perturb it slightlyto the maximum entropy reconstruction, with favorable effects on noise and artifacts. (c) Finally, the effects of finite sampling within a view need to be studied. R. A. Niland REFERENCES 1. B. R. Meyers and M. A. Levine, "Two dimensional line emission reconstruction as a plasma diagnostic," Rev. Sci. Instrum. 49, 610-614 (1978). 2. L. A. Shepp and J. B. Kruskal, "Computerized tomography: the new medical x-ray technology," Am. Math. Mon. 85, 420-438 (1978). 3. A. M. Cormack, "Representation of a function by its line integrals, with some radiological applications," J. Appl. Phys. 34,2722-2727 (1964); 35, 2908-2913 (1964). 4. R. A. Niland, "Some aspects of plasma tomography," Lawrence Berkeley Lab. Rep. LBL-9373 (1982). 5. A. Klug and R. A. Crowther, "Three-dimensional image reconstruction from the viewpoint of information theory," Nature 238, 435-440 (1972). 6. G. Minerbo, "MENT, a maximum entropy algorithm for reconstructing a source from projection data," Computer Graphics and Image Proc. 10, 48-68 (1979). 7. J. A.Shohat and J. D.Tamarkin, The Problem of Moments, Vol. 1 of Math Surveys (American Mathematical Society, Providence, R.I., 1943). 8. J. E. Shore and R. W. Johnson, "Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy," IEEE Trans. Inf. Theory IT-26, 26-37 (1980). 9. L. M. Bregman, "The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming," USSR Comp. Math. Math. Phys. 7, 200-217 (1967). 10. J. Stoer and R. Bulirsch, Introduction to Numerical Analysis (Springer-Verlag, ACKNOWLEDGMENT This work was done under U.S. Department Berlin, 1980). 11. A. H. Stroud and D. Secrest, Gaussian Quadrature Formulas of Energy con- tract no. W-7405-ENG-48 at the Lawrence Berkeley Laboratory. (Prentice-Hall, Englewood, Cliffs, N.J., 1966). 12. G. E. Forsythe, M. A. Malcolm, and C. B. Moler, Computer Methods for Mathematical Computations (Prentice-Hall, Englewood Cliffs, N.J., 1977).