JP6080588B2 - 画像処理装置、画像処理方法、及びコンピュータプログラム - Google Patents

Description

本発明は、画像処理装置、画像処理方法、及びコンピュータプログラムに関し、特に、断層画像を生成するために用いて好適なものである。

Ｘ線を用いた断層画像による診断装置は1970年代に製品化されてから４０年以上経つが、今なお診断技術の中心として発展活躍をしている。近年では、コンピュータのハードウェアの進化に伴って、処理の高速化とメモリ量の増加とが実現されることにより、従来では製品化が難しいと考えられていた逐次近似再構成法により画像を作成する診断装置が製品化されている。

逐次近似再構成法は、ＦＢＰ（Filtered Back Projection）法とは異なり、代数的に画像の再構成を行う方法である。逐次近似再構成法は、再投影と逆投影とを繰り返すことで断層画像の投影画像が実際に撮影された投影画像と一致するように逐次的に断層画像を求めていく方法である。

逐次近似再構成法では、ＦＢＰ（Filtered Back Projection）法よりもはるかに多くのメモリと再構成時間とを費やす。しかしながら、再構成に様々な数学モデルを組み込むことができるため、散乱線補正、アーチファクト補正、検出器分解能補正、及びノイズ低減等を考慮することができる。これにより、従来よりもＳＮ比が高くアーチファクトが抑制された断層画像を生成することができる。

ところで、逐次近似再構成法では、逆投影と再投影とを逐次的に繰り返して断層画像の再構成を行う。このため、逆投影と再投影とが正確に行われないと、誤差が蓄積し画像が破壊されてしまう。特に、断層画像の端ではＦＯＶ（Field Of View）の欠けにより、全角度の投影画像が用いられない場合がある。このため、逐次近似の反復の度に、逆投影時に用いた投影画像の枚数で画素値を規格化しなければ、逐次近似再構成の解が収束せずに画素値が発散してしまう。また、断層画像の端において再投影を行うと、断層の再構成の範囲によっては、全高さの断層の画像が用いられない場合がある。この場合には、逐次近似の反復の度に、再投影時に用いた断層画像の枚数で画素値を規格化しなければ、前述したのと同様に、逐次近似再構成の解が収束せずに画素値が発散してしまう。

特許文献１においては、画素値１（すなわち単位ベクトル）の投影画像・断層画像をそれぞれ逆投影・再投影した画像を用いて、逐次近似の反復の度に、断層画像・再投影画像を規格化する方法が開示されている。

国際公開第２０１０／１００９９６号

しかしながら、特許文献１の方法では、予め、画素値１の投影画像・断層画像をそれぞれ逆投影・再投影した画像を生成しておく必要がある。この場合、投影画像と断層画像のメモリ容量を、それらの画像を生成しない場合に比べ、倍確保しなくてはならない。例えば、断層画像と投影画像の画像サイズが１０００×１０００、枚数が５００枚とすると４ＧＢのメモリを余分に確保する必要がある。

ハイエンドのＣＴ装置等では、このような大きなメモリ容量が必要になっても大きな問題にはならないかもしれない。しかしながら、ローエンドのＣＴ装置やトモシンセシスでは、このような大きなメモリ容量は、再構成に用いる演算装置の価格に直接影響するので問題になる。特に最近では、ＧＰＵ（Graphics Processing Unit）による画像処理が医用画像にも用いられることが多くなり、ローコストな市販のＧＰＵでも十分な速度の再構成演算が可能になっている。こうした市販のＧＰＵを用いてシステムを構築するためにはメモリ容量を仕様範囲内に収めることが必要である。
本発明は、このような問題点に鑑みてなされたものであり、断層画像を再構成する際に必要なメモリ容量を削減することを目的とする。

本発明の画像処理装置は、寝台の面方向の一つの方向と前記寝台の面に対する法線方向とにより定まる平面上を回動する放射線発生源から複数の照射角度で照射された放射線に基づいて平面検出器により検出された投影画像を取得する取得手段と、前記取得手段により取得された投影画像に基づいて逆投影を行って断層画像を生成する逆投影手段と、前記投影画像のうち、前記逆投影手段により前記断層画像が生成された際に用いた前記投影画像の割合を表す逆投影利用率であって、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に値が依存しない逆投影利用率を、前記断層画像の所定の画素のそれぞれについて計算する逆投影利用率計算手段と、前記逆投影手段により生成された断層画像の画素値を前記逆投影利用率で規格化する逆投影利用率演算手段と、を有し、前記逆投影利用率計算手段は、前記断層画像の画素値のうち、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に沿う位置にある複数の画素の画素値については１つの前記逆投影利用率を計算することを特徴とする。

本発明によれば、断層画像を再構成する際に必要なメモリ容量を削減することができる。

画像処理装置の構成を示す図である。 再投影処理の方法を説明する図である。 逐次近似再構成の流れを説明するフローチャートである。 逆投影利用率の計算方法の概念を説明する図である。 再投影利用率の計算方法の概念を説明する図である。 逆投影利用率を説明する図である。 投影利用率を説明する図である。

以下、図面を参照して、本発明の一実施形態を詳細に説明する。尚、以下の説明では、画像処理装置がトモシンセシスである場合を例に挙げて説明する。しかしながら、画像処理装置はトモシンセシスに限定されるものではない。例えば、ＣＴ装置に画像処理装置を適用することが可能である。また、以下の説明では、ＳＩＲＴ法による逐次近似再構成を例に挙げて説明する。しかしながら、逐次近似再構成法はＳＩＲＴ法に限定されるものではない。例えば、ＡＲＴ法やＭＢＩＲ法、ＭＬＥＭ法等の逐次近似再構成法を採用することが可能である。

