1993年[ 技術開発研究助成 ] 成果報告 ： 年報第07号
符号化開口CTを用いた生体組織内RI分布の3次元計測
概要
1.はじめに
符号化開口CTはそれ自身が放射を出す放射源の3次元分布を非接触・非破壊で計測する放射型CTの一手法である')。放射型CTは放射源によって様々な応用が考えられる。医用への応用では,X線CTがX線の吸収係数の分布により組織の形状を与えるのに対して,放射型CTでは放射源として体内に注入した放射性同位元素(RI)が特定の臓器や病変部に集積することを利用して,その分布から組織の機能の情報を得ることができる点で極めて有用である。しかし,実用化にあたっては分解能の向上,装置の低価格化,操作の単純化,測定時間の短縮などの問題点を解決しなければならない。
符号化開口CTは,放射源の断面ごとの投影像を得て,順次に3次元分布を求めるのではなく,放射源をごく少数の方向から観測して得た投影像を用いて一度に3次元の断層像を再構成できる点に特色がある。したがって,測定装置が簡単で,測定時間が短いという利点をもつ。
再構成においては開ロパターンの自己相関関数が近似的にδ関数で表されるという性質を利用する。筆者らはこのような性質をもつパターンとしてM配列2)を符号化開口に用い,RI分布を再構成する方法を示し,放射源分布を深さの異なる断層像として再構成できることを示した3)。しかし,3次元分布の再構成を行うには深さ方向の分解能が十分でない。これを改善するために開発した複数投影を用いる符号化開口CTの手法について報告する。
2.符号化開ロCTの原理
符号化開口CTの測定系は符号化開口,エリアセンサ,データ処理系から構成される(Fig.1)。ここでは符号化開口の面を開口面,エリアセンサの面を検出面,開口面に平行にとった測定対象の任意の断層と開口面の距離をその断層の深さと呼ぶことにする。
Fig.2に示すような座標系を設定し,断層,開口面,検出面内の位置をそれぞれベクトル_yx=(xl,x2),u=(u1、,u2)_,,r=(y1,yz)で表す。測定対象内の位置は断層内の位置.rと深さzの組(yzx,)で表すことができる。また,開口面の透過率分布をf(u)で表す。f(7)の値は実現が容易な0と1の二値とする。点(マ,z)に相対強度Sの点放射源があるとき,これによって生じる投影像p(Y)は次式で与えられる。
投影像p(一・Y)は,開ロパターンを(D+z)/z倍に拡大し,さらに一(D/z刀だけ平行移動したものである。
放射源が広がりをもつ場合,放射源の相対強度分布をS(x,z)とすると,投影像p(Y)は,
で与えられる。ここでdx,dx2を省略してdx2とした。また,積分範囲objは放射源が存在しうる範囲を表す。
放射源分布の断層像を得ることは,観測した投影像p(Y)を用いて(2)式をS(x,z)について解くことに帰着する。(2)式を解析的に解くことは難しいが,開ロパターンがランダムで自己相関関数が近似的にδ関数であれば,マッチトフィルタの考え方を用いた逆投影と呼ばれる簡単な演算でS(x,z)の推定値S(x,z)を得ることができる。
ここで,積分範囲detは投影像の検出された範囲,Aはその面積を表す。
3.実験
核医学の分野でよく用いられるRI, ssmTcを放射源とし,鉛板で作成したM配列符号化開口とガンマカメラを用いてRIから放射されるγ線による投影像を観測し,逆投影法を用いてRI分布を再構成する実験を行った。
用いたM配列符号化開口は3mm厚の鉛板に径2mmのピンホールを3mmピッチであけて製作したものである。放射源分布(ファントム)には,文字CとTの形に切った吸い取り紙にRI溶液を染み込ませたものを用いた。これらをそれぞれ深さ29cm,43cmの位置に開口からみて重なるように置いて放射源分布とした。
観測した投影像(128×128画素)から,(3)式を計算して放射源分布の断層像を得た。Fig.3にいくつかの深さを仮定して再構成した断層像を示す。再構成する深さがファントムを置いた深さにほぼ一致するときには,ほぼ正確な断層像が得られた。
しかし,その他の深さでは,実際にはファントムが存在していないのに,ファントムのぼけた像が現れ,誤った再構成像が得られた。これは,深さ方向の分解能が十分でないためと考えられる。
4.結像系の性能に関する考察4)
結像系の分解能を点放射源の再構成像PSF(Point Spread Function)で評価するものとすると,符号化開口CTでは,開ロバターンの自己相関関数の深さ方向での変化でPSFが表され,これは必ずしも鋭いピークをもたない。このように分解能が低い原因は,僅かに倍率の異なる開ロパターンの相関関数値が倍率の同じ開ロパターンに対する値と類似していることによっている。深さ方向の分解能を上げるには,深さ(=倍率)が異なる投影からの再構成像間の差異を大きくする必要がある。
いま,開口が(1+m)倍に拡大されて投影される場合を考える。開口の符号のピッチをpとする。このとき,原点からkピッチだけ離れた大きさα=p(1一ε)[0<ε〈1]の孔がc(m)の重なりを生ずるものとする。これが小さいほど再構成像間の差異が大きくなる。簡単な考察から,
が得られる。mが一定であるとして,これを小さくするには,p→小,ε→大,k→大とすればよい。つまり,①ピッチpを小さくするか,②符号の単位となる孔の大きさ(1一ε)を小さくするか,③観測範囲kを広くすれば良いことが分かる。一方,mが可変であれば,mは大きいほどよい。これは,④深さの違いがmにできるだけ大きい影響を与えるということで,
から,z=0(z>0)となる。
ここで,これらの条件について考察する。まず①については,開口を細かいピッチで製作して用いることを意味するが,製作上の困難は別にしても,S/Nが低下し,さらに,面検出器の分解能が高くなければ意味がない。
②については,これを直接実行すると,①同様,開口を通過するγ線の量が減り,S/Nの低下をもたらす。
③に関しては,相関関数を計算する範囲を大きくとることによって分解能の向上が図れるが,それを実現するには,面積の大きな面検出器を使用することが必要になり,また,投影の周辺部と中央部で倍率が変わる。
④については,z=0では意味がない。配置を工夫して,z
5.複数投影を用いた再構成アルゴリズム
5.1差動法
複数投影を用いることにより上記②を近似的に実現することができる。これには,開口あるいは検出器を深さ方向に僅かに移動して得られる2枚の投影の差(差動形)を用いればよい。つまり,倍率の僅かに異なる投影の差を得ることで,S/Nを低下させずに,孔を小さくした効果を実現することができる。逆投影法と差動型のPSFをFig.8(d)に示した。差動型のPSFは深さ方向にピークが鋭くなり,深さ方向の分解能が向上させられることが分かる。
差動形のアルゴリズムを用いて3次元の放射源分布を再構成した。放射源はFig.4に示すような形のアクリル容器にssmTcコロイドを満たしたものを用いた。この放射源の断層像を上から3個所順に再構成した結果がFig.5である。Fig,5(a)では上側の円柱の断面,Fig.5(b)では円柱と角柱の断面,Fig.5(c)では下側の円柱の断面が再構成されている。しかし,得られる断層像は断面の輪郭を再構成したものになっている。これは差動形のPSFが2階微分形になっているためである。このように,差動形では深さ方向の分解能が改善できるものの,放射源分布のエッジを強調した断層像が得られるという問題がある。
5.2立体視法5)・6)・7)
複数投影を用いるもう一つの方法は,上記③を実現することに相当して,奥行きを知るのによく用いられる両眼視を利用するものである。検出器をその面内で僅かに移動し,異なる方向で投影を観察し,その視差を利用して深さ方向の情報を得る方法である(立体視形)。
得られる複数の投影像から放射源分布を再構成するために,ここでは,放射源分布と投影像の関係式を行列で表現し,これを解く方法を用いた。つまり,放射源をn点でサンプリングしてn次元ベクトルで表したものをS,同様に(複数の)投影像をm点でサンプリングしてm次元ベクトルで表したものをpとおくと,Sとpの関係はm×n次元の行列Fを用いて,p=Fsと表すことができる。これは(1)式を行列を用いて書き表したものに他ならず,複数投影の場合も形式的に同じ形で書ける。放射源の再構成はこの式をSについて解くことである。しかし,これを直接解くことは行列Fが大きいと困難である。解Sの近似値を漸近的に解に近付ける繰り返し法を用いた。
Fig.6(a)は立体視と繰り返し法を用いたときの深さ方向のPSFを計算で求めたものである。また,放射性分布についてのシミュレーション結果をFig.6(b)(逆投影法),(c)(立体視法)に示す。繰り返し回数が0回のものは逆投影法に相当する。したがって,この方法は逆投影法に較べて深さ方向の分解能が向上できることがわかる。しかし,実際的な状況では行列が非常に大きくなり,計算時間や記憶容量の制約のため,計算が困難となる。
5.3最適フィルタ法8)・9)・10)
エリアセンサを深さ方向に動かして距離1)k(k=1,...,m)で得た複数投影に対する最適ブイルタを考える。投影像pk(Y)は,
と書ける。ここで,vk(Y)は観測雑音である。
複数投影を入力とし,インパルス応答をgk(xとする線形システムによる深さCの断層像の推定値S(x)は,
となる。元の分布Sの推定は推定値SがSにできるだけ近くなるよう,すなわち,推定の2乗平均誤差が最小となるようにgkを定めて行う。伝達関数を用いて表現するとFig.7のような系で再構成を行うことになる。最適フィルタ(一.Skは解析的に定めることができる。
投影像を1枚および複数枚用いたときの各種手法のPSFを計算により求め,逆投影法,Helstromフィルタ11),差動法,最適フィルタの性能の評価を行った。
投影像が1枚の場合には,Helstromフィルタと最適フィルタの分解能は逆投影法のそれより僅かに改善されるものの,両者であまり変わらないが(Fig.8(b)(c)),後者で2乗平均誤差が約30%低減される。
投影像を2枚用いた場合,差動法(Fig.8(d))では分解能は改善されるが,ピークのまわりに負の部分が生ずるため,2乗平均誤差が大きくなる。これに対して,最小2乗ブイルタ(Fig.10(e))では分解能が向上するとともに2乗平均誤差も低減している。
投影像を5枚用いた場合の最適ブイルタの結果(Fig.10(f))から,投影像の枚数を増やすことによって,分解能と精度がよくなったことがわかる。
6.おわりに
符号化開口CTの深さ方向の分解能を改善するために複数投影を用いる方法を開発した。複数投影を入力とする差動法,立体視法,最適ブイルタ法を考案し,数値シミュレーションや実験により,深さ方向の分解能が改善されることを示した。
実際の観測系について性能を評価すること,さらに,開ロパターンや複数投影像の観測方法など,観測系についての最適化を実現し,実際の測定に適用することが今後の課題である。
符号化開口CTはそれ自身が放射を出す放射源の3次元分布を非接触・非破壊で計測する放射型CTの一手法である')。放射型CTは放射源によって様々な応用が考えられる。医用への応用では,X線CTがX線の吸収係数の分布により組織の形状を与えるのに対して,放射型CTでは放射源として体内に注入した放射性同位元素(RI)が特定の臓器や病変部に集積することを利用して,その分布から組織の機能の情報を得ることができる点で極めて有用である。しかし,実用化にあたっては分解能の向上,装置の低価格化,操作の単純化,測定時間の短縮などの問題点を解決しなければならない。
符号化開口CTは,放射源の断面ごとの投影像を得て,順次に3次元分布を求めるのではなく,放射源をごく少数の方向から観測して得た投影像を用いて一度に3次元の断層像を再構成できる点に特色がある。したがって,測定装置が簡単で,測定時間が短いという利点をもつ。
再構成においては開ロパターンの自己相関関数が近似的にδ関数で表されるという性質を利用する。筆者らはこのような性質をもつパターンとしてM配列2)を符号化開口に用い,RI分布を再構成する方法を示し,放射源分布を深さの異なる断層像として再構成できることを示した3)。しかし,3次元分布の再構成を行うには深さ方向の分解能が十分でない。これを改善するために開発した複数投影を用いる符号化開口CTの手法について報告する。
2.符号化開ロCTの原理
符号化開口CTの測定系は符号化開口,エリアセンサ,データ処理系から構成される(Fig.1)。ここでは符号化開口の面を開口面,エリアセンサの面を検出面,開口面に平行にとった測定対象の任意の断層と開口面の距離をその断層の深さと呼ぶことにする。
Fig.2に示すような座標系を設定し,断層,開口面,検出面内の位置をそれぞれベクトル_yx=(xl,x2),u=(u1、,u2)_,,r=(y1,yz)で表す。測定対象内の位置は断層内の位置.rと深さzの組(yzx,)で表すことができる。また,開口面の透過率分布をf(u)で表す。f(7)の値は実現が容易な0と1の二値とする。点(マ,z)に相対強度Sの点放射源があるとき,これによって生じる投影像p(Y)は次式で与えられる。
投影像p(一・Y)は,開ロパターンを(D+z)/z倍に拡大し,さらに一(D/z刀だけ平行移動したものである。
放射源が広がりをもつ場合,放射源の相対強度分布をS(x,z)とすると,投影像p(Y)は,
で与えられる。ここでdx,dx2を省略してdx2とした。また,積分範囲objは放射源が存在しうる範囲を表す。
放射源分布の断層像を得ることは,観測した投影像p(Y)を用いて(2)式をS(x,z)について解くことに帰着する。(2)式を解析的に解くことは難しいが,開ロパターンがランダムで自己相関関数が近似的にδ関数であれば,マッチトフィルタの考え方を用いた逆投影と呼ばれる簡単な演算でS(x,z)の推定値S(x,z)を得ることができる。
ここで,積分範囲detは投影像の検出された範囲,Aはその面積を表す。
3.実験
核医学の分野でよく用いられるRI, ssmTcを放射源とし,鉛板で作成したM配列符号化開口とガンマカメラを用いてRIから放射されるγ線による投影像を観測し,逆投影法を用いてRI分布を再構成する実験を行った。
用いたM配列符号化開口は3mm厚の鉛板に径2mmのピンホールを3mmピッチであけて製作したものである。放射源分布(ファントム)には,文字CとTの形に切った吸い取り紙にRI溶液を染み込ませたものを用いた。これらをそれぞれ深さ29cm,43cmの位置に開口からみて重なるように置いて放射源分布とした。
観測した投影像(128×128画素)から,(3)式を計算して放射源分布の断層像を得た。Fig.3にいくつかの深さを仮定して再構成した断層像を示す。再構成する深さがファントムを置いた深さにほぼ一致するときには,ほぼ正確な断層像が得られた。
しかし,その他の深さでは,実際にはファントムが存在していないのに,ファントムのぼけた像が現れ,誤った再構成像が得られた。これは,深さ方向の分解能が十分でないためと考えられる。
4.結像系の性能に関する考察4)
結像系の分解能を点放射源の再構成像PSF(Point Spread Function)で評価するものとすると,符号化開口CTでは,開ロバターンの自己相関関数の深さ方向での変化でPSFが表され,これは必ずしも鋭いピークをもたない。このように分解能が低い原因は,僅かに倍率の異なる開ロパターンの相関関数値が倍率の同じ開ロパターンに対する値と類似していることによっている。深さ方向の分解能を上げるには,深さ(=倍率)が異なる投影からの再構成像間の差異を大きくする必要がある。
いま,開口が(1+m)倍に拡大されて投影される場合を考える。開口の符号のピッチをpとする。このとき,原点からkピッチだけ離れた大きさα=p(1一ε)[0<ε〈1]の孔がc(m)の重なりを生ずるものとする。これが小さいほど再構成像間の差異が大きくなる。簡単な考察から,
が得られる。mが一定であるとして,これを小さくするには,p→小,ε→大,k→大とすればよい。つまり,①ピッチpを小さくするか,②符号の単位となる孔の大きさ(1一ε)を小さくするか,③観測範囲kを広くすれば良いことが分かる。一方,mが可変であれば,mは大きいほどよい。これは,④深さの違いがmにできるだけ大きい影響を与えるということで,
から,z=0(z>0)となる。
ここで,これらの条件について考察する。まず①については,開口を細かいピッチで製作して用いることを意味するが,製作上の困難は別にしても,S/Nが低下し,さらに,面検出器の分解能が高くなければ意味がない。
②については,これを直接実行すると,①同様,開口を通過するγ線の量が減り,S/Nの低下をもたらす。
③に関しては,相関関数を計算する範囲を大きくとることによって分解能の向上が図れるが,それを実現するには,面積の大きな面検出器を使用することが必要になり,また,投影の周辺部と中央部で倍率が変わる。
④については,z=0では意味がない。配置を工夫して,z
5.複数投影を用いた再構成アルゴリズム
5.1差動法
複数投影を用いることにより上記②を近似的に実現することができる。これには,開口あるいは検出器を深さ方向に僅かに移動して得られる2枚の投影の差(差動形)を用いればよい。つまり,倍率の僅かに異なる投影の差を得ることで,S/Nを低下させずに,孔を小さくした効果を実現することができる。逆投影法と差動型のPSFをFig.8(d)に示した。差動型のPSFは深さ方向にピークが鋭くなり,深さ方向の分解能が向上させられることが分かる。
差動形のアルゴリズムを用いて3次元の放射源分布を再構成した。放射源はFig.4に示すような形のアクリル容器にssmTcコロイドを満たしたものを用いた。この放射源の断層像を上から3個所順に再構成した結果がFig.5である。Fig,5(a)では上側の円柱の断面,Fig.5(b)では円柱と角柱の断面,Fig.5(c)では下側の円柱の断面が再構成されている。しかし,得られる断層像は断面の輪郭を再構成したものになっている。これは差動形のPSFが2階微分形になっているためである。このように,差動形では深さ方向の分解能が改善できるものの,放射源分布のエッジを強調した断層像が得られるという問題がある。
5.2立体視法5)・6)・7)
複数投影を用いるもう一つの方法は,上記③を実現することに相当して,奥行きを知るのによく用いられる両眼視を利用するものである。検出器をその面内で僅かに移動し,異なる方向で投影を観察し,その視差を利用して深さ方向の情報を得る方法である(立体視形)。
得られる複数の投影像から放射源分布を再構成するために,ここでは,放射源分布と投影像の関係式を行列で表現し,これを解く方法を用いた。つまり,放射源をn点でサンプリングしてn次元ベクトルで表したものをS,同様に(複数の)投影像をm点でサンプリングしてm次元ベクトルで表したものをpとおくと,Sとpの関係はm×n次元の行列Fを用いて,p=Fsと表すことができる。これは(1)式を行列を用いて書き表したものに他ならず,複数投影の場合も形式的に同じ形で書ける。放射源の再構成はこの式をSについて解くことである。しかし,これを直接解くことは行列Fが大きいと困難である。解Sの近似値を漸近的に解に近付ける繰り返し法を用いた。
Fig.6(a)は立体視と繰り返し法を用いたときの深さ方向のPSFを計算で求めたものである。また,放射性分布についてのシミュレーション結果をFig.6(b)(逆投影法),(c)(立体視法)に示す。繰り返し回数が0回のものは逆投影法に相当する。したがって,この方法は逆投影法に較べて深さ方向の分解能が向上できることがわかる。しかし,実際的な状況では行列が非常に大きくなり,計算時間や記憶容量の制約のため,計算が困難となる。
5.3最適フィルタ法8)・9)・10)
エリアセンサを深さ方向に動かして距離1)k(k=1,...,m)で得た複数投影に対する最適ブイルタを考える。投影像pk(Y)は,
と書ける。ここで,vk(Y)は観測雑音である。
複数投影を入力とし,インパルス応答をgk(xとする線形システムによる深さCの断層像の推定値S(x)は,
となる。元の分布Sの推定は推定値SがSにできるだけ近くなるよう,すなわち,推定の2乗平均誤差が最小となるようにgkを定めて行う。伝達関数を用いて表現するとFig.7のような系で再構成を行うことになる。最適フィルタ(一.Skは解析的に定めることができる。
投影像を1枚および複数枚用いたときの各種手法のPSFを計算により求め,逆投影法,Helstromフィルタ11),差動法,最適フィルタの性能の評価を行った。
投影像が1枚の場合には,Helstromフィルタと最適フィルタの分解能は逆投影法のそれより僅かに改善されるものの,両者であまり変わらないが(Fig.8(b)(c)),後者で2乗平均誤差が約30%低減される。
投影像を2枚用いた場合,差動法(Fig.8(d))では分解能は改善されるが,ピークのまわりに負の部分が生ずるため,2乗平均誤差が大きくなる。これに対して,最小2乗ブイルタ(Fig.10(e))では分解能が向上するとともに2乗平均誤差も低減している。
投影像を5枚用いた場合の最適ブイルタの結果(Fig.10(f))から,投影像の枚数を増やすことによって,分解能と精度がよくなったことがわかる。
6.おわりに
符号化開口CTの深さ方向の分解能を改善するために複数投影を用いる方法を開発した。複数投影を入力とする差動法,立体視法,最適ブイルタ法を考案し,数値シミュレーションや実験により,深さ方向の分解能が改善されることを示した。
実際の観測系について性能を評価すること,さらに,開ロパターンや複数投影像の観測方法など,観測系についての最適化を実現し,実際の測定に適用することが今後の課題である。