JP2004537731A - 線形方程式系の解を向上させる方法およびシステム - Google Patents

線形方程式系の解を向上させる方法およびシステム Download PDF

Info

Publication number
JP2004537731A
JP2004537731A JP2003517836A JP2003517836A JP2004537731A JP 2004537731 A JP2004537731 A JP 2004537731A JP 2003517836 A JP2003517836 A JP 2003517836A JP 2003517836 A JP2003517836 A JP 2003517836A JP 2004537731 A JP2004537731 A JP 2004537731A
Authority
JP
Japan
Prior art keywords
target medium
energy
image
fesf
applying
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
JP2003517836A
Other languages
English (en)
Inventor
ランドール エル. バーボア、
ハリー エル. グレーバー、
ヤリング ペイ、
Original Assignee
ザ リサーチ ファウンデーション オブ ステイト ユニヴァーシティ オブ ニューヨーク
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by ザ リサーチ ファウンデーション オブ ステイト ユニヴァーシティ オブ ニューヨーク filed Critical ザ リサーチ ファウンデーション オブ ステイト ユニヴァーシティ オブ ニューヨーク
Publication of JP2004537731A publication Critical patent/JP2004537731A/ja
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/008Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/20Image enhancement or restoration using local operators
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/73Deblurring; Sharpening
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10101Optical tomography; Optical coherence tomography [OCT]

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)
  • Complex Calculations (AREA)
  • Image Processing (AREA)

Abstract

線形方程式系の解の品質を大きく向上させる方法である。線形方程式系の解を向上させるため、(1)対象媒体を複数の要素にモデル化し、前記対象媒体に少なくとも1つの局部的変動を加え、(2)少なくとも1つの局部的変動に起因する出力を測定し、(3)測定した出力を処理して結果を再構築し、補正フィルタを決定し、当該補正フィルタを前記結果に適用する。
【選択図】図1

Description