図１は、画像処理装置１００の機能構成の一例を示す図である。
Ｘ線管１０１は、放射線を発生する放射線発生源の一例であり、複数の照射角度でＸ線を照射する。寝台１０３は被写体１０２を寝かせる台である。Ｘ線検出器１０６は、平面検出器の一例であり、Ｘ線管１０１から寝台１０３に向かって照射されたＸ線を受信してＸ線画像を取得する。機構制御部１０５は、Ｘ線管１０１の位置とＸ線検出器１０６の位置との制御を行う。撮影制御部１０４は、Ｘ線検出器１０６を電気的に制御して、Ｘ線画像を取得する。Ｘ線発生装置１０７は、Ｘ線管１０１の動作を電気的に制御して、所定の条件でＸ線を発生させる。Ｘ線撮影システム制御部１０８は、機構制御部１０５と撮影制御部１０４とＸ線発生装置１０７とを制御して、複数の照射角度における複数のＸ線画像を取得する。

Ｘ線撮影システム制御部１０８は、画像処理部１０９及び画像保存部１１２を含む。Ｘ線撮影システム制御部１０８には、１又は複数のコンピュータが内蔵される。コンピュータには、例えば、ＣＰＵ等の主制御手段、ＲＯＭ（Read Only Memory）、及びＲＡＭ（Random Access Memory）等の記憶手段が具備される。また、コンピュータには、ＧＰＵ（Graphics Processing Unit）等のグラフィック制御手段と、ネットワークカード等の通信手段と、キーボード、ディスプレイ又はタッチパネル等の入出力手段等が具備されていてもよい。尚、これら各構成手段は、バス等により相互に接続される。また、これら各構成手段は、主制御手段が、記憶手段に記憶されたプログラムを実行することで制御される。また、Ｘ線撮影システム制御部１０８には、撮影された断層画像等を表示するモニタ１１０と、ユーザが操作する操作部１１１とが備えられている。

Ｘ線撮影システム制御部１０８は、撮影制御部１０４を介して、Ｘ線検出器１０６から取得した様々な照射角度における複数のＸ線画像を画像処理部１０９に入力する。以下の説明では、「Ｘ線画像」を必要に応じて「投影画像」と称する。投影画像に対しては、Ｘ線検出器１０６の欠陥画素や暗電流、Ｘ線管１０１に起因する照射ムラの補正、対数変換等、Ｘ線検出器で一般的に行われる前処理が、画像処理部１０９に入力される前に予め行われる。
画像処理部１０９は、Ｘ線撮影システム制御部１０８の指示に従い、取得した投影画像を再構成し、断層画像を生成する。このために、画像処理部１０９は、差分処理部１１３、逆投影部１１４、加算処理部１１５、再投影部１１６、逆投影利用率計算部１１７、再投影利用率計算部１１８、逆投影利用率メモリ１１９、再投影利用率メモリ１２０、逆投影利用率演算部１２１、及び再投影利用率演算部１２２を有する。

差分処理部１１３は、２つの画像の差分画像を生成する。加算処理部１１５は、２つ以上の画像の加算操作を行う。
再投影部１１６は、Ｘ線による撮影を数値的に模擬した再投影処理（投影処理または順投影処理）を行う。図２は、再投影処理の方法の一例を説明する図である。図２を参照しながら、再投影処理の具体的な方法の一例を説明する。図２において、格子上の点が画素を示し、格子の間隔が画素ピッチに相当する。ある照射角度βにおけるＸ線２０２と交わるＸ線検出器１０６の画素の画素値をｆ_iとし、断層画像２０１の画素の画素値をｇ_jとする。また、断層画像２０１の画素の画素値ｇ_jの、投影画像２０３の画素の画素値ｆ_iへの寄与率をＣ_ijとする。このとき、再投影処理により得られる投影画像２０３の画素値ｆ_iは、以下の式（１）で表される。尚、ここでは、図２に示すように、Ｘ線管１０１の焦点が鉛直下方を向くときに照射角度βが０度になるものとする。

ここで、ｉは投影画像２０３のｉ番目の画素を示し、ｊは、断層画像２０１のｊ番目の画素を示す。また、Ｊは断層画像２０１の全画素数である。前述したように、Ｃ_ijは、断層画像２０１の画素の画素値ｇ_jの、投影画像２０３の画素の画素値ｆ_iへの寄与率を表し、以下の説明では、「Ｃ_ij」を必要に応じて「寄与率」と略称する。ここで、Ｘ線管の焦点と、画素値ｆ_iを有する投影画像２０３の画素とを相互に結ぶ直線が、例えば、画素値ｇ_jを有する断層画像２０１の画素の座標（x，z）の近傍の座標（x，z＋0.2）、（x−0.4，z）を通るとする。この場合、線形補間を用いると、寄与率Ｃ_ijは０．４８（＝（1‐0.2）×（1‐0.4））と計算される。寄与率Ｃ_ijが０．４８であることは、概念的には、断層画像２０１の画素の画素値ｇ_jが、投影画像２０３の画素の画素値ｆ_iに４８％寄与していることを表す。寄与率Ｃ_ijの計算方法は、断層画像２０１の画素の画素値ｇ_jの、投影画像２０３の画素の画素値ｆ_iへの寄与の度合いに比例する計算を行う方法であれば、前述した方法に限定されない。例えば、断層画像の画素を横切る長さから寄与率Ｃ_ijを計算する方法や、バイキュービックやさらに高次の補間を利用して寄与率Ｃ_ijを計算する方法や、断層画像の１画素を投影画像に投影したときの形状の面積から寄与率Ｃ_ijを計算する方法等を採用できる。記述を簡便にするために、以降では式（１）を、以下の式（２）の形で表す。

