JP2016021998A - 画像処理装置、画像処理方法及びプログラム - Google Patents

Abstract

【課題】ユーザが入力した対応点に関する評価を提示することができるようにする。【解決手段】まず、既入力対応点群の情報に基づき、第１の画像を変形させて第２の画像に対応付ける変位場（第１の変位場）を算出する。続いて、対応点候補の組を入力すると、対応点候補の組をも考慮して第１の画像を変形させて第２の画像に対応付ける変位場（第２の変位場）を算出する。そして、二つの変位場の差異に基づいて整合性を算出し、整合か不整合かの判断結果を表示部に表示する。【選択図】図１

Description

本発明は、特に、画像間で位置を対応付けるために用いて好適な画像処理装置、画像処理方法及びプログラムに関する。

医療の分野において、複数のモダリティで撮像された医用画像や、異なる日時、異なる体位、異なる撮影モードで撮像された医用画像を診断に利用するためには、被検体内の注目点を夫々の医用画像間で対応付けることが重要である。この注目点には、画像処理により自動的に対応付けが可能なものもあるが、モダリティの違いや被検体の変形等の影響により自動的に対応付けが困難なものも存在する。そのため、一般的に、医師等は画像を見ながら対応付け作業を目視で行っている。

また、同種の目的のために、複数の画像を位置合わせして表示することが提案されている。例えば、非特許文献１には、画像間における被検体の変形を推定する変形位置合わせ技術が開示されている。ここで、より高精度な位置合わせを実現するためには、ユーザが目視で対応付けした対応点を取得して、この対応点を位置合わせ処理に用いることがしばしばなされている。したがって、対応付けの精度によって位置合わせの精度が異なってくる。

そこで、一方の画像（以下、参照画像）上で指摘した注目点を他方の画像（以下、対象画像）から探索し、対応付けする作業を支援するための技術が開示されている。例えば、特許文献１には、注目点の推定位置が対象画像上に重畳表示される技術が開示されている。また、特許文献２には、注目点の存在確率の分布が対象画像上に重畳表示される技術が開示されている。

特開２００８−２１２６８０号公報 特開２０１３−１９８７２２号公報

しかしながら、上記特許文献１及び２に開示されている技術では、ユーザが入力した対応点の良し悪しが分からないという課題がある。

本発明は前述の問題点に鑑み、ユーザが入力した対応点に関する評価を提示することができるようにすることを目的としている。

本発明に係る画像処理装置は、第１の画像と第２の画像との間の対応情報を取得する対応情報取得手段と、前記第１の画像上における第１の点の座標と前記第２の画像上における第２の点の座標とからなる対応点候補の組を取得する対応点候補取得手段と、前記対応情報取得手段によって取得された対応情報に対する、前記対応点候補取得手段によって取得された対応点候補の組の整合度を算出する整合度算出手段と、前記整合度算出手段の算出結果を表示部に表示させる表示制御手段と、を有することを特徴とする。

本発明によれば、ユーザが入力した対応点に関する評価を提示することができる。

第１及び第２の実施形態に係る画像処理装置の機能構成例を示すブロック図である。 整合度の算出方法及び判定結果の表示画面を説明するための図である。 第１の実施形態において、画像処理装置が行う全体の処理手順の一例を示すフローチャートである。 渦状の変位場の一例を示す図である。 第２の実施形態において、画像処理装置が行う全体の処理手順の一例を示すフローチャートである。 第２の実施形態において、整合度画像を生成する詳細な処理手順の一例を示すフローチャートである。 整合度マップを重畳表示した表示画面の一例を示す図である。 第３の実施形態に係る画像処理装置の機能構成例を示すブロック図である。 第３の実施形態において、画像処理装置が行う全体の処理手順の一例を示すフローチャートである。

以下、添付図面に従って本発明に係る画像処理装置の好ましい実施形態について説明する。なお、本発明は、添付図面に例示された例に限定されるものではない。

（第１の実施形態）
本実施形態に係る画像処理装置は、ユーザが指定した画像間の対応点候補の組について、その２つの点が対応点の組であることの整合度を導出して表示する。ここで、「対応点の組」とは、「第１の画像上の注目点」と、「該注目点と第２の画像上において対応する点」との２点のことである。また「対応点候補の組」とは、対応点の組の候補としてユーザが指定した２点のことである。

このように対応点候補の組の整合度を導出して表示するため、ユーザは、対応点候補の組を指定するたびに、指定した対応点候補の組の整合度を確認することができる。これにより、ユーザが指定した対応点候補の組が不整合と判定された場合は、対応点候補の指定をやり直し、より適切な位置に対応点候補の組を入力することができる。なお、本実施形態において、対応点候補の組の整合度は、第１の画像と第２の画像との間における既知の対応情報に対する整合度として算出される。また、本実施形態における既知の対応情報とは、既に対応付けされている複数の対応点の組（以下、既入力対応点群と呼ぶ）の情報である。以下、図１乃至図３を用いて本実施形態の構成および処理を説明する。

図１は、本実施形態に係る画像処理装置１００の機能構成例を示すブロック図である。
図１に示すように、本実施形態における画像処理装置１００は、データサーバ１７０、操作部１８０、及び表示部１９０と接続されており、画像処理システムを構成している。

