IEICE TRANS. FUNDAMENTALS, VOL.E104–A, NO.4 APRIL 2021

702

PAPER

Multiclass Dictionary-Based Statistical Iterative Reconstruction for
Low-Dose CT*
Hiryu KAMOSHITA† , Nonmember, Daichi KITAHARA†a) , Ken’ichi FUJIMOTO††b) , Members,
Laurent CONDAT†††c) , Nonmember, and Akira HIRABAYASHI†d) , Member

SUMMARY This paper proposes a high-quality computed tomography(CT) image reconstruction method from low-dose X-ray projection data. A
state-of-the-art method, proposed by Xu et al., exploits dictionary learningfor image patches. This method generates an overcomplete dictionary frompatches of standard-dose CT images and reconstructs low-dose CT imagesby minimizing the sum of a data fidelity and a regularization term based onsparse representations with the dictionary. However, this method does nottake characteristics of each patch, such as textures or edges, into account. In
this paper, we propose to classify all patches into several classes and utilize
an individual dictionary with an individual regularization parameter for eachclass. Furthermore, for fast computation, we introduce the orthogonality tocolumn vectors of each dictionary. Since similar patches are collected in thesame cluster, accuracy degradation by the orthogonality hardly occurs. Our
simulations show that the proposed method outperforms the state-of-the-art
in terms of both accuracy and speed.
key words: low-dose CT, image reconstruction, sparse representation, fastdictionary learning, clustering

1.

Introduction

X-ray computed tomography (CT) [2], [3] is widely used for
diagnosis and detection of various diseases because it scans
the inside of the human body noninvasively in a few seconds.
However, the induction of cancerous and genetic diseases
by X-ray exposure is a concern [4]. Therefore, it is desirable
to decrease the X-ray dose as much as possible. If we apply
the standard reconstruction method, filtered back projection(FBP), to low-dose projection data, then unnatural artifacts
and large noise appear in reconstructed images, which leads
to overlooking diseased tissues in the diagnosis [5].
To reconstruct a high-quality CT image from the lowdose projection data, we exploit the sparsity on the basis of
Manuscript received February 13, 2020.
Manuscript revised August 3, 2020.
Manuscript publicized October 6, 2020.
† The authors are with the Dept. of Information Science and Engineering, Ritsumeikan University, Kusatsu-shi, 525-8577 Japan.
†† The author is with the Dept. of Engineering and Design, Kagawa University, Takamatsu-shi, 761-0396 Japan.
††† The author is with the Visual Computing Center, King Abdullah University of Science and Technology, Thuwal, Saudi Arabia.
* This paper is an extended version of our previous paper presented at the IEEE ICIP, Taipei, Taiwan, September 2019 [1]. The
current paper describes the proposed method more precisely andnumerical simulations are conducted in a more realistic setting.
a) E-mail: d-kita@media.ritsumei.ac.jp (Corresponding author)b) E-mail: fujimoto.kenichi@kagawa-u.ac.jpc) E-mail: laurent.condat@kaust.edu.sad) E-mail: akirahrb@media.ritsumei.ac.jp
DOI: 10.1587/transfun.2020EAP1020

the compressed sensing theory [6]. In the image processing
field, patch-based dictionary learning is often used to acquiresparse representations of target images [4]–[9]. A dictionaryis generated as a matrix from training images, and the target
image patches are supposed to be expressed sparsely as linear combinations of a few column vectors of the dictionary.
In [4], Xu et al. utilized the patch-based dictionary learningwith statistical iterative reconstruction (SIR) [10]–[12] fromthe low-dose projection data. Although this method can wellremove streak artifacts, detailed structures of the target images are lost by over-smoothing. In order to reconstruct thedetailed structures while suppressing the artifacts and noise,
it is necessary to express each image patch more sparsely.
To this end, this paper proposes to extend the reconstruction method of Xu et al. [4] to a multiclass version. In
the proposed method, all training image patches are divided
into several classes, and a dictionary is created for each class.
Then, target image patches are also classified and supposed to
be expressed sparsely by the dictionary for each class. In thiscase, approximation errors of the sparse representations forthe patches would differ depending on the class. Therefore,
we also change the regularization parameter value for eachclass. Moreover, for fast reconstruction, we introduce the orthogonality to each dictionary, i.e., we replace overcomplete
dictionaries with orthogonal matrices. In general, reductionof column vectors of the dictionary causes accuracy degradation of the sparse representations, and thus low image quality.
On the other hand, in the proposed method, since the dictionary is utilized only for similar patches because of clustering,
the accuracy hardly deteriorates. Our simulations based on
actual CT images show that the proposed method can obtainreconstructed results better than those by the state-of-the-art
method [4] in a shorter calculation time.
2.

Basic Principles of X-ray CT

We formulate the X-ray CT image reconstruction in the 𝑥𝑦coordinate system as shown in Fig. 1. Let 𝜇(𝑥, 𝑦) ≥ 0 be the
X-ray attenuation coefficient distribution in a field of view.
X-ray projection is described with the 𝑢𝑣-coordinate system:𝑢 = 𝑥 cos 𝜃 + 𝑦 sin 𝜃,
𝑣 = −𝑥 sin 𝜃 + 𝑦 cos 𝜃,which is given by rotating the 𝑥𝑦-coordinate system throughan angle 𝜃 ∈ [0, 𝜋). The ideal projection data 𝑙 ★𝜃,𝑛 are defined

Copyright © 2021 The Institute of Electronics, Information and Communication Engineers

KAMOSHITA et al.: MULTICLASS DICTIONARY-BASED STATISTICAL ITERATIVE RECONSTRUCTION FOR LOW-DOSE CT

703

𝑣

𝑦

reconstructs 𝝁 by iteratively minimizing a cost function𝑢

X-ray source∥𝑅𝝁 − 𝒍 ∥ 2𝒘 + 𝜆Ψ( 𝝁) :=𝜇(𝑥, 𝑦)𝑥𝛾,,. 𝑣

...

𝑧,,-

𝑧,,/
Detectors𝑧,,.

X-ray projection in a parallel-beam CT system.

by line integrals, that∫are parallel to the 𝑣-axis, of the distribution 𝜇(𝑥, 𝑦) as 𝑙 ★𝜃,𝑛 := 𝜇(𝛾 𝜃 ,𝑛 (𝑣)) d𝑣 (𝑛 = 1, 2, . . . , 𝑁), where𝛾 𝜃 ,𝑛 (𝑣) is a straight path from X-ray source to the 𝑛th detector. These line integrals are called the Radon transform. The
Radon transform is discretized by numerical integration, as𝑙 ★𝜃,𝑛 = 𝒓 T𝜃 ,𝑛 𝝁 + 𝜉 ★𝜃,𝑛 ,(1)where 𝝁 ∈ R+𝐽√is a √nonnegative vector corresponding to a CT
image of size 𝐽 × 𝐽, 𝒓 𝜃 ,𝑛 ∈ R+𝐽 is a nonnegative vector corresponding to weights in the numerical integration, (·) T is thetranspose operator, and 𝜉 ★𝜃,𝑛 ∈ R is the discretization error.
The observed data 𝑧 𝜃 ,𝑛 ∈ Z+ at the 𝑛th X-ray detector isa nonnegative integer following a Poisson distribution of expected value 𝐸 [𝑧 𝜃 ,𝑛 ] = 𝑏 exp(−𝑙 ★𝜃,𝑛 ), where 𝑏 > 0 is the X-ray
radiation intensity of the X-ray source [4]. From this relation,
the ideal projection data is expressed as 𝑙 ★𝜃,𝑛 = ln(𝑏/𝐸 [𝑧 𝜃 ,𝑛 ]).
Actually, since 𝐸 [𝑧 𝜃 ,𝑛 ] is unknown, 𝑙 𝜃 ,𝑛 := ln(𝑏/𝑧 𝜃 ,𝑛 ) is utilized, instead of 𝑙 ★𝜃,𝑛 , as the projection data computed fromthe observed data 𝑧 𝜃 ,𝑛 . If we acquire the projection data with
𝑁 detectors for 𝐺 angles 𝜃 𝑔 ∈ [0, 𝜋) (𝜃 1 < 𝜃 2 < · · · < 𝜃 𝐺 ),the observation models in (1) are combined to a matrix form𝒍 = 𝑅𝝁 + 𝝃,(2)

