ボリュームデータ間の対応付け方法
【課題】処理時間を短縮し、かつ精度よくボリュームデータ間の位置合わせを行う。
【解決手段】参照ボリュームデータ（Ｄ１）から探索ウィンドウ（３次元ブロック）を複数抽出する。抽出した複数の探索ウィンドウに設定される基準点に対応する入力ボリュームデータ（Ｄ２）の探索ウィンドウ内の対応点を、各探索ウィンドウより抽出したボクセルデータの相関（３次元ＰＯＣ）から探索する。探索した対応点と基準点との対応関係から参照ボリュームデータと入力ボリュームデータとの間の剛体変形パラメータ（回転ずれ、位置ずれを示すパラメータ）を推定し、入力ボリュームデータを剛体変形（回転ずれ、位置ずれを補正）する。剛体変形処理後の入力ボリュームデータと参照ボリュームデータについて、同様にして対応点を探索し、探索した対応点の情報を用いてＦＦＤ（Free-Form Deformation）により、入力ボリュームデータをさらに非剛体変形させる。参照ボリュームデータを剛体・非剛体変形させてもよい。
【解決手段】参照ボリュームデータ（Ｄ１）から探索ウィンドウ（３次元ブロック）を複数抽出する。抽出した複数の探索ウィンドウに設定される基準点に対応する入力ボリュームデータ（Ｄ２）の探索ウィンドウ内の対応点を、各探索ウィンドウより抽出したボクセルデータの相関（３次元ＰＯＣ）から探索する。探索した対応点と基準点との対応関係から参照ボリュームデータと入力ボリュームデータとの間の剛体変形パラメータ（回転ずれ、位置ずれを示すパラメータ）を推定し、入力ボリュームデータを剛体変形（回転ずれ、位置ずれを補正）する。剛体変形処理後の入力ボリュームデータと参照ボリュームデータについて、同様にして対応点を探索し、探索した対応点の情報を用いてＦＦＤ（Free-Form Deformation）により、入力ボリュームデータをさらに非剛体変形させる。参照ボリュームデータを剛体・非剛体変形させてもよい。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、同一対象物を撮像した２つの３次元の立体画像のデータ（ボリュームデータ）間の対応付けを行うボリュームデータ間の対応付け方法に関するものである。
【背景技術】
【０００２】
ボリュームデータとは、３次元の小立方体（ボクセル）を単位として、各ボクセルにその部分の透過率などのデータを持たせ、物体の内部構造まで表現したデータのことであり、医療や流体解析、コンピュータグラフィックスなど広範な分野において頻繁に用いられるデータ形式である。
【０００３】
例えば医療分野においては、Ｘ線ＣＴ (Computed Tomography) やＭＲＩ (Magnetic Resonance Imaging) から得られるデータが挙げられる。これらの医用画像撮影装置は、人体内部のスライス像を連続的に取得することが可能であり、取得した２次元画像群を積層することで３次元ボリュームデータとして人体内部を可視化することができる。その有用性はすでに広く知られており、日常臨床でも一般的に用いられている。
【０００４】
ボリュームデータを扱う上では、複数のボリュームデータの比較や合成が要求される場面が頻繁に現れる。特に医療分野においては、異なる時間に撮影した病変部位を観察することで、病状の進行や解剖学的な変化について手掛かりを得ることができる。また、異なる装置を用いて撮影したボリュームデータをフュージョン（重ね合わせ）することで、それぞれの装置が苦手とする情報を補完した観察が可能となる。他にも、例えば多列検出器型ＣＴを用いて心臓などのボリュームデータを高速かつ連続的に取得し、時間方向にデータを比較することで心臓の動きを可視化することも可能である。
【０００５】
ここで、撮影のタイミングが異なるボリュームデータ間には一般に未知の位置ずれやひずみが存在するため、正確に比較や合成を行うためにはまず２つのボリュームデータ間のずれやひずみを表すパラメータを知る必要がある。これは、２つのボリュームデータの位置合わせ（レジストレーション）問題として、これまでにランドマークベース（例えば、非特許文献１参照）や幾何学的形状ベース（例えば、非特許文献２参照）、ボクセル類似度ベース（例えば、非特許文献３，４，５，６，７，特許文献１参照）の手法が提案されている。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２００７−５４６３６号公報
【非特許文献】
【０００７】
【非特許文献１】C. R. Maurer Jr., J. M. Fitzpatrick, M. Y. Wang,R. L. Galloeay Jr., R. J. Maciunas, and G. S.Allen. Registration of head volume images usingimplantable fiducial markers. IEEE Trans. MedicalImaging, Vol. 16, No. 4, pp. 447--462, Aug 1997.
【非特許文献２】P. J. Besl and N. D. McKay. A method for registrationof 3-D shapes. IEEE Trans. Pattern Analysisand Machine Intelligence, Vol. 14, No. 2, pp. 239--256,Feb 1992.
【非特許文献３】A. Roche, X. Pennec, G. Malandain, and N. Ayache.Rigid registration of 3-D ultrasound with MR images:A new apploach combining intensity and gradient information.IEEE Trans. Medical Imaging, Vol. 20,No. 10, pp. 1038--1049, Oct 2001.
【非特許文献４】F. Maes, A. Collignon, D. Vandermeulen, G. Marchal,and P. Suetens. Multimodality image registrationby maximization of mutual information. IEEETrans. Medical Imaging, Vol. 16, No. 2, pp. 187--198,Apr 1997.
【非特許文献５】C. Studholme, D. L. G. Hill, and D. J. Hawkes. Aninvariant entropy measure of 3D medical image alignment.Pattern Recognition, Vol. 32, No. 1, pp. 71--86,Apr 1998.
【非特許文献６】D. L. G. Hill, P. G. Batchelor, M. Holden, andD. J. Hawks. Medical image registration. Physicsin Medicine and Biology, Vol. 46, pp. R1--R45, 2001.
【非特許文献７】D. Rueckert, L. I. Sonoda, C. Hayes, D. L. G. Hill,M. O. Leach, and D. J. Hawkes. Nonrigid registrationusing free-form deformation: Application to breastMR images. IEEE Trans. Medical Imaging, Vol. 18,No. 8, pp. 712--721, Aug 1999.
【非特許文献８】K. Takita, T.Aoki, Y. Sasaki, T. Higuchi, andK. Kobayashi. High-accuracy subpixel image registrationbased on phase-only correlation. IEICETrans. Fundamentals, Vol. E86-A, No. 8, pp. 1925--1934, Aug 2003.
【非特許文献９】K. Takita, M. A. Muquit, T. Aoki, and T. Higuchi.A sub-pixel correspondence search technique for computervision applications. IEICE Trans. Fundamentals,Vol. E28-A, No. 8, pp. 1913--1923, 2004.
【非特許文献１０】S. Nagashima, T. Aoki, T. Higuchi, andK. Kobayashi. A subpixel image matching techniqueusing phase-only correlation. Proc. IEEE2006 Int. Symp. Intelligent Signal Processing andCommunication Systems, pp. 701--704, Dec 2006.
【非特許文献１１】K. S. Arun, T. S. Haung, and S. D. Blostein.Least-squares fitting of two 3-D point sets. IEEETrans. Pattern Analysis and Machine Intelligence,Vol. PAMI-9, No. 5, pp. 698--700, Sep 1987.
【非特許文献１２】S. Lee, G. Wolberg, and S. Y. Shin. Scattered datainterpolation with multilevel B-splines. IEEE Trans.Visualization and Computer Graphics, Vol. 3, No. 3,pp. 228--244, Jul 1997.
【非特許文献１３】D. R. Forsey and R. H. Bartels. Hierarchical B-splinerefinement. ACM Trans. Computer Graphics, Vol. 22,No. 4, pp. 205--212, Aug 1988.
【非特許文献１４】M. A. Fischler and R. C. Bolles. Random sampleconsensus: A paradigm for model fitting with applicationsto image analysis and automated cartgraphy.Comm. of the ACM, Vol. 24, No. 6, pp. 381--395, Jun1981.
【非特許文献１５】J. A. Nelder and R. Mead. A simplex method forfunction minimization. Computer Journal, Vol. 7,No. 4, pp. 308--313, 1965.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、ランドマークベースの手法では、埋め込み型マーカなどの位置情報を用いてレジストレーションを行うため、撮影時に侵襲的なマーカが必要であるという問題がある。
【０００９】
幾何学的形状ベースの手法としては、ＩＣＰ（Iterative Closest Point）アルゴリズムが知られている。これは、３次元形状モデル（ポリゴンデータなど）の位置合わせに広く用いられている手法であり、モデル間の対応付けと運動推定を繰り返すことで位置合わせを行う。そのため、まずボリュームから形状モデルを生成する必要があるが、２つのデータ間でモダリティ（様相）が異なるような場合は同一の形状モデルを得ることが難しく、マルチモーダルレジストレーションに適さない。
【００１０】
なお、異なる撮像機器から得られたデータをレジストレーションすることをマルチモーダルレジストレーションと呼んでいる。例えば医療分野では、ＣＴから得られたデータとＭＲＩから得られたデータをレジストレーションするような場合がこれにあたる。また、同じ撮像機器の場合でも、撮影方法が著しく異なるようなデータの場合は、マルチモーダルレジストレーションと呼ぶ場合がある。例えば、ＭＲＩからはＴ１強調画像とＴ２強調画像など、撮影方法によってコントラストが大きく異なった画像が得られる。
【００１１】
ボクセル類似度ベースの手法は、位置合わせパラメータをボリュームデータ間の類似度が最大となるように非線形最適化によって決定する手法であり、医用データのレジストレーションにおいて現在最もよく用いられている。なかでも、非特許文献６に示されているような正規化された相互情報量（Normalized Mutual Information：ＮＭＩ）を類似度として用いる手法が広く知られており、この手法ではボリュームデータ間で計算した相互情報量を最大化することでレジストレーションを行う。なお、特許文献１には、２つのボリュームデータから特徴点を抽出し、この特徴点の近傍のみで相互情報量を最大化することでレジストレーションを行うようにすることが開示されている。
【００１２】
しかし、このような相互情報量の最大化は非線形最適化問題として解かれるため、その膨大な計算時間を要する。特に、３次元のボリュームデータは２次元画像などに比べて大幅に大きなデータ量を持つため、データの解像度によっては位置合わせに数時間を要する場合もある。また、非線形最適化において正しい解を得るためには適切な初期値を設定しなければならず、この初期値を推定する処理も必要となる。
【００１３】
本発明は、このような課題を解決するためになされたもので、その目的とするところは、ボリュームデータ間の高精度な対応付け方法を提供することにある。２つのボリュームデータの間でボクセル同士の対応関係を知ることは、レジストレーションや対象の動きを可視化する際にきわめて有用である。レジストレーションにおいては、ボクセル同士の対応関係からレジストレーションのためのパラメータを推定することができる。本発明においてボリュームデータの対応付け方法に併せて提供されるボリュームデータの位置合わせ方法は、従来技術から精度を低下させることなく大幅な処理時間の削減を達成している。
【課題を解決するための手段】
【００１４】
このような目的を達成するために本発明は、同一対象物を撮像した２つの３次元の立体画像のデータを第１および第２のボリュームデータとして取り込むデータ取込ステップと、第１のボリュームデータから複数の部分領域を抽出し、この抽出した複数の部分領域にそれぞれ設定された基準点に対応する第２のボリュームデータの対応点を、基準点を含む第１のボリュームデータの部分領域と基準点に対応する対応点を含む第２のボリュームデータの部分領域との相関から探索する対応点探索ステップと、探索された対応点と基準点との対応関係から第１のボリュームデータと第２のボリュームデータとの間の回転ずれおよび位置ずれを表すパラメータを剛体変形パラメータとして推定する剛体変形パラメータ推定ステップとを設けたものである。
【００１５】
本発明では、第１のボリュームデータから部分領域が抽出され、この抽出された部分領域に設定された基準点に対応する第２のボリュームデータの部分領域内の対応点が、その部分領域の相関から探索される。同様にして、第１のボリュームデータから複数の部分領域が抽出され、その抽出された部分領域内の基準点に対応する第２のボリュームデータの部分領域内の対応点が探索される。そして、この探索された対応点と基準点との対応関係から第１のボリュームデータと第２のボリュームデータとの間の剛体変形パラメータ（回転ずれおよび位置ずれを表すパラメータ）が推定される。
【００１６】
本発明では、部分領域の相関から対応点を探索するが、この場合の相関法として位相限定相関法（Phase-Only Correlation：ＰＯＣ）を用いることが考えられる。ＰＯＣは、信号を離散フーリエ変換して得られる位相成分のみに着目して信号間の相関を求める手法であり、さまざまな外乱に対してロバストに信号の位置合わせ、類似度評価を行うことができるという特徴がある（例えば、非特許文献８参照）。
【００１７】
ＰＯＣ関数が示す相関ピークは、その理論モデルを入力信号によらず導出することが可能であり、これを用いたフィッティングによって信号間のずれをサンプリング分解能以上の精度で求めることができる。ＰＯＣは、本出願人らによってさまざまな分野への応用が実現されているが、特にサブピクセル精度での画像対応付けが必要となるステレオビジョンなどにおいてその有効性が確認されている（例えば、非特許文献９参照）。
【００１８】
本発明では、その一形態として、これまで１次元および２次元で使用されていたＰＯＣを３次元へ拡張し、この３次元へ拡張したＰＯＣ（３次元ＰＯＣ）を用いた３次元ブロックマッチングによってサブボクセル精度でボリュームデータ間の対応付けを行う。そして、得られた対応関係から、ボリュームデータ間の剛体変形パラメータを推定する。
【００１９】
また、本発明では、その一形態として、推定された剛体変形パラメータに基づいて第１のボリュームデータと第２のボリュームデータとの間の回転ずれおよび位置ずれを補正する。この場合、第１のボリュームデータの回転ずれおよび位置ずれを補正するようにしてもよいし、第２のボリュームデータの回転ずれおよび位置ずれを補正するようにしてもよい。この第１のボリュームデータと第２のボリュームデータとの間の回転ずれおよび位置ずれの補正を「剛体変形」と呼ぶ。すなわち、本発明では、第１のボリュームデータを剛体変形させてもよいし、第２のボリュームデータを剛体変形させてもよい。
【００２０】
また、本発明では、その一形態として、回転ずれおよび位置ずれが補正された第１および第２のボリュームデータを剛体変形処理後の第１および第２のボリュームデータとし、剛体変形処理後の第１のボリュームデータから複数の部分領域を抽出し、この抽出した複数の部分領域にそれぞれ設定された基準点に対応する剛体変形処理後の第２のボリュームデータの対応点を、基準点を含む剛体変形処理後の第１のボリュームデータの部分領域と基準点に対応する対応点を含む剛体変形処理後の第２のボリュームデータの部分領域との相関から探索するようにし、この探索された対応点の情報を用いて剛体変形処理後の第１および第２のボリュームデータの何れか一方の非剛体部分を変形させて重ね合わせるようにしてもよい。この剛体変形処理後の第１および第２のボリュームデータの何れか一方の非剛体部分を変形させて重ね合わせるような補正を「非剛体変形」と呼ぶ。すなわち、本発明では、剛体変形処理後の第１のボリュームデータを非剛体変形させてもよいし、剛体変形処理後の第２のボリュームデータを非剛体変形させてもよい。
【００２１】
本発明は、ボリュームデータの高精度な対応付けを提供する。併せて提供するボクセルの対応関係を利用したボリュームデータの位置合わせ方法は、従来手法で一般的に用いられる複雑な最適化計算を必要としないため、位置合わせに要する計算時間を大幅に削減することができる。また、最適化計算において必須である変形パラメータの初期値を設定する必要もない。さらに、非剛体変形への対応も容易である。
【発明の効果】
【００２２】
本発明によれば、第１のボリュームデータから部分領域を複数抽出し、この抽出した複数の部分領域内の基準点に対応する第２のボリュームデータの部分領域内の対応点を各部分領域の相関から探索することで、ボリュームデータ間の高精度な対応付けが可能となる。
【００２３】
また、本発明によれば、ボリュームデータの対応関係から推定した剛体変形パラメータに基づいて第１のボリュームデータと第２のボリュームデータとの間の回転ずれおよび位置ずれを補正するようにして、複雑な最適化計算を必要とせず、処理時間を短縮し、かつ精度よくボリュームデータ間の位置合わせを行うことが可能となる。
【図面の簡単な説明】
【００２４】
【図１】本発明の実施に用いるボリュームデータ間位置合わせ装置の一例を示すブロック図である。
【図２】このボリュームデータ間位置合わせ装置で行われる特有の処理動作（剛体変形処理）を示すフローチャートである。
【図３】このボリュームデータ間位置合わせ装置における対応点の探索処理（３次元ブロックマッチング）の過程を示すフローチャートである。
【図４】図３に続くフローチャートである。
【図５】このボリュームデータ間位置合わせ装置における対応点の探索処理（３次元ブロックマッチング）の過程を説明するためのボリュームデータの階層構造を示す図である。
【図６】対応点の探索処理に際して用いる３次元ＰＯＣの流れを示す図である。
【図７】このボリュームデータ間位置合わせ装置で行われる特有の処理動作（非剛体変形処理）を示すフローチャートである。
【図８】ＰＯＣに基づく対応点探索におけるパラメータを例示する図である。
【図９】本発明に係る手法と従来手法を用いてレジストレーションを行い、それらの性能を相関係数により比較した結果を示す図である。
【図１０】実験にて使用したボリュームの一部（（ａ），（ｂ））およびレジストレーションの結果（（ｃ），（ｄ））を示すディスプレイ上の写真である。
【図１１】Ｘ線ＣＴデータおよびＭＲＩデータのアキシャル像（（ａ），（ｂ））およびレジストレーションの結果（（ｃ），（ｄ））を示すディスプレイ上の写真である。
【発明を実施するための形態】
【００２５】
以下、本発明の実施の形態を図面に基づいて詳細に説明する。図１は本発明の実施に用いるボリュームデータ間位置合わせ装置の一例を示すブロック図である。
【００２６】
同図において、１００はボリュームデータ間位置合わせ装置であり、プロセッサや記憶装置からなるハードウェアと、これらのハードウェアと協働して各種機能を実現させるプログラムとによって実現される。具体的には、パーソナルコンピュータ（パソコン）にプログラムがインストールされ、このインストールされたプログラムに従うＣＰＵの処理動作として実現される。
【００２７】
ボリュームデータ間位置合わせ装置１００は、参照ボリュームデータ取込部１と、入力ボリュームデータ取込部２と、撮像条件補正部３と、剛体変形処理部４と、非剛体変形処理部５と、ボリュームデータ記憶部６とを備えている。剛体変形処理部４は、オブジェクト領域抽出部４Ａと、対応点探索部４Ｂと、剛体変形パラメータ推定部４Ｃと、剛体変形実行部４Ｄと、反復処理命令部４Ｅとを備え、非剛体変形処理部５は、オブジェクト領域抽出部５Ａと、対応点探索部５Ｂと、非剛体変形実行部５Ｃとを備えている。
【００２８】
この実施の形態において、ボリュームデータ間位置合わせ装置１００には、一例として、参照ボリュームデータＤ１としてＣＴ（図示せず）によって撮影された対象物のスライス画像を再構成した３次元の立体画像のデータが入力され、入力ボリュームデータＤ２としてＭＲＩ（図示せず）によって撮影された対象物のスライス画像を再構成した３次元の立体画像のデータが入力されるものとする。
【００２９】
ここで、本実施の形態においては、位置合わせの基準となるボリュームデータを「参照ボリュームデータ」と呼ぶこととし、この参照ボリュームデータを基準として入力ボリュームデータの位置合わせを行うものとする。なお、この例において、撮影された対象物は同一人の頭部であるとする。
【００３０】
また、本実施の形態において、剛体変形処理とは、撮影対象物自体は変形しないものとして、２つのボリュームデータ間の移動、回転によるずれなどを修正する処理のことを意味する。また、非剛体変形処理とは、撮像対象物が例えば心臓のように、時間的に変形することによって生じる２つのボリュームデータ間のずれを修正する処理のことを意味する。
【００３１】
図２におよび図７にこのボリュームデータ間位置合わせ装置１００で行われる特有の処理動作のフローチャートを示す。以下、このフローチャートを参照しながら、ボリュームデータ間位置合わせ装置１００における各部の処理動作について説明する。
【００３２】
〔剛体変形処理〕
先ず、図２に示すフローチャートを用いて、ボリュームデータ間位置合わせ装置１００で行われる剛体変形処理について説明する。
【００３３】
〔ボリュームデータの取り込み〕
参照ボリュームデータ取込部１は、ＣＴによって撮影され再構成された３次元の立体画像のデータを参照ボリュームデータＤ１として取り込み、この取り込んだ参照ボリュームデータＤ１を撮像条件補正部３へ送る（ステップＳ１０１）。
【００３４】
入力ボリュームデータ取込部２は、ＭＲＩによって撮影され再構成された３次元の立体画像のデータを入力ボリュームデータＤ２として取り込み、この取り込んだ入力ボリュームデータＤ２を撮像条件補正部３へ送る（ステップＳ１０２）。
【００３５】
〔撮像条件の補正〕
撮像条件補正部３は、参照ボリュームデータ取込部１からの参照ボリュームデータＤ１と入力ボリュームデータ取込部２からの入力ボリュームデータＤ２との間の撮像条件に起因する相違を補正する（ステップＳ１０３）。この例では、撮像条件に起因する相違として、参照ボリュームデータＤ１と入力ボリュームデータＤ２との間の倍率や大まかな回転ずれや位置ずれを補正する。補正された参照ボリュームデータＤ１はボリュームデータ記憶部６に保存される。また、補正された参照ボリュームデータＤ１および入力ボリュームデータＤ２は、剛体変形処理部４へ送られる。
【００３６】
なお、撮像条件補正部３において、倍率については装置が持つスケール（１ボクセルと実座標上の大きさ〔ｍ〕の対応）を利用して合わせるなどする。回転ずれについては、それぞれのデータに主成分を計算して合わせるなどする。位置ずれについては、それぞれのデータの重心を利用して合わせるなどする。
【００３７】
〔オブジェクト領域の抽出〕
剛体変形処理部４において、オブジェクト領域抽出部４Ａは、撮像条件補正部３で補正された参照ボリュームデータＤ１および入力ボリュームデータＤ２を取り込み、この参照ボリュームデータＤ１および入力ボリュームデータＤ２から所定の領域をオブジェクト領域としてそれぞれ抽出する（ステップＳ１０４）。
【００３８】
この例では、参照ボリュームデータＤ１および入力ボリュームデータＤ２に含まれる背景を除去し、対象物に関する情報が含まれる領域のみをオブジェクト領域として抽出する。以下、この抽出されたオブジェクト領域のボリュームデータを参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２とする。
【００３９】
なお、ここでの背景とは、ＣＴやＭＲＩなどから得られる医用ボリュームデータで言えば、対象の周りに存在する空気を指す。空気は、対象とは全く異なるボクセル値を持つので、ボクセル値に閾値処理を施すことで容易に対象と空気を分離することができる。
【００４０】
〔対応点の探索（３次元ブロックマッチング）〕
対応点探索部４Ｂは、オブジェクト領域抽出部４Ａからの参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２を取り込み、この参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２間の対応点の探索を行う（ステップＳ１０５）。図３および図４にこの場合の対応点の探索処理のフローチャートを示す。
【００４１】
対応点探索部４Ｂは、先ず、参照ボリュームデータＯＢ１にＮ個の基準点（ｘｉ，ｙｉ，ｚｉ：ｉ＝１，２，・・・・，Ｎ）を定める（ステップＳ２０１（図３））。また、参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２をリサイズして行くことで、各辺の長さ（ｘ，ｙ，ｚ軸方向の長さ）が元のデータの２−ｌとなる一連の階層ボリュームデータを生成する（ステップＳ２０２）。
【００４２】
なお、本実施の形態のリサイズでは、隣接する８ボクセルを平均化することで１ボクセルとし、ボリュームデータの各辺をそれぞれ２分の１の長さにする。しかし、必ずしもこのような方法をとる必要はなく、他にも、例えば平滑化フィルタをかけた後に単純にダウンサンプリング（１ボクセルおきにサンプリング）する方法などが考えられる。
【００４３】
すなわち、参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２を段階的に縮小し、元の参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２を最下位階層とし、最も縮小された参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２を最上位階層とする複数階層（ｌ＝０，１，２，・・・・ｌｍａｘ）のボリュームデータを生成する。以下、参照ボリュームデータＯＢ１の各階層のボリュームデータを第１群のボリュームデータ、入力ボリュームデータＯＢ２の各階層のボリュームデータを第２群のボリュームデータと呼ぶ。
【００４４】
そして、ｉ＝１とし（ステップＳ２０３）、第１群の最上位階層ｌｍａｘのボリュームデータに初期基準点（floor（ｘｉ／２lmax），floor（ｙｉ／２lmax），floor（ｚｉ／２lmax））をｐｌｍａｘとして設定し（図５参照）、第２群の最上位階層ｌｍａｘのボリュームデータの第１群の最上位階層ｌｍａｘに設定した初期基準点ｐｌｍａｘに対応する点に初期対応点をｑｌｍａｘとして設定する（ステップＳ２０４）。この際、最上位階層ｌｍａｘにおいては基準点と対応点とが同じ座標を持つと仮定し、初期基準点ｐｌｍａｘ＝初期対応点ｑｌｍａｘとする。
【００４５】
次に、ｌ＝ｌｍａｘ−１とし（ステップＳ２０５）、第１群の階層ｌ＝ｌｍａｘ−１のボリュームデータに初期基準点ｐｌｍａｘに対応する基準点（floor（ｘｉ／２lmax-1），floor（ｙｉ／２lmax-1），floor（ｚｉ／２lmax-1））をｐｌｍａｘ−１として定め、この基準点ｐｌｍａｘ−１を中心とする３辺の長さが予め定められた３次元のブロック（部分領域）を探索ウィンドウＷ１として設定する（ステップＳ２０６：図５参照）。
【００４６】
基準点（floor（ｘｉ／２lmax-1），floor（ｙｉ／２lmax-1），floor（ｚｉ／２lmax-1））における「floor」の意味について説明する。Ｂ=floor(Ａ)は、ＡをＡ以下の最も近い整数とする。ボリュームデータに設定した基準点（ｘｉ， ｙｉ， ｚｉ）は、階層探索を行う際に解像度を変更したボリュームデータに合わせて（ｘｉ／２lmax，ｙｉ／２lmax，ｚｉ／２lmax）とする必要があるが、この際に必ずしも座標が整数とはならない。このため、本実施の形態では、整数とするために、floor関数を用いている。
【００４７】
また、第２群の階層ｌ＝ｌｍａｘ−１のボリュームデータに初期対応点ｑｌｍａｘに対応する対応候補点（最上位階層ｌｍａｘに設定した初期対応点ｑｌｍａｘを２倍した座標）をｑｌｍａｘ−１’として定め、この対応候補点ｑｌｍａｘ−１’を中心とする３辺の長さが予め定められた３次元のブロック（部分領域）を探索ウィンドウＷ２として設定する（ステップＳ２０７（図４）：図５参照）。なお、探索ウインドウＷ１，Ｗ２の大きさは同じであり、以下に説明する各階層でも同じ大きさの探索ウィンドウＷ１，Ｗ２を用いるものとする。
【００４８】
そして、対応点探索部４Ｂは、ステップＳ２０６で設定した基準点ｐｌｍａｘ−１を中心とする探索ウィンドウＷ１とステップＳ２０７で設定した対応候補点ｑｌｍａｘ−１’を中心とする探索ウィンドウＷ２のボクセルデータを抽出し、例えば後述する位相限定相関法（３次元ＰＯＣ）を用いてウインドウ間のずれをボクセルレベルで求め、このずれから探索ウインドウＷ２におけるボクセルレベルの対応点ｑｌｍａｘ−１を決定する（ステップＳ２０８）。すなわち、第２群の階層ｌｍａｘ−１のボリュームデータにおける対応点ｑｌｍａｘ−１を決定する。
【００４９】
図６にステップＳ２０８で行われる３次元ＰＯＣの処理の流れを示す。対応点探索部４Ｂは、基準点ｐｌｍａｘ−１を中心とする探索ウィンドウＷ１のボクセルデータを抽出し、この抽出したボクセルデータに対して３次元フーリエ変換を行う（ステップＳ３０１）。また、対応候補点ｑｌｍａｘ−１’を中心とする探索ウィンドウＷ２のボクセルデータを抽出し、この抽出したボクセルデータに対して３次元フーリエ変換を行う（ステップＳ３０２）。
【００５０】
そして、この３次元フーリエ変換が施された探索ウィンドウＷ１の３次元スペクトルと探索ウィンドウＷ２の３次元スペクトルとを合成し（ステップＳ３０３）、この合成した３次元スペクトルの振幅を正規化して位相成分のみとし（ステップＳ３０４）、この位相成分のみとした３次元スペクトルに３次元逆フーリエ変換を施して（ステップＳ３０５）、２つの探索ウィンドウＷ１，Ｗ２の相関を表す位相限定相関関数を得る。
【００５１】
この位相限定相関関数から階層ｌ＝ｌｍａｘ−１のボリュームデータの空間内における探索ウィンドウＷ１とＷ２との間のずれをボクセルレベルで求め、このずれから探索ウィンドウＷ２におけるボクセルレベルの対応点ｑｌｍａｘ−１を決定する。なお、３次元ＰＯＣの詳細については後述する。
【００５２】
次に、対応点探索部４Ｂは、ｌ＝ｌ−１とし（ステップＳ２０９）、すなわちｌ＝ｌｍａｘ−２とし、ステップＳ２０６〜Ｓ２０８の処理を繰り返す。
【００５３】
この場合、対応点探索部４Ｂは、第１群の階層ｌｍａｘ−２のボリュームデータに基準点ｐｌｍａｘ−１に対応する基準点（floor（ｘｉ／２lmax-2），floor（ｙｉ／２lmax-2），floor（ｚｉ／２lmax-2））をｐｌｍａｘ−２として定め、この基準点ｐｌｍａｘ−２を中心とする探索ウィンドウ（３次元ブロック）Ｗ１を設定する（ステップＳ２０６：図５参照）。
【００５４】
また、第２群の階層ｌ＝ｌｍａｘ−２のボリュームデータに対応点ｑｌｍａｘ−１に対応する対応候補点（階層ｌｍａｘ−１で求めた対応点ｑｌｍａｘ−１を２倍した座標）をｑｌｍａｘ−２’として定め、この対応候補点ｑｌｍａｘ−２’を中心とする探索ウィンドウ（３次元ブロック）Ｗ２を設定する（ステップＳ２０７：図５参照）。
【００５５】
そして、基準点ｐｌｍａｘ−２を中心とする探索ウィンドウＷ１と対応候補点ｑｌｍａｘ−２’を中心とする探索ウィンドウＷ２のボクセルデータを抽出し、３次元ＰＯＣを用いてウインドウ間のずれをボクセルレベルで求め、このずれから探索ウインドウＷ２におけるボクセルレベルの対応点ｑｌｍａｘ−２を決定する（ステップＳ２０８）。すなわち、第２群の階層ｌｍａｘ−２のボリュームデータにおける対応点ｑｌｍａｘ−２を決定する。
【００５６】
以下、同様にして、対応点探索部４Ｂは、第１群の階層および第２群の階層を１段階ずつ下げながら、ステップＳ２０６〜２０８の処理を繰り返して行く。そして、最下位階層におけるボクセルレベルの対応付けが終了し、ｌ＝−１となると（ステップＳ２１０のＹＥＳ）、対応点探索部４Ｂは、参照ボリュームデータＯＢ１と入力ボリュームデータＯＢ２において、基準点ｐ０を中心とする探索ウィンドウＷ１のボクセルデータとｌ＝０の階層においてステップＳ２０６〜Ｓ２０８を実行することによって得られたボクセルレベルの対応候補点ｑ０を中心とする探索ウインドウＷ２のボクセルデータを抽出する（ステップＳ２１１）。
【００５７】
そして、３次元ＰＯＣを用いてウインドウ間のずれをサブボクセルレベル（ボクセルレベルよりもさらに細かいレベル）で求め、このずれから探索ウインドウＷ２における対応点ｑ０をサブボクセルレベルへ更新し、この更新したサブボクセルレベルの対応点ｑ０を基準点ｐ０に対応する対応点として記憶する（ステップＳ２１２）。すなわち、参照ボリュームデータＯＢ１に定められるｉ＝１番目の基準点ｐ０に対応する入力ボリュームデータＯＢ２における正規の対応点ｑ０として記憶する。
【００５８】
そして、対応点探索部４Ｂは、ｉ＝ｉ＋１とし（ステップＳ２１４）、すなわちｉ＝２とし、２番目の基準点について、ステップＳ２０４〜Ｓ２１２の処理を繰り返す。これにより、上述と同様にして、参照ボリュームデータＯＢ１に定められる基準点ｐ０に対応する入力ボリュームデータＯＢ２における対応点ｑ０が決定され、この決定された対応点ｑ０がｉ＝２番目の基準点に対応する正規の対応点として記憶される。この処理は、ステップＳ２１３においてｉ＝Ｎとなるまで繰り返される。
【００５９】
これにより、参照ボリュームデータＯＢ１におけるＮ個の基準点について、各基準点に対応する入力ボリュームデータＯＢ２における正規の対応点が次々に決定されて行く。
【００６０】
本実施の形態では、上述したように、ボリュームデータの解像度を変えたものをいくつか生成して階層状にし、低い解像度のデータから高い解像度のデータへ順番に対応点探索を行っていくことで、効率の良い対応付けが可能となり、また、データ間で大きな動きが生じている場合にも対応することができるようになる。
【００６１】
〔剛体変形パラメータの推定〕
剛体変形パラメータ推定部４Ｃは、対応点探索部４Ｂによって探索された対応点と基準点との対応関係から参照ボリュームデータＯＢ１と入力ボリュームデータＯＢ２との間の回転ずれおよび位置ずれを表すパラメータを剛体変形パラメータとして推定する（ステップＳ１０６（図２））。
【００６２】
例えば、回転ずれとしてＸ，Ｙ，Ｚ軸の各軸の回転ずれ量Δθｘ，Δθｙ，Δθｚを求め、位置ずれとして各軸の平行ずれ量Δｘ，Δｙ，Δｚを求める。この例では、回転ずれとして回転行列を求め、位置ずれとして並進ベクトルを求める。剛体変形パラメータの推定は、たとえば後述するように最小二乗法を用いて行うことができる。
【００６３】
なお、この剛体変形パラメータ推定部４Ｃにおいて剛体変形パラメータを推定する場合、例えば非特許文献１４に示されているようなロバスト推定法を用いるようにすれば、対応関係に誤りが含まれていても高精度にパラメータの推定を行うことが可能である。
【００６４】
〔剛体変形の実行〕
剛体変形実行部４Ｄは、剛体変形パラメータ推定部４Ｃによって推定された剛体変形パラメータに基づいて、撮影条件補正部３から送られてくる入力ボリュームデータＤ２の回転ずれおよび位置ずれを補正（剛体変形）する（ステップＳ１０７）。すなわち、入力ボリュームデータＤ２を剛体変形させることによって、ボリュームデータ記憶部６に記憶されている参照ボリュームデータＤ１との間の回転ずれおよび位置ずれを補正する。
【００６５】
そして、この剛体変形させた入力ボリュームデータＤ２をＤ２’として、反復処理命令部４Ｅへ送る。なお、剛体変形実行部４Ｄは、剛体変形後の入力ボリュームデータＤ２’を、次の回の剛体変形前の入力ボリュームデータＤ２として記憶する。
【００６６】
〔反復処理〕
反復処理命令部４Ｅは、剛体変形実行部４Ｄからの剛体変形された入力ボリュームデータＤ２’をオブジェクト領域抽出部４Ａへ送り、オブジェクト領域抽出部４Ａでのオブジェクト領域の抽出（ステップＳ１０４）、対応点探索部４Ｂでの対応点の探索（ステップＳ１０５）、剛体変形パラメータ推定部４Ｃでの剛体変形パラメータの推定（ステップＳ１０６）、剛体変形実行部４Ｄでの剛体変形の実行（ステップＳ１０７）を再度行わせる。
【００６７】
反復処理命令部４Ｅは、このステップＳ１０４〜Ｓ１０７の処理を繰り返させ、その回数が所定回数に達した時点で（ステップＳ１０８のＹＥＳ）、その時の補正された入力ボリュームデータＤ２’を剛体変形処理後の入力ボリュームデータとしてボリュームデータ記憶部６に記憶させる（ステップＳ１０９）。
【００６８】
このような反復処理により、剛体変形パラメータの推定とこの推定された剛体変形パラメータを用いての入力ボリュームデータの剛体変形が繰り返され、入力ボリュームデータと参照ボリュームデータとの間で大きな移動や回転が生じている場合にも対応することができるようになる。
【００６９】
〔非剛体変形処理〕
次に、図７に示すフローチャートを用いて、ボリュームデータ間位置合わせ装置１００で行われる非剛体変形処理について説明する。
【００７０】
〔オブジェクト領域の抽出〕
非剛体変形処理部５において、オブジェクト領域抽出部５Ａは、ボリュームデータ記憶部６に記憶されている参照ボリュームデータＤ１を読み出し（ステップＳ４０１）、またボリュームデータ記憶部６に記憶されている剛体変形処理後の入力ボリュームデータＤ２’を読み出し（ステップＳ４０２）、この参照ボリュームデータＤ１および入力ボリュームデータＤ２’からオブジェクト領域を抽出し、参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２とする（ステップＳ４０３）。
【００７１】
なお、この場合、参照ボリュームデータＤ１に対しては剛体変形処理部４での剛体変形は実施されていないが、入力ボリュームデータＤ２に剛体変形を実施することによって参照ボリュームデータＤ１との間の回転ずれおよび位置ずれが補正されており、両者に対して補正が行われているものと考えれば、参照ボリュームデータＤ１も剛体変形処理後のボリュームデータであると言える。ここでは、参照ボリュームデータＤ１も剛体変形処理後のボリュームデータと考える。
【００７２】
〔対応点の探索〕
対応点探索部５Ｂは、オブジェクト領域抽出部５Ａからの参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２を取り込み、この参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２間の対応点の探索を行う（ステップＳ４０４）。この対応点探索部５Ｂでの対応点の探索処理は、剛体変形処理部４における対応点探索部４Ｂでの対応点の探索処理（図３）と同じであるので、その説明は省略する。
【００７３】
非剛体変形実行部５Ｃは、対応点探索部５Ｂによって探索された対応点の情報を用いて、ＦＦＤ（Free-Form Deformation）により入力ボリュームデータＤ２’を補正（非剛体変形）し、剛体変形処理後の入力ボリュームデータＤ２’と参照ボリュームデータＤ１の非剛体部分を重ね合わせる（ステップＳ４０５）。そして、この補正された入力ボリュームデータＤ２’を非剛体変形処理後の入力ボリュームデータＤ２”として、ボリュームデータ記憶部６に記憶させる（ステップＳ４０６）。ＦＦＤによる非剛体変形処理の詳細については後述する。
【００７４】
なお、ＦＦＤ（Free-Form Deformation）については、例えば非特許文献７にも開示されているように、ボリュームデータの非剛体位置合わせにおいては一般的である。本発明によりあらかじめ得られたボクセルの対応関係を利用したＦＦＤによる非剛体変形は、従来技術であるボクセル類似度ベースの手法を用いたＦＦＤに比べて約１００倍の計算速度向上が見られた。従来技術ではボリュームデータの対応関係が未知の状態で計算を始めるため特に計算量の多い非剛体変形では非線形最適化に膨大な時間がかかるが、あらかじめ本発明により高精度にボリュームデータの対応関係が得られていれば、大幅な計算速度向上が可能である。
【００７５】
なお、上述した実施の形態では、剛体変形処理部４の対応点探索部４Ｂにおいて対応点を探索する場合、参照ボリュームデータＯＢ２の全体に基準点を定めるようにするが、例えばＣＴやＭＲＩなどの医用ボリュームデータを扱う場合には、骨などの形状変化が少ない部分に集中的に定めることで、位置合わせ精度を向上させることが可能である。
【００７６】
例えば、参照ボリュームデータＯＢ２について、骨のある領域を所定の特徴のある領域として検索させ、この検索された特徴のある領域に基準点を定めるようにする。すなわち、検索された特徴のある領域を探索ウィンドウの抽出領域とする。非剛体変形処理部５の対応点探索部５Ｂにおいて対応点を探索する場合も同様である。
【００７７】
また、この実施の形態では、オブジェクト領域抽出部４Ａや５Ａにおいて背景を除去して対象物のある領域のみを抽出するようにしたが、上述した骨のある領域など特徴のある領域をオブジェクト領域として抽出するようにしてもよい。また、撮像条件補正部３やオブジェクト領域抽出部４Ａ，５Ａは必ずしも設けなくてもよい。
【００７８】
また、上述した実施の形態では、対応点探索部４Ｂや５Ｂにおいて、参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２の階層は、複数階層であれば何階層であってもよく、２層としても構わない。例えば、例えば図５に示した例では、最下位階層ｌ０を０層とした場合、ｌｍａｘは３層であるが、ｌｍａｘを２層としてもよい。また、ｌｍａｘを４層、５層と、その層を増やして行ってもよい。
【００７９】
また、上述した実施の形態では、対応点の探索を複数階層で行うようにしたが、必ずしも複数階層を利用する必要はなく、対象にそれほど大きなずれがない場合は１層でも対応付けは可能である。
【００８０】
また、上述した実施の形態では、剛体変形処理部４におけるオブジェクト領域抽出部４Ａと非剛体変形処理部５におけるオブジェクト領域抽出部４Ａとを別としたが、また、剛体変形処理部４における対応点探索部４Ｂと非剛体変形処理部５における対応点探索部４Ｂとを別としたが、共通のブロックとすることも可能である。
【００８１】
また、上述した実施の形態では、対応点探索部４Ｂや５Ｂにおいて、３次元ＰＯＣを用いてウインドウ間のずれを求めるようにしたが、必ずしも３次元ＰＯＣ（位相限定相関法）を用いなくてもよく、他の相関法を用いてもよい。また、位相に限定するのではなく、ｌｏｇ処理やルート処理などによって振幅を抑制する振幅抑制相関法を用いてもよい。
【００８２】
また、上述した実施の形態では、剛体変形処理部４における剛体変形実行部４Ｄにおいて、入力ボリュームデータＤ２に対して剛体変形を実施するようにしたが、参照ボリュームデータＤ１に対して剛体変形を実施するようにしてもよい。また、非剛体変形処理部５における非剛体変形実行部５Ｃにおいて、剛体変形処理後の入力ボリュームデータＤ２’に対して非剛体変形を実施するようにしたが、参照ボリュームデータＤ１（剛体変形処理後の参照ボリュームデータＤ１）に対して非剛体変形を実施するようにしてもよい。
【００８３】
また、上述した実施の形態では、対応点探索部４Ｂや５Ｂにおいて対応点を探索する場合、参照ボリュームデータＯＢ２に基準点を複数定め、その基準点を中心する探索ウィンドウを抽出するようにしたが、探索ウィンドウを１つずつ切り出し、その探索ウィンドウの中心を基準点として、対応点を探索して行くようにしてもよい。
【００８４】
また、上述した実施の形態では、剛体変形処理部４の反復処理命令部４Ｅにおいて、所定回数に達するまでステップＳ１０４〜Ｓ１０７の処理を繰り返させるようにしたが、補正された入力ボリュームデータＤ２’と参照入力ボリュームデータＤ１との回転ずれおよび位置ずれが所定の閾値よりも小さくなった場合に、反復処理を停止し、その時の補正された入力ボリュームデータＤ２’を剛体変形処理後の入力ボリュームデータとしてボリュームデータ記憶部６に記憶させるようにしてもよい。またこの他にも、推定される剛体変形パラメータの変化率が所定の閾値よりも小さくなった場合に反復処理を停止してもよい。
【００８５】
また、上述した実施の形態では、参照ボリュームデータＤ１をＣＴによって撮像されたものとし、入力ボリュームデータＤ２をＭＲＩによって撮像されたものとしたが、参照ボリュームデータＤ１と入力ボリュームデータＤ２をともにＣＴによって撮像されたものとしたり、ＭＲＩによって撮像されたものとしてもよく、さらにこれらとも異なる別の撮像装置で撮像されたデータとしてもかまわない。
【００８６】
次に、３次元ＰＯＣを用いたボリュームデータ間の位置合わせ手法について、さらに詳細に説明する。
【００８７】
〔３次元位相限定相関法〕
ここでは３次元ＰＯＣの基本的な定義と、３次元ＰＯＣを用いたボリュームマッチングの高精度化について述べる。
【００８８】
〔３次元位相限定相関関数〕
Ｎ1×Ｎ2×Ｎ3ボクセルの２つのボリュームデータｆ（ｎ1，ｎ2，ｎ3）およびｇ（ｎ1，ｎ2，ｎ3）が与えられたとする。ここで、ボリュームデータの離散空間インデックス（整数）を、便宜上、ｎ1＝−Ｍ1、・・・、Ｍ1、ｎ2＝−Ｍ2、・・・、Ｍ2、ｎ3＝−Ｍ3、・・・、Ｍ3とする。ただし、Ｍ1、Ｍ2、Ｍ3は正の整数であり、Ｎ1＝２Ｍ1＋１、Ｎ2＝２Ｍ2＋１、Ｎ3＝２Ｍ3＋１である。
【００８９】
なお、ここでは説明を簡単にするために離散空間のインデックスを正負対称にとり、かつボリュームデータの大きさＮ1，Ｎ2，Ｎ3を奇数にしているが、これは必須ではない。すなわち、通常よく用いられるように非負のインデックスを用い、Ｎ１，Ｎ2，Ｎ3を任意の正の整数に設定するように一般化することが可能である。ボリュームｆ（ｎ1，ｎ2，ｎ3）およびｇ（ｎ1，ｎ2，ｎ3）の３次元離散フーリエ変換（Discrete Fourie Transform：ＤＦＴ）Ｆ（ｋ1，ｋ2，ｋ3）およびＧ（ｋ1，ｋ2，ｋ3）をそれぞれ次式で定義する。
【００９０】
【数１】
【００９１】
３次元ＰＯＣ関数ｒ（ｎ1，ｎ2，ｎ3）を、正規化相互パワースペクトルの３次元逆離散フーリエ変換（Inverse Discrete Fourier Transform：ＩＤＦＴ）として定義する。
【００９２】
【数２】
【００９３】
すると、（４）式に示す３次元ＰＯＣ関数ｒ（ｎ1，ｎ2，ｎ3）は、２つのボリュームデータｆ（ｎ1，ｎ2，ｎ3）とｇ（ｎ1，ｎ2，ｎ3）の３次元の相互相関を表すことになり、２つのボリュームが類似している場合、３次元ＰＯＣ関数は、デルタ関数に近いきわめて鋭いピークを有する。この相関ピークの高さは２つのボリュームの位相差スペクトルθF（ｋ1，ｋ2，ｋ3）−θR（ｋ1，ｋ2，ｋ3）の線形性を表しており、位相差スペクトルが完全に線形であれば、相関ピークの高さは１となる。この相関ピークの高さはボリュームの類似度の尺度として有用である。一方、相関ピークの座標は２つのボリュームの相対的な位置ずれに対応している。
【００９４】
以下では、ｇ（ｎ1，ｎ2，ｎ3）がｆ（ｎ1，ｎ2，ｎ3）を（δ1，δ2，δ3）だけ微小に平行移動させたボリュームである場合を考える。ここで、δ1，δ2，δ3は、それぞれｎ1，ｎ2，ｎ3方向のサブボクセルレベルの移動量を表している。このとき、ｆ（ｎ1，ｎ2，ｎ3）とｇ（ｎ1，ｎ2，ｎ3）の３次元ＰＯＣ関数は、次式で与えられる。
【００９５】
【数３】
【００９６】
ここで、α＝１である。上式は、ボリュームが（δ1，δ2，δ3）だけ微小に平行移動した場合の３次元ＰＯＣ関数の一般形を表している。αは、相関ピークの高さを表現するために導入されたパラメータである。ボリュームに対して無相関なノイズが加わるとαの値が減少するため、実際にはα≦１となる。この相関ピークのモデルに基づく関数フィッティングにより、パラメータα，δ1，δ2，δ3を推定することで、ボリュームの類似度（位相差スペクトルの線形性）とサブボクセル精度の移動量を求めることができる。
【００９７】
〔３次元ＰＯＣに基づくボリュームマッチングの高精度化〕
以下（Ａ）〜（Ｃ）では、３次元ＰＯＣに基づくボリュームマッチングの高精度化について述べる。
（Ａ）窓関数の適用
ＤＦＴでは、信号が周期的に循環することを仮定するため、画像端での信号の不連続性が問題となる。この不連続性の影響を軽減するため、ボリュームに対して窓関数を適用することが重要である。本手法では下記のハニング窓を用いる。
【００９８】
【数４】
【００９９】
（Ｂ）スペクトルの重み付け
一般に、自然画像のエネルギーは低周波領域に集中し、高周波成分のエネルギーは相対的に小さいことが知られている。このため、エイリアシング、ぼけ、雑音、歪みなどの外乱が加わると、高周波成分のＳ／Ｎが大幅に劣化する。そこで、信頼性の低い高周波成分の影響を抑制するために、正規化相互パワースペクトルＲ（ｋ1，ｋ2，ｋ3）の計算の際に、低域通過型のスペクトル重み付け関数Ｈ（ｋ1，ｋ2，ｋ3）を適用することにより、大幅な精度向上が可能である。本手法ではＨ（ｋ1，ｋ2，ｋ3）として、次式で与えられるガウス関数を用いる。
【０１００】
【数５】
【０１０１】
この場合、３次元ＰＯＣ関数は、Ｈ（ｋ1，ｋ2，ｋ3）Ｒ（ｋ1，ｋ2，ｋ3）の３次元ＩＤＦＴであり、式（５）に対応する相関ピークのモデルは、次式のように変化する。
【０１０２】
【数６】
【０１０３】
なお、重み付け関数については、これら以外にも用途に応じて、任意の関数を適用することが可能である［非特許文献８参照］。
【０１０４】
（Ｃ）ピークモデルのフィッティング
一般に、移動量（δ1，δ2，δ3）は実数値であり、３次元ＰＯＣ関数のピーク座標がサンプリング格子点の間に存在するため、正確に移動量を推定することが困難である。そこで、相関ピークのモデルが式（５）、（９）で与えられることを考慮し、実際に計算されたＰＯＣ関数の数値データに対して本モデルをフィッティングすることで、ボクセル間に存在するピークの位置を推定する。
【０１０５】
上記（Ｂ）のスペクトル重み付け関数Ｈ（ｋ1，ｋ2，ｋ3）を適用した場合は、式（９）の例のように、Ｈ（ｋ1，ｋ2，ｋ3）に応じた相関ピークモデルをフィッティングする必要がある。このとき、α，δ1，δ2，δ3がフィッティングパラメータとなる。
【０１０６】
なお、式（５）の形の相関ピークモデルの場合は、ＰＯＣ関数の数値データから直接的に相関ピークの座標と値を求めるピーク評価式（Peak Evaluation Formula：ＰＥＦ）が導出されており［非特許文献１０参照］、リアルタイム処理に寄与している。それ以外の一般的な相関ピークモデルの場合は、Levenberg-Marquardt法などの非線形最小二乗法を用いたフィッテイングを行う。
【０１０７】
このように、３次元ＰＯＣ関数の相関ピーク形状は、入力ボリュームによらず、Ｈ（ｋ1，ｋ2，ｋ3）のみによって決定されるという特徴があるため、関数フィッティングによって高精度な移動量推定が可能である。
【０１０８】
〔３次元位相限定相関法に基づくボリュームレジストレーション〕
本手法で提案するボリュームレジストレーション手法の詳細について述べる。提案アルゴリズムは２つのステップで構成され、まず、（i）３次元ＰＯＣに基づいてサブボクセルレベルの対応点探索を行い、次に（ii）ボクセルの対応関係からボリュームデータの変形パラメータを推定する。
【０１０９】
まず、ボリュームの変形モデルとして剛体変形を用い、剛体変形を記述するパラメータを推定する。本手法は非剛体変形への拡張も容易であり、本節ではＦＦＤを用いた非剛体レジストレーションについても述べる。以下では本手法における各処理の詳細を述べる。
【０１１０】
〔サブボクセル対応点探索〕
本手法では、まず２つのボリュームデータから小領域（部分領域（３次元ブロック））を切り出し、３次元ＰＯＣを用いてブロック間の移動量をサブボクセル精度で求めることで対応付けを行う。このとき、粗密戦略に基づく階層探索を行うことで効率的な対応付けを実現している。以下では、ボリュームＩ上に設定した基準点ｐ＝（ｐ1，ｐ2，ｐ3）に対応するＪ上の対応点ｑ＝（ｑ1，ｑ2，ｑ3）を求めるものとし、対応付け手法の詳細について述べる。
【０１１１】
【数７】
【０１１２】
【数８】
【０１１３】
以上に示した処理Ｓｔｅｐ１〜Ｓｔｅｐ６によって、基準点ｐの対応点ｑをサブボクセル精度で求めることができる。
【０１１４】
〔剛体変形パラメータの推定〕
前述のブロックマッチングで得られた対応関係から剛体変形パラメータを推定する。剛体変形を記述するパラメータは、前述の対応点が３点以上あれば推定が可能である。本手法では、行列の特異値分解に基づいて回転行列Ｒおよび並進ベクトルｔの最小二乗解を求めるアルゴリズム［非特許文献１１参照］を用いる。以下にその手順を示す。
【０１１５】
【数９】
【０１１６】
〔反復推定〕
ボリューム間に大きな回転などが含まれる場合、一度の対応点探索及び剛体変形パラメータの推定だけでは高精度な位置合わせが難しいことが多い。そこで、まず前述の方法で推定された剛体変形パラメータを利用し、ボリュームデータ間の変形が小さくなるように剛体変形を施す。
【０１１７】
そして、再び対応付けを行い、剛体変形パラメータを推定する。このように、対応付け、剛体変形パラメータの推定、ボリュームの変形を繰り返すことで、ボリューム間で大きな変形が生じている場合でも正確に位置合わせを行うことが可能となる。実際には、例えばボリュームが一方向に２０°程度回転している場合、約１０回の繰り返し推定で誤差が収束することを確認している。
【０１１８】
〔非剛体レジストレーション〕
医用ボリュームデータにおいて、心臓や肺といった臓器に加え関節や筋肉などの組織は、それ自体が変形するなど、剛体変形では表現が不可能な変形を伴う場合があるため、剛体レジストレーションだけでは不十分であるといえる。したがって、剛体レジストレーションに加え、局所的な変形を補正するために非剛体レジストレーションを行う必要がある。
【０１１９】
非剛体レジストレーションでは、一度剛体レジストレーションを行った後、再度対応点探索を行い、この対応関係を用いて非剛体な変形パラメータを推定する。以下では、各ボクセルの移動をFree-Form Deformation（ＦＦＤ）モデルにより推定する手法について述べる。
【０１２０】
本手法では、ローカルな変形モデルとしてＢ-スプライン補間に基づくＦＦＤ［非特許文献１２参照］を用いる。ＦＦＤは３次元変形物体のモデリングに有用であり、これまで主に人体の運動の追跡や解析に利用されてきた。Ｂ-スプライン補間に基づくＦＦＤでは、ボリュームデータの３次元領域内に格子状に制御点が配置され、この制御点の移動がボリュームの変形を決定付ける。ｘ，ｙ，ｚ方向にそれぞれ間隔δｘ，δｙ，δｚで格子状の制御点φｉ，ｊ，ｋが並ぶとき、任意の座標ｘ＝（ｘ，ｙ，ｚ）における移動ベクトルＴ（ｘ）は次式で表される。
【０１２１】
【数１０】
【０１２２】
ここで入力ボリュームをＩ，Ｊとすると、ｐｉ，ｊ，ｋはＩ上の基準点であり、ｑｉ，ｊ，ｋは３次元ＰＯＣを用いた対応付けによって得られたＪ上の対応点である。しかし、後述する高精度化手法を用いて対応点を得た場合、基準点は格子状とならない。そのような場合にはマルチレベルＢ-スプライン［非特許文献１２参照］を利用してＦＦＤを計算する必要がある。マルチレベルＢ-スプラインは、不規則に並んだデータ点からＣ２−連続な補間を高速に計算するアルゴリズムである。本手法の詳細については非特許文献１２、１３を参考にされたい。
【０１２３】
〔レジストレーションの高精度化〕
ボリューム間の対応関係から剛体変形パラメータを推定する場合、対応関係に誤り（誤対応）が含まれていると位置合わせ精度が低下してしまう。そこで、ここでは、３次元ＰＯＣに基づく対応付けにおける誤対応を削減し、剛体変形パラメータの推定精度を向上させるための手法について述べる。
【０１２４】
〔基準点の選別〕
例えば医用ボリュームデータなどを扱う場合、断層画像に造影する部位は撮影に利用した装置のモダリティ（ＣＴやＭＲＩなど）によって異なる。レジストレーションにより両者に不足する情報を補完し合うことは、診断能の向上やより詳細な人体構造の把握につながるため、モダリティの異なるボリュームをレジストレーションする意義は大きい。しかし、モダリティが異なるボリュームデータに対して３次元ＰＯＣに基づく対応点探索を適用する場合、両者で位相情報が全く異なるような箇所に基準点を設定すると正しく対応点を求めることができず、誤対応となる。
【０１２５】
これを避けるため、２つのデータ間で位相情報が共通するような箇所に基準点を配置する。例えばＸ線ＣＴとＭＲＩで位相情報が共通となるような部位としては、骨や皮膚などが挙げられるが、剛体レジストレーションにおいてはその運動が剛体変形で記述しやすい硬組織周辺を利用するのが望ましい。そこで本手法では、ＣＴ値に基づいてＸ線ＣＴデータ上で骨に相当する箇所に基準点を配置し、対応点探索を行う。
【０１２６】
〔ロバスト推定の利用〕
前述の基準点選別に加え、ＰＯＣ関数のピーク値に対して閾値を設けることでさらに誤対応を減らすことが可能である。３次元ブロック間で計算されるＰＯＣ関数は、両ブロックが類似しているほど高いピーク値を示す（１に近付く）ため、ピーク値が低い対応点は信頼性が低いと考えることができる。そこで、本手法では、このピーク値に対して閾値を設定し、閾値を下回る対応点は誤対応として除外する。なお、経験上、閾値は０．１程度が妥当である。
【０１２７】
しかし、実際には、ピーク値に対する閾値処理だけで全ての誤対応を除去することは困難である。そこで、本手法では、対応関係に誤りが含まれていても高精度に剛体変形パラメータを求めることができるロバスト推定手法として、ＲＡＮＳＡＣ（RANdom SAmple Consensus）［非特許文献１４参照］を用いる。以下に、ＲＡＮＳＡＣを用いた剛体変形パラメータの推定手法を示す。
【０１２８】
Ｓｔｅｐ１：基準点と対応点のペアからランダムに３組選び、その３組から剛体変形パラメータの推定を行う。
Ｓｔｅｐ２：Ｓｔｅｐ１で求めたパラメータに基づいて、全ての基準点に剛体変形を施す。
Ｓｔｅｐ３：剛体変形された基準点とその対応点の距離が閾値以内ならば、これをｉｎｌｉｅｒ（外れ点でない点）に含める。
Ｓｔｅｐ４：ｉｎｌｉｅｒの個数を数える。
Ｓｔｅｐ５：Ｓｔｅｐ１〜４をＮｒ回繰り返す。
Ｓｔｅｐ６：Ｎｒ個のｉｎｌｉｅｒ集合のうち最も要素数が多いものを用いて、変形パラメータを再推定する。
【０１２９】
〔実験・考察〕
医用ボリュームデータを例として、本手法（ボリュームレジストレーション手法）の性能評価を行う。まず、剛体レジストレーションについて従来手法と性能比較を行い、本手法では従来手法に比べて高精度かつ高速なレジストレーションが可能であることを示す。続いて、ＦＦＤを用いた非剛体レジストレーションを適用し、本手法により非剛体変形に対しても高精度な位置合わせが可能であることを示す。最後に、Ｘ線ＣＴとＭＲＩから得られたボリュームに対してレジストレーションを行い、本手法がモダリティの違いにロバストであることを示す。
【０１３０】
〔剛体レジストレーション〕
ここでは、歯科用ＣＴで撮影した顎部ボリュームデータを用いる。一人の被験者から異なるタイミングで取得した４つのデータを利用し、これらをそれぞれボリューム１〜４とする。そして、全ての組み合わせ（４Ｃ２＝６組）に対して位置合わせを行った。再構成したボリュームは５１２×５１２×５１２のボクセル数で構成され、０．２×０．２×０．２ｍｍの分解能をもつ。４つのボリュームのうち１つは意図的に首を傾けて撮影されているため、比較的大きな位置ずれや歪みが加わっている。図１０（ａ）、（ｂ）に実験にて使用したボリュームの一部を示す。
【０１３１】
本実験では、従来手法としてStudholmeらの手法［非特許文献５］を利用した。これは、２つのデータ間の正規化相互情報量（Normalized Mutual Information：ＮＭＩ）が最大となるように剛体変形パラメータを最適化する手法であり、医用データのレジストレーションにおいて一般に広く用いられている。非剛体レジストレーションでは、Rueckertらの手法[非特許文献７]を従来手法として用いた。Studholmeらの手法と同じく正規化相互情報量の最大化に基づく手法であるが、変形モデルとしてFFDを利用しているため非剛体なひずみを補正することが可能である。
【０１３２】
実験にあたり、ＰＯＣに基づく対応点探索に用いたパラメータを図８に示す。Studholmeらの手法における非線形最適化にはNelder-Mead法［非特許文献１５参照］を用い、最適化計算における初期値は適当な値を手動で設定した。
【０１３３】
レジストレーションの精度を評価するために、次の相関係数（Correlation Coefficient：ＣＣ）を指標として用いる。
【０１３４】
【数１１】
【０１３５】
本手法と従来手法を用いてレジストレーションを行い、それらの性能を相関係数により比較した結果を図９に示す。ボリューム１，２，４は自然な姿勢で撮影したものであり、ボリューム３は意図的に首を傾けて撮影したものである。相関係数による評価では１に近いほどレジストレーション精度が高いといえる。
【０１３６】
図９を見ると、本手法は全てのボリュームの組み合わせに対して従来手法と同等またはそれ以上の精度を示しており、高いレジストレーション精度を持っていることがわかる。また、本実験において、ボリュームデータ一組の剛体レジストレーションに要する時間は、Studholmeらの手法が約４分であるのに対し、本手法では約１分であり、本手法によって７５％の計算時間削減が達成されている。非剛体レジストレーションの場合は、Rueckertらの手法が約６６５分、本手法が約６分であり、本手法によって約９９％の計算時間削減が達成されている。
【０１３７】
なお、両手法の実装にはMＡＴＬＡＢ7.4.0を利用し、実験に使用したパーソナルコンピュータ（ＰＣ）はCPU：Xeon(3.0GHz×2)、Memory：16GB、OS:CentOS 4.7である。このように、３次元ＰＯＣに基づくボリュームレジストレーション手法は、精度、計算時間の両面において従来手法よりも優れていることがわかる。
【０１３８】
ボリューム１と３を用いた提案手法による非剛体レジストレーション結果を図１０（ｃ），（ｄ）に示す。ボリューム３に非剛体変形を加えて補正したデータのコロナル像にボリューム１の輪郭線を重ねて表示した。ボリューム１と３の間には非剛体な変形が加わっているため、剛体レジストレーションのみでは図１０（ｃ）のように位置ずれが目立つ。このような場合には非剛体レジストレーションが必要であり、最終的には本手法により図１０（ｄ）に示す通り非剛体なひずみが補正されていることがわかる。
【０１３９】
〔マルチモーダルレジストレーション〕
Ｘ線ＣＴとＭＲＩのそれぞれで同一の被験者の頭部データを取得し、本手法を用いて剛体レジストレーションを行った。ボリュームを再構成する際、ボクセル数を５１２×５１２×５１２，ボクセル分解能を０．４６９×０．４６９×０．４６９ｍｍとしている。
【０１４０】
図１１（ａ），（ｂ）にＸ線ＣＴデータおよびＭＲＩデータのアキシャル像を示す。図１１（ｃ），（ｄ）は、本手法によるレジストレーション結果である。図１１（ｃ）は、レジストレーション後のＭＲＩ画像上にＸ線ＣＴデータの輪郭線を重ねたものであり、図１１（ｄ）ではＭＲＩ画像とＸ線ＣＴ画像を格子状に交互に合成している。これらの結果より、モダリティが異なる場合でも本手法によって高精度なレジストレーションが可能であることがわかる。
【０１４１】
〔むすび〕
本手法では、位相限定相関法を３次元に拡張し、これを用いたボリュームデータの対応付け手法を提案した。本手法では、ＰＯＣを利用した３次元ブロックマッチングと階層探索を組み合わせることでサブボクセル精度でボリュームデータ間の対応付けを行う。ボクセル同士の対応関係はボリュームデータのレジストレーションなどの際にきわめて有用であり、従来技術に比べて精度を低下させることなく大幅に高速なレジストレーションを実現することができる。
【０１４２】
本手法は、モダリティの違いに対してもロバストに対応付けを行うことができ、また、非剛体レジストレーションにも容易に拡張することができる。タイミングを変えて取得した複数のＸ線ＣＴおよびＭＲＩデータを用いて、本手法が従来手法に比べて精度、計算時間の両面で高い性能を示すことを確認した。
【産業上の利用可能性】
【０１４３】
本発明のボリュームデータ間の対応付け方法は、同一対象物のほぼ同一箇所を撮像した３次元の立体画像のデータ（ボリュームデータ）間の対応付けを行う方法として、医療などの分野に限らず、各種の分野において適用することが可能である。
【符号の説明】
【０１４４】
１…参照ボリュームデータ取込部、２…入力ボリュームデータ取込部、３…撮像条件補正部、４…剛体変形処理部、４Ａ…オブジェクト領域抽出部、４Ｂ…対応点探索部、４Ｃ…剛体変形パラメータ推定部、４Ｄ…剛体変形実行部、４Ｅ…反復処理命令部、５…非剛体変形処理部、５Ａ…オブジェクト領域抽出部、５Ｂ…対応点探索部、５Ｃ…非剛体変形実行部、６…ボリュームデータ記憶部、Ｗ１，Ｗ２…探索ウィンドウ（３次元ブロック）、１００…ボリュームデータ間位置合わせ装置。
【技術分野】
【０００１】
この発明は、同一対象物を撮像した２つの３次元の立体画像のデータ（ボリュームデータ）間の対応付けを行うボリュームデータ間の対応付け方法に関するものである。
【背景技術】
【０００２】
ボリュームデータとは、３次元の小立方体（ボクセル）を単位として、各ボクセルにその部分の透過率などのデータを持たせ、物体の内部構造まで表現したデータのことであり、医療や流体解析、コンピュータグラフィックスなど広範な分野において頻繁に用いられるデータ形式である。
【０００３】
例えば医療分野においては、Ｘ線ＣＴ (Computed Tomography) やＭＲＩ (Magnetic Resonance Imaging) から得られるデータが挙げられる。これらの医用画像撮影装置は、人体内部のスライス像を連続的に取得することが可能であり、取得した２次元画像群を積層することで３次元ボリュームデータとして人体内部を可視化することができる。その有用性はすでに広く知られており、日常臨床でも一般的に用いられている。
【０００４】
ボリュームデータを扱う上では、複数のボリュームデータの比較や合成が要求される場面が頻繁に現れる。特に医療分野においては、異なる時間に撮影した病変部位を観察することで、病状の進行や解剖学的な変化について手掛かりを得ることができる。また、異なる装置を用いて撮影したボリュームデータをフュージョン（重ね合わせ）することで、それぞれの装置が苦手とする情報を補完した観察が可能となる。他にも、例えば多列検出器型ＣＴを用いて心臓などのボリュームデータを高速かつ連続的に取得し、時間方向にデータを比較することで心臓の動きを可視化することも可能である。
【０００５】
ここで、撮影のタイミングが異なるボリュームデータ間には一般に未知の位置ずれやひずみが存在するため、正確に比較や合成を行うためにはまず２つのボリュームデータ間のずれやひずみを表すパラメータを知る必要がある。これは、２つのボリュームデータの位置合わせ（レジストレーション）問題として、これまでにランドマークベース（例えば、非特許文献１参照）や幾何学的形状ベース（例えば、非特許文献２参照）、ボクセル類似度ベース（例えば、非特許文献３，４，５，６，７，特許文献１参照）の手法が提案されている。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２００７−５４６３６号公報
【非特許文献】
【０００７】
【非特許文献１】C. R. Maurer Jr., J. M. Fitzpatrick, M. Y. Wang,R. L. Galloeay Jr., R. J. Maciunas, and G. S.Allen. Registration of head volume images usingimplantable fiducial markers. IEEE Trans. MedicalImaging, Vol. 16, No. 4, pp. 447--462, Aug 1997.
【非特許文献２】P. J. Besl and N. D. McKay. A method for registrationof 3-D shapes. IEEE Trans. Pattern Analysisand Machine Intelligence, Vol. 14, No. 2, pp. 239--256,Feb 1992.
【非特許文献３】A. Roche, X. Pennec, G. Malandain, and N. Ayache.Rigid registration of 3-D ultrasound with MR images:A new apploach combining intensity and gradient information.IEEE Trans. Medical Imaging, Vol. 20,No. 10, pp. 1038--1049, Oct 2001.
【非特許文献４】F. Maes, A. Collignon, D. Vandermeulen, G. Marchal,and P. Suetens. Multimodality image registrationby maximization of mutual information. IEEETrans. Medical Imaging, Vol. 16, No. 2, pp. 187--198,Apr 1997.
【非特許文献５】C. Studholme, D. L. G. Hill, and D. J. Hawkes. Aninvariant entropy measure of 3D medical image alignment.Pattern Recognition, Vol. 32, No. 1, pp. 71--86,Apr 1998.
【非特許文献６】D. L. G. Hill, P. G. Batchelor, M. Holden, andD. J. Hawks. Medical image registration. Physicsin Medicine and Biology, Vol. 46, pp. R1--R45, 2001.
【非特許文献７】D. Rueckert, L. I. Sonoda, C. Hayes, D. L. G. Hill,M. O. Leach, and D. J. Hawkes. Nonrigid registrationusing free-form deformation: Application to breastMR images. IEEE Trans. Medical Imaging, Vol. 18,No. 8, pp. 712--721, Aug 1999.
【非特許文献８】K. Takita, T.Aoki, Y. Sasaki, T. Higuchi, andK. Kobayashi. High-accuracy subpixel image registrationbased on phase-only correlation. IEICETrans. Fundamentals, Vol. E86-A, No. 8, pp. 1925--1934, Aug 2003.
【非特許文献９】K. Takita, M. A. Muquit, T. Aoki, and T. Higuchi.A sub-pixel correspondence search technique for computervision applications. IEICE Trans. Fundamentals,Vol. E28-A, No. 8, pp. 1913--1923, 2004.
【非特許文献１０】S. Nagashima, T. Aoki, T. Higuchi, andK. Kobayashi. A subpixel image matching techniqueusing phase-only correlation. Proc. IEEE2006 Int. Symp. Intelligent Signal Processing andCommunication Systems, pp. 701--704, Dec 2006.
【非特許文献１１】K. S. Arun, T. S. Haung, and S. D. Blostein.Least-squares fitting of two 3-D point sets. IEEETrans. Pattern Analysis and Machine Intelligence,Vol. PAMI-9, No. 5, pp. 698--700, Sep 1987.
【非特許文献１２】S. Lee, G. Wolberg, and S. Y. Shin. Scattered datainterpolation with multilevel B-splines. IEEE Trans.Visualization and Computer Graphics, Vol. 3, No. 3,pp. 228--244, Jul 1997.
【非特許文献１３】D. R. Forsey and R. H. Bartels. Hierarchical B-splinerefinement. ACM Trans. Computer Graphics, Vol. 22,No. 4, pp. 205--212, Aug 1988.
【非特許文献１４】M. A. Fischler and R. C. Bolles. Random sampleconsensus: A paradigm for model fitting with applicationsto image analysis and automated cartgraphy.Comm. of the ACM, Vol. 24, No. 6, pp. 381--395, Jun1981.
【非特許文献１５】J. A. Nelder and R. Mead. A simplex method forfunction minimization. Computer Journal, Vol. 7,No. 4, pp. 308--313, 1965.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、ランドマークベースの手法では、埋め込み型マーカなどの位置情報を用いてレジストレーションを行うため、撮影時に侵襲的なマーカが必要であるという問題がある。
【０００９】
幾何学的形状ベースの手法としては、ＩＣＰ（Iterative Closest Point）アルゴリズムが知られている。これは、３次元形状モデル（ポリゴンデータなど）の位置合わせに広く用いられている手法であり、モデル間の対応付けと運動推定を繰り返すことで位置合わせを行う。そのため、まずボリュームから形状モデルを生成する必要があるが、２つのデータ間でモダリティ（様相）が異なるような場合は同一の形状モデルを得ることが難しく、マルチモーダルレジストレーションに適さない。
【００１０】
なお、異なる撮像機器から得られたデータをレジストレーションすることをマルチモーダルレジストレーションと呼んでいる。例えば医療分野では、ＣＴから得られたデータとＭＲＩから得られたデータをレジストレーションするような場合がこれにあたる。また、同じ撮像機器の場合でも、撮影方法が著しく異なるようなデータの場合は、マルチモーダルレジストレーションと呼ぶ場合がある。例えば、ＭＲＩからはＴ１強調画像とＴ２強調画像など、撮影方法によってコントラストが大きく異なった画像が得られる。
【００１１】
ボクセル類似度ベースの手法は、位置合わせパラメータをボリュームデータ間の類似度が最大となるように非線形最適化によって決定する手法であり、医用データのレジストレーションにおいて現在最もよく用いられている。なかでも、非特許文献６に示されているような正規化された相互情報量（Normalized Mutual Information：ＮＭＩ）を類似度として用いる手法が広く知られており、この手法ではボリュームデータ間で計算した相互情報量を最大化することでレジストレーションを行う。なお、特許文献１には、２つのボリュームデータから特徴点を抽出し、この特徴点の近傍のみで相互情報量を最大化することでレジストレーションを行うようにすることが開示されている。
【００１２】
しかし、このような相互情報量の最大化は非線形最適化問題として解かれるため、その膨大な計算時間を要する。特に、３次元のボリュームデータは２次元画像などに比べて大幅に大きなデータ量を持つため、データの解像度によっては位置合わせに数時間を要する場合もある。また、非線形最適化において正しい解を得るためには適切な初期値を設定しなければならず、この初期値を推定する処理も必要となる。
【００１３】
本発明は、このような課題を解決するためになされたもので、その目的とするところは、ボリュームデータ間の高精度な対応付け方法を提供することにある。２つのボリュームデータの間でボクセル同士の対応関係を知ることは、レジストレーションや対象の動きを可視化する際にきわめて有用である。レジストレーションにおいては、ボクセル同士の対応関係からレジストレーションのためのパラメータを推定することができる。本発明においてボリュームデータの対応付け方法に併せて提供されるボリュームデータの位置合わせ方法は、従来技術から精度を低下させることなく大幅な処理時間の削減を達成している。
【課題を解決するための手段】
【００１４】
このような目的を達成するために本発明は、同一対象物を撮像した２つの３次元の立体画像のデータを第１および第２のボリュームデータとして取り込むデータ取込ステップと、第１のボリュームデータから複数の部分領域を抽出し、この抽出した複数の部分領域にそれぞれ設定された基準点に対応する第２のボリュームデータの対応点を、基準点を含む第１のボリュームデータの部分領域と基準点に対応する対応点を含む第２のボリュームデータの部分領域との相関から探索する対応点探索ステップと、探索された対応点と基準点との対応関係から第１のボリュームデータと第２のボリュームデータとの間の回転ずれおよび位置ずれを表すパラメータを剛体変形パラメータとして推定する剛体変形パラメータ推定ステップとを設けたものである。
【００１５】
本発明では、第１のボリュームデータから部分領域が抽出され、この抽出された部分領域に設定された基準点に対応する第２のボリュームデータの部分領域内の対応点が、その部分領域の相関から探索される。同様にして、第１のボリュームデータから複数の部分領域が抽出され、その抽出された部分領域内の基準点に対応する第２のボリュームデータの部分領域内の対応点が探索される。そして、この探索された対応点と基準点との対応関係から第１のボリュームデータと第２のボリュームデータとの間の剛体変形パラメータ（回転ずれおよび位置ずれを表すパラメータ）が推定される。
【００１６】
本発明では、部分領域の相関から対応点を探索するが、この場合の相関法として位相限定相関法（Phase-Only Correlation：ＰＯＣ）を用いることが考えられる。ＰＯＣは、信号を離散フーリエ変換して得られる位相成分のみに着目して信号間の相関を求める手法であり、さまざまな外乱に対してロバストに信号の位置合わせ、類似度評価を行うことができるという特徴がある（例えば、非特許文献８参照）。
【００１７】
ＰＯＣ関数が示す相関ピークは、その理論モデルを入力信号によらず導出することが可能であり、これを用いたフィッティングによって信号間のずれをサンプリング分解能以上の精度で求めることができる。ＰＯＣは、本出願人らによってさまざまな分野への応用が実現されているが、特にサブピクセル精度での画像対応付けが必要となるステレオビジョンなどにおいてその有効性が確認されている（例えば、非特許文献９参照）。
【００１８】
本発明では、その一形態として、これまで１次元および２次元で使用されていたＰＯＣを３次元へ拡張し、この３次元へ拡張したＰＯＣ（３次元ＰＯＣ）を用いた３次元ブロックマッチングによってサブボクセル精度でボリュームデータ間の対応付けを行う。そして、得られた対応関係から、ボリュームデータ間の剛体変形パラメータを推定する。
【００１９】
また、本発明では、その一形態として、推定された剛体変形パラメータに基づいて第１のボリュームデータと第２のボリュームデータとの間の回転ずれおよび位置ずれを補正する。この場合、第１のボリュームデータの回転ずれおよび位置ずれを補正するようにしてもよいし、第２のボリュームデータの回転ずれおよび位置ずれを補正するようにしてもよい。この第１のボリュームデータと第２のボリュームデータとの間の回転ずれおよび位置ずれの補正を「剛体変形」と呼ぶ。すなわち、本発明では、第１のボリュームデータを剛体変形させてもよいし、第２のボリュームデータを剛体変形させてもよい。
【００２０】
また、本発明では、その一形態として、回転ずれおよび位置ずれが補正された第１および第２のボリュームデータを剛体変形処理後の第１および第２のボリュームデータとし、剛体変形処理後の第１のボリュームデータから複数の部分領域を抽出し、この抽出した複数の部分領域にそれぞれ設定された基準点に対応する剛体変形処理後の第２のボリュームデータの対応点を、基準点を含む剛体変形処理後の第１のボリュームデータの部分領域と基準点に対応する対応点を含む剛体変形処理後の第２のボリュームデータの部分領域との相関から探索するようにし、この探索された対応点の情報を用いて剛体変形処理後の第１および第２のボリュームデータの何れか一方の非剛体部分を変形させて重ね合わせるようにしてもよい。この剛体変形処理後の第１および第２のボリュームデータの何れか一方の非剛体部分を変形させて重ね合わせるような補正を「非剛体変形」と呼ぶ。すなわち、本発明では、剛体変形処理後の第１のボリュームデータを非剛体変形させてもよいし、剛体変形処理後の第２のボリュームデータを非剛体変形させてもよい。
【００２１】
本発明は、ボリュームデータの高精度な対応付けを提供する。併せて提供するボクセルの対応関係を利用したボリュームデータの位置合わせ方法は、従来手法で一般的に用いられる複雑な最適化計算を必要としないため、位置合わせに要する計算時間を大幅に削減することができる。また、最適化計算において必須である変形パラメータの初期値を設定する必要もない。さらに、非剛体変形への対応も容易である。
【発明の効果】
【００２２】
本発明によれば、第１のボリュームデータから部分領域を複数抽出し、この抽出した複数の部分領域内の基準点に対応する第２のボリュームデータの部分領域内の対応点を各部分領域の相関から探索することで、ボリュームデータ間の高精度な対応付けが可能となる。
【００２３】
また、本発明によれば、ボリュームデータの対応関係から推定した剛体変形パラメータに基づいて第１のボリュームデータと第２のボリュームデータとの間の回転ずれおよび位置ずれを補正するようにして、複雑な最適化計算を必要とせず、処理時間を短縮し、かつ精度よくボリュームデータ間の位置合わせを行うことが可能となる。
【図面の簡単な説明】
【００２４】
【図１】本発明の実施に用いるボリュームデータ間位置合わせ装置の一例を示すブロック図である。
【図２】このボリュームデータ間位置合わせ装置で行われる特有の処理動作（剛体変形処理）を示すフローチャートである。
【図３】このボリュームデータ間位置合わせ装置における対応点の探索処理（３次元ブロックマッチング）の過程を示すフローチャートである。
【図４】図３に続くフローチャートである。
【図５】このボリュームデータ間位置合わせ装置における対応点の探索処理（３次元ブロックマッチング）の過程を説明するためのボリュームデータの階層構造を示す図である。
【図６】対応点の探索処理に際して用いる３次元ＰＯＣの流れを示す図である。
【図７】このボリュームデータ間位置合わせ装置で行われる特有の処理動作（非剛体変形処理）を示すフローチャートである。
【図８】ＰＯＣに基づく対応点探索におけるパラメータを例示する図である。
【図９】本発明に係る手法と従来手法を用いてレジストレーションを行い、それらの性能を相関係数により比較した結果を示す図である。
【図１０】実験にて使用したボリュームの一部（（ａ），（ｂ））およびレジストレーションの結果（（ｃ），（ｄ））を示すディスプレイ上の写真である。
【図１１】Ｘ線ＣＴデータおよびＭＲＩデータのアキシャル像（（ａ），（ｂ））およびレジストレーションの結果（（ｃ），（ｄ））を示すディスプレイ上の写真である。
【発明を実施するための形態】
【００２５】
以下、本発明の実施の形態を図面に基づいて詳細に説明する。図１は本発明の実施に用いるボリュームデータ間位置合わせ装置の一例を示すブロック図である。
【００２６】
同図において、１００はボリュームデータ間位置合わせ装置であり、プロセッサや記憶装置からなるハードウェアと、これらのハードウェアと協働して各種機能を実現させるプログラムとによって実現される。具体的には、パーソナルコンピュータ（パソコン）にプログラムがインストールされ、このインストールされたプログラムに従うＣＰＵの処理動作として実現される。
【００２７】
ボリュームデータ間位置合わせ装置１００は、参照ボリュームデータ取込部１と、入力ボリュームデータ取込部２と、撮像条件補正部３と、剛体変形処理部４と、非剛体変形処理部５と、ボリュームデータ記憶部６とを備えている。剛体変形処理部４は、オブジェクト領域抽出部４Ａと、対応点探索部４Ｂと、剛体変形パラメータ推定部４Ｃと、剛体変形実行部４Ｄと、反復処理命令部４Ｅとを備え、非剛体変形処理部５は、オブジェクト領域抽出部５Ａと、対応点探索部５Ｂと、非剛体変形実行部５Ｃとを備えている。
【００２８】
この実施の形態において、ボリュームデータ間位置合わせ装置１００には、一例として、参照ボリュームデータＤ１としてＣＴ（図示せず）によって撮影された対象物のスライス画像を再構成した３次元の立体画像のデータが入力され、入力ボリュームデータＤ２としてＭＲＩ（図示せず）によって撮影された対象物のスライス画像を再構成した３次元の立体画像のデータが入力されるものとする。
【００２９】
ここで、本実施の形態においては、位置合わせの基準となるボリュームデータを「参照ボリュームデータ」と呼ぶこととし、この参照ボリュームデータを基準として入力ボリュームデータの位置合わせを行うものとする。なお、この例において、撮影された対象物は同一人の頭部であるとする。
【００３０】
また、本実施の形態において、剛体変形処理とは、撮影対象物自体は変形しないものとして、２つのボリュームデータ間の移動、回転によるずれなどを修正する処理のことを意味する。また、非剛体変形処理とは、撮像対象物が例えば心臓のように、時間的に変形することによって生じる２つのボリュームデータ間のずれを修正する処理のことを意味する。
【００３１】
図２におよび図７にこのボリュームデータ間位置合わせ装置１００で行われる特有の処理動作のフローチャートを示す。以下、このフローチャートを参照しながら、ボリュームデータ間位置合わせ装置１００における各部の処理動作について説明する。
【００３２】
〔剛体変形処理〕
先ず、図２に示すフローチャートを用いて、ボリュームデータ間位置合わせ装置１００で行われる剛体変形処理について説明する。
【００３３】
〔ボリュームデータの取り込み〕
参照ボリュームデータ取込部１は、ＣＴによって撮影され再構成された３次元の立体画像のデータを参照ボリュームデータＤ１として取り込み、この取り込んだ参照ボリュームデータＤ１を撮像条件補正部３へ送る（ステップＳ１０１）。
【００３４】
入力ボリュームデータ取込部２は、ＭＲＩによって撮影され再構成された３次元の立体画像のデータを入力ボリュームデータＤ２として取り込み、この取り込んだ入力ボリュームデータＤ２を撮像条件補正部３へ送る（ステップＳ１０２）。
【００３５】
〔撮像条件の補正〕
撮像条件補正部３は、参照ボリュームデータ取込部１からの参照ボリュームデータＤ１と入力ボリュームデータ取込部２からの入力ボリュームデータＤ２との間の撮像条件に起因する相違を補正する（ステップＳ１０３）。この例では、撮像条件に起因する相違として、参照ボリュームデータＤ１と入力ボリュームデータＤ２との間の倍率や大まかな回転ずれや位置ずれを補正する。補正された参照ボリュームデータＤ１はボリュームデータ記憶部６に保存される。また、補正された参照ボリュームデータＤ１および入力ボリュームデータＤ２は、剛体変形処理部４へ送られる。
【００３６】
なお、撮像条件補正部３において、倍率については装置が持つスケール（１ボクセルと実座標上の大きさ〔ｍ〕の対応）を利用して合わせるなどする。回転ずれについては、それぞれのデータに主成分を計算して合わせるなどする。位置ずれについては、それぞれのデータの重心を利用して合わせるなどする。
【００３７】
〔オブジェクト領域の抽出〕
剛体変形処理部４において、オブジェクト領域抽出部４Ａは、撮像条件補正部３で補正された参照ボリュームデータＤ１および入力ボリュームデータＤ２を取り込み、この参照ボリュームデータＤ１および入力ボリュームデータＤ２から所定の領域をオブジェクト領域としてそれぞれ抽出する（ステップＳ１０４）。
【００３８】
この例では、参照ボリュームデータＤ１および入力ボリュームデータＤ２に含まれる背景を除去し、対象物に関する情報が含まれる領域のみをオブジェクト領域として抽出する。以下、この抽出されたオブジェクト領域のボリュームデータを参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２とする。
【００３９】
なお、ここでの背景とは、ＣＴやＭＲＩなどから得られる医用ボリュームデータで言えば、対象の周りに存在する空気を指す。空気は、対象とは全く異なるボクセル値を持つので、ボクセル値に閾値処理を施すことで容易に対象と空気を分離することができる。
【００４０】
〔対応点の探索（３次元ブロックマッチング）〕
対応点探索部４Ｂは、オブジェクト領域抽出部４Ａからの参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２を取り込み、この参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２間の対応点の探索を行う（ステップＳ１０５）。図３および図４にこの場合の対応点の探索処理のフローチャートを示す。
【００４１】
対応点探索部４Ｂは、先ず、参照ボリュームデータＯＢ１にＮ個の基準点（ｘｉ，ｙｉ，ｚｉ：ｉ＝１，２，・・・・，Ｎ）を定める（ステップＳ２０１（図３））。また、参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２をリサイズして行くことで、各辺の長さ（ｘ，ｙ，ｚ軸方向の長さ）が元のデータの２−ｌとなる一連の階層ボリュームデータを生成する（ステップＳ２０２）。
【００４２】
なお、本実施の形態のリサイズでは、隣接する８ボクセルを平均化することで１ボクセルとし、ボリュームデータの各辺をそれぞれ２分の１の長さにする。しかし、必ずしもこのような方法をとる必要はなく、他にも、例えば平滑化フィルタをかけた後に単純にダウンサンプリング（１ボクセルおきにサンプリング）する方法などが考えられる。
【００４３】
すなわち、参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２を段階的に縮小し、元の参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２を最下位階層とし、最も縮小された参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２を最上位階層とする複数階層（ｌ＝０，１，２，・・・・ｌｍａｘ）のボリュームデータを生成する。以下、参照ボリュームデータＯＢ１の各階層のボリュームデータを第１群のボリュームデータ、入力ボリュームデータＯＢ２の各階層のボリュームデータを第２群のボリュームデータと呼ぶ。
【００４４】
そして、ｉ＝１とし（ステップＳ２０３）、第１群の最上位階層ｌｍａｘのボリュームデータに初期基準点（floor（ｘｉ／２lmax），floor（ｙｉ／２lmax），floor（ｚｉ／２lmax））をｐｌｍａｘとして設定し（図５参照）、第２群の最上位階層ｌｍａｘのボリュームデータの第１群の最上位階層ｌｍａｘに設定した初期基準点ｐｌｍａｘに対応する点に初期対応点をｑｌｍａｘとして設定する（ステップＳ２０４）。この際、最上位階層ｌｍａｘにおいては基準点と対応点とが同じ座標を持つと仮定し、初期基準点ｐｌｍａｘ＝初期対応点ｑｌｍａｘとする。
【００４５】
次に、ｌ＝ｌｍａｘ−１とし（ステップＳ２０５）、第１群の階層ｌ＝ｌｍａｘ−１のボリュームデータに初期基準点ｐｌｍａｘに対応する基準点（floor（ｘｉ／２lmax-1），floor（ｙｉ／２lmax-1），floor（ｚｉ／２lmax-1））をｐｌｍａｘ−１として定め、この基準点ｐｌｍａｘ−１を中心とする３辺の長さが予め定められた３次元のブロック（部分領域）を探索ウィンドウＷ１として設定する（ステップＳ２０６：図５参照）。
【００４６】
基準点（floor（ｘｉ／２lmax-1），floor（ｙｉ／２lmax-1），floor（ｚｉ／２lmax-1））における「floor」の意味について説明する。Ｂ=floor(Ａ)は、ＡをＡ以下の最も近い整数とする。ボリュームデータに設定した基準点（ｘｉ， ｙｉ， ｚｉ）は、階層探索を行う際に解像度を変更したボリュームデータに合わせて（ｘｉ／２lmax，ｙｉ／２lmax，ｚｉ／２lmax）とする必要があるが、この際に必ずしも座標が整数とはならない。このため、本実施の形態では、整数とするために、floor関数を用いている。
【００４７】
また、第２群の階層ｌ＝ｌｍａｘ−１のボリュームデータに初期対応点ｑｌｍａｘに対応する対応候補点（最上位階層ｌｍａｘに設定した初期対応点ｑｌｍａｘを２倍した座標）をｑｌｍａｘ−１’として定め、この対応候補点ｑｌｍａｘ−１’を中心とする３辺の長さが予め定められた３次元のブロック（部分領域）を探索ウィンドウＷ２として設定する（ステップＳ２０７（図４）：図５参照）。なお、探索ウインドウＷ１，Ｗ２の大きさは同じであり、以下に説明する各階層でも同じ大きさの探索ウィンドウＷ１，Ｗ２を用いるものとする。
【００４８】
そして、対応点探索部４Ｂは、ステップＳ２０６で設定した基準点ｐｌｍａｘ−１を中心とする探索ウィンドウＷ１とステップＳ２０７で設定した対応候補点ｑｌｍａｘ−１’を中心とする探索ウィンドウＷ２のボクセルデータを抽出し、例えば後述する位相限定相関法（３次元ＰＯＣ）を用いてウインドウ間のずれをボクセルレベルで求め、このずれから探索ウインドウＷ２におけるボクセルレベルの対応点ｑｌｍａｘ−１を決定する（ステップＳ２０８）。すなわち、第２群の階層ｌｍａｘ−１のボリュームデータにおける対応点ｑｌｍａｘ−１を決定する。
【００４９】
図６にステップＳ２０８で行われる３次元ＰＯＣの処理の流れを示す。対応点探索部４Ｂは、基準点ｐｌｍａｘ−１を中心とする探索ウィンドウＷ１のボクセルデータを抽出し、この抽出したボクセルデータに対して３次元フーリエ変換を行う（ステップＳ３０１）。また、対応候補点ｑｌｍａｘ−１’を中心とする探索ウィンドウＷ２のボクセルデータを抽出し、この抽出したボクセルデータに対して３次元フーリエ変換を行う（ステップＳ３０２）。
【００５０】
そして、この３次元フーリエ変換が施された探索ウィンドウＷ１の３次元スペクトルと探索ウィンドウＷ２の３次元スペクトルとを合成し（ステップＳ３０３）、この合成した３次元スペクトルの振幅を正規化して位相成分のみとし（ステップＳ３０４）、この位相成分のみとした３次元スペクトルに３次元逆フーリエ変換を施して（ステップＳ３０５）、２つの探索ウィンドウＷ１，Ｗ２の相関を表す位相限定相関関数を得る。
【００５１】
この位相限定相関関数から階層ｌ＝ｌｍａｘ−１のボリュームデータの空間内における探索ウィンドウＷ１とＷ２との間のずれをボクセルレベルで求め、このずれから探索ウィンドウＷ２におけるボクセルレベルの対応点ｑｌｍａｘ−１を決定する。なお、３次元ＰＯＣの詳細については後述する。
【００５２】
次に、対応点探索部４Ｂは、ｌ＝ｌ−１とし（ステップＳ２０９）、すなわちｌ＝ｌｍａｘ−２とし、ステップＳ２０６〜Ｓ２０８の処理を繰り返す。
【００５３】
この場合、対応点探索部４Ｂは、第１群の階層ｌｍａｘ−２のボリュームデータに基準点ｐｌｍａｘ−１に対応する基準点（floor（ｘｉ／２lmax-2），floor（ｙｉ／２lmax-2），floor（ｚｉ／２lmax-2））をｐｌｍａｘ−２として定め、この基準点ｐｌｍａｘ−２を中心とする探索ウィンドウ（３次元ブロック）Ｗ１を設定する（ステップＳ２０６：図５参照）。
【００５４】
また、第２群の階層ｌ＝ｌｍａｘ−２のボリュームデータに対応点ｑｌｍａｘ−１に対応する対応候補点（階層ｌｍａｘ−１で求めた対応点ｑｌｍａｘ−１を２倍した座標）をｑｌｍａｘ−２’として定め、この対応候補点ｑｌｍａｘ−２’を中心とする探索ウィンドウ（３次元ブロック）Ｗ２を設定する（ステップＳ２０７：図５参照）。
【００５５】
そして、基準点ｐｌｍａｘ−２を中心とする探索ウィンドウＷ１と対応候補点ｑｌｍａｘ−２’を中心とする探索ウィンドウＷ２のボクセルデータを抽出し、３次元ＰＯＣを用いてウインドウ間のずれをボクセルレベルで求め、このずれから探索ウインドウＷ２におけるボクセルレベルの対応点ｑｌｍａｘ−２を決定する（ステップＳ２０８）。すなわち、第２群の階層ｌｍａｘ−２のボリュームデータにおける対応点ｑｌｍａｘ−２を決定する。
【００５６】
以下、同様にして、対応点探索部４Ｂは、第１群の階層および第２群の階層を１段階ずつ下げながら、ステップＳ２０６〜２０８の処理を繰り返して行く。そして、最下位階層におけるボクセルレベルの対応付けが終了し、ｌ＝−１となると（ステップＳ２１０のＹＥＳ）、対応点探索部４Ｂは、参照ボリュームデータＯＢ１と入力ボリュームデータＯＢ２において、基準点ｐ０を中心とする探索ウィンドウＷ１のボクセルデータとｌ＝０の階層においてステップＳ２０６〜Ｓ２０８を実行することによって得られたボクセルレベルの対応候補点ｑ０を中心とする探索ウインドウＷ２のボクセルデータを抽出する（ステップＳ２１１）。
【００５７】
そして、３次元ＰＯＣを用いてウインドウ間のずれをサブボクセルレベル（ボクセルレベルよりもさらに細かいレベル）で求め、このずれから探索ウインドウＷ２における対応点ｑ０をサブボクセルレベルへ更新し、この更新したサブボクセルレベルの対応点ｑ０を基準点ｐ０に対応する対応点として記憶する（ステップＳ２１２）。すなわち、参照ボリュームデータＯＢ１に定められるｉ＝１番目の基準点ｐ０に対応する入力ボリュームデータＯＢ２における正規の対応点ｑ０として記憶する。
【００５８】
そして、対応点探索部４Ｂは、ｉ＝ｉ＋１とし（ステップＳ２１４）、すなわちｉ＝２とし、２番目の基準点について、ステップＳ２０４〜Ｓ２１２の処理を繰り返す。これにより、上述と同様にして、参照ボリュームデータＯＢ１に定められる基準点ｐ０に対応する入力ボリュームデータＯＢ２における対応点ｑ０が決定され、この決定された対応点ｑ０がｉ＝２番目の基準点に対応する正規の対応点として記憶される。この処理は、ステップＳ２１３においてｉ＝Ｎとなるまで繰り返される。
【００５９】
これにより、参照ボリュームデータＯＢ１におけるＮ個の基準点について、各基準点に対応する入力ボリュームデータＯＢ２における正規の対応点が次々に決定されて行く。
【００６０】
本実施の形態では、上述したように、ボリュームデータの解像度を変えたものをいくつか生成して階層状にし、低い解像度のデータから高い解像度のデータへ順番に対応点探索を行っていくことで、効率の良い対応付けが可能となり、また、データ間で大きな動きが生じている場合にも対応することができるようになる。
【００６１】
〔剛体変形パラメータの推定〕
剛体変形パラメータ推定部４Ｃは、対応点探索部４Ｂによって探索された対応点と基準点との対応関係から参照ボリュームデータＯＢ１と入力ボリュームデータＯＢ２との間の回転ずれおよび位置ずれを表すパラメータを剛体変形パラメータとして推定する（ステップＳ１０６（図２））。
【００６２】
例えば、回転ずれとしてＸ，Ｙ，Ｚ軸の各軸の回転ずれ量Δθｘ，Δθｙ，Δθｚを求め、位置ずれとして各軸の平行ずれ量Δｘ，Δｙ，Δｚを求める。この例では、回転ずれとして回転行列を求め、位置ずれとして並進ベクトルを求める。剛体変形パラメータの推定は、たとえば後述するように最小二乗法を用いて行うことができる。
【００６３】
なお、この剛体変形パラメータ推定部４Ｃにおいて剛体変形パラメータを推定する場合、例えば非特許文献１４に示されているようなロバスト推定法を用いるようにすれば、対応関係に誤りが含まれていても高精度にパラメータの推定を行うことが可能である。
【００６４】
〔剛体変形の実行〕
剛体変形実行部４Ｄは、剛体変形パラメータ推定部４Ｃによって推定された剛体変形パラメータに基づいて、撮影条件補正部３から送られてくる入力ボリュームデータＤ２の回転ずれおよび位置ずれを補正（剛体変形）する（ステップＳ１０７）。すなわち、入力ボリュームデータＤ２を剛体変形させることによって、ボリュームデータ記憶部６に記憶されている参照ボリュームデータＤ１との間の回転ずれおよび位置ずれを補正する。
【００６５】
そして、この剛体変形させた入力ボリュームデータＤ２をＤ２’として、反復処理命令部４Ｅへ送る。なお、剛体変形実行部４Ｄは、剛体変形後の入力ボリュームデータＤ２’を、次の回の剛体変形前の入力ボリュームデータＤ２として記憶する。
【００６６】
〔反復処理〕
反復処理命令部４Ｅは、剛体変形実行部４Ｄからの剛体変形された入力ボリュームデータＤ２’をオブジェクト領域抽出部４Ａへ送り、オブジェクト領域抽出部４Ａでのオブジェクト領域の抽出（ステップＳ１０４）、対応点探索部４Ｂでの対応点の探索（ステップＳ１０５）、剛体変形パラメータ推定部４Ｃでの剛体変形パラメータの推定（ステップＳ１０６）、剛体変形実行部４Ｄでの剛体変形の実行（ステップＳ１０７）を再度行わせる。
【００６７】
反復処理命令部４Ｅは、このステップＳ１０４〜Ｓ１０７の処理を繰り返させ、その回数が所定回数に達した時点で（ステップＳ１０８のＹＥＳ）、その時の補正された入力ボリュームデータＤ２’を剛体変形処理後の入力ボリュームデータとしてボリュームデータ記憶部６に記憶させる（ステップＳ１０９）。
【００６８】
このような反復処理により、剛体変形パラメータの推定とこの推定された剛体変形パラメータを用いての入力ボリュームデータの剛体変形が繰り返され、入力ボリュームデータと参照ボリュームデータとの間で大きな移動や回転が生じている場合にも対応することができるようになる。
【００６９】
〔非剛体変形処理〕
次に、図７に示すフローチャートを用いて、ボリュームデータ間位置合わせ装置１００で行われる非剛体変形処理について説明する。
【００７０】
〔オブジェクト領域の抽出〕
非剛体変形処理部５において、オブジェクト領域抽出部５Ａは、ボリュームデータ記憶部６に記憶されている参照ボリュームデータＤ１を読み出し（ステップＳ４０１）、またボリュームデータ記憶部６に記憶されている剛体変形処理後の入力ボリュームデータＤ２’を読み出し（ステップＳ４０２）、この参照ボリュームデータＤ１および入力ボリュームデータＤ２’からオブジェクト領域を抽出し、参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２とする（ステップＳ４０３）。
【００７１】
なお、この場合、参照ボリュームデータＤ１に対しては剛体変形処理部４での剛体変形は実施されていないが、入力ボリュームデータＤ２に剛体変形を実施することによって参照ボリュームデータＤ１との間の回転ずれおよび位置ずれが補正されており、両者に対して補正が行われているものと考えれば、参照ボリュームデータＤ１も剛体変形処理後のボリュームデータであると言える。ここでは、参照ボリュームデータＤ１も剛体変形処理後のボリュームデータと考える。
【００７２】
〔対応点の探索〕
対応点探索部５Ｂは、オブジェクト領域抽出部５Ａからの参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２を取り込み、この参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２間の対応点の探索を行う（ステップＳ４０４）。この対応点探索部５Ｂでの対応点の探索処理は、剛体変形処理部４における対応点探索部４Ｂでの対応点の探索処理（図３）と同じであるので、その説明は省略する。
【００７３】
非剛体変形実行部５Ｃは、対応点探索部５Ｂによって探索された対応点の情報を用いて、ＦＦＤ（Free-Form Deformation）により入力ボリュームデータＤ２’を補正（非剛体変形）し、剛体変形処理後の入力ボリュームデータＤ２’と参照ボリュームデータＤ１の非剛体部分を重ね合わせる（ステップＳ４０５）。そして、この補正された入力ボリュームデータＤ２’を非剛体変形処理後の入力ボリュームデータＤ２”として、ボリュームデータ記憶部６に記憶させる（ステップＳ４０６）。ＦＦＤによる非剛体変形処理の詳細については後述する。
【００７４】
なお、ＦＦＤ（Free-Form Deformation）については、例えば非特許文献７にも開示されているように、ボリュームデータの非剛体位置合わせにおいては一般的である。本発明によりあらかじめ得られたボクセルの対応関係を利用したＦＦＤによる非剛体変形は、従来技術であるボクセル類似度ベースの手法を用いたＦＦＤに比べて約１００倍の計算速度向上が見られた。従来技術ではボリュームデータの対応関係が未知の状態で計算を始めるため特に計算量の多い非剛体変形では非線形最適化に膨大な時間がかかるが、あらかじめ本発明により高精度にボリュームデータの対応関係が得られていれば、大幅な計算速度向上が可能である。
【００７５】
なお、上述した実施の形態では、剛体変形処理部４の対応点探索部４Ｂにおいて対応点を探索する場合、参照ボリュームデータＯＢ２の全体に基準点を定めるようにするが、例えばＣＴやＭＲＩなどの医用ボリュームデータを扱う場合には、骨などの形状変化が少ない部分に集中的に定めることで、位置合わせ精度を向上させることが可能である。
【００７６】
例えば、参照ボリュームデータＯＢ２について、骨のある領域を所定の特徴のある領域として検索させ、この検索された特徴のある領域に基準点を定めるようにする。すなわち、検索された特徴のある領域を探索ウィンドウの抽出領域とする。非剛体変形処理部５の対応点探索部５Ｂにおいて対応点を探索する場合も同様である。
【００７７】
また、この実施の形態では、オブジェクト領域抽出部４Ａや５Ａにおいて背景を除去して対象物のある領域のみを抽出するようにしたが、上述した骨のある領域など特徴のある領域をオブジェクト領域として抽出するようにしてもよい。また、撮像条件補正部３やオブジェクト領域抽出部４Ａ，５Ａは必ずしも設けなくてもよい。
【００７８】
また、上述した実施の形態では、対応点探索部４Ｂや５Ｂにおいて、参照ボリュームデータＯＢ１および入力ボリュームデータＯＢ２の階層は、複数階層であれば何階層であってもよく、２層としても構わない。例えば、例えば図５に示した例では、最下位階層ｌ０を０層とした場合、ｌｍａｘは３層であるが、ｌｍａｘを２層としてもよい。また、ｌｍａｘを４層、５層と、その層を増やして行ってもよい。
【００７９】
また、上述した実施の形態では、対応点の探索を複数階層で行うようにしたが、必ずしも複数階層を利用する必要はなく、対象にそれほど大きなずれがない場合は１層でも対応付けは可能である。
【００８０】
また、上述した実施の形態では、剛体変形処理部４におけるオブジェクト領域抽出部４Ａと非剛体変形処理部５におけるオブジェクト領域抽出部４Ａとを別としたが、また、剛体変形処理部４における対応点探索部４Ｂと非剛体変形処理部５における対応点探索部４Ｂとを別としたが、共通のブロックとすることも可能である。
【００８１】
また、上述した実施の形態では、対応点探索部４Ｂや５Ｂにおいて、３次元ＰＯＣを用いてウインドウ間のずれを求めるようにしたが、必ずしも３次元ＰＯＣ（位相限定相関法）を用いなくてもよく、他の相関法を用いてもよい。また、位相に限定するのではなく、ｌｏｇ処理やルート処理などによって振幅を抑制する振幅抑制相関法を用いてもよい。
【００８２】
また、上述した実施の形態では、剛体変形処理部４における剛体変形実行部４Ｄにおいて、入力ボリュームデータＤ２に対して剛体変形を実施するようにしたが、参照ボリュームデータＤ１に対して剛体変形を実施するようにしてもよい。また、非剛体変形処理部５における非剛体変形実行部５Ｃにおいて、剛体変形処理後の入力ボリュームデータＤ２’に対して非剛体変形を実施するようにしたが、参照ボリュームデータＤ１（剛体変形処理後の参照ボリュームデータＤ１）に対して非剛体変形を実施するようにしてもよい。
【００８３】
また、上述した実施の形態では、対応点探索部４Ｂや５Ｂにおいて対応点を探索する場合、参照ボリュームデータＯＢ２に基準点を複数定め、その基準点を中心する探索ウィンドウを抽出するようにしたが、探索ウィンドウを１つずつ切り出し、その探索ウィンドウの中心を基準点として、対応点を探索して行くようにしてもよい。
【００８４】
また、上述した実施の形態では、剛体変形処理部４の反復処理命令部４Ｅにおいて、所定回数に達するまでステップＳ１０４〜Ｓ１０７の処理を繰り返させるようにしたが、補正された入力ボリュームデータＤ２’と参照入力ボリュームデータＤ１との回転ずれおよび位置ずれが所定の閾値よりも小さくなった場合に、反復処理を停止し、その時の補正された入力ボリュームデータＤ２’を剛体変形処理後の入力ボリュームデータとしてボリュームデータ記憶部６に記憶させるようにしてもよい。またこの他にも、推定される剛体変形パラメータの変化率が所定の閾値よりも小さくなった場合に反復処理を停止してもよい。
【００８５】
また、上述した実施の形態では、参照ボリュームデータＤ１をＣＴによって撮像されたものとし、入力ボリュームデータＤ２をＭＲＩによって撮像されたものとしたが、参照ボリュームデータＤ１と入力ボリュームデータＤ２をともにＣＴによって撮像されたものとしたり、ＭＲＩによって撮像されたものとしてもよく、さらにこれらとも異なる別の撮像装置で撮像されたデータとしてもかまわない。
【００８６】
次に、３次元ＰＯＣを用いたボリュームデータ間の位置合わせ手法について、さらに詳細に説明する。
【００８７】
〔３次元位相限定相関法〕
ここでは３次元ＰＯＣの基本的な定義と、３次元ＰＯＣを用いたボリュームマッチングの高精度化について述べる。
【００８８】
〔３次元位相限定相関関数〕
Ｎ1×Ｎ2×Ｎ3ボクセルの２つのボリュームデータｆ（ｎ1，ｎ2，ｎ3）およびｇ（ｎ1，ｎ2，ｎ3）が与えられたとする。ここで、ボリュームデータの離散空間インデックス（整数）を、便宜上、ｎ1＝−Ｍ1、・・・、Ｍ1、ｎ2＝−Ｍ2、・・・、Ｍ2、ｎ3＝−Ｍ3、・・・、Ｍ3とする。ただし、Ｍ1、Ｍ2、Ｍ3は正の整数であり、Ｎ1＝２Ｍ1＋１、Ｎ2＝２Ｍ2＋１、Ｎ3＝２Ｍ3＋１である。
【００８９】
なお、ここでは説明を簡単にするために離散空間のインデックスを正負対称にとり、かつボリュームデータの大きさＮ1，Ｎ2，Ｎ3を奇数にしているが、これは必須ではない。すなわち、通常よく用いられるように非負のインデックスを用い、Ｎ１，Ｎ2，Ｎ3を任意の正の整数に設定するように一般化することが可能である。ボリュームｆ（ｎ1，ｎ2，ｎ3）およびｇ（ｎ1，ｎ2，ｎ3）の３次元離散フーリエ変換（Discrete Fourie Transform：ＤＦＴ）Ｆ（ｋ1，ｋ2，ｋ3）およびＧ（ｋ1，ｋ2，ｋ3）をそれぞれ次式で定義する。
【００９０】
【数１】
【００９１】
３次元ＰＯＣ関数ｒ（ｎ1，ｎ2，ｎ3）を、正規化相互パワースペクトルの３次元逆離散フーリエ変換（Inverse Discrete Fourier Transform：ＩＤＦＴ）として定義する。
【００９２】
【数２】
【００９３】
すると、（４）式に示す３次元ＰＯＣ関数ｒ（ｎ1，ｎ2，ｎ3）は、２つのボリュームデータｆ（ｎ1，ｎ2，ｎ3）とｇ（ｎ1，ｎ2，ｎ3）の３次元の相互相関を表すことになり、２つのボリュームが類似している場合、３次元ＰＯＣ関数は、デルタ関数に近いきわめて鋭いピークを有する。この相関ピークの高さは２つのボリュームの位相差スペクトルθF（ｋ1，ｋ2，ｋ3）−θR（ｋ1，ｋ2，ｋ3）の線形性を表しており、位相差スペクトルが完全に線形であれば、相関ピークの高さは１となる。この相関ピークの高さはボリュームの類似度の尺度として有用である。一方、相関ピークの座標は２つのボリュームの相対的な位置ずれに対応している。
【００９４】
以下では、ｇ（ｎ1，ｎ2，ｎ3）がｆ（ｎ1，ｎ2，ｎ3）を（δ1，δ2，δ3）だけ微小に平行移動させたボリュームである場合を考える。ここで、δ1，δ2，δ3は、それぞれｎ1，ｎ2，ｎ3方向のサブボクセルレベルの移動量を表している。このとき、ｆ（ｎ1，ｎ2，ｎ3）とｇ（ｎ1，ｎ2，ｎ3）の３次元ＰＯＣ関数は、次式で与えられる。
【００９５】
【数３】
【００９６】
ここで、α＝１である。上式は、ボリュームが（δ1，δ2，δ3）だけ微小に平行移動した場合の３次元ＰＯＣ関数の一般形を表している。αは、相関ピークの高さを表現するために導入されたパラメータである。ボリュームに対して無相関なノイズが加わるとαの値が減少するため、実際にはα≦１となる。この相関ピークのモデルに基づく関数フィッティングにより、パラメータα，δ1，δ2，δ3を推定することで、ボリュームの類似度（位相差スペクトルの線形性）とサブボクセル精度の移動量を求めることができる。
【００９７】
〔３次元ＰＯＣに基づくボリュームマッチングの高精度化〕
以下（Ａ）〜（Ｃ）では、３次元ＰＯＣに基づくボリュームマッチングの高精度化について述べる。
（Ａ）窓関数の適用
ＤＦＴでは、信号が周期的に循環することを仮定するため、画像端での信号の不連続性が問題となる。この不連続性の影響を軽減するため、ボリュームに対して窓関数を適用することが重要である。本手法では下記のハニング窓を用いる。
【００９８】
【数４】
【００９９】
（Ｂ）スペクトルの重み付け
一般に、自然画像のエネルギーは低周波領域に集中し、高周波成分のエネルギーは相対的に小さいことが知られている。このため、エイリアシング、ぼけ、雑音、歪みなどの外乱が加わると、高周波成分のＳ／Ｎが大幅に劣化する。そこで、信頼性の低い高周波成分の影響を抑制するために、正規化相互パワースペクトルＲ（ｋ1，ｋ2，ｋ3）の計算の際に、低域通過型のスペクトル重み付け関数Ｈ（ｋ1，ｋ2，ｋ3）を適用することにより、大幅な精度向上が可能である。本手法ではＨ（ｋ1，ｋ2，ｋ3）として、次式で与えられるガウス関数を用いる。
【０１００】
【数５】
【０１０１】
この場合、３次元ＰＯＣ関数は、Ｈ（ｋ1，ｋ2，ｋ3）Ｒ（ｋ1，ｋ2，ｋ3）の３次元ＩＤＦＴであり、式（５）に対応する相関ピークのモデルは、次式のように変化する。
【０１０２】
【数６】
【０１０３】
なお、重み付け関数については、これら以外にも用途に応じて、任意の関数を適用することが可能である［非特許文献８参照］。
【０１０４】
（Ｃ）ピークモデルのフィッティング
一般に、移動量（δ1，δ2，δ3）は実数値であり、３次元ＰＯＣ関数のピーク座標がサンプリング格子点の間に存在するため、正確に移動量を推定することが困難である。そこで、相関ピークのモデルが式（５）、（９）で与えられることを考慮し、実際に計算されたＰＯＣ関数の数値データに対して本モデルをフィッティングすることで、ボクセル間に存在するピークの位置を推定する。
【０１０５】
上記（Ｂ）のスペクトル重み付け関数Ｈ（ｋ1，ｋ2，ｋ3）を適用した場合は、式（９）の例のように、Ｈ（ｋ1，ｋ2，ｋ3）に応じた相関ピークモデルをフィッティングする必要がある。このとき、α，δ1，δ2，δ3がフィッティングパラメータとなる。
【０１０６】
なお、式（５）の形の相関ピークモデルの場合は、ＰＯＣ関数の数値データから直接的に相関ピークの座標と値を求めるピーク評価式（Peak Evaluation Formula：ＰＥＦ）が導出されており［非特許文献１０参照］、リアルタイム処理に寄与している。それ以外の一般的な相関ピークモデルの場合は、Levenberg-Marquardt法などの非線形最小二乗法を用いたフィッテイングを行う。
【０１０７】
このように、３次元ＰＯＣ関数の相関ピーク形状は、入力ボリュームによらず、Ｈ（ｋ1，ｋ2，ｋ3）のみによって決定されるという特徴があるため、関数フィッティングによって高精度な移動量推定が可能である。
【０１０８】
〔３次元位相限定相関法に基づくボリュームレジストレーション〕
本手法で提案するボリュームレジストレーション手法の詳細について述べる。提案アルゴリズムは２つのステップで構成され、まず、（i）３次元ＰＯＣに基づいてサブボクセルレベルの対応点探索を行い、次に（ii）ボクセルの対応関係からボリュームデータの変形パラメータを推定する。
【０１０９】
まず、ボリュームの変形モデルとして剛体変形を用い、剛体変形を記述するパラメータを推定する。本手法は非剛体変形への拡張も容易であり、本節ではＦＦＤを用いた非剛体レジストレーションについても述べる。以下では本手法における各処理の詳細を述べる。
【０１１０】
〔サブボクセル対応点探索〕
本手法では、まず２つのボリュームデータから小領域（部分領域（３次元ブロック））を切り出し、３次元ＰＯＣを用いてブロック間の移動量をサブボクセル精度で求めることで対応付けを行う。このとき、粗密戦略に基づく階層探索を行うことで効率的な対応付けを実現している。以下では、ボリュームＩ上に設定した基準点ｐ＝（ｐ1，ｐ2，ｐ3）に対応するＪ上の対応点ｑ＝（ｑ1，ｑ2，ｑ3）を求めるものとし、対応付け手法の詳細について述べる。
【０１１１】
【数７】
【０１１２】
【数８】
【０１１３】
以上に示した処理Ｓｔｅｐ１〜Ｓｔｅｐ６によって、基準点ｐの対応点ｑをサブボクセル精度で求めることができる。
【０１１４】
〔剛体変形パラメータの推定〕
前述のブロックマッチングで得られた対応関係から剛体変形パラメータを推定する。剛体変形を記述するパラメータは、前述の対応点が３点以上あれば推定が可能である。本手法では、行列の特異値分解に基づいて回転行列Ｒおよび並進ベクトルｔの最小二乗解を求めるアルゴリズム［非特許文献１１参照］を用いる。以下にその手順を示す。
【０１１５】
【数９】
【０１１６】
〔反復推定〕
ボリューム間に大きな回転などが含まれる場合、一度の対応点探索及び剛体変形パラメータの推定だけでは高精度な位置合わせが難しいことが多い。そこで、まず前述の方法で推定された剛体変形パラメータを利用し、ボリュームデータ間の変形が小さくなるように剛体変形を施す。
【０１１７】
そして、再び対応付けを行い、剛体変形パラメータを推定する。このように、対応付け、剛体変形パラメータの推定、ボリュームの変形を繰り返すことで、ボリューム間で大きな変形が生じている場合でも正確に位置合わせを行うことが可能となる。実際には、例えばボリュームが一方向に２０°程度回転している場合、約１０回の繰り返し推定で誤差が収束することを確認している。
【０１１８】
〔非剛体レジストレーション〕
医用ボリュームデータにおいて、心臓や肺といった臓器に加え関節や筋肉などの組織は、それ自体が変形するなど、剛体変形では表現が不可能な変形を伴う場合があるため、剛体レジストレーションだけでは不十分であるといえる。したがって、剛体レジストレーションに加え、局所的な変形を補正するために非剛体レジストレーションを行う必要がある。
【０１１９】
非剛体レジストレーションでは、一度剛体レジストレーションを行った後、再度対応点探索を行い、この対応関係を用いて非剛体な変形パラメータを推定する。以下では、各ボクセルの移動をFree-Form Deformation（ＦＦＤ）モデルにより推定する手法について述べる。
【０１２０】
本手法では、ローカルな変形モデルとしてＢ-スプライン補間に基づくＦＦＤ［非特許文献１２参照］を用いる。ＦＦＤは３次元変形物体のモデリングに有用であり、これまで主に人体の運動の追跡や解析に利用されてきた。Ｂ-スプライン補間に基づくＦＦＤでは、ボリュームデータの３次元領域内に格子状に制御点が配置され、この制御点の移動がボリュームの変形を決定付ける。ｘ，ｙ，ｚ方向にそれぞれ間隔δｘ，δｙ，δｚで格子状の制御点φｉ，ｊ，ｋが並ぶとき、任意の座標ｘ＝（ｘ，ｙ，ｚ）における移動ベクトルＴ（ｘ）は次式で表される。
【０１２１】
【数１０】
【０１２２】
ここで入力ボリュームをＩ，Ｊとすると、ｐｉ，ｊ，ｋはＩ上の基準点であり、ｑｉ，ｊ，ｋは３次元ＰＯＣを用いた対応付けによって得られたＪ上の対応点である。しかし、後述する高精度化手法を用いて対応点を得た場合、基準点は格子状とならない。そのような場合にはマルチレベルＢ-スプライン［非特許文献１２参照］を利用してＦＦＤを計算する必要がある。マルチレベルＢ-スプラインは、不規則に並んだデータ点からＣ２−連続な補間を高速に計算するアルゴリズムである。本手法の詳細については非特許文献１２、１３を参考にされたい。
【０１２３】
〔レジストレーションの高精度化〕
ボリューム間の対応関係から剛体変形パラメータを推定する場合、対応関係に誤り（誤対応）が含まれていると位置合わせ精度が低下してしまう。そこで、ここでは、３次元ＰＯＣに基づく対応付けにおける誤対応を削減し、剛体変形パラメータの推定精度を向上させるための手法について述べる。
【０１２４】
〔基準点の選別〕
例えば医用ボリュームデータなどを扱う場合、断層画像に造影する部位は撮影に利用した装置のモダリティ（ＣＴやＭＲＩなど）によって異なる。レジストレーションにより両者に不足する情報を補完し合うことは、診断能の向上やより詳細な人体構造の把握につながるため、モダリティの異なるボリュームをレジストレーションする意義は大きい。しかし、モダリティが異なるボリュームデータに対して３次元ＰＯＣに基づく対応点探索を適用する場合、両者で位相情報が全く異なるような箇所に基準点を設定すると正しく対応点を求めることができず、誤対応となる。
【０１２５】
これを避けるため、２つのデータ間で位相情報が共通するような箇所に基準点を配置する。例えばＸ線ＣＴとＭＲＩで位相情報が共通となるような部位としては、骨や皮膚などが挙げられるが、剛体レジストレーションにおいてはその運動が剛体変形で記述しやすい硬組織周辺を利用するのが望ましい。そこで本手法では、ＣＴ値に基づいてＸ線ＣＴデータ上で骨に相当する箇所に基準点を配置し、対応点探索を行う。
【０１２６】
〔ロバスト推定の利用〕
前述の基準点選別に加え、ＰＯＣ関数のピーク値に対して閾値を設けることでさらに誤対応を減らすことが可能である。３次元ブロック間で計算されるＰＯＣ関数は、両ブロックが類似しているほど高いピーク値を示す（１に近付く）ため、ピーク値が低い対応点は信頼性が低いと考えることができる。そこで、本手法では、このピーク値に対して閾値を設定し、閾値を下回る対応点は誤対応として除外する。なお、経験上、閾値は０．１程度が妥当である。
【０１２７】
しかし、実際には、ピーク値に対する閾値処理だけで全ての誤対応を除去することは困難である。そこで、本手法では、対応関係に誤りが含まれていても高精度に剛体変形パラメータを求めることができるロバスト推定手法として、ＲＡＮＳＡＣ（RANdom SAmple Consensus）［非特許文献１４参照］を用いる。以下に、ＲＡＮＳＡＣを用いた剛体変形パラメータの推定手法を示す。
【０１２８】
Ｓｔｅｐ１：基準点と対応点のペアからランダムに３組選び、その３組から剛体変形パラメータの推定を行う。
Ｓｔｅｐ２：Ｓｔｅｐ１で求めたパラメータに基づいて、全ての基準点に剛体変形を施す。
Ｓｔｅｐ３：剛体変形された基準点とその対応点の距離が閾値以内ならば、これをｉｎｌｉｅｒ（外れ点でない点）に含める。
Ｓｔｅｐ４：ｉｎｌｉｅｒの個数を数える。
Ｓｔｅｐ５：Ｓｔｅｐ１〜４をＮｒ回繰り返す。
Ｓｔｅｐ６：Ｎｒ個のｉｎｌｉｅｒ集合のうち最も要素数が多いものを用いて、変形パラメータを再推定する。
【０１２９】
〔実験・考察〕
医用ボリュームデータを例として、本手法（ボリュームレジストレーション手法）の性能評価を行う。まず、剛体レジストレーションについて従来手法と性能比較を行い、本手法では従来手法に比べて高精度かつ高速なレジストレーションが可能であることを示す。続いて、ＦＦＤを用いた非剛体レジストレーションを適用し、本手法により非剛体変形に対しても高精度な位置合わせが可能であることを示す。最後に、Ｘ線ＣＴとＭＲＩから得られたボリュームに対してレジストレーションを行い、本手法がモダリティの違いにロバストであることを示す。
【０１３０】
〔剛体レジストレーション〕
ここでは、歯科用ＣＴで撮影した顎部ボリュームデータを用いる。一人の被験者から異なるタイミングで取得した４つのデータを利用し、これらをそれぞれボリューム１〜４とする。そして、全ての組み合わせ（４Ｃ２＝６組）に対して位置合わせを行った。再構成したボリュームは５１２×５１２×５１２のボクセル数で構成され、０．２×０．２×０．２ｍｍの分解能をもつ。４つのボリュームのうち１つは意図的に首を傾けて撮影されているため、比較的大きな位置ずれや歪みが加わっている。図１０（ａ）、（ｂ）に実験にて使用したボリュームの一部を示す。
【０１３１】
本実験では、従来手法としてStudholmeらの手法［非特許文献５］を利用した。これは、２つのデータ間の正規化相互情報量（Normalized Mutual Information：ＮＭＩ）が最大となるように剛体変形パラメータを最適化する手法であり、医用データのレジストレーションにおいて一般に広く用いられている。非剛体レジストレーションでは、Rueckertらの手法[非特許文献７]を従来手法として用いた。Studholmeらの手法と同じく正規化相互情報量の最大化に基づく手法であるが、変形モデルとしてFFDを利用しているため非剛体なひずみを補正することが可能である。
【０１３２】
実験にあたり、ＰＯＣに基づく対応点探索に用いたパラメータを図８に示す。Studholmeらの手法における非線形最適化にはNelder-Mead法［非特許文献１５参照］を用い、最適化計算における初期値は適当な値を手動で設定した。
【０１３３】
レジストレーションの精度を評価するために、次の相関係数（Correlation Coefficient：ＣＣ）を指標として用いる。
【０１３４】
【数１１】
【０１３５】
本手法と従来手法を用いてレジストレーションを行い、それらの性能を相関係数により比較した結果を図９に示す。ボリューム１，２，４は自然な姿勢で撮影したものであり、ボリューム３は意図的に首を傾けて撮影したものである。相関係数による評価では１に近いほどレジストレーション精度が高いといえる。
【０１３６】
図９を見ると、本手法は全てのボリュームの組み合わせに対して従来手法と同等またはそれ以上の精度を示しており、高いレジストレーション精度を持っていることがわかる。また、本実験において、ボリュームデータ一組の剛体レジストレーションに要する時間は、Studholmeらの手法が約４分であるのに対し、本手法では約１分であり、本手法によって７５％の計算時間削減が達成されている。非剛体レジストレーションの場合は、Rueckertらの手法が約６６５分、本手法が約６分であり、本手法によって約９９％の計算時間削減が達成されている。
【０１３７】
なお、両手法の実装にはMＡＴＬＡＢ7.4.0を利用し、実験に使用したパーソナルコンピュータ（ＰＣ）はCPU：Xeon(3.0GHz×2)、Memory：16GB、OS:CentOS 4.7である。このように、３次元ＰＯＣに基づくボリュームレジストレーション手法は、精度、計算時間の両面において従来手法よりも優れていることがわかる。
【０１３８】
ボリューム１と３を用いた提案手法による非剛体レジストレーション結果を図１０（ｃ），（ｄ）に示す。ボリューム３に非剛体変形を加えて補正したデータのコロナル像にボリューム１の輪郭線を重ねて表示した。ボリューム１と３の間には非剛体な変形が加わっているため、剛体レジストレーションのみでは図１０（ｃ）のように位置ずれが目立つ。このような場合には非剛体レジストレーションが必要であり、最終的には本手法により図１０（ｄ）に示す通り非剛体なひずみが補正されていることがわかる。
【０１３９】
〔マルチモーダルレジストレーション〕
Ｘ線ＣＴとＭＲＩのそれぞれで同一の被験者の頭部データを取得し、本手法を用いて剛体レジストレーションを行った。ボリュームを再構成する際、ボクセル数を５１２×５１２×５１２，ボクセル分解能を０．４６９×０．４６９×０．４６９ｍｍとしている。
【０１４０】
図１１（ａ），（ｂ）にＸ線ＣＴデータおよびＭＲＩデータのアキシャル像を示す。図１１（ｃ），（ｄ）は、本手法によるレジストレーション結果である。図１１（ｃ）は、レジストレーション後のＭＲＩ画像上にＸ線ＣＴデータの輪郭線を重ねたものであり、図１１（ｄ）ではＭＲＩ画像とＸ線ＣＴ画像を格子状に交互に合成している。これらの結果より、モダリティが異なる場合でも本手法によって高精度なレジストレーションが可能であることがわかる。
【０１４１】
〔むすび〕
本手法では、位相限定相関法を３次元に拡張し、これを用いたボリュームデータの対応付け手法を提案した。本手法では、ＰＯＣを利用した３次元ブロックマッチングと階層探索を組み合わせることでサブボクセル精度でボリュームデータ間の対応付けを行う。ボクセル同士の対応関係はボリュームデータのレジストレーションなどの際にきわめて有用であり、従来技術に比べて精度を低下させることなく大幅に高速なレジストレーションを実現することができる。
【０１４２】
本手法は、モダリティの違いに対してもロバストに対応付けを行うことができ、また、非剛体レジストレーションにも容易に拡張することができる。タイミングを変えて取得した複数のＸ線ＣＴおよびＭＲＩデータを用いて、本手法が従来手法に比べて精度、計算時間の両面で高い性能を示すことを確認した。
【産業上の利用可能性】
【０１４３】
本発明のボリュームデータ間の対応付け方法は、同一対象物のほぼ同一箇所を撮像した３次元の立体画像のデータ（ボリュームデータ）間の対応付けを行う方法として、医療などの分野に限らず、各種の分野において適用することが可能である。
【符号の説明】
【０１４４】
１…参照ボリュームデータ取込部、２…入力ボリュームデータ取込部、３…撮像条件補正部、４…剛体変形処理部、４Ａ…オブジェクト領域抽出部、４Ｂ…対応点探索部、４Ｃ…剛体変形パラメータ推定部、４Ｄ…剛体変形実行部、４Ｅ…反復処理命令部、５…非剛体変形処理部、５Ａ…オブジェクト領域抽出部、５Ｂ…対応点探索部、５Ｃ…非剛体変形実行部、６…ボリュームデータ記憶部、Ｗ１，Ｗ２…探索ウィンドウ（３次元ブロック）、１００…ボリュームデータ間位置合わせ装置。
【特許請求の範囲】
【請求項１】
同一対象物を撮像した２つの３次元の立体画像のデータを第１および第２のボリュームデータとして取り込むデータ取込ステップと、
前記第１のボリュームデータから複数の部分領域を抽出し、この抽出した複数の部分領域にそれぞれ設定された基準点に対応する前記第２のボリュームデータの対応点を、前記基準点を含む前記第１のボリュームデータの部分領域と前記基準点に対応する前記対応点を含む前記第２のボリュームデータの部分領域との相関から探索する対応点探索ステップと、
前記探索された対応点と基準点との対応関係から前記第１のボリュームデータと前記第２のボリュームデータとの間の回転ずれおよび位置ずれを表すパラメータを剛体変形パラメータとして推定する剛体変形パラメータ推定ステップと
を備えることを特徴とするボリュームデータ間の対応付け方法。
【請求項２】
請求項１に記載されたボリュームデータ間の対応付け方法において、
前記第１のボリュームデータと前記第２のボリュームデータとの間の撮像条件に起因する相違を補正する補正ステップ
を備えることを特徴とするボリュームデータ間の対応付け方法。
【請求項３】
請求項１又は２に記載されたボリュームデータ間の対応付け方法において、
前記対応点探索ステップは、
前記第１のボリュームデータにおける所定の特徴のある領域を検索する特徴領域検索ステップを備え、
前記特徴領域検索ステップにより検索された特徴のある領域を前記部分領域の抽出領域とする
ことを特徴とするボリュームデータ間の対応付け方法。
【請求項４】
請求項１〜３の何れか１項に記載されたボリュームデータ間の対応付け方法において、
前記対応点探索ステップは、
前記第１および第２のボリュームデータを段階的に縮小し、それぞれ当該第１および第２のボリュームデータを最下位階層とし最も縮小された第１および第２のボリュームデータを最上位階層とする第１群および第２群のボリュームデータを生成する第１ステップと、
前記第１群の最上位階層のボリュームデータに初期基準点を設定する第２ステップと、
前記第２群の最上位階層のボリュームデータの前記初期基準点に対応する点に初期対応点を設定する第３ステップと、
前記第１群のボリュームデータのうち前記初期基準点が設定された階層より１階層下の階層のボリュームデータに前記初期基準点に対応する基準点を定め、この基準点を中心とする部分領域を設定する第４ステップと、
前記第２群のボリュームデータのうち前記初期対応点が設定された階層より１階層下の階層のボリュームデータに前記初期対応点に対応する対応候補点を定め、この対応候補点を中心とする部分領域を設定する第５ステップと、
前記第４ステップで設定した基準点を中心とする部分領域と前記第５ステップで設定した対応候補点を中心とする部分領域との相関から、前記基準点に対応する前記第２群の前記対応候補点が定められたボリュームデータにおける対応点を決定する第６ステップと、
前記第４ステップにおいて設定された基準点を初期基準点に設定し、前記第６ステップにおいて決定された対応点を初期対応点に設定し、前記第１群の階層および前記第２群の階層を順次下げて、前記第４ステップ〜第６ステップを実行することを繰り返すことにより、最終的に前記第１のボリュームデータに定められる基準点に対応する前記第２のボリュームデータにおける正規の対応点を決定する第７ステップと、
前記第２ステップで前記最上位階層のボリュームデータに設定する初期基準点を変更し、前記第２ステップ〜第７ステップを実行することを繰り返すことによって、前記第１のボリュームデータに定められる複数の基準点に対応する前記第２のボリュームデータにおける正規の対応点を決定し、この決定された対応点を探索された対応点として記憶する第８ステップと
を備えることを特徴とするボリュームデータ間の対応付け方法。
【請求項５】
請求項１乃至４の何れか１項に記載されたボリュームデータ間の対応付け方法において、
前記対応点探索ステップで用いられる相関は位相限定相関である
ことを特徴とするボリュームデータ間の対応付け方法。
【請求項６】
請求項１乃至５の何れか１項に記載されたボリュームデータ間の対応付け方法において、
前記推定した剛体変形パラメータに基づいて前記第１のボリュームデータと前記第２のボリュームデータとの間の回転ずれおよび位置ずれを補正する剛体変形ステップ
を備えることを特徴とするボリュームデータ間の対応付け方法。
【請求項７】
請求項６に記載されたボリュームデータ間の対応付け方法において、
所定の条件を満たすまで、前記対応点探索ステップ、前記剛体変形パラメータ推定ステップおよび剛体変形ステップを繰り返す
ことを特徴とするボリュームデータ間の対応付け方法。
【請求項８】
請求項６又は７に記載されたボリュームデータ間の対応付け方法において、
前記剛体変形ステップによって回転ずれおよび位置ずれが補正された前記第１および第２のボリュームデータを剛体変形処理後の第１および第２のボリュームデータとし、前記剛体変形処理後の第１のボリュームデータから複数の部分領域を抽出し、この抽出した複数の部分領域にそれぞれ設定された基準点に対応する前記剛体変形処理後の第２のボリュームデータの対応点を、前記基準点を含む前記剛体変形処理後の第１のボリュームデータの部分領域と前記基準点に対応する前記対応点を含む前記剛体変形処理後の第２のボリュームデータの部分領域との相関から探索する剛体変形処理後対応点探索ステップと、
この剛体変形処理後対応点探索ステップによって探索された対応点の情報を用いて前記剛体変形処理後の第１および第２のボリュームデータの何れか一方の非剛体部分を変形させて重ね合わせる非剛体変形ステップと
をさらに備えることを特徴とするボリュームデータ間の対応付け方法。
【請求項１】
同一対象物を撮像した２つの３次元の立体画像のデータを第１および第２のボリュームデータとして取り込むデータ取込ステップと、
前記第１のボリュームデータから複数の部分領域を抽出し、この抽出した複数の部分領域にそれぞれ設定された基準点に対応する前記第２のボリュームデータの対応点を、前記基準点を含む前記第１のボリュームデータの部分領域と前記基準点に対応する前記対応点を含む前記第２のボリュームデータの部分領域との相関から探索する対応点探索ステップと、
前記探索された対応点と基準点との対応関係から前記第１のボリュームデータと前記第２のボリュームデータとの間の回転ずれおよび位置ずれを表すパラメータを剛体変形パラメータとして推定する剛体変形パラメータ推定ステップと
を備えることを特徴とするボリュームデータ間の対応付け方法。
【請求項２】
請求項１に記載されたボリュームデータ間の対応付け方法において、
前記第１のボリュームデータと前記第２のボリュームデータとの間の撮像条件に起因する相違を補正する補正ステップ
を備えることを特徴とするボリュームデータ間の対応付け方法。
【請求項３】
請求項１又は２に記載されたボリュームデータ間の対応付け方法において、
前記対応点探索ステップは、
前記第１のボリュームデータにおける所定の特徴のある領域を検索する特徴領域検索ステップを備え、
前記特徴領域検索ステップにより検索された特徴のある領域を前記部分領域の抽出領域とする
ことを特徴とするボリュームデータ間の対応付け方法。
【請求項４】
請求項１〜３の何れか１項に記載されたボリュームデータ間の対応付け方法において、
前記対応点探索ステップは、
前記第１および第２のボリュームデータを段階的に縮小し、それぞれ当該第１および第２のボリュームデータを最下位階層とし最も縮小された第１および第２のボリュームデータを最上位階層とする第１群および第２群のボリュームデータを生成する第１ステップと、
前記第１群の最上位階層のボリュームデータに初期基準点を設定する第２ステップと、
前記第２群の最上位階層のボリュームデータの前記初期基準点に対応する点に初期対応点を設定する第３ステップと、
前記第１群のボリュームデータのうち前記初期基準点が設定された階層より１階層下の階層のボリュームデータに前記初期基準点に対応する基準点を定め、この基準点を中心とする部分領域を設定する第４ステップと、
前記第２群のボリュームデータのうち前記初期対応点が設定された階層より１階層下の階層のボリュームデータに前記初期対応点に対応する対応候補点を定め、この対応候補点を中心とする部分領域を設定する第５ステップと、
前記第４ステップで設定した基準点を中心とする部分領域と前記第５ステップで設定した対応候補点を中心とする部分領域との相関から、前記基準点に対応する前記第２群の前記対応候補点が定められたボリュームデータにおける対応点を決定する第６ステップと、
前記第４ステップにおいて設定された基準点を初期基準点に設定し、前記第６ステップにおいて決定された対応点を初期対応点に設定し、前記第１群の階層および前記第２群の階層を順次下げて、前記第４ステップ〜第６ステップを実行することを繰り返すことにより、最終的に前記第１のボリュームデータに定められる基準点に対応する前記第２のボリュームデータにおける正規の対応点を決定する第７ステップと、
前記第２ステップで前記最上位階層のボリュームデータに設定する初期基準点を変更し、前記第２ステップ〜第７ステップを実行することを繰り返すことによって、前記第１のボリュームデータに定められる複数の基準点に対応する前記第２のボリュームデータにおける正規の対応点を決定し、この決定された対応点を探索された対応点として記憶する第８ステップと
を備えることを特徴とするボリュームデータ間の対応付け方法。
【請求項５】
請求項１乃至４の何れか１項に記載されたボリュームデータ間の対応付け方法において、
前記対応点探索ステップで用いられる相関は位相限定相関である
ことを特徴とするボリュームデータ間の対応付け方法。
【請求項６】
請求項１乃至５の何れか１項に記載されたボリュームデータ間の対応付け方法において、
前記推定した剛体変形パラメータに基づいて前記第１のボリュームデータと前記第２のボリュームデータとの間の回転ずれおよび位置ずれを補正する剛体変形ステップ
を備えることを特徴とするボリュームデータ間の対応付け方法。
【請求項７】
請求項６に記載されたボリュームデータ間の対応付け方法において、
所定の条件を満たすまで、前記対応点探索ステップ、前記剛体変形パラメータ推定ステップおよび剛体変形ステップを繰り返す
ことを特徴とするボリュームデータ間の対応付け方法。
【請求項８】
請求項６又は７に記載されたボリュームデータ間の対応付け方法において、
前記剛体変形ステップによって回転ずれおよび位置ずれが補正された前記第１および第２のボリュームデータを剛体変形処理後の第１および第２のボリュームデータとし、前記剛体変形処理後の第１のボリュームデータから複数の部分領域を抽出し、この抽出した複数の部分領域にそれぞれ設定された基準点に対応する前記剛体変形処理後の第２のボリュームデータの対応点を、前記基準点を含む前記剛体変形処理後の第１のボリュームデータの部分領域と前記基準点に対応する前記対応点を含む前記剛体変形処理後の第２のボリュームデータの部分領域との相関から探索する剛体変形処理後対応点探索ステップと、
この剛体変形処理後対応点探索ステップによって探索された対応点の情報を用いて前記剛体変形処理後の第１および第２のボリュームデータの何れか一方の非剛体部分を変形させて重ね合わせる非剛体変形ステップと
をさらに備えることを特徴とするボリュームデータ間の対応付け方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【公開番号】特開２０１１−４１６５６（Ｐ２０１１−４１６５６Ａ）
【公開日】平成２３年３月３日（２０１１．３．３）
【国際特許分類】
【出願番号】特願２００９−１９１３１５（Ｐ２００９−１９１３１５）
【出願日】平成２１年８月２０日（２００９．８．２０）
【出願人】（０００００６６６６）株式会社山武 (1,808)
【出願人】（５０４１５７０２４）国立大学法人東北大学 (2,297)
【Ｆターム（参考）】
【公開日】平成２３年３月３日（２０１１．３．３）
【国際特許分類】
【出願日】平成２１年８月２０日（２００９．８．２０）
【出願人】（０００００６６６６）株式会社山武 (1,808)
【出願人】（５０４１５７０２４）国立大学法人東北大学 (2,297)
【Ｆターム（参考）】
[ Back to top ]