計測自動制御学会東北支部 第 228 回研究集会 (2006.5.17)
資料番号 228-8

長方形スキャンによる画像再構成方式の提案
A Proposal of Image Reconstructing Method
Based on Rectangular Scanning○齋藤昭宏 ∗ ，田山典男 ∗ ，光安康人 ∗ ，加藤弘典 ∗∗ ，多田明広 ∗ , 栗田宏明 ∗ ，佐藤宏明 ∗○ Akihiro Saito∗ , Norio Tayama∗ , Yasuto Mitsuyasu∗ ,Hironori Kato∗∗ , Akihiro Tada∗ , Hiroaki Kurita∗ , Hiroaki Sato∗*岩手大学工学部, **富士通ゼネラル*Faculty of Engineering, Iwate University, **FUJITSU GENERAL LIMITED
キーワード : CT (Computer Tomography)，画像再構成 (Image reconstruction)，ウェーブレット (Wavelet)，
一様透過性 (Nature to transmit uniformly)，長方形スキャン (Rectangular scanning)

連絡先 : 〒 020-8551 岩手県盛岡市上田４丁目３番５号 岩手大学 工学部 電気電子工学科 田山研究室
田山典男，Tel. & Fax.: (019)621-6382 E-mail: tayama@pipe.elc.iwate-u.ac.jp

1.

はじめにによるファンビーム状の投影データから平行ビーム状の投影データへ変換する方法を示し，正しく

医療分野1−3) で大いなる貢献を果たしているX線変換され画像再構成されることを検証する．

CT 装置では専ら FBP 原理が採用されており，400〜1000の多数方向からの投影データを用いている．
これは高精細な画像を再構成行うためであり，装2.

置が大がかりで高価になっている．

FMR 法

まず画像再構成の基礎となる FMR(Fast Model

一方産業分野では，工場や空港でのベルトコンReconstruct) 法について説明する．

ベアなどの流れ工程において，対象物体を非破壊画像を対象とする 2次元対象空間において，遮断的検査するために 3 次元物体として可視化したい

空間角周波数ωm よりも高い周波数成分を制限したという要求が高まっている．そのためには，導入標本化定理を満たす濃度値分布モデルを導入する．

するコストが低く，高速に処理できるコンパクト

このような条件を満たす対象領域の連続空間では，な産業用 CT 装置の開発が求められている．

任意点の濃度値を f (x, y) とすると，t = π/ωm 間隔

本研究では，対象空間に標本化モデルを導入し，以下で標本化された各標本点の濃度値 f (xi , yj ) と，ウェーブレット標本化関数を用いることで局所領標本化関数 S(x) を用いて記述することができる．

域による画像再構成法を検討している4−9) ．X線源

対象空間領域を(n×n)の有限領域とすると，f (x, y)が長方形の周辺上を移動して投影データを採取する長方形スキャン方式を提案する．長方形スキャン–1–

は標本化定理から次のように近似できる．

f (x, y) =n−1
∑ n−1∑f (xi , yj )·S(x−xi )·S(y −yj ) (1)i=0 j=0

ここで透過ビーム上の各点での濃度値 f (xi , yj ) を透過ビームに沿って線積分した値は，投影値 p に等

Fig. 1

しくなり，式 (2) が成り立つ．

Sinc 関数 S(x)∫ Lpf (x, y) dl∫0 L n−1∑ n−1∑=
f (xi , yj ) · S(x − xi ) · S(y − yj ) dl=

0

i=0 j=0(2)ここで透過ビームが対象空間へ侵入する点から脱出する点までの長さを L として，侵入点から透過

Fig. 2Sinc 関数のフーリエ変換 Ŝ(ω)

ビーム上のある点までの長さを l とする．
透過ビームごとに式 (2) が成り立つので，これめる濃度値ベクトル f は式 (6) になる．

らの連立線形方程式を行列の形で表現することに

C + = V Λ+ U t(5)p = (pm ) は，影響の度合いを表す線積分影響係数

f = C+ p(6)行列C = (Cmn )と，標本点濃度値ベクトルf = (fn )このとき最小二乗解が複数あるので，濃度値ノルにより式 (3) のように表される．e = (em ) は誤差ベ

ム kf k が最小となるものを選ぶ．f は最小二乗・最クトルである．

