JPWO2021024420A1 - 屈折率分布推定システム - Google Patents

Abstract

光軸と直交する面内における分解能が高く、標本以外の情報の影響が少ない屈折率分布推定システムを提供する。屈折率分布推定システムは、標本を照明する照明光学系と、標本の光学像を形成する結像光学系と、標本の光学像から撮影画像を取得する撮像素子と、撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、プロセッサは、推定標本を推定するステップと、モデル化した複数の光源から射出された複数の第１波面から結像位置における複数の第１強度分布を計算し、複数の第１強度分布から推定標本の像を算出するステップと、複数の第１波面が推定標本を通過した後の複数の第２波面、撮影画像及び推定標本の像から推定標本の屈折率分布を最適化するステップと、推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、推定標本の構造を再構成して出力するステップと、を含む処理を行う。

Description

本発明は、屈折率分布推定システムに関する。

実際の物体を計算機上の物体モデルで再現する再構成手法がある。この再構成手法では、測定した物体の画像と計算した物体モデルの画像が一致するように、最適化手法で計算機上の物体モデルを変更していく。最終的に、物体の画像と物体モデルの画像が一致した時に、計算機上の物体モデルは実際の物体を再現している。

物体の画像は、測定光学系で取得される。物体モデルの画像がどのような画像になるかは、像演算手法で計算される。そのため、この再構成手法では、測定光学系と、像演算手法の２つが重要である。

測定光学系として、例えば、顕微鏡の光学系を用いることができる。顕微鏡の光学系では、ハロゲンランプ又はＬＥＤを用いて標本の画像を取得する。ハロゲンランプとＬＥＤは、インコヒーレント光源である。

インコヒーレント光源を用いた照明は、照明条件によって、インコヒーレント照明、コヒーレント照明、及びパーシャルコヒーレント照明に分類される。これらの照明について説明する。

顕微鏡では、ケーラー照明が用いられる。ケーラー照明では、光源がコンデンサレンズの焦点面上に配置されるか、又は光源の像がコンデンサレンズの焦点面上に形成される。光源の各点から出た光は、コンデンサレンズで平行光線に変換される。よって、標本は平行光束で照明される。

光源の大きさを変えると、標本面での照明光の空間的コヒーレンスが変化する。照明光の空間的コヒーレンスが変化することで、結像特性が変化する。

インコヒーレント光源であっても、光源の大きさを非常に小さくすると、光源は点光源とみなせるようになる。点光源からの光を標本に照射する照明は、コヒーレント照明と呼ばれる。

光源の大きさを非常に大きくすると、光源は面光源とみなせるようになる。特に、対物レンズの瞳の全域を照明できるように光源の大きさを大きくすると、標本に様々な方向から光源からの光を照射できる。このような照明は、インコヒーレント照明と呼ばれる。

面光源において、対物レンズの瞳の一部を照明できるように、光源の大きさを設定する。このような照明は、パーシャルコヒーレント照明（部分的コヒーレント照明）と呼ばれる。パーシャルコヒーレント照明は、コヒーレント照明とインコヒーレント照明の中間の照明である。

像演算手法は、線形演算と非線形演算に分類できる。線形演算では、標本内の１回の散乱のみが考慮される。非線形演算では、１回の散乱だけでなく、多重散乱も考慮される。

線形演算では、１次ボルン近似を利用している。ボルン近似では、標本で２回以上の散乱が生じても、２回以上の散乱は無視される。線形演算では、標本情報と出力情報が一対一の関係に決まる。そのため、解析的に出力情報を計算できる。出力情報は、例えば、標本の像である。

線形演算について、顕微鏡での結像を例に説明する。標本情報（標本の透過率分布）Ｏと出力情報（像強度分布）Ｉを光学系の点像強度分布ＰＳＦのコンボリューションの線形システムとみなせるとき、出力情報Ｉは、以下の式で表される。
Ｉ＝ＰＳＦ＊Ｏ
ここで、＊は、コンボリューションを表す。

線形演算では、計算時間は短いが、２回以上の散乱は無視しているため計算精度は低くなる。線形演算を利用して再構成した物体モデルの像は、測定した標本の像を点像強度分布でデコンボリューションして求める。

非線形演算は、標本で複数回の散乱が生じることを考慮した演算手法である。非線形演算の１つに、ビーム伝搬法がある。ビーム伝搬法では、物体モデルを複数の薄い層に置き換える。そして、光が各層を通過する際の波面変化を逐次計算して、物体モデルの像を計算する。

ビーム伝搬法は、線形演算に比べて、より正確に、物体モデルの像を計算できる。

非特許文献１に、測定光学系の照明に顕微鏡のコヒーレント照明を用い、像演算手法にビーム伝搬法を用いて物体の再構成を行う手法が提案されている。

この手法では、ＬＥＤアレイが配置された顕微鏡が用いられている。ＬＥＤアレイは、照明光学系の瞳位置に配置されている。ＬＥＤの点灯位置を変えることで、標本に対して様々な角度で照明が照射される。各照射角度で標本の画像を取得するので、複数の標本の画像が取得される。

取得した標本の画像から、ライトフィールド手法を用いて、標本の屈折率分布の初期値を求める。ライトフィールド手法は、再構成手法の１つである。偏射照明では、ピント位置からの距離に比例して、標本の像の位置がずれる。ライトフィールド手法は、この標本の像の位置ずれを利用して標本を再構成手法する。

初期値が求まったら、勾配降下法を用いて、標本の屈折率分布を正しい値に近づける。勾配降下法は、光源側から順方向に標本内を伝搬した波面と、像情報で補正した波面を像側から逆方向に伝播した波面との差が小さくなるように標本のパラメータを更新していく手法である。波面の伝播計算にはビーム伝搬法が用いられる。

"3D intensity and phase imaging from light field measurements in an LED array microscope," L. Tian and L. Waller, Optica 2, 104-111 （2015）.

コヒーレント照明を用いた結像光学像では、光軸と直交する面内における分解能が、パーシャルコヒーレント照明を用いた結像光学系に比べて低い。さらに被写界深度が深く標本以外の情報、例えばカバーガラスに付着したゴミも標本の画像として取得されてしまう。そのため、再構成した像における屈折率分布の精度が悪くなる。

本発明は、このような課題に鑑みてなされたものであって、光軸と直交する面内における分解能が高く、標本以外の情報の影響が少ない屈折率分布推定システムを提供することを目的とする。

上述した課題を解決し、目的を達成するために、本発明の少なくとも幾つかの実施形態に係る屈折率分布推定システムは、
パーシャルコヒーレント照明により標本を照明する照明光学系と、
標本の光学像を形成する結像光学系と、
結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、
撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、
プロセッサは、
標本の屈折率分布を含む推定標本を推定するステップと、
照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、結像光学系の結像位置における複数の第１強度分布を計算し、複数の第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、
複数の第１波面が推定標本を通過した後の複数の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、
推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、
更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする。

本発明の少なくとも幾つかの実施形態に係る屈折率分布推定システムは、
複数の方向から同時に入射する光線により標本を照明する照明光学系と、
標本の光学像を形成する結像光学系と、
結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、
撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、
プロセッサは、
標本の屈折率分布を含む推定標本を推定するステップと、
照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、結像光学系の結像位置における複数の第１強度分布を計算し、複数の第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、
複数の第１波面が推定標本を通過した後の複数の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、
推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、
更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする。

本発明によれば、光軸と直交する面内における分解能が高く、標本以外の情報の影響が少ない屈折率分布推定システムを提供することができる。

本実施形態の屈折率分布推定システムと撮影画像を示す図である。 照明光と照明領域を示す図である。 第１のシミュレーションのフローチャートである。 分割された照明領域を示す図である。 シミュレーションで用いる光学系を示す図である。 推定標本の像を示す図である。 波面の補正を示す図である。 標本の勾配を示す図である。 標本の勾配を示す図である。 本実施形態の屈折率分布推定システムと照明光を示す図である。 撮影画像を示す図である。 第２のシミュレーションのフローチャートである。 第２のシミュレーションのフローチャートである。 シミュレーションで用いる光学系を示す図である。 各層における波面を示す図である。 撮影画像の取得位置における波面と結像面における波面を示す図である。 推定標本の像を示す図である。 波面の補正を示す図である。 標本の勾配と波面の伝搬を示す図である。 標本の勾配と波面の伝搬を示す図である。 標本の勾配を示す図である。 標本の勾配を示す図である。 本実施形態の屈折率分布推定システムを示す図である。 第１開口部材を示す図である。 本実施形態の屈折率分布推定システムを示す図である。 第３のシミュレーションのフローチャートである。 開口部材と撮影画像を示す図である。 標本、撮影画像、及び初期値の画像を示す図である。 伝達可能な空間周波数帯域を示す図である。 標本を示す図である。 撮影画像、伝達可能な空間周波数を示す図、再生像の全体の画像、及び再生像の一部の画像である。 標本を示す図と再生像の画像である。 標本を示す図、初期値の画像、及び再生像の画像である。 標本を示す図と再生像の画像である。

実施例の説明に先立ち、本発明のある態様にかかる実施形態の作用効果を説明する。なお、本実施形態の作用効果を具体的に説明するに際しては、具体的な例を示して説明することになる。しかし、後述する実施例の場合と同様に、それらの例示される態様はあくまでも本発明に含まれる態様のうちの一部に過ぎず、その態様には数多くのバリエーションが存在する。したがって、本発明は例示される態様に限定されるものではない。

本実施形態の屈折率分布推定システムでは、撮影画像と推定標本の像が用いられる。撮影画像は、光学装置で取得した標本の画像である。推定標本の像は、シミュレーションで取得した推定標本の画像である。

シミュレーションでは、残差が閾値以下になるように、推定標本の更新が行われる。残差は、推定標本の像と撮影画像との差である。残差が閾値以下になると、撮影画像と同一の推定標本の像が得られるか、又は、撮影画像と略同一の推定標本の像が得られる。

撮影画像の取得では、パーシャルコヒーレント照明が用いられている。よって、シミュレーションも、パーシャルコヒーレント照明を前提としている。

第１のシミュレーションでは、標本は厚みが薄い標本（以下、「薄い標本」という）である。第２のシミュレーションと第３のシミュレーションでは、標本は厚みが厚い標本（以下、「厚い標本」という）である。

本実施形態の屈折率分布推定システムは、パーシャルコヒーレント照明により標本を照明する照明光学系と、標本の光学像を形成する結像光学系と、結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、プロセッサは、標本の屈折率分布を含む推定標本を推定するステップと、照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、結像光学系の結像位置における複数の第１強度分布を計算し、複数の第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、複数の第１波面が推定標本を通過した後の複数の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする。

図１は本実施形態の屈折率分布推定システムと撮影画像を示す図である。図１（ａ）は、屈折率分布推定システムを示す図である。図１（ｂ）は、撮影画像を示す図である。

図１（ａ）に示すように、屈折率分布推定システム１は、照明光学系２と、結像光学系３と、撮像素子４と、プロセッサ５と、を有する。照明光学系２と結像光学系３で、測定光学系が形成されている。

照明光学系２と結像光学系３との間に、標本６が位置している。標本６は、薄い標本である。結像光学系３の焦点位置Ｆｏは、標本６の内部に位置している。例えば、焦点位置Ｆｏと標本６の表面６ａとの間隔はΔｚ１である。

標本６には、複数の方向から、光線が同時に入射する。図１（ａ）では、光線Ｌ_ＯＮと光線Ｌ_ＯＦＦが示されている。光線Ｌ_ＯＮは、光軸ＡＸと平行な光線である。光線Ｌ_ＯＦＦは、光軸ＡＸと交差する光線である。

このように、複数の方向から同時に入射する光線によって、標本６が照明されている。このような照明としては、インコヒーレント照明とパーシャルコヒーレント照明がある。屈折率分布推定システム１では、パーシャルコヒーレント照明が用いられている。

光のコヒーレンスには、時間的コヒーレンスと空間的コヒーレンスが含まれる。ここでのコヒーレンスは、空間的コヒーレンスを指している。

図２は、照明光と照明領域を示す図である。図２（ａ）は、照明光と照明領域の第１例を示す図である。図２（ｂ）は、照明光と照明領域の第２例を示す図である。各図では、標本に入射する照明光と、結像光学系の瞳位置における照明領域が示されている。

第１例では、光線Ｌ_ＯＮと光線Ｌ_ＯＦＦが照明に用いられている。結像光学系３の瞳位置Ｐｕでは、照明領域ＩＬＬの形状は円である。照明領域ＩＬＬは、結像光学系３の瞳Ｐ_ＯＢよりも狭い。よって、第１例における照明は、パーシャルコヒーレント照明である。

第２例では、光線Ｌ_ＯＦＦだけが照明に用いられている。瞳位置Ｐｕでは、照明領域ＩＬＬの形状は円環である。照明領域ＩＬＬは、瞳Ｐ_ＯＢよりも狭い。よって、第２例における照明は、パーシャルコヒーレント照明である。

図１に戻って説明する。標本６から射出した光は、結像光学系３によって、結像面ＩＰに集光される。結像面ＩＰに、光学像６’が形成される。光学像６’は、標本６の光学像である。

結像面ＩＰには、撮像素子４の撮結像面が位置している。撮像素子４によって、光学像６’の画像の取得が行われる。その結果、図１（ｂ）に示す撮影画像Ｉ_ｍｅａ（ｒ）が得られる。ｒは、（ｘ，ｙ）の２次元座標を示す。

標本６は、薄い標本なので、１つの撮影画像が取得される。よって、結像光学系３と撮像素子４は、光軸方向に移動しない。また、標本６も、光軸方向に移動しない。

撮影画像Ｉ_ｍｅａ（ｒ）は、プロセッサ５に入力される。プロセッサ５では、撮影画像Ｉ_ｍｅａ（ｒ）を用いて、推定標本の像のシミュレーションが行われる。

図３は、第１のシミュレーションのフローチャートである。図１に示す測定光学系では、照明領域は無数の点光源で形成されている。シミュレーションでは、無数の点光源を扱うことはできない。そこで、照明領域を、複数の微小領域に分割する。

微小領域の各々は、点光源とみなすことができる。照明領域を複数の微小領域に分割することで、点光源の数を少なくすることができる。

図４は、分割された照明領域を示す図である。図４（ａ）は、分割された照明領域の第１例を示す図である。図４（ｂ）は、分割された照明領域の第２例を示す図である。

第１例では、円形の照明領域が、複数の微小領域に分割されている。第２例では円環形の照明領域が、複数の微小領域に分割されている。

上述のように、微小領域の各々は、点光源とみなすことができる。よって、第１例では、複数の点光源によって、円形の照明領域が形成されている。第２例では、複数の点光源によって、円環形の照明領域が形成されている。

図５は、シミュレーションで用いる光学系を示す図である。シミュレーションで用いる光学系は、撮影画像Ｉ_ｍｅａ（ｒ）を取得した測定光学系と同一である。シミュレーションでは、標本６の代わりに、推定標本１０が用いられる。

図５には、推定標本１０、波面ｆ_ｉｎ ^ｍ（ｒ）、振幅透過率Ｔ_ｓ（ｒ）、波面ｇ_ｏｕｔ ^ｍ（ｒ）、撮影画像の取得位置における波面ｕ^ｍ（ｒ）、及び結像面における波面ｕ_ｉｍｇ ^ｍ（ｒ）が図示されている。

図５では、１番目の光源からＮ_ＬＳ番目までの光源が図示されている。光源は、照明光学系２の瞳位置に配置することができる。このようにすることで、照明光学系２の瞳の強度分布が、複数の光源でモデル化することができる。