【技術分野】
【0001】
本発明は、NIH援助番号RO1−CA66184の元で米国政府の支援によりなされた。米国政府は本発明に幾分の権利を有する。
【0002】
この出願は、米国特許法120条に基づき、2001年8月2日に先行出願した米国予備特許出願第60/309,572号「散乱媒体の画像化品質を向上させるための周波数符号化空間フィルタリング方法」の利益を主張する。
【0003】
技術分野
本発明は、線形方程式の分野に関し、特に線形方程式系、例えば散乱媒体の画像化に使用する類の方程式について、その解を向上させるため補正フィルタを使用することに関する。
【0004】
背景技術
一般に散乱媒体の画像化は、散乱媒体内部の特性(吸収係数または散乱係数など)の空間分布画像を作成する様式に関する。これは、当該媒体にエネルギを導入し、その媒体が放出する散乱エネルギを検出することによって行う。この種のシステムおよび方法と対照的であるのは、エックス線などの投影画像化システムである。例えばエックス線システムは、エックス線エネルギ源と検出器との間の直線路を移動するエネルギの減衰または吸収を測定し画像化するものであり、散乱エネルギではない。エネルギが基本的に高分散であるか、あるいは基本的に直線路を移動するかは、エネルギの波長およびそれが移動する媒体の関数である。
【0005】
散乱技術に基づく画像化は、人体、地層、大気などの特性を画像化することに関し、新しいエネルギ波長の使用を可能にする。これら特性は、投影技法および投影波長によっては画像化できない。例えばエックス線投影技術は、骨構造等の密度の高い対象を画像化することに優れているが、血液酸化レベルの識別および画像化には比較的効果がない。これは、エックス線波長において、血液の吸収係数が血液酸化に対して大きく変化しないためである。しかしながら、赤外線エネルギは、血液体積における空間変動と血液酸化レベルとを識別できる。これは、これら波長における吸収係数がヘモグロビン状態の関数だからである。他の構造および機能も、赤外線エネルギを照射した組織の散乱係数における変動または変化によって識別できる。例えば、収縮中の筋肉組織、活動中の神経などである。これら構造は、投影技術では画像化できなかった。なぜなら投影技術は、分散係数における変動の測定に有効ではないからである。光学断層撮影法などの分散技術に基づく画像化で得られる測定値は、広範囲の疾病過程を診断する上で極めて大きな可能性を持っている。
【0006】
散乱エネルギ測定に基づく画像化の代表的システムは、媒体に印加するための少なくとも1つのエネルギ源と、放出エネルギを検出するための少なくとも1つの検出器とを含む。エネルギ源は、画像化する媒体内で高度に散乱するように選択する。エネルギ源は、対象散乱媒体にエネルギを印加し、検出器は媒体表面が放出する散乱エネルギを測定する。この測定に基づき、媒体の内部特性の画像を再構築する。
【0007】
この再構築には、代表的に「摂動法」を使用する。この方法は、基本的に対象散乱媒体から得た測定と、既知の参照散乱媒体とを比較する。この参照媒体は、物理媒体または仮想媒体であり、画像化する媒体の特性とできるだけ近い特性を有するように選択する。参照媒体の選択は、基本的に、対象の特性を初期想定することである。再構築の第1段階において、「順モデル」を使用し、参照媒体の既知内部特性に基づき、特定のエネルギ源位置に関する検出器の読みがどのようであるかを予測する。順モデルは、輸送方程式、またはその近似式である拡散方程式に基づく。この拡散方程式は、散乱媒体におけるフォトンの伝播を説明している。次に、輸送方程式の摂動定式化を使用し、(1)対象および参照からの検出器による測定値と予測値との差と、(2)対象および参照の未知および既知の内部特性間の差とを関係させる。この関係を、対象の未知の散乱および吸収特性について解く。内部特性の最終分布は、画像として表示または印刷する。
【0008】
光学断層撮影システム等、散乱技術に基づく画像化システムおよび方法は、吸収、拡散、散乱係数等、散乱媒体の内部特性を検査し画像化する手段を提供する。しかしながら、前記画像化システムおよび方法は、密度の高い散乱媒体のコントラスト特性を復元するものであり、せいぜい中程度の空間分解能を有する結果を生成するに過ぎない。画質を向上させる方策も知られているが(例えばニュートンタイプ)、これらはいずれも計算集約型であり、初期開始状態に極めて敏感である。
【0009】
磁気共鳴映像法(MRI)における画像形成方法の基本は、測定した誘導電流の周波数と、磁界傾斜の空間方位との間に、1対1の対応が存在することである。磁界傾斜の空間方位は既知であるため、この対応により、測定した応答を空間における信号源に直接割り当てることができる。実際において、磁気共鳴現象の物理は、周波数シグネチャを、対象媒体に対して既知の空間関係を有する測定データ中に符号化する。より一般的に言えば、この種の方法は、「周波数符号化空間フィルタリング」として知られる。
【0010】
前記理由により、散乱媒体の再構築画像など線形方程式系の解の品質を大きく向上させることができ計算効率が良い非線形補正方法が要求されている。
【0011】
発明の概要
本発明は、この要求を満足させるため、計算効率が良く、散乱媒体の再構築画像の品質を向上させる画像再構築および画像補整用方法およびシステムを提供する。
【0012】
本発明のシステムおよび方法の一実施例は、散乱媒体の再構築画像を向上させるため、(1)対象媒体を複数の体積要素に分割し、少なくとも1つの体積要素の光学係数に変調周波数を割り当て、(2)少なくとも1つのエネルギ源から前記対象媒体に1期間にわたりエネルギを印加し、前記対象媒体が放出するエネルギを少なくとも1つの検出器を介して測定し、(3)前記測定した対象媒体からの放出エネルギを処理して少なくとも1つの画像を再構築し、周波数符号化空間フィルタ(FESF)を決定し、当該FESFを少なくとも1つの再構築画像に適用する。
【0013】
本発明のシステムおよび方法の他の実施例は、線形方程式系の解を向上させるため、(1)対象媒体を複数の要素にモデル化し、前記対象媒体に少なくとも1つの局部的変動を加え、(2)少なくとも1つの局部的変動に起因する出力を測定し、(3)測定した出力を処理して結果を再構築し、補正フィルタを決定し、当該補正フィルタを前記結果に適用する。
【0014】
前記利点および特徴は、代表的実施例のものに過ぎず、本発明の理解を助けるためだけを目的として提示した。これらは、請求項記載の本発明への制限、あるいは請求項の等価物への制限と見なすべきではない。例えば、これら利点のいくつかは、単一の実施例において同時実施が不可能であるがために、互いに相反するように見えるかもしれない。また、いくつかの利点は、基本的に本発明の1つの態様に適用可能である。従って、ここに要約した特徴および利点は、等価物を決定する上での方向を指定するものと見なすべきではない。本発明の他の特徴および利点は、以下の説明、図面、および請求項から明らかとなろう。
【0015】
詳細な説明
1.はじめに
本発明のシステムおよび方法を、光学断層撮影の分野への適用に関して説明する。しかしながら、本原理体系は、線形応用を扱う広範囲の問題に適用できる。すなわち、線形摂動理論を応用して解を導く、経済、品質管理、疫学、気象学等の問題に適用できる。
【0016】
画像再構築方法は、計算中心のアルゴリズムを使用する。このアルゴリズムは、標準線形摂動法を変更して画像再構築に適用する。このアルゴリズムの発展を今まで困難にしていた要因の1つは、与えられた画像再構築方法に関する情報広がり関数(ISF)を定量的に特徴付ける方法が無かったことである。ここで使用する用語ISFは、対象媒体の与えられた位置に実際に存在する光学係数を画像の空間領域にマップするための厳密な方法を意味する。
【0017】
ISFに関する情報が無い場合、観察した動作品質に応じて再構築アルゴリズムを変更する処理を体系化する明確な方法はない。再構築アルゴリズムと参照媒体との所定組み合わせに対してISFを特徴付けるため、本発明は、磁気共鳴映像法(MRI)において発見された技術を利用する。この技術は、対象媒体と既知の空間関係を有する測定データ中に周波数応答を符号化する。本発明がMRIと異なる点は、前記方法を画像形成に直接適用するのではなく、その概念を適用して周波数符号化空間フィルタ(FESF)を得、そのフィルタを適用し、他の方法であらかじめ記録した画像の空間畳み込みを向上させることである。
【0018】
これを実現するため、時系列画像から得られる位置依存時間周波数スペクトルを検査すればFESFを得られるという認識をする。前記時系列画像の各要素における光学特性は、異なる時間依存特性が割り当てられている。完全画像化方法の場合、時系列の分析は、各ピクセルにおける時間的挙動を正確に復元する。実際には、空間畳み込みが存在する。この場合、畳み込みコントラスト特性の位置および大きさは、ピクセルデータの周波数スペクトルを検査することによって決定できる。しかしながら、全ピクセル中において一意に識別可能な時間特性を割り当てることにより、ピクセル毎に画像コントラストを正確に割り当てることが可能である。その結果の情報を線形演算子として使用し、テスト媒体から復元した画像のコントラスト特性を再配置(すなわち逆重畳)する。この方法において暗黙に仮定していることは、FESFが定義する空間畳み込みは、テスト媒体の画像に存在する畳み込みと同様であるということである。基本的に、どのような数のFESFでも必要に応じて得ることができ適用することができる。
【0019】
本システムおよび方法は、対象散乱媒体の画像を発生するための光学断層映像システムに関連して下記に詳細を説明する。しかしながら、当業者には明らかな通り、本発明の原理体系は、いかなるエネルギ源(例えば電磁気、音響等)、いかなる散乱媒体(例えば体組織、海洋、霧のかかった大気、地質学的地層、各種物質等)、いかなるエネルギ源状態(例えば時間非依存、時間調和、時間分解)、いかなる物理的画像化領域(例えば断面、体積測定)に基づく測定データからの画像再構築にも適用可能である。従って、本原理体系を拡張することにより、非破壊検査、地球物理学画像化、医療画像化、および調査技術を含む広範囲の応用分野において、新しい画像化手法が可能となる。
【0020】
2.光学断層撮影システム
散乱媒体において画像再構築に使用する測定データを集めるための画像化システムは、数多く知られている。光学断層撮影システムの構成を図1に示す。このシステムは、コンピュータ102と、エネルギ源104および106と、ファイバスイッチャ108と、画像化ヘッド110と、検出器112と、データ取得ボード114と、エネルギ源ファイバ120と、検出器ファイバ122とを含む。
【0021】
エネルギ源104および106は、光エネルギを提供する。この光エネルギは、ビームスプリッタ118を経由し、ファイバスイッチ108へ向かい、そこから1度に1本ずつ順次、複数のエネルギ源ファイバ120へ向かう。エネルギ源ファイバ120は、画像化ヘッド110の周囲に配置し、前記エネルギは対象媒体116の周囲の複数のエネルギ源位置に向かう。
【0022】
前記エネルギは、画像化ヘッド110においてエネルギ源ファイバ120から離れ、画像化ヘッド110の中心にある対象媒体116へ入る。前記エネルギは、対象媒体を伝播するに従って散乱し、複数の箇所において対象媒体から出射する。この出射エネルギは、画像化ヘッド110の周囲に配置した検出器ファイバ122が収集する。検出されたエネルギは、検出器ファイバ122を検出器112まで伝搬する。検出器112はエネルギ測定装置を有し、測定に対応する信号を発生する。データ取得ボード114は、測定信号を受信し、それをコンピュータ102へ送る。
【0023】
この処理を各エネルギ源位置について繰り返し、全検出器および全エネルギ源箇所について測定ベクトルを得る。コンピュータ102または適切な処理装置またはハードウエアを使用し、収集したデータを処理し、以下に詳細を説明する方法によって画像を再構築する。
【0024】
3.方法
図2は、本発明の方法を説明するフローチャートである。本発明に基づく第1ステップは、フィルタ機能を計算する散乱媒体をN個の小領域または体積要素に分割する(ステップ200)。次に、各領域/体積要素における光学パラメータ(例えば吸収係数および/または散乱係数)に正弦波時間変動を割り当てる。この時、各箇所は異なる周波数とする(f,f,...,f)(ステップ210)。ステップ210における振動周波数(すなわち変調周波数)は、各周波数率f/f、n mが無理数であるように選択する。
【0025】
次のステップは、順問題の解の時系列を計算することを含む。すなわちI(j,t)。ここでj=1,2,...,Jであり、これは検出器インデックスである。k=1,2,...,Kであり、これは時間ステップインデックスである。この時間ステップインデックスは、結果としての動的媒体に対するものであり、N>2max(fm)/min(Δfm)である(ステップ220)。ステップ220における周波数エイリアシングを避けるため、媒体の連続する状態間の時間間隔は、媒体における最高周波数の逆数に対して小さくなければならない。さらに、計算した時系列において十分高い周波数分解能を確保するため、測定の合計継続時間は、いずれの2つの割り当て周波数間の最小差の逆数に対しても、長くなければならない。ステップ220で得るデータは、検出器読取り値のJ×K行列を形成する。
【0026】
次のステップは、K逆問題を解くことを含む(すなわち断層撮影画像の時系列を再構築する)。この問題は、ステップ220で計算した検出器読取り値の各セットに1つある(ステップ230)。ステップ230において、検出器読取り値のJ×K行列の各列は、方程式
【数1】
Figure 2004537731
【0027】
の左辺を生成するために使用し、対応する
【数2】
Figure 2004537731
【0028】
を計算する。このステップで得られるデータは、再構築した光学パラメータのN×K行列を形成する。第nn行は、第nピクセルまたはボクセルにおける光学パラメータ用時系列である。
【0029】
画像時系列が完了したら、断層撮影画像の各ピクセルについて時間離散フーリエ変換(DFT)を計算する(ステップ240)。その結果、各変調周波数におけるDFT振幅の空間マップが作成される(ステップ250)。ステップ240で計算したDFTを連結して配列または行列にすることによりFESFを決定する。ここで各行は1つの画像ピクセルにおけるDFT振幅に対応し、各列は1つの画像ピクセルにおける振幅に対応し、各列は特定周波数における振幅の空間マップに対応する(ステップ260)。ステップ260の結果は、単一の線形系、例えば
【数3】
Figure 2004537731
【0030】
である。ここで
【数4】
Figure 2004537731
【0031】
およびμaは、再構築のN×1ベクトルおよび真吸収係数であり、
【数5】
Figure 2004537731
【0032】
はN×N行列である。この行列は、後述する通り、画像時系列から計算したDFT振幅の行列と、既知の理想DFT振幅行列とを比較することによって決定する。
【数6】
Figure 2004537731
【0033】
の決定は、すなわちFESFは、直接LU分解(すなわちガウス消去法)を介して行う。特異値分解(SVD)を使用することも可能である。時系列の各再構築画像にFESFを適用することは、単純後退代入(すなわち再構築画像に対する空間逆重畳補正)を行うことである(ステップ270)。
【0034】
3.1 順モデル
順モデル(すなわち順問題)に関する以下の説明は、再構築の第1ステップを明らかにする。この第1ステップは、参照媒体の既知の内部特性に基づき、特定エネルギ源箇所の検出器読取り値がどのようであるかを予測するために使用する。
【0035】
前記した通り、代表的な再構築技術は、摂動法に基づく。基本的にこの方法は、参照媒体から予測した検出器測定と対象からの検出器測定との間の差を、対象の未知特性と参照の既知特性との差に関係づけ、それを解く。従って、再構築における最初のステップの1つは、参照媒体を選択し、モデル化によってあるいは物理測定によって検出器読取り値を予測することである。散乱媒体におけるエネルギ伝播をモデル化することは、輸送方程式またはその近似式である拡散方程式を用いて行う。これら方程式は、散乱媒体におけるフォトンの伝播を説明する。境界∂Λを有する領域について、これは次の式で表せる。
【数7】
Figure 2004537731
【0036】
ここで
【数8】
Figure 2004537731
【0037】
は位置
【数9】
Figure 2004537731
【0038】
におけるフォトン強度であり、
【数10】
Figure 2004537731
【0039】
はDC点エネルギ源の位置であり、
【数11】
Figure 2004537731
【0040】
および
【数12】
Figure 2004537731
【0041】
は位置依存拡散および吸収係数である。ここで拡散係数は、
【数13】
Figure 2004537731
【0042】
と定義し、
【数14】
Figure 2004537731
【0043】
は減少散乱係数である。この式を使用し、各エネルギ源位置に対する各検出器位置における参照媒体からのエネルギ放出を予測する。輸送または拡散方程式も、再構築で使用する摂動定式化または逆定式化のための基礎である。
【0044】
3.2 逆定式化
逆定式化(すなわち逆問題)に関する以下の説明は、時系列の断層撮影画像を再構築する第2ステップを明らかにする。
【0045】
前記した通り、対象媒体の吸収および/または散乱特性の断面画像を再構築することは、放射輸送または拡散方程式の摂動または逆定式化を解くことに基づく。摂動法は、未知の対象媒体の組成は既知の参照媒体から少しだけずれていると仮定する。これは、高度非線形問題を、検査対象媒体と参照媒体間の吸収および散乱特性における差に関して線形な問題に縮小する。結果としての光学逆または摂動定式化は、正規化差分法に基づき、次の形式を有する。
【数15】
Figure 2004537731
【0046】
ここで
【数16】
Figure 2004537731
【0047】
および
【数17】
Figure 2004537731
【0048】
は、対象(測定)媒体と、初期予測を発生するために使用した参照(計算または測定)媒体との、光学特性(吸収および拡散係数)間の断面差のベクトルである。
【数18】
Figure 2004537731
【0049】
および
【数19】
Figure 2004537731
【0050】
は、選択した参照媒体の吸収および拡散係数における局部接合が表面検出器に及ぼす影響を説明する重み行列である。
【数20】
Figure 2004537731
【0051】
は、2セットの検出器読取り値間の正規化差分を表し、次のように定義する。
【数21】
Figure 2004537731
【0052】
ここで
【数22】
Figure 2004537731
【0053】
は、選択した参照媒体に対応する計算した検出器読取り値である。
【数23】
Figure 2004537731
【0054】
は、2セットの測定データを表す(例えば背景対対象、時間平均対特定時間点など)。Mは、測定セットにおけるエネルギ源と検出器とのペア数である。
【0055】
3.3 重み行列スケーリング
重み行列スケーリングに関する以下の説明は、逆問題において得られる重み行列のスケーリングを明らかにする。
【0056】
重み行列をスケーリングすることの効果は、重み行列をより均一にし、良い状態にできることが多いことである。
【数24】
Figure 2004537731
【0057】
および
【数25】
Figure 2004537731
【0058】
の各列をその列ベクトルの平均値にスケーリングするスケーリング方法を使用する。しかしながら、どのような既知スケーリング方法を採用しても構わない。結果としての新しい重み行列の形式は次の通りである。
【数26】
Figure 2004537731
【0059】
ここでkはμまたはDで良く、
【数27】
Figure 2004537731
【0060】
は正規化行列でありその要素は次の通りである。
【数28】
Figure 2004537731
【0061】
ここでNは、領域Λを離散化する時に使用する要素の数である。結果としての系の方程式は次の通りである。
【数29】
Figure 2004537731
【0062】
ここで
【数30】
Figure 2004537731
【0063】
および
【数31】
Figure 2004537731
【0064】
である。なお
【数32】
Figure 2004537731
【0065】
は対角行列(式5)であり、その主対角要素は非ゼロである(式5)。従って、良く定義された逆を有し、その計算は自明のことである。
【0066】
4.周波数符号化空間フィルタ(FESF)
FESFに関する以下の説明は、補正画像を再構築するために使用するFESFの決定を明らかにする。
【0067】
画像再構築により生成する振幅の「空間マップ」および現実におけるDFTの計算は、数字列である。各数字は、再構築アルゴリズムがFEMノードの1つに割り当てた、ある特定周波数における振幅(すなわち3つまたはそれ以上の要素が一緒になる頂点または点の数)である。振幅マップの全セットは、行列に連結できる。
【数33】
Figure 2004537731
【0068】
ここでNは周波数の数(=有限要素の数)であり、NはFEMノードの数である。これはFESFの決定を理解する上で重要である。なぜなら、実際において、ノードの数は常に要素の数より小さいからである。すると
【数34】
Figure 2004537731
【0069】
(iは画像を意味する)は正方行列ではなく、その行数は列数の約半分である。
【0070】
このようにして第2の行列が書ける。これは、各変調周波数が実際に媒体のどこに存在するかを正確に知らせる。
【数35】
Figure 2004537731
【0071】
行列
【数36】
Figure 2004537731
【0072】
は、まばらであることは間違いない(すなわちその要素のほとんどはゼロである)。なぜなら各fは、媒体の有限要素のただ1つのみに割り当てられるからである。実際、
【数37】
Figure 2004537731
【0073】
の各行は、合計で数百から数千(すなわちN=10〜10)の要素を含んでおり、そのうちゼロでないのは正確に3要素(2次元媒体の場合)または4要素(3次元媒体)である。この数字が3または4であるのは、三角または四面体要素を使用するからであり、各要素が3または4ノードで拘束されるからである。
【0074】
前記した再構築処理が完璧であれば、
【数38】
Figure 2004537731
【0075】
である。しかしながら、この理想結果は実際には起こらない。従って、
【数39】
Figure 2004537731
【0076】

