３次元ずれ量計測方法
【課題】少ない演算量で、高速かつ正確に、レジストレーションを行う。
【解決手段】基準立体画像および比較対象立体画像の原画像のデータサイズを縮小し、第１ステージの基準立体画像および比較対象立体画像とする。これら立体画像から概略的な回転軸および回転角度として比較確定回転軸ＪＡおよび比較確定回転角度θＡを求める。次に、第２ステージへ進み、第１ステージでの比較確定回転軸ＪＡおよび比較確定回転角度θＡから、より正確な比較確定回転軸ＪＢおよび比較確定回転角度θＢを求める。そして、第３ステージへ進み、第１ステージよりもデータサイズを大きくした基準立体画像および比較対象立体画像を使用し、第２ステージでの比較確定回転軸ＪＢおよび比較確定回転角度θＢから、さらに正確な比較確定回転軸ＪＣおよび比較確定回転角度θＣを求める。このようにして、最終ステージまで進み、正確な比較確定回転軸ＪENDおよび比較確定回転角度θENDを求める。
【解決手段】基準立体画像および比較対象立体画像の原画像のデータサイズを縮小し、第１ステージの基準立体画像および比較対象立体画像とする。これら立体画像から概略的な回転軸および回転角度として比較確定回転軸ＪＡおよび比較確定回転角度θＡを求める。次に、第２ステージへ進み、第１ステージでの比較確定回転軸ＪＡおよび比較確定回転角度θＡから、より正確な比較確定回転軸ＪＢおよび比較確定回転角度θＢを求める。そして、第３ステージへ進み、第１ステージよりもデータサイズを大きくした基準立体画像および比較対象立体画像を使用し、第２ステージでの比較確定回転軸ＪＢおよび比較確定回転角度θＢから、さらに正確な比較確定回転軸ＪＣおよび比較確定回転角度θＣを求める。このようにして、最終ステージまで進み、正確な比較確定回転軸ＪENDおよび比較確定回転角度θENDを求める。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、基準となる３次元の立体画像を基準立体画像、この基準立体画像に対して比較対象とする３次元の立体画像を比較対象立体画像とし、基準立体画像と比較対象立体画像との間のずれ量を計測する３次元ずれ量計測方法に関するものである。
【背景技術】
【０００２】
従来より、医学分野などにおいては、多様な画像化方式によって立体画像を得ている。例えば、磁気共鳴断層撮影法（ＭＲＩ）、Ｘ線コンピュータ断層撮影法（ＣＴ）、陽電子放射断層撮影法（ＰＥＴ）、単光子射出断層撮影法（ＳＰＥＣＴ）、超音波断層撮影法などによって、３次元の立体画像を得ている。
【０００３】
このような方法で得られた３次元の立体画像は、ボクセル（voxel）値の配列として、通常、コンピュータのメモリなどにデジタル的に保存される。各ボクセルは、互いに直交するＸ，Ｙ，Ｚの３軸で表される３次元空間の座標位置（ｘ，ｙ，ｚ）と結びついており、ボクセル値として色値（典型的にはグレイスケール）が割り当てられる。
【０００４】
図１８にＸ線コンピュータ断層撮影法による３次元の立体画像の取得例を示す。この例では、スキャナ装置としてコーンビームＣＴを用いている。コーンビームＣＴは、中間点Ｃ０を挾んで対向してＸＹ平面内に設けられたＸ線源１と２次元撮像装置２とを備え、このＸ線源１と２次元撮像装置２とを中間点Ｃ０を通るＺ軸を中心として水平方向に定速度で回転させながら、被写体（人体）３の撮影を行う。この撮影に際して、被写体３は、その頭部の中心を中間点Ｃ０に位置させる。
【０００５】
被写体３の撮影は、Ｘ線源１と２次元撮像装置２とが１回転する間、所定の時間間隔で行われる。この結果、１回転する間に任意の異なる角度（以下、この角度を投影角度と呼ぶ）から撮影されたｎ個（例えば、２５６個）の２次元の投影画像を取得して、このｎ個の投影画像を再構成処理することで３次元の立体画像を得ることができる。この３次元の立体画像は、互いに直交するＸ，Ｙ，Ｚの３軸で表される３次元空間の座標位置（ｘ，ｙ，ｚ）と結びつけたボクセルデータ群として表現される。
【０００６】
この３次元の立体画像は、医療や歯科などの分野において、各種の診断を行うために使用される。しかし、情報量が多いために見落としが発生する虞があり、正確な診断を行うために過大な時間を必要とする問題がある。このため、コンピュータ上で、画像処理を行うことで、この問題を解決しようと試みられている。
【０００７】
この場合の画像処理の１つに、レジストレーションと呼ばれる位置合わせを行った後、２つの立体画像を比較するという方法がある。例えば、前回撮影した被写体３（被写体Ａ）の３次元の立体画像を基準立体画像とし、今回撮影した被写体３（被写体Ｂ）の３次元の立体画像を比較対象立体画像とし、基準立体画像と比較対象立体画像との位置合わせを行い、この位置合わせが行われた２つの立体画像の差分画像をディスプレイ上に表示する。これにより、前回撮影時と今撮影時との間の変化部分のみを抽出し、短時間でかつ正確な診断を行うことができるようになる。
【０００８】
この方法において、基準立体画像と比較対象立体画像とのレジストレーションは、６つのパラメータを求めることで実現される。６つのパラメータは、Ｘ，Ｙ，Ｚの各軸の平行ずれ量Δｘ，Δｙ，Δｚと、Ｘ，Ｙ，Ｚの各軸の回転ずれ量Δθ，Δφ，ΔΨである。平行ずれ量Δｘ，Δｙ，Δｚについては、回転ずれ量Δθ，Δφ，ΔΨを求めた後に、３次元相関処理を行うことで、容易に求めることができる。回転ずれ量Δθ，Δφ，ΔΨは、Ｘ，Ｙ，Ｚの各軸の回転角度を３つの自由度とし、この３つの自由度を変化させながら３次元相関処理を行うことで求めることが可能である（例えば、特許文献１参照）。
【０００９】
【特許文献１】特開平９−２２４０６号公報
【特許文献２】特開２００１−２１２１２６号公報
【発明の開示】
【発明が解決しようとする課題】
【００１０】
しかしながら、上述した基準立体画像と比較対象立体画像とのレジストレーションでは、回転ずれ量Δθ，Δφ，ΔΨを求める際、Ｘ，Ｙ，Ｚの各軸の回転角度を自由度として変化させるため、３乗の演算回数が必要となる。例えば、３つの自由度それぞれに１０種類の回転角度を設定した場合、１０×１０×１０＝１０００回の３次元の相関処理が必要となる。このため、膨大な演算を必要とし、高速処理が困難となる。
【００１１】
なお、相関処理を用いないレジストレーションとして、立体画像に現れている被写体の構造に詳しい専門家が、位置合わせする立体画像のそれぞれに目印をつけ、一方の立体画像を移動回転させることによって位置合わせする方法がある（例えば、特許文献２参照）。
【００１２】
しかしながら、この方法では、位置合わせ精度が選択された目印の数と位置によって制限されてしまう。また、目印の数が少な過ぎると、位置合わせが不正確になる。また、目印が多くても必ずしも位置合わせが正確となるわけではなく、位置合わせに要する計算を複雑にするだけである。また、オペレータが手作業で対応点を選択する必要があり、時間がかかるとともに熟練が必要となる。さらに、すべての立体画像で適当な構造上の目印を見出せるわけでもない。
【００１３】
本発明は、このような課題を解決するためになされたもので、その目的とするところは、少ない演算量で、高速かつ正確に、レジストレーションを行うことが可能な３次元ずれ量計測方法を提供することにある。
【課題を解決するための手段】
【００１４】
このような目的を達成するために本発明は、基準となる３次元の立体画像を基準立体画像、この基準立体画像に対して比較対象とする３次元の立体画像を比較対象立体画像とし、基準立体画像と比較対象立体画像との間のずれ量を計測する３次元ずれ量計測方法において、比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ（Ｋ≧２）個の回転軸を比較候補回転軸として設定し（第１ステップ）、Ｋ個の各比較候補回転軸を中心として比較対象立体画像を所定角度範囲内でＬ（Ｌ≧２）通り回転させたＫ×Ｌ個の画像を比較候補立体画像とし、このＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値を各個に算出し（第２ステップ）、この算出されたＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして求め（第３ステップ）、算出されたＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から比較確定回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出し（第４ステップ）、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル（原点からの単位長さ）に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像の第１ステージでの比較確定回転軸ＪＡとする（第５ステップ）、ようにしたものである。
【００１５】
本発明を適用しようとする３次元空間において、その３次元空間における回転角度のずれには、様々な表現方法が考えられる。最も簡単に表現できるのは、互いに直交するＸ，Ｙ，Ｚの３軸のそれぞれを回転軸とした３つの回転角度により表現することである。すなわち、回転軸を３つとし、これらの回転軸を中心としてそれぞれある角度にて回転したと考え、合計３回の回転が発生していると考えることである。本発明では、この考えをさらに進め、原点を中心として複数回の回転を行うことは、ある回転軸を中心として１回の回転を行ったに過ぎないと考え、回転軸と回転角度の２つのパラメータを求めることにする。
【００１６】
この考えに基づき、本発明では、比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ（Ｋ≧２）個の回転軸を比較候補回転軸として設定する（第１ステップ）。例えば、原点を中心とする正二十面体を考え、この正二十面体の中心と各頂点とを結ぶ１２個の軸を比較候補回転軸として設定する。そして、このＫ個の各比較候補回転軸を中心として比較対象立体画像を所定角度範囲内でＬ（Ｌ≧２）通り回転させたＫ×Ｌ個の画像（例えば、２〜３２゜の角度範囲内で２゜刻みで回転させた１２×１６＝１９２個の画像）を比較候補立体画像とし、このＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値を各個に算出する（第２ステップ）。
【００１７】
そして、この算出したＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして求める（第３ステップ）。また、この算出したＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から比較確定回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する（第４ステップ）。そして、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル（原点からの単位長さ）に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像の第１ステージでの比較確定回転軸ＪＡとする（第５ステップ）。
【００１８】
このようにして、本発明では、第１ステージでの比較確定回転軸ＪＡと比較確定回転角度θＡを求める。この場合、第１ステージで求められた比較確定回転軸ＪＡと比較確定回転角度θＡを用いて基準立体画像と比較対象立体画像との位置合わせ（レジストレーション）を行うようにしてもよいが、レジストレーションの精度をさらに高めるためには、さらに第２ステージに進むことが望ましい。
【００１９】
第２ステージでは、第１ステージでの比較確定回転軸ＪＡを中心とし、第１ステージでの比較確定回転角度θＡを含む所定角度範囲内で比較対象立体画像をＭ（Ｍ≧２）通り回転させたＭ個の画像（例えば、比較確定回転角度θＡを中心とする±１０゜の角度範囲内で０．５゜刻みで回転させた４１個の画像）を比較候補立体画像とし、このＭ個の比較候補立体画像と基準立体画像との相関値を各個に算出する（第６ステップ）。そして、この算出したＭ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとして求める（第７ステップ）。
【００２０】
そして、第１ステージでの比較確定回転軸ＪＡを中心に任意のＮ個（例えば、１１×１１＝１２１個）の比較候補回転軸を設定し、このＮ個の比較候補回転軸のそれぞれについてその回転軸を中心として比較確定回転角度θＢ回転させたＮ個の比較候補立体画像と基準立体画像との相関値を各個に算出する（第８ステップ）。そして、このＮ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとする（第９ステップ）。
【００２１】
なお、上述した第７ステップにおいて、Ｍ個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとして求めるようにしてもよい。また、上述した第９ステップにおいて、Ｎ個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとするようにしてもよい。
【００２２】
また、本発明において、第２ステージに進む場合、第１ステップに入る前に、基準立体画像および比較対象立体画像の原画像のデータサイズを小さくして第１ステージで用いる基準立体画像および比較対象立体画像とし、第１ステージで用いる基準立体画像および比較対象立体画像のデータサイズよりもそのデータサイズを大きくしながら、第６〜第９ステップの処理を予め定められた最終ステージまで繰り返すようにしてもよい。
【００２３】
例えば、基準立体画像および比較対象立体画像の原画像のデータサイズ（画素数）を共に５１２×５１２×５１２とした場合、これら立体画像のデータサイズを縮小して１６×１６×１６の立体画像とする。第１ステージでは、この１６×１６×１６の基準立体画像および比較対象立体画像を用いて第１〜第５ステップの処理を行う。第２ステージでも、第１ステージと同様にして、１６×１６×１６の基準立体画像および比較対象立体画像を用いて第６〜第９ステップの処理を行う。そして、第３ステージから、基準立体画像および比較対象立体画像のデータサイズを大きくし（例えば、６４×６４×６４とする）、第６〜第９ステップの処理を行う。このようにして、最終ステージまで、基準立体画像および比較対象立体画像のデータサイズを大きくしながら、第６〜第９ステップの処理を繰り返す。なお、第２ステージから基準立体画像および比較対象立体画像のデータサイズを大きくするようにしてもよく、第３ステージを最終ステージとしてもよい。
【００２４】
このようにして、本発明では、基準立体画像および比較対象立体画像のデータサイズを最初は小さくし、ステージが進むにつれて大きくして行くことにより、最初は少ないボクセルデータ群で概略的な比較確定回転軸と比較確定回転角度を素速く求め、徐々にボクセルデータ群の数を増やして、より正確な比較確定回転軸と比較確定回転角度を求めて行くようにして、高速度でかつ高精度に、レジストレーションを行うことが可能となる。
【００２５】
本発明では、比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ個の回転軸を比較候補回転軸として設定し、このＫ個の比較候補回転軸を中心として回転させた比較対象立体画像と基準立体画像との相関値を求めるが、このような処理において比較対象立体画像と基準立体画像の回転中心とされる原点は同じでなければならない。しかし、３次元の立体画像において、実際の回転中心は不明であるため、原点を決めることが難しい。
【００２６】
本発明において、正確な原点を決めるために、第１ステップの処理に入る前に、基準立体画像の原画像および比較対象立体画像の原画像を構成するボクセルデータ群に対してそれぞれ３次元フーリエ変換を施し、この３次元フーリエ変換が施された基準立体画像の原画像および比較対象立体画像の原画像を構成するボクセルデータ群の３次元周波数空間（ボクセルデータ群の中心を直流とする周波数空間）の中心を基準立体画像および比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心とするようにしてもよい。
【００２７】
なお、ボクセルデータ群に対する３次元フーリエ変換の処理は、必ずしも必要というわけではない。例えば、被写体の回転角度が比較的小さいことが分かっている場合、回転中心がほぼ中心付近にあることが分かっている場合などは、３次元フーリエ変換の処理を行う必要はなく、基準立体画像および比較対象立体画像の原画像を構成するボクセルデータ群をそのまま使って、相関処理を行うことが可能である。本発明において、第１ステップの処理に入る前に基準立体画像の原画像および比較対象立体画像の原画像のデータサイズを小さくする方式とした場合、多少の回転中心のずれは無視することができる。
【発明の効果】
【００２８】
本発明によれば、比較対象立体画像に設定される原点を中心として任意の方向へＫ個の比較候補回転軸を設定し、このＫ個の各比較候補回転軸を中心として比較対象立体画像を所定角度範囲内でＬ通り回転させたＫ×Ｌ個の画像を比較候補立体画像とし、このＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値を各個に算出し、算出したＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして求め、算出したＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から比較確定回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出し、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル（原点からの単位長さ）に当該相関値に基づく重み付けを施し、この重み付けが施された単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像の第１ステージでの比較確定回転軸ＪＡとするようにしたので、この第１ステージで求められた比較確定回転軸ＪＡおよび比較確定回転角度θＡを用いて基準立体画像と比較対象立体画像との位置合わせを行うことにより、少ない演算量で、高速かつ正確に、レジストレーションを行うことが可能となる。
【００２９】
また、本発明によれば、第１ステージでの比較確定回転軸ＪＡを中心とし、第１ステージでの比較確定回転角度θＡを含む所定角度範囲内で比較対象立体画像をＭ通り回転させたＭ個の画像を比較候補立体画像とし、このＭ個の比較候補立体画像と基準立体画像との相関値を算出し、算出したＭ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとして求め、第１ステージでの比較確定回転軸ＪＡを中心に任意のＮ（Ｎ≧２）個の比較候補回転軸を設定し、このＮ個の比較候補回転軸のそれぞれについてその回転軸を中心として比較確定回転角度θＢ回転させたＮ個の比較候補立体画像と基準立体画像との相関値を算出し、Ｎ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとするようにしたので、この第２ステージで求められた比較確定回転軸ＪＢおよび比較確定回転角度θＢを用いて基準立体画像と比較対象立体画像との位置合わせを行うことにより、さらに正確に、レジストレーションを行うことが可能となる。
【００３０】
また、本発明によれば、第１ステップの処理に入る前に、基準立体画像および比較対象立体画像の原画像のデータサイズを小さくして第１ステージで用いる基準立体画像および比較対象立体画像とし、第１ステージで用いる基準立体画像および比較対象立体画像のデータサイズよりもそのデータサイズを大きくしながら第６〜第９ステップの処理を予め定められた最終ステージまで繰り返すことにより、最初は少ないボクセルデータ群で概略的な比較確定回転軸と比較確定回転角度を素速く求め、徐々にボクセルデータ群の数を増やして、より正確な比較確定回転軸と比較確定回転角度を求めて行くようにして、高速度でかつ高精度に、レジストレーションを行うことが可能となる。
【発明を実施するための最良の形態】
【００３１】
以下、本発明を図面に基づいて詳細に説明する。図１は本発明に係る３次元ずれ量計測方法の実施に用いるＸ線コンピュータ断層撮影システムの一例を示す図である。同図において、図１８と同一符号は図１８を参照して説明した構成要素と同一或いは同等構成要素を示し、その説明は省略する。
【００３２】
この実施の形態において、Ｘ線源１および２次元撮像装置２は、３次元ずれ量計測装置４に接続されている。３次元ずれ量計測装置４は、ＣＰＵ４Ａと、ＲＡＭ４Ｂと、ＲＯＭ４Ｃと、ハードディスクなどの記憶装置４Ｄと、インタフェース４Ｅ〜Ｈと、ディスプレイ４Ｉと、キーボード４Ｊと、マウス４Ｋとを備えている。
【００３３】
ＣＰＵ４Ａは、ＲＡＭ４Ｂにアクセスしながら、ＲＯＭ４Ｃや記憶装置４Ｄに格納されたプログラムに従って動作する。記憶装置４Ｄには、本実施の形態特有のプログラムとして３次元ずれ量計測プログラムが格納されている。この３次元ずれ量計測プログラムは、例えばＣＤ−ＲＯＭなどの記録媒体に記録された状態で提供され、この記録媒体から読み出されて記憶装置４Ｄにインストールされている。
【００３４】
以下、図２〜図７に示すフローチャートを参照して、記憶装置４Ｄに格納された３次元ずれ量計測プログラムに従ってＣＰＵ４Ａが実行する本実施の形態特有の処理動作について説明する。
【００３５】
なお、この３次元ずれ量計測プログラムに従う処理動作には、比較対象立体画像の基準立体画像に対する回転軸と回転角度の２つのパラメータ（ＪEND，ΔθEND）を求めるずれ量計測処理と、この求めた２つのパラメータを用いて基準立体画像と比較対象立体画像との間の位置合わせを行う位置合わせ処理と、この位置合わせされた基準立体画像と比較対象立体画像との照合を行う照合処理とが含まれる。
【００３６】
〔基準立体画像の取得〕
先ず、基準立体画像を取得するために、被写体Ａの撮影を行う（図２：ステップＳ１０１）。この被写体Ａの撮影は、Ｘ線源１と２次元撮像装置２とが１回転する間、所定の時間間隔で行う。この結果、１回転する間の任意の異なる角度（投影角度）から撮影されたｎ個（例えば、２５６個）の２次元の投影画像が得られる（ステップＳ１０２）。なお、この撮影において、被写体Ａは、その頭部の中心を中間点ＣＯ付近にきちんと位置させるものとする。
【００３７】
ＣＰＵ４Ａは、この得られたｎ個の２次元の投影画像を記憶装置４Ｄに保存した後（ステップＳ１０３）、このｎ個の投影画像を再構成処理することで３次元の立体画像（ボクセルデータ群）を得る（ステップＳ１０４）。そして、この再構成処理によって得られた３次元の立体画像を基準立体画像の原画像として記憶装置４Ｄに保存する（ステップＳ１０５）。なお、この場合、記憶装置４Ｄに保存される基準立体画像の原画像のデータサイズ（画素数）は、５１２×５１２×５１２とする。
【００３８】
〔比較対象立体画像の取得〕
次に、比較対象立体画像を取得するために、被写体Ｂの撮影を行う（図３：ステップＳ２０１）。この被写体Ｂの撮影も、Ｘ線源１と２次元撮像装置２とが１回転する間、基準立体画像の取得時と同じ所定の時間間隔で行う。この結果、１回転する間の任意の異なる角度（投影角度）から撮影されたｎ個（例えば、２５６個）の２次元の投影画像が得られる（ステップＳ２０２）。なお、この撮影において、被写体Ｂは、その頭部の中心を中間点ＣＯ付近にきちんと位置させるものとする。
【００３９】
ＣＰＵ４Ａは、この得られたｎ個の２次元の投影画像を記憶装置４Ｄに保存した後（ステップＳ２０３）、このｎ個の投影画像を再構成処理することで３次元の立体画像（ボクセルデータ群）を得る（ステップＳ２０４）。そして、この再構成処理によって得られた３次元の立体画像を比較対象立体画像の原画像として記憶装置４Ｄに保存する（ステップＳ２０５）。なお、この場合、記憶装置４Ｄに保存される比較対象立体画像の原画像のデータサイズ（画素数）は、５１２×５１２×５１２とする。
【００４０】
〔基準立体画像と比較対象立体画像との間のずれ量の計測〕
〔第１ステージ〕
次に、ＣＰＵ４Ａは、記憶装置４Ｄに保存されている５１２×５１２×５１２の基準立体画像（原画像）のデータサイズを縮小し、第１ステージで用いる基準立体画像とする（図４：ステップＳ３０１）。また、記憶装置４Ｄに保存されている５１２×５１２×５１２の比較対象立体画像（原画像）のデータサイズを縮小し、第１ステージで用いる比較対象立体画像とする（ステップＳ３０２）。
【００４１】
この例では、Ｘ，Ｙ，Ｚの３軸方向の画素データを所定画素ずつ足し合わせて行くことによって、第１ステージの基準立体画像および比較対象立体画像として１６×１６×１６のデータサイズの立体画像を得る。なお、画素データを間引いて行くことによって、データサイズを縮小するようにしてもよい。
【００４２】
そして、ＣＰＵ４Ａは、ステップ３０１で得た第１ステージの比較対象立体画像のＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ個の回転軸を比較候補回転軸として設定する（ステップＳ３０３）。この例では、原点を中心とする正二十面体（図８参照）を考え、この正二十面体の中心Ｏと各頂点とを結ぶ１２個の軸Ｊ１〜Ｊ１２を比較候補回転軸として設定する。
【００４３】
そして、このＫ個の比較候補回転軸を中心として第１ステージの比較対象立体画像を所定の角度範囲内でＬ通り回転させたＫ×Ｌ個の画像を比較候補立体画像とし（ステップＳ３０４）、このＫ×Ｌ個の比較候補立体画像と第１ステージの基準立体画像との相関値を各個に算出する（ステップＳ３０５）。
【００４４】
この例では、１２個の比較候補回転軸Ｊ１〜Ｊ１２について、その比較候補回転軸を中心として２〜３２゜の角度範囲内で２゜刻みで回転させた１２×１６＝１９２個の画像を比較候補立体画像とし、この１９２個の比較候補立体画像と第１ステージの基準立体画像との相関値ｍ１〜ｍ１９２（図９参照）を算出する。
【００４５】
そして、この算出した相関値ｍ１〜ｍ１９２の中から最大相関値を抽出し（ステップＳ３０６）、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして求める（ステップＳ３０７）。例えば、相関値ｍ１７が最大相関値であった場合、比較候補回転軸Ｊ２を中心とする回転角度θ１を第１ステージでの比較確定回転角度θＡとして求める。
【００４６】
また、ＣＰＵ４Ａは、ステップＳ３０５で算出した相関値ｍ１〜ｍ１９２の中から、第１ステージでの比較確定回転角度として求めた回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する（ステップＳ３０８）。すなわち、比較候補回転軸Ｊ１〜Ｊ１２について、回転角度θＡとした場合の比較候補立体画像の相関値を抽出する。
【００４７】
そして、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸Ｊ１〜Ｊ１２の単位ベクトル（原点からの単位長さ）に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルＶＡ（図１０参照）を求め、この合成ベクトルＶＡが示す方向を比較対象立体画像の第１ステージでの比較確定回転軸ＪＡとする（ステップＳ３０９）。
【００４８】
〔第２ステージ〕
次に、ＣＰＵ４Ａは、記憶装置４Ｄに保存されている５１２×５１２×５１２の基準立体画像（原画像）のデータサイズを縮小し、第２ステージで用いる基準立体画像とする（図５：ステップＳ３１０）。また、記憶装置４Ｄに保存されている５１２×５１２×５１２の比較対象立体画像（原画像）のデータサイズを縮小し、第２ステージで用いる比較対象立体画像とする（ステップＳ３１１）。
【００４９】
この例では、ステップ３１０，３１１は実行せずに、第１ステージで用いた１６×１６×１６の基準立体画像および比較対象立体画像を第２ステージの基準立体画像および比較対象立体画像として用いることとする。なお、ステップ３１０，３１１を実行する場合には、第１ステージで用いた基準立体画像および比較対象立体画像のデータサイズよりも大きいデータサイズに縮小し、第２ステージの基準立体画像および比較対象立体画像とする。
【００５０】
次に、ＣＰＵ４Ａは、第１ステージでの比較確定回転軸ＪＡを中心とし、第１ステージでの比較確定回転角度θＡを含む所定の角度範囲内で第２ステージの比較対象立体画像をＭ１通り回転させたＭ１個の画像を比較候補立体画像とし、このＭ１個の比較候補立体画像と第２ステージの基準立体画像との相関値を各個に算出する（ステップＳ３１２）。この例では、比較確定回転角度θＡを中心とする±１０゜の角度範囲内で０．５゜刻みで回転させた４１個の画像を比較候補立体画像とし（図１１参照）、この４１個の比較候補立体画像と第２ステージの基準立体画像との相関値を各個に算出する。
【００５１】
そして、この算出したＭ１個の比較候補立体画像と第２ステージの基準立体画像との相関値の中から最大相関値を抽出し（ステップＳ３１３）、この最大相関値が得られた比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとする（ステップＳ３１４）。
【００５２】
なお、この例では、Ｍ１個の比較候補立体画像と第２ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとしたが、Ｍ１個の比較候補立体画像と基準立体画像との相関値から関数フィッティング（例えば、シンク関数）により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとするようにしてもよい。
【００５３】
次に、ＣＰＵ４Ａは、第１ステージでの比較確定回転軸ＪＡを中心に任意のＮ１個の比較候補回転軸を設定する（ステップＳ３１５）。例えば、図１２に示すように、第１ステージでの比較確定回転軸ＪＡを中心軸とし、この比較確定回転軸ＪＡの周囲に縦方向１１個、横方向１１個、合計１２１個の原点Ｏからの比較候補回転軸を３゜刻みで設定する。そして、このＮ１個の比較候補回転軸のそれぞれについて、その回転軸を中心として比較確定回転角度θＢ回転させたＮ１個の比較候補立体画像と第２ステージの基準立体画像との相関値を各個に算出する（ステップＳ３１６）。
【００５４】
そして、このＮ１個の比較候補立体画像と第２ステージの基準立体画像との相関値の中から最大相関値を抽出し（ステップＳ３１７）、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとする（ステップＳ３１８）。
【００５５】
なお、この例では、Ｎ１個の比較候補立体画像と第２ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとするようにしたが、Ｎ１個の比較候補立体画像と第２ステージの基準立体画像との相関値から関数フィッティング（例えば、シンク関数）により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとするようにしてもよい。
【００５６】
〔第３ステージ（最終ステージ）〕
次に、ＣＰＵ４Ａは、記憶装置４Ｄに保存されている５１２×５１２×５１２の基準立体画像（原画像）のデータサイズを縮小し、第３ステージで用いる基準立体画像とする（図６：ステップＳ３１９）。また、記憶装置４Ｄに保存されている５１２×５１２×５１２の比較対象立体画像（原画像）のデータサイズを縮小し、第３ステージで用いる比較対象立体画像とする（ステップＳ３２０）。この例では、第３ステージで用いる基準立体画像および比較対象立体画像のデータサイズを６４×６４×６４とし、第１ステージで用いる基準立体画像および比較対象立体画像のデータサイズよりも大きいサイズとする。
【００５７】
次に、ＣＰＵ４Ａは、第２ステージでの比較確定回転軸ＪＢを中心とし、第２ステージでの比較確定回転角度θＢを含む所定の角度範囲内で第３ステージの比較対象立体画像をＭ２通り回転させたＭ２個の画像を比較候補立体画像とし、このＭ２個の比較候補立体画像と第３ステージの基準立体画像との相関値を各個に算出する（ステップＳ３２１）。この例では、比較確定回転角度θＢを中心とする±１゜の角度範囲内で０．１゜刻みで回転させた２１個の画像を比較候補立体画像とし（図１３参照）、この２１個の比較候補立体画像と第３ステージの基準立体画像との相関値を各個に算出する。
【００５８】
そして、この算出したＭ２個の比較候補立体画像と第３ステージの基準立体画像との相関値の中から最大相関値を抽出し（ステップＳ３２２）、この最大相関値が得られた比較候補立体画像の第２ステージでの比較確定回転軸ＪＢを中心とする回転角度を第３ステージでの比較確定回転角度θＣとする（ステップＳ３２３）。
【００５９】
なお、この例では、Ｍ２個の比較候補立体画像と第３ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の第２ステージでの比較確定回転軸ＪＢを中心とする回転角度を第３ステージでの比較確定回転角度θＣとしたが、Ｍ２個の比較候補立体画像と第３ステージの基準立体画像との相関値から関数フィッティング（例えば、シンク関数）により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の第２ステージでの比較確定回転軸ＪＢを中心とする回転角度を第３ステージでの比較確定回転角度θＣとするようにしてもよい。
【００６０】
次に、ＣＰＵ４Ａは、第２ステージでの比較確定回転軸ＪＢを中心に任意のＮ２個の比較候補回転軸を設定する（ステップＳ３２４）。例えば、図１４に示すように、第２ステージでの比較確定回転軸ＪＢを中心軸とし、この中心軸ＪＢの周囲に縦方向７個、横方向７個、合計４９個の原点Ｏからの比較候補回転軸を３゜刻みで設定する。そして、このＮ２個の比較候補回転軸のそれぞれについて、その回転軸を中心として回転角度θＣ回転させたＮ２個の比較候補立体画像と第３ステージの基準立体画像との相関値を各個に算出する（ステップＳ３２５）。
【００６１】
そして、このＮ２個の比較候補立体画像と第３ステージの基準立体画像との相関値の中から最大相関値を抽出し（ステップＳ３２６）、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第３ステージでの比較確定回転軸ＪＣとする（ステップＳ３２７）。
【００６２】
なお、この例では、Ｎ２個の比較候補立体画像と第３ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の比較候補回転軸を第３ステージでの比較確定回転軸ＪＣとするようにしたが、Ｎ２個の比較候補立体画像と第３ステージの基準立体画像との相関値から関数フィッティング（例えば、シンク関数）により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第３ステージでの比較確定回転軸ＪＣとするようにしてもよい。
【００６３】
本実施の形態では、第３ステージを最終ステージとし、この第３ステージでの比較確定回転軸ＪＣおよび比較確定回転角度θＣを最終ステージでの比較確定回転軸ＪENDおよび比較確定回転角度θENDとする。なお、第３ステージを最終ステージとせず、さらに最終ステージまでのステージ数を増やしてもよい。この場合、各ステージでは、第３ステージと同様にして、そのステージで用いる基準立体画像および比較対象立体画像のデータサイズを大きくしながら、第２ステージと同様の処理を繰り返すものとする。
【００６４】
〔基準立体画像と比較対象立体画像との間の位置合わせ〕
ＣＰＵ４Ａは、このようにして最終ステージでの比較確定回転軸ＪENDおよび比較確定回転角度θENDを求めた後、この求めた最終ステージでの比較確定回転軸ＪENDおよび比較確定回転角度θENDに基づいて、記憶装置４に保存されている基準立体画像（原画像）と比較対象立体画像（原画像）との間の回転ずれを補正する（図７：ステップＳ４０１）。すなわち、比較対象立体画像（原画像）を比較確定回転軸ＪENDを中心として比較確定回転角度θEND回転補正する。これにより、一度の処理で、基準立体画像と比較対象立体画像との間の回転ずれが補正される。
【００６５】
次に、この回転ずれが補正された基準立体画像と比較対象立体画像との間で相関演算を行って、この基準立体画像と比較対象立体画像との間のＸ，Ｙ，Ｚの３軸方向のずれ量Δｘ，Δｙ，Δｚを求め、この求めた３軸方向のずれ量Δｘ，Δｙ，Δｚに基づいて回転ずれが補正された基準立体画像と比較対象立体画像との間の平行ずれを補正する（ステップＳ４０２）。すなわち、Ｘ軸方向へΔｘ、Ｙ軸方向へΔｙ、Ｚ軸方向へΔｚ平行移動し、基準立体画像と比較対象立体画像とのＸ，Ｙ，Ｚの３軸方向の位置を合わせる。
【００６６】
これにより、基準立体画像と比較対象立体画像との間の位置合わせが完了する。なお、この際、比較対象立体画像ではなく、基準立体画像の方を回転補正したり、平行移動するようにしてもよい。
【００６７】
〔位置合わせされた基準立体画像と比較対象立体画像との照合〕
ＣＰＵ４Ａは、このようにして基準立体画像と比較対象立体画像との位置合わせを行った後、すなわち基準立体画像と比較対象立体画像とのレジストレーションを行った後、位置合わせされた基準立体画像と比較対象立体画像とを照合する（ステップＳ４０３）。
【００６８】
この例では、位置合わせされた基準立体画像（ボクセルデータ群）と比較対象立体画像（ボクセルデータ群）との相関演算を行い、この相関演算の結果得られた相関値と予め定められた照合用の閾値とを比較することによって、基準立体画像と比較対象立体画像との照合を行う。
【００６９】
なお、３次元の立体画像間の相関演算については特許文献１などに示されているので、ここでの説明は省略する。この特許文献１において、Ｎ＝３として３次元フーリエ変換、３次元逆フーリエ変換を施し、３次元の立体画像間の相関演算を行うことにより、相関成分エリア内でピークの位置ずれがＸ，Ｙ，Ｚの各軸に現れる。このＸ，Ｙ，Ｚの各軸に現れる位置ずれから基準立体画像と比較対象立体画像との相関が分かる。
【００７０】
ＣＰＵ４Ａは、ステップＳ４０２での照合の結果がＯＫであれば（ステップＳ４０４のＹＥＳ）、照合ＯＫの表示を行う（ステップＳ４０５）。照合の結果がＮＧであれば（ステップＳ４０４のＮＯ）、照合ＮＧの表示を行うとともに（ステップＳ４０６）、位置合わせされた基準立体画像と比較対象立体画像との差分画像をディスプレイ４Ｉ上に表示する（ステップＳ４０７）。
【００７１】
以上の説明から分かるように、本実施の形態によれば、基準立体画像および比較対象立体画像のデータサイズが最初は小さくされ、ステージが進むにつれて大きくされて行く。これにより、最初は少ないボクセルデータ群で概略的な比較確定回転軸と比較確定回転角度が素速く求められ、徐々にボクセルデータ群の数を増やして、より正確な比較確定回転軸と比較確定回転角度が求められて行く。このようにして、本実施の形態では、高速度でかつ高精度に、レジストレーションが行われる。この結果、照合スピードもアップし、照合結果の精度も高まる。
【００７２】
〔基準立体画像と比較対象立体画像の回転中心（原点）について〕
上述した実施の形態では、第１ステージの比較対象立体画像のＸ，Ｙ，Ｚの３軸の中心を原点としてＫ個の回転軸を比較候補回転軸として設定し、このＫ個の比較候補回転軸を中心として回転させた比較対象立体画像と第１ステージの基準立体画像との相関値を求めるが、このような処理において比較対象立体画像と基準立体画像の回転中心とされる原点は同じでなければならない。しかし、３次元の立体画像において、実際の回転中心は不明であるため、原点を決めることが難しい。
【００７３】
正確な原点を定めるために、本実施の形態において、ステップＳ３０１の処理に入る前に、記憶装置４Ｄに保存されている基準立体画像（原画像）および比較対象立体画像（原画像）を構成するボクセルデータ群に対してそれぞれ３次元フーリエ変換を施し、この３次元フーリエ変換が施された基準立体画像（原画像）および比較対象立体画像（原画像）を構成するボクセルデータ群の３次元周波数空間（ボクセルデータ群の中心を直流とする周波数空間）の中心を基準立体画像および比較対象立体画像に設定するＸ，Ｙ，Ｚの３軸の中心とすることが考えられる。
【００７４】
この場合、相関値を求める処理は、以下の手順となる。
(１)両方のボクセルデータ群に対し、３次元フーリエ変換（ＦＦＴ）、振幅化、を行う。この時、ボクセルデータ群の中心が直流になるようにする。また、３次元ＦＦＴした結果の振幅の値は、元のボクセルデータに比べて比較的大きな値をとることがある。これが演算上問題になる場合には、振幅データをログ化することで、ダイナミックレンジを抑えることができる。
【００７５】
(２)(１)で得た（必要に応じてログ化された）振幅データの片方のボクセルデータ群に対して、ボクセルデータ群の中心を原点とした任意の回転軸を設定し、任意の回転角度だけ回転させる。
【００７６】
(３)(１)で得た回転させていない方（基準立体画像）のボクセルデータ群と、(２)で得た回転させた方のボクセルデータ群（比較対象立体画像）を３次元相関処理により、相関値を求める。相関値を求めるにには、３次元相関処理として、３次元位相限定相関法を含む任意の相関処理を用いることができる。
【００７７】
なお、ボクセルデータ群に対する３次元フーリエ変換の処理は、必ずしも必要というわけではない。例えば、被写体の回転角度が比較的小さいことが分かっている場合、回転中心がほぼ中心付近にあることが分かっている場合などは、３次元フーリエ変換の処理を行う必要はなく、基準立体画像および比較対象立体画像の原画像を構成するボクセルデータ群をそのまま使って、相関処理を行うことが可能である。
【００７８】
上述した実施の形態では、基準立体画像の原画像および比較対象立体画像の原画像のデータサイズを小さくする方式としているので、多少の回転中心のずれは無視することができる。また、被写体Ａ，Ｂの撮影に際し、被写体Ａ，Ｂはその頭部の中心を中間点ＣＯ付近にきちんと位置させるものとしているので、回転中心がほぼ中心付近にあることが分かっている。このようなことから、上述した実施の形態では、３次元フーリエ変換の処理を省略している。
【００７９】
〔回転軸と回転角度の設定について〕
初期段階において、比較対象立体画像に対して設定するＫ個の回転軸は全く不明である。このため、比較対象立体画像のボクセルデータ群の中心を原点とした任意の回転軸は、３次元上あらゆる方向に均等に分割されるように回転軸を選ぶ必要がある。回転軸の設定には、正多面体を用いることができる。例として、正二十面体を考えると、正二十面体の各面の中心座標を回転軸とすることで、球面に対して均等な２０個の回転軸を設定することができる。また、正二十面体の頂点を回転軸とすることで、球面に対して均等な１２個の回転軸を設定することができる。上述した実施の形態では、正二十面体の頂点を回転軸とし、原点を中心とする任意の方向へのＫ個の回転軸として１２個の回転軸Ｊ１〜Ｊ１２（図８参照）を設定している。
【００８０】
回転角度については、対象物体が最大どの程度回転することがあるか、によって適宜設定すればよい。例えば、ＣＴスキャナなどにより撮像されたボクセルデータであり、撮像時にある程度の位置合わせを行っているのであれば、回転誤差は目視で確認できない程度の角度であり、±１０程度で十分であると考えられる。これを数度刻みで設定すればよい。また、回転軸を設定する際に、対向する回転軸を設定した場合には、この２つの軸に対する回転は、回転方向が逆になるだけなので、対向する回転軸を１つにする、または、回転角度を一方向にする、ことができる。このようなことを考慮して、上述した実施の形態では、回転軸Ｊ１〜Ｊ１２を中心とする角度範囲を２〜３２゜とし、２゜刻みで回転させるようにしている。
【００８１】
〔処理速度について〕
初期段階で、１２個の回転軸（Ｋ＝１２）、１６の回転角度（Ｌ＝１６）を選択した場合、１９２回の３次元回転と３次元相関処理が必要となる。これには膨大な演算を必要とするが、最初に求めるのは概略の回転軸と回転角度であるため、ボクセルデータ群を小さくしても問題はない。ボクセルデータ群を小さくすることで、演算量は少なくなり、高速処理が期待できる。このため、上述した実施の形態では、５１２×５１２×５１２の原画像のデータサイズを縮小し、１６×１６×１６の縮小したデータサイズとしている。
【００８２】
データサイズの一辺の数を１／２にした場合（ボクセル数では１／８）、通常のパーソナルコンピュータ（ＰＣ）にて３次元位相限定相関法を実行すると、処理時間は１／５から１／８程度になる。したがって、５１２×５１２×５１２のデータサイズのボクセルデータ群に対して３次元位相限定相関法を実行した場合の時間を１とすると、１６×１６×１６のデータサイズのボクセルデータ群では、１／３１２５から１／３２７６８となり、多数回実行しても問題ない処理スピードとなる。
【００８３】
〔合成ベクトルを作成する際の重み付けについて〕
合成ベクトルを作成する際に用いる重み付けは、相関処理の出力する値（相関度合いが高くなるほど高い値となる）をそのまま使用しても構わないが、相関処理の特性に合わせて任意の処理を行ってから使用することもある。例えば、相関処理の特性が緩慢な場合（完全に一致しても、少しずれていても、出力値が大きく変わらない場合）、出力値に相関関数を適用することもできるし、逆に急峻な場合には、平方根を適用することもできる。
【００８４】
〔機能ブロック図〕
図１５、図１６および図１７に上述した３次元ずれ量計測装置４の要部の機能ブロック図を分割して示す。３次元ずれ量計測装置４は、基準立体画像の原画像（データサイズ：５１２×５１２×５１２）を記憶するデータベースＤＢ１と、比較対象立体画像の原画像（データサイズ：５１２×５１２×５１２）を記憶するデータベースＤＢ２とを備えている。
【００８５】
また、第１ステージ用として、第１の基準立体画像データサイズ縮小部４−１と、第１の比較対象立体画像データサイズ縮小部４−２と、第１の比較候補回転軸設定部４−３と、第１の比較候補立体画像作成部４−４と、第１の相関値算出部４−５と、第１の最大相関値抽出部４−６と、第１ステージ比較確定回転角度決定部４−７と、相関値抽出部４−８と、第１ステージ比較回転軸決定部４−９とを備えている。
【００８６】
また、第２ステージ用として、第２の比較候補立体画像作成部４−１０と、第２の相関値算出部４−１１と、第２の最大相関値抽出部４−１２と、第２ステージ比較確定回転角度決定部４−１３と、第２の比較候補回転軸設定部４−１４と、第３の比較候補立体画像作成部４−１５と、第３の相関値算出部４−１６と、第３の最大相関値抽出部４−１７と、第２ステージ比較回転軸決定部４−１８とを備えている。なお、第１ステージ用として設けた第１の基準立体画像データサイズ縮小部４−１と第１の比較対象立体画像データサイズ縮小部４−２は、第２ステージ用としても用いられる。
【００８７】
また、第３ステージ用として、第２の基準立体画像データサイズ縮小部４−１９と、第２の比較対象立体画像データサイズ縮小部４−２０と、第４の比較候補立体画像作成部４−２１と、第４の相関値算出部４−２２と、第４の最大相関値抽出部４−２３と、第３ステージ比較確定回転角度決定部４−２４と、第３の比較候補回転軸設定部４−２５と、第５の比較候補立体画像作成部４−２６と、第５の相関値算出部４−２７と、第５の最大相関値抽出部４−２６と、第３ステージ比較回転軸決定部４−２９とを備えている。
【００８８】
第１の基準立体画像データサイズ縮小部４−１は、データベースＤＢ１に保存されている５１２×５１２×５１２の基準立体画像（原画像）ＧＲのデータサイズを縮小し、第１ステージおよび第２ステージで用いる基準立体画像（データサイズ：１６×１６×１６）ＧＲ１とする。第１の比較対象立体画像データサイズ縮小部４−２は、データベースＤＢ２に保存されている５１２×５１２×５１２の比較対象立体画像（原画像）ＧＣのデータサイズを縮小し、第１ステージおよび第２ステージで用いる比較対象立体画像（データサイズ：１６×１６×１６）ＧＣ１とする。
【００８９】
第１の比較回転軸設定部４−３は、第１ステージの比較対象立体画像ＧＣ１のＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ個の回転軸を比較候補回転軸Ｊ１〜ＪＫとして設定する。第１の比較候補立体画像作成部４−４は、第１の比較回転軸設定部４−３によって作成されたＫ個の比較候補回転軸を中心として第１ステージの比較対象立体画像ＧＣ１を所定の角度範囲内でＬ通り回転させたＫ×Ｌ個の画像を比較候補立体画像として作成する。第１の相関値算出部４−５は、第１の比較候補立体画像作成部４−４によって作成されたＫ×Ｌ個の比較候補立体画像と第１ステージの基準立体画像ＧＲ１との相関値を各個に算出する。
【００９０】
第１の最大相関値抽出部４−６は、第１の相関値算出部４−５によって算出されたＫ×Ｌ個の相関値の中から最大相関値を抽出する。第１ステージ比較確定回転角度決定部４−７は、最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして決定する。
【００９１】
相関値抽出部４−８は、第１の相関値算出部４−５によって算出されたＫ×Ｌ個の相関値の中から、第１ステージでの比較確定回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する。
【００９２】
第１ステージ比較確定回転軸決定部４−９は、相関値抽出部４−８によって抽出された相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル（原点からの単位長さ）に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像ＧＣの第１ステージでの比較確定回転軸ＪＡとする。
【００９３】
第２の比較候補立体画像作成部４−１０は、第１ステージでの比較確定回転軸ＪＡを中心とし、第１ステージでの比較確定回転角度θＡを含む所定の角度範囲内で比較対象立体画像ＧＣ１をＭ１通り回転させたＭ１個の画像を比較候補立体画像として作成する。第２の相関値算出部４−１１は、第２の比較候補立体画像作成部４−１０によって作成されたＭ１個の比較候補立体画像と基準立体画像ＧＲ１との相関値を各個に算出する。
【００９４】
第２の最大相関値抽出部４−１２は、第２の相関値算出部４−１１によって算出されたＭ１個の相関値の中から最大相関値を抽出する。第２ステージ比較確定回転角度決定部４−１３は、最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第２ステージでの比較確定回転角度θＢとして決定する。
【００９５】
第２の比較候補回転軸設定部４−１４は、第１ステージでの比較確定回転軸ＪＡを中心にＮ１個の比較候補回転軸を設定する。第３の比較候補立体画像作成部４−１５は、第２の比較候補回転軸設定部４−１４によって設定されたＮ１個の比較候補回転軸のそれぞれについて、その回転軸を中心として比較確定回転角度θＢ回転させたＮ１個の比較候補立体画像を作成する。第３の相関値算出部４−１６は、第３の比較候補立体画像作成部４−１５によって作成されたＮ１個の比較候補立体画像と基準立体画像ＧＲ１との相関値を各個に算出する。
【００９６】
第３の最大相関値抽出部４−１７は、第３の相関値算出部４−１６によって算出されたＮ１個の相関値の中から最大相関値を抽出する。第２ステージ比較確定回軸決定部４−１８は、最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとして決定する。
【００９７】
第２の基準立体画像データサイズ縮小部４−１９は、データベースＤＢ１に保存されている５１２×５１２×５１２の基準立体画像（原画像）ＧＲのデータサイズを縮小し、第３ステージで用いる基準立体画像（データサイズ：６４×６４×６４）ＧＲ２とする。第１の比較対象立体画像データサイズ縮小部４−２０は、データベースＤＢ２に保存されている５１２×５１２×５１２の比較対象立体画像（原画像）ＧＣのデータサイズを縮小し、第３ステージで用いる比較対象立体画像（データサイズ：６４×６４×６４）ＧＣ２とする。
【００９８】
第４の比較候補立体画像作成部４−２１は、第２ステージでの比較確定回転軸ＪＢを中心とし、第２ステージでの比較確定回転角度θＢを含む所定の角度範囲内で比較対象立体画像ＧＣ２をＭ２通り回転させたＭ２個の画像を比較候補立体画像として作成する。第４の相関値算出部４−２２は、第４の比較候補立体画像作成部４−２１によって作成されたＭ２個の比較候補立体画像と基準立体画像ＧＲ２との相関値を各個に算出する。
【００９９】
第４の最大相関値抽出部４−２３は、第４の相関値算出部４−２２によって算出されたＭ２個の相関値の中から最大相関値を抽出する。第３ステージ比較確定回転角度決定部４−２４は、最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第３ステージでの比較確定回転角度θＣとして決定する。
【０１００】
第３の比較候補回転軸設定部４−２５は、第２ステージでの比較確定回転軸ＪＢを中心にＮ２個の比較候補回転軸を設定する。第５の比較候補立体画像作成部４−２６は、第３の比較候補回転軸設定部４−２５によって設定されたＮ２個の比較候補回転軸のそれぞれについてその回転軸を中心として回転角度θＣ回転させたＮ２個の比較候補立体画像を作成する。、第５の相関値算出部４−２７は、第５の比較候補立体画像作成部４−２６によって作成されたＮ２個の比較候補立体画像と基準立体画像ＧＲ２との相関値を各個に算出する。
【０１０１】
第５の最大相関値抽出部４−２８は、第５の相関値算出部４−２７によって算出されたＮ２個の相関値の中から最大相関値を抽出する。第３ステージ比較確定回軸決定部４−２９は、最大相関値が得られた比較候補立体画像の比較候補回転軸を第３ステージでの比較確定回転軸ＪＣとして決定する。
【０１０２】
なお、上述した実施の形態では、最初に撮影した被写体Ａの３次元の立体画像を基準立体画像とし、次に撮影した被写体Ｂの３次元の立体画像を比較対象立体画像としたが、最初に撮影した被写体Ａの３次元の立体画像を比較対象立体画像とし、次に撮影した被写体Ｂの３次元の立体画像を基準立体画像としてもよい。また、基準立体画像や比較対象立体画像は、ボクセルデータ群として上位装置から与えられるものとしてもよい。
【０１０３】
また、上述した実施の形態では、被写体を人体としたが、被写体は人体に限られるものではない。例えば、製薬などの分野において、マウスを使用した実験などにも利用することができる。また、産業分野において、金属などの非破壊検査に利用することもできる。このように、本発明は、医療や歯科などの分野に限らず、各種の分野において有効利用することが可能である。
【０１０４】
また、上述した３次元ずれ量計測処理装置４において、図１５〜図１７に分割して機能ブロックとして示した各部の処理動作は、プログラム（３次元ずれ量計測プログラム）に従うＣＰＵの処理動作として実現されるが、同様の処理をハードウェアの回路構成で実現するようにしてもよい。
【０１０５】
また、上述した実施の形態では、基準立体画像と比較対象立体画像との間のずれ量の計測処理を第３ステージを最終ステージとして実行するようにしたが、第１ステージで終了するようにしてもよく、第２ステージで終了するようにしてもよい。また、第１ステージや第２ステージで終了するような場合、基準立体画像および比較対象立体画像の原画像のデータサイズは、必ずしも縮小しなくてもよい。
【図面の簡単な説明】
【０１０６】
【図１】本発明に係る３次元ずれ量計測方法の実施に用いるＸ線コンピュータ断層撮影システムの一例を示す図である。
【図２】このシステムの３次元ずれ量計測装置におけるＣＰＵの処理動作（基準立体画像の取得処理）を示すフローチャートである。
【図３】このシステムの３次元ずれ量計測装置におけるＣＰＵの処理動作（比較対象立体画像の取得処理）を示すフローチャートである。
【図４】このシステムの３次元ずれ量計測装置におけるＣＰＵの処理動作（基準立体画像と比較対象立体画像との間のずれ量の計測処理）を示すフローチャートである。
【図５】図４に続くフローチャートである。
【図６】図５に続くフローチャートである。
【図７】このシステムの３次元ずれ量計測装置におけるＣＰＵの処理動作（位置合わせ処理＋照合処理）を示すフローチャートである。
【図８】正二十面体の中心と各頂点とを結ぶ軸を比較候補回転軸として設定する状況を示す図である。
【図９】各比較候補回転軸を中心として所定の角度範囲内で回転させた比較候補立体画像と第１ステージの基準立体画像との相関値との関係を示す図である。
【図１０】各比較候補回転軸の単位ベクトル（原点からの単位長さ）に相関値に基づく重み付けを施して得られる合成ベクトルＶＡおよひ第１ステージの比較確定回転軸ＪＡを説明する図である。
【図１１】第１ステージでの比較確定回転軸ＪＡを中心とし第１ステージでの比較確定回転角度θＡを中心とする所定の角度範囲内（±１０゜）で比較対象立体画像をＭ１通り回転させる状況を説明する図である。
【図１２】第１ステージでの比較確定回転軸ＪＡを中心軸としその周囲に設定されるＮ１個の比較候補回転軸を説明する図である。
【図１３】第２ステージでの比較確定回転軸ＪＢを中心とし第２ステージでの比較確定回転角度θＢを中心とする所定の角度範囲内（±１゜）で比較対象立体画像をＭ２通り回転させる状況を説明する図である。
【図１４】第２ステージでの比較確定回転軸ＪＡを中心軸としそのその周囲に設定されるＮ２個の比較候補回転軸を説明する図である。
【図１５】このシステムにおける３次元ずれ量計測装置の要部の機能ブロック図である。
【図１６】図１５に続く３次元ずれ量計測装置の要部の機能ブロック図である。
【図１７】図１６に続く３次元ずれ量計測装置の要部の機能ブロック図である。
【図１８】Ｘ線コンピュータ断層撮影法による３次元の立体画像の取得例を説明する図である。
【符号の説明】
【０１０７】
１…Ｘ線源、２…２次元撮像装置、３（Ａ，Ｂ）…被写体、４Ａ…ＣＰＵ、４Ｂ…ＲＡＭ、４Ｃ…ＲＯＭ、４Ｄ…記憶装置、４Ｅ〜Ｈ…インタフェース、４Ｉ…ディスプレイ、４Ｊ…キーボード、４Ｋ…マウス、ＤＢ１，ＤＢ２…データベース、４−１…第１の基準立体画像データサイズ縮小部、４−２…第１の比較対象立体画像データサイズ縮小部、４−３…第１の比較候補回転軸設定部、４−４…第１の比較候補立体画像作成部、４−５…第１の相関値算出部、４−６…第１の最大相関値抽出部、４−７…第１ステージ比較確定回転角度決定部、４−８…相関値抽出部、４−９…第１ステージ比較回転軸決定部、４−１０…第２の比較候補立体画像作成部、４−１１…第２の相関値算出部、４−１２…第２の最大相関値抽出部、４−１３…第２ステージ比較確定回転角度決定部、４−１４…第２の比較候補回転軸設定部、４−１５…第３の比較候補立体画像作成部、４−１６…第３の相関値算出部、４−１７…第３の最大相関値抽出部、４−１８…第２ステージ比較回転軸決定部、４−１９…第２の基準立体画像データサイズ縮小部、４−２０…第２の比較対象立体画像データサイズ縮小部、４−２１…第４の比較候補立体画像作成部、４−２２…第４の相関値算出部、４−２３…第４の最大相関値抽出部、４−２４…第３ステージ比較確定回転角度決定部、４−２５…第３の比較候補回転軸設定部、４−２６…第５の比較候補立体画像作成部、４−２７…第５の相関値算出部、４−２６…第５の最大相関値抽出部、４−２９…第３ステージ比較回転軸決定部。
【技術分野】
【０００１】
この発明は、基準となる３次元の立体画像を基準立体画像、この基準立体画像に対して比較対象とする３次元の立体画像を比較対象立体画像とし、基準立体画像と比較対象立体画像との間のずれ量を計測する３次元ずれ量計測方法に関するものである。
【背景技術】
【０００２】
従来より、医学分野などにおいては、多様な画像化方式によって立体画像を得ている。例えば、磁気共鳴断層撮影法（ＭＲＩ）、Ｘ線コンピュータ断層撮影法（ＣＴ）、陽電子放射断層撮影法（ＰＥＴ）、単光子射出断層撮影法（ＳＰＥＣＴ）、超音波断層撮影法などによって、３次元の立体画像を得ている。
【０００３】
このような方法で得られた３次元の立体画像は、ボクセル（voxel）値の配列として、通常、コンピュータのメモリなどにデジタル的に保存される。各ボクセルは、互いに直交するＸ，Ｙ，Ｚの３軸で表される３次元空間の座標位置（ｘ，ｙ，ｚ）と結びついており、ボクセル値として色値（典型的にはグレイスケール）が割り当てられる。
【０００４】
図１８にＸ線コンピュータ断層撮影法による３次元の立体画像の取得例を示す。この例では、スキャナ装置としてコーンビームＣＴを用いている。コーンビームＣＴは、中間点Ｃ０を挾んで対向してＸＹ平面内に設けられたＸ線源１と２次元撮像装置２とを備え、このＸ線源１と２次元撮像装置２とを中間点Ｃ０を通るＺ軸を中心として水平方向に定速度で回転させながら、被写体（人体）３の撮影を行う。この撮影に際して、被写体３は、その頭部の中心を中間点Ｃ０に位置させる。
【０００５】
被写体３の撮影は、Ｘ線源１と２次元撮像装置２とが１回転する間、所定の時間間隔で行われる。この結果、１回転する間に任意の異なる角度（以下、この角度を投影角度と呼ぶ）から撮影されたｎ個（例えば、２５６個）の２次元の投影画像を取得して、このｎ個の投影画像を再構成処理することで３次元の立体画像を得ることができる。この３次元の立体画像は、互いに直交するＸ，Ｙ，Ｚの３軸で表される３次元空間の座標位置（ｘ，ｙ，ｚ）と結びつけたボクセルデータ群として表現される。
【０００６】
この３次元の立体画像は、医療や歯科などの分野において、各種の診断を行うために使用される。しかし、情報量が多いために見落としが発生する虞があり、正確な診断を行うために過大な時間を必要とする問題がある。このため、コンピュータ上で、画像処理を行うことで、この問題を解決しようと試みられている。
【０００７】
この場合の画像処理の１つに、レジストレーションと呼ばれる位置合わせを行った後、２つの立体画像を比較するという方法がある。例えば、前回撮影した被写体３（被写体Ａ）の３次元の立体画像を基準立体画像とし、今回撮影した被写体３（被写体Ｂ）の３次元の立体画像を比較対象立体画像とし、基準立体画像と比較対象立体画像との位置合わせを行い、この位置合わせが行われた２つの立体画像の差分画像をディスプレイ上に表示する。これにより、前回撮影時と今撮影時との間の変化部分のみを抽出し、短時間でかつ正確な診断を行うことができるようになる。
【０００８】
この方法において、基準立体画像と比較対象立体画像とのレジストレーションは、６つのパラメータを求めることで実現される。６つのパラメータは、Ｘ，Ｙ，Ｚの各軸の平行ずれ量Δｘ，Δｙ，Δｚと、Ｘ，Ｙ，Ｚの各軸の回転ずれ量Δθ，Δφ，ΔΨである。平行ずれ量Δｘ，Δｙ，Δｚについては、回転ずれ量Δθ，Δφ，ΔΨを求めた後に、３次元相関処理を行うことで、容易に求めることができる。回転ずれ量Δθ，Δφ，ΔΨは、Ｘ，Ｙ，Ｚの各軸の回転角度を３つの自由度とし、この３つの自由度を変化させながら３次元相関処理を行うことで求めることが可能である（例えば、特許文献１参照）。
【０００９】
【特許文献１】特開平９−２２４０６号公報
【特許文献２】特開２００１−２１２１２６号公報
【発明の開示】
【発明が解決しようとする課題】
【００１０】
しかしながら、上述した基準立体画像と比較対象立体画像とのレジストレーションでは、回転ずれ量Δθ，Δφ，ΔΨを求める際、Ｘ，Ｙ，Ｚの各軸の回転角度を自由度として変化させるため、３乗の演算回数が必要となる。例えば、３つの自由度それぞれに１０種類の回転角度を設定した場合、１０×１０×１０＝１０００回の３次元の相関処理が必要となる。このため、膨大な演算を必要とし、高速処理が困難となる。
【００１１】
なお、相関処理を用いないレジストレーションとして、立体画像に現れている被写体の構造に詳しい専門家が、位置合わせする立体画像のそれぞれに目印をつけ、一方の立体画像を移動回転させることによって位置合わせする方法がある（例えば、特許文献２参照）。
【００１２】
しかしながら、この方法では、位置合わせ精度が選択された目印の数と位置によって制限されてしまう。また、目印の数が少な過ぎると、位置合わせが不正確になる。また、目印が多くても必ずしも位置合わせが正確となるわけではなく、位置合わせに要する計算を複雑にするだけである。また、オペレータが手作業で対応点を選択する必要があり、時間がかかるとともに熟練が必要となる。さらに、すべての立体画像で適当な構造上の目印を見出せるわけでもない。
【００１３】
本発明は、このような課題を解決するためになされたもので、その目的とするところは、少ない演算量で、高速かつ正確に、レジストレーションを行うことが可能な３次元ずれ量計測方法を提供することにある。
【課題を解決するための手段】
【００１４】
このような目的を達成するために本発明は、基準となる３次元の立体画像を基準立体画像、この基準立体画像に対して比較対象とする３次元の立体画像を比較対象立体画像とし、基準立体画像と比較対象立体画像との間のずれ量を計測する３次元ずれ量計測方法において、比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ（Ｋ≧２）個の回転軸を比較候補回転軸として設定し（第１ステップ）、Ｋ個の各比較候補回転軸を中心として比較対象立体画像を所定角度範囲内でＬ（Ｌ≧２）通り回転させたＫ×Ｌ個の画像を比較候補立体画像とし、このＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値を各個に算出し（第２ステップ）、この算出されたＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして求め（第３ステップ）、算出されたＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から比較確定回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出し（第４ステップ）、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル（原点からの単位長さ）に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像の第１ステージでの比較確定回転軸ＪＡとする（第５ステップ）、ようにしたものである。
【００１５】
本発明を適用しようとする３次元空間において、その３次元空間における回転角度のずれには、様々な表現方法が考えられる。最も簡単に表現できるのは、互いに直交するＸ，Ｙ，Ｚの３軸のそれぞれを回転軸とした３つの回転角度により表現することである。すなわち、回転軸を３つとし、これらの回転軸を中心としてそれぞれある角度にて回転したと考え、合計３回の回転が発生していると考えることである。本発明では、この考えをさらに進め、原点を中心として複数回の回転を行うことは、ある回転軸を中心として１回の回転を行ったに過ぎないと考え、回転軸と回転角度の２つのパラメータを求めることにする。
【００１６】
この考えに基づき、本発明では、比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ（Ｋ≧２）個の回転軸を比較候補回転軸として設定する（第１ステップ）。例えば、原点を中心とする正二十面体を考え、この正二十面体の中心と各頂点とを結ぶ１２個の軸を比較候補回転軸として設定する。そして、このＫ個の各比較候補回転軸を中心として比較対象立体画像を所定角度範囲内でＬ（Ｌ≧２）通り回転させたＫ×Ｌ個の画像（例えば、２〜３２゜の角度範囲内で２゜刻みで回転させた１２×１６＝１９２個の画像）を比較候補立体画像とし、このＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値を各個に算出する（第２ステップ）。
【００１７】
そして、この算出したＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして求める（第３ステップ）。また、この算出したＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から比較確定回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する（第４ステップ）。そして、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル（原点からの単位長さ）に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像の第１ステージでの比較確定回転軸ＪＡとする（第５ステップ）。
【００１８】
このようにして、本発明では、第１ステージでの比較確定回転軸ＪＡと比較確定回転角度θＡを求める。この場合、第１ステージで求められた比較確定回転軸ＪＡと比較確定回転角度θＡを用いて基準立体画像と比較対象立体画像との位置合わせ（レジストレーション）を行うようにしてもよいが、レジストレーションの精度をさらに高めるためには、さらに第２ステージに進むことが望ましい。
【００１９】
第２ステージでは、第１ステージでの比較確定回転軸ＪＡを中心とし、第１ステージでの比較確定回転角度θＡを含む所定角度範囲内で比較対象立体画像をＭ（Ｍ≧２）通り回転させたＭ個の画像（例えば、比較確定回転角度θＡを中心とする±１０゜の角度範囲内で０．５゜刻みで回転させた４１個の画像）を比較候補立体画像とし、このＭ個の比較候補立体画像と基準立体画像との相関値を各個に算出する（第６ステップ）。そして、この算出したＭ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとして求める（第７ステップ）。
【００２０】
そして、第１ステージでの比較確定回転軸ＪＡを中心に任意のＮ個（例えば、１１×１１＝１２１個）の比較候補回転軸を設定し、このＮ個の比較候補回転軸のそれぞれについてその回転軸を中心として比較確定回転角度θＢ回転させたＮ個の比較候補立体画像と基準立体画像との相関値を各個に算出する（第８ステップ）。そして、このＮ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとする（第９ステップ）。
【００２１】
なお、上述した第７ステップにおいて、Ｍ個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとして求めるようにしてもよい。また、上述した第９ステップにおいて、Ｎ個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとするようにしてもよい。
【００２２】
また、本発明において、第２ステージに進む場合、第１ステップに入る前に、基準立体画像および比較対象立体画像の原画像のデータサイズを小さくして第１ステージで用いる基準立体画像および比較対象立体画像とし、第１ステージで用いる基準立体画像および比較対象立体画像のデータサイズよりもそのデータサイズを大きくしながら、第６〜第９ステップの処理を予め定められた最終ステージまで繰り返すようにしてもよい。
【００２３】
例えば、基準立体画像および比較対象立体画像の原画像のデータサイズ（画素数）を共に５１２×５１２×５１２とした場合、これら立体画像のデータサイズを縮小して１６×１６×１６の立体画像とする。第１ステージでは、この１６×１６×１６の基準立体画像および比較対象立体画像を用いて第１〜第５ステップの処理を行う。第２ステージでも、第１ステージと同様にして、１６×１６×１６の基準立体画像および比較対象立体画像を用いて第６〜第９ステップの処理を行う。そして、第３ステージから、基準立体画像および比較対象立体画像のデータサイズを大きくし（例えば、６４×６４×６４とする）、第６〜第９ステップの処理を行う。このようにして、最終ステージまで、基準立体画像および比較対象立体画像のデータサイズを大きくしながら、第６〜第９ステップの処理を繰り返す。なお、第２ステージから基準立体画像および比較対象立体画像のデータサイズを大きくするようにしてもよく、第３ステージを最終ステージとしてもよい。
【００２４】
このようにして、本発明では、基準立体画像および比較対象立体画像のデータサイズを最初は小さくし、ステージが進むにつれて大きくして行くことにより、最初は少ないボクセルデータ群で概略的な比較確定回転軸と比較確定回転角度を素速く求め、徐々にボクセルデータ群の数を増やして、より正確な比較確定回転軸と比較確定回転角度を求めて行くようにして、高速度でかつ高精度に、レジストレーションを行うことが可能となる。
【００２５】
本発明では、比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ個の回転軸を比較候補回転軸として設定し、このＫ個の比較候補回転軸を中心として回転させた比較対象立体画像と基準立体画像との相関値を求めるが、このような処理において比較対象立体画像と基準立体画像の回転中心とされる原点は同じでなければならない。しかし、３次元の立体画像において、実際の回転中心は不明であるため、原点を決めることが難しい。
【００２６】
本発明において、正確な原点を決めるために、第１ステップの処理に入る前に、基準立体画像の原画像および比較対象立体画像の原画像を構成するボクセルデータ群に対してそれぞれ３次元フーリエ変換を施し、この３次元フーリエ変換が施された基準立体画像の原画像および比較対象立体画像の原画像を構成するボクセルデータ群の３次元周波数空間（ボクセルデータ群の中心を直流とする周波数空間）の中心を基準立体画像および比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心とするようにしてもよい。
【００２７】
なお、ボクセルデータ群に対する３次元フーリエ変換の処理は、必ずしも必要というわけではない。例えば、被写体の回転角度が比較的小さいことが分かっている場合、回転中心がほぼ中心付近にあることが分かっている場合などは、３次元フーリエ変換の処理を行う必要はなく、基準立体画像および比較対象立体画像の原画像を構成するボクセルデータ群をそのまま使って、相関処理を行うことが可能である。本発明において、第１ステップの処理に入る前に基準立体画像の原画像および比較対象立体画像の原画像のデータサイズを小さくする方式とした場合、多少の回転中心のずれは無視することができる。
【発明の効果】
【００２８】
本発明によれば、比較対象立体画像に設定される原点を中心として任意の方向へＫ個の比較候補回転軸を設定し、このＫ個の各比較候補回転軸を中心として比較対象立体画像を所定角度範囲内でＬ通り回転させたＫ×Ｌ個の画像を比較候補立体画像とし、このＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値を各個に算出し、算出したＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして求め、算出したＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から比較確定回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出し、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル（原点からの単位長さ）に当該相関値に基づく重み付けを施し、この重み付けが施された単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像の第１ステージでの比較確定回転軸ＪＡとするようにしたので、この第１ステージで求められた比較確定回転軸ＪＡおよび比較確定回転角度θＡを用いて基準立体画像と比較対象立体画像との位置合わせを行うことにより、少ない演算量で、高速かつ正確に、レジストレーションを行うことが可能となる。
【００２９】
また、本発明によれば、第１ステージでの比較確定回転軸ＪＡを中心とし、第１ステージでの比較確定回転角度θＡを含む所定角度範囲内で比較対象立体画像をＭ通り回転させたＭ個の画像を比較候補立体画像とし、このＭ個の比較候補立体画像と基準立体画像との相関値を算出し、算出したＭ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとして求め、第１ステージでの比較確定回転軸ＪＡを中心に任意のＮ（Ｎ≧２）個の比較候補回転軸を設定し、このＮ個の比較候補回転軸のそれぞれについてその回転軸を中心として比較確定回転角度θＢ回転させたＮ個の比較候補立体画像と基準立体画像との相関値を算出し、Ｎ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとするようにしたので、この第２ステージで求められた比較確定回転軸ＪＢおよび比較確定回転角度θＢを用いて基準立体画像と比較対象立体画像との位置合わせを行うことにより、さらに正確に、レジストレーションを行うことが可能となる。
【００３０】
また、本発明によれば、第１ステップの処理に入る前に、基準立体画像および比較対象立体画像の原画像のデータサイズを小さくして第１ステージで用いる基準立体画像および比較対象立体画像とし、第１ステージで用いる基準立体画像および比較対象立体画像のデータサイズよりもそのデータサイズを大きくしながら第６〜第９ステップの処理を予め定められた最終ステージまで繰り返すことにより、最初は少ないボクセルデータ群で概略的な比較確定回転軸と比較確定回転角度を素速く求め、徐々にボクセルデータ群の数を増やして、より正確な比較確定回転軸と比較確定回転角度を求めて行くようにして、高速度でかつ高精度に、レジストレーションを行うことが可能となる。
【発明を実施するための最良の形態】
【００３１】
以下、本発明を図面に基づいて詳細に説明する。図１は本発明に係る３次元ずれ量計測方法の実施に用いるＸ線コンピュータ断層撮影システムの一例を示す図である。同図において、図１８と同一符号は図１８を参照して説明した構成要素と同一或いは同等構成要素を示し、その説明は省略する。
【００３２】
この実施の形態において、Ｘ線源１および２次元撮像装置２は、３次元ずれ量計測装置４に接続されている。３次元ずれ量計測装置４は、ＣＰＵ４Ａと、ＲＡＭ４Ｂと、ＲＯＭ４Ｃと、ハードディスクなどの記憶装置４Ｄと、インタフェース４Ｅ〜Ｈと、ディスプレイ４Ｉと、キーボード４Ｊと、マウス４Ｋとを備えている。
【００３３】
ＣＰＵ４Ａは、ＲＡＭ４Ｂにアクセスしながら、ＲＯＭ４Ｃや記憶装置４Ｄに格納されたプログラムに従って動作する。記憶装置４Ｄには、本実施の形態特有のプログラムとして３次元ずれ量計測プログラムが格納されている。この３次元ずれ量計測プログラムは、例えばＣＤ−ＲＯＭなどの記録媒体に記録された状態で提供され、この記録媒体から読み出されて記憶装置４Ｄにインストールされている。
【００３４】
以下、図２〜図７に示すフローチャートを参照して、記憶装置４Ｄに格納された３次元ずれ量計測プログラムに従ってＣＰＵ４Ａが実行する本実施の形態特有の処理動作について説明する。
【００３５】
なお、この３次元ずれ量計測プログラムに従う処理動作には、比較対象立体画像の基準立体画像に対する回転軸と回転角度の２つのパラメータ（ＪEND，ΔθEND）を求めるずれ量計測処理と、この求めた２つのパラメータを用いて基準立体画像と比較対象立体画像との間の位置合わせを行う位置合わせ処理と、この位置合わせされた基準立体画像と比較対象立体画像との照合を行う照合処理とが含まれる。
【００３６】
〔基準立体画像の取得〕
先ず、基準立体画像を取得するために、被写体Ａの撮影を行う（図２：ステップＳ１０１）。この被写体Ａの撮影は、Ｘ線源１と２次元撮像装置２とが１回転する間、所定の時間間隔で行う。この結果、１回転する間の任意の異なる角度（投影角度）から撮影されたｎ個（例えば、２５６個）の２次元の投影画像が得られる（ステップＳ１０２）。なお、この撮影において、被写体Ａは、その頭部の中心を中間点ＣＯ付近にきちんと位置させるものとする。
【００３７】
ＣＰＵ４Ａは、この得られたｎ個の２次元の投影画像を記憶装置４Ｄに保存した後（ステップＳ１０３）、このｎ個の投影画像を再構成処理することで３次元の立体画像（ボクセルデータ群）を得る（ステップＳ１０４）。そして、この再構成処理によって得られた３次元の立体画像を基準立体画像の原画像として記憶装置４Ｄに保存する（ステップＳ１０５）。なお、この場合、記憶装置４Ｄに保存される基準立体画像の原画像のデータサイズ（画素数）は、５１２×５１２×５１２とする。
【００３８】
〔比較対象立体画像の取得〕
次に、比較対象立体画像を取得するために、被写体Ｂの撮影を行う（図３：ステップＳ２０１）。この被写体Ｂの撮影も、Ｘ線源１と２次元撮像装置２とが１回転する間、基準立体画像の取得時と同じ所定の時間間隔で行う。この結果、１回転する間の任意の異なる角度（投影角度）から撮影されたｎ個（例えば、２５６個）の２次元の投影画像が得られる（ステップＳ２０２）。なお、この撮影において、被写体Ｂは、その頭部の中心を中間点ＣＯ付近にきちんと位置させるものとする。
【００３９】
ＣＰＵ４Ａは、この得られたｎ個の２次元の投影画像を記憶装置４Ｄに保存した後（ステップＳ２０３）、このｎ個の投影画像を再構成処理することで３次元の立体画像（ボクセルデータ群）を得る（ステップＳ２０４）。そして、この再構成処理によって得られた３次元の立体画像を比較対象立体画像の原画像として記憶装置４Ｄに保存する（ステップＳ２０５）。なお、この場合、記憶装置４Ｄに保存される比較対象立体画像の原画像のデータサイズ（画素数）は、５１２×５１２×５１２とする。
【００４０】
〔基準立体画像と比較対象立体画像との間のずれ量の計測〕
〔第１ステージ〕
次に、ＣＰＵ４Ａは、記憶装置４Ｄに保存されている５１２×５１２×５１２の基準立体画像（原画像）のデータサイズを縮小し、第１ステージで用いる基準立体画像とする（図４：ステップＳ３０１）。また、記憶装置４Ｄに保存されている５１２×５１２×５１２の比較対象立体画像（原画像）のデータサイズを縮小し、第１ステージで用いる比較対象立体画像とする（ステップＳ３０２）。
【００４１】
この例では、Ｘ，Ｙ，Ｚの３軸方向の画素データを所定画素ずつ足し合わせて行くことによって、第１ステージの基準立体画像および比較対象立体画像として１６×１６×１６のデータサイズの立体画像を得る。なお、画素データを間引いて行くことによって、データサイズを縮小するようにしてもよい。
【００４２】
そして、ＣＰＵ４Ａは、ステップ３０１で得た第１ステージの比較対象立体画像のＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ個の回転軸を比較候補回転軸として設定する（ステップＳ３０３）。この例では、原点を中心とする正二十面体（図８参照）を考え、この正二十面体の中心Ｏと各頂点とを結ぶ１２個の軸Ｊ１〜Ｊ１２を比較候補回転軸として設定する。
【００４３】
そして、このＫ個の比較候補回転軸を中心として第１ステージの比較対象立体画像を所定の角度範囲内でＬ通り回転させたＫ×Ｌ個の画像を比較候補立体画像とし（ステップＳ３０４）、このＫ×Ｌ個の比較候補立体画像と第１ステージの基準立体画像との相関値を各個に算出する（ステップＳ３０５）。
【００４４】
この例では、１２個の比較候補回転軸Ｊ１〜Ｊ１２について、その比較候補回転軸を中心として２〜３２゜の角度範囲内で２゜刻みで回転させた１２×１６＝１９２個の画像を比較候補立体画像とし、この１９２個の比較候補立体画像と第１ステージの基準立体画像との相関値ｍ１〜ｍ１９２（図９参照）を算出する。
【００４５】
そして、この算出した相関値ｍ１〜ｍ１９２の中から最大相関値を抽出し（ステップＳ３０６）、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして求める（ステップＳ３０７）。例えば、相関値ｍ１７が最大相関値であった場合、比較候補回転軸Ｊ２を中心とする回転角度θ１を第１ステージでの比較確定回転角度θＡとして求める。
【００４６】
また、ＣＰＵ４Ａは、ステップＳ３０５で算出した相関値ｍ１〜ｍ１９２の中から、第１ステージでの比較確定回転角度として求めた回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する（ステップＳ３０８）。すなわち、比較候補回転軸Ｊ１〜Ｊ１２について、回転角度θＡとした場合の比較候補立体画像の相関値を抽出する。
【００４７】
そして、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸Ｊ１〜Ｊ１２の単位ベクトル（原点からの単位長さ）に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルＶＡ（図１０参照）を求め、この合成ベクトルＶＡが示す方向を比較対象立体画像の第１ステージでの比較確定回転軸ＪＡとする（ステップＳ３０９）。
【００４８】
〔第２ステージ〕
次に、ＣＰＵ４Ａは、記憶装置４Ｄに保存されている５１２×５１２×５１２の基準立体画像（原画像）のデータサイズを縮小し、第２ステージで用いる基準立体画像とする（図５：ステップＳ３１０）。また、記憶装置４Ｄに保存されている５１２×５１２×５１２の比較対象立体画像（原画像）のデータサイズを縮小し、第２ステージで用いる比較対象立体画像とする（ステップＳ３１１）。
【００４９】
この例では、ステップ３１０，３１１は実行せずに、第１ステージで用いた１６×１６×１６の基準立体画像および比較対象立体画像を第２ステージの基準立体画像および比較対象立体画像として用いることとする。なお、ステップ３１０，３１１を実行する場合には、第１ステージで用いた基準立体画像および比較対象立体画像のデータサイズよりも大きいデータサイズに縮小し、第２ステージの基準立体画像および比較対象立体画像とする。
【００５０】
次に、ＣＰＵ４Ａは、第１ステージでの比較確定回転軸ＪＡを中心とし、第１ステージでの比較確定回転角度θＡを含む所定の角度範囲内で第２ステージの比較対象立体画像をＭ１通り回転させたＭ１個の画像を比較候補立体画像とし、このＭ１個の比較候補立体画像と第２ステージの基準立体画像との相関値を各個に算出する（ステップＳ３１２）。この例では、比較確定回転角度θＡを中心とする±１０゜の角度範囲内で０．５゜刻みで回転させた４１個の画像を比較候補立体画像とし（図１１参照）、この４１個の比較候補立体画像と第２ステージの基準立体画像との相関値を各個に算出する。
【００５１】
そして、この算出したＭ１個の比較候補立体画像と第２ステージの基準立体画像との相関値の中から最大相関値を抽出し（ステップＳ３１３）、この最大相関値が得られた比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとする（ステップＳ３１４）。
【００５２】
なお、この例では、Ｍ１個の比較候補立体画像と第２ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとしたが、Ｍ１個の比較候補立体画像と基準立体画像との相関値から関数フィッティング（例えば、シンク関数）により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとするようにしてもよい。
【００５３】
次に、ＣＰＵ４Ａは、第１ステージでの比較確定回転軸ＪＡを中心に任意のＮ１個の比較候補回転軸を設定する（ステップＳ３１５）。例えば、図１２に示すように、第１ステージでの比較確定回転軸ＪＡを中心軸とし、この比較確定回転軸ＪＡの周囲に縦方向１１個、横方向１１個、合計１２１個の原点Ｏからの比較候補回転軸を３゜刻みで設定する。そして、このＮ１個の比較候補回転軸のそれぞれについて、その回転軸を中心として比較確定回転角度θＢ回転させたＮ１個の比較候補立体画像と第２ステージの基準立体画像との相関値を各個に算出する（ステップＳ３１６）。
【００５４】
そして、このＮ１個の比較候補立体画像と第２ステージの基準立体画像との相関値の中から最大相関値を抽出し（ステップＳ３１７）、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとする（ステップＳ３１８）。
【００５５】
なお、この例では、Ｎ１個の比較候補立体画像と第２ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとするようにしたが、Ｎ１個の比較候補立体画像と第２ステージの基準立体画像との相関値から関数フィッティング（例えば、シンク関数）により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとするようにしてもよい。
【００５６】
〔第３ステージ（最終ステージ）〕
次に、ＣＰＵ４Ａは、記憶装置４Ｄに保存されている５１２×５１２×５１２の基準立体画像（原画像）のデータサイズを縮小し、第３ステージで用いる基準立体画像とする（図６：ステップＳ３１９）。また、記憶装置４Ｄに保存されている５１２×５１２×５１２の比較対象立体画像（原画像）のデータサイズを縮小し、第３ステージで用いる比較対象立体画像とする（ステップＳ３２０）。この例では、第３ステージで用いる基準立体画像および比較対象立体画像のデータサイズを６４×６４×６４とし、第１ステージで用いる基準立体画像および比較対象立体画像のデータサイズよりも大きいサイズとする。
【００５７】
次に、ＣＰＵ４Ａは、第２ステージでの比較確定回転軸ＪＢを中心とし、第２ステージでの比較確定回転角度θＢを含む所定の角度範囲内で第３ステージの比較対象立体画像をＭ２通り回転させたＭ２個の画像を比較候補立体画像とし、このＭ２個の比較候補立体画像と第３ステージの基準立体画像との相関値を各個に算出する（ステップＳ３２１）。この例では、比較確定回転角度θＢを中心とする±１゜の角度範囲内で０．１゜刻みで回転させた２１個の画像を比較候補立体画像とし（図１３参照）、この２１個の比較候補立体画像と第３ステージの基準立体画像との相関値を各個に算出する。
【００５８】
そして、この算出したＭ２個の比較候補立体画像と第３ステージの基準立体画像との相関値の中から最大相関値を抽出し（ステップＳ３２２）、この最大相関値が得られた比較候補立体画像の第２ステージでの比較確定回転軸ＪＢを中心とする回転角度を第３ステージでの比較確定回転角度θＣとする（ステップＳ３２３）。
【００５９】
なお、この例では、Ｍ２個の比較候補立体画像と第３ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の第２ステージでの比較確定回転軸ＪＢを中心とする回転角度を第３ステージでの比較確定回転角度θＣとしたが、Ｍ２個の比較候補立体画像と第３ステージの基準立体画像との相関値から関数フィッティング（例えば、シンク関数）により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の第２ステージでの比較確定回転軸ＪＢを中心とする回転角度を第３ステージでの比較確定回転角度θＣとするようにしてもよい。
【００６０】
次に、ＣＰＵ４Ａは、第２ステージでの比較確定回転軸ＪＢを中心に任意のＮ２個の比較候補回転軸を設定する（ステップＳ３２４）。例えば、図１４に示すように、第２ステージでの比較確定回転軸ＪＢを中心軸とし、この中心軸ＪＢの周囲に縦方向７個、横方向７個、合計４９個の原点Ｏからの比較候補回転軸を３゜刻みで設定する。そして、このＮ２個の比較候補回転軸のそれぞれについて、その回転軸を中心として回転角度θＣ回転させたＮ２個の比較候補立体画像と第３ステージの基準立体画像との相関値を各個に算出する（ステップＳ３２５）。
【００６１】
そして、このＮ２個の比較候補立体画像と第３ステージの基準立体画像との相関値の中から最大相関値を抽出し（ステップＳ３２６）、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第３ステージでの比較確定回転軸ＪＣとする（ステップＳ３２７）。
【００６２】
なお、この例では、Ｎ２個の比較候補立体画像と第３ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の比較候補回転軸を第３ステージでの比較確定回転軸ＪＣとするようにしたが、Ｎ２個の比較候補立体画像と第３ステージの基準立体画像との相関値から関数フィッティング（例えば、シンク関数）により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第３ステージでの比較確定回転軸ＪＣとするようにしてもよい。
【００６３】
本実施の形態では、第３ステージを最終ステージとし、この第３ステージでの比較確定回転軸ＪＣおよび比較確定回転角度θＣを最終ステージでの比較確定回転軸ＪENDおよび比較確定回転角度θENDとする。なお、第３ステージを最終ステージとせず、さらに最終ステージまでのステージ数を増やしてもよい。この場合、各ステージでは、第３ステージと同様にして、そのステージで用いる基準立体画像および比較対象立体画像のデータサイズを大きくしながら、第２ステージと同様の処理を繰り返すものとする。
【００６４】
〔基準立体画像と比較対象立体画像との間の位置合わせ〕
ＣＰＵ４Ａは、このようにして最終ステージでの比較確定回転軸ＪENDおよび比較確定回転角度θENDを求めた後、この求めた最終ステージでの比較確定回転軸ＪENDおよび比較確定回転角度θENDに基づいて、記憶装置４に保存されている基準立体画像（原画像）と比較対象立体画像（原画像）との間の回転ずれを補正する（図７：ステップＳ４０１）。すなわち、比較対象立体画像（原画像）を比較確定回転軸ＪENDを中心として比較確定回転角度θEND回転補正する。これにより、一度の処理で、基準立体画像と比較対象立体画像との間の回転ずれが補正される。
【００６５】
次に、この回転ずれが補正された基準立体画像と比較対象立体画像との間で相関演算を行って、この基準立体画像と比較対象立体画像との間のＸ，Ｙ，Ｚの３軸方向のずれ量Δｘ，Δｙ，Δｚを求め、この求めた３軸方向のずれ量Δｘ，Δｙ，Δｚに基づいて回転ずれが補正された基準立体画像と比較対象立体画像との間の平行ずれを補正する（ステップＳ４０２）。すなわち、Ｘ軸方向へΔｘ、Ｙ軸方向へΔｙ、Ｚ軸方向へΔｚ平行移動し、基準立体画像と比較対象立体画像とのＸ，Ｙ，Ｚの３軸方向の位置を合わせる。
【００６６】
これにより、基準立体画像と比較対象立体画像との間の位置合わせが完了する。なお、この際、比較対象立体画像ではなく、基準立体画像の方を回転補正したり、平行移動するようにしてもよい。
【００６７】
〔位置合わせされた基準立体画像と比較対象立体画像との照合〕
ＣＰＵ４Ａは、このようにして基準立体画像と比較対象立体画像との位置合わせを行った後、すなわち基準立体画像と比較対象立体画像とのレジストレーションを行った後、位置合わせされた基準立体画像と比較対象立体画像とを照合する（ステップＳ４０３）。
【００６８】
この例では、位置合わせされた基準立体画像（ボクセルデータ群）と比較対象立体画像（ボクセルデータ群）との相関演算を行い、この相関演算の結果得られた相関値と予め定められた照合用の閾値とを比較することによって、基準立体画像と比較対象立体画像との照合を行う。
【００６９】
なお、３次元の立体画像間の相関演算については特許文献１などに示されているので、ここでの説明は省略する。この特許文献１において、Ｎ＝３として３次元フーリエ変換、３次元逆フーリエ変換を施し、３次元の立体画像間の相関演算を行うことにより、相関成分エリア内でピークの位置ずれがＸ，Ｙ，Ｚの各軸に現れる。このＸ，Ｙ，Ｚの各軸に現れる位置ずれから基準立体画像と比較対象立体画像との相関が分かる。
【００７０】
ＣＰＵ４Ａは、ステップＳ４０２での照合の結果がＯＫであれば（ステップＳ４０４のＹＥＳ）、照合ＯＫの表示を行う（ステップＳ４０５）。照合の結果がＮＧであれば（ステップＳ４０４のＮＯ）、照合ＮＧの表示を行うとともに（ステップＳ４０６）、位置合わせされた基準立体画像と比較対象立体画像との差分画像をディスプレイ４Ｉ上に表示する（ステップＳ４０７）。
【００７１】
以上の説明から分かるように、本実施の形態によれば、基準立体画像および比較対象立体画像のデータサイズが最初は小さくされ、ステージが進むにつれて大きくされて行く。これにより、最初は少ないボクセルデータ群で概略的な比較確定回転軸と比較確定回転角度が素速く求められ、徐々にボクセルデータ群の数を増やして、より正確な比較確定回転軸と比較確定回転角度が求められて行く。このようにして、本実施の形態では、高速度でかつ高精度に、レジストレーションが行われる。この結果、照合スピードもアップし、照合結果の精度も高まる。
【００７２】
〔基準立体画像と比較対象立体画像の回転中心（原点）について〕
上述した実施の形態では、第１ステージの比較対象立体画像のＸ，Ｙ，Ｚの３軸の中心を原点としてＫ個の回転軸を比較候補回転軸として設定し、このＫ個の比較候補回転軸を中心として回転させた比較対象立体画像と第１ステージの基準立体画像との相関値を求めるが、このような処理において比較対象立体画像と基準立体画像の回転中心とされる原点は同じでなければならない。しかし、３次元の立体画像において、実際の回転中心は不明であるため、原点を決めることが難しい。
【００７３】
正確な原点を定めるために、本実施の形態において、ステップＳ３０１の処理に入る前に、記憶装置４Ｄに保存されている基準立体画像（原画像）および比較対象立体画像（原画像）を構成するボクセルデータ群に対してそれぞれ３次元フーリエ変換を施し、この３次元フーリエ変換が施された基準立体画像（原画像）および比較対象立体画像（原画像）を構成するボクセルデータ群の３次元周波数空間（ボクセルデータ群の中心を直流とする周波数空間）の中心を基準立体画像および比較対象立体画像に設定するＸ，Ｙ，Ｚの３軸の中心とすることが考えられる。
【００７４】
この場合、相関値を求める処理は、以下の手順となる。
(１)両方のボクセルデータ群に対し、３次元フーリエ変換（ＦＦＴ）、振幅化、を行う。この時、ボクセルデータ群の中心が直流になるようにする。また、３次元ＦＦＴした結果の振幅の値は、元のボクセルデータに比べて比較的大きな値をとることがある。これが演算上問題になる場合には、振幅データをログ化することで、ダイナミックレンジを抑えることができる。
【００７５】
(２)(１)で得た（必要に応じてログ化された）振幅データの片方のボクセルデータ群に対して、ボクセルデータ群の中心を原点とした任意の回転軸を設定し、任意の回転角度だけ回転させる。
【００７６】
(３)(１)で得た回転させていない方（基準立体画像）のボクセルデータ群と、(２)で得た回転させた方のボクセルデータ群（比較対象立体画像）を３次元相関処理により、相関値を求める。相関値を求めるにには、３次元相関処理として、３次元位相限定相関法を含む任意の相関処理を用いることができる。
【００７７】
なお、ボクセルデータ群に対する３次元フーリエ変換の処理は、必ずしも必要というわけではない。例えば、被写体の回転角度が比較的小さいことが分かっている場合、回転中心がほぼ中心付近にあることが分かっている場合などは、３次元フーリエ変換の処理を行う必要はなく、基準立体画像および比較対象立体画像の原画像を構成するボクセルデータ群をそのまま使って、相関処理を行うことが可能である。
【００７８】
上述した実施の形態では、基準立体画像の原画像および比較対象立体画像の原画像のデータサイズを小さくする方式としているので、多少の回転中心のずれは無視することができる。また、被写体Ａ，Ｂの撮影に際し、被写体Ａ，Ｂはその頭部の中心を中間点ＣＯ付近にきちんと位置させるものとしているので、回転中心がほぼ中心付近にあることが分かっている。このようなことから、上述した実施の形態では、３次元フーリエ変換の処理を省略している。
【００７９】
〔回転軸と回転角度の設定について〕
初期段階において、比較対象立体画像に対して設定するＫ個の回転軸は全く不明である。このため、比較対象立体画像のボクセルデータ群の中心を原点とした任意の回転軸は、３次元上あらゆる方向に均等に分割されるように回転軸を選ぶ必要がある。回転軸の設定には、正多面体を用いることができる。例として、正二十面体を考えると、正二十面体の各面の中心座標を回転軸とすることで、球面に対して均等な２０個の回転軸を設定することができる。また、正二十面体の頂点を回転軸とすることで、球面に対して均等な１２個の回転軸を設定することができる。上述した実施の形態では、正二十面体の頂点を回転軸とし、原点を中心とする任意の方向へのＫ個の回転軸として１２個の回転軸Ｊ１〜Ｊ１２（図８参照）を設定している。
【００８０】
回転角度については、対象物体が最大どの程度回転することがあるか、によって適宜設定すればよい。例えば、ＣＴスキャナなどにより撮像されたボクセルデータであり、撮像時にある程度の位置合わせを行っているのであれば、回転誤差は目視で確認できない程度の角度であり、±１０程度で十分であると考えられる。これを数度刻みで設定すればよい。また、回転軸を設定する際に、対向する回転軸を設定した場合には、この２つの軸に対する回転は、回転方向が逆になるだけなので、対向する回転軸を１つにする、または、回転角度を一方向にする、ことができる。このようなことを考慮して、上述した実施の形態では、回転軸Ｊ１〜Ｊ１２を中心とする角度範囲を２〜３２゜とし、２゜刻みで回転させるようにしている。
【００８１】
〔処理速度について〕
初期段階で、１２個の回転軸（Ｋ＝１２）、１６の回転角度（Ｌ＝１６）を選択した場合、１９２回の３次元回転と３次元相関処理が必要となる。これには膨大な演算を必要とするが、最初に求めるのは概略の回転軸と回転角度であるため、ボクセルデータ群を小さくしても問題はない。ボクセルデータ群を小さくすることで、演算量は少なくなり、高速処理が期待できる。このため、上述した実施の形態では、５１２×５１２×５１２の原画像のデータサイズを縮小し、１６×１６×１６の縮小したデータサイズとしている。
【００８２】
データサイズの一辺の数を１／２にした場合（ボクセル数では１／８）、通常のパーソナルコンピュータ（ＰＣ）にて３次元位相限定相関法を実行すると、処理時間は１／５から１／８程度になる。したがって、５１２×５１２×５１２のデータサイズのボクセルデータ群に対して３次元位相限定相関法を実行した場合の時間を１とすると、１６×１６×１６のデータサイズのボクセルデータ群では、１／３１２５から１／３２７６８となり、多数回実行しても問題ない処理スピードとなる。
【００８３】
〔合成ベクトルを作成する際の重み付けについて〕
合成ベクトルを作成する際に用いる重み付けは、相関処理の出力する値（相関度合いが高くなるほど高い値となる）をそのまま使用しても構わないが、相関処理の特性に合わせて任意の処理を行ってから使用することもある。例えば、相関処理の特性が緩慢な場合（完全に一致しても、少しずれていても、出力値が大きく変わらない場合）、出力値に相関関数を適用することもできるし、逆に急峻な場合には、平方根を適用することもできる。
【００８４】
〔機能ブロック図〕
図１５、図１６および図１７に上述した３次元ずれ量計測装置４の要部の機能ブロック図を分割して示す。３次元ずれ量計測装置４は、基準立体画像の原画像（データサイズ：５１２×５１２×５１２）を記憶するデータベースＤＢ１と、比較対象立体画像の原画像（データサイズ：５１２×５１２×５１２）を記憶するデータベースＤＢ２とを備えている。
【００８５】
また、第１ステージ用として、第１の基準立体画像データサイズ縮小部４−１と、第１の比較対象立体画像データサイズ縮小部４−２と、第１の比較候補回転軸設定部４−３と、第１の比較候補立体画像作成部４−４と、第１の相関値算出部４−５と、第１の最大相関値抽出部４−６と、第１ステージ比較確定回転角度決定部４−７と、相関値抽出部４−８と、第１ステージ比較回転軸決定部４−９とを備えている。
【００８６】
また、第２ステージ用として、第２の比較候補立体画像作成部４−１０と、第２の相関値算出部４−１１と、第２の最大相関値抽出部４−１２と、第２ステージ比較確定回転角度決定部４−１３と、第２の比較候補回転軸設定部４−１４と、第３の比較候補立体画像作成部４−１５と、第３の相関値算出部４−１６と、第３の最大相関値抽出部４−１７と、第２ステージ比較回転軸決定部４−１８とを備えている。なお、第１ステージ用として設けた第１の基準立体画像データサイズ縮小部４−１と第１の比較対象立体画像データサイズ縮小部４−２は、第２ステージ用としても用いられる。
【００８７】
また、第３ステージ用として、第２の基準立体画像データサイズ縮小部４−１９と、第２の比較対象立体画像データサイズ縮小部４−２０と、第４の比較候補立体画像作成部４−２１と、第４の相関値算出部４−２２と、第４の最大相関値抽出部４−２３と、第３ステージ比較確定回転角度決定部４−２４と、第３の比較候補回転軸設定部４−２５と、第５の比較候補立体画像作成部４−２６と、第５の相関値算出部４−２７と、第５の最大相関値抽出部４−２６と、第３ステージ比較回転軸決定部４−２９とを備えている。
【００８８】
第１の基準立体画像データサイズ縮小部４−１は、データベースＤＢ１に保存されている５１２×５１２×５１２の基準立体画像（原画像）ＧＲのデータサイズを縮小し、第１ステージおよび第２ステージで用いる基準立体画像（データサイズ：１６×１６×１６）ＧＲ１とする。第１の比較対象立体画像データサイズ縮小部４−２は、データベースＤＢ２に保存されている５１２×５１２×５１２の比較対象立体画像（原画像）ＧＣのデータサイズを縮小し、第１ステージおよび第２ステージで用いる比較対象立体画像（データサイズ：１６×１６×１６）ＧＣ１とする。
【００８９】
第１の比較回転軸設定部４−３は、第１ステージの比較対象立体画像ＧＣ１のＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ個の回転軸を比較候補回転軸Ｊ１〜ＪＫとして設定する。第１の比較候補立体画像作成部４−４は、第１の比較回転軸設定部４−３によって作成されたＫ個の比較候補回転軸を中心として第１ステージの比較対象立体画像ＧＣ１を所定の角度範囲内でＬ通り回転させたＫ×Ｌ個の画像を比較候補立体画像として作成する。第１の相関値算出部４−５は、第１の比較候補立体画像作成部４−４によって作成されたＫ×Ｌ個の比較候補立体画像と第１ステージの基準立体画像ＧＲ１との相関値を各個に算出する。
【００９０】
第１の最大相関値抽出部４−６は、第１の相関値算出部４−５によって算出されたＫ×Ｌ個の相関値の中から最大相関値を抽出する。第１ステージ比較確定回転角度決定部４−７は、最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして決定する。
【００９１】
相関値抽出部４−８は、第１の相関値算出部４−５によって算出されたＫ×Ｌ個の相関値の中から、第１ステージでの比較確定回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する。
【００９２】
第１ステージ比較確定回転軸決定部４−９は、相関値抽出部４−８によって抽出された相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル（原点からの単位長さ）に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像ＧＣの第１ステージでの比較確定回転軸ＪＡとする。
【００９３】
第２の比較候補立体画像作成部４−１０は、第１ステージでの比較確定回転軸ＪＡを中心とし、第１ステージでの比較確定回転角度θＡを含む所定の角度範囲内で比較対象立体画像ＧＣ１をＭ１通り回転させたＭ１個の画像を比較候補立体画像として作成する。第２の相関値算出部４−１１は、第２の比較候補立体画像作成部４−１０によって作成されたＭ１個の比較候補立体画像と基準立体画像ＧＲ１との相関値を各個に算出する。
【００９４】
第２の最大相関値抽出部４−１２は、第２の相関値算出部４−１１によって算出されたＭ１個の相関値の中から最大相関値を抽出する。第２ステージ比較確定回転角度決定部４−１３は、最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第２ステージでの比較確定回転角度θＢとして決定する。
【００９５】
第２の比較候補回転軸設定部４−１４は、第１ステージでの比較確定回転軸ＪＡを中心にＮ１個の比較候補回転軸を設定する。第３の比較候補立体画像作成部４−１５は、第２の比較候補回転軸設定部４−１４によって設定されたＮ１個の比較候補回転軸のそれぞれについて、その回転軸を中心として比較確定回転角度θＢ回転させたＮ１個の比較候補立体画像を作成する。第３の相関値算出部４−１６は、第３の比較候補立体画像作成部４−１５によって作成されたＮ１個の比較候補立体画像と基準立体画像ＧＲ１との相関値を各個に算出する。
【００９６】
第３の最大相関値抽出部４−１７は、第３の相関値算出部４−１６によって算出されたＮ１個の相関値の中から最大相関値を抽出する。第２ステージ比較確定回軸決定部４−１８は、最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとして決定する。
【００９７】
第２の基準立体画像データサイズ縮小部４−１９は、データベースＤＢ１に保存されている５１２×５１２×５１２の基準立体画像（原画像）ＧＲのデータサイズを縮小し、第３ステージで用いる基準立体画像（データサイズ：６４×６４×６４）ＧＲ２とする。第１の比較対象立体画像データサイズ縮小部４−２０は、データベースＤＢ２に保存されている５１２×５１２×５１２の比較対象立体画像（原画像）ＧＣのデータサイズを縮小し、第３ステージで用いる比較対象立体画像（データサイズ：６４×６４×６４）ＧＣ２とする。
【００９８】
第４の比較候補立体画像作成部４−２１は、第２ステージでの比較確定回転軸ＪＢを中心とし、第２ステージでの比較確定回転角度θＢを含む所定の角度範囲内で比較対象立体画像ＧＣ２をＭ２通り回転させたＭ２個の画像を比較候補立体画像として作成する。第４の相関値算出部４−２２は、第４の比較候補立体画像作成部４−２１によって作成されたＭ２個の比較候補立体画像と基準立体画像ＧＲ２との相関値を各個に算出する。
【００９９】
第４の最大相関値抽出部４−２３は、第４の相関値算出部４−２２によって算出されたＭ２個の相関値の中から最大相関値を抽出する。第３ステージ比較確定回転角度決定部４−２４は、最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第３ステージでの比較確定回転角度θＣとして決定する。
【０１００】
第３の比較候補回転軸設定部４−２５は、第２ステージでの比較確定回転軸ＪＢを中心にＮ２個の比較候補回転軸を設定する。第５の比較候補立体画像作成部４−２６は、第３の比較候補回転軸設定部４−２５によって設定されたＮ２個の比較候補回転軸のそれぞれについてその回転軸を中心として回転角度θＣ回転させたＮ２個の比較候補立体画像を作成する。、第５の相関値算出部４−２７は、第５の比較候補立体画像作成部４−２６によって作成されたＮ２個の比較候補立体画像と基準立体画像ＧＲ２との相関値を各個に算出する。
【０１０１】
第５の最大相関値抽出部４−２８は、第５の相関値算出部４−２７によって算出されたＮ２個の相関値の中から最大相関値を抽出する。第３ステージ比較確定回軸決定部４−２９は、最大相関値が得られた比較候補立体画像の比較候補回転軸を第３ステージでの比較確定回転軸ＪＣとして決定する。
【０１０２】
なお、上述した実施の形態では、最初に撮影した被写体Ａの３次元の立体画像を基準立体画像とし、次に撮影した被写体Ｂの３次元の立体画像を比較対象立体画像としたが、最初に撮影した被写体Ａの３次元の立体画像を比較対象立体画像とし、次に撮影した被写体Ｂの３次元の立体画像を基準立体画像としてもよい。また、基準立体画像や比較対象立体画像は、ボクセルデータ群として上位装置から与えられるものとしてもよい。
【０１０３】
また、上述した実施の形態では、被写体を人体としたが、被写体は人体に限られるものではない。例えば、製薬などの分野において、マウスを使用した実験などにも利用することができる。また、産業分野において、金属などの非破壊検査に利用することもできる。このように、本発明は、医療や歯科などの分野に限らず、各種の分野において有効利用することが可能である。
【０１０４】
また、上述した３次元ずれ量計測処理装置４において、図１５〜図１７に分割して機能ブロックとして示した各部の処理動作は、プログラム（３次元ずれ量計測プログラム）に従うＣＰＵの処理動作として実現されるが、同様の処理をハードウェアの回路構成で実現するようにしてもよい。
【０１０５】
また、上述した実施の形態では、基準立体画像と比較対象立体画像との間のずれ量の計測処理を第３ステージを最終ステージとして実行するようにしたが、第１ステージで終了するようにしてもよく、第２ステージで終了するようにしてもよい。また、第１ステージや第２ステージで終了するような場合、基準立体画像および比較対象立体画像の原画像のデータサイズは、必ずしも縮小しなくてもよい。
【図面の簡単な説明】
【０１０６】
【図１】本発明に係る３次元ずれ量計測方法の実施に用いるＸ線コンピュータ断層撮影システムの一例を示す図である。
【図２】このシステムの３次元ずれ量計測装置におけるＣＰＵの処理動作（基準立体画像の取得処理）を示すフローチャートである。
【図３】このシステムの３次元ずれ量計測装置におけるＣＰＵの処理動作（比較対象立体画像の取得処理）を示すフローチャートである。
【図４】このシステムの３次元ずれ量計測装置におけるＣＰＵの処理動作（基準立体画像と比較対象立体画像との間のずれ量の計測処理）を示すフローチャートである。
【図５】図４に続くフローチャートである。
【図６】図５に続くフローチャートである。
【図７】このシステムの３次元ずれ量計測装置におけるＣＰＵの処理動作（位置合わせ処理＋照合処理）を示すフローチャートである。
【図８】正二十面体の中心と各頂点とを結ぶ軸を比較候補回転軸として設定する状況を示す図である。
【図９】各比較候補回転軸を中心として所定の角度範囲内で回転させた比較候補立体画像と第１ステージの基準立体画像との相関値との関係を示す図である。
【図１０】各比較候補回転軸の単位ベクトル（原点からの単位長さ）に相関値に基づく重み付けを施して得られる合成ベクトルＶＡおよひ第１ステージの比較確定回転軸ＪＡを説明する図である。
【図１１】第１ステージでの比較確定回転軸ＪＡを中心とし第１ステージでの比較確定回転角度θＡを中心とする所定の角度範囲内（±１０゜）で比較対象立体画像をＭ１通り回転させる状況を説明する図である。
【図１２】第１ステージでの比較確定回転軸ＪＡを中心軸としその周囲に設定されるＮ１個の比較候補回転軸を説明する図である。
【図１３】第２ステージでの比較確定回転軸ＪＢを中心とし第２ステージでの比較確定回転角度θＢを中心とする所定の角度範囲内（±１゜）で比較対象立体画像をＭ２通り回転させる状況を説明する図である。
【図１４】第２ステージでの比較確定回転軸ＪＡを中心軸としそのその周囲に設定されるＮ２個の比較候補回転軸を説明する図である。
【図１５】このシステムにおける３次元ずれ量計測装置の要部の機能ブロック図である。
【図１６】図１５に続く３次元ずれ量計測装置の要部の機能ブロック図である。
【図１７】図１６に続く３次元ずれ量計測装置の要部の機能ブロック図である。
【図１８】Ｘ線コンピュータ断層撮影法による３次元の立体画像の取得例を説明する図である。
【符号の説明】
【０１０７】
１…Ｘ線源、２…２次元撮像装置、３（Ａ，Ｂ）…被写体、４Ａ…ＣＰＵ、４Ｂ…ＲＡＭ、４Ｃ…ＲＯＭ、４Ｄ…記憶装置、４Ｅ〜Ｈ…インタフェース、４Ｉ…ディスプレイ、４Ｊ…キーボード、４Ｋ…マウス、ＤＢ１，ＤＢ２…データベース、４−１…第１の基準立体画像データサイズ縮小部、４−２…第１の比較対象立体画像データサイズ縮小部、４−３…第１の比較候補回転軸設定部、４−４…第１の比較候補立体画像作成部、４−５…第１の相関値算出部、４−６…第１の最大相関値抽出部、４−７…第１ステージ比較確定回転角度決定部、４−８…相関値抽出部、４−９…第１ステージ比較回転軸決定部、４−１０…第２の比較候補立体画像作成部、４−１１…第２の相関値算出部、４−１２…第２の最大相関値抽出部、４−１３…第２ステージ比較確定回転角度決定部、４−１４…第２の比較候補回転軸設定部、４−１５…第３の比較候補立体画像作成部、４−１６…第３の相関値算出部、４−１７…第３の最大相関値抽出部、４−１８…第２ステージ比較回転軸決定部、４−１９…第２の基準立体画像データサイズ縮小部、４−２０…第２の比較対象立体画像データサイズ縮小部、４−２１…第４の比較候補立体画像作成部、４−２２…第４の相関値算出部、４−２３…第４の最大相関値抽出部、４−２４…第３ステージ比較確定回転角度決定部、４−２５…第３の比較候補回転軸設定部、４−２６…第５の比較候補立体画像作成部、４−２７…第５の相関値算出部、４−２６…第５の最大相関値抽出部、４−２９…第３ステージ比較回転軸決定部。
【特許請求の範囲】
【請求項１】
基準となる３次元の立体画像を基準立体画像、この基準立体画像に対して比較対象とする３次元の立体画像を比較対象立体画像とし、前記基準立体画像と前記比較対象立体画像との間のずれ量を計測する３次元ずれ量計測方法において、
前記比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ（Ｋ≧２）個の回転軸を比較候補回転軸として設定する第１ステップと、
前記Ｋ個の各比較候補回転軸を中心として前記比較対象立体画像を所定角度範囲内でＬ（Ｌ≧２）通り回転させたＫ×Ｌ個の画像を比較候補立体画像とし、このＫ×Ｌ個の比較候補立体画像と前記基準立体画像との相関値を各個に算出する第２ステップと、
前記算出されたＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして求める第３ステップと、
前記算出されたＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から前記比較確定回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する第４ステップと、
前記比較候補回転軸の前記原点からの単位長さを単位ベクトルとし、前記第４ステップで抽出された相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトルに当該相関値に基づく重み付けを施し、この重み付けが施された単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を前記比較対象立体画像の第１ステージでの比較確定回転軸ＪＡとする第５ステップと
を備えることを特徴とする３次元ずれ量計測方法。
【請求項２】
請求項１に記載された３次元ずれ量計測方法において、
前記第１ステージでの比較確定回転軸ＪＡを中心とし、前記第１ステージでの比較確定回転角度θＡを含む所定角度範囲内で前記比較対象立体画像をＭ（Ｍ≧２）通り回転させたＭ個の画像を比較候補立体画像とし、このＭ個の比較候補立体画像と前記基準立体画像との相関値を各個に算出する第６ステップと、
前記算出されたＭ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の前記第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとして求める第７ステップと、
前記第１ステージでの比較確定回転軸ＪＡを中心に任意のＮ（Ｎ≧２）個の比較候補回転軸を設定し、このＮ個の比較候補回転軸のそれぞれについてその回転軸を中心として前記比較確定回転角度θＢ回転させたＮ個の比較候補立体画像と前記基準立体画像との相関値を各個に算出する第８ステップと、
前記Ｎ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとする第９ステップと
を備えることを特徴とする３次元ずれ量計測方法。
【請求項３】
請求項２に記載された３次元ずれ量計測方法において、
前記第７ステップに代えて、
前記算出されたＭ個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の前記第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとして求める第７ステップ
を備えることを特徴とする３次元ずれ量計測方法。
【請求項４】
請求項２に記載された３次元ずれ量計測方法において、
前記第９ステップに代えて、
前記Ｎ個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとする第９ステップ
を備えることを特徴とする３次元ずれ量計測方法。
【請求項５】
請求項２〜４の何れか１項に記載された３次元ずれ量計測方法において、
前記第１ステップの処理に入る前に、前記基準立体画像および比較対象立体画像の原画像のデータサイズを小さくして前記第１ステージで用いる基準立体画像および比較対象立体画像とするステップと、
前記第１ステージで用いる基準立体画像および比較対象立体画像よりもそのデータサイズを大きくしながら前記第６〜第９ステップの処理を予め定められた最終ステージまで繰り返すステップと
を備えることを特徴とする３次元ずれ量計測方法。
【請求項６】
請求項５に記載された３次元ずれ量計測方法において、
前記最終ステージでの比較確定回転角度および比較確定回転軸に基づいて前記基準立体画像の原画像と前記比較対象立体画像の原画像との間の位置合わせを行う位置合わせ処理ステップ
を備えることを特徴とする３次元ずれ量計測方法。
【請求項７】
請求項６に記載された３次元ずれ量計測方法において、
前記位置合わせ処理ステップによって位置合わせされた前記基準立体画像の原画像と前記比較対象立体画像の原画像とを相関演算によって照合する照合ステップ
を備えることを特徴とする３次元ずれ量計測方法。
【請求項８】
請求項６又は７に記載された３次元ずれ量計測方法において、
前記位置合わせ処理ステップは、
前記最終ステージでの比較確定回転角度および比較確定回転軸に基づいて前記基準立体画像の原画像と前記比較対象立体画像の原画像との間の回転ずれを補正し、この回転ずれが補正された基準立体画像の原画像と比較対象立体画像の原画像との間で相関演算を行って当該基準立体画像の原画像と比較対象立体画像の原画像との間の前記Ｘ，Ｙ，Ｚの３軸方向のずれ量を求め、この求めた３軸方向のずれ量に基づいて前記回転ずれが補正された基準立体画像の原画像と比較対象立体画像の原画像との間の平行ずれを補正する
ことを特徴とする３次元ずれ量計測方法。
【請求項９】
請求項１〜８の何れか１項に記載された３次元ずれ量計測方法において、
前記第１ステップの処理に入る前に、前記基準立体画像の原画像および前記比較対象立体画像の原画像を構成するボクセルデータ群に対してそれぞれ３次元フーリエ変換を施し、この３次元フーリエ変換が施された基準立体画像の原画像および比較対象立体画像の原画像を構成するボクセルデータ群の３次元周波数空間の中心を前記基準立体画像および比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心とするステップ
を備えることを特徴とする３次元ずれ量計測方法。
【請求項１】
基準となる３次元の立体画像を基準立体画像、この基準立体画像に対して比較対象とする３次元の立体画像を比較対象立体画像とし、前記基準立体画像と前記比較対象立体画像との間のずれ量を計測する３次元ずれ量計測方法において、
前記比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心を原点とし、この原点を中心として任意の方向へＫ（Ｋ≧２）個の回転軸を比較候補回転軸として設定する第１ステップと、
前記Ｋ個の各比較候補回転軸を中心として前記比較対象立体画像を所定角度範囲内でＬ（Ｌ≧２）通り回転させたＫ×Ｌ個の画像を比較候補立体画像とし、このＫ×Ｌ個の比較候補立体画像と前記基準立体画像との相関値を各個に算出する第２ステップと、
前記算出されたＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第１ステージでの比較確定回転角度θＡとして求める第３ステップと、
前記算出されたＫ×Ｌ個の比較候補立体画像と基準立体画像との相関値の中から前記比較確定回転角度θＡと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する第４ステップと、
前記比較候補回転軸の前記原点からの単位長さを単位ベクトルとし、前記第４ステップで抽出された相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトルに当該相関値に基づく重み付けを施し、この重み付けが施された単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を前記比較対象立体画像の第１ステージでの比較確定回転軸ＪＡとする第５ステップと
を備えることを特徴とする３次元ずれ量計測方法。
【請求項２】
請求項１に記載された３次元ずれ量計測方法において、
前記第１ステージでの比較確定回転軸ＪＡを中心とし、前記第１ステージでの比較確定回転角度θＡを含む所定角度範囲内で前記比較対象立体画像をＭ（Ｍ≧２）通り回転させたＭ個の画像を比較候補立体画像とし、このＭ個の比較候補立体画像と前記基準立体画像との相関値を各個に算出する第６ステップと、
前記算出されたＭ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の前記第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとして求める第７ステップと、
前記第１ステージでの比較確定回転軸ＪＡを中心に任意のＮ（Ｎ≧２）個の比較候補回転軸を設定し、このＮ個の比較候補回転軸のそれぞれについてその回転軸を中心として前記比較確定回転角度θＢ回転させたＮ個の比較候補立体画像と前記基準立体画像との相関値を各個に算出する第８ステップと、
前記Ｎ個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとする第９ステップと
を備えることを特徴とする３次元ずれ量計測方法。
【請求項３】
請求項２に記載された３次元ずれ量計測方法において、
前記第７ステップに代えて、
前記算出されたＭ個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の前記第１ステージでの比較確定回転軸ＪＡを中心とする回転角度を第２ステージでの比較確定回転角度θＢとして求める第７ステップ
を備えることを特徴とする３次元ずれ量計測方法。
【請求項４】
請求項２に記載された３次元ずれ量計測方法において、
前記第９ステップに代えて、
前記Ｎ個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第２ステージでの比較確定回転軸ＪＢとする第９ステップ
を備えることを特徴とする３次元ずれ量計測方法。
【請求項５】
請求項２〜４の何れか１項に記載された３次元ずれ量計測方法において、
前記第１ステップの処理に入る前に、前記基準立体画像および比較対象立体画像の原画像のデータサイズを小さくして前記第１ステージで用いる基準立体画像および比較対象立体画像とするステップと、
前記第１ステージで用いる基準立体画像および比較対象立体画像よりもそのデータサイズを大きくしながら前記第６〜第９ステップの処理を予め定められた最終ステージまで繰り返すステップと
を備えることを特徴とする３次元ずれ量計測方法。
【請求項６】
請求項５に記載された３次元ずれ量計測方法において、
前記最終ステージでの比較確定回転角度および比較確定回転軸に基づいて前記基準立体画像の原画像と前記比較対象立体画像の原画像との間の位置合わせを行う位置合わせ処理ステップ
を備えることを特徴とする３次元ずれ量計測方法。
【請求項７】
請求項６に記載された３次元ずれ量計測方法において、
前記位置合わせ処理ステップによって位置合わせされた前記基準立体画像の原画像と前記比較対象立体画像の原画像とを相関演算によって照合する照合ステップ
を備えることを特徴とする３次元ずれ量計測方法。
【請求項８】
請求項６又は７に記載された３次元ずれ量計測方法において、
前記位置合わせ処理ステップは、
前記最終ステージでの比較確定回転角度および比較確定回転軸に基づいて前記基準立体画像の原画像と前記比較対象立体画像の原画像との間の回転ずれを補正し、この回転ずれが補正された基準立体画像の原画像と比較対象立体画像の原画像との間で相関演算を行って当該基準立体画像の原画像と比較対象立体画像の原画像との間の前記Ｘ，Ｙ，Ｚの３軸方向のずれ量を求め、この求めた３軸方向のずれ量に基づいて前記回転ずれが補正された基準立体画像の原画像と比較対象立体画像の原画像との間の平行ずれを補正する
ことを特徴とする３次元ずれ量計測方法。
【請求項９】
請求項１〜８の何れか１項に記載された３次元ずれ量計測方法において、
前記第１ステップの処理に入る前に、前記基準立体画像の原画像および前記比較対象立体画像の原画像を構成するボクセルデータ群に対してそれぞれ３次元フーリエ変換を施し、この３次元フーリエ変換が施された基準立体画像の原画像および比較対象立体画像の原画像を構成するボクセルデータ群の３次元周波数空間の中心を前記基準立体画像および比較対象立体画像に設定される互いに直交するＸ，Ｙ，Ｚの３軸の中心とするステップ
を備えることを特徴とする３次元ずれ量計測方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【公開番号】特開２０１０−１２５１１７（Ｐ２０１０−１２５１１７Ａ）
【公開日】平成２２年６月１０日（２０１０．６．１０）
【国際特許分類】
【出願番号】特願２００８−３０３８６７（Ｐ２００８−３０３８６７）
【出願日】平成２０年１１月２８日（２００８．１１．２８）
【出願人】（０００００６６６６）株式会社山武 (1,808)
【出願人】（５０４１５７０２４）国立大学法人東北大学 (2,297)
【Ｆターム（参考）】
【公開日】平成２２年６月１０日（２０１０．６．１０）
【国際特許分類】
【出願日】平成２０年１１月２８日（２００８．１１．２８）
【出願人】（０００００６６６６）株式会社山武 (1,808)
【出願人】（５０４１５７０２４）国立大学法人東北大学 (2,297)
【Ｆターム（参考）】
[ Back to top ]