図３に戻って、シミュレーションについて説明する。シミュレーションは、推定標本を推定するステップと、推定標本の像を算出するステップと、推定標本の屈折率分布を最適化するステップと、推定標本を更新するステップと、推定標本の構造を再構成して出力するステップと、を含む。

ステップＳ１０では、光源の数Ｎ_ＬＳが設定される。光源の数Ｎ_ＬＳは、照明領域を微小領域に分割したときの分割数である。シミュレーションでは、図４（ａ）に示す円形の照明領域が用いられる。

ステップＳ２０は、推定標本を推定するステップである。標本６については、１枚の撮影画像が取得されている。推定標本１０は薄い標本なので、１つの薄い層とみなすことができる。よって、振幅透過率の初期値の設定は１回行われる。

ステップＳ２０では、推定標本１０における振幅透過率Ｔ_ｓ（ｒ）に初期値が設定される。

推定標本１０の像を算出するためには、推定標本１０の情報、例えば、屈折率分布が必要である。推定標本１０は、標本６をモデル化した標本である。よって、推定標本１０の屈折率分布に標本６の屈折率分布を用いることができると良い。

しかしながら、標本６の屈折率分布は、撮影画像Ｉ_ｍｅａ（ｒ）から正確に得られない。よって、推定標本１０の屈折率分布は、推定せざるを得ない。

式（１）に示すように、結像面における屈折率分布ｎ_ｓ（ｒ）は、振幅透過率Ｔ_ｓ（ｒ）に変換することができる。よって、ステップＳ２０では、推定標本１０における振幅透過率Ｔ_ｓ（ｒ）の初期値を設定する。

ここで、
λは、照明光の波長、
ｄｚは、標本の厚み、
である。

振幅透過率Ｔ_ｓ（ｒ）の値を撮影画像Ｉ_ｍｅａ１（ｒ）から推定できる場合、推定した値を初期値に用いても良い。また、他の方法で振幅透過率Ｔ_ｓ（ｒ）の値を推定できる場合、推定した値を初期値に設定することができる。初期値が推定できない場合、例えば、Ｔ_ｓ（ｒ）＝１とする。

ステップＳ３０では、変数ｍの値が初期化される。後述のステップＳ４１、ステップＳ４２、ステップＳ４３、ステップＳ４４、及びステップＳ４５は、全ての光源に対して実行される。変数ｍは、これらのステップが実行された回数を表している。

ステップＳ４０とステップＳ５０は、推定標本の像を算出するステップである。推定標本の像の数は、撮影画像の数と同数である。撮影画像の数は１なので、推定標本の像の数も１である。

ステップＳ４０は、ステップＳ４１、ステップＳ４２、ステップＳ４３、ステップＳ４４、ステップＳ４５、ステップＳ４６、及びステップＳ４７を有する。

ステップＳ４１では、推定標本１０へ入射する波面ｆ_ｉｎ ^ｍ（ｒ）が算出される。ｆ_ｉｎ ^ｍ（ｒ）は、１番目の光源からＮ_ＬＳ番目までの光源から射出された光の波面を表している。

上述のように、照明光学系の瞳の強度分布は、複数の光源でモデル化することができる。１番目の光源からＮ_ＬＳ番目までの光源は、モデル化された複数の光源である。第１波面を、モデル化された複数の光源から射出された波面とすると、波面ｆ_ｉｎ ^ｍ（ｒ）は、第１波面を表している。

上述のように、微小領域の各々は、点光源とみなすことができる。図５では、ｍ番目の光源から照明光Ｌｍが出射されている。照明光Ｌｍは、推定標本１０に入射する。

この場合、波面ｆ_ｉｎ ^ｍ（ｒ）は、式（２）と式（３）で表される。

ここで、
θ_ｘ，ｍ、θ_ｙ，ｍは、推定標本への入射角、
である。

ステップＳ４２では、推定標本１０から射出される波面ｇ_ｏｕｔ ^ｍ（ｒ）が算出される。薄い標本の場合、波面ｇ_ｏｕｔ ^ｍ（ｒ）は、式（４）で表される。

ここで、
Ｔ_ｓ（ｒ）は、推定標本の振幅透過率、
である。

波面ｇ_ｏｕｔ ^ｍ（ｒ）は、波面ｆ_ｉｎ ^ｍ（ｒ）が推定標本１０を通過した後の波面である。波面ｆ_ｉｎ ^ｍ（ｒ）は第１波面を表しているので、波面ｇ_ｏｕｔ ^ｍ（ｒ）は第２波面を表している。

推定標本１０は薄い標本なので、式（４）に示すように、波面ｆ_ｉｎ ^ｍ（ｒ）から、波面ｇ_ｏｕｔ ^ｍ（ｒ）を直接算出することができる。

ステップＳ４３では、撮影画像の取得位置における波面ｕ^ｍ（ｒ）が算出される。撮影画像の取得位置は、撮影画像の取得が行われたときの、標本側における結像光学系３の焦点位置Ｆｏである。

波面ｕ^ｍ（ｒ）は、式（５）で表される。

ここで、
Ｐ_−Δｚ１｛ ｝は、式（６）で表され、

Δｚ１は、推定標本の表面から撮影画像の取得位置までの距離、
λは、波長、
ｕは、瞳面座標（ξ，η）の２次元表記、
Ｆ_２Ｄは、２次元フーリエ変換、
Ｆ_２Ｄ ^−１は、２次元フーリエ逆変換、
である。

後述のステップＳ６０では、残差が算出される。残差の算出は、撮影画像と推定標本の像が用いられる。推定標本の像を算出するためには、撮影画像の取得位置における波面を求める必要がある。

上述のように、焦点位置Ｆｏと表面６ａとの間隔はΔｚ１である。光の進行方向に向かって測定したときの距離の符号をプラスとすると、撮影画像の取得位置は、表面６ａから−Δｚ１だけ離れた位置になる。

よって、シミュレーションで用いる光学系では、撮影画像の取得位置は、推定標本１０の表面１０ａから−Δｚ１だけ離れた位置になる。この場合、撮影画像の取得位置における波面は、表面１０ａから−Δｚ１だけ離れた位置における波面になる。

式（５）における波面ｕ^ｍ（ｒ）は、波面ｇ_ｏｕｔ ^ｍ（ｒ）が、光の進行方向とは逆の方向にΔｚ１だけ伝搬した波面である。この波面は、表面１０ａから−Δｚ１だけ離れた位置における波面である。よって、式（５）における波面ｕ^ｍ（ｒ）は、撮影画像の取得位置における波面を表している。

なお、撮影画像の取得位置と表面６ａの位置は、厳密には異なる。しかしながら、標本６は薄い標本なので、Δｚ１の値は非常に小さい。そのため、撮影画像の取得位置と表面６ａの位置は、ほぼ同じと見なすことができる。

推定標本１０も薄い標本である。そのため、表面１０ａの位置と、表面１０ａから−Δｚ１だけ離れた位置は、ほぼ同じと見なすことができる。すなわち、波面ｇ_ｏｕｔ ^ｍ（ｒ）の位置と波面ｕ^ｍ（ｒ）の位置は、ほぼ同じと見なすことができる。この場合、波面ｕ^ｍ（ｒ）の代わりに、波面ｇ_ｏｕｔ ^ｍ（ｒ）を用いることもできる。

ステップＳ４４では、結像面における波面ｕ_ｉｍｇ ^ｍ（ｒ）が算出される。波面ｕ^ｍ（ｒ）は、結像面ＩＰに伝搬される。このとき、結像光学系３を通過する。結像光学系３は、フーリエ光学系を形成している。よって、式（７）に示すように、結像面ＩＰにおける波面ｕ_ｉｍｇ ^ｍ（ｒ）は、波面ｕ^ｍ（ｒ）と結像光学系の瞳関数Ｐ（ｕ）を用いて算出することができる。

ステップＳ４５では、波面ｕ_ｉｍｇ ^ｍ（ｒ）が２乗される。波面ｕ_ｉｍｇ ^ｍ（ｒ）は、光の振幅を表している。よって、波面ｕ_ｉｍｇ ^ｍ（ｒ）を２乗することで、光強度が算出される。

｜ｕ_ｉｍｇ ^ｍ（ｒ）｜^２は、結像面ＩＰに光強度分布を表している。第１強度分布を、結像光学系の結像位置における光強度分布とすると、｜ｕ_ｉｍｇ ^ｍ（ｒ）｜^２は、結像光学系の結像位置における第１光強度分布を表している。

波面ｆ_ｉｎ ^ｍ（ｒ）、波面ｇ_ｏｕｔ ^ｍ（ｒ）、波面ｕ^ｍ（ｒ）、及び波面ｕ_ｉｍｇ ^ｍ（ｒ）は、ｍ番目の光源から射出された照明光、すなわち、１つの光源から射出された照明光によって生成される波面を表している。

推定標本の像Ｉ_ｅｓｔ（ｒ）は、全ての光源から射出された照明光によって生成される。よって、全ての光源について、波面ｆ_ｉｎ ^ｍ（ｒ）、波面ｇ_ｏｕｔ ^ｍ（ｒ）、波面ｕ^ｍ（ｒ）、及び波面ｕ_ｉｍｇ ^ｍ（ｒ）を求めなくてはならない。

ステップＳ４６では、変数ｍの値が光源の数Ｎ_ＬＳと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ４７が実行される。判断結果がＹＥＳの場合は、ステップＳ５０が実行される。

（判断結果がＮＯの場合：ｍ≠Ｎ_ＬＳ）
判断結果がＮＯの場合、ステップＳ４７で、変数ｍの値に１が加算される。ステップＳ４７が終ると、ステップＳ４１に戻る。

ステップＳ４７で、変数ｍの値が１つ増えている。そのため、別の光源について、ステップＳ４１で波面ｆ_ｉｎ ^ｍ（ｒ）が算出され、ステップＳ４２で波面ｇ_ｏｕｔ ^ｍ（ｒ）が算出され、ステップＳ４３で波面ｕ^ｍ（ｒ）が算出され、ステップＳ４４で波面ｕ_ｉｍｇ ^ｍ（ｒ）が算出され、ステップＳ４５で｜ｕ_ｉｍｇ ^ｍ（ｒ）｜^２が算出される。

ステップＳ４１、ステップＳ４２、ステップＳ４３、ステップＳ４４、及びステップＳ４５は、全ての光源について｜ｕ_ｉｍｇ ^ｍ（ｒ）｜^２が求まるまで、繰り返し行われる。

（判断結果がＹＥＳの場合：ｍ＝Ｎ_ＬＳ）
判断結果がＹＥＳの場合、ステップＳ５０で、｜ｕ_ｉｍｇ ^ｍ（ｒ）｜^２の足し合わせが行われる。その結果、推定標本の像Ｉ_ｅｓｔ（ｒ）が算出される。推定標本の像Ｉ_ｅｓｔ（ｒ）は、式（８）で表される。

図６は、推定標本の像を示す図である。推定標本の像Ｉ_ｅｓｔ（ｒ）は、全ての光源について波面ｕ_ｉｍｇ ^ｍ（ｒ）が求められた場合の像である。図６に示すように、各光源について波面ｕ_ｉｍｇ ^ｍ（ｒ）が算出され、波面ｕ_ｉｍｇ ^ｍ（ｒ）から｜ｕ_ｉｍｇ ^ｍ（ｒ）｜^２が算出され｜ｕ_ｉｍｇ ^ｍ（ｒ）｜^２が全て足し合わされている。その結果、推定標本の像Ｉ_ｅｓｔ（ｒ）が算出される。

ステップＳ６０で、残差が算出される。残差は、式（９）で表される。式（９）に示すように、残差は、撮影画像Ｉ_ｍｅａ（ｒ）と推定標本の像Ｉ_ｅｓｔ（ｒ）とから算出される。

式（９）は、行列のノルムを表している。ノルムは、式（１０）で表される。

ステップＳ７０では、残差と閾値との比較が行われる。判断結果がＮＯの場合は、ステップＳ８０が実行される。判断結果がＹＥＳの場合は、ステップＳ１１０が実行される。

（判断結果がＮＯの場合：残差≧閾値）

ステップＳ８０では、変数ｍの値が初期化される。後述のステップＳ９１とステップＳ９２は、全ての光源に対して実行される。変数ｍは、これらのステップが実行された回数を表している。

ステップＳ９０は、推定標本の屈折率分布を最適化するステップである。

ステップＳ９０は、ステップＳ９１、ステップＳ９２、ステップＳ９３、及びステップＳ９４を有する。

ステップＳ９１では、波面ｕ’^ｍ（ｒ）が算出される。波面ｕ’^ｍ（ｒ）の算出では、撮影画像Ｉ_ｍｅａ（ｒ）と推定標本の像Ｉ_ｅｓｔ（ｒ）が用いられる。また、波面ｕ’^ｍ（ｒ）は、撮影画像の取得位置における波面である。

波面ｕ’^ｍ（ｒ）は式（１１）で表される。

図７は、波面の補正を示す図である。図７（ａ）は、推定標本から射出する補正前の波面を示す図である。図７（ｂ）は、撮影画像の取得位置における補正前の波面を示す図である。図７（ｃ）は、撮影画像の取得位置における補正後の波面を示す図である。図７（ｄ）は、推定標本から射出する補正後の波面を示す図である。

図６に示すように、推定標本の像Ｉ_ｅｓｔ（ｒ）は、波面ｕ_ｉｍｇ ^ｍ（ｒ）に基づいて算出される。また、図６と図７（ｂ）に示すように、波面ｕ_ｉｍｇ ^ｍ（ｒ）は、波面ｕ^ｍ（ｒ）に基づいて算出される。

図７（ａ）に示すように、波面ｕ^ｍ（ｒ）の算出には、振幅透過率Ｔ_ｓ（ｒ）が用いられている。振幅透過率Ｔ_ｓ（ｒ）は、推定された振幅透過率である。ステップＳ９０の１回目の実行時、この振幅透過率Ｔ_ｓ（ｒ）は、標本６の振幅透過率と異なる。

振幅透過率Ｔ_ｓ（ｒ）と標本６の振幅透過率との差が大きくなるほど、推定標本の像Ｉ_ｅｓｔ（ｒ）と撮影画像Ｉ_ｍｅａ（ｒ）との差も大きくなる。よって、推定標本の像Ｉ_ｅｓｔ（ｒ）と撮影画像Ｉ_ｍｅａ（ｒ）との差は、振幅透過率Ｔ_ｓ（ｒ）と標本６の振幅透過率との差を反映していると見なすことができる。

そこで、式（１１）に示すように、推定標本の像Ｉ_ｅｓｔ（ｒ）と撮影画像Ｉ_ｍｅａ（ｒ）と用いて、波面ｕ^ｍ（ｒ）を補正する。その結果、図７（ｃ）に示すように、補正後の波面、すなわち、波面ｕ’^ｍ（ｒ）が得られる。

波面ｕ’^ｍ（ｒ）を用いることで、新たな振幅透過率Ｔ_ｓ（ｒ）を算出できる。波面ｕ’^ｍ（ｒ）は、波面ｕ^ｍ（ｒ）と異なる。よって、新たな振幅透過率Ｔ_ｓ（ｒ）は、波面ｕ^ｍ（ｒ）を算出したときの振幅透過率とは異なる。

このように、波面ｕ’^ｍ（ｒ）を用いて、振幅透過率Ｔ_ｓ（ｒ）を算出することができる。ただし、図７（ａ）に示すように、振幅透過率Ｔ_ｓ（ｒ）の算出には、波面ｇ_ｏｕｔ ^ｍ（ｒ）が必要である。

