JPH10512053A - 散乱放射の直接再生を利用した拡散画像化システムおよび方法 - Google Patents

Abstract

(57)【要約】 対象物７１０に時間領域ソース７２０を照射することによって行われる散乱放射の透過強度の測定から、拡散画像を直接に再生するためのシステムおよび方法が開示されている。検出器７３０によって検出された透過強度は積分演算子によって拡散係数と関係づけられている。画像は積分演算子に関連して決定された、所定の数学的アルゴリズムを透過強度について実行することによりプロセッサ７５０によって直接的に再生される。

Description

【発明の詳細な説明】 散乱放射の直接再生を利用した 拡散画像化システムおよび方法発明の分野 本発明は一般的には、対象物の拡散画像を生成するシステムおよびこれに付随 する方法論に関し、より具体的には、対象物に時間領域ソースを照射することに よって検出された散乱放射の測定から画像が直接的に再生されることを目的とし た、かかるシステムおよび方法論に関する。発明の背景 本発明は、物理的原理、および多重散乱状況(multiple scattering regime)に おける光学的画像化(optical imaging：光学的に画像を得ること）のための直接 再生方法の基礎となる関連の数学的公式とを取り扱うことを主題にしている。そ の結果は画像再生問題の解を直接的に求める方法論である。さらに、この方法は 、一般的には、拡散多重散乱状況において任意のスカラ波(scalar wave)を使用 して画像生成することに応用可能であり、光学的画像化に限定されるものではな い。しかし、本発明の重要な派生的問題を解明するために、この方法の応用分野 を１つ選択して、説明にある程度の明確性と具体性をもたせるようにすることが 非常に役立つ。従って、多くの生物系は本発明の原理、特に、光量子拡散画像化 (photon diffusion imaging)原理を応用するための物理的要件を満たしているの で、本発明の主題がもつ基本的側面は、医学的画像化を本発明の方法の応用例と して説明されている。 医学的画像化（medical imaging）は、この２０年間に３つの主要な開発が行 われ、多数の医学的状態の診断と治療に役立っており、特に、人間の解剖に応用 されている。その主要開発とは、（1）コンピュータ（Ｘ線体軸）断層撮影法 (Computer-Assisted Tomography - ＣＡＴ)、（2）核磁気共鳴映像法(Magnetic Resonance Imaging - ＭＲＩ)、および（3）陽電子放射断層撮影法（Positron E mission Tomogrphy -ＰＥＴ）である。 ＣＡＴスキャナによれば、Ｘ線が、例えば、人間の脳を通り抜けるように透過 され、コンピュータは人間の頭部の外部で検出されたＸ線を使用して、一連の画 像、基本的には、人間の脳の断面の画像を生成し、表示している。画像化される のは、脳内の未散乱の、ハードＸ線のＸ線吸収係数である。ＣＡＴスキャンは、 例えば、脳卒中(stroke)、腫瘍、および癌を検出することができる。ＭＲＩデバ イスでは、コンピュータは脳に突き当たった無線信号からのデータを処理し、生 き写しの３次元画像を組み立てている。ＣＡＴスキャンによる場合と同じように 、腫瘍、凝血塊、萎縮部位などの奇形を検出することができる。ＰＥＴスキャナ では、注入された放射性物質の位置が脳がその物質を使用したとき検出され、画 像化されている。画像化されるのは、ガンマ線ソースの位置である。これらの医 学的画像化手法の各々は、多くの異常な医学的状態を検出し、診断する上で計り 知れないほど貴重であることが実証されている。しかし、多くの点で、これらの 手法は、いずれも、以下の説明で述べている理由により完全には満足すべきもの ではない。 医学画像化手法の最適設計パラメータを設定するに当たり、最も重要な規格と して次の４つの規格がある。これらの規格はその概要を簡単に説明してから、詳 しく説明することにする。さらに、従来の手法の各々がもつ短所についても、そ の概要を説明する。第１は、非イオン化放射ソースを使用することが好ましいこ とである。第２は、メリメートルのオーダで空間解像度を達成すると、診断が容 易化されるという利点があることである。第３は、新陳代謝(metabolic)の情報 を取得することが望ましいことである。第４は、画像化情報をほぼリアルタイム で（ミリ秒のオーダで）得るようにすると、動画に似た画像を見ることができる という利点があることである。３つの従来の画像化手法で、４つの規格をすべて 同時に達成できるものは１つもない。例えば、ＣＡＴスキャナは高解像度の能力 をもっているが、イオン化放射を使用している。このスキャナは新陳代謝の画像 化の能力がなく、その空間解像度は許容ボーダライン上にある。また、ＭＲＩは 非イオン化放射を使用し、許容し得る解像度をもっているが、新陳代謝情報が得 られず、特に速度が遅くなっている。最後に、ＰＥＴスキャナは新陳代謝情報が 得られるが、イオン化放射を使用し、速度が遅く、空間解像度も許容ボーダライ ン上にある。さらに、ＰＥＴ手法は物質が注入されるために、健康な組織を冒す (invasive)。 以下では、これらの４つの規格について詳しく考察することにする。イオン化 放射に関しては、現在、イオン化放射が人体に及ぼす影響について医学界で大い に議論されている。放射レベルが現在許容限界と信じられている範囲内に収まる ようにするために、ＰＥＴスキャンは時間インターバルを狭めて行うことができ ず（多くの場合、スキャンしてから次にスキャンするまでに、少なくとも６か月 待つ必要がある）、照射量を調整する必要がある。さらに、ＰＥＴは、陽電子放 射アイソトープを作るためにサイクロトロンが必要であるため、まだ研究ツール の段階にある。空間解像度に関しては、異なる構造だけでなく、凝血塊や腫瘍な どの望ましくない状態も区別するだけの必要な粒状度（granularity）がないと 、診断が困難であることは説明するまでもない。新陳代謝情報に関しては、望ま しいことは、例えば、人間の頭部内の酸素濃度の空間的マップを作成すること、 あるいは脳内のぶどう糖濃度の空間的マップを作成することである。このような マップを作成できるようにすると、医学関係者は病気だけでなく、正常な機能に ついても知ることができる。残念ながら、ＣＡＴとＰＥＴは濃度測定を報告する が（Ｘ線スキャナでは電子、ＭＲＩでは陽電子）、新陳代謝情報を確認するため のコントラストが充分でない。つまり、ある化学物質（例えば、ぶどう糖など） を他の化学物質と区別することはほとんど不可能である。ＰＥＴスキャナは新陳 代謝情報を得る能力をもっており、この手法が近年人気があるのはそのためであ る。最後に、画像化はかなりの処理時間の後ではじめて行われるために、リアル タイム画像化は従来の手法ではほとんど不可能である。 上述した困難性と制約があるために、現在非常に関心がもたれているのは、生 体組織の陽電子拡散と吸収係数の分布の画像であって、上述した４つの望ましい 規格を満足する画像を生成する手法の開発である。従って、低強度の陽電子を使 用する手法が安全である。この手法は、光学的事象が１０ナノ秒の範囲内で発生 することから高速である必要がある。この速度であれば、多数の測定を行って平 均値をとれば、測定ノイズが減少し、他方では、１ミリ秒の速度でリアルタイム の画像化も達成することが可能になる。さらに、この手法のためのソース（発生 源）と検出機器は、適当に選択した空間パラメータを利用する再生手順に必要な 測定データを得て、これにより望ましい１ミリメータの空間解像度が得られるよ うに構成することが可能である。最後に、この手法による新陳代謝画像化は、局 所化分光法(localized spectroscopy)としての画像化を意図している場合には、 画像内の各点には拡散係数のスペクトルが割り当てられるという意味で、実現可 能になっている必要がある。 光学的画像化のための第１の提案では、コンピュータ（Ｘ線体軸）断層撮影画 像を得るために使用されているものと類似している数学的アプローチ（例えば、 後方投写アルゴリズム）が提案されている。パルス化レーザからの光はソース位 置で試料に入射され、透過光量子(transmitted photon)を受信する個所に戦略的 に置かれた検出器で検出される。最も早く到着した光量子（所謂「弾道光量子(b allistic photon)」）はソースと検出器の間を直線的に走行するものと想定され 、透過強度は数学的再生アルゴリズムで使用されている。実際には、未散乱入射 波だけが試料に埋め込まれた任意の対象物の画像を形成するのに役立つものとみ なされており、従って、「高速ゲーティング時間(fast gating time)」をもつ検 出器を、最も早く到着した光量子だけを処理するように構成するといったように 、検出プロセスからの散乱光を除去するための手法が採用されている。しかるに 、弾道光量子が指数的に減衰されることは公知であるので、試料があらかじめ決 めた値を超える厚さをもつ場合には、画像化は多くの実用分野ではほとんど不可 能である。 光学的画像化のための最近の提案は、散乱放射と拡散放射を使用して試料の内 部表現を再生する画像化システムの傾向にある。この分野の従来技術を代表する ものとして、米国特許第5,070,455 号（Singer他（Singer特許）、１９９１年１ ２月３日特許付与）がある。Singer特許に開示されているシステムは、光量子や 他の粒子などの放射を使用し、この放射は試料の内部構造によって相等の程度ま で散乱させている。このシステムでは、試料に照射し、減衰、散乱された放射 の測定は試料の外部に沿った数個所で行われている。Singer特許によれば、この 測定だけでも、十分に、試料内部の各所の領域の散乱および減衰特性を決定でき ることが明らかにされている。Singer特許の開示内容によれば、試料の内部は体 積要素（volume elements - "voxel"）の配列として表現されている。試料のモ デル内の各体積要素は数値パラメータで表された散乱特性と減衰特性をもち、こ れらは試料の内部の複数の画像を生成するようにマップ作成が可能になっている 。 試料の内部を再生するために Singer 特許によって使用されている特定の手法 を最もよく表している特徴は「反復的(iterative)」手順であることである。以 下では、その短所と欠点を指摘するために、この手順について若干詳しく説明す る。画像化データを収集した後、体積要素の散乱係数と減衰係数に初期値が割り 当てられる。このようにすると、計算プロセスの短縮に役立つが、これは数学的 最小化問題（mathematical minimization problem）の解を反復的にまたは非直 接的に求めることの特徴ともなっている。次に、このシステムは、対象物の内部 が散乱係数と減衰係数の現割当て値によって特徴づけられていた場合、試料から 放出された光の強度を計算する。次に、測定された光強度と計算された光強度と の差は、再生誤差の大きさに関係する「誤差関数(error function)」を計算する ために使用される。この誤差関数（最小化方法では「コスト関数」とも呼ばれる ）は、次に、多次元勾配降下法(multi-dimensional gradient descent metholog y)（Fletcher-Powell 最小化法など）を用いて最小化される。つまり、係数は誤 差関数の値を小さくするように修正される。 放出光強度を散乱係数と減衰係数の現割当て値に基づいて計算し、計算値と測 定値との差を比較して試料の内部の散乱特性と減衰特性の新しい近似値を得るた めのプロセスは、誤差関数が所定しきい値以下になるまで続けられる。このプロ セスからの散乱係数と減衰係数の最終値は試料内部の一連の画像を生成するよう にマップ作成され、試料内部の減衰特性と散乱特性が示される。これらは正常状 態と異常状態の両方を示しているはずである。 以上のように、Singer特許が開示している手法は、反復的最小化方法を用いて 反転により画像を再生するものである。このようなアプローチは、「アルゴリズ ム」というよりも、「ヒューリスティック（発見的）」と特徴づけた方が適して いる。なぜならば、このアプローチによれば、解の存在さえも検証または実証さ れないからである。このような状況の下では、散乱係数と減衰係数の個数はほぼ 無限であり、反復的手法を用いて決定された特定の係数が試料内部の実際の係数 であるとの保証はまったくない。さらに、このヒューリスティック法は計算量（ computational complexity：計算の複雑性）が非常に多く、これは体積要素の数 に比べて指数的であり、多くの局所的最小値をもつ困難な最適化問題の特徴とも なっている。このアプローチの計算量は再生法を画像化に使用するのに役立たな いものにしている。 従来技術において提示されている他のアプローチはSinger特許に提示されてい るものと密接な関係がある。これらのアプローチも、順方向散乱問題の間接的反 転を反復的手法によって行っており、物理的洞察力が得られるとしても、わずか である。発明の概要 従来手法の上記制約およびその他の短所と欠点は、本発明によれば、直接再生 法とこれに付随するシステムを利用して調査対象の対象物の画像を生成すること によって解消している。直接再生の公式化は画像化手法の存在と一意性（単独性 ）の両方を保証する。さらに、この直接的再生法は計算量を大幅に減少する。 本発明の広義の側面によれば、調査対象の対象物は時間領域ソースによって照 射され、拡散的に散乱する放射に主に起因する透過強度は対象物に近接する、選 択した位置で測定され、そこでは、透過強度は積分演算子(integral operator) によって拡散係数と関係づけられる。対象物を表す画像は、積分演算子に関連し て判断された、あらかじめ決めた数学的アルゴリズムを透過強度の測定値に対し て実行することにより直接的に再生される。さらに、異なる波長の放射は画像化 を局所化分光法(localized spectroscopy)として行う。 本発明の構成と動作の理解を容易にするために、以下において、添付図面を参 照して実施例を詳しく説明する。図面の簡単な説明 図１は弾道限界内で対象物を収めている試料を光が通り抜ける様子を示す図で ある。 図２は拡散限界内で対象物を収めている試料を光が通り抜ける様子を示す図で ある。 図３は時間領域源をもつ複数の検出位置の場合の、光が試料を通り抜ける様子 を示す図である。 図４は拡散が一定の場合、吸収が一定の場合、および拡散と吸収が変動する場 合のそれぞれの、媒質内に埋め込まれた対象物を示す図である。 図５は複数の時間的瞬時における拡散核を示すプロット図である。 図６は試料内の対象物を示すと共に、本発明の直接再生法を使用して対象物を 再生する様子をシミュレーションで示す図である。 図７は本発明による画像化システムの一実施例を示すハイレベル・ブロック図 である。 図８は本発明の方法を示すハイレベル・フロー図である。 図９は調査対象の対象物の拡散を計算する一方法を示すフロー図である。 図１０は再生された対象物の例を示す図である。 なお、これらの図において同一のエレメントは、同一の参照符号を付けて示し てある。詳細な説明 本発明の詳細な説明を全体的な観点から展望し、本明細書に開示され、請求の 範囲に記載されている技術からの出発点に焦点を当てるために、理解に役立ち、 参考になることは、まず、本発明による主題に関係する、いくつかの基礎原理を 説明することによって、本発明が作用する画像化環境の基本的理解を得ることで ある。従って、本明細書の最初の部分では、本発明の主題に関係する画像化シス テムを中心にハイレベルに説明することにする。この方法によると、本発明の種 々の具体的側面の理解を助ける概念と用語を理解できるという利点がある。この 概要説明に続いて、本発明をシステムから見た側面と、それに付随する方法論を 具体的に説明することにする。 本発明の概要説明 光の多重散乱(multiple scattering)は、光学的画像化(optical imaging)の基 本的な物理的妨げとなっている。本発明はこの現象を取り扱うことを主題としい るが、拡散光には、高散乱媒質の光学的拡散を画像化するのに十分な情報が含ま れているという、驚くべき結果が得られている。この結論は、拡散限界内の多重 散乱に応用可能である、逆散乱理論を積分方程式で公式化することによって得ら れたものである。高散乱媒質の光学的拡散係数を画像化するために考案されてい た第１直接再生方法は、この表現方法を用いると理解が容易になる。画像形成の ために未散乱（弾道）光量子を利用する手法とは対照的に、本発明の方法による と、平均散乱平均自由行程(average scattering mean free path)に比べてサイ ズが大きい対象物の画像化が可能になる。 衝突光をもつ多くの対象物の、普通の不透明または曇った外観は、多重光散乱 の現象によって説明することができる。（なお、ここで注意すべきことは「対象 物(object)」とは調査対象となっているものを物理的に表わしたものであるとい ったように、本明細書で用いられている用語は概念化されたものである。例えば 、対象物は単独で存在することもあれば、試料またはサンプルに埋め込まれてい ることもある。いずれの場合も、対象物に関する説明内容のコンテキストは、そ のコンテキストにおける「対象物」という概念用語の意味を明らかにして説明さ れている。本発明の開示事項と教示事項は、高散乱媒質に埋め込まれている広が りをもつ対象物(extended object)を画像化する問題を取り扱っている。この問 題の解を求める中心は、上述した積分方程式を公式化することである。拡散的に 透過される光は直接画像再生のための十分な情報を含んでいるので、この問題は 、基本的に閉形体の解を求め易くする、扱いやすい形体で表すことができる。 つまり、このことは、すべての短所と落とし穴のある、反復的／最小化タイプの 再生に頼らないで済むことを意味する。 最も基本レベルの直接再生プロセスを理解し易くするために、まず、直接再生 が応用可能である単純化されたシステムについて説明する。つまり、このシステ ムでは、波長λの平面光波（光量子）は位置に依存する光学的拡散係数によって 特徴づけられる空間拡張対象物を含んでいる線形寸法Ｌのサンプルに入射する。 幅Ｌは衝突入射波と位置合わせされる。さらに、λに比べてそのサイズが大であ る粒子によって光量子が散乱されると想定すると、その散乱は１^*と名づけた走 行平均自由行程で表される。平均自由行程はその方向がランダム化される前に光 量子が走行する平均距離を表している。単一散乱状況では、つまり、１^*＞＞Ｌ のときは、入射波の大部分がサンプルから出るとき未散乱であることが観察され ているので、拡散対象物の投写画像を形成するために使用できる。この効果は図 １に示されている。図１において、波長λの光線１０１は拡散対象物１２０を含 んでいるサンプル１１０の前面１０５に突き当たり、そこでは、サンプル１００ を通り抜けてサンプル１００の背面１０６から出た光線はトレース１３０で表さ れた投写画像を形成する。トレース１３０で表された透過強度は未散乱波の伝播 方向上の光学的拡散係数の線積分と関係づけられている。これにより、拡散係数 の所謂Radon 変換が得られる。Radon 変換を反転すると、拡散係数を回復するこ とが可能であるので、拡散対象物１２０の画像が再生される。すでに上で言及し たように、商用化されている画像化手法はいずれも、この単純な物理原理を基礎 にしている。 多重散乱状態では、つまり、１^*＜＜Ｌのときは、波はサンプルを通り抜ける 間になん度も散乱する。この状態にあるとき、λ＜＜１^*のときは、単一光量子 の通路は拡散ランダムウォーク(diffusive random walk)で表すことができ、こ ｃ／ｎはサンプルの媒質内の光の速度である。未散乱、つまり、弾道光量子は静 的透過係数Ｔ_ball〜ｅｘｐ（−Ｌ／１^*）と共に指数的に減衰する。透過強度に 主に寄与するものは、拡散透過係数Ｔ_diff〜１^*／Ｌをもつ拡散光量子によるも のであり、これは、コヒーレントな照射によっても、サンプルの単純な画像を含 んでいない複雑な干渉パターンを形成する。このようなパターンは図２に示され ている（この図は、図２の物理システムが１^*＜＜Ｌであり、図１の１^*＞＞Ｌと 異なることを除けば、図１とほぼ同じ絵図になっている）。図２において、波長 λの光線２０１はサンプル２１０の前面２０５に突き当たり、最終的には背面２ ０６からサンプル２１０を出ていく。吸収対象物２２０からトレース２３０が得 られ、これは背面２０６から出た複雑な透過光パターンを表している。本発明に よれば、トレース２３０で表されたこの複雑パターンの情報を利用して対象物の 位置を突き止めることによって、多重散乱状態で光学的画像化を行う閉形体手順 が考案されている。 事実、従来技術で頻繁に指摘されているように、弾道光量子は歪みが最も少な い画像情報を透過するのに対し、拡散光量子は画像情報の大部分を失っている。 この理由から、いくつかのエレガントな実験的手法は、光学的ゲーティング(opt ical gating)、ホログラフィ(holography)、または光学的拡散による拡散光量子 のフィルタリングによって、弾道光量子の寄与を選択するように設計されている 。しかし、弾道光量子だけに頼るどの手法にも、内在する物理的制約がある。こ のことは、弾道光量子の指数的減衰を拡散光量子の緩慢な代数的減衰と対比して 考えれば、容易に理解されることである。具体的には、サンプルのサイズＬが１^* に比べて十分に大きければ、Ｔ_ballは、指数的に測定可能なしきい値以下にな る（例えば、１^*が約０．５ｍｍならば、減衰は２ｃｍだけｅ^-40に比例する）。 従って、これまでは再生がほとんど不可能であると信じられていた重要で貴重 な画像を再生する可能性は、画像再生のために多重散乱拡散光量子を採用するこ とによって弾道画像化の制約を解消する強力な動機となっている。基本的な物理 原理から、拡散透過光の干渉パターンからのかかる再生が達成可能になったのは 、かかる再生が２つのパラメータ、つまり、高散乱システムの吸収係数と拡散係 数によって特有な方法で決定されるからである。この最も一般的な問題に対する 解は、ある意味では、拡散画像化に対するRadon 変換を概念化したものである。 本発明によれば、本明細書に説明されているように、拡散透過係数は拡散係数の 積分と関係づけられている。これにより、画像は、この積分演算子を参照す る適当なアルゴリズムを使用すると再生することができる。弾道方法とは対照的 に、その結果の再生アルゴリズムは、サイズＬが１^*に比べて大であるサンプル を画像化するために使用できる。 拡散画像化のための関数理論の基礎 弾道効果と拡散効果を分離することは、最も自然には、時間分解パルス伝播ア プローチ(time-resolved pulse propagation approach)を使用して行われる。こ のアプローチでは、高散乱媒質内を伝播する光学パルスの透過係数の時間依存性 が観測される。短時間スケール上の透過は単一散乱状況における光量子の弾道走 行を測る尺度である。これに対して、長時間スケール上の透過（拡散時間Ｔ_D＝ Ｌ²／Ｄに相当する）は、多重散乱状況における拡散走行を測る尺度である。こ れらの短時間スケールと長時間スケールは図３に示されているが、これは検出器 位置が複数の場合である。図３において、検出された光量子の透過強度はその量 が縦軸に、時間が横軸に示されている。ソースと検出器のペアは、図２に示すよ うに、サンプルに近接する３つの異なる位置に置かれ、各位置からは、透過強度 と時間との対比関係が得られるようになっている。従って、例えば、曲線３１０ はソースと検出器のペアの第１位置に対応し、曲線３２０は第２位置に対応し、 最後に、曲線３３０は第３位置に対応している。ここで指摘すべきことは、３つ の曲線はいずれも、長時間インターバルの後の勾配がほぼ同じになっていること である。 ｒ₁の個所のソースから検出器ｒ₂へ時間ｔで伝播するパルスの、時間依存拡散 強度透過係数Ｔ（ｒ₁，ｒ₂，ｔ）は次の指数法則で与えられる。 上式において、Ｄ（ｒ）は均一バックグラウンドから離れた拡散の変動（Ｄ_o と呼ぶ）である。図４の図はこれらの関係を示している。つまり、図４(ii)に示 すように、対象物４００は吸収α_oと拡散定数Ｄ_oをもつ媒質に埋め込まれており 、他方、対象物４００の方は拡散Ｄ（ｒ）≠０と吸収α_oもっている（図には、 対象物４００を取り囲むｊ番目のソース４０２とｊ番目の検出器４０３も示 されているが、これについては以下で詳しく説明する）。積分方程式の拡散核は 次式で与えられる。 上式において、Ｇ（ｒ₁，ｒ₂；ｔ）は物理的に適当な境界条件を満足する拡散伝 播係数である。自由空間における拡散核の解析式を求めると、次の結果が得られ る。 ただし、 方程式（１）は弱い拡散または拡散の小さな空間的変動の限界内で有効である。 拡散核の輪郭プロットは図５に示されているが、種々の時間瞬時は、約０．００ １γ_D、０．１γ_D、および１．０γ_Dに対応している。図示のように、非常に初 期の時間（図５(i)）では、核５０１はソースと検出器を結ぶライン上に大部分 が集中している。これは、ほぼ弾道的な光量子からの核に支配的に貢献している ことを示している。長い時間（図５(ii)と図５(iii)）では、多重散乱状況にお いて、核（それぞれ５０２と５０３）は長い通路上の光量子からの寄与を含んで いる。従って、これらの核から、拡散限界におけるパルス伝播アプローチで光量 子走行の物理的画像が得られる。 積分方程式（１）は光量子拡散画像化の基本的積分方程式と呼ばれる。この積 分方程式はパルス伝播アプローチにおける透過係数を、吸収係数および拡散係数 の両方と関係づけている。 光量子拡散画像化の中心となる問題は、パルス伝播実験においてソースと検出 器のペア群の透過測定から画像を再生することである。適当な再生手法を説明す るためには、基本的積分方程式（１）の解を求める必要がある。この積分方程式 は第一種のFredholm方程式である。この種の方程式は特異(ill-posed)であるの が代表的であり、その解を求めるためには、正規化法を導入する必要があること は周知である。式（１）のこのような正規化された解は、特異値分解によって求 めることができ、次式によって与えられる。 ただし、 は、ν_D(ｒ；ｒ₁ｒ₂ｔ)）の正規化、一般化された反転であり、積分は対象物を 取り巻く測定表面上にわたって行われる。ここで、σ_n，ｆ_n，ｇ_nは特異値と、 ν，ν^*νｆ_n＝σ_n ²ｆ_n，νｆ_n＝σ_nｇ_n対応する単数関数を示し、βは正規化パ ラメータ、およびＲ_βは適当の正規化定数である。代表例として、Ｒ_βはＲ_β( σ)＝σ／（β＋σ²）となるようにとられ、小さな特異値がスムーズにカットオ フされるようになっている。式（４）と（５）から光量子拡散画像化における画 像再生方法の正式の解が絵得られる。ここで注目すべきことは、上記再生方法の 空間解像度は１^*によって、つまり、拡散近似化が有効である長さスケールによ って制限されていることである。 式（４）から得られる明示的反転公式の存在は画像再生アルゴリズムを開発す る上で重要であることは明らかである。しかし、反転公式は、有限数のソースと 検出器のペアだけからの透過測定が使用できるように適応させる必要がある。こ の問題を解決する１つの方法は、式（５）における正規化された特異値分解を直 接的に数値で実現することを考慮することである。ここでは、積分方程式（１） は、不連続の定数関数をもつコロケーションなどの、適当な離散化方法によって 線形方程式の系に変換される。この方法では、複数のソースと検出器のペアから 透過係数の測定値を得る必要がある。各ペアは複数の時間点にも寄与する。従っ て、ソースと検出器のペア／時間点の組合せの数は、少なくとも再生画像内のピ クセル数と同数にする必要がある。重要なことは、この種の実空間再生アルゴリ ズムの計算量はＯ（Ｎ³）であることを認識することである。ここで、Ｎは再生 画像内のピクセルの数である。なお、これは、関連する数値特異値分解の複雑性 にすぎない。 シミューレートした透過データに対するコンピュータ処理を使用した２次元対 象物の直接再生は図６に示されている。図６(i)に示すように、対象物６２０は サンプル６１０に埋め込まれている。対象物６２０の形状は既知であるので、前 面６０５に突き当たる光量子に起因する背面６０６からの光粒子の放出を数学的 に記述することにより、それを計算することが可能である。つまり、拡散画像化 における所謂順方向問題を解くことが可能である。背面６０６の近くで検出され た光量子の透過強度が与えられているとき、所謂逆方向問題を解くと、対象物６ ２０の画像を直接に再生することができる。このように直接に再生された画像は 図６(ii)に対象物６２１で示されている。ここで強調すべきことは、透過データ はシミューレートされているが、画像を再生するためにコンピュータ処理で使用 されるアルゴリズムは、実際のサンプルから検出された透過強度の実際の測定値 を処理するために使用されるアルゴリズムそのものであることである。このよう なシミュレーションによると、例えば、ノイズが透過強度データに及ぼす影響お よび再生手法がソースと検出器のペアの各位置に対してもつ感度を調べることが できる。 発明の詳細な説明 Ｉ．システム 図７にハイレベル・ブロック図で示すように、システム７００は、対象物に突 き当たる光量子に応じて対象物から放出される透過光量子の測定を使用して、対 象物の画像を生成する直接再生画像化システムである。具体的には、対象物７１ ０は調査対象として示されている。システム７００は対象物７１０に光量子（フ ォトン）を照射する時間領域ソース７２０、対象物７１０に近接する１つまたは 複数の戦略的位置で対象物７１０から放出される光量子の透過強度を測定するデ ータ収集検出器７３０、検出器７３０の位置と時間領域ソース７２０との相対位 置を制御する位置コントローラ７４０、およびコンピュータ・プロセッサ７５０ から構成され、コンピュータ・プロセッサは関連の入力デバイス７６０（例えば 、キーボード）と出力デバイス７７０（例えば、グラフィカル・ディスプレイ端 末）を備えている。コンピュータ・プロセッサ７５０はコントローラ７４０から 位置情報を、検出器７３０から測定した透過強度を入力として受け取る。 システム７００の図示実施例によれば、光量子ソース７２０はCoherent Corp. （コヒーレント社）から提供されているモデルMIRA-900P チューナブル・レーザ を使用している。（このレーザは、実際には、２つの他の補助デバイスが付随し ている。(1)７８ＭＨｚパルスレートを５ＭＨｚにスローダウンさせる音響光学 パルス・ピッカ -- このデバイス例としては、Coherent Corp.から提供されてい るモデル９００がある -- および(2)MIRA-900P をポンピングする別のレーザ -- このポンプ・レーザの例としては、Coherent社から提供されているモデルINNOV A-45 がある -- である。） データ収集検出器７３０は、例えば、Hamamatsu Corp.(浜松社)から提供され ているストリーク・スコープ・モデル64-334-02 などの光量子検出器を使用して いる。 位置コントローラ７４０は光量子検出器７２０および／またはデータ収集検出 器７３０が複数のレーザまたは光量子検出器から構成されているときに使用され 、複数のレーザのどれをある時間期間の間励起させるか、あるいは複数の光量 子検出器のどれをある時間インターバル期間で起動させるかを制御する。以下で 詳しく説明するように、実際に実行された直接再生画像化手法では、対象物７１ ０を取り巻く複数のソース・検出器位置によって透過光量子強度を測定すること が必要になることがよくある。便宜上、必要な透過強度データの生成は、Ｐ個の レーザ・ソースとＱ個の光量子検出器の配列をもつことにより迅速に行われるも のとする。従って、光量子ソース７２０は、最も一般的な実現例では、対象物７ １０の周囲を取り巻くように戦略的に配置されたＰ個のレーザ・ソースまたは類 似ソースで構成することができる。同様に、データ収集検出器は、最も一般的な 実現例では、対象物７１０の周囲を取り巻くように戦略的に配置され、Ｐ個のソ ースと協働関係にあるＱ個の光量子検出器で構成することができる。 コンピュータ７５０は直接再生アルゴリズムを実現するコンピュータ・プログ ラムを格納している。具体的には、その格納されたプログラムは測定した透過デ ータを処理し、積分演算子を参照する、あらかじめ決めた数学的アルゴリズムを 使用して調査対象の対象物の画像を生成する。この積分演算子は、拡散係数Ｄ（ ｒ）が係りをもっている。コンピュータ７５０によって行われる処理は、本明細 書の次のセクションの方法論の説明の中心になっている。 II．方法論 計算モデル 方程式（１）で表された基本的積分方程式は、以下に再度示すように、 第１種のFredholm方程式の形体になっている（具体的には、本明細書では、Shor tland 第２時間領域積分方程式(Schotland Second Time-Domain Integral Equat ion)と呼ぶことにする）。このような方程式はＫｆ＝ｇ、または次式でかかれて いるのが代表的である。 上記において、ｆ，ｇは適切に選択された関数空間の要素である。式（６）は、 (a)解が求められない、(b)固有の解が存在しない、または(c)解がデータに連続 的に依存しないとき、特異(ill-posed)といわれる。最後のケース(c)が安定の特 異な問題を数値的に調べる上で最も関心があるのは、数値不安定になるからであ る。これは、データがあいまいに既知であったり、あるいは測定の不正確性やノ イズのように、画像化のための測定を行うとき起こる統計的不確実性が問題であ る場合に特に重要である。これら特異な問題を調整する方法がいくつかある。ま ず第１に、解が存在しない場合、‖Ｋｆ−ｇ‖の最小値が解と定義される。非一 意性（非単独性）は最小ノルムをもつ最小値を選択することによって処理される 。最後に、連続性は解を求める手順に「正規化」を導入することによって復元さ れる。 最小ノルムをもつ最小値の解を求めると、式（６）に関する「正規方程式」が 得られる。正規方程式は次のような形式になっている。 上式において、Ｋ^*はＫの随伴であり、Ｋ^*Ｋが自己随伴であるとのプロパティ（ 性質）が採用されている。従って、式（６）におけるｆの解は次のような形式に なっている。 式（８）から、 はＫの「一般化反転（generalized inverse）」と呼ばれる。 特異値分解 ＫがＨ₁からＨ₂へのマッピング（写像）が行われるようになっていれば（ただ し、Ｈ₁とＨ₂はHilbert(ヒルベルト)空間である）、Ｋ^*Ｋは自己随伴の正の作用 素(positive operator)である。Ｋ^*Ｋの固有関数(eigenfunction)と固有値(eige nvalue)をそれぞれ｛ｆ_n｝と｛σ_n ²｝で表わすと、次の関係式が得ら れる。 ｛σ_n｝はＫの特異値(singlar value)である。また、｛ｆ_n｝はＨ₁の基底を形成 する。特異値はσ₁ ²≧σ₂ ²≧...≧０と配列され、ここでは、重複度がカウント され、無限大重複度で０が現れることができる。 ｛ｇ_n｝が次式によって定義されていれば、 ｛ｇ_n｝はHilbert 空間Ｈ₂の基底である。さらに、次式が得られる。 Ｋの特異値分解を求めるために、Ｋを次の形式で置き、 次の恒等式を使用する。 および と、次式が導かれる。 式（１６）はＫの「特異値分解(singular value decomposition)」と呼ばれる。 ここで、式（１６）の特異値分解を使用すると、式（９）の一般化反転Ｋ^*の 形式が得られる。式（１６）の結果として、 および が得られ、式（１７）と（１８）を式（９）に代入すると、次のようになる。 次に、式（１９）を使用すると、ｋｆ＝ｇの解はｆ＝ｋ⁺ｇとなり、その形式は 次の通りである。 σ_nの一部がゼロであれば、Ｋ⁺は不明確に定義され、具体的には、連続していな い。この異常を解決するために、正規化法が導入されている。 正規化 特異値分解を条件づけるために、次式が定義されている。 上式において、正規化値Ｒ_β(σ)は次のプロパティ（性質）をもっている。 （i） Ｒ_β(σ)＝１／σ β→０⁺の場合； （ii） Ｒ_β(σ)〜１／σ σ＞＞０に対して（β＞０のとき）； （iii）Ｒ_β(σ)→０ σ→０の場合（β＞０のとき）。 例えば、２つの自然選択（他の選択も可能）には、次のものがある。 （a）Ｒ_β(σ)＝１／σに対してσ＞β；さもなければＲ_β(σ)＝０； （b）Ｒ_β(σ)＝σ／（β＋σ²）。 （１つの代表的なヒューリスティック基準では、β〜０（σ₁）をセットしてい る）。 従って、式（８）の解は次のように書くことができる。 ただし、 Schotland の第２方程式の数値解 第１種の一般的 Fredholm 方程式の公式の解を求めるための上記開発は、特異 値分解と正規化の手法を含めて、式（１）の数値解を実現するために適用するこ とができる。 正式の解は、要約すると、次のように式（４）と（５）で与えられる。 ただし、 Ωで示した３次元対象物については、Ｐ個のソースとＱ個の検出器が対象物を プローブ（検査）するために使用されるものとする。これらのソースは対象物の 周囲を取り巻くように間隔が置かれ、適切な個所に置かれた検出器はソースと一 緒に動作するようになっている。説明を簡単にするために、以下の説明では単一 の時間点が考慮されている。一般的に、その結果は複数の時間点のケースにも容 易に拡張が可能である。ｉをｉ＝１，２，... ，Ｐおよびjをｊ＝ＱとするＰ個 のソースとＱ個の検出器に対応するインデックス（指数）であるとすると、ある 時間のとき、式（１）は次のようになる。 ここで、ν_DとＤ（ｒ）はΩを、対象物をカバーしている"boxel"(つまり、ほ ぼ等辺をもつ体積要素)Ｂ_m，ｍ＝Ｍに分解することによって離散化される。次に 、粒状度は、Ｄ（ｒ）とν_Dが各ボックスで一定になっているものとする。式（ ２７）を標準形式に書き直すと、次のようになる。 および 上記において、｜Ｂ_m｜はboxel の体積であり、ａ_mはｍ番目のboxel の中間のα （ｒ_n）の強度の強度である。次に、これらの定義を使用すると、式（２７）は 次のようになる。 ｍ＝１，２，... ，Ｍ；ｉ＝１，２，... ，Ｐ；およびj＝１，２，... ，Ｑと する。 行列形式では、式（３１）はＡ_a＝ｂと表される。ただし、Ａは（ＰＱ×Ｍ）行 列であるので、式（３１）からＭ個の未知方程式とＰＱ個の方程式が得られる。 従って、少なくともboxel(M)と同数のソースと検出器のペア（ＰＱ）が必要であ る。好ましいことはＰＱ＞Ｍにするか、あるいは各ソースと検出器のペアごとに 複数の時間点を使用して、式（３１）を「確定的」にすることである。Ｋ個の時 間点があれば、行列Ａは（ＫＰＱ×Ｍ）行列になる。実際には、ＫＰＱ＝３Ｍが 代表的である。 行列公式化に適用される特異値分解の解は周知の手法である。例えば、特異値 分解の手法は、文献「数値解法(Numerical Recipes)」（Press、Flannery、Teuk olsky およびVittering 共著、1986、Cambridge University Press，Cambridge ，England）に記載されている。特異値分解を実現する商用化ソフトウェア・パ ッケージはInteractive Data Language(IDL)という名称で、Research Systems I nc.(Denver，Colorado)から提供されており、実用化されている。ＩＤＬは科学 計算、特に画像処理アプリケーションを目的に設計されている。ＩＤＬでは、上 記のＳ行列の観点から見た"SVD［Matrix］"形式(例えば、上記Ａマトリックスに 換算してSVD[A])のサブルーチンに似たコールから特異値が戻され、この値を使 用して線形方程式の系を解くことが可能になっている。 特異値分解が行われた後、式（２６）による正規化を行えば、行列Ａの正規化 された一般化反転（Ａ⁺と呼ぶ）を簡単に得ることができる。離散化されたSchot land の第３積分方程式の解はａ＝Ａ⁺ｂとなる。 フロー図 前のセクションで説明したきた方法論は、図７に示すシステムの図示実施例の 観点から見た図８のハイレベル・フロー図８００に示されている。図８を参照し て説明すると、制御ブロック８１０によって処理が行われると、時間領域光量子 ソース７２０とデータ収集システム７３０が作動して、光量子ソース７２０によ って対象物７１０から放出されたエネルギを測定する。これらの測定値はバス７ ３１を経由して収集システム７３０からコンピュータ・プロセッサ７５０に渡さ れる。次に、処理ブロック８２０が呼び出されて、事前に計算され、ストアされ ていた、式（２）で表された拡散核を取り出す。これを受けて、処置ブロック８ ３０が動作し、式（６）〜（３２）に関連して説明してきた直接再生アルゴリズ ムを実行することにより、拡散係数Ｄ（ｒ）を決定する。最後に、処理ブロック ８４０に示すように、拡散係数Ｄ（ｒ）に対応する再生画像はユーザが決めた形 式で出力デバイス７７０へ送られる。出力デバイス７７０は、例えば、ディスプ レイ・モニタまたはより精細な３次元ビデオ・ディスプレイ・デバイスにする ことが可能である。 ブロック８２０に示す直接再生を行う方法の一例は、図９のハイレベル・フロ ー図９００に示されている。具体的には、処理ブロック９１０に示すように、最 初のステップでは、拡散核行列Ａ（離散化によって求められたＡ行列、式（２９ ）のｍ＝１，２，... ，Ｍ；ｉ＝Ｐ；ｊ＝１，２，... ，ＱのときのＡ_ij ^m）が 作られる。次に、処理ブロック９２０が呼び出され、拡散行列Ａの特異値分解を 計算する。そのあと、処理ブロック９３０が実行され、正規化され、一般化され た反転Ａ⁺を生成する。最後に、ブロック９４０が呼び出され、解ａ＝Ａ⁺ｂを求 める。ただし、ａは拡散係数の離散化値を表している。 図１０を参照して説明すると、図は実例の対象物の再生を示している。具体的 には、図１０（ａ）は１０ｃｍ×１０ｃｍの試料（標本）１００１に埋め込まれ たオリジナルの２次元対象物１００２を示している。図１０（ｂ）は直接再生と ノイズとの相対的不感度を示す、０．１％の付加的ガウス・ノイズが存在すると きの、画像１００２１の直接再生を示している。この例では、α_o＝１．０ｃｍ^{- 1} およびＤ_o＝１．０ｃｍ²ｎｓ^-1であり、輪郭レベル１００３、１００４、１０ ０５、および１０１３、１０１４、１０１５はそれぞれｃｍ^-1の単位で測定され た吸収係数と、ｃｍ²ｎｓ^-1の単位で測定された拡散係数の変動を示している。 透過係数はモンテカルロ・シミュレーションを使用して得られている。 上述してきたシステムおよび方法論は拡散核（式（３））の自由空間モデルを 使用し、拡散核は事前に計算され、コンピュータ・プロセッサ７５０にストアさ れていて、再生プロセス期間にリコールされるようになっている。これが適して いるのは、物質（例えば、Intralipidの名称で市販されている医療製品）が満た された薄く、ラバー状の容器などの構成物によって対象物７１０が取り囲まれて いるときであり、このようにすると、対象物の外部に空間的広がりが得られ、対 象物を取り巻く自由空間条件が実効的に得られる。対象物の実際の境界（例えば 、脳の画像化のときの人間の頭蓋骨）は、直接再生法によって決定された別の形 状となる。Intralipidが有用であるのは、これがコロイド物質であり、そこでは 、０．５ミクロンから２ミクロンの範囲の粒子は懸濁され、この物質はパッ ケージに入っているので、急速に劣化しないからである。さらに、この物質の１^* が測定可能である。 上述してきた実施例は本発明による原理の応用の単なる例示であり、本発明の 精神と範囲内から逸脱しない限り、他の実施例も可能であることはもちろんであ る。従って、これまでに説明してきた方法論は例を示して説明してきた具体的ケ ースに限定されるものではなく、請求の範囲の記載に属するかぎり、他の実施例 も可能である。

Claims (1)

1. 【特許請求の範囲】 １．対象物の拡散画像を生成する方法であって、 対象物に放射型の時間領域ソースを照射するステップと、 拡散的に散乱された放射に主に起因する透過強度を測定するステップであって 、そこでは、前記透過強度は積分演算子によって拡散係数および拡散核と関係づ けられているステップと、 前記積分演算子に関連して決定された、所定の数学的アルゴリズムを該透過強 度に対して実行することにより画像を直接的に再生するステップとを有すること を特徴とする方法。 ２．請求項１に記載の方法において、対象物に照射するステップは、対象物に異 なる波長で連続的に照射するステップを含むことを特徴とする方法。 ３．請求項１に記載の方法において、画像を直接に再生する前記ステップは、拡 散核を計算するステップを含むことを特徴とする方法。 ４．対象物の拡散画像を生成するシステムであって、 対象物に照射する放射型の時間領域ソース手段と、 拡散的に散乱された放射に主に起因する透過強度を測定する検出手段であって 、そこでは、前記透過強度は積分演算子によって拡散係数および拡散核と関係づ けられている検出手段と、 前記積分演算子に関連して決定された、所定の数学的アルゴリズムを該透過強 度に対して実行することにより画像を直接に再生する手段とを具備したことを特 徴とするシステム。 ５．請求項４に記載のシステムにおいて、対象物に照射する前記放射ソース手段 は、対象物に異なる波長で連続的に照射する手段を含むことを特徴とするシステ ム。 ６．請求項４に記載のシステムにおいて、画像を直接に再生する前記手段は拡散 核を計算する手段を含むことを特徴とするシステム。