計算数理工学論文集 Vol. 21 (2021 年 12 月), 論文 No. 05-211218

JASCOME

部分観測による X 線計算機断層撮影法の数値的実現
NUMERICAL COMPUTATION OF X-RAY COMPUTERIZED TOMOGRAPHYFROM PARTIAL MEASUREMENT藤原宏志 1) ，大石直也 2) , SADIQ, Kamran3) , TAMASAN, Alexandru4)

Hiroshi FUJIWARA, Naoya OISHI, Kamran SADIQ, and Alexandru TAMASAN
1) 京都大学大学院情報学研究科(〒 606-8501

京都市左京区吉田本町,

E-mail: fujiwara@acs.i.kyoto-u.ac.jp)

2) 京都大学大学院医学研究科(〒 606-8507

京都市左京区聖護院川原町,

E-mail: noishi@kuhp.kyoto-u.ac.jp)

3)Johann Radon Institute (RICAM)(〒 4040Linz, Austria,

E-mail: kamran.sadiq@ricam.oeaw.ac.at)4)University of Central Florida(〒 32816

Orlando, Florida, USA,

E-mail: tamasan@math.ucf.edu)

The present paper proposes a novel numerical scheme to X-ray Computerized Tomography from partial measurement data. In order to reduce radiation exposure, it is desirable
to irradiate X-ray only around region of interest, while the conventional reconstruction
methods such as ﬁltered back projection could not work due to its intrinsic limitation of
dependency on whole measurement data. The proposed method gives a direct numericalreconstruction employing a Cauchy type boundary integration in A-analytic theory and asingular integral equation which maps boundary measurement to interior data. Numerical examples using experimental data are also exhibited to show validity of the proposed
numerical procedure.
Key Words : Inverse Problems, X-ray Computerized Tomography, Partial Measurement,
Direct Reconstruction, Boundary Integration, Singular Integral Equations
本論文では，部分観測での X 線 CT に対し，A-analytic

1. 緒言
本論文では，部分観測からの X 線計算機断層撮影法 (Com-理論に基く直接解法の数値計算アルゴリズムを提案する．本puterized Tomography; CT) において，近年提案された Radon

手法は，X 線 CT の減衰係数の再構成を，著者らが輸送方程逆変換とは独立な直接再構成手法に対する数値計算スキーム

式を対象として近年提案した，境界の一部での観測から非斉を提案する．また実測データをもちいた再構成により，その次項を再構成する逆問題 (inverse source problem) (7) に帰着

妥当性を示す．

させる．この逆問題解法は，逆 Radon 変換とは独立であり，一般的な X 線 CT では，物体の断層面を通過する全てのA-analytic 函数 (2) の境界値の特徴付け (14) および Cauchy

直線に沿って X 線の投影データ (すなわち，外部から入射さの積分公式に相当する境界積分 (2, 3) が中心的な役割を果れ，物体で減衰を受けて出射される X 線の強度) を測定し，たす．数値計算の視点からは，この境界積分が領域の内部に断層面全体の X 線の減衰率を得る．その数理的再構成手法

メッシュを必要とせず，減衰率が各点で独立に計算できるたに基く Filtered Back Projection

め並列計算に適するという特徴をもつ (9) ．さらにこの手法(FBP) や Algebraic Reconstruction Techniques (ART) が確

は，回転角度が限定された撮影から断層面を再構成するトモとしては逆 Radon 変換
立されており(10)(13)，医用および産業用の CT 実機に実装されシンセシス (15) と呼ばれる近年登場した X 線撮影技術にも適用可能と考えられる．提案手法では，結果の精度と安定性て，広く利用されている．
一方，X 線の被曝量の低減と検査時間の短縮のため，物体を同時に制御するパラメータが現れ，その選択が重要とな中の関心領域 (Region of Interest; ROI) の近傍を通過する Xる．しかし対象の構造，観測誤差，もちいる離散スキームや線のみの照射が望ましい．この部分的な観測からの再構成と離散化数にも依存し，事前に解析的に最適値を決定するこ

