JP5458413B2 - 画像処理装置、画像処理方法、及び、画像処理プログラム - Google Patents

Description

本発明は、２次元又は３次元の複数枚の画像の位置を合わせる画像処理装置に関する。

２次元又は３次元の複数枚の画像の位置合わせ技術は、様々な分野で重要な技術である。

例えば、医用画像の分野では、ＣＴ画像、ＭＲ画像、ＰＥＴ画像、超音波画像など、様々な種類の３次元画像の位置を合わせて、重ね合わせて表示するために、画像の位置合わせ技術が利用される。このような表示方法はフュージョン画像表示と呼ばれており、画像の特徴を生かした表示が可能となる。たとえば、ＣＴ画像は詳細な形状の表示に適しており、ＰＥＴ画像は代謝や血流などの体の機能の表示に適している（例えば、特許文献１参照。）。

さらに、同一患者の経過を観察するためには、時系列に取得した医用画像の位置を合わせることによって、病変を容易に観察することが可能である。

医用画像の分野では、３次元画像だけではなく、２次元画像の位置合わせも重要である。例えば、エネルギーを変更した２種類の胸部放射線画像を用いることによって、骨部を除去した軟部画像を生成することができ、軟部を詳細に観察することが可能となる。この場合にも、２つの画像は撮像時刻が異なるため、位置合わせが必要となる（例えば、特許文献２参照。）。

一方、リモートセンシングの分野では、様々なセンサ及び様々な方法で取得された画像データを総合的に検証することによって、多くの情報が取得できる。

２次元又は３次元の複数枚の画像の位置合わせ、すなわち、レジストレーションは、前述したように幅広い分野で応用されている。レジストレーションでは、固定される画像を参照画像と呼び、位置合わせのために座標変換される画像を浮動画像と呼ぶ。レジストレーションでは、浮動画像に対して座標変換を行い、変換された浮動画像と参照画像との類似度を求め、類似度が最大となる座標変換を繰り返し計算によって求める。レジストレーションでは、二つの画像が正しく位置合わせができたかの指標として、差の二乗や相互相関、相互情報量（Mutual Information）などが用いられている。これらの指標のなかでも、特に、相互情報量は、対象となる画像の種類が異なっていても良好なレジストレーションが可能となるという特徴がある（例えば、非特許文献１参照）。

特開２００９−１１２４６８号公報 特開２００９−１９５４７１号公報 特開２００８−１４２１３７号公報

F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, and P. Suetens, "Multi-modality image registration maximization of mutual information", In Proceedings of MMBIA, pp. 14-22, 1996. R. Shams and N. Barnes: "Speeding up mutual information computation using NVIDIA CUDA hardware", Digital Image Computing Techniques and Applications, pp. 555-560,2007.

前述した従来の相互情報量を用いたレジストレーションによると、計算量が多く、処理時間が長いという問題がある。このため、相互情報量を用いたレジストレーションの高速化が重要な課題となっている。高速化のために、画素データのメモリアクセスが連続となる方法が提案されている（例えば、特許文献３参照。）。

一方、近年のマルチプロセッサやＧＰＵ（グラフィックプロセッシングユニット）の進歩に伴い、並列処理によるレジストレーションを高速化の期待が高まっている。また、ＧＰＵを用いた高速化方法も提案されている（例えば、非特許文献２参照。）。しかし、並列処理を行うとしても、相互情報量を計算するために必要な結合ヒストグラムの算出の計算量が多いという問題がある。

すなわち、相互情報量を計算するためには、位置合わせの対象となる二つの画像の結合ヒストグラムを計算する必要がある。結合ヒストグラムは２次元のヒストグラムであり、ヒストグラムに含まれるビン数が大きくなる。例えば、一つの画像のヒストグラムのビン数を６４とすると、結合ヒストグラム全体のビン数は６４×６４＝４０９６となる。

一般に、ＧＰＵは、高速にアクセスできるローカルメモリの容量が小さい。ＧＰＵを用いて結合ヒストグラムを計算する場合、このような大きなビン数のヒストグラムをローカルメモリに格納することができず、アクセスが低速なグローバルメモリにデータを格納する必要がある。このため、結合ヒストグラムを計算する処理に時間がかかるという問題がある。この問題は、高速にアクセスできるローカルメモリが少ない並列計算機においても同様である。

本発明は、画像レジストレーションの処理速度を向上する画像処理装置を提供することを目的とする。

本発明の代表的な一例を示せば以下の通りである。すなわち、２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせる画像処理システムであって、前記画像処理システムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、前記画像処理システムは、前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手段と、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２の手段と、前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手段と、前記生成された結合ヒストグラムから相互情報量を計算する第４の手段と、を備える。

本発明の実施の形態によると、各プロセッシングユニットに保持されるヒストグラムのビン数を削減でき、使用されるローカルメモリの容量を削減することができる。

