JP2013146540A - 画像処理装置および画像処理方法、並びに、画像処理プログラム - Google Patents

Abstract

【課題】同じ被写体を異なる種類のモダリティにより撮影して取得した２つの画像を被写体の空間的位置が一致するように精度よく位置合わせすることが出来る画像処理装置を提供する。
【解決手段】同じ被写体を異なる種類のモダリティにより撮影して取得した第１画像Ｖ１および第２画像Ｖ２を取得し、第１画像Ｖ１を変形して、変形した第１画像Ｖ１と第２画像Ｖ２との類似度を、両画像の対応する画素値の分布の相関性を評価する評価関数により評価して第１画像Ｖ１の画像変形量を推定し、推定された画像変形量に基づいて、第１画像Ｖ１を変形させた画像を生成する。評価関数は、変形した第１画像Ｖ１の画素値と第２画像Ｖ２の対応する画素値の相関性の尺度を表す項目を有し、この項目が第１画像Ｖ１および第２画像Ｖ２のそれぞれの対応する画素値の組合せの確からしさを表す蓋然性情報に基づいて相関性を評価する。
【選択図】図１

Description

本発明は、同一の被写体を撮影して取得した２つの画像において、両画像中の被写体の空間的な位置が一致するように、一方の画像の画像空間を変形させた画像を生成する画像処理装置および画像処理方法、並びに、画像処理プログラムに関する。

同一の被写体を同一または異なる撮像装置（モダリティ）を用いて異なるタイミングで撮影した２つの３次元画像を用いた画像診断において、両画像を重ねたとき、両画像中の被写体の空間的な位置が一致する変換関数を推定し、推定した変換関数を用いて一方の画像を変形させて２種類の画像を位置合わせする非剛体レジストレーション技術が注目されている。この非剛体レジストレーション技術においては、画像空間を所定間隔で区切る制御点が設定され、この制御点の位置を変位させて変形した一方の画像と他方の画像のそれぞれの画素値の類似性を評価する評価関数が最大となる制御点の画像変形量を決定し、このときの制御点の画像変形量に基づいて変換関数が推定される。

非特許文献１は、種類の異なるモダリティであるＰＥＴ装置（Positron emission tomography）およびＣＴ装置(Computed tomography)により同一の被写体をそれぞれ撮影して得られた画像間に非剛体レジストレーション処理を適用するものであり、評価関数に、２種類の画像の画素値の類似性の尺度として相互相関情報量（Mutual information）を用いている。

David Mattes, David R. Haynor, Hubert Vesselle, Thomas K. Lewellen, William Eubank, "Nonrigid multimodality image registration", Proceedings of the SPIE, volume 4322, p.1609-1620, 2001

ここで、異なる種類のモダリティによって同一の被写体を撮影した２つの画像間においては、モダリティの種類やモダリティの撮影方法の種類により、一方の画像の画素値に対応する他方の画像の画素値がある程度決まった範囲に属するものとなる。しかしながら、非特許文献１に記載の技術では、２つの画像間の画素値の分布の相関性のみを２つの画像の類似性の尺度に用いているため、例えば変形した一方の画像の画素値に対して、実際には対応することがあり得ない他方の画像の画素値を対応付けたものに高い類似度が算出される場合が起こりうる。もし誤った類似性の評価が行われ、誤った評価に基づいて一方の画像を変形させる変換パラメータが決定されると、一方の画像を変形しても他方の画像に一致しない結果となってしまう可能性がある。このため、両画像の画素値の組合せが妥当なものであるか否か区別して、より精度良く類似度の評価を行いたいという要望がある。

本発明は、かかる問題点に鑑みて、異なる種類のモダリティによって同一の被写体を撮影した２つの画像間において、画素値の組合せの確からしさを類似度の評価に反映させて、より正確に類似度の評価を行うことにより、精度良く２つの画像の被写体が一致するように一方の画像を変形させる画像変形量を推定することができる画像処理装置、画像処理方法、及び画像処理プログラムを提供することを目的とする。

上記目的を達成するために、本発明による画像処理装置は、被写体を第１のモダリティで撮像して得られた第１画像と、被写体を第１画像とは異なる第２のモダリティで撮像して得られた第２画像とを取得する画像取得部と、第１画像を変形して、変形した第１画像と第２画像との類似度を、変形した第１画像の画素値と第２画像の対応する画素値の分布の相関性を評価する評価関数により評価して変形した第１画像と第２画像が類似する第１画像の画像変形量を推定する画像変形量推定部と、推定された画像変形量に基づいて、第１画像を変形させた画像を生成する画像生成部を備え、画像変形量推定部が、第１画像および第２画像のそれぞれの対応する画素値の組合せの確からしさを表す蓋然性情報を取得する蓋然性情報取得部を備えるものであり、評価関数が、変形した第１画像の画素値と第２画像の対応する画素値の相関性の尺度を表す項目を備え、この項目が取得した蓋然性情報に基づいて相関性を評価することを特徴とする。

また、本発明による画像処理方法は、上記装置に実行させる画像処理方法であって、被写体を第１のモダリティで撮像して得られた第１画像と、被写体を第１画像とは異なる第２のモダリティで撮像して得られた第２画像とを取得する工程と、第１画像を変形して、変形した第１画像と第２画像との類似度を、変形した第１画像の画素値と第２画像の対応する画素値の分布の相関性を評価する評価関数により評価して変形した第１画像と第２画像が類似する第１画像の画像変形量を推定する工程と、推定された画像変形量に基づいて、第１画像を変形させた画像を生成する工程を備え、画像変形量を推定する工程が、第１画像および第２画像のそれぞれの対応する画素値の組合せの確からしさを表す蓋然性情報を取得する蓋然性情報取得工程を備える画像処理方法であり、評価関数が、変形した第１画像の画素値と第２画像の対応する画素値の相関性の尺度を表す項目を備え、この項目が取得した蓋然性情報に基づいて相関性を評価することを特徴とする。

また、本発明の画像処理プログラムは、コンピュータに上記方法を実行させることを特徴とする。

本発明における第１画像および第２画像は、異なるタイミングで異なるモダリティにより撮影された同一の被写体を表す画像であればよく、例えば、本発明の第１または第２モダリティとして適用可能なモダリティとしてＣＴ、ＭＲＩ、ＰＥＴ、ＳＰＥＣＴ、超音波画像などがあげられる。

また、本発明において、「評価関数が、変形した第１画像の画素値と第２画像の対応する画素値の相関性の尺度を表す項目を備え、この項目が取得した蓋然性情報に基づいて相関性を評価する」とは、蓋然性情報に基づいて、変形した第１画像および第２画像のそれぞれの対応する画素値の組合せの蓋然性が低いと判断できる場合に類似度を低く評価する種々の評価方法を適用できる。

例えば、上記本発明に係る画像処理装置の第１の態様において、評価関数における相関性の尺度を表す項目が、蓋然性情報を用いて重み付けされていることが好ましい。

また、本発明に係る画像処理装置の第１の態様において、蓋然性情報は、第１画像の画素値と第２画像の画素値の組合せの確からしさを表すものであれば、いかなる情報でもよい。例えば、蓋然性情報は、第２画像の各画素値が得られるという事象のもとで第１画像の各画素値が得られるという事象の発生する条件付き確率を表すものとすることができる。

また、本発明に係る画像処理装置の第１の態様において、評価関数は、変形した第１画像の画素値と第２画像の画素値の相関性の尺度を表す項目を備え、この相関性の尺度を表す項目を蓋然性情報に基づいて重み付けしたものに基づいて変形した第１画像と第２画像の類似度を算出するものであればいかなるものでもよい。例えば、評価関数を、両画像が類似しているほど評価値が大きくなるように定義してもよく、両画像が類似しているほど評価値が小さくなるように定義してもよい。なお、両画像が類似しているほど評価値が大きくなるように評価関数を定義する場合には、蓋然性情報は、両画像の画素値の組合せが確からしいほど、評価値が大きくなるように重み付けをするものであり、両画像が類似するほど評価値が小さくなるように評価関数を定義する場合には、蓋然性情報は、両画像の画素値の組合せが確からしいほど、評価値が小さくなるように重み付けをするものであるとする。

また、本発明に係る画像処理装置の第１の態様において、評価関数における相関性の尺度を表す項目が、変形した第１画像の画素値および第２画像の画素値を離散確率変数とした相互相関情報量または二乗損失相互情報量を表すものであることが好ましい。

また、本発明に係る画像処理装置の第１の態様において、蓋然性情報が、第２のモダリティでさらなる被写体を撮影して得られた第２参照画像の各画素値が得られるという事象のもとで第１のモダリティでさらなる被写体を撮影して得られた第１参照画像の各画素値が得られるという事象の発生する条件付き確率を表すものであることが好ましい。また、この場合に、さらなる被写体が被写体と同じ種類の他の被写体であることが好ましい。

また、本発明に係る画像処理装置の第１の態様において、蓋然性情報が、条件付き確率による重み付けを調整するために、第１参照画像における画素値分布を一様分布と近似した確率密度関数の逆数によりさらに重み付けされたものであることが好ましい。

また、本発明に係る画像処理装置の第１の態様において、蓋然性情報が、被写体の種類ごとに、所定の種類の被写体を第１のモダリティで撮影した場合に得られる第１の画素値範囲と、所定の種類の被写体を第２のモダリティで撮影した場合に得られる第２の画素値範囲をそれぞれ対応付けたものであり、評価関数が、蓋然性情報に基づいて、変形した第１画像と第２画像の各画素値が、第１の画素値範囲と第１の画素値範囲に対応付けられた第２の画素値範囲の両方を満たすものでない場合には、変形した第１画像と第２画像の類似度を低くするように重み付けされたものであってもよい。

また、上記場合に、第１または第２の画素値範囲が、それぞれ第１または第２のモダリティの撮像原理に基づいて、被写体の種類ごとに、第１および第２の画素値範囲を推定して算出されたものであることが好ましい。

なお、上記「被写体の種類」とは、診断画像において識別可能な被写体を表すものであれば任意に設定してよい。例えば、画素値としてそれぞれ所定の範囲を表す組成物、解剖学的構造物、解剖学的構造物を構成する各要素、解剖学的構造物の組織構造などに基づいて種類分けすることができ、一例として、水、空気、各解剖学的構造物をそれぞれ異なる種類として定義することができる。

また、本発明に係る画像処理装置の第２の態様において、蓋然性情報が、第２のモダリティでさらなる被写体を撮影して得られた第２参照画像の各画素値が得られるという事象のもとで第１のモダリティでさらなる被写体を撮影して得られた第１参照画像の各画素値が得られるという事象の発生する条件付き確率を表しており、評価関数における相関性の尺度を表す項目が、変形した第１画像の画素値と第２画像の画素値とを離散確率変数とした同時確率分布と、蓋然性情報との差異を表していることが好ましい。