して，逆 Radon 変換に基く直接解法や，圧縮センシング理とは困難と考えられるため，数値計算等を利用する事後 (a論に基く反復解法が提案されている(12)

posteriori) 解析による選択基準を提案する．

．

次節で，部分観測からの X 線 CT を輸送方程式に帰着さ2021 年 9 月 30 日受付，2021 年 11 月 1 日受理− 37 −

力値 I0 で正規化した

yu(z, ξ) = − log

Λ
D

−s

+

ROI

X-ray

•
z

(1)

x
ξ · ∇z u(z, ξ) = µ(z),

λ0

c•
θ

(z, ξ) ∈ D × S 1

は次の輸送方程式を満たすことが知られている．

sSξI(z, ξ),I0u(ζ, ξ) = 0,

ζ − icΓ−(z, ξ) ∈ D × S 1 ,(2)(ζ, ξ) ∈ Γ− .

(3)

た だ し ∇z = (∂x , ∂y ), Γ± = {(ζ, ξ) ∈ ∂D × S 1 ; ±(ζ −•ic) · ξ > 0} である．ここで D × S 1 の集合としての境界は，

Γ+ζ

Γ0 = {(ζ, ξ) ∈ ∂D × S 1 ; (ζ − ic) · ξ = 0} として ∂D × S 1 =Γ− ∪ Γ+ ∪ Γ0 (disjoint union) であることに注意する．ただFig. 1: Introduced coordinate system of R2 for

し Γ0 に沿う X 線の軌跡は D と共通部分をもたない．また

Γ− , Γ+ はそれぞれ X 線の入射と出射に対応する．Γ− でのdescribing the proposed reconstruction proce-

境界条件 (3) は入射に対応しており，u は I0 で正規化してdure. Region of interest (ROI) is surrounded

いるため 0 となる．一方，X 線 CT における代表的な平行

by arc Λ and chord S. The measurement is

ビーム投影において，観測値 (投影データ) では，D の中心

performed only on the arc Λ. Note that the

対する物体からの出射強度を，方向ベクトル ω = ξ ⊥ ∈ S 1 のlength of ξ ∈ S 1 is depicted at reduced scales

投影面への正射影として観測した値を Radon 変換 Rµ(ω, s)

ic から距離 s で ξ 方向に直進する X 線を考え，I0 の入射にとして与える．この投影データは出射強度の観測であるため

for clarity.

Γ+ で与えることになり，u(ζ, ξ) = Rµ(ξ ⊥ , (ζ − ic) · ξ ⊥ ),せて A-analytic 理論を適用する数学解析的手法を述べ，3 節
でその数値計算手法を述べる．4 節では実測データをもちいた数値計算例によって提案手法の妥当性を示し，5 節で精度と安定性を制御するパラメータの選択を論じる．

(ζ, ξ) ∈ Γ+

(4)が成立する．
順問題の設定 では，D で µ を既知として ，境界値 問題(2),(3) の適切性を論じる．また一般的な X 線 CT は，D を
2. 再構成公式

通る全ての直線に沿う投影データを観測するため，∂D 全体

本節では著者らの先行研究 (7) の inverse source problem

での出射，すなわち Γ+ 全体で (4) の Rµ を与えて (2), (3),に基き，X 線 CT の減衰係数の再構成の数学解析的な手順(4) を満たす µ を D で求める逆問題となる (2, 3, 9) ．一方，を示す．

本論文では，観測値が Γ+ 全体ではなく，Λ ⊂ ∂D でのみ与対象物体の断層画像を取得する領域 (ROI) を含む断面をえられる部分観測からの µ の再構成を考える．すなわち µ

D とする．対象物体の外部では X 線の吸収は生じないとし，

