JP5706933B2

JP5706933B2 - 処理装置および処理方法、プログラム

Info

Publication number: JP5706933B2
Application number: JP2013159711A
Authority: JP
Inventors: 亮石川; 池内　克史; 克史池内; 淳高松; 岳史大石; 波鄭
Original assignee: Canon Inc
Current assignee: Canon Inc
Priority date: 2013-07-31
Filing date: 2013-07-31
Publication date: 2015-04-22
Anticipated expiration: 2028-05-13
Also published as: JP2013248466A

Description

本発明は、複数の画像データをコンピュータ処理し、画像の位置合わせをする処理装置に関する。

医療の分野において、医師は、患者を撮影した医用画像をモニタに表示し、表示された医用画像を読影して、病変部の状態や経時変化を観察する。この種の医用画像を生成する装置としては、単純Ｘ線撮影装置、Ｘ線コンピュータ断層撮影装置（X線CT）、核磁気共鳴映像装置（MRI）、核医学診断装置（SPECT、PETなど）、超音波画像診断装置（US）等が挙げられる。

また、より精度の良い診断を目的として複数の装置から得た複数の画像を利用して、診断を行う場合もある。例えば同一の被検査者に対して複数の撮像装置を用いて撮影し、得られた画像を重ね合わせることで、診断や、レーザ照射等の治療に有効な情報を得ることができる。

例えば、超音波画像診断装置から得られる画像とＣＴから得られる三次元画像とを重ね合わせ、治療用のレーザ照射のより正確な位置を得ている。そのため、画像の重ね合わせの高精度化が要求されている。しかし、医師は目視で双方の画像の位置合わせを行っているのが現状である。

一方、画像間の位置合わせを高速に行う方法が研究されつつある（非特許文献１）。
この方法では、以下の１)〜４）を収束するまで計算することにより両者の相対位置を推定する。
１）二次元画像の輪郭線を抽出し、輪郭線からの距離に応じた値を示す距離場を生成する。
２）それに対して、三次元幾何モデルのシルエット画像の輪郭を求め、輪郭上の点に距離場に応じて計算される力を加える。
３）そして、全ての輪郭線上に加えられる力の和と、三次元幾何モデル重心周りのモーメントを求める。
４）得られた力とモーメントに応じて、三次元幾何モデルの位置姿勢を更新する。

2次元距離場を用いた2D-3Dレジストレーション岩下、倉爪、原、油谷、中本、小西、橋爪、長谷川、画像の認識・理解シンポジウム(MIRU3005), 3005 Michael M. Blane, Zhibin Lei, HakanCE ivi, and David B. Cooper, "The 3L Algorithm for Fitting Implicit Polynomial Curves and Surfaces to Data," IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 22, No. 3, pp. 298-313, 3000 Tolga Tasdizen, Jean-PhilippeTarel, David B. Cooper, "Improving the Stability of Algebraic Curves for Applications," IEEE Transactions on Image Processing, Vol. 9, No. 3, pp. 405-416, 3000 Amir Helzer, Meir Barzohar, and David Malah, "Stable Fitting of 2D Curves and 3D Surfaces by Implicit Polynomials," IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 26, No. 10, pp. 1283-1294, 3004 BoZheng, Jun Takamatsu, and Katsushi Ikeuchi, "Adaptively Determining Degrees of Implicit Polynomial Curves and Surfaces," Proc. 8th Asian Conference on Computer Vision, 3007 G. Taubin and D. Cooper, "Symbolic and Numerical Computation for Artificial Intelligence, chapter 6," Computational Mathematics and Applications, Academic Press, 1992

非特許文献１に開示される技術は、輪郭線からの距離に応じた値をしめす距離場を数式で表す思想がなく、距離場を示す値をメモリに記憶しておく必要がある。そのため、上記３）の力などを計算する場合には、三次元幾何モデルのシルエット画像の輪郭座標に対応するメモリ上の座標を求める処理が必要となる。従って、メモリとのアクセスなしに数値計算として解を求めることが出来ない。また、位置合わせの効率化を計るために、メモリのアクセスが制限となり計算の簡略化を計ることも困難である。これは、位置合わせの精度へも影響するものである。

本発明は、上記の課題に鑑みてなされたものであり、複数の画像を精度よく位置合わせ処理する仕組みを提供することを目的とする。

上記の目的を達成するための、本発明の一態様によるによる処理装置は以下の構成を備える。すなわち、
第一の画像の第一の領域と第二の画像の第二の領域との位置合わせを行う処理装置であって、
前記第一の画像における位置合わせの対象となる前記第一の領域からの距離に応じた値を算出する多項式を生成する生成手段と、
前記第二の領域の前記第二の画像における位置座標を前記第一の画像上に座標変換した位置座標と前記多項式に基づいて前記位置合わせを行うための補正値を得る取得手段と、を備える

本発明の構成により、異なる撮像装置で撮像された二つの異なる画像間の位置合わせを高精度に行う仕組みを提供できる。