【数40】
Figure 2004537731
【0077】
に変換する関数の性質に関するある種の仮定をしなければならない。本発明において行う仮定は、画像中のいずれか1つのノード(すなわちピクセル)において存在する周波数スペクトルは、媒体中の全ノードにおいて存在する周波数の線形関数であるということである。これは数学的に
【数41】
Figure 2004537731
【0078】
と書ける。ここで
【数42】
Figure 2004537731
【0079】
はN×N(すなわち正方)行列である。実際には、他の方向へ進む変換が望まれる。すなわち
【数43】
Figure 2004537731
【0080】
から開始し、真の
【数44】
Figure 2004537731
【0081】
にできるだけ近いものを生成する変換である。このため、実際に使用するフィルタの計算は、行列式
【数45】
Figure 2004537731
【0082】
を解くことによって完成する。ここで
【数46】
Figure 2004537731
【0083】
もN×N行列である。
【数47】
Figure 2004537731
【0084】
および
【数48】
Figure 2004537731
【0085】
は、互いに逆である。
【数49】
Figure 2004537731
【0086】
における要素の総数は、
【数50】
Figure 2004537731
【0087】
における数より小さい。なぜならN<Nだからである。これは、この方法を適用する場合、完全な補正は起きないであろうことを意味する。なぜなら利用できる十分な補正項が無いからである。これは、
【数51】
Figure 2004537731
【0088】