は D で未知であり，(4) の Rµ は ζ ∈ Λ, (ζ − ic) · ξ ⊥ > 0 で適当なスケーリングによって D は半径 1 の円の内部としてのみ既知とし，(2), (3), (4) を満たす µ を (D+ 上で) 再構成

一般性を失わない．提案手法では ∂D の空でない連結部分する逆問題を考える．

集合 (弧) Λ からのみ出射される X 線強度を観測し，Λ およ

さて ∂ = (∂x − i∂y )/2, ∂ = (∂x + i∂y )/2 とすると ξ · ∇z = Pe ∂ + eiθ ∂ であり，L2 で u z, ξ(θ) = m um (z)eimθ とすびその端点を結ぶ線分 S で囲まれる領域 D+における X 線の減衰率 µ を再構成する．ここで µ ∈ C01,α (D), α > 1/2 とし，ROI が D+−iθると，(2) は {eimθ ; m ∈ Z} の直交性と µ の等方性よりに含まれるように Λ をとる．

線分 S の中点を原点とし，S に沿って x 軸をとる．また，∂u1 + ∂u−1 = µ,

Λ が R2 の上半平面に位置するように y 軸をとる．線分 S の∂u−m + ∂u−m−2 = 0,

長さを 2s として，D の中心 (0, c) を原点とする極座標にお√ける S の端点 (s, 0) の偏角を λ0 とする．以下では i = −1
を虚数単位として (x, y) ∈ R2 と z = x + iy ∈ C を同一視し，ζ は主として D の境界 ∂D もしくはその一部 Λ の点を表すものとする．また X 線の進む方向 ξ ∈ S 1 を極座標ξ = (cos θ, sin θ) によって θ ∈ [0, 2π) と同一視する (Fig. 1)．
入射する強度を一定値 I0 > 0 とし，位置 z ∈ D において

ξ ∈ S 1 方向に進む X 線の強度を I(z, ξ) とするとき，I を入− 38 −(5)−m = 0, −1, −2, . . .

(6)となる．ここで区分的に滑らかな閉曲線 γ ⊂ D に対し，γ の内部に含まれる点 z に対して，u−1 (z) は γ 上での複素積分

Z1u−1 (ζ)dζ
2πi γ ζ − z

 X
Z 
∞
dζ
1
dζ
+
−
u−1−2j (ζ)
2πi γ ζ − zζ̄ − z̄ j=1
u−1 (z) =

ζ −z
ζ −z

!j 


(7)

で与えられる (2, 3) ．全周 ∂D で観測する場合，(7) で γ = ∂Dで固定し，m は −M ≤ −m ≤ 0 を満たす全ての整数，またとすると，右辺に現れる u−m |∂D は (3) および Γ+ での観測

−m − 2j ≥ −M を満たす最小の整数 j を J = J(m) で表す．

値 (4) から求まり，(5) より D 上での µ が求まる．
しかし本論文で扱う部分観測の場合，u−m (ζ) は ∂D の一部 Λ 上でしか与えられず，(7) を直接利用することはできない．そこで S が実軸に含まれることから，(6) を満たす

u−m |S が特異積分方程式(7, 14)(−m)

Uk．ただし p.v. を Cauchy の

1
p.v.
π

1
πi

Z s

u−m (t)dt,−s x − t

1 Xu(ζk , ξn )eimθn ,N n=0

1 ≤ k ≤ K.

p−m の第 2 項の複素積分は

Z 
dζ
dζ
−
F (ζ)
ζ −x
ζ −x
Λ
!
Z π−λ0
ieiλ dλ−ie−iλ dλ=−F ζ(λ)
ζ(λ) − xζ(λ) − xλ0Z π−λ0
iλ

