JP2005169077A - 第１データセットのデータポイントを第２データセットのデータポイントと関係付ける変換を見つける方法、およびその方法を実行するコンピュータプログラム - Google Patents

Abstract

【課題】異なる時点で取得された画像あるいは異なる画像を重ね合わせて表示することが望まれる場合が多い。その実現のために複数のデータ間を関連付ける変換を見つける方法を提供する。
【解決手段】異なる画像等のデータセットを、その２つの画像を関係付ける変換を推定することによって照合する方法。２つのデータセット自体から導出した確率密度関数に基づいて計算される類似性測度を用いることによって、様々な候補変換をスコアリングする。確率密度関数は局所尤度密度推定（ＬＬＤＥ）に基づく。本技法は、画像の位置合わせに、またデータの融合に適用することができる。
【選択図】 図１

Description

［発明の詳細な説明］
本発明は、データセット(data set)の互いに対する一致度を測定するために用いられる類似性測度に関する。

あるデータセットと別のデータセットの類似性または一致度を測定できることが望ましい分野は多い。例えば、２つの画像を比較して、この２つの画像を関係付ける空間変換を求める状況は多い。例えば、画像の時間シーケンスを圧縮する一方法は、第１フレームに見られる画像特徴が第２フレームのどこに配置されているかを記述する移動ベクトルを検出することである。別法として、医用画像分野では、同じ物体（例えば患者の身体の一部）の２つの画像を重ね合わせて表示することが望ましい場合が多い。この２つの画像は、異なる時点で撮影された画像、または異なるモダリティで撮影された画像である可能性がある。重ね合わせた２つの画像の表示は、それらの画像が位置合わせ(register)されていることを要し、この位置合わせには、それぞれの画像において対応する画像特徴の位置を見つけることが伴う。同様の技法は、例えば可視または赤外線画像を合成開口レーダ等の異なる種類の画像と位置合わせする、いわゆるデータの融合を伴う他の分野でも有用である。この必要性は、位置合わせすることが有用である可能性があるＥＣＧおよび血中酸素濃度等の１次元信号、および３Ｄ画像が時間の経過とともに変化するコントラスト増強ＭＲＩ等の４次元画像においても生じる可能性がある。

このような方法でデータセットを位置合わせする際には通常、様々な候補変換を画像に適用して、２つの画像をどの程度正確に位置合わせするかに基づいてスコアリングする。スコアは、類似性測度を用いて自動的に計算することができ、多数の様々な類似性測度が提案されており（例えば最小二乗誤差（differences）、相互相関等）、異なるタイプのデータセットには概して異なる類似性測度が適している。

現在まで、類似性測度は傾向として、照合する特定のデータセットについて試行錯誤法で選択されるか、あるいは、データ生成の基礎を成す機構、例えば画像生成の背後にある物理プロセスの理解から導出されてきた。しかし、所与の類似性測度が特定のタイプのデータに適していることを予測することは困難である可能性があり、新たな類似性測度の設計は時間のかかる困難なプロセスである。

本発明によれば、所与の変換が正しいものである尤度を測定する類似性測度の設計は、統計的な枠組みにおいて２つのデータセット自体の特性に基づいて、最尤推定法を用いて行われる。

より詳細には、本発明は、第１データセットのデータポイント(data point)を第２データセットのデータポイントと関係付ける変換を見つける方法を提供する。本方法は、複数の候補変換のそれぞれについて、候補変換が２つのデータセット間の関係を正確に表す尤度を表す尤度測度（likelihood measure）を計算することを含み、この尤度測度は、候補変換において対応すると考えられる２つのデータセットのデータポイント対について計算される局所パラメータ結合(local parametric joint)および周辺(marginal)確率密度関数（ｐｄｆ）に基づき、関数のパラメータは、データポイント対の局所尤度基準を最大化することによって設定される。

好ましくは、確率密度関数は、データポイント対に対して局所的なデータポイントの特徴に基づく。したがって、ポイント毎に異なる可能性があるパラメータ関数は、データポイントの近傍に有効であると仮定されるが、互いから離れたポイントに関しては異なる可能性がある。

したがって、本発明では、確率密度関数は、検討中の特定の候補変換によって生成されるデータポイント対について計算される。したがって尤度測度は、照合中の２つのデータセットに適合する。データセットの生成機構に関する事前知識は必要ない。

確率密度関数は、二次多項式関数等の単純なパラメータ関数であってよい。関数のパラメータは、データポイント対の局所尤度基準を最大化することによって設定され、この局所尤度基準は、確率密度関数と、好ましくは窓関数およびスケールパラメータとに基づく。窓関数およびスケールパラメータは、本方法の性能を調整するために発見的に(heuristically)調節することができる。パラメータは、それぞれのデータポイント対について、または少なくとも画像の異なる部分について再計算してもよい。

