JP6934734B2 - 画像処理装置、画像処理装置の制御方法およびプログラム - Google Patents

Description

本発明は、画像処理装置、画像処理装置の制御方法およびプログラムに関し、特に、核磁気共鳴映像装置（ＭＲＩ）、Ｘ線コンピュータ断層撮影装置（Ｘ線ＣＴ）、超音波画像診断装置（ＵＳ）など、種々の撮像装置（モダリティ）で撮像した３次元画像を処理する装置に関する。ＭＲＩはＭａｇｎｅｔｉｃ Ｒｅｓｏｎａｎｃｅ Ｉｍａｇｉｎｇ、ＣＴはＣｏｍｐｕｔｅｄ Ｔｏｍｏｇｒａｐｈｙ、ＵＳはＵｌｔｒａＳｏｎｏｇｒａｐｈｙの略である。

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

また、同種の目的のために、複数の画像を位置合わせして表示することや、位置合わせした画像間の差分画像を表示することが試みられている。このとき、高精度な位置合わせを実現するために、ユーザが目視で同定した対応点を取得して位置合わせ処理に用いることがしばしばなされている。

非特許文献１には、位置合わせ対象のデータだけでなく、ユーザが目視で同定した特徴点の組（対応点）を位置合わせ処理に用いている。より具体的には、医用画像などに含まれる解剖構造の形状間において、夫々の形状の一致度とユーザが同定した対応点の一致度の両方の拘束条件を考慮したコスト関数を、非線形最適化手法で最適化することで位置合わせする手法が開示されている。このとき、物理的に不自然な変形になることを抑制するための正則化項をコスト関数に付加することが、一般的に行われている。

Ｘｉａｏｌｅｉ Ｈｕａｎｇ，Ｎｉｋｏｓ Ｐａｒａｇｉｏｓ，ａｎｄ Ｄｉｍｉｔｒｉｓ Ｎ．Ｍｅｔａｘａｓ，"Ｓｈａｐｅ Ｒｅｇｉｓｔｒａｔｉｏｎ ｉｎ Ｉｍｐｌｉｃｉｔ Ｓｐａｃｅｓ Ｕｓｉｎｇ Ｉｎｆｏｒｍａｔｉｏｎ Ｔｈｅｏｒｙ ａｎｄ Ｆｒｅｅ Ｆｏｒｍ Ｄｅｆｏｒｍａｔｉｏｎｓ"，ＩＥＥＥ ＴＲＡＮＳＡＣＴＩＯＮＳ ＯＮ ＰＡＴＴＥＲＮ ＡＮＡＬＹＳＩＳ ＡＮＤ ＭＡＣＨＩＮＥ ＩＮＴＥＬＬＩＧＥＮＣＥ，ＶＯＬ．２８，ＮＯ．８，ＡＵＧＵＳＴ ２００６．

しかし、非特許文献１の技術では、対応点の一致度以外の拘束条件の働きにより対応点位置における位置合わせ精度が不十分であるという課題があった。

本発明は、上記の課題に鑑みてなされたものであり、２つの画像データの間で、対応点を高精度に一致させつつ、全体として所望の位置合わせを行うための技術を提供することを目的とする。

上記の目的を達成する本発明に係る画像処理装置は、以下の構成を備える。即ち、
第１の画像と第２の画像との間で被検体の特徴点を対応付けた対応点を取得する対応点取得手段と、
前記対応点に基づく第１の変形を取得する第１の変形取得手段と、
前記第１の画像と前記第２の画像との間で画像の輝度情報を一致させることを含む第２の変形を取得する第２の変形取得手段と、
前記第１の画像と前記第２の画像との少なくとも一方の画像上の各位置について、前記対応点から該位置への距離に基づいて、前記第１の変形と前記第２の変形とを合成する合成手段と、
を備えることを特徴とする。

本発明によれば、２つの画像データの間で、対応点を高精度に一致させつつ、全体として所望の位置合わせを行うことが可能となる。

第１の実施形態に係る画像処理装置の機器構成を示す図である。 第１の実施形態における全体の処理手順を示すフローチャートである。 位置合わせ対象である第１の画像と第２の画像の一例を示す図である。 第１の画像と第２の画像の間の対応点の一例を示す図である。 画像上の所定位置における第１の変位場の合成比率のグラフを示す図である。 第２の実施形態における全体の処理手順を示すフローチャートである。 第３の実施形態における全体の処理手順を示すフローチャートである。 第４の実施形態に係る画像処理装置の機器構成を示す図である。 第４の実施形態における全体の処理手順を示すフローチャートである。

以下、図面を参照しながら実施形態を説明する。なお、以下の実施形態において示す構成は一例に過ぎず、本発明は図示された構成に限定されるものではない。

（第１の実施形態）
＜概要＞
本実施形態に係る画像処理装置は、３次元画像である第１の画像と第２の画像の変形位置合わせを行う。このとき、画像間で対応付けられた対応点を高精度に一致させる第１の変形と、所望の拘束条件に基づく位置合わせを行う第２の変形とを、対応点からの距離に応じた割合で合成する。より具体的には、対応点に近い位置ほど第１の変形の変位場の合成比率が高くなり（特に対応点上では合成比率＝１．０）、対応点からの距離が離れるにつれ、第２の変形の変位場の合成比率が高くなるように、２つの変位場を合成する。これにより、対応点に近いほど、対応点を高精度に一致させる第１の変形がより大きく影響し、対応点から離れるにつれて所望の拘束条件の第２の変形がより大きく影響する変形結果を取得できる。そのため、画像間で対応点は高精度に一致させつつ、全体として所望の位置合わせを行うことができる。なお、第１の変形は必ずしも対応点を一致させる変形である必要はなく、第２の変形よりも相対的に対応点の一致度が高い変形であればよい。また、第２の変形は、全体としての精度の期待値が第１の変形よりも高い変形であれば、いずれの変形であってもよい。

＜構成＞
以下、図１乃至図６を用いて、本実施形態の構成及び処理を説明する。

図１は、本実施形態に係る画像処理システム及び画像処理装置の構成例を示す。同図に示すように、本実施形態における画像処理システムは、画像処理装置１０００、データサーバ１１０、および表示部１２０を備えている。画像処理装置１０００は、データサーバ１１０および表示部１２０と接続されている。

データサーバ１１０は第１および第２の画像を保持する。ここで、第１および第２の画像は、異なる条件（異なるモダリティ、撮像モード、日時、体位等）で被検体を予め撮像して得られた３次元断層画像（ボリュームデータ）である。３次元断層画像を撮像するモダリティは、ＭＲＩ装置、Ｘ線ＣＴ装置、３次元超音波撮影装置、光音響トモグラフィ装置、ＰＥＴ／ＳＰＥＣＴ、ＯＣＴ装置などであってもよい。ＰＥＴは、Ｐｏｓｉｔｒｏｎ Ｅｍｉｓｓｉｏｎ Ｔｏｍｏｇｒａｐｈｙの略であり、ＳＰＥＣＴはＳｉｎｇｌｅ Ｐｈｏｔｏｎ Ｅｍｉｓｓｉｏn Ｃｏｍｐｕｔｅｄ Ｔｏｍｏｇｒａｐｈｙの略である。ＯＣＴはＯｐｔｉｃａｌ Ｃｏｈｅｒｅｎｃｅ Ｔｏｍｏｇｒａｐｈｙの略である。第１および第２の画像は、例えば、異なるモダリティや異なる撮像モードで同時期に撮像されたものであってもよい。また、経過観察のために同一患者を同一モダリティ、同一体位で異なる日時に撮像した画像であってもよい。第１および第２の画像は、データ取得部１０１０を介して画像処理装置１０００に入力される。

表示部１２０は液晶モニタ等であり、画像処理装置１０００が生成する表示画像等の各種の情報を表示する。また、表示部１２０には、ユーザからの指示を取得するためのＧＵＩ（Ｇｒａｐｈｉｃａｌ Ｕｓｅｒ Ｉｎｔｅｒｆａｃｅ）も配置されている。

画像処理装置１０００は、データ取得部１０１０と、対応点取得部１０２０と、第１の変形取得部１０３０と、第２の変形取得部１０４０と、変形合成部１０５０と、画像生成部１０６０と、表示制御部１０７０とを備えている。

データ取得部１０１０は、画像処理装置１０００へと入力される第１および第２の画像を取得する。対応点取得部１０２０は、第１および第２の画像上で対応する対応点の位置情報を取得する。第１の変形取得部１０３０は、第１の画像と第２の画像との間で対応点を高精度に一致させる第１の変形位置合わせを行い、第１の画像と第２の画像との間の変形を表す第１の変位場を取得する。第２の変形取得部１０４０は、所定の拘束条件に基づく第２の変形位置合わせを行い、第１の画像と第２の画像との間の変形を表す第２の変位場を取得する。変形合成部１０５０は、第１の変位場と第２の変位場とを対応点からの距離に応じた比率で合成した合成変位場を生成する。画像生成部１０６０は、第１の画像を合成変位場に基づいて変形させた変形画像や、第２の画像と変形画像との差分画像などを生成する。表示制御部１０７０は、第１の画像、第２の画像、変形画像、差分画像の断層画像などの情報を表示部１２０に表示させる制御を行う。