本発明の実施の形態の画像処理装置の論理的な構成を示す図である。 本発明の実施の形態の画像処理装置のハードウェア構成を示すブロック図である。 本発明の実施の形態の画像処理装置の構成の一例を示すブロック図である。 本発明の実施の形態の処理を示すフローチャートである。 レジストレーションの対象となる参照画像を説明する図である。 レジストレーションの対象となる浮動画像を説明する図である。 レジストレーション処理の概要を説明する図である。 レジストレーション処理の概要を説明する図である。 従来技術によって結合ヒストグラムを計算する一般的な方法を説明する図である。 本発明によって結合ヒストグラムを計算する方法を説明する図である。 本発明の実施の形態の参照画像ヒストグラム作成処理の手順の詳細を示すフローチャートである。 本実施の形態の座標リスト群の詳細を説明する図である。 本発明の実施の形態の相互情報量計算処理の手順の詳細を示すフローチャートである。 本発明の実施の形態の浮動画像のヒストグラムを説明する図である。 本発明の実施の形態の結合ヒストグラムを説明する図である。 本発明の実施の形態の画像処理装置で実行される処理を示すシーケンス図である。 本発明の実施の形態の結合ヒストグラムの具体例を説明する図である。 本発明の実施の形態の結合ヒストグラムの具体例を説明する図である。 本発明の実施の形態の画像処理装置の構成の一例を示すブロック図である。

＜構成＞
図１は、本発明の実施の形態の画像処理装置の論理的な構成を示す図である。

本発明の実施の形態の画像処理装置は、前処理部１６及びレジストレーション部１７を備える。本実施の形態の画像処理装置は、レジストレーションの対象となる参照画像１及び浮動画像２が入力されると、レジストレーション済の浮動画像１５を生成して、生成された画像を出力する。

前処理部１６に参照画像１が入力されると、ヒストグラム生成部４は１次元のヒストグラムを生成し、参照画像ヒストグラムの各ビンの座標リストである座標リスト群５を生成する。

レジストレーション部１７に前記座標リスト群５及び浮動画像２が入力されると、相互情報量計算部１８は、幾何変換情報３を参照して、相互情報量１２を求める。なお、相互情報量計算部１８の処理内容の詳細は後述する。相互情報量最大化部１３は、求められた相互情報量１２の幾何変換情報３を更新し、相互情報量計算部１８が、相互情報量１２を再計算し、相互情報量１２を最大化する。浮動画像幾何変換部１４は、相互情報量が最大となる幾何変換情報３を用いて浮動画像２を幾何変換することによって、レジストレーション済の浮動画像１５を生成し、生成された画像を出力する。以上の処理によってレジストレーションが実行される。

次に、相互情報量計算部１８によって実行される処理の概要について説明する。相互情報量計算部１８は、幾何変換情報３を参照して、座標リスト群５の各画素を座標変換し、変換後の座標を求める。浮動画像参照・量子化部７は、変換後の座標に基づいて浮動画像２を参照し、輝度値を量子化し、参照画像ヒストグラムのビンごとに浮動画像のヒストグラム８を生成する。結合ヒストグラム生成部９は、参照画像ヒストグラムのビンごとの浮動画像のヒストグラム８から結合ヒストグラム１０を求める。相互情報量計算部１１は、結合ヒストグラム１０を用いて相互情報量１２を求める。

図２は、本発明の実施の形態の画像処理装置のハードウェア構成を示すブロック図である。

本実施の形態の画像処理装置は、ＧＰＵ及び／又は並列プロセッサシステムを備える一般的な計算機上に実装可能である。画像処理装置は、ＣＰＵ４０、ＲＯＭ２１、ＲＡＭ２２、記憶装置２３、ＧＰＵ４１、画像入力部３０、媒体入力部２６、入力制御部２８及び画像生成部２９を備える。ＣＰＵ４０、ＲＯＭ２１、ＲＡＭ２２、ハードディスク２３、ＧＰＵ４１、画像入力部３０、媒体入力部２６、入力制御部２８及び画像生成部２９は、データバス２４によって相互に接続される。

ＣＰＵ４０は、ＲＯＭ２１又はＲＡＭ２２に記憶されたプログラムを実行することによって、各種処理を実行する。

ＲＯＭ２１及びＲＡＭ２２は、各種処理を実行するために必要なプログラム及びデータを記憶する。ＲＯＭ２１は、データの読み出し専用の記憶媒体である。ＲＡＭ２２は、データの読み書きが可能な記憶媒体である。

記憶装置２３は、入力画像などを格納する磁気記憶装置である。記憶装置２３は、不揮発性半導体記憶媒体（例えば、フラッシュメモリ）を備えてもよい。また、ネットワークなどを介して接続された外部記憶装置を利用してもよい。

ＧＰＵ４１は、画像データの処理に適する集積回路であり、その構成は図３を用いて詳述する。なお、ＧＰＵ４１は並列プロセッサシステムでもよい。

本実施の形態の画像処理装置において実行されるプログラムは、記憶媒体２５（例えば、光ディスク）に格納されており、媒体入力部２６（例えば、光ディスクドライブ）によって読み込まれ、読み込まれたプログラムがＲＡＭ２２に格納される。なお、記憶装置２３に当該プログラムを格納し、記憶装置２３からプログラムをＲＡＭ２２にロードしてもよい。また、ＲＯＭ２１にあらかじめプログラムを記憶させてもよい。

