画像処理装置、画像処理方法、及び、画像処理プログラム
【課題】画像レジストレーションの処理速度を向上する画像処理装置を提供する。
【解決手段】参照画像と浮動画像との位置を合わせる画像処理システムであって、前記画像処理システムは、前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手段と、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、ローカルメモリ上に生成する処理を実行する第２の手段と、前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手段と、を備える。
【解決手段】参照画像と浮動画像との位置を合わせる画像処理システムであって、前記画像処理システムは、前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手段と、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、ローカルメモリ上に生成する処理を実行する第２の手段と、前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手段と、を備える。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、２次元又は３次元の複数枚の画像の位置を合わせる画像処理装置に関する。
【背景技術】
【０００２】
２次元又は３次元の複数枚の画像の位置合わせ技術は、様々な分野で重要な技術である。
【０００３】
例えば、医用画像の分野では、ＣＴ画像、ＭＲ画像、ＰＥＴ画像、超音波画像など、様々な種類の３次元画像の位置を合わせて、重ね合わせて表示するために、画像の位置合わせ技術が利用される。このような表示方法はフュージョン画像表示と呼ばれており、画像の特徴を生かした表示が可能となる。たとえば、ＣＴ画像は詳細な形状の表示に適しており、ＰＥＴ画像は代謝や血流などの体の機能の表示に適している（例えば、特許文献１参照。）。
【０００４】
さらに、同一患者の経過を観察するためには、時系列に取得した医用画像の位置を合わせることによって、病変を容易に観察することが可能である。
【０００５】
医用画像の分野では、３次元画像だけではなく、２次元画像の位置合わせも重要である。例えば、エネルギーを変更した２種類の胸部放射線画像を用いることによって、骨部を除去した軟部画像を生成することができ、軟部を詳細に観察することが可能となる。この場合にも、２つの画像は撮像時刻が異なるため、位置合わせが必要となる（例えば、特許文献２参照。）。
【０００６】
一方、リモートセンシングの分野では、様々なセンサ及び様々な方法で取得された画像データを総合的に検証することによって、多くの情報が取得できる。
【０００７】
２次元又は３次元の複数枚の画像の位置合わせ、すなわち、レジストレーションは、前述したように幅広い分野で応用されている。レジストレーションでは、固定される画像を参照画像と呼び、位置合わせのために座標変換される画像を浮動画像と呼ぶ。レジストレーションでは、浮動画像に対して座標変換を行い、変換された浮動画像と参照画像との類似度を求め、類似度が最大となる座標変換を繰り返し計算によって求める。レジストレーションでは、二つの画像が正しく位置合わせができたかの指標として、差の二乗や相互相関、相互情報量（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.
【発明の概要】
【発明が解決しようとする課題】
【００１０】
前述した従来の相互情報量を用いたレジストレーションによると、計算量が多く、処理時間が長いという問題がある。このため、相互情報量を用いたレジストレーションの高速化が重要な課題となっている。高速化のために、画素データのメモリアクセスが連続となる方法が提案されている（例えば、特許文献３参照。）。
【００１１】
一方、近年のマルチプロセッサやＧＰＵ（グラフィックプロセッシングユニット）の進歩に伴い、並列処理によるレジストレーションを高速化の期待が高まっている。また、ＧＰＵを用いた高速化方法も提案されている（例えば、非特許文献２参照。）。しかし、並列処理を行うとしても、相互情報量を計算するために必要な結合ヒストグラムの算出の計算量が多いという問題がある。
【００１２】
すなわち、相互情報量を計算するためには、位置合わせの対象となる二つの画像の結合ヒストグラムを計算する必要がある。結合ヒストグラムは２次元のヒストグラムであり、ヒストグラムに含まれるビン数が大きくなる。例えば、一つの画像のヒストグラムのビン数を６４とすると、結合ヒストグラム全体のビン数は６４×６４＝４０９６となる。
【００１３】
一般に、ＧＰＵは、高速にアクセスできるローカルメモリの容量が小さい。ＧＰＵを用いて結合ヒストグラムを計算する場合、このような大きなビン数のヒストグラムをローカルメモリに格納することができず、アクセスが低速なグローバルメモリにデータを格納する必要がある。このため、結合ヒストグラムを計算する処理に時間がかかるという問題がある。この問題は、高速にアクセスできるローカルメモリが少ない並列計算機においても同様である。
【００１４】
本発明は、画像レジストレーションの処理速度を向上する画像処理装置を提供することを目的とする。
【課題を解決するための手段】
【００１５】
本発明の代表的な一例を示せば以下の通りである。すなわち、２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせる画像処理システムであって、前記画像処理システムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、前記画像処理システムは、前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手段と、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２の手段と、前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手段と、前記生成された結合ヒストグラムから相互情報量を計算する第４の手段と、を備える。
【発明の効果】
【００１６】
本発明の実施の形態によると、各プロセッシングユニットに保持されるヒストグラムのビン数を削減でき、使用されるローカルメモリの容量を削減することができる。
【図面の簡単な説明】
【００１７】
【図１】本発明の実施の形態の画像処理装置の論理的な構成を示す図である。
【図２】本発明の実施の形態の画像処理装置のハードウェア構成を示すブロック図である。
【図３】本発明の実施の形態の画像処理装置の構成の一例を示すブロック図である。
【図４】本発明の実施の形態の処理を示すフローチャートである。
【図５Ａ】レジストレーションの対象となる参照画像を説明する図である。
【図５Ｂ】レジストレーションの対象となる浮動画像を説明する図である。
【図６Ａ】レジストレーション処理の概要を説明する図である。
【図６Ｂ】レジストレーション処理の概要を説明する図である。
【図７Ａ】従来技術によって結合ヒストグラムを計算する一般的な方法を説明する図である。
【図７Ｂ】本発明によって結合ヒストグラムを計算する方法を説明する図である。
【図８】本発明の実施の形態の参照画像ヒストグラム作成処理の手順の詳細を示すフローチャートである。
【図９】本実施の形態の座標リスト群の詳細を説明する図である。
【図１０】本発明の実施の形態の相互情報量計算処理の手順の詳細を示すフローチャートである。
【図１１】本発明の実施の形態の浮動画像のヒストグラムを説明する図である。
【図１２】本発明の実施の形態の結合ヒストグラムを説明する図である。
【図１３】本発明の実施の形態の画像処理装置で実行される処理を示すシーケンス図である。
【図１４Ａ】本発明の実施の形態の結合ヒストグラムの具体例を説明する図である。
【図１４Ｂ】本発明の実施の形態の結合ヒストグラムの具体例を説明する図である。
【図１５】本発明の実施の形態の画像処理装置の構成の一例を示すブロック図である。
【発明を実施するための形態】
【００１８】
＜構成＞
図１は、本発明の実施の形態の画像処理装置の論理的な構成を示す図である。
【００１９】
本発明の実施の形態の画像処理装置は、前処理部１６及びレジストレーション部１７を備える。本実施の形態の画像処理装置は、レジストレーションの対象となる参照画像１及び浮動画像２が入力されると、レジストレーション済の浮動画像１５を生成して、生成された画像を出力する。
【００２０】
前処理部１６に参照画像１が入力されると、ヒストグラム生成部４は１次元のヒストグラムを生成し、参照画像ヒストグラムの各ビンの座標リストである座標リスト群５を生成する。
【００２１】
レジストレーション部１７に前記座標リスト群５及び浮動画像２が入力されると、相互情報量計算部１８は、幾何変換情報３を参照して、相互情報量１２を求める。なお、相互情報量計算部１８の処理内容の詳細は後述する。相互情報量最大化部１３は、求められた相互情報量１２の幾何変換情報３を更新し、相互情報量計算部１８が、相互情報量１２を再計算し、相互情報量１２を最大化する。浮動画像幾何変換部１４は、相互情報量が最大となる幾何変換情報３を用いて浮動画像２を幾何変換することによって、レジストレーション済の浮動画像１５を生成し、生成された画像を出力する。以上の処理によってレジストレーションが実行される。
【００２２】
次に、相互情報量計算部１８によって実行される処理の概要について説明する。相互情報量計算部１８は、幾何変換情報３を参照して、座標リスト群５の各画素を座標変換し、変換後の座標を求める。浮動画像参照・量子化部７は、変換後の座標に基づいて浮動画像２を参照し、輝度値を量子化し、参照画像ヒストグラムのビンごとに浮動画像のヒストグラム８を生成する。結合ヒストグラム生成部９は、参照画像ヒストグラムのビンごとの浮動画像のヒストグラム８から結合ヒストグラム１０を求める。相互情報量計算部１１は、結合ヒストグラム１０を用いて相互情報量１２を求める。
【００２３】
図２は、本発明の実施の形態の画像処理装置のハードウェア構成を示すブロック図である。
【００２４】
本実施の形態の画像処理装置は、ＧＰＵ及び／又は並列プロセッサシステムを備える一般的な計算機上に実装可能である。画像処理装置は、ＣＰＵ４０、ＲＯＭ２１、ＲＡＭ２２、記憶装置２３、ＧＰＵ４１、画像入力部３０、媒体入力部２６、入力制御部２８及び画像生成部２９を備える。ＣＰＵ４０、ＲＯＭ２１、ＲＡＭ２２、ハードディスク２３、ＧＰＵ４１、画像入力部３０、媒体入力部２６、入力制御部２８及び画像生成部２９は、データバス２４によって相互に接続される。
【００２５】
ＣＰＵ４０は、ＲＯＭ２１又はＲＡＭ２２に記憶されたプログラムを実行することによって、各種処理を実行する。
【００２６】
ＲＯＭ２１及びＲＡＭ２２は、各種処理を実行するために必要なプログラム及びデータを記憶する。ＲＯＭ２１は、データの読み出し専用の記憶媒体である。ＲＡＭ２２は、データの読み書きが可能な記憶媒体である。
【００２７】
記憶装置２３は、入力画像などを格納する磁気記憶装置である。記憶装置２３は、不揮発性半導体記憶媒体（例えば、フラッシュメモリ）を備えてもよい。また、ネットワークなどを介して接続された外部記憶装置を利用してもよい。
【００２８】
ＧＰＵ４１は、画像データの処理に適する集積回路であり、その構成は図３を用いて詳述する。なお、ＧＰＵ４１は並列プロセッサシステムでもよい。
【００２９】
本実施の形態の画像処理装置において実行されるプログラムは、記憶媒体２５（例えば、光ディスク）に格納されており、媒体入力部２６（例えば、光ディスクドライブ）によって読み込まれ、読み込まれたプログラムがＲＡＭ２２に格納される。なお、記憶装置２３に当該プログラムを格納し、記憶装置２３からプログラムをＲＡＭ２２にロードしてもよい。また、ＲＯＭ２１にあらかじめプログラムを記憶させてもよい。
【００３０】
画像処理装置においてプログラムが実行されることによって、図１に示す前処理部１６及びレジストレーション部１７が機能する。ＲＡＭ２２は、図１に示す座標リスト群５、幾何変換情報３、浮動画像のヒストグラム８、結合ヒストグラム１０、相互情報量１２及びレジストレーション済浮動画像１５を格納する。
【００３１】
なお、ＧＰＵ４１がレジストレーション処理を実行する場合、ＧＰＵ４１上のグローバルメモリが、座標リスト群５、幾何変換情報３、浮動画像のヒストグラム８、結合ヒストグラム１０、相互情報量１２及びレジストレーション済浮動画像１５を格納してもよい。
【００３２】
画像入力部３０は、画像撮像装置２０が接続され、画像撮像装置２０によって撮影された画像が入力されるインターフェースである。画像撮像装置２０から入力された画像は、ＣＰＵ４０によって各種処理がされる。
【００３３】
媒体入力部２６は、記憶媒体２５に記憶されたデータ及びプログラムを読み出す。記憶媒体２５から読み出されたデータ及びプログラムは、ＣＰＵ４０によって、ＲＡＭ２２又はハードディスク２３に格納される。
【００３４】
入力制御部２８は、入力装置２７（例えば、キーボード）から、ユーザによって入力されたデータを受け付けるインターフェースである。入力制御部２８が受けたデータは、ＣＰＵ４０によって処理される。
【００３５】
画像生成部２９は、図１に示すレジストレーション部１７によってレジストレーション処理がされた浮動画像１５から画像データを生成し、生成された画像データをディスプレイ１４に送る。
【００３６】
図３は、ＧＰＵ４１の構成の一例を示すブロック図である。
【００３７】
ＧＰＵ４１は、グローバルメモリ４０３及び複数のマルチプロセッシングユニット４０５を備える。各マルチプロセッシングユニット４０５は、ローカル共有メモリ４０２及び複数のプロセッシングユニット４０１を備える。プロセッシングユニット４０１は、同じマルチプロセッシングユニット４０５内のローカル共有メモリ４０２を読み書き可能である。しかし、プロセッシングユニット４０１は、別のマルチプロセッシングユニット４０５のローカル共有メモリ４０２にはアクセスできない。各プロセッシングユニット４０１は、グローバルメモリ４０３を読み書き可能であるので、グローバルメモリ４０３には、複数のプロセッシングユニット４０１で共有されるデータが格納される。
【００３８】
なお、一般的にＧＰＵのローカル共有メモリ４０２は高速に読み書き可能であるが、グローバルメモリ４０３へのアクセスは、ローカル共有メモリ４０２のアクセスと比較すると低速である。このため、ＧＰＵを用いて効率よく処理を実行するためには、グローバルメモリ４０３のアクセスを削減し、ローカル共有メモリ４０２を用いて処理を実行することが重要である。ただし、ローカル共有メモリ４０２は、一般にその容量に制限があるので、ローカル共有メモリの容量内で効率のよいアルゴリズムを用いることがＧＰＵを用いた処理の高速化において重要である。
【００３９】
図１５は、図２の構成において、ＧＰＵ４１の代わりにマルチプロセッサシステムを用いる場合のマルチプロセッサシステムの構成の一例を示すブロック図である。
【００４０】
マルチプロセッサシステムは、グローバルメモリ１５０５及び複数のマルチプロセッシングエレメント１５０３を備える。各マルチプロセッシングエレメント１５０３は、プロセッシングユニット１５０１とローカルメモリ１５０２を備える。ローカルメモリはバス１５０４を介してグローバルメモリ１５０５と結合される。またバス１５０４はホストインターフェース１５０６を介して、図２のバス２４と接続される。
図１５に示すマルチプロセッサシステムを用いる場合には，ローカルメモリ１５０２は高速に読み書き可能であるが、グローバルメモリ１５０５へのアクセスは、ローカルメモリ１５０２のアクセスと比較すると低速である。このため、マルチプロセッサシステムを用いて効率よく処理を実行するためには、グローバルメモリ１５０５のアクセスを削減し、ローカルメモリ１５０２を用いて処理を実行することが重要である。ただし、ローカルメモリ１５０２は、一般にその容量に制限があるので、ローカルメモリの容量内１５０２で効率のよいアルゴリズムを用いることがＧＰＵを用いた処理の高速化において重要である。
【００４１】
本願発明は、以上の点を考慮した効率的なレジストレーション方法を提案するものである。
【００４２】
＜動作＞
次に、本発明の実施の形態の動作について説明する。
【００４３】
図４は、本発明の実施の形態の処理を示すフローチャートである。
【００４４】
図４に示す処理のうち、ステップＳ１０１及びＳ１０２の処理が、前処理部１６によって実行される処理であり、ステップＳ１０３〜Ｓ１０８の処理が、レジストレーション部１７によって実行される処理である。
【００４５】
まず、画像入力部３０が、参照画像１の入力を受け付ける（Ｓ１０１）。
【００４６】
次に、ヒストグラム生成部４が、入力された参照画像１から１次元のヒストグラムを作成する（Ｓ１０２）。さらに、ヒストグラムの各ビンに属する参照画像の画素の座標を座標リスト群５に格納する。なお、Ｓ１０２の処理内容の詳細は、図８を用いて後述する。
【００４７】
次に、画像入力部３０が、浮動画像２の入力を受け付ける（Ｓ１０３）。その後、幾何変換情報３の初期値を設定する（Ｓ１０４）。レジストレーションでは、二つの画像の相互情報量が最大となる幾何変換を求めるために収束計算を行う。このため、幾何変換情報３の初期値が必要となる。幾何変換情報の初期値は、ユーザが与えてもよく、まったく幾何変換しない、すなわち恒等変換を初期値として与えてもよい。
【００４８】
次に、相互情報量計算部１８が、相互情報量１２を求め、メモリに格納する（Ｓ１０５）。その後、求められた相互情報量が収束したかどうかを判定する（Ｓ１０６）。判定の結果、相互情報量が収束していない場合、より高い相互情報量を得るために、幾何変換情報３を更新する（Ｓ１０７）。そして、更新された幾何変換情報３を用いて、相互情報量１２を再計算し（Ｓ１０５）、再度収束を判定する（Ｓ１０６）。
【００４９】
一方、相互情報量が収束している場合、求められた幾何変換情報３を用いて、浮動画像２を幾何変換し、レジストレーションされた浮動画像１５にＲＡＭ２２に格納する（Ｓ１０８）。以上の処理によって、画像のレジストレーションが完了する。
【００５０】
図５Ａは、レジストレーションの対象となる参照画像１を説明する図であり、図５Ｂは、レジストレーションの対象となる浮動画像２を説明する図である。
【００５１】
図５Ａ、図５Ｂでは、簡単のため画像を２次元として図示しているが、参照画像１及び浮動画像２は３次元画像である。参照画像１は、図５Ａに示すように、Ｘ軸方向にＲｘ画素、Ｙ軸方向にＲｙ画素、図示は省略したがＺ軸方向にＲｚ画素、合計Ｒｘ×Ｒｙ×Ｒｚ個の画素３０１から構成される３次元画像である。同様に、浮動画像２は、Ｘ軸方向にＦｘ画素、Ｙ軸方向にＦｙ画素、図示は省略したがＺ軸方向にＦｚ画素、合計Ｆｘ×Ｆｙ×Ｆｚ個の画素３０２から構成される３次元画像である。
【００５２】
図６Ａ、図６Ｂは、レジストレーション処理の概要を説明する図である。
【００５３】
図６Ａに示すように、参照画像１と浮動画像２は類似した画像ではあるが、その位置又は方向が違っている。例えば、異なる視点から同一の患者のＣＴ画像、ＰＥＴ画像などを撮影した場合、類似した複数の画像が取得される。また、ＣＴ画像とＰＥＴ画像のように、輝度値の性質が異なる場合も、類似した複数の画像が取得される。
【００５４】
そして、図６Ｂに示すように、浮動画像２にレジストレーション処理を行うことによって、浮動画像２を幾何変換し、幾何変換された浮動画像１３と参照画像１とを位置合わせすることができる。
【００５５】
レジストレーション処理における幾何変換は、３次元の移動及び回転から構成される剛体変換を用いることが可能である。また、一般的なアフィン変換や、Ｂスプライン関数を用いた非剛体変換を適用することもできる。
【００５６】
剛体変換を用いた幾何変換には、式（１）を用いることができる。
【００５７】
【数１】
【００５８】
式（１）において，（ｘ，ｙ，ｚ）は幾何変換前の座標、（ｘ’，ｙ’，ｚ’）は変換後の座標、Ｒｘ（θｘ）、Ｒｙ（θｙ）、Ｒｚ（θｚ）はそれぞれ，Ｘ軸周りの回転行列、Ｙ軸周りの回転行列、Ｚ軸周りの回転行列である。また、（ｔｘ，ｔｙ，ｔｚ）は移動ベクトルである、すなわち、３次元画像の剛体変換は６自由度を有する。なお、２次元画像の剛体変換では、ｘ方向への移動、ｙ方向への移動及び回転の３自由度を有する。
【００５９】
アフィン変換を用いた幾何変換には、式（２）を用いることができる。
【００６０】
【数２】
【００６１】
式（２）において、（ｘ，ｙ，ｚ）は幾何変換前の座標、（ｘ’，ｙ’，ｚ’）は変換後の座標、ａ、ｂ、ｃ、ｄ、ｅ、ｆ、ｇ、ｈ、ｉはアフィン変換のパラメータ、（ｔｘ，ｔｙ，ｔｚ）は移動ベクトルである。
【００６２】
Ｂスプライン関数を用いた非剛体変換には、式（３）を用いることができる。
【００６３】
【数３】
【００６４】
Ｂスプライン関数を用いた非剛体変換では、変換パラメータはＢスプライン立体の格子点位置となる。式（３）において、φｉ，ｊ，ｋは、ｎｘ×ｎｙ×ｎｚ個の制御点である。
【００６５】
ここで、本発明との対比のため従来技術における結合ヒストグラムの計算について説明する。図７Ａは、従来技術によって結合ヒストグラムを計算する一般的な方法を説明する図である。
【００６６】
従来の結合ヒストグラムの計算方法では、図７Ａに示すように、参照画像１の画素値と、幾何変換された浮動画像の対応する位置における画素値とを求める。画素値には、例えば、各画素の輝度値を用いることができる。なお、特定波長の輝度値を用いてもよい。
【００６７】
そして、この求められた二つの輝度値に基づいて、結合ヒストグラムを求める。具体的には、参照画像のヒストグラムのビン数Ｎｒを参照して、参照画像の輝度値が属する参照画像ヒストグラムのビンｉを求める。同様に、幾何変換された浮動画像のヒストグラムのビン数Ｎｆを参照して、浮動画像の輝度値が属する浮動画像のヒストグラムのビンｊを求める。求められた（ｉ，ｊ）から対応する結合ヒストグラムのビンの頻度を更新する。
【００６８】
以上説明したように、従来技術によると、参照画像１と幾何変換された浮動画像１３とに基づいて結合ヒストグラム１０を求めていた。レジストレーション処理における画像の類似度の指標である相互情報量を求めるためには、前述した結合ヒストグラムの算出が必要である。しかし、結合ヒストグラムは２次元のヒストグラムであり、その容量が大きい性質があることから、結合ヒストグラムを記憶するために大きなメモリ容量が必要である。相互情報量の算出では、ヒストグラムのビン数は６４〜２５６程度である。ビン数が６４の場合では、結合ヒストグラムの総ビン数は６４×６４＝４０９６ビンとなる。
【００６９】
一般に、ＧＰＵ及び並列プロセッサでは、各プロセッサユニットが局所的に利用できるメモリ容量に制限がある場合が多い。これは、例えば、図３に示すＧＰＵ装置において、ローカル共有メモリ４０２の容量が制限されているためである。このため、共有メモリの容量の不足が、ＧＰＵや並列プロセッサで結合ヒストグラムを計算する場合の問題点となっていた。
【００７０】
一方、本実施の形態では、図７Ｂに示すように、参照画像ヒストグラムの座標リスト５を作成し、その後、参照画像ヒストグラムのビンごとに浮動画像の１次元ヒストグラム８を生成して、結合ヒストグラムを生成する。このように、各プロセッサユニットが２次元の結合ヒストグラムを保持する必要がなく、１次元のヒストグラムを保持するだけでよい計算方法を提供することによって、前述した問題を解決することができる。
【００７１】
相互情報量の計算には、まず、参照画像１と浮動画像２の同じ位置の画素値から生成される２次元ヒストグラムを計算する。その後、計算された２次元ヒストグラムを用いて計算される確率密度関数を用いて、相互情報量を計算することができる。
【００７２】
相互情報量とは、２つの事象ＡとＢについて、事象Ａが有している事象Ｂに関する情報量を定量化した尺度である。レジストレーションで用いられる正規化相互情報量ＮＭＩ（Ａ，Ｂ）は事象Ａと事象Ｂの２次元の結合ヒストグラムＨｉｓｔ（Ａ，Ｂ）から下式（４）によって求められる。
【００７３】
【数４】
【００７４】
ここで、Ｈ（Ａ）は事象Ａのエントロピー、Ｈ（Ｂ）は事象Ｂのエントロピー、Ｈ（Ａ，Ｂ）は事象Ａ，Ｂの結合エントロピーである。ｐ（ａ）はａの確率密度分布、ｐ（ｂ）はｂの確率密度分布、ｐ（ａ，ｂ）はａ，ｂの同時確率分布であり、式（４）のように、Ｈｉｓｔ（Ａ，Ｂ）から求められる。
【００７５】
事象Ａと事象Ｂが完全に独立である場合には、Ｈ（Ａ，Ｂ）＝Ｈ（Ａ）＋Ｈ（Ｂ）となる。また、事象Ａと事象Ｂが完全に従属である場合には、Ｈ（Ａ，Ｂ）＝Ｈ（Ａ）＝Ｈ（Ｂ）となる。以上から，ＮＭＩ（Ａ，Ｂ）のとりうる範囲は、１以上２以下となる。
【００７６】
図８は、図４のフローチャートにおける参照画像ヒストグラム作成（Ｓ１０２）の処理手順の詳細を示すフローチャートである。
【００７７】
まず、変数Ｚを１から参照画像のＺ軸方向の画素数Ｒｚまで増加させながら、ステップＳ２０２からＳ２０８の処理を繰り返し実行する（Ｓ２０１）。次に、変数Ｙを１から参照画像のＹ軸方向の画素数Ｒｙまで増加させながら、ステップＳ２０３からＳ２０７の処理を繰り返し実行する（Ｓ２０２）。さらに、変数Ｘを１から参照画像のＸ軸方向の画素数Ｒｘまで増加させながら、ステップＳ２０４からＳ２０６の処理を繰り返し実行する（Ｓ２０３）。以上の処理によって、参照画像のすべての画素（Ｘ，Ｙ，Ｚ）について、ステップＳ２０４からＳ２０６の処理が実行される。
【００７８】
ステップＳ２０４では、参照画像の画素値Ｒ（Ｘ，Ｙ，Ｚ）を変数Ｉに代入する。次に、画素値Ｉを量子化して、量子化された画素値Ｉを用いて参照画像のヒストグラムのビンを求める（Ｓ２０５）。この処理は、例えば、参照画像の輝度値の最大値と最小値との間をＮｒ個の領域に分割し、輝度値が属する領域の番号からビンを求める。
【００７９】
その後、ヒストグラムのビンと、座標値（Ｘ，Ｙ，Ｚ）とを用いて、座標リスト群５にデータを格納する（Ｓ２０６）。座標リスト群５へのデータの格納方法の詳細は、図９を用いて後述する。
【００８０】
なお、前述した説明では、参照画像の全画素について、ヒストグラム及びヒストグラムリストを求めたが、ダウンサンプリングを行って少ないサンプル数の画素についてヒストグラム及びヒストグラムリストを求めてもよい。また、スーパーサンプリングを行って、サンプル数を増やしてもよい。さらに、等間隔にサンプリングするのではなく、ランダムにサンプリングしてもよい。
【００８１】
図９は、本実施の形態の座標リスト群５の詳細を説明する図である。
【００８２】
ヒストグラムリスト５は、Ｎｒ個の可変長のリストを含む。各リストの構成要素は、画素の座標である。
【００８３】
図８のステップＳ２０６において、ビン番号ｂに座標（Ｘ，Ｙ，Ｚ）を登録する場合、ｂ番目のリストの最後尾に座標（Ｘ，Ｙ，Ｚ）を追加する。
【００８４】
図９に示すように、参照画像１のすべての画素に対して、座標（Ｘ，Ｙ，Ｚ）を追加した場合の、ｂ番目のリストの長さをＢ（ｂ）とする。Ｂ（ｂ）は参照画像の１次元のヒストグラムと等価である。
【００８５】
図１０は、図４のフローチャートにおける相互情報量計算（Ｓ１０５）の処理手順の詳細を示すフローチャートである。
【００８６】
まず、変数ｂを１から参照画像のヒストグラムのビン数Ｎｒまで増加させながら、ステップＳ３０２からステップＳ３０８を繰り返し実行する（Ｓ３０１）。
【００８７】
次に、変数ｉを１からヒストグラムの長さＢ（ｂ）まで増加させながら、ステップＳ３０３からステップＳ３０７を繰り返し実行する（Ｓ３０２）。
【００８８】
ループ内では、座標リスト群５を参照して、対応する座標（Ｘ，Ｙ，Ｚ）を取り出し（Ｓ３０３）、幾何変換情報３を用いて座標（Ｘ，Ｙ，Ｚ）を座標変換し、対応する浮動画像の座標系に変換する（Ｓ３０４）。
【００８９】
次に、浮動画像の対応する画素値Ｒ（Ｘ，Ｙ，Ｚ）を変数Ｉに代入し（Ｓ３０５）、画素値Ｉから浮動画像のヒストグラムのビンを求める。この処理は、浮動画像の最大値と最小値とをＮｆ個の領域に分割し、輝度値が属する領域の番号からビンを求めることで実現できる（Ｓ３０６）。
【００９０】
その後、量子化された値に対応するビン番号ｂ’の浮動画像ヒストグラム８を更新する（Ｓ３０７）。処理の詳細は後述する。
【００９１】
ステップＳ３０１からステップＳ３０９のループによって、参照画像の全画素に対応する参照画像ヒストグラムのビンごとの浮動画像のヒストグラム８が求められる。ヒストグラム８は、図１１に示すように、ビン数Ｎｆの１次元ヒストグラムをＮｒ個含む。
【００９２】
その後、参照画像ヒストグラムのビンごとの浮動画像のヒストグラム８を合成して、２次元の結合ヒストグラム１０を求める（Ｓ３１０）。結合ヒストグラム１０は、図１２に示すように、ヒストグラム８に含まれる一次元ヒストグラムを結合したもので、Ｎｒ×Ｎｆのマトリックスで構成されている。
【００９３】
結合ヒストグラム１０の具体例を図１４Ａ、図１４Ｂに示す。図１４Ａに示すように、参照画像１と浮動画像２との位置を合わせない場合、結合ヒストグラム１０には複数のピークが現れる。しかし、図１４Ｂに示すように、参照画像１と浮動画像２との位置が合っている状態では、結合ヒストグラム１０には一つのピークが現れる。
【００９４】
その後、結合ヒストグラム１０から相互情報量１２を求める（Ｓ３１１）。
【００９５】
次に、図１０のステップＳ３０７の処理の詳細について説明する。
【００９６】
ステップＳ３０７では、参照画像のビンがｂ、浮動画像のビンがｂ’である場合、ｂ番目のヒストグラムのｂ’のビンの頻度に１を加算する。すなわち、２次元の結合ヒストグラム１０ではなく、１次元のヒストグラム８が更新される。
【００９７】
このことから、ステップＳ３０２からのループであるステップＳ３０３からステップＳ３０８の処理をＧＰＵ又は並列プロセッサで実行する場合、各プロセッシングユニット４０１は２次元のヒストグラムをローカル共有メモリ４０２に保持する必要がなく、１次元のヒストグラムをローカル共有メモリ４０２に保持するだけでよい。
【００９８】
さらに、ステップＳ３０１からのループであるステップＳ３０２からステップ３０９の処理をＧＰＵ又は並列プロセッサで実行することも可能である。この場合にも、各プロセッシングユニット４０１は２次元のヒストグラムをローカル共有メモリ４０２に保持する必要がなく、１次元のヒストグラムをローカル共有メモリ４０２に保持するだけでよい。
【００９９】
３次元画像のレジストレーション実験によれば、本発明は非特許文献２の方法と比較して、より高速である。具体的な実験は、ファントムのＣＴ画像を用いてレジストレーションの性能を評価した。使用したファントム画像は５１２×５１２×２００ボクセルの３次元画像である。ファントム画像に対して、Ｘ軸、Ｙ軸及びＺ軸に対して１０度回転したボリュームデータを作成しレジストレーションを行った。レジストレーション実験によればもとの位置に正しく座標変換された。
【０１００】
レジストレーションに使用するサンプル数を種々に変更して計算時間を計測した。サンプル数が６４×６４×６４の場合の計算時間は、ＣＰＵ実装が１４．２秒、非特許文献２の方法が１．８２秒、本発明による方法が１．１６秒であった。サンプル数が２５６×２５６×２５６の場合にはＣＰＵ実装の計算時間は７８０秒、非特許文献２の方法が１４．９秒、本発明による方法が５．７秒であった。以上のように本発明による方法はＣＰＵ実装と比較して数１０倍以上、非特許文献２の方法と比較して１．５倍から２．６倍程度の高速化が実現されている。
【０１０１】
図１３は、本実施の形態の画像処理装置で実行される処理を示すシーケンス図である。
【０１０２】
まず、レジストレーション命令が発行されると、入力装置２７は、ユーザからの参照画像及び浮動画像の指定を受けて、画像の指定を入力制御部２８を介してＣＰＵ４０に伝える（１３０１）。ＣＰＵ４０は、指定された参照画像及び浮動画像を取得し、取得した参照画像及び浮動画像を記憶装置２３に格納する（１３０２）。ＧＰＵ４１は、記憶装置２３に格納された参照画像を取得し、参照画像ヒストグラムを生成し（１３０３）、生成された参照画像ヒストグラムを記憶装置に格納する。
【０１０３】
さらに、入力装置２７は、ユーザからの幾何変換情報の指定を受けて、入力された幾何変換情報を入力制御部２８を介してＣＰＵ４０に伝える（１３０４）。この幾何変換情報は、幾何変換の初期条件を含み、さらに、幾何変換情報を更新する条件を含んでもよい。ＣＰＵ４０は、入力された幾何変換情報をＧＰＵ４１に転送する（１３０５）。
【０１０４】
ＧＰＵ４１は、幾何変換を受けると、浮動画像を幾何変換し、幾何変換された浮動画像の画素値をサンプリングし、参照画像ヒストグラムのビンごとに浮動画像の１次元ヒストグラム８を生成する（１３０６）。このため、ＧＰＵ４１は、浮動画像の１次元ヒストグラム８をローカル共有メモリ４０２に格納すればよく。小容量で高速アクセス可能なメモリのみを用いて（すなわち、大容量だがアクセス速度が低いグローバルメモリ４０３を用いることなく）、浮動画像の１次元ヒストグラム８を生成することができる。その後、ＧＰＵ４１は、生成された浮動画像の１次元ヒストグラム８を、記憶装置２３に格納する。
【０１０５】
ＧＰＵ４１は、全ての１次元ヒストグラムの計算を終了した後、記憶装置２３から１次元ヒストグラム８を読み出して、読み出された１次元ヒストグラム８を結合する。そして結合されたヒストグラム１０から相互情報量を計算する（１３０７）。計算された相互情報量はＣＰＵ４０に転送される（１３０８）。
【０１０６】
ＣＰＵ４０は、相互情報量が最大となるように幾何変換情報を更新し（１３０９）、更新された幾何変換情報をＧＰＵに転送する（１３１０）。
【０１０７】
ＧＰＵ４１は、転送された幾何変換情報に基づいて、１３０６〜１３０９の処理と同様に、幾何変換された浮動画像の１次元ヒストグラム８を生成し、結合ヒストグラムを生成し、相互情報量を計算し（１３１１）、計算された相互情報量をＣＰＵ４０に転送する（１３１２）。
【０１０８】
ＣＰＵ４０は、相互情報量が最大となる幾何変換情報で浮動画像を変換し（１３１３）、この条件で幾何変換された浮動画像（レジストレーション済画像）をディスプレイ３１に表示する（１３１４）。
【０１０９】
以上に説明したように、本発明の実施の形態では、２次元の結合ヒストグラムの計算をＧＰＵ又は並列プロセッサで実行する場合に、各プロセッシングユニットは１次元のヒストグラムにアクセスするので、必要なメモリ容量を削減することができる。このため、結合ヒストグラムを効率よく計算でき、レジストレーション計算を効率的に処理することができる。
【符号の説明】
【０１１０】
１ 参照画像
２ 浮動画像
３ 幾何変換情報
５ 参照画像ヒストグラム（座標リスト群）
８ 浮動画像ヒストグラム
１０ 結合ヒストグラム
１２ 相互情報量
１５ レジストレーション済み浮動画像
１６ 前処理部
１７ レジストレーション部
４１ ＧＰＵ
【技術分野】
【０００１】
本発明は、２次元又は３次元の複数枚の画像の位置を合わせる画像処理装置に関する。
【背景技術】
【０００２】
２次元又は３次元の複数枚の画像の位置合わせ技術は、様々な分野で重要な技術である。
【０００３】
例えば、医用画像の分野では、ＣＴ画像、ＭＲ画像、ＰＥＴ画像、超音波画像など、様々な種類の３次元画像の位置を合わせて、重ね合わせて表示するために、画像の位置合わせ技術が利用される。このような表示方法はフュージョン画像表示と呼ばれており、画像の特徴を生かした表示が可能となる。たとえば、ＣＴ画像は詳細な形状の表示に適しており、ＰＥＴ画像は代謝や血流などの体の機能の表示に適している（例えば、特許文献１参照。）。
【０００４】
さらに、同一患者の経過を観察するためには、時系列に取得した医用画像の位置を合わせることによって、病変を容易に観察することが可能である。
【０００５】
医用画像の分野では、３次元画像だけではなく、２次元画像の位置合わせも重要である。例えば、エネルギーを変更した２種類の胸部放射線画像を用いることによって、骨部を除去した軟部画像を生成することができ、軟部を詳細に観察することが可能となる。この場合にも、２つの画像は撮像時刻が異なるため、位置合わせが必要となる（例えば、特許文献２参照。）。
【０００６】
一方、リモートセンシングの分野では、様々なセンサ及び様々な方法で取得された画像データを総合的に検証することによって、多くの情報が取得できる。
【０００７】
２次元又は３次元の複数枚の画像の位置合わせ、すなわち、レジストレーションは、前述したように幅広い分野で応用されている。レジストレーションでは、固定される画像を参照画像と呼び、位置合わせのために座標変換される画像を浮動画像と呼ぶ。レジストレーションでは、浮動画像に対して座標変換を行い、変換された浮動画像と参照画像との類似度を求め、類似度が最大となる座標変換を繰り返し計算によって求める。レジストレーションでは、二つの画像が正しく位置合わせができたかの指標として、差の二乗や相互相関、相互情報量（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.
【発明の概要】
【発明が解決しようとする課題】
【００１０】
前述した従来の相互情報量を用いたレジストレーションによると、計算量が多く、処理時間が長いという問題がある。このため、相互情報量を用いたレジストレーションの高速化が重要な課題となっている。高速化のために、画素データのメモリアクセスが連続となる方法が提案されている（例えば、特許文献３参照。）。
【００１１】
一方、近年のマルチプロセッサやＧＰＵ（グラフィックプロセッシングユニット）の進歩に伴い、並列処理によるレジストレーションを高速化の期待が高まっている。また、ＧＰＵを用いた高速化方法も提案されている（例えば、非特許文献２参照。）。しかし、並列処理を行うとしても、相互情報量を計算するために必要な結合ヒストグラムの算出の計算量が多いという問題がある。
【００１２】
すなわち、相互情報量を計算するためには、位置合わせの対象となる二つの画像の結合ヒストグラムを計算する必要がある。結合ヒストグラムは２次元のヒストグラムであり、ヒストグラムに含まれるビン数が大きくなる。例えば、一つの画像のヒストグラムのビン数を６４とすると、結合ヒストグラム全体のビン数は６４×６４＝４０９６となる。
【００１３】
一般に、ＧＰＵは、高速にアクセスできるローカルメモリの容量が小さい。ＧＰＵを用いて結合ヒストグラムを計算する場合、このような大きなビン数のヒストグラムをローカルメモリに格納することができず、アクセスが低速なグローバルメモリにデータを格納する必要がある。このため、結合ヒストグラムを計算する処理に時間がかかるという問題がある。この問題は、高速にアクセスできるローカルメモリが少ない並列計算機においても同様である。
【００１４】
本発明は、画像レジストレーションの処理速度を向上する画像処理装置を提供することを目的とする。
【課題を解決するための手段】
【００１５】
本発明の代表的な一例を示せば以下の通りである。すなわち、２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせる画像処理システムであって、前記画像処理システムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、前記画像処理システムは、前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手段と、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２の手段と、前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手段と、前記生成された結合ヒストグラムから相互情報量を計算する第４の手段と、を備える。
【発明の効果】
【００１６】
本発明の実施の形態によると、各プロセッシングユニットに保持されるヒストグラムのビン数を削減でき、使用されるローカルメモリの容量を削減することができる。
【図面の簡単な説明】
【００１７】
【図１】本発明の実施の形態の画像処理装置の論理的な構成を示す図である。
【図２】本発明の実施の形態の画像処理装置のハードウェア構成を示すブロック図である。
【図３】本発明の実施の形態の画像処理装置の構成の一例を示すブロック図である。
【図４】本発明の実施の形態の処理を示すフローチャートである。
【図５Ａ】レジストレーションの対象となる参照画像を説明する図である。
【図５Ｂ】レジストレーションの対象となる浮動画像を説明する図である。
【図６Ａ】レジストレーション処理の概要を説明する図である。
【図６Ｂ】レジストレーション処理の概要を説明する図である。
【図７Ａ】従来技術によって結合ヒストグラムを計算する一般的な方法を説明する図である。
【図７Ｂ】本発明によって結合ヒストグラムを計算する方法を説明する図である。
【図８】本発明の実施の形態の参照画像ヒストグラム作成処理の手順の詳細を示すフローチャートである。
【図９】本実施の形態の座標リスト群の詳細を説明する図である。
【図１０】本発明の実施の形態の相互情報量計算処理の手順の詳細を示すフローチャートである。
【図１１】本発明の実施の形態の浮動画像のヒストグラムを説明する図である。
【図１２】本発明の実施の形態の結合ヒストグラムを説明する図である。
【図１３】本発明の実施の形態の画像処理装置で実行される処理を示すシーケンス図である。
【図１４Ａ】本発明の実施の形態の結合ヒストグラムの具体例を説明する図である。
【図１４Ｂ】本発明の実施の形態の結合ヒストグラムの具体例を説明する図である。
【図１５】本発明の実施の形態の画像処理装置の構成の一例を示すブロック図である。
【発明を実施するための形態】
【００１８】
＜構成＞
図１は、本発明の実施の形態の画像処理装置の論理的な構成を示す図である。
【００１９】
本発明の実施の形態の画像処理装置は、前処理部１６及びレジストレーション部１７を備える。本実施の形態の画像処理装置は、レジストレーションの対象となる参照画像１及び浮動画像２が入力されると、レジストレーション済の浮動画像１５を生成して、生成された画像を出力する。
【００２０】
前処理部１６に参照画像１が入力されると、ヒストグラム生成部４は１次元のヒストグラムを生成し、参照画像ヒストグラムの各ビンの座標リストである座標リスト群５を生成する。
【００２１】
レジストレーション部１７に前記座標リスト群５及び浮動画像２が入力されると、相互情報量計算部１８は、幾何変換情報３を参照して、相互情報量１２を求める。なお、相互情報量計算部１８の処理内容の詳細は後述する。相互情報量最大化部１３は、求められた相互情報量１２の幾何変換情報３を更新し、相互情報量計算部１８が、相互情報量１２を再計算し、相互情報量１２を最大化する。浮動画像幾何変換部１４は、相互情報量が最大となる幾何変換情報３を用いて浮動画像２を幾何変換することによって、レジストレーション済の浮動画像１５を生成し、生成された画像を出力する。以上の処理によってレジストレーションが実行される。
【００２２】
次に、相互情報量計算部１８によって実行される処理の概要について説明する。相互情報量計算部１８は、幾何変換情報３を参照して、座標リスト群５の各画素を座標変換し、変換後の座標を求める。浮動画像参照・量子化部７は、変換後の座標に基づいて浮動画像２を参照し、輝度値を量子化し、参照画像ヒストグラムのビンごとに浮動画像のヒストグラム８を生成する。結合ヒストグラム生成部９は、参照画像ヒストグラムのビンごとの浮動画像のヒストグラム８から結合ヒストグラム１０を求める。相互情報量計算部１１は、結合ヒストグラム１０を用いて相互情報量１２を求める。
【００２３】
図２は、本発明の実施の形態の画像処理装置のハードウェア構成を示すブロック図である。
【００２４】
本実施の形態の画像処理装置は、ＧＰＵ及び／又は並列プロセッサシステムを備える一般的な計算機上に実装可能である。画像処理装置は、ＣＰＵ４０、ＲＯＭ２１、ＲＡＭ２２、記憶装置２３、ＧＰＵ４１、画像入力部３０、媒体入力部２６、入力制御部２８及び画像生成部２９を備える。ＣＰＵ４０、ＲＯＭ２１、ＲＡＭ２２、ハードディスク２３、ＧＰＵ４１、画像入力部３０、媒体入力部２６、入力制御部２８及び画像生成部２９は、データバス２４によって相互に接続される。
【００２５】
ＣＰＵ４０は、ＲＯＭ２１又はＲＡＭ２２に記憶されたプログラムを実行することによって、各種処理を実行する。
【００２６】
ＲＯＭ２１及びＲＡＭ２２は、各種処理を実行するために必要なプログラム及びデータを記憶する。ＲＯＭ２１は、データの読み出し専用の記憶媒体である。ＲＡＭ２２は、データの読み書きが可能な記憶媒体である。
【００２７】
記憶装置２３は、入力画像などを格納する磁気記憶装置である。記憶装置２３は、不揮発性半導体記憶媒体（例えば、フラッシュメモリ）を備えてもよい。また、ネットワークなどを介して接続された外部記憶装置を利用してもよい。
【００２８】
ＧＰＵ４１は、画像データの処理に適する集積回路であり、その構成は図３を用いて詳述する。なお、ＧＰＵ４１は並列プロセッサシステムでもよい。
【００２９】
本実施の形態の画像処理装置において実行されるプログラムは、記憶媒体２５（例えば、光ディスク）に格納されており、媒体入力部２６（例えば、光ディスクドライブ）によって読み込まれ、読み込まれたプログラムがＲＡＭ２２に格納される。なお、記憶装置２３に当該プログラムを格納し、記憶装置２３からプログラムをＲＡＭ２２にロードしてもよい。また、ＲＯＭ２１にあらかじめプログラムを記憶させてもよい。
【００３０】
画像処理装置においてプログラムが実行されることによって、図１に示す前処理部１６及びレジストレーション部１７が機能する。ＲＡＭ２２は、図１に示す座標リスト群５、幾何変換情報３、浮動画像のヒストグラム８、結合ヒストグラム１０、相互情報量１２及びレジストレーション済浮動画像１５を格納する。
【００３１】
なお、ＧＰＵ４１がレジストレーション処理を実行する場合、ＧＰＵ４１上のグローバルメモリが、座標リスト群５、幾何変換情報３、浮動画像のヒストグラム８、結合ヒストグラム１０、相互情報量１２及びレジストレーション済浮動画像１５を格納してもよい。
【００３２】
画像入力部３０は、画像撮像装置２０が接続され、画像撮像装置２０によって撮影された画像が入力されるインターフェースである。画像撮像装置２０から入力された画像は、ＣＰＵ４０によって各種処理がされる。
【００３３】
媒体入力部２６は、記憶媒体２５に記憶されたデータ及びプログラムを読み出す。記憶媒体２５から読み出されたデータ及びプログラムは、ＣＰＵ４０によって、ＲＡＭ２２又はハードディスク２３に格納される。
【００３４】
入力制御部２８は、入力装置２７（例えば、キーボード）から、ユーザによって入力されたデータを受け付けるインターフェースである。入力制御部２８が受けたデータは、ＣＰＵ４０によって処理される。
【００３５】
画像生成部２９は、図１に示すレジストレーション部１７によってレジストレーション処理がされた浮動画像１５から画像データを生成し、生成された画像データをディスプレイ１４に送る。
【００３６】
図３は、ＧＰＵ４１の構成の一例を示すブロック図である。
【００３７】
ＧＰＵ４１は、グローバルメモリ４０３及び複数のマルチプロセッシングユニット４０５を備える。各マルチプロセッシングユニット４０５は、ローカル共有メモリ４０２及び複数のプロセッシングユニット４０１を備える。プロセッシングユニット４０１は、同じマルチプロセッシングユニット４０５内のローカル共有メモリ４０２を読み書き可能である。しかし、プロセッシングユニット４０１は、別のマルチプロセッシングユニット４０５のローカル共有メモリ４０２にはアクセスできない。各プロセッシングユニット４０１は、グローバルメモリ４０３を読み書き可能であるので、グローバルメモリ４０３には、複数のプロセッシングユニット４０１で共有されるデータが格納される。
【００３８】
なお、一般的にＧＰＵのローカル共有メモリ４０２は高速に読み書き可能であるが、グローバルメモリ４０３へのアクセスは、ローカル共有メモリ４０２のアクセスと比較すると低速である。このため、ＧＰＵを用いて効率よく処理を実行するためには、グローバルメモリ４０３のアクセスを削減し、ローカル共有メモリ４０２を用いて処理を実行することが重要である。ただし、ローカル共有メモリ４０２は、一般にその容量に制限があるので、ローカル共有メモリの容量内で効率のよいアルゴリズムを用いることがＧＰＵを用いた処理の高速化において重要である。
【００３９】
図１５は、図２の構成において、ＧＰＵ４１の代わりにマルチプロセッサシステムを用いる場合のマルチプロセッサシステムの構成の一例を示すブロック図である。
【００４０】
マルチプロセッサシステムは、グローバルメモリ１５０５及び複数のマルチプロセッシングエレメント１５０３を備える。各マルチプロセッシングエレメント１５０３は、プロセッシングユニット１５０１とローカルメモリ１５０２を備える。ローカルメモリはバス１５０４を介してグローバルメモリ１５０５と結合される。またバス１５０４はホストインターフェース１５０６を介して、図２のバス２４と接続される。
図１５に示すマルチプロセッサシステムを用いる場合には，ローカルメモリ１５０２は高速に読み書き可能であるが、グローバルメモリ１５０５へのアクセスは、ローカルメモリ１５０２のアクセスと比較すると低速である。このため、マルチプロセッサシステムを用いて効率よく処理を実行するためには、グローバルメモリ１５０５のアクセスを削減し、ローカルメモリ１５０２を用いて処理を実行することが重要である。ただし、ローカルメモリ１５０２は、一般にその容量に制限があるので、ローカルメモリの容量内１５０２で効率のよいアルゴリズムを用いることがＧＰＵを用いた処理の高速化において重要である。
【００４１】
本願発明は、以上の点を考慮した効率的なレジストレーション方法を提案するものである。
【００４２】
＜動作＞
次に、本発明の実施の形態の動作について説明する。
【００４３】
図４は、本発明の実施の形態の処理を示すフローチャートである。
【００４４】
図４に示す処理のうち、ステップＳ１０１及びＳ１０２の処理が、前処理部１６によって実行される処理であり、ステップＳ１０３〜Ｓ１０８の処理が、レジストレーション部１７によって実行される処理である。
【００４５】
まず、画像入力部３０が、参照画像１の入力を受け付ける（Ｓ１０１）。
【００４６】
次に、ヒストグラム生成部４が、入力された参照画像１から１次元のヒストグラムを作成する（Ｓ１０２）。さらに、ヒストグラムの各ビンに属する参照画像の画素の座標を座標リスト群５に格納する。なお、Ｓ１０２の処理内容の詳細は、図８を用いて後述する。
【００４７】
次に、画像入力部３０が、浮動画像２の入力を受け付ける（Ｓ１０３）。その後、幾何変換情報３の初期値を設定する（Ｓ１０４）。レジストレーションでは、二つの画像の相互情報量が最大となる幾何変換を求めるために収束計算を行う。このため、幾何変換情報３の初期値が必要となる。幾何変換情報の初期値は、ユーザが与えてもよく、まったく幾何変換しない、すなわち恒等変換を初期値として与えてもよい。
【００４８】
次に、相互情報量計算部１８が、相互情報量１２を求め、メモリに格納する（Ｓ１０５）。その後、求められた相互情報量が収束したかどうかを判定する（Ｓ１０６）。判定の結果、相互情報量が収束していない場合、より高い相互情報量を得るために、幾何変換情報３を更新する（Ｓ１０７）。そして、更新された幾何変換情報３を用いて、相互情報量１２を再計算し（Ｓ１０５）、再度収束を判定する（Ｓ１０６）。
【００４９】
一方、相互情報量が収束している場合、求められた幾何変換情報３を用いて、浮動画像２を幾何変換し、レジストレーションされた浮動画像１５にＲＡＭ２２に格納する（Ｓ１０８）。以上の処理によって、画像のレジストレーションが完了する。
【００５０】
図５Ａは、レジストレーションの対象となる参照画像１を説明する図であり、図５Ｂは、レジストレーションの対象となる浮動画像２を説明する図である。
【００５１】
図５Ａ、図５Ｂでは、簡単のため画像を２次元として図示しているが、参照画像１及び浮動画像２は３次元画像である。参照画像１は、図５Ａに示すように、Ｘ軸方向にＲｘ画素、Ｙ軸方向にＲｙ画素、図示は省略したがＺ軸方向にＲｚ画素、合計Ｒｘ×Ｒｙ×Ｒｚ個の画素３０１から構成される３次元画像である。同様に、浮動画像２は、Ｘ軸方向にＦｘ画素、Ｙ軸方向にＦｙ画素、図示は省略したがＺ軸方向にＦｚ画素、合計Ｆｘ×Ｆｙ×Ｆｚ個の画素３０２から構成される３次元画像である。
【００５２】
図６Ａ、図６Ｂは、レジストレーション処理の概要を説明する図である。
【００５３】
図６Ａに示すように、参照画像１と浮動画像２は類似した画像ではあるが、その位置又は方向が違っている。例えば、異なる視点から同一の患者のＣＴ画像、ＰＥＴ画像などを撮影した場合、類似した複数の画像が取得される。また、ＣＴ画像とＰＥＴ画像のように、輝度値の性質が異なる場合も、類似した複数の画像が取得される。
【００５４】
そして、図６Ｂに示すように、浮動画像２にレジストレーション処理を行うことによって、浮動画像２を幾何変換し、幾何変換された浮動画像１３と参照画像１とを位置合わせすることができる。
【００５５】
レジストレーション処理における幾何変換は、３次元の移動及び回転から構成される剛体変換を用いることが可能である。また、一般的なアフィン変換や、Ｂスプライン関数を用いた非剛体変換を適用することもできる。
【００５６】
剛体変換を用いた幾何変換には、式（１）を用いることができる。
【００５７】
【数１】
【００５８】
式（１）において，（ｘ，ｙ，ｚ）は幾何変換前の座標、（ｘ’，ｙ’，ｚ’）は変換後の座標、Ｒｘ（θｘ）、Ｒｙ（θｙ）、Ｒｚ（θｚ）はそれぞれ，Ｘ軸周りの回転行列、Ｙ軸周りの回転行列、Ｚ軸周りの回転行列である。また、（ｔｘ，ｔｙ，ｔｚ）は移動ベクトルである、すなわち、３次元画像の剛体変換は６自由度を有する。なお、２次元画像の剛体変換では、ｘ方向への移動、ｙ方向への移動及び回転の３自由度を有する。
【００５９】
アフィン変換を用いた幾何変換には、式（２）を用いることができる。
【００６０】
【数２】
【００６１】
式（２）において、（ｘ，ｙ，ｚ）は幾何変換前の座標、（ｘ’，ｙ’，ｚ’）は変換後の座標、ａ、ｂ、ｃ、ｄ、ｅ、ｆ、ｇ、ｈ、ｉはアフィン変換のパラメータ、（ｔｘ，ｔｙ，ｔｚ）は移動ベクトルである。
【００６２】
Ｂスプライン関数を用いた非剛体変換には、式（３）を用いることができる。
【００６３】
【数３】
【００６４】
Ｂスプライン関数を用いた非剛体変換では、変換パラメータはＢスプライン立体の格子点位置となる。式（３）において、φｉ，ｊ，ｋは、ｎｘ×ｎｙ×ｎｚ個の制御点である。
【００６５】
ここで、本発明との対比のため従来技術における結合ヒストグラムの計算について説明する。図７Ａは、従来技術によって結合ヒストグラムを計算する一般的な方法を説明する図である。
【００６６】
従来の結合ヒストグラムの計算方法では、図７Ａに示すように、参照画像１の画素値と、幾何変換された浮動画像の対応する位置における画素値とを求める。画素値には、例えば、各画素の輝度値を用いることができる。なお、特定波長の輝度値を用いてもよい。
【００６７】
そして、この求められた二つの輝度値に基づいて、結合ヒストグラムを求める。具体的には、参照画像のヒストグラムのビン数Ｎｒを参照して、参照画像の輝度値が属する参照画像ヒストグラムのビンｉを求める。同様に、幾何変換された浮動画像のヒストグラムのビン数Ｎｆを参照して、浮動画像の輝度値が属する浮動画像のヒストグラムのビンｊを求める。求められた（ｉ，ｊ）から対応する結合ヒストグラムのビンの頻度を更新する。
【００６８】
以上説明したように、従来技術によると、参照画像１と幾何変換された浮動画像１３とに基づいて結合ヒストグラム１０を求めていた。レジストレーション処理における画像の類似度の指標である相互情報量を求めるためには、前述した結合ヒストグラムの算出が必要である。しかし、結合ヒストグラムは２次元のヒストグラムであり、その容量が大きい性質があることから、結合ヒストグラムを記憶するために大きなメモリ容量が必要である。相互情報量の算出では、ヒストグラムのビン数は６４〜２５６程度である。ビン数が６４の場合では、結合ヒストグラムの総ビン数は６４×６４＝４０９６ビンとなる。
【００６９】
一般に、ＧＰＵ及び並列プロセッサでは、各プロセッサユニットが局所的に利用できるメモリ容量に制限がある場合が多い。これは、例えば、図３に示すＧＰＵ装置において、ローカル共有メモリ４０２の容量が制限されているためである。このため、共有メモリの容量の不足が、ＧＰＵや並列プロセッサで結合ヒストグラムを計算する場合の問題点となっていた。
【００７０】
一方、本実施の形態では、図７Ｂに示すように、参照画像ヒストグラムの座標リスト５を作成し、その後、参照画像ヒストグラムのビンごとに浮動画像の１次元ヒストグラム８を生成して、結合ヒストグラムを生成する。このように、各プロセッサユニットが２次元の結合ヒストグラムを保持する必要がなく、１次元のヒストグラムを保持するだけでよい計算方法を提供することによって、前述した問題を解決することができる。
【００７１】
相互情報量の計算には、まず、参照画像１と浮動画像２の同じ位置の画素値から生成される２次元ヒストグラムを計算する。その後、計算された２次元ヒストグラムを用いて計算される確率密度関数を用いて、相互情報量を計算することができる。
【００７２】
相互情報量とは、２つの事象ＡとＢについて、事象Ａが有している事象Ｂに関する情報量を定量化した尺度である。レジストレーションで用いられる正規化相互情報量ＮＭＩ（Ａ，Ｂ）は事象Ａと事象Ｂの２次元の結合ヒストグラムＨｉｓｔ（Ａ，Ｂ）から下式（４）によって求められる。
【００７３】
【数４】
【００７４】
ここで、Ｈ（Ａ）は事象Ａのエントロピー、Ｈ（Ｂ）は事象Ｂのエントロピー、Ｈ（Ａ，Ｂ）は事象Ａ，Ｂの結合エントロピーである。ｐ（ａ）はａの確率密度分布、ｐ（ｂ）はｂの確率密度分布、ｐ（ａ，ｂ）はａ，ｂの同時確率分布であり、式（４）のように、Ｈｉｓｔ（Ａ，Ｂ）から求められる。
【００７５】
事象Ａと事象Ｂが完全に独立である場合には、Ｈ（Ａ，Ｂ）＝Ｈ（Ａ）＋Ｈ（Ｂ）となる。また、事象Ａと事象Ｂが完全に従属である場合には、Ｈ（Ａ，Ｂ）＝Ｈ（Ａ）＝Ｈ（Ｂ）となる。以上から，ＮＭＩ（Ａ，Ｂ）のとりうる範囲は、１以上２以下となる。
【００７６】
図８は、図４のフローチャートにおける参照画像ヒストグラム作成（Ｓ１０２）の処理手順の詳細を示すフローチャートである。
【００７７】
まず、変数Ｚを１から参照画像のＺ軸方向の画素数Ｒｚまで増加させながら、ステップＳ２０２からＳ２０８の処理を繰り返し実行する（Ｓ２０１）。次に、変数Ｙを１から参照画像のＹ軸方向の画素数Ｒｙまで増加させながら、ステップＳ２０３からＳ２０７の処理を繰り返し実行する（Ｓ２０２）。さらに、変数Ｘを１から参照画像のＸ軸方向の画素数Ｒｘまで増加させながら、ステップＳ２０４からＳ２０６の処理を繰り返し実行する（Ｓ２０３）。以上の処理によって、参照画像のすべての画素（Ｘ，Ｙ，Ｚ）について、ステップＳ２０４からＳ２０６の処理が実行される。
【００７８】
ステップＳ２０４では、参照画像の画素値Ｒ（Ｘ，Ｙ，Ｚ）を変数Ｉに代入する。次に、画素値Ｉを量子化して、量子化された画素値Ｉを用いて参照画像のヒストグラムのビンを求める（Ｓ２０５）。この処理は、例えば、参照画像の輝度値の最大値と最小値との間をＮｒ個の領域に分割し、輝度値が属する領域の番号からビンを求める。
【００７９】
その後、ヒストグラムのビンと、座標値（Ｘ，Ｙ，Ｚ）とを用いて、座標リスト群５にデータを格納する（Ｓ２０６）。座標リスト群５へのデータの格納方法の詳細は、図９を用いて後述する。
【００８０】
なお、前述した説明では、参照画像の全画素について、ヒストグラム及びヒストグラムリストを求めたが、ダウンサンプリングを行って少ないサンプル数の画素についてヒストグラム及びヒストグラムリストを求めてもよい。また、スーパーサンプリングを行って、サンプル数を増やしてもよい。さらに、等間隔にサンプリングするのではなく、ランダムにサンプリングしてもよい。
【００８１】
図９は、本実施の形態の座標リスト群５の詳細を説明する図である。
【００８２】
ヒストグラムリスト５は、Ｎｒ個の可変長のリストを含む。各リストの構成要素は、画素の座標である。
【００８３】
図８のステップＳ２０６において、ビン番号ｂに座標（Ｘ，Ｙ，Ｚ）を登録する場合、ｂ番目のリストの最後尾に座標（Ｘ，Ｙ，Ｚ）を追加する。
【００８４】
図９に示すように、参照画像１のすべての画素に対して、座標（Ｘ，Ｙ，Ｚ）を追加した場合の、ｂ番目のリストの長さをＢ（ｂ）とする。Ｂ（ｂ）は参照画像の１次元のヒストグラムと等価である。
【００８５】
図１０は、図４のフローチャートにおける相互情報量計算（Ｓ１０５）の処理手順の詳細を示すフローチャートである。
【００８６】
まず、変数ｂを１から参照画像のヒストグラムのビン数Ｎｒまで増加させながら、ステップＳ３０２からステップＳ３０８を繰り返し実行する（Ｓ３０１）。
【００８７】
次に、変数ｉを１からヒストグラムの長さＢ（ｂ）まで増加させながら、ステップＳ３０３からステップＳ３０７を繰り返し実行する（Ｓ３０２）。
【００８８】
ループ内では、座標リスト群５を参照して、対応する座標（Ｘ，Ｙ，Ｚ）を取り出し（Ｓ３０３）、幾何変換情報３を用いて座標（Ｘ，Ｙ，Ｚ）を座標変換し、対応する浮動画像の座標系に変換する（Ｓ３０４）。
【００８９】
次に、浮動画像の対応する画素値Ｒ（Ｘ，Ｙ，Ｚ）を変数Ｉに代入し（Ｓ３０５）、画素値Ｉから浮動画像のヒストグラムのビンを求める。この処理は、浮動画像の最大値と最小値とをＮｆ個の領域に分割し、輝度値が属する領域の番号からビンを求めることで実現できる（Ｓ３０６）。
【００９０】
その後、量子化された値に対応するビン番号ｂ’の浮動画像ヒストグラム８を更新する（Ｓ３０７）。処理の詳細は後述する。
【００９１】
ステップＳ３０１からステップＳ３０９のループによって、参照画像の全画素に対応する参照画像ヒストグラムのビンごとの浮動画像のヒストグラム８が求められる。ヒストグラム８は、図１１に示すように、ビン数Ｎｆの１次元ヒストグラムをＮｒ個含む。
【００９２】
その後、参照画像ヒストグラムのビンごとの浮動画像のヒストグラム８を合成して、２次元の結合ヒストグラム１０を求める（Ｓ３１０）。結合ヒストグラム１０は、図１２に示すように、ヒストグラム８に含まれる一次元ヒストグラムを結合したもので、Ｎｒ×Ｎｆのマトリックスで構成されている。
【００９３】
結合ヒストグラム１０の具体例を図１４Ａ、図１４Ｂに示す。図１４Ａに示すように、参照画像１と浮動画像２との位置を合わせない場合、結合ヒストグラム１０には複数のピークが現れる。しかし、図１４Ｂに示すように、参照画像１と浮動画像２との位置が合っている状態では、結合ヒストグラム１０には一つのピークが現れる。
【００９４】
その後、結合ヒストグラム１０から相互情報量１２を求める（Ｓ３１１）。
【００９５】
次に、図１０のステップＳ３０７の処理の詳細について説明する。
【００９６】
ステップＳ３０７では、参照画像のビンがｂ、浮動画像のビンがｂ’である場合、ｂ番目のヒストグラムのｂ’のビンの頻度に１を加算する。すなわち、２次元の結合ヒストグラム１０ではなく、１次元のヒストグラム８が更新される。
【００９７】
このことから、ステップＳ３０２からのループであるステップＳ３０３からステップＳ３０８の処理をＧＰＵ又は並列プロセッサで実行する場合、各プロセッシングユニット４０１は２次元のヒストグラムをローカル共有メモリ４０２に保持する必要がなく、１次元のヒストグラムをローカル共有メモリ４０２に保持するだけでよい。
【００９８】
さらに、ステップＳ３０１からのループであるステップＳ３０２からステップ３０９の処理をＧＰＵ又は並列プロセッサで実行することも可能である。この場合にも、各プロセッシングユニット４０１は２次元のヒストグラムをローカル共有メモリ４０２に保持する必要がなく、１次元のヒストグラムをローカル共有メモリ４０２に保持するだけでよい。
【００９９】
３次元画像のレジストレーション実験によれば、本発明は非特許文献２の方法と比較して、より高速である。具体的な実験は、ファントムのＣＴ画像を用いてレジストレーションの性能を評価した。使用したファントム画像は５１２×５１２×２００ボクセルの３次元画像である。ファントム画像に対して、Ｘ軸、Ｙ軸及びＺ軸に対して１０度回転したボリュームデータを作成しレジストレーションを行った。レジストレーション実験によればもとの位置に正しく座標変換された。
【０１００】
レジストレーションに使用するサンプル数を種々に変更して計算時間を計測した。サンプル数が６４×６４×６４の場合の計算時間は、ＣＰＵ実装が１４．２秒、非特許文献２の方法が１．８２秒、本発明による方法が１．１６秒であった。サンプル数が２５６×２５６×２５６の場合にはＣＰＵ実装の計算時間は７８０秒、非特許文献２の方法が１４．９秒、本発明による方法が５．７秒であった。以上のように本発明による方法はＣＰＵ実装と比較して数１０倍以上、非特許文献２の方法と比較して１．５倍から２．６倍程度の高速化が実現されている。
【０１０１】
図１３は、本実施の形態の画像処理装置で実行される処理を示すシーケンス図である。
【０１０２】
まず、レジストレーション命令が発行されると、入力装置２７は、ユーザからの参照画像及び浮動画像の指定を受けて、画像の指定を入力制御部２８を介してＣＰＵ４０に伝える（１３０１）。ＣＰＵ４０は、指定された参照画像及び浮動画像を取得し、取得した参照画像及び浮動画像を記憶装置２３に格納する（１３０２）。ＧＰＵ４１は、記憶装置２３に格納された参照画像を取得し、参照画像ヒストグラムを生成し（１３０３）、生成された参照画像ヒストグラムを記憶装置に格納する。
【０１０３】
さらに、入力装置２７は、ユーザからの幾何変換情報の指定を受けて、入力された幾何変換情報を入力制御部２８を介してＣＰＵ４０に伝える（１３０４）。この幾何変換情報は、幾何変換の初期条件を含み、さらに、幾何変換情報を更新する条件を含んでもよい。ＣＰＵ４０は、入力された幾何変換情報をＧＰＵ４１に転送する（１３０５）。
【０１０４】
ＧＰＵ４１は、幾何変換を受けると、浮動画像を幾何変換し、幾何変換された浮動画像の画素値をサンプリングし、参照画像ヒストグラムのビンごとに浮動画像の１次元ヒストグラム８を生成する（１３０６）。このため、ＧＰＵ４１は、浮動画像の１次元ヒストグラム８をローカル共有メモリ４０２に格納すればよく。小容量で高速アクセス可能なメモリのみを用いて（すなわち、大容量だがアクセス速度が低いグローバルメモリ４０３を用いることなく）、浮動画像の１次元ヒストグラム８を生成することができる。その後、ＧＰＵ４１は、生成された浮動画像の１次元ヒストグラム８を、記憶装置２３に格納する。
【０１０５】
ＧＰＵ４１は、全ての１次元ヒストグラムの計算を終了した後、記憶装置２３から１次元ヒストグラム８を読み出して、読み出された１次元ヒストグラム８を結合する。そして結合されたヒストグラム１０から相互情報量を計算する（１３０７）。計算された相互情報量はＣＰＵ４０に転送される（１３０８）。
【０１０６】
ＣＰＵ４０は、相互情報量が最大となるように幾何変換情報を更新し（１３０９）、更新された幾何変換情報をＧＰＵに転送する（１３１０）。
【０１０７】
ＧＰＵ４１は、転送された幾何変換情報に基づいて、１３０６〜１３０９の処理と同様に、幾何変換された浮動画像の１次元ヒストグラム８を生成し、結合ヒストグラムを生成し、相互情報量を計算し（１３１１）、計算された相互情報量をＣＰＵ４０に転送する（１３１２）。
【０１０８】
ＣＰＵ４０は、相互情報量が最大となる幾何変換情報で浮動画像を変換し（１３１３）、この条件で幾何変換された浮動画像（レジストレーション済画像）をディスプレイ３１に表示する（１３１４）。
【０１０９】
以上に説明したように、本発明の実施の形態では、２次元の結合ヒストグラムの計算をＧＰＵ又は並列プロセッサで実行する場合に、各プロセッシングユニットは１次元のヒストグラムにアクセスするので、必要なメモリ容量を削減することができる。このため、結合ヒストグラムを効率よく計算でき、レジストレーション計算を効率的に処理することができる。
【符号の説明】
【０１１０】
１ 参照画像
２ 浮動画像
３ 幾何変換情報
５ 参照画像ヒストグラム（座標リスト群）
８ 浮動画像ヒストグラム
１０ 結合ヒストグラム
１２ 相互情報量
１５ レジストレーション済み浮動画像
１６ 前処理部
１７ レジストレーション部
４１ ＧＰＵ
【特許請求の範囲】
【請求項１】
２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせる画像処理システムであって、
前記画像処理システムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、
前記画像処理システムは、
前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手段と、
所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２の手段と、
前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手段と、
前記生成された結合ヒストグラムから相互情報量を計算する第４の手段と、を備えることを特徴とする画像処理システム。
【請求項２】
前記画像処理システムは、並列して演算処理が可能な複数の前記プロセッシングエレメントと、前記ローカルメモリと、前記グローバルメモリとを備えるメニーコアプロセッサシステムによって構成され、
前記各プロセシングエレメントが所定のプログラムを実行することによって、前記ローカルメモリを用いて、参照画像のヒストグラムのビン毎に並列して、前記第２の手段を実現することを特徴とする請求項１に記載の画像処理システム。
【請求項３】
前記メニーコアプロセッサシステムは、グラフィックプロセッシングユニットであることを特徴とする請求項２に記載の画像処理システム。
【請求項４】
前記第２の手段は、前記生成された参照画像のヒストグラムの一つのビンに属する画素の位置を取得し、前記取得された画素の位置に対応する前記浮動画像の画素の画素値を取得し、前記取得された浮動画像の画素の画素値を量子化し、前記量子化された画素値に画素数を積算することによって、浮動画像の１次元ヒストグラムを生成する処理を、前記参照画像のヒストグラムのビン数だけ実行することを特徴とする請求項１に記載の画像処理システム。
【請求項５】
前記第１の手段は、前記参照画像の画素を表すパラメータを変化させ、当該パラメータによって特定される画素の画素値を量子化し、前記量子化された画素値のビンに前記特定された画素の位置を記録することによって、前記参照画像のヒストグラムを生成することを特徴とする請求項１に記載の画像処理システム。
【請求項６】
前記第２の手段は、
前記生成された参照画像のヒストグラムの一つのビンに属する画素に対応する位置の前記変換された浮動画像の画素値を取得し、前記取得された浮動画像の画素値を量子化し、前記量子化された画素値を前記浮動画像のヒストグラムに加える処理を、前記参照画像のヒストグラムの一つのビンに属する画素の数だけ実行し、
さらに、前記の処理を、前記生成された１次元ヒストグラムの一つのビンの数だけ実行することを特徴とする請求項１に記載の画像処理システム。
【請求項７】
２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせるためのパラメータを計算する画像処理方法であって、
前記画像処理方法が実行されるシステムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、
前記画像処理方法は、
前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１のステップと、
前記各プロセッシングエレメントが、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２のステップと、
前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３のステップと、
前記生成された結合ヒストグラムから相互情報量を計算する第４のステップと、を含むことを特徴とする画像処理方法。
【請求項８】
前記画像処理システムは、複数の前記プロセッシングエレメントと、前記ローカルメモリと、前記グローバルメモリとを備え、各プロセシングエレメントが並列して演算処理が可能なメニーコアプロセッサシステムによって構成され、
前記画像処理システムは、並列して演算処理が可能な複数の前記プロセッシングエレメントと、前記ローカルメモリと、前記グローバルメモリとを備えるメニーコアプロセッサシステムによって構成され、
前記各プロセシングエレメントは、前記ローカルメモリを用いて、参照画像のヒストグラムのビン毎に並列して、前記第２のステップの処理を実行することを特徴とする請求項７に記載の画像処理方法。
【請求項９】
前記第２ステップでは、前記生成された参照画像のヒストグラムの一つのビンに属する画素の位置を取得し、前記取得された画素の位置に対応する前記浮動画像の画素の画素値を取得し、前記取得された浮動画像の画素の画素値を量子化し、前記量子化された画素値に画素数を積算することによって、浮動画像の１次元ヒストグラムを生成する処理を、前記参照画像のヒストグラムのビン数だけ実行することを特徴とする請求項７に記載の画像処理方法。
【請求項１０】
前記第１ステップでは、前記参照画像の画素を表すパラメータを変化させ、当該パラメータによって特定される画素の画素値を量子化し、前記量子化された画素値のビンに前記特定された画素の位置を記録することによって、前記参照画像のヒストグラムを生成することを特徴とする請求項７に記載の画像処理方法。
【請求項１１】
前記第２ステップでは、
前記生成された参照画像のヒストグラムの一つのビンに属する画素に対応する位置の前記変換された浮動画像の画素値を取得し、前記取得された浮動画像の画素値を量子化し、前記量子化された画素値を前記浮動画像のヒストグラムに加える処理を、前記参照画像のヒストグラムの一つのビンに属する画素の数だけ実行し、
さらに、前記の処理を、前記生成された１次元ヒストグラムの一つのビンの数だけ実行することを特徴とする請求項７に記載の画像処理方法。
【請求項１２】
２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせるためのパラメータを画像処理システムに計算させるための画像処理プログラムであって、
前記画像処理方法が実行されるシステムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、
前記プログラムは、
前記各プロセッシングエレメントが、前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手順と、
前記各プロセッシングエレメントが、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２の手順と、
前記各プロセッシングエレメントが、前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手順と、
前記各プロセッシングエレメントが、前記生成された結合ヒストグラムから相互情報量を計算する第４の手順と、を含むことを特徴とする画像処理プログラム。
【請求項１】
２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせる画像処理システムであって、
前記画像処理システムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、
前記画像処理システムは、
前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手段と、
所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２の手段と、
前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手段と、
前記生成された結合ヒストグラムから相互情報量を計算する第４の手段と、を備えることを特徴とする画像処理システム。
【請求項２】
前記画像処理システムは、並列して演算処理が可能な複数の前記プロセッシングエレメントと、前記ローカルメモリと、前記グローバルメモリとを備えるメニーコアプロセッサシステムによって構成され、
前記各プロセシングエレメントが所定のプログラムを実行することによって、前記ローカルメモリを用いて、参照画像のヒストグラムのビン毎に並列して、前記第２の手段を実現することを特徴とする請求項１に記載の画像処理システム。
【請求項３】
前記メニーコアプロセッサシステムは、グラフィックプロセッシングユニットであることを特徴とする請求項２に記載の画像処理システム。
【請求項４】
前記第２の手段は、前記生成された参照画像のヒストグラムの一つのビンに属する画素の位置を取得し、前記取得された画素の位置に対応する前記浮動画像の画素の画素値を取得し、前記取得された浮動画像の画素の画素値を量子化し、前記量子化された画素値に画素数を積算することによって、浮動画像の１次元ヒストグラムを生成する処理を、前記参照画像のヒストグラムのビン数だけ実行することを特徴とする請求項１に記載の画像処理システム。
【請求項５】
前記第１の手段は、前記参照画像の画素を表すパラメータを変化させ、当該パラメータによって特定される画素の画素値を量子化し、前記量子化された画素値のビンに前記特定された画素の位置を記録することによって、前記参照画像のヒストグラムを生成することを特徴とする請求項１に記載の画像処理システム。
【請求項６】
前記第２の手段は、
前記生成された参照画像のヒストグラムの一つのビンに属する画素に対応する位置の前記変換された浮動画像の画素値を取得し、前記取得された浮動画像の画素値を量子化し、前記量子化された画素値を前記浮動画像のヒストグラムに加える処理を、前記参照画像のヒストグラムの一つのビンに属する画素の数だけ実行し、
さらに、前記の処理を、前記生成された１次元ヒストグラムの一つのビンの数だけ実行することを特徴とする請求項１に記載の画像処理システム。
【請求項７】
２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせるためのパラメータを計算する画像処理方法であって、
前記画像処理方法が実行されるシステムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、
前記画像処理方法は、
前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１のステップと、
前記各プロセッシングエレメントが、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２のステップと、
前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３のステップと、
前記生成された結合ヒストグラムから相互情報量を計算する第４のステップと、を含むことを特徴とする画像処理方法。
【請求項８】
前記画像処理システムは、複数の前記プロセッシングエレメントと、前記ローカルメモリと、前記グローバルメモリとを備え、各プロセシングエレメントが並列して演算処理が可能なメニーコアプロセッサシステムによって構成され、
前記画像処理システムは、並列して演算処理が可能な複数の前記プロセッシングエレメントと、前記ローカルメモリと、前記グローバルメモリとを備えるメニーコアプロセッサシステムによって構成され、
前記各プロセシングエレメントは、前記ローカルメモリを用いて、参照画像のヒストグラムのビン毎に並列して、前記第２のステップの処理を実行することを特徴とする請求項７に記載の画像処理方法。
【請求項９】
前記第２ステップでは、前記生成された参照画像のヒストグラムの一つのビンに属する画素の位置を取得し、前記取得された画素の位置に対応する前記浮動画像の画素の画素値を取得し、前記取得された浮動画像の画素の画素値を量子化し、前記量子化された画素値に画素数を積算することによって、浮動画像の１次元ヒストグラムを生成する処理を、前記参照画像のヒストグラムのビン数だけ実行することを特徴とする請求項７に記載の画像処理方法。
【請求項１０】
前記第１ステップでは、前記参照画像の画素を表すパラメータを変化させ、当該パラメータによって特定される画素の画素値を量子化し、前記量子化された画素値のビンに前記特定された画素の位置を記録することによって、前記参照画像のヒストグラムを生成することを特徴とする請求項７に記載の画像処理方法。
【請求項１１】
前記第２ステップでは、
前記生成された参照画像のヒストグラムの一つのビンに属する画素に対応する位置の前記変換された浮動画像の画素値を取得し、前記取得された浮動画像の画素値を量子化し、前記量子化された画素値を前記浮動画像のヒストグラムに加える処理を、前記参照画像のヒストグラムの一つのビンに属する画素の数だけ実行し、
さらに、前記の処理を、前記生成された１次元ヒストグラムの一つのビンの数だけ実行することを特徴とする請求項７に記載の画像処理方法。
【請求項１２】
２次元又は３次元画像の参照画像と浮動画像との間で相互情報量が大きくなる条件で、両画像の位置を合わせるためのパラメータを画像処理システムに計算させるための画像処理プログラムであって、
前記画像処理方法が実行されるシステムは、１以上のプロセッシングエレメントと、前記各プロセシングエレメントから高速にアクセス可能なローカルメモリと、前記各プロセシングエレメントに共有されるデータを格納するグローバルメモリと、を備え、
前記プログラムは、
前記各プロセッシングエレメントが、前記参照画像の画素値の１次元のヒストグラムを生成し、前記生成された参照画像のヒストグラムの各ビンに属する画素の位置を記憶装置に格納する第１の手順と、
前記各プロセッシングエレメントが、所定の幾何変換情報に従って浮動画像の座標を変換し、前記生成された参照画像のヒストグラムの一つのビンに属する参照画像の画素の位置に対応する前記変換された浮動画像の画素の画素値の頻度を示す、浮動画像の１次元ヒストグラムを、前記ローカルメモリ上に生成する処理を実行する第２の手順と、
前記各プロセッシングエレメントが、前記生成された浮動画像の１次元ヒストグラムを合成することによって、前記浮動画像の２次元の結合ヒストグラムを生成する第３の手順と、
前記各プロセッシングエレメントが、前記生成された結合ヒストグラムから相互情報量を計算する第４の手順と、を含むことを特徴とする画像処理プログラム。
【図１】
【図２】
【図３】
【図４】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１５】
【図５Ａ】
【図５Ｂ】
【図６Ａ】
【図６Ｂ】
【図７Ａ】
【図７Ｂ】
【図１４Ａ】
【図１４Ｂ】
【図２】
【図３】
【図４】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１５】
【図５Ａ】
【図５Ｂ】
【図６Ａ】
【図６Ｂ】
【図７Ａ】
【図７Ｂ】
【図１４Ａ】
【図１４Ｂ】
【公開番号】特開２０１１−２３９８１２（Ｐ２０１１−２３９８１２Ａ）
【公開日】平成２３年１２月１日（２０１１．１２．１）
【国際特許分類】
【出願番号】特願２０１０−１１１８８２（Ｐ２０１０−１１１８８２）
【出願日】平成２２年５月１４日（２０１０．５．１４）
【出願人】（０００００５１０８）株式会社日立製作所 (27,607)
【Ｆターム（参考）】
【公開日】平成２３年１２月１日（２０１１．１２．１）
【国際特許分類】
【出願日】平成２２年５月１４日（２０１０．５．１４）
【出願人】（０００００５１０８）株式会社日立製作所 (27,607)
【Ｆターム（参考）】
[ Back to top ]