Title
反復解法による3次元画像再構成の高精度化
Author(s)
門馬, 寛明; 陳, 延偉; 仲尾, 善勝
Citation
琉球大学工学部紀要(56): 43-49
Issue Date
1998-09
URL
http://hdl.handle.net/20.500.12000/1967
Rights
琉球大学工学部紀要第56号11998年 4３
反復解法による３次元画像再構成の高精度化
門馬寛明． 陳延偉↑ 仲尾善勝↑ ＡｎＡ(:(:111.atcA1g()l･itl〕ｍＲ)1.1tGI･ativeThree-DimensionalReconstl･uction HiI･oal(ｉＭ()ｎｍａ，Ｙ・-ＷＣｈｅｎ，andZenshoNaka(） Ａｂｓｔｒａｃｔ ＡＲＴ法のような反復解法による画像再梢成では,反復毎に投影を計算する必要があるた め,投影を計算機で高精度かつ高速に計算しなければならない.本論文では,高精度かつ高速 な再投影計算法を提案し，従来法との比較を行う． 要があるため,投影を計算機で高精度かつ高速に計算し なければならない．本論文では,高精度かつ高速な再投 影計算法を提案する． １はじめに 現在,快)||Ｘ線ＣＴ画像再ji稚成にはフィルター補正逆 投影法が使われている．異なる原理に基づく再栂成法， 例えば,ＡＲＴ(AIgcbrai(:RBconstrucbion亜chniques） 椎も知られている．この方法は,反復再ｷｨi成怯の一つとし て１９３７年にＫａ(:zm勘rz[1]によって提案され,Tanabe がその収束性について論じた[２１初期の頃は射影法 (MeCllo(IofProjecCi〔)､)とl平ばれていた．ＡＲＴ法によ ゐｌｖ榊成li1li像の画質はフィルター補正逆投影法による ものより劣るので，あまり注目されなかった．しかし， フィルター補正逆投影法を投影データが部分的に欠損 している場合の再構成やコーンビーム投影データから の３次j6mii像再榊成などに適用す為のは困難である．こ れらの場合には，ＡＲＴ法が有効であるといわれている [31．それ以後Kaczmarzによって提案されたＡＲＴ法 は再び脚光をあびるようになってきた[4,51.この方法 は医用画像だけではなく，地理断層の解析やレーザー核 融合で応用されており[6,71,応用技術の確立が望まれ る．１９７０年代になってGordon,HerｍａｎなどによりＡ ＲＴの近似解法が提案された[８１ ＡＲＴ法のアルゴリズムは,(1)加法的ＡＲＴ(Addi‐ tiveART),(2)乗法的ＡＲＴ(MulbiplicaCiveART)の ２つに大別される．多くのアルゴリズムが提案され,そ の性能比較もされてきた．一方,ＡＲＴ法のような反復 解法による画像再構成では,反復毎に投影を計算する必 ２反復解法による画像再構成 反復解法による画像再構成では，まず適当な初期推定 画像を設定し，次のような処理を繰り返す．Stepｌ推定画像に投影を行い,投影データＲを算出する
Step２原画像からの投影データＰと推定画像からの投
影データＲを比較 Step３投影データＰとＲが近似されるように推定画像 に補正を加える このように反復毎に投影を計算する必要があるため， 高精度かつ高速な投影計算を行わなければならない．さらに,離散画像における投影アルゴリズムを大別す
ると次のようになる．Stcpl離散画像を回転させた回転画像を計算する
Step２回転画像の直交座標上において零度方向から投
影を行う この投影アルゴリズムから,高精度な投影データを得るには,高糀度な回転画像の計算を必要とすることがわ
かる．従来の方法では,座標回転時において離散的に座 標を配置する際,新しい座標軸に対して各格子間の画素値は,近傍の格子点からの線形補間によって決定される．
受理：1998年５月２５日 ･大学院理工学lijf究科電気電子工学専攻 （CriLdu塾tcStwlent,Depも.ofElectricalan〔lBIectronicEng.） ↑工学部勉気勉子工学科 （[”1兆`)｢Electrical孔UdEIcctronicEng.,【泡c､、『Ellg.）門馬・陳・仲尾：反復解法による３次元画像再綱成の高精度化 4４ 3.3線形重み補間法 芯線形重み補間法 3.1線形補間法 線形重みｷIlillM法は，格子点におけあ重みの綿|がｌと なるように線形ﾄﾞﾙ間法に重み補正を加えたﾉﾉ法であり， 高精度な回転ijHi像を得る 図２は，格f･点（w,妙）に樹:\する重みの分hjlxlであ る．笑線はlijl転後のＵＶ畷標糸であり，点線はＩｉＩｌｉＩ唖ＩｊｉｉＵ） ＸＹ座標系を表している．
｛
エcos8＋y6il'０ －.BHirl0＋WcosO 皿Ｕ (1) 上式は,座標(ｗ）を角度０逆回転した回転座標系 における座標を表しており，皿,Ｕは実数値をとるこの とき,座標(１J,Ｕ)は回転座標上の格子点（整数値）と一 致しない場合が多い． 図２:格子点（'`,'ノル対する砿みの分布 例として,格子点(秘,T')の画素値を求めJXLまず線形 補間法を用い】画素値/(皿,U)を求めると ノ(皿，ひ）＝池ｗ,).８１＋/(鯉2,z/2）８２ ＋/(⑪3'"3).SIB＋（zｲ'"4)･SII（２）次に式(3)のように,格子点(u,U)に対する重みの合
計Ｑ(u,U)を計算する． Ｑ(皿，U)＝β'＋３２＋ｓ;＋ｓ４（３）Ｑ(u，〃)を用い,式(2)で得られた画素値/(?いりの補
正を行う． 図１：２次元画像の線形補間 線形補間法では，図１で示されるような面積 81,82,s３，８４を配分割合（以下，重みとよぶ）として，座標(u,Ｕ)がもつ画素値ノ(皿,U）を近傍の格子点
Al'A２，A､3,ｋJに配分する．例えば,Ｓ１＝０．２１の場合,画
素値ノ(u,ひ)の21％が格子点ｈ,の画素値ノ(k,）に加
算される． 3.2線形補間法の問題点線形補間法の問題点として,回転画像のある格子点に
対する重みの合計がｌとならないことがある．この問題
点に対する解決策として，ファインサンプリング法の適
用が提案された[9）
この方法は，まず各座標間を細分化し,細分化された
座標に対して線形補間を行うものである．これにより，
線形補間のみを用いた大まかな画素値の配分に比べ，よ
り高精度な画素値配分が行われる．しかし,この方法では各座標間をより細分化するほど
計算時間を必要とする計算時間は，２次元においては
ファインサンプリング数の２乗に比例し，３次元におい
てはファインサンプリング数の３乗に比例する．
．この従来法を用いて高精度な投影データを得るため
には，ファインサンプリング数を増加させる必要があり，
その結果として多大な計算時間を要する[9｝そこで本
研究では,高精度かつ高速で投影データを計算できる線
形重み補間法を提案する(1０１．
小咋佇:::'二④
このように式(4)から得られる画素値／(…)が,線
形重み補間法により計算された格子点(⑪,｡)の画素値
である．３次元離散画像に対する線形重み補間法は’２次元の
場合と同様な考え方によって行われる．３次元では重み
として体積が用いられ，回転画像の格子点に寄与する亜
みが8個となる[1０１．
４投影の離散的表現 ４１２次元まず２次元画像を対象とする投影の離散的表現を示
す.平面上に直交座標系（)-.Ｗを設け,その':の画像を
琉球大学工学部紀要第56号，1998年 4５
/(｡:’1ﾉ）で表す．左上隅を原点として適当な大きさの疋
万領域を考え，それをさらに小正方形（これを画素とい う）に分割し，ｊ＝１，…,ハノと番号付ける．第ｉ番目の llii素におけるiiW像ノ(｡:,I/)の代表値を画素値といい,ノｉで衣す．（ﾉｶﾞ)’'１＝1,…,八'を離散画像という．本欄で
は，あれllUi素ﾄﾆおけも画像の代表値として，画素の中心 にｵjけるl11Ii像の値をとる． ０，<b方向の投影における投影データＰＣ,Ｃｌは次式を計 算することにより得られる． (9) Ｐ０，゛＝、ハノＦＤ,｡ ＦＩＬのはＦの回転画像を表し，Ｗは０．方向からの３ 次元離散IiHi像に対する投影の重みベクトルである．０° ﾉﾌﾞ向からの２次元離散画像ｆｊに対する投影の重みベク トルをｗｊとすると，重みベクトルＷは次式のように 表される．Ｗ＝（ｗ,，ｗ２，…，ｗＬ)Ｔ
（10） このときⅣ＝Ｌ２とすると，ｗｊはＬｘＮ行列,Ｗは /Ｖｘ/Ｖ行列である． 図３：離散画像 4.3回転画像の計算 ここで,線形重み補間法を用いた回転画像の計算を示 す．まず行列Ａが次のように定義されるとき ''１１i素ｸﾘをベクトルで表し，『＝(ん/E'…'/〕ｖ)7．
（Ｔは転置を表す） とおく． このとき次式によりノ(211")を八に対応づける． ｊ＝ＬｘＵ＋。；＋ｌ ただしＬ：画像の幅(ﾉＶ＝Ｌ２） (5)ﾊｲ焉篭:::i費iT､jiWl川
られろ． 回転後の座標(､ｗ’7．)は次式を計算することにより得 (6)IiH三ﾙﾛ
このような定義の下で,例えばini像３×３[pixels1の /'１度０にｶﾅする投影データは，式(7)のように表される． Ｐ(】＝ｗｆＤ (12）NＩＩ
ただし’
１００ １００ １００ ０１０ ０１０ ０１０ ００１ ００Ｉ ｚｃ：画像の中心点の⑰座標"噂：画像の中心点のｙ座標
ｚｃ：画像の中心点のｚ座標 有効範囲１≦u≦L’’三U≦Ｌ’１≦7,≦Ｌ (7) このとき１，＝３，Ｎ＝９ ここで，ｆｌ〕はDITi像ｆを,〃向に逆回転させたもので， ｆと'１１じ／ＶｘｌのiiIIi素ベクトルである．１，０はＬｘＩの 投影ベクトルであり，ｗは敏みベクトルからなゐＬ･×Ⅳ 行列である．次にｕ,U,ｒの小数部を切り捨てた整数[?↓ＭＵ１ｔ１を
求め,式(13)を計算する．'三二鍋|菫冒|ミル)
ここで,回転画像の画素値を/解として次式を計算
する．(に踊柳,＃`,土|("）
ノー1,2,…,Ｅ，Ａ-1,2,…,Ｅ（Ｅ＝Ｌ３）
４．２３次元 ２次元噸像は式(5)のように表されるとし，ＬｘＬ×Ｌ[Ⅵ)xq11H1V)３次元imi像を考える．２次元断層画像を
f:j,ｊ＝1,2,…,1`とすｵしば３次元画像Ｆは次式のよう
に表される．Ｆ＝（fi，Ｌ，…，ｆＬ)】．
(8)門馬・陳・仲尾：反復解法による３次元画像再構成の高精度化 4６ ５シミュレーション
錦,に対する（の重みを要素にもつ行列ｓを用意
し,すべての要素は初期値０を持っている.式(13),(14） で得られた値を用いて,次式を計算する． （１）ファントムファントムとして有名なSI1epljnndLol1Ianファントム を用いて，２次元画像におけるシミニルーションを行った． ファントムの画像サイズは,64ｘ64{piXCl蘭1であるま た投影方向数を１０方向とし,投影角８Ｍ。＜０＜180｡） を18.の等間隔にとり，投影を行っている．ただし、投 影データは再投影と同じ方法で得られてい’(jため、再投 影による誤差は考慮していない。再生誤差はアルゴリズ ム自身の誤差となる。 ３次元1mi像におけるシミュレーションでは（】八Wｉｓフア ントムを)l1いた．ファントムのimi像サイズは,３２ｘ３２ｘ;l2 1vnx畷15]であるまた画像の中心点(３ぃ,加毒,）におけ為 ｌＵ１ｉ紫仙（I[蜆分布関数のピーク値）は255.(),lmi紫の最 小値が０である．投影方向数を４方向として,投影角は 適当な値を与えており，投影ベクトルＰには,理論値を 与えている．本研究における理論値とは,１１．方向からの 投影データである． (2)誤差 原画像Ｆと推定画像Ｆの誤差評価式を示す． ｛了】７】ｒ一丁７７γｒ △△△△△△△△ Ｑの●。’〃。？■己■■□Ｐ■■〃■か０□．■Ｂ Ｇ■■ヴ▲』■■ヴ●も■■ＰＩＦＪ。Ｌｑ０ログ、曰０■▽■■Ｄｐ△ｑⅡＪＦ ＩＣリ》｜⑰几〉β、》５－卯一⑪｜旬Ｉ』】⑥ ロロェ、コル、ＯＬ氏ｎ匹『〃血（○ｍ一ｍりゆへ０小Ｌ ｌＪ■１日ｎＦ６已已ｊｒ■可〃Ⅱ■、凸■■二■。夕■日日▽〃■■も少一■■Ｕ ’（』『】〈配『、〉〈』『｜⑫）知』「四】銅》｜「四）（回『叩｝知。』加）ハ一一ｍ〉 △△△△△△△△ －Ⅷ弘一皿皿一⑭皿一弘皿 ××ｘ×ｘ×ｘｘ △⑪× △Ｕ× △Ｕ× △7）ｘ △⑪ｘ △７１× △Ｔｌ× △Tjx Ａｌ＝ん 182＝ん＋l As8l＝ん＋Ｌ Ｍ＝AP＋Ｌ＋１ k月＝化1＋Ⅳ A(i＝Ａ２＋/V AP７畠ごA艤;〕＋/Ｖ Ａｌ８:-:１Ｍ＋Ⅳ ただしSKj：Ｓのｋ行ｉダリ|]の妾莱 以上の計算をすべての画素に行い,行列Ｑを求めるイ±■!’
(15） ＥＱｋ＝ESX:ｉ
ｆ＝１ Ｑ上：第１６番目の画素ﾄﾆ対する重みの合計最後に次式を計算すると，回転画像ＦＤ,‘を得る．
Ｅ７７℃「 （19） また,原画像Ｆの投影データＰと推定画像Ｆの投 影データＲの誤差評価式を示す． Ｆ9,小＝ＱＳＦ Ｕ６） また次式を計算すると，従来の線形補間法による回転 画像が得られる． (20） Ｅ７．T･ｏ７･＝ これらは３次元画像における評劒価式であり，２次元画 像における評価式はz=ルー０とすることにより得 られる． (3)画像再構成法反復解法の一種であるＭＡＲＴ法[11]とＥＭ－ＭＬ法
[12]を用いて,画像再構成を行った．以下にＭＡＲＴ法 とＥＭ－ＭＬ法を示す． Ｆ０，ｄ＝ＳＦ (17） 4.4投影データの計算 式(16)を式(9)に代入すると，次式が成立する． ＰＣ,ウーｗＱｓＦ ＝ＨＦ （18） このとき Ｈ＝ＷＱＳ i)ＭＡＲＴ法 lＴｊｄ鮒Ⅶ(篝）
(21）ここでｗは一定解をもち,Ｑ,ｓは回転角０，のに対
して一様に求まるので,Ｈは一定解をもつ．したがって，
反復毎にＨを計算する必要はなく，回転画像を計算せ
ずに投影データを得る， ‘＝１１２，…,Ⅳ，ノー１，２，…,A‘ ここで,画素数をⅣ,投影数Ｍ,反復回数Aとした． Ｐは原画像Ｆから得られる投影ベクトル!Ｒは推定画琉球大学工学部紀要第56号，1998年 4７ 表ｌ：誤差と実行時間の比較 図４８再構成結果 15.4 5０ 4５ 階』 4０ ２ ０ ● ５ 時 ｌ ｈＣ自由目 』●』白回 ３５ 。0 2５ 2０ Ｍ』 “､ロ６０００ 額00。OOO Iterntl⑨ⅡⅡ ０１０曲 4０００６０００ 20,0３０DD Iteratlon ＤＩＯＯＯ 図６：ＥＭ－ＭＬ法による薗潔値誤差 図５:ＭＡＲＴ法による薗素値誤差 ＭＡＲＴ ＥＭ－ＭＬ
mli素値の誤差
17.645 ⑬871） 1４．９２３ (75）実行時間[秒］
657 253Pｈ制ntUm MART法 EM-ML法
SIUeppmUdLog5m E｢】⑪r＝17.65[9'６］ E】･ﾛ･ｏｒ＝14.,2[恥】
－－－■０
Ｐ００■■０Ⅲ０■０，０，■■■■■
Ｆﾖ祠
～－－－－－－－－－￣－￣
門馬・陳・仲尾：反復解法による３次元画像再構成の高精度化 4８ 図７；Ｇ;Lussフアントム 図８:Gaussフアントムの再橘成結果 図９:Gaussフアントムの再構成結果 1０ 1２
=澤轆鱗闘瀧
８ ５４Ⅱ０ 白●角自国 角。白白回 ６ ４ ２ 0 ０ Ⅱ｡ 2０ ｡｡ ＄０ 6０ 0 1０ 2０３０ 4０５０ IteratioIn Iterntion 図１０:画素値の誤差 図１１：投影の誤差 表２：誤差と実行時間の比較 線形重み補間法 Ｚ＝１５ Ｚ＝1３ Ｚ＝１１ 線形補間法 １１回目 ５０回目 線形重み補間法 １４回目 ５０回目 画素値の誤差 4.243 5.410 1.723 2.002 投影誤差 4.263 4.098 0.769 0.636 実行時間 10[秒１ 11[秒]琉球大学工学部紀要第56号，1998年 4９