画像処理装置においてプログラムが実行されることによって、図１に示す前処理部１６及びレジストレーション部１７が機能する。ＲＡＭ２２は、図１に示す座標リスト群５、幾何変換情報３、浮動画像のヒストグラム８、結合ヒストグラム１０、相互情報量１２及びレジストレーション済浮動画像１５を格納する。

なお、ＧＰＵ４１がレジストレーション処理を実行する場合、ＧＰＵ４１上のグローバルメモリが、座標リスト群５、幾何変換情報３、浮動画像のヒストグラム８、結合ヒストグラム１０、相互情報量１２及びレジストレーション済浮動画像１５を格納してもよい。

画像入力部３０は、画像撮像装置２０が接続され、画像撮像装置２０によって撮影された画像が入力されるインターフェースである。画像撮像装置２０から入力された画像は、ＣＰＵ４０によって各種処理がされる。

媒体入力部２６は、記憶媒体２５に記憶されたデータ及びプログラムを読み出す。記憶媒体２５から読み出されたデータ及びプログラムは、ＣＰＵ４０によって、ＲＡＭ２２又はハードディスク２３に格納される。

入力制御部２８は、入力装置２７（例えば、キーボード）から、ユーザによって入力されたデータを受け付けるインターフェースである。入力制御部２８が受けたデータは、ＣＰＵ４０によって処理される。

画像生成部２９は、図１に示すレジストレーション部１７によってレジストレーション処理がされた浮動画像１５から画像データを生成し、生成された画像データをディスプレイ１４に送る。

図３は、ＧＰＵ４１の構成の一例を示すブロック図である。

ＧＰＵ４１は、グローバルメモリ４０３及び複数のマルチプロセッシングユニット４０５を備える。各マルチプロセッシングユニット４０５は、ローカル共有メモリ４０２及び複数のプロセッシングユニット４０１を備える。プロセッシングユニット４０１は、同じマルチプロセッシングユニット４０５内のローカル共有メモリ４０２を読み書き可能である。しかし、プロセッシングユニット４０１は、別のマルチプロセッシングユニット４０５のローカル共有メモリ４０２にはアクセスできない。各プロセッシングユニット４０１は、グローバルメモリ４０３を読み書き可能であるので、グローバルメモリ４０３には、複数のプロセッシングユニット４０１で共有されるデータが格納される。

なお、一般的にＧＰＵのローカル共有メモリ４０２は高速に読み書き可能であるが、グローバルメモリ４０３へのアクセスは、ローカル共有メモリ４０２のアクセスと比較すると低速である。このため、ＧＰＵを用いて効率よく処理を実行するためには、グローバルメモリ４０３のアクセスを削減し、ローカル共有メモリ４０２を用いて処理を実行することが重要である。ただし、ローカル共有メモリ４０２は、一般にその容量に制限があるので、ローカル共有メモリの容量内で効率のよいアルゴリズムを用いることがＧＰＵを用いた処理の高速化において重要である。

図１５は、図２の構成において、ＧＰＵ４１の代わりにマルチプロセッサシステムを用いる場合のマルチプロセッサシステムの構成の一例を示すブロック図である。

マルチプロセッサシステムは、グローバルメモリ１５０５及び複数のマルチプロセッシングエレメント１５０３を備える。各マルチプロセッシングエレメント１５０３は、プロセッシングユニット１５０１とローカルメモリ１５０２を備える。ローカルメモリはバス１５０４を介してグローバルメモリ１５０５と結合される。またバス１５０４はホストインターフェース１５０６を介して、図２のバス２４と接続される。
図１５に示すマルチプロセッサシステムを用いる場合には，ローカルメモリ１５０２は高速に読み書き可能であるが、グローバルメモリ１５０５へのアクセスは、ローカルメモリ１５０２のアクセスと比較すると低速である。このため、マルチプロセッサシステムを用いて効率よく処理を実行するためには、グローバルメモリ１５０５のアクセスを削減し、ローカルメモリ１５０２を用いて処理を実行することが重要である。ただし、ローカルメモリ１５０２は、一般にその容量に制限があるので、ローカルメモリの容量内１５０２で効率のよいアルゴリズムを用いることがＧＰＵを用いた処理の高速化において重要である。

本願発明は、以上の点を考慮した効率的なレジストレーション方法を提案するものである。

＜動作＞
次に、本発明の実施の形態の動作について説明する。

図４は、本発明の実施の形態の処理を示すフローチャートである。

図４に示す処理のうち、ステップＳ１０１及びＳ１０２の処理が、前処理部１６によって実行される処理であり、ステップＳ１０３〜Ｓ１０８の処理が、レジストレーション部１７によって実行される処理である。

まず、画像入力部３０が、参照画像１の入力を受け付ける（Ｓ１０１）。