where 𝒍 ∈ R𝐼 is a vector concatenating all 𝑙 𝜃𝑔 ,𝑛 , 𝑅 ∈ R+𝐼 × 𝐽 is amatrix whose row vectors are 𝒓 T𝜃𝑔 ,𝑛 , 𝝃 ∈ R𝐼 is an error vector,
and the number of the observations is 𝐼 := 𝐺 𝑁. In standarddose CT, since 𝐼 ≫ 𝐽 and the noise level ∥𝝃 ∥ 2 is small, 𝝁 isreconstructed with high accuracy by the FBP method [2], [3].
3.

𝑤𝑖 (𝒓 𝑖T 𝝁 − 𝑙 𝑖 ) 2 + 𝜆Ψ( 𝝁), (3)𝑖=1𝜃Fig. 1𝐼∑︁Dictionary-Based SIR from Low-Dose Projection

There exist two approaches for low-dose CT: (i) reducing thenumber 𝐺 of projection views and (ii) lowering the X-ray radiation intensity 𝑏. In the first approach, the number 𝐼 of row
vectors of 𝑅 decreases, implying rank(𝑅) < 𝐽. In the second
approach, the noise level ∥𝝃 ∥ 2 in (2) becomes larger. As a result, in either case, the standard FBP method generates streak
artifacts and large noise in the reconstructed images [5].
Sauer et al. proposed a method that reconstructs 𝝁 by amaximum a posteriori (MAP) estimation [10]. This method

