顕微鏡画像を強調する方法及びシステム
或る照明方向でサンプルに光を当てることによりサンプルを照明し散乱光を捕捉することによって形成された顕微鏡画像を用いて強調画像を作成する。これは、画像の構成部分の強度を、サンプルの複数のそれぞれの要素における散乱パラメータのそれぞれの値に関連させる式を用いて行われる。散乱パラメータは、放射係数ρｅｍとすることもできるし、そうでない場合は吸収係数ρａｂに等しくすることもできる。この式を解いて、散乱パラメータの値を見つける。散乱パラメータを用いて、強調画像、例えば、散乱パラメータ自体の変動をマッピングする画像を構成する。散乱パラメータが正確に見つけられることを条件として、強調画像は、不均一な光の減衰及び散乱に起因した劣化を原画像よりも受けにくくなる。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、顕微鏡を用いて取得された画像である顕微鏡画像を強調する方法及びシステムに関する。特に、強調画像は、不均一な光の減衰及び散乱による劣化の影響がより少ない画像である。
【背景技術】
【０００２】
顕微鏡法［１］は、生物学の重要な光学撮像技法である。２光子励起顕微鏡法及び単一平面照明顕微鏡法等の多くの顕微鏡技法があるが、共焦点顕微鏡法［１］は、バイオイメージングにとって最も重要なツールの１つとなってきている。共焦点顕微鏡法では、焦点ずれの光が、ピンホールの使用により取り除かれる。入射照明光は、ピンホールを通過し、サンプル内の小さな領域上に合焦され、この領域で、入射照明光はサンプルによって散乱される。入射照明光と同じ経路に沿って進む散乱光のみが、戻ってピンホールを通過し、そのような光は、光電子増倍管等の光検出器において再び合焦され、これによって、画像が生成される。共焦点顕微鏡を通じて取得された画像は、従来の広視野顕微鏡によって作成される画像よりも鮮明である。しかしながら、光減衰効果による劣化は、共焦点顕微鏡法では深刻である。共焦点顕微鏡法における基本的な問題は光透過問題である。入射光は、散乱されるので減衰し、したがって、厚いサンプルを透過することができない。その結果、サンプル内の深い領域から取得された画像は、サンプルの表面に近い領域から取得された画像よりも指数関数的に暗く見える。光透過の難題は、共焦点顕微鏡法に限ったものではない。単一平面照明顕微鏡法及び広視野顕微鏡法等の他の光学顕微鏡技法も同じ問題に悩まされている。古典的な空間不変デコンボリューション手法［２］、［３］、［４］は、顕微鏡法撮像のこの問題にうまく対処することができない。
【０００３】
上述した問題を解決する試みは、レーザー出力を増加させるか又は光電子増倍管の感度を増大させるかのいずれかによって行われてきた［２０］。双方の技法は、不十分であり、レーザー出力の増加は光退色効果を加速させるのに対して、光電子増倍管の感度の増大は雑音を追加するという欠点を有する。Umesh Adiga及びB. B. Chaudhury［２１］は、画像スタックの深度に沿った光減衰を考慮に入れて画像を復元するために、単純なしきい値処理方法を用いて前景から背景を分離することを議論している。この技法は、軸方向解像度を改善するために、画像ボクセルが、画像スタックに画像スライスを仮想的に挿入するＸＯＲ輪郭形成（XOR contouring）及びモーフィングに基づいて等方性である（これは共焦点顕微鏡法には当てはまらない）と仮定する。
【０００４】
表面的には無関係な技術分野が屋外撮像の分野である。その分野内では、大気エアロゾルに起因した劣化画像の復元の問題が、監視、ナビゲーション、追跡、及びリモートセンシング［５］、［１０］等のその重要な用途のために広範囲にわたって研究されてきた［５］、［６］、［７］、［８］、［９］、［１０］、［１１］、［１２］、［１３］。これらの劣化画像の復元に用いられる復元技法と類似の復元技法は、水中ビジョン［８］、［９］において、具体的には水中ケーブル及び水中パイプラインの監視等について、新しい用途も創出している。さまざまな復元アルゴリズムが、均一な媒体を通る光減衰及び光散乱（大気散乱光）の物理モデルに基づいて提案されている。そのような画像復元アルゴリズムに関するこれまでの研究の１つ［５］は、シーンの奥行きの正確な情報を必要とする。その後の研究は、シーンの奥行きの必要性を巧みに回避したが、必要とされる情報を回復するのに複数の画像を必要とする［８］、［９］、［１０］、［１４］。Narasimhan及びNayar［１０］、［１１］、［１２］、［１３］は、１つの劣化画像のみから全ての必要な情報を抽出する対話型アルゴリズムを開発した。この方法は、大気散乱光の色（airlight color）及び「良好な色の領域」の手動選択を必要とする。これらの復元技法に関する基本的な問題は、雑音の増幅である。この基本的な問題を取り扱う試みは、Kaftory他［１４］によって提案された変分手法（variational approach）における正則化項を用いることで行われる。
【発明の概要】
【発明が解決しようとする課題】
【０００５】
本発明は、上記問題を克服することができる、画像を復元する方法を提供することを目的とする。
【課題を解決するための手段】
【０００６】
概括的に言えば、本発明は、或る照明方向でサンプルに光を当てることによりサンプルを照明し散乱光を捕捉することによって形成された顕微鏡画像を用いて強調画像を作成することを提案する。これは、顕微鏡画像の構成部分の強度を、サンプルの複数のそれぞれの要素における散乱パラメータのそれぞれの値に関連させる式を用いて行われる。散乱パラメータは、放射係数ρｅｍとすることもできるし、そうでない場合は吸収係数ρａｂに等しくすることもできる。この式を解いて、散乱パラメータの値を見つける。散乱パラメータの値を用いて、強調画像、例えば、散乱パラメータ自体の変動をマッピングする画像を構成する。
【０００７】
散乱パラメータの値が正確に見つけられることを条件として、強調画像は、不均一な光の減衰及び散乱に起因した劣化を原画像よりも受けにくくなる。
【０００８】
上記式は、各要素の散乱パラメータの値を、入射光の方向に沿った要素の散乱パラメータの値の関数として与えることができる。この場合、散乱パラメータの値は、ロケーションについて連続的に、さらに照明方向でも連続的に見つけることができる。
【０００９】
上記式は、照明方向に平行な一組の要素を取り囲む領域にわたって規定された散乱パラメータの平均値を用いることができる。
【００１０】
本発明は、代替的に、そのような方法を実行するコンピュータシステムとして表現することもできる。このコンピュータシステムは、画像を取得するデバイス、例えば顕微鏡と統合することができる。本発明は、上記方法のステップを実行する、コンピュータシステムによって動作可能なプログラム命令を含む、有形のコンピュータ媒体に記録されるようなコンピュータプログラム製品として表現することもできる。
【００１１】
次に、単なる例示の目的で、以下の図面を参照して本発明の一実施形態を示すことにする。
【図面の簡単な説明】
【００１２】
【図１】本発明の一実施形態による顕微鏡画像を強調する方法のフロー図である。
【図２】サンプルの要素に入射する光の減衰を示す図である。
【図３】微小体積によって放射される光の散乱を示す図である。
【図４】共焦点顕微鏡法のジオメトリを示す図である。
【図５】顕微鏡の焦点レンズとサンプルとの間の複数のｚスタックを生成するプロセスを示す図である。
【図６】単一平面照明顕微鏡法の側方散乱ジオメトリを示す図である。
【図７】図７（ａ）〜図７（ｄ）は、図１の方法を用いて得られた一組の強調結果を示す図であり、入力画像は、フルオレセイン及び液体ゲルを用いて調製されたサンプルの画像であり、共焦点顕微鏡法を用いて取得されている。
【図８】図８（ａ）〜図８（ｃ）は、図１の方法を用いて得られた第１の一組の強調結果を示す図であり、入力画像は、ニューロステムセル（neuro stem cell）の画像であり、共焦点顕微鏡法を用いて取得されている。
【図９】図９（ａ）〜図９（ｃ）は、１の方法を用いて得られた第２の一組の強調結果を示す図であり、入力画像は、ニューロステムセルの画像であり、共焦点顕微鏡法を用いて取得されている。
【図１０】図１の方法を用いて得られた一組の強調結果を示す図であり、入力画像は、合成劣化画像（synthetically degraded image）であり、単一平面照明顕微鏡法を用いて取得されている。
【図１１】図１１（ａ）および図１１（ｂ）は、図１の方法を用いて得られた強調結果に対して１／α’ｎ０の値を変動させた効果を示す図であり、入力画像は、合成劣化画像であり、単一平面照明顕微鏡法を用いて取得されている。
【発明を実施するための形態】
【００１３】
図１を参照すると、本発明の一実施形態であって、顕微鏡画像を強調する方法である方法１００のステップが示されている。
【００１４】
方法１００への入力は、サンプルを照明しかつサンプルの要素により吸収されて散乱された光を収集する顕微鏡を用いて取得されたサンプルの画像である。画像のピクセルはサンプルの要素に対応する。サンプル内の各点における強度は、（サンプルを通過するにつれて次第に減衰される）入射光の成分と散乱に起因した散乱成分との和である。サンプルの所与の要素による入射光の散乱は、散乱パラメータの値によって記述される。散乱パラメータは、通常、放射係数ρｅｍであるか、そうでない場合は吸収係数ρａｂに等しい。ステップ１０２において、画像ピクセルごとに、対応する要素の散乱パラメータの値が、散乱パラメータの値と画像の強度とを関連させる数式を用いて計算される。ステップ１０４において、強調画像が、散乱パラメータの計算された値を用いて形成される。ステップ１０２の前に入力画像を線形に正規化することができる。さらに、方法１００は、強調画像内のピクセルの強度を０〜（２ｎ−１）の範囲に線形にスケーリングするステップをさらに含みうる。ここで、ｎは、ピクセルを表すのに用いられるビット数である。
【００１５】
数式の微分を以下で検討する。以下の検討は、共焦点顕微鏡法の場合と側方散乱ジオメトリの場合との式を解く方法に関する議論をさらに含んでいる。
【００１６】
１．場の理論的定式化
サンプル、場合によっては減衰媒体、光源、及び検出器（例えばカメラ）を含む撮像システム全体を含んだ関心領域Ω∈Ｒ３を仮定する。実際には、光源は無限遠から光を発することもできるが、一例では、光源は、関心領域∂Ωの境界から光を発するものとみなされる。しかしながら、これは本発明の実施形態の要件ではないことに留意されたい。ｒｓ∈Ωｓは、光源内の点の集合であり、ｒｐは検出器内のボクセルのロケーションである。
【００１７】
１．１ 光子密度及び光強度
光子（光）密度及び光強度場の数学モデルは次のように記述される。図２は、サンプルの要素に入射する光の減衰を示している。図２において、ｄｌ及びｄＡは、それぞれ要素の微小長さ及び微小面積を示す。式（１）において、ｃは媒体内での光速度であり、ｆ（ｒ）ｄｌｄＡは、微小体積ｄＶ＝ｄｌｄＡ（図２参照）内の光子の総数である。式（１）は、単位体積あたりの光子数ｆ（ｒ）と光強度ｎ（ｒ）との間の関係を与える。ｎ（ｒ）は、単位時間あたりに単位面積を通過する光子の数である。
【数１】
【００１８】
１．２ 減衰係数及び吸収係数
媒体を通る光の減衰度は、媒体の不透明度及び媒体内を進む距離に依存する。図２を参照して、光は、（図２にｎ（ｒ）とラベル付けされた矢印によって示されるように）点ｒにおいてｘ軸に沿って要素に入射すると仮定すると、微小厚さｄｌを有する媒体を通る強度の微分変化は、式（２）によって与えられる。式（２）において、ρａｂ（ｒ）は点ｒにおける光の吸収係数である。いくつかの論文では、ρａｂ（ｒ）は、吸光係数（extinction coefficient）としても知られている［５］、［８］、［１４］。ρａｂ（ｒ）は、一般的には光の波長の関数、すなわち、
【数２】
であるが、簡単にするために、下付き文字は以下の検討では省略されていることに留意されたい。複数の波長を取り扱うために、これらの式の一般化は簡単に行うことができる。
【数３】
【００１９】
ｒｓの光源から点ｒへの全減衰効果を計算するために、式（２）は、式（３）に示すようにｒｓからｒまで積分される。ここで、γ（ｒｓ：ｒ）は、ｒｓとｒとを結ぶ光線を示す。
【数４】
【００２０】
次に、ｒｓからｒへの全減衰効果を記述する式（４）は、全ての光源から点ｒへの光線にわたって式（３）のｎ（ｒ）を総和することによって導出される。式（４）において、下付き文字Ａを有するｎＡ（ｒ）は、減衰成分の影響を受けた光強度を示し、全ての光源から発して光源と点ｒ（ｒｓ∈Ωｓ）との間の光線に沿って減衰された全光強度である。式（４）は、光強度が一般に指数関数的に減衰し、指数関数的減衰率は、点が異なれば変動することを表している。
【数５】
【００２１】
１．３ 光子吸収率及び光子放射率
光子吸収率及び光子放射率は、式（５）に示すような連続の式（continuity equation）を用いて関係付けられる［１８］。
【数６】
ここで、
【数７】
は入射光の方向である。
【００２２】
式（１）及び式（５）を結合すると、単位体積あたり単位時間あたりの吸収光子の数が式（６）に与えられる。
【数８】
【００２３】
連続の式（すなわち式（５））を式（２）に関係付けると、単位体積あたり単位時間あたりの光子数を記述する式（７）が導出される。
【数９】
【００２４】
１．４ 散乱係数及び放射係数
媒体は全ての方向に光を散乱するので、散乱光は、媒体の他の部分の粒子によって再び吸収及び散乱される可能性がある。吸収される光子の数が、放射される光子の数に等しい場合、単位体積あたり単位時間あたりに放射される光子の数は、ｄｎ（ｒ）／ｄｌであり、式（７）によって、ｒ’における微小体積ｄｒ’によって放射される光の率は、ｎ（ｒ’）ρａｂ（ｒ’）ｄｒ’によって与えられる。しかしながら、一部の光エネルギーは、熱を通じて又は他の或る手段によって散逸される場合がある。散逸性媒体の場合、一部の光エネルギーは散逸される。この明細書では、ρｅｍは、放射係数を表すのに用いられ、放射率は、ｎ（ｒ’）ρｅｍ（ｒ’）ｄｒ’によって与えられ、ここで、ｎ（ｒ’）ρｅｍ（ｒ’）ｄｒ’≦ｎ（ｒ’）ρａｂ（ｒ’）ｄｒ’である。
【００２５】
図３は、微小体積によって放射された光の散乱を示している。図３に示すように、微小体積ｄｒ’によって放射された光子の総数は、ｎ（ｒ’）ρｅｍ（ｒ’）によって与えられる。図３を参照すると、これらの光子の一部分が点ｒに到達し、微小体積ｄｒ’からの光の散乱に起因して点ｒで受け取られる微小入射光強度は、式（８）に与えられる。
【数１０】
【００２６】
式（８）において、下付き文字Ｓは散乱成分を表し、γ（ｒ’：ｒ）はｒ’からｒへの光線である。第１項の分母は、３Ｄ空間のジオメトリを反映するジオメトリックファクターである。分子は、体積要素ｄｒ’によって単位時間あたりに放射される光子の数である。第２項は、ｒ’からｒへの光の減衰を表す。全てのｒ’∈Ω，ｒ’≠ｒにわたって積分すると、点ｒで受け取られる全散乱光は、式（９）に与えられる。
【数１１】
【００２７】
１．５ 画像形成
点ｒ∈Ωにおける全光強度は、式（１０）に示すように減衰成分及び散乱成分の和として記述することができる。
【数１２】
【００２８】
式（１０）に記述されるような物理モデルは、観察された画像に関係付けることができる。微小体積ｄｒによって単位時間あたりに放射される全光量は、ρｅｍ（ｒ）ｎ（ｒ）ｄｒである。検出器が（例えば、共焦点顕微鏡法により）この光の一部を検出して３次元観察画像ｕ０のピクセルｒｐを形成すると仮定すると、ｒｐにおけるピクセル強度は、式（１１）によってｕ０（ｒｐ）として与えられる。
【数１３】
【００２９】
式（１１）の積分は、全ての点ｒ∈Ωから点ｒｐへの全ての光線にわたって行われる。光は、媒体ロケーションｒからピクセルロケーションｒｐに進むときに減衰されるので、減衰項はこの式にも再び現れる。αγ＝αγ（ｒ，ｒｐ）は、検出器のレンズ系に依存する関数である。下付き文字γは、αγが光の経路に依存することを示すのに用いられる。
【００３０】
撮像の目的は、どのような物体が関心領域Ωに存在するのかを見つけ出すことである。換言すれば、Ω内の物質の光学特性を見つけ出すことが必要である。これらの特性はρａｂ（ｒ）及びρｅｍ（ｒ）によって与えられる。観察された画像ｕ０（ｒｐ）が与えられると、ρａｂ（ｒ）及びρｅｍ（ｒ）は、これらのパラメータについて式（１１）を解くことによって推定される。
【００３１】
以下の観察結果が上記式から作成される。
１．ジオメトリ：全てのジオメトリ情報は、点ｒから点ｒ’への光線を表す経路γ（ｒ’，ｒ）に組み込まれる。
２．光源：光源情報は、式（４）のΩＳ及びγ（ｒＳ：ｒ）にわたる総和によって与えられる。
３．大気散乱光：大気散乱光効果［１０］は、屋外の撮像の分野で知られている。屋外の撮像では、大気中の水粒子が太陽光を観察者に向けて反射し、これによって、光源として動作する。類似の効果は、本顕微鏡法の分野でも生じ、散乱成分に含まれる（式（１０）参照）。
４．非一意解：式（１１）の解は、一般に非一意である。例えば、Ωが不透明な箱を含み、この箱の画像が撮られる場合を考える。箱は不透明であるので、箱内のρａｂ及びρｅｍの値は確定されない。
【００３２】
１．６ 離散化
行列方程式が、各点ｒにおける全光強度（式（１０））を先ず離散化してＮ個の有限要素を形成することにより導出される。Ｎ個の有限要素はｒｉで示される。ここで、ｉ＝１，...，Ｎは、有限要素をラベル付けする整数インデックスである。有限要素は、以下では「ボクセル」と呼ばれ、離散化は、画像データ内の各それぞれのボクセルｒｐがボクセルｒｉのうちの１つに対応するように実行される。式（１２）において、ｒｋ∈γ（ｒｊ：ｒｉ）にわたる総和は、光線γ（ｒｊ：ｒｉ）が通過する全ての有限要素にわたる総和である。ρａｂ（ｒｋ）は、ボクセルｋにおける吸収係数であり、Δｒｋは、ボクセルｋ内に存在するγ（ｒｊ：ｒｉ）の線分の長さである。ΔＶはボクセルの体積である。式（１０）は、この場合、次のように書き換えることができる。
【数１４】
【００３３】
各ｉ＝１，...，Ｎについて、ｂｉ（ｒｉ）をｂｉ（ｒｉ）＝ｎＡ（ｒｉ）によって定義し、ｕｉ（ｒｉ）（又は略してｕｉ）をｕｉ／ρｅｍ（ｒｉ）＝ｎ（ｒｉ）によって定義する。その結果、式（１２）は次のように書き換えることができる。
【数１５】
【００３４】
Ｇ（ｒｉ，ｒｊ）を式（１４）によって定義する。
【数１６】
【００３５】
式（１３）は次のように書き換えることができる。
【数１７】
【００３６】
Ｇｉｊ＝Ｇ（ｒｉ，ｒｊ）並びにｕ＝（ｕ１，…，ｕＮ）及びｂ＝（ｂ１，…，ｂＮ）と定義すると、式（１５）は、式（１６）に示すように行列形式に書き換えることができる。
【数１８】
【００３７】
式（１６）は、以下に示すような式（１６ａ）を用いて数値的に解くことができ、それによって、ρ＝ρａｂ（ｒ）であり、これは、次にｑ−１ρｅｍ（ｒ）又はρｅｍ（ｒ）に等しくすることができる。ρ＞０であるので、カルーシュ・キューン・タッカー（Karush-Kuhn-Tucker）条件を回避するために、絶対符号が式（１６ａ）に必要とされる。
【数１９】
【００３８】
２．共焦点顕微鏡法
図４は、共焦点顕微鏡法の構成のジオメトリを示している。入射光は焦点レンズを通過し、点ｒｆに合焦される。式（４）の全ての光線にわたる総和は、焦点レンズからの全ての光線γ（ｒｓ：ｒｆ）にわたり総和する。レンズの面積は、入射光の発する点の集合、換言すれば入射光源Ωｓの点の集合となるように取ることができる。検出された光は、焦点レンズを通る同じ経路を介して進む。したがって、式（１１）における全ての光線にわたる総和は、同じ経路γ（ｒｓ：ｒｆ）上であるが反対方向に実行される。一例では、放射係数ρｅｍは、蛍光体の量子収量によって吸収係数ρａｂに関係付けられ［１９］、すなわち量子収量ｑによって吸収係数ρａｂに関係付けられる。ρａｂ（ｒ）をより簡単にρ（ｒ）と記述すると、ρ（ｒ）＝ｑ−１ρｅｍ（ｒ）が得られる。例えば、フルオレセインが蛍光体として用いられるとき、ｑは０．９２の値をとり、Ｈｏｅｓｃｈｔ ３３３４２が蛍光体として用いられるとき、ｑは０．８３の値をとる。代替的に、蛍光顕微鏡法では、蛍光体が光子を吸収し、ほぼ即時に光子を再放射することが可能な場合がある。したがって、別の例では、吸収された光子の総数は、放射された光子の総数に等しい（すなわちｑ＝１）と仮定され、ρ（ｒ）は、ρ（ｒ）＝ρｅｍ（ｒ）＝ρａｂ（ｒ）と設定される。
【００３９】
図５は、サンプルが離散ロケーションでどのように走査されて、画像取得プロセスでｚスタック（グレイの陰影）を生成するのかを示している。一例では、焦点レンズは、最初に、第１のスライスｚ＝０から特定の焦点距離離れた箇所に配置される。この第１のスライスｚ＝０の画像を捕捉した後、焦点レンズは、第２のスライスｚ＝１から同じ焦点距離離れた第２の位置にシフトされ、第２のスライスｚ＝１の画像が捕捉される。このプロセスは、全てのスライスの画像が捕捉されるまで繰り返される。
【００４０】
式（１７）に示すように、ｚスタックごとの光円錐の円盤面積にわたる平均ρ（ｒ）（すなわち〈ρ〉（ｚ））を計算することによる近似が用いられて、式（４）及び式（１１）が簡単化される。
【数２０】
【００４１】
式（１７）を用い、焦点レンズの入射光強度は一定の入射光強度ｎ０であると仮定すると、式（４）は式（１８）として記述することができる。
【数２１】
ここで、
【数２２】
である。βは光路の複雑な関数であるが、焦点レンズの焦点距離が変化しない限り定数である。
【００４２】
共焦点顕微鏡法の場合、図４に示すように、ｒｆで放射された光のみが光電子増倍管によって収集される。したがって、式（１１）のρｅｍ及びρａｂを、式（１７）で定義した〈ρ〉を用いて置き換えることができ、式（１９）が得られる。式（１９）において、ΔＶｆは、共焦点体積であり、
【数２３】
は定数であり、ｕ（ｒ）＝ρ（ｒ）ｎ（ｒ）である。以下に示すような式（１９）は、ρ（ｒ）をρ（ｒ）＝ρａｂ（ｒ）＝ｑ−１ρｅｍ（ｒ）と設定することによって導出される。ρａｂ（ｒ）＝ρｅｍ（ｒ）であると仮定する場合、項ｑは式（１９）から省略される。
【数２４】
【００４３】
一例では、散乱項が無視できる（すなわちｎｓ≪ｎＡ，ｎ（ｒ）≒ｎＡ（ｒ）、換言すれば、ｉ≠ｊの場合にＧｉｊ＝０及びＧｉｉ＝１／ρ（ｒｆ））と仮定することによって、式（１６）の解析解が得られる。言い換えると、この仮定は、サンプルの各要素における光強度が、他の要素からの散乱に起因した無視できる成分のみを含むということである。式（１８）及び式（１９）を式（１６）に代入すると、式（２０）が得られる。
【数２５】
ここで、ｕ０ｉ＝ｕ０（ｒｉ）であり、式（２０）のρＡは、散乱項が無視される場合に、真の光放射係数となる。画像は、方法１００において、画像ピクセルｒｉごとに計算された放射係数ρＡ（ｒｉ）を用いて強調される。
【００４４】
一例では、ρＡは、第１のスライスから開始して、ｚスタックを通じてスライスごとに観察画像から計算される。
【００４５】
１．第１のスライスｚ＝０の場合、式（２０）の積分は、ゼロの値、すなわちρＡ（ｒｉ，ｚ＝０）＝ｕ０ｉ／α’βｎ０を与える。これは、ρＡが、観察画像内の強度に比例することを暗に意味している。一例では、α’βは、照明を最も均一にするように較正することができるチューニングパラメータである。
【００４６】
２．第２のスライスのρＡは、式（２１）に従って第１のスライスのρＡに依存する。ここで、Δｚは離散化されたｚスタックの厚さであり、〈ρＡ（ｚ＝０）〉は、第１のスライスのρＡの値を用いて計算された平均値である。
【数２６】
【００４７】
３．第ｋのスライスのρＡは、式（２２）によって与えられる。ρＡの値は、第１のスライスから開始してスライスごとに計算されるので、第ｋのスライスのρＡを計算する時点では、ρＡの値は、第１のスライスから第（ｋ−１）のスライスまでの全てのスライスについてすでに得られている。したがって、項
【数２７】
の値を容易に得ることができる。強調画像全体を得るために、第１のスライスから最後のスライスまでのρＡが順に計算される。
【数２８】
【００４８】
代替的に、式（１６）を解くときに散乱項を含めることができる。散乱項を含めた結果、非解析解がもたらされ、これは、上記に示したような式（１６ａ）によって数値的に解くことができる。式（１６ａ）は、式（１８））によるｂｉ＝ｂ（ｒｉ）＝ｎＡ（ｒｉ）及び式（１９）による
【数２９】
と共に用いることができる。ベクトルｂ、ｕ、及び行列Ｇは、画像のボクセルごとに形成され、行列Ｇは式（１４）に従って形成される。ｉ≠ｊの場合、Ｇｉｊの計算は、点ｒｉと点ｒｊとの間の光線（すなわち直線）に沿ってρａｂ（ｒｋ）Δｒｋを総和することを伴う。計算時間を削減するために、サンプリングは、各ｚスタックの光円錐の円盤面積にわたる平均ρ（ｒ）を計算するときに実行することができる。
【００４９】
代替的に、∂Ｊ／∂ρｋ，∀ｋの値は数値的に求めることができるので、式（１６ａ）は勾配降下法を用いて数値的に解くことができる。一例では、ρＡ（式（２０））は、勾配降下法においてρの初期推定値に用いられる。本発明の実施形態を用いて実行される数値シミュレーションを通じて、ρＡ（式（２０））がρの良好な近似値であることが分かる。ρＡ（式（２０））をρの初期推定値として用いると、勾配降下法におけるローカルミニマム問題が低減される。
【００５０】
３．側方散乱ジオメトリ
側方散乱ジオメトリは、実際には、単一平面照明顕微鏡法（ＳＰＩＭ）のジオメトリである［２２］、［２３］、［２４］、［２５］、［２６］。図６は、光源が側方から発してサンプル、すなわち図６の「サンプル」の１つの平面を照明する側方散乱のジオメトリ配置を示している。図６に示すように、散乱光は、ＣＣＤカメラによって直交方向で収集される。このジオメトリ配置では、入射光線は一定でかつ平行である。式（４）は、点ｒ＝（ｘ，ｙ）における一定の入射強度をｎ０と表記することによって以下の式（２４）に簡約することができる。式（２４）では、積分は、図６に示すように水平のｘ方向にわたる。共焦点顕微鏡法の場合のように、ρは、ρ＝ｑ−１ρｅｍ＝ρａｂとして設定することができる。
【数３０】
【００５１】
散乱光は、減衰することなくＣＣＤカメラに直接進むものと仮定される。ほとんどのカメラ設定と同様に、（ＣＣＤカメラにおける）ピクセル点ｒｐとサンプルロケーションｒとの間には１対１の対応関係がある。したがって、式（１１）は、（ρ＝ｑ−１ρｅｍ＝ρａｂが用いられるとき）以下に示すような式（２５）として記述することができ、ここで、α’は、全ての光線にわたる総和等を含む量子収量及び検出器の積分された効果を表す。
【数３１】
【００５２】
式（２５）を用いると、式（１６）の行列方程式は、式（２６）として記述することができる。
【数３２】
【００５３】
一例では、解析解は、散乱項が無視できると仮定することによって、式（２７）に示すように得られる。式（２７）では、下付き文字Ａは、近似解が減衰項のみを用いて得られることを示すのに用いられる。この近似を用いると、式（２７）を数値的により容易に解くことができる。式（２７）の総和は、光源ｒｓからそれぞれの点ｒｉへのレーザービームの方向の光線（すなわち直線）に沿って実行される。
【数３３】
【００５４】
本発明の実施形態の数値実験は、式（２６）を用いた強調画像及び式（１６ａ）を用いた強調画像が約１％しか異ならないことを示し、ｎｓ≪ｎＡの近似が有効であることを示していた。
【００５５】
４．結果
数値計算が実行され、その結果はグラウンド・トルース（ground truth）と比較される。他の物理学に基づく復元方法［５］、［６］、［７］、［８］、［９］、［１０］、［１１］、［１２］、［１３］と比較することは可能ではない。その理由は、これらの方法を顕微鏡画像に適用することはできないからである。第１に、他の物理学に基づく方法は、３次元画像を強調するように設計されていない。第２に、これらの方法は、一定減衰媒体を仮定しているが、この仮定は顕微鏡画像では大きく破られる。
【００５６】
４．１ 有効性確認及び較正
方法１００は、グラウンド・トルースが実験的設計によって判明している特別に調製されたサンプルに関して有効性が確認される。この場合、画像強調が方法１００を用いて実行され、得られた結果がグラウンド・トルースと比較される。実験では、サンプルが、フルオレセイン及び液体ゲルを、ゲルが硬化するまでオービタルシェーカー上で混合することによって作製される。このようにして、サンプルは、３Ｄの体積全体にわたって均一になる。しかしながら、取得された画像の強度プロファイルは、減衰のために均一とはならない。その代わり、強度プロファイルは深度と共に減少する。図７（ａ）に示すように、方法１００を用いて強調された画像のそれぞれは、原入力画像よりも均一な強度プロファイル（最大強度投影）を有する。強調画像は、図７では「復元（restored）」と表記されている。パラメータα’βｎ０は、顕微鏡に対して較正される。図７（ａ）に示すように、α’βｎ０＝１８１．２７が最良の結果を与える。他方、それよりも低い値１２１．５１及び９０．０２は過補償を招く。パラメータα’βｎ０の値をＣと表記すると、２つのパラメータ（Ｃ１，Ｃ２）と異なるレーザー強度（ｎ１，ｎ２）との間の関係は、単純にＣ１／Ｃ２＝ｎ１／ｎ２である。したがって、図７（ａ）の較正されたパラメータ値は、図７（ｃ）に示すように、異なるレーザー強度、例えば１．５ｎ０及び２．０ｎ０で撮った画像にも用いることができる。ここで、ｎ０は、図７（ａ）で用いられたレーザー強度である。図７（ｂ）には、原画像７０２及び（レーザー強度が１．５ｎ０であるときの）強調画像７０４の２Ｄ投影が示されているのに対して、図７（ｄ）には、原画像７０６及び（レーザー強度が２．０ｎ０であるときの）強調画像７０８の２Ｄ投影が示されている。強調画像はより均一に照明されることが分かる。
【００５７】
４．２ 共焦点顕微鏡法
方法１００の有効性を実証するために、核がＨｏｅｓｃｈｔ ３３３４２で染色された、マウスの胚からのニューロステムセルの３Ｄ画像が強調される。オリンパス社製の点走査共焦点（Point Scanning Confocal）ＦＶ１０００システムを用いて画像を取得した。開口数１．２を有する６０×の水レンズで撮像を行った。ダイオードレーザー４０ｎｍを用いて、Ｈｏｅｓｃｈｔで染色されたニューロスフェア（neurosphere）を励起させた。サンプリング速度は、２μｍ／ピクセルに設定した。原顕微鏡画像は、サイズがｘ方向及びｙ方向に０．１３７μｍの解像度並びにｚ方向に０．２μｍの解像度を有する５１２×５１２×ｎｚ個のボクセルである。ここで、ｎｚは、画像内のｚスタックの数である。
【００５８】
計算時間を低減するために、原画像は、ｚ方向の解像度を維持しながら、ｘ方向及びｙ方向でボクセルを平均することによって２５６×２５６×ｎｚ個のボクセルにダウンサンプリングされる。図８及び図９は、式（２０）を用いて強調された２５６×２５６×ｎｚ個のボクセルの画像の強調結果を示している。図８及び図９に示す共焦点顕微鏡画像は、それぞれ１５５個のｚスタック及び１６３個のｚスタック（すなわちｎｚ＝１５５及びｎｚ＝１６３）を有する。より具体的には、図８（ａ）及び図９（ａ）は、それぞれの原画像のｙｚ平面上への最大強度投影を示している（図８及図９に原画像と表記されている）のに対して、図８（ｂ）及び図９（ｂ）は、それぞれの強調画像のｙｚ平面上への最大強度投影を示している（図８及び図９に復元図と表記されている）。図８（ｂ）及び図９（ｂ）において、調整されたチューニングパラメータは、１／α’βｎ０＝０．０１４９９５として設定される。これは、最適な強調結果を与えることが分かる。図８（ｃ）及び図９（ｃ）は、原画像（実線）及び強調画像（破線）の双方のｚ軸上への（ｘｙ平面における最も明るい０．１％のボクセルにわたって平均された）最大強度投影を示している。照明レーザーは底部から発するので、原画像（図８（ａ）又は図９（ａ））について、ボクセルが画像の底部ではるかに明るく、画像の強度が画像の上部に向かって低下する（換言すれば、照明は均一でない）ことが図８及び図９から容易に観察することができる。しかしながら、図８（ｂ）及び図９（ｂ）に示すように、照明は、強調後、均一になる。さらに、図８（ｃ）及び図９（ｃ）は、原画像（実線）の強度プロファイルと強調画像（破線）の強度プロファイルとの間の相違を明らかに示している。加えて、図８及び図９に示すように、底部のｚスタックの露出過度のエリアも、強調方法１００によって正しく補償される。したがって、強調画像の均一な照明を達成することができるので、方法１００を用いて画像を復元することは有利であることが図８及び図９から分かる。その後、この強調画像に対してヒストグラム等化等の他の画像強調方法を用いることができる。強調画像は、平均してより暗くもなるが、多くの画像処理技法は、平均ボクセル強度に対してロバストである。
【００５９】
４．３ 側方散乱顕微鏡法
合成劣化画像に対して側方散乱ジオメトリの計算を行った。不均一な照明（すなわち、光源が左から発すると仮定して最大強度投影が指数関数的に低下する）による合成画像を均一な照明の画像から生成した。図１０は、２５６×２５６ピクセル画像の強調結果を示している（強調画像は、図１０では「復元（Restored）」とラベル付けされている）。図１０において、式（２７）は、画像を強調するのに用いられる。ｎ０は、入射光強度であり、α’は、通例未知であるジオメトリックファクターである。チューニングパラメータ１／α’ｎ０は、最適な結果を得るように調整することができる。図１１（ａ）は、小さな１／α’ｎ０（０．００１１）が用いられたときの強調画像を示しているのに対して、図１１（ｂ）は、大きな１／α’ｎ０（０．０１１１）が用いられたときの強調画像を示している。図１１（ａ）及び図１１（ｂ）から分かるように、小さな１／α’ｎ０が用いられたとき、画像はほとんど強調されず、大きな１／α’ｎ０が用いられたとき、減衰効果が過補償される。１／α’ｎ０の最適値は、図１０の強調画像を得るのに用いられた０．００９５である。図１０に示すように、１／α’ｎ０が０．００９５に設定されたとき、強調画像はほぼ完全に（均一に）照明される。
【００６０】
上述したように、方法１００は、均一な照明を有する強調画像を得ることができるので有利である。換言すれば、方法１００を用いて画像を強調することによって、光学顕微鏡法の基本的な光減衰及び散乱の問題を軽減することができる。
【００６１】
方法１００で用いられる式ｂ＝Ｇ・ｕの微分は、強い理論的基礎に基づいて定式化されており、連続の式によって表わされる保存則等の物理学の基本法則に基づいている。さらに、場の理論的手法（field theoretical approach）がこの微分に用いられている。
【００６２】
方法１００は、物理学に基づく復元方法のタイプであり、物理学に基づく復元方法は、モデルに基づくコントラスト強調方法（例えば、ヒストグラム等化）を上回る多くの利点を有する。モデルに基づく方法［１５］、［１６］、［１７］は、一般に、画像特性が画像全体にわたって一定であると仮定し、この仮定は、天候で劣化した画像では無視される。さらに、物理モデルは物理学の法則に基づいて構築され、物理学は、否定できない真実である可能性が最も高い。物理学に基づく復元技法は、多くの用途で用いることができる。そのような復元技法の１つの態様は、その有効性が数桁の物理的な長さスケールにわたることである。航空機による監視では、物理的な長さスケールはほぼ≒１０ｋｍ程度であり、水中監視では、物理的な長さスケールはほぼ≒１０ｍ程度である。
【００６３】
物理学に基づく復元方法が天候により劣化した画像の復元に用いられてきたが、それらの方法は、（≒１００μｍの物理的な長さスケールを有する）光学顕微鏡法の画像強調の分野で適切に探求されていなかった。方法１００は、顕微鏡画像に対して用いられる物理学に基づく復元技法のタイプであり、物理学に基づく復元の長さスケールを８桁に拡大する。
【００６４】
方法１００は、物理学に基づく復元技法のタイプであっても、全ての既存の物理学に基づく復元技法とは根本的に異なっている。既存の物理学に基づく復元技法では、減衰媒体の一定の吸収係数が前提とされるのに対して、このことは、方法１００では前提とされない。さらに、方法１００では、サンプルと減衰媒体とが区別されない。一般的な一組の式が導出され、画像取得におけるあらゆるジオメトリ設定を取り扱うために方法１００で用いられる。方法１００を用いるには、光源及びカメラ等の検出機器の詳細を指定することだけが必要である。他方、既存の物理学に基づく方法［５］、［６］、［７］、［８］、［９］、［１０］、［１１］、［１２］、［１３］は、以下の理由によって３次元顕微鏡画像に適用することさえできない。第１に、既存の物理学に基づく方法は、減衰媒体を「除去」して、２次元シーンをリトリーブする。これとは対照的に、方法１００では、減衰媒体も画像情報を含む。これは、媒体を除去するのではなく媒体の真の信号を復元することが必要であるので有利である。第２に、既存の方法は均一な減衰媒体を仮定しているが、この仮定は、顕微鏡画像では大きく破られる。これとは対照的に、そのような仮定は方法１００ではなされない。
【００６５】
［参考文献］
【表１】
【表２】
【表３】
【技術分野】
【０００１】
本発明は、顕微鏡を用いて取得された画像である顕微鏡画像を強調する方法及びシステムに関する。特に、強調画像は、不均一な光の減衰及び散乱による劣化の影響がより少ない画像である。
【背景技術】
【０００２】
顕微鏡法［１］は、生物学の重要な光学撮像技法である。２光子励起顕微鏡法及び単一平面照明顕微鏡法等の多くの顕微鏡技法があるが、共焦点顕微鏡法［１］は、バイオイメージングにとって最も重要なツールの１つとなってきている。共焦点顕微鏡法では、焦点ずれの光が、ピンホールの使用により取り除かれる。入射照明光は、ピンホールを通過し、サンプル内の小さな領域上に合焦され、この領域で、入射照明光はサンプルによって散乱される。入射照明光と同じ経路に沿って進む散乱光のみが、戻ってピンホールを通過し、そのような光は、光電子増倍管等の光検出器において再び合焦され、これによって、画像が生成される。共焦点顕微鏡を通じて取得された画像は、従来の広視野顕微鏡によって作成される画像よりも鮮明である。しかしながら、光減衰効果による劣化は、共焦点顕微鏡法では深刻である。共焦点顕微鏡法における基本的な問題は光透過問題である。入射光は、散乱されるので減衰し、したがって、厚いサンプルを透過することができない。その結果、サンプル内の深い領域から取得された画像は、サンプルの表面に近い領域から取得された画像よりも指数関数的に暗く見える。光透過の難題は、共焦点顕微鏡法に限ったものではない。単一平面照明顕微鏡法及び広視野顕微鏡法等の他の光学顕微鏡技法も同じ問題に悩まされている。古典的な空間不変デコンボリューション手法［２］、［３］、［４］は、顕微鏡法撮像のこの問題にうまく対処することができない。
【０００３】
上述した問題を解決する試みは、レーザー出力を増加させるか又は光電子増倍管の感度を増大させるかのいずれかによって行われてきた［２０］。双方の技法は、不十分であり、レーザー出力の増加は光退色効果を加速させるのに対して、光電子増倍管の感度の増大は雑音を追加するという欠点を有する。Umesh Adiga及びB. B. Chaudhury［２１］は、画像スタックの深度に沿った光減衰を考慮に入れて画像を復元するために、単純なしきい値処理方法を用いて前景から背景を分離することを議論している。この技法は、軸方向解像度を改善するために、画像ボクセルが、画像スタックに画像スライスを仮想的に挿入するＸＯＲ輪郭形成（XOR contouring）及びモーフィングに基づいて等方性である（これは共焦点顕微鏡法には当てはまらない）と仮定する。
【０００４】
表面的には無関係な技術分野が屋外撮像の分野である。その分野内では、大気エアロゾルに起因した劣化画像の復元の問題が、監視、ナビゲーション、追跡、及びリモートセンシング［５］、［１０］等のその重要な用途のために広範囲にわたって研究されてきた［５］、［６］、［７］、［８］、［９］、［１０］、［１１］、［１２］、［１３］。これらの劣化画像の復元に用いられる復元技法と類似の復元技法は、水中ビジョン［８］、［９］において、具体的には水中ケーブル及び水中パイプラインの監視等について、新しい用途も創出している。さまざまな復元アルゴリズムが、均一な媒体を通る光減衰及び光散乱（大気散乱光）の物理モデルに基づいて提案されている。そのような画像復元アルゴリズムに関するこれまでの研究の１つ［５］は、シーンの奥行きの正確な情報を必要とする。その後の研究は、シーンの奥行きの必要性を巧みに回避したが、必要とされる情報を回復するのに複数の画像を必要とする［８］、［９］、［１０］、［１４］。Narasimhan及びNayar［１０］、［１１］、［１２］、［１３］は、１つの劣化画像のみから全ての必要な情報を抽出する対話型アルゴリズムを開発した。この方法は、大気散乱光の色（airlight color）及び「良好な色の領域」の手動選択を必要とする。これらの復元技法に関する基本的な問題は、雑音の増幅である。この基本的な問題を取り扱う試みは、Kaftory他［１４］によって提案された変分手法（variational approach）における正則化項を用いることで行われる。
【発明の概要】
【発明が解決しようとする課題】
【０００５】
本発明は、上記問題を克服することができる、画像を復元する方法を提供することを目的とする。
【課題を解決するための手段】
【０００６】
概括的に言えば、本発明は、或る照明方向でサンプルに光を当てることによりサンプルを照明し散乱光を捕捉することによって形成された顕微鏡画像を用いて強調画像を作成することを提案する。これは、顕微鏡画像の構成部分の強度を、サンプルの複数のそれぞれの要素における散乱パラメータのそれぞれの値に関連させる式を用いて行われる。散乱パラメータは、放射係数ρｅｍとすることもできるし、そうでない場合は吸収係数ρａｂに等しくすることもできる。この式を解いて、散乱パラメータの値を見つける。散乱パラメータの値を用いて、強調画像、例えば、散乱パラメータ自体の変動をマッピングする画像を構成する。
【０００７】
散乱パラメータの値が正確に見つけられることを条件として、強調画像は、不均一な光の減衰及び散乱に起因した劣化を原画像よりも受けにくくなる。
【０００８】
上記式は、各要素の散乱パラメータの値を、入射光の方向に沿った要素の散乱パラメータの値の関数として与えることができる。この場合、散乱パラメータの値は、ロケーションについて連続的に、さらに照明方向でも連続的に見つけることができる。
【０００９】
上記式は、照明方向に平行な一組の要素を取り囲む領域にわたって規定された散乱パラメータの平均値を用いることができる。
【００１０】
本発明は、代替的に、そのような方法を実行するコンピュータシステムとして表現することもできる。このコンピュータシステムは、画像を取得するデバイス、例えば顕微鏡と統合することができる。本発明は、上記方法のステップを実行する、コンピュータシステムによって動作可能なプログラム命令を含む、有形のコンピュータ媒体に記録されるようなコンピュータプログラム製品として表現することもできる。
【００１１】
次に、単なる例示の目的で、以下の図面を参照して本発明の一実施形態を示すことにする。
【図面の簡単な説明】
【００１２】
【図１】本発明の一実施形態による顕微鏡画像を強調する方法のフロー図である。
【図２】サンプルの要素に入射する光の減衰を示す図である。
【図３】微小体積によって放射される光の散乱を示す図である。
【図４】共焦点顕微鏡法のジオメトリを示す図である。
【図５】顕微鏡の焦点レンズとサンプルとの間の複数のｚスタックを生成するプロセスを示す図である。
【図６】単一平面照明顕微鏡法の側方散乱ジオメトリを示す図である。
【図７】図７（ａ）〜図７（ｄ）は、図１の方法を用いて得られた一組の強調結果を示す図であり、入力画像は、フルオレセイン及び液体ゲルを用いて調製されたサンプルの画像であり、共焦点顕微鏡法を用いて取得されている。
【図８】図８（ａ）〜図８（ｃ）は、図１の方法を用いて得られた第１の一組の強調結果を示す図であり、入力画像は、ニューロステムセル（neuro stem cell）の画像であり、共焦点顕微鏡法を用いて取得されている。
【図９】図９（ａ）〜図９（ｃ）は、１の方法を用いて得られた第２の一組の強調結果を示す図であり、入力画像は、ニューロステムセルの画像であり、共焦点顕微鏡法を用いて取得されている。
【図１０】図１の方法を用いて得られた一組の強調結果を示す図であり、入力画像は、合成劣化画像（synthetically degraded image）であり、単一平面照明顕微鏡法を用いて取得されている。
【図１１】図１１（ａ）および図１１（ｂ）は、図１の方法を用いて得られた強調結果に対して１／α’ｎ０の値を変動させた効果を示す図であり、入力画像は、合成劣化画像であり、単一平面照明顕微鏡法を用いて取得されている。
【発明を実施するための形態】
【００１３】
図１を参照すると、本発明の一実施形態であって、顕微鏡画像を強調する方法である方法１００のステップが示されている。
【００１４】
方法１００への入力は、サンプルを照明しかつサンプルの要素により吸収されて散乱された光を収集する顕微鏡を用いて取得されたサンプルの画像である。画像のピクセルはサンプルの要素に対応する。サンプル内の各点における強度は、（サンプルを通過するにつれて次第に減衰される）入射光の成分と散乱に起因した散乱成分との和である。サンプルの所与の要素による入射光の散乱は、散乱パラメータの値によって記述される。散乱パラメータは、通常、放射係数ρｅｍであるか、そうでない場合は吸収係数ρａｂに等しい。ステップ１０２において、画像ピクセルごとに、対応する要素の散乱パラメータの値が、散乱パラメータの値と画像の強度とを関連させる数式を用いて計算される。ステップ１０４において、強調画像が、散乱パラメータの計算された値を用いて形成される。ステップ１０２の前に入力画像を線形に正規化することができる。さらに、方法１００は、強調画像内のピクセルの強度を０〜（２ｎ−１）の範囲に線形にスケーリングするステップをさらに含みうる。ここで、ｎは、ピクセルを表すのに用いられるビット数である。
【００１５】
数式の微分を以下で検討する。以下の検討は、共焦点顕微鏡法の場合と側方散乱ジオメトリの場合との式を解く方法に関する議論をさらに含んでいる。
【００１６】
１．場の理論的定式化
サンプル、場合によっては減衰媒体、光源、及び検出器（例えばカメラ）を含む撮像システム全体を含んだ関心領域Ω∈Ｒ３を仮定する。実際には、光源は無限遠から光を発することもできるが、一例では、光源は、関心領域∂Ωの境界から光を発するものとみなされる。しかしながら、これは本発明の実施形態の要件ではないことに留意されたい。ｒｓ∈Ωｓは、光源内の点の集合であり、ｒｐは検出器内のボクセルのロケーションである。
【００１７】
１．１ 光子密度及び光強度
光子（光）密度及び光強度場の数学モデルは次のように記述される。図２は、サンプルの要素に入射する光の減衰を示している。図２において、ｄｌ及びｄＡは、それぞれ要素の微小長さ及び微小面積を示す。式（１）において、ｃは媒体内での光速度であり、ｆ（ｒ）ｄｌｄＡは、微小体積ｄＶ＝ｄｌｄＡ（図２参照）内の光子の総数である。式（１）は、単位体積あたりの光子数ｆ（ｒ）と光強度ｎ（ｒ）との間の関係を与える。ｎ（ｒ）は、単位時間あたりに単位面積を通過する光子の数である。
【数１】
【００１８】
１．２ 減衰係数及び吸収係数
媒体を通る光の減衰度は、媒体の不透明度及び媒体内を進む距離に依存する。図２を参照して、光は、（図２にｎ（ｒ）とラベル付けされた矢印によって示されるように）点ｒにおいてｘ軸に沿って要素に入射すると仮定すると、微小厚さｄｌを有する媒体を通る強度の微分変化は、式（２）によって与えられる。式（２）において、ρａｂ（ｒ）は点ｒにおける光の吸収係数である。いくつかの論文では、ρａｂ（ｒ）は、吸光係数（extinction coefficient）としても知られている［５］、［８］、［１４］。ρａｂ（ｒ）は、一般的には光の波長の関数、すなわち、
【数２】
であるが、簡単にするために、下付き文字は以下の検討では省略されていることに留意されたい。複数の波長を取り扱うために、これらの式の一般化は簡単に行うことができる。
【数３】
【００１９】
ｒｓの光源から点ｒへの全減衰効果を計算するために、式（２）は、式（３）に示すようにｒｓからｒまで積分される。ここで、γ（ｒｓ：ｒ）は、ｒｓとｒとを結ぶ光線を示す。
【数４】
【００２０】
次に、ｒｓからｒへの全減衰効果を記述する式（４）は、全ての光源から点ｒへの光線にわたって式（３）のｎ（ｒ）を総和することによって導出される。式（４）において、下付き文字Ａを有するｎＡ（ｒ）は、減衰成分の影響を受けた光強度を示し、全ての光源から発して光源と点ｒ（ｒｓ∈Ωｓ）との間の光線に沿って減衰された全光強度である。式（４）は、光強度が一般に指数関数的に減衰し、指数関数的減衰率は、点が異なれば変動することを表している。
【数５】
【００２１】
１．３ 光子吸収率及び光子放射率
光子吸収率及び光子放射率は、式（５）に示すような連続の式（continuity equation）を用いて関係付けられる［１８］。
【数６】
ここで、
【数７】
は入射光の方向である。
【００２２】
式（１）及び式（５）を結合すると、単位体積あたり単位時間あたりの吸収光子の数が式（６）に与えられる。
【数８】
【００２３】
連続の式（すなわち式（５））を式（２）に関係付けると、単位体積あたり単位時間あたりの光子数を記述する式（７）が導出される。
【数９】
【００２４】
１．４ 散乱係数及び放射係数
媒体は全ての方向に光を散乱するので、散乱光は、媒体の他の部分の粒子によって再び吸収及び散乱される可能性がある。吸収される光子の数が、放射される光子の数に等しい場合、単位体積あたり単位時間あたりに放射される光子の数は、ｄｎ（ｒ）／ｄｌであり、式（７）によって、ｒ’における微小体積ｄｒ’によって放射される光の率は、ｎ（ｒ’）ρａｂ（ｒ’）ｄｒ’によって与えられる。しかしながら、一部の光エネルギーは、熱を通じて又は他の或る手段によって散逸される場合がある。散逸性媒体の場合、一部の光エネルギーは散逸される。この明細書では、ρｅｍは、放射係数を表すのに用いられ、放射率は、ｎ（ｒ’）ρｅｍ（ｒ’）ｄｒ’によって与えられ、ここで、ｎ（ｒ’）ρｅｍ（ｒ’）ｄｒ’≦ｎ（ｒ’）ρａｂ（ｒ’）ｄｒ’である。
【００２５】
図３は、微小体積によって放射された光の散乱を示している。図３に示すように、微小体積ｄｒ’によって放射された光子の総数は、ｎ（ｒ’）ρｅｍ（ｒ’）によって与えられる。図３を参照すると、これらの光子の一部分が点ｒに到達し、微小体積ｄｒ’からの光の散乱に起因して点ｒで受け取られる微小入射光強度は、式（８）に与えられる。
【数１０】
【００２６】
式（８）において、下付き文字Ｓは散乱成分を表し、γ（ｒ’：ｒ）はｒ’からｒへの光線である。第１項の分母は、３Ｄ空間のジオメトリを反映するジオメトリックファクターである。分子は、体積要素ｄｒ’によって単位時間あたりに放射される光子の数である。第２項は、ｒ’からｒへの光の減衰を表す。全てのｒ’∈Ω，ｒ’≠ｒにわたって積分すると、点ｒで受け取られる全散乱光は、式（９）に与えられる。
【数１１】
【００２７】
１．５ 画像形成
点ｒ∈Ωにおける全光強度は、式（１０）に示すように減衰成分及び散乱成分の和として記述することができる。
【数１２】
【００２８】
式（１０）に記述されるような物理モデルは、観察された画像に関係付けることができる。微小体積ｄｒによって単位時間あたりに放射される全光量は、ρｅｍ（ｒ）ｎ（ｒ）ｄｒである。検出器が（例えば、共焦点顕微鏡法により）この光の一部を検出して３次元観察画像ｕ０のピクセルｒｐを形成すると仮定すると、ｒｐにおけるピクセル強度は、式（１１）によってｕ０（ｒｐ）として与えられる。
【数１３】
【００２９】
式（１１）の積分は、全ての点ｒ∈Ωから点ｒｐへの全ての光線にわたって行われる。光は、媒体ロケーションｒからピクセルロケーションｒｐに進むときに減衰されるので、減衰項はこの式にも再び現れる。αγ＝αγ（ｒ，ｒｐ）は、検出器のレンズ系に依存する関数である。下付き文字γは、αγが光の経路に依存することを示すのに用いられる。
【００３０】
撮像の目的は、どのような物体が関心領域Ωに存在するのかを見つけ出すことである。換言すれば、Ω内の物質の光学特性を見つけ出すことが必要である。これらの特性はρａｂ（ｒ）及びρｅｍ（ｒ）によって与えられる。観察された画像ｕ０（ｒｐ）が与えられると、ρａｂ（ｒ）及びρｅｍ（ｒ）は、これらのパラメータについて式（１１）を解くことによって推定される。
【００３１】
以下の観察結果が上記式から作成される。
１．ジオメトリ：全てのジオメトリ情報は、点ｒから点ｒ’への光線を表す経路γ（ｒ’，ｒ）に組み込まれる。
２．光源：光源情報は、式（４）のΩＳ及びγ（ｒＳ：ｒ）にわたる総和によって与えられる。
３．大気散乱光：大気散乱光効果［１０］は、屋外の撮像の分野で知られている。屋外の撮像では、大気中の水粒子が太陽光を観察者に向けて反射し、これによって、光源として動作する。類似の効果は、本顕微鏡法の分野でも生じ、散乱成分に含まれる（式（１０）参照）。
４．非一意解：式（１１）の解は、一般に非一意である。例えば、Ωが不透明な箱を含み、この箱の画像が撮られる場合を考える。箱は不透明であるので、箱内のρａｂ及びρｅｍの値は確定されない。
【００３２】
１．６ 離散化
行列方程式が、各点ｒにおける全光強度（式（１０））を先ず離散化してＮ個の有限要素を形成することにより導出される。Ｎ個の有限要素はｒｉで示される。ここで、ｉ＝１，...，Ｎは、有限要素をラベル付けする整数インデックスである。有限要素は、以下では「ボクセル」と呼ばれ、離散化は、画像データ内の各それぞれのボクセルｒｐがボクセルｒｉのうちの１つに対応するように実行される。式（１２）において、ｒｋ∈γ（ｒｊ：ｒｉ）にわたる総和は、光線γ（ｒｊ：ｒｉ）が通過する全ての有限要素にわたる総和である。ρａｂ（ｒｋ）は、ボクセルｋにおける吸収係数であり、Δｒｋは、ボクセルｋ内に存在するγ（ｒｊ：ｒｉ）の線分の長さである。ΔＶはボクセルの体積である。式（１０）は、この場合、次のように書き換えることができる。
【数１４】
【００３３】
各ｉ＝１，...，Ｎについて、ｂｉ（ｒｉ）をｂｉ（ｒｉ）＝ｎＡ（ｒｉ）によって定義し、ｕｉ（ｒｉ）（又は略してｕｉ）をｕｉ／ρｅｍ（ｒｉ）＝ｎ（ｒｉ）によって定義する。その結果、式（１２）は次のように書き換えることができる。
【数１５】
【００３４】
Ｇ（ｒｉ，ｒｊ）を式（１４）によって定義する。
【数１６】
【００３５】
式（１３）は次のように書き換えることができる。
【数１７】
【００３６】
Ｇｉｊ＝Ｇ（ｒｉ，ｒｊ）並びにｕ＝（ｕ１，…，ｕＮ）及びｂ＝（ｂ１，…，ｂＮ）と定義すると、式（１５）は、式（１６）に示すように行列形式に書き換えることができる。
【数１８】
【００３７】
式（１６）は、以下に示すような式（１６ａ）を用いて数値的に解くことができ、それによって、ρ＝ρａｂ（ｒ）であり、これは、次にｑ−１ρｅｍ（ｒ）又はρｅｍ（ｒ）に等しくすることができる。ρ＞０であるので、カルーシュ・キューン・タッカー（Karush-Kuhn-Tucker）条件を回避するために、絶対符号が式（１６ａ）に必要とされる。
【数１９】
【００３８】
２．共焦点顕微鏡法
図４は、共焦点顕微鏡法の構成のジオメトリを示している。入射光は焦点レンズを通過し、点ｒｆに合焦される。式（４）の全ての光線にわたる総和は、焦点レンズからの全ての光線γ（ｒｓ：ｒｆ）にわたり総和する。レンズの面積は、入射光の発する点の集合、換言すれば入射光源Ωｓの点の集合となるように取ることができる。検出された光は、焦点レンズを通る同じ経路を介して進む。したがって、式（１１）における全ての光線にわたる総和は、同じ経路γ（ｒｓ：ｒｆ）上であるが反対方向に実行される。一例では、放射係数ρｅｍは、蛍光体の量子収量によって吸収係数ρａｂに関係付けられ［１９］、すなわち量子収量ｑによって吸収係数ρａｂに関係付けられる。ρａｂ（ｒ）をより簡単にρ（ｒ）と記述すると、ρ（ｒ）＝ｑ−１ρｅｍ（ｒ）が得られる。例えば、フルオレセインが蛍光体として用いられるとき、ｑは０．９２の値をとり、Ｈｏｅｓｃｈｔ ３３３４２が蛍光体として用いられるとき、ｑは０．８３の値をとる。代替的に、蛍光顕微鏡法では、蛍光体が光子を吸収し、ほぼ即時に光子を再放射することが可能な場合がある。したがって、別の例では、吸収された光子の総数は、放射された光子の総数に等しい（すなわちｑ＝１）と仮定され、ρ（ｒ）は、ρ（ｒ）＝ρｅｍ（ｒ）＝ρａｂ（ｒ）と設定される。
【００３９】
図５は、サンプルが離散ロケーションでどのように走査されて、画像取得プロセスでｚスタック（グレイの陰影）を生成するのかを示している。一例では、焦点レンズは、最初に、第１のスライスｚ＝０から特定の焦点距離離れた箇所に配置される。この第１のスライスｚ＝０の画像を捕捉した後、焦点レンズは、第２のスライスｚ＝１から同じ焦点距離離れた第２の位置にシフトされ、第２のスライスｚ＝１の画像が捕捉される。このプロセスは、全てのスライスの画像が捕捉されるまで繰り返される。
【００４０】
式（１７）に示すように、ｚスタックごとの光円錐の円盤面積にわたる平均ρ（ｒ）（すなわち〈ρ〉（ｚ））を計算することによる近似が用いられて、式（４）及び式（１１）が簡単化される。
【数２０】
【００４１】
式（１７）を用い、焦点レンズの入射光強度は一定の入射光強度ｎ０であると仮定すると、式（４）は式（１８）として記述することができる。
【数２１】
ここで、
【数２２】
である。βは光路の複雑な関数であるが、焦点レンズの焦点距離が変化しない限り定数である。
【００４２】
共焦点顕微鏡法の場合、図４に示すように、ｒｆで放射された光のみが光電子増倍管によって収集される。したがって、式（１１）のρｅｍ及びρａｂを、式（１７）で定義した〈ρ〉を用いて置き換えることができ、式（１９）が得られる。式（１９）において、ΔＶｆは、共焦点体積であり、
【数２３】
は定数であり、ｕ（ｒ）＝ρ（ｒ）ｎ（ｒ）である。以下に示すような式（１９）は、ρ（ｒ）をρ（ｒ）＝ρａｂ（ｒ）＝ｑ−１ρｅｍ（ｒ）と設定することによって導出される。ρａｂ（ｒ）＝ρｅｍ（ｒ）であると仮定する場合、項ｑは式（１９）から省略される。
【数２４】
【００４３】
一例では、散乱項が無視できる（すなわちｎｓ≪ｎＡ，ｎ（ｒ）≒ｎＡ（ｒ）、換言すれば、ｉ≠ｊの場合にＧｉｊ＝０及びＧｉｉ＝１／ρ（ｒｆ））と仮定することによって、式（１６）の解析解が得られる。言い換えると、この仮定は、サンプルの各要素における光強度が、他の要素からの散乱に起因した無視できる成分のみを含むということである。式（１８）及び式（１９）を式（１６）に代入すると、式（２０）が得られる。
【数２５】
ここで、ｕ０ｉ＝ｕ０（ｒｉ）であり、式（２０）のρＡは、散乱項が無視される場合に、真の光放射係数となる。画像は、方法１００において、画像ピクセルｒｉごとに計算された放射係数ρＡ（ｒｉ）を用いて強調される。
【００４４】
一例では、ρＡは、第１のスライスから開始して、ｚスタックを通じてスライスごとに観察画像から計算される。
【００４５】
１．第１のスライスｚ＝０の場合、式（２０）の積分は、ゼロの値、すなわちρＡ（ｒｉ，ｚ＝０）＝ｕ０ｉ／α’βｎ０を与える。これは、ρＡが、観察画像内の強度に比例することを暗に意味している。一例では、α’βは、照明を最も均一にするように較正することができるチューニングパラメータである。
【００４６】
２．第２のスライスのρＡは、式（２１）に従って第１のスライスのρＡに依存する。ここで、Δｚは離散化されたｚスタックの厚さであり、〈ρＡ（ｚ＝０）〉は、第１のスライスのρＡの値を用いて計算された平均値である。
【数２６】
【００４７】
３．第ｋのスライスのρＡは、式（２２）によって与えられる。ρＡの値は、第１のスライスから開始してスライスごとに計算されるので、第ｋのスライスのρＡを計算する時点では、ρＡの値は、第１のスライスから第（ｋ−１）のスライスまでの全てのスライスについてすでに得られている。したがって、項
【数２７】
の値を容易に得ることができる。強調画像全体を得るために、第１のスライスから最後のスライスまでのρＡが順に計算される。
【数２８】
【００４８】
代替的に、式（１６）を解くときに散乱項を含めることができる。散乱項を含めた結果、非解析解がもたらされ、これは、上記に示したような式（１６ａ）によって数値的に解くことができる。式（１６ａ）は、式（１８））によるｂｉ＝ｂ（ｒｉ）＝ｎＡ（ｒｉ）及び式（１９）による
【数２９】
と共に用いることができる。ベクトルｂ、ｕ、及び行列Ｇは、画像のボクセルごとに形成され、行列Ｇは式（１４）に従って形成される。ｉ≠ｊの場合、Ｇｉｊの計算は、点ｒｉと点ｒｊとの間の光線（すなわち直線）に沿ってρａｂ（ｒｋ）Δｒｋを総和することを伴う。計算時間を削減するために、サンプリングは、各ｚスタックの光円錐の円盤面積にわたる平均ρ（ｒ）を計算するときに実行することができる。
【００４９】
代替的に、∂Ｊ／∂ρｋ，∀ｋの値は数値的に求めることができるので、式（１６ａ）は勾配降下法を用いて数値的に解くことができる。一例では、ρＡ（式（２０））は、勾配降下法においてρの初期推定値に用いられる。本発明の実施形態を用いて実行される数値シミュレーションを通じて、ρＡ（式（２０））がρの良好な近似値であることが分かる。ρＡ（式（２０））をρの初期推定値として用いると、勾配降下法におけるローカルミニマム問題が低減される。
【００５０】
３．側方散乱ジオメトリ
側方散乱ジオメトリは、実際には、単一平面照明顕微鏡法（ＳＰＩＭ）のジオメトリである［２２］、［２３］、［２４］、［２５］、［２６］。図６は、光源が側方から発してサンプル、すなわち図６の「サンプル」の１つの平面を照明する側方散乱のジオメトリ配置を示している。図６に示すように、散乱光は、ＣＣＤカメラによって直交方向で収集される。このジオメトリ配置では、入射光線は一定でかつ平行である。式（４）は、点ｒ＝（ｘ，ｙ）における一定の入射強度をｎ０と表記することによって以下の式（２４）に簡約することができる。式（２４）では、積分は、図６に示すように水平のｘ方向にわたる。共焦点顕微鏡法の場合のように、ρは、ρ＝ｑ−１ρｅｍ＝ρａｂとして設定することができる。
【数３０】
【００５１】
散乱光は、減衰することなくＣＣＤカメラに直接進むものと仮定される。ほとんどのカメラ設定と同様に、（ＣＣＤカメラにおける）ピクセル点ｒｐとサンプルロケーションｒとの間には１対１の対応関係がある。したがって、式（１１）は、（ρ＝ｑ−１ρｅｍ＝ρａｂが用いられるとき）以下に示すような式（２５）として記述することができ、ここで、α’は、全ての光線にわたる総和等を含む量子収量及び検出器の積分された効果を表す。
【数３１】
【００５２】
式（２５）を用いると、式（１６）の行列方程式は、式（２６）として記述することができる。
【数３２】
【００５３】
一例では、解析解は、散乱項が無視できると仮定することによって、式（２７）に示すように得られる。式（２７）では、下付き文字Ａは、近似解が減衰項のみを用いて得られることを示すのに用いられる。この近似を用いると、式（２７）を数値的により容易に解くことができる。式（２７）の総和は、光源ｒｓからそれぞれの点ｒｉへのレーザービームの方向の光線（すなわち直線）に沿って実行される。
【数３３】
【００５４】
本発明の実施形態の数値実験は、式（２６）を用いた強調画像及び式（１６ａ）を用いた強調画像が約１％しか異ならないことを示し、ｎｓ≪ｎＡの近似が有効であることを示していた。
【００５５】
４．結果
数値計算が実行され、その結果はグラウンド・トルース（ground truth）と比較される。他の物理学に基づく復元方法［５］、［６］、［７］、［８］、［９］、［１０］、［１１］、［１２］、［１３］と比較することは可能ではない。その理由は、これらの方法を顕微鏡画像に適用することはできないからである。第１に、他の物理学に基づく方法は、３次元画像を強調するように設計されていない。第２に、これらの方法は、一定減衰媒体を仮定しているが、この仮定は顕微鏡画像では大きく破られる。
【００５６】
４．１ 有効性確認及び較正
方法１００は、グラウンド・トルースが実験的設計によって判明している特別に調製されたサンプルに関して有効性が確認される。この場合、画像強調が方法１００を用いて実行され、得られた結果がグラウンド・トルースと比較される。実験では、サンプルが、フルオレセイン及び液体ゲルを、ゲルが硬化するまでオービタルシェーカー上で混合することによって作製される。このようにして、サンプルは、３Ｄの体積全体にわたって均一になる。しかしながら、取得された画像の強度プロファイルは、減衰のために均一とはならない。その代わり、強度プロファイルは深度と共に減少する。図７（ａ）に示すように、方法１００を用いて強調された画像のそれぞれは、原入力画像よりも均一な強度プロファイル（最大強度投影）を有する。強調画像は、図７では「復元（restored）」と表記されている。パラメータα’βｎ０は、顕微鏡に対して較正される。図７（ａ）に示すように、α’βｎ０＝１８１．２７が最良の結果を与える。他方、それよりも低い値１２１．５１及び９０．０２は過補償を招く。パラメータα’βｎ０の値をＣと表記すると、２つのパラメータ（Ｃ１，Ｃ２）と異なるレーザー強度（ｎ１，ｎ２）との間の関係は、単純にＣ１／Ｃ２＝ｎ１／ｎ２である。したがって、図７（ａ）の較正されたパラメータ値は、図７（ｃ）に示すように、異なるレーザー強度、例えば１．５ｎ０及び２．０ｎ０で撮った画像にも用いることができる。ここで、ｎ０は、図７（ａ）で用いられたレーザー強度である。図７（ｂ）には、原画像７０２及び（レーザー強度が１．５ｎ０であるときの）強調画像７０４の２Ｄ投影が示されているのに対して、図７（ｄ）には、原画像７０６及び（レーザー強度が２．０ｎ０であるときの）強調画像７０８の２Ｄ投影が示されている。強調画像はより均一に照明されることが分かる。
【００５７】
４．２ 共焦点顕微鏡法
方法１００の有効性を実証するために、核がＨｏｅｓｃｈｔ ３３３４２で染色された、マウスの胚からのニューロステムセルの３Ｄ画像が強調される。オリンパス社製の点走査共焦点（Point Scanning Confocal）ＦＶ１０００システムを用いて画像を取得した。開口数１．２を有する６０×の水レンズで撮像を行った。ダイオードレーザー４０ｎｍを用いて、Ｈｏｅｓｃｈｔで染色されたニューロスフェア（neurosphere）を励起させた。サンプリング速度は、２μｍ／ピクセルに設定した。原顕微鏡画像は、サイズがｘ方向及びｙ方向に０．１３７μｍの解像度並びにｚ方向に０．２μｍの解像度を有する５１２×５１２×ｎｚ個のボクセルである。ここで、ｎｚは、画像内のｚスタックの数である。
【００５８】
計算時間を低減するために、原画像は、ｚ方向の解像度を維持しながら、ｘ方向及びｙ方向でボクセルを平均することによって２５６×２５６×ｎｚ個のボクセルにダウンサンプリングされる。図８及び図９は、式（２０）を用いて強調された２５６×２５６×ｎｚ個のボクセルの画像の強調結果を示している。図８及び図９に示す共焦点顕微鏡画像は、それぞれ１５５個のｚスタック及び１６３個のｚスタック（すなわちｎｚ＝１５５及びｎｚ＝１６３）を有する。より具体的には、図８（ａ）及び図９（ａ）は、それぞれの原画像のｙｚ平面上への最大強度投影を示している（図８及図９に原画像と表記されている）のに対して、図８（ｂ）及び図９（ｂ）は、それぞれの強調画像のｙｚ平面上への最大強度投影を示している（図８及び図９に復元図と表記されている）。図８（ｂ）及び図９（ｂ）において、調整されたチューニングパラメータは、１／α’βｎ０＝０．０１４９９５として設定される。これは、最適な強調結果を与えることが分かる。図８（ｃ）及び図９（ｃ）は、原画像（実線）及び強調画像（破線）の双方のｚ軸上への（ｘｙ平面における最も明るい０．１％のボクセルにわたって平均された）最大強度投影を示している。照明レーザーは底部から発するので、原画像（図８（ａ）又は図９（ａ））について、ボクセルが画像の底部ではるかに明るく、画像の強度が画像の上部に向かって低下する（換言すれば、照明は均一でない）ことが図８及び図９から容易に観察することができる。しかしながら、図８（ｂ）及び図９（ｂ）に示すように、照明は、強調後、均一になる。さらに、図８（ｃ）及び図９（ｃ）は、原画像（実線）の強度プロファイルと強調画像（破線）の強度プロファイルとの間の相違を明らかに示している。加えて、図８及び図９に示すように、底部のｚスタックの露出過度のエリアも、強調方法１００によって正しく補償される。したがって、強調画像の均一な照明を達成することができるので、方法１００を用いて画像を復元することは有利であることが図８及び図９から分かる。その後、この強調画像に対してヒストグラム等化等の他の画像強調方法を用いることができる。強調画像は、平均してより暗くもなるが、多くの画像処理技法は、平均ボクセル強度に対してロバストである。
【００５９】
４．３ 側方散乱顕微鏡法
合成劣化画像に対して側方散乱ジオメトリの計算を行った。不均一な照明（すなわち、光源が左から発すると仮定して最大強度投影が指数関数的に低下する）による合成画像を均一な照明の画像から生成した。図１０は、２５６×２５６ピクセル画像の強調結果を示している（強調画像は、図１０では「復元（Restored）」とラベル付けされている）。図１０において、式（２７）は、画像を強調するのに用いられる。ｎ０は、入射光強度であり、α’は、通例未知であるジオメトリックファクターである。チューニングパラメータ１／α’ｎ０は、最適な結果を得るように調整することができる。図１１（ａ）は、小さな１／α’ｎ０（０．００１１）が用いられたときの強調画像を示しているのに対して、図１１（ｂ）は、大きな１／α’ｎ０（０．０１１１）が用いられたときの強調画像を示している。図１１（ａ）及び図１１（ｂ）から分かるように、小さな１／α’ｎ０が用いられたとき、画像はほとんど強調されず、大きな１／α’ｎ０が用いられたとき、減衰効果が過補償される。１／α’ｎ０の最適値は、図１０の強調画像を得るのに用いられた０．００９５である。図１０に示すように、１／α’ｎ０が０．００９５に設定されたとき、強調画像はほぼ完全に（均一に）照明される。
【００６０】
上述したように、方法１００は、均一な照明を有する強調画像を得ることができるので有利である。換言すれば、方法１００を用いて画像を強調することによって、光学顕微鏡法の基本的な光減衰及び散乱の問題を軽減することができる。
【００６１】
方法１００で用いられる式ｂ＝Ｇ・ｕの微分は、強い理論的基礎に基づいて定式化されており、連続の式によって表わされる保存則等の物理学の基本法則に基づいている。さらに、場の理論的手法（field theoretical approach）がこの微分に用いられている。
【００６２】
方法１００は、物理学に基づく復元方法のタイプであり、物理学に基づく復元方法は、モデルに基づくコントラスト強調方法（例えば、ヒストグラム等化）を上回る多くの利点を有する。モデルに基づく方法［１５］、［１６］、［１７］は、一般に、画像特性が画像全体にわたって一定であると仮定し、この仮定は、天候で劣化した画像では無視される。さらに、物理モデルは物理学の法則に基づいて構築され、物理学は、否定できない真実である可能性が最も高い。物理学に基づく復元技法は、多くの用途で用いることができる。そのような復元技法の１つの態様は、その有効性が数桁の物理的な長さスケールにわたることである。航空機による監視では、物理的な長さスケールはほぼ≒１０ｋｍ程度であり、水中監視では、物理的な長さスケールはほぼ≒１０ｍ程度である。
【００６３】
物理学に基づく復元方法が天候により劣化した画像の復元に用いられてきたが、それらの方法は、（≒１００μｍの物理的な長さスケールを有する）光学顕微鏡法の画像強調の分野で適切に探求されていなかった。方法１００は、顕微鏡画像に対して用いられる物理学に基づく復元技法のタイプであり、物理学に基づく復元の長さスケールを８桁に拡大する。
【００６４】
方法１００は、物理学に基づく復元技法のタイプであっても、全ての既存の物理学に基づく復元技法とは根本的に異なっている。既存の物理学に基づく復元技法では、減衰媒体の一定の吸収係数が前提とされるのに対して、このことは、方法１００では前提とされない。さらに、方法１００では、サンプルと減衰媒体とが区別されない。一般的な一組の式が導出され、画像取得におけるあらゆるジオメトリ設定を取り扱うために方法１００で用いられる。方法１００を用いるには、光源及びカメラ等の検出機器の詳細を指定することだけが必要である。他方、既存の物理学に基づく方法［５］、［６］、［７］、［８］、［９］、［１０］、［１１］、［１２］、［１３］は、以下の理由によって３次元顕微鏡画像に適用することさえできない。第１に、既存の物理学に基づく方法は、減衰媒体を「除去」して、２次元シーンをリトリーブする。これとは対照的に、方法１００では、減衰媒体も画像情報を含む。これは、媒体を除去するのではなく媒体の真の信号を復元することが必要であるので有利である。第２に、既存の方法は均一な減衰媒体を仮定しているが、この仮定は、顕微鏡画像では大きく破られる。これとは対照的に、そのような仮定は方法１００ではなされない。
【００６５】
［参考文献］
【表１】
【表２】
【表３】
【特許請求の範囲】
【請求項１】
サンプルの顕微鏡画像を強調する方法であって、該顕微鏡画像は、該サンプルを或る照明方向で照明しカメラを用いて該サンプルにより散乱された光を捕捉することによって形成されたものであり、前記顕微鏡画像のそれぞれの構成部分は、該サンプルのそれぞれの要素から捕捉された散乱光の量を表しており、
該方法は、
（ｉ）前記顕微鏡画像の構成要素と、前記サンプルの複数のそれぞれの要素の散乱パラメータの値とを関連させる数式を用いて該散乱パラメータの前記値を取得するステップであって、前記サンプルの各要素の前記散乱パラメータのそれぞれの前記値は、前記サンプルの当該各要素についての、入射光を散乱する傾向を示すものである、取得するステップと、
（ｉｉ）前記散乱パラメータの前記取得された値を用いて、前記サンプルの強調画像を形成するステップと
を含んでなる方法。
【請求項２】
前記サンプルの所与の前記要素ごとに、前記数式は、該所与の前記要素の前記散乱パラメータの前記値を、該所与の前記要素として前記照明方向にある要素の前記散乱パラメータの前記値に関して表しており、
上記ステップ（ｉ）は、要素について連続的に、さらに前記照明方向でも連続的に前記散乱パラメータの前記値を取得することを含むものである、請求項１に記載の方法。
【請求項３】
前記顕微鏡画像は、共焦点顕微鏡法又は単一平面照明顕微鏡法を用いて取得される、請求項１に記載の方法。
【請求項４】
前記強調画像は、各構成部分が前記サンプルのそれぞれの要素に対応する画像であり、その要素の前記散乱パラメータの前記取得された値に対応する強度を有する、請求項１〜３のいずれか１項に記載の方法。
【請求項５】
前記画像は、前記散乱パラメータの前記値を取得するステップの前に線形に正規化される、請求項１〜４のいずれか１項に記載の方法。
【請求項６】
前記画像は、前記散乱パラメータの前記値を取得するステップの前にダウンサンプリングされる、請求項１〜５のいずれか１項に記載の方法。
【請求項７】
前記強調画像のピクセルの強度を０〜（２ｎ−１）の範囲に再スケーリングするステップをさらに含み、ここで、ｎは、前記ピクセルを表すのに用いられるビット数である、請求項１〜６のいずれか１項に記載の方法。
【請求項８】
前記数式はチューナブルパラメータを含み、前記方法は、前記強調画像に実質的に一定の平均強度を与える前記チューナブルパラメータの値を選択することを含む、請求項１〜７のいずれか１項に記載の方法。
【請求項９】
前記数式は、前記サンプルの各要素における前記光強度が、前記サンプルの他の要素からの散乱による無視できる成分のみを含むとの仮定と矛盾しない、請求項１〜８のいずれか１項に記載の方法。
【請求項１０】
前記数式は、ｂ＝Ｇ・ｕの形式であり、ここで、ｂ及びｕは、前記サンプルを含む３次元空間内の複数のそれぞれの点の各点についての成分を有するベクトルであり、ｂは、減衰後に残っている入射光の振幅を表すデータ値を含み、ｕは、各点が散乱光を生成する度合いを表すデータ値を含み、Ｇは、前記散乱パラメータを組み込んだ行列である、請求項１〜９のいずれか１項に記載の方法。
【請求項１１】
前記数式は、１つ以上の平均パラメータを用いることによって、前記サンプルの所与の前記要素の前記散乱パラメータの前記値を表し、各平均パラメータは、前記照明方向に平行に前記サンプルの前記所与の要素まで延びる線を取り囲む対応する領域にわたる前記散乱パラメータの前記値の平均を示す、請求項１〜１０のいずれか１項に記載の方法。
【請求項１２】
前記照明は、光をレンズに透過させて同じレンズを通る前記散乱光を収集することによって行われ、前記数式は、前記レンズと前記サンプルの前記所与の要素との間の円盤である複数の前記領域のそれぞれの前記平均パラメータを用い、該円盤のそれぞれは前記レンズに平行である、請求項１１に記載の方法。
【請求項１３】
前記サンプルは平面サンプルであり、該平面サンプルは、該平面サンプルの平面を或る照明方向で照明され、前記散乱光は、前記サンプルの前記平面を横切る方向で該サンプルから離間したカメラによって収集される、請求項１〜１０のいずれか１項に記載の方法。
【請求項１４】
請求項１〜１３のいずれか１項に記載の方法を実行するように構成されたプロセッサを有するコンピュータシステム。
【請求項１５】
有形のデータストレージデバイス等のコンピュータプログラム製品であって、コンピュータが読み取ることができるものであり、請求項１〜１３のいずれか１項に記載の方法を実行させるために、コンピュータシステムのプロセッサによって動作可能な命令を含む、コンピュータプログラム製品。
【請求項１】
サンプルの顕微鏡画像を強調する方法であって、該顕微鏡画像は、該サンプルを或る照明方向で照明しカメラを用いて該サンプルにより散乱された光を捕捉することによって形成されたものであり、前記顕微鏡画像のそれぞれの構成部分は、該サンプルのそれぞれの要素から捕捉された散乱光の量を表しており、
該方法は、
（ｉ）前記顕微鏡画像の構成要素と、前記サンプルの複数のそれぞれの要素の散乱パラメータの値とを関連させる数式を用いて該散乱パラメータの前記値を取得するステップであって、前記サンプルの各要素の前記散乱パラメータのそれぞれの前記値は、前記サンプルの当該各要素についての、入射光を散乱する傾向を示すものである、取得するステップと、
（ｉｉ）前記散乱パラメータの前記取得された値を用いて、前記サンプルの強調画像を形成するステップと
を含んでなる方法。
【請求項２】
前記サンプルの所与の前記要素ごとに、前記数式は、該所与の前記要素の前記散乱パラメータの前記値を、該所与の前記要素として前記照明方向にある要素の前記散乱パラメータの前記値に関して表しており、
上記ステップ（ｉ）は、要素について連続的に、さらに前記照明方向でも連続的に前記散乱パラメータの前記値を取得することを含むものである、請求項１に記載の方法。
【請求項３】
前記顕微鏡画像は、共焦点顕微鏡法又は単一平面照明顕微鏡法を用いて取得される、請求項１に記載の方法。
【請求項４】
前記強調画像は、各構成部分が前記サンプルのそれぞれの要素に対応する画像であり、その要素の前記散乱パラメータの前記取得された値に対応する強度を有する、請求項１〜３のいずれか１項に記載の方法。
【請求項５】
前記画像は、前記散乱パラメータの前記値を取得するステップの前に線形に正規化される、請求項１〜４のいずれか１項に記載の方法。
【請求項６】
前記画像は、前記散乱パラメータの前記値を取得するステップの前にダウンサンプリングされる、請求項１〜５のいずれか１項に記載の方法。
【請求項７】
前記強調画像のピクセルの強度を０〜（２ｎ−１）の範囲に再スケーリングするステップをさらに含み、ここで、ｎは、前記ピクセルを表すのに用いられるビット数である、請求項１〜６のいずれか１項に記載の方法。
【請求項８】
前記数式はチューナブルパラメータを含み、前記方法は、前記強調画像に実質的に一定の平均強度を与える前記チューナブルパラメータの値を選択することを含む、請求項１〜７のいずれか１項に記載の方法。
【請求項９】
前記数式は、前記サンプルの各要素における前記光強度が、前記サンプルの他の要素からの散乱による無視できる成分のみを含むとの仮定と矛盾しない、請求項１〜８のいずれか１項に記載の方法。
【請求項１０】
前記数式は、ｂ＝Ｇ・ｕの形式であり、ここで、ｂ及びｕは、前記サンプルを含む３次元空間内の複数のそれぞれの点の各点についての成分を有するベクトルであり、ｂは、減衰後に残っている入射光の振幅を表すデータ値を含み、ｕは、各点が散乱光を生成する度合いを表すデータ値を含み、Ｇは、前記散乱パラメータを組み込んだ行列である、請求項１〜９のいずれか１項に記載の方法。
【請求項１１】
前記数式は、１つ以上の平均パラメータを用いることによって、前記サンプルの所与の前記要素の前記散乱パラメータの前記値を表し、各平均パラメータは、前記照明方向に平行に前記サンプルの前記所与の要素まで延びる線を取り囲む対応する領域にわたる前記散乱パラメータの前記値の平均を示す、請求項１〜１０のいずれか１項に記載の方法。
【請求項１２】
前記照明は、光をレンズに透過させて同じレンズを通る前記散乱光を収集することによって行われ、前記数式は、前記レンズと前記サンプルの前記所与の要素との間の円盤である複数の前記領域のそれぞれの前記平均パラメータを用い、該円盤のそれぞれは前記レンズに平行である、請求項１１に記載の方法。
【請求項１３】
前記サンプルは平面サンプルであり、該平面サンプルは、該平面サンプルの平面を或る照明方向で照明され、前記散乱光は、前記サンプルの前記平面を横切る方向で該サンプルから離間したカメラによって収集される、請求項１〜１０のいずれか１項に記載の方法。
【請求項１４】
請求項１〜１３のいずれか１項に記載の方法を実行するように構成されたプロセッサを有するコンピュータシステム。
【請求項１５】
有形のデータストレージデバイス等のコンピュータプログラム製品であって、コンピュータが読み取ることができるものであり、請求項１〜１３のいずれか１項に記載の方法を実行させるために、コンピュータシステムのプロセッサによって動作可能な命令を含む、コンピュータプログラム製品。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【公表番号】特表２０１２−５１９８７６（Ｐ２０１２−５１９８７６Ａ）
【公表日】平成２４年８月３０日（２０１２．８．３０）
【国際特許分類】
【出願番号】特願２０１１−５５２９１３（Ｐ２０１１−５５２９１３）
【出願日】平成２２年２月１１日（２０１０．２．１１）
【国際出願番号】ＰＣＴ／ＳＧ２０１０／００００５１
【国際公開番号】ＷＯ２０１０／１０１５２５
【国際公開日】平成２２年９月１０日（２０１０．９．１０）
【出願人】（５０４３５４１１７）エイジェンシー・フォー・サイエンス，テクノロジー・アンド・リサーチ (17)
【Ｆターム（参考）】
【公表日】平成２４年８月３０日（２０１２．８．３０）
【国際特許分類】
【出願日】平成２２年２月１１日（２０１０．２．１１）
【国際出願番号】ＰＣＴ／ＳＧ２０１０／００００５１
【国際公開番号】ＷＯ２０１０／１０１５２５
【国際公開日】平成２２年９月１０日（２０１０．９．１０）
【出願人】（５０４３５４１１７）エイジェンシー・フォー・サイエンス，テクノロジー・アンド・リサーチ (17)
【Ｆターム（参考）】
[ Back to top ]