次に、ヒストグラム生成部４が、入力された参照画像１から１次元のヒストグラムを作成する（Ｓ１０２）。さらに、ヒストグラムの各ビンに属する参照画像の画素の座標を座標リスト群５に格納する。なお、Ｓ１０２の処理内容の詳細は、図８を用いて後述する。

次に、画像入力部３０が、浮動画像２の入力を受け付ける（Ｓ１０３）。その後、幾何変換情報３の初期値を設定する（Ｓ１０４）。レジストレーションでは、二つの画像の相互情報量が最大となる幾何変換を求めるために収束計算を行う。このため、幾何変換情報３の初期値が必要となる。幾何変換情報の初期値は、ユーザが与えてもよく、まったく幾何変換しない、すなわち恒等変換を初期値として与えてもよい。

次に、相互情報量計算部１８が、相互情報量１２を求め、メモリに格納する（Ｓ１０５）。その後、求められた相互情報量が収束したかどうかを判定する（Ｓ１０６）。判定の結果、相互情報量が収束していない場合、より高い相互情報量を得るために、幾何変換情報３を更新する（Ｓ１０７）。そして、更新された幾何変換情報３を用いて、相互情報量１２を再計算し（Ｓ１０５）、再度収束を判定する（Ｓ１０６）。

一方、相互情報量が収束している場合、求められた幾何変換情報３を用いて、浮動画像２を幾何変換し、レジストレーションされた浮動画像１５にＲＡＭ２２に格納する（Ｓ１０８）。以上の処理によって、画像のレジストレーションが完了する。

図５Ａは、レジストレーションの対象となる参照画像１を説明する図であり、図５Ｂは、レジストレーションの対象となる浮動画像２を説明する図である。

図５Ａ、図５Ｂでは、簡単のため画像を２次元として図示しているが、参照画像１及び浮動画像２は３次元画像である。参照画像１は、図５Ａに示すように、Ｘ軸方向にＲｘ画素、Ｙ軸方向にＲｙ画素、図示は省略したがＺ軸方向にＲｚ画素、合計Ｒｘ×Ｒｙ×Ｒｚ個の画素３０１から構成される３次元画像である。同様に、浮動画像２は、Ｘ軸方向にＦｘ画素、Ｙ軸方向にＦｙ画素、図示は省略したがＺ軸方向にＦｚ画素、合計Ｆｘ×Ｆｙ×Ｆｚ個の画素３０２から構成される３次元画像である。

図６Ａ、図６Ｂは、レジストレーション処理の概要を説明する図である。

図６Ａに示すように、参照画像１と浮動画像２は類似した画像ではあるが、その位置又は方向が違っている。例えば、異なる視点から同一の患者のＣＴ画像、ＰＥＴ画像などを撮影した場合、類似した複数の画像が取得される。また、ＣＴ画像とＰＥＴ画像のように、輝度値の性質が異なる場合も、類似した複数の画像が取得される。

そして、図６Ｂに示すように、浮動画像２にレジストレーション処理を行うことによって、浮動画像２を幾何変換し、幾何変換された浮動画像１３と参照画像１とを位置合わせすることができる。

レジストレーション処理における幾何変換は、３次元の移動及び回転から構成される剛体変換を用いることが可能である。また、一般的なアフィン変換や、Ｂスプライン関数を用いた非剛体変換を適用することもできる。

剛体変換を用いた幾何変換には、式（１）を用いることができる。

式（１）において，（ｘ，ｙ，ｚ）は幾何変換前の座標、（ｘ’，ｙ’，ｚ’）は変換後の座標、Ｒｘ（θｘ）、Ｒｙ（θｙ）、Ｒｚ（θｚ）はそれぞれ，Ｘ軸周りの回転行列、Ｙ軸周りの回転行列、Ｚ軸周りの回転行列である。また、（ｔｘ，ｔｙ，ｔｚ）は移動ベクトルである、すなわち、３次元画像の剛体変換は６自由度を有する。なお、２次元画像の剛体変換では、ｘ方向への移動、ｙ方向への移動及び回転の３自由度を有する。

アフィン変換を用いた幾何変換には、式（２）を用いることができる。

式（２）において、（ｘ，ｙ，ｚ）は幾何変換前の座標、（ｘ’，ｙ’，ｚ’）は変換後の座標、ａ、ｂ、ｃ、ｄ、ｅ、ｆ、ｇ、ｈ、ｉはアフィン変換のパラメータ、（ｔｘ，ｔｙ，ｔｚ）は移動ベクトルである。

Ｂスプライン関数を用いた非剛体変換には、式（３）を用いることができる。

Ｂスプライン関数を用いた非剛体変換では、変換パラメータはＢスプライン立体の格子点位置となる。式（３）において、φｉ，ｊ，ｋは、ｎｘ×ｎｙ×ｎｚ個の制御点である。

ここで、本発明との対比のため従来技術における結合ヒストグラムの計算について説明する。図７Ａは、従来技術によって結合ヒストグラムを計算する一般的な方法を説明する図である。