本方法は、データセットが例えば画素の輝度(intensity)である画像の位置合わせに適用することができ、この場合、変換は、２つの画像を位置合わせする空間変換である。画像は例えば、異なるモダリティおよび／または時間間隔をおいた画像である可能性があり、本方法は特に、医用画像の位置合わせに適している。

しかし、本方法は、信号または測定値を表すもの等、他のデータセットにも適用することができる。

本発明は、コンピュータソフトウェアにおいて実施することができ、よって適切にプログラムされたコンピュータシステム上で本方法を実行するコンピュータプログラムを提供することができることが理解されるだろう。別法として、当然ながら、本発明は、ファームウェアまたは専用ハードウェアにおいて実施してもよい。

本発明を、例として、添付図面を参照してさらに説明する。

本発明の１実施形態を、画像である２つのデータセットを参照して説明する。たとえば、図１に示すように、２つの画像ＩおよびＪが与えられ、ＪがＩを変換Ｔによって変換したものであり、ＴがＩおよびＪそれぞれの座標系を関係付けることになっている場合、目的は、変換Ｔがどの程度「正しい」変換に近いかを正確にスコアリングする類似性測度を開発することである。言い換えれば、Ｊが変換Ｔによって変換された画像Ｉである尤度はどれほどか、である。変換Ｔは、Ｉの画素とＪの画素の間の空間的対応のセットを定める。

Ｉが元画像であり、Ｊがターゲット画像であると仮定すると、これは、画像Ｉの座標系において座標ｘを有する画素が画像Ｊの座標系において座標Ｔ（ｘ）を有する画素と一致することを意味する。｛ｘ_ｋ｝（ｋ＝１，．．．，ｎ）は、画像Ｉの画素セットを示すものとする。ｘ_ｋに対応する輝度値をｉ_ｋで示す。同様に、Ｊにおける対応画素の輝度値、すなわち座標Ｔ（ｘ_ｋ）を有する画素をｊ_ｋで示す。

所与の変換Ｔからこのように進んでいくと、輝度対（intensity pair）Ａ（Ｔ）＝｛（ｉ_ｋ，ｊ_ｋ），ｋ＝１，．．．，ｎ｝のセットが得られることに留意されたい。また、輝度対のセットは、推定変換Ｔの変化に伴って変化することが予期される、すなわち輝度対のセットは数学的にＴの関数として考えることができることにも留意されたい。Ｔの尤度、すなわちＴが正しい位置合わせ変換である確率は、次式によって得られる。

ここで、ｐは結合輝度ｐｄｆであり、ｐ_１およびｐ_２はそれぞれＩおよびＪに対応する周辺ｐｄｆである。しかし、式（１）は実際には、これら３つのｐｄｆが指定された場合にのみ利用することができる。

本実施形態において用いられる基本原理は、局所尤度密度推定(Local Likelihood Density Estimation)（ＬＬＤＥ）を用いて、一般にそうであるように、ｐｄｆ全体の単純なパラメータ形式がない場合にも、ｐｄｆが局所的に単純なパラメータ形式によって近似できると仮定することによって、データセットからｐｄｆを推定することである。本実施形態ではさらに、ｐｄｆの対数は局所的に二次多項式であると仮定するが、任意の他のパラメータ形式、例えばより高次の多項式を用いることもできることを理解することができる。

したがって、任意の観察された輝度対（ｉ，ｊ）について、ｉおよびｊにそれぞれ十分近い

および

に関してｐｄｆの多項式近似が（ｉ，ｊ）の近傍において有効であるように、係数セット（ａ_０，ａ_１，．．．，ａ_５）を計算する。

本方法の１実施形態において、多項式の係数は、次の局所尤度基準を最大化するように設定される。

ここで、

は所与の窓関数であり、ｓ_１ｓ_２はスケールパラメータである。

の通常の選択肢は、１）箱形関数、すなわち

ならびに２）ガウス関数、すなわち

である。