ここで、本実施形態において画像処理装置１０００で行われる画像の変形とその際に用いる変位場について説明する。本実施形態では、画像間で位置合わせする際の変形をΦと表記する。さらに、本実施形態では、変形Φによる位置合わせの適切さを評価する評価関数（コスト関数）を定義し、該評価関数を最小化するような変形Φを推定する。このとき推定された変形Φは変位場によって表現される。ここで、変位場の形態は、画像上の各位置を変位させる変位ベクトルからなる変位ベクトル場であってもよいし、空間的な場で表現されておらずとも、画像上の任意の位置を変位可能な所定の変換モデルが有する変換パラメータの形で表わされていてもよい。また、本実施形態では、変形Φとして、第１の画像を変形させて第２の画像に位置合わせするための変位場を推定するものとする。

このとき、第１の画像を変形させる変位場として、第１の画像の各座標において第２の画像へ向かう方向の変位量を有する変位場（この方向の変位場をΦとする）を算出すると、必ずしも変形先である第２の画像上の全ての座標に対応する変位量が取得されない。従って、この変位場Φを用いて第１の画像上の各画素値を第２の画像上にマッピングすると、穴あきの変形画像が生成されてしまう。そこで、本実施形態では、第１の画像を変形させる変位場として、第２の画像上の各座標において第２の画像から第１の画像へ向かう方向の変位量を有する変位場（上述の逆の方向の変位場Φ^−１）を算出する。これにより、逆の変位場Φ^−１を用いて第２の画像上の各座標に対応する第１の画像上の位置の画素値を参照して持ってきて、第２の画像上にマッピングすることで、穴のない第１の画像の変形画像を生成できる。すなわち、本実施形態では、第１の画像を第２の画像に位置合わせするための変形Φを表す変位場として、逆方向の変位場Φ^−１を求めるものとする。

＜処理＞
図２は、本実施形態に係る画像処理装置１０００が行う全体の処理手順を示すフローチャートである。

（Ｓ２０１０：データの取得）
ステップＳ２０１０において、データ取得部１０１０は、データサーバ１１０から、第１の画像と第２の画像を取得する。そして、第１および第２の画像を、対応点取得部１０２０、第２の変形取得部１０４０、画像生成部１０６０、および表示制御部１０７０へと出力する。

図３は、位置合わせ対象である第１の画像と第２の画像の一例を示す図である。図３の３０１０は第１の画像、３０２０は第２の画像を示している。それぞれの画像は例えばＣＴやＭＲＩ、３Ｄ超音波画像などの３次元のボリュームデータを表すが、３０１０、３０２０は、紙面が２次元である都合上、ボリュームデータからｘ−ｙ平面に平行に切り出された断面画像の形式で示されている。また、３０３０、３０４０は、それぞれ、第１の画像上の被検体の領域、第２の画像上の被検体の領域を示している。図３の例では、第１の画像と第２の画像は異なる種類のモダリティの撮像画像であり、被検体が人体における乳房である場合を示している。図３の例では、第１の画像と第２の画像の間で、被検体の領域（乳房）は大きく変形している。

（Ｓ２０２０：対応点の取得）
ステップＳ２０１０において、対応点取得部１０２０は、第１の画像と第２の画像の間における解剖学的に対応する特徴点（対応点）の位置情報を取得する。そして、取得した対応点の位置情報を第１の変形取得部１０３０、第２の変形取得部１０４０、そして、変形合成部１０５０へと出力する。

図４は、第１の画像と第２の画像の間の対応点の一例を示す図である。図４において、第１の画像３０１０上の点Ａ１〜Ａ３、および第２の画像３０２０上の点Ｂ１〜Ｂ３は、２つの画像間の対応点（乳頭や血管分岐部）であり、点Ａ１〜Ａ３は夫々点Ｂ１〜Ｂ３と対応している。なお、図４では、図示を簡潔にするため夫々の対応点が同一平面（断面）上に存在しているように図示しているが、実際の対応点は夫々の画像（第１の画像および第２の画像）が定義する３次元空間中に分布している。このような対応点の取得方法として、例えば、まず第１の画像、第２の画像に、それぞれインタレストオペレータ等の公知の画像特徴点抽出処理を施して、それぞれの画像から乳頭や血管分岐等の解剖学的な特徴点を抽出する。その後、抽出されたそれぞれの特徴点の間で、１対１に対応付けすることによって対応点の組を生成する処理を行う。これは、各画像上の特徴点ごとに特徴点近傍に関心領域を設定し、画像間で関心領域の画像類似度が高くなる特徴点同士を対応付けすることで行う。

なお、対応点の取得方法は、上記の方法に限られるものではなく、表示部１２０上に第１の画像と第２の画像を表示させ、ユーザが不図示の操作部を用いて手動で対応点を入力することで取得するようにしてもよい。また、第１の画像および第２の画像の付帯情報としてデータサーバ１１０が画像間の対応点群の座標情報を予め保持している場合には、これをデータサーバ１１０から入力する構成であってもよい。

（Ｓ２０３０：対応点を一致させる第１の変形の取得）
ステップＳ２０３０において、第１の変形取得部１０３０は、対応点取得部１０２０から取得した対応点を用いて、これらの対応点を高精度に一致させる第１の変形位置合わせ（第２の画像から第１の画像への変形を表す第１の変位場の推定）を実行する。このとき、第１の変形位置合わせによって推定された変形Φを、第１の変位場Φ_１ ^−１と定義する。そして、生成した第１の変位場Φ_１ ^−１を変形合成部１０５０へと出力する。

ここで、本ステップで定義する評価関数を、第１の評価関数Ｅ１（Φ）と称する。本実施形態では、画像間で対応点を一致させるために、第１の評価関数を構成する指標として、次式で表わされる、第１の画像と第２の画像の間における対応点の位置の乖離度項Ｌ（Φ）を適用する。

ここで、Φ_１ ^−１（ｘ）は、第１の変形Φ_１を表現する変位場Φ_１ ^−１によって第２の画像上の座標ｘを第１の画像に向けて変位させた座標を表す。また、（ｘ_Ａｉ，ｘ_Ｂｉ）は、第１の画像と第２の画像上におけるｉ番目の対応点ＡｉとＢｉの座標の組を表す。また、Ｎは対応点の数を表す。

第１の変形取得部１０３０は、第１の評価関数Ｅ１を最小化するように、変形Φ_１を表現する所定の変換モデルの変換パラメータを最適化する。これにより、画像内の被検体全体の位置合わせとして、対応点を一致させるという点で最適な変形が推定される。本実施形態では、変形Φ_１を表現する所定の変換モデルとして、評価関数Ｅ１の最適化により対応点を一致させる変形を導出可能な変換モデルを用いる。例えば、第２の画像上の夫々の対応点位置（すなわち、ｘ_Ｂｉ）を基準点とした放射基底関数で構成される変換モデルを適用する。ここで、放射基底関数とは、所定の基準点からの距離のみに依存して値が定まる関数を意味しており、Ｒａｄｉａｌ Ｂａｓｉｓ Ｆｕｎｃｔｉｏｎ（ＲＢＦ）とも表記する。放射基底関数に基づく変形推定手法は、（１）式の評価関数Ｅ１を最適化することで、理論上は画像間で対応点を完全に一致させられる（Ｅ１＝０にできる）という特徴を有している。実際には、計算上の誤差により完全には一致しないが、対応点間の位置誤差を限りなく０に近づけられるため、高精度に対応点を一致させることができる。本実施形態では、放射基底関数として、変形の曲げエネルギーを最小化できる、つまり滑らかな変形が得られるＴｈｉｎ Ｐｌａｔｅ Ｓｐｌｉｎｅ（ＴＰＳ）を適用する。但し、放射基底関数はこれに限らず、ガウス関数など、どのような放射基底関数であってもよい。さらに、適用する変換モデルは、対応点を一致させる評価関数Ｅ１を最適化可能であり、かつ高精度に対応点を一致させられる手法であれば、放射基底関数に限らず、どのような手法を用いてもよい。例えば、Ｆｒｅｅ Ｆｏｒｍ Ｄｅｆｏｒｍａｔｉｏｎ（ＦＦＤ）や、Ｌａｒｇｅ Ｄｅｆｏｒｍａｔｉｏｎ Ｄｉｆｆｅｏｍｏｒｐｈｉｃ Ｍｅｔｒｉｃ Ｍａｐｐｉｎｇ（ＬＤＤＭＭ）など、変形を表現できる公知の何れの変換モデルであってもよい。なお、第１の変形は必ずしも対応点を一致させる変形である必要はなく、後述する第２の変形よりも相対的に対応点の一致度が高い変形であればよい。例えば、第２の変形よりも正則化項の重みを小さくした第２の変形と同様な変形モデルによる変形を用いてもよい。