図７（ａ）、（ｃ）に示すように、波面ｕ’^ｍ（ｒ）の位置と、波面ｇ_ｏｕｔ ^ｍ（ｒ）の位置は異なる。よって、振幅透過率Ｔ_ｓ（ｒ）を算出するためには、図７（ｄ）に示すように、波面ｇ’_ｏｕｔ ^ｍ（ｒ）が必要である。

波面ｇ’_ｏｕｔ ^ｍ（ｒ）は、式（１２）で表される。波面ｕ’^ｍ（ｒ）は補正後の波面なので、波面ｇ’_ｏｕｔ ^ｍ（ｒ）も補正後の波面である。

上述のように、撮影画像の取得位置は、表面１０ａから−Δｚ１だけ離れた位置になる。言い換えると、表面１０ａの位置は、撮影画像の取得位置からΔｚ１だけ離れた位置になる。よって、表面１０ａの位置における波面は、撮影画像の取得位置からΔｚ１だけ離れた位置における波面になる。

式（１２）における波面ｇ’_ｏｕｔ ^ｍ（ｒ）は、波面ｕ’^ｍ（ｒ）が、光の進行方向にΔｚ１だけ伝搬した波面である。この波面は、影画像の取得位置からΔｚ１だけ離れた位置における波面である。よって、式（１２）における波面ｇ’_ｏｕｔ ^ｍ（ｒ）は、表面１０ａの位置における波面を表している。

表面１０ａの位置における波面は、ｆ_ｉｎ ^ｍ（ｒ）が推定標本１０を通過した後の波面である。上述のように、ｆ_ｉｎ ^ｍ（ｒ）は、第１波面を表している。第２波面を、第１波面が推定標本を通過した後の波面とすると、波面ｇ’_ｏｕｔ ^ｍ（ｒ）は、第２波面を表している。

上述のように、Δｚ１の値は非常に小さい。また、推定標本１０は薄い標本である。そのため、撮影画像の取得位置と、撮影画像の取得位置からΔｚ１だけ離れた位置は、ほぼ同じと見なすことができる。すなわち、波面ｕ’^ｍ（ｒ）の位置と、波面ｇ_ｏｕｔ ^ｍ（ｒ）の位置は、ほぼ同じと見なすことができる。この場合、波面ｇ’_ｏｕｔ ^ｍ（ｒ）の代わりに、波面ｕ’^ｍ（ｒ）を用いることもできる。

ステップＳ９２では、標本の勾配ΔＴ_ｓ ^ｍ（ｒ）が算出される。標本の勾配ΔＴ_ｓ ^ｍは式（１３）で表される。標本の勾配ΔＴ_ｓ ^ｍ（ｒ）の算出には、例えば、勾配降下法を用いることができる。

ここで、
ｆ^＊は、ｆの複素共役、
δは、ゼロ除算を防ぐための正規化定数、
である。

図７（ａ）に示すように、波面ｇ_ｏｕｔ ^ｍ（ｒ）の算出には、振幅透過率Ｔ_ｓ（ｒ）が用いられている。振幅透過率Ｔ_ｓ（ｒ）は、推定された振幅透過率である。よって、この振幅透過率Ｔ_ｓ（ｒ）は、標本６の振幅透過率と異なる。

振幅透過率Ｔ_ｓ（ｒ）と標本６の振幅透過率との差が大きくなるほど、波面ｇ_ｏｕｔ ^ｍ（ｒ）と波面ｇ’_ｏｕｔ ^ｍ（ｒ）との差も大きくなる。よって、波面ｇ_ｏｕｔ ^ｍ（ｒ）と波面ｇ’_ｏｕｔ ^ｍ（ｒ）との差は、振幅透過率Ｔ_ｓ（ｒ）と標本６の振幅透過率との差を反映していると見なすことができる。

波面ｆ_ｉｎ ^ｍ（ｒ）、振幅透過率Ｔ_ｓ（ｒ）、波面ｇ_ｏｕｔ ^ｍ（ｒ）、及び波面ｇ’_ｏｕｔ ^ｍ（ｒ）は、既知である。そこで、式（１２）に示すように、波面ｆ_ｉｎ ^ｍ（ｒ）、振幅透過率Ｔ_ｓ（ｒ）、波面ｇ_ｏｕｔ ^ｍ（ｒ）、及び波面ｇ’_ｏｕｔ ^ｍ（ｒ）用いて、標本の勾配ΔＴ_ｓ ^ｍ（ｒ）を算出することができる。

図８は、標本の勾配を示す図である。

ステップＳ９２で得られる標本の勾配ΔＴ_ｓ ^ｍ（ｒ）は、１つの光源から射出された照明光における標本の勾配を表している。標本の勾配ΔＴ_ｓ ^ｍ（ｒ）は、全ての光源から射出された照明光によって決まる。よって、全ての光源について、標本の勾配ΔＴ_ｓ ^ｍ（ｒ）を求めなくてはならない。

ステップＳ９３では、変数ｍの値が光源の数Ｎ_ＬＳと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ９４が実行される。判断結果がＹＥＳの場合は、ステップＳ１００が実行される。

（判断結果がＮＯの場合：ｍ≠Ｎ_ＬＳ）
判断結果がＮＯの場合、ステップＳ９４で、変数ｍの値に１が加算される。ステップＳ９４が終ると、ステップＳ９１に戻る。

ステップＳ９４で、変数ｍの値が１つ増えている。そのため、別の光源について、ステップＳ９１で波面ｕ’^ｍ（ｒ）が算出され、ステップＳ９２で標本の勾配ΔＴ_ｓ ^ｍ（ｒ）が算出される。

ステップＳ９１とステップＳ９２は、全ての光源について標本の勾配ΔＴ_ｓ ^ｍ（ｒ）が求まるまで、繰り返し行われる。

図９は、標本の勾配を示す図である。図９では、全ての光源について、標本の勾配ΔＴ_ｓ ^ｍ（ｒ）が求められている。

（判断結果がＹＥＳの場合：ｍ＝Ｎ_ＬＳ）
判断結果がＹＥＳの場合、ステップＳ１００で、振幅透過率Ｔ_ｓ（ｒ）が更新される。ステップＳ１００は、推定標本を更新するステップするステップである。

更新された振幅透過率Ｔ_ｓ（ｒ）は、式（１４）で表される。

ここで、
αは、標本の勾配の補正係数、
である。

さらに、標本６を吸収の無い完全位相物体として考える場合、式（１５）を用いて、振幅透過率Ｔ_ｓ（ｒ）をさらに更新することができる。

ステップＳ１００が終ると、ステップＳ３０に戻る。更新された振幅透過率Ｔ_ｓ（ｒ）で、ステップＳ３０からステップＳ１００までが実行される。

ステップＳ３０からステップＳ１００までが繰り返し実行されることで、更新された振幅透過率Ｔ_ｓ（ｒ）は、徐々に、標本６の振幅透過率に近づく。すなわち、残差が小さくなる。やがて、残差は閾値よりも小さくなる。

（判断結果がＹＥＳの場合：残差＜閾値）
ステップ１１０では、推定標本の屈折率分布が算出される。得られた振幅透過率Ｔ_ｓ（ｒ）は、標本６の振幅透過率と同一か、又は、略同一である。得られた振幅透過率Ｔ_ｓ（ｒ）と式（１）から、屈折率分布ｎ（ｒ）が求まる。

ステップＳ１１０で得られた屈折率分布ｎ（ｒ）を用いることで、推定標本の構造を再構成することができる。再構成された推定標本の構造は、例えば、表示装置に出力することができる。推定標本１０は、薄い標本である。第１のシミュレーションでは、薄い標本の構造を再構成することができる。

上述のように、ステップＳ１１０で得られた振幅透過率Ｔ_ｓ（ｒ）は、標本６の振幅透過率と同一か、又は、略同一である。この場合、屈折率分布ｎ（ｒ）も、標本６の屈折率分布と同一か、又は、略同一と見なすことができる。よって、再構成された推定標本１０の構造は、標本６の構造と同一か、又は、略同一と見なすことができる。

第１のシミュレーションでは、ステップＳ４０、ステップＳ５０、及びステップＳ９０が、繰り返し実行される。その結果、振幅透過率Ｔ_ｓ（ｒ）が更新される。上述のように、ステップＳ４０とステップＳ５０は、推定標本の像を算出するステップである。ステップＳ９０は、推定標本の屈折率分布を最適化するステップである。

振幅透過率Ｔ_ｓ（ｒ）は、推定標本を表している。よって、推定標本の像を算出するステップと推定標本の屈折率分布を最適化するステップが繰り返し実行されることで、推定標本が更新される。

本実施形態の屈折率分布推定システムは、結像光学系の焦点位置と標本の位置との間隔を、結像光学系の光軸方向に変化させる駆動機構を更に備え、駆動機構により間隔を変化させて、複数の間隔のそれぞれに対応する複数の標本の撮影画像を取得し、複数の間隔のそれぞれに対応する複数の推定標本の像を算出し、複数の間隔のそれぞれにおいて、推定標本の屈折率分布を最適化することが好ましい。

図１０は、本実施形態の屈折率分布推定システムと照明光を示す図である。図１０（ａ）は、屈折率分布推定システムを示す図である。図１０（ｂ）は、照明光を示す図である。図１（ａ）と同じ構成については同じ番号を付し、説明は省略する。

標本２１は、厚い標本である。標本２１には、複数の方向から、光線が同時に入射する。図１０（ａ）では、光線Ｌ_ＯＮと光線Ｌ_ＯＦＦが示されている。図１０（ｂ）に示すように、光線Ｌ_ＯＦＦだけを照明に用いても良い。

標本２１から射出した光は、結像光学系３によって、結像面ＩＰに集光される。結像面ＩＰに、光学像２１’が形成される。光学像２１’は、標本２１の光学像である。

屈折率分布推定システム２０は、移動ステージ２２を有する。移動ステージ２２は、光軸ＡＸの方向に移動する。

上述のように、推定標本の屈折率分布の最適化には、撮影画像が用いられる。標本２１は、厚い標本なので、複数の撮像画像が取得される。複数の撮像画像を取得するために、標本２１を固定し、移動ステージ２２で、結像光学系３の焦点位置を移動させている。

結像光学系３は、例えば、無限遠補正対物レンズと結像レンズを有する。この場合、対物レンズを移動させることで、結像光学系３の焦点位置を移動させることができる。結像光学系３と撮像素子４を固定し、標本２１を移動させても良い。

以下、取得された撮像画像が４枚の場合について説明する。

図１１は、撮影画像を示す図である。図１１（ａ）は、第１の位置における撮影画像を示す図である。図１１（ｂ）は、第２の位置における撮影画像を示す図である。図１１（ｃ）は、第３の位置における撮影画像を示す図である。図１１（ｄ）は、第４の位置における撮影画像を示す図である。

結像光学系３と標本２１との間隔を変化させることで、標本２１に対する焦点位置Ｆｏが変化する。ここでは、標本２１に対する焦点位置Ｆｏを、４回変化させている。これにより、以下の４つの撮影画像が取得される。
撮影画像Ｉ_ｍｅａ１（ｒ）：表面２１ａからの距離が３×Δｚの位置の画像。
撮影画像Ｉ_ｍｅａ２（ｒ）：表面２１ａからの距離が２×Δｚの位置の画像。
撮影画像Ｉ_ｍｅａ３（ｒ）：表面２１ａからの距離がΔｚの位置の画像。
撮影画像Ｉ_ｍｅａ４（ｒ）：表面２１ａの画像。

撮影画像Ｉ_ｍｅａ１（ｒ）、撮影画像Ｉ_ｍｅａ２（ｒ）、撮影画像Ｉ_ｍｅａ３（ｒ）、及び撮影画像Ｉ_ｍｅａ４（ｒ）は、プロセッサ５に入力される。プロセッサ５では、４つの撮影画像を用いて、推定標本の像のシミュレーションが行われる。

図１２と図１３は、第２のシミュレーションのフローチャートである。第１のフローチャートにおける処置と同じ処理につては、同じ番号を付し、説明は省略する。

図１４は、シミュレーションで用いる光学系を示す図である。図５と同じ構成については同じ番号を付し、説明は省略する。

シミュレーションで用いる光学系は、撮影画像Ｉ_ｍｅａ１（ｒ）、撮影画像Ｉ_ｍｅａ２（ｒ）、撮影画像Ｉ_ｍｅａ３（ｒ）、及び撮影画像Ｉ_ｍｅａ４（ｒ）を取得した測定光学系と同一である。シミュレーションでは、標本２１の代わりに、推定標本３０が用いられる。

図１４には、推定標本３０、波面ｆ_ｉｎ ^ｍ（ｒ）、振幅透過率Ｔ_ｚ（ｒ）、及び波面ｇ_ｏｕｔ ^ｍ（ｒ）が図示されている。

推定標本が薄い標本の場合、式（４）に示すように、波面ｆ_ｉｎ ^ｍ（ｒ）から、波面ｇ_ｏｕｔ ^ｍ（ｒ）を直接算出することができる。しかしながら、推定標本が厚い標本の場合、波面ｆ_ｉｎ ^ｍ（ｒ）から波面ｇ_ｏｕｔ ^ｍ（ｒ）を直接算出することは困難である。

推定標本３０は、厚い標本である。そこで、光軸方向に沿って、推定標本３０を複数の薄い層に置き換える。そして、薄い層の各々について、層の両側の波面を算出する。

図１５は、各層における波面を示す図である。波面の算出については後述する。層の数は、取得画像の数と同じにすることができる。ただし、層の数は、取得画像の数よりも多くしても良い。図１５では、層の数は、取得画像の数と同じである。

図１５では、Ｚ＝１の位置が第１層の位置、Ｚ＝２の位置が第２層の位置、Ｚ＝３の位置が第３層の位置、Ｚ＝４の位置が第４層の位置である。

図１２、図１３を用いて、シミュレーションについて説明する。

ステップＳ１０では、光源の数Ｎ_ＬＳが設定される。

ステップＳ２００では、層の数Ｎ_ＩＭが設定される。推定標本３０は、厚い標本である。よって、上述のように、推定標本３０を複数の薄い層に置き換える。層の数Ｎ_ＩＭは、薄い層の数を表している。

標本２１では、複数の位置で撮影画像が取得される。層の数Ｎ_ＩＭは、撮影画像が取得された位置の数と同じにすることができる。標本２１に対する焦点位置Ｆｏを４回変化させている場合、Ｎ_ＩＭ＝４になる。

１からＮ_ＩＭまでの数字は、薄い層の位置を表している。例えば、Ｎ_ＩＭ＝４の場合、数字の１は第１層の位置を表し、数字の２は第２層の位置を表し、数字の３は第３層の位置を表し、数字の４は第４層の位置を表している。

推定標本の像の算出は、シミュレーションで行われる。そのため、層の数Ｎ_ＩＭは、自由に設定できる。例えば、層の数Ｎ_ＩＭは、撮影画像が取得された位置の数よりも多くすることができる。

例えば、Ｎ_ＩＭ＝７の場合、薄い層の数は７つである。この場合、７つの推定標本の像が算出される。シミュレーションでは、後述のように、撮影画像と薄い層における推定標本の像が用いられる。よって、推定標本の像が算出される７つ位置には、撮影画像が取得された４つの位置が含まれている。

７つ位置と撮影画像の関係は、例えば、以下のようにすることができる。

数字の１は、第１層の位置を表している。この位置では、撮影画像Ｉ_ｍｅａ１（ｒ）が取得されている。また、この位置では、第１層における推定標本の像が算出される。よって、第１層における推定標本の像と撮影画像Ｉ_ｍｅａ１（ｒ）が、後述のステップで用いられる。