が局所重み関数の形状を決める一方、スケールパラメータｓ_１およびｓ_２は、窓のサポート（window support）によって定められる近傍の拡張を調節し、よって局所サイズを制御する。ｓ１およびｓ２はともに本方法のパラメータであり、これらのパラメータを調整する方法は、本発明とは完全に独立している。実験に基づく実施では、ｓ１およびｓ２は、検討中の特定の画像タイプに応じて、同様の画像を用いた以前の実験において良好な位置合わせ結果をもたらすことが立証された値に設定することができる。例として、実際に有用であることが分かった規則の例は、ｓ１を画像Ｉ中の最大輝度の所与の割合に設定することであり、この割合は、ＣＴ画像の場合に１％、ＭＲ画像の場合に２％、ＰＥＴまたはＳＰＥＣＴ画像の場合に５％である。同一の規則を、画像Ｊのタイプに応じてｓ２に当てはめてもよい。

理論的な観点で見れば、ｓ１およびｓ２が大きくなるほど、局所多項式モデルを当てはめるために用いられる輝度近傍は大きくなり、よって結果として得られるｐｄｆの推定値は滑らかになる。一方、ｓ１およびｓ２に小さな値を選択すると、ｐｄｆ推定値により急激な変化が許容されるため、より大きな柔軟性が得られるが、その代償として、位置合わせは画像のアーチファクトおよび雑音により敏感になる傾向がある。このよく知られたトレードオフのスケールパラメータは、統計的推定では一般的なものであり、「帯域幅選択」として知られる。

周辺ｐｄｆ ｐ_１およびｐ_２の推定は同じ原理から得られる。局所的な二次多項式の形式を次のように仮定し、

ｐ_２についても同様に行う。

局所多項式は、式（３）を１次元にしたものを最大化することによって当てはめ、

ｐ_２についても同様に行う。結合ｐｄｆの推定との一貫性を保つために、窓関数

および

は、二次元窓関数

から、

として定義する。

ＬＬＤＥの主な利点は、ｐｄｆの形状に関する非常にわずかな仮定に頼る一方で正確な密度推定値を与えることができることにある。これは、類似性測度（１）に伴う密度ｐ、ｐ_１およびｐ_２の指定誤りが実質的な位置合わせ誤差に反映される場合があるため、重要である。

パラメータ手法とは異なり、選択に画像取得プロセスをモデル化する退屈で困難なタスクを必要とする結合ｐｄｆのグローバルモデルは不要である。さらに、パラメータ手法は通常、その実際の実施を困難にする非凸形の（non-convex）最適化問題をもたらす。

図２は、ステップを例示的なアルゴリズムで行われることのできる順番で示すフロー図である。ステップ２００において２つのデータセット、例えば２つの画像を取得し、次にステップ２０１では、上述のようにスケールパラメータｓ_１およびｓ_２ならびに窓関数

を選択する必要がある。

ステップ２０２において、式４および式５を用いて、それぞれのデータセットの周辺ｐｄｆを計算する、つまり、式５を最大化することにより式４の局所多項式を当てはめることによってａ_０ ａ_１ ａ_２を計算する。（このステップは、ａ_０ ａ_１ ａ_２が必要となる尤度の計算以前の任意の時点で実行することができる。）

次にステップ２０３において、所与の候補変換Ｔについて、輝度対Ａ（Ｔ）を見つける。図１を参照すると、図示される変換Ｔは、右側へ１画素分の単純なシフトである。したがって、変換Ｔにおいて、画素ｊ２は画素ｉ１に対応し、画素ｊ３は画素ｉ２に対応する、等となる。このように、輝度対Ａ（Ｔ）のセットが２つの画像フレームの下に示される。

次にステップ２０４において輝度対を取得し、ステップ２０５において、式（３）の局所尤度基準を最大化することによって多項式係数ａ_１、ａ_２、ａ_３、ａ_４、ａ_５を計算する。（係数ａ_１、ａ_２、ａ_３、ａ_４、ａ_５は対毎に異なることに留意されたい。）

多項式係数を取得したら、ステップ２０６において式（２）を用いて結合確率密度関数ｐを計算することができる。

ステップ２０７では、全ての輝度対を処理したかどうかを調べ、全ての輝度対を処理していない場合、全ての対が完了するまでプロセスを繰り返す。全ての対が完了したら、ステップ２０９において、式（１）を用いてその変換の尤度測度を計算することができる。

次にステップ２１０および２１１において、他の候補変換を取得し、プロセスを繰り返す。

すると、各候補変換に尤度が関連付けられ、ステップ２１２において、尤度の最も高い変換を選択することができる。これは、第２画像フレームが第１画像フレームに関係する方法の最良推定値を表す。ステップ２１３に示すように、ｓ_１、ｓ_２および

に異なる値を用いてプロセスを繰り返し、結果を改善することが可能である。