小ノルム解と呼ばれる．式 (5) の Λ+ は，Λ の転置する．このように定式化すると，投影値ベクトル

p=Cf +e(3)行列の対角要素を逆数に置き換え，さらにランク

この式 (3) の影響係数行列 C は M × N 行列であ落ちの起こった対角要素を 0 にしたものである．

り，m 行 n 列の要素が線積分影響係数 Cmn である．

これら C + は定数行列であり，あらかじめ計算C は透過ビームの形状や投影方向数，投影方法なして蓄えておくことができる．このようにして得どから定まるものである．この影響係数行列 C を，られた定数行列C + と計測した投影値pとを掛け合

M × M 直交行列 U と M × N 対角行列 Λ と，N × Nわせるだけで，最小二乗・最小ノルム解である濃直交行列V の転置行列V t とに分解することを特異度値ベクトル f を求めることができる．

値分解という．

C = U ΛV t(4)3.

式 (4) の行列 Λ の対角要素は非負値であり，降順に並んでいる．この各要素を特異値という．

ウェーブレット標本化関数

これまで標本化画像再構成においては，標本化関数として Shannon の標本化関数 (Sinc 関数) を用

2

式 (3) において，再構成の誤差二乗和 (kek ) を最いてきた．しかし Fig.1 に示すように影響領域が局小にする解 f を求める再構成問題は，種々の一般在化しておらず，有限領域で打ち切って使用する逆行列のうちで Moore-Penrose 型行列 10) を求めると誤差が大きくなるという問題があった．

問題として定式化できる．従って，行列 C がラン

これに対して，Meyer スケーリング関数 φ(x)11)

ク落ちしているときには式 (5) として表される．求では振動の幅が狭く局在性がよいことが分かって–2–

Fig. 5
Fig. 3Fig. 6Fig.

4

ウェーブレット関数 W (x)Meyer スケーリング関数 φ(x)ウェーブレット関数のフーリエ変換 Ŵ (ω)

Meyer ス ケ ー リ ン グ 関 数 の フ ー リ エ 変スケーリング関数のフーリエ変換 φ̂(ω) は，Fig.4

換 φ̂(ω)に示すように，角がなく滑らかになり，局在性がいる．そこでこの Meyer スケーリング関数に対しよい．フーリエ逆変換を施したスケーリング関数て，標本化関数の性質を持たせたウェーブレットφ(x) の形状は，Fig.3 のようにすぐ減衰する局在性標本化関数を定式化していくことにする．

がよいものになる．

まず Meyerウェーブレットでは，スケーリング関

このウェーブレットのスケーリング関数をもと数のフーリエ変換 φ̂(ω) が，次の性質を持つと定義に標本化定理を満たす標本化関数を求める方法にされている．

ついてはすでに知られている．多重解像度解析空間において，有限角周波数の任意信号 g(x) は，ス1) φ̂(ω) は無限回微分可能である．

ケーリング関数 φ(x) を用いて一意に表現できる．

2) φ̂(ω)

1,|ω| ≤ 2π
33=
cos[ π2 ν( 2π|ω| − 1)], 2π≤|ω|≤ 4π33
4π0,|ω| ≥ 3

スケーリング関数 φ(x) のフーリエ変換を φ̂(ω) として，離散フーリエ変換を φ̂∗ (ω)とすると，標本化定理を満たすウェーブレット標本化関数のフーリエ(7)変換 Ŵ (ω) が次式で求められる．

Ŵ (ω) = φ̂(ω)/φ̂∗ (ω)

