(57)【要約】線形方程式系の解の品質を大きく向上させる方法である。線形方程式系の解を向上させるため、（１）対象媒体を複数の要素にモデル化し、前記対象媒体に少なくとも１つの局部的変動を加え、（２）少なくとも１つの局部的変動に起因する出力を測定し、（３）測定した出力を処理して結果を再構築し、補正フィルタを決定し、当該補正フィルタを前記結果に適用する。【選択図】図１
2004537731.tif 000002
【特許請求の範囲】
【請求項１】
第１対象媒体を複数の体積要素に分割するステップと、
前記体積要素の光学係数の少なくとも１つに変調周波数を割り当てるステップと、
少なくとも１つのエネルギ源から１期間にわたり前記第１対象媒体にエネルギを印加するステップと、
前記第１対象媒体が放出するエネルギを少なくとも１つの検出器を介して測定するステップと、
前記第１対象媒体の測定放出エネルギを処理して少なくとも１つの画像を再構築するステップと、
周波数符号化空間フィルタ（ＦＥＳＦ）を決定するステップと、
前記第１対象媒体の少なくとも１つの再構築画像に前記ＦＥＳＦを適用するステップと、
を備える散乱媒体の再構築画像を向上させる方法。
【請求項２】
前記変調周波数は、体積要素の吸収係数に割り当てる請求項１記載の方法。
【請求項３】
前記変調周波数は、体積要素の散乱係数に割り当てる請求項１記載の方法。
【請求項４】
前記第１対象媒体の測定放出エネルギは、摂動法を使用して処理する請求項１記載の方法。
【請求項５】
前記使用する摂動法は、順問題の解を用いて断層撮影画像を再構築する請求項４記載の方法。
【請求項６】
前記使用する摂動法は、逆問題を用いて断層撮影画像を再構築する請求項４記載の方法。
【請求項７】
ＦＥＳＦを決定するステップは、
前記再構築した断層撮影画像の時間離散フーリエ変換を計算するステップと、
前記計算した時間離散フーリエ変換を処理し、対応する体積要素に関する変調周波数における振幅を決定するステップと、
を備える請求項１記載の方法。
【請求項８】
前記ＦＥＳＦは、単純行列乗算を実行することにより、少なくとも１つの再構築画像に適用する請求項１記載の方法。
【請求項９】
少なくとも１つのエネルギ源から１期間にわたり第２対象媒体にエネルギを印加するステップと、
前記第２対象媒体が放出するエネルギを少なくとも１つの検出器を介して測定するステップと、
前記第２対象媒体の測定放出エネルギを処理して少なくとも１つの画像を再構築するステップと、
前記第１対象媒体から決定した前記ＦＥＳＦを前記第２対象媒体の少なくとも１つの再構築画像に適用するステップと、
をさらに備える請求項１記載の方法。
【請求項１０】
第１対象媒体を複数の体積要素に分割する手段と、
前記体積要素の光学係数の少なくとも１つに変調周波数を割り当てる手段と、
少なくとも１つのエネルギ源から１期間にわたり前記第１対象媒体にエネルギを印加する手段と、
前記第１対象媒体が放出するエネルギを少なくとも１つの検出器を介して測定する手段と、
前記第１対象媒体の測定放出エネルギを処理して少なくとも１つの画像を再構築する手段と、
ＦＥＳＦを決定する手段と、
前記第１対象媒体の少なくとも１つの再構築画像に前記ＦＥＳＦを適用する手段と、
を備える散乱媒体の再構築画像を向上させるシステム。
【請求項１１】
前記変調周波数は、体積要素の吸収係数に割り当てる請求項１０記載の方法。
【請求項１２】
前記変調周波数は、体積要素の散乱係数に割り当てる請求項１０記載の方法。
【請求項１３】
前記第１対象媒体の測定放出エネルギは、摂動法を使用して処理する請求項１０記載の方法。
【請求項１４】
前記使用する摂動法は、順問題の解を用いて断層撮影画像を再構築する請求項１３記載の方法。
【請求項１５】
前記使用する摂動法は、逆問題を用いて断層撮影画像を再構築する請求項１３記載の方法。
【請求項１６】
ＦＥＳＦを決定するステップは、
前記再構築した断層撮影画像の時間離散フーリエ変換を計算する手段と、
前記計算した時間離散フーリエ変換を処理し、対応する体積要素に関する変調周波数における振幅を決定する手段と、
を備える請求項１０記載の方法。
【請求項１７】
前記ＦＥＳＦは、単純行列乗算を実行することにより、少なくとも１つの再構築画像に適用する請求項１０記載の方法。
【請求項１８】
少なくとも１つのエネルギ源から１期間にわたり第２対象媒体にエネルギを印加する手段と、
前記第２対象媒体が放出するエネルギを少なくとも１つの検出器を介して測定する手段と、
前記第２対象媒体の測定放出エネルギを処理して少なくとも１つの画像を再構築する手段と、
前記第１対象媒体から決定した前記ＦＥＳＦを前記第２対象媒体の少なくとも１つの再構築画像に適用する手段と、
をさらに備える請求項１０記載の方法。
【請求項１９】
プロセッサにより実行されると、第１対象媒体を複数の体積要素に分割する命令コードと、
プロセッサにより実行されると、前記体積要素の光学係数の少なくとも１つに変調周波数を割り当てる命令コードと、
プロセッサにより実行されると、少なくとも１つのエネルギ源から１期間にわたり前記第１対象媒体にエネルギを印加する命令コードと、
プロセッサにより実行されると、前記第１対象媒体が放出するエネルギを少なくとも１つの検出器を介して測定する命令コードと、
プロセッサにより実行されると、前記第１対象媒体の測定放出エネルギを処理して少なくとも１つの画像を再構築する命令コードと、
プロセッサにより実行されると、ＦＥＳＦを決定する命令コードと、
プロセッサにより実行されると、前記第１対象媒体の少なくとも１つの再構築画像に前記ＦＥＳＦを適用する命令コードと、
を備えるコンピュータ読み取り可能媒体に格納されプロセッサによって実行可能なプログラム。
【請求項２０】
プロセッサにより実行されると、少なくとも１つのエネルギ源から１期間にわたり第２対象媒体にエネルギを印加する命令コードと、
プロセッサにより実行されると、前記第２対象媒体が放出するエネルギを少なくとも１つの検出器を介して測定する命令コードと、
プロセッサにより実行されると、前記第２対象媒体の測定放出エネルギを処理して少なくとも１つの画像を再構築する命令コードと、
プロセッサにより実行されると、前記第１対象媒体から決定した前記ＦＥＳＦを前記第２対象媒体の少なくとも１つの再構築画像に適用する命令コードと、
をさらに備える請求項１９記載のプログラム。
【請求項２１】
第１対象媒体を複数の要素にモデル化するステップと、
前記対象媒体に少なくとも１つの局部変動を加えるステップと、
少なくとも１つの局部変動に起因する出力を測定するステップと、
前記測定した出力を処理して結果を再構築するステップと、
補正フィルタを決定するステップと、
前記補正フィルタを前記結果に適用するステップと、
を備える線形方程式系の解を向上させる方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ＮＩＨ援助番号ＲＯ１−ＣＡ６６１８４の元で米国政府の支援によりなされた。米国政府は本発明に幾分の権利を有する。
【０００２】
この出願は、米国特許法１２０条に基づき、２００１年８月２日に先行出願した米国予備特許出願第６０／３０９，５７２号「散乱媒体の画像化品質を向上させるための周波数符号化空間フィルタリング方法」の利益を主張する。
【０００３】
技術分野
本発明は、線形方程式の分野に関し、特に線形方程式系、例えば散乱媒体の画像化に使用する類の方程式について、その解を向上させるため補正フィルタを使用することに関する。
【０００４】
背景技術
一般に散乱媒体の画像化は、散乱媒体内部の特性（吸収係数または散乱係数など）の空間分布画像を作成する様式に関する。これは、当該媒体にエネルギを導入し、その媒体が放出する散乱エネルギを検出することによって行う。この種のシステムおよび方法と対照的であるのは、エックス線などの投影画像化システムである。例えばエックス線システムは、エックス線エネルギ源と検出器との間の直線路を移動するエネルギの減衰または吸収を測定し画像化するものであり、散乱エネルギではない。エネルギが基本的に高分散であるか、あるいは基本的に直線路を移動するかは、エネルギの波長およびそれが移動する媒体の関数である。
【０００５】
散乱技術に基づく画像化は、人体、地層、大気などの特性を画像化することに関し、新しいエネルギ波長の使用を可能にする。これら特性は、投影技法および投影波長によっては画像化できない。例えばエックス線投影技術は、骨構造等の密度の高い対象を画像化することに優れているが、血液酸化レベルの識別および画像化には比較的効果がない。これは、エックス線波長において、血液の吸収係数が血液酸化に対して大きく変化しないためである。しかしながら、赤外線エネルギは、血液体積における空間変動と血液酸化レベルとを識別できる。これは、これら波長における吸収係数がヘモグロビン状態の関数だからである。他の構造および機能も、赤外線エネルギを照射した組織の散乱係数における変動または変化によって識別できる。例えば、収縮中の筋肉組織、活動中の神経などである。これら構造は、投影技術では画像化できなかった。なぜなら投影技術は、分散係数における変動の測定に有効ではないからである。光学断層撮影法などの分散技術に基づく画像化で得られる測定値は、広範囲の疾病過程を診断する上で極めて大きな可能性を持っている。
【０００６】
散乱エネルギ測定に基づく画像化の代表的システムは、媒体に印加するための少なくとも１つのエネルギ源と、放出エネルギを検出するための少なくとも１つの検出器とを含む。エネルギ源は、画像化する媒体内で高度に散乱するように選択する。エネルギ源は、対象散乱媒体にエネルギを印加し、検出器は媒体表面が放出する散乱エネルギを測定する。この測定に基づき、媒体の内部特性の画像を再構築する。
【０００７】
この再構築には、代表的に「摂動法」を使用する。この方法は、基本的に対象散乱媒体から得た測定と、既知の参照散乱媒体とを比較する。この参照媒体は、物理媒体または仮想媒体であり、画像化する媒体の特性とできるだけ近い特性を有するように選択する。参照媒体の選択は、基本的に、対象の特性を初期想定することである。再構築の第１段階において、「順モデル」を使用し、参照媒体の既知内部特性に基づき、特定のエネルギ源位置に関する検出器の読みがどのようであるかを予測する。順モデルは、輸送方程式、またはその近似式である拡散方程式に基づく。この拡散方程式は、散乱媒体におけるフォトンの伝播を説明している。次に、輸送方程式の摂動定式化を使用し、（１）対象および参照からの検出器による測定値と予測値との差と、（２）対象および参照の未知および既知の内部特性間の差とを関係させる。この関係を、対象の未知の散乱および吸収特性について解く。内部特性の最終分布は、画像として表示または印刷する。
【０００８】
光学断層撮影システム等、散乱技術に基づく画像化システムおよび方法は、吸収、拡散、散乱係数等、散乱媒体の内部特性を検査し画像化する手段を提供する。しかしながら、前記画像化システムおよび方法は、密度の高い散乱媒体のコントラスト特性を復元するものであり、せいぜい中程度の空間分解能を有する結果を生成するに過ぎない。画質を向上させる方策も知られているが（例えばニュートンタイプ）、これらはいずれも計算集約型であり、初期開始状態に極めて敏感である。
【０００９】
磁気共鳴映像法（ＭＲＩ）における画像形成方法の基本は、測定した誘導電流の周波数と、磁界傾斜の空間方位との間に、１対１の対応が存在することである。磁界傾斜の空間方位は既知であるため、この対応により、測定した応答を空間における信号源に直接割り当てることができる。実際において、磁気共鳴現象の物理は、周波数シグネチャを、対象媒体に対して既知の空間関係を有する測定データ中に符号化する。より一般的に言えば、この種の方法は、「周波数符号化空間フィルタリング」として知られる。
【００１０】
前記理由により、散乱媒体の再構築画像など線形方程式系の解の品質を大きく向上させることができ計算効率が良い非線形補正方法が要求されている。
【００１１】
発明の概要
本発明は、この要求を満足させるため、計算効率が良く、散乱媒体の再構築画像の品質を向上させる画像再構築および画像補整用方法およびシステムを提供する。
【００１２】
本発明のシステムおよび方法の一実施例は、散乱媒体の再構築画像を向上させるため、（１）対象媒体を複数の体積要素に分割し、少なくとも１つの体積要素の光学係数に変調周波数を割り当て、（２）少なくとも１つのエネルギ源から前記対象媒体に１期間にわたりエネルギを印加し、前記対象媒体が放出するエネルギを少なくとも１つの検出器を介して測定し、（３）前記測定した対象媒体からの放出エネルギを処理して少なくとも１つの画像を再構築し、周波数符号化空間フィルタ（ＦＥＳＦ）を決定し、当該ＦＥＳＦを少なくとも１つの再構築画像に適用する。
【００１３】
本発明のシステムおよび方法の他の実施例は、線形方程式系の解を向上させるため、（１）対象媒体を複数の要素にモデル化し、前記対象媒体に少なくとも１つの局部的変動を加え、（２）少なくとも１つの局部的変動に起因する出力を測定し、（３）測定した出力を処理して結果を再構築し、補正フィルタを決定し、当該補正フィルタを前記結果に適用する。
【００１４】
前記利点および特徴は、代表的実施例のものに過ぎず、本発明の理解を助けるためだけを目的として提示した。これらは、請求項記載の本発明への制限、あるいは請求項の等価物への制限と見なすべきではない。例えば、これら利点のいくつかは、単一の実施例において同時実施が不可能であるがために、互いに相反するように見えるかもしれない。また、いくつかの利点は、基本的に本発明の１つの態様に適用可能である。従って、ここに要約した特徴および利点は、等価物を決定する上での方向を指定するものと見なすべきではない。本発明の他の特徴および利点は、以下の説明、図面、および請求項から明らかとなろう。
【００１５】
詳細な説明
１．はじめに
本発明のシステムおよび方法を、光学断層撮影の分野への適用に関して説明する。しかしながら、本原理体系は、線形応用を扱う広範囲の問題に適用できる。すなわち、線形摂動理論を応用して解を導く、経済、品質管理、疫学、気象学等の問題に適用できる。
【００１６】
画像再構築方法は、計算中心のアルゴリズムを使用する。このアルゴリズムは、標準線形摂動法を変更して画像再構築に適用する。このアルゴリズムの発展を今まで困難にしていた要因の１つは、与えられた画像再構築方法に関する情報広がり関数（ＩＳＦ）を定量的に特徴付ける方法が無かったことである。ここで使用する用語ＩＳＦは、対象媒体の与えられた位置に実際に存在する光学係数を画像の空間領域にマップするための厳密な方法を意味する。
【００１７】
ＩＳＦに関する情報が無い場合、観察した動作品質に応じて再構築アルゴリズムを変更する処理を体系化する明確な方法はない。再構築アルゴリズムと参照媒体との所定組み合わせに対してＩＳＦを特徴付けるため、本発明は、磁気共鳴映像法（ＭＲＩ）において発見された技術を利用する。この技術は、対象媒体と既知の空間関係を有する測定データ中に周波数応答を符号化する。本発明がＭＲＩと異なる点は、前記方法を画像形成に直接適用するのではなく、その概念を適用して周波数符号化空間フィルタ（ＦＥＳＦ）を得、そのフィルタを適用し、他の方法であらかじめ記録した画像の空間畳み込みを向上させることである。
【００１８】
これを実現するため、時系列画像から得られる位置依存時間周波数スペクトルを検査すればＦＥＳＦを得られるという認識をする。前記時系列画像の各要素における光学特性は、異なる時間依存特性が割り当てられている。完全画像化方法の場合、時系列の分析は、各ピクセルにおける時間的挙動を正確に復元する。実際には、空間畳み込みが存在する。この場合、畳み込みコントラスト特性の位置および大きさは、ピクセルデータの周波数スペクトルを検査することによって決定できる。しかしながら、全ピクセル中において一意に識別可能な時間特性を割り当てることにより、ピクセル毎に画像コントラストを正確に割り当てることが可能である。その結果の情報を線形演算子として使用し、テスト媒体から復元した画像のコントラスト特性を再配置（すなわち逆重畳）する。この方法において暗黙に仮定していることは、ＦＥＳＦが定義する空間畳み込みは、テスト媒体の画像に存在する畳み込みと同様であるということである。基本的に、どのような数のＦＥＳＦでも必要に応じて得ることができ適用することができる。
【００１９】
本システムおよび方法は、対象散乱媒体の画像を発生するための光学断層映像システムに関連して下記に詳細を説明する。しかしながら、当業者には明らかな通り、本発明の原理体系は、いかなるエネルギ源（例えば電磁気、音響等）、いかなる散乱媒体（例えば体組織、海洋、霧のかかった大気、地質学的地層、各種物質等）、いかなるエネルギ源状態（例えば時間非依存、時間調和、時間分解）、いかなる物理的画像化領域（例えば断面、体積測定）に基づく測定データからの画像再構築にも適用可能である。従って、本原理体系を拡張することにより、非破壊検査、地球物理学画像化、医療画像化、および調査技術を含む広範囲の応用分野において、新しい画像化手法が可能となる。
【００２０】
２．光学断層撮影システム
散乱媒体において画像再構築に使用する測定データを集めるための画像化システムは、数多く知られている。光学断層撮影システムの構成を図１に示す。このシステムは、コンピュータ１０２と、エネルギ源１０４および１０６と、ファイバスイッチャ１０８と、画像化ヘッド１１０と、検出器１１２と、データ取得ボード１１４と、エネルギ源ファイバ１２０と、検出器ファイバ１２２とを含む。
【００２１】
エネルギ源１０４および１０６は、光エネルギを提供する。この光エネルギは、ビームスプリッタ１１８を経由し、ファイバスイッチ１０８へ向かい、そこから１度に１本ずつ順次、複数のエネルギ源ファイバ１２０へ向かう。エネルギ源ファイバ１２０は、画像化ヘッド１１０の周囲に配置し、前記エネルギは対象媒体１１６の周囲の複数のエネルギ源位置に向かう。
【００２２】
前記エネルギは、画像化ヘッド１１０においてエネルギ源ファイバ１２０から離れ、画像化ヘッド１１０の中心にある対象媒体１１６へ入る。前記エネルギは、対象媒体を伝播するに従って散乱し、複数の箇所において対象媒体から出射する。この出射エネルギは、画像化ヘッド１１０の周囲に配置した検出器ファイバ１２２が収集する。検出されたエネルギは、検出器ファイバ１２２を検出器１１２まで伝搬する。検出器１１２はエネルギ測定装置を有し、測定に対応する信号を発生する。データ取得ボード１１４は、測定信号を受信し、それをコンピュータ１０２へ送る。
【００２３】
この処理を各エネルギ源位置について繰り返し、全検出器および全エネルギ源箇所について測定ベクトルを得る。コンピュータ１０２または適切な処理装置またはハードウエアを使用し、収集したデータを処理し、以下に詳細を説明する方法によって画像を再構築する。
【００２４】
３．方法
図２は、本発明の方法を説明するフローチャートである。本発明に基づく第１ステップは、フィルタ機能を計算する散乱媒体をＮ個の小領域または体積要素に分割する（ステップ２００）。次に、各領域／体積要素における光学パラメータ（例えば吸収係数および／または散乱係数）に正弦波時間変動を割り当てる。この時、各箇所は異なる周波数とする（ｆ１，ｆ２，．．．，ｆＮ）（ステップ２１０）。ステップ２１０における振動周波数（すなわち変調周波数）は、各周波数率ｆｍ／ｆｎ、ｎ  ｍが無理数であるように選択する。
【００２５】
次のステップは、順問題の解の時系列を計算することを含む。すなわちＩ（ｊ，ｔｋ）。ここでｊ＝１，２，．．．，Ｊであり、これは検出器インデックスである。ｋ＝１，２，．．．，Ｋであり、これは時間ステップインデックスである。この時間ステップインデックスは、結果としての動的媒体に対するものであり、Ｎ＞２ｍａｘ（ｆｍ）／ｍｉｎ（Δｆｍ）である（ステップ２２０）。ステップ２２０における周波数エイリアシングを避けるため、媒体の連続する状態間の時間間隔は、媒体における最高周波数の逆数に対して小さくなければならない。さらに、計算した時系列において十分高い周波数分解能を確保するため、測定の合計継続時間は、いずれの２つの割り当て周波数間の最小差の逆数に対しても、長くなければならない。ステップ２２０で得るデータは、検出器読取り値のＪ×Ｋ行列を形成する。
【００２６】
次のステップは、Ｋ逆問題を解くことを含む（すなわち断層撮影画像の時系列を再構築する）。この問題は、ステップ２２０で計算した検出器読取り値の各セットに１つある（ステップ２３０）。ステップ２３０において、検出器読取り値のＪ×Ｋ行列の各列は、方程式
【数１】
2004537731.tif 000003
【００２７】
の左辺を生成するために使用し、対応する
【数２】
2004537731.tif 000004
【００２８】
を計算する。このステップで得られるデータは、再構築した光学パラメータのＮ×Ｋ行列を形成する。第ｎｎ行は、第ｎピクセルまたはボクセルにおける光学パラメータ用時系列である。
【００２９】
画像時系列が完了したら、断層撮影画像の各ピクセルについて時間離散フーリエ変換（ＤＦＴ）を計算する（ステップ２４０）。その結果、各変調周波数におけるＤＦＴ振幅の空間マップが作成される（ステップ２５０）。ステップ２４０で計算したＤＦＴを連結して配列または行列にすることによりＦＥＳＦを決定する。ここで各行は１つの画像ピクセルにおけるＤＦＴ振幅に対応し、各列は１つの画像ピクセルにおける振幅に対応し、各列は特定周波数における振幅の空間マップに対応する（ステップ２６０）。ステップ２６０の結果は、単一の線形系、例えば
【数３】
2004537731.tif 000005
【００３０】
である。ここで
【数４】
2004537731.tif 000006
【００３１】
およびμａは、再構築のＮ×１ベクトルおよび真吸収係数であり、
【数５】
2004537731.tif 000007
【００３２】
はＮ×Ｎ行列である。この行列は、後述する通り、画像時系列から計算したＤＦＴ振幅の行列と、既知の理想ＤＦＴ振幅行列とを比較することによって決定する。
【数６】
2004537731.tif 000008
【００３３】
の決定は、すなわちＦＥＳＦは、直接ＬＵ分解（すなわちガウス消去法）を介して行う。特異値分解（ＳＶＤ）を使用することも可能である。時系列の各再構築画像にＦＥＳＦを適用することは、単純後退代入（すなわち再構築画像に対する空間逆重畳補正）を行うことである（ステップ２７０）。
【００３４】
３．１  順モデル
順モデル（すなわち順問題）に関する以下の説明は、再構築の第１ステップを明らかにする。この第１ステップは、参照媒体の既知の内部特性に基づき、特定エネルギ源箇所の検出器読取り値がどのようであるかを予測するために使用する。
【００３５】
前記した通り、代表的な再構築技術は、摂動法に基づく。基本的にこの方法は、参照媒体から予測した検出器測定と対象からの検出器測定との間の差を、対象の未知特性と参照の既知特性との差に関係づけ、それを解く。従って、再構築における最初のステップの１つは、参照媒体を選択し、モデル化によってあるいは物理測定によって検出器読取り値を予測することである。散乱媒体におけるエネルギ伝播をモデル化することは、輸送方程式またはその近似式である拡散方程式を用いて行う。これら方程式は、散乱媒体におけるフォトンの伝播を説明する。境界∂Λを有する領域について、これは次の式で表せる。
【数７】
2004537731.tif 000009
【００３６】
ここで
【数８】
2004537731.tif 000010
【００３７】
は位置
【数９】
2004537731.tif 000011
【００３８】
におけるフォトン強度であり、
【数１０】
2004537731.tif 000012
【００３９】
はＤＣ点エネルギ源の位置であり、
【数１１】
2004537731.tif 000013
【００４０】
および
【数１２】
2004537731.tif 000014
【００４１】
は位置依存拡散および吸収係数である。ここで拡散係数は、
【数１３】
2004537731.tif 000015
【００４２】
と定義し、
【数１４】
2004537731.tif 000016
【００４３】
は減少散乱係数である。この式を使用し、各エネルギ源位置に対する各検出器位置における参照媒体からのエネルギ放出を予測する。輸送または拡散方程式も、再構築で使用する摂動定式化または逆定式化のための基礎である。
【００４４】
３．２  逆定式化
逆定式化（すなわち逆問題）に関する以下の説明は、時系列の断層撮影画像を再構築する第２ステップを明らかにする。
【００４５】
前記した通り、対象媒体の吸収および／または散乱特性の断面画像を再構築することは、放射輸送または拡散方程式の摂動または逆定式化を解くことに基づく。摂動法は、未知の対象媒体の組成は既知の参照媒体から少しだけずれていると仮定する。これは、高度非線形問題を、検査対象媒体と参照媒体間の吸収および散乱特性における差に関して線形な問題に縮小する。結果としての光学逆または摂動定式化は、正規化差分法に基づき、次の形式を有する。
【数１５】
2004537731.tif 000017
【００４６】
ここで
【数１６】
2004537731.tif 000018
【００４７】
および
【数１７】
2004537731.tif 000019
【００４８】
は、対象（測定）媒体と、初期予測を発生するために使用した参照（計算または測定）媒体との、光学特性（吸収および拡散係数）間の断面差のベクトルである。
【数１８】
2004537731.tif 000020
【００４９】
および
【数１９】
2004537731.tif 000021
【００５０】
は、選択した参照媒体の吸収および拡散係数における局部接合が表面検出器に及ぼす影響を説明する重み行列である。
【数２０】
2004537731.tif 000022
【００５１】
は、２セットの検出器読取り値間の正規化差分を表し、次のように定義する。
【数２１】
2004537731.tif 000023
【００５２】
ここで
【数２２】
2004537731.tif 000024
【００５３】
は、選択した参照媒体に対応する計算した検出器読取り値である。
【数２３】
2004537731.tif 000025
【００５４】
は、２セットの測定データを表す（例えば背景対対象、時間平均対特定時間点など）。Ｍは、測定セットにおけるエネルギ源と検出器とのペア数である。
【００５５】
３．３  重み行列スケーリング
重み行列スケーリングに関する以下の説明は、逆問題において得られる重み行列のスケーリングを明らかにする。
【００５６】
重み行列をスケーリングすることの効果は、重み行列をより均一にし、良い状態にできることが多いことである。
【数２４】
2004537731.tif 000026
【００５７】
および
【数２５】
2004537731.tif 000027
【００５８】
の各列をその列ベクトルの平均値にスケーリングするスケーリング方法を使用する。しかしながら、どのような既知スケーリング方法を採用しても構わない。結果としての新しい重み行列の形式は次の通りである。
【数２６】
2004537731.tif 000028
【００５９】
ここでｋはμａまたはＤで良く、
【数２７】
2004537731.tif 000029
【００６０】
は正規化行列でありその要素は次の通りである。
【数２８】
2004537731.tif 000030
【００６１】
ここでＮは、領域Λを離散化する時に使用する要素の数である。結果としての系の方程式は次の通りである。
【数２９】
2004537731.tif 000031
【００６２】
ここで
【数３０】
2004537731.tif 000032
【００６３】
および
【数３１】
2004537731.tif 000033
【００６４】
である。なお
【数３２】
2004537731.tif 000034
【００６５】
は対角行列（式５）であり、その主対角要素は非ゼロである（式５）。従って、良く定義された逆を有し、その計算は自明のことである。
【００６６】
４．周波数符号化空間フィルタ（ＦＥＳＦ）
ＦＥＳＦに関する以下の説明は、補正画像を再構築するために使用するＦＥＳＦの決定を明らかにする。
【００６７】
画像再構築により生成する振幅の「空間マップ」および現実におけるＤＦＴの計算は、数字列である。各数字は、再構築アルゴリズムがＦＥＭノードの１つに割り当てた、ある特定周波数における振幅（すなわち３つまたはそれ以上の要素が一緒になる頂点または点の数）である。振幅マップの全セットは、行列に連結できる。
【数３３】
2004537731.tif 000035
【００６８】
ここでＮｆは周波数の数（＝有限要素の数）であり、ＮｎはＦＥＭノードの数である。これはＦＥＳＦの決定を理解する上で重要である。なぜなら、実際において、ノードの数は常に要素の数より小さいからである。すると
【数３４】
2004537731.tif 000036
【００６９】
（ｉは画像を意味する）は正方行列ではなく、その行数は列数の約半分である。
【００７０】
このようにして第２の行列が書ける。これは、各変調周波数が実際に媒体のどこに存在するかを正確に知らせる。
【数３５】
2004537731.tif 000037
【００７１】
行列
【数３６】
2004537731.tif 000038
【００７２】
は、まばらであることは間違いない（すなわちその要素のほとんどはゼロである）。なぜなら各ｆｍは、媒体の有限要素のただ１つのみに割り当てられるからである。実際、
【数３７】
2004537731.tif 000039
【００７３】
の各行は、合計で数百から数千（すなわちＮｎ＝１０３〜１０４）の要素を含んでおり、そのうちゼロでないのは正確に３要素（２次元媒体の場合）または４要素（３次元媒体）である。この数字が３または４であるのは、三角または四面体要素を使用するからであり、各要素が３または４ノードで拘束されるからである。
【００７４】
前記した再構築処理が完璧であれば、
【数３８】
2004537731.tif 000040
【００７５】
である。しかしながら、この理想結果は実際には起こらない。従って、
【数３９】
2004537731.tif 000041
【００７６】
を
【数４０】
2004537731.tif 000042
【００７７】
に変換する関数の性質に関するある種の仮定をしなければならない。本発明において行う仮定は、画像中のいずれか１つのノード（すなわちピクセル）において存在する周波数スペクトルは、媒体中の全ノードにおいて存在する周波数の線形関数であるということである。これは数学的に
【数４１】
2004537731.tif 000043
【００７８】
と書ける。ここで
【数４２】
2004537731.tif 000044
【００７９】
はＮｎ×Ｎｎ（すなわち正方）行列である。実際には、他の方向へ進む変換が望まれる。すなわち
【数４３】
2004537731.tif 000045
【００８０】
から開始し、真の
【数４４】
2004537731.tif 000046
【００８１】
にできるだけ近いものを生成する変換である。このため、実際に使用するフィルタの計算は、行列式
【数４５】
2004537731.tif 000047
【００８２】
を解くことによって完成する。ここで
【数４６】
2004537731.tif 000048
【００８３】
もＮｎ×Ｎｎ行列である。
【数４７】
2004537731.tif 000049
【００８４】
および
【数４８】
2004537731.tif 000050
【００８５】
は、互いに逆である。
【数４９】
2004537731.tif 000051
【００８６】
における要素の総数は、
【数５０】
2004537731.tif 000052
【００８７】
における数より小さい。なぜならＮｎ＜Ｎｆだからである。これは、この方法を適用する場合、完全な補正は起きないであろうことを意味する。なぜなら利用できる十分な補正項が無いからである。これは、
【数５１】
2004537731.tif 000053
【００８８】
と
【数５２】
2004537731.tif 000054
【００８９】
間に仮定した線形関係のためではなく、Ｎｆ要素における周波数をより少ない数のノードＮｎにマップすることに関して避け難い情報の損失があるために発生する。
【００９０】
このように計算されるＦＥＳＦは、品質制御において、再構築アルゴリズムの精度を定量化する方法に利用できる。さらにＦＥＳＦは、フィルタを作成するために使用した媒体とは異なる実験媒体のデータに使用する場合、画像向上ツールとして使用できる。例えば画像
【数５３】
2004537731.tif 000055
【００９１】
等のセットを再構築する場合、フィルタ機能が媒体特性に強く依存しない限り、
【数５４】
2004537731.tif 000056
【００９２】
等を計算することにより再構築の空間精度を向上できる。
【００９３】
５．実証結果
以下の例は、本発明の特徴および特性を示すため、本発明の実例を説明することだけを目的としており、本発明を限定するものではない。
【００９４】
ＦＥＳＦの効用は前記に例示した通りである。以下の説明において、ＦＥＳＦを２つの異なる時系列画像に適用した。いずれも同一セットの検出器読取り値から得たものであるが、異なる再構築アルゴリズムを使用した。使用した再構築方法は、原則として同一結果を生成すべきものであった。なぜなら、いずれか一方の使用を選択すべき自明の妥当な理由が無いからである。しかしながら、ＦＥＳＦ法の適用は、一方の再構築法はモデル化した媒体のいずれの箇所においても空間的に正確な摂動の画像を生成できるが、他方はできないことを示している。従って、両アルゴリズム用に計算したＩＳＦは、いずれも再構築画像に対する空間逆重畳補正の適用を可能にする。
【００９５】
図３は、規則的形状の２次元媒体（すなわち対象媒体）を示す。図３に示す通り、この媒体は直径８ｃｍの均一な円板であり、光学係数の値は  ａ＝０．０６ｃｍ−１、  ｓ＝１０ｃｍ−１である。順および逆問題のより好都合な解のため、図３に示す通り、円板の数学的境界を「物理」媒体の境界の外側へ０．５ｃｍ延ばした。延長領域内の係数の値は、「物理」媒体のものと同じとした。等間隔の単位強さを有する均一な１６の点エネルギ源を、物理境界上に示した位置に置いた。
【００９６】
図示のメッシュにおける有限要素およびノードの数は、それぞれ１６０４と８５０であり、最小要素面積と最大要素面積は、それぞれ０．０２５ｃｍ２と０．０７３ｃｍ２（平均±標準偏差＝０．０４０±０．００６ｃｍ２）である。正弦波変調を各要素における吸収係数に加えた。一意のｆｍを各々に割り当てた。一方で、振幅は各箇所において０．００６ｃｍ−１（すなわち平均値の１０％）であった。この予備調査に関し、要素の散乱係数は時間で変調しなかった。
【００９７】
分解能帯域幅がｆｍ間の最小差より小さく、かつナイキスト周波数が最大ｆｍより大きくなることを確実にするため、断層撮影検出器読取り値の１万セットの時系列をｔ−０．０１ｓにおいて計算した。画像再構築は、２つのアルゴリズムによって実行した。いずれのアルゴリズムも、画像演算子行列のＳＶＤに基づく。使用した第１のアルゴリズムは、先に説明した重み変換ＳＶＤＷＴ法である。第２の再構築方法は、ＳＶＤＷＴＷＲＳ、すなわちＳＶＤＷＴと追加の行列予備演算との結合である。後者において、各方程式をスケーリングし、重み行列の全行が同一合計を持つようにした。
【００９８】
画像再構築中に生成した１６０４のＤＦＴ振幅マップに２種類の分析を行った。第１に、各振幅マップと対象媒体における対応周波数の既知空間分布との間の全体空間相関を計算した（理想結果は、全ての周波数において正確に１．０に等しい相関である）。第２に、各マップの質量中心の座標を計算し、それを元に有限要素の幾何学的質量中心からの変位（理想結果は、全ての周波数において変位は正確に０．０に等しい）を容易に決定した。この有限要素の  ａは対応周波数において変調されている。これら２つの量をｆｍ（または、その等価として、対象媒体中の位置）の関数としてグラフにすると、ＳＶＤＷＴアルゴリズムについては図４Ａおよび４Ｂとなり、ＳＶＤＷＴＷＲＳアルゴリズムについては図５Ａおよび５Ｂとなる。図４Ａ、４Ｂおよび図５Ａ、５Ｂにおける薄い色の曲線は、フィルタ無しのＦＴ振幅空間分布から得た。濃い曲線は、前記した理論モデルに沿った可能な最良の補正を行った後、計算を繰り返して得た結果である。前記理論モデルによれば、再構築画像から得た振幅マップは、対象媒体に存在する真の空間分布の単純な線形変換である。
【００９９】
図４Ａ、４Ｂ、５Ａ、５Ｂを調べると、描かれた各関数は、第４００番ｆｍの後、挙動に質的変化を示している。この変化は、最初の４００の有限要素の全ては、物理境界と拡張境界との間、すなわちエネルギ源と検出器との輪の外側に位置する領域（図３参照）にあるという事実の単なる結果である。図４Ａおよび４Ｂをさらに調べて分かることは、いずれの空間精度測定値も、ほぼ第８００番目から第１１００番目のｆｍに対応する有限要素の理想値から特に離れていることである。これら要素は、対象媒体の中央領域に位置している。すなわち、ＳＶＤＷＴアルゴリズムで再構築した画像は、空間的に強くひずんでおり、中央領域の吸収計数値は、表面方向に著しく変位し、より周辺領域のものは、より高い精度で再構築されている。これに対し、（予備調整した）ＳＶＤＷＴＷＲＳアルゴリズムによって再構築した画像から得た振幅マップの場合、空間相関および質量中心変位は、空間的により均一である。これは重要な観察である。なぜなら、２つの再構築タイプは、所定セットの検出器データで動作する場合、理論的に同一解を生成すべきだからである。最後に、図４Ａ、４Ｂ、５Ａ、５Ｂの各々において、濃い（補正画像）曲線上のほとんどの点は、薄い（非補正画像）曲線上のものより、理想値に近いことが分かる。これは、再構築後に画像の空間精度を向上させるため、ＩＳＦ中の情報を使用できる可能性を示している。
【０１００】
以上の記述は、説明のための実施例を示したに過ぎない。以上の記述は、本発明の理解を容易にするため、可能な実施例の代表例、および本発明原理を明らかにするための例に焦点を当てた。本明細書は、可能な全ての変更形態を網羅するものではない。本発明の特定部分に関する変更実施例を提示していなくとも、あるいは本発明の一部について未説明の変更実施例が可能であっても、それら変更実施例を放棄するものではない。当業者であれば、本発明の要旨および範囲を逸脱することなく、他の応用例や実施例を推察できるであろう。従って本発明は、開示した実施例に限定されるものではなく、請求項に基づいて定義されるべきものである。未説明の実施例の多くは、請求項の範囲に入るものであり、その他も同等であると理解すべきである。
【図面の簡単な説明】
【０１０１】
本発明の様々な特徴および利点は、好適実施例の詳細説明および添付図面を参照することにより、よりよく理解できるであろう。
【図１】本発明に基づき使用される光学断層撮影システムの構成図である。
【図２】本発明の方法を説明するフローチャートである。
【図３】対象媒体を示す図である。
【図４】図４Ａは対象媒体における各振幅マップと対応周波数の既知空間分布との全体空間相関を示す図であり、重み変換特異値分解（ＳＶＤＷＴ）アルゴリズム用の変調周波数（ｆｍ）の関数として描いた図である。図４Ｂは各振幅マップの質量中心を示す図であり、ＳＶＤＷＴアルゴリズム用ｆｍの関数として描いた図である。
【図５】図５Ａは対象媒体における各振幅マップと対応周波数の既知空間分布との全体空間相関を示す図であり、追加の行列予備調整演算を伴う結合ＳＶＤＷＴ（ＳＶＤＷＴＷＲＳ）アルゴリズム用ｆｍの関数として描いた図である。図５Ｂは各振幅マップの質量中心を示す図であり、ＳＶＤＷＴＷＲＳアルゴリズム用ｆｍの関数として描いた図である。
【図１】
2004537731.tif 000057
【図２】
2004537731.tif 000058
【図４】
2004537731.tif 000059
【図５】
2004537731.tif 000060
【国際公開パンフレット】
2004537731.tif 000061
2004537731.tif 000062
2004537731.tif 000063
2004537731.tif 000064
2004537731.tif 000065
2004537731.tif 000066
2004537731.tif 000067
2004537731.tif 000068
2004537731.tif 000069
2004537731.tif 000070
2004537731.tif 000071
2004537731.tif 000072
2004537731.tif 000073
2004537731.tif 000074
2004537731.tif 000075
2004537731.tif 000076
2004537731.tif 000077
2004537731.tif 000078
2004537731.tif 000079
2004537731.tif 000080
2004537731.tif 000081
2004537731.tif 000082
2004537731.tif 000083
2004537731.tif 000084
2004537731.tif 000085
2004537731.tif 000086
2004537731.tif 000087
2004537731.tif 000088
2004537731.tif 000089
2004537731.tif 000090
【国際公開パンフレット（コレクトバージョン）】
2004537731.tif 000091
2004537731.tif 000092
【国際調査報告】
2004537731.tif 000093