Reconstruction of CT Images by
Iterative Least Squares Methods
with Nonnegative Constraint
Hiromasa Kono, Yuichi Tanji, Ken’ichi Fujimoto,
Hiroyuki Kitajima, Yo Horikawa (Kagawa Univ.),
Norikazu Takahashi(Okayama Univ.)
Background
• Filtered Back-Projection(FBP) and iterative methods have been
used to reconstruct CT images
• Advantage of these methods
Image reconstruction time is short for the FBP method
Reconstructed image is high quality for the Iterative Method
• Disadvantage of these methods
Artifacts are appeared in images obtained by the FBP method
The huge amount of memory and time are required for the iterative
method
• The linear system represented by
Basic problem in CT system


Projection operator

(known)
Projection dataset
(known)
Pixel value of the true image
(unknown)
P is projections, N is the pixel size
Due to the noise and systematic measurement errors,
solving Eq.(1) is often treated as an ill-posed problem


Nonnegative matrix factorization
• NMF(Nonnegative Matrix Factorization）is factor analysis method
• For given nonnegative matrix

,
approximate by the product of two nonnegative matrices

,


NMF for CT image reconstruction
• For given nonnegative matrix

and nonegative
vector
the solution vector
is a pixel value
• Evaluating both sides with Frobenius norm

 


• The projection operator

is represented by 
and the solution vector x is represented by
when updating
the corresponding residual
must be evaluated
By replacing the calculation of
with the above formula, we can speed up
the calculation
HALS






 


CG(Conjugate gradient) method
• CG is effective for solving the linear equations with sparse
matrix

 
If the matrix is symmetric positive definite
matrix
 represents the inner product
Starting from the initial matrix
, search for that
minimizes the cost function
•
is the searching direction vector for solution vector
If the n searching direction vector

are linearly independent,
the solution of linear system is represented by linear combination


By satisfying an appropriate orthogonal relationship for

the algorithm becomes effective
Considering that 


 

CG(Conjugate gradient) method
CG(Conjugate gradient) method
• Apply CG to CT image reconstruction
For given projection operator R

(P<N)


 

By representing rectangular matrix as

Algorithm of CG
forend
for
1-nnn
1-n1-nn
prp
rr
rr
prr
pxx
pAp
rr
rpAxyArx
n
n
T
n
n
T
n
n
n
nnnn
n
TT
n
n
T
n
n
T
n
),(
),(
),(
),(
...,3,2,1
),(,0
11
11
11
11
0000
Starting from the initial vector
,
search for that minimizes the
cost function
is selected so that

and
are orthogonal
is selected to minimize
the cost function f(x)
PCG（preconditioned conjugate gradient）
• PCG(preconditioned conjugate gradient)
We use the diagonal scaling as preconditioner
Extract only diagonal elements of
and create the inverse matrix
can be regard as an approximation of inverse matrix



 




nearly approaches an identity matrix
forend
if end
thenif
forend
for
for
000
)(
)(
11
11
11
11
0000
),(
0
0
~
)
~
,(
)
~
,(
~
),(
),(
...,,3,2,1
...,3,2,1
),(,
frequencyrestart theChoose
MrpAxyAr
x
x
prp
rMr
rMr
Mrr
prr
pxx
pAp
Mrr
MrpAxyAryMAx
1-nnn
nn
1-n1-nn
T
T
N
N
n
n
T
n
n
T
n
n
n
nnnn
n
TT
n
n
T
n
n
T
mn
k
m
Algorithm of PCG with
nonnegative constraint
Since the solution in CT image
reconstruction is a pixel value, it
should be taken nonnegative values
The Solution by PCG is not
guaranteed to be positive
When an element of

is negative
at restart, it is forcibly replaced with
zero
NC : Nonnegative Constraint
Parameter Setting
• Setting parameter
Phantom image in MATLAB is used as the original image.
The size of image is 256×256 (N=65,536)
We set the scan range 180 and 100 projections,
assuming irradiation every 1.8 degrees (P=36,700)
For PCG with restart, set the number of restarts every 5, 10,
20 iterations, respectively
Number of iterations : 200
Evaluation criteria
• PSNR(Peak Signal-to-Noise Ratio) [dB]
Evaluation for quality of original image and
reconstructed image
• Residual Norm
Distance between the original projection dataset
and back-projection 
• Reconstruction time [sec]
HALS
• PSNR and residual norm obtained by HALS
0 100 200
10
0
10
2
10
4
Residual norm
Iterations
PSNR(dB) Residual norm Processing time
HALS 29.587 1.1792 836.36
0 100 200
10
20
30
PSNR (dB)
Iterations
CG and PCG without NC
• PSNR and residual norm obtained by PCG and CG without NC
0 100 200
10
0
10
2
10
4
Iterations
PCG without NC
CG without NC
Residual norm
0 100 200
0
20
CG without NC
PCG without NC
PSNR (dB)
Iterations
PSNR(dB) Residual norm Processing time
CG 14.030 475.1 87.844
PCG 27.783 0.225 88.755
PCG with NC
0 100 200
10
20
30
Restart every 5 iterations
Restart every 10 iterations
Restart every 20 iterations
Iterations
PSNR (dB)
• PSNR and residual norm obtained by PCG with NC
0 100 200
10
1
10
2
10
3
10
4
Iterations
Residual norm
Restart every 5 iterations
Restart every 10 iterations
Restart every 20 iterations
PSNR(dB) Residual norm Processing time
PCG with NC (every 5 iterations) 31.998 22.757 93.198
PCG（w/o NC）
0 100 200
10
0
10
2
10
4
Iterations
Residual norm
Restart every 5 iterations
Restart every 10 iterations
Restart every 20 iterations
Without NC and restart
• The difference between PSNR and residual norm （w/o NC）
Without NC : good residual norm
NC
With NC : good PSNR value
HALS,PCG with NC
0 100 200
10
20
30
Restart every 5 iterations
Restart every 10 iterations
Restart every 20 iterations
Iterations
PSNR (dB)
HALS
• PSNR and residual norm obtained by HALS and PCG with NC
0 100 200
10
0
10
2
10
4
Iterations
Residual norm
Restart every 5 iterations
Restart every 10 iterations
Restart every 20 iterations
HALS
PCG with NC (every 5 iterations) HALS
Processing time 93.198 836.36
Conclusion
• We presented CT image reconstruction simulation by HALS,
CG and PCG
By using CG and PCG, the CT image reconstruction time can
be shorter than that of HALS
There was a difference in evaluation criteria depending on NC
We will conduct the image reconstruction using GMRES
and Lanczos method
We will perform image reconstruction with the projection
operator and the projection dataset added noise
Future Works