また、本発明に係る画像処理装置の第２の態様において、評価関数として、下記式（９）により定義される評価関数を用いることが好ましい。
ただし、
ｆ：第２画像の画素値
ｍ：第１画像の画素値
Ｆ：第２画像の全画素値からなる集合
Ｍ：第１画像の全画素値からなる集合
ｍ；μ：第１画像の変形量をμとしたときの変形した第１画像の画素値。
ｐ（ｆ,ｍ；μ）：変形した第１画像の画素値と第２画像の画素値とを離散確率変数とした同時確率分布
ｐ（ｆ）: 第２画像の画素値を離散確率変数とした周辺確率分布
ｐ_Ｌ（ｍ；μ｜ｆ）：蓋然性情報
とする。

また、本発明に係る第２の態様の評価関数が、第１画像を変形して変形した第１画像と第２画像との類似度を、変形した第１画像の画素値と第２画像の対応する画素値の分布の相関性を評価するさらなる評価関数を備え、さらなる評価関数は、変形した第１画像の画素値と第２画像の対応する画素値の相関性の尺度を表す項目を備え、この項目が変形した第１画像の画素値および第２画像の画素値を離散確率変数とした相互相関情報量または二乗損失相互情報量を表していることが好ましい。

また、本発明に係る画像処理装置の第２の態様において、さらなる評価関数は、変形した第１画像の画素値と第２画像の対応する画素値の相関性の尺度を表す項目を取得した蓋然性情報に基づいて重み付けされていてもよい。

なお、本明細書では、変形した第１画像と第２画像が類似している場合を類似度が高い、両画像が類似していない場合を類似度が低いと称する。上記「類似度を低くするように重み付け」するとは、両画像が類似するほど評価値が大きくなるように評価関数を定義する場合には、評価値が小さくなるように重み付けをすることを意味し、両画像が類似するほど評価値が小さくなるように評価関数を定義する場合には、評価値が大きくなるように重み付けをすることを意味する。

本発明によれば、被写体を第１のモダリティで撮像して得られた第１画像と、被写体を第１画像とは異なる第２のモダリティで撮像して得られた第２画像とを取得する画像取得部と、第１画像を変形して、変形した第１画像と第２画像との類似度を、変形した第１画像の画素値と第２画像の対応する画素値の分布の相関性を評価する評価関数により評価して変形した第１画像と第２画像が類似する第１画像の画像変形量を推定する画像変形量推定部と、推定された画像変形量に基づいて、第１画像を変形させた画像を生成する画像生成部を備え、画像変形量推定部が、第１画像および第２画像のそれぞれの対応する画素値の組合せの確からしさを表す蓋然性情報を取得する蓋然性情報取得部を備え、評価関数が、変形した第１画像の画素値と第２画像の対応する画素値の相関性の尺度を表す項目を備え、この項目が取得した蓋然性情報に基づいて相関性を評価するようにしたので、実際にはあり得ないような画素値の組合せに対して、蓋然性情報に応じて類似度を低くするようにすることができる。このことにより、従来よりも正確に類似度を評価することができ、結果として好適に第１画像を第２画像に一致させた画像を生成することができる。

本発明の第１の実施形態における画像処理装置の電気的な概略ブロック図である。 本発明の第１の実施形態における画像処理装置の動作を示すフローチャートである。 位置合わせ前の第１画像（ＭＲ画像）と第２画像（ＣＴ画像）の一例を示す図である。 位置合わせ後の変形した第１画像（ＭＲ画像）と第２画像（ＣＴ画像）の一例を示す図である。 本発明の第４の実施形態において、変形した第１画像と第２画像を基準点からの距離に基づいて分割した図を表す図である。 本発明の第４の実施形態において、変形した第１画像と第２画像をｘｙ平面においてｘ軸からの角度に基づいて分割した図を表す図である。 本発明の第４の実施形態における画像処理装置の動作を示すフローチャートである。 本発明の第２の実施形態における画像処理装置の動作を示すフローチャートである。 本発明の第３の実施形態における画像処理装置の動作を示すフローチャートである。

以下、本発明の画像処理装置、画像処理プログラムおよび画像処理方法の実施形態について、図面を参照しながら詳細に説明する。本発明は、異なるモダリティにより同一の被写体を異なるタイミングで撮影して取得した２つの画像を位置合わせする処理を行う様々な分野で応用可能であるが、ここでは、医療分野における画像診断に本発明を適用した例をもとに説明を行う。

図１に、医師が使用するワークステーションに、画像処理プログラムをインストールすることにより実現された画像処理装置の概略構成を示す。画像処理装置１は、標準的なワークステーションの構成として、プロセッサおよびメモリ（いずれも図示せず）を備え、さらに、ＨＤＤ（Hard Disk Drive）等のストレージ２を備えている。また、画像処理装置１には、ディスプレイ３と、マウス、キーボード等の入力装置４が接続されている。

画像処理プログラムと画像処理プログラムが参照するデータは、インストール時にストレージ２に記憶され、起動時にメモリにロードされる。画像処理プログラムは、ＣＰＵに実行させる処理として、画像取得処理と、画像変形量推定処理と、画像生成処理と、表示制御処理とを規定している。

そして、プログラムの規定にしたがって、ＣＰＵが上記各処理を実行することにより、汎用のワークステーションは、後述の画像取得部１１と、画像変形量推定部１２と、画像生成部１４と、表示制御部１５として機能する。

ストレージ２には、撮影を担当する検査部門から転送された、第１画像および第２画像、もしくはデータベース検索により取得された第１画像および第２画像が記憶される。本実施の形態では、ある患者の検査において、同じ日の異なる時間に異なるモダリティにより一人の患者の頭部を撮影して得られた、第１画像Ｖ１（ＭＲ画像）および第２画像Ｖ２（ＣＴ画像）がそれぞれ、検査部門から転送され、両画像Ｖ１、Ｖ２のヘッダ情報などの任意の情報から周知の方法で取得した画素サイズ（Pixel Spacing）とスライス間隔(Slice Spacing)に基づいて両画像Ｖ１、Ｖ２のいずれか一方に両画像Ｖ１、Ｖ２のスケールを統一する処理が施され、ストレージ２に記憶されている。ここでは、図３Ａに第１画像Ｖ１および第２画像Ｖ２のイメージ画像を示す。第１画像Ｖ１および第２画像Ｖ２は、第１画像内の頭蓋骨などの解剖学的構造物が、やや第２画像の対応する解剖学的構造物より大きく、また、左右の大脳半球の境界部分の位置など解剖学的構造物の特徴的な位置もずれている。

画像取得部１１は、第１画像Ｖ１および第２画像Ｖ２をストレージ２からメモリに取得する。本実施形態では、画像処理装置１は、選択メニューにおいて所定の位置合わせ機能が選択されたことを検出すると、ユーザに、第１および第２画像の特定に必要な情報の選択または入力を促す。そして、ユーザの入力装置４の操作により、第１および第２画像が特定されると、画像取得部１１は、第１画像Ｖ１および第２画像Ｖ２をストレージ２からメモリに取得する。

画像変形量推定部１２は、第１画像Ｖ１を変形させて、変形した第１画像Ｖ１ａと第２画像Ｖ２との類似度を、変形した第１画像Ｖ１ａの画素値と第２画像Ｖ２の対応する画素値の分布の相関性を評価する評価関数Ｓ（μ）により評価して変形した第１画像と第２画像の被写体が一致する第１画像の画像変形量を推定する。

さらに詳細には、画像変形量推定部１２は、変形した第１画像Ｖ１ａおよび第２画像Ｖ２において、画像空間を所定間隔で区切る制御点ｘ１、ｘ２、…、ｘｎからなる集合Ｘをそれぞれ設定する。以下、制御点ｘ１、ｘ２、…、ｘｎの集合を制御点Ｘと記載する。なお、画像変形量推定部１２は、第１画像Ｖ１の制御点Ｘを変換関数ｇにより画像変形量μだけ変位させることにより、第１画像Ｖ１を変形する。以下、第１画像Ｖ１の制御点Ｘを変換関数ｇにより画像変形量μだけ変位させた制御点をｇ(Ｘ,μ)と記載し、第１画像Ｖ１の制御点Ｘを変換関数ｇにより画像変形量μだけ変位させることにより第１画像Ｖ１を変形した画像を、変形した第１画像Ｖ１ａと記載する。なお、本実施形態における画像変形量推定部１２は、非特許文献１に開示された方法により制御点Ｘや変画像形量μを設定し、非特許文献１に記載された変換関数ｇを変換関数ｇとして用いる。

次いで、画像変形量推定部１２は、変形した第１画像Ｖ１ａの制御点ｇ(Ｘ,μ)における画素値Ｍ(ｇ(Ｘ,μ))を取得するとともに第２画像Ｖ２の制御点Ｘの画素値Ｆ(Ｘ)を取得する。そして、変形した第１画像Ｖ１ａの各制御点ｇ(Ｘ,μ)における画素値Ｍ(ｇ(Ｘ,μ))と第２画像Ｖ２の各制御点Ｘの画素値Ｆ(Ｘ)の類似性を評価する評価関数Ｓ（μ）（レジストレーション関数）が最大となる制御点Ｘの画像変形量μを決定し、このときの制御点Ｘの画像変形量μに基づいて第１画像Ｖ１に対する変換関数を推定する。

本第１の実施形態では、画像変形量推定部１２は、後述の式（３）で表す評価関数Ｓ（μ）を用いて、変形した第１画像Ｖ１ａと第２画像Ｖ２の類似度を評価する。本実施形態における評価関数Ｓ（μ）によれば、変形した第１画像Ｖ１ａと第２画像Ｖ２の画素値分布が類似するほど評価関数Ｓ（μ）の値が大きくなるため、画像変形量推定部１２は、画像変形量μを変化させながら、評価関数Ｓ（μ）の変化量（または、μに対する偏導関数▽Ｓ（μ）の値の絶対値）が所定のしきい値以下となる場合の画像変形量μを、両画像の類似度が最大値をとる画像変形量（両画像が最も類似する画像変形量）と判断する。そして、この画像変形量μに基づいて、第１画像Ｖ１を変形する変換関数を決定する。なお、所定のしきい値は、式（３）に示す評価関数Ｓ（μ）の変化量が十分小さい値であると見なせる任意の値を設定してよい。また、この画像変形量μに基づいて、第１画像Ｖ１を変形する変換関数を決定する方法は周知の種々の方法を適用でき、ここでは非特許文献１に記載された方法を適用するものとする。

なお、ここでは評価関数を、類似度（評価値）が大きいほど類似度が大きくなるように定義しているが、評価値が小さいほど類似度が大きくなるように定義してもよい。また、評価関数の評価値（類似度）が最大（または最小）となる画像変形量を特定可能な方法であれば、非剛体レジストレーション技術の評価関数の最大値（または最小値）を算出するいかなる方法を用いて類似度が最大となる画像変形量μを特定してもよい。また、複数の異なる画像変形量μに対して評価関数Ｓ（μ）による評価値をそれぞれ算出して、算出した評価値のうち評価値が最大（または最小）であるものを特定し、特定された最大評価値（または最小評価値）に対応する画像変形量μを特定してもよい。