データサーバ１７０は、ネットワークを介して画像処理装置１００と接続されており、以下に説明する第１の画像及び第２の画像を保持している。データサーバ１７０が保持する第１の画像及び第２の画像は、それぞれ異なる条件（異なるモダリティ、撮影モード、日時、体位等）で被検体を予め撮像して得られた３次元断層画像である。３次元断層画像を生成するモダリティは、ＭＲＩ装置、Ｘ線ＣＴ装置、３次元超音波撮影装置、光音響トモグラフィ装置、ＰＥＴ／ＳＰＥＣＴ、ＯＣＴ装置などである。第１の画像及び第２の画像は、例えば、異なるモダリティや異なる撮影モードで同時期に撮像されたものであってもよく、経過観察のために同一患者を同一モダリティ、同一体位で異なる日時に撮像した画像であってもよい。なお、第１の画像及び第２の画像は２次元断層画像の集合として構成されており、各２次元断層画像の位置および姿勢は、基準座標系（被検体を基準とした空間中の座標系）に変換した上でデータサーバ１７０に保持されているものとする。基準座標系で表現された第１の画像及び第２の画像は、データ取得部１０２により画像処理装置１００に入力される。

さらに、データサーバ１７０は、第１の画像及び第２の画像における既入力対応点群の情報として、夫々の対応点の組の３次元座標値の情報を保持している。これらのデータも、データ取得部１０２により画像処理装置１００に入力される。既入力対応点群の中で、互いに対応付けられている第１の画像上の点と第２の画像上の点とには、共通するＩＤ番号が付与されているものとする。すなわち、既入力対応点群の各点は、点のＩＤ番号と、基準座標系における３次元座標（３次元空間におけるＸ，Ｙ，Ｚの３つの位置成分）とによって表されているものとする。

操作部１８０は、ユーザによるマウスやキーボードの操作を受け付ける。この操作により、ユーザが指定する対応点候補の組が、対応点候補取得部１０４により画像処理装置１００に入力される。また、ユーザは、操作部１８０から対応点の指定処理を終了するか否かを指示することができる。この対応点の指定処理を終了するか否かの判定結果は、終了判定部１１２により画像処理装置１００に入力される。

表示部１９０は、画像処理装置１００が生成する表示画像や整合度情報等の各種の情報を表示する。また、表示部１９０には、ユーザからの指示を取得するためのＧＵＩも配置されている。

画像処理装置１００は、以下に説明する構成要素により構成されている。
データ取得部１０２（対応情報取得手段）は、データサーバ１７０から第１の画像及び第２の画像と、既入力対応点群の情報とを取得する。そして、既入力対応点群の情報を変位場算出部１０６へ出力し、第１の画像及び第２の画像を表示制御部１１０へ出力する。記憶部１１１は、表示部１９０に表示するための対応点の組の情報などの情報を記憶する。

対応点候補取得部１０４（対応点候補取得手段）は、ユーザによる操作部１８０の操作に従って、第１の画像上と第２の画像上とにおける対応点候補の組（３次元座標の組）を取得して、変位場算出部１０６へ出力する。変位場算出部１０６（変位場算出手段）は、第１の画像を変形させて第２の画像に位置合わせするための変位場として、既入力対応点群のみに基づく第１の変位場と、対応点候補の組をも考慮した第２の変位場とを算出する。そして、算出した変位場を整合度算出部１０８へ出力する。

整合度算出部１０８（整合度算出手段）は、変位場算出部１０６から第１の変位場及び第２の変位場を取得し、これらの情報に基づいて既知の対応情報に対する対応点候補の組の整合度を算出する。そして、算出した整合度を表示制御部１１０へ出力する。表示制御部１１０は、第１の画像及び第２の画像の断面画像と、整合度に関する情報とを、表示部１９０に表示させる制御を行う。終了判定部１１２は、ユーザによる操作部１８０の操作に従って、対応付けを終了するか否かを判定する。

図３は、本実施形態において、画像処理装置１００が行う全体の処理手順の一例を示すフローチャートである。以下、本実施形態における処理の詳細について説明する。

（ステップＳ３００：データの取得）
まず、ステップＳ３００において、データ取得部１０２は、データサーバ１７０から、第１の画像、第２の画像、および既入力対応点群の情報を取得する。

（ステップＳ３０１：第１の変位場を算出）
次に、ステップＳ３０１において、変位場算出部１０６は、ステップＳ３００で得た既入力対応点群の情報に基づき、第１の画像を変形させて第２の画像に対応付ける変位場（第１の変位場）を算出する。この第１の変位場は、第１の画像上の既入力対応点群の位置を、対応する第２の画像上の既入力点群の位置へと極力近づくように変位させ、かつ空間的に滑らかに変位を補間した場である。なお、変位場の算出は、具体的にはＢスプライン関数の補間によって変形を表すＦｒｅｅ−ｆｏｒｍ ｄｅｆｏｒｍａｔｉｏｎ（ＦＦＤ）等の公知の手法を用いて行うものとする。本実施形態では、ＦＦＤが表す変位の場を第１の変位場とする。

（ステップＳ３０２：第１の点を取得）
次に、ステップＳ３０２において、対応点候補取得部１０４は、ユーザによる操作部１８０の操作に従って、第１の画像における注目点（第１の点）の位置（３次元座標）を取得する。すなわち、対応点候補取得部１０４は、ユーザによる操作部１８０の操作がなされるまで待機し、ユーザによる操作部１８０の操作がなされると、第１の画像における注目点の位置を取得する。具体的には、まず、表示制御部１１０は、第１の画像の断面画像を表示部１９０に表示させる。ユーザはマウスやキーボードを操作して、着目する断面画像を選択して表示させる。そして、表示した断面画像において、注目位置を注目点として入力する。対応点候補取得部１０４は、断面画像上で指定された注目位置を、第１の画像における断層画像の位置および姿勢を利用して、基準座標系における注目点の位置に変換する。なお、インタレストオペレータ等の画像特徴点抽出処理によって第１の画像から注目点候補を抽出し、その中からユーザが注目点を選択する構成であってもよい。

（ステップＳ３０３：第２の点を取得）
次に、ステップＳ３０３において、対応点候補取得部１０４は、ユーザによる操作部１８０の操作に従って、ステップＳ３０２で取得した注目点（第１の点）と第２の画像上において対応する点の候補（第２の点）の位置（３次元座標）を取得する。具体的な取得方法は、ステップＳ３０２と同様である。そして、取得した第１の点と第２の点とを対応点候補の組とする。