where the statistical weight 𝑤𝑖 := 𝑧𝑖 = 𝑏 exp(−𝑙𝑖 ) for the 𝑖thprojection data 𝑙𝑖 is derived by the second-order Taylor series
approximation for the log-likelihood of the Poisson distribution, and Ψ : R 𝐽 → R is some regularization term based on aprior distribution of 𝝁. The regularization parameter 𝜆 > 0controls the balance between the data fidelity and regularization terms. The optimization algorithm of Sauer et al. and its
accelerated version which uses a majorization-minimization(MM) technique [12] are called SIR methods [11].
Xu et al. combined (3) and patch-based dictionary learning [4]. In this method, a target image 𝝁 is decomposed into𝑃
small√images. . . , 𝑆), called patches, of√ 𝐻𝑠 𝝁 ∈ R+ (𝑠 = 1, 2, 𝑃×
size 𝑃 × 𝑃, where 𝐻𝑠 ∈ {0, 1} 𝐽 denotes the 𝑠th patch
extraction matrix and 𝑆 is the number of the patches. Then, 𝝁
is reconstructed under an assumption that each patch can be
expressed by a linear combination of a few column vectors of
an appropriate dictionary 𝐷 ∈ R 𝑃×𝐾 , where 𝐾 is the numberof column vectors 𝒅 𝑘 ∈ R 𝑃 s.t. ∥ 𝒅 𝑘 ∥ 2 = 1 (𝑘 = 1, 2, . . . , 𝐾).
Let 𝒄 𝑠 ∈ R𝐾 be sparse representations (coefficients) for 𝐻𝑠 𝝁,and this method reconstructs the CT image 𝝁 by solving
minimize ∥𝑅𝝁 − 𝒍 ∥ 2𝒘 + 𝜆𝑆∑︁𝝁, (𝐷),𝐶(∥𝐻𝑠 𝝁 − 𝐷𝒄 𝑠 ∥ 22 + 𝜈𝑠 ∥𝒄 𝑠 ∥ 0 ), (4)𝑠=1where 𝐶 := (𝒄 1 , 𝒄 2 , . . . , 𝒄 𝑆 ) ∈ R𝐾 ×𝑆 is a coefficient matrixand 𝜈𝑠 > 0 is a Lagrangian multiplier for each patch. There
are two setups for the problem in (4): (i) 𝐷 is fixed to a matrix pre-learned from training images, and (ii) 𝐷 is updated
adaptively from a current estimate of the target image 𝝁. The
former is called global dictionary-based SIR (GDSIR), andthe latter is called adaptive dictionary-based SIR (ADSIR).
Since it is difficult to directly solve the problem in (4),the coefficient matrix 𝐶 (and the dictionary 𝐷 in ADSIR) andthe CT image 𝝁 are alternately updated from the initial image
𝝁 (0) . In this method, the initial image 𝝁 (0) is set to the FBP
result, in which negative pixel values are replaced with zeros.
The matrix 𝐶 (and 𝐷 in ADSIR) is updated by solving
minimize

𝑆
∑︁

(𝐷),𝐶

(∥𝐻𝑠 𝝁 (𝑚−1) − 𝐷𝒄 𝑠 ∥ 22 + 𝜈𝑠 ∥𝒄 𝑠 ∥ 0 ),

(5)

𝑠=1

where 𝝁 is fixed to the current estimate 𝝁 (𝑚−1) . In GDSIR, 𝐶
is updated by orthogonal matching pursuit (OMP) [13] (seeSect. 4.2 in detail). In ADSIR, after 𝐷 is updated by an online
dictionary learning method e.g. [14], 𝐶 is updated by OMP.
Thus ADSIR requires computational cost more than GDSIR.
The CT image 𝝁 is updated by solving
minimize ∥𝑅𝝁 − 𝒍 ∥ 2𝒘 + 𝜆
𝝁∈R+𝐽

𝑆
∑︁

∥𝐻𝑠 𝝁 − 𝐷𝒄 𝑠 ∥ 22

(6)

𝑠=1

with fixed 𝐶 and 𝐷. The optimal solution 𝝁∗ to (6) is given by

IEICE TRANS. FUNDAMENTALS, VOL.E104–A, NO.4 APRIL 2021

704

∗

T

𝝁 = 𝑅 𝑊𝑅 + 𝜆

𝑆
∑︁

! −1
𝐻𝑠T 𝐻𝑠T

𝑅 𝑊𝒍 + 𝜆

𝑠=1

𝑆
∑︁

!
𝐻𝑠T 𝐷𝒄 𝑠

4.2

Significance of Multiple Regularization Parameters 𝜆 𝑞

𝑠=1with a diagonal matrix 𝑊 = diag(𝒘) ∈ R+𝐼 ×𝐼 . However, if theimage size is large, this linear system is difficult to solve. To
avoid this, Xu et al. used the same MM technique as in [12].
Specifically, 𝝁 is updated from 𝝁 (𝑚−1) by (7), where 0 and 1
denote vectors with all components 0 and 1, respectively, ⊘ isthe component-wise division, and max returns a larger value
for each pair of corresponding components of two vectors.

In this part, we explain the merit of using multiple regularization parameters 𝜆 𝑞 . When updating coefficient matrix 𝐶 by
OMP in a similar manner to (5) of GDSIR, there are mainly
three types of update-rules [17]. Type 1 update-rule fixes the
Lagrangian multipliers as 𝜈𝑠 = 𝜈 (𝑠 = 1, 2, . . . , 𝑆) and solves
minimize𝑆∑︁𝐶(∥𝐻𝑠 𝝁 (𝑚−1) − 𝐷 𝑞 (𝑠) 𝒄 𝑠 ∥ 22 + 𝜈∥𝒄 𝑠 ∥ 0 ),

𝑠=1

4.

SIR with Multiclass Dictionary Learning

4.1

Multiclass GDSIR (MGDSIR)

where 𝑞(𝑠) denotes the class of the 𝑠th patch. Type 2 updaterule restricts the sparse level of each 𝒄 𝑠 under 𝑇 and solves

Recently, in high-speed MRI and low-dose CT fields, image
reconstruction methods using multiclass dictionary learning
have been proposed for higher compression ratio, noise reduction, and artifact reduction [5], [15], [16]. In these methods, all image patches are divided into multiple classes using
geometric directions or pixel values. Since an individual dictionary is generated for each class, sparse representations forthe patches become more accurate, and better reconstruction
results are obtained than when we utilize only one dictionary.
The methods in [15], [16] are proposed for high-speed MRI
and thus cannot be directly applied to low-dose CT images.
The method in [5], which is proposed for low-dose CT, uses
∥ 𝝁 − 𝝁 (0) ∥ 22 as the data fidelity term, and hence the reconstruction results largely depend on 𝝁 (0) obtained by the FBP
method. On the other hand, since GDSIR uses ∥𝑅𝝁 − 𝒍 ∥ 2𝒘
in (4) as the data fidelity term, the reconstruction results arerobust against the initial image 𝝁 (0) . Thus, in this paper, we
propose to extend GDSIR to a multiclass version in order toreconstruct higher-quality images from low-dose projection.
Let 𝑄 be the number of classes and 𝑞 = 1, 2, . . . , 𝑄 be
the index of the class. A dictionary 𝐷 𝑞 ∈ R 𝑃×𝐾𝑞 for the 𝑞thclass is generated in advance from standard-dose CT images,
where 𝐾𝑞 is the number of column vectors of 𝐷 𝑞 . Each patch
𝐻𝑠 𝝁 is classified into one of 𝑄 classes, and a patch index set
Ðfor the 𝑞th class is denoted by S𝑞 s.t. 𝑄
𝑞=1 S𝑞 = {1, 2, . . . , 𝑆}
and S𝑞 ∩S𝑞 ′ = ∅ (𝑞 ≠ 𝑞 ′ ). We suppose S𝑞 is fixed hereafter.
We reconstruct 𝝁 by solving a multiclass version of (4):minimize ∥𝑅𝝁−𝒍 ∥ 2𝒘 +

𝑄
∑︁

𝝁,𝐶

𝑞=1

𝜆𝑞

∑︁

(∥𝐻𝑠 𝝁−𝐷 𝑞 𝒄 𝑠 ∥ 22 +𝜈𝑠 ∥𝒄 𝑠 ∥ 0 ),

𝑠∈ S𝑞

(8)
where 𝜆 𝑞 > 0 is a regularization parameter for the 𝑞th class.
The problem in (8) can be solved by a similar update to that of
GDSIR. We call this method multiclass GDSIR (MGDSIR).

(
𝝁

(𝑚)

= max 0, 𝝁

"
(𝑚−1)

T

− 𝑅 𝑊 (𝑅𝝁(𝑚−1)− 𝒍) + 𝜆𝑆∑︁
𝑠=1minimize𝑆∑︁𝐶∥𝐻𝑠 𝝁 (𝑚−1) − 𝐷 𝑞 (𝑠) 𝒄 𝑠 ∥ 22

s.t. ∀𝑠 ∥𝒄 𝑠 ∥ 0 ≤ 𝑇.

𝑠=1

Type 3 update-rule restricts each error under 𝜀 and solves
minimize𝑆∑︁∥𝒄 𝑠 ∥ 0𝐶

s.t. ∀𝑠 ∥𝐻𝑠 𝝁 (𝑚−1) − 𝐷 𝑞 (𝑠) 𝒄 𝑠 ∥ 22 ≤ 𝜀.

𝑠=1Note that Type 2 and 3 problems can be expressed in Type 1
form as in (8) by using different Lagrangian multipliers 𝜈𝑠 .
In multiclass dictionary learning, the number of patches
is different for each class. Further, there is a possibility thatthe ratio of the number of patches in each class is quite different between training images and a target image. As a result,the approximation errors of the sparse representations, withthe same sparse level, differ depending on the class. For highquality image reconstruction, 𝐾𝑞 , 𝜈, 𝑇, and 𝜀 should be set todifferent values for each class, and then we should choose thebest update-rule. However, such a parameter tuning becomes
more burdensome as the number 𝑄 of the classes increases.
To avoid the above burden in the multiclass dictionary
learning, in the proposed method, we also extended the regularization parameter 𝜆 in (4) to the multiclass ones 𝜆 𝑞 in (8).
This extension allows us to reconstruct high-quality imagesby an appropriate tuning of each 𝜆 𝑞 even if we set 𝐾𝑞 , 𝜈, 𝑇,and 𝜀 to the same values for all classes. Specifically, we set
𝜆 𝑞 to a large value for a class of high representation accuracy,
and to a small value for a class of low accuracy. Indeed, 𝜆 𝑞removes the burdensome tuning in the dictionary learning.
5.

Acceleration and Update Formulae of MGDSIR

5.1

Fast MGDSIR (FMGDSIR)

To accelerate MGDSIR, we further propose to introduce theorthogonality to column vectors of each dictionary in (8) as#
𝐻𝑠T (𝐻𝑠 𝝁 (𝑚−1) − 𝐷𝒄 𝑠 )"⊘T

𝑅 𝑊𝑅 + 𝜆

𝑆
∑︁
𝑠=1

! #)
𝐻𝑠T 𝐻𝑠1(7)KAMOSHITA et al.: MULTICLASS DICTIONARY-BASED STATISTICAL ITERATIVE RECONSTRUCTION FOR LOW-DOSE CT

705

minimize ∥𝑅𝝁−𝒍 ∥ 2𝒘 +

𝑄
∑︁

𝝁,𝐶

𝑞=1

𝜆𝑞

∑︁

b𝑞 𝒄 𝑠 ∥ 2 +𝜈𝑠 ∥𝒄 𝑠 ∥ 0 ),(∥𝐻𝑠 𝝁− 𝐷
2

𝑠∈ S𝑞

(9)
b𝑞 ∈ R 𝑃× 𝑃 is an orthogonal dictionary s.t. 𝐷
bT𝑞 𝐷
b𝑞 =
where 𝐷𝐸, and 𝐸 denotes the identity matrix of order 𝑃. In general,
if we replace an overcomplete dictionary with an orthogonal
dictionary, then dictionary learning can be accelerated while
the accuracy of the sparse representation deteriorates. On theother hand, in the proposed method, since each dictionary has
only to deal with similar patches thanks to clustering, the accuracy hardly deteriorates even if we use the orthogonal dictionaries. We call this method Fast MGDSIR (FMGDSIR).
MGDSIR and FMGDSIR reconstruct 𝝁 from the initial
estimate 𝝁 (0) by alternately updating the sparse coefficientsand the image 𝑀 times after generating the global dictionaries. Because of the orthogonality, both dictionary generationand image reconstruction can be quickly done in FMGDSIR.
5.2

First, patches are extracted from training images and divided
into 𝑄 classes by the K-means algorithm in the same manneras in [16]. Note that the clustering centers are also used in thereconstruction step to classify target patches. In FMGDSIR,b𝑞 is generated† for each class asan orthogonal dictionary 𝐷an approximate solution to the optimization problem∑︁b𝑞 𝒄tr𝑠 ∥ 2 +𝜈𝑠 ∥𝒄tr𝑠 ∥ 0 ) s.t. 𝐷
bT𝑞 𝐷
b𝑞 = 𝐸,
minimize(∥𝐻𝑠 𝝁tr − 𝐷2𝑠∈ S𝑞tr(10)
where 𝝁tr is a standard-dose training image, 𝒄tr𝑠 is a sparse coefficient vector for the 𝑠th training patch 𝐻𝑠 𝝁tr , S𝑞tr is a patch
index set for training patches of class 𝑞, and all 𝒄tr𝑠 for 𝑠 ∈ S𝑞tr
are combined to a matrix 𝐶𝑞tr . We solve the problem in (10)b𝑞 is updated byaccording to [8], and thus while fixing 𝐶𝑞tr , 𝐷b∗𝑞 = argmin ∥Φtr𝑞 − 𝐷
b𝑞 𝐶𝑞tr ∥ 2F𝐷

bT𝑞 𝑈 is an orthogonal matrix and every component
Since 𝑉 T 𝐷
is smaller than 1, the cost function in (12) is maximized whenbT𝑞 𝑈 = 𝐸. As a result, the optimal solution is given by
𝑉T𝐷b∗𝑞 = 𝑈𝑉 T .
𝐷

b𝑞 = 𝐸,bT𝑞 𝐷
s.t. 𝐷

(11)

b𝑞
𝐷
tr

where Φtr𝑞 ∈ R 𝑃× | S𝑞 | is a matrix whose column vectors are
𝐻𝑠 𝝁tr s.t. 𝑠 ∈ S𝑞tr , and ∥·∥ F denotes the Frobenius norm.
The cost function in (11) is expressed asbT𝑞 Φtr𝑞
b𝑞 𝐶 tr ∥ 2F = ∥Φtr𝑞 ∥ 2F + ∥ 𝐷
b𝑞 𝐶 tr ∥ 2F − 2 Trace 𝐶 trT 𝐷∥Φtr𝑞 − 𝐷bT𝑞 ,= ∥Φtr𝑞 ∥ 2F + ∥𝐶 tr ∥ 2F − 2 Trace Φtr𝑞 𝐶 trT 𝐷
where Trace returns the sum of all diagonal components of asquare matrix. In this equation, the first and second terms on
b𝑞 . Thus, (11) is expressed as
the right side are constant for 𝐷b∗𝑞 = argmax Trace Φtr𝑞 𝐶 trT 𝐷
bT𝑞 s.t. 𝐷
bT𝑞 𝐷
b𝑞 = 𝐸. (12)𝐷b𝑞
𝐷
† In MGDSIR, an over-complete dictionary 𝐷

𝑞 is generated bysolving the problem in (10) without the orthogonality constraint.

(13)

After the dictionary update by (13), we update a coeffib𝑞 = 𝐷b∗𝑞 in (10). In FMGDSIR,
cient matrix 𝐶 tr while fixing 𝐷
we utilize the property that the ℓ2 -norm is unitarily invariant,and the coefficient matrix is updated†† by∑︁
b𝑞 𝒄tr𝑠 ∥ 2 + 𝜈𝑠 ∥𝒄tr𝑠 ∥ 0 )𝐶𝑞tr∗ = argmin(∥𝐻𝑠 𝝁tr − 𝐷
2
𝐶𝑞tr

Global Dictionary Generation

b𝑞 ,𝐶𝑞tr
𝐷

This problem is known as the orthogonal Procrustes problem[6] and can be solved as follows. By using singular value decomposition (SVD) Φtr𝑞 𝐶 trT = 𝑈Σ𝑉 T with one diagonal matrix Σ ∈ R+𝑃× 𝑃 and two orthogonal matrices 𝑈 ∈ R 𝑃× 𝑃 and𝑉 ∈ R 𝑃× 𝑃 , the cost function in (12) yieldsbT𝑞 = Trace 𝑈Σ𝑉 T 𝐷
bT𝑞 = Trace Σ𝑉 T 𝐷
bT𝑞 𝑈 .
Trace Φtr𝑞 𝐶 trT 𝐷= argmin𝐶𝑞tr𝑠∈ S𝑞tr∑︁bT𝑞 𝐻𝑠 𝝁tr − 𝒄tr𝑠 ∥ 2 + 𝜈𝑠 ∥𝒄tr𝑠 ∥ 0 ).
(∥ 𝐷
2

(14)

𝑠∈ S𝑞tr

In the orthogonal dictionary learning, Type 1 update-rule isoften used [8] and (14) is solved for each column vector 𝒄 𝑠 as
tr
tr 2
tr
bT
𝒄tr∗
𝑠 = argmin ∥ 𝐷 𝑞 𝐻 𝑠 𝝁 − 𝒄 𝑠 ∥ 2 + 𝜈∥𝒄 𝑠 ∥ 0 ,

(15)

𝒄 𝑠tr

where the Lagrangian multipliers are fixed as 𝜈𝑠 = 𝜈 for all𝑠 ∈ S𝑞tr . The optimal solution to (15) is easily given by
tr
√
bT
𝒄tr∗
𝑠 = hard_thresh 𝜈 ( 𝐷 𝑞 𝐻 𝑠 𝝁 )(16)with the hard thresholding operator:(√
𝒄[𝑖] if |𝒄[𝑖] | ≥ 𝜈,
hard_thresh√𝜈 (𝒄) [𝑖] :=√0
if |𝒄[𝑖] | < 𝜈,
where [𝑖] stands for the 𝑖th component of a vector. Note thatType 2 and 3 update-rules are also easily given. By repeatingb𝑞 are obtained.
(13) and (16) for each class, the dictionaries 𝐷
5.3

Image Reconstructionb𝑞 ,To solve (9) with the multiclass orthogonal dictionaries 𝐷(0)first, the initial estimate 𝝁 is obtained by applying the FBP
method to the low-dose projection data 𝒍. Then, by using thestored clustering centers in the dictionary generation step, theinitial patches 𝐻𝑠 𝝁 (0) are classified into one of 𝑄 classes, andhereafter the patch index sets S𝑞 are fixed. We reconstruct animage by alternately updating 𝐶 and 𝝁 for 𝑚 = 1, 2, . . . , 𝑀.
At the 𝑚th update of 𝐶, if we choose Type 1 update-rule,then we solve, as shown in Sect. 4.2, the following problem†† In MGDSIR, coefficients are updated by OMP as in Sect. 4.2.

Note that 𝐷 𝑞 can be updated by the same rule for any OMP type.

IEICE TRANS. FUNDAMENTALS, VOL.E104–A, NO.4 APRIL 2021

706

Class 1 (38,529)(a) Training image (f_1601).

𝑆
∑︁

𝐶

b𝑞 (𝑠) 𝒄 𝑠 ∥ 2 + 𝜈∥𝒄 𝑠 ∥ 0 ), (17)(∥𝐻𝑠 𝝁 (𝑚−1) − 𝐷
2where 𝝁 is fixed to the current estimate 𝝁 (𝑚−1) . In the samemanner as (14), the problem in (17) is solved for each 𝒄 𝑠 asbT 𝐻𝑠 𝝁 (𝑚−1) ).
𝒄∗𝑠 = hard_thresh√𝜈 ( 𝐷
𝑞 (𝑠)

(18)

Differently from the problem in (5), the orthogonal dictionaries allow us to quickly compute the exact closed-form solution to (18) without using a greedy algorithm such as OMP.
At the 𝑚th update of 𝝁, we solve a least squares problem
minimize ∥𝑅𝝁 − 𝒍 ∥ 2𝒘 +
𝝁∈R+𝐽

𝑄
∑︁

𝜆𝑞

𝑞=1

∑︁

b𝑞 𝒄 𝑠 ∥ 2 (19)
∥𝐻𝑠 𝝁 − 𝐷2𝑠∈ S𝑞

for a fixed 𝐶. For the problem in (19), to avoid the computation of a large inverse matrix, we use the same MM technique
as in [12], and 𝝁 is updated from 𝝁 (𝑚−1) by (20). Note thatb𝑞 is replaced with 𝐷 𝑞 in the update by (20).
in MGDSIR, 𝐷
Numerical Simulations and Discussions

In this section, we compare the proposed methods, MGDSIRand FMGDSIR, with the conventional methods, ADSIR and
GDSIR, by simulations based on actual CT images in [18].
6.1

Simulation Data and Environment

We used a female head image ‘f_1601’ as training data, andfemale head images ‘f_1610’ and ‘f_1615’ and male head
images ‘m_1114’ and ‘m_1132’ were targets. The original
sizes of these images were 512 × 512. By cropping the outer
black pixels, we made training and target images composedof the central 256 × 256 pixels (𝐽 = 65,536). Thus, the pixel
length in the training and target images was the same as the

(
𝝁

(𝑚)

Class 4 (2,246)

Class 5 (1,571)(b) Clustering results of 𝑄 = 5 for the training image and the five clustering centers.

𝑠=1

6.

Class 3 (3,400)

Training image used for dictionary learning and the clustering results of the K-means algorithm in the proposed method.

Fig. 2minimizeClass 2 (18,263)= max 0, 𝝁

"
(𝑚−1)

T

− 𝑅 𝑊 (𝑅𝝁(𝑚−1)− 𝒍) +

𝑄∑︁
𝑞=1

original pixel length, i.e., 1 mm. CT values 𝒗CT in the images
were transformed to the X-ray attenuation coefficients by 𝝁 =𝜇water (1+𝒗CT /1000), where 𝜇water = 0.2059 cm−1 . In (2), theRadon transform 𝑅 was created on the basis of parallel beams
to 579 X-ray detectors (𝑁 = 579) of width 0.625 mm. Initial
images 𝝁 (0) were obtained by FBP with the Ram–Lak filter.
We reconstructed the target images by ADSIR, GDSIR,MGDSIR, and FMGDSIR with 𝑀 = 1,000 iterations. We set
𝐾𝑞 = 𝐾 = 256 and 𝜈𝑠 = 𝜈 = 0.001 for ADSIR, GDSIR and
MGDSIR, and set 𝜈𝑠 = 𝜈 = 0.0007 for FMGDSIR. In overcomplete dictionary leaning and OMP, we used the programin [17] for fast implementation. Dictionaries were generated
with 100 iterations in ADSIR, with 2,000 iterations in GDSIRand MGDSIR, and with 1,000 iterations in FMGDSIR. Numerical simulations were conducted with MATLAB R2018a64-bit on a MacBook Pro, OS Mojave ver. 10.14.6, CPU 3.1
GHz Intel Core i5, and memory 8 GB 2133 MHz LPDDR3.
6.2

Clustering and Regularization Parameters

We extracted 𝑆 = 64,009 image patches of size 4×4 (𝑃 = 16)
by shifting the areas pixel by pixel. After clustering by the Kmeans algorithm, the mean value, called a DC component, ofeach patch is subtracted right before dictionary learning [17].
The dictionaries of GDSIR, MGDSIR, and FMGDSIR weregenerated from the training image in Fig. 2(a), where we divided the training patches into 𝑄 = 5 classes in the proposed
methods. The clustering results of the K-means algorithm
and each clustering center are shown in Fig. 2(b). The cluster
indices were ordered by the number of patches in descending order. The class 1 is mainly composed of dark patches,e.g., the background and the maxillary sinus. The class 2 ismainly composed of gray patches of soft tissues. The class 4
is mainly composed of middle-level patches in the boundary
between the classes 1 and 2. The classes 3 and 5 are composed of bright patches of bones and detailed structures.

#
𝜆𝑞

∑︁
𝑠∈ S𝑞

b𝑞 𝒄 𝑠 )𝐻𝑠T (𝐻𝑠 𝝁 (𝑚−1) − 𝐷

"
⊘

T

𝑅 𝑊𝑅+

𝑄
∑︁
𝑞=1

! #)
𝜆𝑞

∑︁
𝑠∈ S𝑞

𝐻𝑠T 𝐻𝑠1(20)

KAMOSHITA et al.: MULTICLASS DICTIONARY-BASED STATISTICAL ITERATIVE RECONSTRUCTION FOR LOW-DOSE CT

707(a) Ground truth (m_1132).

(b) FBP [41.03 / 0.950].

(c) ADSIR [47.07 / 0.995].

(d) GDSIR [47.00 / 0.995].
(e) MGDSIR [47.20 / 0.995].
(f) FMGDSIR [47.27 / 0.995].
Fig. 3 Reconstruction result of m_1132 from its standard-dose projection data by each method [PSNR / SSIM].

Table 1
SIR \ Image
ADSIR
GDSIR
MGDSIRFMGDSIRPSNR [dB] and calculation time [sec] in 60/300 views.
f_1610
37.43
37.27
38.13
38.14

f_1615
37.81
37.66
38.34
38.41

m_1114
38.16
38.04
38.78
38.92

m_1132
36.76
36.67
37.24
37.20

calc. time
0.240/iter.
0.146/iter.
0.168/iter.
0.090/iter.

Table 2PSNR [dB] and calculation time [sec] in 1/40 X-ray intensity.

SIR \ Image
ADSIR
GDSIR
MGDSIRFMGDSIR

f_1610
38.15
38.46
39.12
38.50

f_1615
39.07
39.36
39.97
39.36

m_1114
39.00
39.33
39.73
39.24

m_1132
38.00
38.20
38.68
38.18

calc. time
0.345/iter.
0.244/iter.
0.266/iter.
0.185/iter.

To determine the values of the multiple regularization
parameters 𝜆 𝑞 in MGDSIR and FMGDSIR, we use that of 𝜆in GDSIR as a reference. For patches in the classes 1 and 2,
DC components are dominant, and the patches are well approximated by both the overcomplete dictionaries 𝐷 𝑞 and theb𝑞 . Hence, 𝜆1 and 𝜆2 can be larger
orthogonal dictionaries 𝐷
than 𝜆. On the other hand, patches in the classes 3 and 5 arevery complex and not well approximated by the dictionaries,
and hence 𝜆3 and 𝜆5 can be smaller than 𝜆. The class 4 is difficult to handle. Patches in the class 4 mainly correspond tothe boundary between the classes 1 and 2, which are better
approximated by the overcomplete dictionary than by the orthogonal dictionary. Based on these considerations, we firstadjusted 𝜆 in GDSIR, and then adjusted 𝜆 𝑞 in MGDSIR and
FMGDSIR manually by the simple change-and-test process.

larization parameters to 𝜆 = 2,000 in ADSIR and GDSIR, to
𝜆1 = 2,560, 𝜆2 = 2,560, 𝜆3 = 1,280, 𝜆4 = 2,880, 𝜆5 = 1,280
in MGDSIR, and to 𝜆1 = 4,000, 𝜆2 = 2,600, 𝜆3 = 1,600,
𝜆4 = 2,200, 𝜆5 = 1,600 in FMGDSIR. Figures 3(a), (b), (c),(d), (e), and (f) show the original image, the initial image 𝝁 (0)given by FBP, the reconstructed images by ADSIR, GDSIR,MGDSIR, and FMGDSIR, respectively. In this case, since
a lot of projection data (𝐼 = 173,700 ≫ 𝐽) is available withsmall noise, we can see that FBP reconstructs a high-quality
image. Therefore, we can confirm that 𝐺 = 300 and 𝑏 = 106
are proper values as the setting of a standard-dose situation.
We can also see that all SIR methods reconstruct even better
results than that by FBP with mostly the same quality.

6.3

Low-dose projection data were generated in two situations:(i) reducing projection views from 𝐺 = 300 to 𝐺 = 60, and(ii) lowering X-ray intensity from 𝑏 = 106 to 𝑏 = 2.5 × 104 .
Note that there is a significant difference between these two
approaches. The first approach reduces the number of obser-

Standard-Dose Simulation Results

We first reconstructed ‘m_1132’ from standard-dose projection data in 𝐺 = 300 views, where we set 𝜃 𝑔 := 0.6(𝑔−1) deg
and the X-ray radiation intensity to 𝑏 = 106 . We set the regu-

6.4

Low-Dose Simulation Results

IEICE TRANS. FUNDAMENTALS, VOL.E104–A, NO.4 APRIL 2021

708(a) Ground truth (f_1615).

(b) FBP [32.92 / 0.893].

(c) ADSIR [37.81 / 0.983].

(d) GDSIR [37.66 / 0.983].
(e) MGDSIR [38.34 / 0.984].
(f) FMGDSIR [38.41 / 0.985].
Fig. 4 Reconstruction result of f_1615 from its low-dose projection data (60/300 views) by each method [PSNR / SSIM].

(a) Ground truth (m_1132).

(b) FBP [31.66 / 0.858].

(c) ADSIR [36.76 / 0.978].

(d) GDSIR [36.67 / 0.978].
(e) MGDSIR [37.24 / 0.981].
(f) FMGDSIR [37.20 / 0.981].
Fig. 5 Reconstruction result of m_1132 from its low-dose projection data (60/300 views) by each method [PSNR / SSIM].

KAMOSHITA et al.: MULTICLASS DICTIONARY-BASED STATISTICAL ITERATIVE RECONSTRUCTION FOR LOW-DOSE CT

709(a) Ground truth (f_1615).

(b) FBP [30.48 / 0.525].

(c) ADSIR [39.07 / 0.985].

(d) GDSIR [39.36 / 0.986].
(e) MGDSIR [39.97 / 0.987].
(f) FMGDSIR [39.36 / 0.987].
Fig. 6 Reconstruction result of f_1615 from its low-dose projection data (1/40 intensity) by each method [PSNR / SSIM].

(a) Ground truth (m_1132).

(b) FBP [29.92 / 0.525].

(c) ADSIR [38.00 / 0.980].

(d) GDSIR [38.20 / 0.981].
(e) MGDSIR [38.68 / 0.982].
(f) FMGDSIR [38.18 / 0.981].
Fig. 7 Reconstruction result of m_1132 from its low-dose projection data (1/40 intensity) by each method [PSNR / SSIM].

IEICE TRANS. FUNDAMENTALS, VOL.E104–A, NO.4 APRIL 2021

710(a) f_1615.

(a) f_1615.

(b) FBP.

(c) ADSIR.
(d) GDSIR.
(e) MGDSIR.
Enlarged images of each reconstruction result in Fig. 5.

(f) FMGDSIR.

Fig. 8(c) ADSIR.
(d) GDSIR.
(e) MGDSIR.
Enlarged images of each reconstruction result in Fig. 7.

(f) FMGDSIR.

Fig. 9(b) FBP.

vations (𝐼 = 34,740 < 𝐽), while the second approach does
not reduce the number of observations (𝐼 = 173,700 ≫ 𝐽)
but enlarges the variance, 40 times, of 𝑧 𝜃𝑔 ,𝑛 /𝑏 = exp(−𝑙 𝜃𝑔 ,𝑛 ).
The reduction of the number 𝐼 of the observations in the first
approach makes the reconstruction problem in (2) more difficult than the increase of the noise level ∥𝝃 ∥ 2 in the second
approach. If we set a same reduction ratio of the X-ray dosefor both approaches, then the reconstruction quality is low in
the first approach or difference of the quality is not visible inthe second one. Hence, we used the different reduction ratios
1/5 and 1/40 in the first and second approaches, respectively.
1) 𝐺 = 60 and 𝑏 = 106 : We set 𝜃 𝑔 := 3(𝑔 − 1) deg andthe regularization parameters to 𝜆 = 3,800 in ADSIR and
GDSIR, to 𝜆1 = 7,500, 𝜆2 = 3,800, 𝜆 3 = 1,000, 𝜆4 = 2,500,
𝜆5 = 1,000 in MGDSIR, and to 𝜆1 = 7,500, 𝜆2 = 6,000, 𝜆 3 =1,000, 𝜆4 = 1,500, 𝜆5 = 1,000 in FMGDSIR. When 𝝁 (0) was
created by FBP, the projection data were linearly interpolated
to 300 views. Table 1 summarizes PSNR of the reconstructed
images and the calculation time per iteration for each method.
From Table 1, we see that MGDSIR and FMGDSIR achieved
higher reconstruction accuracy than by ADSIR and GDSIR
with the same regularization parameters for all the images.

In addition, FMGDSIR reduced the calculation time approximately 46% from MGDSIR without accuracy degradation.
2) 𝐺 = 300 and 𝑏 = 2.5×104 : We set the regularization
parameters to 𝜆 = 800 in ADSIR and GDSIR, to 𝜆1 = 800,
𝜆2 = 800, 𝜆3 = 400, 𝜆 4 = 900, 𝜆5 = 400 in MGDSIR, and to
𝜆1 = 2,000, 𝜆2 = 1,300, 𝜆3 = 800, 𝜆4 = 1,100, 𝜆5 = 800 in
FMGDSIR, which are constant multiples of 𝜆 𝑞 in Sect. 6.3
since the matrix 𝑅 is the same and the noise level is differentin (2). Table 2 shows PSNR and the calculation time for eachmethod. From Table 2, we see that MGDSIR achieved thehighest reconstruction accuracy. FMGDSIR reduced the calculation time by about 30% from MGDSIR with almost thesame accuracy as with the conventional GDSIR. In our simulations, FMGDSIR performed more effectively in the view
reduced situation than in the weakened intensity situation.
In Figs. 4–7, (a), (b), (c), (d), (e), and (f) show the original image, the initial image 𝝁 (0) by FBP, the reconstructed
images by ADSIR, GDSIR, MGDSIR, and FMGDSIR, respectively, in the two low-dose simulations for ‘f_1615’ and‘m_1132’, where Figs. 4 and 5 show the results in the view reduced situation, and Figs. 6 and 7 show ones in the weakened
intensity situation. From these figures, we see that each ini-

KAMOSHITA et al.: MULTICLASS DICTIONARY-BASED STATISTICAL ITERATIVE RECONSTRUCTION FOR LOW-DOSE CT

711
Table 3
𝑃
4×4
6×6
8×8
10 × 10
12 × 12

4×4

Settings of 𝑃, 𝑄 and 𝜆𝑞 in FMGDSIR and the reconstruction accuracy for m_1132 in the view reduced situation.
(𝜆1 , 𝜆2 , . . . , 𝜆𝑄 )𝑄

537.20 / 0.981
37.05 / 0.980
36.91 / 0.980
36.78 / 0.979
36.62 / 0.978

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

4200
5000 × (1.5, 0.8)5000 × (1.5, 1.0, 0.2)
5000 × (1.5, 1.2, 0.2, 0.3)
5000 × (1.5, 1.2, 0.2, 0.3, 0.2)
5000 × (1.5, 1.2, 0.2, 0.2, 0.3, 0.3)
5000 × (1.5, 1.2, 0.3, 0.2, 0.2, 0.2, 0.2)
5000 × (1.5, 1.2, 0.3, 0.2, 0.2, 0.2, 0.2, 0.2)
5000 × (1.5, 1.2, 0.3, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2)
5000 × (1.5, 1.2, 0.3, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2)
5000 × (1.5, 1.2, 0.3, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2)
5000 × (1.5, 1.2, 0.3, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2)
5000 × (1.5, 1.2, 0.2, 0.2, 0.2, 0.2, 0.3, 0.2, 0.3, 0.3, 0.2, 0.3, 0.2)
5000 × (1.5, 1.2, 0.2, 0.2, 0.2, 0.3, 0.2, 0.3, 0.3, 0.2, 0.3, 0.2, 0.3, 0.2)
5000 × (1.5, 1.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3)
5000 × (1.5, 1.2, 0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.2, 0.2, 0.3, 0.3, 0.2, 0.2, 0.2)

36.18 / 0.976
36.23 / 0.976
36.72 / 0.979
37.08 / 0.980
37.20 / 0.981
37.22 / 0.981
37.25 / 0.981
37.23 / 0.981
37.27 / 0.981
37.26 / 0.981
37.23 / 0.981
37.27 / 0.981
37.24 / 0.980
37.24 / 0.981
37.26 / 0.981
37.25 / 0.981(a) Relation between the accuracy and the patch size.
Fig. 10(b) Relation between the accuracy and the number of classes.

Values of PSNR (⃝) and SSIM (∗) in terms of the patch size 𝑃 and the number 𝑄 of classes.

tial estimate has serious artifacts or noise while all SIR methods reconstruct high-quality images. Figures 8 and 9 show
enlarged images of Figs. 5 and 7, and the red arrows indicatethat detailed structures are lost in ADSIR and GDSIR while
the proposed methods can reconstruct them more precisely.
For example, a light-gray vertical line is visible clearly at the
tip of the left arrow in the upper image of Fig. 9(e) while thelines in Figs. 9(c) and 9(d) are blurred. Such detailed structures are often used for medical diagnosis, and a radiographer
in a hospital supported the results of the proposed methods.

the pixel length is smaller than 1 mm, the optimal 𝑃 will belarger. If the target is not a head, we have to adjust 𝑃 again.
2) Number of classes: We used 𝑄 = 1, 2, . . . , 16. Figure 10(b) shows that PSNR and SSIM largely increased until
𝑄 = 5. After that, the accuracy vibrated and the image quality hardly changed. Hence, we have chosen 𝑄 = 5 considering the difficulty in tuning of the regularization parameters
𝜆 𝑞 , although we can see that, in Table 3, FMGDSIR keeps
the reconstruction accuracy with the rough parameter tuning.
7.

6.5

PSNR / SSIM

5000 × (1.5, 1.2, 0.2, 0.3, 0.2)
2500 × (1.5, 1.2, 0.2, 0.3, 0.2)
1500 × (1.5, 1.2, 0.2, 0.3, 0.2)
800 × (1.5, 1.2, 0.3, 0.2, 0.2)
600 × (1.5, 1.2, 0.3, 0.2, 0.2)

Conclusion

On the Patch Size 𝑃 and the Number 𝑄 of Classes

To analyze the effects of 𝑃 and 𝑄, we reconstructed ‘m_1132’
by FMGDSIR in the view reduced situation, using several 𝑃
and 𝑄 with a rough parameter tuning for 𝜆 𝑞 as in Table 3.
√
1) Patch size: We used 𝑃 = 4, 6, . . . , 12. Figure 10(a)shows PSNR and SSIM in terms of the patch size. We foundthat the reconstruction accuracy decreased as the patch size
increased. Hence, we have chosen 𝑃 = 16 as the patch size.
Note that the width and height of each pixel were 1 mm and
𝑃 = 16 was optimal for this pixel length. We expect that if

