(57)【要約】【課題】視覚的に相似していない局所領域間の局所位置合わせを良好な精度と効率性とをもって行うことができるようにする。【解決手段】画像対の視覚的に相似していない局所領域をリンクする方法において、対象物の少なくとも第１及び第２のデジタル画像を取得するステップと、３Ｄランドマーク検出法を使用して第１及び第２のデジタル画像の一部との間の対応を決定するステップとを有し、前記第１の画像の一部はユーザ定義されたクリックポイントとしてユーザにより選択された関心点を含んでおり、前記第１の画像の一部は、その位置合わせの精度が第１の画像の他の部分における位置合わせの精度よりも高くなくてはならないような第１の画像内の箇所であり、前記対応は第１の画像内のユーザ定義されたクリックポイントと第２の画像内の所定の点との間の幾何学的関係を用いて決定される。【選択図】なし【特許請求の範囲】
【請求項１】
画像対の視覚的に相似していない局所領域をリンクする方法において、対象物の少なくとも第１及び第２のデジタル画像を取得するステップと、３Ｄランドマーク検出法を使用して第１のデジタル画像の一部と第２のデジタル画像の一部との間の対応を決定するステップとを有しており、前記第１の画像の一部はユーザ定義されたクリックポイントとしてユーザにより選択された関心点を含んでおり、前記第１の画像の一部は、その位置合わせの精度が第１の画像の他の部分における位置合わせの精度よりも高くなくてはならないような第１の画像内の箇所であり、前記対応は第１の画像内のユーザ定義されたクリックポイントと第２の画像内の所定の点との間の幾何学的関係を用いて決定される、ことを特徴とする方法。
【請求項２】
前記３Ｄランドマーク検出法は学習をベースにした検出法であり、訓練データとして全身ボリュームをデータベースに収集するステップと、各訓練ボリュームについて関心ランドマークを手作業でラベル付けするステップと、所与の訓練ボリュームの関心ランドマークに他のランドマークよりも大きな重みを付与する重み付けアラインメントアルゴリズムを用いて訓練ボリュームをアラインするステップと、所与の訓練ボリュームの関心ランドマークの周囲のデータをクロッピングして前記関心ランドマークを含んだ正事例を形成するステップと、前記訓練ボリュームの他の箇所のデータをランダムにクロッピングして負事例を形成するステップと、前記正及び負の訓練事例から成るデータベースを用いて判別分類子を開発するステップと、ランドマークを検出するために対象物に対して新たな全身ボリューム走査を行うステップとを有しており、前記分類子により正として分類された窓が可能性のある検出結果として返され、検出された窓はその応答に基づいてランク付けされ、その中で最も応答の大きな窓は新たな全身ボリュームにおいて検出されたランドマークである尤度が最も高い、請求項１記載の方法。
【請求項３】
前記３Ｄランドマーク検出法は、第１のデジタル画像についてスケール不変の顕著特徴の集合を計算するステップと、第２のデジタル画像についてスケール不変の顕著特徴の集合を計算するステップと、要素分布の平均と幅とを有するガウシアン混合分布から成る幾何構成コンテクスト（ＧＣＣ）を第１のデジタル画像の顕著特徴と第２のデジタル画像の顕著特徴との間の対応に応じてモデル化するステップとから成る顕著領域検出法を含んでおり、前記ＧＣＣは第１の画像の中のユーザ定義されたクリックポイントに対応する第２のデジタル画像内の点の空間尤度を表すものである、請求項１記載の方法。
【請求項４】
純粋な並進、スケーリング及び並進、ならびに相似変換から成る３つの変換クラスのうちの少なくとも１つについて、前記ＧＣＣモデルを導出する、請求項３記載の方法。
【請求項５】
第１及び第２のデジタル画像の少なくとも一部の間の前記対応は、前記第１のデジタル画像内の第１の点と前記第２のデジタル画像内の第２の点との間の点ごとの対応を確立するステップを含む、請求項１記載の方法。
【請求項６】
前記第１のデジタル画像内の第１の点と前記第２のデジタル画像内の第２の点との間の点ごとの対応を確立するステップを画像内の異なる複数の点について順次行う、請求項１記載の方法。
【請求項７】
前記第１及び第２のデジタル画像を（ａ）コンピュータ断層撮影（ＣＴ）ボリューム、（ｂ）磁気共鳴画像法（ＭＲＩ）、又は（ｃ）陽電子放射断層撮影法（ＰＥＴ）のうちの１つを用いて取得する、請求項１記載の方法。
【請求項８】
一連のステップを実行する機械により実行可能な一連の命令を有体的に実現した機械可読記憶装置において、前記ステップは、対象物の少なくとも第１及び第２のデジタル画像を取得するステップと、３Ｄランドマーク検出法を使用して第１のデジタル画像の一部と第２のデジタル画像の一部との間の対応を決定するステップとを有しており、前記第１の画像の一部はユーザ定義されたクリックポイントとしてユーザにより選択された関心点を含んでおり、前記第１の画像の一部は、その位置合わせの精度が第１の画像の他の部分における位置合わせの精度よりも高くなくてはならないような第１の画像内の箇所であり、前記対応は第１の画像内のユーザ定義されたクリックポイントと第２の画像内の所定の点との間の幾何学的関係を用いて決定される、ことを特徴とする機械可読記憶装置。
【請求項９】
前記３Ｄランドマーク検出法は学習をベースにした検出法であり、訓練データとして全身ボリュームをデータベースに収集するステップと、各訓練ボリュームについて関心ランドマークを手作業でラベル付けするステップと、所与の訓練ボリュームの関心ランドマークに他のランドマークよりも大きな重みを付与する重み付けアラインメントアルゴリズムを用いて訓練ボリュームをアラインするステップと、所与の訓練ボリュームの関心ランドマークの周囲のデータをクロッピングして前記関心ランドマークを含んだ正事例を形成するステップと、前記訓練ボリュームの他の箇所のデータをランダムにクロッピングして負事例を形成するステップと、前記正及び負の訓練事例から成るデータベースを用いて判別分類子を開発するステップと、ランドマークを検出するために対象物に対して新たな全身ボリューム走査を行うステップとを有しており、前記分類子により正として分類された窓が可能性のある検出結果として返され、検出された窓はその応答に基づいてランク付けされ、その中で最も応答の大きな窓は新たな全身ボリュームにおいて検出されたランドマークである尤度が最も高い、請求項８記載の機械可読記憶装置。
【請求項１０】
前記３Ｄランドマーク検出法は、第１のデジタル画像についてスケール不変の顕著特徴の集合を計算するステップと、第２のデジタル画像についてスケール不変の顕著特徴の集合を計算するステップと、要素分布の平均と幅とを有するガウシアン混合分布から成る幾何構成コンテクスト（ＧＣＣ）を第１のデジタル画像の顕著特徴と第２のデジタル画像の顕著特徴との間の対応に応じてモデル化するステップとから成る顕著領域検出法を含んでおり、前記ＧＣＣは第１の画像の中のユーザ定義されたクリックポイントに対応する第２のデジタル画像内の点の空間尤度を表すものである、請求項８記載の機械可読記憶装置。
【請求項１１】
前記ＧＣＣモデルは、純粋な並進、スケーリング及び並進、ならびに相似変換から成る３つの変換クラスのうちの少なくとも１つについて導出される、請求項１０記載の機械可読記憶装置。
【請求項１２】
第１及び第２のデジタル画像の少なくとも一部の間の前記対応は、前記第１のデジタル画像内の第１の点と前記第２のデジタル画像内の第２の点との間の点ごとの対応の確立を含む、請求項１０記載の機械可読記憶装置。
【請求項１３】
前記第１のデジタル画像内の第１の点と前記第２のデジタル画像内の第２の点との間の点ごとの対応の確立は画像内の異なる複数の点について順次行われる、請求項１２記載の機械可読記憶装置。
【請求項１４】
前記第１及び第２のデジタル画像は、（ａ）コンピュータ断層撮影（ＣＴ）ボリューム、（ｂ）磁気共鳴画像法（ＭＲＩ）、又は（ｃ）陽電子放射断層撮影法（ＰＥＴ）のうちの１つを用いて取得される、請求項１０記載の機械可読記憶装置。
【請求項１５】
医用画像データ対の間の空間的対応に高い自由度を実現するために医用画像データ対を位置合わせするシステムにおいて、該システムは、対象物の少なくとも第１及び第２のデジタル画像を取得する手段と、第１及び第２の医用画像を表すデータを記憶する手段と、３Ｄランドマーク検出法を使用して第１のデジタル画像の一部と第２のデジタル画像の一部との間の対応を決定する手段とを有しており、前記第１の画像の一部はユーザ定義されたクリックポイントとしてユーザにより選択された関心点を含んでおり、前記第１の画像の一部は、その位置合わせの精度が第１の画像の他の部分における位置合わせの精度よりも高くなくてはならないような第１の画像内の箇所であり、前記対応は第１の画像内のユーザ定義されたクリックポイントと第２の画像内の所定の点との間の幾何学的関係を用いて決定される、ことを特徴とする医用画像データ対を位置あわせするシステム。
【請求項１６】
３Ｄランドマーク検出法を実行する手段は学習をベースにした検出法を使用するものであり、該検出法は、訓練データとして全身ボリュームをデータベースに収集するステップと、各訓練ボリュームについて関心ランドマークを手作業でラベル付けするステップと、所与の訓練ボリュームの関心ランドマークに他のランドマークよりも大きな重みを付与する重み付けアラインメントアルゴリズムを用いて訓練ボリュームをアラインするステップと、所与の訓練ボリュームの関心ランドマークの周囲のデータをクロッピングして前記関心ランドマークを含んだ正事例を形成するステップと、前記訓練ボリュームの他の箇所のデータをランダムにクロッピングして負事例を形成するステップと、前記正及び負の訓練事例から成るデータベースを用いて判別分類子を開発するステップと、ランドマークを検出するために対象物に対して新たな全身ボリューム走査を行うステップとを有しており、前記分類子により正として分類された窓が可能性のある検出結果として返され、検出された窓はその応答に基づいてランク付けされ、その中で最も応答の大きな窓は新たな全身ボリュームにおいて検出されたランドマークである尤度が最も高い、請求項１５記載のシステム。
【請求項１７】
３Ｄランドマーク検出法を使用する前記手段は、第１のデジタル画像についてスケール不変の顕著特徴の集合を計算するステップと、第２のデジタル画像についてスケール不変の顕著特徴の集合を計算するステップと、要素分布の平均と幅とを有するガウシアン混合分布から成る幾何構成コンテクスト（ＧＣＣ）を第１のデジタル画像の顕著特徴と第２のデジタル画像の顕著特徴との間の対応に応じてモデル化するステップとから成る顕著領域検出法の使用を含んでおり、前記ＧＣＣは第１の画像の中のユーザ定義されたクリックポイントに対応する第２のデジタル画像内の点の空間尤度を表すものである、請求項１記載の方法。
【請求項１８】
前記ＧＣＣモデルは純粋な並進、スケーリング及び並進、ならびに相似変換から成る３つの変換クラスのうちの少なくとも１つについて導出される、請求項１７記載のシステム。
【請求項１９】
第１及び第２のデジタル画像の少なくとも一部の間の前記対応は、前記第１のデジタル画像内の第１の点と前記第２のデジタル画像内の第２の点との間の点ごとの対応を確立することを含む、請求項１６記載のシステム。
【請求項２０】
前記第１及び第２のデジタル画像は（ａ）コンピュータ断層撮影（ＣＴ）ボリューム、（ｂ）磁気共鳴画像法（ＭＲＩ）、又は（ｃ）陽電子放射断層撮影法（ＰＥＴ）のうちの１つを用いて取得される、請求項１６記載のシステム。
【発明の詳細な説明】
【技術分野】
【０００１】
関連出願の相互参照
本願は２００６年４月１７日に提出されたOkada他による同時係属中の米国特許仮出願通し番号第６０／７９２，５０７号の本出願である。
【０００２】
発明の分野
本発明は画像の位置合わせを容易にするシステム及び方法、より詳しくは、画像間の高度な空間的対応を達成するために医用画像データ対をアラインするシステム及び方法に関するものである。
【背景技術】
【０００３】
位置合わせは医用画像分析における主要な問題の１つである。単一モード位置合わせの場合、幾つかの部分で異なる画像データ対に最も高い空間的一致が見られるように、この対をアラインする幾何学的変換を修復することを目的としている。変化分析やデータ融合のような臨床上重要な多くのタスクは、このようなデータ対の精確な空間的アラインメントを必要とする。
【０００４】
単一モード位置合わせの従来の解決法は、全体的なミスマッチ誤差を最小化する特定タイプのドメインマッピングを見つけ出すことを目的としている。厳格な位置合わせの場合、このような全体的誤差はドメイン全体にわたって平均された大域的な性質のものである。非厳密位置合わせの場合には誤差を局所的に扱うことができるが、問題を扱いやすくするために幾つかの大域的な正則化条件が用いられることが多い。これらの大域的ファクタを用いれば、相似してないデータ部分の間に一致を確立することが可能である。しかし、同時に、これらの大域的ファクタは、相似していない部分が相似している部分においてさえ全体的な位置合わせに影響を及ぼすことを許容してしまう。さらに、これらの解決法において、位置合わせがどこで正確であり、どこで不正確であるかは、データと費用関数の特定の選択により決定されてしまい、利用可能な臨床上の意味関係や要求がいずれも無視されてしまう。
【非特許文献１】T.Kadir and M.Brady, "Saliency, scale and imgage description", International Journal of Computer Vision, vol.45, no.2, pp.83-105, 2001
【非特許文献２】X.Huang, Y.Sun, D.Metaxas, F.Sauer; C.Xu, "Hybrid image registration based on configural matching of scale-invariant salient region features", in Second IEEE Workshop on Image and Video Registration, in conjunction with CVPR '04, 2004
【非特許文献３】D.Hahn, Y.Sun, J.Homegger, C.Xu, G.Wolz, and T.Kuwert, "A practical salient region feature based 3D multimodality regtistration method for medical images", in SPIE Med.Imag., 2006
【非特許文献４】D.Comaniciu, "An algorithm for data-driven bandwidth selection", IEEE Trans.Pat.Anal. Mach.Intell, vol.25, no.2, pp.281-288, 2003
【非特許文献５】R.Fergus, P.Perona, and A.Zisserman, "Object class recognition by unsuperviesed scale-invariant learining", in IEEE Conf.on Computer Vision and Pattern Recognition, vol.2, 2003, pp.264-271
【非特許文献６】B.Epshtein and S.Ulman, "Identifying semantically equivalent object fragments", in IEEE Conf. on Computer Vision and Pattern Recognition, vol.1,2005, pp.2-9
【非特許文献７】X.Pennec, N.Ayache, and J.Thirion, "Landmark-based registration using features identified through differnetial geometry", in Handbook of Medical Imaging, Academic Press, 2000, pp.499-513
【非特許文献８】K.Okada, D.Cornaniciu, and A.Krishnan, "Robust anisotropic Gaussian fitting for volumetric characterization of pulmonary nodules in multislice CT", IEEE Trans.Med.Imag., vol.24, no.3, pp.409-423, 2005
【非特許文献９】C.Novak, H.Shen, B.Odry, J.Ko, and D.Naidich, "System for automatic detection of lung nodules exhibiting growth", in SPIE Med.Imag., 2004
【非特許文献１０】P.J.Besl and N.D.McKay, "A method for registration of 3-d shapes", IEEE Trans. Pattern Anal.Mach.Intell, vol.14, no.2, pp.239-256, 1992
【非特許文献１１】J.Sill and C.Tomasi, "Good features to track", in IEEE Conf.On Comptuter Vision and Pattern Recognition, 1994, pp.593-600
【非特許文献１２】P.Viola and M.Jones, "Rapid object detection using a boosted cascade of simple features", in Proc. of IEEE Int'l Cont on Comptuter Vision and Pattern Recognition, 2001, pp.511-518
【非特許文献１３】B.Georgescu, X.S.Zhou, D.Comaniciu, and A.Gupta, "Database-guided segmentation of anatomical structures with complex appearances", in Proc. of IEEE Int'l Cont on Computer Vision and Pattern Recognition, 2005, pp.429-436
【発明の開示】
【発明が解決しようとする課題】
【０００５】
本発明の課題は、視覚的に相似していない局所領域間の局所位置合わせを良好な精度と効率性とをもって行うことができるようにすることである。
【課題を解決するための手段】
【０００６】
上記課題は、対象物の少なくとも第１及び第２のデジタル画像を取得するステップと、３Ｄランドマーク検出法を使用して第１のデジタル画像の一部と第２のデジタル画像の一部との間の対応を決定するステップとを有しており、前記第１の画像の一部はユーザ定義されたクリックポイントとしてユーザにより選択された関心点を含んでおり、前記第１の画像の一部は、その位置合わせの精度が第１の画像の他の部分における位置合わせの精度よりも高くなくてはならないような第１の画像内の箇所であり、前記対応は第１の画像内のユーザ定義されたクリックポイントと第２の画像内の所定の点との間の幾何学的関係を用いて決定されることを特徴とする、画像対の視覚的に相似していない局所領域をリンクする方法により解決される。
【０００７】
上記課題はまた、上記方法ステップを実行する機械により実行可能な一連の命令を有体的に実現した機械可読記憶装置、ならびに、上記方法ステップを実行する手段を備えた医用画像データ対を位置合わせするシステムによっても解決される。
【発明を実施するための最良の形態】
【０００８】
本願はロバストクリックポイントリンキングの技法を開示する。クリックポイントリンキングは局所的位置合わせの新しいフレームワークであり、これによりユーザは高い精度が必要とされる箇所を対話的に指定することができる。発明者は、データ対のうちの一方の関心領域の近傍に位置する点位置をユーザ（又は自律エージェント）が指定すると仮定している。ユーザにより提供されたこのような点位置は関心点又は"ＰＯＩ"と呼ばれる。したがって、対話的局所位置合わせのタスクは原データにおける所与のＰＯＩに対応する他データ内の点を見つけ出すことである。本願で発明者が考察する適用シナリオは、同一患者の一連の追跡調査を分析の対象とする長期的な３Ｄデータ調査である。ＰＯＩはユーザのマウスクリックにより与えられると仮定する。したがって、クリックポイントリンキング問題は一方のデータにおけるマウスクリック位置を他方のデータにおける対応する点に自動的にリンクすることにより解決される。このフレームワークは、上で述べた一般的な位置合わせ問題を、（１）特定の点における点ごとの対応を確立すること、及び、（２）これを異なる点に対して逐次実行することとして解釈することを支持する。
【０００９】
対話的局所位置合わせ方式の主な利点の１つは、位置合わせ結果の実際の用いられ方に忠実であることである。臨床上の多くの状況において、医用画像は局所的にしか評価されない。特定の病変や解剖学的構造を評価する際、目標位置における位置合わせの精度は高くなくてはならない。しかし、現場では、他の目標領域が同様に正しく位置合わせされていることにも、データ対が最小平均誤差でアラインされていることにも関心がないことが多い。このような局所的な関心の集中は、ＰＯＩから離れた非目標領域からの影響を無視し、非目標領域の計算を避けることにより、より良い精度と効率性の達成を容易にする。
【００１０】
他方で、このフレームワークの主要な問題は、変化しつつある又は本質的に異なる対応領域をどのようにしてリンクするかということである。治療の前後に撮像された肝臓癌を内容とする追跡データ対を調査するとしよう。治療の効果を定量化するには、データ対の位置合わせと、それに続く変化分析とが必要とされる。これは典型的な循環問題である。といのも、関心対象である時間的変化を分析するには位置合わせが必要であるが、これらの変化そのものが位置合わせを困難にするからである。局所位置合わせは視覚的に非常に異なる局所領域間に対応を見つけ出すというより困難なタスクを要するので、問題をさらに悪くしてしまう。
【００１１】
上記の問題に対処するため、発明者は幾何構成コンテクストを用いた局所位置合わせの解決法を開発した。点ごとの対応は他の安定したデータ点を基準とした所与のＰＯＩの幾何構成を利用することにより確立される。より重要なことに、これは潜在的に信頼のできない局所的な外観／輝度情報を使用せずに行われる。この解決法は各データについて最初に検出されたスケール不変の顕著特徴点（非特許文献１，２，及び３参照）の集合を利用している。したがって、任意のＰＯＩは顕著特徴位置の集合を基準として幾何学的に表現される。幾何構成コンテクスト（ＧＣＣ)は、与えられた特徴の集合に対してＰＯＩの空間尤度をモデル化するガウシアン混合分布として定義される。ＰＯＩの与えられた１つのデータドメインにおいて定義されたＧＣＣは、特徴点集合の粗い対応があれば、他のデータドメインに変換することができる。新しいドメインに変換されたＧＣＣの最大尤度推定値は所望の対応点を提供し、また可変帯域平均シフト法（非特許文献４参照）を用いて効率的に解かれる。
【００１２】
本明細書で開示される位置合わせのフレームワークは部分ベース物体認識の研究における最近の発展に触発されたものである（非特許文献５及び６参照）。最近、EpshteinとUllmanは、我々のモデルに似た尤度モデルを利用して、意味論的には等しいが視覚的には異なる物体の部分を検出する自動アルゴリズムを提案している。我々の仕事は彼らのフレームワークの柔軟なオンライン学習バージョンとして解釈しうるものであり、ブートストラップ尤度を互いに素なデータ集合からでなくテストデータインスタンスから推定する。ランドマークベースの位置合わせ（例えば、非特許文献７参照）も、位置合わせが正確でなければならない点位置がユーザにより与えられることを前提としている点で、上で提案されたフレームワークに関連している。しかし、それらは技術的にも使用上でも異なる目的をもっている。ロバストクリックポイントリンキングが与えられたＰＯＩに関してボリューム対の間の対応を確立しようとするのに対して、ランドマークベースの位置合わせはユーザにより提供された対応の集合に関してドメイン変換を推定することを目的としている。ランドマークベースの位置合わせが所与の対応と交わるスムーズなドメインマッピングの全体を得ようとするのに対して、クリックポイントリンキングは局所的である。本明細書で開示される技術は、平均シフト法アルゴリズム（例えば、非特許文献４参照）を用いて、上記の顕著特徴に基づいた位置合わせ法（例えば、非特許文献２及び３参照）を改良することも目的としている。平均シフト法はコンピュータビジョンにおける追跡とセグメンテーションのための一般的な解決法である。平均シフト法アルゴリズムは医用画像のセグメンテーション問題に適用され成功を収めてきたが（例えば、非特許文献８参照）、本願の発明者は本願こそが平均シフト法アルゴリズムの医用画像位置合わせ問題への初めての適用であると考える。
【００１３】
本願では、（１）純粋な並進、（２）スケーリングと並進、（３）相似変換（スケーリング、並進、及び回転）を考慮した、上記ＧＣＣベースのフレームワークの３つの例が提示される。このシステムの性能は手作業で注釈をつけた１６個の全身ＣＴ追跡データを用いてすでに評価されている。このクリックポイントリンキングの構想は以前に肺小結節検出に関連してすでに提案されたものである（例えば、非特許文献９参照）。しかしながら、本願は、肺小結節のコンテクストを超え、一般的な状況のもとで、視覚的に異なる領域の処理に重きをおいてこの問題を解決することを目指している。
【実施例】
【００１４】
本願は次のように組織されている。以下のセクションでは、ＧＣＣと、純粋な並進のもとでの３Ｄロバストクリックポイントリンキングの問題に対するＧＣＣを用いた解決法とが説明される。セクションＩＩ−Ａでは、ロバストクリックポイントリンキング問題が形式的に定義され、ＧＣＣを用いたその解決法が概観される。本願全体を通して使用される用語及び記号もこのセクションで定義される。セクションＩＩ−Ｂ及びＩＩ−Ｃでは、顕著特徴を抽出し、全数探索法をベースとした戦略を用いてボリューム対間の特徴の対応を推定するために使用されるアルゴリズムが説明される。セクションＩＩ−Ｄでは、所与のＰＯＩに対応する点の空間尤度関数を表現するＧＣＣが提案される。セクションＩＩ−Ｅでは、ＧＣＣの最大尤度推定を解くために使用される最適化法が導入される。次の２つのセクションでは、上記アプローチを陰伏的なドメイン変換のより一般的なクラスへ拡張することが検討される。セクションＩＩＩでは、このような一般化されたフレームワークが説明され、セクションＩＶではさらに、ドメイン変換を陽に推定せずに、相似変換までの複数の変換のもとでどのようにしてこのような解を導き出すことができるのかが詳述される。セクションＶでは、幾何学的アンカーとして使用される顕著特徴に対する代案が説明される。セクションＶＩでは、提案された方法の実行可能性が評価される。
【００１５】
ＩＩ．  ロバストクリックポイントリンキング
Ａ．  定義と概観
位置合わせすべき画像関数の対が与えられているとしよう。一般性を失うことなく、一方は基準画像と呼ばれ、Ｉr（ｘr）で表され、他方は浮動画像と呼ばれ、Ｉf（ｘf）で表される。ここで、ｘr∈R3とｘf∈R3は各々の連続ドメインにおける座標変数を表している。ドメインの対はθをパラメータとする未知の変換Ｔθによって陰伏的に関係づけられているものと仮定される。
【００１６】
【数１】
2007283108.tif 000002
【００１７】
クリックポイントリンキングのタスクは、基準画像Ｉr(ｘr）内の所与のクリックポイント又はＰＯＩ ｃrに対応する浮動画像Ｉf（ｘf）内の点ｃｆの推定として定義される。真のドメイン変換Ｔθが知られていれば、真の解ｃfを定義することができる。
【００１８】
Ｃf＝Ｔθ（ｃr）                                                  （２）
次に、Ｐで表される３Ｄ座標をもつ顕著特徴を導入する。ここで、それぞれ、基準画像についてはＮr個の特徴から成る集合Ｃrが計算され、浮動画像についてはＮf個の特徴から成る集合Ｃfが計算されるものと仮定する。
【００１９】
Ｃr＝｛Ｐｒ１，．．，ＰｒＮr｝                              （３）
Ｃf＝｛Ｐｆ１，．．，ＰｆＮf｝                              （４）
ＣrとＣfから構成されたＭ個の対応する特徴の対Ｑから成る集合をＱで表せば、
Ｑ＝｛（ｑｒ１，ｑｆ１），．．，（ｑｒM，ｑｆM）｝      （５）
ここで、ｑｒi∈Ｃr、ｑｆi∈Ｃf、Ｍ＜ｍｉｎ（Ｎr,Ｎｆ）。
【００２０】
標準的な位置合わせの解法はエネルギー最小化問題
2007283108.tif 000003
を解くことによりドメイン変換
2007283108.tif 000004
を推定することを狙いとしている。例えば、特徴ベースの位置合わせは、
2007283108.tif 000005
となるように
2007283108.tif 000006
とＱを同時に推定する反復最近傍点（ＩＣＰ）アルゴリズムを使用して解くことができる（非特許文献１０参照）。いったんドメイン変換が正しく推定されれば、クリックポイントリンキングは
2007283108.tif 000007
とトリビアルになる。しかし、雑音のあるデータから変換を推定することはトリビアルと言うには程遠い。推定精度は対応の誤差に非常に敏感である。ＩＣＰのような反復的解法はまた計算コストが高くなりがちでもある。
【００２１】
本発明のアプローチによれば、リンキング問題はドメイン変換を陽に推定せずに局所変数の上の空間尤度関数を直に最適化することにより解かれる。
【００２２】
【数２】
2007283108.tif 000008
ここで、Ｌ（ｘf｜ｃr，Ｑ）は、基準画像内のＰＯＩ ｃrと対応する特徴Ｑの集合とを条件とした浮動画像のドメインにおける空間尤度関数を表している。平均シフト法アルゴリズムは良好な計算効率と誤った対応に対する所望のロバストネスとを可能にするが、最大尤度のこの一般的定式化は平均シフト法アルゴリズムの利用を可能にする。以下では、解法の詳細をステップを追って説明する。
【００２３】
Ｂ．スケール不変顕著特徴抽出
特徴をベースにしたマッチング法の研究の道筋は長いこと「どんな特徴を使用すべきか」という問題によって制限されてきた。隠蔽化及び非隠蔽化のもとでの物体追跡に関して、興味深い特徴選択の基準が非特許文献１１において提案されている。この基準は、追跡のための正しい特徴とはまさに追跡装置を最も良好に動作させる特徴に他ならないと述べる。この基準に含まれる直観は使用目的が方法論の選択を決定するというものである。同様の理屈を適用すれば、信頼できる対応を形成するという脈絡で使用されるべき良い特徴とは「一意の」又は「まれな」特徴であるべきである。すなわち、基準画像から特徴が与えられたときに、この特徴がマッチング画像上に対応する複数の特徴を有する尤度が低ければ、この特徴は「良い」と見なされる。この基準を満たす特徴はますます研究されるようになってきており、画像位置合わせや物体追跡の目的で使用されている。特に、スケール不変顕著「領域」特徴検出器は非特許文献１において提案されている。顕著領域は空間的空間とスケール空間の両方において最も高い顕著度を有する局所画像領域として選択される。局所領域の顕著度と最良スケールはエントロピーをベースとした基準に基づいて決定される。顕著領域特徴の画像位置合わせへの適用は、２Ｄ（例えば、非特許文献２参照）と３Ｄ（例えば、非特許文献３参照）の両方において研究されており、領域特徴の主な利点の１つは回転、並進、及びスケールに関する不変性であることが示されている（例えば、非特許文献２参照）。本明細書において、発明者は顕著領域特徴検出器を３Ｄで使用してＣＴボリュームから特徴を抽出している。
【００２４】
まず、ＣＴボリュームＩの各ボクセルｘについて、半径ｓと中心ｘとにより記述されるあるスケールの球領域Ｒ(s,x)における輝度値iから確率密度関数（ｐ．ｄ．ｆ）ｐ（ｉ｜Ｒ(s,x)）が計算される。
【００２５】
【数３】
2007283108.tif 000009
ここで、iは可能なすべての輝度値の集合の中に値をとり、Ｖ（Ｒ(s,x)）は局所領域Ｒ(s,x)の体積を表し、ｙは領域Ｒ(s,x)内のボクセルを表し、σは上記したカーネルベースのノンパラメトリックｐ．ｄ．ｆ推定におけるガウシアンカーネルの幅を指定する定数である（σはすべてのＣＴボリュームに対して例えば１０の一定値に設定してもよい）。
【００２６】
領域の輝度ｐ．ｄ．ｆが与えられれば、その輝度分布の微分エントロピーが下式により定義される。
【００２７】
【数４】
2007283108.tif 000010
ここで、ｉ（ｓ，ｘ）は領域Ｒ(s,x)内の輝度値の範囲を表している。
【００２８】
ｘを中心とした周囲の領域に対する最良スケールＳｘは局所エントロピーを最大化するスケールとして選択される：Ｓｘ＝ａｒｇｍａｘsＨ（Ｒ(s,x)）。したがって、最良スケールをもつ領域の顕著度Ａ（Ｒ(s,x)）は最良スケールＳｘで重み付けされたエントロピー極値とスケール空間内の微分自己相似測度とにより定義される。
【００２９】
【数５】
2007283108.tif 000011
【００３０】
上記の顕著度測度は空間的空間とスケール空間の両方に適用可能であるから、異なる位置及びスケールにおける領域特徴の顕著度の値を比較することができる。
【００３１】
次に、大域的に最も顕著な少数Ｎ（Ｎ＜１００）の領域特徴（各々その中心と最良スケールとにより定義される）を選ぶために、以下の処理ステップが導入される。
【００３２】
特徴抽出
Ａ１  各ボクセル位置ｘについて、ｘを中心とした周囲の領域の最良スケールＳｘと顕著度Ａ（Ｒ(s,x)）を計算する。
【００３３】
Ａ２  顕著度が局所極大値であるボクセルを同定する。ここで、関心顕著領域とはこれらのボクセルを中心とし、かつ最良スケールを有する関心領域である。
【００３４】
Ａ３  局所極大な顕著領域の中から、ＣＴボリュームの領域特徴として最も顕著なＮ個の領域｛ｐi｝（最高の顕著度を有する）を選ぶ。
【００３５】
１つのＣＴボリュームの顕著領域特徴を計算するにはおよそ２．５分かかる。この性能は我々の用途においては許容しうるものである。というのも、特徴はポイントをクリックする前にオフラインで計算されるからである。
【００３６】
臨床上の実務では、一般に１２ビットの正値ＣＴボリュームデータがＣＴ数値の標準的な物理単位であるハウンスフィールド単位（ＨＵ）に変換される。本研究で使用するデータは０〜４０９５の範囲であり、−１０２４のオフセットと１の傾斜をもつＨＵに関係している。特定のタイプの組織を標準的な８ビットのグレイスケールで視覚化するために、一般にＨＵ値に窓が当てはめられる。本研究では、脂肪（−１００〜−５０ＨＵ）のようなある種の柔軟な組織や水分（０ＨＵ）を抑制するために、３０〜２８５ＨＵの特定の窓を使用する。この種の柔軟な組織の除去は特徴抽出プロセスの安定化を助け、異時点間の対応の確立を可能にするより硬い構造への集中が可能になる。
【００３７】
Ｃ．  特徴対応マッチング
基準画像Ｉr（ｘr）と浮動画像Ｉf（ｘf）は双方独立して、前のセクションＩＩ−Ｂで説明したスケール不変特徴を抽出するための手続きを受ける。その結果、式（３）で定義されたＩrのＮr個の特徴から成るＣrと呼ばれる集合と式（４）で定義されたＩfのＮf個の特徴から成るＣfと呼ばれる集合との対が得られる。
【００３８】
基準ドメインｘrのＰＯＩ  ｃrを所与として、式（５）で定義されたＭ個の対応する特徴から成る集合Ｑを、以下の全数探索戦略を用いて見つけ出す。
【００３９】
全数探索による特徴マッチング
Ｂ１  ユークリッド距離の観点でｃrに最も近いＣrからＭ＜Ｎr個の特徴｛ｑr1，．．，ｑrM｝を選択する。
【００４０】
Ｂ２  各基準特徴ｑriについて、
Ｂ２ａ  ＣfからのＮf個の浮動ドメイン特徴｛ｐfj｝について、外観に基づいた相似性をしらみつぶしに計算する。
外観の相似性は、それぞれｑriとｐfjを中心とする円形領域Ｒ（ｓｑri，ｑri）とＲ（ｓｐfj，ｐfj）から導出された輝度ヒストグラムの対の間のカイ二乗距離によって測られる。
【００４１】
Ｂ２ｂ  最も相似したｐfjを選択し、それをｑfiに設定する。
【００４２】
注意しなければならないのは、これは非常に単純なマッチングアルゴリズムであり、粗い結果のみを出すことを意図したものであるということである。上記の全数探索はＩＣＰのようなより複雑な方式に比べれば正確さに劣ることもある。したがって、Ｑは無視できない数の誤った対応を含んでいることがありうる。しかし、その計算の複雑さは他の複雑な方式よりも著しく低くなることが予想され、より効率的な解法の実現が可能になる。
【００４３】
Ｄ．  空間尤度モデリング
このサブセクションでは、どのように空間尤度関数Ｌ（ｘf｜ｃr，Ｑ）をモデル化するかを説明する。このようなモデルがどのように構成されるかは、Ｔθでどのような変換のクラスを考えるのかに依存する。説明のために、まず、極度に単純化された並進のみの変換に対するモデル構成を具体的に示す。より一般的な変換クラスへの拡張はセクションＩＩＩで論じられる。
【００４４】
基準ドメインと浮動ドメインの両方に対して、まず、それぞれの顕著特徴位置に原点の設定された局所フレームを導入する。
ｘr＝ｘri＋ｐri                                            （１０）
ｘf＝ｘfi＋ｐfi                                            （１１）
ここで、ｘriとｘfiはそれぞれ基準ドメイン内の特徴ｐriを中心とする局所フレーム内の座標値と浮動ドメイン内の特徴ｐfiを中心とする局所フレーム内の座標値を表している。
【００４５】
我々は基準ドメイン内の幾何構成コンテクスト、すなわち、ＧＣＣをＭ個の要素分布を有するガウシアン混合分布として定義する。この混合分布モデルは、ＱからのＭ個の顕著特徴｛ｑri｝を基準としたＰＯＩ位置ｃrの空間尤度を表す。
【００４６】
【数６】
2007283108.tif 000012
ここで、Ｉは３次元単位行列であり、ｍriとσriはそれぞれi番目のガウシアン要素分布の平均と幅を表しており、ｃriはｑriを中心とする局所フレームにおけるＰＯＩの３次元座標であり、Sｑriは顕著特徴ｑriに関連した等方性スケールであり、ｆrは特徴のスケールをガウシアン幅に関連付ける関数である。ここで注意すべきことは、基準ドメイン内でＰＯＩがどこにあるかは正確に知られているため、基準ドメイン内のＧＣＣがM個の同心ガウシアン要素分布から成る混合分布を形成するということである。したがって、結果として得られる尤度は凸であり、ＰＯＩ位置において一意の極大値を有する。その一方で、i番目の局所フレーム内の３次元ベクトルｃrは顕著特徴位置ｐriに対するＰＯＩの幾何学的関係を定める。
【００４７】
最後に、ｘrにおける所与のＰＯＩに対応する浮動ドメインｘfにおける点の空間尤度を表す所望のＧＣＣを定義する。
【００４８】
【数７】
2007283108.tif 000013
ここで、ｇtはｃr及びQを特徴ｑfiに関するi番目のガウシアン要素分布の平均に関係付ける関数を表しており、ｆtはSｑri及びSqfiをi番目のガウシアン要素分布の幅に関係付ける関数を表している。
式（１０）及び（１１）で定義されたｑriとｑfiを中心とする局所フレームに関する基盤的変換Ｔθはここでは直線的な剛体並進のみから構成されるものと仮定されていることを想起されたい。ベクトルは純粋な並進のもとでは不変であるから、ｘfiはｘriに等しくなければならない。結果として、式（１８）の導出は、上で述べた純粋な並進の仮定のもとでｑriとｑfiとの間の対応を用いることにより、基準ドメインｘrにおいて式（１３）で定義されたガウシアン要素分布を浮動ドメインｘfへ移すことと解釈しうる。まず始めに、ｚ番目の基準局所フレームにおけるＰＯＩを特定するために式（１０）を使用する。
ｃr＝ｃri＋ｑri                                                  （２０）
【００４９】
次に、平均ｍfiが式（１１）のようにi番目の浮動局所フレーム内で表現されたｃfに等しいものとして定義される。続いて、並進不変性ｃfi＝ｃriを適用して式（２０）への代入が行われ、結果として式（１８）に示されている形が得られる。
ｍfi＝ｃf＝ｃfi＋ｑfi＝ｃri＋ｑfi＝ｃr−ｑri＋ｑfi                （２１）
ガウシアン幅のモデル化のために、我々は前に述べた最大エントロピーアルゴリズムにより導出された各特徴のスケールを点位置同定の統計的不確実性として解釈する。平均推定量（２１）は両方のドメインにおける特徴ｑriとｑfiとの間で均衡しているので、ガウシアン幅は式（１９）に示されているように推定されたスケールの対の不偏平均としてモデル化される。
【００５０】
空間尤度L(ｘf｜ｃr，Ｑ）の混合分布モデル（１６）は、（１２）の場合とは異なり、平均と幅が可変のガウシアン要素分布から成っている。これは抽出された特徴位置に異なる時点にわたって分散を生じさせる測定誤差に因るものである。さらに、Ｑにおける誤った特徴の対応が平均推定値を真の平均から大幅に外れさせることもありうる。スケール推定値Ｓｑri及びＳｑfiが空間的に可変であるため、ガウス幅も可変である。
【００５１】
Ｅ．  平均シフトに基づいたロバスト最大尤度推定
このセクションでは、（１６）〜（１９）において定義された尤度モデルを用いて、（６）における最大尤度推定問題に対する我々のロバストかつ効率的な解法が説明される。セクションＩＩ−Ｃ及びＩＩ−Ｄで論じた特徴マッチングの誤差のため、尤度関数は外れ値モードを成す誤った対応によってマルチモーダルとなる。したがって、解決すべき課題は正しく見つけ出された対応のみによって混合分布を推定することである。言い換えれば、任意の分布における外れ値の中から正しいモードが選び出されなければならない。この課題は非特許文献４において提案されている可変帯域平均シフト（ＶＢＭＳ）法を用いて解決される。平均シフト法は適応的な反復ステップサイズを有する効率的かつ収束の保証された勾配上昇アルゴリズムである。本来のアルゴリズムは与えられたデータ点集合に対してカーネル密度推定値のモード構造を分析するために設計された。この設定では、カーネル帯域幅は空間的に一定であると見なされることが多い。ＶＢＭＳは上記平均シフト法のフレームワークを帯域幅が空間的に可変のケースに拡張したものであり、このケースでは、データ点が異なれば、有意性も異なる。この拡張により、一般的なガウシアン混合分布モデルのモード構造の分析と一般的な情報融合問題の解決への適用が可能になる。ここで、情報融合問題とは、与えられた仮説の集合に対して、最も妥当な解を推定することを課題とする問題のことである。以下では、ＶＢＭＳフレームワークを手短に要約する。
【００５２】
ｘi∈Ｒ3、ｉ＝１，．．，Ｍが３Ｄデータ点の集合を表し、Ｈiが点ｘiに関する不確実性又は有意性を表す３次元行列であるとする。点ｘにおける３Ｄ正規カーネル点密度推定量は次式により与えられる：
【数８】
2007283108.tif 000014
【００５３】
ＶＢＭＳベクトルｍv（ｘ）は（２３）により定義される。
【００５４】
【数９】
2007283108.tif 000015
ここで、Ｈh（ｘ）はｘにおける帯域幅行列のデータ加重した調和平均を表しており、
【数１０】
2007283108.tif 000016
重みｗi（ｘ）はすべての要素分布にわたって正規化されたｘにおけるi番目の要素分布からの影響を表している。
【数１１】
2007283108.tif 000017
【００５５】
ＶＢＭＳベクトルが基礎となる密度の正規化された勾配の適応推定量であることは示すことができる：
【数１２】
2007283108.tif 000018
ＶＢＭＳベクトルを用いた以下の反復アルゴリズムは、面倒なパラメータ調整
ｙ0＝ｘinit
yn+1＝ｍv(yn）＋ｙn                                              （２６）
なしに、初期値ｘinitの近傍における密度推定値のモードに勾配上昇方向で収束することが保証されている。
【００５６】
イテレータの収束値はｙ*で表される。
【００５７】
上記ＶＢＭＳアルゴリズム（２６）を用いてマルチモード尤度モデル（１６）の最大尤度を推定するために、２つのロバストなアルゴリズムが使用される。（１６）における特定のモデルへのアルゴリズムの適用は直接的であり、それぞれ（１８）と（１９）において定義されているように単にｘi＝ｍfi，Ｈi＝σfi2Iと設定するだけである。第１の解法であるシングルＶＢＭＳは純粋な並進の仮定のもとでCrとCｆとから推定された初期値ｘinitからの単一のＶＢＭＳ反復を必要とする。その一方で、第２の解法であるマルチプルＶＢＭＳは各要素分布の平均ｍfiにおいて初期化された複数のＶＢＭＳ反復からの収束値の集合からの投票を必要とする。
【００５８】
シングルＶＢＭＳ：
Ｃ１  ＣrとＣfにおける顕著特徴点の平均ｚrとｚfをそれぞれ計算する。
【００５９】
Ｃ２  ＣrとＣfの間の平均バイアスｚ＝ｚf−ｚrを計算する。
【００６０】
Ｃ３  浮動ドメイン内の平均バイアス補正されたＰＯＩによりＶＢＭＳイテレータの初期値を設定する：ｘinit＝ｃr＋ｚ。
【００６１】
Ｃ４  （２６）におけるＶＢＭＳアルゴリズムを実行し、収束値ｙ*を得る。
【００６２】
Ｃ５  リンキング推定値
2007283108.tif 000019
を得る。
【００６３】
マルチプルＶＢＭＳ：
Ｄ１  Ｍ個のＶＢＭＳイテレータを（１８）におけるＭ個の要素分布の平均｛ｍfi｝で初期化する：ｘinit,i＝ｍfi。
【００６４】
Ｄ２  Ｍ個の独立したＶＢＭＳ反復（２６）を実行し、収束値の集合｛ｙi*｝を得る。
【００６５】
Ｄ３  点対間ユークリッド距離に従って｛ｙi*｝をグループ分けし、最も多くの要素を含むクラスタを形成する部分集合｛ｙk*｝を選ぶ。
【００６６】
Ｄ４  部分集合の平均によりリンキング推定値を得る：
【数１３】
2007283108.tif 000020
【００６７】
前者の解法は単一のＶＢＭＳ反復しか必要としないため効率性を重視している。しかしながら、推定精度は推定された初期値に大きく依存している。この解法は、ボリューム対の間に顕著特徴の大きな変化と純粋な並進を越えた変形とが存在する場合には有効でない。後者の解法は複数のＶＢＭＳ反復を行うことによってロバストネスを改善する。この解法は前者の解法よりも効率性の点で劣るものの、特定のドメイン変換を仮定しないため、より一般性がある。
【００６８】
ＩＩＩ．  一般の変換クラスへの拡張
セクションＩＩ−Ｄで述べたように、浮動ドメインにおける空間尤度関数（１６）は基礎となるドメイン変換としてどのようなタイプのドメイン変換を考えるのかに依存している。例えば、（１８）及び（１９）におけるガウシアン要素分布の平均と幅の式は純粋な並進を越えた変換のもとでは成り立たない。これは（２１）に示されているような（１８）の導出が変換クラスに固有の不変性を利用しているからである。このセクションでは、尤度のモデル化に使用される戦略をどのようにしてより一般的な変換クラスに拡張するかを論じる。
【００６９】
ＱはセクションＩＩ−Ｃで説明した手続きを用いて計算された基準ドメインと浮動ドメインとにおける顕著特徴の間のM個の対応から成る集合を表していたことを想起されたい。各対応ｑiは２つ組（ｑri，ｑfi）により表現される。Pは
【数１４】
2007283108.tif 000021
であるようなＱのすべての部分集合から成る集合を表すものと仮定する。
ここで、ＬはＰのカーディナル数であり、ＰlはＱのＫ部分集合である。一般性を失うことなく、Ｐlから任意に１つの対応ｑkを選び、それをアンカー対応ｑl＝（ｑrl，ｑfl）と呼ぶことができる。これにより、各Ｋ部分集合｛ｑl｜ｌ＝１，．．，Ｌ｝に対してアンカーの集合が得られる。
【００７０】
上記設定を所与として、（１６−１９）で定義された空間尤度関数をより一般的な線形変換のクラスに拡張することを検討する。一般化された空間尤度はセクションＩＩ−ＤIで説明した純粋な並進のケースに似てＬ個のガウシアン要素分布の混合分布として定義される。
【００７１】
【数１５】
2007283108.tif 000022
【００７２】
これらは特定の変換クラスに依存する関数ｇgとｆgを除いて（１６−１９）に似た形をとっていることに注意されたい。尤度関数も集合Ｑに代わって集合Ｐに従属している。したがって、各ガウシアン要素分布はＫ個の特徴の対応から成る集合Ｐlに依存する。
【００７３】
ここで、基礎となるドメイン変換Ｔgはある種の一般的なパラメトリックモデルによりモデル化することができ、ドメイン変換を一意に決定するにはＰlにおいて与えられたＫ個の対応から成る集合で十分であることを前提とする。すると、下記によって変換クラスに対してｇg及びｆgを決定する一般的な手続きが提供される。
【００７４】
所与のＰｌに対するl番目のガウシアン要素分布の平均と幅：
Ｅ１  Ｐlを用いて変換パラメータ
2007283108.tif 000023
を推定する。
【００７５】
Ｅ２  ｇｇ：
2007283108.tif 000024
によりl番目の要素分布の平均を計算する。
【００７６】
Ｅ３  Ｐl内のすべての特徴点のスケールの関数としてl番目の要素分布の幅を計算する。
【００７７】
上記の尤度の定式化は、適切なＫ値を選ぶことによりドメイン変換の異なるクラスを扱う一般的なフレームワークを提供する。自由度（ＤＯＦ)の大きな複雑な変換を扱うには高いＫ値が必要である。例えば、Ｒ3では、Ｋ＝３の尤度関数は相似変換までの変換を扱うことができる。Ｋ＝２の尤度はスケーリングと並進による変換しか決定することができない。Ｋ＝１の場合には、純粋な並進しか扱うことができない。（１６−１９）の空間尤度関数が一般的な尤度（２８−３１）のＫ＝１のときの具体例であることは容易に示される。このような状態においては、（２７）で定義された集合ＰはＱに還元されるので、PＰ＝ＱかつＬ＝Ｍである。（１８）及び（１９）における元々の平均と幅の式はＥ１−Ｅ３の手続きを満たし、Ｋ＝１の場合には、ｇg−ｇtとｆg−ｆtに等しくなる。
【００７８】
関数ｆg（ｃr，Ｐl）はl番目のガウシアン要素分布の幅を決定する。セクションＩＩ−Ｄで論じたように、Ｐlにおける顕著特徴のスケールＳrkとＳfkは特徴点の位置を同定する際の統計的不確実性として解釈できる。特徴位置におけるこのような不確実性は変換Ｔｇを介して伝播させられることにより要素分布の推定された平均ｍflの位置へ移されなければならない。しかし、ドメイン変換による変形が過度に大きくないことを仮定して、伝播による変化を無視できるものとして扱えるようにしてもよい。それゆえ、特徴における不確実性は要素分布の推定された平均として扱うこともできる。本研究では、不確実性の平均を不確実性集合の融合の簡単な例として論じる。
【００７９】
【数１６】
2007283108.tif 000025
【００８０】
ＩＶ．  不変性に基づいた陰伏変換の推定
前セクションで説明した手続きＥ１−Ｅ３はＱにおける元の対応の各Ｋ部分集合Ｐlから変換パラメータθlを陽に推定することを要する。ｃrとｃfならびにｍrlとｍflの間にリンクを確立するという我々の目的にとって、推定問題のより直観的な幾何学的解釈が存在しており、この解釈がより効率的な解法へと導いてくれる。このセクションでは、Ｅ１及びＥ２のような陽的なパラメータ推定をせずに幾何学的不変性に基づいてl番目のガウシアン要素分布の平均を特定する代替的な解法を提示する。
【００８１】
あるＫ値とこれに対応する許容される変換のクラスに対して、特定の変換クラスのもとで真である幾何学的不変性を用いてｍflの推定量を構成する。以下がその設定である。（１０）及び（１１）で定義されているように、基準局所フレームと浮動局所フレームの対を考える。これら２つの局所フレームの原点を相互に対応するアンカー点ｑrlとｑflに設定する。セクションＩＩＩに記載されている定義により、基礎となるドメイン変換Ｔθlを決定するには、Ｐlにおいて与えられたK個の対応から成る集合で十分である。変換Ｔθｌは所与のＰlに関して陰伏的に推定されると仮定する。すると、l番目のガウシアン要素分布の平均ｍflは、陰伏的な推定値において与えられた制約条件のもとで所与のＰＯＩ ｃrに対して
2007283108.tif 000026
を与える関数ｇgにより定義される。：
2007283108.tif 000027
【００８２】
２つの重要な問題が存在している。１つは、ここで考察している変換
2007283108.tif 000028
が局所フレーム内で定義されているため、手続きＥ１−Ｅ３の場合とは異なり、各部分集合Ｐlごとに変換の各インスタンスが互いに異なることである。しかし、このことは何ら余分な計算コストを付加するものではない。というのも、手続きＥ１−Ｅ３もまた各Ｐlごとに大域フレームｘfで変換を再計算するからである。この組合せ的計算はＱにおける推定された対応における誤差を処理するためにも必要である。２つめは、大域フレームにおけるｘrとｘfとの間の並進はｑlにおける対応点をこれらのフレームの原点として設定することにより局所フレームの対に吸収されてしまうということである。したがって、局所フレームにおける
2007283108.tif 000029
のＤＯＦは相応して調整されなければならない。
【００８３】
本明細書は相似変換までのドメイン変換をカバーするＫ＝１，２，３の場合の解法を提供する。Ｋ＝４の場合はＲ3におけるアフィン変換をカバーするが、その解法は将来の仕事として残しておく。Ｋ＝１の場合については、セクションＩＩ−Ｄにおける導出において既に１つの解法が提供されている。Ｋ＝１という条件はＲdにおける純粋な並進しか許容しない。さらに、ベクトルはＲdにおけるこの並進という変換クラスのもとで不変である。局所フレームの対は等しく、したがってｃfl＝ｃrlであることを理解するのは容易である。以下では、Ｋ＝１，２の場合について検討する。
【００８４】
Ａ．  スケーリング及び並進変換
Ｒ3では、Ｋ＝２の場合のＫ部分集合Ｐlは２つの対応を生じさせ、６つの制約条件を提供する。スケーリングと推進とによる変換（４自由度（ＤＯＦ））と純粋な並進（３ＤＯＦ）を決定するには、これらの制約条件で十分である。
【００８５】
まず、式（１０）と（１１）を用いて、それぞれアンカーｑrl、ｑflを中心とするl番目の基準局所フレームと浮動局所フレームにおける所与のＰＯＩとその対応する点を特定する。
ｃr＝ｃrl＋ｑrl  ｃf＝ｃfl＋ｑfl                                  （３４）
ここで、ｃflは未知数であり、既知数ｃrl，ｑl及びＰlの関数として決定されなければならない。
【００８６】
Ｋ＝１の場合と同じ論法を用いるならば、並進要素をキャンセルしてしまえば、スケーリングが局所フレームｃrlとｃflとの間の座標マッピングにおける唯一の可変要素として残る。Ｐlは２つの対応しか含んでいないので、アンカーを選んだ後には残りの１つの対応しか存在しない。ｑla＝（ｑrla，ｑfla）がこの残りを表すものとしよう。また、arlとaflがそれぞれの局所フレームにおける相対ベクトルarl＝ｑrla−ｑrlとafl＝ｑfla−ｑflを表すとしよう。
【００８７】
すると、Ｐlにおける対応の対と、ｃrlとｃflとの対は、一対の相似三角形（０，ａrl，ｃrl）と（０，ａfl，ｃfl）として解釈することができる。なおここで、これらの相似三角形は回転なしの異なる大きさの相似三角形である。したがって、この解釈はスケーリングのもとでの下記の２つの不変性を提供する：正規化されたベクトル
【数１７】
2007283108.tif 000030
及びベクトルノルムの比
【数１８】
2007283108.tif 000031
ここで、||・||はベクトルノルムを表している。（３５）、（３６）及び（３３）を組み合わせれば、K＝２の場合のl番目のガウシアン要素分布の平均を推定する所望の関数が得られる。
【００８８】
【数１９】
2007283108.tif 000032
【００８９】
Ｂ．  相似及びユークリッド変換
Ｒ3では、Ｋ＝３の場合のＫ部分集合Ｐlは３つの対応を生じさせ、９つの制約条件を提供する。相似変換（７ＤＯＦ)とユークリッド変換（６ＤＯＦ)までの射影変換を決定するには、これらの制約条件で十分である。局所フレームにおける部分集合内の２つの点を結ぶベクトルの周りの３Ｄ回転の多義性のため、２部分集合で与えられる６つの制約条件では６DOFのユークリッド変換を一意に決定するには不十分である。
【００９０】
ｑla＝（ｑrla，ｑfla）とｑlb＝（ｑrlb，ｑflb）がPlからアンカーｑlを選んだ後の残りの２つを表すものとしよう。また、arlとaflがそれぞれの局所フレームにおける相対ベクトルarl＝ｑrla−ｑrlとafl＝ｑfla−ｑflを表すとしよう。同様に、ｂrlとｂflがそれぞれ相対ベクトルｂrl＝ｑrlb−ｑrlとｂfl＝ｑflb−ｑflを表すとしよう。
【００９１】
すると、Ｐlにおける対応とｃrlとｃｆｌとの対応の３つの対応を一対の相似四面体（０，ａrl，ｂrl，ｃrl）と（０，ａfl，ｂfl，ｃfl）として解釈することができる。定義により、幾何学的相似性は、対応するそれぞれの角が等しく、かつ、各辺が同じ係数でスケーリングされていることを保証する。以下では、上記四面体の幾何学的相似性を利用して、未知数ｃflの閉形式の式を他の既知数の関数として導出する。
【００９２】
基準四面体の表面を含んだ３Ｄ平面Ｌrlを考えてみよう。ただし、基準四面体とは基準局所フレーム内の２つのベクトルａrlとｂrlを含んでいるものである。同様に、表面（０，ａfl，ｂfl）を含んだ浮動局所フレーム内の３Ｄ平面をＬflで表すことにしよう。次に、ベクトルｃrlとｃflをそれぞれ平面ＬrlとＬflへ直交射影する。このような直交射影をｕrlとｕflで表すことにする。その結果、ベクトルｃrlとｃflは直交射影と上記平面に対して垂直なベクトルとのベクトル和に線形分解することができる。
ｃrl＝ｕrl＋ｖrl                                                  （３８）
ｃfl＝ｕfl＋ｖfl                                                  （３９）
ｖrl＝ｋrｎｒl    ｋr∈Ｒ                                        （４０）
ｖfl＝ｋfｎfl    ｋf∈Ｒ                                          （４１）
ここで、ｎrlとｎflはそれぞれＬrlとＬflの単位法線を表している。直交射影により、ベクトルｖrlとｖflはそれぞれ（ａrl，ｂrl）と（ａfl，ｂfl）に直交し、ｕrlとｕflはそれぞれＬrlとＬflに含まれることが保証される。
【００９３】
ｕrlはＬrl内にあるのであるから、ａrlとｂrlの線形結合で書き表すことができる。重みは上記の直交性制約条件を用いて陽に解くことができる。
【００９４】
【数２０】
2007283108.tif 000033
【００９５】
ｃrlは既知であるから、ｕrlが与えられれば、ｖrlは（３８）から求められる。
ｖrl＝ｃrl−ｕrl＝ｃrl−ｗaａrl−ｗbｂrl                          （４５）
浮動四面体の場合、ｕflは（４２）のｕrlのようにａflとｂflの線形結合としても書き表すことができる。さらに、幾何学的相似性のため、浮動局所フレーム内のｕflも（４３）及び（４４）と同じ重みで決定される。
ｕfl＝ｗaａfl＋ｗbｂfl                                            （４６）
したがって、Ｌflの単位法線もａflとｂflから導き出すことができる。
【００９６】
【数２１】
2007283108.tif 000034
ここで、ａfl＝（ａfl1，ａfl2，ａfl3）Tであり、、ｂfl＝（ｂfl1，ｂfl2，ｂfl3）Tである。以下のベクトルノルムの比が相似変換のもとで不変であることは明らかである。
【００９７】
【数２２】
2007283108.tif 000035
【００９８】
（４１）と（４９）を組み合わせれば、陽形式のサイズ係数ｋfが得られる。
【００９９】
【数２３】
2007283108.tif 000036
【０１００】
最後に、（４６）と（４１）を（３９）に挿入すれば、K＝３の場合のl番目のガウシアン要素分布の平均を推定する所望の関数が得られる。
ｍfl＝ｇg，Ｋ＝３（ｃr，Ｐl）＝ｗaａfl＋ｗbｂfl＋ｋfnfl          （５１）
ここで、ｗa，ｗb及びｋfはそれぞれ（４３）、（４４）及び（５０）において与えられている。そして、ｎは（４７）及び（４８）において与えられている。
【０１０１】
Ｖ．  学習に基づいた解剖学的ランドマーク検出
ロバストな融合のための幾何構成コンテクスト（ＧＣＣ）方式において、クリックポイントの対応点はＮ個の最近傍顕著領域特徴の対応に基づいて推論される。これらの最近傍顕著領域特徴はクリックポイントのコンテクスト特徴と考えることができる。したがって、クリックポイントに関する良好な対応を導き出すためには、そのＮ個のコンテクスト特徴対応の主要部分が正確であることが重要である。顕著領域に基づいた特徴は並進、回転、及びスケーリングに対して不変であることが理論的に証明されているが（非特許文献１参照）、構造の出現／消滅がある場合には、変化した構造の近傍で検出された特徴は劇的に変化することがありうる。例えば、大きな腫瘍のある患者を異なる３時点にイメージングした場合、画像の外観は腫瘍の変化により変化し、腫瘍の近傍で検出された顕著領域特徴は３つの画像において安定していない。ユーザが時間の経過とともに急速に変化及び進展する傾向のある患部近傍にある関心点をクリックすることもよくあることを考慮するならば、患部近傍で検出された顕著領域特徴はロバストな融合アルゴリズムにとってさえ厄介となりうる。
【０１０２】
事前学習した判別分類子を用いて検出された解剖学的ランドマークは顕著領域特徴の代案となりうる。この学習フレームワークは非特許文献１２及び非特許文献１３によるリアルタイム検出アルゴリズムに基づいている。
基本的なアイデアは、画像のデータベースを使用し、このデータベースから、学習に基づいたアルゴリズムが特に解剖学的構造と構造間の大域的な解剖学的関係とにおける変化を表すものを抽出することができるようにするというものである。この方法は全身ボリュームを訓練データとしてデータベースに収集することから始まる。次に、各訓練ボリュームについて、関心ランドマークに手作業でラベル付けする。特定のランドマークについての判別分類子を学習するために、まず、関心ランドマークに他のランドマークよりも大きな重みを与える重み付けアラインメントアルゴリズムを用いて訓練ボリュームのアラインメントを行う。ボリュームのアラインメントが済むと、ランドマークを含んだ正事例を形成するために関心ランドマークの周囲のデータがクロッピングされ、負事例を形成するために訓練ボリュームの他の箇所がランダムにクロッピングされる。つぎに、正及び負の訓練事例のデータベースを使用することにより、カスケード型ＡｄａＢｏｏｓｔ法を用いた判別分類子が学習される。実行時には、新たな全身ボリュームが走査され、分類子によって正として分類された窓が可能性のある検出結果として返される。検出された窓はまたそれらの応答に基づいてランク付けされ、最も応答値の高い窓が新たな全身ボリュームにおける検出されたランドマークである尤度が最も高い。
学習をベースとした方法を用いて検出されたランドマークは、特に構造の出現又は消滅がある場合には、顕著領域特徴よりも安定している。なぜならば、これらのランドマークは事前定義されたものであり、訓練プロセスがこれらのランドマークの近傍にある構造の変化をある程度吸収するからである。検出されたランドマークのロバストネスも保証され、また正負検出は粗密検出の実施とランドマーク間の幾何学的制約条件の利用とにより除くことができる。我々のこの研究では、訓練データとして４６個の全身ＣＴボリュームにおけるランドマークを手作業でラベル付けした。次に、重み付けアラインメントを用いてボリュームのアラインメントが行われ、正及び負の事例が収集され、カスケード型ＡｄａＢｏｏｓｔ学習器に入力として与えられる。実行時には、１４個のランドマークについて８７％の平均検出率が得られる。そのあと、検出されたランドマークをＧＣＣロバスト融合のためのコンテクスト特徴として使用することもできる。
【０１０３】
ＶＩ．  実験研究
Ａ．  方法
１６個の全身ＣＴボリューム対を用いて上記アルゴリズムの３Ｄ実施をテストすることにより、提案されたフレームワークの実行可能性の評価を行った。各対の２つのボリュームは異なる時点において同じ患者に対して行われたスキャンである。各対の間で同じスキャナプロトコルを使用した。５１２×５１２の体軸スライスのスタックによる元々のボリュームを１２８×１２８のスライスにダウンサンプリングした。各対のうちの一方を任意に選んで基準ボリュームとし、他方を浮動ボリュームとした。
【０１０４】
実際に実施してテストしたのは、提案されたアルゴリズムの以下の特定の構成である。各ボリュームについて、５０個の顕著特徴が事前計算される：Ｎr＝Ｎf＝５０。各対に対して、セクションＩＩ−Ｃで説明した特徴マッチングアルゴリズムが各クリックポイントｃrに最も近い１０個の特徴を基準特徴として実行される：Ｍ＝１０。
幾何学的ユークリッド距離と輝度ヒストグラムのカイ二乗距離の２つの相似関数が検討される。（１）セクションＩＩ−ＤのＫ＝１のときの純粋な並進と、（２）セクションＩＶ−ＡのＫ＝２のときのスケーリング及び並進の場合について、２つのＧＣＣ解がテストされる。セクションＩＩ−Ｅに記載されているように両方のＧＣＣを最適化するために、シングルＶＢＭＳ法を使用する。
【０１０５】
テストのために、セクションＶに記載されているように、手作業でラベル付けした判別ランドマークを使用した。それぞれの人について、骨盤、肺、腎臓、及び鎖骨を含む重要な解剖学的ランドマークに散らばった１４個のランドマークが存在した。各対について、基準ボリューム内のこれら１４個の点をクリックポイントとして使用し、推定されたリンクｃfと浮動ドメイン内のグラウンドトゥルースランドマークとの間のユークリッド誤差を計算した。その結果、全部で２２４（１４個のランドマークに関して１６人の患者）のテスト事例が得られた。ＧＣＣを用いてクリックポイントリンキングを行った後、テンプレートマッチングに基づいた精緻化により浮動ドメイン内の推定されたクリックポイントを精緻化することも考慮した。球形テンプレートのサイズは最大エントロピー基準を用いて自動的に評価した（例えば、非特許文献２参照）。
【０１０６】
Ａ．  結果と考察
図１は本発明の実験の結果を示したものである。誤差は対応するグラウンドトゥルースとＲ3における推定値との間のユークリッド距離として計算されている。左側の例のプロットでは、異なる複数のランドマークにわたって平均した１６人の患者の関数として誤差が示されている。一方、右側のプロットでは、異なる複数の患者にわたって平均した１４個のランドマークの関数として誤差が示されている。特徴対応マッチングのために、２つの相似関数を検討した。一方は平均バイアス調整付きの幾何学的ユークリッド距離であり、他方は幾何学的距離と輝度ヒストグラムのカイ二乗距離を用いた外観をベースとした距離との線形結合である。一番上の行には幾何学的距離に関する結果が示されている。トータルの平均誤差と中央誤差はそれぞれ４．２３ボクセル、３．５０ボクセルであった。次に、自動スケール選択を備えたテンプレートマッチングを使用してポストリファインメントを行った。その結果は真ん中の行に示されている。平均誤差と中央誤差はそれぞれ４．３９、３．２４であった。最後に、外観をベースにした相似性とポストリファインメントの結果が一番下の行に示されている。平均誤差と中央誤差はそれぞれ４．６８、３．１０であった。
【０１０７】
全体的に、平均誤差は３〜５ボクセルの範囲内にあり、提案された方法の実行可能性を実証するものであった。この結果はまた、精度は患者に強く依存するが、ランドマークにはそれほど強く依存しないということも示している。目視での観察により、比較的大きな誤差（例えば、患者７及び１４）は主に大量のミスマッチした特徴に起因する外れ値によって引き起こされていることが明らかになった。外観に基づいた相似とポストリファインメントの使用により精度は少しばかり改善された。しかし、この改善は小さく、実際には外れ値誤差をより悪化させたため、平均誤差は実際には悪化した。インライアー（inliers）については、ポストリファインメントを行った場合、平均誤差は３ボクセルよりも小さかった。
【０１０８】
ＶＩＩ．  結論
本願はロバストクリックポイントリンキングの新しいフレームワークを開示している。クリックポイントリンキングは所与の関心点に対応する点を推定するための局所化された位置合わせとして解釈できる。変化する腫瘍のような視覚的に相似していない局所領域をリンクするロバストな解を導き出すために、我々は相似変換までのドメインマッピングのクラスのもとでリンク点の空間尤度を表すガウシアン混合分布に基づいた幾何構成コンテクストを提案した。可変帯域平均シフト法に基づいた最適解はマルチモーダル混合分布の１つのモードにおいてロバストかつ効率的にリンク点を見つけ出すために開示されている。我々の実験研究により、手作業でラベル付けした全身ＣＴデータセットを用いる提案された手法の有望性が実証された。
【０１０９】
開示されたロバストクリックポイントリンキングのフレームワークを実施すれば、前述の癌治療モニタリングのための長期３Ｄデータ調査における放射線医の典型的な作業の流れに忠実な直観的な図表的視覚化を用いた医用画像のブラウジングと分析のためのコンピュータ支援診断ツールが提供される。
【０１１０】
典型的な医用画像／ボリュームブラウジングシステムは医用画像のうちの少なくとも一対を並べて視覚化する。このような視覚化方式により、さまざまな画像にわたって画像の任意の３Ｄ位置を位置指示子（例えば、マーカー、十字線）によって独立に指定することが可能になる。最終的に、ユーザはユーザ制御デバイス（例えば、マウス及びキーパッド）により制御される位置ロケータ（例えば、カーソル）を用いていつでも自由に選択した任意の画像の任意の箇所を探索及びブラウズすることができる。
【０１１１】
ロバストクリックポイントリンキングはモニタに表示された複数の画像にわたって位置マーカー／十字線を結合する自動的方法を提供する。ユーザが特定の画像をアクティブにして、ブラウジング後に関心点を見つけた場合、ユーザはマウス制御されたカーソルを用いてその箇所をクリックすることができる。これが開示されたリンキングアルゴリズムを開始させ、このリンキングアルゴリズムが他の画像におけるマーカーを自動的に対応する箇所に定位させ、「粘着マーカー」と呼ばれうる作用を実現する。これにより、ユーザは手動でブラウズしなくても、リンクされたマーカーによって指定された対応する病変に直接赴き、分析することができる。このことが、対応する領域を見つけ出すために続く他の画像を手作業で再ブラウズする必要性をなくすことにより、ユーザ／放射線医の作業の流れを改善する。これらの粘着マーカーは異なる関心位置において次々とマウスをシングルクリックすること（すなわち、クリックアンドルック）により対話的な作業の流れを容易にする。
【０１１２】
本明細書に記載された発明は、例えば、命令のプログラムをこの命令を実行することのできるマシーンによって読むことのできるコンピュータ可読プログラム記憶媒体上に有体的に実現することにより自動化してもよい。汎用コンピュータはこのようなマシーンの一例である。適切な記憶媒体の例は当該分野では周知であり、読み書き可能なＣＤ、フラッシュメモリチップ（例えば、サムドライブ）、さまざまな磁気記憶媒体などのようなデバイスを含む。
【０１１３】
以上、本発明の特徴を開示したが、さらなる変更形態は当業者には明らかである。このような変更形態のすべては添付した請求項の範囲内にあると考えられる。本発明の真の範囲を示すものは以上の明細書ではなく、添付した請求項が参照されなければならない。
【図面の簡単な説明】
【０１１４】
【図１】本発明の技術を用いた医用画像ボリューム対の位置合わせの実験結果を図式的に示す。
【図１】
2007283108.tif 000037