JP5171575B2 - 3次元ずれ量計測方法 - Google Patents
3次元ずれ量計測方法 Download PDFInfo
- Publication number
- JP5171575B2 JP5171575B2 JP2008303867A JP2008303867A JP5171575B2 JP 5171575 B2 JP5171575 B2 JP 5171575B2 JP 2008303867 A JP2008303867 A JP 2008303867A JP 2008303867 A JP2008303867 A JP 2008303867A JP 5171575 B2 JP5171575 B2 JP 5171575B2
- Authority
- JP
- Japan
- Prior art keywords
- comparison
- stereoscopic image
- stage
- candidate
- image
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000691 measurement method Methods 0.000 title claims 3
- 238000006073 displacement reaction Methods 0.000 title description 7
- 238000000034 method Methods 0.000 claims description 54
- 238000012545 processing Methods 0.000 claims description 44
- 239000013598 vector Substances 0.000 claims description 29
- 238000012790 confirmation Methods 0.000 claims description 13
- 230000000052 comparative effect Effects 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 description 33
- 238000000605 extraction Methods 0.000 description 19
- 238000005549 size reduction Methods 0.000 description 14
- 238000002591 computed tomography Methods 0.000 description 9
- 238000003384 imaging method Methods 0.000 description 9
- 230000029777 axis specification Effects 0.000 description 8
- 238000010586 diagram Methods 0.000 description 8
- 239000000284 extract Substances 0.000 description 6
- 238000005259 measurement Methods 0.000 description 6
- 238000003325 tomography Methods 0.000 description 3
- 241000699666 Mus <mouse, genus> Species 0.000 description 2
- 239000002131 composite material Substances 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 239000003814 drug Substances 0.000 description 2
- 241000699670 Mus sp. Species 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 238000005314 correlation function Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 229910052751 metal Inorganic materials 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 150000002739 metals Chemical class 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 238000002600 positron emission tomography Methods 0.000 description 1
- 238000002603 single-photon emission computed tomography Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Landscapes
- Apparatus For Radiation Diagnosis (AREA)
- Image Processing (AREA)
- Image Generation (AREA)
- Image Analysis (AREA)
Description
この発明は、基準となる3次元の立体画像を基準立体画像、この基準立体画像に対して比較対象とする3次元の立体画像を比較対象立体画像とし、基準立体画像と比較対象立体画像との間のずれ量を計測する3次元ずれ量計測方法に関するものである。
従来より、医学分野などにおいては、多様な画像化方式によって立体画像を得ている。例えば、磁気共鳴断層撮影法(MRI)、X線コンピュータ断層撮影法(CT)、陽電子放射断層撮影法(PET)、単光子射出断層撮影法(SPECT)、超音波断層撮影法などによって、3次元の立体画像を得ている。
このような方法で得られた3次元の立体画像は、ボクセル(voxel)値の配列として、通常、コンピュータのメモリなどにデジタル的に保存される。各ボクセルは、互いに直交するX,Y,Zの3軸で表される3次元空間の座標位置(x,y,z)と結びついており、ボクセル値として色値(典型的にはグレイスケール)が割り当てられる。
図18にX線コンピュータ断層撮影法による3次元の立体画像の取得例を示す。この例では、スキャナ装置としてコーンビームCTを用いている。コーンビームCTは、中間点C0を挾んで対向してXY平面内に設けられたX線源1と2次元撮像装置2とを備え、このX線源1と2次元撮像装置2とを中間点C0を通るZ軸を中心として水平方向に定速度で回転させながら、被写体(人体)3の撮影を行う。この撮影に際して、被写体3は、その頭部の中心を中間点C0に位置させる。
被写体3の撮影は、X線源1と2次元撮像装置2とが1回転する間、所定の時間間隔で行われる。この結果、1回転する間に任意の異なる角度(以下、この角度を投影角度と呼ぶ)から撮影されたn個(例えば、256個)の2次元の投影画像を取得して、このn個の投影画像を再構成処理することで3次元の立体画像を得ることができる。この3次元の立体画像は、互いに直交するX,Y,Zの3軸で表される3次元空間の座標位置(x,y,z)と結びつけたボクセルデータ群として表現される。
この3次元の立体画像は、医療や歯科などの分野において、各種の診断を行うために使用される。しかし、情報量が多いために見落としが発生する虞があり、正確な診断を行うために過大な時間を必要とする問題がある。このため、コンピュータ上で、画像処理を行うことで、この問題を解決しようと試みられている。
この場合の画像処理の1つに、レジストレーションと呼ばれる位置合わせを行った後、2つの立体画像を比較するという方法がある。例えば、前回撮影した被写体3(被写体A)の3次元の立体画像を基準立体画像とし、今回撮影した被写体3(被写体B)の3次元の立体画像を比較対象立体画像とし、基準立体画像と比較対象立体画像との位置合わせを行い、この位置合わせが行われた2つの立体画像の差分画像をディスプレイ上に表示する。これにより、前回撮影時と今撮影時との間の変化部分のみを抽出し、短時間でかつ正確な診断を行うことができるようになる。
この方法において、基準立体画像と比較対象立体画像とのレジストレーションは、6つのパラメータを求めることで実現される。6つのパラメータは、X,Y,Zの各軸の平行ずれ量Δx,Δy,Δzと、X,Y,Zの各軸の回転ずれ量Δθ,Δφ,ΔΨである。平行ずれ量Δx,Δy,Δzについては、回転ずれ量Δθ,Δφ,ΔΨを求めた後に、3次元相関処理を行うことで、容易に求めることができる。回転ずれ量Δθ,Δφ,ΔΨは、X,Y,Zの各軸の回転角度を3つの自由度とし、この3つの自由度を変化させながら3次元相関処理を行うことで求めることが可能である(例えば、特許文献1参照)。
しかしながら、上述した基準立体画像と比較対象立体画像とのレジストレーションでは、回転ずれ量Δθ,Δφ,ΔΨを求める際、X,Y,Zの各軸の回転角度を自由度として変化させるため、3乗の演算回数が必要となる。例えば、3つの自由度それぞれに10種類の回転角度を設定した場合、10×10×10=1000回の3次元の相関処理が必要となる。このため、膨大な演算を必要とし、高速処理が困難となる。
なお、相関処理を用いないレジストレーションとして、立体画像に現れている被写体の構造に詳しい専門家が、位置合わせする立体画像のそれぞれに目印をつけ、一方の立体画像を移動回転させることによって位置合わせする方法がある(例えば、特許文献2参照)。
しかしながら、この方法では、位置合わせ精度が選択された目印の数と位置によって制限されてしまう。また、目印の数が少な過ぎると、位置合わせが不正確になる。また、目印が多くても必ずしも位置合わせが正確となるわけではなく、位置合わせに要する計算を複雑にするだけである。また、オペレータが手作業で対応点を選択する必要があり、時間がかかるとともに熟練が必要となる。さらに、すべての立体画像で適当な構造上の目印を見出せるわけでもない。
本発明は、このような課題を解決するためになされたもので、その目的とするところは、少ない演算量で、高速かつ正確に、レジストレーションを行うことが可能な3次元ずれ量計測方法を提供することにある。
このような目的を達成するために本発明は、基準となる3次元の立体画像を基準立体画像、この基準立体画像に対して比較対象とする3次元の立体画像を比較対象立体画像とし、基準立体画像と比較対象立体画像との間のずれ量を計測する3次元ずれ量計測方法において、比較対象立体画像に設定される互いに直交するX,Y,Zの3軸の中心を原点とし、この原点を中心として任意の方向へK(K≧2)個の回転軸を比較候補回転軸として設定し(第1ステップ)、K個の各比較候補回転軸を中心として比較対象立体画像を所定角度範囲内でL(L≧2)通り回転させたK×L個の画像を比較候補立体画像とし、このK×L個の比較候補立体画像と基準立体画像との相関値を各個に算出し(第2ステップ)、この算出されたK×L個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第1ステージでの比較確定回転角度θAとして求め(第3ステップ)、算出されたK×L個の比較候補立体画像と基準立体画像との相関値の中から比較確定回転角度θAと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出し(第4ステップ)、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル(原点からの単位長さ)に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像の第1ステージでの比較確定回転軸JAとする(第5ステップ)、ようにしたものである。
本発明を適用しようとする3次元空間において、その3次元空間における回転角度のずれには、様々な表現方法が考えられる。最も簡単に表現できるのは、互いに直交するX,Y,Zの3軸のそれぞれを回転軸とした3つの回転角度により表現することである。すなわち、回転軸を3つとし、これらの回転軸を中心としてそれぞれある角度にて回転したと考え、合計3回の回転が発生していると考えることである。本発明では、この考えをさらに進め、原点を中心として複数回の回転を行うことは、ある回転軸を中心として1回の回転を行ったに過ぎないと考え、回転軸と回転角度の2つのパラメータを求めることにする。
この考えに基づき、本発明では、比較対象立体画像に設定される互いに直交するX,Y,Zの3軸の中心を原点とし、この原点を中心として任意の方向へK(K≧2)個の回転軸を比較候補回転軸として設定する(第1ステップ)。例えば、原点を中心とする正二十面体を考え、この正二十面体の中心と各頂点とを結ぶ12個の軸を比較候補回転軸として設定する。そして、このK個の各比較候補回転軸を中心として比較対象立体画像を所定角度範囲内でL(L≧2)通り回転させたK×L個の画像(例えば、2〜32゜の角度範囲内で2゜刻みで回転させた12×16=192個の画像)を比較候補立体画像とし、このK×L個の比較候補立体画像と基準立体画像との相関値を各個に算出する(第2ステップ)。
そして、この算出したK×L個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第1ステージでの比較確定回転角度θAとして求める(第3ステップ)。また、この算出したK×L個の比較候補立体画像と基準立体画像との相関値の中から比較確定回転角度θAと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する(第4ステップ)。そして、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル(原点からの単位長さ)に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像の第1ステージでの比較確定回転軸JAとする(第5ステップ)。
このようにして、本発明では、第1ステージでの比較確定回転軸JAと比較確定回転角度θAを求める。この場合、第1ステージで求められた比較確定回転軸JAと比較確定回転角度θAを用いて基準立体画像と比較対象立体画像との位置合わせ(レジストレーション)を行うようにしてもよいが、レジストレーションの精度をさらに高めるためには、さらに第2ステージに進むことが望ましい。
第2ステージでは、第1ステージでの比較確定回転軸JAを中心とし、第1ステージでの比較確定回転角度θAを含む所定角度範囲内で比較対象立体画像をM(M≧2)通り回転させたM個の画像(例えば、比較確定回転角度θAを中心とする±10゜の角度範囲内で0.5゜刻みで回転させた41個の画像)を比較候補立体画像とし、このM個の比較候補立体画像と基準立体画像との相関値を各個に算出する(第6ステップ)。そして、この算出したM個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の第1ステージでの比較確定回転軸JAを中心とする回転角度を第2ステージでの比較確定回転角度θBとして求める(第7ステップ)。
そして、第1ステージでの比較確定回転軸JAを中心に任意のN個(例えば、11×11=121個)の比較候補回転軸を設定し、このN個の比較候補回転軸のそれぞれについてその回転軸を中心として比較確定回転角度θB回転させたN個の比較候補立体画像と基準立体画像との相関値を各個に算出する(第8ステップ)。そして、このN個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第2ステージでの比較確定回転軸JBとする(第9ステップ)。
なお、上述した第7ステップにおいて、M個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の第1ステージでの比較確定回転軸JAを中心とする回転角度を第2ステージでの比較確定回転角度θBとして求めるようにしてもよい。また、上述した第9ステップにおいて、N個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第2ステージでの比較確定回転軸JBとするようにしてもよい。
また、本発明において、第2ステージに進む場合、第1ステップに入る前に、基準立体画像および比較対象立体画像の原画像のデータサイズを小さくして第1ステージで用いる基準立体画像および比較対象立体画像とし、第1ステージで用いる基準立体画像および比較対象立体画像のデータサイズよりもそのデータサイズを大きくしながら、第6〜第9ステップの処理を予め定められた最終ステージまで繰り返すようにしてもよい。
例えば、基準立体画像および比較対象立体画像の原画像のデータサイズ(画素数)を共に512×512×512とした場合、これら立体画像のデータサイズを縮小して16×16×16の立体画像とする。第1ステージでは、この16×16×16の基準立体画像および比較対象立体画像を用いて第1〜第5ステップの処理を行う。第2ステージでも、第1ステージと同様にして、16×16×16の基準立体画像および比較対象立体画像を用いて第6〜第9ステップの処理を行う。そして、第3ステージから、基準立体画像および比較対象立体画像のデータサイズを大きくし(例えば、64×64×64とする)、第6〜第9ステップの処理を行う。このようにして、最終ステージまで、基準立体画像および比較対象立体画像のデータサイズを大きくしながら、第6〜第9ステップの処理を繰り返す。なお、第2ステージから基準立体画像および比較対象立体画像のデータサイズを大きくするようにしてもよく、第3ステージを最終ステージとしてもよい。
このようにして、本発明では、基準立体画像および比較対象立体画像のデータサイズを最初は小さくし、ステージが進むにつれて大きくして行くことにより、最初は少ないボクセルデータ群で概略的な比較確定回転軸と比較確定回転角度を素速く求め、徐々にボクセルデータ群の数を増やして、より正確な比較確定回転軸と比較確定回転角度を求めて行くようにして、高速度でかつ高精度に、レジストレーションを行うことが可能となる。
本発明では、比較対象立体画像に設定される互いに直交するX,Y,Zの3軸の中心を原点とし、この原点を中心として任意の方向へK個の回転軸を比較候補回転軸として設定し、このK個の比較候補回転軸を中心として回転させた比較対象立体画像と基準立体画像との相関値を求めるが、このような処理において比較対象立体画像と基準立体画像の回転中心とされる原点は同じでなければならない。しかし、3次元の立体画像において、実際の回転中心は不明であるため、原点を決めることが難しい。
本発明において、正確な原点を決めるために、第1ステップの処理に入る前に、基準立体画像の原画像および比較対象立体画像の原画像を構成するボクセルデータ群に対してそれぞれ3次元フーリエ変換を施し、この3次元フーリエ変換が施された基準立体画像の原画像および比較対象立体画像の原画像を構成するボクセルデータ群の3次元周波数空間(ボクセルデータ群の中心を直流とする周波数空間)の中心を基準立体画像および比較対象立体画像に設定される互いに直交するX,Y,Zの3軸の中心とするようにしてもよい。
なお、ボクセルデータ群に対する3次元フーリエ変換の処理は、必ずしも必要というわけではない。例えば、被写体の回転角度が比較的小さいことが分かっている場合、回転中心がほぼ中心付近にあることが分かっている場合などは、3次元フーリエ変換の処理を行う必要はなく、基準立体画像および比較対象立体画像の原画像を構成するボクセルデータ群をそのまま使って、相関処理を行うことが可能である。本発明において、第1ステップの処理に入る前に基準立体画像の原画像および比較対象立体画像の原画像のデータサイズを小さくする方式とした場合、多少の回転中心のずれは無視することができる。
本発明によれば、比較対象立体画像に設定される原点を中心として任意の方向へK個の比較候補回転軸を設定し、このK個の各比較候補回転軸を中心として比較対象立体画像を所定角度範囲内でL通り回転させたK×L個の画像を比較候補立体画像とし、このK×L個の比較候補立体画像と基準立体画像との相関値を各個に算出し、算出したK×L個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第1ステージでの比較確定回転角度θAとして求め、算出したK×L個の比較候補立体画像と基準立体画像との相関値の中から比較確定回転角度θAと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出し、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル(原点からの単位長さ)に当該相関値に基づく重み付けを施し、この重み付けが施された単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像の第1ステージでの比較確定回転軸JAとするようにしたので、この第1ステージで求められた比較確定回転軸JAおよび比較確定回転角度θAを用いて基準立体画像と比較対象立体画像との位置合わせを行うことにより、少ない演算量で、高速かつ正確に、レジストレーションを行うことが可能となる。
また、本発明によれば、第1ステージでの比較確定回転軸JAを中心とし、第1ステージでの比較確定回転角度θAを含む所定角度範囲内で比較対象立体画像をM通り回転させたM個の画像を比較候補立体画像とし、このM個の比較候補立体画像と基準立体画像との相関値を算出し、算出したM個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の第1ステージでの比較確定回転軸JAを中心とする回転角度を第2ステージでの比較確定回転角度θBとして求め、第1ステージでの比較確定回転軸JAを中心に任意のN(N≧2)個の比較候補回転軸を設定し、このN個の比較候補回転軸のそれぞれについてその回転軸を中心として比較確定回転角度θB回転させたN個の比較候補立体画像と基準立体画像との相関値を算出し、N個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第2ステージでの比較確定回転軸JBとするようにしたので、この第2ステージで求められた比較確定回転軸JBおよび比較確定回転角度θBを用いて基準立体画像と比較対象立体画像との位置合わせを行うことにより、さらに正確に、レジストレーションを行うことが可能となる。
また、本発明によれば、第1ステップの処理に入る前に、基準立体画像および比較対象立体画像の原画像のデータサイズを小さくして第1ステージで用いる基準立体画像および比較対象立体画像とし、第1ステージで用いる基準立体画像および比較対象立体画像のデータサイズよりもそのデータサイズを大きくしながら第6〜第9ステップの処理を予め定められた最終ステージまで繰り返すことにより、最初は少ないボクセルデータ群で概略的な比較確定回転軸と比較確定回転角度を素速く求め、徐々にボクセルデータ群の数を増やして、より正確な比較確定回転軸と比較確定回転角度を求めて行くようにして、高速度でかつ高精度に、レジストレーションを行うことが可能となる。
以下、本発明を図面に基づいて詳細に説明する。図1は本発明に係る3次元ずれ量計測方法の実施に用いるX線コンピュータ断層撮影システムの一例を示す図である。同図において、図18と同一符号は図18を参照して説明した構成要素と同一或いは同等構成要素を示し、その説明は省略する。
この実施の形態において、X線源1および2次元撮像装置2は、3次元ずれ量計測装置4に接続されている。3次元ずれ量計測装置4は、CPU4Aと、RAM4Bと、ROM4Cと、ハードディスクなどの記憶装置4Dと、インタフェース4E〜Hと、ディスプレイ4Iと、キーボード4Jと、マウス4Kとを備えている。
CPU4Aは、RAM4Bにアクセスしながら、ROM4Cや記憶装置4Dに格納されたプログラムに従って動作する。記憶装置4Dには、本実施の形態特有のプログラムとして3次元ずれ量計測プログラムが格納されている。この3次元ずれ量計測プログラムは、例えばCD−ROMなどの記録媒体に記録された状態で提供され、この記録媒体から読み出されて記憶装置4Dにインストールされている。
以下、図2〜図7に示すフローチャートを参照して、記憶装置4Dに格納された3次元ずれ量計測プログラムに従ってCPU4Aが実行する本実施の形態特有の処理動作について説明する。
なお、この3次元ずれ量計測プログラムに従う処理動作には、比較対象立体画像の基準立体画像に対する回転軸と回転角度の2つのパラメータ(JEND,ΔθEND)を求めるずれ量計測処理と、この求めた2つのパラメータを用いて基準立体画像と比較対象立体画像との間の位置合わせを行う位置合わせ処理と、この位置合わせされた基準立体画像と比較対象立体画像との照合を行う照合処理とが含まれる。
〔基準立体画像の取得〕
先ず、基準立体画像を取得するために、被写体Aの撮影を行う(図2:ステップS101)。この被写体Aの撮影は、X線源1と2次元撮像装置2とが1回転する間、所定の時間間隔で行う。この結果、1回転する間の任意の異なる角度(投影角度)から撮影されたn個(例えば、256個)の2次元の投影画像が得られる(ステップS102)。なお、この撮影において、被写体Aは、その頭部の中心を中間点CO付近にきちんと位置させるものとする。
先ず、基準立体画像を取得するために、被写体Aの撮影を行う(図2:ステップS101)。この被写体Aの撮影は、X線源1と2次元撮像装置2とが1回転する間、所定の時間間隔で行う。この結果、1回転する間の任意の異なる角度(投影角度)から撮影されたn個(例えば、256個)の2次元の投影画像が得られる(ステップS102)。なお、この撮影において、被写体Aは、その頭部の中心を中間点CO付近にきちんと位置させるものとする。
CPU4Aは、この得られたn個の2次元の投影画像を記憶装置4Dに保存した後(ステップS103)、このn個の投影画像を再構成処理することで3次元の立体画像(ボクセルデータ群)を得る(ステップS104)。そして、この再構成処理によって得られた3次元の立体画像を基準立体画像の原画像として記憶装置4Dに保存する(ステップS105)。なお、この場合、記憶装置4Dに保存される基準立体画像の原画像のデータサイズ(画素数)は、512×512×512とする。
〔比較対象立体画像の取得〕
次に、比較対象立体画像を取得するために、被写体Bの撮影を行う(図3:ステップS201)。この被写体Bの撮影も、X線源1と2次元撮像装置2とが1回転する間、基準立体画像の取得時と同じ所定の時間間隔で行う。この結果、1回転する間の任意の異なる角度(投影角度)から撮影されたn個(例えば、256個)の2次元の投影画像が得られる(ステップS202)。なお、この撮影において、被写体Bは、その頭部の中心を中間点CO付近にきちんと位置させるものとする。
次に、比較対象立体画像を取得するために、被写体Bの撮影を行う(図3:ステップS201)。この被写体Bの撮影も、X線源1と2次元撮像装置2とが1回転する間、基準立体画像の取得時と同じ所定の時間間隔で行う。この結果、1回転する間の任意の異なる角度(投影角度)から撮影されたn個(例えば、256個)の2次元の投影画像が得られる(ステップS202)。なお、この撮影において、被写体Bは、その頭部の中心を中間点CO付近にきちんと位置させるものとする。
CPU4Aは、この得られたn個の2次元の投影画像を記憶装置4Dに保存した後(ステップS203)、このn個の投影画像を再構成処理することで3次元の立体画像(ボクセルデータ群)を得る(ステップS204)。そして、この再構成処理によって得られた3次元の立体画像を比較対象立体画像の原画像として記憶装置4Dに保存する(ステップS205)。なお、この場合、記憶装置4Dに保存される比較対象立体画像の原画像のデータサイズ(画素数)は、512×512×512とする。
〔基準立体画像と比較対象立体画像との間のずれ量の計測〕
〔第1ステージ〕
次に、CPU4Aは、記憶装置4Dに保存されている512×512×512の基準立体画像(原画像)のデータサイズを縮小し、第1ステージで用いる基準立体画像とする(図4:ステップS301)。また、記憶装置4Dに保存されている512×512×512の比較対象立体画像(原画像)のデータサイズを縮小し、第1ステージで用いる比較対象立体画像とする(ステップS302)。
〔第1ステージ〕
次に、CPU4Aは、記憶装置4Dに保存されている512×512×512の基準立体画像(原画像)のデータサイズを縮小し、第1ステージで用いる基準立体画像とする(図4:ステップS301)。また、記憶装置4Dに保存されている512×512×512の比較対象立体画像(原画像)のデータサイズを縮小し、第1ステージで用いる比較対象立体画像とする(ステップS302)。
この例では、X,Y,Zの3軸方向の画素データを所定画素ずつ足し合わせて行くことによって、第1ステージの基準立体画像および比較対象立体画像として16×16×16のデータサイズの立体画像を得る。なお、画素データを間引いて行くことによって、データサイズを縮小するようにしてもよい。
そして、CPU4Aは、ステップ301で得た第1ステージの比較対象立体画像のX,Y,Zの3軸の中心を原点とし、この原点を中心として任意の方向へK個の回転軸を比較候補回転軸として設定する(ステップS303)。この例では、原点を中心とする正二十面体(図8参照)を考え、この正二十面体の中心Oと各頂点とを結ぶ12個の軸J1〜J12を比較候補回転軸として設定する。
そして、このK個の比較候補回転軸を中心として第1ステージの比較対象立体画像を所定の角度範囲内でL通り回転させたK×L個の画像を比較候補立体画像とし(ステップS304)、このK×L個の比較候補立体画像と第1ステージの基準立体画像との相関値を各個に算出する(ステップS305)。
この例では、12個の比較候補回転軸J1〜J12について、その比較候補回転軸を中心として2〜32゜の角度範囲内で2゜刻みで回転させた12×16=192個の画像を比較候補立体画像とし、この192個の比較候補立体画像と第1ステージの基準立体画像との相関値m1〜m192(図9参照)を算出する。
そして、この算出した相関値m1〜m192の中から最大相関値を抽出し(ステップS306)、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第1ステージでの比較確定回転角度θAとして求める(ステップS307)。例えば、相関値m17が最大相関値であった場合、比較候補回転軸J2を中心とする回転角度θ1を第1ステージでの比較確定回転角度θAとして求める。
また、CPU4Aは、ステップS305で算出した相関値m1〜m192の中から、第1ステージでの比較確定回転角度として求めた回転角度θAと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する(ステップS308)。すなわち、比較候補回転軸J1〜J12について、回転角度θAとした場合の比較候補立体画像の相関値を抽出する。
そして、この抽出した相関値の抽出元の比較候補立体画像の比較候補回転軸J1〜J12の単位ベクトル(原点からの単位長さ)に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルVA(図10参照)を求め、この合成ベクトルVAが示す方向を比較対象立体画像の第1ステージでの比較確定回転軸JAとする(ステップS309)。
〔第2ステージ〕
次に、CPU4Aは、記憶装置4Dに保存されている512×512×512の基準立体画像(原画像)のデータサイズを縮小し、第2ステージで用いる基準立体画像とする(図5:ステップS310)。また、記憶装置4Dに保存されている512×512×512の比較対象立体画像(原画像)のデータサイズを縮小し、第2ステージで用いる比較対象立体画像とする(ステップS311)。
次に、CPU4Aは、記憶装置4Dに保存されている512×512×512の基準立体画像(原画像)のデータサイズを縮小し、第2ステージで用いる基準立体画像とする(図5:ステップS310)。また、記憶装置4Dに保存されている512×512×512の比較対象立体画像(原画像)のデータサイズを縮小し、第2ステージで用いる比較対象立体画像とする(ステップS311)。
この例では、ステップ310,311は実行せずに、第1ステージで用いた16×16×16の基準立体画像および比較対象立体画像を第2ステージの基準立体画像および比較対象立体画像として用いることとする。なお、ステップ310,311を実行する場合には、第1ステージで用いた基準立体画像および比較対象立体画像のデータサイズよりも大きいデータサイズに縮小し、第2ステージの基準立体画像および比較対象立体画像とする。
次に、CPU4Aは、第1ステージでの比較確定回転軸JAを中心とし、第1ステージでの比較確定回転角度θAを含む所定の角度範囲内で第2ステージの比較対象立体画像をM1通り回転させたM1個の画像を比較候補立体画像とし、このM1個の比較候補立体画像と第2ステージの基準立体画像との相関値を各個に算出する(ステップS312)。この例では、比較確定回転角度θAを中心とする±10゜の角度範囲内で0.5゜刻みで回転させた41個の画像を比較候補立体画像とし(図11参照)、この41個の比較候補立体画像と第2ステージの基準立体画像との相関値を各個に算出する。
そして、この算出したM1個の比較候補立体画像と第2ステージの基準立体画像との相関値の中から最大相関値を抽出し(ステップS313)、この最大相関値が得られた比較候補立体画像の第1ステージでの比較確定回転軸JAを中心とする回転角度を第2ステージでの比較確定回転角度θBとする(ステップS314)。
なお、この例では、M1個の比較候補立体画像と第2ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の第1ステージでの比較確定回転軸JAを中心とする回転角度を第2ステージでの比較確定回転角度θBとしたが、M1個の比較候補立体画像と基準立体画像との相関値から関数フィッティング(例えば、シンク関数)により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の第1ステージでの比較確定回転軸JAを中心とする回転角度を第2ステージでの比較確定回転角度θBとするようにしてもよい。
次に、CPU4Aは、第1ステージでの比較確定回転軸JAを中心に任意のN1個の比較候補回転軸を設定する(ステップS315)。例えば、図12に示すように、第1ステージでの比較確定回転軸JAを中心軸とし、この比較確定回転軸JAの周囲に縦方向11個、横方向11個、合計121個の原点Oからの比較候補回転軸を3゜刻みで設定する。そして、このN1個の比較候補回転軸のそれぞれについて、その回転軸を中心として比較確定回転角度θB回転させたN1個の比較候補立体画像と第2ステージの基準立体画像との相関値を各個に算出する(ステップS316)。
そして、このN1個の比較候補立体画像と第2ステージの基準立体画像との相関値の中から最大相関値を抽出し(ステップS317)、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第2ステージでの比較確定回転軸JBとする(ステップS318)。
なお、この例では、N1個の比較候補立体画像と第2ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の比較候補回転軸を第2ステージでの比較確定回転軸JBとするようにしたが、N1個の比較候補立体画像と第2ステージの基準立体画像との相関値から関数フィッティング(例えば、シンク関数)により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第2ステージでの比較確定回転軸JBとするようにしてもよい。
〔第3ステージ(最終ステージ)〕
次に、CPU4Aは、記憶装置4Dに保存されている512×512×512の基準立体画像(原画像)のデータサイズを縮小し、第3ステージで用いる基準立体画像とする(図6:ステップS319)。また、記憶装置4Dに保存されている512×512×512の比較対象立体画像(原画像)のデータサイズを縮小し、第3ステージで用いる比較対象立体画像とする(ステップS320)。この例では、第3ステージで用いる基準立体画像および比較対象立体画像のデータサイズを64×64×64とし、第1ステージで用いる基準立体画像および比較対象立体画像のデータサイズよりも大きいサイズとする。
次に、CPU4Aは、記憶装置4Dに保存されている512×512×512の基準立体画像(原画像)のデータサイズを縮小し、第3ステージで用いる基準立体画像とする(図6:ステップS319)。また、記憶装置4Dに保存されている512×512×512の比較対象立体画像(原画像)のデータサイズを縮小し、第3ステージで用いる比較対象立体画像とする(ステップS320)。この例では、第3ステージで用いる基準立体画像および比較対象立体画像のデータサイズを64×64×64とし、第1ステージで用いる基準立体画像および比較対象立体画像のデータサイズよりも大きいサイズとする。
次に、CPU4Aは、第2ステージでの比較確定回転軸JBを中心とし、第2ステージでの比較確定回転角度θBを含む所定の角度範囲内で第3ステージの比較対象立体画像をM2通り回転させたM2個の画像を比較候補立体画像とし、このM2個の比較候補立体画像と第3ステージの基準立体画像との相関値を各個に算出する(ステップS321)。この例では、比較確定回転角度θBを中心とする±1゜の角度範囲内で0.1゜刻みで回転させた21個の画像を比較候補立体画像とし(図13参照)、この21個の比較候補立体画像と第3ステージの基準立体画像との相関値を各個に算出する。
そして、この算出したM2個の比較候補立体画像と第3ステージの基準立体画像との相関値の中から最大相関値を抽出し(ステップS322)、この最大相関値が得られた比較候補立体画像の第2ステージでの比較確定回転軸JBを中心とする回転角度を第3ステージでの比較確定回転角度θCとする(ステップS323)。
なお、この例では、M2個の比較候補立体画像と第3ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の第2ステージでの比較確定回転軸JBを中心とする回転角度を第3ステージでの比較確定回転角度θCとしたが、M2個の比較候補立体画像と第3ステージの基準立体画像との相関値から関数フィッティング(例えば、シンク関数)により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の第2ステージでの比較確定回転軸JBを中心とする回転角度を第3ステージでの比較確定回転角度θCとするようにしてもよい。
次に、CPU4Aは、第2ステージでの比較確定回転軸JBを中心に任意のN2個の比較候補回転軸を設定する(ステップS324)。例えば、図14に示すように、第2ステージでの比較確定回転軸JBを中心軸とし、この中心軸JBの周囲に縦方向7個、横方向7個、合計49個の原点Oからの比較候補回転軸を3゜刻みで設定する。そして、このN2個の比較候補回転軸のそれぞれについて、その回転軸を中心として回転角度θC回転させたN2個の比較候補立体画像と第3ステージの基準立体画像との相関値を各個に算出する(ステップS325)。
そして、このN2個の比較候補立体画像と第3ステージの基準立体画像との相関値の中から最大相関値を抽出し(ステップS326)、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第3ステージでの比較確定回転軸JCとする(ステップS327)。
なお、この例では、N2個の比較候補立体画像と第3ステージの基準立体画像との相関値中、最大相関値が得られた比較候補立体画像の比較候補回転軸を第3ステージでの比較確定回転軸JCとするようにしたが、N2個の比較候補立体画像と第3ステージの基準立体画像との相関値から関数フィッティング(例えば、シンク関数)により最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第3ステージでの比較確定回転軸JCとするようにしてもよい。
本実施の形態では、第3ステージを最終ステージとし、この第3ステージでの比較確定回転軸JCおよび比較確定回転角度θCを最終ステージでの比較確定回転軸JENDおよび比較確定回転角度θENDとする。なお、第3ステージを最終ステージとせず、さらに最終ステージまでのステージ数を増やしてもよい。この場合、各ステージでは、第3ステージと同様にして、そのステージで用いる基準立体画像および比較対象立体画像のデータサイズを大きくしながら、第2ステージと同様の処理を繰り返すものとする。
〔基準立体画像と比較対象立体画像との間の位置合わせ〕
CPU4Aは、このようにして最終ステージでの比較確定回転軸JENDおよび比較確定回転角度θENDを求めた後、この求めた最終ステージでの比較確定回転軸JENDおよび比較確定回転角度θENDに基づいて、記憶装置4に保存されている基準立体画像(原画像)と比較対象立体画像(原画像)との間の回転ずれを補正する(図7:ステップS401)。すなわち、比較対象立体画像(原画像)を比較確定回転軸JENDを中心として比較確定回転角度θEND回転補正する。これにより、一度の処理で、基準立体画像と比較対象立体画像との間の回転ずれが補正される。
CPU4Aは、このようにして最終ステージでの比較確定回転軸JENDおよび比較確定回転角度θENDを求めた後、この求めた最終ステージでの比較確定回転軸JENDおよび比較確定回転角度θENDに基づいて、記憶装置4に保存されている基準立体画像(原画像)と比較対象立体画像(原画像)との間の回転ずれを補正する(図7:ステップS401)。すなわち、比較対象立体画像(原画像)を比較確定回転軸JENDを中心として比較確定回転角度θEND回転補正する。これにより、一度の処理で、基準立体画像と比較対象立体画像との間の回転ずれが補正される。
次に、この回転ずれが補正された基準立体画像と比較対象立体画像との間で相関演算を行って、この基準立体画像と比較対象立体画像との間のX,Y,Zの3軸方向のずれ量Δx,Δy,Δzを求め、この求めた3軸方向のずれ量Δx,Δy,Δzに基づいて回転ずれが補正された基準立体画像と比較対象立体画像との間の平行ずれを補正する(ステップS402)。すなわち、X軸方向へΔx、Y軸方向へΔy、Z軸方向へΔz平行移動し、基準立体画像と比較対象立体画像とのX,Y,Zの3軸方向の位置を合わせる。
これにより、基準立体画像と比較対象立体画像との間の位置合わせが完了する。なお、この際、比較対象立体画像ではなく、基準立体画像の方を回転補正したり、平行移動するようにしてもよい。
〔位置合わせされた基準立体画像と比較対象立体画像との照合〕
CPU4Aは、このようにして基準立体画像と比較対象立体画像との位置合わせを行った後、すなわち基準立体画像と比較対象立体画像とのレジストレーションを行った後、位置合わせされた基準立体画像と比較対象立体画像とを照合する(ステップS403)。
CPU4Aは、このようにして基準立体画像と比較対象立体画像との位置合わせを行った後、すなわち基準立体画像と比較対象立体画像とのレジストレーションを行った後、位置合わせされた基準立体画像と比較対象立体画像とを照合する(ステップS403)。
この例では、位置合わせされた基準立体画像(ボクセルデータ群)と比較対象立体画像(ボクセルデータ群)との相関演算を行い、この相関演算の結果得られた相関値と予め定められた照合用の閾値とを比較することによって、基準立体画像と比較対象立体画像との照合を行う。
なお、3次元の立体画像間の相関演算については特許文献1などに示されているので、ここでの説明は省略する。この特許文献1において、N=3として3次元フーリエ変換、3次元逆フーリエ変換を施し、3次元の立体画像間の相関演算を行うことにより、相関成分エリア内でピークの位置ずれがX,Y,Zの各軸に現れる。このX,Y,Zの各軸に現れる位置ずれから基準立体画像と比較対象立体画像との相関が分かる。
CPU4Aは、ステップS402での照合の結果がOKであれば(ステップS404のYES)、照合OKの表示を行う(ステップS405)。照合の結果がNGであれば(ステップS404のNO)、照合NGの表示を行うとともに(ステップS406)、位置合わせされた基準立体画像と比較対象立体画像との差分画像をディスプレイ4I上に表示する(ステップS407)。
以上の説明から分かるように、本実施の形態によれば、基準立体画像および比較対象立体画像のデータサイズが最初は小さくされ、ステージが進むにつれて大きくされて行く。これにより、最初は少ないボクセルデータ群で概略的な比較確定回転軸と比較確定回転角度が素速く求められ、徐々にボクセルデータ群の数を増やして、より正確な比較確定回転軸と比較確定回転角度が求められて行く。このようにして、本実施の形態では、高速度でかつ高精度に、レジストレーションが行われる。この結果、照合スピードもアップし、照合結果の精度も高まる。
〔基準立体画像と比較対象立体画像の回転中心(原点)について〕
上述した実施の形態では、第1ステージの比較対象立体画像のX,Y,Zの3軸の中心を原点としてK個の回転軸を比較候補回転軸として設定し、このK個の比較候補回転軸を中心として回転させた比較対象立体画像と第1ステージの基準立体画像との相関値を求めるが、このような処理において比較対象立体画像と基準立体画像の回転中心とされる原点は同じでなければならない。しかし、3次元の立体画像において、実際の回転中心は不明であるため、原点を決めることが難しい。
上述した実施の形態では、第1ステージの比較対象立体画像のX,Y,Zの3軸の中心を原点としてK個の回転軸を比較候補回転軸として設定し、このK個の比較候補回転軸を中心として回転させた比較対象立体画像と第1ステージの基準立体画像との相関値を求めるが、このような処理において比較対象立体画像と基準立体画像の回転中心とされる原点は同じでなければならない。しかし、3次元の立体画像において、実際の回転中心は不明であるため、原点を決めることが難しい。
正確な原点を定めるために、本実施の形態において、ステップS301の処理に入る前に、記憶装置4Dに保存されている基準立体画像(原画像)および比較対象立体画像(原画像)を構成するボクセルデータ群に対してそれぞれ3次元フーリエ変換を施し、この3次元フーリエ変換が施された基準立体画像(原画像)および比較対象立体画像(原画像)を構成するボクセルデータ群の3次元周波数空間(ボクセルデータ群の中心を直流とする周波数空間)の中心を基準立体画像および比較対象立体画像に設定するX,Y,Zの3軸の中心とすることが考えられる。
この場合、相関値を求める処理は、以下の手順となる。
(1)両方のボクセルデータ群に対し、3次元フーリエ変換(FFT)、振幅化、を行う。この時、ボクセルデータ群の中心が直流になるようにする。また、3次元FFTした結果の振幅の値は、元のボクセルデータに比べて比較的大きな値をとることがある。これが演算上問題になる場合には、振幅データをログ化することで、ダイナミックレンジを抑えることができる。
(1)両方のボクセルデータ群に対し、3次元フーリエ変換(FFT)、振幅化、を行う。この時、ボクセルデータ群の中心が直流になるようにする。また、3次元FFTした結果の振幅の値は、元のボクセルデータに比べて比較的大きな値をとることがある。これが演算上問題になる場合には、振幅データをログ化することで、ダイナミックレンジを抑えることができる。
(2)(1)で得た(必要に応じてログ化された)振幅データの片方のボクセルデータ群に対して、ボクセルデータ群の中心を原点とした任意の回転軸を設定し、任意の回転角度だけ回転させる。
(3)(1)で得た回転させていない方(基準立体画像)のボクセルデータ群と、(2)で得た回転させた方のボクセルデータ群(比較対象立体画像)を3次元相関処理により、相関値を求める。相関値を求めるにには、3次元相関処理として、3次元位相限定相関法を含む任意の相関処理を用いることができる。
なお、ボクセルデータ群に対する3次元フーリエ変換の処理は、必ずしも必要というわけではない。例えば、被写体の回転角度が比較的小さいことが分かっている場合、回転中心がほぼ中心付近にあることが分かっている場合などは、3次元フーリエ変換の処理を行う必要はなく、基準立体画像および比較対象立体画像の原画像を構成するボクセルデータ群をそのまま使って、相関処理を行うことが可能である。
上述した実施の形態では、基準立体画像の原画像および比較対象立体画像の原画像のデータサイズを小さくする方式としているので、多少の回転中心のずれは無視することができる。また、被写体A,Bの撮影に際し、被写体A,Bはその頭部の中心を中間点CO付近にきちんと位置させるものとしているので、回転中心がほぼ中心付近にあることが分かっている。このようなことから、上述した実施の形態では、3次元フーリエ変換の処理を省略している。
〔回転軸と回転角度の設定について〕
初期段階において、比較対象立体画像に対して設定するK個の回転軸は全く不明である。このため、比較対象立体画像のボクセルデータ群の中心を原点とした任意の回転軸は、3次元上あらゆる方向に均等に分割されるように回転軸を選ぶ必要がある。回転軸の設定には、正多面体を用いることができる。例として、正二十面体を考えると、正二十面体の各面の中心座標を回転軸とすることで、球面に対して均等な20個の回転軸を設定することができる。また、正二十面体の頂点を回転軸とすることで、球面に対して均等な12個の回転軸を設定することができる。上述した実施の形態では、正二十面体の頂点を回転軸とし、原点を中心とする任意の方向へのK個の回転軸として12個の回転軸J1〜J12(図8参照)を設定している。
初期段階において、比較対象立体画像に対して設定するK個の回転軸は全く不明である。このため、比較対象立体画像のボクセルデータ群の中心を原点とした任意の回転軸は、3次元上あらゆる方向に均等に分割されるように回転軸を選ぶ必要がある。回転軸の設定には、正多面体を用いることができる。例として、正二十面体を考えると、正二十面体の各面の中心座標を回転軸とすることで、球面に対して均等な20個の回転軸を設定することができる。また、正二十面体の頂点を回転軸とすることで、球面に対して均等な12個の回転軸を設定することができる。上述した実施の形態では、正二十面体の頂点を回転軸とし、原点を中心とする任意の方向へのK個の回転軸として12個の回転軸J1〜J12(図8参照)を設定している。
回転角度については、対象物体が最大どの程度回転することがあるか、によって適宜設定すればよい。例えば、CTスキャナなどにより撮像されたボクセルデータであり、撮像時にある程度の位置合わせを行っているのであれば、回転誤差は目視で確認できない程度の角度であり、±10程度で十分であると考えられる。これを数度刻みで設定すればよい。また、回転軸を設定する際に、対向する回転軸を設定した場合には、この2つの軸に対する回転は、回転方向が逆になるだけなので、対向する回転軸を1つにする、または、回転角度を一方向にする、ことができる。このようなことを考慮して、上述した実施の形態では、回転軸J1〜J12を中心とする角度範囲を2〜32゜とし、2゜刻みで回転させるようにしている。
〔処理速度について〕
初期段階で、12個の回転軸(K=12)、16の回転角度(L=16)を選択した場合、192回の3次元回転と3次元相関処理が必要となる。これには膨大な演算を必要とするが、最初に求めるのは概略の回転軸と回転角度であるため、ボクセルデータ群を小さくしても問題はない。ボクセルデータ群を小さくすることで、演算量は少なくなり、高速処理が期待できる。このため、上述した実施の形態では、512×512×512の原画像のデータサイズを縮小し、16×16×16の縮小したデータサイズとしている。
初期段階で、12個の回転軸(K=12)、16の回転角度(L=16)を選択した場合、192回の3次元回転と3次元相関処理が必要となる。これには膨大な演算を必要とするが、最初に求めるのは概略の回転軸と回転角度であるため、ボクセルデータ群を小さくしても問題はない。ボクセルデータ群を小さくすることで、演算量は少なくなり、高速処理が期待できる。このため、上述した実施の形態では、512×512×512の原画像のデータサイズを縮小し、16×16×16の縮小したデータサイズとしている。
データサイズの一辺の数を1/2にした場合(ボクセル数では1/8)、通常のパーソナルコンピュータ(PC)にて3次元位相限定相関法を実行すると、処理時間は1/5から1/8程度になる。したがって、512×512×512のデータサイズのボクセルデータ群に対して3次元位相限定相関法を実行した場合の時間を1とすると、16×16×16のデータサイズのボクセルデータ群では、1/3125から1/32768となり、多数回実行しても問題ない処理スピードとなる。
〔合成ベクトルを作成する際の重み付けについて〕
合成ベクトルを作成する際に用いる重み付けは、相関処理の出力する値(相関度合いが高くなるほど高い値となる)をそのまま使用しても構わないが、相関処理の特性に合わせて任意の処理を行ってから使用することもある。例えば、相関処理の特性が緩慢な場合(完全に一致しても、少しずれていても、出力値が大きく変わらない場合)、出力値に相関関数を適用することもできるし、逆に急峻な場合には、平方根を適用することもできる。
合成ベクトルを作成する際に用いる重み付けは、相関処理の出力する値(相関度合いが高くなるほど高い値となる)をそのまま使用しても構わないが、相関処理の特性に合わせて任意の処理を行ってから使用することもある。例えば、相関処理の特性が緩慢な場合(完全に一致しても、少しずれていても、出力値が大きく変わらない場合)、出力値に相関関数を適用することもできるし、逆に急峻な場合には、平方根を適用することもできる。
〔機能ブロック図〕
図15、図16および図17に上述した3次元ずれ量計測装置4の要部の機能ブロック図を分割して示す。3次元ずれ量計測装置4は、基準立体画像の原画像(データサイズ:512×512×512)を記憶するデータベースDB1と、比較対象立体画像の原画像(データサイズ:512×512×512)を記憶するデータベースDB2とを備えている。
図15、図16および図17に上述した3次元ずれ量計測装置4の要部の機能ブロック図を分割して示す。3次元ずれ量計測装置4は、基準立体画像の原画像(データサイズ:512×512×512)を記憶するデータベースDB1と、比較対象立体画像の原画像(データサイズ:512×512×512)を記憶するデータベースDB2とを備えている。
また、第1ステージ用として、第1の基準立体画像データサイズ縮小部4−1と、第1の比較対象立体画像データサイズ縮小部4−2と、第1の比較候補回転軸設定部4−3と、第1の比較候補立体画像作成部4−4と、第1の相関値算出部4−5と、第1の最大相関値抽出部4−6と、第1ステージ比較確定回転角度決定部4−7と、相関値抽出部4−8と、第1ステージ比較回転軸決定部4−9とを備えている。
また、第2ステージ用として、第2の比較候補立体画像作成部4−10と、第2の相関値算出部4−11と、第2の最大相関値抽出部4−12と、第2ステージ比較確定回転角度決定部4−13と、第2の比較候補回転軸設定部4−14と、第3の比較候補立体画像作成部4−15と、第3の相関値算出部4−16と、第3の最大相関値抽出部4−17と、第2ステージ比較回転軸決定部4−18とを備えている。なお、第1ステージ用として設けた第1の基準立体画像データサイズ縮小部4−1と第1の比較対象立体画像データサイズ縮小部4−2は、第2ステージ用としても用いられる。
また、第3ステージ用として、第2の基準立体画像データサイズ縮小部4−19と、第2の比較対象立体画像データサイズ縮小部4−20と、第4の比較候補立体画像作成部4−21と、第4の相関値算出部4−22と、第4の最大相関値抽出部4−23と、第3ステージ比較確定回転角度決定部4−24と、第3の比較候補回転軸設定部4−25と、第5の比較候補立体画像作成部4−26と、第5の相関値算出部4−27と、第5の最大相関値抽出部4−26と、第3ステージ比較回転軸決定部4−29とを備えている。
第1の基準立体画像データサイズ縮小部4−1は、データベースDB1に保存されている512×512×512の基準立体画像(原画像)GRのデータサイズを縮小し、第1ステージおよび第2ステージで用いる基準立体画像(データサイズ:16×16×16)GR1とする。第1の比較対象立体画像データサイズ縮小部4−2は、データベースDB2に保存されている512×512×512の比較対象立体画像(原画像)GCのデータサイズを縮小し、第1ステージおよび第2ステージで用いる比較対象立体画像(データサイズ:16×16×16)GC1とする。
第1の比較回転軸設定部4−3は、第1ステージの比較対象立体画像GC1のX,Y,Zの3軸の中心を原点とし、この原点を中心として任意の方向へK個の回転軸を比較候補回転軸J1〜JKとして設定する。第1の比較候補立体画像作成部4−4は、第1の比較回転軸設定部4−3によって作成されたK個の比較候補回転軸を中心として第1ステージの比較対象立体画像GC1を所定の角度範囲内でL通り回転させたK×L個の画像を比較候補立体画像として作成する。第1の相関値算出部4−5は、第1の比較候補立体画像作成部4−4によって作成されたK×L個の比較候補立体画像と第1ステージの基準立体画像GR1との相関値を各個に算出する。
第1の最大相関値抽出部4−6は、第1の相関値算出部4−5によって算出されたK×L個の相関値の中から最大相関値を抽出する。第1ステージ比較確定回転角度決定部4−7は、最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第1ステージでの比較確定回転角度θAとして決定する。
相関値抽出部4−8は、第1の相関値算出部4−5によって算出されたK×L個の相関値の中から、第1ステージでの比較確定回転角度θAと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する。
第1ステージ比較確定回転軸決定部4−9は、相関値抽出部4−8によって抽出された相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトル(原点からの単位長さ)に当該相関値に基づく重み付けを施し、この重み付けを施した単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を比較対象立体画像GCの第1ステージでの比較確定回転軸JAとする。
第2の比較候補立体画像作成部4−10は、第1ステージでの比較確定回転軸JAを中心とし、第1ステージでの比較確定回転角度θAを含む所定の角度範囲内で比較対象立体画像GC1をM1通り回転させたM1個の画像を比較候補立体画像として作成する。第2の相関値算出部4−11は、第2の比較候補立体画像作成部4−10によって作成されたM1個の比較候補立体画像と基準立体画像GR1との相関値を各個に算出する。
第2の最大相関値抽出部4−12は、第2の相関値算出部4−11によって算出されたM1個の相関値の中から最大相関値を抽出する。第2ステージ比較確定回転角度決定部4−13は、最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第2ステージでの比較確定回転角度θBとして決定する。
第2の比較候補回転軸設定部4−14は、第1ステージでの比較確定回転軸JAを中心にN1個の比較候補回転軸を設定する。第3の比較候補立体画像作成部4−15は、第2の比較候補回転軸設定部4−14によって設定されたN1個の比較候補回転軸のそれぞれについて、その回転軸を中心として比較確定回転角度θB回転させたN1個の比較候補立体画像を作成する。第3の相関値算出部4−16は、第3の比較候補立体画像作成部4−15によって作成されたN1個の比較候補立体画像と基準立体画像GR1との相関値を各個に算出する。
第3の最大相関値抽出部4−17は、第3の相関値算出部4−16によって算出されたN1個の相関値の中から最大相関値を抽出する。第2ステージ比較確定回軸決定部4−18は、最大相関値が得られた比較候補立体画像の比較候補回転軸を第2ステージでの比較確定回転軸JBとして決定する。
第2の基準立体画像データサイズ縮小部4−19は、データベースDB1に保存されている512×512×512の基準立体画像(原画像)GRのデータサイズを縮小し、第3ステージで用いる基準立体画像(データサイズ:64×64×64)GR2とする。第1の比較対象立体画像データサイズ縮小部4−20は、データベースDB2に保存されている512×512×512の比較対象立体画像(原画像)GCのデータサイズを縮小し、第3ステージで用いる比較対象立体画像(データサイズ:64×64×64)GC2とする。
第4の比較候補立体画像作成部4−21は、第2ステージでの比較確定回転軸JBを中心とし、第2ステージでの比較確定回転角度θBを含む所定の角度範囲内で比較対象立体画像GC2をM2通り回転させたM2個の画像を比較候補立体画像として作成する。第4の相関値算出部4−22は、第4の比較候補立体画像作成部4−21によって作成されたM2個の比較候補立体画像と基準立体画像GR2との相関値を各個に算出する。
第4の最大相関値抽出部4−23は、第4の相関値算出部4−22によって算出されたM2個の相関値の中から最大相関値を抽出する。第3ステージ比較確定回転角度決定部4−24は、最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第3ステージでの比較確定回転角度θCとして決定する。
第3の比較候補回転軸設定部4−25は、第2ステージでの比較確定回転軸JBを中心にN2個の比較候補回転軸を設定する。第5の比較候補立体画像作成部4−26は、第3の比較候補回転軸設定部4−25によって設定されたN2個の比較候補回転軸のそれぞれについてその回転軸を中心として回転角度θC回転させたN2個の比較候補立体画像を作成する。、第5の相関値算出部4−27は、第5の比較候補立体画像作成部4−26によって作成されたN2個の比較候補立体画像と基準立体画像GR2との相関値を各個に算出する。
第5の最大相関値抽出部4−28は、第5の相関値算出部4−27によって算出されたN2個の相関値の中から最大相関値を抽出する。第3ステージ比較確定回軸決定部4−29は、最大相関値が得られた比較候補立体画像の比較候補回転軸を第3ステージでの比較確定回転軸JCとして決定する。
なお、上述した実施の形態では、最初に撮影した被写体Aの3次元の立体画像を基準立体画像とし、次に撮影した被写体Bの3次元の立体画像を比較対象立体画像としたが、最初に撮影した被写体Aの3次元の立体画像を比較対象立体画像とし、次に撮影した被写体Bの3次元の立体画像を基準立体画像としてもよい。また、基準立体画像や比較対象立体画像は、ボクセルデータ群として上位装置から与えられるものとしてもよい。
また、上述した実施の形態では、被写体を人体としたが、被写体は人体に限られるものではない。例えば、製薬などの分野において、マウスを使用した実験などにも利用することができる。また、産業分野において、金属などの非破壊検査に利用することもできる。このように、本発明は、医療や歯科などの分野に限らず、各種の分野において有効利用することが可能である。
また、上述した3次元ずれ量計測処理装置4において、図15〜図17に分割して機能ブロックとして示した各部の処理動作は、プログラム(3次元ずれ量計測プログラム)に従うCPUの処理動作として実現されるが、同様の処理をハードウェアの回路構成で実現するようにしてもよい。
また、上述した実施の形態では、基準立体画像と比較対象立体画像との間のずれ量の計測処理を第3ステージを最終ステージとして実行するようにしたが、第1ステージで終了するようにしてもよく、第2ステージで終了するようにしてもよい。また、第1ステージや第2ステージで終了するような場合、基準立体画像および比較対象立体画像の原画像のデータサイズは、必ずしも縮小しなくてもよい。
1…X線源、2…2次元撮像装置、3(A,B)…被写体、4A…CPU、4B…RAM、4C…ROM、4D…記憶装置、4E〜H…インタフェース、4I…ディスプレイ、4J…キーボード、4K…マウス、DB1,DB2…データベース、4−1…第1の基準立体画像データサイズ縮小部、4−2…第1の比較対象立体画像データサイズ縮小部、4−3…第1の比較候補回転軸設定部、4−4…第1の比較候補立体画像作成部、4−5…第1の相関値算出部、4−6…第1の最大相関値抽出部、4−7…第1ステージ比較確定回転角度決定部、4−8…相関値抽出部、4−9…第1ステージ比較回転軸決定部、4−10…第2の比較候補立体画像作成部、4−11…第2の相関値算出部、4−12…第2の最大相関値抽出部、4−13…第2ステージ比較確定回転角度決定部、4−14…第2の比較候補回転軸設定部、4−15…第3の比較候補立体画像作成部、4−16…第3の相関値算出部、4−17…第3の最大相関値抽出部、4−18…第2ステージ比較回転軸決定部、4−19…第2の基準立体画像データサイズ縮小部、4−20…第2の比較対象立体画像データサイズ縮小部、4−21…第4の比較候補立体画像作成部、4−22…第4の相関値算出部、4−23…第4の最大相関値抽出部、4−24…第3ステージ比較確定回転角度決定部、4−25…第3の比較候補回転軸設定部、4−26…第5の比較候補立体画像作成部、4−27…第5の相関値算出部、4−26…第5の最大相関値抽出部、4−29…第3ステージ比較回転軸決定部。
Claims (9)
- 基準となる3次元の立体画像を基準立体画像、この基準立体画像に対して比較対象とする3次元の立体画像を比較対象立体画像とし、前記基準立体画像と前記比較対象立体画像との間のずれ量を計測する3次元ずれ量計測方法において、
前記比較対象立体画像に設定される互いに直交するX,Y,Zの3軸の中心を原点とし、この原点を中心として任意の方向へK(K≧2)個の回転軸を比較候補回転軸として設定する第1ステップと、
前記K個の各比較候補回転軸を中心として前記比較対象立体画像を所定角度範囲内でL(L≧2)通り回転させたK×L個の画像を比較候補立体画像とし、このK×L個の比較候補立体画像と前記基準立体画像との相関値を各個に算出する第2ステップと、
前記算出されたK×L個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を中心とする回転角度を第1ステージでの比較確定回転角度θAとして求める第3ステップと、
前記算出されたK×L個の比較候補立体画像と基準立体画像との相関値の中から前記比較確定回転角度θAと同じ回転角度でその比較候補回転軸を回転させた比較候補立体画像の相関値を抽出する第4ステップと、
前記比較候補回転軸の前記原点からの単位長さを単位ベクトルとし、前記第4ステップで抽出された相関値の抽出元の比較候補立体画像の比較候補回転軸の単位ベクトルに当該相関値に基づく重み付けを施し、この重み付けが施された単位ベクトルを合成して合成ベクトルを求め、この合成ベクトルが示す方向を前記比較対象立体画像の第1ステージでの比較確定回転軸JAとする第5ステップと
を備えることを特徴とする3次元ずれ量計測方法。 - 請求項1に記載された3次元ずれ量計測方法において、
前記第1ステージでの比較確定回転軸JAを中心とし、前記第1ステージでの比較確定回転角度θAを含む所定角度範囲内で前記比較対象立体画像をM(M≧2)通り回転させたM個の画像を比較候補立体画像とし、このM個の比較候補立体画像と前記基準立体画像との相関値を各個に算出する第6ステップと、
前記算出されたM個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の前記第1ステージでの比較確定回転軸JAを中心とする回転角度を第2ステージでの比較確定回転角度θBとして求める第7ステップと、
前記第1ステージでの比較確定回転軸JAを中心に任意のN(N≧2)個の比較候補回転軸を設定し、このN個の比較候補回転軸のそれぞれについてその回転軸を中心として前記比較確定回転角度θB回転させたN個の比較候補立体画像と前記基準立体画像との相関値を各個に算出する第8ステップと、
前記N個の比較候補立体画像と基準立体画像との相関値の中から最大相関値を抽出し、この最大相関値が得られた比較候補立体画像の比較候補回転軸を第2ステージでの比較確定回転軸JBとする第9ステップと
を備えることを特徴とする3次元ずれ量計測方法。 - 請求項2に記載された3次元ずれ量計測方法において、
前記第7ステップに代えて、
前記算出されたM個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の前記第1ステージでの比較確定回転軸JAを中心とする回転角度を第2ステージでの比較確定回転角度θBとして求める第7ステップ
を備えることを特徴とする3次元ずれ量計測方法。 - 請求項2に記載された3次元ずれ量計測方法において、
前記第9ステップに代えて、
前記N個の比較候補立体画像と基準立体画像との相関値から最大相関値を有する比較候補立体画像を推測し、この推測した比較候補立体画像の比較候補回転軸を第2ステージでの比較確定回転軸JBとする第9ステップ
を備えることを特徴とする3次元ずれ量計測方法。 - 請求項2〜4の何れか1項に記載された3次元ずれ量計測方法において、
前記第1ステップの処理に入る前に、前記基準立体画像および比較対象立体画像の原画像のデータサイズを小さくして前記第1ステージで用いる基準立体画像および比較対象立体画像とするステップと、
前記第1ステージで用いる基準立体画像および比較対象立体画像よりもそのデータサイズを大きくしながら前記第6〜第9ステップの処理を予め定められた最終ステージまで繰り返すステップと
を備えることを特徴とする3次元ずれ量計測方法。 - 請求項5に記載された3次元ずれ量計測方法において、
前記最終ステージでの比較確定回転角度および比較確定回転軸に基づいて前記基準立体画像の原画像と前記比較対象立体画像の原画像との間の位置合わせを行う位置合わせ処理ステップ
を備えることを特徴とする3次元ずれ量計測方法。 - 請求項6に記載された3次元ずれ量計測方法において、
前記位置合わせ処理ステップによって位置合わせされた前記基準立体画像の原画像と前記比較対象立体画像の原画像とを相関演算によって照合する照合ステップ
を備えることを特徴とする3次元ずれ量計測方法。 - 請求項6又は7に記載された3次元ずれ量計測方法において、
前記位置合わせ処理ステップは、
前記最終ステージでの比較確定回転角度および比較確定回転軸に基づいて前記基準立体画像の原画像と前記比較対象立体画像の原画像との間の回転ずれを補正し、この回転ずれが補正された基準立体画像の原画像と比較対象立体画像の原画像との間で相関演算を行って当該基準立体画像の原画像と比較対象立体画像の原画像との間の前記X,Y,Zの3軸方向のずれ量を求め、この求めた3軸方向のずれ量に基づいて前記回転ずれが補正された基準立体画像の原画像と比較対象立体画像の原画像との間の平行ずれを補正する
ことを特徴とする3次元ずれ量計測方法。 - 請求項1〜8の何れか1項に記載された3次元ずれ量計測方法において、
前記第1ステップの処理に入る前に、前記基準立体画像の原画像および前記比較対象立体画像の原画像を構成するボクセルデータ群に対してそれぞれ3次元フーリエ変換を施し、この3次元フーリエ変換が施された基準立体画像の原画像および比較対象立体画像の原画像を構成するボクセルデータ群の3次元周波数空間の中心を前記基準立体画像および比較対象立体画像に設定される互いに直交するX,Y,Zの3軸の中心とするステップ
を備えることを特徴とする3次元ずれ量計測方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2008303867A JP5171575B2 (ja) | 2008-11-28 | 2008-11-28 | 3次元ずれ量計測方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2008303867A JP5171575B2 (ja) | 2008-11-28 | 2008-11-28 | 3次元ずれ量計測方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
JP2010125117A JP2010125117A (ja) | 2010-06-10 |
JP5171575B2 true JP5171575B2 (ja) | 2013-03-27 |
Family
ID=42325908
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2008303867A Expired - Fee Related JP5171575B2 (ja) | 2008-11-28 | 2008-11-28 | 3次元ずれ量計測方法 |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP5171575B2 (ja) |
Family Cites Families (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP3773111B2 (ja) * | 2002-07-16 | 2006-05-10 | 独立行政法人科学技術振興機構 | 3次元物体の投影画像位置合わせシステム |
JP4828920B2 (ja) * | 2004-11-24 | 2011-11-30 | 株式会社東芝 | 3次元画像処理装置 |
JP4991181B2 (ja) * | 2005-05-10 | 2012-08-01 | 株式会社東芝 | 3次元画像処理装置、3次元画像処理方法及び3次元画像処理装置で使用される制御プログラム |
US8340374B2 (en) * | 2007-01-11 | 2012-12-25 | Kabushiki Kaisha Toshiba | 3-dimensional diagnostic imaging system |
-
2008
- 2008-11-28 JP JP2008303867A patent/JP5171575B2/ja not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
JP2010125117A (ja) | 2010-06-10 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US9955951B2 (en) | Sensor coordinate calibration in an ultrasound system | |
JP2019526124A (ja) | 三次元表面の画像を再構築するための方法、装置およびシステム | |
US20160007970A1 (en) | Segmentation of large objects from multiple three-dimensional views | |
JP2010154920A (ja) | 3次元形状データ処理装置、3次元形状データ処理方法及びプログラム | |
US10278663B2 (en) | Sensor coordinate calibration in an ultrasound system | |
WO2021104153A1 (en) | Systems and methods for automated calibration | |
CN108113700B (zh) | 一种应用于三维超声成像数据采集中的位置标定方法 | |
US10078906B2 (en) | Device and method for image registration, and non-transitory recording medium | |
US9020215B2 (en) | Systems and methods for detecting and visualizing correspondence corridors on two-dimensional and volumetric medical images | |
US9633433B1 (en) | Scanning system and display for aligning 3D images with each other and/or for detecting and quantifying similarities or differences between scanned images | |
CN110133014A (zh) | 一种芯片内部缺陷检测方法及系统 | |
JP5526401B2 (ja) | 心室壁情報抽出装置 | |
JP2011120825A (ja) | 医用画像表示装置および方法並びにプログラム | |
JP6429958B2 (ja) | 画像処理装置、画像処理方法、及びプログラム | |
JP5171575B2 (ja) | 3次元ずれ量計測方法 | |
JP5219759B2 (ja) | 3次元ずれ量計測方法 | |
CN111904450B (zh) | 左心室中心和感兴趣区域的提取方法、装置及系统 | |
JP7083427B2 (ja) | 修正指示領域表示装置、方法およびプログラム | |
CN104700419A (zh) | 一种放射科x光片的图像处理方法 | |
Kalmykova et al. | An approach to point-to-point reconstruction of 3D structure of coronary arteries from 2D X-ray angiography, based on epipolar constraints | |
JP6391544B2 (ja) | 医用画像処理装置、医用画像処理方法、及びプログラム | |
US11664116B2 (en) | Medical image data | |
Kim et al. | Dense femur reconstruction from two x-ray images using generic 3D model with twist correction | |
CN117541493A (zh) | 一种基于改进分数阶累积量的三维特征医学图像融合方法 | |
CN117058197A (zh) | 用于神经调控的空间注册方法及装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20110928 |
|
A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20121218 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20121225 |
|
LAPS | Cancellation because of no payment of annual fees |