本第１の実施形態の画像変形量推定部１２は第１画像および第２画像のそれぞれの対応する画素値の組合せの確からしさを蓋然性情報として取得する蓋然性情報取得部１３を備え、評価関数における類似度の尺度を表す項目に対しこの蓋然性情報による重み付けを行う。

ここで、従来の評価関数について説明し、その後本第１の実施形態の評価関数および蓋然性情報について詳細に説明する。

非特許文献１に示されるように、異なるモダリティにより同一の被写体をそれぞれ撮影して取得した第１画像と第２画像の位置合わせを行う非剛体レジストレーションの技術において、相互相関情報量（mutual information、相互情報量）に基づく評価関数を用いて、第１画像の画素値と第２画像の対応する画素値の分布の類似度を評価することができる。

相互相関情報量とは、２つの確率変数ｆ、ｍの同時確率分布関数ｐ（ｆ，ｍ）と２つの確率変数ｆ、ｍの周辺確率分布関数ｐ（ｆ）、ｐ（ｍ）に基づいて２つの確率変数ｆ、ｍの相関性の尺度を表すもので、代表的には下記式（１）により定義されるものである。式（１）において、ｆは集合Ｆに属する離散確率変数であり、ｍは集合Ｍに属する離散確率変数である。また、相互相関情報量は、２つの確率変数ｆ、ｍが相関しているほど情報量が大きくなるものであり、言い換えると、一方の変数が与えられたときに他方の変数を高い確率で推測できるほど情報量が大きくなるものである。なお、２つの確率変数が完全に独立であれば、相互相関情報量は０となる。

異なるモダリティにより同一の被写体を撮影した画像間においては、各画像の撮影の物理的原理が異なるため、同じ種類の被写体であっても画像ごとに異なる画素値（信号値）を示す場合があり、単に画素値を比較しても両画像の類似性を判断することができない。例えば、ＣＴ画像においては、放射線の吸収率（透過度）が高いほど画素値（ＣＴ値）が大きくなり、空気、水、筋肉や肝臓や心臓などの臓器組織、骨の順に画素値が大きくなる。一方、ＭＲ画像においては、撮影対象物に含まれる水素原子の核磁気共鳴に応じて画素値が決定され、Ｔ１、Ｔ２強調画像などの撮影方法に応じて物質の信号値の大きさが変化する。例えば、ＭＲ画像におけるＴ１強調画像では、脂肪、筋肉、水は画素値の大きさがこの順に小さくなる。しかしながら、ＣＴ画像では、脂肪、筋肉、水は画素値の大きさの順はＭＲ画像のＴ１強調画像とは異なっているので、両者の画素値を単純に比較しても両画像の類似を判断することができない。

しかし、このような場合であっても、両画像の同一の解剖学的構造物を表している部分は、同一の解剖学的構造物に基づく共通の特徴に従って画素値が分布しているため、両者の画素値分布は互いに相関するものとなる。このことを利用して、非特許文献１では、異なるモダリティにより同一の被写体を撮影した２つの画像を、両者の画素値をそれぞれ確率変数とする相互相関情報量に基づいて、両画像の画素値の分布の相関が大きいほど、両画像が類似していると判断している。詳細には、式（２）で示すように、評価関数として、第１のモダリティにより撮影して得られた画像を変形した画像の画素値ｍ;μと、第２のモダリティにより撮影して得られた画像の各画素値ｆをそれぞれ確率変数とする相互相関情報量を用いている。（以下、第１画像を画像変形量μだけ変形した画像の各画素値を、第１画像からの画像変形量μを用いてｍ;μと表す。）式（２）において、相互相関情報量は、第２画像と変形した第１画像のそれぞれの画素値ｆ、ｍ;μの分布の相関を表しているため、第２画像と変形した第１画像の類似性の尺度として機能している。ここで、集合Ｆは変形第１画像の全画素値からなる集合であり、集合Ｍは第２画像の全画素値からなる集合である。なお、以下、式（２）のような相互相関情報量に基づく評価関数を、第１の評価関数と場合により称する。

しかしながら、上記式（２）のような従来の評価関数では、各モダリティの種類などにより得られる画素値の組合せの確からしさを全く評価していないため、第１画像を変形した画像の画素値に対応する第２画像の画素値がある程度決まった範囲となることが予想される場合などであっても、類似度算出の評価に反映させることができないという問題がある（第１の問題）。

本発明者は上記第１の問題に鑑み、評価関数に、第１画像を変形した画像と第２画像の対応する画素値の組合せの確からしさを表す蓋然性情報を表す項目を組み込むことが有効であることを見出した。そして、その第１の態様として、式（２）の評価関数の変形第１画像と第２画像の類似性の尺度を表す項目に式（３）のように、蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）を適用した。つまり、第１の評価関数における、第１画像を変形した画像と第２画像の類似性の尺度を表す項目を、第１画像を変形した画像と第２画像の画素値の組合せごとに、両画像の画素値の組合せの確からしさに応じて、両画像の組合せの確からしさを表す蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）で重み付けするものとした。

蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）は、所定の被写体を第１モダリティにより撮影して得られた画像と同じ種類の所定の被写体を第２モダリティにより撮影された画像の対応する画素値の組合せの確からしさを規定するものであればいかなる情報でもよい。ここでいう所定の被写体は、第１画像Ｖ１および第２画像Ｖ２に表す被写体と同じであってもよく、異なっていてもよい。

本実施形態では、第１モダリティにより撮影して得られた画像と第２モダリティにより撮影して得られた画像の対応する画素値の組合せの確からしさを、第１画像の各画素の画素値がｆとなるという事象をもとに第２画像の対応する各画素の画素値がｍ；μとなる事象が生じる条件付き確率として規定する。なお、下記式（４）に表されるように、条件付き確率分布関数は、同時確率分布関数ｐ_Ｌ（ｆ，ｍ；μ）と確率変数ｆの周辺確率分布関数ｐ_Ｌ（ｆ）により算出することができる。

本第１の実施形態においては、蓋然性情報取得部１３は、第１画像を画像変形量μだけ変形した画像の全ての制御点の画素値ｆ（ｆ∈Ｆ）と第２画像の対応する各制御点の画素値ｍ（ｍ∈Ｍ）を取得する。そして、周知の方法により、各制御点の画素値ｆ（ｆ∈Ｆ）ごとに対応する各制御点の各画素値ｍ（ｍ∈Ｍ）の分布をそれぞれ調べることにより同時確率分布関数ｐ_Ｌ（ｆ，ｍ；μ）を取得し、各制御点の画素値ｆについて分布を調べることにより周辺確率分布関数ｐ_Ｌ（ｆ）を取得する。そして、式（４）に基づいて、蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）を取得する。

なお、上記の条件付き確率の算出には、所定の被写体を第１モダリティにより撮影して得られた第１参照画像および同じ所定の被写体を第２モダリティにより撮影して得られた第２参照画像を用いることができる。なお、所定の被写体は必ずしも第１画像および第２画像に表される被写体と同じ種類の被写体である必要はないが、蓋然性情報をより精度よく算出するために、第１参照画像及び第２参照画像は、第１画像および第２画像に表す被写体と同じ種類の被写体を撮影したものであることが好ましい。また、第１参照画像と第２参照画像のペアは、１つであってもよく、複数であってもよい。第１参照画像と第２参照画像のペアを複数用いて、条件付き確率を算出した場合には、画素値の組合せの確からしさをより精度よく推定できると考えられるため好ましい。

また、蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）の算出または取得は、評価関数Ｓ（μ）を最大とする画像変形量を算出する処理以前に行うものであれば任意のタイミングで行ってよい。

画像生成部１４は、画像変形量推定部１２により決定された変換関数により第１画像Ｖ１を変換した画像Ｖ１Ａを生成する。

表示制御部１５は、画像生成部１４により生成された画像Ｖ１Ａおよび第２画像Ｖ２を比較可能にディスプレイ３に表示させる。また、ユーザの入力などの必要に応じて、取得した第１画像Ｖ１および第２画像Ｖ２および／または本実施形態の画像処理プログラムの実行過程で生成した各画像をディスプレイ３に表示させる。

図２は本発明の画像処理方法の好ましい実施形態を示すフローチャートである。図２を参照して、本実施形態の画像処理について説明する。

まず、画像取得部１１は、被写体を撮像して得られた第１画像（第１画像データ）Ｖ１および第２画像（第２画像データ）Ｖ２を取得する（Ｓ０１）。

次に、画像変形量推定部１２は、画像変形量μを設定して（Ｓ０２）、第１画像Ｖ１を変形する（Ｓ０３）。

次いで、画像変形量推定部１２は、変形した第１画像Ｖ１ａと第２画像Ｖ２各画素値の各画素値を取得して同時確率分布関数ｐ（ｆ,ｍ；μ）や周辺確率分布関数ｐ（ｆ）、ｐ（ｍ；μ）を算出し、蓋然性情報取得部１３は、算出された同時確率分布関数ｐ（ｆ,ｍ；μ）、周辺確率分布関数ｐ（ｆ）、ｐ（ｍ；μ）および蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）に基づいて、先述のように条件付き確率分布関数を算出することにより蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）を取得する（Ｓ０４）。そして、画像変形量推定部１２は、評価値として式（３）であらわされる評価関数Ｓ（μ）の変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜を算出する（Ｓ０５）。そして、算出された評価関数Ｓの変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜が所定のしきい値より大きい場合には（Ｓ０６のＮ）、画像変形量μを所定量Δμ増加させたμ＋Δμを新しい画像変形量μとして設定し（Ｓ０９）、再度Ｓ０４からＳ０５の処理を繰り返す。なお、画像変形量推定部１２は、評価値として、評価関数Ｓの変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜に換えて、評価関数Ｓ（μ）の偏導関数の絶対値｜▽Ｓ（μ）｜を用いてもよい。

一方算出された評価値の変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜が所定のしきい値以下である場合には（Ｓ０６のＹ）、画像変形量推定部１２は、この場合の画像変形量μを評価関数Ｓ（μ）が最大となる場合の画像変形量μとして取得し、この画像変形量μに基づいて、第１画像を変形する変換関数を決定する。画像生成部１４は、決定された変換関数を用いて第１画像を変換して再構成し、画像Ｖ１Ａを生成する（Ｓ０７）。そして、表示制御部１５は、生成された画像Ｖ１Ａと第２画像Ｖ２を比較可能に並列表示する（Ｓ０８）。図３Ｂに、生成された画像Ｖ１Ａと第２画像Ｖ２の比較表示の一例を示す。図３Ｂでは、画像Ｖ１Ａは第２画像Ｖ２に位置合わせされたものとなっており、画像Ｖ１Ａと第２画像Ｖ２は頭蓋骨などの解剖学的構造物のサイズが一致し、左右の大脳半球の境界などの位置も一致している。

上記本実施形態によれば、評価関数において、変形第１画像と第２画像の各画素値の類似性を表す項目を、変形第１画像および第２画像の各制御点の画素値の組合せの確からしさを表す蓋然性情報により重み付けしたため、第１モダリティにより撮影された画像の各制御点の画素値と第２モダリティにより撮影された画像の画素値とが、実際にはあり得ないような画素値の組合せに対して、蓋然性条件に応じて類似度を低くするように重み付けすることができる。このことにより、従来よりも正確に類似度を評価することができ、結果として好適に第１画像を第２画像に一致させた画像を生成することができる。