ここで、ｆは、投影画像２０３の総画素数であるＩ個の要素からなる投影画像ベクトルである。投影画像ベクトルｆの各成分はｆ_i（画素値）である。Ｈは、Ｉ行Ｊ列の投影行列である。投影行列Ｈの各成分はＣ_ijである。ｇは、断層画像２０１の総画素数であるＪ個の要素からなる断層画像ベクトルである。断層画像ベクトルｇの各成分はｇ_j（画素値）である。
図１の説明に戻り、逆投影部１１４は、再投影部１１６とは逆に、投影画像２０３を、Ｘ線管１０１の焦点とＸ線検出器１０６とを相互に結ぶ直線上で、断層画像２０１に再配置する逆投影処理を行う。逆投影処理により得られる断層画像２０１の画素値ｇ_jは、以下の式（３）で表される。

ここで、Ｉは、投影画像２０３の全画素数である。式（２）と同様に、式（３）をベクトル表記すると、以下の式（４）のようになる。

ここで、逆投影行列Ｈ^Tは投影行列Ｈの転置行列である。
再投影利用率計算部１１８は、断層画像２０１から投影画像２０３へ再投影を行うときに、何枚の断層画像を用いたかを表す再投影利用率を計算する。再投影利用率の具体的な計算方法は後述する。
逆投影利用率計算部１１７は、投影画像２０３からの断層画像２０１へ逆投影を行うときに、何枚の投影画像を用いたかを表す逆投影利用率を計算する。逆投影利用率の具体的な計算方法は後述する。

再投影利用率メモリ１２０は、再投影利用率計算部１１８で計算された再投影利用率を記憶する。
逆投影利用率メモリ１１９は、逆投影利用率計算部１１７で計算された逆投影利用率を記憶する。

再投影利用率演算部１２２は、再投影利用率メモリ１２０に記憶された再投影利用率で投影画像２０３の各画素の画素値を除算する演算処理を行う。
逆投影利用率演算部１２１は、逆投影利用率メモリ１１９に記憶された逆投影利用率で断層画像２０１の各画素の画素値を除算する演算処理を行う。

次に、図３のフローチャートを参照しながら、図１に示す画像処理装置１００における逐次近似再構成の流れの一例について説明する。
まず、ステップＳ３０１で、Ｘ線検出器１０６は、投影画像を取得する。これは、Ｘ線管１０１から照射されるＸ線の照射角度βを変えながら、被写体１０２をＸ線で撮影することにより行われる。照射角度βと撮影枚数は、装置の性能や仕様によるが、例えば、照射角度βを−４０度から＋４０度まで１度ずつ変えながら、８０枚の投影画像を１５ＦＰＳで撮影すると６秒程度で投影画像の収集ができる。Ｘ線の撮影条件も任意の値を設定可能である。例えば、胸部等の撮影では１００ｋＶ、１ｍＡｓ程度の撮影条件で撮影画像の収集を行えばよい。また、Ｘ線検出器１０６とＸ線管１０１との間の距離は透視撮影装置や一般撮影装置の設定範囲である１００ｃｍ〜１５０ｃｍ程度に設定される。

Ｘ線管１０１から照射されるＸ線（ベクトル）の水平成分の方向とは反対方向に、Ｘ線検出器１０６を（水平面に沿って）平行移動させる。このときの移動量は、照射角度βが０度のときの「Ｘ線管１０１の回転中心とＸ線検出器１０６の中心との距離」をＰとすれば、Ｐtanβで与えられる。このようにしてＸ線検出器１０６を平行移動させれば、Ｘ線管１０１からのＸ線の照射方向が変わっても、Ｘ線管１０１の基準軸はＸ線検出器１０６の中心を常に通る。

得られた一連の投影画像は、前述した前処理が行われた後、画像処理部１０９に入力される。前処理（対数変換）により、投影画像２０３の画素値はＸ線減弱係数を線積分したものになる。このＸ線減弱係数の加法性に基づいて投影画像２０３の再構成が行われる。
ステップＳ３０９では、逆投影利用率計算部１１７は、再構成に先立ち、逆投影利用率を予め計算する。
図６は、逆投影利用率を概念的に説明する図である。図６において、６０１は各投影角度でのＸ線検出器、６０２は断層画像、６０３は各照射角度でのＸ線管、６０４は断層画像中央への逆投影、６０５は断層画像端への逆投影である。逆投影操作を行う場合、断層画像中央への逆投影６０４の場合、全投影画像（ここでは簡単のため３つの投影）が寄与するため、逆投影利用率は１００％（＝（３／３）×１００）である。これに対して、断層画像端への逆投影５０５では最も大きな照射角度では、断層画像の画素とＸ線管とを結ぶ線がＸ線検出器５０１からはみ出してしまう。この場合、投影画像２枚分しか寄与しないため、逆投影利用率は６７％（＝（２／３）×１００）となる。