（ステップＳ３０４：第２の変位場を算出）
次に、ステップＳ３０４において、変位場算出部１０６は、第１の画像を変形させて第２の画像に対応付ける変位場（第２の変位場）を算出する。第２の変位場の算出は、ステップＳ３０１で実行した第１の変位場の算出と同様の処理により実行する。ただし、第２の変位場の算出では、既入力対応点群に加え、ステップＳ３０２及びステップＳ３０３で取得した対応点候補の組も拘束条件に含める。

（ステップＳ３０５：整合度を算出）
次に、ステップＳ３０５において、整合度算出部１０８は、ステップＳ３０１で算出された第１の変位場とステップＳ３０４で算出された第２の変位場とを用いて、既知の対応情報に対する対応点候補の組の整合度を算出する。ここで、整合度とは、対応付けによって生じる変位場の局所的な歪みを定量化したものである。例えば、二つの変位場の差異に基づいて整合性を算出する。すなわち、対応付けを行ったことで、その位置における変位場の向きが大きく変化した場合に、整合度は低いと評価する。その反対に、変位場の向きの変化が小さかった場合には、整合度は高いと評価する。本実施形態では、対応付けを行った位置（第２の点の位置）における変位場の角度変化に基づいて、整合度を定式化する。

ここで、図２を参照しながら整合度の算出方法について説明する。図２（ａ）には、第１の断面画像２０１と第２の断面画像２０２との一例を示す。第１の断面画像２０１は、３次元画像データである第１の画像を、任意の平面（第１の切断面）で切断し、３次元画像データをその切断面上に切り出した２次元画像である。同様に第２の断面画像２０２は、第２の画像を第２の切断面で切断し、３次元画像データをその切断面上に切り出した２次元画像である。さらに、各断面画像において、既入力対応点群および対応点候補の組（第１の点及び第２の点）が示されている。図２（ａ）において、○印で示されている点が既入力対応点群であり、正確には、既入力対応点群のうちの断面画像上に存在する点である。また、第１の断面画像２０１上の×印で示されている点が第１の点である。また、第２の断面画像２０２上の×印で示されている点が第２の点である。

図２（ｂ）には、第２の断面画像２０２上における第１の変位場２０４および第２の変位場２０５を示す。ここで、第１の変位場をＴ₁と表し、第２の変位場をＴ₂と表す。第１の変位場Ｔ₁および第２の変位場Ｔ₂は、第２の画像の座標系（基準座標系）が表す空間上で定義される３次元の値を持つ関数である。各関数は、第２の画像の座標系から第１の画像の座標系へ変換するための関数を意味しており、第２の画像の座標系における任意の位置を引数とし、第１の画像の座標系における対応する位置への変位（３次元の変位ベクトル）を返す関数である。

図２（ｃ）には、第２の点の位置ｐにおける第１の変位場Ｔ₁と第２の変位場Ｔ₂との変位ベクトル（それぞれｔ_p,1、ｔ_p,2とする）の例を示す。このとき、整合度算出部１０８は、対応点候補の組の整合度Ｍを、以下の式（１）を用いて算出する。
Ｍ ＝ ｔ_p,1・ｔ_p,2 ／｜ｔ_p,1｜｜ｔ_p,2｜ ・・・（１）

式（１）に示すように整合度Ｍの計算は、２つのベクトルｔ_p,1、ｔ_p,2がなす角度の余弦を計算することを意味する。したがって、２つのベクトルがなす角度が小さい場合には整合度Ｍの値は高くなり、２つのベクトルのなす角度が大きい場合には整合度Ｍの値は低くなる。

なお、第１の変位場Ｔ₁と第２の変位場Ｔ₂との差異に基づいて他の方法により整合度を算出してもよい。例えば、位置ｐにおける第１の変位場Ｔ₁と第２の変位場Ｔ₂とにおけるヤコビ行列の行列式（ヤコビアン）を夫々求め、その差分に基づいて整合度を算出するようにしてもよい。この場合、局所領域の体積が変化した場合には整合度が低く判定するようにできる。なお、位置ｐの近傍の各点においても同様に、式（１）の値またはヤコビアンの値を求め、その積分値を整合度としてもよい。他にも、非特許文献（Ｄ．Ｒｕｅｃｋｅｒｔ, Ｌ．Ｓｏｎｏｄａ, Ｃ．Ｈａｙｅｓ, Ｄ．Ｈｉｌｌ, Ｍ．Ｌｅａｃｈ，ａｎｄ Ｄ．Ｈａｗｋｅｓ，"Ｎｏｎｒｉｇｉｄ ｒｅｇｉｓｔｒａｔｉｏｎ ｕｓｉｎｇ ｆｒｅｅ−ｆｏｒｍ ｄｅｆｏｒｍａｔｉｏｎｓ： ａｐｐｌｉｃａｔｉｏｎ ｔｏ ｂｒｅａｓｔ ＭＲ ｉｍａｇｅｓ"，ＩＥＥＥ ｍｅｄ． ｉｍａｇ．，ｖｏｌ．１８（８），ｐｐ．７１２−７２１，１９９９．）に開示されている変位場の滑らかさの算出方法を用いて、位置ｐの近傍における第１の変位場と第２の変位場との夫々の変位場の滑らかさを求め、それらの差分（すなわち、滑らかさの悪化度）に基づいて整合度を算出してもよい。なお、これらの方法は一例に過ぎず、変位場を入力として計算される定量値であれば、いかなる方法で整合度を算出しても良い。