（Ｓ２０４０：所望の拘束条件に基づく第２の変形の取得）
ステップＳ２０４０において、第２の変形取得部１０４０は、対応点取得部１０２０から取得した対応点を用いて、所望の拘束条件に基づく第２の変形位置合わせ（第２の画像から第１の画像への変形を表す第２の変位場の推定）を実行する。このとき、第２の変形位置合わせによって推定された変形Φを、第２の変位場Φ_２ ^−１と定義する。そして、生成した第２の変位場Φ_２ ^−１を変形合成部１０５０へと出力する。

ここで、本ステップで定義する評価関数を、第２の評価関数Ｅ２（Φ）と称する。本実施形態では、変形における所望の拘束条件を、ステップＳ２０３０における対応点を一致させる拘束条件（評価関数Ｅ１（Φ））に加え、画像間で画像全体の輝度情報を一致させ、かつ物理的に不自然な変形にならないこととする。このような拘束条件を満たすために、第２の評価関数を構成する指標として、次式で表わされる、対応点の乖離度項Ｌ（Φ）（＝Ｅ１（Φ））、画像類似度項Ｓ（Φ）および正則化項Ｒ（Φ）からなる式を適用する。

ここで、画像類似度項をＳ（Φ）＝Ｓ（Φ（Ｉ_１），Ｉ_２）と定義する。ただし、Ｉ_１は第１の画像、Ｉ_２は第２の画像を表す。また、Φ（Ｉ_１）は、変位場Φ^−１を用いて画像Ｉ_１を変形させた変形画像を表す。画像類似度を計算する尺度としては、一般的に用いられているＳｕｍ ｏｆ Ｓｑｕａｒｅｄ Ｄｉｆｆｅｒｅｎｃｅ（ＳＳＤ）や相互情報量、相互相関係数などの公知の手法を用いるものとする。本実施形態では、例えば、相互情報量を適用する。また、Ｒ（Φ）は、変形Φが物理的に不自然にならないように、変形が所定の物理モデルから逸脱することに対してペナルティを与える正則化項である。正則化項としては、一般的に用いられている曲げエネルギー正則化や、体積保存正則化を用いるものとする。本実施形態では、例えば、体積保存正則化を適用する。また、λ_Ｓとλ_Ｒは夫々、類似度項Ｓ（Φ）と正則化項Ｒ（Φ）の重みを表し、本実施形態では、例えば、λ_Ｓ＝１．０、λ_Ｒ＝１．０を適用する。

なお、本実施形態における評価関数Ｅ２は、（２）式の形に限定されるものではなく、Ｌ（Φ）・Ｓ（Φ）・Ｒ（Φ）の何れかの項がなくてもよい。また、Ｒ（Φ）は例えば曲げエネルギー正則化と体積保存正則化など複数の正則化からなるものであってもよい。また、評価関数Ｅ２はユーザが所望する拘束条件を表現できるものであれば、どのような評価関数であってもよい。例えば、被検体の表面形状を一致させる形状項などを含むようにしてもよい。但し、後述するステップＳ２０５０において２種類の変位場を合成する性質上、評価関数Ｅ２は、評価関数Ｅ１すなわち対応点の乖離度項Ｌ（Φ）を含むものであることが望ましい（詳細な理由は後述する）。

第２の変形取得部１０３０は、第２の評価関数Ｅ２を最小化するように、変形Φ_２を表現する所定の変換モデルの変換パラメータを最適化する。これにより、画像内の被検体全体の位置合わせとして、所望の拘束条件を満たすという点で最適な変形が推定される。本実施形態では、任意の評価関数に基づいて変形推定が可能な手法であるＦＦＤを変換モデルとして適用し、非線形最適化処理によって評価関数の最適化を行う。

なお，適用する変換モデルやその最適化手法は、所望の拘束条件を付加した評価関数に基づいて変形推定が可能なモデルや手法であれば、何れのモデルや手法であってもよい。また、第２の変形は、全体としての精度の期待値が第１の変形よりも高い変形であれば、いずれの変形であってもよい。

（Ｓ２０５０：対応点からの距離に基づく第１と第２の変形の合成）
ステップＳ２０５０において、変形合成部１０５０は、対応点取得部１０２０から取得した対応点からの距離に基づき、第１の変位場Φ_１ ^−１と第２の変位場Φ_２ ^−１を合成した合成変位場Φ_１＋２ ^−１を生成する。そして、生成した合成変位場Φ_１＋２ ^−１を画像生成部１０６０へと出力する。

以下に、具体的な処理を説明する。変形合成部１０５０は、第２の画像上の各座標において、第２の画像上の対応点からの距離に基づいて、第１の変位場Φ_１ ^−１と第２の変位場Φ_２ ^−１を合成する。このとき、第１の変位場Φ_１ ^−１と第２の変位場Φ_２ ^−１は夫々、上述のとおり第２の画像上の各座標において、第１の画像へと向かう逆方向の変位場として取得されている。従って、以下の式のように、第２の画像上の各座標ｘにおいて、単純に第１の変位場Φ_１ ^−１と第２の変位場Φ_２ ^−１の変位量を所定の合成比率で重み付けして足し合わせるだけで、第２の画像上の各座標において合成された変位量を有する合成変位場を取得できる。

ここで、ｗ_１（ｘ）およびｗ_２（ｘ）は夫々、第２の画像上の座標ｘにおいて第１の変位場および第２の変位場を合成する合成比率を表している。両者はｗ_１（ｘ）＋ｗ_２（ｘ）＝１．０の関係を満たしている。なお、合成比率ｗ_１（ｘ）およびｗ_２（ｘ）は、次式のように、注目する座標ｘに対して、第２の画像上の各対応点Ｂｉ（１≦ｉ≦Ｎ）からの距離に基づいて算出される。

ここで、ｘ_Ｂｉは第２の画像上のｉ番目の対応点の座標を表わす。（４）式におけるｇ（ｘ｜μ，σ^２）は、座標μを平均、σ^２を分散とする正規化係数を除いたガウス関数であり、座標μにおいて最大値１をとり、μから離れるにつれμを中心に放射的に値が小さくなり０に近づく関数である。このガウス関数は、放射基底関数の一種である。ここで、放射基底関数は、ステップＳ２０３０に記載した通り所定の基準点からの距離のみに依存して値が定まる関数である。

従って、（４）式のｇ（ｘ｜ｘ_Ｂｉ，σ^２）ように、平均μが対応点座標ｘ_Ｂｉの場合は、対応点座標ｘ_Ｂｉで最大値１をとり、対応点座標ｘ_Ｂｉから離れるにつれ値が小さくなり０に近づいていく放射基底関数となる。本実施形態では、例えばσ＝１．０を適用する。さらに、１−ｇ（ｘ｜ｘ_Ｂｉ，σ^２）のように、この放射基底関数を１から減算することで、対応点座標ｘ_Ｂｉで最小値０をとり、点ｘ_Ｂｉから離れるにつれ値が大きくなり最大値１に近づいていく放射基底関数となる。これを放射基底関数ｒ（ｘ｜ｘ_Ｂｉ）と定義する。そして、（４）式のように合成比率ｗ_２（ｘ）として、放射基底関数ｒ（ｘ｜ｘ_Ｂｉ）を対応点の数Ｎだけ乗算することによって統合すると、全ての対応点の位置において必ず最小値０をとり、何れの対応点からも離れるにつれて最大値１に近づくような関数（合成比率関数）が得られる。これは、ある座標ｘにおいて複数の関数を掛け合わせるとき、何れか一つでも関数の値が０であれば結果も必ず０になり、複数の関数全ての最大値が１であれば、それらを掛け合わせた関数の最大値も１になるという性質があるためである。このようにして得られた合成比率ｗ_２（ｘ）は、対応点からの距離のみに応じて値が決まる夫々の放射基底関数ｒ（ｘ｜ｘ_Ｂｉ）を乗算した関数なので、ｗ_２（ｘ）も夫々の対応点からの距離に応じて値が決まる関数といえる。

なお、本実施形態では、（４）式において関数ｇとしてガウス関数を用いたが、これに限られるものではなく、対応点上で最大値１をとり、対応点からの距離が離れるにつれ値が減少する放射基底関数であれば何であってもよい。また、ｗ_２（ｘ）は、全ての対応点の位置において必ず０をとり、何れの対応点からも離れるにつれて１に近づく関数であれば何でもよく、必ずしも（４）式で表現される必要はない。最後に、（５）式のように合成比率ｗ_１（ｘ）として、１から合成比率ｗ_２（ｘ）を減算することで、全ての対応点の位置において必ず１をとり、何れの対応点からも離れるにつれて０に近づくような関数が得られる。従って、ｗ_１（ｘ）も放射基底関数であり、夫々の対応点からの距離に応じて値が決まる関数である。