実施形態に係る位置合わせ処理装置の機能構成を示す図である。実施形態に係る位置合わせ処理システムの装置構成を示す図である。実施形態に係る位置合わせ処理装置の処理手順を示すフローチャートである。実施形態に係る図２のステップＳ３０２の処理を説明する図である。実施形態に係る図２のステップＳ３０３の処理を説明する図である。実施形態に係る図２のステップＳ３０５の処理を説明する図である。実施形態に係る位置合わせ処理を説明する図である。

以下、添付図面に従って本発明に係る位置合わせ処理装置及び方法の好ましい実施形態について詳説する。ただし、発明の範囲は図示例に限定されるものではない。

図１は、本実施形態に係る位置合わせ処理装置１０００の機能構成を示している。位置合わせ処理装置１０００は第一の撮影装置１１００および、第二の撮影装置１１１０と接続されている。

これらの撮影装置として、単純X線撮影装置、X線コンピュータ断層撮影装置（X線CT）、核磁気共鳴映像装置（MRI）、核医学診断装置（SPECT、PETなど）、超音波画像診断装置（US）等の撮影装置が挙げられる。また、医療用の画像のみでなく、一般のカメラで撮影した画像にも適用できるものである。

第一の画像入力部１０１０は、第一の撮影装置１１００が撮影した対象物体の画像を第一の画像として入力する。

生成部１０２０は、入力された第一の画像を処理して、位置合わせの対象となる領域から複数の位置座標を取得し、この複数の位置座標に基づく多項式を生成する。ここで位置合わせの対象となる領域とは、注目領域のことであり、例えば、医療画像であれば肺、胃、心臓等の人体内部の臓器や、頭、手等が相当する。一般のカメラ画像であれば、人の顔などが相当する。位置合わせの対象となる領域からの距離とは、例えば肺ならば肺と他の人体内部組織との境界領域からの距離を意味する。すなわち、位置合わせの対象となる領域を三次元構成する外面からの距離を意味する。

また、二次元画像であれば、位置合わせの対象となる領域からの距離とは、位置合わせの対象となる領域と外部領域の境界線である輪郭線からの距離を意味する。

第二の画像入力部１０３０は、第二の撮影装置１１１０が撮影した対象物体の画像を第二の画像として入力する。

取得部１０４０は、第二の画像における位置合わせの対象となる領域から複数の位置座標を取得する。

なお、本実施形態では三次元画像をI_3D(x,y,z)と表記する。I_3D(x,y,z)は三次元空間の位置座標(x,y,z)における三次元画像の画素値を表現したものである。二次元画像であれば、Ｚ＝０として処理される。また、この画素値は撮影装置毎に物理的な意味が異なる。また、上記座標系は直交座標系を用いて表記しているがこれにかぎらず極座標系などを用いて表現してもよい。

例えば、単純Ｘ線撮影装置は、Ｘ線を人体に対して照射し、その透過Ｘ線を記録することで、人体内部のX線吸収率の平行投影像を得る。すなわち、Ｘ線吸収率に関する情報が画素値となる。

また、Ｘ線ＣＴはＸ線を人体に対して様々な方向から照射し、多数の透視像を得て、それらの像を解析することで、人体内部の三次元的な情報を得る。このような三次元的な情報をＣＴ撮影ではボクセルデータとよぶ。すなわち、Ｘ線ＣＴから得られた画像では、ボクセル値に関する情報が画素値となる。

また、ＭＲＩは、ＣＴと同様に人体内部の３次元的な情報を得ることができる装置であるが、ＭＲＩは磁気共鳴現象を利用して画像化する装置であるため、Ｘ線吸収率を画像化するＣＴとは異なる物理情報が得られる。

また、超音波画像診断装置は人体に超音波を照射して、人体の内部から反射してくる超音波を検出することで、人体内部の情報を取得する。超音波画像診断装置は通常、Ｂモードと称される方法で人体内部の断層画像を撮影する。超音波画像診断装置は、Ｘ線による被爆などといった人体に対する浸襲が無く、また撮影と観察を同時に行えるといった特徴がある。
以上のように、撮影装置毎に、物理的な意味が異なる情報が画素値として取得される。

算出部１０５０は、第二の画像から取得した複数の位置座標を座標変換して得た位置座標に対して、第一の画像における位置合わせの対象となる領域からの距離に応じた値を、生成部１０２０で生成された多項式を用いて夫々算出する。
決定部１０６０は、算出部１０５０で算出された値に基づいて、前記複数の位置座標の座標変換方法を決定する。座標変換方法に関しては後述する。

座標変換部１０７０は、決定部１０６０で決定した座標変換方法で、前記第二の画像の位置座標を座標変換する。
また、画像合成部１０８０は、第一の画像と座標変換した第二の画像を合成した合成画像を生成し、画像表示装置１１３０に出力する。画像表示装置１１３０は、入力された画像を表示する。