数字の２は、第２層の位置を表している。この位置で取得された撮影画像は無い。

数字の３は、第３層の位置を表している。この位置では、撮影画像Ｉ_ｍｅａ２（ｒ）が取得されている。また、この位置では、第３層における推定標本の像が算出される。よって、第３層における推定標本の像と撮影画像Ｉ_ｍｅａ２（ｒ）が、後述のステップで用いられる。

数字の４は、第４層の位置を表している。この位置で取得された撮影画像は無い。

数字の５は、第５層の位置を表している。この位置では、撮影画像Ｉ_ｍｅａ３（ｒ）が取得されている。また、この位置では、第５層における推定標本の像が算出される。よって、第５層における推定標本の像と撮影画像Ｉ_ｍｅａ３（ｒ）が、後述のステップで用いられる。

数字の６は、第６層の位置を表している。この位置で取得された撮影画像は無い。

数字の７は、第７層の位置を表している。この位置では、撮影画像Ｉ_ｍｅａ４（ｒ）が取得されている。また、この位置では、第７層における推定標本の像が算出される。よって、第７層における推定標本の像と撮影画像Ｉ_ｍｅａ４（ｒ）が、後述のステップで用いられる。

ステップＳ２１０では、補正の回数Ｎ_ＣＲが設定される。

ステップＳ２２０では、変数ｚの値が初期化される。後述のステップＳ２３１は、全ての取得位置に対して実行される。変数ｚは、ステップＳ２３１が実行された回数を表している。

ステップＳ２３０は、推定標本を推定するステップである。標本２１では、４つの撮影画像が取得されている。上述のように、推定標本３０は４つの薄い層に置き換えられている。よって、振幅透過率の初期値の設定は４回行われる。

ステップＳ２３０は、ステップＳ２３１、ステップＳ２３２、及びステップＳ２３３を有する。

ステップＳ２３１では、推定標本３０における振幅透過率Ｔ_ｚ（ｒ）に初期値が設定される。

初期値の設定では、強度輸送方程式を用いても良い。強度輸送方程式は、例えば、以下の文献に開示されている。
M. R. Teague, “Deterministic phase retrieval: a Greens function solution,” J. Opt. Soc. Am. 73, 1434-1441（1983）

焦点位置Ｚ０での強度輸送方程式は、式（１６）で表される。

ここで、
∇^２は、２次のラプラシアン、
ｋは、波数、
φ_Ｚ０（ｒ）は，結像面での標本の位相分布、
Ｉ_Ｚ０は，光学像の平均光強度、
δＩ_{ｍｅａＺ０}（ｒ）／δ_Ｚは、結像面から±△zだけ離れた２つのデフォーカス像の差分像、
である。

式（１６）を用いると、フォーカス像と２つのデフォーカス像から、標本の位相分布φ_Ｚ０（ｒ）を簡単に求めることができる。

ただし、２つのデフォーカス像の同一点における光強度の差が０か、又は、非常に小さいと、位相を計測できない。パーシャルコヒーレント照明であっても、照明光の開口数が小さい場合は、この光強度の差が０になるか、又は、非常に小さくなる。そのため、このような場合は、強度輸送方程式を用いて初期値を設定することは難しい。

上述のように、位相分布φ_Ｚ０（ｒ）は、フォーカス像と２つのデフォーカス像から算出される。フォーカス像、例えば、対物レンズを一定の間隔で、光軸方向に移動させることで取得される。この場合、光軸に沿って、複数のフォーカス像が離散的に取得される。よって、２つのデフォーカス像も離散的に取得される。

式（１６）で表される位相分布φ_Ｚ０（ｒ）は、光軸と直交する面内における位相分布である。フォーカス像と２つのデフォーカス像は離散的に取得されるので、位相分布φ_Ｚ０（ｒ）を表す面も、光軸に沿って離散的に位置している。

式（１７）に示すように、位相分布φ_ｚ（ｒ）は、振幅透過率Ｔ_ｓ（ｒ）に変換することができる。このようにして、振幅透過率Ｔ_ｚ（ｒ）に初期値を設定することができる。

強度輸送方程式で得られた位相分布φ_Ｚ０は、位相分布φ_ｚ（ｒ）に用いることができる。強度輸送方程式を用いて、初期値を設定することができる。なお、初期値の推定が難しい場合、例えば、Ｔ_ｚ（ｒ）＝１としても構わない。

ステップＳ２３２では、変数ｚの値が取得位置の数Ｎ_ＩＭと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ２３３が実行される。判断結果がＹＥＳの場合は、ステップＳ３０が実行される。

（判断結果がＮＯの場合：ｚ≠Ｎ_ＩＭ）
判断結果がＮＯの場合、ステップＳ２３３で、変数ｚの値に１が加算される。ステップＳ２３３が終ると、ステップＳ２３１に戻る。

ステップＳ２３３で、変数ｚの値が１つ増えている。そのため、別の取得位置について、ステップＳ２３１で振幅透過率Ｔ_ｚ（ｒ）に初期値が設定される。

ステップＳ２３１は、全ての取得位置について初期値が設定されるまで、繰り返し行われる。

（判断結果がＹＥＳの場合：ｚ＝Ｎ_ＩＭ）
ステップＳ３０で、変数ｍの値が初期化される。後述のステップＳ２４０、ステップＳ４１、ステップＳ４２、ステップＳ２５１、及びステップＳ２６０は、全ての光源に対して実行される。変数ｍは、これらのステップが実行された回数を表している。

ステップＳ２４０で、関数Ｉ_ｅｓｔｚ（ｒ）の値が初期化される。Ｉ_ｅｓｔｚ（ｒ）は、推定標本３０の像を表している。上述のように、推定標本３０の像は、４つの薄い層に置き換えられている。よって、Ｉ_ｅｓｔｚ（ｒ）は、薄い層の像を表している。

ステップＳ２５０とステップＳ２７０は、推定標本の像を算出するステップである。推定標本の像の数は、撮影画像の数と同数である。撮影画像の数は４なので、推定標本の像の数も４である。

ステップＳ２５０は、ステップＳ４１、ステップＳ４２、ステップＳ２５１、ステップＳ２５２、ステップＳ２５３、及びステップＳ２６０を有する。

ステップＳ４１では、推定標本３０へ入射する波面ｆ_ｉｎ ^ｍ（ｒ）が算出される。波面ｆ_ｉｎ ^ｍ（ｒ）は、上述の式（２）、（３）で表される。

ステップＳ４２では、推定標本３０から射出する波面ｇ_ｏｕｔ ^ｍ（ｒ）が算出される。波面ｇ_ｏｕｔ ^ｍ（ｒ）は、波面ｆ_ｉｎ ^ｍ（ｒ）に基づいて算出される。推定標本３０は４つの薄い層に置き換えられている。よって、薄い層の各々で波面を算出する。

図１５では、Ｚ＝１の位置が第１層の位置、Ｚ＝２の位置が第２層の位置、Ｚ＝３の位置が第３層の位置、Ｚ＝４の位置が第４層の位置である。

４つの薄い層は、等間隔で並んでいる。隣り合う２つの層の間隔はΔｚである。波面は、２つの層の間を伝搬する。よって、Δｚは、伝搬距離を表している。

第１層における波面ｆ_１ ^ｍ（ｒ）は、式（１８）と式（３）で表される。

第１層の位置は、推定標本３０の表面３０ｂの位置と一致している。表面３０ｂには、波面ｆ_ｉｎ ^ｍ（ｒ）が入射する。よって、波面ｆ_１ ^ｍ（ｒ）は、波面ｆ_ｉｎ ^ｍ（ｒ）を表している。図１５では、波面ｆ_１ ^ｍ（ｒ）の代わりに、波面ｆ_ｉｎ ^ｍ（ｒ）が示されている。

第１層における波面ｇ_１ ^ｍ（ｒ）は、式（１９）で表される。

ここで、
Ｔ_１（ｒ）は、第１層における振幅透過率、
である。

第２層における波面ｆ_２ ^ｍ（ｒ）は、波面ｇ_１ ^ｍ（ｒ）が、Δｚだけ伝搬したときの波面である。波面ｆ_２ ^ｍ（ｒ）は、式（２０）で表される。式（２０）で、ΔＤ＝Δｚとすることで、波面ｆ_２ ^ｍ（ｒ）を算出することができる。

ここで、Ｐ_ΔＤ｛ ｝は、式（２１）で表され、

ΔＤは、隣り合う２つの層の間隔、
λは、波長、
ｕは、瞳面座標（ξ，η）の２次元表記、
Ｆ_２Ｄは、２次元フーリエ変換、
Ｆ_２Ｄ ^−１は、２次元フーリエ逆変換、
である。

第２層における波面ｇ_２ ^ｍ（ｒ）は、式（２２）で表される。

ここで、
Ｔ_２（ｒ）は、第２層における振幅透過率、
である。

第３層における波面ｆ_３ ^ｍ（ｒ）は、波面ｇ_２ ^ｍ（ｒ）が、Δｚだけ伝搬したときの波面である。第３層における波面ｆ_３ ^ｍ（ｒ）は、式（２３）で表される。式（２１）で、ΔＤ＝Δｚとすることで、波面ｆ_３ ^ｍ（ｒ）を算出することができる。

第３層における波面ｇ_３ ^ｍ（ｒ）は、式（２４）で表される。

ここで、
Ｔ_３（ｒ）は、第３層における振幅透過率、
である。

第４層における波面ｆ_４ ^ｍ（ｒ）は、波面ｇ_３ ^ｍ（ｒ）が、Δｚだけ伝搬したときの波面である。第４層における波面ｆ_４ ^ｍ（ｒ）は、式（２５）で表される。式（２１）で、ΔＤ＝Δｚとすることで、波面ｆ_４ ^ｍ（ｒ）を算出することができる。

第４層における波面ｇ_４ ^ｍ（ｒ）は、式（２６）で表される。

ここで、
Ｔ_４（ｒ）は、第４層における振幅透過率、
である。

第４層の位置は、推定標本３０の表面３０ａの位置と一致している。表面３０ａからは、波面ｇ_ｏｕｔ ^ｍ（ｒ）が射出される。よって、波面ｇ_４ ^ｍ（ｒ）は、波面ｇ_ｏｕｔ ^ｍ（ｒ）を表している。図１５では、波面ｇ_４ ^ｍ（ｒ）の代わりに、波面ｇ_ｏｕｔ ^ｍ（ｒ）が示されている。

以上のように、推定標本が厚い標本の場合は、複数の薄い層に置き換えると共に、２つの層の間を伝搬する波面を求めることで、波面ｇ_ｏｕｔ ^ｍ（ｒ）を算出することができる。

ステップＳ２５１では、変数ｚの値が初期化される。後述のステップＳ２６１、ステップＳ２６２、及びステップＳ２６３は、全ての取得位置に対して実行される。変数ｚは、これらのステップが実行された回数を表している。

ステップＳ２６０は、ステップＳ２６１、ステップＳ２６２、ステップＳ２６３、ステップＳ２６４、及びステップＳ２６５を有する。

ステップＳ２６１では、撮影画像の取得位置における波面ｕ_ｚ ^ｍ（ｒ）が算出される。波面ｕ_ｚ ^ｍ（ｒ）は、式（２７）で表される。

ここで、
△Ｄは推定標本の表面から薄い層までの距離、
である。
Ｐ_ΔＤ｛ ｝は、式（２８）で表される。

ステップＳ２６２では、結像面における波面ｕ_ｉｍｇｚ ^ｍ（ｒ）が算出される。波面ｕ_ｉｍｇｚ ^ｍ（ｒ）は、式（２９）で表される。

ステップＳ２６３では、波面ｕ_ｉｍｇｚ ^ｍ（ｒ）が２乗される。波面ｕ_ｉｍｇｚ ^ｍ（ｒ）は、光の振幅を表している。よって、波面ｕ_ｉｍｇｚ ^ｍ（ｒ）を２乗することで、光強度が算出される。

｜ｕ_ｉｍｇｚ ^ｍ（ｒ）｜^２は、結像面ＩＰに光強度分布を表している。第１強度分布を、結像光学系の結像位置における光強度分布とすると、｜ｕ_ｉｍｇｚ ^ｍ（ｒ）｜^２は、結像光学系の結像位置における第１光強度分布を表している。

ステップＳ２６４では、変数ｚの値が取得位置の数Ｎ_ＩＭと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ２６５が実行される。判断結果がＹＥＳの場合は、ステップＳ２５２が実行される。

（判断結果がＮＯの場合：ｚ≠Ｎ_ＩＭ）
判断結果がＮＯの場合、ステップＳ２６５で、変数ｚの値に１が加算される。ステップＳ２６５が終ると、ステップＳ２６１に戻る。

ステップＳ２６５で、変数ｚの値が１つ増えている。そのため、別の取得位置について、ステップＳ２６１、ステップＳ２６２、及びステップＳ２６３が実行される。

ステップＳ２６１、ステップＳ２６２、及びステップＳ２６３は、全ての取得位置について初期値が設定されるまで、繰り返し行われる。

ステップＳ２５０における処理を、第１層と第４層を用いて説明する。第２層と第３層については、第１層と同じように考えれば良い。

図１６は、撮影画像の取得位置における波面と結像面における波面を示す図である。図１６（ａ）は、第１層におけるにおける２つの波面を示す図である。図１６（ｂ）は、第４層における２つの波面を示す図である。

ｚ＝１における撮影画像は、撮影画像Ｉ_ｍｅａ１（ｒ）である。撮影画像Ｉ_ｍｅａ１（ｒ）は、表面２１ａからの距離が３×Δｚの位置の画像である。第１層は、表面３０ａから３×Δｚだけ離れている。よって、第１層の位置が、撮影画像Ｉ_ｍｅａ１（ｒ）の取得位置に対応する。

波面ｇ_ｏｕｔ ^ｍ（ｒ）の射出位置は、表面３０ａと一致している。図１６（ａ）に示すように、波面ｇ_ｏｕｔ ^ｍ（ｒ）の射出位置は、第１層の位置と異なる。第１層は、波面ｇ_ｏｕｔ ^ｍ（ｒ）の射出位置から３×Δｚだけ離れている。

第１層における波面ｕ_１ ^ｍ（ｒ）は、波面ｇ_ｏｕｔ ^ｍ（ｒ）が光の進行方向とは反対に３×Δｚだけ伝搬したときの波面である。よって、ステップＳ２６１においてΔＤ＝−３×Δｚとすることで、式（２７）と式（２８）から、波面ｕ_１ ^ｍ（ｒ）を算出することができる。

波面ｕ_１ ^ｍ（ｒ）が算出されると、ステップＳ２６２で、式（２９）から、結像面における波面ｕ_ｉｍｇ１ ^ｍ（ｒ）が算出される。

更に、ステップＳ２６３で、第１層の像の光強度｜ｕ_ｉｍｇ１（ｒ）｜^２が算出される。

ｚ＝２における撮影画像は、撮影画像Ｉ_ｍｅａ２（ｒ）である。撮影画像Ｉ_ｍｅａ２（ｒ）は、表面２１ａからの距離が２×Δｚの位置の画像である。第２層は、表面３０ａから２×Δｚだけ離れている。よって、第２層の位置が、撮影画像Ｉ_ｍｅａ２（ｒ）の取得位置に対応する。

波面ｇ_ｏｕｔ ^ｍ（ｒ）の射出位置は、第２層の位置と異なる。第２層は、波面ｇ_ｏｕｔ ^ｍ（ｒ）の射出位置から２×Δｚだけ離れている。