In this paper, we proposed low-dose CT reconstruction methods exploiting multiclass dictionary learning. We used multiclass dictionaries and multiclass regularization parameters.
The proposed methods, MGDSIR and FMGDSIR, achieved
higher reconstruction accuracy than the conventional singleclass dictionary-based methods, ADSIR and GDSIR. In addition, since MGDSIR performs clustering for target image
patches only once, computational cost hardly increases fromGDSIR. Moreover, by introducing orthogonal dictionaries,
FMGDSIR achieved fast reconstruction. Simulation results

IEICE TRANS. FUNDAMENTALS, VOL.E104–A, NO.4 APRIL 2021

712

showed that the proposed methods reconstruct detailed structures of the targets while suppressing artifacts and noise.
Acknowledgments
This work was supported by JSPS KAKENHI (17H02069).
The authors are very grateful to Mr. Hiroki Okada in Kusatsu
General Hospital for his valuable advice as a radiographer.

[17] J. Mairal, J.P. Chieze, Y. Chen, G. Durif, R. Jenatton, F. Bach,
J. Ponce, G. Obozinski, B. Yu, G. Sapiro, and Z. Harchaoui, Spams: a
SPArse Modeling Software ver. 2.6, http://spams-devel.gforge.inria.
fr/doc/html, accessed Jan. 20, 2020.
[18] Magnetic Resonance Research Facility, Visible Human Project CT
Datasets, https://medicine.uiowa.edu/mri/facility-resources/images/
visible-human-project-ct-datasets, accessed Jan. 20, 2020.