図２は、実施形態に係る位置合わせ処理装置１０００を用いたシステムの装置構成を示す図である。本実施形態の位置合わせ処理システムは、位置合わせ処理装置１０００、第一の撮影装置１１００、画像記録装置３、ローカルエリアネットワーク（ＬＡＮ）４、第二の撮影装置１１１０、外部位置計測センサ６を有する。

位置合わせ処理装置１０００は、例えばパーソナルコンピュータ（ＰＣ）などで実現することができる。即ち、位置合わせ処理装置１０００は、中央演算処理装置（ＣＰＵ）１００、主メモリ１０１、磁気ディスク１０２、表示メモリ１０３、モニタ１０４、マウス１０５、キーボード１０６を有する。

ＣＰＵ１００は、主として位置合わせ処理装置１０００の各構成要素の動作を制御する。主メモリ１０１は、ＣＰＵ１００が実行する制御プログラムを格納したり、ＣＰＵ１００によるプログラム実行時の作業領域を提供したりする。磁気ディスク１０２は、オペレーティングシステム（ＯＳ）、周辺機器のデバイスドライブ、後述する位置合わせ処理等を行うためのプログラムを含む各種アプリケーションソフト等を格納する。表示メモリ１０３は、モニタ１０４のための表示用データを一時記憶する。モニタ１０４は、例えばＣＲＴモニタや液晶モニタ等であり、表示メモリ１０３からのデータに基づいて画像を表示する。表示メモリ１０３とモニタ１０４は表示部を構成する。マウス１０５及びキーボード１０６はユーザによるポインティング入力及び文字等の入力をそれぞれ行う。上記各構成要素は共通バス１０７により互いに通信可能に接続されている。

本実実施形態において、位置合わせ処理装置１０００は、ＬＡＮ４を介して画像記録装置３から三次元画像データ等を読み出すことができる。また、ＬＡＮ４を経由して第一の撮影装置１１００から直接に三次元画像、二次元画像のいずれかを取得できるようにしてもよい。なお、本発明の形態はこれに限定されず、例えば位置合わせ処理装置１０００に記憶装置、例えばＦＤＤ、ＣＤ−ＲＷドライブ、ＭＯドライブ、ＺＩＰドライブ等を接続し、それらのドライブから画像データ等を読み込むようにしても良い。また、第一の撮影装置１１００、第二の撮影装置１１１０としては、例えばＸ線ＣＴ装置、ＭＲＩ装置、超音波画像診断装置、核医学装置、一般的なデジタルカメラなどが挙げられる。

また位置合わせ処理装置１０００は、第二の撮影装置１１１０と接続されていて、撮影した画像を入力することができる。本実施形態では超音波画像診断装置を第二の撮影装置１１１０とした例に関して説明する。

位置合わせ処理装置１０００と第二の撮影装置１１１０とは直接に接続されても良いし、ＬＡＮ４を介して接続されても良い。また、第二の撮像装置５から画像記録装置３に画像を記録し、位置合わせ処理装置１０００が画像記録装置３から画像を読み出せるようにしても良い。

ここで、I_3D(x,y,z)で示される画像は、装置内の処理では画像データとして処理され、表示部における表示では可視画像として表示される。それゆえ、本実施形態では、I_3D(x,y,z)で示される画像データを画像データ又は画像と呼ぶものとする。

外部位置計測センサ６は、第二の撮影装置１１１０の不図示の撮影プローブに装着され、その位置・姿勢などを計測する。また、外部位置計測センサ６と第二の撮影装置１１１０の撮影範囲、三次元画像との間の関係は事前に校正が行われているものとする。位置合わせ処理装置１０００は外部位置計測センサ６の出力信号を計測結果として読み出すことで、第一の画像を基準とした座標系で、第二の撮影装置１１１０の撮影範囲をおおよその精度で取得できるものである。

次に、図３のフローチャートを用いて、位置合わせ処理装置１０００の動作に関して説明する。

ステップＳ３０１において、第一の画像入力部１０１０は、第一の撮影装置１１００によって撮影した、被検査体の三次元画像データ又は二次元画像データを第一の画像として入力する。ここで第一の画像は第一の撮影装置１１００から位置合わせ処理装置１０００に直接入力しても良い。あるいは、第一の撮影装置１１００が撮影した画像を、図２に記載の画像記録装置３に記録し、第一の画像入力部１０１０は画像記録装置３から所望の三次元画像データを読み出して入力してもよい。これら装置間の情報伝達は、図２に示すように例えばＬＡＮ４を介して行うことができ、その際の情報伝達のプロトコルとしてＤＩＣＯＭ形式を利用する事ができる。
入力された三次元画像I_3D(x,y,z)は、生成部１０２０に送信される。

ステップＳ３０２において、生成部１０２０は、ステップＳ３０１で入力された三次元画像から位置合わせの対象となる領域の輪郭点を抽出する。ここで位置合わせの対象となる領域とは、注目領域のことであり、例えば、医療画像であれば肺、胃、心臓、頭、手等の部位等が相当する。

この処理について図４を用いて説明する。図４の（ａ）はステップＳ３０１で入力した三次元画像の断層像（第１の画像）を、二次元画像として表した模式図である。