（ステップＳ３０６：整合度を判定）
次に、ステップＳ３０６において、整合度算出部１０８は、ステップＳ３０５で算出した整合度の値が一定値以上か否かの判定を行う。例えば式（１）によって算出される整合度Ｍを用いる場合に、その値が０．５以上の場合、すなわち、２つのベクトルｔ_p,1、ｔ_p,2のなす角度が６０度以下の場合は整合と判定する。また、整合度Ｍの値が０．５未満の場合は不整合と判定する。

（ステップＳ３０７：追加指示の判定）
次に、ステップＳ３０７において、表示制御部１１０は、対応点候補の組の整合度の算出結果の情報として、ステップＳ３０５で算出された整合度の値、及びステップＳ３０６の判定結果を表示部１９０に表示させる。そして、対応点候補取得部１０４は、ユーザにより操作部１８０が操作され、対応点候補の組を対応点の組として採用する指示を受けたか否か判定する。なお、この表示は、ステップＳ３０６の処理で不整合と判定された場合にのみ実施するようにしてもよく、整合と判定された場合にはユーザからの指示を待たずにステップＳ３０８に処理を進めるようにしてもよい。

図２（ｄ）には、判定結果が表示された状態で採用／不採用を選択させるＧＵＩの例を示す。図２（ｄ）に示す例では、整合度の値、判定結果、及び採用か不採用かを選択させる２つのボタンが配置されたダイアログボックス２１０を表示する。そして、対応点候補取得部１０４は、ユーザの操作によりいずれかのボタンをマウスでクリックするなどして入力された指示を取得する。ステップＳ３０７の判定の結果、採用が指示された場合は、ステップＳ３０８へと処理を進め、不採用が指示された場合は、ステップＳ３０９へと処理を進める。なお、ダイアログボックス２１０には、判定結果を表示せずに整合度の値のみを表示させるようにしてもよいし、整合度の値は表示せずに判定結果のみを表示させるようにしてもよい。

（ステップＳ３０８：既入力対応点群に追加）
ステップＳ３０８においては、対応点候補取得部１０４は、対応点候補の組を対応点の組として既入力対応点群に追加する。

（ステップＳ３０９：第２の点を削除）
ステップＳ３０９において、対応点候補取得部１０４は、ステップＳ３０３で取得した第２の点の情報を削除する。これにより、第１の点のみが取得された状態で、第２の点は未取得の状態となる。そして、ステップＳ３０３以降の処理を再度実行する。

（ステップＳ３１０：終了判定）
ステップＳ３１０においては、終了判定部１１２は、対応点候補の指定処理を終了する指示を受けたか否か判定する。例えば、表示制御部１１０が表示画面に終了するか否かを選択させる終了ボタンを配置し、終了判定部１１２は、ユーザの操作により終了ボタンをマウスでクリックするなどして終了の指示を受けたか否かを判定する。ステップＳ３１０の判定の結果、終了する指示を受けた場合にはステップＳ３１１へと処理を進める。一方、終了する指示を受けていない場合には、ステップＳ３０１へ戻り、更新された既入力対応点群を用いてステップＳ３０１以降の処理を再度実行する。ここで、再指定の入力がない場合には、終了の入力があるまで断面画像が表示された状態で待機状態となり、入力を受け付ける。

（ステップＳ３１１：対応点群の出力）
ステップＳ３１１においては、データ取得部１０２は、以上の処理で新たに取得した２つの３次元画像データ間の対応点の組の情報を、記憶部１１１に格納する。さらに、表示制御部１１０は、表示部１９０に対してそれらの対応点の情報を出力し、表示部１９０にその旨の情報を表示させる。このとき、データサーバ１７０に格納されている既入力対応点群の情報を更新するよう、更新された既入力対応点群の情報をデータサーバ１７０に不図示の通信部から送信するようにしてもよい。もちろん、最終的に取得した既入力対応点群に基づいて推定した変位場を記憶部１１１やデータサーバ１７０に格納するようにしてもよい。また、該変位場を用いて、第２の画像に対応付くように第１の画像を変形させた変形画像を生成して、これを記憶部１１１やデータサーバ１７０に格納するようにしてもよい。

以上のように本実施形態によれば、第１の画像と第２の画像との間の対応情報を取得する対応情報取得手段（データ取得部１０２）と、第１の画像上における第１の点の座標と第２の画像上における第２の点の座標とからなる対応点候補の組を取得する対応点候補取得手段（対応点候補取得部１０４）と、対応情報取得手段（データ取得部１０２）によって取得された対応情報に対する、対応点候補取得手段（対応点候補取得部１０４）によって取得された対応点候補の組の整合度を算出する整合度算出手段（整合度算出部１０８）と、整合度算出手段の算出結果を表示部に表示させる表示制御手段（表示制御部１１０）を備える。よって、ユーザに先入観を与えることなく、指定した対応点候補の組が既入力対応点群に整合しているかどうかを確認する機能が提供される。このため、不整合と判定された場合は対応点候補となる点の入力をやり直し、より適切な位置に指定点を入力することができる。

＜変形例１＞
本実施形態では、既知の対応情報として既入力対応点群を用いる例について説明したが、他の既知の対応情報を用いる構成であってもよい。例えば、他の手段で取得した変位場そのものを既知の対応情報として用いてもよい。ここで、他の変位場は、例えば、有限要素モデルを用いた物理変形シミュレーション結果や、解剖学的知識に基づいて得ることができる。

この場合、ステップＳ３００において、データ取得部１０２は、既知の対応情報として前述の他の変位場を取得する。そして、ステップＳ３０１において、変位場算出部１０６は、ステップＳ３００で取得した変位場を第１の変位場とする。また、ステップＳ３０４において、変位場算出部１０６は、既入力対応点群に相当する情報を第１の変位場から導出し、この導出した既入力対応点群に相当する情報に基づいて第２の変位場を算出する。具体的にはまず、第２の画像の定義範囲を一定間隔でサンプリングし、第２の画像上に格子状に配置された点群（第２の画像上の点群）を得る。次に、第１の変位場を用いて、前記第２の画像上の点群に対応する第１の画像上の位置を全ての点について求め、第１の画像上の点群とする。そして、この第１の画像及び第２の画像上の点群の組を、既入力対応点群とする。

