ROBUST REGULARIZED TOMOGRAPHIC IMAGING WITH CONVEX PROJECTIONS Farzad Kamalabadi and Behzad Sharif Department of Electrical and Computer Engineering and Coordinated Science Laboratory University of Illinois at Urbana-Champaign, IL, USA Email: [email protected] Fax: (217) 333-4303 ABSTRACT A robust method for tomographic image reconstruction from limited-angle noisy measurements is proposed which builds upon a combination of regularization theory and the method of projections onto convex sets (POCS). Two specific formulations of the proposed method, namely, Tikhonov-POCS and TV-POCS, are introduced and investigated. A statistical framework is developed that provides insight into the behavior of the two algorithms. The inclusion of a reference image is approached by either a coarse reconstruction or a model generated background image. The method is validated in the context of simulations for the reconstruction of highly structured images from partial projections. Results demonstrate significant improvement over conventional regularization methods in situations where the conventional techniques are inadequate. 1. INTRODUCTION The need to reconstruct images where the data is limited in its angular range occurs in many tomographic imaging applications. Data acquisition may be limited by physical or practical constraints, e.g., space imaging, or by time constraints, e.g., medical imaging. In principle, the problem of image reconstruction from projections can be formulated by the Fredholm integral equation of the first kind [1] . In practice, the observations are often a discrete sequence of measured data. Discretizing the integral equation and appropriate lexicographic ordering leads to the following observation model in form of a matrix equation: y = Ax + w where y is the vector of measurements, x is the discretized unknown field and w is the additive measurement noise that is assumed to be independent of the signal. The discretized tomographic kernel is denoted by A ∈ Rm×n . In limiteddata imaging applications A is highly ill-conditioned and in many cases it is also a singular matrix. Consequently, the pseudo-inverse estimate is unstable. A formal framework for addressing this problem is provided by regularization theory which attempts at replacing a given ill-conditioned problem with a well-conditioned one whose solution is an acceptable approximation to the solution of the original problem. This is usually done by adding a penalty term to the least squares functional. This added term is typically a roughness penalty which indirectly incorporates a smoothness constraint on the underlying field. However, standard regularization methods, such as Tikhonov or total variation (TV) regularization, prove to be inadequate in case of highly structured images and noisy partial data. Increasing the regularization parameter can resolve this issue only at the expense of smoothing out the edges in the image, hence decreasing its information content. The key to overcoming this limitation of a reconstruction technique is to formulate and include additional types of available a priori knowledge in the form of constraints on the reconstructed image. A large class of constraints in image processing are convex. Examples include knowledge about spatial or frequency domain support of the image, amplitude range, energy constraint, the reference (prototype) image constraint, and similarity to a reference image (see Chapters 2, 8, and 11 of [2]). Unfortunately, there is no straightforward way of formulating such constraints in a standard regularized reconstruction method. A class of image reconstruction techniques that specifically deals with such constraints is the method of POCS (see [3] and references therein). In POCS one summarizes all the prior information and any available data about the unknown image in a sequence of closed, convex constraint sets C1 , · · · , Cn . T Then, a feasible solution must n lie in the intersection C = i=1 Ci . Such a solution satisfies, by definition, all the a priori and data-driven constraints, and in the absence of further knowledge or other criteria is the best that one can hope to achieve. The fundamental theorem of POCS states that a point in C can be found iteratively through xk+1 = PC1 PC2 ...PCn xk (1) where xk+1 is a new estimate, xk is the previous estimate, and PCi is a projection operator (projector) that projects onto Ci for i ∈ {1, 2, ..., n}. The projection PCi x is the element in Ci closest to x. Under mild technical conditions on the constraint sets, the POCS algorithm (equation (1)) is guaranteed to converge to a point in C [3]. The main drawback of POCS as a set-theoretic method is its inability to incorporate both signal and noise statistics. In contrast, methods such as TV regularization incorporate such information and have superior edge preserving characteristics [4]. The goal of this work is to overcome the limitations of each approach through incorporation of a priori information in the form of combined regularization and convex projections. Section 2 describes the proposed method in detail. In Section 3, the proposed technique is studied in the context of constrained maximum likelihood estimation. In Section 4, through numerical simulations, the proposed method is validated for the application of tomographic imaging of the ionosphere. Fig. 1. Left panel: flowchart for an iterative reconstruction algorithm. Right panel: flowchart for the proposed algorithm (Reg-POCS). 2. THE PROPOSED METHOD For computational efficiency, most reconstruction methods (e.g., Tikhonov regularization) are formulated to utilize an iterative algorithm such as conjugate gradient (CG) for optimizing the corresponding objective functional. The algorithm iteratively updates the current estimate of the image until it converges to a solution (see left panel of Figure 1). The form of the update equation varies with the choice of the optimization method used. In this work, for clarification purposes, it is assumed that the optimization method used is CG and the only available convex constraint is the reference constraint defined as: CR = {x ∈ Rm×n : kx − xR k ≤ ρ} (2) which defines a sphere whose center is at xR , the reference image, and has a radius of ρ > 0. The corresponding projection is given by ½ x if kx − xR k ≤ ρ (3) PCR x = x−xR if kx − xR k > ρ xR + ρ kx−x Rk The proposed method, named Reg-POCS, projects the output of the CG loop onto the reference convex set using equation (3) (see right panel of Figure 1). The projection result is then fed back to the next iteration of the Reg-POCS algorithm (the outer loop). The stopping criteria for CG loop and the Reg-POCS iterations are important factors in the behavior of the algorithm. The convergence condition for the CG loop is when the ℓ2 norm of the difference of the previous estimate and the updated estimate is less than a predesigned threshold. The convergence condition for the Reg-POCS algorithm is when all of the projections are inactive (projection P is inactive on x if P x = x). In Tikhonov regularization, the update equation is the CG update for the following linear inverse problem [1]: (AT A + λLT L)ˆ xTikhonov = AT y (4) where λ is the regularization parameter and L is the positive semidefinite regularization matrix usually taken to be the first order approximation to the gradient operator. Taking the Tikhonov regularization method as the original algorithm, the proposed method, named Tikhonov-POCS, has one extra step which is the convex projections and one extra loop (see right panel of Figure 1). The following pseudocode describes the Tikhonov-POCS algorithm using CG and reference constraint: 1. Q := AT A + λLT L 2. % CG starts here 3. x′0 := xinit , r0 := y − Qx′0 , p0 := r0 , n := 0 4. While kx′n+1 − x′n k ≤ tolerance Do: 5. αn := pTn rn /pTn Qpn 6. x′n+1 := x′n + αn pn 7. rn+1 := y − Qx′n+1 8. βn := pTn Qrn+1 /pTn Qpn 9. pn+1 := rn+1 − βn pn 10. n := n + 1 11. % CG ends here 12. xcg := x′n+1 13. if kxcg − xR k > ρ % reference projection x −x 14. cg R x := xR + ρ kxcg −xR k 15. xinit := x, Go to 2 16. else x := xcg % happens when x satisfies the reference constraint 17. Return x Also, the TV solution is obtained by finding the fixed point solution to the following linear inverse problem [4]: (n) (n+1) (AT A + λLT W(ˆ xT V )L)ˆ xT V = AT y (5) where W(x) is a weight matrix defined as # " 1 1 W(x) = diag p 2 |[Lx]i |2 + φ and φ is a real constant usually chosen to be very small [4]. The update equation for TV regularization is the CG update equation for solving (5). Taking the TV regularization method as the original algorithm, the proposed method, named TV-POCS, has one extra step which is the convex projections and one extra loop (Figure 1). The pseudo-code for TV-POCS method is similar to the Tikhonov-POCS case but there is another loop over the CG loop that implements the fixed point iterations of equation (5). Also the Q matrix (n) is defined as Q = AT A + λLT W(ˆ xT V )L. As will be presented in the results (Section 4), the reference image can be chosen adaptively through a coarse initial reconstruction or through a physics based model. 3. ANALYSIS OF THE ALGORITHM This section develops a statistical framework which provides insight into the Tikhonov-POCS and TV-POCS algorithms. Typically the intersection of the available convex constraint sets, denoted by C in Section 1, contains more than one image. The image that belongs to C and has the maximum likelihood of generating the observed data is a reasonable estimate under modeling uncertainty. This provides a signal estimation framework which is referred to as the constrained maximum likelihood (CML) [5]. However, for the case of Tikhonov-POCS and TV-POCS, the distribution of the unknown image, pX (x) is also available. Hence, the problem is in fact a constrained maximum a posteriori (CMAP) estimation problem: ˆ CMAP (y, C) = arg max pX|Y (x|y) x x∈C £ ¤ = arg min − log pY|X (y|x) − log pX (x) x∈ C 3.1. Statistical Analysis of CMAP Tikhonov-POCS Tikhonov regularization is equivalent to assuming Gaussian statistics for both the signal and noise [1]. Here we assume that the only convex constraint available is a ℓ2 -norm reference constraint as stated in equation (2). According to the CMAP formulation and a fundamental theorem in convex optimization (Theorem 1, p. 217 of [6]), we have: i h ˆ CMAP (y) = arg min ky − Axk2Σ−1 + kx − xR k2Σ−1 x w x x∈CR h ky − Axk2Σ−1 + kx − xR k2Σ−1 = arg min w x x∈Rm×n ¤ + ν(ρ, y)kx − xR k22 h i = arg min ky − Axk2Σ−1 + kx − xR k2(Σ−1 +ν(ρ,y)I) x∈Rm×n w x where ν(ρ, y) ≥ 0 is the data-dependent Lagrange multiplier associated with the CMAP problem and Σx and Σw are the covariance matrices of the signal and noise respectively (in equation (4), it was assumed that noise is i.i.d., xR = 0 and Σx is proportional to LT L). In general, the Lagrange multiplier has a complicated dependence on ρ and y and it can be shown that this dependence is not one-to-one [5]. It is evident from the foregoing analysis that the CMAP estimation is equivalent to a new unconstrained Tikhonovregularized reconstruction. The new signal covariance ma˜ x = (Σ−1 +ν(ρ, y)I)−1 and the new regularization trix is Σ x 1 2 ˜ =Σ ˜− ˜ −1 , i.e., L matrix is the matrix square root of Σ x . It x is worth mentioning that computing the Lagrange multiplier ν(ρ, y) involves numerical evaluation of high-dimensional integrals. A larger value of ρ models greater uncertainty in the reference constraint and as ρ → +∞, ν(ρ, y) → 0. The Tikhonov-POCS algorithm is an approximate implementation of the equivalent MAP estimation derived above. Also note that this analysis was done assuming a single reference constraint. In case of multiple convex constraints the statistical analysis may become very difficult. 3.2. Statistical Analysis of CMAP TV-POCS TV regularization is equivalent to assuming Gaussian distribution for the noise and Laplacian distribution for the signal [1]. Here it is assumed that the only convex constraint available is a ℓ1 -norm reference constraint: C1 = {x ∈ Rm×n : |x − xR | ≤ ρ1 }, ρ1 ≥ 0 where |.| denotes the ℓ1 norm. According to the CMAP formulation and the theorem mentioned in Section 3.1, we have: h i ˆ CMAP (y) = arg min ky − Axk2Σ−1 + λ|Lx| = x w x∈C1 h i arg min ky − Axk2Σ−1 + λ|Lx| + ν1 (ρ1 , y)|x − xR | x∈Rm×n w where ν1 (ρ1 , y) is the associated Lagrange multiplier. Here it can be seen that, unlike the case in Section 3.1, the CMAP estimation does not reduce to unconstrained TV regularization. Nevertheless, the terms inside the sum can be interpreted as a new penalty term. The regularization parameter λ and the Lagrange multiplier ν1 (ρ1 , y) control the balance between how much sharp gradients are penalized versus how far the reconstruction is from the reference image xR . The TV-POCS algorithm is an approximate implementation of the equivalent MAP estimation. 4. NUMERICAL RESULTS Figure 2(a) corresponds to a two-dimensional section of the electron density in the ionosphere containing a plasma depletion. The reconstruction of such images from projection 0.5 −0.5 (a) 5 0 0.5 1 2 300 2 300 1 200 1 200 0 100 −1 Altitude(Km) 4 400 3 300 200 0.5 1 3 5 (cmx) 10 6 5 −0.5 600 500 400 2 300 1 200 100 −1 700 3 2 1 −0.5 0 0.5 5 4 400 3 2 3 300 2 300 200 1 200 (d) Fig. 2. (a) Original ionosphere (b) Best possible ART reconstruction result (c) Best possible Tikhonov reconstruction (d) Best possible TV reconstruction. data is subject of many experiments, and poses grand challenges due to the fact that the imaging scenario is inherently data-limited and the underlying image has a non-smooth structure. The reconstructed images were obtained from such limited-angle noisy data using algebraic reconstruction technique (ART), Tikhonov regularization, and TV regularization and are shown in Figure 2. For each technique, the corresponding image represents the best possible reconstruction with respect to its tunable parameters. As can be seen, these reconstructions bear little resemblance to the original image. The Tikhonov regularized reconstruction is better than ART, but the best result is achieved by the TV reconstruction shown in Figure 2(d). This is expected as TV regularization gives superior results in presence of sharp edges which is the case at hand [4]. But even the TV reconstruction is not satisfactory since it is physically impossible to have high electron density at low (≤ 150km) or high (≥ 650km) altitudes. In order to overcome this limitation, we apply the proposed TV-POCS technique method. The reference image in Figure 3(a) was obtained by an initial coarse reconstruction. In other words, to obtain a reference image the number of unknowns is reduced, hence reducing the ill-conditioning. The TV-POCS reconstruction using this reference image is shown in Figure 3(b). It can be seen that the background and the bubble structure are better reconstructed. To have a quantitative measure of the improvement, the normalized mean square error defined as −xtrue k was computed. The MSE for TV reconMSE = kˆxkx true k struction is 40.85% and for the TV-POCS reconstruction is 29.27%. Next we study the proposed method using a reference image that is not adaptive to the data but instead it is taken to be the quiescent background generated from a physics-based model. The reference image is shown in Figure 3(c). The output of TV-POCS using this reference is shown in Figure 3(d). The MSE as defined above for TV 0 0.5 1 6 500 400 Longitude 3 (cm x) 600 4 −0.5 el/ TV−POCS Output 700 5 100 −1 1 (b) 5 Reference Image (PIM generated) el/(cm3) x 10 6 500 1 0.5 Longitude 600 1 0 100 −1 1 4 2 Longitude 0.5 5 3 −0.5 0 Longitude 4 400 (a) el/(cmx3)105 6 Total Variation Output 700 Altitude(Km) x 10 6 5 500 (c) 3 (b) el/(cm3) 600 0 400 3 Longitude Tikhonov Output Longitude 500 400 100 −1 1 700 −0.5 4 4 Altitude(Km) 0 500 500 200 1 Longitude 600 300 2 −0.5 5 Altitude(Km) Altitude(Km) Altitude(Km) 3 200 el/ TV−POCS Output 700 600 5 4 300 5 Reference Image(Coarse Reconstruction) el/(cm3)x 10 700 6 5 400 100 −1 el/(cm3) x 105 ART Reconstruction 700 600 500 100 −1 5 10 6 Altitude(Km) 3 (cm )x 600 Altitude(Km) el/ Original Ionosphere 700 100 −1 5 10 1 −0.5 0 Longitude (c) 0.5 1 (d) Fig. 3. (a) Reference image taken to be a data-adaptive coarse reconstruction. (b) TV-POCS reconstruction using the image in panel (a) as the reference. (c) Reference image taken to be a model-generated quiescent background. (d) TV-POCS reconstruction using the image in panel (c) as the reference. reconstruction is 40.85% and for the TV-POCS reconstruction is 20.3%. The reconstructed background is much improved and the MSE is significantly lowered. Also, the depleted structure is more pronounced. The result in this case is better than the case where the reference was obtained by a coarse reconstruction. The reason is that by incorporating the knowledge of the background a large amount of a priori knowledge is included which was not incorporated for the former case. 5. REFERENCES [1] W. C. Karl, “Regularization in image restoration and reconstruction,” in Handbook of Image and Video Processing, Al Bovik, Ed., chapter 3.6, pp. 141–160. San Diego: Academic Press, 2000. [2] H. Stark, Ed., Image Recovery, Theory and Application, New York: Academic Press, 1987. [3] P. L. Combettes, “The foundations of set theoretic estimation,” Proc. IEEE, vol. 81, no. 2, pp. 182–208, Feb 1993. [4] C. R. Vogel and M. E. Oman, “Fast, robust total variation-based reconstruction of noisy, blurred images,” IEEE Trans. Image Processing, vol. 7, no. 6, pp. 813–824, 1998. [5] P. Ishwar and P. Moulin, “On the equivalence of settheoretic and maxent map estimation,” IEEE Trans. Image Processing, vol. 51, no. 3, pp. 698–713, 2003. [6] D. G. Luenberger, Optimization by Vector Space Methods, New York: Wiley, 1969.
© Copyright 2024