References[1] H. Kamoshita, T. Shibata, D. Kitahara, K. Fujimoto, and A. Hirabayashi, “Low-dose CT reconstruction with multiclass orthogonal
dictionaries,” Proc. IEEE Int. Conf. Image Process. (ICIP), Taipei,
Taiwan, pp.2055–2059, Sept. 2019.
[2] A.C. Kak and M. Slaney, Principles of Computerized Tomographic
Imaging, SIAM, Philadelphia, PA, 2001.
[3] J. Hsieh, Computed Tomography: Principles, Design, Artifacts, andRecent Advances, 3rd ed., SPIE Press, Bellingham, WA, 2015.
[4] Q. Xu, H. Yu, X. Mou, L. Zhang, J. Hsieh, and G. Wang, “Low-dose
X-ray CT reconstruction via dictionary learning,” IEEE Trans. Med.
Imag., vol.31, no.9, pp.1682–1697, Sept. 2012.
[5] Y. Chen, L. Shi, Q. Feng, J. Yang, H. Shu, L. Luo, J.L. Coatrieux, andW. Chen, “Artifact suppressed dictionary learning for low-dose CT
image processing,” IEEE Trans. Med. Imag., vol.33, no.12, pp.2271–
2292, Dec. 2014.
[6] M. Elad, Sparse and Redundant Representations: From Theory to
Applications in Signal and Image Processing, Springer, New York,NY, 2010.
[7] Y. Chen, X. Yin, L. Shi, H. Shu, L. Luo, J.L. Coatrieux, and C. Toumoulin, “Improving abdomen tumor low-dose CT images using afast dictionary learning based processing,” Phys. Med. Biol., vol.58,
no.16, pp.5803–5820, Aug. 2013.
[8] C. Bao, J.F. Cai, and H. Ji, “Fast sparsity-based orthogonal dictionary
learning for image restoration,” Proc. IEEE Int. Conf. Comput. Vis.
(ICCV), Sydney, Australia, pp.3384–3391, Dec. 2013.
[9] C. Zhang, T. Zhang, M. Li, C. Peng, Z. Liu, and J. Zheng, “Low-dose
CT reconstruction via ℓ1 dictionary learning regularization using
iteratively reweighted least-squares,” Biomed. Eng. Online, vol.15,
no.66, 21 pages, June 2016.
[10] K. Sauer and C. Bouman, “A local update strategy for iterative reconstruction from projections,” IEEE Trans. Signal Process., vol.41,
no.2, pp.534–548, Feb. 1993.
[11] J.A. Fessler, “Statistical image reconstruction methods for transmission tomography,” Handbook of Medical Imaging, Volume 2:
Medical Image Processing and Analysis, eds. M. Sonka and J.M.
Fitzpatrick, pp.1–70, SPIE Press, Bellingham, WA, 2000.
[12] I.A. Elbakri and J.A. Fessler, “Statistical image reconstruction forpolyenergetic X-ray computed tomography,” IEEE Trans. Med.
Imag., vol.21, no.2, pp.89–99, Feb. 2002.
[13] Y.C. Pati, R. Rezaiifar, and P.S. Krishnaprasad, “Orthogonal matching pursuit: Recursive function approximation with applications towavelet decomposition,” Proc. Asilomar Conf. Signal Syst. Comput.
(ACSSC), Pacific Grove, CA, pp.40–44, Nov. 1993.
[14] J. Mairal, F. Bach, J. Ponce, and G. Sapiro, “Online learning for
matrix factorization and sparse coding,” J. Mach. Learn. Res., vol.11,
pp.19–60, Jan. 2010.
[15] Z. Zhan, J.F. Cai, D. Guo, Y. Liu, Z. Chen, and X. Qu, “Fast multiclass dictionaries learning with geometrical directions in MRI reconstruction,” IEEE Trans. Biomed. Eng., vol.63, no.9, pp.1850–1861,
Sept. 2016.
[16] N. Inamuro, M. Shibata, C. Tang, T. Ijiri, and A. Hirabayashi, “Highquality MR imaging based on dictionary learning using training
images and observed signals,” IEICE Technical Report, vol.116,
no.395, pp.123–128, Jan. 2017 (in Japanese).