また、評価関数の第１モダリティにより撮影された画像と第２モダリティにより撮影された画像の各画素値の類似性を表す項目を、第１モダリティにより撮影された画像の各制御点の画素値がｆとなるという事象をもとに第２モダリティにより撮影された画像の対応する制御点の画素値がｍ；μとなるという事象が生じる条件付き確率により重み付けしたものとしたため、両画像の画素値の組合せの確からしさに応じてより精度良く両画像の類似度を算出することができる。

また、第１画像及び第２画像の各画素値に基づいて条件付き確率を算出して蓋然性情報として用いる場合には、両画像の画素値分布の情報さえあれば蓋然性情報を設定することができるため、蓋然性情報の設定のための情報収集や解析が不要であり、評価関数への適用が容易である。また、第１参照画像及び第２参照画像の各画素値に基づいて条件付き確率を算出して蓋然性情報として用いる場合には、両画像の画素値分布の情報さえあれば蓋然性情報を設定することができるため、評価関数への適用が容易である。また、同時確率分布関数ｐ（ｆ,ｍ；μ）や周辺確率分布関数ｐ（ｆ）に基づいて、相互相関情報量および条件付き確率の両方を算出できるため、計算効率がよい。

また、第１の実施形態の変形例として、蓋然性情報が、被写体の種類ごとに、所定の種類の被写体を第１のモダリティで撮影した場合に得られる第１の画素値範囲と、所定の種類の被写体を第２のモダリティで撮影した場合に得られる第２の画素値範囲をそれぞれ対応付けたものであり、評価関数が、蓋然性情報に基づいて、変形した第１画像と第２画像の各画素値が、第１の画素値範囲と第１の画素値範囲に対応付けられた第２の画素値範囲の両方を満たすものでない場合には、変形した第１画像と第２画像の類似度を低くするように重み付けされたものであってもよい。なお、本明細書において、類似度を低くするように重み付けするとは、類似度が大きいほど第１変形画像と第２画像が類似していると評価する評価関数を用いる場合には、類似度を小さくするように重み付けすることを意味し、類似度が小さいほど第１変形画像と第２画像が類似していると評価する評価関数を用いる場合には、類似度を大きくするように重み付けすることを意味する。

一例として、上記第１の実施形態の変形例の場合の評価関数を、以下の式（５）のように定義できる。以下の例では、第２モダリティにより撮影された画像の画素値ｆ（ｆ∈Ｆ）ごとに、第１モダリティにより撮影された画像の取り得る画素値の範囲ｍ_ｍｉｎ ^ｆ≦ｍ（μ）≦ｍ_ｍａｘ ^ｆを対応づけて記憶する。そして、第１画像の画素値ｆ（第１の範囲）に対して第２画像の画素値ｍ（μ）が対応づけられた範囲ｍ_ｍｉｎ ^ｆ≦ｍ（μ）≦ｍ_ｍａｘ ^ｆ（第２の範囲）に属する場合には、評価関数の第１画像と第２画像の相関性の尺度を示す項目に１により重み付けし、それ以外の場合には、実際にはあり得ない画素値の組合せと判断して、評価関数の第１画像と第２画像の相関性の尺度を示す項目に０により類似度が低くなるように重み付けする。

上記の第１の実施形態の変形例において、第１画像と第２画像の画素値の対応する画素値の範囲を、予め撮影した試験画像を解析した既知の情報により規定してもよい。例えば、第１のモダリティにより撮影された画像における、水、空気、各解剖学的構造物（または解剖学的構造物を構成する各要素）などの被写体の種類ごとに、各種類の被写体の示す画素値の範囲を、第２のモダリティにより撮影された画像における、水、空気、各解剖学的構造物などのそれぞれ対応する対象物の示す画素値の範囲に対応づけて蓋然性情報として記憶する。そして、画像変形量推定部１２が、蓋然性情報に対応づけられた画素値の範囲の組合せのいずれにも該当しない画素値の組合せに対して、類似度を低くすることが考えられる。この場合には、蓋然性情報において、同一の被写体を表していることが既知である２つの画像の画素値の情報に基づいて対応する画素値範囲を設定しているため、より正確に蓋然性情報を定義することができる。

また、第１画像と第２画像の画素値の対応する画素値の範囲を、各モダリティの撮像原理に基づいて推定される理論的な画素値から決定してもよい。例えば、被写体となりうる対象物の種類毎に、各モダリティによりその対象物を撮影した場合の画素値を推定して対応付けて記憶しておく。例えば、水、空気、各解剖学的構造物（または解剖学的構造物を構成する各要素）などの被写体の種類ごとに、各被写体を撮影したときに得られる第１のモダリティによる画像の画素値の範囲と、第２のモダリティによる画像における画素値の範囲を、撮影原理に基づいて推定する。そして、対象物ごとに、推定された第１のモダリティによる画像の画素値の範囲および第２のモダリティによる画像の画素値の範囲を対応付けて記憶する。そして、蓋然性情報に対応づけられた画素値の範囲の組合せのいずれにも該当しない画素値の組合せに対して、類似度の重み付けを低くするようにすることが考えられる。この場合には、蓋然性情報において、撮像原理から推定される２つの画像の画素値の情報に基づいて対応する画素値範囲を設定しているため、より正確に蓋然性情報を定義することができる。

また、上述の各蓋然性情報を、さらに重み付けを調整するようにスケーリング（正規化）したものとしてもよい。この場合には、条件付き確率が非常に小さい値になる場合など、第１および第２参照画像の相関性を表す項目に対する重み付けが大きくなりすぎない（または小さくなりすぎない）ように調整することができる。また、重み付けの調整をする周知の種々の方法を適用することができる。

例えば、条件付き確率を蓋然性情報に用いる場合、画素値の分布の形状による影響を受けにくくなるように、蓋然性情報を正規化することが好ましい。条件付き確率の算出に用いる要素である、画素値を確率変数とした確率密度関数は、画素値の分布が平坦である場合には値が小さくなり、画素値の分布が集中する部分を有する急峻なものである場合には値が大きくなるように、画素値の分布の形状に応じて値が変動する性質を有するからである。このため、まず、第２モダリティにより撮影された画像の画素値ｆ（ｆ∈Ｆ）ごとに、第１モダリティにより撮影された画像の取り得る画素値の範囲ｍ_ｍｉｎ ^ｆ≦ｍ（μ）≦ｍ_ｍａｘ ^ｆを対応づけて記憶する。そして、第２モダリティにより撮影された画像の画素値ｆ（ｆ∈Ｆ）ごとに、この画素値ｆに対応付けられた第１参照画像の画素値分布を一様分布と近似し、第１参照画像の確率密度関数Ｐ_Ｕ（ｍ）（＝ｍ_ｍａｘ ^ｆ−ｍ_ｍｉｎ ^ｆ）を算出する。そして、下記式（６）に示すように、蓋然性情報を確率密度関数Ｐ_Ｕ（ｍ）の逆数Ａ_ｍによりさらに重み付けしたものとする。

この場合には、第２モダリティにより撮影された画像の画素値ｆに対応づけられた第１モダリティにより撮影された画像の画素値ｍ_ｍａｘ ^ｆ＜ｍ（μ）＜ｍ_ｍｉｎ ^ｆの分布形状の影響を抑制することができ、第１モダリティにより撮影された画像の画素値と第２モダリティに寄り撮影された画像の画素値の組合せの確かしさをより好適に反映して、より適切に重み付けを行うことができる。

また、本第１の実施形態においては、異なるモダリティにより同一の被写体を撮影して得られた２つの画像の相関性の尺度を表す項目として、式（２）で表す相互相関情報量を用いるが、２つの画像の相関性の尺度を表す項目は、周知の種々の方法により定義された相互相関情報量を用いることができる。また、評価関数は、２つの画像の相関性の尺度を表す項目が蓋然性情報により重み付けされたものであれば、変形の滑らかさの限界条件を規定するなどの他の項目を更に加えたものであってもよく、周知の種々の変形を行ったものであってもよい。

例えば、第１の実施形態において第１の評価関数における、２つの画像の相関性の尺度を表す項目を、相互相関情報量に換えて二乗損失相互情報量を表すものとしてもよく、この場合にも、上記第１の実施形態と同様の効果が得られる。異なるモダリティにより同一の被写体をそれぞれ撮影して取得した第１画像と第２画像の位置合わせを行う非剛体レジストレーションの技術において、二乗損失相互情報量に基づく評価関数を用いて、第１画像の画素値と第２画像の対応する画素値の分布の類似度を評価する場合にも、相互相関情報量に基づく評価関数を用いた場合と同様に、各モダリティの種類などにより得られる画素値の組合せの確からしさを全く評価しておらず、第１の問題が共通して発生するためである。

二乗損失相互情報量は、２つの確率変数ｆ、ｍの同時確率分布関数ｐ（ｆ，ｍ；μ）と２つの確率変数ｆ、ｍ；μの周辺確率分布関数ｐ（ｆ）、ｐ（ｍ；μ）に基づいて２つの確率変数ｆ、ｍの相関性の尺度を表すものである。例えば、第１実施形態における第１画像と第２画像の相関性を評価する評価関数として、二乗損失相互情報量を下記式（８）のように適用することができる。式（８）において、ｆは集合Ｆに属する離散確率変数であり、ｍ；μは集合Ｍに属する離散確率変数である。また、二乗損失相互情報量は、２つの確率変数ｆ、ｍ；μが相関しているほど情報量が大きくなるものであり、言い換えると、一方の変数が与えられたときに他方の変数を高い確率で推測できるほど情報量が大きくなるものである。なお、２つの確率変数が完全に独立であれば、二乗損失相互情報量は０となる。

なお、式（８）においては、二乗損失相互情報量を第１のモダリティにより撮影して得られた画像を変形した画像の画素値ｍ;μと、第２のモダリティにより撮影して得られた画像の各画素値ｆをそれぞれ確率変数としており、式（８）において、二乗損失相互情報量は、第２画像と変形した第１画像のそれぞれの画素値ｆ、ｍ;μの分布の相関を表しているため、第２画像と変形した第１画像の類似性の尺度として機能している。ここで、集合Ｆは変形第１画像の全画素値からなる集合であり、集合Ｍは第２画像の全画素値からなる集合である。

また、本発明者は、第１の問題に対し、蓋然性情報を別の態様（第２の態様）で評価関数に組み込む方法も有効であることを見出した。以下、この方法を第２の実施形態として説明する。なお、第２の実施形態においては、蓋然性情報の取得、評価関数Ｓ（μ）を用いた第１画像Ｖ１を変形した画像Ｖ１ａと第２画像Ｖ２の類似度の評価以外の処理については、第１の実施形態と同じであり、各機能ブロックの機能も共通している。以下、第１の実施形態と異なる点を中心に説明し、第１の実施形態と同じ点については説明を省略する。