なお、本実施形態では、合成比率関数は、各放射基底関数の乗算で表される関数としたが、これに限られるものではなく、放射基底関数を統合して得られるものであれば何であってもよい。例えば、各対応点上で全ての放射基底関数を重み付けした上で加算したときの和の値が０をとるような関数を適用することができる。これは、以下のようにして求めることができる。まず、各放射基底関数を加算する際の重みを未知の値、全ての放射基底関数の重み付け和によって表される合成比率関数の対応点上の値を既知の値０とする。そして、（対応点における全ての放射基底関数の重み付け和）＝０となるような連立方程式を対応点の数だけ立てて解くことで、未知数である全ての放射基底関数の重みを算出することができる。そして、合成比率関数は、このようにして算出された重みを有する放射基底関数の和で表すことができ、対応点上で値０をとる。なお、このまま合成比率関数を採用すると、放射基底関数の和で表されるため、対応点から離れた場所では値が１を超えてしまう。そこで最後に、合成比率関数全体を放射基底関数の数（対応点の数）Ｎで割って正規化することで、合成比率関数の最大値を１とすることができる。なお、対応点上では合成比率関数の値は０であるため、正規化後も値は０となることは保障される。

図５は、画像上の所定位置における第１の変位場の合成比率のグラフを示す図である。図５において、５０１０は、第２の画像３０２０上のあるｙ，ｚ座標におけるｘ方向の第１の変位場の合成比率ｗ_１（ｘ）を表わすグラフである。５０２０は、グラフ５０１０に対応する、あるｙ，ｚ座標におけるｘ方向のラインを表わす。さらに、グラフ５０１０上の曲線５０３０は、合成比率ｗ_１（ｘ）の関数の形状を表わす。そして、曲線５０３０上の点ｂ１・ｂ２・ｂ３は夫々、ｘ座標が対応点Ｂ１・Ｂ２・Ｂ３に対応する合成比率ｗ_１（ｘ）の値を表わす。この図において、曲線５０３０上の点ｂ２に関して、対応する点Ｂ２がちょうどライン５０２０上に位置するため、点Ｂ２に対応する放射基底関数ｒ（ｘ｜ｘ_Ｂ２）が値０となり、ｗ_２（ｘ_Ｂ２）も０となる。従って、点ｂ２における値がｗ_１（ｘ_Ｂ２）＝１．０をとっていることが分かる。また、点ｂ３に関しては、対応する点Ｂ３がライン５０２０上にはないが近くに位置する。そのため、ライン５０２０上の位置が点Ｂ３に対応する放射基底関数ｒ（ｘ｜ｘ_Ｂ３）から受ける影響が比較的大きくなり、点ｂ３における値ｗ_１（ｘ_Ｂ３）も大きな値をとっていることが分かる。一方、点ｂ１に関しては、対応する点Ｂ１がライン５０２０から遠く離れている。そのため、ライン５０２０上の位置が点Ｂ１に対応する放射基底関数ｒ（ｘ｜ｘ_Ｂ１）から受ける影響は多少はあるが小さくなり、点ｂ１における値ｗ_１（ｘ_Ｂ１）も小さな値をとっていることが分かる。このように、ライン５０２０上の位置において、何れかの対応点に近い位置ほど（点Ｂ２付近など）第１の変位場の合成割合が大きくなり、何れの対応点からも離れた位置ほど（ライン５０２０の中央付近など）第２の変位場の合成割合が大きくなることが分かる。実際には、通常対応点は３次元の画像空間中に疎に分布しているため、対応点近傍は対応点を高精度に一致させる第１の変位場の影響が大きく、それ以外の大部分の領域は所望の拘束条件に基づく第２の変位場の影響が大きくなる。従って、対応点は高精度に一致させ、画像全体としては所望の位置合わせを実現することができる。

なお、第２の変位場は、所望の拘束条件に基づくものであるため、位置合わせのための評価関数Ｅ２はステップＳ２０２０の（２）式となるように必ずしも対応点の乖離度項Ｌ（Φ）を含まなくともよい。しかしながら、下記の理由から、Ｌ（Φ）を含むようにした方がより適切な位置合わせを実現できる。すなわち、評価関数Ｅ２が対応点の乖離度項Ｌ（Φ）を含まない場合、第２の変位場による変形は対応点の残差が大きく発生する可能性がある。従って、第２の変位場と対応点を高精度で一致させる第１の変位場との間で、対応点の周辺における変位場の相違が大きく発生してしまう。その結果、両者を合成した合成変位場は、対応点の周辺での変位場の変化が激しくなり、不自然な変形結果になってしまう可能性がある。評価関数Ｅ２にもＬ（Φ）を含むようにすることで、第２の変位場は、対応点間の残差は発生するものの、対応点が概ね一致するような変位場にすることができる。従って、対応点近傍における第１の変位場との変位場の違いがあまり発生しないため、２つの変位場を合成した結果、より自然な位置合わせを行うことができる。

（Ｓ２０６０：画像の生成）
ステップＳ２０６０において、画像生成部１０６０は、変形合成部１０５０から取得した合成変位場に基づいて、データ取得部１０１０から取得した第１の画像を変形させた変形画像を生成する。また、データ取得部１０１０から取得した第２の画像と生成した変形画像との間で減算を行った差分画像を生成する。本実施形態では、差分画像は、第２の画像から変形画像を減算することで生成する。なお、必要に応じて、第１の画像の形状に一致するように第２の画像を変形させた変形画像（逆変形画像）を、合成変位場に基づいて生成することもできる。また、減算する方向はその逆であってもよい。また、変形画像と差分画像の何れか一方を生成しなくてもよい。そして、画像生成部１０６０は、生成した変形画像と差分画像とを表示制御部１０７０へと出力する。また、生成した画像をデータサーバ１１０へと出力して保存する。

（Ｓ２０７０：画像の表示）
ステップＳ２０７０において、表示制御部１０７０は、ユーザによる不図示のＵＩの操作に応じて、取得した画像を表示部１２０へと表示する制御を行う。具体的には、データ取得部１０１０から取得した第１の画像と第２の画像、および画像生成部１０６０から取得した変形画像と差分画像の断面画像を表示する制御を行う。このとき、第１の画像を変形させた変形画像と、第２の画像とを並べて表示するように構成することができる。さらに、第１の画像を加えて、第２の画像および変形画像とともに並べて表示する構成であってもよい。また、変形画像と差分画像は、ユーザの操作に応じて、１つの表示画面領域内で切り替えて表示できるような構成であってもよい。なお、変形画像や差分画像の表示を行わない構成であってもよい。以上によって、画像処理装置１０００の処理が実施される。

以上説明したように、本実施形態によれば、対応点に近いほど、対応点を高精度に一致させる第１の変形がより大きく影響し、対応点から離れるにつれて所望の拘束条件の第２の変形がより大きく影響した変形結果を取得できる。そのため、画像間で対応点は高精度に一致させつつ、画像全体として所望の拘束条件に従った位置合わせを行うことができる。

（変形例１−１：第１の画像上の対応点からの距離に基づいて変位場を合成してもよい）
第１の実施形態では、ステップＳ２０５０において、第２の画像上の対応点からの距離に基づいて、第１の変位場と第２の変位場を合成していた。しかし、第１の画像上の対応点からの距離に基づいて、２つの変位場を合成するようにしてもよい。ここで、本実施形態における変形Φは、ステップＳ２０２０に記載したとおり第２の画像上の各座標について第２の画像から第１の画像へと向かう方向の変位量を有する逆方向の変位場Φ^−１で表わされる。そこで、この変位場Φ^−１を用いて第２の画像上の各座標に対応する第１の画像上の位置を参照することで、合成比率を算出する。すなわち、以下の式のように変位場Φ^−１を介して第２の画像上の座標に対応する第１の画像上の各参照位置における第１の画像上の各対応点Ａｉ（１≦ｉ≦Ｎ）からの距離に基づいた比率を算出できる。

ここで、ｗ_２（ｘ）は、第１実施形態の（４）式の代わりに、第２の画像上の座標ｘにおいて第２の変位場を合成する合成比率を表している。また、Φ^−１（ｘ）は、変位場Φ^−１を用いて第２の画像上の座標ｘを第１の画像に向けて変位させた座標を表す。ここで、座標ｘを第１の画像上の座標Φ^−１（ｘ）に変換するために、変位場Φ^−１として何らかの変位場を与える必要があるが、当然、合成比率を算出する前の段階では、目的とする合成変位場は取得されていない。従って、ここでは、Φ^−１（ｘ）として、第１の変位場Φ_１ ^−１（ｘ）を適用する。但し、Φ^−１（ｘ）はこれに限られるものではなく、第２の変位場Φ_２ ^−１（ｘ）であってもよい。さらに、第１の変位場Φ_１ ^−１（ｘ）と第２の変位場Φ_２ ^−１（ｘ）とを決まった比率（例：同じ比率）で合成した変位場であってもよい。