従来の結合ヒストグラムの計算方法では、図７Ａに示すように、参照画像１の画素値と、幾何変換された浮動画像の対応する位置における画素値とを求める。画素値には、例えば、各画素の輝度値を用いることができる。なお、特定波長の輝度値を用いてもよい。

そして、この求められた二つの輝度値に基づいて、結合ヒストグラムを求める。具体的には、参照画像のヒストグラムのビン数Ｎｒを参照して、参照画像の輝度値が属する参照画像ヒストグラムのビンｉを求める。同様に、幾何変換された浮動画像のヒストグラムのビン数Ｎｆを参照して、浮動画像の輝度値が属する浮動画像のヒストグラムのビンｊを求める。求められた（ｉ，ｊ）から対応する結合ヒストグラムのビンの頻度を更新する。

以上説明したように、従来技術によると、参照画像１と幾何変換された浮動画像１３とに基づいて結合ヒストグラム１０を求めていた。レジストレーション処理における画像の類似度の指標である相互情報量を求めるためには、前述した結合ヒストグラムの算出が必要である。しかし、結合ヒストグラムは２次元のヒストグラムであり、その容量が大きい性質があることから、結合ヒストグラムを記憶するために大きなメモリ容量が必要である。相互情報量の算出では、ヒストグラムのビン数は６４〜２５６程度である。ビン数が６４の場合では、結合ヒストグラムの総ビン数は６４×６４＝４０９６ビンとなる。

一般に、ＧＰＵ及び並列プロセッサでは、各プロセッサユニットが局所的に利用できるメモリ容量に制限がある場合が多い。これは、例えば、図３に示すＧＰＵ装置において、ローカル共有メモリ４０２の容量が制限されているためである。このため、共有メモリの容量の不足が、ＧＰＵや並列プロセッサで結合ヒストグラムを計算する場合の問題点となっていた。

一方、本実施の形態では、図７Ｂに示すように、参照画像ヒストグラムの座標リスト５を作成し、その後、参照画像ヒストグラムのビンごとに浮動画像の１次元ヒストグラム８を生成して、結合ヒストグラムを生成する。このように、各プロセッサユニットが２次元の結合ヒストグラムを保持する必要がなく、１次元のヒストグラムを保持するだけでよい計算方法を提供することによって、前述した問題を解決することができる。

相互情報量の計算には、まず、参照画像１と浮動画像２の同じ位置の画素値から生成される２次元ヒストグラムを計算する。その後、計算された２次元ヒストグラムを用いて計算される確率密度関数を用いて、相互情報量を計算することができる。

相互情報量とは、２つの事象ＡとＢについて、事象Ａが有している事象Ｂに関する情報量を定量化した尺度である。レジストレーションで用いられる正規化相互情報量ＮＭＩ（Ａ，Ｂ）は事象Ａと事象Ｂの２次元の結合ヒストグラムＨｉｓｔ（Ａ，Ｂ）から下式（４）によって求められる。

ここで、Ｈ（Ａ）は事象Ａのエントロピー、Ｈ（Ｂ）は事象Ｂのエントロピー、Ｈ（Ａ，Ｂ）は事象Ａ，Ｂの結合エントロピーである。ｐ（ａ）はａの確率密度分布、ｐ（ｂ）はｂの確率密度分布、ｐ（ａ，ｂ）はａ，ｂの同時確率分布であり、式（４）のように、Ｈｉｓｔ（Ａ，Ｂ）から求められる。

事象Ａと事象Ｂが完全に独立である場合には、Ｈ（Ａ，Ｂ）＝Ｈ（Ａ）＋Ｈ（Ｂ）となる。また、事象Ａと事象Ｂが完全に従属である場合には、Ｈ（Ａ，Ｂ）＝Ｈ（Ａ）＝Ｈ（Ｂ）となる。以上から，ＮＭＩ（Ａ，Ｂ）のとりうる範囲は、１以上２以下となる。

図８は、図４のフローチャートにおける参照画像ヒストグラム作成（Ｓ１０２）の処理手順の詳細を示すフローチャートである。

まず、変数Ｚを１から参照画像のＺ軸方向の画素数Ｒｚまで増加させながら、ステップＳ２０２からＳ２０８の処理を繰り返し実行する（Ｓ２０１）。次に、変数Ｙを１から参照画像のＹ軸方向の画素数Ｒｙまで増加させながら、ステップＳ２０３からＳ２０７の処理を繰り返し実行する（Ｓ２０２）。さらに、変数Ｘを１から参照画像のＸ軸方向の画素数Ｒｘまで増加させながら、ステップＳ２０４からＳ２０６の処理を繰り返し実行する（Ｓ２０３）。以上の処理によって、参照画像のすべての画素（Ｘ，Ｙ，Ｚ）について、ステップＳ２０４からＳ２０６の処理が実行される。

ステップＳ２０４では、参照画像の画素値Ｒ（Ｘ，Ｙ，Ｚ）を変数Ｉに代入する。次に、画素値Ｉを量子化して、量子化された画素値Ｉを用いて参照画像のヒストグラムのビンを求める（Ｓ２０５）。この処理は、例えば、参照画像の輝度値の最大値と最小値との間をＮｒ個の領域に分割し、輝度値が属する領域の番号からビンを求める。