図６は、第２の実施形態における画像処理の流れを表すフローチャートである。以下、第２の実施例について図６に従って説明する。まず、画像取得部１１は、第１の実施形態と同様に、被写体を撮像して得られた第１画像（第１画像データ）Ｖ１および第２画像（第２画像データ）Ｖ２を取得する（Ｓ３１）。また、画像変形量推定部１２は、第１の実施形態と同様に、画像変形量μを設定して（Ｓ３２）、第１画像Ｖ１を変形する（Ｓ３３）。なお、図６におけるＳ３１−Ｓ３３の処理は、図２におけるＳ０１−Ｓ０３の処理にそれぞれ相当する。

続いて、蓋然性情報取得部１３は、変形した第１画像Ｖ１ａと第２画像Ｖ２各画素値の各画素値を取得して同時確率分布関数ｐ（ｆ,ｍ；μ）や周辺確率分布関数ｐ（ｆ）、ｐ（ｍ；μ）を算出する。また、蓋然性情報を第２のモダリティ（ＣＴ装置）でさらなる被写体を撮影して得られた第２参照画像（ＣＴ画像）の各画素値が得られるという事象のもとで第１のモダリティ（ＭＲ装置）でさらなる被写体を撮影して得られた第１参照画像（ＭＲ画像）の各画素値が得られるという事象の発生する条件付き確率を算出する（Ｓ３４）。

ここでは、蓋然性情報取得部１３は、Ｎ個の被写体について、Ｎ個の被写体に含まれる１つの被写体を第１モダリティにより撮影して得られた第１参照画像とその被写体を第２モダリティにより撮影して得られた第２参照画像のペア（参照画像のペア）を用意し、各参照画像のペアごとに、第１および第２参照画像に含まれる全ての画素について、第１参照画像の画素の画素値と、第１参照画像のその画素が表す被写体部分を表す第２参照画像の画素の画素値を対応付けて取得する。そして、周知の方法により、Ｎ個の参照画像のペアから取得した第１参照画像の画素値ｍ;μ_refと、その画素値ｍ;μ_refに対応付けられた第２参照画像の画素の画素値ｆ_refに基づいて、各第１参照画像の画素値ｍ;μ_ref（または第２参照画像の画素値）ごとに、第２参照画像の各画素値ｆ_refの分布を算出して、条件付き確率分布関数Ｐ_Ｌ（ｍ；μ_ref｜ｆ_ref）を取得する。なお、蓋然性情報取得部１３は、条件付き確率分布関数Ｐ_Ｌ（ｍ；μ_ref｜ｆ_ref）が算出可能なものであれば、周知の種々の方法を用いてよい。ここでは、Ｎ個の被写体をそれぞれ撮影した、Ｎ個の第１参照画像と第２参照画像のペアが用いられたものとする。より正確に蓋然性情報を算出するために、参照画像のペアの数Ｎが多いことが好ましく、Ｎ個の被写体が、第１画像Ｖ１および第２画像Ｖ２と同じ種類の被写体であることが好ましい。

次いで、画像変形量推定部１２は、算出された条件付き確率分布関数Ｐ_Ｌ（ｍ；μ_ref｜ｆ_ref）を蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）として後述の式（９）で表す評価関数Ｓ_Ｄ（μ）に適用し、評価関数Ｓ_Ｄ（μ）の変化量｜Ｓ_Ｄ（μ）−Ｓ_Ｄ（μ−Δμ）｜を評価値として算出する（Ｓ３５）。なお、Ｓ３５およびＳ３６の処理において、画像変形量推定部１２は、評価関数Ｓ_Ｄ（μ）の変化量｜Ｓ_Ｄ（μ）−Ｓ_Ｄ（μ−Δμ）｜に換えて、評価関数Ｓ_Ｄ（μ）の偏導関数▽Ｓ_Ｄ（μ）の絶対値｜▽Ｓ_Ｄ（μ）｜を用いてもよい。

第２の実施形態における画像変形量推定部１２は、評価関数Ｓ_Ｄにおける、変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の画素値ｆとの相関性の尺度を表す項目（相関性の尺度を表す項目）として、変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の画素値ｆとを離散確率変数とした同時確率分布関数ｐ（ｆ，ｍ；μ）と、蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）との差異（以下、差異Ａと場合により記載する。）を表す項目を用いる。以下、本明細書において、この評価関数（相関性の尺度を表す項目として、差異Ａを表す項目を用いた評価関数）を第２の評価関数と場合により称する。

本実施形態においては、上記差異Ａを、上記同時確率分布関数ｐ（ｆ，ｍ；μ）と、上記蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）とのカルバック・ライブラー距離（ＫＬ距離）として下記式（９）のように定義する。

式（９）において、蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）と周辺確率分布関数ｐ（ｆ）の乗算は、同時確率分布（以下、蓋然性同時確率分布と称す。）を表している。式（９）において、相関性の尺度を表す項目は、同時確率分布関数ｐ（ｆ，ｍ；μ）と蓋然性同時確率分布Ｐ_Ｌ（ｍ；μ｜ｆ）＊ｐ（ｆ）との関係に着目すると、変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の画素値ｆとを離散確率変数とした同時確率分布関数ｐ（ｆ，ｍ；μ）と、蓋然性同時確率分布Ｐ_Ｌ（ｍ；μ｜ｆ）＊ｐ（ｆ）との差異（差異Ａ）を表している。式（９）によれば、同時確率分布関数ｐ（ｆ，ｍ；μ）が蓋然性同時確率分布Ｐ_Ｌ（ｍ；μ｜ｆ）＊ｐ（ｆ）に近くなるほど小さい評価値が得られるため、評価関数Ｓ_Ｄ（μ）の評価値を最小化することにより、変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の画素値ｆの同時確率分布関数を、参照画像のペアの画素値の組合せに基づく蓋然性同時確率分布Ｐ_Ｌ（ｍ；μ｜ｆ）＊ｐ（ｆ）に近づけることができる。

このため、式（９）によれば、差異Ａに基づいて、同時確率分布関数ｐ（ｆ，ｍ；μ）と、第１参照画像と第２参照画像のペアから得られる蓋然性同時確率分布との差異を評価して、参照画像のペアの画素値の組合せに基づく蓋然性同時確率分布に近い（類似する）ほど類似度を高く評価する（低い評価値を算出する）ことにより、変形した第１画像Ｖ１ａの各制御点の画素値ｍ；μと第２画像Ｖ２の画素値ｆとが、実際にはあり得ないような画素値の組合せを取る場合に対して、蓋然性同時確率分布に基づいて類似度を低く（評価値を高く）評価することができる。このことにより、従来よりも正確に類似度を評価することができ、結果として好適に第１画像Ｖ１を第２画像Ｖ２に一致させた画像を生成することができる。

以上により、評価関数を式（９）のように定義することにより、先述の差異Ａに基づいて、変形した第１画像Ｖ１ａの各制御点の画素値ｍ；μと第２画像Ｖ２の画素値ｆとが、実際にはあり得ないような画素値の組合せを取る場合に類似度を低く評価しつつ、変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の画素値ｆの同時確率分布関数が第２画像Ｖ２の周辺確率分布関数に近づくほど類似度が高くなるように評価することができる。

続いて、第２の実施形態において、画像変形量推定部１２は、式（９）に示す評価関数Ｓ_Ｄ（μ）の変化量｜Ｓ_Ｄ（μ）−Ｓ_Ｄ（μ−Δμ）｜が所定のしきい値より大きい場合には（Ｓ３６のＮ）、画像変形量μを所定量Δμ増加させたμ＋Δμを新しい画像変形量μとして設定し（Ｓ３９）、再度Ｓ３４からＳ３５の処理を繰り返す。一方、算出された評価関数Ｓ_Ｄ（μ）の変化量｜Ｓ_Ｄ（μ）−Ｓ_Ｄ（μ−Δμ）｜が所定のしきい値以下である場合には（Ｓ３６のＹ）、画像変形量推定部１２は、この場合の画像変形量μを評価関数Ｓ_Ｄ（μ）が最小となる場合の画像変形量μとして取得し、この画像変形量μに基づいて、第１画像を変形する変換関数を決定する。

画像生成部１４は、決定された変換関数を用いて第１画像を変換して再構成し、画像Ｖ１Ａを生成する（Ｓ３７）。そして、表示制御部１５は、生成された画像Ｖ１Ａと第２画像Ｖ２を比較可能に並列表示して処理を終了する（Ｓ３８）。なお、図６におけるＳ３６−Ｓ３９の処理は、図２におけるＳ０６−Ｓ０９の処理にそれぞれ相当する。

なお、第２の実施形態において、変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の画素値ｆとを離散確率変数とした同時確率分布関数ｐ（ｆ，ｍ；μ）と、蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）との差異を表せるものであれば、評価関数Ｓ_Ｄ（μ）における変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の画素値ｆとの相関性の尺度を表す項目を種々の方法により定義してよい。

また、第３の実施形態として、第１の実施形態に示すような第１の評価関数と第２の実施形態に示すような第２の評価関数を組み合わせたものを評価関数として用いてもよい。

第３の実施形態においては、第１および第２の実施形態とは、蓋然性情報の取得、評価関数Ｓ_Ｄ（μ）−λＳ_Ｄ（μ）を用いた第１画像Ｖ１を変形した画像Ｖ１ａと第２画像Ｖ２の類似度の評価以外の処理については、第１の実施形態と同じであり、各機能ブロックの機能も共通している。以下、第１の実施形態と異なる点を中心に説明し、第１の実施形態と同じ点については説明を省略する。

第３の実施形態における評価関数は、第２の実施形態における第２の評価関数Ｓ_Ｄ（μ）と第１の評価関数Ｓ（μ）（さらなる評価関数）とを組み合わせたものである。

第３の実施形態における第１の評価関数は、第１の実施形態における第１の評価関数に相当するものであって、上記式（２）で定義される。第３の実施形態における第１の評価関数は、第１画像Ｖ１を変形して、変形した第１画像Ｖ１ａと第２画像Ｖ２との類似度を、変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の対応する画素値ｆの分布の相関性を評価するものであって、変形した第１画像の画素値と第２画像の対応する画素値の相関性の尺度を表す項目を備えたものであればよい。例えば、項目が変形した第１画像の画素値および第２画像の画素値を離散確率変数とした相互相関情報量または二乗損失相互情報量を表すものとすることができる。

また、第３の実施形態における第２の評価関数は、第２の実施形態における第２の評価関数に相当するものであって、上記式（９）で定義される。また、式（９）における蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）として、第２の実施形態と共通する蓋然性情報が用いられる。第３の実施形態における第２の評価関数は、第１画像Ｖ１を変形して、変形した第１画像Ｖ１ａと第２画像Ｖ２との類似度を、変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の対応する画素値ｆの分布の相関性を評価するものであって、変形した第１画像の画素値と第２画像の対応する画素値の相関性の尺度を表す項目を備え、この項目が変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の画素値ｆとを離散確率変数とした同時確率分布関数ｐ（ｆ，ｍ；μ）と、蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）との差異（差異Ａ）を表すものであればよい。