図４は、逆投影利用率の計算方法の概念を説明する図である。ステップＳ３０１で取得された投影画像２０３に基づき、断層画像２０１を再構成するには、断層画像２０１を生成する範囲を決めなくてはならない。この範囲を再構成エリア４０１として図３に図示する（矩形状の太い実線を参照）。

前述したようにトモシンセシスでは、Ｘ線検出器１０６とＸ線管１０１を移動させながら投影画像２０３を取得する。図２において、実線で示す１０６ａは、照射角度βが０度のときのＸ線検出器１０６の位置を示し、破線で示す１０６ｂは、照射角度βが２０度のときのＸ線検出器１０６の位置を示す。
図４に示すように再構成エリア４０１を設定した場合、照射角度βが０度のときには再構成エリア４０１の全てを覆う投影画像２０３がＸ線検出器１０６によって撮影される。しかしながら、照射角度βが２０度のときには、再構成エリア４０１の一部の領域４０２に対応するＸ線検出器１０６上の画素が存在しないため、この領域４０２については逆投影することができない。
単純にこの領域４０２の画素値を０として、再構成の反復を進めていくと、領域４０２の部分の断層画像２０１の画素値がオーバーフロー（発散）してしまい、診断価値が無い部分となってしまう。

これを防ぐには、特許文献１に示されるように、逆投影するたびに、画素の逆投影に用いた投影画像２０３の枚数（逆投影利用率）で断層画像２０１の画素値を規格化すればよい。しかしながら、特許文献１に記載の逆投影利用率は、全ての照射角度βで全ての断層画像２０１を逆投影しないと分からないため、予め計算してメモリに記憶しておく必要がある。この逆投影利用率は、逆投影の対象となる画素に対する投影画像の寄与率を示す値であり、また逆投影に用いられる投影画像の枚数の割合を示す値である。
特許文献１に記載のように、断層画像２０１の全ての画素毎の逆投影利用率を全てメモリに記憶すると、断層画像２０１と同じサイズのメモリ容量を確保する必要がある。そこで本実施形態では、次のようにして逆投影利用率を計算する。

図２に示すように、Ｘ線管１０１の焦点の座標を（Xs，Ys，Zs）とする。また、被写体１０２が寝ている寝台１０３の長手方向と平行な軸をＸ軸とし、寝台１０３の面に対する法線方向の軸をＺ軸とする。また、寝台１０３の幅方向（紙面に垂直な方向）の軸をＹ軸とする。さらに、断層画像２０１の座標を（Xr，Yr，Zr）とし、Ｘ線検出器１０６の画素、すなわち投影画像２０３の座標を（Xd，Yd，Zd）とする。

逆投影では、Ｘ線管１０１の焦点と断層画像２０１の座標とを相互に結ぶ直線がＸ線検出器１０６と交わる点（投影画像２０３の座標）を参照する。Ｘ線管１０１の焦点のＺ軸の座標がZsであり、断層画像２０１のＺ軸の座標がZrであり、Ｘ線検出器１０６のＺ軸の座標がZdであるとき、逆投影の際に参照する投影画像２０３のＸ軸の座標Xdは、以下の式（５）で表される。

なお、本実施形態の画像処理装置１００をトモシンセシスに適用した場合には投影画像２０３のＺ軸の座標Zdは常に一定であるので、この場合には式（５）の計算が容易になる。
座標（Xr，Yr，Zr）の断層画像２０１を逆投影により再構成する際に参照する投影画像２０３のＸ軸の座標Xdを全ての照射角度β（全投影画像）で計算する。そして、当該座標Xdの位置にＸ線検出器１０６の画素が存在する数をカウントすれば逆投影利用率を計算できる。尚、以下の説明では、「座標（Xr，Yr，Zr）の断層画像２０１を逆投影により再構成する際に参照する投影画像２０３のＸ軸の座標Xd」を必要に応じて「逆投影時参照Ｘ軸座標Xd」と称する。

式（５）において、逆投影時参照Ｘ軸座標Xdは、断層画像２０１のＹ軸の座標Yrに依存しない。すなわち、ある照射角度βにおいて、Ｘ軸の座標XrとＺ軸の座標Zrとが同じ画素の断層画像２０１を逆投影により再構成する際に参照する逆投影時参照Ｘ軸座標Xdは全て同じである。したがって、逆投影利用率も、断層画像２０１のＹ軸の座標Yrに依存しない。
そこで、本実施形態では、逆投影により再構成する断層画像２０１のＹ軸の所定の一軸上（Ｙ軸に平行な所定の一軸上）でのみ逆投影利用率を計算する。
逆投影利用率ＢＵ（Xr，Zr）は、以下の式（６）により計算される。

ここで、ＢＵ（Xr，Zr）は逆投影利用率であり、ＢＣ（Xr，Zr）は式（５）で計算した逆投影時参照Ｘ軸座標Xdの位置にＸ線検出器１０６の画素が存在する数である。また、ViewNumberは全投影画像の枚数である。

ステップＳ３０９において逆投影利用率計算部１１７が式（６）の計算を行って導出した逆投影利用率ＢＵ（Xr，Zr）は、ステップＳ３１０において、逆投影利用率メモリ１１９に記憶される。この逆投影利用率メモリ１１９のメモリ容量は、断層画像２０１のＹ軸の画素数をＷとした場合、断層画像２０１のメモリ容量の１／Ｗになる。例えば、断層画像２０１が１０００×１０００画素であり、断層画像２０１の枚数が５００枚の場合、断層画像２０１の全ての逆投影利用率を記憶すると４ＧＢのメモリ容量が必要なところ、その１０００分の１である４ＭＢまでメモリ容量を削減することができる。また、同様に逆投影利用率ＢＵ（Xr，Zr）の計算コストも大幅に削減できる。