図３（ａ）は、脳のＣＴ−ＭＲ画像対の位置合わせ後の結合ヒストグラムを示す。このように、ヒストグラム中の各ポイントは、位置合わせされた画像において所与のＣＴ輝度および所与のＭＲ輝度を有する画素数を表す。良好な位置合わせは結果として、明確な構造および少数のピークをヒストグラム中に生じる。不良な位置合わせは結果として、より多くの散乱したピークを生じる。図３（ａ）の未処理のヒストグラムは、特にＣＴ骨−ＭＲ軟組織領域で雑音が多い。この未処理のヒストグラム（すなわちｐ（ｉ，ｊ）＝ｈ（ｉ，ｊ））のみからｐ（ｉ，ｊ）を直接推定すると、変換の小さな変化がヒストグラムひいては類似性スコアに大きな変化を生じるため、結果として位置合わせ変換の推定の頑強性が限定される。図３（ｂ）は、単純で均一なぼかし技法、または単純なパルゼン(Parzen)窓掛け処理を用いてヒストグラムを平滑化したものを示す。ＣＴ軟組織−ＭＲ軟組織領域は平滑化され、よって変換に対する類似性測度の感受性を有効に低減することが分かる。しかし、分布の鮮明さは必然的に低くなるため、結合ヒストグラムのそのような単純すぎる平滑化は、分布の重要な特性を失い、正確なｐ（ｉ，ｊ）の推定を偏らせる可能性がある。この結果、２つの画像の位置合わせ精度は低くなる（変換に対する類似性スコアの感受性を低くしたことによる必然的な結果である）。図３（ｃ）および図３（ｄ）は、本発明に従ってｓ_１およびｓ_２に異なる値を用いて推定されたｐ（ｉ，ｊ）を示す。図示のように、雑音は相殺されるが、分布の特性、例えば高いピークの鮮明さは保持される。雑音の低減により頑強性が増す一方で、鮮明なピークの保持により正確な位置合わせを行うことができる。

図４（ａ）および図４（ｂ）はそれぞれＭＲ画像およびＣＴ画像のスライスの例を示し、図４（ｃ）は本発明の技法を用いてＭＲ画像とＣＴ画像を融合した結果を示す。頭蓋構造がＣＴ（明るい楕円形）とＭＲ（脳と頭皮の間）の両方の画像で良好に整列していることが分かる。

本発明に従って処理することができる２つのデータセットを概略的に示す。 本発明の一例のフロー図である。 （ａ）（ｂ）は脳のＣＴ−ＭＲ画像対の平滑化された結合ヒストグラムを示す。（ｃ）（ｄ）は本発明の技法を用いた対応する推定ｐｄｆである。 （ａ）（ｂ）（ｃ）はそれぞれ、ＭＲ画像のスライスの例、ＣＴ画像のスライスの例、本発明の技法を用いてＭＲ画像とＣＴ画像を融合した結果を示す。

Claims (11)

1. 第１データセットのデータポイントを第２データセットのデータポイントと関係付ける変換を見つける方法であって、
   複数の候補変換のそれぞれについて、前記候補変換が２つの前記データセット間の関係を正確に表す尤度を表す尤度測度を計算することを含み、
   前記尤度測度は、前記候補変換において対応すると考えられる２つの前記データセットのデータポイント対について計算される局所パラメータ結合および周辺確率密度関数に基づき、
   前記関数のパラメータは、前記データポイント対の局所尤度基準を最大化することによって設定される、
   第１データセットのデータポイントを第２データセットのデータポイントと関係付ける変換を見つける方法。
2. 各確率密度関数の前記局所尤度基準は、
   前記確率密度関数と、
   前記局所データポイントの特徴の窓関数およびスケールパラメータのうちの少なくとも１つと
   に基づいて計算される請求項１に記載の方法。
3. 前記窓関数および前記スケールパラメータは発見的に調節される請求項２に記載の方法。
4. 前記確率密度関数は多項式関数である請求項１〜３のいずれか一項に記載の方法。
5. 前記確率密度関数は二次多項式関数である請求項４に記載の方法。
6. 前記データセットは信号または測定値を表す請求項１〜５のいずれか一項に記載の方法。
7. 前記データセットはそれぞれ画像を表し、
   前記変換は前記画像を位置合わせする空間変換である請求項１〜６のいずれか一項に記載の方法。
8. 前記データセットは異なるモダリティ画像を表す請求項７に記載の方法。
9. 前記データセットは医用画像を表す請求項７または８に記載の方法。
10. 前記データセットは時間的に分離される請求項１〜９のいずれか一項に記載の方法。
11. プログラムされたコンピュータシステム上で請求項１〜１０のいずれか一項に記載の方法を実行するプログラムコード手段を含むコンピュータプログラム。
    

Images