その後、ヒストグラムのビンと、座標値（Ｘ，Ｙ，Ｚ）とを用いて、座標リスト群５にデータを格納する（Ｓ２０６）。座標リスト群５へのデータの格納方法の詳細は、図９を用いて後述する。

なお、前述した説明では、参照画像の全画素について、ヒストグラム及びヒストグラムリストを求めたが、ダウンサンプリングを行って少ないサンプル数の画素についてヒストグラム及びヒストグラムリストを求めてもよい。また、スーパーサンプリングを行って、サンプル数を増やしてもよい。さらに、等間隔にサンプリングするのではなく、ランダムにサンプリングしてもよい。

図９は、本実施の形態の座標リスト群５の詳細を説明する図である。

ヒストグラムリスト５は、Ｎｒ個の可変長のリストを含む。各リストの構成要素は、画素の座標である。

図８のステップＳ２０６において、ビン番号ｂに座標（Ｘ，Ｙ，Ｚ）を登録する場合、ｂ番目のリストの最後尾に座標（Ｘ，Ｙ，Ｚ）を追加する。

図９に示すように、参照画像１のすべての画素に対して、座標（Ｘ，Ｙ，Ｚ）を追加した場合の、ｂ番目のリストの長さをＢ（ｂ）とする。Ｂ（ｂ）は参照画像の１次元のヒストグラムと等価である。

図１０は、図４のフローチャートにおける相互情報量計算（Ｓ１０５）の処理手順の詳細を示すフローチャートである。

まず、変数ｂを１から参照画像のヒストグラムのビン数Ｎｒまで増加させながら、ステップＳ３０２からステップＳ３０８を繰り返し実行する（Ｓ３０１）。

次に、変数ｉを１からヒストグラムの長さＢ（ｂ）まで増加させながら、ステップＳ３０３からステップＳ３０７を繰り返し実行する（Ｓ３０２）。

ループ内では、座標リスト群５を参照して、対応する座標（Ｘ，Ｙ，Ｚ）を取り出し（Ｓ３０３）、幾何変換情報３を用いて座標（Ｘ，Ｙ，Ｚ）を座標変換し、対応する浮動画像の座標系に変換する（Ｓ３０４）。

次に、浮動画像の対応する画素値Ｒ（Ｘ，Ｙ，Ｚ）を変数Ｉに代入し（Ｓ３０５）、画素値Ｉから浮動画像のヒストグラムのビンを求める。この処理は、浮動画像の最大値と最小値とをＮｆ個の領域に分割し、輝度値が属する領域の番号からビンを求めることで実現できる（Ｓ３０６）。

その後、量子化された値に対応するビン番号ｂ’の浮動画像ヒストグラム８を更新する（Ｓ３０７）。処理の詳細は後述する。

ステップＳ３０１からステップＳ３０９のループによって、参照画像の全画素に対応する参照画像ヒストグラムのビンごとの浮動画像のヒストグラム８が求められる。ヒストグラム８は、図１１に示すように、ビン数Ｎｆの１次元ヒストグラムをＮｒ個含む。

その後、参照画像ヒストグラムのビンごとの浮動画像のヒストグラム８を合成して、２次元の結合ヒストグラム１０を求める（Ｓ３１０）。結合ヒストグラム１０は、図１２に示すように、ヒストグラム８に含まれる一次元ヒストグラムを結合したもので、Ｎｒ×Ｎｆのマトリックスで構成されている。

結合ヒストグラム１０の具体例を図１４Ａ、図１４Ｂに示す。図１４Ａに示すように、参照画像１と浮動画像２との位置を合わせない場合、結合ヒストグラム１０には複数のピークが現れる。しかし、図１４Ｂに示すように、参照画像１と浮動画像２との位置が合っている状態では、結合ヒストグラム１０には一つのピークが現れる。

その後、結合ヒストグラム１０から相互情報量１２を求める（Ｓ３１１）。

次に、図１０のステップＳ３０７の処理の詳細について説明する。

ステップＳ３０７では、参照画像のビンがｂ、浮動画像のビンがｂ’である場合、ｂ番目のヒストグラムのｂ’のビンの頻度に１を加算する。すなわち、２次元の結合ヒストグラム１０ではなく、１次元のヒストグラム８が更新される。

このことから、ステップＳ３０２からのループであるステップＳ３０３からステップＳ３０８の処理をＧＰＵ又は並列プロセッサで実行する場合、各プロセッシングユニット４０１は２次元のヒストグラムをローカル共有メモリ４０２に保持する必要がなく、１次元のヒストグラムをローカル共有メモリ４０２に保持するだけでよい。

さらに、ステップＳ３０１からのループであるステップＳ３０２からステップ３０９の処理をＧＰＵ又は並列プロセッサで実行することも可能である。この場合にも、各プロセッシングユニット４０１は２次元のヒストグラムをローカル共有メモリ４０２に保持する必要がなく、１次元のヒストグラムをローカル共有メモリ４０２に保持するだけでよい。

