(57)【特許請求の範囲】
【請求項１】
組織の少なくとも１つの弾性率に関連するデータを決定するためのシステムであって、
所定の技術を実行する際に、
ａ）前記組織に関する情報を含む少なくとも１つの干渉計測信号を受信し、
ｂ）前記組織の前記弾性率を、前記受信した少なくとも１つの干渉計測信号及び少なくとも前記組織の速度分布の関数として決定するように構成された処理装置を備えるシステム。
【請求項２】
前記情報は、前記組織の構造及び前記組織の組成のうち、少なくとも１つを含む請求項１に記載のシステム。
【請求項３】
前記少なくとも１つの干渉計測信号は、第１の応力レベルにおける前記組織に関する第１の情報を含む第１の干渉計測信号と、第２の応力レベルにおける前記組織に関する第２の情報を含む第２の干渉計測信号とを含み、第２の応力が第１の応力とは異なる請求項１に記載のシステム。
【請求項４】
前記処理装置は、比較情報を生成するために、前記第１の情報と前記第２の情報とを比較するように更に構成され、前記弾性率が比較情報の関数として決定される請求項３に記載のシステム。
【請求項５】
前記処理装置は、前記第１の情報又は前記第２の情報の１つ以上の関数として、数値モデルを生成するように更に構成される請求項３に記載のシステム。
【請求項６】
前記数値モデルが動的な数値モデルである請求項５に記載のシステム。
【請求項７】
前記動的な数値モデルは、前記第１の情報又は前記第２の情報の１つ以上の関数として修正可能な、制限、モデルの複雑さ又はモデルの次数のうちの少なくとも１つを含む請求項６に記載のシステム。
【請求項８】
前記モデルの複雑さ又はモデルの次数は、前記第１の情報又は前記第２の情報の少なくとも１つの関数として修正可能である請求項７に記載のシステム。
【請求項９】
前記動的な数値モデルは、詳細情報を生成するために実行され、前記詳細情報は、前記制限、前記モデルの複雑さ又は前記のモデル次数のうちの１つ以上を修正するために、前記動的な数値モデルに対して提供される請求項７に記載のシステム。
【請求項１０】
前記モデルの複雑さは、複数のモデル要素を有し、該モデル要素のうちの少なくとも第１の要素は、前記モデル要素のうちの少なくとも第２の要素と関連し、該関連は前記モデル要素の前記第１の要素及び前記第２の要素の重みに基づいている請求項７に記載のシステム。
【請求項１１】
前記処理装置は、前記数値モデルを使用して前記組織に関する詳細情報を生成するように更に構成され、前記詳細情報は、前記組織に加える応力に対する前記組織の応答に関連している請求項５に記載のシステム。
【請求項１２】
前記処理装置は、比較情報と前記詳細情報の関数として、更に別のデータを生成するように更に構成される請求項１１に記載のシステム。
【請求項１３】
前記処理装置は、前記更に別のデータの関数として、前記数値モデルを修正するように更に構成される請求項１２に記載のシステム。
【請求項１４】
前記弾性率が前記数値モデルに基づいて決定される請求項１３に記載のシステム。
【請求項１５】
前記処理装置は、前記数値モデルに基づいて前記組織の歪み情報を決定するように更に構成される請求項１３に記載のシステム。
【請求項１６】
前記比較情報は更に、（ｉ）前記組織又は組織内の造影剤の少なくとも一方についての１つ以上の既知の特性又は（ｉｉ）干渉計システムの特性のうち、１つ以上に関連する詳細情報に依存する請求項４に記載のシステム。
【請求項１７】
前記詳細情報は、前記組織の速度分布、前記組織の圧縮的特性又は前記組織の弾性的特性の少なくとも１つ、組織型、及び前記組織の構造のうち、１つ以上を含む請求項１６に記載のシステム。
【請求項１８】
前記処理装置は、前記組織から得られるドップラ信号に基づいて前記組織の前記速度分布を決定するように更に構成される請求項１７に記載のシステム。
【請求項１９】
組織の少なくとも１つの弾性率に関連するデータを決定するための方法であって、
ａ）前記組織に関する情報を含む少なくとも１つの干渉計測信号を受信するステップと、
ｂ）前記組織の弾性率を、前記受信した少なくとも１つの干渉計測信号及び少なくとも前記組織の速度分布の関数として決定するステップと、を有する方法。
【請求項２０】
前記情報は、前記組織の構造及び前記組織の組成のうち、少なくとも１つを含む請求項１９に記載の方法。
【請求項２１】
前記少なくとも１つの干渉計測信号は、第１の応力レベルにおける前記組織に関する第１の情報を含む第１の干渉計測信号と、第２の応力レベルにおける前記組織に関する第２の情報を含む第２の干渉計測信号と、を含み、第２の応力が第１の応力とは異なる請求項１９に記載の方法。
【請求項２２】
比較情報を生成するために、前記第１の情報と前記第２の情報とを比較するステップを更に有し、前記弾性率が比較情報の関数として決定される請求項２１に記載の方法。
【請求項２３】
前記第１の情報又は前記第２の情報の１つ以上の関数として、数値モデルを生成するステップを更に有する請求項２１に記載の方法。
【請求項２４】
前記数値モデルを使用して前記組織に関する詳細情報を生成するステップを更に有し、前記詳細情報は、前記組織に加える応力に対する前記組織の応答に関連している請求項２３に記載の方法。
【請求項２５】
比較情報と前記詳細情報の関数として、更に別のデータを生成するステップを更に有する請求項２４に記載の方法。
【請求項２６】
前記更に別のデータの関数として、前記数値モデルを修正するステップを更に有する請求項２５に記載の方法。
【請求項２７】
前記弾性率が前記数値モデルに基づいて決定される請求項２６に記載の方法。
【請求項２８】
前記数値モデルに基づいて前記組織の歪み情報を決定するステップを更に有する請求項２６に記載の方法。
【請求項２９】
前記比較情報は更に、（ｉ）前記組織又は前記組織内の造影剤の少なくとも一方についての１つ以上の既知の特性及び（ｉｉ）光コヒーレンス断層撮影システムの特性のうち、１つ以上の詳細情報に依存する請求項２２に記載の方法。
【請求項３０】
前記詳細情報は、前記組織の速度分布、前記組織の圧縮的特性又は前記組織の弾性的特性の少なくとも１つ、組織型、及び前記組織の構造のうち、１つ以上を含む請求項２９に記載の方法。
【請求項３１】
前記組織から得られるドップラ信号に基づいて前記組織の速度分布を決定するステップを更に有する請求項３０に記載の方法。
【請求項３２】
組織の少なくとも１つの弾性率に関連するデータを決定するためのソフトウェアであって、
処理装置によって実行される際に、前記組織に関する情報を含む少なくとも１つの干渉計測信号を受信するように前記処理装置を設定する第１の命令群と、
前記処理装置によって実行される際に、前記組織の弾性率を、前記受信した少なくとも１つの干渉計測信号及び少なくとも前記組織の速度分布の関数として決定するように前記処理装置を設定する第２の命令群と、を有するソフトウェア。
【請求項３３】
前記少なくとも１つの干渉計測信号は、第１の応力レベルにおける前記組織に関する第１の情報を含む第１の干渉計測信号と、第２の応力レベルにおける前記組織に関する第２の情報を含む第２の干渉計測信号とを含み、
前記処理装置によって実行される際に、比較情報を生成するために、前記第１の情報と前記第２の情報とを比較するように前記処理装置を設定する第３の命令群と、
前記処理装置によって実行される際に、組織の前記少なくとも１つの構造的変化に関連する前記データを、前記比較情報と、（ｉ）前記組織の少なくとも１つの既知の特性又は（ｉｉ）ＯＣＴシステムの特性のうちの少なくとも１つに関連する詳細情報と、の関数として決定するように前記処理装置を設定する第４の命令群と、をさらに含む請求項３２に記載のソフトウェア。
【請求項３４】
前記少なくとも１つの構造的変化は、前記組織の歪みである請求項３３に記載のソフトウェア。
【請求項３５】
第２の応力が第１の応力とは異なる請求項３４に記載のソフトウェア。
【請求項３６】
前記処理装置は、前記干渉計測信号の少なくとも１つのドップラ周波数偏移に基づいて前記速度分布を決定するように構成される請求項１に記載のシステム。
【請求項３７】
前記処理装置は、前記干渉計測信号におけるスペックルに基づいて前記速度分布を決定するように構成される請求項１に記載のシステム。
【請求項３８】
前記干渉計測信号の少なくとも１つのドップラ周波数偏移に基づいて前記速度分布を決定するステップをさらに含む請求項１９に記載の方法。
【請求項３９】
前記干渉計測信号におけるスペックルに基づいて前記速度分布を決定するステップをさらに含む請求項１９に記載の方法。
【請求項４０】
前記処理装置は、前記干渉計測信号の少なくとも１つのドップラ周波数偏移に基づいて前記速度分布を決定するように構成される請求項３２に記載のソフトウェア。
【請求項４１】
前記処理装置は、前記干渉計測信号におけるスペックルに基づいて前記速度分布を決定するように構成される請求項３２に記載のソフトウェア。
【発明の詳細な説明】
【技術分野】
【０００１】
関連出願への相互参照
本発明は、２００４年８月２４日出願の米国特許出願番号第６０／６０４，１３７号の利益を請求し、その全開示が参照によりここに援用される。
【０００２】
本発明は一般に、試料の機械的歪み及び弾性的性質の測定に関し、より詳細には、生体試料の機械的歪み及び弾性的性質の空間的分布を非侵襲的に測定して決定するためのプロセス、システム及びソフトウェアに関する。
【背景技術】
【０００３】
心筋梗塞は、先進工業国の主要な死因である。脆弱な動脈硬化性プラークの破裂は、多くの場合に突然死の原因となる急性心筋梗塞の重要なメカニズムとして認識されている。心臓血管系の研究における最近の進歩により、動脈硬化性プラークに破裂傾向をもたらす該プラークの構造及び組成上の特徴が識別された。大多数の脆弱なプラークには、ａ）プラークの端部や肩部における活性化マクロファージの存在、ｂ）薄い不安定な線維性の被膜、ｃ）弾力的な脂質プールといった、３つの特徴がある。これらの３つの特徴によって表される生化学的に惹き起こされる弱化と、機械的応力の増大との組み合わせは、特にリスクの高い状況を表す虞がある。
【０００４】
破裂するほどに脆弱なプラークを探知可能な技術は、リスクのある患者の管理を指導するための有益なツールとなり、急性疾患の最終的な予防を支援することができる。脆弱なプラークを検出するための数多くの異なる方法が研究されている。これらの方法には、血管内超音波診断法（「ＩＶＵＳ」）、光コヒーレンス断層撮影法（「ＯＣＴ」）、蛍光分光法、磁気共鳴映像法（「ＭＲＩ」）、コンピュータ断層撮影法（「ＣＴ」）、陽電子放出型断層撮影法（「ＰＥＴ」）、赤外分光法がある。
【０００５】
ＯＣＴは、参照光線と試料から反射された検出ビームとの間の干渉を測定することができる画像化技術である。従来の時間領域ＯＣＴのシステムに関する詳細な説明は、非特許文献１に記載されている。スペクトル領域でのＯＣＴの改良は、スペクトル領域光コヒーレンス断層撮影法（「ＳＤ−ＯＣＴ」）と呼ばれており、極めて高解像度での眼科の画像診断に適した方法である。この方法は、非特許文献２及び特許文献１に記載されている。更に、２００２年１０月１６日出願の特許文献２、非特許文献３、非特許文献４も、この主題に関連している。更に（非特許文献５、特許文献１、及び２００４年１０月２７日出願の特許文献３に記載されているように）、光周波数領域干渉分光法（「ＯＦＤＩ」）の装置も本発明の対象に関連する。
【０００６】
ＳＤ−ＯＣＴ法とＯＦＤＩ法は、組織の高解像断面画像を提供するという点で、ＯＣＴ法と似ている。このような例示的な方法も、組織組成の正確な特徴化を可能とし、画像取得レートの大幅な向上をもたらす。これらの例示的な変形技術を、ここではまとめてＯＣＴと呼ぶことにする。提案されている上記の技術のうち、ＯＣＴ法は、プラーク破裂の直接的原因であると考えられている構造及び組成の特徴を空間的に解明する能力をもつことがわかっている。しかし、構造及び組成の特徴をただ知っているだけでは、プラーク破裂を詳細に理解し、高精度で予測するには不十分である。構造的及び組成的な情報を、歪み及び弾性率の測定と結び付ける技術が望ましい。
【０００７】
ある種の数値的方法（例えば有限要素分析法）は、機械的応力及び歪みと、プラーク破裂におけるこれらの役割を理解するために使用されてきた。現在の各種分析は、組織学及びＩＶＵＳに漫然と基づく血管断面モデルを頼りにしており、組織の弾性的性質の推定値又は間接的な測定値が得られている。これらの数値的方法により、プラーク破裂についてある程度の知見が得られてはいるが、これらの方法は不利である。その理由は、例えば、ａ）その精度について弾性的性質及びその分布についての知識が不正確であるために制限されること、ｂ）これらの方法が、遡及的データに基づいており、生きた患者の血管構造の診断に直接適用できないためである。
【０００８】
ＩＶＵＳ弾性率計測法は、生体内で血管構造内の歪みを測定するための方法として開発されてきた。この例示的な方法は、血管内の圧力が変化する間に、多数の断面画像を得ることによって実行可能である。これらの画像を相関させることによって、圧力変化に対する血管の機械的応答が調べられ、この結果、歪み、局所的変位、変形又は空間分解速度の断面マップが得られる。この方法は生体内で実施可能であるものの、動脈硬化性プラークの典型的な組織成分同士の間での空間分解能が低く、コントラストが低い。更に、このような方法は、歪みと独立して応力を決定することができず、このため、弾性率分布を決定できない可能性がある。ＯＣＴ弾性率計測法は、ＩＶＵＳ弾性率計測法に使用される技術に関連する技術に基づく。このＯＣＴ弾性率計測法は、原則的には、ＩＶＵＳ弾性率計測法よりも高い解像度と、相対弾性率分布を提供することができる。また、動脈内腔における圧力負荷に関する知識と組み合わせれば、絶対弾性率の高解像度での評価が可能である。
【０００９】
ＩＶＵＳ及びＯＣＴと併用したドップラ画像化法は、画像化プローブに近づくか又はこれから離れる試料の、深さ分解速度を決定するために使用されてきた。これらの方法の変形例がいくつか知られてはいるが、これらに共通する基盤は、試料中の散乱体を移動させることによって、プローブビーム、つまり、ＩＶＵＳの超音波やＯＣＴの光に与えられるドップラ周波数偏移を測定することにある。
【先行技術文献】
【特許文献】
【００１０】
【特許文献１】国際公開第０３／０６２８０２号パンフレット
【特許文献２】米国特許出願第１０／２７２，１７１号明細書
【特許文献３】米国仮特許出願第６０／５１４，７６９号明細書
【特許文献４】米国特許出願第１０／１３６，８１３号明細書
【特許文献５】国際特許出願第ＰＣＴ／ＵＳ０３／０２３４９号明細書
【非特許文献】
【００１１】
【非特許文献１】Ｈｕａｎｇ他、“Ｏｐｔｉｃａｌ  ｃｏｈｅｒｅｎｃｅ  ｔｏｍｏｇｒａｐｈｙ，”  Ｓｃｉｅｎｃｅ  ２５４  （５０３５），  １１７８−８１  （１９９１）
【非特許文献２】Ｃｅｎｓｅ  Ｂ．他、“Ｕｌｔｒａｈｉｇｈ−ｒｅｓｏｌｕｔｉｏｎ  ｈｉｇｈ−ｓｐｅｅｄ  ｒｅｔｉｎａｌ  ｉｍａｇｉｎｇ  ｕｓｉｎｇ  ｓｐｅｃｔｒａｌ−ｄｏｍａｉｎ  ｏｐｔｉｃａｌ  ｃｏｈｅｒｅｎｃｅ  ｔｏｍｏｇｒａｐｈｙ”，  Ｏｐｔｉｃｓ  Ｅｘｐｒｅｓｓ，  ２００４
【非特許文献３】Ｗｏｊｔｋｏｗｓｋｉ他、“Ｉｎ  Ｖｉｖｏ  Ｈｕｍａｎ  Ｒｅｔｉｎａｌ  Ｉｍａｇｉｎｇ  ｂｙ  Ｆｏｕｒｉｅｒ  Ｄｏｍａｉｎ  Ｏｐｔｉｃａｌ  Ｃｏｈｅｒｅｎｃｅ  Ｔｏｍｏｇｒａｐｈｙ”，  Ｊｏｕｒｎａｌ  ｏｆ  Ｂｉｏｍｅｄｉｃａｌ  Ｏｐｔｉｃｓ，  ２００２，  ７（３），  ｐｐ．  ４５７−４６３
【非特許文献４】Ｎａｓｓｉｆ，  Ｎ．他、“Ｉｎ  Ｖｉｖｏ  Ｈｕｍａｎ  Ｒｅｔｉｎａｌ  Ｉｍａｇｉｎｇ  ｂｙ  Ｕｌｔｒａｈｉｇｈ−Ｓｐｅｅｄ  Ｓｐｅｃｔｒａｌ  Ｄｏｍａｉｎ  Ｏｐｔｉｃａｌ  Ｃｏｈｅｒｅｎｃｅ  Ｔｏｍｏｇｒａｐｈｙ”，  Ｏｐｔｉｃｓ  Ｌｅｔｔｅｒｓ，  ２００４，  ２９（５），  ｐｐ．  ４８０−４８２
【非特許文献５】Ｙｕｎ，  Ｓ．Ｈ．他、“Ｈｉｇｈ−Ｓｐｅｅｄ  Ｏｐｔｉｃａｌ  Ｆｒｅｑｕｅｎｃｙ−Ｄｏｍａｉｎ  Ｉｍａｇｉｎｇ”，  Ｏｐｔｉｃｓ  Ｅｘｐｒｅｓｓ，  ２００３，  １１（２２），  ｐｐ．  ２９５３−２９６３
【非特許文献６】Ａ．Ｒ．  Ｓｋｏｖｏｒｏｄａ他、“Ｔｉｓｓｕｅ  ｅｌａｓｔｉｃｉｔｙ  ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ  ｂａｓｅｄ  ｏｎ  ｕｌｔｒａｓｏｎｉｃ  ｄｉｓｐｌａｃｅｍｅｎｔ  ａｎｄ  ｓｔｒａｉｎ  ｉｍａｇｅｓ”、ＩＥＥＥ  Ｔｒａｎｓ  Ｕｌｔｒａｓｏｎ  Ｆｅｒｒｏｅｌｅｃｔｒ  Ｆｒｅｑ  Ｃｏｎｔｒｏｌ，  Ｃｏｌ．  ４２，１９９５，  ｐｐ．  ７４７−７６５
【非特許文献７】Ｆ．  Ｋａｌｌｅｌ他、“Ｔｉｓｓｕｅ  ｅｌａｓｔｉｃｉｔｙ  ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ  ｕｓｉｎｇ  ｌｉｎｅａｒ  ｐｅｒｔｕｒｂａｔｉｏｎ  ｍｅｔｈｏｄ”，  ＩＥＥＥ  Ｔｒａｎｓ  Ｍｅｄ  Ｉｍａｇｉｎｇ，  Ｖｏｌ．  １５，  １９９６，  ｐｐ．  ２９９−３１３
【非特許文献８】Ｒ．Ｃ．  Ｃｈａｎ他、“ＯＣＴ−ｂａｓｅｄ  ａｒｔｅｒｉａｌ  ｅｌａｓｔｏｇｒａｐｈｙ：  ｒｏｂｕｓｔ  ｅｓｔｉｍａｔｉｏｎ  ｅｘｐｌｏｉｔｉｎｇ  ｔｉｓｓｕｅ  ｂｉｏｍｅｃｈａｎｉｃｓ”，  Ｏｐｔｉｃｓ  Ｅｘｐｒｅｓｓ，  Ｖｏｌ．  １２（１９），  ２００４，  ｐｐ．  ４５５８−４５７２
【発明の概要】
【発明が解決しようとする課題】
【００１２】
しかしながら、試料の構造、組成及び生体力学的特性を同時に決定する方法が利用可能でない。このような能力は、生物医学において広い用途をもつが、特に、脆弱なプラークを見つけだして、急性心筋梗塞との関連性を理解する上で有効となろう。
【００１３】
更に、弾性率計測法及び弾性率画像化法は一般に、多数の独立した有限要素、つまり対象領域で空間的に分布している画像ピクセルに亘る未知の歪み又は弾性率パラメータの推定を使用する。歪み又は弾性率の画像化に使用する空間分解能が高くなるほど、推定すべき独立の未知数の数が増える。パラメータ空間が大きくなるにつれて、所望の目的関数を満たすパラメータの推定についての探索が、困難な劣決定問題（ｕｎｄｅｒｄｅｔｅｒｍｉｎｅｄ  ｐｒｏｂｌｅｍ）となってしまう虞がある。通常は、未知数の数が、基礎となる画像化データから一意に決定することができる数を大幅に超えると、目的関数を満たす多くの解が得られることになる。また、高次元（１００次元を超えるオーダー）のパラメータ空間を探索するためには、一般に高い計算コストと長い計算時間がかかる。
【００１４】
生物組織の弾性率計測法及び弾性率画像化のための従来の方法は、有限要素毎に又は対象のピクセルでの歪み又は弾性率を独立の未知数として処理する。その際、通常は、上記非特許文献６、非特許文献７に記載されるように、目的関数の最適化のためのＬｅｖｅｎｂｕｒｇ−Ｍａｒｑｕａｒｄｔ又は同様のアルゴリズムを使用する。極小値に対するロバスト性を実現するために、多重解像度法が用いられ、この方法では、少ない未知数を用いて低解像度のグリッド（格子）で推定値が得られ、次に、最大解像度の領域におけるパラメータ最適化の初期化のために、これらの低解像度の推定値が高解像度のグリッドにマッピングされる。これらの従来の方法は時間がかかり、対象領域が大きいと数分の処理を要する場合がある。
【課題を解決するための手段】
【００１５】
従来の方法とは対照的に、本発明によるシステム、プロセス、及びソフトウェア構成の例示的な実施形態は、少なくとも１つの試料における歪み及び弾性率の空間分布を、高い空間分解能及び高感度で決定することができるとともに、構造及び組成の高解像度画像を同時に提供できる。本発明のシステム、プロセス、及びソフトウェア構成は広く適用可能であり、その能力は特に生物組織及び血管組織に関連している。
【００１６】
本発明の１つの例示的な実施形態では、血管の構造及び組織の組成を決定するためにＯＣＴが使用される。そして、この情報が、血管を表現する数値モデルを構築するために使用され、その後、任意の応力負荷に対する血管の機械的応答を予測するために、弾性率の推定を用いた有限要素モデリングが使用される。この例示的な計算とは別に、本発明による例示的なＯＣＴ弾性率計測法は、血管の機械的応答を測定するために実行される。モデリングと画像化という２つの経路は、ａ）想定された弾性率分布に基づく予測と、ｂ）測定を、それぞれに表している。この２つの結果の差については、弾性率分布の初期推定を修正することによって最小化される誤差関数とみなすことができる。本発明によるこの例示的な方法を反復することによって、弾性率の分布及び大きさを決定することができる。この情報は、弾性率の断面画像又は三次元画像として表示される。更に、誤差関数を最小化することによって、改善された、歪みの弾性率計測法による画像を生成することができる。その結果、本発明によるシステム、プロセス、及びソフトウェア構成の例示的な実施形態では、現在の診断技術の制約、すなわち、構造及び／又は歪みは測定されるが、組織の生体力学的特性が不明であるといった制約を解消することができる。更に、本発明は、弾性率計測法のための従来法の解像度及び感度を改善する。
【００１７】
要約すると、本発明によるシステム、プロセス、及びソフトウェア構成の例示的な実施形態によって、医学的用途や非医学的な用途のために、試料の構造、組成、歪み及び弾性率を同時に決定できるようになる。
【００１８】
本発明の１つの例示的な実施形態では、組織の少なくとも１つの構造的変化に関連するデータを決定するためのシステム、プロセス及びソフトウェア構成が提供される。より詳細には、第１の応力レベルでの組織に関する第１の情報を含む第１の光コヒーレンス断層撮影法（「ＯＣＴ」）信号と、第２の応力レベルでの組織に関する第２の情報を含む第２のＯＣＴ信号と、が受信される。第１の情報と第２の情報は、比較情報を生成するために比較される。少なくとも１つの構造的変化に関連するデータが、比較情報と、（ｉ）組織についての１つ以上の既知特性及び／又は（ｉｉ）ＯＣＴシステムの特性に関連する詳細情報との関数として決定される。
【００１９】
例えば、前記構造的変化は、組織の歪みである。また、第２の応力は第１の応力と異なる。前記詳細情報には、組織の速度分布、組織の機械的特性（例えば、圧縮可能性（ｃｏｍｐｒｅｓｓａｂｉｌｉｔｙ）及び／又は弾性的性質）、組織型、組織中の造影剤の光学特性、及び／又は組織の構造がある。更に、組織の速度分布は、組織から得られるドップラ信号に基づいて決定される。
【００２０】
本発明の別の例示的な実施形態によれば、組織の少なくとも１つの弾性率に関連するデータを決定するための方法、システム、及びソフトウェア構成が提供される。例えば、組織に関する情報を含む、少なくとも１つの光コヒーレンス断層撮影法（「ＯＣＴ」）信号が受信される。そして、組織の弾性率が、受信された１つ以上のＯＣＴ信号の関数として決定される。
【００２１】
例えば、この情報には、組織の構造及び／又は組織の組成を含むことができる。ＯＣＴ信号は、第１の応力レベルでの組織に関する第１の情報を含む第１のＯＣＴ信号と、第２の応力レベルでの組織に関する第２の情報を含む第２のＯＣＴ信号とを含むことができ、その際、第２の応力が第１の応力とは異なってもよい。第１の情報と第２の情報については、比較情報を生成するために比較することができ、弾性率は比較情報の関数として決定される。また、数値モデルが、第１の情報又は第２の情報のうちの１つ以上の関数として生成される。組織に関する詳細情報を、数値モデルを使用して生成することができ、この詳細情報は、組織に加える応力に対する組織の応答に関連している。
【００２２】
前記数値モデルは動的な数値モデルとすることが可能であり、この数値モデルは、（ｉ）制限、（ｉｉ）モデルの複雑さ、及び／又は（ｉｉｉ）モデルの次数を含むことができ、これらについては、第１の情報及び／又は第２の情報の関数として修正可能である。モデルの複雑さ及び／又はモデルの次数は、第１の情報及び／又は第２の情報の関数として修正可能である。動的な数値モデルは詳細情報を生成するために実行され、この詳細情報は、制限、モデルの複雑さ及び／又はモデルの次数を修正するために動的な数値モデルに与えられる。モデルの複雑さには、複数のモデル要素が含まれ、該モデル要素のうち、少なくとも第１の要素は、該モデル要素における第１の要素及び／又は第２の要素の重みに基づいて、複数の要素と関連をもつ。
【００２３】
また、前記比較情報及び詳細情報の関数として、更に別のデータが生成される。数値モデルは、この更なるデータの関数として修正される。そして、弾性率が数値モデルに基づいて決定される。組織の歪み情報は数値モデルに基づいて得られる。前記比較情報は更に、詳細情報、つまり、（ｉ）組織の少なくとも１つの既知の特性、及び／又は、（ｉｉ）ＯＣＴシステムの特性についての情報に依存する。この詳細情報には、組織の速度分布、組織の圧縮性／弾性的性質、組織型、組織内の造影剤の光学特性、及び／又は組織の構造が含まれる。組織の速度分布は、組織から得られるドップラ信号に基づいて決定される。
【００２４】
更に、生体力学的画像化のための従来の方法及び技術とは対照的に、本発明の別の例示的な実施形態では、同じ型の組織が同様の、そしておそらくはほぼ同一の機械的な特性を有することや、対象領域を、区別できる組織構成要素（例えば、繊維性、脂質性、石灰化など）に区分化するために、組織の高コントラストがＯＣＴ法に利用可能であることを考慮している。本発明によるプロセスの、このような例示的な実施形態は、組織間に存在する境界を保持したままで、パラメータの探索空間を低減させることができる。例えば、低解像度グリッドでの推定とは異なり、この技術によって、各要素中の組織型の部分的なボリューミングが最小化され、これにより歪み又は弾性率の急峻な空間的勾配が保たれる。また、生体力学モデルがデータに十分良好に適合していない要素の適応的メッシュの精緻化は、弾性率計測法及び弾性率画像化法に対して有益に寄与することとなる。
【００２５】
本発明の上記の目的及び他の目的、特徴及び利点は、添付の特許請求の範囲と併せて、以下に記載する本発明の詳細な説明を読む際に、更に明白となるであろう。
【００２６】
本発明の更に別の目的、特徴及び利点については、本発明の例示的な実施形態を示す図面を参照して、以下の発明の詳細な説明から明らとなるであろう。
【図面の簡単な説明】
【００２７】
【図１】病変血管の断面を示す例示的な概略図である。
【図２】本発明に従って弾性率分布を決定するために、ＯＣＴ弾性率計測法と例示的な有限要素モデリング技術とを組み合わせた、一般的なフローチャートである。
【図３】弾性率分布を決定するために速度分布を用いる、本発明による例示的な方法のフローチャートである。
【図４】弾性率分布及び歪み分布を決定するために構造及び構造的変形を用いる、本発明による別の例示的な方法のフローチャートである。
【図５】本発明による多重解像度の速度場推定のための、例示的な方法のブロック図である。
【図６】本発明の例示的な実施形態の結果を検証するために、実験で使用される例示的な有限要素形状と有限要素メッシュをそれぞれに示す図である。
【図７】フリンジ分解測定によりシミュレーションされたＯＣＴ点広がり関数、及び組織ブロック内の内包物のシミュレーションされたＯＣＴ画像を例示的に示す図である。
【図８】本発明の例示的な実施形態による、弾力的な内包物の軸方向速度場を例示的に示す図であり、フレーム（４００）は有限要素モデリングによる例示的な真の軸方向速度を示し、フレーム（４０５）は、従来の運動追跡による例示的な軸方向速度推定を示し、フレーム（４１０）は、例示的な多重解像度変分法による例示的な軸方向速度推定を示す。
【図９】図８と同様のフレームの画像表現に対応する、フレーム（４５０）、（４５５）及び（４６０）を用いて、固い内包物の軸方向速度場を例示的に示す図である。
【図１０】図８、９と同様のフレームの画像表現に対応する、フレーム（５００）、（５０５）及び（５１０）を用いて、弾力的な内包物の軸方向の歪場を例示的に示す図である。
【図１１】図８、９及び１０と同様のフレームの画像表現に対応する、フレーム（５５０）、（５５５）及び（５６０）を用いて、固い内包物の軸方向の歪場を例示的に示す図である。
【図１２】大動脈検体の例示的なＯＣＴ画像である。
【図１３】横方向の伸張を受けている、図１２の大動脈検体について例示的な横方向速度分布を示す図である。
【図１４】横方向の伸張を受けている、図１２の大動脈検体について例示的な横方向歪み分布を別に示す図である。
【図１５】パラメータの効果的な低減に使用可能な、ＯＣＴに基づく生体力学特性の推定のために、本発明による方法の別の例示的な実施形態について一例を示すフローチャートである。
【発明を実施するための形態】
【００２８】
本発明の特定の例示的な実施形態では、試料の構造、組成、歪み及び／又は弾性率を決定するために、ＯＣＴ法を有限要素モデリング法と組み合わせたハイブリッド技術を利用することができる。
【００２９】
図１は、正常な血管壁４内に埋設された脂質プール３を有する、罹患動脈の断面を例示的に示す。内腔５内の血圧変動によって、血管及びプラーク形状の変形が生じる。例えば、図１では、腔内の圧力が上昇すると、点線の外形線１が実線の外形線２の位置へと変形する。ここに記載するＯＣＴ弾性率計測法の例示的な実施形態では、血管及びプラークの境界の変位を追跡することはもとより、組織自体内で発生する生体力学的な歪みをも推定することができる。
【００３０】
本発明によるＯＣＴ弾性率計測法の１つの例示的な実施形態によれば、試料内の散乱体の速度（例えば大きさ及び方向）を決定することができる。また、ＯＣＴの構造的な画像と組織の弾性率推定に基づいて、ある例示的な有限要素モデリング法を用いることで、対応する速度分布を予測することが可能になる。これら２つの例示的な速度分布間の差は、弾性率の初期推定について最適化を繰り返すことによって最小化される誤差関数とみなせる。そして、結果的に得られた弾性率分布が、画像形態で視覚化される。更に、歪みを画像形態で視覚的に表現するために、ＯＣＴ弾性率計測法のデータを使用することができる。
【００３１】
図２は、本発明に従い、弾性率分布を決定するために、ＯＣＴ弾性率計測法と有限要素の例示的なモデリング法とを組み合わせて、別々のステップとして図示的に示す一般的なフローチャートである。処理の最初には、例えば、動脈が心臓周期に亘って動的な変形を受ける際に、ＯＣＴの画像化及び取得がステップ５０で行われる。同時に、ステップ５５において、腔内の圧力がデジタル化されるか、対応するＯＣＴフレームで記録されるか、この両方が行われる。ステップ６５において、１つの圧力レベルから取得されたＯＣＴ画像が、罹患血管のための幾何形状モデルの基礎となり、これは、ステップ７０において、有限要素モデリング（ＦＥＭ）による数値シミュレーションのためにメッシュ化される。ステップ６０において組織の速度場を得るために、時間の関数としてのＯＣＴ画像データの変化が、運動推定のための例示的な方法を用いて追跡される。次に、ステップ９０において、測定された速度場から、これに対応する歪みの固有値及び固有ベクトルが計算され又は決定される。ステップ９５において、得られた組織の歪みの画像が、画像表示される。
【００３２】
また、ステップ６０で生成された、推定済みの組織の速度は、モデルに基づいた弾性率７０の決定に対して、その基礎となる。ステップ７５では、測定された圧力負荷（ステップ５５）及び既知の境界条件を、弾性率の値についての既定の分布でもって初期化した有限要素メッシュに適用することによって、数値モデルが実行され、予測の速度分布が得られる。ステップ８０において、２乗誤差測定（ｓｑｕａｒｅｄ−ｅｒｒｏｒ−ｍｅａｓｕｒｅ）法を使用することによって、ステップ７５で得られたモデル予測による速度が、測定された速度６０と比較される。ステップ８０での比較によって生じる誤差に基づいて、弾性率の値が更新されて、ステップ６５にてモデルが再構築され、次にステップ７０において再シミュレーションが行われて、新しい予測速度の組が得られる。このプロセスは、弾性率の推定値が指定された許容レベルに収束するまで続けられる。収束の後に、ステップ８５では、最終的に得られた例示的な弾性率分布が画像表示される。
【００３３】
図３は、本発明による例示的な方法、つまり、例示的なＯＣＴ法を使用して弾性率分布を決定するために速度分布を使用する方法についてフローチャートを示す。より詳細には、（ステップ１００の）ＯＣＴ画像取得と（ステップ１０５の）腔内圧力の記録が同時に実行される。本発明のこの例示的な実施形態では、ステップ１１０において、多重解像度変分法を使用して、動的なＯＣＴデータセットが処理され、その結果は、２つの画像化時点間の組織変位のロバストな推定となる。組織歪みの固有値と固有ベクトルは、ステップ１４０において速度推定から決定され、そして、ステップ１４５において視覚的に画像表示される。
【００３４】
例示的な弾性率推定では、ステップ１１５において、血管及びプラークの表面を抽出するために、参照時点でのＯＣＴデータが区分化される。また、これらの表面については、動脈特有の形状を画定するために、三次元で再構築し得る。そして、この血管形状がメッシュ化され、境界条件が適用されて、ステップ１２０においてメッシュ要素が初期の弾性率の値を割り当てられるが、これは、ステップ１２５の有限要素モデリング法で更に使用するためである。この例示的な処理／方法から、２乗誤差測定値を決定するために、測定された組織速度と共に使用される組織の速度予測の組がステップ１３０で得られ、この２乗誤差測定は、数値モデルで使用される弾性率値及び境界条件の更新を推し進める。弾性率の更新と数値シミュレーションの例示的な方法は、２乗誤差測定値を最小化して、指定の許容レベルに収束し得る弾性率の推定値を生成するために、反復的に続行される。最終的な弾性率分布は、ステップ１３５において視覚的に画像表示される。
【００３５】
本発明の別の例示的な実施形態では、加える応力を変えて測定しながら、試料の構造及び組成を決定するためにＯＣＴ法を実行することができる。ＯＣＴ画像取得レートについては、個々の画像内における運動への顕著な人為的影響を回避できるように十分大きい。加えた応力の１つの値について決定される構造及び組成は、組織を表わす数値モデルを生成するために使用され、そして、数値モデリングは、測定された応力変化及び弾性率分布の初期推定を組み込むことで、第２の応力に対する構造を予測するために使用される。対応する第２の応力で取得されたＯＣＴ画像は予測された構造と比較され、予測された構造と測定された構造との差が、弾性率分布の初期推定を反復的に最適化することによって最小化される。この例示的な実施形態では、数値モデリング、例えば、最適化された弾性率分布に基づくモデリングが、弾性率分布と歪みを統合した最終的な決定のために使用される。この結果が、画像形態で視覚的に表示される。
【００３６】
本発明の更に別の例示的な実施形態では、試料の構造及び組成を決定するためにＯＣＴ法が実行され、この情報からは、組織を表現する数値モデルが生成される。数値モデリングは、試料の弾性率分布の推定に基づき、速度分布を予測するために使用され、該速度分布は、加えた応力に対する応答として試料内に生じるものである。そして、試料中を移動する散乱体からのＯＣＴビームの反射によって生じるドップラ周波数偏移は、画像強度データに加えて、試料中の深さ分解速度分布を決定するために使用される。モデルによる速度予測値と、ＯＣＴドップラデータ及び画像強度データからの速度測定値との差は、初期弾性率分布を反復的に最適化することによって最小化される。そして、得られた弾性率分布が、画像形態で可視化される。これとは独立に、試料中の歪みを視覚的に表現するために、ドップラＯＣＴデータを使用することができる。
【００３７】
図４は、本発明による別の例示的な方法のフローチャートを示し、この方法では弾性率分布及び歪み分布を決定するために構造及び構造的変形を使用する。ステップ２００において、例示的なＯＣＴ画像化法が実行され、同時にステップ２０５で腔内圧力の記録が実行される。得られたデータは、ステップ２２０における画像及びドップラ情報を有する探索データセットと、ステップ２１０における数値モデル構築のために処理される参照データセットと、に分けられる。そして、ステップ２１５では、参照データセットと探索データセットとの間で、有限要素法シミュレーションとモデルの剛体運動推定との統合のために、参照形状と圧力負荷が使用される。ステップ２２５では、推定された剛体モデルの変換が、参照強度データの再サンプリングを行うために、モデル予測によるメッシュ変形（例えば、参照強度データを、参照する探索データセットフレームへと効果的に歪ませること）と組み合わされる。ステップ２３０では、歪められた参照データと測定された探索データが、ＯＣＴで指定される目的関数内で結合される。未知の弾性率値は、目的関数を最大化するために、モデル構築ステップ２１０にて反復的に更新される。弾性率の推定値が収束すると、ステップ２３５、２４０において、数値モデルからの対応する弾性率及び歪みの分布が出力され、視覚的に画像表示される。
【００３８】
光コヒーレンス弾性率計測法
光コヒーレンス弾性率計測法は、好ましくは、超音波弾性率計測法の基礎をなす原理と同じ原理に基づく。例えば、組織が機械的負荷の下で画像化されると、その変位は、巨視的な構造（例えば組織界面）に対応する画像特徴に現れる。また、運動はコヒーレント画像のスペックルに現れるが、これは微視的な組織散乱体の空間分布が負荷の下で変化するためである。巨視的構造及び微視的スペックルからの運動推定では、連続した画像同士の間で画像特徴が良好に保たれると仮定している。このため、所望の速度推定は、異なる負荷条件下で取得した、参照画像内のブロックと、探索画像内のブロックとの間の類似性測度を最大化し得る。ＯＣＴで得られる干渉画像は、ビーム減衰をモデル化する指数関数的な減衰項と、空間的な畳み込みとの積によって近似することができる。
【００３９】
【数１】
0005524904.tif 000002
【００４０】
座標ｘは側方の走査方向に対応し、そして座標ｙは軸上の走査方向に対応しており、前者はサンプルビーム軸に直交する方向を向き、後者はサンプルビーム軸に対して平行な方向を向いている。上付きバー記号を付したパラメータμｓは、試料の散乱による平均減衰量であり、σｂ（ｘ，ｙ）は、様々な後方散乱断面積を有する、点の分布として試料内の後方散乱をモデル化しており、ｈ（ｘ，ｙ）は、ＯＣＴシステムの点広がり関数（ＰＳＦ）を表している。ＯＣＴのＰＳＦは、ソース検出光学系の自己相関関数Γ（ｙ）と、瞳関数ｐ（ｘ）との間での、分離可能な空間的不変積として近似できる。
【００４１】
【数２】
0005524904.tif 000003
【００４２】
ガウス型ビームの場合、以下のようになる。
【００４３】
【数３】
0005524904.tif 000004
【００４４】
【数４】
0005524904.tif 000005
【００４５】
ここで、ＥＳは光源の電界ベクトルの振幅、０は光源の自由空間中心波長、Δは光源のＦＷＨＭ（半値全幅）でのスペクトルバンド幅、ｆは対物レンズの焦点距離、Ｄは対物レンズの入射瞳におけるサンプルビームの１／ｅ２強度直径である。
【００４６】
この画像形成モデルに基づいて、参照画像における位置（ｘ０，ｙ０）から探索画像内の新位置（ｘ０＋ｕ，ｙ０＋ｖ）へと変位を受ける１つの点状散乱体は、これに対応する以下の参照及び探索の干渉像を有することになる。
【００４７】
【数５】
0005524904.tif 000006
【００４８】
従来の速度測定を用いて、組織の運動は、式（５）及び（６）のエンベロープ又は複素振幅のサブブロック間での相関係数を最大化することによって概ね推定される。各画像は、事前に決められたサイズのブロックへと更に分割できる。各参照ブロックについて、全ての探索画像ブロックとの相互相関が計算されて、相関係数が、参照位置と探索位置との間の相対変位の関数として得られる。例えば、本発明の１つの例示的な実施形態によれば、探索画像内で最もマッチングするブロックは、正規化された相互相関関数を最大にし、このブロックと参照ブロックとの間の相対的なオフセットが所望の速度推定を与える。この手順は、参照位置（ｘ，ｙ）、平均強度μＲ及びμＳを有するＭ×Ｎ個のサブブロックに関して、式（７）及び（８）で数学的に表現される。重複するサブブロックは、参照画像内のより細密な格子で速度を推定するために使用される。
【００４９】
【数６】
0005524904.tif 000007
【００５０】
平行移動するインパルス散乱体の場合に、式（５）〜（８）による速度推定は、
【数７】
0005524904.tif 000008
のときに、インパルス応答のシフトを追跡することができる。非剛性の変形を受ける散乱体の集団を含む現実の組織の場合には、各散乱体からのコヒーレントインパルス応答が後方散乱信号の干渉パターンを生成し、この信号は、時間的に連続する画像間で単純に平行移動しそうにない。速度推定は、干渉「ノイズ」、つまり式（８）で最大に達成可能とされる相関係数を低減させるノイズに敏感となる。上記のような非相関効果は、相関ウィンドウ（窓）のサイズが対象となる変形構造と比べて広いときは常に生じ、また負荷によって生起される歪みが大きい場合に生じることがある。いずれの場合にも、機械的負荷の影響は、単純なスペックルの平行移動ではモデル化できず、その理由は、空間歪みが根本的な散乱体分布で生じるためである。特定の現実的な状況では、画像ノイズと非相関性は、相関表面ｘ，ｙ（ｕ，ｖ）内の真の変位における相関値を低減させるだけではなく、真の変位における相関値を超える値をもった複数の極大値又は偽ピークに加えて、相関ピークの位置をシフトさせるジッタをもたらす。
【００５１】
Ａライン間の一次元相関は、理想的には真の変位において単一の明確なピークとなるべきである。しかし、スペックル非相関性とノイズのために、一般に複数のピークが相関関数に存在し、そのうちの最も高いピークは真の変位よりも非常に小さい速度に存在する。また、二次元の運動推定の場合には、理想的な相関関数が、単一の明確なピークを示さねばならないが、複数の極大値が存在することがある。相関ウィンドウ全体に及ぶ境界等の特徴を有する画像からの二次元推定の場合には、境界の接線方向の速度成分を決定することが難しくなる。このような場合は、相関関数は明確なピークを含まず、相関値は、境界の接線方向を向いて幅広い範囲の変位に亘って高くなる。一次元及び二次元のいずれの場合にも、得られる速度推定は、血管ＯＣＥでの使用にとって極度にノイズが多い歪み推定をもたらす虞がある。
【００５２】
ロバストなコヒーレンス弾性率計測法
速度推定を改善するための方策には、ノイズ抑圧のための画像シーケンス（列）のブラーリング、広い相関ウィンドウの採用、及び相関最大化による推定後の速度場のスムージングがある。これらの方策は、ある観察に基づくと、速度及び歪みの推定をある程度改善するが、弾性率計測法に対するＯＣＴの空間分解能の利点を損なうこともある。例えば、画像シーケンスのブラーリングは、ノイズだけでなく、運動追跡に有用とされる精緻な画像特徴をも除去する場合がある。広い相関ウィンドウは、速度場の微細な変化を追跡する能力を低下させる虞があり、また、式（７）及び（８）に仮定した並進スペックルモデルの破綻につながる可能性がある。メジアンフィルタ又は他のスムージングカーネルを用いた速度又は歪みのフィルタリングは、測定が行われた後の測定値に作用する。このため、このような方法は、速度及び歪みの推定を改善するために、基礎をなす相関関数内に存在する情報を使用できない可能性がある。より好ましい推定法は、相関を最大化する手順自体の間に、データ駆動型の速度フィルタリングを可能にすることである。そのような例示的な方法は、下記に記載する変分法である。
【００５３】
速度推定問題は、変分エネルギーの最小化として提示され、これは、拍動性動脈壁内の速度場に関する従来の知識を用いて、推定にロバスト性を与える一方で、相関関数に存在する速度情報を有効に利用することを目的とする。この手法では、最大解像度データからの利用可能な全情報を維持するために、画像のスムージングを使わない。全体的な変分エネルギー関数は、以下のようになる。
【００５４】
【数８】
0005524904.tif 000009
【００５５】
このエネルギーは、未知の速度場Ｖ＝［ｕ，ｖ］に依存し、データ忠実度ＥＤ（Ｖ）、歪み場の平滑度ＥＳ（Ｖ）、及び動脈壁の非圧縮性ＥＩ（Ｖ）を調整する３つの項の重み付けされた組み合わせである。これら各項の関数形は、以下の通りである。
【００５６】
【数９】
0005524904.tif 000010
【００５７】
ここで、式（１０）の表式「ｘ，ｙ（Ｖ）」は、式（８）に示した相関係数関数である。歪み平滑度項及び組織非圧縮性項が存在しない場合において、データ忠実度項の最小化は、相関関数の最大化と同じであり、これは従来の速度測定からの推定と同一の速度推定をもたらす。歪み平滑度項及び組織非圧縮性項は、動脈の組織生体力学に関する従来の知識から逸脱しないように、速度推定を制限する。隣接する参照位置からの相関関数の情報は、非相関性、偽ピーク、及び高い相関係数値の明確に定義されない領域について、ロバスト性を与えるために、効果的に組み合わされる。歪み平滑度項は、動脈の速度場の二次導関数を、壁に亘って滑らかに変化させるのに対し、非圧縮性モデルは、ｕ及びｖの速度場の挙動と結び付くことで、壁内部の点が非圧縮性から大きく逸脱しないようにする。所望の速度場推定は、全体的な変分エネルギーを最小化する。
【００５８】
【数１０】
0005524904.tif 000011
【００５９】
エネルギー最小化問題に対する数値解を得るためには、式（１３）の連続的表現を以下のように離散化する。
【００６０】
【数１１】
0005524904.tif 000012
【００６１】
ここで、列（ｘ）方向及び行（ｙ）方向の離散的な速度成分は、順序付けられた以下の縦ベクトルとしてそれぞれに表現される。
【００６２】
【数１２】
0005524904.tif 000013
【００６３】
ここで、ｋは、対象のｋ番目の参照位置についての、順序付けの添え字であり、［ｉｋ，ｊｋ］は、参照画像行列ＩＲ［ｉ，ｊ］内のこの位置での行座標及び列座標である。
【００６４】
離散化されたデータ忠実度項は、以下の式で表される。
【００６５】
【数１３】
0005524904.tif 000014
【００６６】
ここで、離散的な領域における相関係数関数は、参照画像行列をＩＲ［ｉ，ｊ］、探索画像行列をＩＳ［ｉ，ｊ］とした場合に、以下の通りである。
【００６７】
【数１４】
0005524904.tif 000015
【００６８】
参照画像行列ＩＲ［ｉ，ｊ］及び探索画像行列ＩＳ［ｉ，ｊ］については、直線的な格子上で定義された、一定間隔を置いた点集合からサンプリングされる。実際には、全画像について秒オーダーの迅速な計算を実現するために、二次元ＦＦＴを使用する高速な正規化された相互相関近似を使用して、式（１７）の分子と、式（１７）の分母についての事前に計算される移動和を計算する。離散化された歪み平滑度項及び非圧縮性項は、それぞれ下式のようになる。
【００６９】
【数１５】
0005524904.tif 000016
【００７０】
【数１６】
0005524904.tif 000017
【００７１】
ここで、Ｄ２ｒ及びＤ２ｃは、列ベクトルｕ及びｖにおける近隣位置からの速度に作用する第二次の行差分行列及び列差分行列である。行列Ｄ１ｒ及びＤ１ｃは、これらに対応する第一次の行差分演算子及び列差分演算子である。Ｍ×Ｎの矩形状領域で定義された二次元の速度場から生成される、順序付けられた速度ベクトルの場合に、第一次の行差演算子Ｄ１ｒ及び第一次の列差演算子Ｄ１ｃは、以下のように定義される。
【００７２】
【数１７】
0005524904.tif 000018
【００７３】
ここで、Ｄ１（Ｍ−１）Ｎは（Ｍ−１）×Ｎの第一次の差分行列であり、ＩＭはＭ×Ｍの単位行列である。これに対応する第二次の行差分演算子及び列差分演算子は、それぞれ以下のようになる。
【００７４】
【数１８】
0005524904.tif 000019
【００７５】
ここで、Ｄ２（Ｍ−２）Ｎは（Ｍ−２）×Ｎの第二次の差分行列であり、ＩＭはＭ×Ｍの単位行列である。
【００７６】
式（１４）を最小化するために、第一次の変分を導出して、オイラー方程式を得る。
【００７７】
【数１９】
0005524904.tif 000020
【００７８】
ここで、データ忠実度項の第一次変分は以下のように定義される。
【００７９】
【数２０】
0005524904.tif 000021
【００８０】
式（２４）のオイラー方程式は、以下の発展方程式を作ることで反復的に解くことができる。
【００８１】
【数２１】
0005524904.tif 000022
【００８２】
ここで、Ａは以下の通りである。
【数２２】
0005524904.tif 000023
【００８３】
ｔは、反復毎にとった時間ステップである。時間ｔにおける更新された速度推定を解くためにこれを整理すると、以下の行列及びベクトルの方程式を得る。
【００８４】
【数２３】
0005524904.tif 000024
【００８５】
定常状態では時間導関数が消え、これにより得られる速度推定は式（２４）を満たす。実際には、速度場の初期推定から開始して、式（２７）の各反復においてＬＵ分解法を使って、更新された速度推定を解くことが可能である。非整数の速度推定では、双三次補間を用いて、データ忠実度項で必要な勾配を計算する。この例示的な方法は、速度場の大きさについての最大変化が０．０１％未満となるまで続行される。
【００８６】
大域的最小値への収束のための多重解像度推定
式（２７）の解は、変分エネルギー関数の極小値へと収束することとなる。よって、良好な大域的収束性を保証するために、未知の速度場は大域的最小値の近くで初期化する必要がある。これを実現するために、図５のブロック図として示した、本発明による例示的な多重解像度技術が使用可能である。例えば、入力される最大解像度の参照画像及び探索画像は先ず、ステップ２５５において１０分の１にダウンサンプリングされて、低解像度シーケンスが得られ、ステップ２６０では、それから、式（７）、（８）の相関最大化によって、速度場の低解像度初期推定が得られる。この推定は、ステップ２６５において、低解像度領域に適用される変分法の初期化に使用される。そして、ステップ２７０において、速度についてのロバストな低解像度推定が、高解像度領域にマップ化され、高解像度の探索領域を定義するために使用されるが、これは、ステップ２７０において、対象の各参照位置における最大解像度相関関数を計算するためである。また、変分法から得た低解像度推定は、最大解像度相関関数による速度場の反復的な推定にとって、良好な初期推測として役立つ。そして、ステップ２７５で得られた最大解像度の速度推定は、ステップ２８０での表示や後続の歪み計算に使用される。
【００８７】
弾性率の決定及び改良された歪みマッピング
ＯＣＴに基づく組織弾性画像化については、組織弾性率及び歪み分布を併せて推定するために、（図４に示すような）一元化された計算フレームワークを利用するが、これは、以下の基本手順から構成される。
【００８８】
１．参照ＯＣＴデータセットによる数値的再構築（ステップ２１０，２１５）。
ａ．ＯＣＴ区分化及び組織分類。
ｂ．弾性率の値及び境界条件の初期化。
ｃ．測定の圧力負荷の適用。この負荷は、研究対象の系に固有（例えば、心臓周期中の血管内圧の正常な変動）であってもよく、また外部的に制御されてもよい。
ｄ．血管変形の予測を得るための、例えば有限要素モデリング（ＦＥＭ）などの数値シミュレーション。
２．参照ＯＣＴ画像データと探索ＯＣＴ画像データとの間での、モデルの剛体並進及び回転の推定（ステップ２１５）。
３．剛体変換と、参照ＯＣＴデータセットを歪ませるための数値予測された変形場との組み合わせ（ステップ２２５）。
４．ＯＣＴ固有のデータ忠実度項の計算（ステップ２３０）。
５．弾性率の更新と、ＯＣＴデータ忠実度を最大にするための剛体変換、及び更新された数値モデルのシミュレーション（ステップ２１０）。
６．弾性率ベクトル推定値の収束後における、最終的な弾性率及び歪みの表示（ステップ２３５，２４０）。
【００８９】
例えば、取得された参照ＯＣＴデータは最初、血管壁の形状を画定するために区分化される（ステップ２１０）。勾配に基づく能動輪郭を用いて内腔境界が抽出され、これは、血管のＯＣＴ画像化において、大きな強度勾配の値を示す。好ましい実施形態では、レベル設定に基づく能動輪郭が使用される。このレベル設定法では、動脈表面が、高次元の埋め込み空間のゼロレベルセットとしてモデル化される。閉じた平面曲線Ｃ（ｐ）：［０，１］→Ｒ２の場合には、埋め込み表面（ｅｍｂｅｄｄｉｎｇ  ｓｕｒｆａｃｅ）、ｕ：→Ｒ２→Ｒは、曲線に対する符号付きの距離関数によって表される。ｕの各点に割り当てられる値は、曲線Ｃ上での最も近い点までの符号付き距離である（曲線の内側で負、外側で正である）。画像内の物体境界に最も適合する曲線を特定するために、以下の曲線の発展方程式が離散化されて、勾配降下法によって反復的に解かれる。
【００９０】
【数２４】
0005524904.tif 000025
【００９１】
ｇは画像勾配値の関数であり、κはユークリッド曲率であり、
0005524904.tif 000026
は、単位法線ベクトルである。曲率項により、大きな画像勾配が存在する場合以外は、曲線が滑らかとなり、また、画像から導出される曲線伝播項は、物体境界において曲線を大きな勾配に向けて引きつける。ｕはＣの間接的表現であるため、式（２８）を解くことは以下の式を解き、曲線Ｃを抽出するためにゼロレベルセットを選択することと同じである。
【００９２】
【数２５】
0005524904.tif 000027
【００９３】
式（２９）のレベルセット発展方程式は位相幾何学的に柔軟であり、血管分岐の複雑な変化を自動的に受け入れる。エントロピーを満たす上流有限差分法は、数値の安定性及び速さをもって（２）を解くために、上述した狭帯域の更新法と共に使用される。
【００９４】
ＯＣＴデータにおける血管壁の最も外側の境界は、強度閾値処理及び計算幾何学に基づいて定義することができる。画像化データには、ＯＣＴシステムの測定されたノイズフロアを超える信号強度を有する全ての点を探すために、最初に閾値が設定される。そして当業者には理解されるように、これらの点から幾何学的な凸包が形成されて、得られた表面が外側の血管境界を画定するために使用される。
【００９５】
次に、内側及び外側の血管表面に含まれるＯＣＴボクセルが、管理された最大事後（Ｍａｘｉｍｕｍ  Ａ  Ｐｏｓｔｅｒｉｏｒｉ：ＭＡＰ）分類、及び実験的に導出されたクラス条件付き（ｃｌａｓｓｏ−ｃｏｎｄｉｔｉｏｎａｌ）の強度確率密度関数に基づいて、脂質に富む組織、繊維性組織又は、石灰化組織の領域に割り当てられる。本発明の例示的な一実施形態では、ＯＣＴ及び組織学からのトレーニングデータが記録され、脂質に富む（Ｌ）組織、繊維性（Ｆ）組織及び石灰化（Ｃ）組織の領域が、ＯＣＴデータにおいて、これに対応する組織学データに基づいて位置される。各組織クラスに対して、ＯＣＴ画像強度値が抽出されて、クラス条件付きの強度確率分布を近似する頻度ヒストグラムが生成される。この手順により、ＯＣＴ強度Ｉに対して、任意のボクセルＰ（Ｉ｜Ｆ）、Ｐ（Ｉ｜Ｌ）及びＰ（Ｉ｜Ｃ）における３つのクラス条件付き確率分布が得られる。ベイズ式を使用すると、ボクセルが組織クラスＦ、Ｌ又はＣに属する確率は、その強度値をｇとして、下式で与えられる。
【００９６】
【数２６】
0005524904.tif 000028
【００９７】
先の組織クラスの確率Ｐ（Ｆ）、Ｐ（Ｌ）、Ｐ（Ｃ）は、各病変型の頻度に関する生体内の観察から導出することができ、これらの観察は各ボクセルについてもあてはまると仮定される。クラス条件付きの確率値から、任意の組織クラスを、あるボクセルに割り当てることに関連する分類誤差の確率が計算される。分類誤差の確率を最小化する組織割当は、ボクセル組織クラスに対して選択される。
【００９８】
区分化と分類を行った後で、血管形状、及び、同様の弾性を有するプラーク内ボクセル間の境界面を規定する表面を得ることができる。そして技術的に理解されるように、これらの表面は、市販されるか又はカスタムコード化された有限要素分析プログラムにおいて、有限要素モデルを作成するために使用される。このプログラム又は当業界で公知の代替のメッシュ生成手順を用いて、血管及びプラーク内成分の有限要素メッシュが生成される。境界条件については、動脈入口における内腔輪郭の固定点と、径方向に自由に移動できる内腔重心の反対側の点を定義することによって適用される。測定された圧力負荷は、血管モデルの内腔表面に適用されて、（対応する画像領域における）ＯＣＴ組織分類の結果と、血管の生体力学に関する文献から得た脂質組織、繊維性組織及び石灰化組織の平均弾性率値に基づいて、各メッシュ要素に対して初期弾性率推定値が割り当てられる。
【００９９】
本発明による技術の１つの例示的な実施において、各メッシュは、入力ＯＣＴデータが二次元画像であるか、三次元ボリュームであるかに応じて、二次元の四辺形又は三次元のアイソパラメトリックな線形弾性有限要素のいずれかより構成される。変位（ｕ）に基づく要素又は変位／圧力（ｕ／ｐ）に基づく要素は、ポアソン比０．５の手法での「ロッキング」の回避に使用できる。しかし、ＯＣＴによる組織弾性画像化のための一元的なフレームワークについては、このような特定の構造要素又は構成的材料モデルの使用に限定されない。例えば、非線形Ｍｏｏｎｅｙ−  Ｒｉｖｌｉｎ歪みエネルギー関数を、本実施形態において使用する線形材料モデルに置き換えことができる。モデル構築後に、有限要素シミュレーションを行って、探索ＯＣＴデータセットにおける変形メッシュのノード位置が予測される。
【０１００】
剛体の並進及び回転は、参照画像と探索画像との間で血管内に生じ、上記したＦＥＭモデルは、適用された境界条件のために、これを反映しない。よって、当業界で公知のように、多次元画像相関法のための標準的な方法を用いて、参照データセットと探索データセットとの間で、有限要素モデルにおける剛体の並進及び回転が計算される。そして、推定された剛体の並進がＦＥＭで予測されるメッシュ変形に適用されて、参照データセットと探索データセットとの間で有限要素モデルの全体的な挙動が得られる（ステップ２１５を参照）。
【０１０１】
ＦＥＭのノード変位は、剛体運動を受けるときにメッシュ生成し、変形は、参照データの空間座標を、探索データセットの予測位置に対してマッピングにより位置付ける歪み場を規定する。この歪み場は、多次元補間を用いてＯＣＴ参照画像化データを空間的に再サンプリングする（２２５）ために使用される。好ましい実施形態では、線形補間カーネルがその単純性及び計算速度を理由に使用される。しかし、原理上は、これに代えて、スプライン補間又は三次補間等、任意の多次元補間法が使用可能である。
【０１０２】
結果としての歪められた参照データは、探索データ及びＯＣＴドップラ測定値と共に、以下の項の線形結合によって定義されるＯＣＴ特有の目的関数を計算するために使用される（ステップ２３０を参照）。
【０１０３】
１．歪められた参照画像と探索画像中の測定されたＯＣＴデータの強度間での二乗差についての負の和。
２．歪められた参照画像と探索画像中の測定されたＯＣＴデータの強度間での正規化された相互情報量。
３．歪められた参照画像と探索画像中の測定されたＯＣＴデータの強度間での相関係数の要素に亘る総和。
４．加えられた圧力／変位負荷による、組織又は境界の造影剤の光学特性についてモデル予測の変化と測定された変化との間の二乗誤差のノードに亘る総和。
５．画像化ビーム方向における、光学的ドップラ速度測定値と、ＦＥＭ予測の変位との間の内積のノードに亘る総和。
【０１０４】
この例示的な目的関数については、最初に弾性率及び剛体変換パラメータを更新し、次に、参照フレームと探索フレームとの間の血管変形の改善された予測を得るために、各反復においてＦＥＭモデルを繰り返し実行することによって、反復的に最大化される。当業者には理解されるように、全体的な目的関数を最大化にするためには、多次元の制限付き共役勾配法を使用することができる。この例示的な方法による制限は、各要素のとりうる弾性率値の範囲を制限するために課され、これは、生体力学的挙動（生体力学試験から先験的に測定される）、つまり、ＯＣＴ画像強度に基づいて要素に割り当てられる組織クラスでの典型的な挙動を基礎とする。制限された関数を最大化するためのこの手順は、代用可能な多くの可能な技術のうちの１つである。本発明の１つの例示的な実装では、最大弾性率値の絶対変化が０．０００１を下回る場合に、推定手順が完了したとみなせる。収束時には、最終的な弾性率分布と、これに対応する有限要素法シミュレーションからの、対応する歪み及び応力の場について、入力データの次元に応じて、色分けによりマップ化された画像又はボリュームレンダリングとして表示することができる（２３５，２４０）。
【０１０５】
自動レンジ調節
本発明はまた、自動レンジ調節技術を使用することができ、これには、２００２年４月３０日に出願された同時係属中の特許文献４に説明する処理法が含まれ、この出願の全ての開示が参照によりここに援用される。
【０１０６】
１つの例示的な実施形態では、自動レンジ調節メカニズムにより、本発明の技術を血管画像化に適用できるようになり、これにより、画像カテーテルを血管内腔の中心に位置させる必要がなくなる。この例示的な実施形態では、自動レンジ調節メカニズムのフィードバック信号を本発明の画像化メカニズムに組み込むことが望ましく、これは、例えば、血管構造の精度の高い決定を容易にし、かつ維持するためである。
【０１０７】
実施例
以下の説明は、本発明による方法の例示的な実施形態について、実験的な検査に関する詳説である。具体的には、例示的な多重解像度変分法が、円形の内包物を含む組織ブロックの軸方向における圧縮中に、シミュレーションによるＯＣＴ画像化において実行された。図６は、組織ブロック及び円形の内包物のために使用した、有限要素モデルの形状３００、及び、これに対応する有限要素メッシュ３０５を示す。連続する干渉画像は、式（１）〜（６）で説明したように、指数関数的な減衰項と、コヒーレントＯＣＴ点広がり関数と試料内を移動する点状散乱体から生じる後方散乱分布との間の畳み込みと、の積を計算することによって生成された。
【０１０８】
組織ブロック中の離散的な点における後方散乱値については、ブロック中の散乱体に対する自由度１０と、円形の内包物中の散乱体に対する自由度２でもって、独立した一様な乱数変数としてシミュレーションが行われた。これらの値は、内包物に対するブロック内の平均後方散乱値を大きくするために経験的に選ばれた。得られたコントラスト差は、ＯＣＴ画像内の脂質と正常な動脈壁との間で観察されたコントラスト差を模倣する。
【０１０９】
圧縮中における組織の散乱体の運動については、組織形状の有限要素モデリングから、変位場を使用してシミュレーションが行われた。円形の内包物３１５が組み込まれた二次元の矩形状をした幾何形状３１０が定義された。全シミュレーションにおいて内包物の直径は５００μｍ（マイクロメートル）であった。
【０１１０】
不動点境界条件がブロック３２０の下端中央に課され、緩い境界条件が下端３２５における他の全ての点に与えられた。軸方向（下流）の変位負荷が上面３３０に加えられ、５回の時間ステップでブロックが０．１５ｍｍだけ圧縮されて、全体的な歪みが４．３％となった。各時間ステップでは、負荷により、ブロック内に約０．９％の歪みが生起された。
【０１１１】
ブロックと内包物はともに、ほぼ非圧縮の線形弾性材料（ν＝０．４９５）としてモデル化された。全てのシミュレーションについて、ブロックには単位弾性率が割り当てられる一方で、内包物の弾性率は、線維組織に埋め込まれた脂質に富む病巣又は石灰化の病巣を表現するように変更された。線維プラークに対する脂質の係数比は約０．０００１で、線維プラークに対するカルシウムの係数比は約５である。有限要素モデリングは、ＡＤＩＮＡ  ８．０（マサチューセッツ州ウォータータウン）を使用して、９ノードの四辺形での二次元平面歪要素からなるメッシュ３３５，３４０を用いて実施された。メッシュ密度は、各要素のエッジの長さによって定義され、内包物３４０の内外では０．０２５ｍｍ、周辺ブロック３３５では０．１ｍｍであった。各シミュレーションモデルは、約３２００の要素及び１３０００のノードから構成された。各時間ステップで計算される変位場は、連続するＯＣＴフレーム間での、組織の散乱体の速度ｕ（ｘ，ｙ）及びｖ（ｘ，ｙ）を表現するために使用された。後方散乱場σｂ（ｘ，ｙ，ｔ＋１）は、時刻「ｔ＋１」での移動する組織散乱体に関連しており、下式を使用して計算された。
σｂ（ｘ，ｙ，ｔ＋１）＝σｂ（ｘ−ｕ（ｘ，ｙ），ｙ−ｖ（ｘ，ｙ），ｔ）  （２８）
【０１１２】
実際には、第１フレームにおける組織散乱体の場が、最初にアップサンプリングされて、次に式（２８）の線形補間を用いて不均一に再サンプリングされることで、後続フレーム中の組織散乱体の場が得られる。
【０１１３】
式（２）〜（４）におけるコヒーレントＰＳＦ近似を用いることにより、本ＯＣＴシステムの軸方向のフリンジパターンの測定に基づいて、ＰＳＦのシミュレーションが可能である。図７のグラフ３５０は、シミュレーションを行ったＯＣＴ点広がり関数（ＰＳＦ）３６５に対する、フリンジ分解ＯＣＴ  ＰＳＦ３６０の測定を示す。シミュレーションと測定は、ある一致した特徴を有する。このコヒーレントＰＳＦについて、式（２８）からの組織散乱体の場との畳み込みを行った後で、軸方向の指数関数的な減衰との乗算が行われ、ダウンサンプリングが行われることで、ピクセルサイズが１μｍ×２５μｍとされる一連のシミュレーション画像Ｉ（ｘ，ｙ，ｔ）が得られる。このようなシミュレーションに存在するスペックル又は乗算ノイズの度合いは、以下のノイズモデルを使用することより、更に制御された。
Ｉｎ（ｘ，ｙ，ｔ）＝Ｉ（ｘ，ｙ，ｔ）＋ｎＩ（ｘ，ｙ，ｔ）    （２９）
【０１１４】
ここで、ｎは、平均値がゼロで分散がσｎ２とされた一様分布の乱数変数である。図７において、３５５は、図６の３００、及びＯＣＴ  ＰＳＦ３６５におけるＦＥＭ形状から生成されたシミュレーションＯＣＴ画像の復調結果を示す。内包物３７０が、組織ブロック３７５内に見え、正常な動脈壁中の脂質プールに近い外観を有する。
【０１１５】
上記の方法は、迅速な処理のためにＣ＋＋プログラミング言語を使用して実装することができ、可視化及び分析はＭＡＴＬＡＢ（米国マサチューセッツ州ナティック所在ＭａｔｈＷｏｒｋｓ）で実行することができる。数秒程度という高速で全画像の計算を実現するために、式（８）の相関係数関数に対して高速な正規化された相互相関近似を使用したが、これには、この式の分子、及び分母について事前に計算される移動和を計算するために二次元ＦＦＴを使用する。
【０１１６】
相関の最大化に基づく従来の速度推定技術の場合、各相関関数は、ピーク周辺で双三次補間を用いて５０倍のアップサンプリングが行われたが、これは軸方向で０．０２μｍ、横方向で０．５μｍのサブピクセルシフトを検出するためである。従来の方法では、相関関数を計算するために、８１×７ピクセル（８１×１７５μｍ）の参照ブロックサイズと３６１×２１ピクセル（３６１×５２５μｍ）の探索領域が使用された。これらのパラメータは、速度推定の誤差を最小化する要求に対して、速度の空間変化の検出感度に対する要求を均衡させるために経験的に決定された。従来の運動追跡からの速度及び歪み推定についてのメジアンフィルタリングは、本発明による例示的な変分法から得た結果と比較する前に、５×５のカーネルで実行された。
【０１１７】
変分法の例示的な実施形態における低解像度ステップでは、各相関関数を計算するために、１５×１１ピクセル（７５×１３７５μｍ）の参照ブロックサイズと、６１×４１ピクセル（３０５×５１２５μｍ）の探索領域が使用された。最大解像度では、２５×７ピクセル（２５×１７５μｍ）の参照ブロックサイズと、１０１×２１ピクセル（１０１×５２５μｍ）の探索領域が使用された。変分法の重み付けパラメータは、シミュレーションされたＯＣＴシーケンスのサブセットに基づいて、経験的にａ＝１、ｂ＝２０、ｃ＝０．１に決定された。そして、これらの値は後述するＯＣＴ実験の速度推定に使用された。
【０１１８】
速度推定の後、変形行列Ｆが、以下の表式に対する第一次の有限差分近似を使用して速度場の全ての点で計算された。
【０１１９】
【数２７】
0005524904.tif 000029
【０１２０】
歪みを微小とする近似を仮定した場合に、変形行列は、以下の関係によって歪み行列Ｅ及び単位行列Ｉに関連付けられる。
【０１２１】
【数２８】
0005524904.tif 000030
【０１２２】
軸方向圧縮のシミュレーションＯＣＴ画像化では、歪み行列の軸方向における歪み成分ｙｙの結果を示し、また横方向の伸張における画像化実験では、横方向の歪み成分ｘｘの結果を示す。速度及び歪み場の推定の誤差は、正規化された二乗平均平方根の誤差測定に基づいて評価された。
【数２９】
0005524904.tif 000031
【０１２３】
ここでＮは推定値の合計数、ｖｋは軸方向又は横方向の推定速度、ｖｋ，ｒｅａｌは軸方向又は横方向の実際の速度、ｋは速度場中のｋ番目の対象点での軸方向又は横方向の推定歪み、ｋ，ｒｅａｌは、速度場中のｋ番目の対象点での軸方向又は横方向の実際の歪みである。
【０１２４】
シミュレーション実験
５００μｍの内包物の検出能力が、スペックル非相関性の関数として検査され、ロバストな推定から得た値に対して、従来の運動推定から得た軸方向の速度及び歪みの推定におけるＲＭＳ（二乗平均平方根）誤差が比較された。ｖｋ，ｒｅａｌ及びｋ，ｒｅａｌの値は、有限要素モデリングから直接得た。速度は、下方及び右方への変位にそれぞれ対応した、正の軸方向速度及び正の横方向速度についてピクセル単位で報告された。
【０１２５】
図８、図９はそれぞれ、弾力的な内包物と固い内包物の運動シーケンスにそれぞれ対応するＦＥＭによる軸方向速度場４００，４５０を示す。また、これらは、従来の運動追跡４０５，４５５と、変分法フレームワーク４１０，４６０内でのロバストな推定から得た速度測定結果の軸方向成分を示す。これらの例では、内包物の位置は、符号４１５，４２０，４２５，４６５，４７０，及び４７５に対応している。この２つの方法から推定された速度場は、ＦＥＭから得た真の速度と定性的に同じである。変分法による推定値は、従来の追跡よりも著しく滑らかとなる。更に、軸方向速度のＲＭＳ誤差は、ロバストな推定よりも従来の速度測定による方が大きい。固い内包物の場合、従来の速度測定による軸方向速度場のＲＭＳ誤差が１．６０％であったのに対し、ロバストな推定のＲＭＳ誤差は１．０４％であった。同様の結果は、弾力的な内包物についても得られた（従来法の１．８３％に対し、変分法では１．４０％）。
【０１２６】
図１０、図１１は、弾力的な内包物と堅い内包物のそれぞれのシミュレーションに対応する軸方向の歪み場を示す。フレーム５００，５５０は、有限要素モデリングによる真の軸歪みを示し、フレーム５０５，５５５は、従来の運動追跡による軸方向歪みの推定を示し、フレーム５１０，５６０は、多重解像度変分法による軸方向歪みの推定を示す。これらの画像における内包物の位置は、符号５１５，５２０，５２５，５６５，５７０，及び５７５に対応している。
【０１２７】
歪み場測定に関する速度ノイズの影響は、これらの例から明らかである。従来のロバストな速度測定による推定は、ＦＥＭで導出される変位と定性的に同様であるが、歪み推定は、真の歪に対して同様ではない。軸方向歪みは軸方向速度場の導関数であり、このため、速度推定に存在するノイズは、導関数演算子のハイパス（高域通過）特性のため、歪み画像内で強調される。メジアンフィルタリングを用いた場合であっても、従来の運動追跡から導出される歪みについて、内包物内の歪みの差の程度と大きさを充分に評価することは困難である。これとは対照的に、ロバストな歪み推定を用いれば、内包物の位置及びサイズを決定することが可能である。更に、変分法によって、弾力的な内包物と固い内包物とを視覚的に区別することが可能である。定量的には、従来法と変分法との間のＲＭＳ歪み誤差の差は顕著である。例えば、固い内包物の場合、従来の手法ではＲＭＳ歪み誤差が１０９．１％であるが、変分法ではＲＭＳ歪み誤差が２７．５％である。従来の追跡から得られる歪みは、解釈が困難であるのに対し、ロバストな推定では、病変型の間での比較について、歪み測定をより簡単にかつ高精度をもって解釈可能となる。
【０１２８】
画像化実験
中心波長１３１０ｎｍ及び帯域幅７０ｎｍの広帯域光源の光が、干渉計内で参照場とサンプル場に分かれる。サンプル場は、参照アームの光路長に対応する深さで組織を調べるために、走査光学装置によって合焦される。サンプルアームからの戻り後方散乱光が参照場と混合されて干渉信号が生成され、該信号は、２５μｍ（横方向）×１μｍ（軸方向）のピクセルを生成するためにデジタル化される。干渉信号の振幅は、参照アームに規定される走査深さにおける組織構造及び光学特性に関する情報を含んでいる。ＸＹ断面における組織構造については、参照アームの光路長を変更することにより軸（Ｙ）方向で探索され、検体を横切るサンプルビームを掃引することによって横（Ｘ）方向で探索される。この例示的な方法では、横方向５００ピクセル×軸方向２５００ピクセルからなる画像フレームが、２５０ミリ秒で取得される。
【０１２９】
検視解剖で得たヒトの大動脈の正常断片は、リン酸緩衝生理食塩水中で３７℃に暖められた。画像化は、断片を取り出して２４時間以内に実施された。円筒状の大動脈の断片が、縦方向に切断されて開かれ、管腔の表面を検体光束に露出させる矩形の組織検体が得られた。横方向の機械的負荷が、無処置の大動脈断片の周縁部における伸張に近くなるように、縦方向に切断された端が検体ホルダに固定された。検体は検体ホルダ内に水平に取り付けられ、これにより、その一端がしっかりと固定され、他端がマイクロマニピュレータ、つまり、水平軸に沿って一次元の平行移動ができる極微操作装置に固定された。試料内の画像化位置は、検体光束と同期する照準ビーム（レーザダイオード、６３５ｎｍ）を視覚化することによって監視された。走査位置は検体中心に設定され、サンプルアーム内の走査方向は、散乱体の変位が画像化面内に制限されるように位置合わせされた。
【０１３０】
この画像化実験において、速度及び歪みの推定が、静止又は横方向の伸張のいずれかの状態下にある正常な大動脈の検体で調べられた。静止状態の場合、グランドトゥルースの速度場及び歪み場はともにゼロとすべきである。このため、速度及び歪みの推定の、ゼロからの偏差が１つの性能評価基準として決定された。横方向の伸張の場合、予測された歪み分布は均一であり、推定された歪みの標準偏差を、試料内の歪みの均一性についての評価基準として報告する。
【０１３１】
静止状態では、従来法で推定された平均横方向速度が、フレーム当たり０．１７μｍであり、これは、名目上用いられる並進量、フレーム当たり０ピクセルと比較される。場についての横方向速度推定の標準偏差は、フレーム当たり３．９４μｍであった。軸方向では、平均偏差及び標準偏差がそれぞれ０．６６μｍ及び４．２８μｍであった。変分法による速度推定は、横方向でフレーム当たり０．００１７μｍ（軸方向ではフレーム当たり０．５４μｍ）の平均速度を示し、これは、従来法の結果よりも、実際の並進量に近かった。また、横方向でフレーム当たり０．１６μｍ（軸方向ではフレーム当たり０．５８μｍ）の速度標準偏差は従来の方法の結果よりも小さかった。
【０１３２】
場の歪み推定の平均偏差及び標準偏差は、従来法の場合に横方向ではそれぞれ−０．０４３％と５．５％（軸方向では−０．２４％と１１．２６％）であったが、これは場の名目上用いられる歪み０％と比較される。これに対して、ロバストな推定では、一般に、より良好な歪み推定性能を示す。歪み推定の平均偏差及び標準偏差は、横方向ではそれぞれ−０．００３８％と０．１２％（軸方向では−０．０２４％と０．１２％）であった。
【０１３３】
わずかな横方向の伸張による横方向の移動を受けた大動脈断片についての例示的なＯＣＴ画像を、図１２に６００で示す。符号６０５で示した領域がマスクされており、これは大動脈の組織６１０の表面より上の空気に対応している。図１３は、速度場測定の横方向成分を示す。フレーム６５０は、従来の速度推定手法によって得られ、またフレーム６５５は、本発明による変分法の例示的な実施形態を使用している。組織より上の空気に対応する領域６６０，６６５がマスクされ、領域６７０，６７５は大動脈組織に対応している。図１４は、推定された歪み場の横方向成分を示す。フレーム７００は従来法によって得られ、またフレーム７０５は、本発明による変分法の例示的な実施形態から得られた。組織より上の空気に対応する領域７１０，７１５がマスクされ、領域７２０，７２５は大動脈組織に対応する。
【０１３４】
ロバストな推定から得た速度は、従来の推定から得た推定値よりも、横方向及び軸方向において滑らかに変化することが認められた。従来法では、歪みの平均偏差及び歪みの標準偏差は、横方向で０．２２％、４７．２６％であり、軸方向で−３．４６％と８５．０８％であった。これに対して、ロバストな推定による歪み場は、平均偏差と標準偏差が、横方向では−０．０２６％と１．３２％であり、軸方向では−０．０２％と０．１４％であることを示した。これらの結果は、ロバストな歪み計測から得た推定値が、従来の歪み計測から得た推定値よりも、予測される均一な歪み場に対して、より合致していることを示している。
【０１３５】
本発明の更に別の例示的な実施形態によれば、例示的なフレームワークが、効率的なパラメータ削減策に使用可能な生体力学特性のＯＣＴベースの推定に提供される。このような方策は、基礎データに存在し得る機械的特性の急峻な勾配を維持しつつ、迅速な手術中での推定を支援することができる。通常、パラメータを減らすことと空間分解能とはトレードオフ（二律背反）の関係にある。組織領域同士を正しく判別するためにＯＣＴの軟組織のコントラストを使用するとともに、同種の組織が類似した力学的挙動を示すという知見によって、このトレードオフを低減させるか又は最小にし、計算性能を改善することが好ましい。これを実現するための本発明による方法の１つの例示的な実施形態が、図１５に示されており、以下の通りである。
【０１３６】
ａ．取得したＯＣＴ参照フレームを主要な組織領域に区分化する（ステップ８０５）。区分化は、同様の平均強度、強度分散、Ａラインの深さの関数としての勾配又は他の候補となる強度計量を用いたクラスタリング領域に基づいて自動化することができる。異なる領域の数は、推定すべき未知パラメータの数（例えば、約２〜６個の未知パラメータ）であり得る。
【０１３７】
ｂ．ステップ８１０において、特定された組織の形状及び境界が、粗い有限要素モデルメッシュを生成するために使用される（ステップ８１０）。組織境界は、メッシュ生成ステップへの入力として使用されるため、要素が１つ以上の組織型に跨ってはならない。各要素は１つの組織型に属し、該要素はその組織についての対応する未知の弾性率に割り当てられる。
【０１３８】
ｃ．メッシュの変位境界条件を推定する（ステップ８１５）。
【０１３９】
ｄ．ステップ８１５で得た境界条件を使用することによって、所望の目的関数を最小化する値を特定するために、上記非特許文献６乃至８に記載されるように、低減された弾性率パラメータの組に亘って反復される最適化法で前進型有限要素法シミュレーションを拡張する（ステップ８２０）。
【０１４０】
ｅ．各要素において目的関数の値を決定する（ステップ８２５）。目的関数の値が所定の許容閾値を超える要素は、有限要素モデルがデータに十分に適合していない位置に対応する。このような要素は、更なる精緻化のためのリストを形成することができる。
【０１４１】
ｆ．精緻化リストの要素を更に分割して、元の要素が占めていた空間に及ぶ、より小さな新要素を形成する（ステップ８３０）。このような精緻化された要素はそれぞれ、元の要素組織クラスについて推定された弾性率に初期化される、独立した弾性率の値をもつことができる。
【０１４２】
ｇ．元の組織クラスの割当に特有とされる限定範囲に亘って値が変化するように、精緻化された要素に対して制限条件付き最適化法を用い、ステップ８２０〜８２５を繰り返す（ステップ８３５）。
【０１４３】
ｈ．収束時に、最終的な歪み及び弾性率の分布を表示する（ステップ８４０）。
【０１４４】
メッシュの精緻化のための閾値を適切に調整することによって、未知パラメータの数を、従来法の最大次元数、つまり全要素が独立の未知数として処理される次元数から大幅に減らすことができる。
【０１４５】
上記説明は、本発明の原理を示したに過ぎない。記載した実施形態を様々に変更及び変形することは、ここに記載した教示を鑑みれば、当業者にとって明らかであろう。例えば、ここに記載した本発明は、２００３年１０月２７日出願の特許文献３及び２００３年１月２４日出願の特許文献５に記載した例示的な方法、システム及び装置において有用である。これら出願の開示は、その全体をここに参照することで援用される。よって、当業者には分かるように、ここに明確に図示又は記載されない、数多くのシステム、構成及び方法を工夫し、本発明の原理を具体化できるが、これらは本発明の趣旨及び範囲に含まれる。また、上記に示した全ての刊行物、特許及び特許出願は、その全体が参照によりここに援用される。
【図１】
0005524904.tif 000032
【図２】
0005524904.tif 000033
【図３】
0005524904.tif 000034
【図４】
0005524904.tif 000035
【図５】
0005524904.tif 000036
【図６】
0005524904.tif 000037
【図７】
0005524904.tif 000038
【図８】
0005524904.tif 000039
【図９】
0005524904.tif 000040
【図１０】
0005524904.tif 000041
【図１１】
0005524904.tif 000042
【図１２】
0005524904.tif 000043
【図１３】
0005524904.tif 000044
【図１４】
0005524904.tif 000045
【図１５】
0005524904.tif 000046