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.