本実施形態では三次元画像で説明するが、二次元画像にも同様に適用できるものである。

図４の（ｂ）の画像は図４の（ａ）の画像から位置合わせの対象となる領域の輪郭点を検出した結果である。輪郭点の検出方法としては様々な方法が考えられるが、例えば三次元画像の画素値の空間勾配を求め、その空間画素勾配の大きさに対して閾値処理を施すことで実現することができる。

三次元画像からの輪郭点検出の方法は、必ずしも自動的な処理の必要は無く、例えば、マウス１０５やキーボード１０６などからの入力を用いて、ユーザが手動で輪郭を入力できるようにしてもよい。また、自動で抽出した複数の輪郭点の中からユーザの指定に基づいて、後段の処理に用いる輪郭点を選択できるようにしても良い。

ここではステップＳ３０２の処理により三次元画像の表面形状からＮ個の輪郭点が抽出されたものとし、その位置座標をベクトルｘ_3di=(x_3di,y_3di,z_3di)^T，1≦i≦Nと表記する。

ステップＳ３０３において、生成部１０２０は、ステップＳ３０２で抽出した輪郭点ｘ_3diに対して陰多項式（以下Implicit Polynomialの略記号として「ＩＰ」と呼ぶ）の当てはめを行う。それにより、位置合わせの対象となる領域の形状を陰多項式であらわす。

ここで、陰多項式とは、陰関数を多項式の形式で定義したものである。より具体的には、ステップＳ３０２で抽出した位置合わせの対象となる領域の外面上のＩＰの値が、ゼロ等値面となるように近似する。これにより、位置合わせの対象となる領域からの距離に応じた値を算出する。この位置合わせの対象となる領域からの距離に応じた値の分布を距離場と呼ぶことがある。
また、陰関数を多項式の形式で表現し、その多項式の係数を求める処理を対象物体のモデル化と呼ぶ場合もある。

この処理について図５を用いて説明する。図５の輪郭点４００はステップＳ３０２で求めた輪郭点であり、ＩＰゼロ等値面４０１は、このステップで求める陰関数のゼロ等値面である。求めたい陰関数Ψ_3D(x)は、これらの点の位置において以下の制約を満たす。

ここで、ベクトルｘ_3diはステップＳ３０２で求めた複数の輪郭点のうちｉ番目の輪郭点の位置座標ベクトル(x_3di,y_3di,z_3di)^Tである。

一方、陰多項式は、陰関数Ψ_3D(x)を多項式で表現するため、以下のように定式化される。

ここでｎは多項式の次数を表している。例えばｎ＝３の場合には、

となる。

式２の右辺の第一項と第二項をそれぞれ横ベクトルa^T、縦ベクトルm(x)と表記すると、

と書き直すことができる。

さらに、陰多項式Ψⁿ _3D(x)はＮ個の輪郭点の全ての位置で式４が成り立つことから、それらの輪郭点の全てを連結した行列Ｍ=(m(x_3d0)m(x_3d1)・・・ m(x_3dN))^Tを定義することにより、次式を得る。

ここで、０は要素が全て０となるＮ次元ベクトルである。

式５から直接的にａを求めるのは解の不安定性を招くため、幾つかの制約を加えた安定的な解法が非特許文献２、非特許文献３、非特許文献４に開示されている。本実施形態においても、これらの技術を利用して実現することができる。

また、次数ｎは、表現する対象物体の表面形状の複雑さに応じて適切に選択する必要があり、例えば非特許文献５には、次数ｎを対象物体の形状に対して適応的に求める技術が開示されている。また、対象物体の形状に適した次数を予め設定しておくことも本発明の１実施形態になりうる。

以上の処理を行うことで、陰多項式の係数行列ａが求まり、ステップＳ３０１で入力した三次元画像における撮影対象物体の表面形状がモデル化される。

ここで、次数ｎは、位置合わせの精度に依存するものであり、ｎが高次になるほど位置合わせの精度があがる。一方、ｎが高次になるに従い計算処理の負荷は増すものである。

本実施形態では、複数の次元ｎに対して多項式を求めておき、使い分けるものとする。

計算の初期の段階では、荒い位置合わせで高速化するとともに、ローカルミニマムと呼ばれる極所解に入らないようにｎの次元を小さくしておく、そして計算回数が増すに従いｎの値を上げた式に変更する。これは、後述する算出部１０５０の値に応じて変更しても良いし、予め設定された計算回数に応じて使用する多項式の次数を変更しても良い。
すなわち、高精度な位置合わせをする場合にはｎの次元を上げる。

医療画像においては、撮影条件と画質に強い相関があり、例えば、Ｘ線の撮影では線量が増すほど画質があがる。そのため、高画質な画像に対しては高次のｎが要求される。つまり、撮影条件によりｎの次元を定めることが好ましい。