このように変形例１によれば、既入力対応点群が得られない場合であっても整合度を算出することができる。例えば画像が不明瞭であり、第１の画像と第２の画像とで確信を持って対応していることが分かる点がない場合でも、対応付けを支援することができる。

＜変形例２＞
本実施形態では、整合度算出部１０８の処理として、第１の変位場と第２の変位場とに基づいて整合度を算出する例について説明した。一方、第１の変位場と対応点候補の組の情報とから、第２の変位場を算出せずに整合度を算出してもよい。この場合、例えば整合度Ｍを、以下の式（２）を用いて算出する。
Ｍ ＝ ｔ_p,1・ｔ_p,q ／｜ｔ_p,1｜｜ｔ_p,q｜ ・・・（２）

ここで、ｔ_p,qは、第２の点の位置ｐと第１の点の位置ｑとを結ぶベクトルを表す。式（２）に示すようにこの場合は、第１の変位場における変位ベクトルｔ_p,1と、対応点候補の組を結ぶベクトルｔ_p,qとの角度差の余弦を算出することを意味する。他にも、２つのベクトルの差分のノルム、すなわち、以下の式（３）により整合度Ｍを算出してもよい。この場合、第１の変位場によって推定される位置と、実際に入力された第１の点の位置ｑとの距離が整合度として算出される。
Ｍ ＝ ｜ｔ_p,1−ｔ_p,q ｜ ・・・（３）

以上のように変形例２によれば、３次元画像データ領域全体にわたって第２の変位場を算出することなく整合度を算出することが可能であるため、計算量を減少させることができる。

＜変形例３＞
本実施形態では、整合度算出部１０８の処理として、第１の変位場と第２の変位場との差異に基づいて整合度を算出する例について説明した。一方、第１の変位場を用いずに第２の変位場から整合度を算出してもよい。例えば、第２の変位場によって既入力対応点を一方から他方へ変位させた際の残差を求め、残差の平均値の逆数を整合度としてもよい。また、第２の変位場における変位場の乱れを計量し、この変位場の乱れを整合度としてもよい。例えば、非特許文献１に開示されている変位場の滑らかさの算出手法を用いて、位置ｐの近傍における第２の変位場の滑らかさを求め、その絶対値を整合度としてもよい。他にも、変位場のヤコビ行列の行列式（ヤコビアン）を整合度としてもよい。この場合、変位場の中でヤコビアンの値の絶対値が大きい位置は、局所的に大きな変位を表示しているため、整合度が低いと判定される。特に、ヤコビアンの値が負の場合、空間が折りたたまれるような変位を生じていることを意味するため、不整合の可能性が高い。また、第２の変位場が渦状になっている度合いなどを計量し、これを対応点候補の組の整合度としてもよい。

ここで、渦状の変位場が発生する場合の例について、図４を用いて説明する。図４（ａ）における格子の角は、第２の画像におけるボクセル位置４０１を表す。また、第１の画像上の既入力対応点４０２及び第２の画像上の対応する既入力対応点４０３が与えられている状況である。そして、対応点候補取得部１０４により、既入力対応点の近傍において、既入力対応点とは逆向きの変位を持った第２の点４０４と第１の点４０５とからなる対応点候補の組を取得した場合を考える。図４（ｂ）には、このような状況下で算出される、第２の点４０４の近傍領域における第２の変位場４０６の例を示す。図４（ｂ）に示すように、近傍に大きく方向が異なる２つの変位ベクトルが存在するため、算出される変位場は、第２の点４０４付近を中心とした回転を表す渦状の形となる。このような変位場が生成される場合は、整合度は低いと考えられる。

以下に、図４（ｂ）に示した渦状の変位場を検出するための具体的な方法を示す。まず、第２の点４０４の位置ｐから距離値が一定範囲内の３次元領域を、変位場の渦を検出する範囲４０７として定義する。そして、この範囲に含まれるボクセル位置の数をＮとする。また、第２の点４０４の位置ｐにおける変位ベクトルをｔ₁とし、範囲４０７に含まれる第２の点４０４以外のボクセル位置における変位ベクトルを、それぞれｔ₂，ｔ₃，…ｔ_Nとする。このとき、範囲４０７に含まれる変位場の渦度Ｕは、以下の式（４）で定義される。

式（４）に示すようにこの場合は、範囲４０７に含まれるすべてのボクセル位置における変位ベクトルについて、変位ベクトルｔ₁との角度差の余弦を計算し、その角度差の余弦の平均値を求めることを意味する。渦度Ｕが１に近い場合、すなわち、範囲４０７内の変位場の多くが変位ベクトルｔ₁に近い方向である場合は整合度が高いと考えられる。一方、渦度Ｕが−１に近づくほど整合度が低くなると考えられる。そして、例えば、渦度Ｕが０．５以上であれば整合と判定し、渦度Ｕが０．５未満であれば不整合と判定することができる。

以上のように変形例３によれば、整合度が低くなる場合に変位場に現れるパターンに対する知識を整合度の算出に用いることができる。

（第２の実施形態）
本実施形態に係る画像処理装置は、第１の画像で入力された注目点との整合度の分布を表示することにより、対応点の入力を支援する。本実施形態に係る画像処理システムの構成は図１と同様である。ただし、以下に説明するように、各部の処理内容の一部は第１の実施形態とは異なっている。以下、本実施形態に係る画像処理装置１００について、第１の実施形態との相違部分についてのみ説明する。