第２層における波面ｕ_２ ^ｍ（ｒ）は、波面ｇ_ｏｕｔ ^ｍ（ｒ）が光の進行方向とは反対に２×Δｚだけ伝搬したときの波面である。よって、ステップＳ２６１においてΔＤ＝−２×Δｚとすることで、波面ｕ_２ ^ｍ（ｒ）を算出することができる。

波面ｕ_２ ^ｍ（ｒ）が算出されると、ステップＳ２６２で、結像面における波面ｕ_ｉｍｇ２ ^ｍ（ｒ）が算出される。

更に、ステップＳ２６３で、第２層の像の光強度｜ｕ_ｉｍｇ２（ｒ）｜^２が算出される。

ｚ＝３における撮影画像は、撮影画像Ｉ_ｍｅａ３（ｒ）である。撮影画像Ｉ_ｍｅａ３（ｒ）は、表面２１ａからの距離がΔｚの位置の画像である。第３層は、表面３０ａからΔｚだけ離れている。よって、第３層の位置が、撮影画像Ｉ_ｍｅａ３（ｒ）の取得位置に対応する。

波面ｇ_ｏｕｔ ^ｍ（ｒ）の射出位置は、第３層の位置と異なる。第３層は、波面ｇ_ｏｕｔ ^ｍ（ｒ）の射出位置からΔｚだけ離れている。

第３層における波面ｕ_３ ^ｍ（ｒ）は、波面ｇ_ｏｕｔ ^ｍ（ｒ）が光の進行方向とは反対にΔｚだけ伝搬したときの波面である。よって、ステップＳ２６１においてΔＤ＝Δｚとすることで、波面ｕ_３ ^ｍ（ｒ）を算出することができる。

波面ｕ_３ ^ｍ（ｒ）が算出されると、ステップＳ２６２で、結像面における波面ｕ_ｉｍｇ３ ^ｍ（ｒ）が算出される。

更に、ステップＳ２６３で、第３層の像の光強度｜ｕ_ｉｍｇ３（ｒ）｜^２が算出される。

ｚ＝４における撮影画像は、撮影画像Ｉ_ｍｅａ４（ｒ）である。撮影画像Ｉ_ｍｅａ４（ｒ）は、表面２１ａの画像である。第４層は、表面３０ａと一致している。よって、第４層の位置が、撮影画像Ｉ_ｍｅａ４（ｒ）の取得位置に対応する。

波面ｇ_ｏｕｔ ^ｍ（ｒ）の射出位置は、表面３０ａである。図１６（ｂ）に示すように、波面ｇ_ｏｕｔ ^ｍ（ｒ）の射出位置は、第４層の位置と同じである。

第４層における波面ｕ_４ ^ｍ（ｒ）は、波面ｇ_ｏｕｔ ^ｍ（ｒ）同じである。波面ｇ_ｏｕｔ ^ｍ（ｒ）を波面ｕ_４ ^ｍ（ｒ）に置き換えることができる。

波面ｕ_４ ^ｍ（ｒ）が算出されると、ステップＳ２６２で、結像面における波面ｕ_ｉｍｇ４ ^ｍ（ｒ）が算出される。

更に、ステップＳ２６３で、第４層の像の光強度｜ｕ_ｉｍｇ４（ｒ）｜^２が算出される。

波面ｕ_ｚ ^ｍ（ｒ）と波面ｕ_ｉｍｇｚ ^ｍ（ｒ）は、ｍ番目の光源から射出された照明光、すなわち、１つの光源から射出された照明光によって生成される波面を表している。

推定標本の像Ｉ_ｅｓｔｚ（ｒ）は、取得位置において、全ての光源から射出された照明光によって生成される。よって、全ての光源について、波面を求めなくてはならない。

（判断結果がＹＥＳの場合：ｚ＝Ｎ_ＩＭ）
ステップＳ２４２が実行される。

波面ｆ_ｉｎ ^ｍ（ｒ）、波面ｇ_ｏｕｔ ^ｍ（ｒ）、波面ｕ_ｚ ^ｍ（ｒ）、及び波面ｕ_ｉｍｇｚ ^ｍ（ｒ）は、ｍ番目の光源から射出された照明光、すなわち、１つの光源から射出された照明光によって生成される波面を表している。

推定標本の像Ｉ_ｅｓｔｚ（ｒ）は、全ての光源から射出された照明光によって生成される。よって、全ての光源について、波面ｆ_ｉｎ ^ｍ（ｒ）、波面ｇ_ｏｕｔ ^ｍ（ｒ）、波面ｕ_ｚ ^ｍ（ｒ）、及び波面ｕ_ｉｍｇｚ ^ｍ（ｒ）を求めなくてはならない。

ステップＳ２５２では、変数ｍの値が光源の数Ｎ_ＬＳと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ２５３が実行される。判断結果がＹＥＳの場合は、ステップＳ２７０が実行される。

（判断結果がＮＯの場合：ｍ≠Ｎ_ＬＳ）
判断結果がＮＯの場合、ステップＳ２５３で、変数ｍの値に１が加算される。ステップＳ２５３が終ると、ステップＳ４１に戻る。

ステップＳ２５３で、変数ｍの値が１つ増えている。そのため、別の光源について、ステップＳ４１で波面ｆ_ｉｎ ^ｍ（ｒ）が算出され、ステップＳ４２で波面ｇ_ｏｕｔ ^ｍ（ｒ）が算出され、ステップＳ２６１で波面ｕ_ｚ ^ｍ（ｒ）が算出され、ステップＳ２６２で波面ｕ_ｉｍｇｚ ^ｍ（ｒ）が算出され，ステップＳ２６３で｜ｕ_ｉｍｇｚ ^ｍ（ｒ）｜^２が算出される。

ステップＳ４１、ステップＳ４２、ステップＳ２５１、及びステップＳ２６０は、全ての光源について｜ｕ_ｉｍｇｚ ^ｍ（ｒ）｜^２が求まるまで、繰り返し行われる。

（判断結果がＹＥＳの場合：ｍ＝Ｎ_ＬＳ）
判断結果がＹＥＳの場合、ステップＳ２７０で、｜ｕ_ｉｍｇｚ ^ｍ（ｒ）｜^２の足し合わせが行われる。その結果、推定標本の像Ｉ_ｅｓｔｚ（ｒ）が算出される。推定標本の像Ｉ_ｅｓｔｚ（ｒ）は、式（３０）で表される。

図１７は、推定標本の像を示す図である。図１７（ａ）は、第１層における推定標本の像を示す図である。図１７（ｂ）は、第４層における推定標本の像を示す図である。

推定標本の像Ｉ_ｅｓｔ１（ｒ）は、全ての光源について波面ｕ_ｉｍｇ１ ^ｍ（ｒ）が求められた場合の像である。推定標本の像Ｉ_ｅｓｔ４（ｒ）は、全ての光源について波面ｕ_ｉｍｇ４ ^ｍ（ｒ）が求められた場合の像である。

図１７（ａ）に示すように、各光源について波面ｕ_ｉｍｇ１ ^ｍ（ｒ）が算出され、波面ｕ_ｉｍｇ１ ^ｍ（ｒ）から｜ｕ_ｉｍｇ１ ^ｍ（ｒ）｜^２が算出され、｜ｕ_ｉｍｇ１ ^ｍ（ｒ）｜^２が全て足し合わされている。その結果、第１層における推定標本の像Ｉ_ｅｓｔ１（ｒ）が算出される。

図１７（ｂ）に示すように、各光源について波面ｕ_ｉｍｇ４ ^ｍ（ｒ）が算出され、波面ｕ_ｉｍｇ４ ^ｍ（ｒ）から｜ｕ_ｉｍｇ４ ^ｍ（ｒ）｜^２が算出され、｜ｕ_ｉｍｇ４ ^ｍ（ｒ）｜^２が全て足し合わされている。その結果、第４層における推定標本の像Ｉ_ｅｓｔ４（ｒ）が算出される。

（判断結果がＹＥＳの場合：ｍ＝Ｎ_ＬＳ）
ステップＳ２８０で、残差が算出される。残差は、式（３１）で表される。式（３１）に示すように、残差は、撮影画像Ｉ_ｍｅａｚ（ｒ）と推定標本の像Ｉ_ｅｓｔｚ（ｒ）とから算出される。

上述のように、撮影画像の数は４で、推定標本の像の数も４である。よって、Ｉ_ｍｅａ１（ｒ）とＩ_ｅｓｔ１（ｒ）から、第１層における残差が算出される。Ｉ_ｍｅａ２（ｒ）とＩ_ｅｓｔ２（ｒ）から、第２層における残差が算出される。Ｉ_ｍｅａ３（ｒ）とＩ_ｅｓｔ３（ｒ）から、第３層における残差が算出される。Ｉ_ｍｅａ４（ｒ）とＩ_ｅｓｔ４（ｒ）から、第４層における残差が算出される。

第１層における残差、第２層における残差、第３層における残差、及び第４層における残差から、ステップ７０で用いられる残差が算出される。

ステップＳ７０では、残差と閾値との比較が行われる。判断結果がＮＯの場合は、ステップＳ２９０が実行される。判断結果がＹＥＳの場合は、ステップＳ１１０が実行される。

（判断結果がＮＯの場合：残差≧閾値）
ステップＳ２９０では、変数Ｌの値が初期化される。後述のステップＳ３０１、ステップＳ３０２、ステップＳ３０３、ステップＳ３０４、及びステップＳ３１０は、ステップＳ２１０で設定した回数だけ実行される。変数Ｌは、これらのステップが実行された回数を表している。

ステップＳ３００は、ステップＳ３０１、ステップＳ３０２、ステップＳ３０３、ステップＳ３０４、ステップＳ３０５、ステップＳ３０６、及びステップＳ３１０を有する。

ステップＳ３０１では、１からＮ_ＩＭまでのなかから１つがランダムに選択される。後述のステップＳ３１１では、補正後の波面が算出される。補正後の波面の算出では、１つの撮影画像と１つの推定標本の像が用いられる。

上述のように、ステップＳ２７０では、複数の推定標本の像が算出される。補正後の波面の算出に用いる推定標本の像は、１つである。よって補正後の波面の算出に用いる推定標本の像を、複数の推定標本の像のなかから選択する。

Ｎ_ＩＭは、層の数である。Ｎ_ＩＭ＝４の場合、ステップＳ３０１では、１から４までの数字のなかから、１つの数字がランダムに選択される。

例えば、選択された数字が１の場合、数字の１は第１層を表している。第１層における推定標本の像には、１番目の取得位置における撮影画像が対応する。よって、補正後の波面の算出には、１番目の取得位置における撮影画像と、第１層における推定標本の像が用いられる。

また、例えば、選択された数字が４の場合、選択された数字は第４層を表している。第４層における推定標本の像には、４番目の取得位置における撮影画像が対応する。よって、補正後の波面の算出には、４番目の取得位置における撮影画像と、第４層における推定標本の像が用いられる。

ステップＳ３０２では、ステップＳ３０１で選択した値が変数ｚＬに入力される。上述のように、ステップＳ３０１では、１からＮ_ＩＭまでの数字のなかから、１つの数字がランダムに選択される。例えば、選択された数字が１の場合、ステップＳ３０２で、変数ｚＬに１が入力されることになる。

ステップＳ３０３では、変数ｍの値が初期化される。後述のステップＳ３１１、ステップＳ３１２、及びステップＳ３１３は、全ての光源に対して実行される。変数ｍは、これらのステップが実行された回数を表している。

ステップＳ３１０は、推定標本の屈折率分布を最適化するステップである。

ステップＳ３１０は、ステップＳ３１１、ステップＳ３１２、ステップＳ３１３、ステップＳ３１４、及びステップＳ３１５を有する。

ステップＳ３１１では、波面ｕ’_ｚＬ ^ｍ（ｒ）が算出される。波面ｕ’_ｚＬ ^ｍ（ｒ）は、変数ｚＬの値で示される層の位置における波面である。

波面ｕ’_ｚＬ ^ｍ（ｒ）の算出では、撮影画像Ｉ_{ｍｅａｚＬ}（ｒ）と推定標本の像Ｉ_{ｅｓｔｚＬ}（ｒ）が用いられる。撮影画像Ｉ_{ｍｅａｚＬ}（ｒ）は、撮影画像Ｉ_ｍｅａｚのうち、変数ｚＬの値で示される位置の撮像画像である。推定標本の像Ｉ_{ｅｓｔｚＬ}（ｒ）は、推定標本の像Ｉ_ｅｓｔｚのうち、変数ｚＬの値で示される位置の推定標本の像である。

波面ｕ’_ｚＬ ^ｍ（ｒ）は式（３２）で表される。

図１８は、波面の補正を示す図である。図１８（ａ）は、推定標本から射出する補正前の波面を示す図である。図１８（ｂ）は、撮影画像の取得位置における補正前の波面を示す図である。図１８（ｃ）は、撮影画像の取得位置における補正後の波面を示す図である。図１８（ｄ）は、推定標本から射出する補正後の波面を示す図である。

ステップＳ３０１で選択された数字が１の場合、すなわち、ｚＬ＝１の場合について説明する。

図１７（ａ）に示すように、推定標本の像Ｉ_ｅｓｔ１（ｒ）は、波面ｕ_ｉｍｇ１ ^ｍ（ｒ）に基づいて算出される。また、図１７（ａ）と図１８（ｂ）に示すように、波面ｕ_ｉｍｇ１ ^ｍ（ｒ）は、波面ｕ_１ ^ｍ（ｒ）に基づいて算出される。

図１８（ａ）に示すように、波面ｕ_１ ^ｍ（ｒ）の算出には、振幅透過率Ｔ_ｚ（ｒ）が用いられている。振幅透過率Ｔ_ｚ（ｒ）は、推定された振幅透過率である。ステップＳ３００の１回目の実行時、この振幅透過率Ｔ_ｚ（ｒ）は、標本２１の振幅透過率と異なる。

振幅透過率Ｔ_ｚ（ｒ）と標本２１の振幅透過率との差が大きくなるほど、推定標本の像Ｉ_ｅｓｔｚ（ｒ）と撮影画像Ｉ_ｍｅａｚ（ｒ）との差も大きくなる。よって、推定標本の像Ｉ_ｅｓｔｚ（ｒ）と撮影画像Ｉ_ｍｅａｚ（ｒ）との差は、振幅透過率Ｔ_ｚ（ｒ）と標本２１の振幅透過率との差を反映していると見なすことができる。

上述のように、ｚＬ＝１である。そこで、式（３２）でｚＬ＝１として、推定標本の像Ｉ_ｅｓｔ１（ｒ）と撮影画像Ｉ_ｍｅａ１（ｒ）と用いて、波面ｕ_１ ^ｍ（ｒ）を補正する。その結果、図１８（ｃ）に示すように、補正された波面、すなわち、波面ｕ’_１ ^ｍ（ｒ）が得られる。

波面ｕ’_１ ^ｍ（ｒ）を用いることで、新たな振幅透過率を算出できる。波面ｕ’_１ ^ｍ（ｒ）は、波面ｕ_１ ^ｍ（ｒ）と異なる。よって、新たな振幅透過率は、波面ｕ_１ ^ｍ（ｒ）を算出したときの振幅透過率とは異なる。

ステップＳ３１２では、補正後の波面ｇ’_ｏｕｔ ^ｍ，ｚＬ（ｒ）が算出される。波面ｇ’_ｏｕｔ ^ｍ，ｚＬ（ｒ）は、波面ｕ’_ｚＬ ^ｍ（ｒ）がΔＤだけ伝搬したときの波面である。波面ｇ’_ｏｕｔ ^ｍ，ｚＬ（ｒ）は、式（３３）で表される。