【数52】
Figure 2004537731
【0089】
間に仮定した線形関係のためではなく、N要素における周波数をより少ない数のノードNにマップすることに関して避け難い情報の損失があるために発生する。
【0090】
このように計算されるFESFは、品質制御において、再構築アルゴリズムの精度を定量化する方法に利用できる。さらにFESFは、フィルタを作成するために使用した媒体とは異なる実験媒体のデータに使用する場合、画像向上ツールとして使用できる。例えば画像
【数53】
Figure 2004537731
【0091】
等のセットを再構築する場合、フィルタ機能が媒体特性に強く依存しない限り、
【数54】
Figure 2004537731
【0092】
等を計算することにより再構築の空間精度を向上できる。
【0093】
5.実証結果
以下の例は、本発明の特徴および特性を示すため、本発明の実例を説明することだけを目的としており、本発明を限定するものではない。
【0094】
FESFの効用は前記に例示した通りである。以下の説明において、FESFを2つの異なる時系列画像に適用した。いずれも同一セットの検出器読取り値から得たものであるが、異なる再構築アルゴリズムを使用した。使用した再構築方法は、原則として同一結果を生成すべきものであった。なぜなら、いずれか一方の使用を選択すべき自明の妥当な理由が無いからである。しかしながら、FESF法の適用は、一方の再構築法はモデル化した媒体のいずれの箇所においても空間的に正確な摂動の画像を生成できるが、他方はできないことを示している。従って、両アルゴリズム用に計算したISFは、いずれも再構築画像に対する空間逆重畳補正の適用を可能にする。
【0095】
図3は、規則的形状の2次元媒体(すなわち対象媒体)を示す。図3に示す通り、この媒体は直径8cmの均一な円板であり、光学係数の値は =0.06cm−1=10cm−1である。順および逆問題のより好都合な解のため、図3に示す通り、円板の数学的境界を「物理」媒体の境界の外側へ0.5cm延ばした。延長領域内の係数の値は、「物理」媒体のものと同じとした。等間隔の単位強さを有する均一な16の点エネルギ源を、物理境界上に示した位置に置いた。
【0096】
図示のメッシュにおける有限要素およびノードの数は、それぞれ1604と850であり、最小要素面積と最大要素面積は、それぞれ0.025cmと0.073cm(平均±標準偏差=0.040±0.006cm)である。正弦波変調を各要素における吸収係数に加えた。一意のfを各々に割り当てた。一方で、振幅は各箇所において0.006cm−1(すなわち平均値の10%)であった。この予備調査に関し、要素の散乱係数は時間で変調しなかった。
【0097】
分解能帯域幅がf間の最小差より小さく、かつナイキスト周波数が最大fより大きくなることを確実にするため、断層撮影検出器読取り値の1万セットの時系列をt−0.01sにおいて計算した。画像再構築は、2つのアルゴリズムによって実行した。いずれのアルゴリズムも、画像演算子行列のSVDに基づく。使用した第1のアルゴリズムは、先に説明した重み変換SVDWT法である。第2の再構築方法は、SVDWTWRS、すなわちSVDWTと追加の行列予備演算との結合である。後者において、各方程式をスケーリングし、重み行列の全行が同一合計を持つようにした。
【0098】
画像再構築中に生成した1604のDFT振幅マップに2種類の分析を行った。第1に、各振幅マップと対象媒体における対応周波数の既知空間分布との間の全体空間相関を計算した(理想結果は、全ての周波数において正確に1.0に等しい相関である)。第2に、各マップの質量中心の座標を計算し、それを元に有限要素の幾何学的質量中心からの変位(理想結果は、全ての周波数において変位は正確に0.0に等しい)を容易に決定した。この有限要素の は対応周波数において変調されている。これら2つの量をf(または、その等価として、対象媒体中の位置)の関数としてグラフにすると、SVDWTアルゴリズムについては図4Aおよび4Bとなり、SVDWTWRSアルゴリズムについては図5Aおよび5Bとなる。図4A、4Bおよび図5A、5Bにおける薄い色の曲線は、フィルタ無しのFT振幅空間分布から得た。濃い曲線は、前記した理論モデルに沿った可能な最良の補正を行った後、計算を繰り返して得た結果である。前記理論モデルによれば、再構築画像から得た振幅マップは、対象媒体に存在する真の空間分布の単純な線形変換である。
【0099】
図4A、4B、5A、5Bを調べると、描かれた各関数は、第400番fの後、挙動に質的変化を示している。この変化は、最初の400の有限要素の全ては、物理境界と拡張境界との間、すなわちエネルギ源と検出器との輪の外側に位置する領域(図3参照)にあるという事実の単なる結果である。図4Aおよび4Bをさらに調べて分かることは、いずれの空間精度測定値も、ほぼ第800番目から第1100番目のfに対応する有限要素の理想値から特に離れていることである。これら要素は、対象媒体の中央領域に位置している。すなわち、SVDWTアルゴリズムで再構築した画像は、空間的に強くひずんでおり、中央領域の吸収計数値は、表面方向に著しく変位し、より周辺領域のものは、より高い精度で再構築されている。これに対し、(予備調整した)SVDWTWRSアルゴリズムによって再構築した画像から得た振幅マップの場合、空間相関および質量中心変位は、空間的により均一である。これは重要な観察である。なぜなら、2つの再構築タイプは、所定セットの検出器データで動作する場合、理論的に同一解を生成すべきだからである。最後に、図4A、4B、5A、5Bの各々において、濃い(補正画像)曲線上のほとんどの点は、薄い(非補正画像)曲線上のものより、理想値に近いことが分かる。これは、再構築後に画像の空間精度を向上させるため、ISF中の情報を使用できる可能性を示している。
【0100】
以上の記述は、説明のための実施例を示したに過ぎない。以上の記述は、本発明の理解を容易にするため、可能な実施例の代表例、および本発明原理を明らかにするための例に焦点を当てた。本明細書は、可能な全ての変更形態を網羅するものではない。本発明の特定部分に関する変更実施例を提示していなくとも、あるいは本発明の一部について未説明の変更実施例が可能であっても、それら変更実施例を放棄するものではない。当業者であれば、本発明の要旨および範囲を逸脱することなく、他の応用例や実施例を推察できるであろう。従って本発明は、開示した実施例に限定されるものではなく、請求項に基づいて定義されるべきものである。未説明の実施例の多くは、請求項の範囲に入るものであり、その他も同等であると理解すべきである。
【図面の簡単な説明】
【0101】
本発明の様々な特徴および利点は、好適実施例の詳細説明および添付図面を参照することにより、よりよく理解できるであろう。
【図1】本発明に基づき使用される光学断層撮影システムの構成図である。
【図2】本発明の方法を説明するフローチャートである。
【図3】対象媒体を示す図である。
【図4】図4Aは対象媒体における各振幅マップと対応周波数の既知空間分布との全体空間相関を示す図であり、重み変換特異値分解(SVDWT)アルゴリズム用の変調周波数(f)の関数として描いた図である。図4Bは各振幅マップの質量中心を示す図であり、SVDWTアルゴリズム用fの関数として描いた図である。
【図5】図5Aは対象媒体における各振幅マップと対応周波数の既知空間分布との全体空間相関を示す図であり、追加の行列予備調整演算を伴う結合SVDWT(SVDWTWRS)アルゴリズム用fの関数として描いた図である。図5Bは各振幅マップの質量中心を示す図であり、SVDWTWRSアルゴリズム用fの関数として描いた図である。