図５は、本実施形態において、画像処理装置１００が行う全体の処理手順の一例を示すフローチャートである。なお、ステップＳ５００からＳ５０２までの処理は、それぞれ図３のステップＳ３００からＳ３０２の処理と同様である。また、ステップＳ５０５からＳ５０７までの処理は、それぞれ図３のステップＳ３０８、Ｓ３１０及びＳ３１１の処理と同様であるため、説明を省略する。

（ステップＳ５０３：整合度画像を生成）
ステップＳ５０３において、画像処理装置１００は、第２の画像上の各ボクセル位置について、当該位置を注目点に対応する点の候補と仮定した際の、既知の対応情報との整合度を算出する。すなわち、それぞれのボクセル位置を、第１の点に対応する第２の点と仮定した場合（対応点候補の組と仮定した場合）の既知の対応情報との整合度を算出する。そして、各ボクセル位置における整合度をボクセル値として有するボリューム画像（以下、整合度画像）を生成する。以下、図６のフローチャートを用いて、ステップＳ５０３の処理の詳細について説明する。

図６は、図５のステップＳ５０３の詳細な処理手順の一例を示すフローチャートである。
（ステップＳ６００：整合度画像を初期化）
まず、ステップＳ６００において、整合度算出部１０８は、第２の画像と同サイズの整合度画像を生成して初期化する。このとき、全てのボクセル位置に０を与える。

（ステップＳ６０１：画像中のボクセル位置を取得）
次に、ステップＳ６０１において、対応点候補取得部１０４は、第２の画像中のある１つのボクセル位置の３次元座標値を取得し、データ取得部１０２は、取得した３次元座標値のインデックスを記憶部１１１に格納する。ボクセル位置の３次元座標値を取得する際には、記憶部１１１に格納されたインデックスを検索し、まだ取得されていないボクセル位置のみを対象とする。取得された３次元座標値は、第１の実施形態における第２の点と同じ役割を持つ。すなわち、ステップＳ５０２で取得した第１の点と組み合わせたものを、対応点候補の組とする。

（ステップＳ６０２：第２の変位場を算出）
次に、ステップＳ６０２において、変位場算出部１０６は、ステップＳ３０４と同様の処理により、既知の対応情報と対応点候補の組とに基づいて第２の変位場を算出する。

（ステップＳ６０３：整合度を算出）
次に、ステップＳ６０３において、整合度算出部１０８は、ステップＳ６０２で取得した第２の変位場に基づいて、既知の対応情報に対する対応点候補の組の整合度を算出する。本処理は、図３のステップＳ３０５と同様の処理によって実施できる。そして、算出した値を、ステップＳ６０１で取得したボクセル位置に対応する整合度画像のボクセルに、該位置における整合度の値として設定する。

（ステップＳ６０４：整合度を判定）
次に、ステップＳ６０４において、整合度算出部１０８は、ステップＳ３０６と同様の処理により、算出した整合度の値が一定値以上か否かを判定する。そして、不整合と判定した場合には、ステップＳ６０１で取得したボクセル位置に対応する整合度画像のボクセルに、「不整合」を表す特殊な値（例えば、負の最大値）を設定する。

（ステップＳ６０５：全ボクセル位置で整合度の算出の完了の判定）
次に、ステップＳ６０５において、整合度算出部１０８は、記憶部１１１に格納されたインデックスを参照し、整合度が算出されていないボクセル位置が存在するか否かを判定する。この判定の結果、全てのボクセル位置について整合度の算出が終了している場合はステップＳ５０３の処理を終了する。一方、まだ整合度が算出されていないボクセル位置が存在する場合は、ステップＳ６０１へ戻り、残されたボクセル位置に対して同様の処理を実行する。以上の手順により、整合度画像が生成される。

（ステップＳ５０４及びＳ５０５：第２の点を取得、及び対応点の組の追加）
図５の説明に戻り、続いてステップＳ５０４において、対応点候補取得部１０４は、ユーザによる操作部１８０の操作に従って、第２の点の位置を取得する。そして、ステップＳ５０５において、対応点候補取得部１０４は、ステップＳ５０２で取得した第１の点とステップＳ５０４で取得した第２の点とを新たな対応点の組として既入力対応点群に追加する。

なお、ステップＳ５０４において第２の点を取得する具体的な方法は、基本的には図３のステップＳ３０３の処理と同様であるが、以下の点が異なっている。まず、表示制御部１１０は、第２の点の位置を指定するための第２の画像の断面画像を表示するとともに、整合度画像において該断面画像と対応する断面の整合度を得て、その分布の情報を表示部１９０に表示させる。この表示は、例えば、濃淡で表現した第２の画像の断面画像上に整合度の分布をカラーマップ（整合度マップ）として重畳表示するような形態とすることができる。例えば、整合度の値と表示色との対応を予めテーブルで定義しておき、そのテーブルに従って決定した表示色を第２の画像の対応するボクセル位置に重畳する。この時、「不整合」と判定された位置には、何も重畳しないようにする。なお、整合度の分布の情報の表示・非表示は、ユーザの指示に応じて選択可能であることが望ましい。

また、整合度の分布の表示形態は、整合度を表すカラーマップの重畳表示に限らず、例えば、整合度の等値線を重畳表示してもよい。また、第２の画像の断面画像上でユーザの操作によるカーソルの動きに連動して、カーソルの示す座標における整合度を整合度画像から取得し、その整合度を文字情報として画面内の所定の位置やカーソル付近に表示する形態でもよい。また、整合度の分布の表示形態は、第２の画像へ重畳表示する形態に限定されるものではなく、第２の画像と対比しやすいように整合度画像の同一断面を並べて表示する構成であってもよい。これら整合度の分布の表示形態は、ユーザの指示に応じて選択可能であることが望ましい。