上述のように、波面ｕ’_１ ^ｍ（ｒ）を用いて、振幅透過率Ｔ_ｚ（ｒ）を算出することができる。ただし、図１７（ａ）に示すように、振幅透過率Ｔ_ｚ（ｒ）の算出には、波面ｇ_ｏｕｔ ^ｍ（ｒ）の位置における波面が必要である。

図１８（ａ）、（ｃ）に示すように、波面ｕ’_１ ^ｍ（ｒ）の位置と、波面ｇ_ｏｕｔ ^ｍ（ｒ）の位置は異なる。よって、振幅透過率Ｔ_ｚ（ｒ）を算出するためには、図１８（ｄ）に示すように、波面ｇ’_ｏｕｔ ^ｍ，１（ｒ）が必要である。

波面ｇ’_ｏｕｔ ^ｍ，１（ｒ）は、波面ｕ’_１ ^ｍ（ｒ）が、３×Δｚだけ伝搬したときの波面である。式（３３）で、ΔＤ＝３×Δｚ、ｚＬ＝１とすることで、波面ｇ’_ｏｕｔ ^ｍ，１（ｒ）を算出することができる。

ステップＳ３１３では、標本の勾配ΔＴ_ｚ ^ｍ，ｚＬ（ｒ）が算出される。ΔＴ_ｚ ^ｍ，ｚＬ（ｒ）は、ｍ番目の光源で照明し、且つ、変数ｚＬの値で示される層の位置における撮影画像と推定標本の像で補正した時の標本の勾配である。

標本の勾配ΔＴ_ｚ ^ｍ，ｚＬは式（３４）で表される。標本の勾配ΔＴ_ｚ ^ｍ，ｚＬ（ｒ）の算出には、例えば、勾配降下法を用いることができる。

ここで、
ｆ^＊は、ｆの複素共役
δは、ゼロ除算を防ぐための正規化定数、
である。

上述のように、推定標本３０は、複数の薄い層に置き換えられている。よって、薄い層の各々について、標本の勾配ΔＴ_ｚ ^ｍ，ｚＬ（ｒ）を算出する必要がある。

図１９は、標本の勾配と波面の伝搬を示す図である。推定標本３０が４つの薄い層に置き換えられた場合、つまりＮ_ＩＭ＝４の場合について説明する。また、ｚＬ＝１とする。図１９（ａ）は、標本の勾配を示す図である。図１９（ｂ）は、波面の伝搬を示す図である。

波面ｇ_ｏｕｔ ^ｍ（ｒ）の算出には、振幅透過率Ｔ_４（ｒ）が用いられている。振幅透過率Ｔ_４（ｒ）は、推定された振幅透過率である。よって、この振幅透過率Ｔ_４（ｒ）は、標本２１の振幅透過率と異なる。

振幅透過率Ｔ_４（ｒ）と標本２１の振幅透過率との差が大きくなるほど、波面ｇ_ｏｕｔ ^ｍ（ｒ）と波面ｇ’_ｏｕｔ ^ｍ，１（ｒ）との差も大きくなる。よって、波面ｇ_ｏｕｔ ^ｍ（ｒ）と波面ｇ’_ｏｕｔ ^ｍ，１（ｒ）との差は、振幅透過率Ｔ_４（ｒ）と標本２１の振幅透過率との差を反映していると見なすことができる。

波面ｆ_４ ^ｍ（ｒ）、振幅透過率Ｔ_４（ｒ）、波面ｇ_ｏｕｔ ^ｍ（ｒ）、及び波面ｇ’_ｏｕｔ ^ｍ，１（ｒ）は、既知である。そこで、式（３４）でｚ＝４、ｚＬ＝１とすることで、図１９（ａ）に示すように、標本の勾配ΔＴ_４ ^ｍ，１（ｒ）を算出することができる。

ｇ_４ ^ｍ（ｒ）と波面ｇ_ｏｕｔ ^ｍ（ｒ）は同じなので、ｇ_４ ^ｍ（ｒ）の代わりに、波面ｇ_ｏｕｔ ^ｍ（ｒ）を用いれば良い。また、ｇ’_４ ^ｍ，１（ｒ）はｇ’_ｏｕｔ ^ｍ，１（ｒ）と同じなので、ｇ’_４ ^ｍ，１（ｒ）の代わりに、ｇ’_ｏｕｔ ^ｍ，１（ｒ）を用いれば良い。

次に、標本の勾配ΔＴ_３ ^ｍ，１（ｒ）を算出する。標本の勾配ΔＴ_３ ^ｍ，１（ｒ）の算出には、波面ｇ_３ ^ｍ（ｒ）の位置における波面が必要である。この波面を算出するためには、図１９（ｂ）に示すように、波面ｆ’_４ ^ｍ，１（ｒ）が必要である。

波面ｆ’_４ ^ｍ，１（ｒ）は、式（３５）でｚ＝４、ｚＬ＝１とすることで算出できる。

次に、算出された波面ｆ’_４ ^ｍ，１（ｒ）を用いて、波面ｇ_３ ^ｍ，１（ｒ）の位置における波面を算出する。

図２０は、標本の勾配と波面の伝搬を示す図である。図２０（ａ）は、波面の伝搬と標本の勾配を示す図である。図２０（ｂ）は、波面の伝搬を示す図である。

図２０（ａ）に示すように、波面ｇ’_３ ^ｍ，１（ｒ）は、波面ｆ’_４ ^ｍ，１（ｒ）がΔｚだけ伝搬したときの波面である。これは、第４層から第３層への波面の伝搬である。

上述のように、第３層から第４層への波面の伝搬は、式（２５）で表わされる。よって、式（２５）において、以下のようにすることで、波面ｇ’_３ ^ｍ，１（ｒ）を算出できる。
波面ｆ_４ ^ｍ（ｒ）を、波面ｇ’_３ ^ｍ，１（ｒ）に置き換える。
波面ｇ_３ ^ｍ（ｒ）を、波面ｆ’_４ ^ｍ，１（ｒ）に置き換える。
ΔＤ＝−Δｚとする。

波面ｆ_３ ^ｍ（ｒ）、振幅透過率Ｔ_３（ｒ）、波面ｇ_３ ^ｍ（ｒ）、及び波面ｇ’_３ ^ｍ，１（ｒ）は、既知である。そこで、式（３４）でｚ＝３、ｚＬ＝１とすることで、図２０（ｂ）に示すように、標本の勾配ΔＴ_３ ^ｍ，１（ｒ）を算出することができる。

波面ｆ’_３ ^ｍ，１（ｒ）は、式（３５）でｚ＝３、ｚＬ＝１とすることで算出できる。

第２層と第１層についても、第３層と同様にして、標本の勾配の算出を行えばよい。

図２１は、標本の勾配を示す図である。図２１では、第１層における標本の勾配ΔＴ_１ ^ｍ，１（ｒ）、第２層における標本の勾配ΔＴ_２ ^ｍ，１（ｒ）、第３層における標本の勾配ΔＴ_３ ^ｍ，１（ｒ）、及び第４層における標本の勾配ΔＴ_４ ^ｍ，１（ｒ）、が算出されている。

ステップＳ３１３で得られる標本の勾配ΔＴ_ｚ ^ｍ，１（ｒ）は、ｍ番目の光源で照明し、且つ、第１層の位置における撮影画像と第１層の位置における推定標本の像で補正した時の標本の勾配である。標本の勾配ΔＴ_ｚ ^ｍ，１（ｒ）は、全ての光源から射出された照明光によって決まる。よって、全ての光源について、標本の勾配ΔＴ_ｚ ^ｍ，１（ｒ）を求めなくてはならない。

ステップＳ３１４では、変数ｍの値が光源の数Ｎ_ＬＳと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ３１５が実行される。判断結果がＹＥＳの場合は、ステップＳ３０４が実行される。

（判断結果がＮＯの場合：ｍ≠Ｎ_ＬＳ）
判断結果がＮＯの場合、ステップＳ３１５で、変数ｍの値に１が加算される。ステップＳ３１５が終ると、ステップＳ３１１に戻る。

ステップＳ３１５で、変数ｍの値が１つ増えている。そのため、別の光源について、ステップＳ３１１で波面ｕ_ｚ’^ｍ，１（ｒ）が算出され、ステップＳ３１２で波面ｇ_ｏｕｔ’^ｍ，１（ｒ）が算出され、ステップＳ３１３で標本の勾配ΔＴ_ｚ ^ｍ，１（ｒ）が算出される。

ステップＳ３１１、ステップＳ３１２、及びステップＳ３１３は、全ての光源について標本の勾配ΔＴ_ｚ ^ｍ，１（ｒ）が求まるまで、繰り返し行われる。

図２２は、標本の勾配を示す図である。図２２では、全ての光源について、標本の勾配ΔＴ_ｚ ^ｍ．１（ｒ）が求められている。

（判断結果がＹＥＳの場合：ｍ＝Ｎ_ＬＳ）
判断結果がＹＥＳの場合、ステップＳ３０４で、振幅透過率Ｔ_ｚ（ｒ）が更新される。ステップＳ３０４は、推定標本を更新するステップするステップである。

更新された振幅透過率Ｔ_ｚ（ｒ）は、式（３６）で表される。

ここで、
αは、標本の勾配の補正係数、
である。

ステップＳ３０５では、変数Ｌの値が補正の回数Ｎ_ＣＲと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ３０６が実行される。判断結果がＹＥＳの場合は、ステップＳ３０が実行される。

（判断結果がＮＯの場合：ｍ≠Ｎ_ＣＲ）
判断結果がＮＯの場合、ステップＳ３０６で、変数Ｌの値に１が加算される。ステップＳ３０６が終ると、ステップＳ３０１に戻る。

ステップＳ３０１で、１からＮ_ＩＭまでのなかから１つがランダムに選択される。選択された数字に基づいて、補正に用いる推定標本の像と取得位置が決まる。

そして、ステップＳ３１１で波面ｕ_ｚ’^ｍ，１（ｒ）が算出され、ステップＳ３１２で波面ｇ_ｏｕｔ’^ｍ，１（ｒ）が算出され、ステップＳ３１３で標本の勾配ΔＴ_ｚ ^ｍ，１（ｒ）が算出され、ステップＳ３０４で振幅透過率Ｔ_ｚ（ｒ）が更新される。

ステップＳ３０１、ステップＳ３０２、ステップＳ３０３、ステップＳ３０４、及びステップＳ３１０は、設定された回数の補正が終るまで、繰り返し行われる。

（判断結果がＹＥＳの場合：ｍ＝Ｎ_ＣＲ）
判断結果がＹＥＳの場合、ステップＳ３０に戻る。更新された振幅透過率Ｔ_ｚ（ｒ）で、ステップＳ３０からステップＳ３００までが実行される。

ステップＳ３０からステップＳ３００までが繰り返し実行されることで、更新された振幅透過率Ｔ_ｓ（ｒ）は、徐々に、標本２１の振幅透過率に近づく。すなわち、残差が小さくなる。やがて、残差は閾値よりも小さくなる。

（判断結果がＹＥＳの場合：残差＜閾値）
ステップＳ１１０では、推定標本の屈折率分布が算出される。得られた振幅透過率Ｔ_ｚ（ｒ）は、標本２１の振幅透過率と同一か、又は、略同一である。得られた振幅透過率Ｔ_ｚ（ｒ）と式（１）から、屈折率分布ｎ_ｚ（ｒ）が求まる。

ステップＳ１１０で得られた屈折率分布ｎ_ｚ（ｒ）を用いることで、推定標本の構造を再構成することができる。再構成された推定標本の構造は、例えば、表示装置に出力することができる。推定標本３０は、厚い標本である。第２のシミュレーションでは、厚い標本の構造は、推定標本の３次元構成を再構成することができる。

上述のように、ステップＳ１１０で得られた振幅透過率Ｔ_ｚ（ｒ）は、標本２１の振幅透過率と同一か、又は、略同一である。この場合、屈折率分布ｎ_ｚ（ｒ）も、標本２１の屈折率分布と同一か、又は、略同一と見なすことができる。よって、再構成された推定標本３０の構造は、標本６の構造と同一か、又は、略同一と見なすことができる。

第２のシミュレーションでは、ステップＳ２５０、ステップＳ２７０、及びステップＳ３１０が、繰り返し実行される。その結果、振幅透過率Ｔ_ｚ（ｒ）が更新される。上述のように、ステップＳ２５０とステップＳ２７０は、推定標本の像を算出するステップである。ステップＳ３１０は、推定標本の屈折率分布を最適化するステップである。

振幅透過率Ｔ_ｚ（ｒ）は、推定標本を表している。よって、推定標本の像を算出するステップと推定標本の屈折率分布を最適化するステップが繰り返し実行されることで、推定標本が更新される。

本実施形態の屈折率分布推定システムは、光源を有し、照明光学系は、コンデンサレンズと、第１開口部材と、を有し、結像光学系は、対物レンズと、結像レンズと、を有し、対物レンズの瞳位置に、第１開口部材の像が結像されることが好ましい。

図２３は、本実施形態の屈折率分布推定システムを示す図である。図１（ａ）と同じ構成については同じ番号を付し、説明は省略する。

屈折率分布推定システム４０は、光源４１を有する。照明光学系２は、コンデンサレンズ４２と、第１開口部材４３と、レンズ４４と、を有する。結像光学系３は、対物レンズ４５と、結像レンズ４６と、を有する。

光源４１から射出された光は、レンズ４４で集光される。集光位置に、第１開口部材４３が配置されている。第１開口部材４３は、コンデンサレンズ４２の前側焦点位置に配置されている。よって、コンデンサレンズ４２から、平行光束が標本６に射出される。

標本６は、対物レンズ４５の前側焦点位置に配置されている。標本６から対物レンズ４５に入射した光は、対物レンズ４５から平行光束となって射出される。対物レンズ４５から射出された平行光束は、結像レンズ４６で集光される。

対物レンズ４５の前側焦点位置は、コンデンサレンズ４２の後側焦点位置と一致している。よって、対物レンズ４５の瞳位置Ｐｕに、第１開口部材４３の像が結像成される。照明光学系２と結像光学系３には、例えば、顕微鏡の光学系を用いることができる。

本実施形態の屈折率分布推定システムでは、第１開口部材は、輪帯形状の透過部又は減光部を有することが好ましい。

図２４は、第１開口部材を示す図である。図２４（ａ）は、第１例を示す図である。図２４（ｂ）は、第２例を示す図である。

第１例では、第１開口部材５０は、輪帯形状の透過部５１と、円形の遮光部５２と、輪帯形状の遮光部５３と、を有する。透過部５１は、減光部にしても良い。

第２例では、第２開口部材６０は、輪帯形状の透過部６１と、円形の遮光部６２と、輪帯形状の遮光部６３と、を有する。透過部６１は、減光部にしても良い。

第１開口部材５０と第２開口部材６０では、透過部の位置と、透過部の幅が異なる。透過部５１は、透過部６１よりも外側に位置している。透過部５１の幅は、透過部６１の幅よりも狭い。

本実施形態の屈折率分布推定システムは、第１開口部材とは異なる第２開口部材と、第１開口部材と第２開口部材とを切り替える移動機構と、を備えることが好ましい。

図２５は、本実施形態の屈折率分布推定システムを示す図である。図９（ａ）と同じ構成については同じ番号を付し、説明は省略する。

屈折率分布推定システム７０は、複数の開口部材を有する。例えば、第１開口部材５０と、第２開口部材６０と、を使用することができる。第１開口部材５０と第２開口部材６０は、移動機構７１に配置することができる。

移動機構７１として、例えば、スライダーを用いることができる。この場合、スライダーを移動することで、第１開口部材５０と第２開口部材６０のどちらか一方を、光軸ＡＸ上に位置させることができる。