ただし，ν は次式を満たす平滑関数である．
{0, x ≤ 0ν(x) =, ν(x) + ν(1 − x) = 1
1, x ≥ 1(9)そこで，式 (9) よりフーリエ逆変換を施して，ウェーブレット標本化関数W (x)を次式で算出できる．

またウェーブレットのフーリエ変換 ψ̂(ω) は，式1
W (x) =2π(8) で表現される．
ψ̂(ω)
34πiω/2|ω| − 1)], 2πsin[ π2 ν( 2π e3 ≤ |ω| ≤ 338π=eiω/2 cos[ π2 ν( 4π|ω| − 1)], 4π3 ≤ |ω| ≤ 30,|ω| ≥ 8π3∫ ∞φ̂(ω)∗−∞ φ̂ (ω)· ejωx dω(10)ウェーブレット標本化関数のフーリエ変換 Ŵ (ω)
をFig.6に，ウェーブレット標本化関数W (x)をFig.5に示した． このウェーブレット標本化関数 W (x)の(8)形状は，Fig.1 の Sinc 関数 S(x) に比べて，局在性が–3–改善されている様子が分かる．標本化定理の条件は，角周波数の帯域幅については −π < ω/T < +πとなり，近似的に満たすことが分かる．今回の実験ではこのウェーブレット標本化関数 W (x) を採用している．
このように局在性があるウェーブレット標本化関数 W (x) を用いることで，FMR 法と組み合わせると高速化が期待できる．それは従来までの FMR

Fig. 7法では Sinc 関数 S(x) を採用していたために，ひと一様透過性

つの点を再構成する際にも領域全体からのデータが必要であったが，局在性のあるウェーブレット標本化関数 W (x) を用いることで局所領域で再構成が可能になる．つまり再構成に必要な部分だけを読み出すだけでよくなり，計算量が大幅に削減可能となった．

4.

長方形スキャン方式
Fig. 8

長方形スキャンVLSI による高速な画像再構成を行うために，以下ではあらかじめ蓄えておく一般逆行列 C + の係

4.2

長方形による投影角度の選定数値の個数を大幅削減する方法を提案する．
一様透過性を満たすことで，センサーに対して
そこで本稿では，最初にX線を平行状ビームとし水平方向では再構成係数並びを大幅削減することて，ビームの入射パターンに規則性を見いだす．次ができた．さらに，垂直方向にある規則性を持たにX線を物体に近づけた状態で投影計測したファンせることにより，垂直方向でも再構成係数並びを状ビームから平行状ビームへ変換する方法を示す．
大幅削減する方法を示す．

4.1

Fig.8 のように，X 線源が長方形の周辺を移動し一様透過性ながら，相対応する辺のセンサで投影値を計測す

標本空間においてセンサーを固定配置する．する方法を「長方形スキャン方式」と呼ぶ．ここで，るとセンサーに対して水平方向の標本点の各層に平行ビームの通り方が何層おきかに同じようになおいては，すべて同じ入射パターンを示すことがるように投影方向の選定を工夫する．

分かる．これを「一様透過性」という．Fig.7 にそつまり，Fig.9 のように長方形の中心画素とそのの様子を示す．この一様透過性を常に満たすよう周辺画素とを結ぶビームの方向を投影方向としてに，標本空間にビームを入射させることで，同一選定する．この例では 7 × 7 の正方形であるが，セの層に対して 1 パターンの再構成係数並びを用意ンサーに対して垂直方向に 3 層おきに平行ビーム+しておくだけでよくなり，一般逆行列 C の再構が入射することが分かる．このように投影方向を成係数並びを垂直方向の層数分のみ蓄えておくだ

選定すると，垂直方向に対しては 3 パターンの再けでよいことになる．

–4–

Fig. 10

Fig. 9

原画像 (256 × 256)

投影方向の取り方構成係数並びを蓄えておくだけでよい．大きな再構成領域でも，再構成係数並びを繰り返し使用することができるため，情報量を大幅に削減することができ，VLSI の内部に再構成係数並びを格納すFig. 11ることが可能になる．例えば，1000 × 1000 の画像

再構成画像 (256 × 256)では，一様透過性により水平層の 1000 画素が 1 画案する．ファンビーム照射で計測された投影デー素分の再構成係数並びでよくなるので 1000 分の 1タから平行ビーム投影で得られる投影データへ変に削減できる．さらに投影角度の選定により，垂換することができるならば，これまでの平行ビー直方向の 1000 画素が 3 画素分でよくなるので 1000

ム投影での種々の利点を活かすことができる．

分の 3 に削減できる．

ここで平行ビームとは，センサーを基準としたこのような長方形スキャン方式を採用すること水平方向の標本点層に対して，空間的に平行に入で，装置全体の大きさをコンパクトにすることが射するビームのことである．本研究では，ファン期待でき，一様透過性や垂直方向に現れる規則性ビーム照射をセンサー上の標本点間隔で時間的にを活用できるので，専用 VLSI のハードウエアの設移動させるならば，その計測データの中から平行計が可能になり，高速で安価な画像再構成の可能ビーム投影データを抽出生成することができるこ性が期待でき，実用化に向けて大きな前進となる．とに留意する．
本研究ではこの性質を利用し，長方形スキャン

5.

投影データの変換によるファンビーム投影データから平行ビーム投

これまで X 線を対象物体から遠く離れた状態で投影計測する平行ビーム投影を想定してきた．装影データへ変換することで，これまでの一般逆行列 C + に基づく画像再構成を行っている．

置を小型化する観点からは，X 線源を対象物体に近づけた状態で投影計測することが望ましい．そ6.

実験

こで本研究では，X 線源が長方形の周辺上を移動

Fig.10 に示す画像サイズが 256 × 256 の画像を原して，対象物体の近くから扇形ビームを照射する，画像として，シミュレーション実験を行った． 平行

ファンビーム投影による長方形スキャン方式を提ビームに変換した際の投影方向数を 40 方向，局所–5–領域サイズを 89 とした場合の再構成結果を Fig.11

参考文献に示す．このように比較的良好な画像が得られた．

1) 尾上 守夫：医用画像処理, 朝倉書店 (1982).

両者の正規化誤差分散をとってみると 0.6370 であ2) 江尻 正員: 画像処理産業応用総覧, フジ・テクノシステム, (1995).