上記第３の実施形態における第１の評価関数と第２の評価関数は、変形した第１画像と第２画像が類似しているほど類似度が高くなるように、評価関数に組み込まれているものであればよく、例えば、第１の評価関数と第２の評価関数を重み付け加算したものを評価関数とすることができる。

ここでは、評価関数を、第１の評価関数Ｓ（μ）と第２の評価関数Ｓ_Ｄ（μ）を重み付け加算した下記式（１０）により定義する。なお、式（１０）においてλは、条件に応じて適宜設定される重み付け係数である。
Ｓ（μ）−λＳ_Ｄ（μ） （１０）

図７は、第２の実施形態における画像処理の流れを表すフローチャートである。以下、第２の実施例について図７に従って説明する。

まず、画像取得部１１は、第１の実施形態と同様に、被写体を撮像して得られた第１画像（第１画像データ）Ｖ１および第２画像（第２画像データ）Ｖ２を取得する（Ｓ４１）。次に、画像変形量推定部１２は、第１の実施形態と同様に、画像変形量μを設定して（Ｓ４２）、第１画像Ｖ１を変形する（Ｓ４３）。なお、図７におけるＳ４１−Ｓ４３の処理は、図２におけるＳ０１−Ｓ０３の処理にそれぞれ相当する。

次いで、画像変形量推定部１２は、変形した第１画像Ｖ１ａと第２画像Ｖ２各画素値の各画素値を取得して同時確率分布関数ｐ（ｆ,ｍ；μ）や周辺確率分布関数ｐ（ｆ）、ｐ（ｍ；μ）を算出し、蓋然性情報取得部１３は、第２の実施形態と同様に条件付き確率を算出することにより蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）を取得する（Ｓ４４）。

次に、画像変形量推定部１２は、第１の実施形態と同様に、式（２）であらわされる第１の評価関数Ｓ（μ）の変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜を評価値として算出する（Ｓ４５）。

次に、画像変形量推定部１２は、第２の実施形態と同様に、式（９）であらわされる第２の評価関数Ｓ_Ｄ（μ）の変化量｜Ｓ_Ｄ（μ）−Ｓ_Ｄ（μ−Δμ）｜を評価値として算出する（Ｓ４６）。

そして、画像変形量推定部１２は、式（１０）で示す評価式の変化量を評価値として算出する。つまり、式（１０）に基づいて、第１の評価関数の変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜と第２の評価関数の変化量｜Ｓ_Ｄ（μ）−Ｓ_Ｄ（μ−Δμ）｜を重み付け加算して式（１０）で示す評価式の変化量（｜｜Ｓ（μ）−Ｓ（μ−Δμ）｜−λ｜Ｓ_Ｄ（μ）−Ｓ_Ｄ（μ−Δμ）｜｜）を算出する（Ｓ４７）。

なお、Ｓ４７およびＳ４８の処理において、画像変形量推定部１２は、上記の評価値の変化量（｜｜Ｓ（μ）−Ｓ（μ−Δμ）｜−λ｜Ｓ_Ｄ（μ）−Ｓ_Ｄ（μ−Δμ）｜｜）に換えて、式（１０）の偏導関数（｜▽Ｓ（μ）｜−λ｜▽Ｓ_Ｄ（μ）｜｜）を用いてもよい。この場合、Ｓ４５の処理において、画像変形量推定部１２は、第１の評価関数Ｓ（μ）の変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜に換えて、第１の評価関数Ｓ（μ）の偏導関数▽Ｓ（μ）の絶対値｜▽Ｓ（μ）｜を算出し、Ｓ４６の処理において、画像変形量推定部１２は、第２の評価関数Ｓ_Ｄ（μ）の変化量｜Ｓ_Ｄ（μ）−Ｓ_Ｄ（μ−Δμ）｜に換えて、第２の評価関数Ｓ_Ｄ（μ）の偏導関数▽Ｓ_Ｄ（μ）の絶対値｜▽Ｓ_Ｄ（μ）｜を算出する。

そして、算出された評価値の絶対値が所定のしきい値より大きい場合には（Ｓ４８のＮ）、画像変形量μを所定量Δμ増加させたμ＋Δμを新しい画像変形量μとして設定し（Ｓ５１）、再度Ｓ４３からＳ４８の処理を繰り返す。

一方算出された評価関数Ｓ（μ）−λＳ_Ｄ（μ）の変化量（｜｜Ｓ（μ）−Ｓ（μ−Δμ）｜−λ｜Ｓ_Ｄ（μ）−Ｓ_Ｄ（μ−Δμ）｜｜）が所定のしきい値以下である場合には（Ｓ４８のＹ）、画像変形量推定部１２は、この場合の画像変形量μを式（１０）で表される評価関数が最大となる場合の画像変形量μとして取得し、この画像変形量μに基づいて、第１画像を変形する変換関数を決定する。画像生成部１４は、決定された変換関数を用いて第１画像を変換して再構成し、画像Ｖ１Ａを生成する（Ｓ４９）。そして、表示制御部１５は、生成された画像Ｖ１Ａと第２画像Ｖ２を比較可能に並列表示して処理を終了する（Ｓ５０）。なお、図６におけるＳ４８−Ｓ５１の処理は、図２におけるＳ０６−Ｓ０９の処理にそれぞれ相当する。

また、第３の実施形態の変形例として、第１の評価関数において、第１の実施形態のように、変形した第１画像の画素値と第２画像の対応する画素値の相関性の尺度を表す項目を取得した蓋然性情報に基づいて重み付けされたものとしてもよい。例えば、下記式（１１）に示すように、上記式（１０）のＳ（μ）に換えて、式（３）に示す評価関数を用いてよい。なお、式（１１）においては、式（３）に示す評価関数をＳ_Ｌ（μ）と表している。なお、式（１１）においてλは、条件に応じて適宜設定される重み付け係数である。
Ｓ_Ｌ（μ）−λＳ_Ｄ（μ） （１１）

第３の実施形態においては、第１の評価関数と第２の評価関数を組み合わせて用いることにより、第1の評価関数によって、変形した第１画像Ｖ１ａと第２画像Ｖ２との相関性があるほど類似度を高く評価しつつ、第２の評価関数によって、変形した第１画像Ｖ１ａの画素値ｍ；μと第２画像Ｖ２の画素値ｆとを離散確率変数とした同時確率分布関数ｐ（ｆ，ｍ；μ）と、第１参照画像と第２参照画像のペアから予想される蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）との差異が小さくなるほど類似度を高く評価することにより、異なる２つの観点から類似度を評価して、精度よく変形した第１画像と第２画像との相関性を評価することができる。

式（１１）のように、第１の評価関数を取得した蓋然性情報に基づいて重み付けされたものとした場合には、第１の実施形態における評価関数の利点をさらに備えることができる。つまり、第１の評価関数において、変形第１画像と第２画像の各画素値の類似性を表す項目を、変形第１画像および第２画像の各制御点の画素値の組合せの確からしさを表す蓋然性情報により重み付けしたため、第１モダリティにより撮影された画像の各制御点の画素値と第２モダリティにより撮影された画像の画素値とが、実際にはあり得ないような画素値の組合せに対して、蓋然性条件に応じて類似度を低くするようにより重み付けすることができる。このことにより、さらに好適に第１画像を第２画像に一致させた画像を生成することができる。

なお、第３の実施形態における第１の評価関数に換えて、変形した第１画像と第２画像との相関性を評価できる他の評価関数を用いてもよい。また、変形した第１画像と第２画像との相関性を評価することができるものであれば、第１の評価関数と第２の評価関数を別の方法で組み合わせてもよい。

次いで、本発明の第４の実施形態について以下に説明する。第４の実施形態では、画像変形量推定部１２が、所定の分割条件に従って変形した第１画像を複数の第１分割画像に分割するとともに第２画像を第１分割画像に対応する複数の第２分割画像に分割し、第１分割画像およびこれに対応する第２分割画像の一対ごとに、一対の分割画像間の画素値の分布の類似性を表すように定義された分割画像類似度に基づいて、第１画像を変形した画像と第２画像の類似性を評価する点が第１の実施形態と異なる。

第４の実施形態において、本発明者は、非特許文献１に記載された式（２）で表されるような従来の２つの画像の画素値の分布の相関性のみに基づいて２つの画像の類似性を判断している評価関数に関し、画像の空間的な特徴を判別できず、類似性の判定を誤る場合があるというさらなる問題（第２の問題）に着目した。

例えば、非特許文献１に記載された方法では、同じ画素値範囲に属する複数の被写体を撮影した２つの画像間において被写体の数や空間的位置が異なっていても、２つの画像間で同じ画素値範囲に属する画素の総数が同じ場合には、類似していると判断してしまう。具体的には、患者の胸部を撮影した一つの画像において、膵臓と肝臓を表す各画素の画素値は、所定の同じ画素値範囲に属するものとなる。もし、同一の患者の胸部を撮影した他方の画像に、この所定の範囲の画素値範囲を示す１つの陰影があり、この陰影の体積が、先述の一つの画像の肝臓と膵臓の体積の総和と等しい場合、非特許文献１の方法では両画像が類似していると判断してしまう可能性がある。

本発明者は上記第２の問題に鑑み、第１変形画像と第２画像を互いに対応する領域ごとに所定の分割条件に従って分割し、この分割した領域ごとに類似度を算出し、分割した領域ごとの類似度に基づいて、第１画像と第２画像の類似度を評価することが有効であることを見出した。この場合には、離れて位置する同じ画素値範囲に属する複数の被写体が別々の分割された領域に含まれる可能性が高いため、互いに異なる分割領域に位置する複数の被写体を表す画素値間で相関性が評価されることが抑制され、上記第２の問題の発生を低減することができるからである。

なお、分割条件は、所定の決まりに基づいて第１変形画像と第２画像を互いに対応する領域ごとに分割できるものであれば、いかなる方法でもよい。

例えば、分割条件として、変形した第１画像を所定の図形に対する所定の第１の空間パラメータに基づいて複数の第１分割画像に分割し、第２画像を第１の空間パラメータに対応する第２の空間パラメータに基づいて第１分割画像に対応する複数の第２分割画像に分割することを規定することが考えられる。なお、ここでいう図形とは、一定の決まりによって定められる形状を意味し、点、直線、平面、曲線、球などの立体またはその一部を含む。第１画像および第２画像の対応する図形は、点、直線、曲線、平面、曲面や球などの立体やこれらの一部など任意に定義してよいが、両図形は第１画像と第２画像の対応する位置に位置することが必要であり、さらに、計算の便宜上、同一の種類の図形であることが望ましい。また、第１および第２の空間パラメータは、変形した第１画像Ｖ１ａおよび第２画像間の対応する図形に対して同様の決まりにより規定される同じ種類のパラメータであればよく、１つであってもよく複数であってもよい。例えば、第１および第２の空間パラメータとして距離や角度などを用いることができる。