続いてステップＳ３１１において、再投影利用率計算部１１８は、再構成に先立ち、再投影利用率を予め計算する。
図７は、投影利用率を概念的に示す図である。図７において、７０１はＸ線検出器、７０２は断層画像、７０３はＸ線管、７０４はＸ線検出器中央への投影線、７０５はＸ線検出器端への投影線である。投影操作を行う場合、Ｘ線検出器中央への投影７０４の場合、全断層画像（ここでは簡単のため５枚）が寄与するため、投影利用率は１００％（＝（５／５）×１００）である。これに対して、Ｘ線検出器端への投影７０５では断層画像２枚分しか寄与しないため、投影利用率は４０％（＝（２／５）×１００）となる。
図５は、再投影利用率の計算方法の概念を説明する図である。
図５において、実線で示す１０６ｃは、照射角度βが０度のときのＸ線検出器１０６の位置を示し、破線で示す１０６ｄは、照射角度βが２０度のときのＸ線検出器１０６の位置を示す。

図５に示すように再構成エリア５０１を設定した場合、例えば、照射角度βが２０度のときには、再投影時に、領域５０２で示される部分の断層画像２０１が存在しないため、この領域５０２については再投影することができない。
この問題を解決するために、再構成エリアを充分に広くとれば解決できるように思えるが、再構成エリアを広く取るとメモリ容量を多く必要とするだけでなく、図４に示した領域４０２が大きくなってしまい、逆投影時の欠損部分が大きくなってしまう。

本来は、全ての照射角度βにおけるＦＯＶが重なる範囲に再構成エリアを設定するのが計算上は好ましい。しかし、ＦＯＶが重なる範囲に再構成エリアを限定すると、トモシンセシスのような撮影では変則的な撮影動作を行うのでＦＯＶが複雑な形になり、再構成エリアが小さくなってしまう。また、トモシンセシスでは自由な被写体の配置が可能なため、図４や図５に示すように被写体１０２がＦＯＶをはみ出してしまうことが頻繁にある。

そこで、本実施形態では、逆投影時と同様に再投影時にも、再投影により得られた投影画像（再投影画像）の画素値を再投影利用率で再投影画像の画素値を規格化する。この再投影利用率は、投影の対象となる画素に対する再構成画像の寄与率を示す値であり、また投影に用いられる再構成画像の枚数の割合を示す値である。
ただし、逆投影の場合と同様に、特許文献１に記載のようにして、再投影画像の全ての画素毎の再投影利用率を全てメモリに記憶すると、投影画像と同じサイズのメモリ容量を確保する必要がある。そこで本実施形態では、次のようにして再投影利用率を計算する。

図５に示すように、Ｘ線管１０１の焦点の座標を（Xs，Ys，Zs）とする。また、被写体１０２が寝ている寝台１０３の長手方向と平行な軸をＸ軸とし、寝台１０３の面に対しする法線方向の軸をＺ軸とする。また、寝台１０３の幅方向（紙面に垂直な方向）の軸をＹ軸とする。断層画像２０１の座標を（Xｒ，Yr，Zr）とし、Ｘ線検出器１０６の画素、すなわち投影画像２０３の座標を（Xd，Yd，Zd）とする。

再投影では、Ｘ線管１０１の焦点とＸ線検出器１０６の座標とを相互に結ぶ直線上で断層画像２０１の画素値を線積分する。実際には断層画像２０１は離散的なため、直線上で断層画像２０１の座標を参照して画素値を積算する。Ｘ線管１０１の焦点のＺ軸の座標がZsであり、Ｘ線検出器１０６上の再投影画像のＺ軸の座標がZdであるとき、再投影の際に参照する断層画像２０１のＸ軸の座標Xrは、以下の式（７）で表される。

座標（Xd，Yd，Zd）の投影画像２０３（再投影画像）を再投影により再構成する際に参照すべき断層画像２０１のＸ軸の座標Xrを全ての照射角度β（全断層画像）で計算する。そして、当該座標Xrの位置に断層画像２０１の画素が存在する数をカウントすれば再投影利用率を計算できる。尚、以下の説明では、「座標（Xd，Yd，Zd）の投影画像２０３（再投影画像）を再投影により再構成する際に参照する断層画像２０１のＸ軸の座標Xr」を必要に応じて「再投影時参照Ｘ軸座標Xr」と称する。

式（７）において、再投影時参照Ｘ軸座標Xrは、再投影画像のＹ軸の座標Ydに依存しない。すなわち、ある照射角度βにおいて、Ｘ軸の座標XdとＺ軸の座標Zdとが同じ画素の投影画像２０３（再投影画像）を再投影により再構成する際に参照する再投影時参照Ｘ軸座標Xrは全て同じである。したがって、再投影利用率も、再投影画像のＹ軸の座標Ydに依存しない。
そこで、本実施形態では、再投影により再構成する投影画像２０３（再投影画像）のＹ軸の所定の一軸上（Ｙ軸に平行な所定の一軸上）でのみ再投影利用率を計算する。
再投影利用率ＦＵ（Xd，Zd）は、以下の式（８）により計算される。