Hiryu Kamoshita
received his B.E. degreeand M.E. degree in computer science from Ritsumeikan University, Shiga, Japan, in 2018 and2020, respectively. His research interest includes
compressed sensing and its applications. He received the Outstanding Research Award from the
Graduate School of Information Science and Engineering, Ritsumeikan University in 2020.

Daichi Kitahara
received his B.E. degree
in computer science in 2012, and his M.E. and
Ph.D. degrees in communications and computerengineering in 2014 and 2017 from the Tokyo Institute of Technology, Tokyo, Japan, respectively.
Since 2017, he has been an Assistant Professor
with the College of Information Science and Engineering, Ritsumeikan University, Shiga, Japan.
His current research interest includes signal processing and its applications, inverse problem, optimization theory, and multivariate spline theory.

Ken’ichi Fujimoto
received his M.E. degree from the Tokyo Institute of Technology in1997 and his Ph.D. degree from the University ofTokushima in 2000. He is currently an Associate
Professor at Kagawa University. His research interest spans interdisciplinary basic research on
the dynamics of nonlinear dynamical systems
and its applications to medical engineering as tomographic image reconstruction. He is a member of IEEE, IEICE, IEEJ, SICE, and ISCIE.

Laurent Condat
received his Ph.D. degree
in applied mathematics from Grenoble Instituteof Technology in 2006. From 2006 to 2008, he
was a postdoc in the Helmholtz Zentrum München in Munich, Germany. He was hired as a permanent researcher by the French National Center
for Scientific Research (CNRS) in 2008. From2016 to 2019, he was a member of the French National Committee for Scientific Research. Since
Nov. 2019, he is a Research Scientist in the Visual Computing Center of King Abdullah University of Science and Technology (KAUST), Saudi Arabia. His area of interest includes convex optimization models and algorithms, and signal andimage processing. He received a best student paper award at the IEEE ICIP
2005 and a best Ph.D. award from Grenoble Institute of Technology in 2007.

KAMOSHITA et al.: MULTICLASS DICTIONARY-BASED STATISTICAL ITERATIVE RECONSTRUCTION FOR LOW-DOSE CT

713

Akira Hirabayashi
received his D.E. degree
in computer science in 1999, from the Tokyo Institute of Technology, Tokyo, Japan. Currently,
he is a Professor of the College of Information
Science and Engineering, Ritsumeikan University, Shiga, Japan. His research interests include
compressed sensing, sampling theory, and signal and image processing. Prof. Hirabayashi isalso a member of the Institute of Electrical andElectronics Engineers (IEEE) and the AcousticalSociety of Japan (ASJ).