図７には、整合度の分布の表示例を示す。図７（ａ）は、第１の画像上の注目点７０１及び第２の画像上で定義される整合度マップ７０２を示す図である。図７（ｂ）は、整合度マップ７０２を第２の断面画像２０２上に重畳表示した結果７０３を示す図である。また、図７（ｃ）は、第２の断面画像２０２上にマウスカーソル７０４が移動している際に、その座標を第２の点としたときの整合度を示す文字情報７０５をカーソル右下の位置に表示している例を示す図である。これらの提示によって、ユーザは、表示部１９０に表示されている第２の画像と整合度の分布との表示を確認しながら、整合度が高く、かつ、対応すると思われる位置に第２の点を入力することができる。

以上のように本実施形態によれば、ユーザは、整合度の分布を確認することができるため、整合度と画像上の対応関係のバランスとを考慮しながら対応付けを容易に行うことができる。また、対応点候補の組が変位場に与える影響に基づいて整合性を求めているため、信頼性が高い整合度が得られるという効果がある。

＜変形例＞
本実施形態では、ユーザの操作により第２の点が入力される前に整合度の分布を表示した。一方、ステップＳ５０４において、最初は整合度の分布を表示せず、第２の点を取得した後でユーザの指示に応じて整合度の分布を表示する構成であってもよい。このような構成とすることにより、対応点に関する先入観をユーザに与えることをより回避することができる。

また、第２の点を取得した後に、取得した対応点候補の組が不整合と判定された場合にのみ整合度の分布を表示する構成であってもよい。この場合、ステップＳ５０４の後に、整合度算出部１０８は、取得した第２の点の位置における整合度を整合度画像から取得し、その整合度に基づいて整合度の分布を表示するか否かを判定する。すなわち、不整合と判定された場合に、表示制御部１１０は、整合度の判定結果と共に整合度マップを第２の画像に重畳表示する。そして、対応点候補取得部１０４は、ユーザがこの第２の点を採用するか否かの指示を待つことになる。ユーザの操作により不採用と指示された場合には、図３のステップＳ３０９と同様の処理により第２の点を削除して、ステップＳ５０４以降の処理を再度実行する。このように、不整合と判定されたときにのみ整合度の分布が表示されるので、ユーザを煩わさずに、必要な場合には警告できるという効果がある。

（第３の実施形態）
第１及び第２の実施形態では、ユーザが対話的操作によって１点ずつ整合度を確認しながら対応付けを行う例について説明したが、これらの処理を繰り返すことによって、既入力対応点群の整合度が変化する可能性がある。本実施形態では、既入力対応点群の整合度を算出する例について説明する。以下、本実施形態に係る画像処理装置について、第１の実施形態との相違部分についてのみ説明する。

図８は、本実施形態に係る画像処理装置８００の機能構成例を示すブロック図である。なお、図１と同じ構成については、同じ符号を付けており、その説明は省略する。

対応点候補取得部８０４は、データ取得部１０２から既入力対応点群を取得し、その中のある１組を整合度の評価対象である対応点候補の組として順次選択する。また、前記既入力対応点群から前記対応点候補の組を除外した点群を、変位場算出点群として取得する。そして、変位場算出点群を、変位場算出部１０６へ出力する。

図９は、画像処理装置８００が行う全体の処理手順を示すフローチャートである。なお、ステップＳ９００の処理は、図３のステップＳ３００の処理と同様であるため、説明は省略する。

（ステップＳ９０１：第３の変位場を算出）
次に、ステップＳ９０１において、変位場算出部１０６は、データ取得部１０２により取得された既入力対応点群の全ての点を用いて、第３の変位場を算出する。ここで第３の変位場は、第１の実施形態で説明した第１の変位場と同様であり、図３のステップＳ３０１と同様の方法により算出することができる。

（ステップＳ９０２：対応点候補の組を選択）
次に、ステップＳ９０２において、対応点候補取得部８０４は、既入力対応点群の中のある１組の点を対応点候補の組として選択する。そして、前記対応点候補の組を除外した残りの点群を変位場算出点群として取得する。そして、データ取得部１０２は、取得した対応点候補の組のＩＤ番号を記憶部１１１に格納する。なお、対応点候補の組を取得する際には、記憶部１１１に格納されたＩＤ番号を検索し、まだ取得されていない点のみを対象とする。

（ステップＳ９０３：第４の変位場を算出）
ステップＳ９０３において、変位場算出部１０６は、ステップＳ９０２で生成された変位場算出点群を用いて、第４の変位場を算出する。この第４の変位場についても、図３のステップＳ３０１と同様の方法により算出することができる。

（ステップＳ９０４：整合度を算出）
次に、ステップＳ９０４において、整合度算出部１０８は、ステップＳ９０１で取得した第３の変位場と、ステップＳ９０３で算出した第４の変位場とを用いて、対応点候補の組の整合度を算出する。ここで、整合度の評価方法は、図３のステップＳ３０５の処理と同様である。第４の変位場は、既入力対応点群から対応点候補の組を除外した点群を用いて算出されているのに対し、第３の変位場は、既入力対応点群の全ての点を用いて算出される。そのため、第４の変位場に対する第３の変位場の整合度は、対応点候補の組がそれ以外の既入力対応点群とどの程度整合しているかを表す。

（ステップＳ９０５：全点の整合度算出完了？）
次に、ステップＳ９０５において、データ取得部１０２は、記憶部１１１に格納された取得済み対応点のＩＤ番号を参照し、既入力対応点群の中で対応点候補の組として整合度が算出されていない点が存在するか否かを判定する。この判定の結果、全ての点について整合度の算出が終了している場合は、ステップＳ９０６へ進み、まだ整合度が算出されていない点が存在する場合はステップＳ９０２へ戻る。