３次元画像のレジストレーション実験によれば、本発明は非特許文献２の方法と比較して、より高速である。具体的な実験は、ファントムのＣＴ画像を用いてレジストレーションの性能を評価した。使用したファントム画像は５１２×５１２×２００ボクセルの３次元画像である。ファントム画像に対して、Ｘ軸、Ｙ軸及びＺ軸に対して１０度回転したボリュームデータを作成しレジストレーションを行った。レジストレーション実験によればもとの位置に正しく座標変換された。

レジストレーションに使用するサンプル数を種々に変更して計算時間を計測した。サンプル数が６４×６４×６４の場合の計算時間は、ＣＰＵ実装が１４．２秒、非特許文献２の方法が１．８２秒、本発明による方法が１．１６秒であった。サンプル数が２５６×２５６×２５６の場合にはＣＰＵ実装の計算時間は７８０秒、非特許文献２の方法が１４．９秒、本発明による方法が５．７秒であった。以上のように本発明による方法はＣＰＵ実装と比較して数１０倍以上、非特許文献２の方法と比較して１．５倍から２．６倍程度の高速化が実現されている。

図１３は、本実施の形態の画像処理装置で実行される処理を示すシーケンス図である。

まず、レジストレーション命令が発行されると、入力装置２７は、ユーザからの参照画像及び浮動画像の指定を受けて、画像の指定を入力制御部２８を介してＣＰＵ４０に伝える（１３０１）。ＣＰＵ４０は、指定された参照画像及び浮動画像を取得し、取得した参照画像及び浮動画像を記憶装置２３に格納する（１３０２）。ＧＰＵ４１は、記憶装置２３に格納された参照画像を取得し、参照画像ヒストグラムを生成し（１３０３）、生成された参照画像ヒストグラムを記憶装置に格納する。

さらに、入力装置２７は、ユーザからの幾何変換情報の指定を受けて、入力された幾何変換情報を入力制御部２８を介してＣＰＵ４０に伝える（１３０４）。この幾何変換情報は、幾何変換の初期条件を含み、さらに、幾何変換情報を更新する条件を含んでもよい。ＣＰＵ４０は、入力された幾何変換情報をＧＰＵ４１に転送する（１３０５）。

ＧＰＵ４１は、幾何変換を受けると、浮動画像を幾何変換し、幾何変換された浮動画像の画素値をサンプリングし、参照画像ヒストグラムのビンごとに浮動画像の１次元ヒストグラム８を生成する（１３０６）。このため、ＧＰＵ４１は、浮動画像の１次元ヒストグラム８をローカル共有メモリ４０２に格納すればよく。小容量で高速アクセス可能なメモリのみを用いて（すなわち、大容量だがアクセス速度が低いグローバルメモリ４０３を用いることなく）、浮動画像の１次元ヒストグラム８を生成することができる。その後、ＧＰＵ４１は、生成された浮動画像の１次元ヒストグラム８を、記憶装置２３に格納する。

ＧＰＵ４１は、全ての１次元ヒストグラムの計算を終了した後、記憶装置２３から１次元ヒストグラム８を読み出して、読み出された１次元ヒストグラム８を結合する。そして結合されたヒストグラム１０から相互情報量を計算する（１３０７）。計算された相互情報量はＣＰＵ４０に転送される（１３０８）。

ＣＰＵ４０は、相互情報量が最大となるように幾何変換情報を更新し（１３０９）、更新された幾何変換情報をＧＰＵに転送する（１３１０）。

ＧＰＵ４１は、転送された幾何変換情報に基づいて、１３０６〜１３０９の処理と同様に、幾何変換された浮動画像の１次元ヒストグラム８を生成し、結合ヒストグラムを生成し、相互情報量を計算し（１３１１）、計算された相互情報量をＣＰＵ４０に転送する（１３１２）。

ＣＰＵ４０は、相互情報量が最大となる幾何変換情報で浮動画像を変換し（１３１３）、この条件で幾何変換された浮動画像（レジストレーション済画像）をディスプレイ３１に表示する（１３１４）。

以上に説明したように、本発明の実施の形態では、２次元の結合ヒストグラムの計算をＧＰＵ又は並列プロセッサで実行する場合に、各プロセッシングユニットは１次元のヒストグラムにアクセスするので、必要なメモリ容量を削減することができる。このため、結合ヒストグラムを効率よく計算でき、レジストレーション計算を効率的に処理することができる。

１ 参照画像
２ 浮動画像
３ 幾何変換情報
５ 参照画像ヒストグラム（座標リスト群）
８ 浮動画像ヒストグラム
１０ 結合ヒストグラム
１２ 相互情報量
１５ レジストレーション済み浮動画像
１６ 前処理部
１７ レジストレーション部
４１ ＧＰＵ

Claims (12)