Claims (21)

  1. 第1対象媒体を複数の体積要素に分割するステップと、
    前記体積要素の光学係数の少なくとも1つに変調周波数を割り当てるステップと、
    少なくとも1つのエネルギ源から1期間にわたり前記第1対象媒体にエネルギを印加するステップと、
    前記第1対象媒体が放出するエネルギを少なくとも1つの検出器を介して測定するステップと、
    前記第1対象媒体の測定放出エネルギを処理して少なくとも1つの画像を再構築するステップと、
    周波数符号化空間フィルタ(FESF)を決定するステップと、
    前記第1対象媒体の少なくとも1つの再構築画像に前記FESFを適用するステップと、
    を備える散乱媒体の再構築画像を向上させる方法。
  2. 前記変調周波数は、体積要素の吸収係数に割り当てる請求項1記載の方法。
  3. 前記変調周波数は、体積要素の散乱係数に割り当てる請求項1記載の方法。
  4. 前記第1対象媒体の測定放出エネルギは、摂動法を使用して処理する請求項1記載の方法。
  5. 前記使用する摂動法は、順問題の解を用いて断層撮影画像を再構築する請求項4記載の方法。
  6. 前記使用する摂動法は、逆問題を用いて断層撮影画像を再構築する請求項4記載の方法。
  7. FESFを決定するステップは、
    前記再構築した断層撮影画像の時間離散フーリエ変換を計算するステップと、
    前記計算した時間離散フーリエ変換を処理し、対応する体積要素に関する変調周波数における振幅を決定するステップと、
    を備える請求項1記載の方法。
  8. 前記FESFは、単純行列乗算を実行することにより、少なくとも1つの再構築画像に適用する請求項1記載の方法。
  9. 少なくとも1つのエネルギ源から1期間にわたり第2対象媒体にエネルギを印加するステップと、
    前記第2対象媒体が放出するエネルギを少なくとも1つの検出器を介して測定するステップと、
    前記第2対象媒体の測定放出エネルギを処理して少なくとも1つの画像を再構築するステップと、
    前記第1対象媒体から決定した前記FESFを前記第2対象媒体の少なくとも1つの再構築画像に適用するステップと、
    をさらに備える請求項1記載の方法。
  10. 第1対象媒体を複数の体積要素に分割する手段と、
    前記体積要素の光学係数の少なくとも1つに変調周波数を割り当てる手段と、
    少なくとも1つのエネルギ源から1期間にわたり前記第1対象媒体にエネルギを印加する手段と、
    前記第1対象媒体が放出するエネルギを少なくとも1つの検出器を介して測定する手段と、
    前記第1対象媒体の測定放出エネルギを処理して少なくとも1つの画像を再構築する手段と、
    FESFを決定する手段と、
    前記第1対象媒体の少なくとも1つの再構築画像に前記FESFを適用する手段と、
    を備える散乱媒体の再構築画像を向上させるシステム。
  11. 前記変調周波数は、体積要素の吸収係数に割り当てる請求項10記載の方法。
  12. 前記変調周波数は、体積要素の散乱係数に割り当てる請求項10記載の方法。
  13. 前記第1対象媒体の測定放出エネルギは、摂動法を使用して処理する請求項10記載の方法。
  14. 前記使用する摂動法は、順問題の解を用いて断層撮影画像を再構築する請求項13記載の方法。
  15. 前記使用する摂動法は、逆問題を用いて断層撮影画像を再構築する請求項13記載の方法。
  16. FESFを決定するステップは、
    前記再構築した断層撮影画像の時間離散フーリエ変換を計算する手段と、
    前記計算した時間離散フーリエ変換を処理し、対応する体積要素に関する変調周波数における振幅を決定する手段と、
    を備える請求項10記載の方法。
  17. 前記FESFは、単純行列乗算を実行することにより、少なくとも1つの再構築画像に適用する請求項10記載の方法。
  18. 少なくとも1つのエネルギ源から1期間にわたり第2対象媒体にエネルギを印加する手段と、
    前記第2対象媒体が放出するエネルギを少なくとも1つの検出器を介して測定する手段と、
    前記第2対象媒体の測定放出エネルギを処理して少なくとも1つの画像を再構築する手段と、
    前記第1対象媒体から決定した前記FESFを前記第2対象媒体の少なくとも1つの再構築画像に適用する手段と、
    をさらに備える請求項10記載の方法。
  19. プロセッサにより実行されると、第1対象媒体を複数の体積要素に分割する命令コードと、
    プロセッサにより実行されると、前記体積要素の光学係数の少なくとも1つに変調周波数を割り当てる命令コードと、
    プロセッサにより実行されると、少なくとも1つのエネルギ源から1期間にわたり前記第1対象媒体にエネルギを印加する命令コードと、
    プロセッサにより実行されると、前記第1対象媒体が放出するエネルギを少なくとも1つの検出器を介して測定する命令コードと、
    プロセッサにより実行されると、前記第1対象媒体の測定放出エネルギを処理して少なくとも1つの画像を再構築する命令コードと、
    プロセッサにより実行されると、FESFを決定する命令コードと、
    プロセッサにより実行されると、前記第1対象媒体の少なくとも1つの再構築画像に前記FESFを適用する命令コードと、
    を備えるコンピュータ読み取り可能媒体に格納されプロセッサによって実行可能なプログラム。
  20. プロセッサにより実行されると、少なくとも1つのエネルギ源から1期間にわたり第2対象媒体にエネルギを印加する命令コードと、
    プロセッサにより実行されると、前記第2対象媒体が放出するエネルギを少なくとも1つの検出器を介して測定する命令コードと、
    プロセッサにより実行されると、前記第2対象媒体の測定放出エネルギを処理して少なくとも1つの画像を再構築する命令コードと、
    プロセッサにより実行されると、前記第1対象媒体から決定した前記FESFを前記第2対象媒体の少なくとも1つの再構築画像に適用する命令コードと、
    をさらに備える請求項19記載のプログラム。
  21. 第1対象媒体を複数の要素にモデル化するステップと、
    前記対象媒体に少なくとも1つの局部変動を加えるステップと、
    少なくとも1つの局部変動に起因する出力を測定するステップと、
    前記測定した出力を処理して結果を再構築するステップと、
    補正フィルタを決定するステップと、
    前記補正フィルタを前記結果に適用するステップと、
    を備える線形方程式系の解を向上させる方法。