（ステップＳ９０６：判定結果を表示）
ステップＳ９０６においては、表示制御部１１０は判定結果を表示する。判定結果の表示は、例えば以下の方法によって行う。まず、整合度算出部１０８は、既入力対応点群における夫々の対応点の組について、図３のステップＳ３０６と同様の方法によって整合か不整合かを判定する。そして、表示制御部１１０は、不整合と判定された点のＩＤ番号を表示部１９０に表示させる。

以上のように本実施形態によれば、既に対応付けが行われた既入力対応点群における夫々の対応点の組について、他の対応点群に対して整合しているか否かを判定することができる。そのため、正しい対応点の組が後から追加されることによって不整合が生じるような、不適切な対応点の組を検出することが可能となる。

本実施形態で説明した既入力対応点群を対象とした整合性の判定処理は、第１及び第２の実施形態と組み合わせて用いることにより一層の効果を発揮する。例えば、既入力対応点群がある程度追加された時点でユーザが既入力対応点群の整合度を算出する処理を要求した場合に、それまでに取得した既入力対応点群を対象とした整合性の評価を実施するような形態が望ましい。あるいは、図３のステップＳ３１０や図５のステップＳ５０６における終了指示後の確認工程として、既入力対応点群を対象とした整合性の評価を実施するような形態が望ましい。

（その他の実施形態）
また、本発明は、以下の処理を実行することによっても実現される。即ち、上述した実施形態の機能を実現するソフトウェア（プログラム）を、ネットワーク又は各種記憶媒体を介してシステム或いは装置に供給し、そのシステム或いは装置のコンピュータ（またはＣＰＵやＭＰＵ等）がプログラムを読み出して実行する処理である。

１０２ データ取得部
１０４ 対応点候補取得部
１０８ 整合度算出部
１１０ 表示制御部

Claims (10)

1. 第１の画像と第２の画像との間の対応情報を取得する対応情報取得手段と、
   前記第１の画像上における第１の点の座標と前記第２の画像上における第２の点の座標とからなる対応点候補の組を取得する対応点候補取得手段と、
   前記対応情報取得手段によって取得された対応情報に対する、前記対応点候補取得手段によって取得された対応点候補の組の整合度を算出する整合度算出手段と、
   前記整合度算出手段の算出結果を表示部に表示させる表示制御手段と、
   を有することを特徴とする画像処理装置。
2. 前記対応情報を用いて前記第１の画像と前記第２の画像との間の第１の変位場を算出するとともに、前記対応情報と前記対応点候補の組とを用いて前記第１の画像と前記第２の画像との間の第２の変位場を算出する変位場算出手段をさらに備え、
   前記整合度算出手段は、前記第１の変位場と前記第２の変位場との差異に基づいて前記整合度を算出することを特徴とする請求項１に記載の画像処理装置。
3. 前記対応情報を用いて前記第１の画像と前記第２の画像との間の第１の変位場を算出する変位場算出手段をさらに備え、
   前記整合度算出手段は、前記変位場算出手段によって算出された第１の変位場と前記対応点候補の組とに基づいて整合度を算出することを特徴とする請求項１に記載の画像処理装置。
4. 前記対応情報と前記対応点候補の組とを用いて前記第１の画像と前記第２の画像との間の第２の変位場を算出する変位場算出手段をさらに備え、
   前記整合度算出手段は、前記変位場算出手段によって算出された第２の変位場に基づいて整合度を算出することを特徴とする請求項１に記載の画像処理装置。
5. 前記対応点候補取得手段は、前記第２の画像上における複数の点の夫々を前記第２の点と仮定した複数の対応点候補の組を取得し、
   前記整合度算出手段は、前記対応点候補の組の夫々に関して整合度を算出し、
   前記表示制御手段は、前記複数の対応点候補の組に対応する複数の整合度の分布を前記表示部に表示させることを特徴とする請求項１〜４の何れか１項に記載の画像処理装置。
6. 前記表示制御手段は、ユーザの操作によって指定された対応点候補の組を前記対応点候補取得手段が取得した後に、前記複数の整合度の分布を前記表示部に表示させることを特徴とする請求項５に記載の画像処理装置。
7. 前記対応情報は、前記第１の画像と前記第２の画像との間で既に対応付けされている複数の対応点の組の情報であることを特徴とする請求項１〜６の何れか１項に記載の画像処理装置。
8. 前記整合度算出手段によって算出された整合度が所定の条件を満たしている場合に、前記対応点候補の組を前記対応点の組として追加する追加手段をさらに備え、
   前記整合度算出手段は、前記追加手段によって前記対応点候補の組が追加された場合に、対応点の組のそれぞれを対応点候補の組として、それぞれの対応点候補の組について、該対応点候補の組を除いた対応点の組に対する整合度を算出することを特徴とする請求項７に記載の画像処理装置。
9. 第１の画像と第２の画像との間の対応情報を取得する対応情報取得工程と、
   前記第１の画像上における第１の点の座標と前記第２の画像上における第２の点の座標とからなる対応点候補の組を取得する対応点候補取得工程と、
   前記対応情報取得工程において取得された対応情報に対する、前記対応点候補取得工程において取得された対応点候補の組の整合度を算出する算出工程と、
   前記算出工程の算出結果を表示部に表示させる表示制御工程と、
   を有することを特徴とする画像処理方法。
10. 第１の画像と第２の画像との間の対応情報を取得する対応情報取得工程と、
    前記第１の画像上における第１の点の座標と前記第２の画像上における第２の点の座標とからなる対応点候補の組を取得する対応点候補取得工程と、
    前記対応情報取得工程において取得された対応情報に対する、前記対応点候補取得工程において取得された対応点候補の組の整合度を算出する算出工程と、
    前記算出工程の算出結果を表示部に表示させる表示制御工程と、
    をコンピュータに実行させることを特徴とするプログラム。