移動機構７１として、例えば、ターレットを用いることができる。この場合、ターレットを回転することで、第１開口部材５０と第２開口部材６０のどちらか一方を、光軸ＡＸ上に位置させることができる。

図２６は、第３のシミュレーションのフローチャートである。第２のフローチャートにおける処置と同じ処理につては、同じ番号を付し、説明は省略する。

ステップＳ４００では、開口部材の数Ｎ_ＡＰが設定される。上述のように、屈折率分布推定システム７０では、複数の開口部材が使用される。

ステップＳ４１０では、変数ｎの値が初期化される。後述のステップＳ２３１、ステップＳ５０１、ステップＳ５０２、ステップＳ７０、ステップＳ２９０、及びステップＳ６００は、少なくとも１つの開口部材に対して実行される。変数ｎは、これらのステップが実行された回数を表している。

ステップＳ５００は、ステップＳ２３１、ステップＳ５０１、ステップＳ５０２、ステップＳ７０、ステップＳ２９０、ステップＳ５０３、ステップＳ５０４、ステップＳ５０５、及びステップＳ６００を有する。

ステップＳ２３１では、推定標本３０における振幅透過率Ｔ_ｚ（ｒ）に初期値が設定される。振幅透過率Ｔ_ｚ（ｒ）は、ｚ番目の層における振幅透過率である。

図２６では、ステップＳ２３１しか図示されていない。しかしながら、振幅透過率Ｔ_ｚ（ｒ）を算出する具体的な処理は、第２のシミュレーションにおける振幅透過率Ｔ_ｚ（ｒ）の算出処理と同じである。よって、振幅透過率Ｔ_ｚ（ｒ）の数は、層の数と同じである。

ステップＳ５０１では、推定標本の像Ｉ_ｅｓｔｚ ^ｎ（ｒ）が算出される。推定標本の像Ｉ_ｅｓｔｚ ^ｎ（ｒ）は、ｎ番目の開口部材を用いたときのｚ番目の層における推定標本の像である。

図２６では、ステップＳ５０１しか図示されていない。しかしながら、推定標本の像Ｉ_ｅｓｔｚ ^ｎ（ｒ）を算出する具体的な処理は、第２のシミュレーションにおける推定標本の像Ｉ_ｅｓｔｚ（ｒ）の算出処理と同じである。よって、推定標本の像Ｉ_ｅｓｔｚ ^ｎ（ｒ）の数は、層の数と同じである。

推定標本の像Ｉ_ｅｓｔｚ ^ｎ（ｒ）の算出では、以下の式（３７）〜式（４４）が用いられる。

ここで、
ΔＤは、隣り合う２つの層の間隔、
である。

ここで、
△Ｄは、推定標本の表面から撮影画像の取得位置までの距離、
である。

ステップＳ７０では、残差と閾値との比較が行われる。判断結果がＮＯの場合は、ステップＳ２９０が実行される。判断結果がＹＥＳの場合は、ステップＳ１１０が実行される。

（判断結果がＮＯの場合：残差≧閾値）
ステップＳ２９０で変数Ｌの値が初期化された後、ステップＳ６００が実行される。

ステップＳ６０１で、標本の勾配ΔＴ_ｚ ^{ｎ，ｍ，ｚＬ}（ｒ）が算出される。標本の勾配ΔＴ_ｚ ^{ｎ，ｍ，ｚＬ}（ｒ）は、ｎ番目の開口部材を用いたときのｚ番目の層における標本の勾配である。

標本の勾配ΔＴ_ｚ ^{ｎ，ｍ，ｚＬ}（ｒ）では、ｍ番目の光源で標本が照明されている。また、ｚＬ番目の取得位置における撮影画像と、ｚＬ番目の層における推定標本の像を用いて、補正が行われている。

標本の勾配ΔＴ_ｚ ^{ｎ，ｍ，ｚＬ}（ｒ）の算出では、以下の式（４５）〜式（４９）が用いられる。

ここで、
△Ｄは、ＺＬ番目の層から推定標本の表面までの距離、
である。

ステップＳ６０１で得られる標本の勾配ΔＴ_ｚ ^{ｎ，ｍ，ＺＬ}（ｒ）は、ｎ番目の開口部材を用い、ｍ番目の光源で照明し、且つ、ｚＬ番目の位置における撮影画像とｚＬ番目の位置における推定標本の像で補正した時の標本の勾配である。標本の勾配ΔＴ_ｚ ^{ｎ，ｍ，ＺＬ}（ｒ）は、全ての光源から射出された照明光によって決まる。よって、全ての光源について、標本の勾配ΔＴ_ｚ ^{ｎ，ｍ，ｚＬ}（ｒ）を求めなくてはならない。

ステップＳ６０２では、変数ｍの値が光源の数Ｎ_ＬＳと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ６０３が実行される。判断結果がＹＥＳの場合は、ステップＳ６０４が実行される。

（判断結果がＮＯの場合：ｍ≠Ｎ_ＬＳ）
判断結果がＮＯの場合、ステップＳ６０３で、変数ｍの値に１が加算される。ステップＳ６０３が終ると、ステップＳ３０１に戻る。

ステップＳ６０３で、変数ｍの値が１つ増えている。そのため、別の光源について、ステップＳ６０１で標本の勾配Ｔ_ｚ ^{ｎ，ｍ，ｚＬ}（ｒ）が算出される。

（判断結果がＹＥＳの場合：ｍ＝Ｎ_ＬＳ）
判断結果がＹＥＳの場合、ステップＳ６０４で、振幅透過率Ｔ_ｚ（ｒ）が更新される。ステップＳ６０４は、推定標本を更新するステップするステップである。

更新された振幅透過率Ｔ_ｚ（ｒ）は、式（５０）で表される。

ここで、
αは、標本の勾配の補正係数、
である。

ステップＳ５０３では、変数ｎの値が開口部材の数Ｎ_ＡＰと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ５０４が実行される。判断結果がＹＥＳの場合は、ステップＳ５０５が実行される。

（判断結果がＮＯの場合：ｎ≠Ｎ_ＡＰ）
判断結果がＮＯの場合、ステップＳ５０４で、変数ｎの値に１が加算される。ステップＳ５０４が終ると、ステップＳ２３１に戻る。

ステップＳ５０４で、変数ｎの値が１つ増えている。そのため、別の開口部材について、ステップＳ５００が実行される。

（判断結果がＹＥＳの場合：ｎ＝Ｎ_ＡＰ）
判断結果がＹＥＳの場合、ステップＳ５０５で、変数ｎの値に１が設定される。ステップＳ５０５が終ると、ステップＳ２３１に戻る。

全ての開口部材についてステップＳ５００を実行しても、残差が閾値よりも大きい場合がある。このような場合は、再び、最初の開口部材を用いて、ステップＳ５００が実行される。

（判断結果がＹＥＳの場合：残差＜閾値）
ステップＳ１１０では、推定標本の屈折率分布が算出される。得られた振幅透過率Ｔ_ｚ（ｒ）は、標本２１の振幅透過率と同一か、又は、略同一である。得られた振幅透過率Ｔ_ｚ（ｒ）と式（１５）から、屈折率分布ｎ_ｚ（ｒ）が求まる。

図２７は、開口部材と撮影画像を示す図である。図２７（ａ）は、第１の開口部材と撮影画像を示す図である。図２７（ｂ）は、第２の開口部材と撮影画像を示す図である。図２７（ｃ）は、第３の開口部材と撮影画像を示す図である。

第１の開口部材、第２の開口部材及び第３の開口部材は、輪帯形状の透過部を有する。透過部の位置、透過部の幅は、各々の開口部材で異なる。

撮影画像は、ＸＺ面における光学像の画像である。Ｘ方向は、光軸と直交する方向である。Ｚ方向は光軸方向である。よって、撮影画像では、線状の画像が光軸方向に積み重なっている。

上述のシミュレーションでは、推定標本における振幅透過率Ｔ_ｚ（ｒ）に、初期値が設定されている。第２のシミュレーションで説明したように、初期値の設定に、強度輸送方程式を用いることができる。

図２８は、標本、撮影画像、及び初期値の画像を示す図である。図２８（ａ）は、標本を示す図である。図２８（ｂ）は、撮影画像を示す図である。図２８（ｃ）は、初期値の画像を示す図である。

標本は、フォトニッククリスタルファイバー（以下、「ＰＣＦ」という）である。撮影画像は、ＸＺ面における光学像の画像である。初期値を示す画像は、式（１５）で算出した値に基づく画像である。

図２９は、伝達可能な空間周波数帯域を示す図である。図２９（ａ）は、透過部の形状が円形の場合の空間周波数帯域を示す図である。図２９（ｂ）は、透過部の形状が輪帯形状の場合の空間周波数帯域を示す図である。

２つの図において、横軸はＺ方向、縦軸はＸ方向、又はＹ方向である。Ｙ方向は、Ｚ方向とＸ方向の両方と直交する方向である。また、空間周波数帯域の算出では、開口数が１．４の対物レンズを用い、照明光学系の開口数を変化させている。

標本の像が再生されるためには、標本の持つ空間周波数が結像面まで伝達されなくてはならない。伝達可能な空間周波数が決まる要因として、照明光学系の開口数と、結像光学系の開口数がある。

透過部の形状が円形の場合、図２９（ａ）に示すように、破線で示す範囲の内側には、空間周波数帯域が存在しない。これに対して、透過部の形状が輪帯形状の場合、図２９（ｂ）に示すように、破線で示す範囲の内側には、空間周波数帯域が存在する。

破線で示す範囲は、空間周波数が低い範囲である。よって、透過部の形状が輪帯形状の場合、低い空間周波数から高い空間周波数までを、伝達することができる。その結果、標本の像を、より正確に再生することができる。

輪帯形状の位置、輪帯形状の幅によって、伝達可能な空間周波数は変化する。

図３０は、標本を示す図である。図３０（ａ）は、標本の全体を示す断面図である。図３０（ｂ）は、標本の一部を示す断面図である。

標本は、ＰＣＦである。ＰＣＦでは、クラッドの直径は１８０．０μｍで、穴の直径は６．０μｍで、穴の間隔は１２．９μｍである。クラッドの屈折率は１．４５６である。

ＰＣＦは液体に浸されている。よって、穴は液体で満たされている。また、ＰＣＦの外側も、同じ液体で覆われている。液体の屈折率は１．４３６である。クラッドの屈折率は、クラッド以外の部分の屈折率よりも０．０２だけ高い。

照明光の波長は１５００ｎｍである。また、対物レンズの開口数は１．４である。

図３１は、撮影画像、伝達可能な空間周波数を示す図、再生像の全体の画像、及び再生像の一部の画像である。図３１（ａ）は、第１輪帯における撮影画像である。図３１（ｂ）は、第２輪帯における撮影画像である。図３１（ｃ）は、第３輪帯における撮影画像である。図３１（ｄ）は、第４輪帯における撮影画像である。

輪帯の位置と輪帯の幅は、照明光学系の開口数で表すことができる。上述のように、４つの輪帯では、以下のようになっている。数値は、照明光学系の開口数である。左側の数値は、輪帯の内側を表している。右側の数値は、輪帯の外側を表している。

第１輪帯：０．１４−０．５４。
第２輪帯：０．７−０．９。
第３輪帯：１．０−１．１。
第４輪帯：１．１５−１．２５。

照明光学系の開口数の数値が大きくなるほど、輪帯は外側に位置する。よって、第１輪帯が最も内側に位置し、第４輪帯が最も外側に位置している。

図３１（ａ）、図３１（ｂ）、図３１（ｃ）、図３１（ｄ）は、ＸＺ面における光学像の画像である。

図３１（ａ）、図３１（ｂ）、図３１（ｃ）、及び図３１（ｄ）を比較すると、第１輪帯、第２輪帯、第３輪帯、第４輪帯の順で、画像が鮮明になっている。すなわち、輪帯が外側に位置するほど、画像が鮮明になっている。

画像は、伝達可能な空間周波数帯域に含まれる空間周波数が高くなるほど鮮明になる。上述のように、輪帯が外側に位置するほど、画像が鮮明になっている。よって、輪帯が外側に位置するほど、伝達可能な空間周波数帯域に含まれる空間周波数が高くなる。

このように、輪帯の位置と伝達可能な空間周波数には、相関がある。そこで、複数の輪帯を組み合わせることで、伝達可能な空間周波数帯域に含まれる空間周波数を変えることができる。

図３１（ｅ）、図３１（ｆ）、図３１（ｇ）、及び図３１（ｈ）では、伝達可能な空間周波数を示し、輪帯の数と組み合わせは、以下のようになっている。
図３１（ｅ）：第１輪帯、第２輪帯、第３輪帯、第４輪帯。
図３１（ｆ）：第１輪帯、第４輪帯。
図３１（ｇ）：第４輪帯。
図３１（ｈ）：第１輪帯。

図３１（ｅ）、図３１（ｆ）、図３１（ｇ）、及び図３１（ｈ）を比較すると、輪帯が外側に位置するほど、伝達可能な空間周波数に含まれる空間周波数は高くなる。また、輪帯の組み合わせ数が多くなるほど、伝達可能な空間周波数に含まれる空間周波数が多くなる。

標本の像の再生では、撮影画像が用いられる。標本の像をより正確に再生するためには、撮影画像に、低い空間周波数から高い空間周波数までが含まれていると良い。

図３１（ｉ）、図３１（ｊ）、図３１（ｋ）、図３１（ｌ）、図３１（ｍ）、図３１（ｎ）、図３１（ｏ）、及び図３１（ｐ）では、再生像を示し、輪帯の数と組み合わせは、以下のようになっている。
図３１（ｉ）、（ｍ）：第１輪帯、第２輪帯、第３輪帯、第４輪帯。
図３１（ｊ）、（ｎ）：第１輪帯、第４輪帯。
図３１（ｋ）、（ｏ）：第４輪帯。
図３１（ｌ）、（ｐ）：第１輪帯。

図３１（ｉ）、図３１（ｊ）、図３１（ｋ）、及び図３１（ｌ）を比較すると、輪帯が外側に位置するほど、再生像が鮮明になっている。また、輪帯の組み合わせ数が多くなるほど、再生像が鮮明になっている。

４つの輪帯を用いた場合、２つの輪帯を用いた場合に比べて、ＰＣＦの外側のアーチファクトが少ない。そのため、図３１（ｍ）と図３１（ｎ）から分かるように、４つの輪帯を用いた方が、実際の穴の形状に近くなっている。

穴の径は、ＰＣＦの外径よりも小さい、穴の方が、ＰＣＦの外形（輪郭）よりも、高い空間周波数を有する。よって、第４輪帯を用いると、再生像において、高い空間周波数の部分を、鮮明に再生することができる。そのため、図３１（ｋ）、（ｏ）に示すように、穴の再構成はできている。

しかしながら、第４輪帯を用いると、再生像において、低い空間周波数の部分については情報が欠落している。そのため、図３１（ｋ）に示すように、ＰＣＦの外形が上下に伸びていている。

また、第１輪帯を用いると、図３１（ｌ）に示すようにＰＣＦの外形が、概ね再生できている。しかしながら、第１輪帯を用いると、再生像において、高い低い空間周波数の部分については情報が欠落している。そのため、図３１（ｐ）に示すように、穴の再構成ができていない。

このように、複数の輪帯を組み合わせることで、像の再構成に利用できる空間周波数帯域を広げることができる。その結果、像の再構成の性能を向上させることができる。

図３２は、標本を示す図と再生像の画像である。図３２（ａ）は、第１標本を示す図である。図３２（ｂ）は、第２標本を示す図である。図３２（ｃ）は、第１標本の再生像の画像である。図３２（ｄ）は、第２標本の再生像の画像である。