ステップＳ３０４において、第二の画像入力部１０３０は、第二の撮影装置１１１０が撮影した第二の画像データを入力する。この入力は第二の撮影装置１１１０の撮影と同期して直接入力することもできるし、あるいは第二の撮影装置１１１０が過去に撮影した画像を、図２に記載の画像記録装置３に記録し、その画像を読み出して入力することもできる。いずれにしても、この第二の画像はステップＳ３０１で入力された三次元画像、又は二次元画像に写る対象物体の少なくとも一部を撮影した画像である。ここで、説明を簡略化するために入力する画像データは撮影対象物体の二次元断層像であるものとして説明を行う。

ステップＳ３０５において、取得部１０４０は、ステップＳ３０４で入力された第二の画像から撮影対象物の輪郭点を検出する。ここで撮影対象とは、注目領域のことであり、例えば、医療画像であれば肺、胃、心臓、頭、手等の部位等が相当する。第二の画像が三次元画像であれば、撮影対象物の表面の輪郭点を検出する。

輪郭点は、例えばエッジ検出処理により抽出される。エッジ検出処理は、画像上の画素値とが大きく変化する位置を検出する処理である。例えば画素値の空間勾配の計算や、ラプラシアンフィルタ、ソーベルフィルタ、キャニーオペレータなどによって求めることができる。

図６は、画像に対して空間勾配の計算を行うことでエッジ検出を行った結果を模式的に表している。図６の（ａ）は入力された第二の画像を表している。（ｂ）は（ａ）の画像に対して空間勾配の絶対値を値とする画像を示している。この図において、（ａ）で撮影された撮影対象の輪郭付近において（ｂ）の値が大きくなる。以上のようにして求めたエッジに対して、予め設定した閾値処理を施すことにより、エッジが検出された画素と、エッジでない画素とを分別する。そして、検出されたＮ個の点群の夫々の点の画像座標を、ｘ_Usi=(x_USi,y_USi,z_USi)^T（iは検出された点群におけるi番目の点を表す識別子）として記憶する。

なお、本実施形態における第二の画像は二次元の断層像であるので、画像の横方向の座標値をx、縦方向の座標値をyとし、画像の厚み方向の座標値をz=0とする。
この場合、点群は三次元空間中の同一平面上に存在することになる。

ここで、第二の画像にノイズの混入がある場合には、上記エッジ検出処理の前にガウシアンフィルタやメディアンフィルタなどの平滑化フィルタを施し、ノイズの低減を行うことが望ましい。また、画像領域に被写体以外の領域が含まれる場合には、（ｃ）のように画像領域に限定して処理を施すことが望ましい。例えば、画像に対してマスク処理を施したり、エッジ検出結果にマスク処理を施したりするなどして、撮影対象以外に起因するエッジを極力排除することが望ましい。
このステップで第二の画像の表面形状上における複数の位置座標が得られる。

ここで、第一の画像の表面形状上における複数の位置座標であるｘ_3D=(x_3D,y_3D,z_3D)^Tと、第二の画像の表面形状上における複数の位置座標であるｘ_Us=(x_US,y_US,z_US)^Tとの関係を考える。これらの位置座標は基準とする座標系が異なるため、例えば同一の被検査者を撮影した夫々の画像中で人体の同一の部位を示す画素の位置座標は異なる。図７は、異なる座標系で撮影された画像同士の位置合わせを一般的に説明する図である。図７の（ａ）は、座標系６１０で対象物体を撮影した画像であり、本実施形態の第二の画像であるものとする。図７の（ｂ）は、座標系６１０とは異なる座標系６１１の元で対象物を撮影した画像であり、本実施形態の第一の画像であるものとする。ここでは書面での説明の都合で三次元画像を二次元平面で提示している。また、図７の（ｃ）は、（ａ）および（ｂ）の画像の座標系が一致していると仮定して単純に合成した時の合成画像である。実際には座標系６１０と座標系６１１は異なるので（ｃ）の合成画像では、撮影対象物体の輪郭がずれて表示されてしまう。一方、図７の（ｄ）は、座標系６１０と座標系６１１との関係を正しく反映させて、両方の画像を合成した図である。以下に説明する本実施形態のステップＳ３０６〜ステップＳ３１０の処理は、このように異なる座標系の元で撮影された複数の画像の対応関係を求める処理（位置合わせ処理）であり、この処理は算出部１０５０で行われる。

このとき、第一の画像と第二の画像における同一部位を示す画素ｘ_3Dとｘ_USの関係は、

と書く事ができる。ここで、Ｒ_3D→USは正規直行する３×３の回転行列であり、ｔ_3D→USは各軸方向への並進ベクトルでる。ここで、座標値の表記を拡張ベクトル(x,y,z)^Tにすることで、式６は

と書き直せる。ここで、Ｔ_US→3Dは次式に示す４×４の行列である。

本実施形態における位置合わせ処理とは、第二の画像の対象物体の表面上における複数の輪郭点の位置座標の座標変換方法としての変換行列Ｔ_US→3Dを決定することである。