なお、上記第１分割画像と第２分割画像を対応する範囲を表すものとするために、第1分割画像および第２分割画像のスケールを統一する必要がある。このために、第１分割画像と第２分割画像のスケールを統一する処理を、ヘッダ情報などの任意の情報から周知の方法で取得した画素サイズ（Pixel Spacing）とスライス間隔(Slice Spacing)に基づいて、両画像の第１分割画像と第２分割画像の設定の前に、第１画像と第２画像に対して実施してもよく、第１分割画像と第２分割画像のそれぞれに対して実施してもよい。

なお、対応する図形の初期位置は周知の任意の方法で特定してよい。例えば第１画像と第２画像の対応する位置をユーザが入力した位置により特定してもよく、周知の自動認識技術により得られた解剖学構造物の特徴的な位置から特定してもよい。

第４の実施形態における分割条件は、変形した第１画像を所定位置からの距離に応じて複数の第１分割画像に分割し、第２画像を第１画像の所定位置に対応する位置からの距離に基づいて第１分割画像に対応する複数の第２分割画像に分割することを規定する。そして、第４の実施形態における評価関数Ｓ（μ）は、分割条件に従って、第１分割画像とこれに対応する第２分割画像とからなる一対の分割画像ごとに、一対の分割画像のそれぞれの画素値の分布の類似性を表す分割画像類似度を定義し、各分割画像類似度に基づいて、第１画像を変形した画像Ｖ１ａと第２画像Ｖ２の類似度を評価する。

ここで、評価関数Ｓ（μ）は、複数の分割画像類似度に基づいて変形した第１画像Ｖ１ａと第２画像Ｖ２との類似度を算出するものであれば、種々の方法を用いることができる。一例として、評価関数Ｓ（μ）は、複数の分割画像類似度の総和により定義できる。

第４の実施形態においては、評価関数Ｓ（μ）を用いた第１画像を変形した画像と第２画像の類似度の評価以外の処理については、第１の実施形態と同じであり、各機能ブロックの機能も共通している。以下、第１の実施形態と異なる点を中心に説明し、第１の実施形態と同じ点については説明を省略する。

図４Ａは第４の実施形態における第１および第２画像の分割方法を説明する図である。図４Ａにおける変形した第１画像Ｖ１ａおよび第２画像Ｖ２は、共に３次元画像の中心を原点としたｘｙ平面を表しており、変形した第１画像Ｖ１ａの各制御点Ｘを説明のために、ｘｙ平面上でのみ変位させて表している。図４Ａに示すように、本実施形態における分割条件は、変形した第１画像Ｖ１ａおよび第２画像Ｖ２のそれぞれの基準点Ｐ１、Ｐ２を設定し、各基準点Ｐ１、Ｐ２から距離に応じて、各基準点Ｐ１、Ｐ２からの距離の範囲に応じてＶ１ａおよび第２画像Ｖ２をそれぞれ分割することを規定している。詳細には、基準点からの距離に応じて、０≦ｄ＜ｄ_０、ｄ_０≦ｄ＜ｄ_１、ｄ_１≦ｄ＜ｄ_２、…、ｄ_ｋ−１≦ｄ＜ｄ_ｋを満たすｋ個の球状または中空の球状の領域に，変形した第１画像Ｖ１ａおよび第２画像Ｖ２をそれぞれ分割する。そして、評価関数Ｓ（μ）を、各第１分割画像と第２分割画像のペア（Ａ_０、Ｂ_０）、（Ａ_１、Ｂ_１）、…、（Ａ_ｋ、Ｂ_ｋ）ごとの分割画像類似度の総和によって規定する。

また、評価関数として、以下の式（７）を用いる。
式（７）において、ｄは、各距離の範囲を意味し、集合Ｄは第１画像の所定点からの距離の範囲ｄ_０、ｄ_１、ｄ_２、…、ｄ_ｋ（ｋは正の整数）からなる集合を表す。なお、場合により、距離の範囲０≦ｄ＜ｄ_０を距離の範囲ｄ_０と、記載する。また、各基準点の位置は、ユーザのマニュアル操作により入力装置４から入力されたものとする。また、式（７）に示す評価関数Ｓ（μ）においては、第１分割画像と第２分割画像のペア（Ａ_０、Ｂ_０）、（Ａ_１、Ｂ_１）、…、（Ａ_ｋ、Ｂ_ｋ）の類似度（分割画像類似度）を、各ペア画像の相互相関情報に先述の蓋然性情報により重みづけた項目により定義している。

図５は、第４の実施形態の画像処理装置の処理の流れを示すフローチャートである。図５に従って、処理の流れを説明する。まず、画像取得部１１は、被写体を撮像して得られた第１画像（第１画像データ）Ｖ１および第２画像（第２画像データ）Ｖ２を取得する（Ｓ１１）。

次に、画像変形量推定部１２は、予め設定された分割条件に従って、第２画像を分割した第２分割画像を設定する。具体的には、ユーザの入力装置の入力に基づいて第２画像Ｖ２の基準点Ｐ２を特定し、基準点Ｐ２と第２画像における各制御点Ｘとの距離をそれぞれ算出する。そして、基準点Ｐ２からの距離に応じて、所定の距離の範囲ｄ_０、ｄ_１、ｄ_２、…、ｄ_ｋごとに、第２画像Ｖ２を複数の第２分割画像として設定し、各分割画像を特定する情報をメモリに記憶する（Ｓ１２）。

次に、画像変形量推定部１２は、画像変形量μを設定して（Ｓ１３）、変換関数ｇを用いて第１画像Ｖ１を変形する（Ｓ１４）。

また、画像変形量推定部１２は、予め設定された分割条件に従って、第１画像を分割した第１分割画像を設定する。具体的には、ユーザの入力装置の入力に基づいて変形した第１画像Ｖ１ａの基準点Ｐ１を特定し、基準点Ｐ１と変形した第１画像Ｖ１ａにおける各制御点Ｘとの距離をそれぞれ算出する（Ｓ１５）。そして、基準点Ｐ１からの距離に応じて、所定の距離の範囲ｄ_０、ｄ_１、ｄ_２、…、ｄ_ｋごとに、変形した第１画像Ｖ１ａを複数の第１分割画像として設定し、各分割画像を特定する情報をメモリに記憶する。

次いで、画像変形量推定部１２は、変形した第１画像Ｖ１ａと第２画像Ｖ２各画素値の各画素値を取得して同時確率分布関数ｐ（ｆ,ｍ；μ）や周辺確率分布関数ｐ（ｆ）、ｐ（ｍ；μ）を算出し、蓋然性情報取得部１３は、算出された同時確率分布関数ｐ（ｆ,ｍ；μ）、周辺確率分布関数ｐ（ｆ）、ｐ（ｍ；μ）および蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）に基づいて、先述のように条件付き確率を算出することにより蓋然性情報Ｐ_Ｌ（ｍ；μ｜ｆ）を取得する（Ｓ１６）。そして、画像変形量推定部１２は、評価値として式（７）で表される評価関数Ｓ（μ）の変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜を算出する（Ｓ１７）。そして、算出された評価値Ｓの変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜が所定のしきい値より大きい場合には（Ｓ１８のＮ）、画像変形量μを所定量Δμ増加させたμ＋Δμを新しい画像変形量μとして設定し（Ｓ２１）、再度Ｓ１３からＳ１８の処理を繰り返す。なお、画像変形量推定部１２は、評価値として、評価関数Ｓ（μ）の変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜に換えて、式（７）であらわされる評価関数Ｓ（μ）の偏導関数の絶対値｜▽Ｓ（μ）｜を用いてもよい。

一方、評価関数Ｓ（μ）の変化量｜Ｓ（μ）−Ｓ（μ−Δμ）｜が所定のしきい値以下である場合には（Ｓ１８のＹ）、画像変形量推定部１２は、この場合の画像変形量μを評価関数Ｓ（μ）が最大となる場合の画像変形量μとして取得し、この画像変形量μに基づいて、第１画像を変形する変換関数を決定する。画像生成部１４は、決定された変換関数を用いて第１画像を変換して再構成し、画像Ｖ１Ａを生成する（Ｓ１９）。そして、表示制御部１５は、生成された画像Ｖ１Ａと第２画像Ｖ２を比較可能に並列表示する（Ｓ２０）。

第４の実施形態によれば、変形した第１画像Ｖ１ａと第２画像Ｖ２をそれぞれ対応する複数の分割画像に分割し、評価関数Ｓ（μ）が複数の第１分割画像と第２分割画像のペアごとに第１および第２分割画像間の画素値の分布の相関性をそれぞれ定義した複数の分割画像類似度に基づいて類似度を評価することにより、互いに異なる分割領域に位置する複数の被写体を表す画素値間で相関性が評価されることが抑制され、上記第２の問題の発生を低減することができる。結果として、さらに精度良く類似度を評価することができるため、より好適に第１画像を第２画像と一致するように変換した画像を生成することができる。

また、分割条件を基準点からの距離という１つの空間パラメータに応じて変形した第１画像および第２画像を分割するものとし、基準点Ｐ１、Ｐ２からの距離に応じて、所定の距離の範囲ｄ_０、ｄ_１、ｄ_２、…、ｄ_ｋごとに、対応する各第１分割画像および各第２分割画像を特定する情報をメモリに記憶した場合には、変形した第１画像および第２画像の類似度算出のための計算負荷を必要以上に増加させることがない。また、第１の空間パラメータを、基準となる点からの距離としたため、第１画像および第２画像において基準となる図形の設定が容易である。

第４の実施形態の分割条件のさらなる変形例として、距離に換えて角度により変形した第１画像Ｖ１ａおよび第２画像Ｖ２を分割してもよい。このように、分割条件を基準となる図形に対する角度θに応じて所定角度範囲ごとに、変形した第１画像Ｖ１ａおよび第２画像Ｖ２をそれぞれ分割することが考えられる。なお、この場合、式（７）における距離ｄを角度θに置換したものを評価関数として用いることができる。

図４Ｂに、第４の実施形態における第１および第２画像の分割方法の他の例として、変形した第１画像Ｖ１ａおよび第２画像Ｖ２をｘｙ平面（所定の平面）において、ｘ軸（所定の軸）からの角度θに基づいて分割した例を示す。図４Ｂにおける変形した第１画像Ｖ１ａおよび第２画像Ｖ２は、共に３次元画像の原点を通るｘｙ平面を表しており、変形した第１画像Ｖ１ａの各制御点Ｘを説明のために、ｘｙ平面上でのみ変位させて表している。なお、角度θは、両画像Ｖ１ａ、Ｖ２をそれぞれ円柱座標系で表した場合の角度θとして取得できる。

図４Ｂでは、０≦θ＜９０°、９０°≦θ＜１８０°、１８０°≦θ＜２７０°、２７０°≦θ＜３６０°の４つの範囲に、角度に応じて変形した第１画像Ｖ１ａを第１分割画像Ａ_０、Ａ_１、Ａ_２、Ａ_３に分割し、第２画像Ｖ２を第２分割画像Ｂ_０、Ｂ_１、Ｂ_２、Ｂ_３に分割している。θの範囲を任意に設定することにより、分割条件を基準となる図形に対する任意の角度θに応じて、０≦θ＜θ_０、θ_０≦θ＜θ_１、θ_１≦θ＜θ_２、…、θ_ｋ−１≦θ＜θ_ｋのｋ個に変形した第１画像Ｖ１ａおよび第２画像Ｖ２を分割するものとしてもよい。