第１標本は、５０μｍの細胞の塊である。第２標本は、１００μｍの細胞の塊である。照明光の波長は、１．５μｍである。１つの細胞の大きさは、１０μｍである。核の大きさは３μｍである。細胞質の屈折率は、１．３６である。 核の屈折率は、１．４である。溶液の屈折率１．３３である。

図３２（ｃ）と図３２（ｄ）に示すように、第１標本と第２標本のいずれについても、鮮明な再生像が得られている。

図３３は、標本を示す図、初期値の画像、及び再生像の画像である。図３３（ａ）は、標本を示す図である。図３３（ｂ）は、初期値の画像である。図３３（ｃ）は、標本の再生像の画像である。

標本は、５０μｍの細胞塊である。図３３（ｃ）に示すように、３次元の標本でも鮮明で正確な再生像をえることができている。

本実施形態の屈折率分布推定システムでは、推定標本像を算出するステップは、照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を計算し、複数の第１波面が推定標本を通過した後の複数の第２波面を計算し、複数の第２波面を用いて結像光学系の結像位置における複数の第１強度分布を計算し、複数の第１強度分布を足し合わせることにより、推定標本の像を算出し、標本の屈折率分布を最適化するステップは、複数の光源のそれぞれについて、複数の第２波面を撮影画像及び推定標本の像を用いて補正した複数の第２補正波面を算出し、複数の第２波面と複数の第２補正波面との誤差から推定標本の屈折率分布の勾配を算出し、屈折率分布の勾配を用いて推定標本の屈折率分布を最適化することが好ましい。

第１波面の計算は、ステップＳ４１で行われる。第２波面の計算は、ステップＳ４２で行われる。第１強度分布の計算は、ステップＳ４５、又はステップＳ２６３で行われる。推定標本の像の算出は、ステップＳ５０、又はステップＳ２７０で行われる。

標本の屈折率分布を最適化するステップは、ステップＳ９０、又はステップＳ３１０である。第２補正波面の算出は、ステップＳ９１、又はステップＳ３１１で行われる。推定標本の屈折率分布の勾配の算出は、ステップＳ９２、又はステップＳ３１２で行われる。

本実施形態の屈折率分布推定システムでは、推定標本像を算出するステップは、照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を計算するステップと、複数の第１波面が推定標本を通過した後の複数の第２波面を計算するステップと、複数の第２波面から結像光学系の標本側の焦点位置における複数の第３波面を計算するステップと、複数の第３波面と結像光学系の瞳関数とを用いて結像光学系の結像位置における複数の第４波面を計算し、複数の第４波面をそれぞれ２乗して複数の第１強度分布を計算するステップと、複数の第１強度分布を足し合わせることにより、推定標本像を算出するステップと、を有することが好ましい。

第１波面を計算するステップは、ステップＳ４１である。第２波面を計算するステップは、ステップＳ４２である。第３波面を計算するステップは、ステップＳ４３、又はステップＳ２６２である。

第１強度分布を計算するステップは、ステップＳ４４とステップＳ４５、又はステップＳ２６２とステップＳ２６３である。推定標本像を算出するステップは、ステップＳ５０、又はステップＳ２７０である。

本実施形態の屈折率分布推定システムでは、標本の屈折率分布を最適化するステップは、複数の光源のそれぞれについて、複数の第２波面を撮影画像及び推定標本の像を用いて補正した複数の第２補正波面を算出し、複数の第２波面と複数の第２補正波面との誤差から推定標本の屈折率分布の勾配を算出し、屈折率分布の勾配を用いて推定標本の屈折率分布を最適化することが好ましい。

本実施形態の屈折率分布推定システムでは、ＴＶ正則化を用いて振幅透過率を更新することが好ましい。

ＴＶ正則化は、エッジを保持しつつ、振動成分を抑制できる。そのため、ノイズ除去やぼけ画像の修正などの画像処理で用いられる。

振幅透過率Ｔ_ｓ（ｒ）又は振幅透過率Ｔ_ｚ（ｒ）を、式（５１）に示される最小化問題によって、Ｔ’_ｓ（ｒ）又は振幅透過率Ｔ’_ｚ（ｒ）に更新する。

右辺の第１項は、推定残差のＬ２ノルムを示す。右辺の第２項は、正則化項と呼ばれ、振幅透過率Ｔ_ｓ（ｒ）又は振幅透過率Ｔ_ｚ（ｒ）の局所変化が少なければ少ないほど、小さい値をとる性質の関数が一般的に用いられる。τは正則化パラメータと呼ばれる定数である。

この正則化項として、式（５２）に示す隣接する画素間の推定値の差分の絶対値和を意味するＴＶ正則化項を加えると、エッジを保持したまま平滑化することができる。

図３４は、標本を示す図と再生像の画像である。図３４（ａ）は、標本を示す図である。図３４（ｂ）、（ｃ）は撮影画像である。図３４（ｄ）、（ｅ）、（ｆ）、（ｇ）、（ｈ）は、標本の再生像の画像である。

標本は円筒である。円筒の直径は１０μｍである。円筒の屈折率は１．３６である。円筒は液体に浸されている。液体の屈折率は１．３３である。つまり、円筒の屈折率は、円筒以外の部分の屈折率よりも０．０３だけ高い。

照明光の波長は１５００ｎｍである。また、対物レンズの開口数は１．４である。

図３４（ｂ）、（ｃ）は、ＸＺ面における光学像の画像である。図３４（ｂ）の撮影画像では、照明光学系の開口数は１．０である。図３４（ｃ）の撮影画像では、照明光学系の開口数は０．５−１．０である。

照明条件、初期値の有無、ＴＶ正則化の有無と、各図の関係を以下に示す。

図３４（ｄ）と図３４（ｆ）の比較から、初期値を設定すると、より正確に像を再生することができることが分かる。

図３４（ｄ）と図３４（ｅ）の比較から、ＴＶ正則化を行うと、エッジを保持しつつ振動成分を抑制できるため、ノイズ除去され自然な画像に再構成できることが分かる。

図３４（ｇ）と図３４（ｈ）の比較から、輪帯を用いた照明だと、より正確に像を再生することができることが分かる。

本発明には、以下の発明が含まれる。
（付記項１）
複数の方向から同時に入射する光線により標本を照明する照明光学系と、
標本の光学像を形成する結像光学系と、
結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、
撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、
プロセッサは、
標本の屈折率分布を含む推定標本を推定するステップと、
照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、結像光学系の結像位置における複数の第１強度分布を計算し、複数の第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、
複数の第１波面が推定標本を通過した後の複数の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、
推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、
更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする屈折率分布推定システム。

（付記項２）
照明光学系は、空間的パーシャルコヒーレントな光で標本を照明することを特徴とする付記項１に記載の屈折率分布推定システム。

以上のように、本発明は、光軸と直交する面内における分解能高く、標本以外の情報の影響が少ない屈折率分布推定システムに適している。

１ 屈折率分布推定システム
２ 照明光学系
３ 結像光学系
４ 撮像素子
５ プロセッサ
６ 標本
６ａ 表面
６’ 光学像
１０ 推定標本
１０ａ 表面
２０ 屈折率分布推定システム
２１ 標本
２１ａ 表面
２１’ 光学像
２２ 移動ステージ
３０ 推定標本
３０ａ、３０ｂ 表面
４０ 屈折率分布推定システム
４１ 光源
４２ コンデンサレンズ
４３ 第１開口部材
４４ レンズ
４５ 対物レンズ
４６ 結像レンズ
５０ 第１開口部材
６０ 第２開口部材
５１、６１ 透過部
５２、５３、６２、６３ 遮光部
７０ 屈折率分布推定システム
７１ 移動機構
ＡＸ 光軸
Ｌ_ＯＮ、Ｌ_ＯＦＦ 光線
Ｆｏ 結像光学系の焦点位置
Ｐｕ 結像光学系の瞳位置（対物レンズの瞳位置）
ＩＬＬ 照明領域
Ｐ_ＯＢ 結像光学系の瞳
ＩＰ 結像面

上述した課題を解決し、目的を達成するために、本発明の少なくとも幾つかの実施形態に係る屈折率分布推定システムは、
パーシャルコヒーレント照明により標本を照明する照明光学系と、
標本の光学像を形成する結像光学系と、
結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、
撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、
プロセッサは、
標本の屈折率分布を含む推定標本を推定するステップと、
照明光学系の瞳の強度分布をモデル化した複数の光源のそれぞれに対応する第１波面を用いて、複数の光源のそれぞれに対応する結像光学系の結像位置における第１強度分布を計算し、複数の光源のそれぞれに対応する第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、
複数の光源のそれぞれに対応する第１波面のそれぞれが推定標本を通過した後の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、
推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、
更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする。

本発明の少なくとも幾つかの実施形態に係る屈折率分布推定システムは、
複数の方向から同時に入射する光線により標本を照明する照明光学系と、
標本の光学像を形成する結像光学系と、
結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、
撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、
プロセッサは、
標本の屈折率分布を含む推定標本を推定するステップと、
照明光学系の瞳の強度分布をモデル化した複数の光源のそれぞれに対応する第１波面を用いて、複数の光源のそれぞれに対応する結像光学系の結像位置における第１強度分布を計算し、複数の光源のそれぞれに対応する第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、
複数の光源のそれぞれに対応する第１波面のそれぞれが推定標本を通過した後の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、
推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、
更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする。

結像面ＩＰには、撮像素子４の撮像面が位置している。撮像素子４によって、光学像６’の画像の取得が行われる。その結果、図１（ｂ）に示す撮影画像Ｉ_ｍｅａ（ｒ）が得られる。ｒは、（ｘ，ｙ）の２次元座標を示す。

図５では、１番目の光源からＮ_ＬＳ番目までの光源が図示されている。光源は、照明光学系２の瞳位置に配置することができる。このようにすることで、照明光学系２の瞳の強度分布を、複数の光源でモデル化することができる。

上述のように、波面ｕ’_１ ^ｍ（ｒ）を用いて、振幅透過率Ｔ_ｚ（ｒ）を算出することができる。ただし、図１８（ａ）に示すように、振幅透過率Ｔ_ｚ（ｒ）の算出には、波面ｇ_ｏｕｔ ^ｍ（ｒ）の位置における波面が必要である。

Claims (10)

1. パーシャルコヒーレント照明により標本を照明する照明光学系と、
   前記標本の光学像を形成する結像光学系と、
   前記結像光学系により形成された前記標本の光学像から撮影画像を取得する撮像素子と、
   前記撮影画像から前記標本の屈折率分布を再構成するプロセッサと、を備え、
   前記プロセッサは、
   前記標本の屈折率分布を含む推定標本を推定するステップと、
   前記照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、前記結像光学系の結像位置における複数の第１強度分布を計算し、前記複数の第１強度分布を足し合わせることにより、前記推定標本の像を算出するステップと、
   前記複数の第１波面が前記推定標本を通過した後の複数の第２波面、前記撮影画像及び前記推定標本の像を用いて前記推定標本の屈折率分布を最適化するステップと、
   前記推定標本の像の算出と前記推定標本の屈折率分布の最適化とを繰り返して、前記推定標本を更新するステップと、
   更新された前記推定標本の屈折率分布を用いて前記推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする屈折率分布推定システム。
2. 前記結像光学系の焦点位置と前記標本の位置との間隔を、前記結像光学系の光軸方向に変化させる駆動機構を更に備え、前記駆動機構により前記間隔を変化させて、複数の間隔のそれぞれに対応する複数の標本の撮影画像を取得し、
   前記複数の間隔のそれぞれに対応する複数の推定標本の像を算出し、
   前記複数の間隔のそれぞれにおいて、前記推定標本の屈折率分布を最適化することを特徴とする請求項１に記載の屈折率分布推定システム。
3. 光源を有し、
   前記照明光学系は、コンデンサレンズと、第１開口部材と、を有し、
   前記結像光学系は、対物レンズと、結像レンズと、を有し、
   前記対物レンズの瞳位置に、前記第１開口部材の像が結像されることを特徴とする請求項１に記載の屈折率分布推定システム。
4. 前記第１開口部材は、輪帯形状の透過部又は減光部を有することを特徴とする請求項３に記載の屈折率分布推定システム。
5. 前記第１開口部材とは異なる第２開口部材と、
   前記第１開口部材と前記第２開口部材とを切り替える移動機構と、を備えることを特徴とする請求項３に記載の屈折率分布推定システム。
6. 前記推定標本像を算出するステップは、
   前記照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を計算し、
   前記複数の第１波面が推定標本を通過した後の複数の第２波面を計算し、
   前記複数の第２波面を用いて前記結像光学系の結像位置における複数の第１強度分布を計算し、前記複数の第１強度分布を足し合わせることにより、前記推定標本の像を算出し、
   前記標本の屈折率分布を最適化するステップは、
   前記複数の光源のそれぞれについて、前記複数の第２波面を前記撮影画像及び前記推定標本の像を用いて補正した複数の第２補正波面を算出し、
   前記複数の第２波面と前記複数の第２補正波面との誤差から前記推定標本の屈折率分布の勾配を算出し、前記屈折率分布の勾配を用いて前記推定標本の屈折率分布を最適化することを特徴とする請求項１に記載の屈折率分布推定システム。
7. 前記推定標本像を算出するステップは、
   前記照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を計算するステップと、
   前記複数の第１波面が推定標本を通過した後の複数の第２波面を計算するステップと、
   前記複数の第２波面から前記結像光学系の前記標本側の焦点位置における複数の第３波面を計算するステップと、
   前記複数の第３波面と前記結像光学系の瞳関数とを用いて前記結像光学系の結像位置における複数の第４波面を計算し、前記複数の第４波面をそれぞれ２乗して複数の第１強度分布を計算するステップと、
   前記複数の第１強度分布を足し合わせることにより、前記推定標本像を算出するステップと、を有することを特徴とする請求項１に記載の屈折率分布推定システム。
8. 前記標本の屈折率分布を最適化するステップは、
   前記複数の光源のそれぞれについて、前記複数の第２波面を前記撮影画像及び前記推定標本の像を用いて補正した複数の第２補正波面を算出し、
   前記複数の第２波面と前記複数の第２補正波面との誤差から前記推定標本の屈折率分布の勾配を算出し、前記屈折率分布の勾配を用いて前記推定標本の屈折率分布を最適化することを特徴とする請求項７に記載の屈折率分布推定システム。
9. 複数の方向から同時に入射する光線により標本を照明する照明光学系と、
   前記標本の光学像を形成する結像光学系と、
   前記結像光学系により形成された前記標本の光学像から撮影画像を取得する撮像素子と、
   前記撮影画像から前記標本の屈折率分布を再構成するプロセッサと、を備え、
   前記プロセッサは、
   前記標本の屈折率分布を含む推定標本を推定するステップと、
   前記照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、前記結像光学系の結像位置における複数の第１強度分布を計算し、前記複数の第１強度分布を足し合わせることにより、前記推定標本の像を算出するステップと、
   前記複数の第１波面が前記推定標本を通過した後の複数の第２波面、前記撮影画像及び前記推定標本の像を用いて前記推定標本の屈折率分布を最適化するステップと、
   前記推定標本の像の算出と前記推定標本の屈折率分布の最適化とを繰り返して、前記推定標本を更新するステップと、
   更新された前記推定標本の屈折率分布を用いて前記推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする屈折率分布推定システム。
10. 前記照明光学系は、空間的パーシャルコヒーレントな光で前記標本を照明することを特徴とする請求項９に記載の屈折率分布推定システム。
    

Images