さらに、第１の画像上の対応点からの距離と、第２の画像上の対応点からの距離との両方に基づいて、変位場を合成してもよい。このとき、合成比率は以下の式で算出できる。

ここで、ｗ_２（ｘ）は、第１実施形態の（４）式や上述の（６）式の代わりに、第２の画像上の座標ｘにおいて第２の変位場を合成する合成比率を表している。また、ａは、第１の画像上の対応点からの距離を利用する比率、ｂは、第２の画像上の対応点からの距離を利用する比率を表しており、ａ＋ｂ＝１．０を満たす。ここでは、例えば、ａ＝０．５、ｂ＝０．５を適用する。

以上説明したように、本変形例によれば、第１の画像上の対応点からの距離に基づいて変位場を合成することで、第１の画像上の対応点の位置関係を考慮した合成変位場を取得できる。

（変形例１−２：対応点の信頼度に基づいて合成比率を算出してもよい）
第１の実施形態では、対応点上で必ずｗ_１（ｘ）＝１となるように合成比率を算出する、つまり、対応点が必ず合致するように２つの変位場を合成していた。しかし、対応点上で必ずしもｗ_１（ｘ）＝１とならないように合成比率を算出してもよい。例えば、対応点ごとに、その対応関係に関する信頼度を設定し、信頼度に応じて合成比率を算出してもよい。より具体的には、対応点の信頼度が高ければ、ｗ_１（ｘ）の値を１に近付けることでより対応点が一致するように２つの変位場を合成する。一方、対応点の信頼度が低ければ、ｗ_１（ｘ）の値を０に近付けることで、より対応点が一致しなくてもいいように２つの変位場を合成する。このとき、合成比率は以下の式で表わされる。

ここで、ｗ_２（ｘ）は、第１実施形態の（４）式の代わりに、第２の画像上の座標ｘにおいて第２の変位場を合成する合成比率を表している。また、Ｃｉは、０≦Ｃｉ≦１．０を満たす比率であり、対応点ごとの信頼度を表わしている（１≦ｉ≦Ｎ）。（８）式において、対応点上で最大値１．０をとるガウス関数ｇ（ｘ｜ｘ_Ｂｉ，σ^２）に信頼度Ｃｉを掛け合わせた関数は、対応点上で最大値Ｃｉをとるような放射基底関数となる。従って、信頼度Ｃｉが低くなると、第２の変形の合成比率ｗ_２（ｘ）のｉ番目の対応点周辺の最小値が０．０よりも増加するようになる、つまり、第１の変形の合成比率ｗ_１（ｘ）のｉ番目の対応点周辺の最大値が１．０よりも減少するようになる。このように、対応点ごとにその信頼度に応じた合成比率を設定することができる。

信頼度Ｃｉは、例えばユーザが指定した対応点ならば高く設定し、ステップＳ２０２０に記載したようにコンピュータが自動検出した対応点ならば低く設定する。例えば、ユーザが指定した点ならばＣｉ＝１．０、コンピュータが自動検出した点ならば０≦Ｃｉ≦１．０を設定する。さらに、コンピュータが自動検出した点は、ステップＳ２０２０で特徴点間を対応点として対応付ける際に算出した、画像間の特徴点近傍の関心領域の画像類似度の大きさに応じて値を決めることができる。例えば、特徴点ごとの画像間の画像類似度をＳｉ（０≦Ｓｉ≦１．０）と定義すると、Ｃｉ＝Ｓｉとして信頼度Ｃｉを算出する。また、信頼度Ｃｉは、例えば、画像間で取得されている対応点の信頼性が被検体の表面形状から離れるにつれて低下するような場合に、被検体の表面形状からの距離に応じて算出してもよい。その場合、信頼度Ｃｉは、非検体の表面形状からの距離に反比例する値にすることで、表面形状に近いほど信頼度Ｃｉが高く、遠いほど信頼度Ｃｉが低くなるように値を設定することができる。また、ユーザが対応点を指定する場合において、夫々の対応点の対応の確信度を不図示のＵＩを介してユーザが指定できるようにして、該確信度に基づきＣｉを定めるようにしてもよい。確信度の与え方は、例えば、０〜１の値を無段階で与える形式を用いることができる。この場合、確信度をそのままＣｉとして用いることができる。このとき、確信度が一定以下の対応点については、式（８）の積の対象から除外するようにしてもよい（Ｃｉ＝０とするのと同値となる）。また、「確実に一致させる」「概ね一致させる」といったラベルのみを不図示のＵＩを介してユーザが夫々の対応点に対して指定し、「確実に一致させる」対応点はＣｉ＝１．０、「概ね一致させる」対応点は０≦Ｃｉ≦１．０の所定の値を与えるようにしてもよい。

以上説明したように、本変形例によれば、対応点の信頼度に応じて合成比率を算出することで、信頼性が高い対応点の周辺は対応点を一致させる変形に近づけ、信頼性が低い対応点の周辺は所望の拘束条件に基づく変形に近づけることができ、適切な位置あわせが実現できる。

（第２の実施形態）
＜概要＞
第１の実施形態では、ステップＳ２０５０において、第２の画像上の各座標において、第２の画像上の全ての対応点からの距離に基づいて変位場の合成比率を算出していた。このような場合、以下のような問題が生じることがある。すなわち、画像間で対応点の数が多数ある場合には、（３）式における放射基底関数ｒ（ｘ｜ｘ_Ｂｉ）の数が増えることで、それを乗算した関数ｗ_２（ｘ）が、各対応点の近傍において値０を中心に非常に急峻に１まで変化してしまう。つまり、各対応点のごく近傍で合成比率が急激に変化するため、合成される変位場も急激に変化してしまい、各対応点の近傍で不自然な変形が生じる可能性がある。

そこで、本実施形態では、第２の画像上の各座標において、その座標の近傍に存在する第２の画像上の対応点からの距離のみに基づいて変位場の合成比率を算出する。これにより、画像間で対応点の数が多くとも、画像上の各位置において合成比率の算出に利用される対応点はその近傍に存在する点に限られるため、合成比率の関数は急峻になることなく、２つの変位場が滑らかに合成されて自然な変形位置合わせが実現できる。

＜構成＞
以下、本実施形態に係る画像処理装置について、第１の実施形態との相違部分についてのみ説明する。

本実施形態に係る画像処理装置１０００の構成は、第１の実施形態と同様であるため、説明を省略する。しかし、変形合成部１０５０に関しては、第１の実施形態と機能が異なるため、以下でその機能を説明する。その他の構成については第１の実施形態と機能が同じであるため、同じ機能に関しては説明を省略する。

変形合成部１０５０は、第２の画像上の各座標において、第１の変位場と第２の変位場をその座標の近傍に存在する対応点からの距離のみに応じた比率で合成した合成変位場を生成する。

＜処理＞
図６は、本実施形態に係る画像処理装置１０００が行う全体の処理手順を示すフローチャートである。但し、図６において、ステップＳ６０１０からＳ６０４０、およびＳ６０６０からＳ６０７０の処理はそれぞれ、図２のフローチャートにおけるステップＳ２０１０からＳ２０４０、およびＳ２０６０からＳ２０７０の処理と同じであるため説明は省略する。以下、図２のフローチャートとの相違部分についてのみ説明する。

（Ｓ６０５０：各座標の近傍の対応点からの距離のみに基づく第１と第２の変形の合成）
ステップＳ６０５０において、変形合成部１０５０は、対応点取得部１０２０から取得した対応点からの距離に基づき、第１の変位場Φ_１ ^−１と第２の変位場Φ_２ ^−１を合成した合成変位場Φ_１＋２ ^−１を生成する。但し、本ステップでは、以下の点で、第１の実施形態のステップＳ２０５０とは異なる。ステップＳ２０５０では、変位場を合成する座標である第２の画像上の各座標において、全ての対応点からの距離に基づいて変位場の合成比率を算出していた。一方、本ステップでは、第２の画像上の各座標において、その座標の近傍に存在する第２の画像上の対応点を選択し、当該対応点からの距離のみに基づいて変位場の合成比率を算出する。そして、生成した合成変位場Φ_１＋２ ^−１を画像生成部１０６０へと出力する。

以下、第１の実施形態のステップＳ２０５０と相違する点について、具体的な処理を説明する。本ステップでは、第２の画像上の各座標ｘにおける近傍Ｋ個の（Ｋ近傍の）対応点からの距離のみに基づいて合成比率を算出する。具体的には、以下の式によって、合成比率を算出する。

ここで、ｗ_２（ｘ）は、第１実施形態の（４）式の代わりに、第２の画像上の座標ｘにおいて第２の変位場を合成する合成比率を表している。また、座標ｘのＫ近傍の対応点を点Ｂｋ（ｘ）（１≦ｋ≦Ｋ）と表記し、ｘ_Ｂｋ（ｘ）は点Ｂｋ（ｘ）の座標を表す。本実施形態では例えば、Ｋ＝５を適用する。

