(57)【特許請求の範囲】
【請求項１】
生の像データと物体の構成要素の性質のデータベースであって、前記構成要素の内部運動特性、及び隣接する前記構成要素間の相対的運動特性を具備するデータベースを入力する入力手段と、
（ａ）（i）前記像中のそれぞれの前記像要素中の像の局所統計的性質を推算し、（ii）前記像要素の局所統計的性質における類似性に従い、前記像要素を結合して領域を形成し、前記像の一つの像を、それぞれ前記物体の構成要素に対応する複数の前記領域に分解し、（ｂ）（i）前記データベースに従い、それぞれの前記領域の運動を推算し、（iii）前記運動に従って、すべての像における前記領域の位置を決定して、すべての前記像における前記領域を定義するプロセス手段と、
すべての前記像に配置された前記領域により、前記像を可視化し、仮想現実表示をするための図式表現手段
とを具備し、時間的に互いに分離し、且つそれぞれ複数の像要素を含む、物体の像の一連のシーケンスからなる生の像データから時間的空間的像を再構築するシステム。
【請求項２】
前記相対的運動特性は、バネ定数（ｓｐｒｉｎｇ  ｃｏｎｓｔａｎｔ）で表現されることを特徴とする請求項１記載のシステム。
【請求項３】
前記プロセス手段は、前記一連のシーケンス中の連続した像のそれぞれのエネルギーを最小にする手段を含むことを特徴とする請求項１記載のシステム。
【請求項４】
前記推算される運動は、剛体運動及び変形運動を含み、前記エネルギーは、前記剛体運動のエネルギー及び変形運動のエネルギーを含むことを特徴とする請求項３記載のシステム。
【請求項５】
前記プロセス手段は、前記最小化されたエネルギーに従って、前記像要素のそれぞれの速度を推算する手段を含むことを特徴とする請求項３記載のシステム。
【請求項６】
前記プロセス手段は、
（Ａ）少なくとも一つの閾値に従い、像を結合し、繰り返し領域を形成する手段と、
（Ｂ）前記繰り返し領域からコスト関数（cost  functional）を導出する手段と、
（Ｃ）該コスト関数を最適化する手段
とを含むことを特徴とする請求項１記載のシステム。
【請求項７】
前記コスト関数が減少しなくなるまで、前記少なくとも一つの閾値を調整することにより、前記コスト関数を減少させる操作を繰り返して最適化されることを含むことを特徴とする請求項６記載のシステム。
【発明の詳細な説明】
【０００１】
（発明の背景）
本発明は、生のイメージデータを解析するために対象物の既に知られている特徴を利用して、対象物の４次元（空間の３次元と時間の１次元）の生のイメージデータから対象物を再構築及び可視化するためのシステム及び方法に関する。より詳しくは、本発明は、対象物の時間系列のビットマップイメージから、対象物を再構築及び可視化するためのシステム及び方法に関する。１つの実施形態では、本発明は、より具体的に、人体のパーツのＭＲＩ（磁気共鳴映像法）データのイメージングに用いられるシステム及び方法に関する。動的な人体の”ディジタルクローン”即ち仮想現実レプリカを創造して、仮想（ディジタル）空間に４次元の明瞭な器官、組織及び骨が再構築される。
【０００２】
骨及び軟組織構造体の両方を含む、人間の関節の動作パターンを評価する能力は、骨関節症などの症状の発生及びステージングの理解を大いに増すであろう。人間の関節の解剖学的構造のバイオメカニクス及び運動力学に関する多くの情報が存在するが、その情報は主に、死体の研究や２次元イメージに基づく研究に基づいている。今までのところ、生体内の３次元の関節運動を表示することに成功したシステムはなかった。
【０００３】
３次元空間でハイコントラストの解剖学的構造を再構築及び可視化する技術は知られているが、その結果である再構築したものは静的であり、時間をベースにした情報は何ら提供しない。診断され治療されるべき多くの問題は、関節のような複雑な構造体におけるアブノーマルな動作運動学的事象に関連する；それゆえ４次元（時間と空間の両方）でのイメージングが望ましい。
【０００４】
複雑な軟組織や骨組織の観察に用いる生のイメージデータの最も良い現在の情報源は、磁気共鳴映像法（ＭＲＩ）である。しかしながら、ＭＲＩはしばしば下記の技術的な困難を持ち込む。なぜならイメージ化されるべき筋骨格構造の多くは、小さくそして込み入っているので、可視化されるべき解剖学的構造の多くは高分解能を必要とし、低いコントラストを提示する。ＭＲＩは、電磁場を発生させるために局所電場コイル（ｌｏｃａｌ  ｆｉｅｌｄ  ｃｏｉｌｓ）を使用する；そのような局所電場コイルは不均一なイルミネーションフィールドを形成する。ＭＲＩイメージはまた、ノイジーである可能性がある。ＭＲＩイメージデータを用いる再構築及び可視化のための技術はいずれもそのような困難を乗り切らなければならない。
【０００５】
変形動作の追尾及び解析の領域ではいつくかのリサーチがなされているが、そのようなリサーチの多くは、シンプルパラメトリック変形モデルに基づくイメージにおける単一組織の時間的発展に集中している。関節の運動力学的な解析は、多くの組織の運動と絡み合っている、従って従来のアプローチは適さない。
【０００６】
更に、３Ｄ再構築及び可視化のための既知の技術は、下記のような不都合を呈する。第一に、そのような既知の技術は、時系列的イメージに実際に用いるには計算的に強烈すぎる。第二に、３Ｄ再構築及び可視化のための最も知られているアルゴリズムは、監督され（つまりオペレータの干渉を必要とし）そして放射線医師などの専門的知識に頼る。そのような監督は、そのようなアルゴリズムの適用範囲を、骨や心臓や海馬などのいくつかの特定の解剖学的にシンプルな組織に限定する。他方、従来技術であって監督を必要としないでうまくおこなう手法は、一つの器官に特定され、且つそれらが処理可能な解剖学的な組織の範囲に限定される。
【０００７】
監督される技法の例は、下記の参考文献に教示されている。
【０００８】
ダブリュー．イー．ヒギンス他（Ｗ．Ｅ．Ｈｉｇｇｉｎｓ  ｅｔ  ａｌ．）“Ｅｘｔｒａｃｔｉｏｎ  ｏｆ  ｌｅｆｔ−ｖｅｎｔｒｉｃｕｌａｒ  ｃｈａｍｂｅｒ  ｆｒｏｍ  ３Ｄ  ＣＴ  ｉｍａｇｅｓ  ｏｆ  ｔｈｅ  ｈｅａｒｔ，”Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ｍｅｄｉｃａｌ  Ｉｍａｇｉｎｇ，Ｖｏｌ．９，ｎｏ．４，ｐｐ．３８４−３９４，１９９０；  イー．エイ．アシュトン他（Ｅ．Ａ．Ａｓｈｔｏｎ  ｅｔ  ａｌ）“Ａ  ｎｏｖｅｌ  ｖｏｌｕｍｅｔｒｉｃ  ｆｅａｔｕｒｅ  ｅｘｔｒａｃｔｉｏｎ  ｔｅｃｈｎｉｑｕｅ  ｗｉｔｈ  ａｐｐｌｉｃａｔｉｏｎｓ  ｔｏ  ＭＲＩ  ｉｍａｇｅｓ，”ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ｍｅｄｉｃａｌ  Ｉｍａｇｉｎｇ，Ｖｏｌ．１６，ｐｐ．３６５−３７１，１９９７；ａｎｄ
エム．ダブリュー．ハンセン他（Ｍ．Ｗ．Ｈａｎｓｅｎ  ｅｔ  ａｌ，）“Ｒｅｌａｘａｔｉｏｎ  ｍｅｔｈｏｄｓ  ｆｏｒ  ｓｕｐｅｒｖｉｓｅｄ  ｉｍａｇｅ  ｓｅｇｍｅｎｔａｔｉｏｎ，”ＩＥＥＥ  Ｔｒａｎｓ．Ｐａｔｔ．Ａｎａｌ．Ｍａｒｃｈ．Ｉｎｔｅｌ．，Ｖｏｌ．１９，ｐｐ．９４９−９６２，１９９７．
監督されない技法の例は、下記の参考文献に教示されている。
【０００９】
ダブリュー．エム．ウェルズＩＩＩ他（Ｗ．Ｍ．Ｗｅｌｌｓ  ＩＩＩ  ｅｔ  ａｌ，）“Ａｄａｐｔｉｖｅ  ｓｅｇｍｅｎｔａｔｉｏｎ  ｏｆ  ＭＲ  ｄａｔａ，”ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ｍｅｄｉｃａｌ  Ｉｍａｇｉｎｇ，ｐｐ．４２９−４４０，１９９６；ａｎｄ
ジェイ．ラジャパルス他（Ｊ．Ｒａｊａｐａｌｓｅ  ｅｔ  ａｌ，）“Ｓｔａｔｉｓｔｉｃａｌ  ａｐｐｒｏａｃｈ  ｔｏ  ｓｅｇｍｅｎｔａｔｉｏｎ  ｏｆ  ｓｉｎｇｌｅ−ｃｈａｎｎｅｌ  ｃｅｒｅｂｒａｌ  ＭＲ  ｉｍａｇｅｓ，”ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ｍｅｄｉｃａｌ  Ｉｍａｇｉｎｇ，Ｖｏｌ．１６，ｎｏ．２，ｐｐ．１７６−１８６，１９９７．
発明者等は、タメズ−ペナ他（Ｔａｍｅｚ−Ｐｅｎａ  ｅｔ  ａｌ，）“Ａｕｔｏｍａｔｉｃ  Ｓｔａｔｉｓｔｉｃａｌ  Ｓｅｇｍｅｎｔａｔｉｏｎ  ｏｆ  Ｍｅｄｉｃａｌ  Ｖｏｌｕｍｅｔｒｉｃ  Ｉｍａｇｅｓ，”ＩＥＥＥ  Ｃｏｍｐｕｔｅｒ  Ｖｉｓｉｏｎ  ａｎｄ  Ｐａｔｔｅｒｎ  Ｒｅｃｏｇｎｉｔｉｏｎ  ９８中の関連する概念を提示する。上記全参考文献は全内容が、本開示中に参照によって取り込まれる。
【００１０】
（発明の要約及び目的）
本発明の目的は、身体組織等の４Ｄ（空間と時間）運動の再現化（ｋｉｎｅｔｉｃ  ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ）と視覚化のシステムと方法を提供することである。
【００１１】
本発明の別の目的は、運動の経過（ｍｏｔｉｏｎ  ｔｒａｃｋｉｎｇ）を組み入れた再現化と視覚化のシステムと方法を提供することである。
【００１２】
本発明の更に別の目的は、同時に無限の組織セグメント数を与える再現化と視覚化のシステムと方法を提供することである。
【００１３】
本発明の更に別の目的は、再現され視覚化された組織の生物力学的特徴を評価（ａｃｃｏｕｎｔ）に取り入れる再現化と視覚化のためのシステムと方法を提供することである。
【００１４】
本発明の更に別の目的は、そのような生物力学的特徴のデータベースを組み入れる再現化と視覚化のためのシステムと方法を提供することである。
【００１５】
本発明の他の目的は、有限なエレメント生物力学モデルを作成する再現化と視覚化のシステムと方法を提供することである。
【００１６】
本発明の更に他の目的は、バーチャル空間で動くことができ、身体の各部が与えられた力に対しバーチャル空間でリアルな動作で反応する、「ディジタルクローン」もしくは人体のバーチャルリアリティレプリカを作製する再現化と視覚化のシステムと方法を提供することである。
【００１７】
これら及びその他の目的を達成するため、本発明は次に挙げる手法を実施するシステムと方法に向けられたものである。画像シーケンス中のある画像、典型的には最初のものが、その画像中に示された組織と構造を特定するために、統計的手法によってセグメント化される。それから、すべての画像に対し同じことを行う代わりに、各組織と構造の運動を推算（ｅｓｔｉｍａｔｅ）することによってセグメント化を持ち越す。この推算（ｅｓｔｉｍａｔｉｏｎ）は、既知の組織と構造との特徴を評価（ａｃｃｏｕｎｔ）に加える、なおこの既知の特徴はバネ弾性係数で表現されるものである。
【００１８】
このセグメント化は、次のようなやり方で行うことができる。まず、粗画像データを得るために、ＭＲＩによる場合と同様に、３Ｄ画像の時間シーケンスを作成する。その最初の画像を、複数の領域に分割化（ｄｉｖｉｄｅｄ）もしくはセグメント化（ｓｅｇｍｅｎｔｅｄ）する。なお、この各領域は画像化された身体部分の一組織もしくは一構造に対応する。セグメント化のためのこの統計的手法は、各ボクセル（３次元ピクセル）中の画像に対するの局所的平均値（ａ  ｌｏｃａｌ  ｍｅａｎ  ｖａｌｕｅ）と局所的分散（ａ  ｌｏｃａｌ  ｖａｒｉａｎｃｅ）の推算を行うこと（ｅｓｔｉｍａｔｉｎｇ）を含んでいる。推算された局所的平均値と推定された局所的分散が十分近似していれば、このボクセルは領域に加えられる。こうして、各領域は、様々な組織或いは構造の応じて境界が決められる。
【００１９】
一旦、最初の画像についてのセグメント化が終了すると、同じ方法ですべての画像のセグメント化を行う場合に使用されるであろうコンピュータの電力を節約するため、他の画像に対する処理を持ち越す。この持ち越しは、種々の部品の既知の運動特性を使用し、それらをバネ定数として表現する。例えば、腱と骨の接続の特性はその二つを結ぶバネをイメージしそのバネ定数を与えることで表現できる。これらの既知の運動特性は、各組織或いは各構造が動く場合の速度（速さと方向）の推算（ｅｓｔｉｍａｔｅ）に使用できる。推算された速度は、次の画像において各組織或いは各構造がどこに存在するはずであるかを決定するために使用することができ、こうすることで、このシーケンス中のある画像からのセグメンテーションが次の段階に持ち越される。
【００２０】
このセグメーテーションは、正確なものというより統計的なものであるから、どの二つのボクセルを接続すべきかを決定するために、多様な統計的テストを使用することが望ましい。間違った接続は、間違った接続不足より大きなダメージを最終的な視覚化の品質に与えるので、この手法の実行は、接続を形成する前の統計的なテストが合致することを必要とすることができる。
【００２１】
また、この統計的テストは、ある境界値もしくはカットオフ値の設定を必要とする。勿論、カットオフ値のまずい選択は、すべての視覚化技術を混乱させることができる。そのような結果を避けるため、本発明の一態様は、カットオフ値を最適化するための反復手法（ｉｔｅｒａｔｉｖｅ  ｔｅｃｈｎｉｑｕｅ）を与えている。
【００２２】
これらのステップの最後に、骨、筋肉、肌及び他の組織の詳細な４次元の画像がバーチャルリアリティで「ディジタルクローン」として利用できる。この「ディジタルクローン」は、運動と生物力学的特徴を包含している。このバーチャルリアリティ画像は臨床医による、例えば特定の組織を検査したりすることによって、或いは力を与えたりすることによって相互に影響される。
【００２３】
（良好な実施の形態）
本実施形態のプロセッサは、シリコングラフィックスインディゴ（Ｓｉｌｉｃｏｎ  Ｇｒａｐｈｉｃｓ  Ｉｎｄｉｇｏ）ワークステーションで実施され、ビジュアル化には、Ｉｒｉｘ  Ｅｘｐｌｏｒｅｒ  ３Ｄグラフィックスエンジンを用いたが、当業者にとっては、本明細書の開示に基いて本発明を任意のハードウェア及びソフトウェアに容易に適用し得るものである。
【００２４】
ＭＲＩ等で形成された一連の画像は、生画像データとして受信され、各画像は一般にグレイスケールビットマップの形式にある。生画像データを使用する場合は、各画像をセグメントに分割し、骨、それを取り巻く筋肉などどいった画像中の各部位（フィーチャー）を個別に取り出す必要がある。連続画像の場合、各画像がそれぞれ個別にセグメント化されるわけであるが、発明者等は、まず第１の３Ｄイメージをセグメント化して種々の部位の運動をトラックし、この第１の画像のセグメントから外挿補間することによって、残りの画像をセグメント化することが計算上非常に効率的であることを発見した。
【００２５】
より具体的には、被写体自体は、連続空間における点ｘｃ＝（ｘｃ，ｙｃ，ｚｃ）に存在するが、そのイメージは離散的空間の点ｘ＝（ｘ，ｙ，ｚ）として表現される。ここで、各離散点ｘはボクセル（ｖｏｘｅｌ）と呼ばれる。本発明では、通常はグレイスケールイメージであるｔｅｍｐｏｒａｌな立体測定ディジタル画像ｇ（ｘ，ｙ，ｚ）を、連続３Ｄ時間イメージｇｃ（ｘｃ，ｙｃ，ｚｃ，ｔ）のサンプリングバージョンとして表現する。連続画像は、時間的に不変な照射フィールドＩ（ｘ）で、３次元の動的な複合物体ｆ（ｘｃ，ｔ）を照射することによって形成され得る。
【００２６】
各領域Ｒｉは、ランダムなフィールドＺ＝｛Ｚ（ｘ，ｙ，ｚ），（ｘ，ｙ，ｚ）∈Ｒｉ｝として表現し得ると考えられる。即ち、各領域は、その領域全体の確率過程であり、ランダムフィールドＺ（ｘ，ｙ，ｚ）は、領域Ｒにおける任意の点のボクセルＸ＝（ｘ，ｙ，ｚ）に関するランダム変数である。
【００２７】
物体ｆは、Ｍ個の異なる領域Ｒ１、Ｒ２、…、ＲＭからなるものとする。これらの領域が別個のものとして記述される場合、任意のｉ≠ｊに対してＲｉ∩Ｒｊ＝φであることを意味する。連続画像ｇｃが独立加法ウィーナー過程ηによって転化される場合、時間ｔ＝ｋＴでのボクセルｘ∈ｇにおける離散画像ｇは、
ｇ（ｘ，ｔ）＝Ｉ（ｘ＋Δｘ，ｙ＋Δｙ，ｚ＋Δｚ）ｆ（ｘ＋Δｘ，ｙ＋Δｙ，ｚ＋Δｚ，ｋＴ）
＋η（ｘ＋Δｘ，ｙ＋Δｙ，ｚ＋Δｚ，ｋＴ）
で与えられる。ここで、（Δｘ，Δｙ，Δｚ）は、サンプリング空間間隔である。
【００２８】
本発明の目的は、物体ｆを構成するＭ個の領域の運動を分析することである。これを行なうために、上述したようにウィーナー過程により転化された測定画像ｇから、これらの領域をセグメント化し、トラッキングする。即ち、立体（３Ｄ）のセグメント化、運動の推算、トラッキングを階層的に実行する。最初に、３Ｄセグメント化について説明する。
【００２９】
図１のフローチャートは、セグメント化のプロセスを示す。図面において、同様の構成要素には同じ符号を付してある。このプロセスの詳細な工程は以下の通りである。
【００３０】
まず、ステップ１０１で、ＭＲＩにより連続する画像を撮像する。これにより、生画像データが得られる。次に、ステップ１０３で、一連の画像のうち、第１の画像の生データを、先に例示したようなプロセッサに入力する。ステップ１０５で、各ボクセルにおいて、画像データの局所平均値と領域変数を求める。ステップ１０７で、ステップ１０５で求めた各画像データの平均値と変数を比較することによって、ボクセル間の連続性を求める。連続性を求めたならば、ステップ１０９で、スプリットを必要とする領域を決定し、これらの領域をスプリットする。ステップ１１１でセグメントを緩和することによって、スプリットされた領域の精度を向上させる。ステップ１１３で、どの領域をマージすべきかを決定し、決定した領域をマージする。ステップ１１５で、再度セグメント緩和を行なう。このようにして、生画像データをセグメント化された画像に変換し、ステップ１１７で処理を終了する。
【００３１】
ステップ１０５の局所平均値と変数の計算、及びステップ１０７の連続性の推算は、立体の各ボクセルをグラフ中のノードとみなす結合領域発展法（ｇｒｏｗｉｎｇ  ｓｃｈｅｍｅ）に基いて行なわれる。まず、３次元空間における隣接する系（ｓｙｓｔｅｍ）とクリーク（ｃｌｉｑｕｅ）を定義する。Ｇｉ（ｘ，ｙ，ｚ）＝Ｇｘ，ｙ，ｚ、（ｘ，ｙ，ｚ）∈Ｒｉ１を、領域Ｒｉ（ｘｋ）
【外１】
0004360757.tif 000002
を満たす領域Ｒｉの部分集合の集まりであることを意味する。したがって、Ｇｉ（ｘ）はｘの隣接点の集合である。Ｃｉ（ｘ）における各ボクセルのペアが隣接するならば、部分集合Ｃｉ（ｘ）⊆Ｒｉはクリークである。領域Ｒｉの境界或いはその近傍のボクセルの隣接点は、領域内部のボクセルが有する隣接点よりも少ない。
【００３２】
特性が十分に近似する隣接ボクセルは、弧で結ばれる。各画像セグメントＲｉは、結合された同一の部分に所属するすべてのボクセルの最大集合である。画像が、ガウス構造（ｔｅｘｔｕｒｅｄ）の領域Ｚｉで構成されると仮定すると、この画像中の各ボクセルに、局所平均値μｉ（ｘ）と、局所領域変数σｉ２（ｘ）という２つの特性パラメータが与えられる。即ち、Ｚｉ＝Ｎ（μｉ，σｉ）で表わされる。２つのボクセルの平均値及び変数のそれぞれが互いに十分近似するならば、これら２つのボクセルは弧で結合される。比較されるボクセルについて２つの特性値がわかれば、２つのボクセルの近似性は、ガウス累積分布関数にしたがってテストされる。一方、ある一定の確さでしか特性値を求めることができない場合は、比較統計テストを行なう。
【００３３】
ステップＳ１０５において、ノイズがある場合の局所平均と変数の求め方は、次の通りである。標準画像復元方法は、領域境界を適正に取り扱うことができない。境界問題は、推算されるパラメータのギブスのランダムフィールドを仮定することによって取り扱われる。即ち、所定のボクセルで求められる値は、結合されたボクセルの画像値の関数になる。ボクセルｘと、これに関連して結合される
【外２】
0004360757.tif 000003
【数１】
0004360757.tif 000004
ここで、
【数２】
0004360757.tif 000005
【外３】
0004360757.tif 000006
階ではまだわかっていない。この問題を取り扱うには、３次元空間の単純なラインセグメントの連続性に関して、一連の前提条件を仮定する必要がある。まず、３次元空間のラインセグメントは、同じ隣接ボクセルを共有する一対の構成クリークとして定義される。したがって、各ボクセルにおいて、複数のラインセグメントがあり、各ラインセグメントは、互いには隣接しないが前記共有するボクセルとは隣接する２つの隣接ボクセルによって形成される。例えば、図３〜５に示す位置関係にあるボクセルを隣接ボクセルとし、図２に示す隣接の系を考えるなら、１９のラインセグメントがあり、その中で直交平面内にある１２のラインセグメントが図５Ａ〜５Ｌに示される。これらのラインセグメントは、データの局所分散を計算することによって求められる。
【００３４】
【数３】
0004360757.tif 000007
ここで、Ｍは図２〜５の隣接点の数である。したがって、ラインセグメント検出は、
【数４】
0004360757.tif 000008
で定義される。ここで、ｐ（ｘ）はボクセルの特性であり、このボクセルにおける平均値又は変数である。Ｔ０はしきい値である。Ｌｐ＝｛ｘ，ｘｉ，ｘｊ｝＝１ならば、点｛ｘ，ｘｉ，ｘｊ｝は同じラインセグメントに属する。この手順をラインセグメント検出と呼ぶ。
【００３５】
求めたボクセルの平均値と変数は、ほとんどすべての医療用画像形成において、ラインセグメント検出におけるボクセル特性として使用することができる。同じラインセグメントに属するボクセルは結合され、パラメータの算出に使用される。このようなラインセグメント検出によって、連続性を求めることができる。求めた連続性情報を用いて、点ｘｉに連結されるボクセルの重みａｎｉｊを強調する。
【００３６】
【数５】
0004360757.tif 000009
ここで、ａｎはステップｎでの平均値と変数パラメータの初期重みである。
【００３７】
【数６】
0004360757.tif 000010
重みＷμ、Ｗσは自明ではなく、求めたパラメータの振るまいを反映したものでなければならない。平均値及び変数の算出においては、局所エネルギーを最小にする重みを選択することが好ましく、例えば次式で表わされる重みが考えられる。
【００３８】
【数７】
0004360757.tif 000011
このように、重みゲインは局所分散と関連する。なお、ｋｉは次式で与えられる正規化定数である。
【００３９】
【数８】
0004360757.tif 000012
ここで、Ｎは隣接ボクセルにおける点の数である。したがって、点ｘｉとｘｊが各パラメータにおけるラインセグメント検出により結合されるべき点であると判断されたならば、重みａｎｉｊは増大する。ライン構造がある場合は、ラインセグメント検出はエンドライン検出を含むことになる。エンドラインの検出は、ライン検出に次式で表わされる条件を追加し、この条件を、重みａｎｉｊ、ａｎｉｉに関連付けて用いることによって行なわれる。
【００４０】
【数９】
0004360757.tif 000013
丁度開示されたアプローチは繰り返し行われるアプローチであり、平均及び分散についての概算値が反復的に改良される。それ故にそのアプローチは収束の基準を満たすまで数回繰り返されるべきである。経験及び数回の試行錯誤は２０回の繰り返しが平均及び分散についての信頼しうる概算値を得られるのに十分であることを示している。それ故、丁度開示されたアプローチは確率的緩和技術（ｓｔｏｃｈａｓｔｉｃ  ｒｅｌａｘａｔｉｏｎ  ｔｅｃｈｎｉｑｕｅｓ）よりも速い。
【００４１】
ボクセル・コネクティビティ（ｖｏｘｅｌ  ｃｏｎｎｅｃｔｉｖｉｔｙ）の概算には分散の基準を用いる。概算された複数の分散が類似するか否かを決定するためにσｊ２（ｘ）のσｉ２（ｘ）に対する比が計算される。
【００４２】
【数１０】
0004360757.tif 000014
（ｓｑｕａｒｅ  ｆｉｅｌｄ  ｉｎｔｅｎｓｉｔｙ  ｆｕｎｃｔｉｏｎ）がノイズ分散よりも非常に大きいかＩ２（ｘ）＞＞ση２であれば、この比はｆ分布（ｆ  ｄｉｓｔｒｉｂｕｔｉｏｎ）に従い、ノイズηから独立であると見なすことができる。もしその考慮のもとで複数のクリーク（ｃｌｉｑｕｅ）が同様の大きさであれば、Ｆｊが与えられたしきい値ｔｈ≧１よりも小さいことを条件としてσｊ２（ｘ）とσｉ２（ｘ）は同じ分散の概算であるということができる。他の言葉で言い換えると、もしＦｊ＜ｔｈならば、ボクセルｘ∈Ｒｉと隣接するボクセル（ａ，ｂ，ｃ）∈Ｃｊ⊂Ｒｉとが結合する弧が存在する。
【００４３】
コネクティビティに対する他の独立したテストを行うことは有用である。一つのそのような独立したテストはボクセル間のマハラノビス（Ｍａｈａｌａｎｏｂｉｓ）距離を基礎とする。マハラノビス距離は
【数１１】
0004360757.tif 000015
として定義される。場の強度が一定でない場合において、小さな移動距離Δｘに対して、場の強度はＩ（ｘ＋Δｘ）＝（１−αｘ）Ｉ（ｘ）となり、もし
【数１２】
0004360757.tif 000016
ならば、しきい値距離ｄｍｉｎを規定して近接したボクセルが連続しているということができる。
【００４４】
そのアプローチはイルミネーション（ｉｌｌｕｍｉｎａｔｉｏｎ）分散に対して確固としたものであり、ボクセル・コネクティビティの計算が画像の１回のスキャンで行うことが可能であるため計算的に有効である。すべての生じ得る弧に対して、２つのコネクティビティの概算値が存在する。間違ったコネクションを形成するコストは正しいコネクションを形成しない場合のコストに比べて非常に大きいため、コネクションは双方の概算値が一致しない限り形成されない。
【００４５】
一旦コネクティビティのグラフが作られると、それぞれの領域は同じコネクテッド・コンポーネンツ（ｃｏｎｎｅｃｔｅｄ  ｃｏｍｐｏｎｅｎｔｓ）に属したすべてのボクセルのセットとなる。しかしながら、雑音と、貧弱な平均及び分散の概算に起因して領域間でいくらかの漏れが存在しうる。多様な技術によってそのような漏れを避けることが可能である。あるそのような技術において、すべての孤立した結合弧を移動させることによって緩やかに結合された領域を分割するために、単純３次元形態クロージングオペレーション（ｓｉｍｐｌｅ  ３Ｄ  ｍｏｒｐｈｏｌｏｇｉｃａｌ  ｃｌｏｓｉｎｇ  ｏｐｅｒａｔｉｏｎ）をコネクティビティのグラフに適用することが可能である。他のそのような技術は２つの異なるしきい値を選んでヒステリシス（ｈｙｓｔｅｒｅｓｉｓ）を使用する。最終的なコネクティビティのマップ（ｍａｐ）はより低いしきい値から、隣接したコンポーネントを加えた高いしきい値のコネクティビティにより形成される。後者の技術はより計算的に有効なものとして好まれている。
【００４６】
勿論、実像は、小さな連続的線構造が常に領域の間に存在するという上で説明した基礎的な画像モデルには従わない。従って、丁度開示されたアプローチの結果である区分は立体的画像の分類の誤りから自由ということではない。更に、しきい値ｔｈの不正確な選択は画像の不正確な区分の原因となる。大きなしきい値はいくつかの領域を１つにまとめる原因となり、一方小さなしきい値は多種多様の破損した領域の原因となる。もし１つの領域がいくつかの隣接した断片へと破壊されたならば、断片を１つにまとめることによりその問題に取り組むことが可能である。画像モデルを維持して、平均が等しい隣接するイメージのの断片は同じ領域に属するものということができる。
【００４７】
ステップ１０９の領域分割は今詳細に述べられる。ヒステリシスの後でさえも、類似した複数の構造が不適当に１つにまとめられ得る。一つの抽出された領域が２つかそれ以上の独立した構造から構成されているのか否かをテストする単純な方法はその領域を平滑な表面へとフィッティングした際のエラーを計算し、その後分散イメージを用いて抽出された領域の平均分散とエラーを比較することである。もしそのエラーが分散と非常に異なっていたならば抽出された領域は単一の構造ではなく、分割されなくてはならない。
【００４８】
領域Ｒｉ対する観測データにフィットする最も平滑な表面ｕ（ｘ）は２次元エラー（ｓｑｕａｒｅ  ｅｒｒｏｒ）を最小化したものである。
【００４９】
【数１３】
0004360757.tif 000017
ｇは観測されたイメージである。もし観測されたイメージが場の分解（ｄｅｇｒａｄａｔｉｏｎｓ）による損傷を受けるならばｕはそのような分解に合わせられるべきである。ＭＲＩのケースにおいてコイル場の分解（ｃｏｉｌ  ｆｉｅｌｄ  ｄｅｇｒａｄａｔｉｏｎｓ）に対する演繹的な知識から、それは
Ｉ（ｘ）＝ｅｘｐ［−（ａ０＋ａ１ｘ＋ａ２ｙ＋ａ３ｚ＋ａ４ｘｙ＋ａ５ｘｚ＋ａ６ｙｚ＋ａ７ｘ２＋ａ８ｙ２＋ａ９ｚ２］
のように仮定される。このケースでは同じ表現でｕｉ（ｘ）のよいモデルが与えられ、ａ１からａ９までの係数は対数イメージであるｌｏｇ  ｇ（ｘ）へのスタンダード・リースト・スクエア・アプローチ（ｓｔａｎｄａｒｄ  ｌｅａｓｔ−ｓｑｕａｒｅ  ａｐｐｒｏａｃｈ）により発見することができる。一旦最も平滑な表面が抽出された領域であることが明らかになり、もしフィッティングのエラーが平均標準偏差と類似するならばその領域に属さないすべての点は移動させることができる。
【００５０】
【数１４】
0004360757.tif 000018
σ（ｘ）は上で導かれた、概算された標準偏差である。その値は表面での点の距離ｄｈｐｌａｎｅと比較される。ｄｈｐｌａｎｅは次のように与えられる。
【００５１】
【数１５】
0004360757.tif 000019
もし与えられたボクセルｘに対する距離が与えられたしきい値ｔｐσｉよりも小さかったならば、そのボクセルはその領域に属すると決定できる。
【００５２】
他方では、上記から引き出されたＭＳＥフィッティングエラーがσｉよりかなり大きい場合、その領域は、おそらく二つ或いはそれ以上の構成要素からなり、分割されなければならない。二つの構成要素のみが合併されることが想定されれば、それらを分割するための一つの簡単な方法は、次の反復操作を用いることである。
【００５３】
（１）ＭＳＥ基準を用いるｕ（ｘ）を見積もる
（２）記述された距離基準（距離としきい値の比較）を使用するポイントメンバーシップのテストをする
（３）領域外のポイントを使用する新しいｕ（ｘ）を見積もる
（４）エラーは減少するか？そうでなければ、（２）へ戻る
その簡単な反復処理は、その領域内の主な構成要素を抽出する。その領域の残りは、新しい領域としてマークされる。そのアプローチは、既知のｋ−方法アルゴリズムの簡単な概念であり、ここではｋ＝２とし、最適なクラス重心の代わりに最適なクラス表面を求める。
【００５４】
その段階で、フィールド・デグラデーションを修正することは可能である。上述したＩ（ｘ）の式が与えられると、ａ１からａ９までのパラメータは、主な抽出された領域の表面関数以上のそれらのパラメータの平均値であり、ａ０は、回復できないと想定され、任意値に設定される。そのようなパラメータがＩ（ｘ）の距離偏差に寄与しないが代わりにべき指数の前の係数にのみ寄与することがリコールされるとき、そのようなａ０の処理が意味をなす。画像モデルからガウス雑音（ｘ）が十分に小さい場合、ｇ（ｘ）＝Ｉ（ｘ）・（ｘ）、或いは言い換えれば、ｆ（ｘ）＝ｇ（ｘ）／Ｉ（ｘ）となる。
【００５５】
ステップ１１１のセグメンテーション緩和に基づくギブス・ランダム・フィールドは、後の確率を最小化することによってセグメンテーションエラーを修正する目的にかなう。その最後に、標識画像のマルコフ・ランダム・フィールド・モデルが呈され、ボクセルラベル確率が隣接したラベルの関数である。ギブス・ランダム・フィールドは、ボクセルの確率の判断で想定され、隣接のラベルの分布が与えられる。後の確率は、画像回復で首尾よく使用されている反復暫定モードアプローチを使用し、最適化され得る。
【００５６】
Ｌ（ｘ）はラベル化された画像で、ｆ１（ｘ）はフィールド修正画像の主な構成要素である場合、その画像のギブス・ランダム・フィールド・モデルは後の確率分布関数を与える。
【００５７】
【数１６】
0004360757.tif 000020
ここで、Ｖｃ（Ｌ）は、容積測定の画像Ｇｘ，ｘ∈Ｌにおける隣接したシステムに関連したクリーク・ポテンシャルであり、ｙｉとｓｉはその平均値と分散であり、そしてγは規格化定数である。Ｖｃは、任意の費用関数ｆやｈの点から言えば、表現され得る。
【００５８】
【数１７】
0004360757.tif 000021
一般的に、費用関数は、共に一定であるが、最良の実施において、ｆ（ｘ）＝αｓｉ／（ｓｉ＋｜ｘ｜２）とｈ（ｘ）＝βは、中央ボクセルからのグレイ値偏差を不利にするために使用される。その不利益の程度は、αとβを設定することによって、決めることができる。
【００５９】
その等式の全体的な最小化は可能であるが、時間を費やす。代わりに、反復暫定モード（ＩＣＭ）が用いられる。全段階でガウスのＰＤＦ（確率密度関数）とその確率を見積もるための以前の解決策を使用することによって、その技術は繰り返して、その関数を最小化する。そのような処理は、与えられた適した最初の条件以外で、全体的な最小化を保証しないが、平均最適なローカル最小値に収束する。
【００６０】
ステップ１１３でマージしている統計領域はここで描かれている。簡単な統計的なｔ−テストは、与えられた信頼度αで、二つの隣接した画像パッチが等しい意味を持つかどうか決定するために適用されることができる。フィールド修正画像からの平均値と分散（ｙｉ、ｓｉ）を持つパッチＰｉと平均値と分散（ｙｉ、ｓｉ）を持つ隣接するパッチＰｊが与えられると、もし、
【数１８】
0004360757.tif 000022
がｔα，νよりも大きく、ここでνは、自由度の数であり、
【数１９】
0004360757.tif 000023
によって与えられ、ｎ値は、ｙ値とｓ値の推算に含まれるポイントの数である場合、その二つのパッチは１−αの信頼度内で同じ平均値を持つ。それらの計算はその分野において知られるているｔ−テストに基づいている。
【００６１】
もし、そのイメージが異質でないならば、単純な基準はほとんどの誤った分類のセグメントを同質の領域に統合することが可能である。イメージの異質はパッチの各ボクセルに対するローカル平均と分散推定値のみの使用をユーザに強いる。したがって、パッチｉとｊとの間の境界での全ボクセルのローカル平均と分散は、両方の領域が等しいかどうであれ、ａｂｏｖｅ  ｔｏ  ｔｅｓｔ  ｉ０に対する方程式で用いられる。平均と分散に対して連結グラフからの推定値を用いることができるが、そのような推定値は次の穴だらけの推定値より信頼性が低い。
【００６２】
【数２０】
0004360757.tif 000024
ここで、係数ｗは、ウィンドウ形状を決定し満たす。
【００６３】
【数２１】
0004360757.tif 000025
ガウシアンウィンドウ（Ｇａｕｓｓｉａｎ  ｗｉｎｄｏｗ）が用いられ、境界ボクセルｘｉ∈Ｐｉで集中される。
【００６４】
一方、平均と分散が正規分布からもたらされると仮定すると、チャーノフ（Ｃｈｅｒｎｏｆｆ）距離が用いられる。
【００６５】
【数２２】
0004360757.tif 000026
ここで、ΣｉとΣｊは共分散行列であり、ｓは正規分布に対するベイズ誤差を最小にする定数であり、
【数２３】
0004360757.tif 000027
チャーノフ（Ｃｈｅｒｎｏｆｆ）距離が最大になる点によって与えられる。もし、最適距離が関係しないならば、ｓ＝１／２とセットすることが正規分布を単純化し、バタチャリア（Ｂｈａｔｔａｃｈａｒｙｙａ）距離を与える。
【００６６】
【数２４】
0004360757.tif 000028
その距離は計算することが容易であり、二つの分布の可分性の良い尺度となる。その距離を用いることによって、ユーザは二つのパッチの点の分布が類似しているか否かを決定づける閾値を決めることができる。
【００６７】
一旦、ステップ１１３の統計的領域結合が実行されると、ステップ１１５で再びセグメンテーション緩和が実行される。その処理はセグメント化されたイメージのステップ１１７で終了する。
【００６８】
一旦、イメージが領域に適当に分割されると、運動の推定と追跡は一連のＭＲＩスキャンに引き受けられる。そのような処理は、図６のフローチャートに示されており、下記に詳細について説明する。もし、イメージが接合部についてであるならば、腱、靱帯、筋肉のような柔らかい組織の堅くない運動が含まれるので運動は複雑となる。そのような堅くない運動は、自動化４Ｄセグメンテーションへの応用を困難にする位相変化をもたらす。
【００６９】
一連の次の３Ｄイメージのセグメンテーションが同様の方法でなされうるにもかかわらず、そのアプローチはある問題を引き起こす。第１に、低コントラストな構造間の境界がはっきり定義されず、それ故、ノイズとサンプリングが一連のイメージ間を変化するセグメンテーションによって形成される境界を引き起こす。第２に、多くの構造をとるため、一連の各イメージに対する一貫した標識を保持することが難しい。第３に、セグメンテーションアルゴリズムは計算的にきつく、一連の各イメージに対してセグメンテーションアルゴリズムを再実行することは好ましくない。それゆえ、代替的なアプローチが一連の残ったイメージのセグメンテーションに対して開発されてきた。
【００７０】
運動の追跡と推定のアルゴリズムは一つのフレームから一連の最初のイメージまでセグメント化されたイメージを通過するのに必要とされる情報を与え、上述したようにそこから引き出された完全にセグメント化されたイメージは、ステップ６０１で入力される。堅い構成要素と堅くない構成要素の両方の存在は、３Ｄの運動の推定において観念的に考慮すべきである。本発明によれば、各ボクセルの運動ベクトルはイメージの選択された特徴点のレジストレーション後に推定される。
【００７１】
接合部に存在する多くの構造の運動を考慮するため、本発明のアプローチはイメージのセグメンテーションで見出された異なった構造の物質的特性である先行の知識を用いることにより柔らかい組織の局所変形を考慮する。そのような知識はステップ６０３で適当なデータベースフォームに入力される。また、種々の方法が堅い構造の運動と柔らかい組織の運動に適用されうる。一旦、選択された点が記憶されると、イメージの全ボクセルの運動ベクトルは、選択された点の運動ベクトルを補間することによって計算される。一旦、各ボクセルの運動ベクトルが推定されると、一連の次のイメージのセグメンテーションは、前のイメージのセグメンテーションの丁度伝達となる。そのような技術は、一連の全イメージが解析されるまで繰り返される。
【００７２】
有限要素モデル（ＦＥＭ）は、イメージの分析法と時間発展分析法として知られている。本発明は同様の方法にアプローチに従い、解剖学的構造の物理的特性のモデルとなる質量とバネの組み合わせの全エネルギーを最小化することによって点対応付けを回復する。本発明において、組み合わせはイメージ中の一つの構造によって強いられないが、その代わり全容量測定のイメージのモデルとはならない。ここで、位相特性は最初のセグメント化されたイメージによって供給され、物理的特性は先験的な特性と最初のセグメント化されたイメージによって供給される。運動の推定のアプローチは、一連の二つの連続したイメージの間のＦＥＭを基礎とした点対応付け回復アルゴリズムである。組み合わせの各ノードは、追跡されるイメージの自動的に選択された特徴点であり、バネ剛性は最初のセグメント化されたイメージとヒトの解剖学的構造の従来の知識と筋肉や骨のような典型的な生物の力学的構造特性から計算される。
【００７３】
多くの変形できるモデルは、スプリングで押さえられた質点を駆動するベクトル力場が、イメージから取り出されうると仮定する。ほとんどのそのようなモデルは、セミ自動的な特徴抽出アルゴリズムを組み立てるアプローチを使う。本発明は同様なアプローチを使用し、ｔ＝ｎでサンプリングされたイメージが３つの動的なスカラー場のセットであると仮定する：Φ（ｘ、ｔ）＝｛ｇｎ（ｘ），｜▽ｇｎ（ｘ）｜，▽２ｇｎ（ｘ）｝、即ち、白黒縮尺イメージの値、そのイメージの値のグラジエント（傾き）の絶対値、及びそのイメージの値のラプラシアン。従って、Φ（ｘ、ｔ）の変更は、スカラー場エネルギーＵΦ（ｘ）∝（△Φ（ｘ））２の二次の変更を起こす。更に、イメージの基になる構造体は、それらのスカラー場によって、平衡の状態において、スプリングで押さえられた質点のメッシュとして作成されると仮定される。平衡は、外力場があると仮定するけれども、力場の形は重要ではない。その質点の分配は、時間内に変わるように仮定されて、ｔ＝ｎ回の後の回周期△ｔの間の全体のエネルギー変更は、以下の式で与えられる。
【００７４】
【数２５】
0004360757.tif 000029
ここで、α、β、ａｎｄγがすべての固有場変更の分担のための重みであり、ηがひずみエネルギーの利得を重み付けし、Ｋが有限要素法固有のマトリクスであり、△Ｘが有限要素法節変位マトリクスである。その方程式の分析は、イメージ場又はメッシュポイント分配のどのような変更でも系統総計エネルギーを増大させることを示す。従って、ｇｎからｇｎ＋１へのポイントコレスポンデンスは、全体のエネルギー変化が最小であるメッシュ配置により与えられます。従って、ポイントコレスポンデンスは、以下の式により与えられる。
【００７５】
【数２６】
0004360757.tif 000030
ここで、
【数２７】
0004360757.tif 000031
その記法では、ｍｉｎｐｑは、ｑを最小化するｐの価値である。
【００７６】
上記の方程式が、多分、イメージのすべてのボクセルの運動（ポイントコレスポンデンス）を推算するために使えるとき、ボクセルの数と、一般にボクセルの数は１００万より多くなるが、及びこの式の複雑な性質とは、包括的に最小化することを困難にしている。問題を簡素化するために、目の粗い有限要素法網目は、ステップ６０５で、選ばれたポイントでイメージから組み立てられる。エネルギー最小化は、選ばれたポイントのポイントコレスポンデンスを与える。
【００７７】
そのようなポイントの選択は些細なことではない。第一に、実用的な目的のために、ポイントの数を非常に小さくする必要がある。一般に  ≒  １０４  ；選ばれたポイントが全体のイメージの運動を記述するように気をつけなければならない。第二に、境界追跡が正確な領域運動記述に十分なので、領域境界は重要な特徴です。第三に、領域境界において、グラジエントの絶対値は高く、ラプラシアンはゼロの交差ポイントにあり、領域境界を、追跡するための簡単な特徴とする。従って、分割された境界点は、有限要素法の作図において選ばれる。
【００７８】
境界点はイメージポイントの小さい部分集合を表しているけれども、より非常に多い境界点が実際上はある。ポイントの数を減らすために、ポイント注出ステップにおいて、境界点の強制された無作為抽出が使われる。強制は、すでに選ばれていたポイントと近すぎるポイントの選択を避けることからなっています。その強制は境界を横切るポイントのより均一な選択を許します。最後に、個々の領域の内部のポイントの運動推算の誤差の縮小のため、イメージの更にいくつかのポイントが、同じ距離の強制を使って無作為に選ばれる。実験的な結果は、５０００から１００００までのポイントが、２５６  ｘ  ２５６  ｘ  ３４個のボクセルの典型的な容積測定イメージの運動を推算し、説明するのに十分であることを示す。選ばれたポイントで、７５％は境界点として恣意的に選ばれる一方、残留２５％は内点である。勿論、適当な場所では、他のパーセンテージは使われえる。
【００７９】
一旦追跡するポイントの集合が選ばれたら、次のステップは、ステップ６０７でそれらのポイントのための有限要素法網目を組み立てることです。メッシュは、それぞれの領域の材料の属性と相互作用の属性を符号化することにより許された運動の種類を強制する。最初のステップは、すべての接点のために近隣の接点を見つけることである。当業者らは、近隣の接点を見つける操作が、メッシュのボロノイ（Ｖｏｒｏｎｏｉ）ダイアグラムを組むことと一致していると理解する。そのデュアル、デラウナイ（Ｄｅｌａｕｎａｙ）三角測量、は、与えられた節の配置のための最もよい可能な四面体の有限の素子を表す。ボロノイダイアグラムは拡張アプローチにより構成される。そのアプローチのもとで、離散的なボリュームの個々の接点が拡張される。そのような拡張は２つの目的を達成する。第一に、ある拡張されたポイントが別のものに接する時は、近隣のポイントが識別できるようにテストされる。第二に、すべてのボクセルはメッシュのポイントと関連できる。
【００８０】
一旦、すべてのポイントｘｉが近隣のポイントｘｊと関連すると、その２つのポイントは、バネ定数ｋｉｌｊｍを有するスプリングによって押さえられていると考えられます。ここで、ｌとｍにより材料を識別する。バネ定数は連結ポイントの物質的な相互作用の性質により定義されます；それらの物質的な相互作用の性質は、材料の知られている性質に従ってユーザによりあらかじめ決められる。もし連結ポイントが同じ領域に属しているならば、バネ定数はｋｉｌｊｍに減り、領域においてマテリアルの弾性特性から引き出される。もし連結ポイントが種々の領域に属しているならば、バネ定数は境界で材料の間の平均的な相互作用の力から引き出される。もしイメージされているオブジェクトが人の肩であるならば、バネ定数は以下のようなテーブルから引き出されうる：
0004360757.tif 000032
理論上では、相互作用は、任意の２つの調整領域間で定義されなければならない。しかし実際は、任意定数として残余を残して、映像中の主要なａｎａｔｏｍｉｃａｌ（解剖学上：構造上）の要素間のみの相互作用を定義し、許容される近似を行っている。このような近似において導入されるエラーは、上述したように実行される推測において導入される他のエラーと比較して、重要ではない。
【００８１】
バネ定数は、近似の大きさ及び通常用いられている（ａ  ｐｒｉｏｒｉ）骨のイメージ強度として、自動的に割り当てられることができる。通常使用される（ａ  ｐｒｉｏｒｉ）期待値（予測値）に整合する分割されたイメージ領域は、相対的に固定された骨のバネ定数に割り当てられる。柔らかい組織は、柔らかいバネ定数が相対的に割り当てられる。
【００８２】
一度メッシュを設定したら、ステップ６０９においてシーケンス中の次のイメージを入力し、ステップ６１１において、シーケンス中の２つの連続したイメージの間のエネルギーを最小限にする。エネルギーＵを最小限にする際の問題は、固定した運動から連想されるエネルギーと、変形可能（ｄｅｆｏｒｍａｂｌｅ）運動から連想されるエネルギーとからなる、２つの分離された問題に分割することができる。
【００８３】
剛体運動の見積もりは、メッシュの一様でないエネルギー（△ＸＴＫ△Ｘ）／２の剛体運動への寄与はゼロにとても近いという事実に従う。区分及び通常使用される解剖（分析）の知識は、堅いボディへ属する点を指し示す。このような点が、個々の固定された領域としてそれぞれ選択されたならば、剛体運動のエネルギーの最小限化は、剛体運動の回転Ｒｉ、及び領域の自己のエネルギーを最小限化する変換Ｔｉを、固定された領域Ｒｉ毎に、見つけることで成し遂げられる。
【００８４】
【数２８】
0004360757.tif 000033
【外４】
0004360757.tif 000034
点とって最適な置き換え母体（ｍａｔｒｉｘ）である。この最小限化の問題は、固定された領域毎の６段階の自由度のみを有する。：回転母体の中に３個あり、変換母体の中に３個ある。したがって、もし、シーケンス中の任意の２つのイメージの間の違いが十分小さいとき、６個の次元的に最も急激な降下技術（ａ  ｓｉｘ−ｄｉｍｅｎｓｉｏｎａｌ  ｓｔｅｅｐｅｓｔ−ｄｅｓｃｅｎｔ  ｔｅｃｈｎｉｑｕｅ）を介して、１２個の構成要素（９個の回転、３個の変換）を見つけることができる。
【００８５】
一度、剛体運動のパラメータが見つかったら、一様でない運動は、全システムのエネルギーＵの最小限化を通して見積もられる。この最小限化は、固定されたエネルギーの最小限化と同程度には単純化することができない。更に考慮が無い場合、変形可能な３Ｄ物体の中の自由度の数は、すべてのメッシュの中の接続点の数の３倍である。問題の本質は、メッシュ中の各接続点に対する単調減少技術（ａ  ｓｉｍｐｌｅ  ｇｒａｄｉｅｎｔ  ｄｅｓｃｅｎｔ  ｔｅｃｈｎｉｑｕｅ）の使用を許可する。ポテンシャルエネルギー及び運動エネルギーから、システムのラグラシアン（Ｌａｇｒａｎｇｉａｎ）（或いは、運動エネルギーからポテンシャルエネルギーを減じる物理学から定義される運動ポテンシャル）を、エネルギー場の傾斜に丁度合う駆動局所的力（ｔｈｅ  ｄｒｉｖｉｎｇ  ｌｏｃａｌ  ｆｏｒｃｅ）が存在するシステム中のすべての接続点について、ヨーラー−ラグランジェ方程式を引き出すために使用することができる。メッシュ中のすべての接続点にとって、局所的エネルギーは、以下の式により与えられる。
【００８６】
【数２９】
0004360757.tif 000035
ここで、Ｇｍは、ボロノイ（Ｖｏｒｏｎｏｉ）図の中の近傍を示す。
【００８７】
このように、すべての接続点において３個の自由度についての問題が存在する。３個の自由度の最小限化は、局所的な接続点のエネルギーを反復して減少させる単調減少技術を使用することで、実行される。接続点の局所的な減少方程式（Ｔｈｅ  ｌｏｃａｌ  ｎｏｄｅ  ｇｒａｄｉｅｎｔ  ｄｅｓｃｅｎｔ  ｅｑｕａｔｉｏｎ）は、以下の式の通りである。
【００８８】
【数３０】
0004360757.tif 000036
ここで、メッシュエネルギーの傾斜は分析的に計算可能である。場のエネルギーの傾斜は、２つの異なる分解により、イメージから数字的に見積もられる。ｘ（ｎ＋１）は次の接続点の位置であり、ｖは傾斜の寄与のための重み付け係数である。
【００８９】
最小限化のすべてのステップにおいて、接続点毎の過程は、隣接する接続点の先の配置を考慮に入れる。この過程は、エネルギー全体が局所的な最小限に達するまで繰り返され、小さな変形とって、全体的な最小限化に近づく或いは等しくなる。このように見つけられた配置ベクトルは、接続点における見積もられた運動を示す。
【００９０】
一度、丁度表わされた最小限化過程は、見本とされた配置場△Ｘ強度を産出する。配置場は、シーケンス中の１つのイメージから次への割り当て（ステップ６１１）をたどる必要がある密集した運動の場を見積もるために使用される。密集した運動は、メッシュ中のすべての隣接する接続点の寄与を重み付けすることにより、見積もられる。一定の速度モデルが、想定される。時間ｔに置けるボクセルｘの見積もられた速度は、ｖ（ｘ，ｔ）＝△ｘ（ｔ）／△ｔである。密集した運動の場は、以下の式から見積もられる。
【００９１】
【数３１】
0004360757.tif 000037
ここで、
【数３２】
0004360757.tif 000038
ｋｌ，ｍはバネ定数或いはボクセルｘ及びｘｊから連想される材料ｌとｍ間の堅さを示し、△ｔはシーケンス中の連続イメージの間の間隔を示し、｜ｘ−ｘｊ｜はボクセル間の単純なユークリジャン（Ｅｕｃｌｉｄｅａｎ）距離を示し、改ざんが、ボクセルｘに最も近い隣接する接続点を使用することで実行される。この改ざんは、材料の設定ｋｌ，ｍによりすべての隣接する接続点の寄与に重みをつける。；このように、見積もられたボクセルモーション（ｖｏｘｅｌ  ｍｏｔｉｏｎ）は、すべての同種の領域だけでなく、領域の境界においてさえも、類似している。
【００９２】
そして、ステップ６１５において、シーケンス中の次のイメージは、区分データにより満たされる。これは、１つのイメージの中で決定された領域が次のイメージの中に引き継がれることを意味する。このようにするために、速度は、次のイメージ中のすべてのボクセルについて見積もられる。これは、以下の式に示す見積もられた運動の逆マッピングにより成し遂げられる。
【００９３】
【数３３】
0004360757.tif 000039
ここで、Ｈは次のイメージ中の同じボクセル空間Ｓ（ｘ）へ落ちてゆく点の数である。このマッピングは時間ｔ＋△ｔにおけるすべての空間を満たすものではなく、調査された（ｍａｐｐｅｄ）隣接するボクセルの間の単純な改ざんは、この空間を膨らませる（ｆｉｌｌ  ｏｕｔ）ために使用することができる。一度、次のイメージ中のすべてのボクセルについて速度が見積もられたら、このイメージの割り当ては単純に、以下の式で表わされる。
【００９４】
【数３４】
0004360757.tif 000040
ここで、Ｌ（ｘ、ｔ）及びＬ（ｘ、ｔ＋△ｔ）は時間ｔ及びｔ＋△ｔのボクセルｘにおける割り当ての札である。
【００９５】
ステップ６１７において、このように展開された割り当ては、ステップ１１１及び１１５において行われるようにして、緩和札付け（ｒｅｌａｘａｔｉｏｎ  ｌａｂｅｌｉｎｇ）を通じて調整される。イメージ中のメッシュ接続点に微妙な調整が施される。そして、ステップ６１９において最後のイメージが割り当てられていない限り、次のイメージがステップ６０９において入力される。最後のイメージが割り当てられている場合、ステップ６２１において動作は終了する。
【００９６】
上記アプローチは、ｔｈ、ｄｍｉｎ、及びバタチャリア距離閾値の少なくとも３つのパラメータの入力を必要とする。セグメンテーション・オペレーションの変数は、初めの２つのパラメータの反復最適化によって得られる管理されない（ｕｎｓｕｐｅｒｖｉｓｅｄ）最適なセグメンテーションにおいて示される。
【００９７】
図７は、初めの２つのパラメータの費用関数を限定し最小にすることによって初めの２つのパラメータが最適化される、上述したような変数のフローチャートを示す。上述したように、ステップ７０１において、画像を入力する工程、ステップ７９３において、平均値と分散を概算し（見積もり）、ステップ７０５において、イニシャル・パラメータを入力し、ステップ７０７において、接続性を概算し、ステップ７０９において、領域を分割し、これらの工程は、上述した最適化は領域を最適化するため、最後に使用される領域を形成するのではなく、反復領域を形成する。
【００９８】
ステップ７１１において、パラメータｔｈ及びｄｍｉｎの値は、次のようにして調整される。  平均二乗誤差（ＭＳＥ）は、ＭＳＥ＝（１／Ｖ）∫（ｕ（ｘ）−μ（ｘ））２ｄｘによって与えられる。このとき、Ｖは、画像の体積（Ｖｏｌｕｍｅ）である。また、画像は、表面積Ｓを有する。各領域Ｒｉは、体積Ｖｉ及び表面積Ｓｉを有する。画像及び各領域の双方は、（表面積）１．５／体積からなるコンパクトネスを有する。２つのパラメータの値を最適化するためには、評価コスト関数（ｃｏｓｔ  ｆｕｎｃｔｉｏｎａｌ）は、２つのパラメータの中の特定の値のどれか２つを用いて費用を想定することを導入する。汎関数（ｆｕｎｃｔｉｏｎａｌ）は、関数の与えられたセットもしくはドメインにおけるすべての関数に、独自の実数（ｕｎｉｑｕｅ  ｒｅａｌ  ｎｕｍｂｅｒ）が割り当てられるというルールによって定義される。評価コスト関数は、
【数３５】
0004360757.tif 000041
によって与えらる。、Ｎ領域が存在し、ｗは重みファクタである。評価（汎）関数（Ｔｈｅ  ｅｖａｌｕａｔｉｏｎ  ｆｕｎｃｔｉｏｎａｌ）  Ｉ  は、ｔｈ及びｄｍｉｎの関数であり、従って、∂Ｉ／∂ｔｈ＝０＝∂Ｉ／∂ｄｍｉｎにセット（設定）することによって最小化することができる。これらの式は、最も勾配のある下降（ｄｅｓｃｅｎｔ）によって解くことができる。即ち、ｔｈ及びｄｍｉｎの値をステップ７１１における評価コスト関数の勾配（変化率）に近づくように調整すること、及びステップ７１３において評価コスト関数が減少しているかを確定することによって解くことができる。評価コスト関数が減少している限り、このプロセスは繰り返される。
【００９９】
マージンする領域について説明する。前述した説明において、ステップ７１５におけるバイアス又はイルミネーション・フィールド推算及び調整（訂正）、ステップ７１７における領域分割、ステップ７１９における主コンポーネント分析、ステップ７２１におけるギブズ・ランダム・フィールド・ベース・セグメンテーション・リラクゼーションについては既に触れている。バタチャリア閾値Ｔを最適化するのが可能であるとき、更に有能な計算技術はマージンを制約する。ステップ７２３に示す、ローカル  バタチャリア距離によって制約された統計領域マージンは、コンパクトネスが減少する場合、好ましくは、合体するべきであり、コンパクトネスが増加する場合、好ましくは、合体するべきではないという原則に基づく。このような原則に基づく理由は、軟骨や腱といった２つの薄い領域のマージンを避けるためである。即ち、別々に考慮される領域ｉ及びｊのコンパクトネスがＣｉ，Ｃｊ，及びＣｉｊで表され、領域のコンパクトネスが組み合わされる場合、バタチャリア距離がＴ（Ｃｉ＋Ｃｊ）／２Ｃｉｊよりも小さくない限り領域はマージンしない。ステップ７２５において、プロセスはセグメンテーション・リラクゼーションに進み、ステップ７２７において終了する。
【０１００】
上述したアプローチはいくつかのＭＲＩ体積像（ｖｏｌｕｍｅｔｒｉｃ  ｉｍａｇｅ）シーケンスで試された。膝の屈伸を例にとってみる。膝の６つの画像が取得される。各画像は、６０ｘ２５６ｘ２５６ボクセルの解像度を有し、各ボクセルの大きさは０．４７ｍｍ  ｘ  ０．４７ｍｍ  ｘ  １．４ｍｍである。
【０１０１】
図８Ａから８Ｄは、シーケンス中の各画像を示したものである。図８Ａは、未加工のデータを示したものである。図８Ｂは、ファイナル・セグメンテーションを示したものである。図８Ｃは、未加工のデータにファイナル・セグメンテーションを重ねたものからの予想領域の境界（ｅｓｔｉｍａｔｅｄ  ｒｅｇｉｏｎ  ｂｏｕｎｄａｒｉｅｓ）を示したものである。図８Ｄは、５番目の画像から６番目の画像への予想動作フローを示したものである。曲げた膝の６つの画像からなる外観表現（ｓｕｒｆａｃｅ  ｒｅｎｄｅｒｉｎｇ）を図９Ａから９Ｆに示す。
【０１０２】
前述した実施例及びその他の実施例は、図１０のブロック図に示すようなシステムにおいて実装されることが可能である。システム１０００は、画像データを入力する入力装置１００２、素材データベース及びそれらに代替するものを有する。入力装置１００２から入力される情報は、ワークステーション１００４で受信される。ワークステーション１００４は、ハードドライブ等の記憶装置１００６、４次元データを提供するために上述した処理を行う処理部１００８、及び外観表現などを表示するための４次元データを用意する画像表現エンジン１０１０とを有する。出力装置１０１２は、表現エンジン１０１０で表現された画像を表示するためのモニター、画像を記憶するためのビデオ・レコーダーのような更に記憶装置、又はその両方を含んでいても構わない。上述したように、実施例のワークステーション１００４及び画像表現エンジン１０１０は、シリコングラフィックス・インディゴ・ワークステーション及びＩｒｉｘエクスプローラ３Ｄグラフィックス・エンジンである。
【０１０３】
更に、図１０に示したシステム１０００及び更に一般的にここで説明した発明は、ディジタルクローンと共に、臨床医とのバーチャル・リアリティ対話を可能にする。図９Ａから図９Ｆに示すアウトプットにおいて、例えば、セグメント化された膝蓋骨（ｐａｔｅｌｌａ）は選択的に診察、回転、測定、又はその他の方法で操作される。更に、仮想膝蓋骨に、力が加えられ、上述したＦＥＭモデルによって、仮想膝蓋骨の反応が診察される。このバーチャル・リアリティ対話は、図１０に示すシステム１０００もしくは簡易的なバーチャル・リアリティ技術と本発明との組み合わせによるその他の同等及び適当なシステムにおいて実施することができる。
【０１０４】
本発明は上記のような最適な実施例を記載したが、この開示から当業者には様々な代替実施例が明らかとなろう。例えば、最適な実施例は、人間の組織のＭＲＩ分析に関する文脈で開示したが、画像化する対象は、人間のものではない組織や生き物でないものであっても構わない。また、画像シーケンスは、いかなる方法で取得されても構わず、例えば、Ｘ線やＣＴスキャン、若しくはその他の方法であっても構わない。更に、画像ｇはグレイ・スケール・画像として開示されているが、カラー若しくはｐｓｅｕｄｏ−ｃｏｌｏｒ（擬似色）画像、マルチメンショナル・画像、又はマルチスペクタル・画像であっても分析することができる。
【図面の簡単な説明】
【図１】  画像のセグメント化に関わるオペレーショナル・ステップのフローチャートを示す図である。
【図２】  図１のセグメント化に用いられるボクセルの近傍方式を示す図である。
【図３】  図１のセグメント化に用いられるボクセルの近傍方式を示す図である。
【図４】  図１のセグメント化に用いられるボクセルの近傍方式を示す図である。
【図５】  図１のセグメント化に用いられるボクセルの近傍方式を示す図である。
【図５Ａ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｂ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｃ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｄ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｅ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｆ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｇ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｈ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｉ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｊ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｋ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図５Ｌ】  図２から図５の近傍方式の２次元平面におけるボクセル間に引くことのできる線分を示す図である。
【図６】  連続画像のトラッキング及び動作見積もりに関わるオペレーショナル・ステップのフローチャートを示す図である。
【図７】  図１の処理の変形を示すフローチャートである。
【図８】  図８Ａから図８Ｄは、４段階でのセグメント化及び動作見積もりにおける屈曲膝画像を示す図である。
【図９】  図９Ａから図９Ｆは、動作見積もり及び外観表現後における、図８Ａから図８Ｄの屈曲膝の６画像を示す図である。
【図１０】  本発明を実行可能なシステムを示すブロック図である。
【図１】
0004360757.tif 000042
【図２】
0004360757.tif 000043
【図３】
0004360757.tif 000044
【図４】
0004360757.tif 000045
【図５】
0004360757.tif 000046
【図５Ａ】
0004360757.tif 000047
【図５Ｂ】
0004360757.tif 000048
【図５Ｃ】
0004360757.tif 000049
【図５Ｄ】
0004360757.tif 000050
【図５Ｅ】
0004360757.tif 000051
【図５Ｆ】
0004360757.tif 000052
【図５Ｇ】
0004360757.tif 000053
【図５Ｈ】
0004360757.tif 000054
【図５Ｉ】
0004360757.tif 000055
【図５Ｊ】
0004360757.tif 000056
【図５Ｋ】
0004360757.tif 000057
【図５Ｌ】
0004360757.tif 000058
【図６】
0004360757.tif 000059
【図７】
0004360757.tif 000060
【図８Ａ】
0004360757.tif 000061
【図８Ｂ】
0004360757.tif 000062
【図８Ｃ】
0004360757.tif 000063
【図８Ｄ】
0004360757.tif 000064
【図９Ａ】
0004360757.tif 000065
【図９Ｂ】
0004360757.tif 000066
【図９Ｃ】
0004360757.tif 000067
【図９Ｄ】
0004360757.tif 000068
【図９Ｅ】
0004360757.tif 000069
【図９Ｆ】
0004360757.tif 000070
【図１０】
0004360757.tif 000071