ここで、ＦＵ（Xd，Zd）は再投影利用率である。ＦＣ（Xd，Zd）は式（７）で計算した再投影時参照Ｘ軸座標Xrの位置に断層画像２０１の画素が（欠落せずに）存在する数である。また、SliceNumberは全断層画像の枚数である。

ステップＳ３１１において再投影利用率計算部１１８が式（８）の計算を行って導出した再投影利用率ＦＵ（Xd，Zd）は、ステップＳ３１２において、再投影利用率メモリ１２０に記憶される。この再投影利用率メモリ１２０のメモリ容量は、投影画像２０３のＹ軸の画素数をＷとした場合、投影画像２０３のメモリ容量の１／Ｗになる。例えば、投影画像２０３が１０００×１０００画素であり、投影画像２０３の枚数が１２５枚の場合、投影画像２０３の全ての再投影利用率を記憶すると１ＧＢのメモリが必要なところ、その１０００分の１である１ＭＢまでメモリ容量を削減することができる。また、同様に再投影利用率ＦＵ（Xd，Zd）の計算コストも大幅に削減できる。

ステップＳ３０７において、再投影部１１６は再投影を行って再投影画像ベクトルを生成する。次に、ステップＳ３０８において、再投影利用率演算部１２２は、ステップＳ３０７で生成された再投影画像ベクトルを、再投影利用率ＦＵ（Xd，Zd）で除算して再投影画像ベクトルを規格化する。規格化された再投影画像ベクトルは「Ｈｇ^k／ＦＵ」となる。前述したように、Ｈは、Ｉ行Ｊ列の投影行列で、各成分はＣ_ijである（Ｉは投影画像２０３の総画素数であり、Ｊは断層画像２０１の総画素数である）。ｇは、断層画像２０１の総画素数であるＪ個の要素からなる断層画像ベクトルであり、各成分はｇ_j（画素値）である。ｋはｋ回目の反復を示す。尚、ｋ＝０での初期値ｇ⁰は負でない適当な値を与えてもよいし、ステップＳ３０１で取得された投影画像ベクトルｆを逆投影したものを初期値として用いてもよい。

ステップＳ３０２において、差分処理部１１３は、ステップＳ３０１で取得された投影画像ベクトルｆから、ステップＳ３０８で規格化された最新の再投影画像ベクトルＨｇ^k／ＦＵを減算して、投影画像ベクトルの差分ｆ−Ｈｇ^k／ＦＵを導出する。前述したように、投影画像ベクトルｆは、Ｉ個の要素からなる投影画像ベクトルであり、その各成分はｆ_iである。
次に、ステップＳ３０３において、逆投影部１１４は、ステップＳ３０２で得られた投影画像ベクトルの差分ｆ−Ｈｇ^k／ＦＵを逆投影する。

次に、ステップＳ３０４において、逆投影利用率演算部１２１は、ステップＳ３０３で得られた投影画像ベクトルの差分ｆ−Ｈｇ^k／ＦＵの逆投影を、逆投影利用率メモリ１１９に記憶された逆投影利用率ＢＵ（Xr，Zr）で除算する。

次に、ステップＳ３０５において、加算処理部１１５は、ステップＳ３０４において逆投影利用率ＢＵ（Xr，Zr）で除算された投影画像ベクトルの差分ｆ−Ｈｇ^k／ＦＵの逆投影を断層画像ベクトルｇ^kに加算する。
ここまでの結果、更新される断層画像ベクトルｇ^k+1は、以下の式（９）で表される。

次に、ステップＳ３０６において、画像処理部１０９は、断層画像ベクトルｇ^k+1、ｇ^kを比較した結果から、所定の収束条件を満足し、反復計算を終了するか否かを判定する。本実施形態では、ｇ^k+1−ｇ^kの絶対値が所定の値より小さくなるという収束条件、又は、ｋが所定の値より大きくなるという収束条件を満足する場合に、反復計算を終了すると判定する。反復計算を終了すると判定されたときの断層画像ベクトルｇ^k+1に基づく断層画像は、例えば、操作部１１１の操作に基づいて、モニタ１１０に表示されたり、画像保存部１１２に保存されたりする。式（９）で表される逐次近似再構成法は一般的にＳＩＲＴと呼ばれる最小二乗型の最も簡単な逐次近似再構成手法である。ＳＩＲＴには、加速係数（式（９）の係数「２」が加速係数に対応する）を最適化しながら収束を加速させる最急勾配法や共役勾配法もあるが、これらの手法により逐次近似再構成法を実行することができる。

以上のように本実施形態では、予め設定された照射角度βの全てにおけるＸ線管１０１の焦点と断層画像２０１の座標（所定の画素）とを相互に結ぶ直線のうち、Ｘ線検出器１０６と交わる直線の割合を逆投影利用率ＢＵとして計算する。また、断層画像２０１の総数のうち、Ｘ線管１０１の焦点とＸ線検出器１０６上の座標（再投影画像の所定の画素）とを相互に結ぶ直線が交わる断層画像２０１の数の割合を再投影利用率ＦＵとして計算する。ここで、逆投影利用率計算ＢＵと再投影利用率ＦＵは、寝台１０３の幅方向（Ｙ軸方向）に依存しない値になる。よって、再投影と逆投影における投影画像と断層画像の規格化用に使用されるメモリのメモリ容量を著しく削減することができる。これにより、ＣＴ装置やトモシンセシスにおける再構成演算装置の演算メモリを大幅に削減しコストを削減することが可能である。また、演算メモリが制限されている市販のＧＰＵによる逐次近似再構成も実現することが可能になる。以上のように、逆投影と再投影とを逐次的に繰り返して断層画像の再構成を行う逐次近似再構成の解が収束しないことの抑制と、断層画像を再構成する際に必要なメモリ容量の削減とを共に実現することができる。