なお、第２の画像上の各座標ｘの近傍の対応点からの距離に基づく合成比率の算出方法は、上記の方法に限らず、例えば、座標ｘを中心とする一定範囲内（例：半径５ｍｍなど）に存在する対応点からの距離のみに基づいて算出するようにしてもよい。あるいは、座標ｘから一定範囲内にあり、かつ、所定の数の近傍対応点を用いて算出してもよい。これにより、合成比率ｗ_２（ｘ）は、画像間の全体の対応点の数に左右されずに、座標ｘの近傍の対応点からの距離のみに基づいて算出することができる。以上によって、本実施形態における画像処理装置１０００の処理が実施される。

以上説明したように、本実施形態によれば、画像間の対応点の数が多くとも、画像上の各位置において合成比率の算出に利用される対応点はその近傍に存在する点に限られる。そのため、合成比率の関数は急峻になることなく、２つの変位場が滑らかに合成されてより自然な変形位置合わせが実現できる。

（第３の実施形態）
＜概要＞
第１の実施形態では、ステップＳ２０５０において、分散が所定の値σに固定されたガウス関数に基づいて変位場の合成比率を算出していた。このような場合、以下のような問題が生じることがある。すなわち、第２の変形を表す第２の変位場において対応点残差が大きく発生した場合、第１の変形では高精度に対応点を一致させるため、対応点の周辺での変位場の相違が大きく発生してしまう。このとき、第１の実施形態では所定の分散σに固定されたガウス関数に基づいて合成比率を算出すると、対応点を中心とする合成比率の関数の広がりが固定されている。従って、固定されたσの値が小さい場合は、相違が大きい２つの変位場を、空間的に狭い範囲内で比率が急激に変化するように合成した結果、対応点の周辺での変位場の変化が大きくなってしまい、不自然な変形結果になる可能性がある。一方、固定されたσの値が大きい場合は、たとえ２つの変位場の相違が小さくても、空間的に広い範囲内で比率が変化するように合成した結果、不必要に広い範囲にまで第１の変位場の影響を大きくしてしまう可能性がある。つまり、所望の拘束条件に基づく第２の変位場の影響が小さくなってしまう可能性がある。

そこで、本実施形態では、第２の変位場における夫々の対応点残差の大きさに応じて、それらの対応点に対応するガウス関数のσの大きさを変えて、合成比率を算出する。これにより、対応点残差が大きい対応点周辺は比較的なだらかに広範囲で比率が変化するように変位場を合成し、対応点残差が小さい対応点周辺は比較的急峻に狭い範囲で比率が変化するように変位場を合成でき、より自然な変形結果を得ることができる。

＜構成＞
以下、本実施形態に係る画像処理装置について、第１の実施形態との相違部分についてのみ説明する。

本実施形態に係る画像処理装置の構成は、第１の実施形態と同様であるため、説明を省略する。しかし、第２の変形取得部１０４０および変形合成部１０５０に関しては、第１の実施形態と機能が異なるため、以下でその機能を説明する。その他の構成については第１の実施形態と機能が同じであるため、同じ機能に関しては説明を省略する。

第２の変形取得部１０４０は、第１の実施形態と同様に、第１の画像を変形させるための第２の変位場を取得する。さらに、本実施形態における第２の変形取得部１０４０は、第２の変位場における対応点残差を取得する（対応点残差取得処理）。変形合成部１０５０は、第１の変位場と第２の変位場を対応点からの距離と対応点残差の大きさに応じた比率で合成した合成変位場を生成する。

＜処理＞
図７は、本実施形態に係る画像処理装置１０００が行う全体の処理手順を示すフローチャートである。但し、図７において、ステップＳ７０１０からＳ７０４０、およびＳ７０６０からＳ７０７０の処理はそれぞれ、図２のフローチャートにおけるステップＳ２０１０からＳ２０４０、およびＳ２０６０からＳ２０７０の処理と同じであるため説明は省略する。以下、図２のフローチャートとの相違部分についてのみ説明する。

（Ｓ７０４５：第２の変形における対応点残差の取得）
ステップＳ７０４５において、第２の変形取得部１０４０は、ステップＳ７０４０で取得した第２の変位場Φ_２ ^−１における対応点残差を算出する。そして、算出した対応点残差を変形合成部１０５０へと出力する。

以下に、具体的な処理を説明する。第２の変位場Φ_２ ^−１における各対応点の対応点残差Ｅ_Ｌ（ｉ）（１≦ｉ≦Ｎ）を、以下の式に基づいて算出する。

ここで、Φ_２ ^−１（ｘ）は、第２の変形Φ_２によって第２の画像上の座標ｘを第１の画像に向けて変位させた座標を表す。従って、各対応点の対応点残差Ｅ_Ｌ（ｉ）は、第２の画像上の各対応点座標ｘ_Ｂｉを第２の変位場Φ_２ ^−１を用いて第１の画像上の座標に変換した座標Φ_２ ^−１（ｘ_Ｂｉ）と、第１の画像上の各対応点座標ｘ_Ａｉとの間の距離を表す。

（Ｓ７０５０：対応点からの距離と対応点残差に基づく 第１と第２の変形の合成 ）
ステップＳ７０５０において、変形合成部１０５０は、対応点取得部１０２０から取得した対応点からの距離と、対応点残差Ｅ_Ｌ（ｉ）とに基づき、第１の変位場Φ_１ ^−１と第２の変位場Φ_２ ^−１を合成した合成変位場Φ_１＋２ ^−１を生成する。そして、生成した合成変位場Φ_１＋２ ^−１を画像生成部１０６０へと出力する。

以下、第１の実施形態のステップＳ２０５０と相違する点について、具体的な処理を説明する。本実施形態では、第２の変位場における夫々の対応点残差Ｅ_Ｌ（ｉ）の大きさに応じて、それらの対応点に対応するガウス関数のσの大きさを変えて、合成比率を算出する。このとき、合成比率は以下の式ように算出される。

ここで、ｗ_２（ｘ）は、第１実施形態の（４）式の代わりに、第２の画像上の座標ｘにおいて第２の変位場を合成する合成比率を表している。（１１）式におけるσ（ｉ）は、対応点残差Ｅ_Ｌ（ｉ）の値に応じて決まる変数であり、（１２）式のように、対応点残差Ｅ_Ｌ（ｉ）に比例する。（１１）式において、αは定数であり、本実施形態では例えばα＝５．０を適用する。これにより、対応点残差Ｅ_Ｌ（ｉ）が大きいときは、ガウス関数ｇの分散σ（ｉ）^２も大きくなり、対応点残差Ｅ_Ｌ（ｉ）が小さいときは、ガウス関数ｇの分散σ（ｉ）^２も小さくなる。従って、対応点残差Ｅ_Ｌ（ｉ）の大きさに応じて、合成比率ｗ_２（ｘ）の各対応点周辺の関数の空間的な広がりが異なるようになる。具体的には、対応点残差Ｅ_Ｌ（ｉ）が大きい対応点周辺では合成比率ｗ_２（ｘ）の空間的広がりの度合いが大きく決定され、対応点残差Ｅ_Ｌ（ｉ）が小さい対応点周辺では合成比率ｗ_２（ｘ）の空間的広がりの度合いが小さく決定される。以上によって、本実施形態における画像処理装置１０００の処理が実施される。

以上説明したように、本実施形態によれば、対応点残差が大きい対応点周辺では、２つの変位場の相違が大きいため、比較的なだらかに広範囲で比率が変化するように変位場を合成する。一方、対応点残差が小さい対応点周辺では、２つの変位場の相違が小さいため、比較的急峻に狭い範囲で比率が変化するように変位場を合成する。これにより、より自然な変形結果を得ることができる。

（第４の実施形態）
＜概要＞
第１の実施形態では、対応点からの距離に応じた比率で、第１の変形と第２の変形を合成していた。しかしながら、距離に応じて第１の変形と第２の変形を合成するための基準は対応点に限られるものではない。本実施形態では、画像上の被検体の表面形状からの距離に応じた比率で第１の変形と第２の変形を合成する。これにより、画像間で取得されている対応点の信頼性が被検体の表面形状から離れるにつれて低下するような場合に、以下のような効果が得られる。すなわち、被検体の表面形状からの距離が大きくなるにつれ、対応点を高精度に一致させる変形よりも所望の拘束条件に基づく変形の比率が大きくなるように合成することで、対応点の信頼性に応じた適切な変形位置合わせを行うことができる。また、対応点が被検体の表面形状の近傍に主に分布する場合には、体表形状からの距離に応じて２つの変形の合成比率を変えることは、対応点からの距離に応じて合成比率を変えることと同様の効果を得ることができる。従って、画像上の各位置における表面形状からの距離さえ分かっていれば、第１の実施形態のような複数の対応点における放射基底関数を統合した（４）式よりも、表面形状からの距離に基づく単純な式に基づいて計算可能である。すなわち、より高速に合成変位場の計算が可能となる。