上記のように、分割条件を基準図形からの角度という１つの空間パラメータに応じて変形した第１画像および第２画像を分割するものとし、基準図形（基準となるｘ軸）からの角度に応じて、所定の角度の範囲θ_０、θ_１、θ_２、…、θ_ｋごとに、対応する各第１分割画像および各第２分割画像を特定する情報をメモリに記憶した場合には、変形した第１画像および第２画像の類似度算出のための計算負荷を必要以上に増加させることがない。

また、ｘ軸に対する任意の角度（０≦θ＜θ_０、θ_０≦θ＜θ_１、θ_１≦θ＜θ_２、…、θ_ｋ−１≦θ＜θ_ｋ）に基づいて、変形した第１画像Ｖ１ａおよび第２画像Ｖ２をそれぞれ分割し、この分割した画像をさらに、ｚ軸との任意の角度（０≦β＜β_０、β_０≦β＜β_１、β_１≦β＜β_２、…、β_ｍ−１≦β＜β_ｍ）（ｍは０以上の整数）に基づいて、さらに分割してもよい。このように２つの軸に対するそれぞれの角度に基づいて変形した第１画像Ｖ１ａおよび第２画像Ｖ２を分割した場合には、変形した第１画像Ｖ１ａおよび第２画像Ｖ２を３次元的に近接した領域ごとに分割できるため、第１の問題をより好適に抑制でき、変形した第１画像および第２画像の類似度を精度よく評価できる。

また、変形した第１画像Ｖ１ａを、距離に基づいて分割するとともに、角度にも基づいてさらに分割して第１分割画像を設定してもよい。

上記第４の実施形態において、各第１分割画像とこれに対応する第２分割画像のペアについて、式（３）に表す評価関数によって分割画像類似度を算出しているが、これに限定されることなく、各第１分割画像とこれに対応する第２分割画像のペアについて、第２の実施形態における第２の評価関数によって分割画像類似度を算出してもよく、第３の実施形態における、第１の評価関数および第２の評価関数を組み合わせた評価関数によって分割画像類似度を算出してもよく、各第１分割画像とこれに対応する第２分割画像のペアの類似度を算出可能な種々の評価関数によって分割画像類似度を算出することができる。

上記の各実施形態はあくまでも例示であり、上記のすべての説明が本発明の技術的範囲を限定的に解釈するために利用されるべきものではない。

この他、上記の実施形態におけるシステム構成、ハードウェア構成、処理フロー、モジュール構成、ユーザインターフェースや具体的処理内容等に対して、本発明の趣旨から逸脱しない範囲で様々な改変を行ったものも、本発明の技術的範囲に含まれる。

また、画像処理装置１は、複数台のコンピュータにより、手段としての機能を分担する構成としてもよい。また、入力装置、ディスプレイ等、システムを構成する装置としては、公知のあらゆる装置を採用することができる。

１ 画像処理装置、 ２ ストレージ、 ３ ディスプレイ（表示装置）
４ 入力装置
１１ 画像取得部
１２ 画像変形量推定部
１３ 蓋然性取得部
１４ 画像生成部
１５ 表示制御部
Ｖ１ ＭＲ画像（第１画像）
Ｖ１ａ 第１画像を変形した画像
Ｖ２ ＣＴ画像（第２画像）
Ｓ 評価関数
Ｘ 制御点
μ 画像変形量

Claims (13)

1. 被写体を第１のモダリティで撮像して得られた第１画像と、前記被写体を前記第１画像とは異なる第２のモダリティで撮像して得られた第２画像とを取得する画像取得部と、
   前記第１画像を変形して該変形した第１画像と前記第２画像との類似度を、前記変形した第１画像の画素値と前記第２画像の対応する画素値の分布の相関性を評価する評価関数により評価して前記変形した第１画像と前記第２画像が類似する前記第１画像の画像変形量を推定する画像変形量推定部と、
   前記推定された画像変形量に基づいて、前記第１画像を変形させた画像を生成する画像生成部を備え、
   前記画像変形量推定部が、前記第１画像および前記第２画像のそれぞれの対応する画素値の組合せの確からしさを表す蓋然性情報を取得する蓋然性情報取得部を備えるものであり、
   前記評価関数が、前記変形した第１画像の画素値と前記第２画像の対応する画素値の相関性の尺度を表す項目を備え、該項目が前記取得した蓋然性情報に基づいて相関性を評価することを特徴とする画像処理装置。
2. 前記評価関数における前記相関性の尺度を表す項目が、前記蓋然性情報を用いて重み付けされていることを特徴とする請求項１記載の画像処理装置。
3. 前記蓋然性情報が、前記第２のモダリティでさらなる被写体を撮影して得られた第２参照画像の各画素値が得られるという事象のもとで前記第１のモダリティで前記さらなる被写体を撮影して得られた第１参照画像の各画素値が得られるという事象の発生する条件付き確率を表すものであることを特徴とする請求項２記載の画像処理装置。
4. 前記蓋然性情報が、前記条件付き確率による重み付けを調整するために、前記第１参照画像における画素値分布を一様分布と近似した確率密度関数の逆数によりさらに重み付けされたものであることを特徴とする請求項３項記載の画像処理装置。
5. 前記蓋然性情報が、被写体の種類ごとに、所定の種類の被写体を前記第１のモダリティで撮影した場合に得られる第１の画素値範囲と、前記所定の種類の被写体を前記第２のモダリティで撮影した場合に得られる第２の画素値範囲をそれぞれ対応付けたものであり、
   前記評価関数が、前記蓋然性情報に基づいて、前記変形した第１画像と前記第２画像の各画素値が、前記第１の画素値範囲と該第１の画素値範囲に対応付けられた前記第２の画素値範囲の両方を満たすものでない場合には、前記変形した第１画像と前記第２画像の類似度を低くするように重み付けされたものであることを特徴とする請求項２記載の画像処理装置。
6. 前記第１または第２の画素値範囲が、それぞれ前記第１または第２のモダリティの撮像原理に基づいて、前記被写体の種類ごとに、前記第１および第２の画素値範囲を推定して算出されたものであることを特徴とする請求項５記載の画像処理装置。
7. 前記評価関数における相関性の尺度を表す項目が、前記変形した第１画像の画素値および前記第２画像の画素値を離散確率変数とした相互相関情報量または二乗損失相互情報量を表すものであることを特徴とする請求項２から６のいずれか１項記載の画像処理装置。
8. 前記蓋然性情報が、前記第２のモダリティでさらなる被写体を撮影して得られた第２参照画像の各画素値が得られるという事象のもとで前記第１のモダリティで前記さらなる被写体を撮影して得られた第１参照画像の各画素値が得られるという事象の発生する条件付き確率を表しており、
   前記評価関数における前記相関性の尺度を表す項目が、前記変形した第１画像の画素値と前記第２画像の画素値とを離散確率変数とした同時確率分布と、前記蓋然性情報との差異を表していることを特徴とする請求項１記載の画像処理装置。
9. 前記評価関数が、下記式（９）により定義されることを特徴とする請求項８記載の画像処理装置。
   ただし、
   ｆ：前記第２画像の画素値
   ｍ：前記第１画像の画素値
   Ｆ：前記第２画像の全画素値からなる集合
   Ｍ：前記第１画像の全画素値からなる集合
   ｍ；μ：前記第１画像の変形量をμとしたときの変形した前記第１画像の画素値。
   ｐ（ｆ,ｍ；μ）：前記変形した第１画像の画素値と前記第２画像の画素値とを離散確率変数とした同時確率分布
   ｐ（ｆ）: 前記第２画像の画素値を離散確率変数とした周辺確率分布
   ｐ_Ｌ（ｍ；μ｜ｆ）：前記蓋然性情報
   とする。
10. 前記評価関数が、前記第１画像を変形して該変形した第１画像と前記第２画像との類似度を、前記変形した第１画像の画素値と前記第２画像の対応する画素値の分布の相関性を評価するさらなる評価関数を備え、
    前記さらなる評価関数は、前記変形した第１画像の画素値と前記第２画像の対応する画素値の相関性の尺度を表す項目を備え、該項目が前記変形した第１画像の画素値および前記第２画像の画素値を離散確率変数とした相互相関情報量または二乗損失相互情報量を表していることを特徴とする請求項８または９項記載の画像処理装置。
11. 前記さらなる評価関数において、前記変形した第１画像の画素値と前記第２画像の対応する画素値の相関性の尺度を表す項目が前記取得した蓋然性情報に基づいて重み付けされていることを特徴とする請求項１０項記載の画像処理装置。
12. 画像処理装置に実行させる画像処理方法であって、
    被写体を第１のモダリティで撮像して得られた第１画像と、前記被写体を前記第１画像とは異なる第２のモダリティで撮像して得られた第２画像とを取得する工程と、
    前記第１画像を変形して該変形した第１画像と前記第２画像との類似度を、前記変形した第１画像の画素値と前記第２画像の対応する画素値の分布の相関性を評価する評価関数により評価して前記変形した第１画像と前記第２画像が類似する前記第１画像の画像変形量を推定する工程と、
    前記推定された画像変形量に基づいて、前記第１画像を変形させた画像を生成する工程を備え、
    前記画像変形量を推定する工程が、前記第１画像および前記第２画像のそれぞれの対応する画素値の組合せの確からしさを表す蓋然性情報を取得する蓋然性情報取得工程を備えるものであり、
    前記評価関数が、前記変形した第１画像の画素値と前記第２画像の対応する画素値の相関性の尺度を表す項目を備え、該項目が前記取得した蓋然性情報に基づいて相関性を評価していることを特徴とする画像処理方法。
13. コンピュータに、
    被写体を第１のモダリティで撮像して得られた第１画像と、前記被写体を前記第１画像とは異なる第２のモダリティで撮像して得られた第２画像とを取得する工程と、
    前記第１画像を変形して該変形した第１画像と前記第２画像との類似度を、前記変形した第１画像の画素値と前記第２画像の対応する画素値の分布の相関性を評価する評価関数により評価して前記変形した第１画像と前記第２画像が類似する前記第１画像の画像変形量を推定する工程と、
    前記推定された画像変形量に基づいて、前記第１画像を変形させた画像を生成する工程とを実行させる画像処理プログラムであって、
    前記画像変形量を推定する工程が、前記第１画像および前記第２画像のそれぞれの対応する画素値の組合せの確からしさを表す蓋然性情報を取得する蓋然性情報取得工程を備えるものであり、
    前記評価関数が、前記変形した第１画像の画素値と前記第２画像の対応する画素値の相関性の尺度を表す項目を備え、該項目が前記取得した蓋然性情報に基づいて相関性を評価していることを特徴とする画像処理プログラム。