JP2003517836A 2001-08-02 2002-08-02 線形方程式系の解を向上させる方法およびシステム Pending JP2004537731A (ja)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US30957201P 2001-08-02 2001-08-02
PCT/US2002/024520 WO2003012736A2 (en) 2001-08-02 2002-08-02 Method and system for enhancing solutions to a system of linear equations

Publications (1)

Publication Number Publication Date
JP2004537731A true JP2004537731A (ja) 2004-12-16

Family

ID=23198755

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2003517836A Pending JP2004537731A (ja) 2001-08-02 2002-08-02 線形方程式系の解を向上させる方法およびシステム

Country Status (6)

Country Link
US (1) US7099519B2 (ja)
EP (1) EP1421545A4 (ja)
JP (1) JP2004537731A (ja)
AU (1) AU2002330963A1 (ja)
CA (1) CA2456112A1 (ja)
WO (1) WO2003012736A2 (ja)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2003088133A1 (en) * 2002-04-06 2003-10-23 Barbour Randall L Modification of the normalized difference method for real-time optical tomography
US7376282B2 (en) * 2003-11-20 2008-05-20 Xerox Corporation Method for designing nearly circularly symmetric descreening filters that can be efficiently implemented in VLIW (very long instruction word) media processors
CA2843525A1 (en) * 2011-07-29 2013-02-07 Shell Internationale Research Maatschappij B.V. Method for increasing broadside sensitivity in seismic sensing system
CN111096761B (zh) * 2018-10-29 2024-03-08 上海西门子医疗器械有限公司 修正楔形滤波器散射的方法、装置和相关设备

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11337476A (ja) * 1998-05-26 1999-12-10 Hamamatsu Photonics Kk 散乱吸収体の内部特性分布の計測方法及び装置
WO2001020305A1 (en) * 1999-09-14 2001-03-22 The Research Foundation Of State University Of New York, Technology Transfer Office Method and system for imaging the dynamics of scattering medium
JP2001324444A (ja) * 2000-05-12 2001-11-22 Shimadzu Corp 光測定装置

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
NL8005311A (nl) 1980-09-24 1982-04-16 Philips Nv Inrichting voor het afvragen van een meetgegeven met voorbewerking van de gedetekteerde signaalwaarden.
US5384573A (en) * 1990-10-29 1995-01-24 Essex Corporation Image synthesis using time sequential holography
US5751243A (en) * 1990-10-29 1998-05-12 Essex Corporation Image synthesis using time sequential holography
IL100530A (en) * 1991-12-26 1996-05-14 Elscint Ltd Evolving image
US5414623A (en) * 1992-05-08 1995-05-09 Iowa State University Research Foundation Optoelectronic system for implementation of iterative computer tomography algorithms
US5546472A (en) * 1992-08-07 1996-08-13 Arch Development Corp. Feature guided method and apparatus for obtaining an image of an object
US5802218A (en) * 1994-11-04 1998-09-01 Motorola, Inc. Method, post-processing filter, and video compression system for suppressing mosquito and blocking atrifacts
US5764307A (en) * 1995-07-24 1998-06-09 Motorola, Inc. Method and apparatus for spatially adaptive filtering for video encoding
US6041248A (en) * 1997-05-07 2000-03-21 The Texas A&M University System Method and apparatus for frequency encoded ultrasound-modulated optical tomography of dense turbid media

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11337476A (ja) * 1998-05-26 1999-12-10 Hamamatsu Photonics Kk 散乱吸収体の内部特性分布の計測方法及び装置
WO2001020305A1 (en) * 1999-09-14 2001-03-22 The Research Foundation Of State University Of New York, Technology Transfer Office Method and system for imaging the dynamics of scattering medium
JP2003527883A (ja) * 1999-09-14 2003-09-24 ザ・リサーチ・ファンデーション・オブ・ステート・ユニバーシティ・オブ・ニューヨーク 散乱媒体の動力学を画像化する方法及びシステム
JP2001324444A (ja) * 2000-05-12 2001-11-22 Shimadzu Corp 光測定装置