＜構成＞
以下、本実施形態に係る画像処理装置について、第１の実施形態との相違部分についてのみ説明する。

図８は、本実施形態に係る画像処理装置８０００の構成を示す図であり、新たに形状取得部８０１０が追加された点以外は、第１の実施形態における画像処理装置１０００と同様である。また、変形合成部１０５０に関しては、第１の実施形態と機能が異なるため以下でその機能を説明する。その他の構成については第１の実施形態と機能が同じであるため、同じ機能に関しては説明を省略する。

形状取得部８０１０は、画像から被検体の表面形状を抽出する。変形合成部１０５０は、第１の変位場と第２の変位場を被検体の表面形状からの距離に応じた比率で合成した合成変位場を生成する。

＜処理＞
図９は、本実施形態に係る画像処理装置８０００が行う全体の処理手順を示すフローチャートである。但し、図９において、ステップＳ９０１０からＳ９０４０、およびＳ９０６０からＳ９０７０の処理はそれぞれ、図２のフローチャートにおけるステップＳ２０１０からＳ２０４０、およびＳ２０６０からＳ２０７０の処理と同じであるため説明は省略する。以下、図２のフローチャートとの相違部分についてのみ説明する。

（Ｓ９０４５：被検体の表面形状の取得）
ステップＳ９０４５において、形状取得部８０１０は、画像から被検体の表面形状を抽出する。そして、抽出した被検体の表面形状のデータを変形合成部１０５０へと出力する。以下に具体的な処理を説明する。本ステップにおいて、形状取得部８０１０は、第２の画像から公知の形状抽出手法を用いて被検体の表面形状を抽出する。本実施形態では、例えば、動的輪郭法の一手法であるＳＮＡＫＥＳを用いて被検体の表面形状を抽出する。但し、表面形状の抽出手法はこれに限られるものではなく、画像の輝度の閾値で領域分割する手法や、Ｒｅｇｉｏｎ Ｇｒｏｕｉｎｇ法、Ｌｅｖｅｌ Ｓｅｔ法、Ｇｒａｐｈ Ｃｕｔ法など公知の何れの手法を用いてもよい。また、上述のような自動抽出手法ではなく、ユーザが不図示のＵＩを操作して手動で抽出してもよい。

（Ｓ９０５０：表面形状からの距離に基づく第１と第２の変形の合成）
ステップＳ９０５０において、変形合成部１０５０は、対応点取得部１０２０から取得した対応点からの距離と、形状取得部８０１０から取得した被検体の表面形状からの距離とに基づき、第１の変位場Φ_１ ^−１と第２の変位場Φ_２ ^−１を合成する。そして、生成した合成変位場Φ_１＋２ ^−１を画像生成部１０６０へと出力する。

以下、第１の実施形態のステップＳ２０５０と相違する点について、具体的な処理を説明する。本ステップでは、第２の画像上の各座標ｘにおける、第２の画像の被検体の表面形状から距離に基づいて合成比率を算出する。具体的には、以下の式によって、合成比率を算出する。

ここで、ｗ_１（ｘ）は、第１実施形態の（５）式の代わりに、第２の画像上の座標ｘにおいて第１の変位場を合成する合成比率を表している。ｗ_１（ｘ）は、第２の画像の被検体の表面形状からの距離に応じて決まる比率（０≦ｗ_１（ｘ）≦１）であり、（１３）式のように、第２の画像上の被検体の表面形状から座標ｘまでの距離をｄ_Ｓ２（ｘ）と定義したとき、距離ｄ_Ｓ２（ｘ）にほぼ反比例する。また、（１３）式において、βは定数であり、本実施形態では、例えばβ＝１．０を適用する。合成比率ｗ_１（ｘ）は、座標ｘが表面形状上に位置する、つまり距離ｄ_Ｓ２（ｘ）＝０．０のときに値１．０をとり、表面形状からの距離が無限遠に近付く、つまり距離ｄ_Ｓ２（ｘ）≒∞のときに、値０．０に限りなく近づく。すなわち、合成比率ｗ_１（ｘ）は、表面形状からの距離が遠ざかるにつれ、値が１．０から減少する。また、（１４）式より、ｗ_２（ｘ）は表面形状からの距離が遠ざかるにつれ、値が０．０から増加する。これにより、被検体の表面形状からの距離が大きくなるにつれ、対応点上における第１の変形の合成比率を減少させ、第２の変形の合成比率を増加させる効果がある。このように、合成比率は、第１の実施形態のような複数の対応点における放射基底関数を統合した（４）式及び（５）式よりも、（１３）式及び（１４）式による表面形状からの距離に基づく単純な式に基づいて計算可能である。以上によって、本実施形態における画像処理装置８０００の処理が実施される。

このような処理を行うことで、以下のような場合に効果を得ることができる。本実施形態における第２の画像が例えば光音響トモグラフィ画像の場合、撮像の際の光音響信号は、被検体の表面形状から体内の深部に向かうにつれ減衰する。そのため、画像化されたときに描出される解剖構造の輝度情報が表面形状から遠ざかるにつれ弱くなる傾向がある。ここで、図４の第２の画像３０２０が、本実施形態の第２の画像である光音響トモグラフィ画像を表しているものとする。図４の対応点Ｂ２は、表面形状の近くに位置することから、周辺の血管像が明瞭に描出されているが、対応点Ｂ３は、点Ｂ２よりも表面形状から離れているため、周辺の血管像がやや不明瞭に描出されていることが分かる。従って、これらの対応点をステップＳ９０２０で取得する際、自動抽出で取得する場合は血管を表す輝度が低下するため、抽出精度が低下する可能性があり、手動で抽出する場合でも、見えづらさから抽出精度が低下する可能性がある。このような場合に、本実施形態の処理を適用することで、対応点の抽出精度の信頼性が高い表面形状に近い領域は、対応点周辺において高精度に対応点を一致させる。一方、抽出精度の信頼性が低い表面形状から遠い領域は、対応点周辺において、所望の拘束条件に基づいて変形させる。これにより、適切な変形位置合わせを行うことができる。

さらに、光音響トモグラフィ画像で描出される血管構造は、体表に近いほど細かい血管（毛細血管）となり、特徴点である血管分岐部の数が多くなるため、対応点の数も多くなる傾向にある。さらに、上述のとおり、光音響信号は、体内の深部に向かうにつれ減衰するため、可視化される血管構造が少なくなり、対応点の数も少なくなる傾向にある。従って、光音響トモグラフィ画像では、対応点は被検体の表面形状の近傍に主に分布することになる。このことから、第１の実施形態における変形の合成方法と同様の効果を得ることができ、かつ第１の実施形態よりも高速に結果を取得できる。

以上説明したように、本実施形態によれば、被検体の表面形状からの距離が大きくなるにつれ、対応点を高精度に一致させる変形よりも所望の拘束条件に基づく変形に近付くようにすることで、対応点の信頼性に応じた適切な変形位置合わせを行うことができる。さらに、対応点が被検体の表面形状の近傍に主に分布する場合には、第１の実施形態の対応点からの距離に応じて２つの変形の合成比率を変えることと同様の効果を得ることができ、かつ第１の実施形態よりも高速に結果を取得できる。

（変形例４−１：第１の画像上の表面形状からの距離に基づいて変位場を合成してもよい）
第４の実施形態では、ステップＳ９０５０において、第２の画像上の被検体の表面形状からの距離に基づいて、第１の変位場と第２の変位場を合成していた。しかし、第１の画像上の被検体の表面形状からの距離に基づいて、２つの変位場を合成するようにしてもよい。この場合、ステップＳ９０４５において、形状取得部８０１０は、被検体の表面形状を第１の画像から抽出する。ここでは、第１の画像上で表面形状からの距離を算出するために、第１の実施形態の変形例１−１と同様、変位場Φ^−１を用いて第２の画像上の各座標に対応する第１の画像上の位置を参照することで、合成比率を算出する。すなわち、以下の式のように変位場Φ^−１を介して第２の画像上の座標に対応する第１の画像上の各参照位置における第１の画像上の表面形状からの距離に基づいた比率を算出できる。

ここで、ｗ_１（ｘ）は、第４実施形態の（１３）式の代わりに、第２の画像上の座標ｘにおいて第１の変位場を合成する合成比率を表している。以下に、第１の実施形態の変形例１−１と相違する点のみ説明する。（１５）式において、Φ^−１（ｘ）は、変位場Φ^−１を用いて第２の画像上の座標ｘを第１の画像に向けて変位させた座標を表す。この変位場の算出方法は、第１の実施形態の変形例１−１と同様であるため、説明を省略する。また、ｗ_１（ｘ）は、第１の画像の被検体の表面形状からの距離に応じて決まる比率（０≦ｗ_１（ｘ）≦１）である。（１５）式のように、第１の画像上の被検体の表面形状から第１の画像上の座標Φ^−１（ｘ）までの距離をｄ_Ｓ１（Φ^−１（ｘ））と定義したとき、距離ｄ_Ｓ１（Φ^−１（ｘ））に反比例する。具体的な説明については、（１３）式において、ｄ_Ｓ２（ｘ）をｄ_Ｓ１（Φ^−１（ｘ））に置き換え、説明における第２の画像を第１の画像に読み替えればよいため、省略する。