1. ２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせる画像処理システムであって、
   前記画像処理システムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、
   前記画像処理システムは、
   前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手段と、
   所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２の手段と、
   前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手段と、
   前記生成された結合ヒストグラムから相互情報量を計算する第４の手段と、を備えることを特徴とする画像処理システム。
2. 前記画像処理システムは、並列して演算処理が可能な複数の前記プロセッシングエレメントと、前記ローカルメモリと、前記グローバルメモリとを備えるメニーコアプロセッサシステムによって構成され、
   前記各プロセシングエレメントが所定のプログラムを実行することによって、前記ローカルメモリを用いて、参照画像のヒストグラムのビン毎に並列して、前記第２の手段を実現することを特徴とする請求項１に記載の画像処理システム。
3. 前記メニーコアプロセッサシステムは、グラフィックプロセッシングユニットであることを特徴とする請求項２に記載の画像処理システム。
4. 前記第２の手段は、前記生成された参照画像のヒストグラムの一つのビンに属する画素の位置を取得し、前記取得された画素の位置に対応する前記浮動画像の画素の画素値を取得し、前記取得された浮動画像の画素の画素値を量子化し、前記量子化された画素値に画素数を積算することによって、浮動画像の１次元ヒストグラムを生成する処理を、前記参照画像のヒストグラムのビン数だけ実行することを特徴とする請求項１に記載の画像処理システム。
5. 前記第１の手段は、前記参照画像の画素を表すパラメータを変化させ、当該パラメータによって特定される画素の画素値を量子化し、前記量子化された画素値のビンに前記特定された画素の位置を記録することによって、前記参照画像のヒストグラムを生成することを特徴とする請求項１に記載の画像処理システム。
6. 前記第２の手段は、
   前記生成された参照画像のヒストグラムの一つのビンに属する画素に対応する位置の前記変換された浮動画像の画素値を取得し、前記取得された浮動画像の画素値を量子化し、前記量子化された画素値を前記浮動画像のヒストグラムに加える処理を、前記参照画像のヒストグラムの一つのビンに属する画素の数だけ実行し、
   さらに、前記の処理を、前記生成された１次元ヒストグラムの一つのビンの数だけ実行することを特徴とする請求項１に記載の画像処理システム。
7. ２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせるためのパラメータを計算する画像処理方法であって、
   前記画像処理方法が実行されるシステムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、
   前記画像処理方法は、
   前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１のステップと、
   前記各プロセッシングエレメントが、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２のステップと、
   前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３のステップと、
   前記生成された結合ヒストグラムから相互情報量を計算する第４のステップと、を含むことを特徴とする画像処理方法。
8. 前記画像処理システムは、複数の前記プロセッシングエレメントと、前記ローカルメモリと、前記グローバルメモリとを備え、各プロセシングエレメントが並列して演算処理が可能なメニーコアプロセッサシステムによって構成され、
   前記画像処理システムは、並列して演算処理が可能な複数の前記プロセッシングエレメントと、前記ローカルメモリと、前記グローバルメモリとを備えるメニーコアプロセッサシステムによって構成され、
   前記各プロセシングエレメントは、前記ローカルメモリを用いて、参照画像のヒストグラムのビン毎に並列して、前記第２のステップの処理を実行することを特徴とする請求項７に記載の画像処理方法。
9. 前記第２ステップでは、前記生成された参照画像のヒストグラムの一つのビンに属する画素の位置を取得し、前記取得された画素の位置に対応する前記浮動画像の画素の画素値を取得し、前記取得された浮動画像の画素の画素値を量子化し、前記量子化された画素値に画素数を積算することによって、浮動画像の１次元ヒストグラムを生成する処理を、前記参照画像のヒストグラムのビン数だけ実行することを特徴とする請求項７に記載の画像処理方法。
10. 前記第１ステップでは、前記参照画像の画素を表すパラメータを変化させ、当該パラメータによって特定される画素の画素値を量子化し、前記量子化された画素値のビンに前記特定された画素の位置を記録することによって、前記参照画像のヒストグラムを生成することを特徴とする請求項７に記載の画像処理方法。
11. 前記第２ステップでは、
    前記生成された参照画像のヒストグラムの一つのビンに属する画素に対応する位置の前記変換された浮動画像の画素値を取得し、前記取得された浮動画像の画素値を量子化し、前記量子化された画素値を前記浮動画像のヒストグラムに加える処理を、前記参照画像のヒストグラムの一つのビンに属する画素の数だけ実行し、
    さらに、前記の処理を、前記生成された１次元ヒストグラムの一つのビンの数だけ実行することを特徴とする請求項７に記載の画像処理方法。
12. ２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせるためのパラメータを画像処理システムに計算させるための画像処理プログラムであって、
    前記画像処理方法が実行されるシステムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、
    前記プログラムは、
    前記各プロセッシングエレメントが、前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手順と、
    前記各プロセッシングエレメントが、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２の手順と、
    前記各プロセッシングエレメントが、前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手順と、
    前記各プロセッシングエレメントが、前記生成された結合ヒストグラムから相互情報量を計算する第４の手順と、を含むことを特徴とする画像処理プログラム。