本実施形態では、逆投影利用率計算ＢＵと再投影利用率ＦＵを両方用いる場合を例に挙げて説明した。しかしながら、両方の利用率を必ずしも用いる必要はなく、発明の利用者の実施形態に応じて片方のみを用いても良い。。また、本実施形態では、逐次近似再構成を例に挙げて説明したが、反復計算を行わない解析的な再構成においても本発明を用いることができる。この場合例えば、断層画像ベクトルｇは、式（９）の替わりに、ＦＢＰ（Filtered Back Projection）法に基づく以下の式（１０）のようになる。

式（１０）において、ｆは投影画像ベクトルであり、ＢＵは逆投影利用率であり、Ｈ^Tは、逆投影操作を表す行列（投影行列Ｈの転置行列）である。また、Ｒはフィルタ操作を表す行列（再構成フィルタ）であり、投影画像ベクトルｆに再構成フィルタを施すためのものである。この再構成フィルタＲとしては、例えば、一般に知られているRampフィルタやShepp＆Loganフィルタを用いることができる。

式（５）において説明したように、本実施形態では、寝台１０３の面が水平方向であり、投影画像２０３のＺ軸の座標Zdが常に一定である場合を例に挙げて説明した。しかしながら、必ずしもこのようにする必要はない。
また、Ｘ線管１０１が、寝台１０３の長手方向（Ｘ軸方向）と寝台１０３の面に対する法線方向（Ｚ軸方向）とにより定まる平面上を回動する場合を例に挙げて説明した。しかしながら、必ずしもこのようにする必要はない。例えば、寝台１０３の長手方向以外の寝台１０３の面方向の一つの方向と、寝台１０３の面に対する法線方向とにより定まる平面上を回動させることもできる。ここで、寝台１０３の面方向の一つの方向（例えばＸ軸方向）と、寝台１０３の面方向の一つの方向及び寝台１０３の面に対する法線方向に垂直な方向（例えばＹ軸方向）が、Ｘ線検出器１０６に配置される画素の行方向又は列方向になるようにするのが好ましい。
また、Ｘ線管１０１の焦点が鉛直下方を向くときに照射角度βが０度になるようにした。しかしながら、照射角度βの基準となる所定の方向は、鉛直下方に限定されるものではない。

尚、前述した実施形態は、何れも本発明を実施するにあたっての具体化の例を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。すなわち、本発明はその技術思想、又はその主要な特徴から逸脱することなく、様々な形で実施することができる。例えば、複数の機器から構成されるシステムに適用してもよいし、また、一つの機器からなる装置に適用してもよい。

（その他の実施例）
本発明は、以下の処理を実行することによっても実現される。即ち、まず、以上の実施形態の機能を実現するソフトウェア（コンピュータプログラム）を、ネットワーク又は各種記憶媒体を介してシステム或いは装置に供給する。そして、そのシステム或いは装置のコンピュータ（又はＣＰＵやＭＰＵ等）が当該コンピュータプログラムを読み出して実行する。

１０１ Ｘ線管、１０３ 寝台、１０９ 画像処理部

Claims (11)

1. 寝台の面方向の一つの方向と前記寝台の面に対する法線方向とにより定まる平面上を回動する放射線発生源から複数の照射角度で照射された放射線に基づいて平面検出器により検出された投影画像を取得する取得手段と、
   前記取得手段により取得された投影画像に基づいて逆投影を行って断層画像を生成する逆投影手段と、
   前記投影画像のうち、前記逆投影手段により前記断層画像が生成された際に用いた前記投影画像の割合を表す逆投影利用率であって、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に値が依存しない逆投影利用率を、前記断層画像の所定の画素のそれぞれについて計算する逆投影利用率計算手段と、
   前記逆投影手段により生成された断層画像の画素値を前記逆投影利用率で規格化する逆投影利用率演算手段と、
   を有し、
   前記逆投影利用率計算手段は、前記断層画像の画素値のうち、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に沿う位置にある複数の画素の画素値については１つの前記逆投影利用率を計算することを特徴とする画像処理装置。
2. 前記逆投影利用率計算手段は、前記複数の照射角度における前記放射線発生源の焦点と前記断層画像の所定の画素とを相互に結ぶ直線のうち、前記平面検出器と交わる直線の割合を表す逆投影利用率であって、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に値が依存しない逆投影利用率を計算することを特徴とする請求項１に記載の画像処理装置。
3. 前記断層画像に基づいて再投影を行って再投影画像を生成する再投影手段と、
   前記断層画像のうち、前記再投影手段により前記再投影画像が生成された際に用いた前記断層画像の割合を表す再投影利用率であって、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に値が依存しない再投影利用率を、前記再投影画像の所定の画素のそれぞれについて計算する再投影利用率計算手段と、
   前記再投影手段により生成された再投影画像の画素値を前記再投影利用率で規格化する再投影利用率演算手段と、
   前記取得手段により取得された投影画像と、前記再投影利用率演算手段により規格化された再投影画像との差分画像に基づく収束条件を含む所定の収束条件を満足するか否かを判定する判定手段と、
   を更に有し、
   前記判定手段により収束条件を満足しないと判定されると、前記再投影利用率演算手段により画素値が規格化された再投影画像が最新の投影画像として更新され、
   前記再投影利用率計算手段は、前記再投影画像の画素値のうち、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に沿う位置にある複数の画素の画素値については１つの前記再投影利用率を計算することを特徴とする請求項１又は２に記載の画像処理装置。