さらに、第１の画像上の表面形状からの距離と、第２の画像上の表面形状からの距離との両方に基づいて、変位場を合成してもよい。この場合、ステップＳ９０４５において、形状取得部８０１０は、第１の画像と第２の画像の両方から夫々被検体の表面形状を抽出する。このとき、合成比率は以下の式で算出できる。

ここで、ｗ_１（ｘ）は、第４の実施形態の（１３）式や上述の（１５）式の代わりに、第２の画像上の座標ｘにおいて第１の変位場を合成する合成比率を表している。具体的な説明は、第１の実施形態の変形例１−１と同様であるため、省略する。

以上説明したように、本変形例によれば、第１の画像上の被検体の表面形状からの距離に基づいて変位場を合成することで、第１の画像上の表面形状を考慮した合成変位場を取得できる。

（その他の実施形態）
本発明は、上述の実施形態の１以上の機能を実現するプログラムを、ネットワーク又は記憶媒体を介してシステム又は装置に供給し、そのシステム又は装置のコンピュータにおける１つ以上のプロセッサがプログラムを読出し実行する処理でも実現可能である。また、１以上の機能を実現する回路（例えば、ＡＳＩＣ）によっても実現可能である。

１１０：データサーバ、１２０：表示部、１０００：画像処理装置、１０１０：データ取得部、１０２０：対応点取得部、１０３０：第１の変形取得部、１０４０：第２の変形取得部、１０５０：変形合成部、１０６０：画像生成部、１０７０：表示制御部、８０００：画像処理装置、８０１０：形状取得部

Claims (18)

1. 第１の画像と第２の画像との間で被検体の特徴点を対応付けた対応点を取得する対応点取得手段と、
   前記対応点に基づく第１の変形を取得する第１の変形取得手段と、
   前記第１の画像と前記第２の画像との間で画像の輝度情報を一致させることを含む第２の変形を取得する第２の変形取得手段と、
   前記第１の画像と前記第２の画像との少なくとも一方の画像上の各位置について、前記対応点から該位置への距離に基づいて、前記第１の変形と前記第２の変形とを合成する合成手段と、
   を備えることを特徴とする画像処理装置。
2. 前記合成手段は、前記対応点に近い位置ほど前記第１の変形の変位場の合成比率が高くなり、前記対応点から遠い位置ほど前記第２の変形の変位場の合成比率が高くなるように、前記第１の変形と前記第２の変形とを合成することを特徴とする請求項１に記載の画像処理装置。
3. 前記合成手段は、対応点ごとに、対応関係に関する信頼度を設定し、各対応点の信頼度に応じた合成比率に基づいて、前記第１の変形と前記第２の変形とを合成することを特徴とする請求項１又は２に記載の画像処理装置。
4. 前記合成手段は、前記画像上の各位置について、該位置の近傍に存在する対応点を前記対応点の中から選択し、該選択した対応点から該位置への距離に基づいて算出した合成比率を用いて、前記第１の変形と前記第２の変形とを合成することを特徴とする請求項１乃至３の何れか１項に記載の画像処理装置。
5. 前記対応点取得手段は、複数の対応点を取得し、
   前記合成手段は、前記画像上の各位置について、前記複数の対応点の夫々に対応する、該対応点からの距離に基づく放射基底関数を統合することで得られる合成比率関数を用いて、前記第１の変形と前記第２の変形とを合成することを特徴とする請求項１乃至３の何れか１項に記載の画像処理装置。
6. 前記放射基底関数は、前記対応点上で最小値０をとり、最大値が１であることを特徴とする請求項５に記載の画像処理装置。
7. 前記合成比率関数は、前記複数の対応点の夫々に対応する前記放射基底関数を乗算によって統合した関数であることを特徴とする請求項５又は６に記載の画像処理装置。
8. 前記合成比率関数は、前記複数の対応点の夫々に対応する前記放射基底関数を加算によって統合した関数であることを特徴とする請求項５又は６に記載の画像処理装置。
9. 前記第２の変形における対応点間の位置のずれを表す対応点残差を取得する対応点残差取得手段を更に備え、
   前記合成手段は、前記対応点残差に基づいて前記放射基底関数の空間的な広がりの度合いを決定することを特徴とする請求項５に記載の画像処理装置。
10. 前記第１の画像を、前記合成手段により合成された変形に基づいて変形させた変形画像を生成する画像生成手段を更に備えることを特徴とする請求項１乃至９の何れか１項に記載の画像処理装置。
11. 前記画像生成手段は、前記第２の画像と前記変形画像との差分画像を更に生成することを特徴とする請求項１０に記載の画像処理装置。
12. 前記第１の変形は、前記第１の画像と前記第２の画像との間で対応点を一致させる変形であることを特徴とする請求項１乃至１１の何れか１項に記載の画像処理装置。
13. 第１の画像と第２の画像の間で被検体の特徴点を対応付けた対応点を取得する対応点取得手段と、
    前記対応点に基づく第１の変形を取得する第１の変形取得手段と、
    前記第１の画像と前記第２の画像との間で画像の輝度情報を一致させることを含む第２の変形を取得する第２の変形取得手段と、
    前記第１の画像と前記第２の画像との少なくとも一方の画像における前記被検体の表面形状を取得する形状取得手段と、
    前記第１の画像と前記第２の画像との少なくとも一方の画像上の各位置について、前記表面形状から該位置への距離に基づいて、前記第１の変形と前記第２の変形とを合成する合成手段と、
    を備えることを特徴とする画像処理装置。
14. 第１の画像と第２の画像との間で被検体の特徴点を対応付けた対応点を取得する対応点取得手段と、
    前記対応点に基づく第１の変形を取得する第１の変形取得手段と、
    前記第１の画像と前記第２の画像の輝度情報に関する画像類似度項と、変形が所定の物理モデルから逸脱することに対してペナルティを与える正則化項との少なくとも一つを含む評価関数を用いて第２の変形を取得する第２の変形取得手段と、
    前記第１の画像と前記第２の画像との少なくとも一方の画像上の各位置について、前記対応点から該位置への距離に基づいて、前記第１の変形と前記第２の変形とを合成する合成手段と、
    を備えることを特徴とする画像処理装置。
15. 画像処理装置の制御方法であって、
    第１の画像と第２の画像との間で被検体の特徴点を対応付けた対応点を取得する対応点取得工程と、
    前記対応点に基づく第１の変形を取得する第１の変形取得工程と、
    前記第１の画像と前記第２の画像との間で画像の輝度情報を一致させることを含む第２の変形を取得する第２の変形取得工程と、
    前記第１の画像と前記第２の画像との少なくとも一方の画像上の各位置について、前記対応点から該位置への距離に基づいて、前記第１の変形と前記第２の変形とを合成する合成工程と、
    を有することを特徴とする画像処理装置の制御方法。
16. 画像処理装置の制御方法であって、
    第１の画像と第２の画像の間で被検体の特徴点を対応付けた対応点を取得する対応点取得工程と、
    前記対応点に基づく第１の変形を取得する第１の変形取得工程と、
    前記第１の画像と前記第２の画像との間で画像の輝度情報を一致させることを含む第２の変形を取得する第２の変形取得工程と、
    前記第１の画像と前記第２の画像との少なくとも一方の画像における前記被検体の表面形状を取得する形状取得工程と、
    前記第１の画像と前記第２の画像との少なくとも一方の画像上の各位置について、前記表面形状から該位置への距離に基づいて、前記第１の変形と前記第２の変形とを合成する合成工程と、
    を有することを特徴とする画像処理装置の制御方法。
17. 画像処理装置の制御方法であって、
    第１の画像と第２の画像との間で被検体の特徴点を対応付けた対応点を取得する対応点取得工程と、
    前記対応点に基づく第１の変形を取得する第１の変形取得工程と、
    前記第１の画像と前記第２の画像の輝度情報に関する画像類似度項と、変形が所定の物理モデルから逸脱することに対してペナルティを与える正則化項との少なくとも一つを含む評価関数を用いて第２の変形を取得する第２の変形取得工程と、
    前記第１の画像と前記第２の画像との少なくとも一方の画像上の各位置について、前記対応点から該位置への距離に基づいて、前記第１の変形と前記第２の変形とを合成する合成工程と、
    を有することを特徴とする画像処理装置の制御方法。
18. コンピュータを、請求項１乃至１４の何れか１項に記載の画像処理装置の各手段として機能させるためのプログラム。

Images