すなわち、ステップＳ３０３で得た第一の画像の撮影対象の表面と、ステップＳ３０５で得た第二の画像の表面上における複数の位置座標との距離がなるべく近くなるような変換行列Ｔ_US→3Dを決定する。この処理を行う方法には様々なものが考えられるが、本実施形態では、比較的安定的に解が求まる繰り返し計算によって解を得る。すなわち、変換行列の推定値T~_US→3Dを逐次更新する処理を繰り返し行い、最終的に真の変換行列Ｔ_US→3Dに近い値を導出する。上述したように、位置合わせの精度に応じて用いる多項式の次数ｎも計算回数に応じて変更する。これにより位置合わせの高精度化が行えると共に、計算回数も短縮できる。また、撮影装置及び／または撮影条件に応じて、多項式の次数の変更方法を変えることが好ましい。

ステップＳ３０６において、算出部１０５０は、繰り返し計算で逐次更新していく変換行列T~_US→3Dの初期値を設定する。本実施形態では、初期値の外部位置計測センサ６の出力信号を計測結果として読み出し、その計測値から得られる変換行列をT~_US→3Dの初期値とする。これによって得られる変換行列は、センサの誤差や被検査者の体動の程度により、真の変換行列とは異なった変換行列となる。

ここで、ステップＳ３０５で検出された第二の画像の対象物の表面上における複数の位置座標がＮ個あり、そのうちｉ番目のエッジ点の位置座標がｘ_USi=(x_USi,y_USi,z_USi)^Tであるとする。ステップＳ３０７において、算出部１０５０は、このエッジ点の位置座標値を第一の画像を基準とした座標系へ座標変換する。この変換は以下の計算によって行うことができる。

ここで、T~_US→3Dは、第二の画像の座標系から第一の画像の座標系への変換行列の推定値を表す。T~_US→3Dの値は、最初はステップＳ３０６で設定した初期値を用いるが、繰り返し計算の過程では順次更新された値を使用する。

ステップＳ３０８において、算出部１０５０はさらに、陰多項式表現された第一の画像の対象物体表面に関する情報を入力し、それと第二の画像上の各エッジ点が近づくように変換行列T~_US→3Dを更新する。

これを行うために、まずステップＳ３０７で、算出部１０５０は、第一の画像の座標系に投影された第二の画像上のエッジ点ｘ_3diに対して、ＩＰにより算出されるベクトルｇ(ｘ_3Di)によって暫定的な移動先座標ｘ^’ _3Diを算出する。

ここでベクトルｇ(ｘ_3Di)は、

である。ここで∇は空間勾配を求める演算子を表している。また、dist(ｘ_3Di)はエッジ点ｘ_3DiとＩＰによって表現された物体表面との近似距離であり、

である。なお、Ψⁿ _3Dは対象物体の形状モデルを表すＩＰであるので、Ψⁿ _3D(x_3Di)は位置座標ｘ_3diにおける距離場の値を表すことになる。

次に、第一の画像の座標系に投影された第二の画像上のエッジ点ｘ_3diと、それを式１０により暫定的に移動させたエッジ点ｘ^’ _3diから、この移動を最小二乗の評価基準で近似する剛体変換のパラメータを求める。そのために、まずｘ_3diおよびｘ^’ _3diについてi=0〜Nの座標ベクトルを結合した行列Ｘと行列Ｘ’を以下のように生成する。

ここで、x-_3di，x-^’ _3di はそれぞれｘ_3di，ｘ^’ _3diの平均値である。そして、

を求め、このＡの特異値分解Ａ=ＵＳＶ^Tを行う。そして、その結果を用いて回転行列Ｒおよび、並進ベクトルｔを次のように求める。

最後に、位置合わせのための変換行列T~_US→3Dの更新を行う。式１６、式１７の結果を用いて更新後の変換行列T~_US→3Dを次のように算出する。

ステップＳ３０９において、決定部１０６０は、ステップＳ３０８で更新された変換行列T~^’ _US→3Dを評価する。そのために、ステップＳ３０５で得た点群の第一の画像の座標系への投影結果と、ステップＳ３０３で得たＩＰとの距離distを計算する。距離の計算は、

として近似的に求めることができる。ここで、x^” _3Di=T^~’ _US→3Dx_USiであり、ステップＳ３０８で求めた更新後の変換行列T~^’ _US→3Dにより第一の画像の座標系へ投影した第二の画像上のエッジ点の位置座標ベクトルである。

ステップＳ３１０において、決定部１０６０はさらに、ステップＳ３０７からステップＳ３０９までの位置合わせ処理を終了するか、さらに繰り返しを行うかの判定を行う。この判定は例えば、ステップＳ３０９で求めた各エッジ点の距離の総和と予め設定した閾値との比較により行うことができる。この比較の結果、各エッジ点の距離の総和が閾値よりも小さい場合には位置合わせ処理を終了し、そうでない場合にはステップＳ３０７に処理を戻して位置合わせ処理を繰り返すように制御される。また、それ以外にも位置合わせ処理によりステップＳ３０９で求めた各エッジ点の距離の総和の減少量が予め設定した閾値以下になった場合に位置合わせ処理の終了と判定するとしても良い。上述のように、本実施形態では各エッジ点の距離の総和に応じて、多項式の次数ｎも変更する。