e
=2i
Re
F ζ(λ) dλζ(λ) − xλ0|x| < sZu−m (ζ)dζΛ ζ −x
j )
 (X

Z 
∞
1
dζ
dζ
ζ −x
+
−
u−m−2j (ζ)πi Λ ζ − x
ζ −x
ζ −x
j=1

であることに注意して，である．
ここで作用素 iHs の L2 (S) 上でのスペクトルは −1 ≤ ξ ≤ 1
であるが=次に，Λ = {ζ(λ) = eiλ + ic ; λ0 < λ < π − λ0 } と表すと，およびp−m (x) =(ζ − ic) · ξ ⊥ > 0

N −1|x| < s主値積分としてHs u−m (x) =(ζ − ic) · ξ ⊥ ≤ 0;から ζk ∈ Λ における Fourier 係数 u−m (ζk ) を次で近似する．

(I − iHs )u−m (x) = 2p−m (x),を満たすことに着目する

数値計算の手順としては，まず0,u(ζ, ξ) =Rµ(ξ ⊥ , (ζ − ic) · ξ ⊥ ),(11)

，端点での I −iHs は単射であり

(7)

−1

，(I −iHs )は自然な定義域で一意に存在することに注意する．そこでこの解 u−m |S をもちい，γ = Λ ∪ S として (7) をもちいることで u−1 |D+ が求まる．さらに u は実数値ゆえ u1 = u−1 で

1 X Ukeiλk ∆λπζk − zk=1 J(−m)
j 

K
X
2X
ζ
−
z
eiλk
k
(−m−2j)+URe∆λπ
ζk − z  j=1 kζk − z K

k=1とし，積分方程式を選点法で離散化して，連立一次方程式あり，(5) より z ∈ D+ において


µ = ∂u−1 + ∂u−1 = Re ∂x u−1 + Im ∂y u−1

(−m)

Vℓ

−

(−m)

3. 提案手法の数値的実現を得る．ここで Vℓ

応用の視点からは，前節で述べた数学解析的な再構成手法の数値的手法が重要となる．そこで本節では，前節の再構成
手法を数値的に実現するため，現れる線積分を (複合) 中点則で近似するアルゴリズムを提案する．
離散化にあたり，K, L, N , M をそれぞれ正整数とする．
そこで Λ を K 等分して

π − 2λ0,K

1
λk = λ0 + k −2(8)は S 上での u−m (xℓ ) 相当値であり，数値解析学の立場からは，計算が停止せず数値解が唯一つ求まるために，次が重要である (8) ．
定理 1. 任意の正整数 L に対して，連立一次方程式 (8) に現れ (−m)る係数行列は (狭義) 正定値である．特に Vℓ;1≤ℓ≤L
は唯一つ存在する．

 (−m)と Vℓ, −M ≤ −m ≤ 0 を

もちいると


∆λ,

1 ≤ k ≤ K,

1 X Vℓ∆x2πixℓ − zℓ=1 J(−1)
L 
X (−1−2j)  xℓ − z j 
1 X
∆x
∆x
+
−
Vℓ
2πi
xℓ − zxℓ − z xℓ − z L(−1)

U (−1) (z) = P−1 (z) +とする．また S を L 等分して
2s
1
, xℓ = −s + ℓ −
∆x,
∆x =
L
2

j=1

ℓ=1

1≤ℓ≤L

が z ∈ D+ で求まる．そこで，充分小さい正数 h1 , h2 を定とする．X 線の進行方向は

θn = n∆θ,1≤ℓ≤L (−m)
以上で得られる Uk

ζk = eiλk + ic ∈ Λ

2π
∆θ =
,
N

i X
∆x
(−m)
V
= 2P−m (xℓ ),πxℓ − xj j
1≤j<L
j̸=ℓ

によって減衰率 µ が求まる．

∆λ =

(−m)

P−m (z) =めて，中心差分近似により

ξn = (cos θn , sin θn ) ∈ S 1

U (−1) (z + h1 ) − U (−1) (z − h1 )2h1
U (−1) (z + ih2 ) − U (−1) (z − ih2 )
+ Im
2h2

µ(z) ≈ Re

で離散化し，各 ζk ∈ Λ で (ζk − ic) · ξn⊥ > 0 の範囲の方向ξn で Rµ が与えられているとする．さらに M をひとつ選ん− 39 −

の，D の中心 ic を通過する X 線との距離を表しており，こ1の図の範囲では [−1, 1] 区間を 1887 分割している．

FBP (10, 13) で得られる再構成を Fig. 3 に示す．FBP で

0.5

-0.5

(ζ-ic)•ξ⊥0

典型的な streak アーチファクトを含むものの，D の画像が得られている．

1
√3/2
0.5

π

-1
0

arg ξ⊥
Fig. 2:-0.5

Measurement data (sinogram) bywhole measurement on Γ+

π/2

0
を得る．

π

(ζ-ic)•ξ⊥

π/2

0

-√3/2
-1

arg ξ⊥

4. 実測値からの再構成例

Fig. 4:

提案手法の妥当性を示すため，Finnish Inverse Problems

Sinogram corresponding to partial

Society が公開している実測値のひとつ (1, 4) をもちいた再measurement data on arc Λ = {π/6 < arg(ζ −構成例を示す．

ic) < 5π/6}.
FBP

本節では部分観測をおこなう Λ として π/6 < arg(ζ −ic) <1

5π/6 の部分を考える．このとき，Fig. 1 のとおり座標軸を√設定すると，c = −0.5, λ0 = π/6, s = 3/2 となる．部分観測に対応する観測値を Fig. 4 に示す．図示していない部

0.5

分は，投影データを取得する必要のない部分を示しており，Fig. 2 と比較して，観測回数が有意に減少していることが0わかる．例えば ξ = 0 (x 軸に平行に X 線を照射) の場合，Fig. 5 の左図に示すとおり，D+ を通過するのは ic との距
離 s が 0.5 < s < 1 の範囲である．また ξ = π/2 (y 軸に平-0.5

行に X 線を照射) の場合，Fig. 5 の右図に示すとおり，D+√を通過するのは |s| < 3/2 の範囲のみとなる．

FBP は，Pω (s) = Rµ(ω, s), s = (ζ − ic) · ξ ⊥ と |r| の逆

-1
-1

-0.5

0

0.5

1

s

Fig. 3: Reconstructed image by conventional

FBP from whole measurement in Fig. 2. TheD+

sub domain y > 0.5 (above the dotted line)

ic •

corresponds to ROI in following numerical ex-

X-ray
ξ=0

periments.

1
0.5
•O

s

√
− 3/2

√

O
•

3/2

D+
ic •
X-ray
ξ = π/2

Fig. 5: Parallel beam projection in partial
Fig. 2 に，公開されている投影データ (Radon 変換像) の実測値の一部を示す．これは部分観測ではなく全周観測のも

measurement. The axis s is the distance be-

ので，公開されているデータから物体近傍のみを取り出した

tween X-ray and ic. (left) X-ray’s direction is

ものである．横軸は X 線の照射 (進行) 方向に垂直な ξ⊥の偏角であり，このデータは，方向の区間 [0, π) を 180 等分方向して提供されている．また縦軸は，平行に照射した X 線− 40 −

ξ = 0; (right) X-ray’s direction is ξ = π/2

Fourier 変換の適当な近似 h の合成積によって
Z 2π
1µ(z) ≈Pω ∗ h(z · ω) d(arg ω)
2π 0

も実数値であることに注意する．さらに式 (7) と同様にで表されるが，本論文で考える部分観測では Pω (s) が s について局所的にのみ与えられるため，任意の z ∈ D+ で合成
積 Pω ∗ h が計算されず，µ が求まらない．

0.25
0
0



して

0.5

-0.5

j=1

!j 


ζ −z
ζ −z

が成立し，前節の離散化アルゴリズムに対応する近似値と

M=309

-1

Z
1
u0 (ζ)dζ
2πi Λ∪S ζ − z


 X
Z
∞
dζ
1
dζ
+
−
u−2j (ζ)
2πi Λ∪S ζ − zζ̄ − z̄ 
u0 (z) =0.5

1 X Vℓ∆x2πixℓ − zℓ=1 M/2L 
X (−2j)  xℓ − z j 
∆x
1 X
∆x
−
+
Vℓ
2πi
xℓ − zxℓ − z xℓ − z L(0)

U (0) (z) = P0 (z) +

j=1

ℓ=1

1

が考えられる．

Fig. 6: Reconstructed image in D + from par-上述のとおり u0 ∈ R であるが，この近似値 U (0) は実数

tial measurement in Fig. 4 on Λ with M =値とは限らない．そこで，虚部 Im U (0) が誤差に対応すると

309.

考え，Im U (0) を最小とする M が打ち切りパラメータの候補と考える (6) ．

一方，提案する手法による再構成の結果を Fig. 6 に示す．

0.0136FBP の再構成 Fig. 3 の y > 0.5 の部分と同様の再構成結
0.0135果が得られており，提案手法の妥当性を示している．ここで Fig. 4 の観測は，Λ を K = 1977 等分することに相当し，となるよう，L = 1634 とした．また h1 = h2 = 1/200 として，格子点 Dh+ = {zj = (j1 h1 , j2 h2 ) ∈ D+ ; j = (j1 , j2 ) ∈ Z2 } 上

0.0134|| Im U(0) ||2(∆λ = (2π/3)/1977 ≈ 2/1887)，L 上の分点数は ∆x ≈ ∆λ0.01330.0132で µ(zj ) を求めた．M = 309 とするとき，再構成に要した計算時間は，Xeon E5-2650 v4 (2.2GHz) での計算で約 970秒であり，本研究の実装で主として計算時間を要したのは{P−m (xℓ ) ; −M ≤ −m ≤ 0, 1 ≤ ℓ ≤ L} の計算に約 688-730.01310.0130
-400

-350

-300

-250
-200
-150
Minimum Mode -M

-100

-50

0

Fig. 7: L2 norm of the imaginary part of u0

秒，{U (−1) (zj ) ; zj ∈ Dh+ } の計算に約 271 秒であった．また{P−m (xℓ )}, {U (−1) (zj )} は xℓ , zj の点ごとに求める自明な並列化が可能であり，この場合， OpenMP による 24 スレッ
ドでの並列計算では約 157 秒であった．そのうち {P−m (xℓ )},{U (−1) (zj )} の時間はそれぞれ約 62 秒，約 92 秒であり，提5. 打ち切りパラメータ M の選択への注意
提案する数値計算手法では，u−1 および p−m に現れる級数を打ち切るために，現れる Fourier モード (周波数) を制限する整数 M を導入した．全周観測の場合には，この Mは正則化パラメータに相当する (5) ．本節ではその選択について，幾つかの注意を与える．

|| Im U(0) ||2

案手法は並列計算に適することがわかる．

0.07

0.01500.06

0.01450.05

0.01400.04

0.01350.0130

0.03

-359

-309

-340

-320

-300

-280

-260

0.02
0.01
-400

-350

本論文では，u の Fourier 係数のひとつである u0 に注目

-300

-250
-200
-150
Minimum Mode -M

-100

-50

する．この u0 は提案する再構成手法では現れないが，X 線

Fig. 8: L2 norm of the imaginary part of u0 inの強度が正値であることから (1) で定めた u は実数値であ

wider ranges

0

り，したがってu0 (z) =

1
2π

Z 2π


u z, ξ(θ) dθ

Fig. 7 および 8 に M と Im U (0)の値を示す．Fig. 7 よ
20り，M がある程度大きな範囲 (例えば M > 10) で Im U (0)2が最小となるのは M = 73 であった．一方で Fig. 7 の縦軸− 41 −

の範囲の幅は 6 × 10−4 であり，数値計

様の傾向であり，多くの点で 20% 以下となっている．また

算で用いたパラメータに比して極めて小さいことに注意す

Dh+ に含まれる格子点の個数は 24725 個であり，(9) の値がる．実際，縦軸の範囲を変えると Fig. 8 を得るが，この場

100 を越えた格子点数は，M = 73, 309, 357 の場合にそれぞ

合，M = 73 や M = 309 付近での Im U (0) の変化に比してれ 7326 個 (30%)，7956 個 (32%)，11090 個 (45%) であった．

に示す

Im U (0)2が有意に増大していることがわ

このことから M = 357 での再構成は FBP との差が大きい

かる．多くの Fourier モードを利用するほど細部の構造が再ことが定量的にもわかる．さらに M = 309 の場合に (9) の現されることを期待すると，極小値を与える M = 309 およ値が 100 以上となった点を Fig. 11 に示す．これより，(9) のび急激な増大前の M = 357 も，再構成にもちいるパラメー

差が大きくなるのは対象の「外部」(µ = 0) に相当しているタの候補である．

こともわかり，物体「内部」における断層画像は FBP と提M = 359 付近で Im U

(0)2

案手法で充分に一致していることも定量的にわかる．

M=73

600

M=73
M=309
M=357

0.5
0.25
0
-1

-0.5

0

0.5

1

Number of Lattices

500

400

300

200

100

M=357
0

0.5

0

0.25

-0.5

0

0.5

40
60
Relative Difference [%] in Equation (8)

80

100

Fig. 10: Number of lattices taking relative dif-ferences deﬁned in (9). The number of total

0
-1

20

lattices is 24725 in Dh+ , and numbers of lat-1

tices taking relative diﬀerences over 100 are
Fig. 9:Image reconstruction by proposed

7326, 2956, and 11090 for M = 73, M = 309,

method from partial measurement in Fig. 4

and M = 357 respectively.

(top) M = 73 achieving the minimum; (bottom) M = 357 including artifacts
これら M = 73, 357 による再構成結果を Fig. 9 に示す．

0.5

M = 309 とした Fig. 6 と同様，M = 73 では良好な再構成が得られる．一方， Im U (0)0.25の増大の直前の M = 357 で2は画像全体にアーチファクトが現れており，この観測データ

0および分割数においては，より小さな M を選択すべきであることを示唆している．すなわち， Im U

(0)も重要であることがわかる．

の微小な変化

おける再構成と，FBP による再構成の差の相対値を，µ(M ) (zij ) − µ(FBP) (zij ),µ(FBP) (zij )

-0.5

0

0.5

1

2

再構成の差を定量的に調べるため，M = 73, 309, 357 にDiﬀ(M ) (zij ) =-1

zij ∈ Dh+Fig. 11: Region where relative diﬀerences (9)is greater than 100 for M = 309.
(9)

なお Fig. 7 および Fig. 8 で示した Im U (0)

の挙動や，
2

で定義する．ただし µ(M ) (zij ), µ(FBP) (zij ) はそれぞれ打ち

それをもとに論じた M の具体的な値は，一般に，再構成す切りパラメータ M の提案手法と (全周観測の) FBP で得る領域によって異なることに注意する．

られる zij での再構成値とする．前節の再構成例において，

Diﬀ(M ) (zij ) を横軸に，その値をとった格子点数を縦軸にとったものを Fig. 10 に示す．この結果より， FBP の結果を基準としたときの提案手法の誤差は，M = 73, 309 の場合で同− 42 −

6. 結言(6) H. Fujiwara, K. Sadiq, and A. Tamasan : Numerical Re-

本論文では，X 線計算機断層撮影において，A-analytic 理

construction of Radiative Sources in an Absorbing and論にもとづく境界積分の手法により，部分観測から再構成画

Non-diﬀusing Scattering Medium in Two Dimensions,

像を取得する数値計算アルゴリズムを提案した．また実測値

SIAM J. Imagin Sci., 13 (2020), pp.535–555.

からの部分的な再構成をおこない，その妥当性を示した．提(7) H. Fujiwara, K. Sadiq, and A. Tamasan : Partial Inver-

案手法には正則化パラメータの役割を果たすパラメータ M

sion of the 2D Attenuated X-ray Transform with Dataも現れるが，結果の精度と安定性を両立させるための選択に

on an Arc, Inverse Probl. Imaging, (2021), DOI :ついての基準を論じた．提案する再構成手法は各点ごとの再10.3934/ipi.2021047

構成を与えることから，本論文で論じた「領域 D+ におい(8) H. Fujiwara, K. Sadiq, and A. Tamasan : Numericalて単一の最適な M 」という条件での選択基準は目安に留め，

Source Reconstruction from Partial Measurement in an

実用化においては，例えば ∂D+の近傍では Cauchy 核の数値積分による精度の劣化も考慮し，部分領域ごとに局所的にM を選択するなど，多様な検討が必要と考えられる．

Absorbing and Planar Scattering Medium, (in preparation).
(9) 藤 原 宏 志 ，A. Tamasan : Cauchy 型 積 分 に よ る メッ

謝辞 本研究の一部は JSPS 科研費 JP20H01821, 21K07593,
Austrian Science Fund (FWF) P31053-N32, お よ び NSF
grant DMS-1907097 の助成を受けたものです．

シュレス X 線計算機断層撮影法, 計算数理工学論文集,19(2019) pp. 1–6．
(10) A. C. Kak and M. Slaney : Principles of Computerized
Tomographic Imaging, (1988), IEEE Press, New York.
(11) W. Koppelman and J. D. Pincus : Spectral Representa-

参考文献

tion for Finite Hilbert Transforms, Math. Z., 71 (1959),(1)

T. A. Bubba, et al.:Tomographic X-ray Dataof a Lotus Root Filled with Attenuation Objects,
arXiv:1609.07299v2 [physics.data-an] (2016)(12) 工藤博幸: インテリア CT における画像再構成, Med.
Imag. Tech., 34 (2016), pp. 186–197.

(2) A. L. Bukhgeim : Inversion Formulas in Inverse Problems,

pp. 399-407.

(13) J. Radon : Über die Bestimmung von Funktionen durch

Linear Operators and Ill-Posed Problems byihre Integralwerte längs gewisser Mannigfaltigkeiten,

M. M. Lavrentiev and L. Ya. Savalev, (1995), Plenum,

Berichte Sächsische Akademie der Wissenschaften zu

New York, pp. 323–378.

Leipzig, Math.-Phys. Kl., 69 (1917), pp. 262–277.

(3) D. V. Finch : The Attenuated X-ray Transform: Recent(English translation : On the Determination of Func-

Developments, in Inside Out: Inverse Problems andtions from Their Integral Values Along Certain Man-

Applications, Math. Sci. Res. Inst. Publ., 47 (2003),

ifolds, in IEEE Trans. Med. Imaging, MI-5(1986),

Cambridge Univ. Press, Cambridge, pp. 47–66.

pp. 170–176.)(4) Finnish Inverse Problems Society : Open X-ray Tomo-(14) K. Sadiq and A. Tamasan : On the Range of the Atten-

graphic Datasets, https://www.fips.fi/dataset.php

uated Radon Transform in Strictly Convex Sets, Trans.

または DOI : 10.5281/zenodo.1254204

Amer. Math. Soc., 367 (2015), pp. 5375–5398.

(5) 藤原宏志，大石直也 : 境界積分による X 線計算機断
層 撮 影 法 の 正 則 化 ， 計 算 数 理 工 学 論 文 集,20(2020)

pp. 75–79．

− 43 −(15) 塩見剛 : トモシンセシスの原理と応用, 医用画像情報
学会雑誌, 24 (2007), pp. 22–27.