Also Published As

Publication number Publication date
EP1421545A2 (en) 2004-05-26
EP1421545A4 (en) 2010-08-04
US7099519B2 (en) 2006-08-29
US20050041882A1 (en) 2005-02-24
WO2003012736A3 (en) 2004-03-04
CA2456112A1 (en) 2003-02-13
AU2002330963A1 (en) 2003-02-17
WO2003012736A2 (en) 2003-02-13

Similar Documents

Publication Publication Date Title
JP4271040B2 (ja) 実時間の光学的トモグラフィーの正規化された差方法の変更
Barbour et al. MRI-guided optical tomography: prospects and computation for a new imaging method
JP5869476B2 (ja) 断層撮影画像の取得および再構成を実行するシステムおよび方法
JP2018529409A (ja) 超高密度の電極に基づく脳撮像システム
JP6214819B1 (ja) 微分位相コントラストx線撮像における暗視野信号の最適なエネルギ加重
JP6884344B2 (ja) 脳内ネットワークの活動推定システム、脳内ネットワークの活動推定方法、脳内ネットワークの活動推定プログラム、および、学習済み脳活動推定モデル
Gross et al. Properties of MEG tomographic maps obtained with spatial filtering
JP6076981B2 (ja) 異なるモダリティのx線画像の周波数依存複合
Bertero et al. Inverse problems in biomedical imaging: modeling and methods of solution
WO2013159250A1 (zh) 一种动态荧光分子断层图像的重建方法
JP2017521167A (ja) 信号処理システム、信号処理方法、コンピュータプログラム及び記憶媒体
US7617080B2 (en) Image enhancement by spatial linear deconvolution
JP2004537731A (ja) 線形方程式系の解を向上させる方法およびシステム
KR20110121536A (ko) 양전자 단층 촬영 영상에서 워블 동작과 psf을 이용한 초해상도 촬영 장치 및 방법
Shifa et al. Improved image reconstruction using multi frequency data for diffuse optical tomography
JP4543774B2 (ja) 生体光計測装置
JP4567318B2 (ja) 生体組織部分内の領域の位置測定方法および装置
Goncharsky et al. Three-dimensional ultrasound tomography: mathematical methods and experimental results
Wong et al. Objective assessment and design improvement of a staring, sparse transducer array by the spatial crosstalk matrix for 3d photoacoustic tomography
Yao et al. Frequency domain optical tomography in human tissue
Lila et al. Representation and reconstruction of covariance operators in linear inverse problems
CN114159021B (zh) 基于双输入-单输出深度学习的切伦可夫激发的荧光扫描断层成像重建方法
Holt et al. Toward ideal imaging geometry for recovery independence in fluorescence molecular tomography
Abidi et al. Radial noise filtering in positron emission tomography
Ding et al. Optical topography guided semi-three-dimensional diffuse optical tomography for a multi-layer model of occipital cortex: a pilot methodological study

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20050720

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20071211

A601 Written request for extension of time

Free format text: JAPANESE INTERMEDIATE CODE: A601

Effective date: 20080310

A602 Written permission of extension of time

Free format text: JAPANESE INTERMEDIATE CODE: A602

Effective date: 20080317

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20080421

A02 Decision of refusal

Free format text: JAPANESE INTERMEDIATE CODE: A02

Effective date: 20081209

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20090309

A911 Transfer to examiner for re-examination before appeal (zenchi)

Free format text: JAPANESE INTERMEDIATE CODE: A911

Effective date: 20090417

A912 Re-examination (zenchi) completed and case transferred to appeal board

Free format text: JAPANESE INTERMEDIATE CODE: A912

Effective date: 20090619

A601 Written request for extension of time

Free format text: JAPANESE INTERMEDIATE CODE: A601

Effective date: 20100716

A602 Written permission of extension of time

Free format text: JAPANESE INTERMEDIATE CODE: A602

Effective date: 20100722

A601 Written request for extension of time

Free format text: JAPANESE INTERMEDIATE CODE: A601

Effective date: 20110331

A602 Written permission of extension of time

Free format text: JAPANESE INTERMEDIATE CODE: A602

Effective date: 20110405