ステップＳ３１０において位置合わせ処理を終了しないと判定した場合には、再びステップＳ３０７へ戻って、位置合わせ処理の繰り返しが継続される。この時、ステップＳ３０８で求めた更新後の変換行列T^~’ _US→3Dを用いてステップＳ３０７以降の処理が行われる。

ステップＳ３１１において、座標変換部１０７０は、ステップＳ３１０までで求められた変換行列T~^’ _US→3Dに基づいて、第二の画像の座標変換を行う。そして第一の画像と第二の画像の統合表示を行う。統合表示の方法には様々なものが考えられるが、第一の画像と第二の画像の画素値を予め定まる比率で加算して表示することができる。

また、第一の画像の一部の領域を、対応する第二の画像領域の面像に置き換えて表示することができる。

また、第一の画像と、座標変換後の第二の画像を、並べて表示することもできる。

以上のような本実施形態における位置合わせ処理装置１０００によれば、複数画像の高精度な位置合わせを、高速に行うことが可能となる。また、ＩＰを用いて対象物体をモデル化することにより、メモリ上の距離場にアクセスする必要がない。このため、移動量の算出を行うための演算（式１１、１２）を非常に高速に行える効果がある。また、メモリ展開した距離場との照合が不用であり、多項式の次数を変更することが可能である。これにより、ローカルミニマムと呼ばれる極所解に陥ることなく、位置合わせ精度をあげることが出来る効果を有する。

また、三次元画像の座標系からから二次元画像の座標系への変換行列を求める事にしても良い。その場合、ステップＳ３０７の処理はステップＳ３０３で求めたＩＰを移動・回転させることになる。ＩＰの移動・回転は、非特許文献６に開示されている方法により、ＩＰの係数行列に対する演算処理で実現できる。

また、外部位置計測センサ６の計測値を用いて変換行列の初期値を設定する方法を例として説明したが、本発明はそれに限られるものではない。
例えば、変換行列の初期値の設定に外部位置計測センサ６の計測値を用いずに、単位行列などの適当な値を与えて繰り返し計算の処理を開始しても良い。

その場合、真の変換行列と初期値との乖離が大きく、繰り返し計算によって真の変換行列とは異なる局所最適値に陥ってしまったり、計算の収束に非常に時間を要してしまったりする可能性がある。
それに対処するために、繰り返し処理を行った回数もカウントし、そのカウントが所定値以上に到達してもステップＳ３１０の終了条件を満たさない場合に、ステップＳ３０６で設定した初期値と個なる初期値を再度設定しなおして処理をやり直すようにしても良い。

また、繰り返し処理の回数をカウントせずとも、ステップＳ３０６で複数の初期値を設定し、その初期値毎にステップＳ３０７以後の処理を行い、得られた複数結果から最終的な結果を導く処理手順にしても良い。

また、複数の超音波画像を対象に連続して位置合わせ処理を行う場合には、直前の超音波画像に対する位置合わせ結果を次の処理の初期値として与えるようにしても良い。この方法によれば、連続する超音波画像同士の撮影位置が近い場合には、最適解に近い初期値から処理を開始できるため、繰り返し処理の回数を省略でき、処理の効率化が期待できる。

また、エッジ検出を行うことで輪郭点を抽出し、その輪郭点に対してＩＰをモデリングする処理を例として説明したが、本発明はこれに限定されるものではない。

例えば、三次元画像から撮影対象とする臓器のテクスチャなどの領域特徴を抽出し、その領域が陰多項式の内部に含まれるように陰多項式のモデリングを行うこともできる。この場合は、対象物体は臓器のテクスチャで定められることになる。

この方法によれば、三次元画像の撮影対象の輪郭部画素においてエッジが無い場合や、位置合わせ処理で容易にエッジ検出をできないような場合にも本発明を適用することができる。

また、ステップＳ３０８において、算出部１０５０が、第二の座標系におけるＩＰの距離場に基づいて各エッジ点の移動を計算し、変換行列を更する場合を例として説明したが、本発明の実施はこれに限定されない。
例えば、式１３のg(x_3Di)は、

（kはスカラーの定数）としても良い。

すなわち、式１４で示したＩＰの距離場（距離の近似）を直接用いずとも、ＩＰの距離場から算出可能な他の値を用いてエッジ点の移動や変換行列の更新を行うものも、本発明の一つの実施形態となりうる。

［その他の実施形態］
また、本発明は、前述した実施形態の機能を実現するソフトウェアのプログラムコードを記録した記憶媒体をシステム或いは装置に供給し、コンピュータがそのプログラムコードを読み出し実行することで、前述の実施形態の機能が達成される場合を含む。この場合、記録媒体から読み出されたプログラムコード自体が前述した実施形態の機能を実現することになり、そのプログラムコードを記録したコンピュータ読み取り可能な記憶媒体は本発明を構成することになる。