4. 前記再投影利用率計算手段は、前記断層画像の総数のうち、前記放射線発生源の焦点と前記平面検出器で得られた前記再投影画像の所定の画素とを相互に結ぶ直線が交わる断層画像の数の割合を表す再投影利用率であって、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に値が依存しない再投影利用率を計算することを特徴とする請求項３に記載の画像処理装置。
5. 前記取得手段により取得された投影画像と、前記再投影利用率演算手段により規格化された再投影画像との差分画像を生成する差分処理手段と、
   前記逆投影手段により生成された断層画像の画素値と、前記逆投影利用率演算手段により規格化された断層画像の画素値とを加算する加算手段と、
   を更に有し、
   前記判定手段は、前記加算手段により加算された断層画像の画素値と、前記逆投影手段により生成された断層画像の画素値とを比較した結果から、収束条件を満足するか否かを判定し、
   前記逆投影手段は、前記取得手段により取得された投影画像の替わりに、前記差分処理手段により生成された差分画像に基づいて逆投影を行って断層画像を生成し、
   前記再投影手段は、前記判定手段により収束条件を満足しないと判定されると、前記加算手段により画素値が加算された断層画像に基づいて再投影を行って再投影画像を生成することを特徴とする請求項４に記載の画像処理装置。
6. 寝台の面方向の一つの方向と前記寝台の面に対する法線方向とにより定まる平面上を回動する放射線発生源から複数の照射角度で照射された放射線に基づいて平面検出器により検出された投影画像を取得する取得工程と、
   前記取得工程により取得された投影画像に基づいて逆投影を行って断層画像を生成する逆投影工程と、
   前記投影画像のうち、前記逆投影工程により前記断層画像が生成された際に用いた前記投影画像の割合を表す逆投影利用率であって、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に値が依存しない逆投影利用率を、前記断層画像の所定の画素のそれぞれについて計算する逆投影利用率計算工程と、
   前記逆投影工程により生成された断層画像の画素値を前記逆投影利用率で規格化する逆投影利用率演算工程と、
   を有し、
   前記逆投影利用率計算工程は、前記断層画像の画素値のうち、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に沿う位置にある複数の画素の画素値については１つの前記逆投影利用率を計算することを特徴とする画像処理方法。
7. 前記逆投影利用率計算工程は、前記複数の照射角度における前記放射線発生源の焦点と前記断層画像の所定の画素とを相互に結ぶ直線のうち、前記平面検出器と交わる直線の割合を表す逆投影利用率であって、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に値が依存しない逆投影利用率を計算することを特徴とする請求項６に記載の画像処理方法。
8. 前記断層画像に基づいて再投影を行って再投影画像を生成する再投影工程と、
   前記断層画像のうち、前記再投影工程により前記再投影画像が生成された際に用いた前記断層画像の割合を表す再投影利用率であって、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に値が依存しない再投影利用率を、前記再投影画像の所定の画素のそれぞれについて計算する再投影利用率計算工程と、
   前記再投影工程により生成された再投影画像の画素値を前記再投影利用率で規格化する再投影利用率演算工程と、
   前記取得工程により取得された投影画像と、前記再投影利用率演算工程により規格化された再投影画像との差分画像に基づく収束条件を含む所定の収束条件を満足するか否かを判定する判定工程と、
   を更に有し、
   前記判定工程により収束条件を満足しないと判定されると、前記再投影利用率演算工程により画素値が規格化された再投影画像が最新の投影画像として更新され、
   前記再投影利用率計算工程は、前記再投影画像の画素値のうち、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に沿う位置にある複数の画素の画素値については１つの前記再投影利用率を計算することを特徴とする請求項６又は７に記載の画像処理方法。
9. 前記再投影利用率計算工程は、前記断層画像の総数のうち、前記複数の照射角度における前記放射線発生源の焦点と前記平面検出器で得られた前記再投影画像の所定の画素とを相互に結ぶ直線が交わる断層画像の数の割合を表す再投影利用率であって、前記寝台の面方向の一つの方向と前記寝台の面に対する法線方向とに垂直な方向に値が依存しない再投影利用率を計算することを特徴とする請求項８に記載の画像処理方法。
10. 前記取得工程により取得された投影画像と、前記再投影利用率演算工程により規格化された再投影画像との差分画像を生成する差分処理工程と、
    前記逆投影工程により生成された断層画像の画素値と、前記逆投影利用率演算工程により規格化された断層画像の画素値とを加算する加算工程と、
    を更に有し、
    前記判定工程は、前記加算工程により加算された断層画像の画素値と、前記逆投影工程により生成された断層画像の画素値とを比較した結果から、収束条件を満足するか否かを判定し、
    前記逆投影工程は、前記取得工程により取得された投影画像の替わりに、前記差分処理工程により生成された差分画像に基づいて逆投影を行って断層画像を生成し、
    前記再投影工程は、前記判定工程により収束条件を満足しないと判定されると、前記加算工程により画素値が加算された断層画像に基づいて再投影を行って再投影画像を生成することを特徴とする請求項９に記載の画像処理方法。
11. 請求項６〜１０の何れか１項に記載の画像処理方法の各工程をコンピュータに実行させることを特徴とするコンピュータプログラム。

Images