る．再構成の計算時間は，通常のワークステーションを用いて 3.97 秒であり，大変高速であることが分かった．

7.

3) H. Hiriyannaiah: X-ray Computed Tomography
for Medical Imaging, IEEE Signal Processing,Vol.114, No.2 45/59 (1997).
4) 田山, 船岡：緩やかな 2 次元濃度分布に対する特異値分解による可視化計測, 可視化情報学会誌, 14 巻
52 号 21/27 (1994).

むすび

本研究では，産業分野向けの X 線 CT 装置の開発を目指して，画像再構成 FMR 法にウェーブレット標本化関数 W (x) を導入して，さらに長方形スキャン方式を採り入れた，新しい画像再構成手法を提案した．
本手法では，長方形の周辺を移動するX線源からのファンビーム投影データを計測して，平行ビーム投影データへ変換し，各方向の投影ビーム間での規則的な透過位置関係を活用する工夫により，再構成係数並びのデータ量を大幅に削減している．
本手法の主な特色を，以下に示す．

1) 長方形スキャン方式に基づきファンビームで投影計測するので，装置を小型にできること．

2) 再構成の計算構造が規則的であり，再構成係数並びのデータ量が大幅削減できるので，VLSI にアルゴリズムを搭載して能率よく実行できる可能性が高いこと．
今後の課題は，具体的検討を進めてX線スキャナ装置を構成し，再構成を確認評価することである．
さらに，長方形スキャンの電子化を検討し，超並列アーキテクチャの VLSI 設計を行うことである．
本手法は長方形電子スキャンと専用 VLSI を導入することにより，ベルトコンベア上でのような高速性が求められる分野でも手ごろな価格で対応できる可能性があり，産業分野の不可視な対象物体の内部をリアルタイムで画像化する簡易装置の構築可能性が期待される．

–6–5) 田山, 加藤, 大坊, 長谷川, 杜, 栗田：傾斜スキャンによる X 線投影データ取り込みと 3 次元画像再構成の実験, 画像電子学会誌, 28 巻 2 号, 152/160 (1999).
6) N. Tayama, H. Du and M. Daibo, H. Kurita, K.
Seki: Development of Computed 3D Imaging System from a Few Projections, Proc, IASTED SIP’99, 119/123 (1999).
7) H. Du, N. Tayama, M. Daibo, T. Hasegawa
and K. Seki: A Computed Imaging System Using Wavelets Sampling Model, Proc. MVA ’98,115/118 (1998).
8) 杜, 田山, 渡邊, 関: VLSI 向きの極少数方向投影によるウェーブレット部分画像再構成, 画像電子学会
誌 30 巻 3 号, 233/241 (2001).
9) 田山, 杜, 大坊: 産業用実時間 3 次元 CT の画像再構成プロセッサ, 計測と制御, 40 巻 12 号, 907/910(2001).
10) A. Albert: Regression and the Moore-Penrose
Pseudo Inverse, Academic Press, New York (1972).
11) Y. Meyer: “Wavelets : Algorithms and Applications”, SIAM, Philadelphia (1993).