また、コンピュータが読み出したプログラムコードを実行することにより、前述した実施形態の機能が実現される構成に限られるものではない。例えば、そのプログラムコードの指示に基づき、コンピュータ上で稼働しているオペレーティングシステム（ＯＳ）などが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。

さらに、記憶媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張カードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれ、前述した実施形態の機能が実現される場合も含まれる。その場合、そのプログラムコードの指示に基づき、その機能拡張カードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現されることになる。

本発明を上記記録媒体に適用する場合、その記録媒体には、先に説明したフローチャートに対応するプログラムコードが格納されることになる。

なお、上述した本実施の形態における記述は、本発明に係る好適な位置合わせ処理装置の一例であり、本発明はこれに限定されるものではない。

Claims

第一の画像の第一の領域と第二の画像の第二の領域との位置合わせを行う処理装置であって、
前記第一の画像における位置合わせの対象となる前記第一の領域からの距離に応じた値を算出する多項式を生成する生成手段と、
前記第二の領域の前記第二の画像における位置座標を前記第一の画像上に座標変換した位置座標と前記多項式に基づいて前記位置合わせを行うための補正値を得る取得手段と、を備えることを特徴する処理装置。
前記第二の画像における位置合わせの対象となる第二の領域の複数の位置座標を前記第一の画像における座標に変換して得た位置座標に対して、前記第一の領域からの距離に応じた値を前記多項式を用いて夫々算出する算出手段を更に備え、
前記取得手段は、前記算出手段で算出された値に基づいて前記補正値を得ることを特徴とする請求項１に記載の処理装置。
第一の領域と第二の領域との位置合わせを行う処理装置であって、
第一の画像における位置合わせの対象となる前記第一の領域からの距離に応じた値を算出する多項式を生成する生成手段と、
第二の画像における位置合わせの対象となる前記第二の領域から得た複数の位置座標を前記第一の画像の座標に変換して得た位置座標に対して、前記第一の画像における前記第一の領域からの距離に応じた値を前記多項式を用いて夫々算出する算出手段と、
前記算出手段で算出された値に基づいて、前記複数の位置座標の座標変換方法を決定する決定手段と、
前記決定手段で決定した座標変換方法で、前記第二の画像の位置座標を座標変換する座標変換手段と、を備えることを特徴する処理装置。
前記生成手段は、計算回数に応じて前記多項式の次数を変更することを特徴とする請求項１乃至３のいずれか１項に記載の処理装置。
前記生成手段は、前記第一の画像の撮影条件で前記多項式の次数を定めることを特徴とする請求項１乃至４のいずれか１項に記載の処理装置。
前記生成手段で生成される多項式は、前記第一の画像における位置合わせの対象となる前記第一の領域の外面上の位置座標に対して０の値を示す陰関数であることを特徴とする請求項１乃至５のいずれか１項に記載の処理装置。
前記第一の画像における前記第一の領域としての対象物体、前記第二の画像における前記第二の領域としての対象物体は、それぞれ、単純Ｘ線撮影装置、Ｘ線コンピュータ断層撮影装置、核磁気共鳴映像装置、核医学診断装置、超音波画像診断装置のいずれかで撮像されたデータから構成された物体であることを特徴とする請求項１乃至６のいずれか１項に記載の処理装置。
前記第一の画像における前記第一の領域としての対象物体と前記座標変換手段で座標変換した前記第二の画像における前記第二の領域としての対象物体を重ねて表示する表示手段を更に備えることを特徴とする請求項３に記載の処理装置。
前記第二の画像は、超音波画像診断装置で撮影されたものであり、
前記超音波画像診断装置は撮影に用いる撮影プローブの位置および姿勢を計測する位置計測センサを備え、
前記座標変換手段は、前記座標変換方法の初期値を該位置計測センサの出力信号に基づいて定めることを特徴とする請求項３に記載の処理装置。
第一の画像における位置合わせの対象となる第一の領域の輪郭からの距離に応じた値を算出する多項式に、前記第一の画像とは別の第二の画像における位置合わせの対象となる第二の領域の輪郭から得た位置座標を入力して計算した値に応じて前記第二の領域を前記第一の領域に位置合わせするための位置変更量を決定する決定手段と、
前記位置変更量に基づいて前記第二の領域と前記第一の領域の位置を合わせる位置合わせ手段と、を備えることを特徴する処理装置。
第一の画像の第一の領域と第二の画像の第二の領域との位置合わせを行う処理方法であって、
前記第一の画像における位置合わせの対象となる前記第一の領域からの距離に応じた値を算出する多項式を生成する生成工程と、
前記第二の領域の前記第二の画像における位置座標を前記第一の画像上に座標変換した位置座標と前記多項式に基づいて前記位置合わせを行うための補正値を得る取得工程と、を有することを特徴する処理方法。
コンピュータを、請求項１乃至１０のいずれか１項に記載の処理装置の各手段として機能させるためのプログラム。