(57)【要約】  光軸と直交する面内における分解能が高く、標本以外の情報の影響が少ない屈折率分布推定システムを提供する。  屈折率分布推定システムは、標本を照明する照明光学系と、標本の光学像を形成する結像光学系と、標本の光学像から撮影画像を取得する撮像素子と、撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、プロセッサは、推定標本を推定するステップと、モデル化した複数の光源から射出された複数の第１波面から結像位置における複数の第１強度分布を計算し、複数の第１強度分布から推定標本の像を算出するステップと、複数の第１波面が推定標本を通過した後の複数の第２波面、撮影画像及び推定標本の像から推定標本の屈折率分布を最適化するステップと、推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、推定標本の構造を再構成して出力するステップと、を含む処理を行う。
2021024420.tif 000002
【特許請求の範囲】
【請求項１】
パーシャルコヒーレント照明により標本を照明する照明光学系と、
前記標本の光学像を形成する結像光学系と、
前記結像光学系により形成された前記標本の光学像から撮影画像を取得する撮像素子と、
前記撮影画像から前記標本の屈折率分布を再構成するプロセッサと、を備え、
前記プロセッサは、
前記標本の屈折率分布を含む推定標本を推定するステップと、
前記照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、前記結像光学系の結像位置における複数の第１強度分布を計算し、前記複数の第１強度分布を足し合わせることにより、前記推定標本の像を算出するステップと、
前記複数の第１波面が前記推定標本を通過した後の複数の第２波面、前記撮影画像及び前記推定標本の像を用いて前記推定標本の屈折率分布を最適化するステップと、
前記推定標本の像の算出と前記推定標本の屈折率分布の最適化とを繰り返して、前記推定標本を更新するステップと、
更新された前記推定標本の屈折率分布を用いて前記推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする屈折率分布推定システム。
【請求項２】
前記結像光学系の焦点位置と前記標本の位置との間隔を、前記結像光学系の光軸方向に変化させる駆動機構を更に備え、前記駆動機構により前記間隔を変化させて、複数の間隔のそれぞれに対応する複数の標本の撮影画像を取得し、
前記複数の間隔のそれぞれに対応する複数の推定標本の像を算出し、
前記複数の間隔のそれぞれにおいて、前記推定標本の屈折率分布を最適化することを特徴とする請求項１に記載の屈折率分布推定システム。
【請求項３】
光源を有し、
前記照明光学系は、コンデンサレンズと、第１開口部材と、を有し、
前記結像光学系は、対物レンズと、結像レンズと、を有し、
前記対物レンズの瞳位置に、前記第１開口部材の像が結像されることを特徴とする請求項１に記載の屈折率分布推定システム。
【請求項４】
前記第１開口部材は、輪帯形状の透過部又は減光部を有することを特徴とする請求項３に記載の屈折率分布推定システム。
【請求項５】
前記第１開口部材とは異なる第２開口部材と、
前記第１開口部材と前記第２開口部材とを切り替える移動機構と、を備えることを特徴とする請求項３に記載の屈折率分布推定システム。
【請求項６】
前記推定標本像を算出するステップは、
前記照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を計算し、
前記複数の第１波面が推定標本を通過した後の複数の第２波面を計算し、
前記複数の第２波面を用いて前記結像光学系の結像位置における複数の第１強度分布を計算し、前記複数の第１強度分布を足し合わせることにより、前記推定標本の像を算出し、
前記標本の屈折率分布を最適化するステップは、
前記複数の光源のそれぞれについて、前記複数の第２波面を前記撮影画像及び前記推定標本の像を用いて補正した複数の第２補正波面を算出し、
前記複数の第２波面と前記複数の第２補正波面との誤差から前記推定標本の屈折率分布の勾配を算出し、前記屈折率分布の勾配を用いて前記推定標本の屈折率分布を最適化することを特徴とする請求項１に記載の屈折率分布推定システム。
【請求項７】
前記推定標本像を算出するステップは、
前記照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を計算するステップと、
前記複数の第１波面が推定標本を通過した後の複数の第２波面を計算するステップと、
前記複数の第２波面から前記結像光学系の前記標本側の焦点位置における複数の第３波面を計算するステップと、
前記複数の第３波面と前記結像光学系の瞳関数とを用いて前記結像光学系の結像位置における複数の第４波面を計算し、前記複数の第４波面をそれぞれ２乗して複数の第１強度分布を計算するステップと、
前記複数の第１強度分布を足し合わせることにより、前記推定標本像を算出するステップと、を有することを特徴とする請求項１に記載の屈折率分布推定システム。
【請求項８】
前記標本の屈折率分布を最適化するステップは、
前記複数の光源のそれぞれについて、前記複数の第２波面を前記撮影画像及び前記推定標本の像を用いて補正した複数の第２補正波面を算出し、
前記複数の第２波面と前記複数の第２補正波面との誤差から前記推定標本の屈折率分布の勾配を算出し、前記屈折率分布の勾配を用いて前記推定標本の屈折率分布を最適化することを特徴とする請求項７に記載の屈折率分布推定システム。
【請求項９】
複数の方向から同時に入射する光線により標本を照明する照明光学系と、
前記標本の光学像を形成する結像光学系と、
前記結像光学系により形成された前記標本の光学像から撮影画像を取得する撮像素子と、
前記撮影画像から前記標本の屈折率分布を再構成するプロセッサと、を備え、
前記プロセッサは、
前記標本の屈折率分布を含む推定標本を推定するステップと、
前記照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、前記結像光学系の結像位置における複数の第１強度分布を計算し、前記複数の第１強度分布を足し合わせることにより、前記推定標本の像を算出するステップと、
前記複数の第１波面が前記推定標本を通過した後の複数の第２波面、前記撮影画像及び前記推定標本の像を用いて前記推定標本の屈折率分布を最適化するステップと、
前記推定標本の像の算出と前記推定標本の屈折率分布の最適化とを繰り返して、前記推定標本を更新するステップと、
更新された前記推定標本の屈折率分布を用いて前記推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする屈折率分布推定システム。
【請求項１０】
前記照明光学系は、空間的パーシャルコヒーレントな光で前記標本を照明することを特徴とする請求項９に記載の屈折率分布推定システム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、屈折率分布推定システムに関する。
【背景技術】
【０００２】
実際の物体を計算機上の物体モデルで再現する再構成手法がある。この再構成手法では、測定した物体の画像と計算した物体モデルの画像が一致するように、最適化手法で計算機上の物体モデルを変更していく。最終的に、物体の画像と物体モデルの画像が一致した時に、計算機上の物体モデルは実際の物体を再現している。【０００３】
物体の画像は、測定光学系で取得される。物体モデルの画像がどのような画像になるかは、像演算手法で計算される。そのため、この再構成手法では、測定光学系と、像演算手法の２つが重要である。【０００４】
測定光学系として、例えば、顕微鏡の光学系を用いることができる。顕微鏡の光学系では、ハロゲンランプ又はＬＥＤを用いて標本の画像を取得する。ハロゲンランプとＬＥＤは、インコヒーレント光源である。【０００５】
インコヒーレント光源を用いた照明は、照明条件によって、インコヒーレント照明、コヒーレント照明、及びパーシャルコヒーレント照明に分類される。これらの照明について説明する。【０００６】
顕微鏡では、ケーラー照明が用いられる。ケーラー照明では、光源がコンデンサレンズの焦点面上に配置されるか、又は光源の像がコンデンサレンズの焦点面上に形成される。光源の各点から出た光は、コンデンサレンズで平行光線に変換される。よって、標本は平行光束で照明される。【０００７】
光源の大きさを変えると、標本面での照明光の空間的コヒーレンスが変化する。照明光の空間的コヒーレンスが変化することで、結像特性が変化する。【０００８】
インコヒーレント光源であっても、光源の大きさを非常に小さくすると、光源は点光源とみなせるようになる。点光源からの光を標本に照射する照明は、コヒーレント照明と呼ばれる。【０００９】
光源の大きさを非常に大きくすると、光源は面光源とみなせるようになる。特に、対物レンズの瞳の全域を照明できるように光源の大きさを大きくすると、標本に様々な方向から光源からの光を照射できる。このような照明は、インコヒーレント照明と呼ばれる。【００１０】
面光源において、対物レンズの瞳の一部を照明できるように、光源の大きさを設定する。このような照明は、パーシャルコヒーレント照明（部分的コヒーレント照明）と呼ばれる。パーシャルコヒーレント照明は、コヒーレント照明とインコヒーレント照明の中間の照明である。【００１１】
像演算手法は、線形演算と非線形演算に分類できる。線形演算では、標本内の１回の散乱のみが考慮される。非線形演算では、１回の散乱だけでなく、多重散乱も考慮される。【００１２】
線形演算では、１次ボルン近似を利用している。ボルン近似では、標本で２回以上の散乱が生じても、２回以上の散乱は無視される。線形演算では、標本情報と出力情報が一対一の関係に決まる。そのため、解析的に出力情報を計算できる。出力情報は、例えば、標本の像である。【００１３】
線形演算について、顕微鏡での結像を例に説明する。標本情報（標本の透過率分布）Ｏと出力情報（像強度分布）Ｉを光学系の点像強度分布ＰＳＦのコンボリューションの線形システムとみなせるとき、出力情報Ｉは、以下の式で表される。
Ｉ＝ＰＳＦ＊Ｏ
ここで、＊は、コンボリューションを表す。【００１４】
線形演算では、計算時間は短いが、２回以上の散乱は無視しているため計算精度は低くなる。線形演算を利用して再構成した物体モデルの像は、測定した標本の像を点像強度分布でデコンボリューションして求める。【００１５】
非線形演算は、標本で複数回の散乱が生じることを考慮した演算手法である。非線形演算の１つに、ビーム伝搬法がある。ビーム伝搬法では、物体モデルを複数の薄い層に置き換える。そして、光が各層を通過する際の波面変化を逐次計算して、物体モデルの像を計算する。【００１６】
ビーム伝搬法は、線形演算に比べて、より正確に、物体モデルの像を計算できる。【００１７】
非特許文献１に、測定光学系の照明に顕微鏡のコヒーレント照明を用い、像演算手法にビーム伝搬法を用いて物体の再構成を行う手法が提案されている。【００１８】
この手法では、ＬＥＤアレイが配置された顕微鏡が用いられている。ＬＥＤアレイは、照明光学系の瞳位置に配置されている。ＬＥＤの点灯位置を変えることで、標本に対して様々な角度で照明が照射される。各照射角度で標本の画像を取得するので、複数の標本の画像が取得される。【００１９】
取得した標本の画像から、ライトフィールド手法を用いて、標本の屈折率分布の初期値を求める。ライトフィールド手法は、再構成手法の１つである。偏射照明では、ピント位置からの距離に比例して、標本の像の位置がずれる。ライトフィールド手法は、この標本の像の位置ずれを利用して標本を再構成手法する。【００２０】
初期値が求まったら、勾配降下法を用いて、標本の屈折率分布を正しい値に近づける。勾配降下法は、光源側から順方向に標本内を伝搬した波面と、像情報で補正した波面を像側から逆方向に伝播した波面との差が小さくなるように標本のパラメータを更新していく手法である。波面の伝播計算にはビーム伝搬法が用いられる。
【先行技術文献】
【非特許文献】
【００２１】
【非特許文献１】“3D intensity and phase imaging from light field measurements in an LED array microscope,” L. Tian and L. Waller, Optica 2, 104-111 （2015）.
【発明の概要】
【発明が解決しようとする課題】
【００２２】
コヒーレント照明を用いた結像光学像では、光軸と直交する面内における分解能が、パーシャルコヒーレント照明を用いた結像光学系に比べて低い。さらに被写界深度が深く標本以外の情報、例えばカバーガラスに付着したゴミも標本の画像として取得されてしまう。そのため、再構成した像における屈折率分布の精度が悪くなる。【００２３】
本発明は、このような課題に鑑みてなされたものであって、光軸と直交する面内における分解能が高く、標本以外の情報の影響が少ない屈折率分布推定システムを提供することを目的とする。
【課題を解決するための手段】
【００２４】
上述した課題を解決し、目的を達成するために、本発明の少なくとも幾つかの実施形態に係る屈折率分布推定システムは、
パーシャルコヒーレント照明により標本を照明する照明光学系と、
標本の光学像を形成する結像光学系と、
結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、
撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、
プロセッサは、
標本の屈折率分布を含む推定標本を推定するステップと、
照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、結像光学系の結像位置における複数の第１強度分布を計算し、複数の第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、
複数の第１波面が推定標本を通過した後の複数の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、
推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、
更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする。【００２５】
本発明の少なくとも幾つかの実施形態に係る屈折率分布推定システムは、
複数の方向から同時に入射する光線により標本を照明する照明光学系と、
標本の光学像を形成する結像光学系と、
結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、
撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、
プロセッサは、
標本の屈折率分布を含む推定標本を推定するステップと、
照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、結像光学系の結像位置における複数の第１強度分布を計算し、複数の第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、
複数の第１波面が推定標本を通過した後の複数の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、
推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、
更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする。
【発明の効果】
【００２６】
本発明によれば、光軸と直交する面内における分解能が高く、標本以外の情報の影響が少ない屈折率分布推定システムを提供することができる。
【図面の簡単な説明】
【００２７】
【図１】本実施形態の屈折率分布推定システムと撮影画像を示す図である。
【図２】照明光と照明領域を示す図である。
【図３】第１のシミュレーションのフローチャートである。
【図４】分割された照明領域を示す図である。
【図５】シミュレーションで用いる光学系を示す図である。
【図６】推定標本の像を示す図である。
【図７】波面の補正を示す図である。
【図８】標本の勾配を示す図である。
【図９】標本の勾配を示す図である。
【図１０】本実施形態の屈折率分布推定システムと照明光を示す図である。
【図１１】撮影画像を示す図である。
【図１２】第２のシミュレーションのフローチャートである。
【図１３】第２のシミュレーションのフローチャートである。
【図１４】シミュレーションで用いる光学系を示す図である。
【図１５】各層における波面を示す図である。
【図１６】撮影画像の取得位置における波面と結像面における波面を示す図である。
【図１７】推定標本の像を示す図である。
【図１８】波面の補正を示す図である。
【図１９】標本の勾配と波面の伝搬を示す図である。
【図２０】標本の勾配と波面の伝搬を示す図である。
【図２１】標本の勾配を示す図である。
【図２２】標本の勾配を示す図である。
【図２３】本実施形態の屈折率分布推定システムを示す図である。
【図２４】第１開口部材を示す図である。
【図２５】本実施形態の屈折率分布推定システムを示す図である。
【図２６】第３のシミュレーションのフローチャートである。
【図２７】開口部材と撮影画像を示す図である。
【図２８】標本、撮影画像、及び初期値の画像を示す図である。
【図２９】伝達可能な空間周波数帯域を示す図である。
【図３０】標本を示す図である。
【図３１】撮影画像、伝達可能な空間周波数を示す図、再生像の全体の画像、及び再生像の一部の画像である。
【図３２】標本を示す図と再生像の画像である。
【図３３】標本を示す図、初期値の画像、及び再生像の画像である。
【図３４】標本を示す図と再生像の画像である。
【発明を実施するための形態】
【００２８】
実施例の説明に先立ち、本発明のある態様にかかる実施形態の作用効果を説明する。なお、本実施形態の作用効果を具体的に説明するに際しては、具体的な例を示して説明することになる。しかし、後述する実施例の場合と同様に、それらの例示される態様はあくまでも本発明に含まれる態様のうちの一部に過ぎず、その態様には数多くのバリエーションが存在する。したがって、本発明は例示される態様に限定されるものではない。
【００２９】
本実施形態の屈折率分布推定システムでは、撮影画像と推定標本の像が用いられる。撮影画像は、光学装置で取得した標本の画像である。推定標本の像は、シミュレーションで取得した推定標本の画像である。
【００３０】
シミュレーションでは、残差が閾値以下になるように、推定標本の更新が行われる。残差は、推定標本の像と撮影画像との差である。残差が閾値以下になると、撮影画像と同一の推定標本の像が得られるか、又は、撮影画像と略同一の推定標本の像が得られる。
【００３１】
撮影画像の取得では、パーシャルコヒーレント照明が用いられている。よって、シミュレーションも、パーシャルコヒーレント照明を前提としている。
【００３２】
第１のシミュレーションでは、標本は厚みが薄い標本（以下、「薄い標本」という）である。第２のシミュレーションと第３のシミュレーションでは、標本は厚みが厚い標本（以下、「厚い標本」という）である。
【００３３】
本実施形態の屈折率分布推定システムは、パーシャルコヒーレント照明により標本を照明する照明光学系と、標本の光学像を形成する結像光学系と、結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、プロセッサは、標本の屈折率分布を含む推定標本を推定するステップと、照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、結像光学系の結像位置における複数の第１強度分布を計算し、複数の第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、複数の第１波面が推定標本を通過した後の複数の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする。
【００３４】
図１は本実施形態の屈折率分布推定システムと撮影画像を示す図である。図１（ａ）は、屈折率分布推定システムを示す図である。図１（ｂ）は、撮影画像を示す図である。
【００３５】
図１（ａ）に示すように、屈折率分布推定システム１は、照明光学系２と、結像光学系３と、撮像素子４と、プロセッサ５と、を有する。照明光学系２と結像光学系３で、測定光学系が形成されている。
【００３６】
照明光学系２と結像光学系３との間に、標本６が位置している。標本６は、薄い標本である。結像光学系３の焦点位置Ｆｏは、標本６の内部に位置している。例えば、焦点位置Ｆｏと標本６の表面６ａとの間隔はΔｚ１である。
【００３７】
標本６には、複数の方向から、光線が同時に入射する。図１（ａ）では、光線ＬＯＮと光線ＬＯＦＦが示されている。光線ＬＯＮは、光軸ＡＸと平行な光線である。光線ＬＯＦＦは、光軸ＡＸと交差する光線である。
【００３８】
このように、複数の方向から同時に入射する光線によって、標本６が照明されている。このような照明としては、インコヒーレント照明とパーシャルコヒーレント照明がある。屈折率分布推定システム１では、パーシャルコヒーレント照明が用いられている。
【００３９】
光のコヒーレンスには、時間的コヒーレンスと空間的コヒーレンスが含まれる。ここでのコヒーレンスは、空間的コヒーレンスを指している。
【００４０】
図２は、照明光と照明領域を示す図である。図２（ａ）は、照明光と照明領域の第１例を示す図である。図２（ｂ）は、照明光と照明領域の第２例を示す図である。各図では、標本に入射する照明光と、結像光学系の瞳位置における照明領域が示されている。
【００４１】
第１例では、光線ＬＯＮと光線ＬＯＦＦが照明に用いられている。結像光学系３の瞳位置Ｐｕでは、照明領域ＩＬＬの形状は円である。照明領域ＩＬＬは、結像光学系３の瞳ＰＯＢよりも狭い。よって、第１例における照明は、パーシャルコヒーレント照明である。
【００４２】
第２例では、光線ＬＯＦＦだけが照明に用いられている。瞳位置Ｐｕでは、照明領域ＩＬＬの形状は円環である。照明領域ＩＬＬは、瞳ＰＯＢよりも狭い。よって、第２例における照明は、パーシャルコヒーレント照明である。
【００４３】
図１に戻って説明する。標本６から射出した光は、結像光学系３によって、結像面ＩＰに集光される。結像面ＩＰに、光学像６’が形成される。光学像６’は、標本６の光学像である。
【００４４】
結像面ＩＰには、撮像素子４の撮結像面が位置している。撮像素子４によって、光学像６’の画像の取得が行われる。その結果、図１（ｂ）に示す撮影画像Ｉｍｅａ（ｒ）が得られる。ｒは、（ｘ，ｙ）の２次元座標を示す。
【００４５】
標本６は、薄い標本なので、１つの撮影画像が取得される。よって、結像光学系３と撮像素子４は、光軸方向に移動しない。また、標本６も、光軸方向に移動しない。
【００４６】
撮影画像Ｉｍｅａ（ｒ）は、プロセッサ５に入力される。プロセッサ５では、撮影画像Ｉｍｅａ（ｒ）を用いて、推定標本の像のシミュレーションが行われる。
【００４７】
図３は、第１のシミュレーションのフローチャートである。図１に示す測定光学系では、照明領域は無数の点光源で形成されている。シミュレーションでは、無数の点光源を扱うことはできない。そこで、照明領域を、複数の微小領域に分割する。
【００４８】
微小領域の各々は、点光源とみなすことができる。照明領域を複数の微小領域に分割することで、点光源の数を少なくすることができる。
【００４９】
図４は、分割された照明領域を示す図である。図４（ａ）は、分割された照明領域の第１例を示す図である。図４（ｂ）は、分割された照明領域の第２例を示す図である。
【００５０】
第１例では、円形の照明領域が、複数の微小領域に分割されている。第２例では円環形の照明領域が、複数の微小領域に分割されている。
【００５１】
上述のように、微小領域の各々は、点光源とみなすことができる。よって、第１例では、複数の点光源によって、円形の照明領域が形成されている。第２例では、複数の点光源によって、円環形の照明領域が形成されている。
【００５２】
図５は、シミュレーションで用いる光学系を示す図である。シミュレーションで用いる光学系は、撮影画像Ｉｍｅａ（ｒ）を取得した測定光学系と同一である。シミュレーションでは、標本６の代わりに、推定標本１０が用いられる。
【００５３】
図５には、推定標本１０、波面ｆｉｎｍ（ｒ）、振幅透過率Ｔｓ（ｒ）、波面ｇｏｕｔｍ（ｒ）、撮影画像の取得位置における波面ｕｍ（ｒ）、及び結像面における波面ｕｉｍｇｍ（ｒ）が図示されている。
【００５４】
図５では、１番目の光源からＮＬＳ番目までの光源が図示されている。光源は、照明光学系２の瞳位置に配置することができる。このようにすることで、照明光学系２の瞳の強度分布が、複数の光源でモデル化することができる。
【００５５】
図３に戻って、シミュレーションについて説明する。シミュレーションは、推定標本を推定するステップと、推定標本の像を算出するステップと、推定標本の屈折率分布を最適化するステップと、推定標本を更新するステップと、推定標本の構造を再構成して出力するステップと、を含む。
【００５６】
ステップＳ１０では、光源の数ＮＬＳが設定される。光源の数ＮＬＳは、照明領域を微小領域に分割したときの分割数である。シミュレーションでは、図４（ａ）に示す円形の照明領域が用いられる。
【００５７】
ステップＳ２０は、推定標本を推定するステップである。標本６については、１枚の撮影画像が取得されている。推定標本１０は薄い標本なので、１つの薄い層とみなすことができる。よって、振幅透過率の初期値の設定は１回行われる。
【００５８】
ステップＳ２０では、推定標本１０における振幅透過率Ｔｓ（ｒ）に初期値が設定される。
【００５９】
推定標本１０の像を算出するためには、推定標本１０の情報、例えば、屈折率分布が必要である。推定標本１０は、標本６をモデル化した標本である。よって、推定標本１０の屈折率分布に標本６の屈折率分布を用いることができると良い。
【００６０】
しかしながら、標本６の屈折率分布は、撮影画像Ｉｍｅａ（ｒ）から正確に得られない。よって、推定標本１０の屈折率分布は、推定せざるを得ない。
【００６１】
式（１）に示すように、結像面における屈折率分布ｎｓ（ｒ）は、振幅透過率Ｔｓ（ｒ）に変換することができる。よって、ステップＳ２０では、推定標本１０における振幅透過率Ｔｓ（ｒ）の初期値を設定する。
【数１】
2021024420.tif 000003
ここで、
λは、照明光の波長、
ｄｚは、標本の厚み、
である。
【００６２】
振幅透過率Ｔｓ（ｒ）の値を撮影画像Ｉｍｅａ１（ｒ）から推定できる場合、推定した値を初期値に用いても良い。また、他の方法で振幅透過率Ｔｓ（ｒ）の値を推定できる場合、推定した値を初期値に設定することができる。初期値が推定できない場合、例えば、Ｔｓ（ｒ）＝１とする。
【００６３】
ステップＳ３０では、変数ｍの値が初期化される。後述のステップＳ４１、ステップＳ４２、ステップＳ４３、ステップＳ４４、及びステップＳ４５は、全ての光源に対して実行される。変数ｍは、これらのステップが実行された回数を表している。
【００６４】
ステップＳ４０とステップＳ５０は、推定標本の像を算出するステップである。推定標本の像の数は、撮影画像の数と同数である。撮影画像の数は１なので、推定標本の像の数も１である。
【００６５】
ステップＳ４０は、ステップＳ４１、ステップＳ４２、ステップＳ４３、ステップＳ４４、ステップＳ４５、ステップＳ４６、及びステップＳ４７を有する。
【００６６】
ステップＳ４１では、推定標本１０へ入射する波面ｆｉｎｍ（ｒ）が算出される。ｆｉｎｍ（ｒ）は、１番目の光源からＮＬＳ番目までの光源から射出された光の波面を表している。
【００６７】
上述のように、照明光学系の瞳の強度分布は、複数の光源でモデル化することができる。１番目の光源からＮＬＳ番目までの光源は、モデル化された複数の光源である。第１波面を、モデル化された複数の光源から射出された波面とすると、波面ｆｉｎｍ（ｒ）は、第１波面を表している。
【００６８】
上述のように、微小領域の各々は、点光源とみなすことができる。図５では、ｍ番目の光源から照明光Ｌｍが出射されている。照明光Ｌｍは、推定標本１０に入射する。
【００６９】
この場合、波面ｆｉｎｍ（ｒ）は、式（２）と式（３）で表される。
【数２】
2021024420.tif 000004
【数３】
2021024420.tif 000005
ここで、
θｘ，ｍ、θｙ，ｍは、推定標本への入射角、
である。
【００７０】
ステップＳ４２では、推定標本１０から射出される波面ｇｏｕｔｍ（ｒ）が算出される。薄い標本の場合、波面ｇｏｕｔｍ（ｒ）は、式（４）で表される。
【数４】
2021024420.tif 000006
ここで、
Ｔｓ（ｒ）は、推定標本の振幅透過率、
である。
【００７１】
波面ｇｏｕｔｍ（ｒ）は、波面ｆｉｎｍ（ｒ）が推定標本１０を通過した後の波面である。波面ｆｉｎｍ（ｒ）は第１波面を表しているので、波面ｇｏｕｔｍ（ｒ）は第２波面を表している。
【００７２】
推定標本１０は薄い標本なので、式（４）に示すように、波面ｆｉｎｍ（ｒ）から、波面ｇｏｕｔｍ（ｒ）を直接算出することができる。
【００７３】
ステップＳ４３では、撮影画像の取得位置における波面ｕｍ（ｒ）が算出される。撮影画像の取得位置は、撮影画像の取得が行われたときの、標本側における結像光学系３の焦点位置Ｆｏである。
【００７４】
波面ｕｍ（ｒ）は、式（５）で表される。
【数５】
2021024420.tif 000007
ここで、
Ｐ−Δｚ１｛  ｝は、式（６）で表され、
【数６】
2021024420.tif 000008
Δｚ１は、推定標本の表面から撮影画像の取得位置までの距離、
λは、波長、
ｕは、瞳面座標（ξ，η）の２次元表記、
Ｆ２Ｄは、２次元フーリエ変換、
Ｆ２Ｄ−１は、２次元フーリエ逆変換、
である。
【００７５】
後述のステップＳ６０では、残差が算出される。残差の算出は、撮影画像と推定標本の像が用いられる。推定標本の像を算出するためには、撮影画像の取得位置における波面を求める必要がある。
【００７６】
上述のように、焦点位置Ｆｏと表面６ａとの間隔はΔｚ１である。光の進行方向に向かって測定したときの距離の符号をプラスとすると、撮影画像の取得位置は、表面６ａから−Δｚ１だけ離れた位置になる。
【００７７】
よって、シミュレーションで用いる光学系では、撮影画像の取得位置は、推定標本１０の表面１０ａから−Δｚ１だけ離れた位置になる。この場合、撮影画像の取得位置における波面は、表面１０ａから−Δｚ１だけ離れた位置における波面になる。
【００７８】
式（５）における波面ｕｍ（ｒ）は、波面ｇｏｕｔｍ（ｒ）が、光の進行方向とは逆の方向にΔｚ１だけ伝搬した波面である。この波面は、表面１０ａから−Δｚ１だけ離れた位置における波面である。よって、式（５）における波面ｕｍ（ｒ）は、撮影画像の取得位置における波面を表している。
【００７９】
なお、撮影画像の取得位置と表面６ａの位置は、厳密には異なる。しかしながら、標本６は薄い標本なので、Δｚ１の値は非常に小さい。そのため、撮影画像の取得位置と表面６ａの位置は、ほぼ同じと見なすことができる。
【００８０】
推定標本１０も薄い標本である。そのため、表面１０ａの位置と、表面１０ａから−Δｚ１だけ離れた位置は、ほぼ同じと見なすことができる。すなわち、波面ｇｏｕｔｍ（ｒ）の位置と波面ｕｍ（ｒ）の位置は、ほぼ同じと見なすことができる。この場合、波面ｕｍ（ｒ）の代わりに、波面ｇｏｕｔｍ（ｒ）を用いることもできる。
【００８１】
ステップＳ４４では、結像面における波面ｕｉｍｇｍ（ｒ）が算出される。波面ｕｍ（ｒ）は、結像面ＩＰに伝搬される。このとき、結像光学系３を通過する。結像光学系３は、フーリエ光学系を形成している。よって、式（７）に示すように、結像面ＩＰにおける波面ｕｉｍｇｍ（ｒ）は、波面ｕｍ（ｒ）と結像光学系の瞳関数Ｐ（ｕ）を用いて算出することができる。
【数７】
2021024420.tif 000009
【００８２】
ステップＳ４５では、波面ｕｉｍｇｍ（ｒ）が２乗される。波面ｕｉｍｇｍ（ｒ）は、光の振幅を表している。よって、波面ｕｉｍｇｍ（ｒ）を２乗することで、光強度が算出される。
【００８３】
｜ｕｉｍｇｍ（ｒ）｜２は、結像面ＩＰに光強度分布を表している。第１強度分布を、結像光学系の結像位置における光強度分布とすると、｜ｕｉｍｇｍ（ｒ）｜２は、結像光学系の結像位置における第１光強度分布を表している。
【００８４】
波面ｆｉｎｍ（ｒ）、波面ｇｏｕｔｍ（ｒ）、波面ｕｍ（ｒ）、及び波面ｕｉｍｇｍ（ｒ）は、ｍ番目の光源から射出された照明光、すなわち、１つの光源から射出された照明光によって生成される波面を表している。
【００８５】
推定標本の像Ｉｅｓｔ（ｒ）は、全ての光源から射出された照明光によって生成される。よって、全ての光源について、波面ｆｉｎｍ（ｒ）、波面ｇｏｕｔｍ（ｒ）、波面ｕｍ（ｒ）、及び波面ｕｉｍｇｍ（ｒ）を求めなくてはならない。
【００８６】
ステップＳ４６では、変数ｍの値が光源の数ＮＬＳと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ４７が実行される。判断結果がＹＥＳの場合は、ステップＳ５０が実行される。
【００８７】
（判断結果がＮＯの場合：ｍ≠ＮＬＳ）
判断結果がＮＯの場合、ステップＳ４７で、変数ｍの値に１が加算される。ステップＳ４７が終ると、ステップＳ４１に戻る。
【００８８】
ステップＳ４７で、変数ｍの値が１つ増えている。そのため、別の光源について、ステップＳ４１で波面ｆｉｎｍ（ｒ）が算出され、ステップＳ４２で波面ｇｏｕｔｍ（ｒ）が算出され、ステップＳ４３で波面ｕｍ（ｒ）が算出され、ステップＳ４４で波面ｕｉｍｇｍ（ｒ）が算出され、ステップＳ４５で｜ｕｉｍｇｍ（ｒ）｜２が算出される。
【００８９】
ステップＳ４１、ステップＳ４２、ステップＳ４３、ステップＳ４４、及びステップＳ４５は、全ての光源について｜ｕｉｍｇｍ（ｒ）｜２が求まるまで、繰り返し行われる。
【００９０】
（判断結果がＹＥＳの場合：ｍ＝ＮＬＳ）
判断結果がＹＥＳの場合、ステップＳ５０で、｜ｕｉｍｇｍ（ｒ）｜２の足し合わせが行われる。その結果、推定標本の像Ｉｅｓｔ（ｒ）が算出される。推定標本の像Ｉｅｓｔ（ｒ）は、式（８）で表される。
【数８】
2021024420.tif 000010
【００９１】
図６は、推定標本の像を示す図である。推定標本の像Ｉｅｓｔ（ｒ）は、全ての光源について波面ｕｉｍｇｍ（ｒ）が求められた場合の像である。図６に示すように、各光源について波面ｕｉｍｇｍ（ｒ）が算出され、波面ｕｉｍｇｍ（ｒ）から｜ｕｉｍｇｍ（ｒ）｜２が算出され｜ｕｉｍｇｍ（ｒ）｜２が全て足し合わされている。その結果、推定標本の像Ｉｅｓｔ（ｒ）が算出される。
【００９２】
ステップＳ６０で、残差が算出される。残差は、式（９）で表される。式（９）に示すように、残差は、撮影画像Ｉｍｅａ（ｒ）と推定標本の像Ｉｅｓｔ（ｒ）とから算出される。
【数９】
2021024420.tif 000011
【００９３】
式（９）は、行列のノルムを表している。ノルムは、式（１０）で表される。
【数１０】
2021024420.tif 000012
【００９４】
ステップＳ７０では、残差と閾値との比較が行われる。判断結果がＮＯの場合は、ステップＳ８０が実行される。判断結果がＹＥＳの場合は、ステップＳ１１０が実行される。
【００９５】
（判断結果がＮＯの場合：残差≧閾値）
【００９６】
ステップＳ８０では、変数ｍの値が初期化される。後述のステップＳ９１とステップＳ９２は、全ての光源に対して実行される。変数ｍは、これらのステップが実行された回数を表している。
【００９７】
ステップＳ９０は、推定標本の屈折率分布を最適化するステップである。
【００９８】
ステップＳ９０は、ステップＳ９１、ステップＳ９２、ステップＳ９３、及びステップＳ９４を有する。
【００９９】
ステップＳ９１では、波面ｕ’ｍ（ｒ）が算出される。波面ｕ’ｍ（ｒ）の算出では、撮影画像Ｉｍｅａ（ｒ）と推定標本の像Ｉｅｓｔ（ｒ）が用いられる。また、波面ｕ’ｍ（ｒ）は、撮影画像の取得位置における波面である。
【０１００】
波面ｕ’ｍ（ｒ）は式（１１）で表される。
【数１１】
2021024420.tif 000013
【０１０１】
図７は、波面の補正を示す図である。図７（ａ）は、推定標本から射出する補正前の波面を示す図である。図７（ｂ）は、撮影画像の取得位置における補正前の波面を示す図である。図７（ｃ）は、撮影画像の取得位置における補正後の波面を示す図である。図７（ｄ）は、推定標本から射出する補正後の波面を示す図である。
【０１０２】
図６に示すように、推定標本の像Ｉｅｓｔ（ｒ）は、波面ｕｉｍｇｍ（ｒ）に基づいて算出される。また、図６と図７（ｂ）に示すように、波面ｕｉｍｇｍ（ｒ）は、波面ｕｍ（ｒ）に基づいて算出される。
【０１０３】
図７（ａ）に示すように、波面ｕｍ（ｒ）の算出には、振幅透過率Ｔｓ（ｒ）が用いられている。振幅透過率Ｔｓ（ｒ）は、推定された振幅透過率である。ステップＳ９０の１回目の実行時、この振幅透過率Ｔｓ（ｒ）は、標本６の振幅透過率と異なる。
【０１０４】
振幅透過率Ｔｓ（ｒ）と標本６の振幅透過率との差が大きくなるほど、推定標本の像Ｉｅｓｔ（ｒ）と撮影画像Ｉｍｅａ（ｒ）との差も大きくなる。よって、推定標本の像Ｉｅｓｔ（ｒ）と撮影画像Ｉｍｅａ（ｒ）との差は、振幅透過率Ｔｓ（ｒ）と標本６の振幅透過率との差を反映していると見なすことができる。
【０１０５】
そこで、式（１１）に示すように、推定標本の像Ｉｅｓｔ（ｒ）と撮影画像Ｉｍｅａ（ｒ）と用いて、波面ｕｍ（ｒ）を補正する。その結果、図７（ｃ）に示すように、補正後の波面、すなわち、波面ｕ’ｍ（ｒ）が得られる。
【０１０６】
波面ｕ’ｍ（ｒ）を用いることで、新たな振幅透過率Ｔｓ（ｒ）を算出できる。波面ｕ’ｍ（ｒ）は、波面ｕｍ（ｒ）と異なる。よって、新たな振幅透過率Ｔｓ（ｒ）は、波面ｕｍ（ｒ）を算出したときの振幅透過率とは異なる。
【０１０７】
このように、波面ｕ’ｍ（ｒ）を用いて、振幅透過率Ｔｓ（ｒ）を算出することができる。ただし、図７（ａ）に示すように、振幅透過率Ｔｓ（ｒ）の算出には、波面ｇｏｕｔｍ（ｒ）が必要である。
【０１０８】
図７（ａ）、（ｃ）に示すように、波面ｕ’ｍ（ｒ）の位置と、波面ｇｏｕｔｍ（ｒ）の位置は異なる。よって、振幅透過率Ｔｓ（ｒ）を算出するためには、図７（ｄ）に示すように、波面ｇ’ｏｕｔｍ（ｒ）が必要である。
【０１０９】
波面ｇ’ｏｕｔｍ（ｒ）は、式（１２）で表される。波面ｕ’ｍ（ｒ）は補正後の波面なので、波面ｇ’ｏｕｔｍ（ｒ）も補正後の波面である。
【数１２】
2021024420.tif 000014
【０１１０】
上述のように、撮影画像の取得位置は、表面１０ａから−Δｚ１だけ離れた位置になる。言い換えると、表面１０ａの位置は、撮影画像の取得位置からΔｚ１だけ離れた位置になる。よって、表面１０ａの位置における波面は、撮影画像の取得位置からΔｚ１だけ離れた位置における波面になる。
【０１１１】
式（１２）における波面ｇ’ｏｕｔｍ（ｒ）は、波面ｕ’ｍ（ｒ）が、光の進行方向にΔｚ１だけ伝搬した波面である。この波面は、影画像の取得位置からΔｚ１だけ離れた位置における波面である。よって、式（１２）における波面ｇ’ｏｕｔｍ（ｒ）は、表面１０ａの位置における波面を表している。
【０１１２】
表面１０ａの位置における波面は、ｆｉｎｍ（ｒ）が推定標本１０を通過した後の波面である。上述のように、ｆｉｎｍ（ｒ）は、第１波面を表している。第２波面を、第１波面が推定標本を通過した後の波面とすると、波面ｇ’ｏｕｔｍ（ｒ）は、第２波面を表している。
【０１１３】
上述のように、Δｚ１の値は非常に小さい。また、推定標本１０は薄い標本である。そのため、撮影画像の取得位置と、撮影画像の取得位置からΔｚ１だけ離れた位置は、ほぼ同じと見なすことができる。すなわち、波面ｕ’ｍ（ｒ）の位置と、波面ｇｏｕｔｍ（ｒ）の位置は、ほぼ同じと見なすことができる。この場合、波面ｇ’ｏｕｔｍ（ｒ）の代わりに、波面ｕ’ｍ（ｒ）を用いることもできる。
【０１１４】
ステップＳ９２では、標本の勾配ΔＴｓｍ（ｒ）が算出される。標本の勾配ΔＴｓｍは式（１３）で表される。標本の勾配ΔＴｓｍ（ｒ）の算出には、例えば、勾配降下法を用いることができる。
【数１３】
2021024420.tif 000015
ここで、
ｆ＊は、ｆの複素共役、
δは、ゼロ除算を防ぐための正規化定数、
である。
【０１１５】
図７（ａ）に示すように、波面ｇｏｕｔｍ（ｒ）の算出には、振幅透過率Ｔｓ（ｒ）が用いられている。振幅透過率Ｔｓ（ｒ）は、推定された振幅透過率である。よって、この振幅透過率Ｔｓ（ｒ）は、標本６の振幅透過率と異なる。
【０１１６】
振幅透過率Ｔｓ（ｒ）と標本６の振幅透過率との差が大きくなるほど、波面ｇｏｕｔｍ（ｒ）と波面ｇ’ｏｕｔｍ（ｒ）との差も大きくなる。よって、波面ｇｏｕｔｍ（ｒ）と波面ｇ’ｏｕｔｍ（ｒ）との差は、振幅透過率Ｔｓ（ｒ）と標本６の振幅透過率との差を反映していると見なすことができる。
【０１１７】
波面ｆｉｎｍ（ｒ）、振幅透過率Ｔｓ（ｒ）、波面ｇｏｕｔｍ（ｒ）、及び波面ｇ’ｏｕｔｍ（ｒ）は、既知である。そこで、式（１２）に示すように、波面ｆｉｎｍ（ｒ）、振幅透過率Ｔｓ（ｒ）、波面ｇｏｕｔｍ（ｒ）、及び波面ｇ’ｏｕｔｍ（ｒ）用いて、標本の勾配ΔＴｓｍ（ｒ）を算出することができる。
【０１１８】
図８は、標本の勾配を示す図である。
【０１１９】
ステップＳ９２で得られる標本の勾配ΔＴｓｍ（ｒ）は、１つの光源から射出された照明光における標本の勾配を表している。標本の勾配ΔＴｓｍ（ｒ）は、全ての光源から射出された照明光によって決まる。よって、全ての光源について、標本の勾配ΔＴｓｍ（ｒ）を求めなくてはならない。
【０１２０】
ステップＳ９３では、変数ｍの値が光源の数ＮＬＳと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ９４が実行される。判断結果がＹＥＳの場合は、ステップＳ１００が実行される。
【０１２１】
（判断結果がＮＯの場合：ｍ≠ＮＬＳ）
判断結果がＮＯの場合、ステップＳ９４で、変数ｍの値に１が加算される。ステップＳ９４が終ると、ステップＳ９１に戻る。
【０１２２】
ステップＳ９４で、変数ｍの値が１つ増えている。そのため、別の光源について、ステップＳ９１で波面ｕ’ｍ（ｒ）が算出され、ステップＳ９２で標本の勾配ΔＴｓｍ（ｒ）が算出される。
【０１２３】
ステップＳ９１とステップＳ９２は、全ての光源について標本の勾配ΔＴｓｍ（ｒ）が求まるまで、繰り返し行われる。
【０１２４】
図９は、標本の勾配を示す図である。図９では、全ての光源について、標本の勾配ΔＴｓｍ（ｒ）が求められている。
【０１２５】
（判断結果がＹＥＳの場合：ｍ＝ＮＬＳ）
判断結果がＹＥＳの場合、ステップＳ１００で、振幅透過率Ｔｓ（ｒ）が更新される。ステップＳ１００は、推定標本を更新するステップするステップである。
【０１２６】
更新された振幅透過率Ｔｓ（ｒ）は、式（１４）で表される。
【数１４】
2021024420.tif 000016
ここで、
αは、標本の勾配の補正係数、
である。
【０１２７】
さらに、標本６を吸収の無い完全位相物体として考える場合、式（１５）を用いて、振幅透過率Ｔｓ（ｒ）をさらに更新することができる。
【数１５】
2021024420.tif 000017
【０１２８】
ステップＳ１００が終ると、ステップＳ３０に戻る。更新された振幅透過率Ｔｓ（ｒ）で、ステップＳ３０からステップＳ１００までが実行される。
【０１２９】
ステップＳ３０からステップＳ１００までが繰り返し実行されることで、更新された振幅透過率Ｔｓ（ｒ）は、徐々に、標本６の振幅透過率に近づく。すなわち、残差が小さくなる。やがて、残差は閾値よりも小さくなる。
【０１３０】
（判断結果がＹＥＳの場合：残差＜閾値）
ステップ１１０では、推定標本の屈折率分布が算出される。得られた振幅透過率Ｔｓ（ｒ）は、標本６の振幅透過率と同一か、又は、略同一である。得られた振幅透過率Ｔｓ（ｒ）と式（１）から、屈折率分布ｎ（ｒ）が求まる。
【０１３１】
ステップＳ１１０で得られた屈折率分布ｎ（ｒ）を用いることで、推定標本の構造を再構成することができる。再構成された推定標本の構造は、例えば、表示装置に出力することができる。推定標本１０は、薄い標本である。第１のシミュレーションでは、薄い標本の構造を再構成することができる。
【０１３２】
上述のように、ステップＳ１１０で得られた振幅透過率Ｔｓ（ｒ）は、標本６の振幅透過率と同一か、又は、略同一である。この場合、屈折率分布ｎ（ｒ）も、標本６の屈折率分布と同一か、又は、略同一と見なすことができる。よって、再構成された推定標本１０の構造は、標本６の構造と同一か、又は、略同一と見なすことができる。
【０１３３】
第１のシミュレーションでは、ステップＳ４０、ステップＳ５０、及びステップＳ９０が、繰り返し実行される。その結果、振幅透過率Ｔｓ（ｒ）が更新される。上述のように、ステップＳ４０とステップＳ５０は、推定標本の像を算出するステップである。ステップＳ９０は、推定標本の屈折率分布を最適化するステップである。
【０１３４】
振幅透過率Ｔｓ（ｒ）は、推定標本を表している。よって、推定標本の像を算出するステップと推定標本の屈折率分布を最適化するステップが繰り返し実行されることで、推定標本が更新される。
【０１３５】
本実施形態の屈折率分布推定システムは、結像光学系の焦点位置と標本の位置との間隔を、結像光学系の光軸方向に変化させる駆動機構を更に備え、駆動機構により間隔を変化させて、複数の間隔のそれぞれに対応する複数の標本の撮影画像を取得し、複数の間隔のそれぞれに対応する複数の推定標本の像を算出し、複数の間隔のそれぞれにおいて、推定標本の屈折率分布を最適化することが好ましい。
【０１３６】
図１０は、本実施形態の屈折率分布推定システムと照明光を示す図である。図１０（ａ）は、屈折率分布推定システムを示す図である。図１０（ｂ）は、照明光を示す図である。図１（ａ）と同じ構成については同じ番号を付し、説明は省略する。
【０１３７】
標本２１は、厚い標本である。標本２１には、複数の方向から、光線が同時に入射する。図１０（ａ）では、光線ＬＯＮと光線ＬＯＦＦが示されている。図１０（ｂ）に示すように、光線ＬＯＦＦだけを照明に用いても良い。
【０１３８】
標本２１から射出した光は、結像光学系３によって、結像面ＩＰに集光される。結像面ＩＰに、光学像２１’が形成される。光学像２１’は、標本２１の光学像である。
【０１３９】
屈折率分布推定システム２０は、移動ステージ２２を有する。移動ステージ２２は、光軸ＡＸの方向に移動する。
【０１４０】
上述のように、推定標本の屈折率分布の最適化には、撮影画像が用いられる。標本２１は、厚い標本なので、複数の撮像画像が取得される。複数の撮像画像を取得するために、標本２１を固定し、移動ステージ２２で、結像光学系３の焦点位置を移動させている。
【０１４１】
結像光学系３は、例えば、無限遠補正対物レンズと結像レンズを有する。この場合、対物レンズを移動させることで、結像光学系３の焦点位置を移動させることができる。結像光学系３と撮像素子４を固定し、標本２１を移動させても良い。
【０１４２】
以下、取得された撮像画像が４枚の場合について説明する。
【０１４３】
図１１は、撮影画像を示す図である。図１１（ａ）は、第１の位置における撮影画像を示す図である。図１１（ｂ）は、第２の位置における撮影画像を示す図である。図１１（ｃ）は、第３の位置における撮影画像を示す図である。図１１（ｄ）は、第４の位置における撮影画像を示す図である。
【０１４４】
結像光学系３と標本２１との間隔を変化させることで、標本２１に対する焦点位置Ｆｏが変化する。ここでは、標本２１に対する焦点位置Ｆｏを、４回変化させている。これにより、以下の４つの撮影画像が取得される。
撮影画像Ｉｍｅａ１（ｒ）：表面２１ａからの距離が３×Δｚの位置の画像。
撮影画像Ｉｍｅａ２（ｒ）：表面２１ａからの距離が２×Δｚの位置の画像。
撮影画像Ｉｍｅａ３（ｒ）：表面２１ａからの距離がΔｚの位置の画像。
撮影画像Ｉｍｅａ４（ｒ）：表面２１ａの画像。
【０１４５】
撮影画像Ｉｍｅａ１（ｒ）、撮影画像Ｉｍｅａ２（ｒ）、撮影画像Ｉｍｅａ３（ｒ）、及び撮影画像Ｉｍｅａ４（ｒ）は、プロセッサ５に入力される。プロセッサ５では、４つの撮影画像を用いて、推定標本の像のシミュレーションが行われる。
【０１４６】
図１２と図１３は、第２のシミュレーションのフローチャートである。第１のフローチャートにおける処置と同じ処理につては、同じ番号を付し、説明は省略する。
【０１４７】
図１４は、シミュレーションで用いる光学系を示す図である。図５と同じ構成については同じ番号を付し、説明は省略する。
【０１４８】
シミュレーションで用いる光学系は、撮影画像Ｉｍｅａ１（ｒ）、撮影画像Ｉｍｅａ２（ｒ）、撮影画像Ｉｍｅａ３（ｒ）、及び撮影画像Ｉｍｅａ４（ｒ）を取得した測定光学系と同一である。シミュレーションでは、標本２１の代わりに、推定標本３０が用いられる。
【０１４９】
図１４には、推定標本３０、波面ｆｉｎｍ（ｒ）、振幅透過率Ｔｚ（ｒ）、及び波面ｇｏｕｔｍ（ｒ）が図示されている。
【０１５０】
推定標本が薄い標本の場合、式（４）に示すように、波面ｆｉｎｍ（ｒ）から、波面ｇｏｕｔｍ（ｒ）を直接算出することができる。しかしながら、推定標本が厚い標本の場合、波面ｆｉｎｍ（ｒ）から波面ｇｏｕｔｍ（ｒ）を直接算出することは困難である。
【０１５１】
推定標本３０は、厚い標本である。そこで、光軸方向に沿って、推定標本３０を複数の薄い層に置き換える。そして、薄い層の各々について、層の両側の波面を算出する。
【０１５２】
図１５は、各層における波面を示す図である。波面の算出については後述する。層の数は、取得画像の数と同じにすることができる。ただし、層の数は、取得画像の数よりも多くしても良い。図１５では、層の数は、取得画像の数と同じである。
【０１５３】
図１５では、Ｚ＝１の位置が第１層の位置、Ｚ＝２の位置が第２層の位置、Ｚ＝３の位置が第３層の位置、Ｚ＝４の位置が第４層の位置である。
【０１５４】
図１２、図１３を用いて、シミュレーションについて説明する。
【０１５５】
ステップＳ１０では、光源の数ＮＬＳが設定される。
【０１５６】
ステップＳ２００では、層の数ＮＩＭが設定される。推定標本３０は、厚い標本である。よって、上述のように、推定標本３０を複数の薄い層に置き換える。層の数ＮＩＭは、薄い層の数を表している。
【０１５７】
標本２１では、複数の位置で撮影画像が取得される。層の数ＮＩＭは、撮影画像が取得された位置の数と同じにすることができる。標本２１に対する焦点位置Ｆｏを４回変化させている場合、ＮＩＭ＝４になる。
【０１５８】
１からＮＩＭまでの数字は、薄い層の位置を表している。例えば、ＮＩＭ＝４の場合、数字の１は第１層の位置を表し、数字の２は第２層の位置を表し、数字の３は第３層の位置を表し、数字の４は第４層の位置を表している。
【０１５９】
推定標本の像の算出は、シミュレーションで行われる。そのため、層の数ＮＩＭは、自由に設定できる。例えば、層の数ＮＩＭは、撮影画像が取得された位置の数よりも多くすることができる。
【０１６０】
例えば、ＮＩＭ＝７の場合、薄い層の数は７つである。この場合、７つの推定標本の像が算出される。シミュレーションでは、後述のように、撮影画像と薄い層における推定標本の像が用いられる。よって、推定標本の像が算出される７つ位置には、撮影画像が取得された４つの位置が含まれている。
【０１６１】
７つ位置と撮影画像の関係は、例えば、以下のようにすることができる。
【０１６２】
数字の１は、第１層の位置を表している。この位置では、撮影画像Ｉｍｅａ１（ｒ）が取得されている。また、この位置では、第１層における推定標本の像が算出される。よって、第１層における推定標本の像と撮影画像Ｉｍｅａ１（ｒ）が、後述のステップで用いられる。
【０１６３】
数字の２は、第２層の位置を表している。この位置で取得された撮影画像は無い。
【０１６４】
数字の３は、第３層の位置を表している。この位置では、撮影画像Ｉｍｅａ２（ｒ）が取得されている。また、この位置では、第３層における推定標本の像が算出される。よって、第３層における推定標本の像と撮影画像Ｉｍｅａ２（ｒ）が、後述のステップで用いられる。
【０１６５】
数字の４は、第４層の位置を表している。この位置で取得された撮影画像は無い。
【０１６６】
数字の５は、第５層の位置を表している。この位置では、撮影画像Ｉｍｅａ３（ｒ）が取得されている。また、この位置では、第５層における推定標本の像が算出される。よって、第５層における推定標本の像と撮影画像Ｉｍｅａ３（ｒ）が、後述のステップで用いられる。
【０１６７】
数字の６は、第６層の位置を表している。この位置で取得された撮影画像は無い。
【０１６８】
数字の７は、第７層の位置を表している。この位置では、撮影画像Ｉｍｅａ４（ｒ）が取得されている。また、この位置では、第７層における推定標本の像が算出される。よって、第７層における推定標本の像と撮影画像Ｉｍｅａ４（ｒ）が、後述のステップで用いられる。
【０１６９】
ステップＳ２１０では、補正の回数ＮＣＲが設定される。
【０１７０】
ステップＳ２２０では、変数ｚの値が初期化される。後述のステップＳ２３１は、全ての取得位置に対して実行される。変数ｚは、ステップＳ２３１が実行された回数を表している。
【０１７１】
ステップＳ２３０は、推定標本を推定するステップである。標本２１では、４つの撮影画像が取得されている。上述のように、推定標本３０は４つの薄い層に置き換えられている。よって、振幅透過率の初期値の設定は４回行われる。
【０１７２】
ステップＳ２３０は、ステップＳ２３１、ステップＳ２３２、及びステップＳ２３３を有する。
【０１７３】
ステップＳ２３１では、推定標本３０における振幅透過率Ｔｚ（ｒ）に初期値が設定される。
【０１７４】
初期値の設定では、強度輸送方程式を用いても良い。強度輸送方程式は、例えば、以下の文献に開示されている。
M. R. Teague, “Deterministic phase retrieval: a Greens function solution,” J. Opt. Soc. Am. 73, 1434-1441（1983）
【０１７５】
焦点位置Ｚ０での強度輸送方程式は、式（１６）で表される。
【数１６】
2021024420.tif 000018
ここで、
∇２は、２次のラプラシアン、
ｋは、波数、
φＺ０（ｒ）は，結像面での標本の位相分布、
ＩＺ０は，光学像の平均光強度、
δＩｍｅａＺ０（ｒ）／δＺは、結像面から±△zだけ離れた２つのデフォーカス像の差分像、
である。
【０１７６】
式（１６）を用いると、フォーカス像と２つのデフォーカス像から、標本の位相分布φＺ０（ｒ）を簡単に求めることができる。
【０１７７】
ただし、２つのデフォーカス像の同一点における光強度の差が０か、又は、非常に小さいと、位相を計測できない。パーシャルコヒーレント照明であっても、照明光の開口数が小さい場合は、この光強度の差が０になるか、又は、非常に小さくなる。そのため、このような場合は、強度輸送方程式を用いて初期値を設定することは難しい。
【０１７８】
上述のように、位相分布φＺ０（ｒ）は、フォーカス像と２つのデフォーカス像から算出される。フォーカス像、例えば、対物レンズを一定の間隔で、光軸方向に移動させることで取得される。この場合、光軸に沿って、複数のフォーカス像が離散的に取得される。よって、２つのデフォーカス像も離散的に取得される。
【０１７９】
式（１６）で表される位相分布φＺ０（ｒ）は、光軸と直交する面内における位相分布である。フォーカス像と２つのデフォーカス像は離散的に取得されるので、位相分布φＺ０（ｒ）を表す面も、光軸に沿って離散的に位置している。
【０１８０】
式（１７）に示すように、位相分布φｚ（ｒ）は、振幅透過率Ｔｓ（ｒ）に変換することができる。このようにして、振幅透過率Ｔｚ（ｒ）に初期値を設定することができる。
【数１７】
2021024420.tif 000019
【０１８１】
強度輸送方程式で得られた位相分布φＺ０は、位相分布φｚ（ｒ）に用いることができる。強度輸送方程式を用いて、初期値を設定することができる。なお、初期値の推定が難しい場合、例えば、Ｔｚ（ｒ）＝１としても構わない。
【０１８２】
ステップＳ２３２では、変数ｚの値が取得位置の数ＮＩＭと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ２３３が実行される。判断結果がＹＥＳの場合は、ステップＳ３０が実行される。
【０１８３】
（判断結果がＮＯの場合：ｚ≠ＮＩＭ）
判断結果がＮＯの場合、ステップＳ２３３で、変数ｚの値に１が加算される。ステップＳ２３３が終ると、ステップＳ２３１に戻る。
【０１８４】
ステップＳ２３３で、変数ｚの値が１つ増えている。そのため、別の取得位置について、ステップＳ２３１で振幅透過率Ｔｚ（ｒ）に初期値が設定される。
【０１８５】
ステップＳ２３１は、全ての取得位置について初期値が設定されるまで、繰り返し行われる。
【０１８６】
（判断結果がＹＥＳの場合：ｚ＝ＮＩＭ）
ステップＳ３０で、変数ｍの値が初期化される。後述のステップＳ２４０、ステップＳ４１、ステップＳ４２、ステップＳ２５１、及びステップＳ２６０は、全ての光源に対して実行される。変数ｍは、これらのステップが実行された回数を表している。
【０１８７】
ステップＳ２４０で、関数Ｉｅｓｔｚ（ｒ）の値が初期化される。Ｉｅｓｔｚ（ｒ）は、推定標本３０の像を表している。上述のように、推定標本３０の像は、４つの薄い層に置き換えられている。よって、Ｉｅｓｔｚ（ｒ）は、薄い層の像を表している。
【０１８８】
ステップＳ２５０とステップＳ２７０は、推定標本の像を算出するステップである。推定標本の像の数は、撮影画像の数と同数である。撮影画像の数は４なので、推定標本の像の数も４である。
【０１８９】
ステップＳ２５０は、ステップＳ４１、ステップＳ４２、ステップＳ２５１、ステップＳ２５２、ステップＳ２５３、及びステップＳ２６０を有する。
【０１９０】
ステップＳ４１では、推定標本３０へ入射する波面ｆｉｎｍ（ｒ）が算出される。波面ｆｉｎｍ（ｒ）は、上述の式（２）、（３）で表される。
【０１９１】
ステップＳ４２では、推定標本３０から射出する波面ｇｏｕｔｍ（ｒ）が算出される。波面ｇｏｕｔｍ（ｒ）は、波面ｆｉｎｍ（ｒ）に基づいて算出される。推定標本３０は４つの薄い層に置き換えられている。よって、薄い層の各々で波面を算出する。
【０１９２】
図１５では、Ｚ＝１の位置が第１層の位置、Ｚ＝２の位置が第２層の位置、Ｚ＝３の位置が第３層の位置、Ｚ＝４の位置が第４層の位置である。
【０１９３】
４つの薄い層は、等間隔で並んでいる。隣り合う２つの層の間隔はΔｚである。波面は、２つの層の間を伝搬する。よって、Δｚは、伝搬距離を表している。
【０１９４】
第１層における波面ｆ１ｍ（ｒ）は、式（１８）と式（３）で表される。
【数１８】
2021024420.tif 000020
【０１９５】
第１層の位置は、推定標本３０の表面３０ｂの位置と一致している。表面３０ｂには、波面ｆｉｎｍ（ｒ）が入射する。よって、波面ｆ１ｍ（ｒ）は、波面ｆｉｎｍ（ｒ）を表している。図１５では、波面ｆ１ｍ（ｒ）の代わりに、波面ｆｉｎｍ（ｒ）が示されている。
【０１９６】
第１層における波面ｇ１ｍ（ｒ）は、式（１９）で表される。
【数１９】
2021024420.tif 000021
ここで、
Ｔ１（ｒ）は、第１層における振幅透過率、
である。
【０１９７】
第２層における波面ｆ２ｍ（ｒ）は、波面ｇ１ｍ（ｒ）が、Δｚだけ伝搬したときの波面である。波面ｆ２ｍ（ｒ）は、式（２０）で表される。式（２０）で、ΔＤ＝Δｚとすることで、波面ｆ２ｍ（ｒ）を算出することができる。
【数２０】
2021024420.tif 000022
ここで、ＰΔＤ｛  ｝は、式（２１）で表され、
【数２１】
2021024420.tif 000023
ΔＤは、隣り合う２つの層の間隔、
λは、波長、
ｕは、瞳面座標（ξ，η）の２次元表記、
Ｆ２Ｄは、２次元フーリエ変換、
Ｆ２Ｄ−１は、２次元フーリエ逆変換、
である。
【０１９８】
第２層における波面ｇ２ｍ（ｒ）は、式（２２）で表される。
【数２２】
2021024420.tif 000024
ここで、
Ｔ２（ｒ）は、第２層における振幅透過率、
である。
【０１９９】
第３層における波面ｆ３ｍ（ｒ）は、波面ｇ２ｍ（ｒ）が、Δｚだけ伝搬したときの波面である。第３層における波面ｆ３ｍ（ｒ）は、式（２３）で表される。式（２１）で、ΔＤ＝Δｚとすることで、波面ｆ３ｍ（ｒ）を算出することができる。
【数２３】
2021024420.tif 000025
【０２００】
第３層における波面ｇ３ｍ（ｒ）は、式（２４）で表される。
【数２４】
2021024420.tif 000026
ここで、
Ｔ３（ｒ）は、第３層における振幅透過率、
である。
【０２０１】
第４層における波面ｆ４ｍ（ｒ）は、波面ｇ３ｍ（ｒ）が、Δｚだけ伝搬したときの波面である。第４層における波面ｆ４ｍ（ｒ）は、式（２５）で表される。式（２１）で、ΔＤ＝Δｚとすることで、波面ｆ４ｍ（ｒ）を算出することができる。
【数２５】
2021024420.tif 000027
【０２０２】
第４層における波面ｇ４ｍ（ｒ）は、式（２６）で表される。
【数２６】
2021024420.tif 000028
ここで、
Ｔ４（ｒ）は、第４層における振幅透過率、
である。
【０２０３】
第４層の位置は、推定標本３０の表面３０ａの位置と一致している。表面３０ａからは、波面ｇｏｕｔｍ（ｒ）が射出される。よって、波面ｇ４ｍ（ｒ）は、波面ｇｏｕｔｍ（ｒ）を表している。図１５では、波面ｇ４ｍ（ｒ）の代わりに、波面ｇｏｕｔｍ（ｒ）が示されている。
【０２０４】
以上のように、推定標本が厚い標本の場合は、複数の薄い層に置き換えると共に、２つの層の間を伝搬する波面を求めることで、波面ｇｏｕｔｍ（ｒ）を算出することができる。
【０２０５】
ステップＳ２５１では、変数ｚの値が初期化される。後述のステップＳ２６１、ステップＳ２６２、及びステップＳ２６３は、全ての取得位置に対して実行される。変数ｚは、これらのステップが実行された回数を表している。
【０２０６】
ステップＳ２６０は、ステップＳ２６１、ステップＳ２６２、ステップＳ２６３、ステップＳ２６４、及びステップＳ２６５を有する。
【０２０７】
ステップＳ２６１では、撮影画像の取得位置における波面ｕｚｍ（ｒ）が算出される。波面ｕｚｍ（ｒ）は、式（２７）で表される。
【数２７】
2021024420.tif 000029
ここで、
△Ｄは推定標本の表面から薄い層までの距離、
である。
ＰΔＤ｛  ｝は、式（２８）で表される。
【数２８】
2021024420.tif 000030
【０２０８】
ステップＳ２６２では、結像面における波面ｕｉｍｇｚｍ（ｒ）が算出される。波面ｕｉｍｇｚｍ（ｒ）は、式（２９）で表される。
【数２９】
2021024420.tif 000031
【０２０９】
ステップＳ２６３では、波面ｕｉｍｇｚｍ（ｒ）が２乗される。波面ｕｉｍｇｚｍ（ｒ）は、光の振幅を表している。よって、波面ｕｉｍｇｚｍ（ｒ）を２乗することで、光強度が算出される。
【０２１０】
｜ｕｉｍｇｚｍ（ｒ）｜２は、結像面ＩＰに光強度分布を表している。第１強度分布を、結像光学系の結像位置における光強度分布とすると、｜ｕｉｍｇｚｍ（ｒ）｜２は、結像光学系の結像位置における第１光強度分布を表している。
【０２１１】
ステップＳ２６４では、変数ｚの値が取得位置の数ＮＩＭと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ２６５が実行される。判断結果がＹＥＳの場合は、ステップＳ２５２が実行される。
【０２１２】
（判断結果がＮＯの場合：ｚ≠ＮＩＭ）
判断結果がＮＯの場合、ステップＳ２６５で、変数ｚの値に１が加算される。ステップＳ２６５が終ると、ステップＳ２６１に戻る。
【０２１３】
ステップＳ２６５で、変数ｚの値が１つ増えている。そのため、別の取得位置について、ステップＳ２６１、ステップＳ２６２、及びステップＳ２６３が実行される。
【０２１４】
ステップＳ２６１、ステップＳ２６２、及びステップＳ２６３は、全ての取得位置について初期値が設定されるまで、繰り返し行われる。
【０２１５】
ステップＳ２５０における処理を、第１層と第４層を用いて説明する。第２層と第３層については、第１層と同じように考えれば良い。
【０２１６】
図１６は、撮影画像の取得位置における波面と結像面における波面を示す図である。図１６（ａ）は、第１層におけるにおける２つの波面を示す図である。図１６（ｂ）は、第４層における２つの波面を示す図である。
【０２１７】
ｚ＝１における撮影画像は、撮影画像Ｉｍｅａ１（ｒ）である。撮影画像Ｉｍｅａ１（ｒ）は、表面２１ａからの距離が３×Δｚの位置の画像である。第１層は、表面３０ａから３×Δｚだけ離れている。よって、第１層の位置が、撮影画像Ｉｍｅａ１（ｒ）の取得位置に対応する。
【０２１８】
波面ｇｏｕｔｍ（ｒ）の射出位置は、表面３０ａと一致している。図１６（ａ）に示すように、波面ｇｏｕｔｍ（ｒ）の射出位置は、第１層の位置と異なる。第１層は、波面ｇｏｕｔｍ（ｒ）の射出位置から３×Δｚだけ離れている。
【０２１９】
第１層における波面ｕ１ｍ（ｒ）は、波面ｇｏｕｔｍ（ｒ）が光の進行方向とは反対に３×Δｚだけ伝搬したときの波面である。よって、ステップＳ２６１においてΔＤ＝−３×Δｚとすることで、式（２７）と式（２８）から、波面ｕ１ｍ（ｒ）を算出することができる。
【０２２０】
波面ｕ１ｍ（ｒ）が算出されると、ステップＳ２６２で、式（２９）から、結像面における波面ｕｉｍｇ１ｍ（ｒ）が算出される。
【０２２１】
更に、ステップＳ２６３で、第１層の像の光強度｜ｕｉｍｇ１（ｒ）｜２が算出される。
【０２２２】
ｚ＝２における撮影画像は、撮影画像Ｉｍｅａ２（ｒ）である。撮影画像Ｉｍｅａ２（ｒ）は、表面２１ａからの距離が２×Δｚの位置の画像である。第２層は、表面３０ａから２×Δｚだけ離れている。よって、第２層の位置が、撮影画像Ｉｍｅａ２（ｒ）の取得位置に対応する。
【０２２３】
波面ｇｏｕｔｍ（ｒ）の射出位置は、第２層の位置と異なる。第２層は、波面ｇｏｕｔｍ（ｒ）の射出位置から２×Δｚだけ離れている。
【０２２４】
第２層における波面ｕ２ｍ（ｒ）は、波面ｇｏｕｔｍ（ｒ）が光の進行方向とは反対に２×Δｚだけ伝搬したときの波面である。よって、ステップＳ２６１においてΔＤ＝−２×Δｚとすることで、波面ｕ２ｍ（ｒ）を算出することができる。
【０２２５】
波面ｕ２ｍ（ｒ）が算出されると、ステップＳ２６２で、結像面における波面ｕｉｍｇ２ｍ（ｒ）が算出される。
【０２２６】
更に、ステップＳ２６３で、第２層の像の光強度｜ｕｉｍｇ２（ｒ）｜２が算出される。
【０２２７】
ｚ＝３における撮影画像は、撮影画像Ｉｍｅａ３（ｒ）である。撮影画像Ｉｍｅａ３（ｒ）は、表面２１ａからの距離がΔｚの位置の画像である。第３層は、表面３０ａからΔｚだけ離れている。よって、第３層の位置が、撮影画像Ｉｍｅａ３（ｒ）の取得位置に対応する。
【０２２８】
波面ｇｏｕｔｍ（ｒ）の射出位置は、第３層の位置と異なる。第３層は、波面ｇｏｕｔｍ（ｒ）の射出位置からΔｚだけ離れている。
【０２２９】
第３層における波面ｕ３ｍ（ｒ）は、波面ｇｏｕｔｍ（ｒ）が光の進行方向とは反対にΔｚだけ伝搬したときの波面である。よって、ステップＳ２６１においてΔＤ＝Δｚとすることで、波面ｕ３ｍ（ｒ）を算出することができる。
【０２３０】
波面ｕ３ｍ（ｒ）が算出されると、ステップＳ２６２で、結像面における波面ｕｉｍｇ３ｍ（ｒ）が算出される。
【０２３１】
更に、ステップＳ２６３で、第３層の像の光強度｜ｕｉｍｇ３（ｒ）｜２が算出される。
【０２３２】
ｚ＝４における撮影画像は、撮影画像Ｉｍｅａ４（ｒ）である。撮影画像Ｉｍｅａ４（ｒ）は、表面２１ａの画像である。第４層は、表面３０ａと一致している。よって、第４層の位置が、撮影画像Ｉｍｅａ４（ｒ）の取得位置に対応する。
【０２３３】
波面ｇｏｕｔｍ（ｒ）の射出位置は、表面３０ａである。図１６（ｂ）に示すように、波面ｇｏｕｔｍ（ｒ）の射出位置は、第４層の位置と同じである。
【０２３４】
第４層における波面ｕ４ｍ（ｒ）は、波面ｇｏｕｔｍ（ｒ）同じである。波面ｇｏｕｔｍ（ｒ）を波面ｕ４ｍ（ｒ）に置き換えることができる。
【０２３５】
波面ｕ４ｍ（ｒ）が算出されると、ステップＳ２６２で、結像面における波面ｕｉｍｇ４ｍ（ｒ）が算出される。
【０２３６】
更に、ステップＳ２６３で、第４層の像の光強度｜ｕｉｍｇ４（ｒ）｜２が算出される。
【０２３７】
波面ｕｚｍ（ｒ）と波面ｕｉｍｇｚｍ（ｒ）は、ｍ番目の光源から射出された照明光、すなわち、１つの光源から射出された照明光によって生成される波面を表している。
【０２３８】
推定標本の像Ｉｅｓｔｚ（ｒ）は、取得位置において、全ての光源から射出された照明光によって生成される。よって、全ての光源について、波面を求めなくてはならない。
【０２３９】
（判断結果がＹＥＳの場合：ｚ＝ＮＩＭ）
ステップＳ２４２が実行される。
【０２４０】
波面ｆｉｎｍ（ｒ）、波面ｇｏｕｔｍ（ｒ）、波面ｕｚｍ（ｒ）、及び波面ｕｉｍｇｚｍ（ｒ）は、ｍ番目の光源から射出された照明光、すなわち、１つの光源から射出された照明光によって生成される波面を表している。
【０２４１】
推定標本の像Ｉｅｓｔｚ（ｒ）は、全ての光源から射出された照明光によって生成される。よって、全ての光源について、波面ｆｉｎｍ（ｒ）、波面ｇｏｕｔｍ（ｒ）、波面ｕｚｍ（ｒ）、及び波面ｕｉｍｇｚｍ（ｒ）を求めなくてはならない。
【０２４２】
ステップＳ２５２では、変数ｍの値が光源の数ＮＬＳと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ２５３が実行される。判断結果がＹＥＳの場合は、ステップＳ２７０が実行される。
【０２４３】
（判断結果がＮＯの場合：ｍ≠ＮＬＳ）
判断結果がＮＯの場合、ステップＳ２５３で、変数ｍの値に１が加算される。ステップＳ２５３が終ると、ステップＳ４１に戻る。
【０２４４】
ステップＳ２５３で、変数ｍの値が１つ増えている。そのため、別の光源について、ステップＳ４１で波面ｆｉｎｍ（ｒ）が算出され、ステップＳ４２で波面ｇｏｕｔｍ（ｒ）が算出され、ステップＳ２６１で波面ｕｚｍ（ｒ）が算出され、ステップＳ２６２で波面ｕｉｍｇｚｍ（ｒ）が算出され，ステップＳ２６３で｜ｕｉｍｇｚｍ（ｒ）｜２が算出される。
【０２４５】
ステップＳ４１、ステップＳ４２、ステップＳ２５１、及びステップＳ２６０は、全ての光源について｜ｕｉｍｇｚｍ（ｒ）｜２が求まるまで、繰り返し行われる。
【０２４６】
（判断結果がＹＥＳの場合：ｍ＝ＮＬＳ）
判断結果がＹＥＳの場合、ステップＳ２７０で、｜ｕｉｍｇｚｍ（ｒ）｜２の足し合わせが行われる。その結果、推定標本の像Ｉｅｓｔｚ（ｒ）が算出される。推定標本の像Ｉｅｓｔｚ（ｒ）は、式（３０）で表される。
【数３０】
2021024420.tif 000032
【０２４７】
図１７は、推定標本の像を示す図である。図１７（ａ）は、第１層における推定標本の像を示す図である。図１７（ｂ）は、第４層における推定標本の像を示す図である。
【０２４８】
推定標本の像Ｉｅｓｔ１（ｒ）は、全ての光源について波面ｕｉｍｇ１ｍ（ｒ）が求められた場合の像である。推定標本の像Ｉｅｓｔ４（ｒ）は、全ての光源について波面ｕｉｍｇ４ｍ（ｒ）が求められた場合の像である。
【０２４９】
図１７（ａ）に示すように、各光源について波面ｕｉｍｇ１ｍ（ｒ）が算出され、波面ｕｉｍｇ１ｍ（ｒ）から｜ｕｉｍｇ１ｍ（ｒ）｜２が算出され、｜ｕｉｍｇ１ｍ（ｒ）｜２が全て足し合わされている。その結果、第１層における推定標本の像Ｉｅｓｔ１（ｒ）が算出される。
【０２５０】
図１７（ｂ）に示すように、各光源について波面ｕｉｍｇ４ｍ（ｒ）が算出され、波面ｕｉｍｇ４ｍ（ｒ）から｜ｕｉｍｇ４ｍ（ｒ）｜２が算出され、｜ｕｉｍｇ４ｍ（ｒ）｜２が全て足し合わされている。その結果、第４層における推定標本の像Ｉｅｓｔ４（ｒ）が算出される。
【０２５１】
（判断結果がＹＥＳの場合：ｍ＝ＮＬＳ）
ステップＳ２８０で、残差が算出される。残差は、式（３１）で表される。式（３１）に示すように、残差は、撮影画像Ｉｍｅａｚ（ｒ）と推定標本の像Ｉｅｓｔｚ（ｒ）とから算出される。
【数３１】
2021024420.tif 000033
【０２５２】
上述のように、撮影画像の数は４で、推定標本の像の数も４である。よって、Ｉｍｅａ１（ｒ）とＩｅｓｔ１（ｒ）から、第１層における残差が算出される。Ｉｍｅａ２（ｒ）とＩｅｓｔ２（ｒ）から、第２層における残差が算出される。Ｉｍｅａ３（ｒ）とＩｅｓｔ３（ｒ）から、第３層における残差が算出される。Ｉｍｅａ４（ｒ）とＩｅｓｔ４（ｒ）から、第４層における残差が算出される。
【０２５３】
第１層における残差、第２層における残差、第３層における残差、及び第４層における残差から、ステップ７０で用いられる残差が算出される。
【０２５４】
ステップＳ７０では、残差と閾値との比較が行われる。判断結果がＮＯの場合は、ステップＳ２９０が実行される。判断結果がＹＥＳの場合は、ステップＳ１１０が実行される。
【０２５５】
（判断結果がＮＯの場合：残差≧閾値）
ステップＳ２９０では、変数Ｌの値が初期化される。後述のステップＳ３０１、ステップＳ３０２、ステップＳ３０３、ステップＳ３０４、及びステップＳ３１０は、ステップＳ２１０で設定した回数だけ実行される。変数Ｌは、これらのステップが実行された回数を表している。
【０２５６】
ステップＳ３００は、ステップＳ３０１、ステップＳ３０２、ステップＳ３０３、ステップＳ３０４、ステップＳ３０５、ステップＳ３０６、及びステップＳ３１０を有する。
【０２５７】
ステップＳ３０１では、１からＮＩＭまでのなかから１つがランダムに選択される。後述のステップＳ３１１では、補正後の波面が算出される。補正後の波面の算出では、１つの撮影画像と１つの推定標本の像が用いられる。
【０２５８】
上述のように、ステップＳ２７０では、複数の推定標本の像が算出される。補正後の波面の算出に用いる推定標本の像は、１つである。よって補正後の波面の算出に用いる推定標本の像を、複数の推定標本の像のなかから選択する。
【０２５９】
ＮＩＭは、層の数である。ＮＩＭ＝４の場合、ステップＳ３０１では、１から４までの数字のなかから、１つの数字がランダムに選択される。
【０２６０】
例えば、選択された数字が１の場合、数字の１は第１層を表している。第１層における推定標本の像には、１番目の取得位置における撮影画像が対応する。よって、補正後の波面の算出には、１番目の取得位置における撮影画像と、第１層における推定標本の像が用いられる。
【０２６１】
また、例えば、選択された数字が４の場合、選択された数字は第４層を表している。第４層における推定標本の像には、４番目の取得位置における撮影画像が対応する。よって、補正後の波面の算出には、４番目の取得位置における撮影画像と、第４層における推定標本の像が用いられる。
【０２６２】
ステップＳ３０２では、ステップＳ３０１で選択した値が変数ｚＬに入力される。上述のように、ステップＳ３０１では、１からＮＩＭまでの数字のなかから、１つの数字がランダムに選択される。例えば、選択された数字が１の場合、ステップＳ３０２で、変数ｚＬに１が入力されることになる。
【０２６３】
ステップＳ３０３では、変数ｍの値が初期化される。後述のステップＳ３１１、ステップＳ３１２、及びステップＳ３１３は、全ての光源に対して実行される。変数ｍは、これらのステップが実行された回数を表している。
【０２６４】
ステップＳ３１０は、推定標本の屈折率分布を最適化するステップである。
【０２６５】
ステップＳ３１０は、ステップＳ３１１、ステップＳ３１２、ステップＳ３１３、ステップＳ３１４、及びステップＳ３１５を有する。
【０２６６】
ステップＳ３１１では、波面ｕ’ｚＬｍ（ｒ）が算出される。波面ｕ’ｚＬｍ（ｒ）は、変数ｚＬの値で示される層の位置における波面である。
【０２６７】
波面ｕ’ｚＬｍ（ｒ）の算出では、撮影画像ＩｍｅａｚＬ（ｒ）と推定標本の像ＩｅｓｔｚＬ（ｒ）が用いられる。撮影画像ＩｍｅａｚＬ（ｒ）は、撮影画像Ｉｍｅａｚのうち、変数ｚＬの値で示される位置の撮像画像である。推定標本の像ＩｅｓｔｚＬ（ｒ）は、推定標本の像Ｉｅｓｔｚのうち、変数ｚＬの値で示される位置の推定標本の像である。
【０２６８】
波面ｕ’ｚＬｍ（ｒ）は式（３２）で表される。
【数３２】
2021024420.tif 000034
【０２６９】
図１８は、波面の補正を示す図である。図１８（ａ）は、推定標本から射出する補正前の波面を示す図である。図１８（ｂ）は、撮影画像の取得位置における補正前の波面を示す図である。図１８（ｃ）は、撮影画像の取得位置における補正後の波面を示す図である。図１８（ｄ）は、推定標本から射出する補正後の波面を示す図である。
【０２７０】
ステップＳ３０１で選択された数字が１の場合、すなわち、ｚＬ＝１の場合について説明する。
【０２７１】
図１７（ａ）に示すように、推定標本の像Ｉｅｓｔ１（ｒ）は、波面ｕｉｍｇ１ｍ（ｒ）に基づいて算出される。また、図１７（ａ）と図１８（ｂ）に示すように、波面ｕｉｍｇ１ｍ（ｒ）は、波面ｕ１ｍ（ｒ）に基づいて算出される。
【０２７２】
図１８（ａ）に示すように、波面ｕ１ｍ（ｒ）の算出には、振幅透過率Ｔｚ（ｒ）が用いられている。振幅透過率Ｔｚ（ｒ）は、推定された振幅透過率である。ステップＳ３００の１回目の実行時、この振幅透過率Ｔｚ（ｒ）は、標本２１の振幅透過率と異なる。
【０２７３】
振幅透過率Ｔｚ（ｒ）と標本２１の振幅透過率との差が大きくなるほど、推定標本の像Ｉｅｓｔｚ（ｒ）と撮影画像Ｉｍｅａｚ（ｒ）との差も大きくなる。よって、推定標本の像Ｉｅｓｔｚ（ｒ）と撮影画像Ｉｍｅａｚ（ｒ）との差は、振幅透過率Ｔｚ（ｒ）と標本２１の振幅透過率との差を反映していると見なすことができる。
【０２７４】
上述のように、ｚＬ＝１である。そこで、式（３２）でｚＬ＝１として、推定標本の像Ｉｅｓｔ１（ｒ）と撮影画像Ｉｍｅａ１（ｒ）と用いて、波面ｕ１ｍ（ｒ）を補正する。その結果、図１８（ｃ）に示すように、補正された波面、すなわち、波面ｕ’１ｍ（ｒ）が得られる。
【０２７５】
波面ｕ’１ｍ（ｒ）を用いることで、新たな振幅透過率を算出できる。波面ｕ’１ｍ（ｒ）は、波面ｕ１ｍ（ｒ）と異なる。よって、新たな振幅透過率は、波面ｕ１ｍ（ｒ）を算出したときの振幅透過率とは異なる。
【０２７６】
ステップＳ３１２では、補正後の波面ｇ’ｏｕｔｍ，ｚＬ（ｒ）が算出される。波面ｇ’ｏｕｔｍ，ｚＬ（ｒ）は、波面ｕ’ｚＬｍ（ｒ）がΔＤだけ伝搬したときの波面である。波面ｇ’ｏｕｔｍ，ｚＬ（ｒ）は、式（３３）で表される。
【数３３】
2021024420.tif 000035
【０２７７】
上述のように、波面ｕ’１ｍ（ｒ）を用いて、振幅透過率Ｔｚ（ｒ）を算出することができる。ただし、図１７（ａ）に示すように、振幅透過率Ｔｚ（ｒ）の算出には、波面ｇｏｕｔｍ（ｒ）の位置における波面が必要である。
【０２７８】
図１８（ａ）、（ｃ）に示すように、波面ｕ’１ｍ（ｒ）の位置と、波面ｇｏｕｔｍ（ｒ）の位置は異なる。よって、振幅透過率Ｔｚ（ｒ）を算出するためには、図１８（ｄ）に示すように、波面ｇ’ｏｕｔｍ，１（ｒ）が必要である。
【０２７９】
波面ｇ’ｏｕｔｍ，１（ｒ）は、波面ｕ’１ｍ（ｒ）が、３×Δｚだけ伝搬したときの波面である。式（３３）で、ΔＤ＝３×Δｚ、ｚＬ＝１とすることで、波面ｇ’ｏｕｔｍ，１（ｒ）を算出することができる。
【０２８０】
ステップＳ３１３では、標本の勾配ΔＴｚｍ，ｚＬ（ｒ）が算出される。ΔＴｚｍ，ｚＬ（ｒ）は、ｍ番目の光源で照明し、且つ、変数ｚＬの値で示される層の位置における撮影画像と推定標本の像で補正した時の標本の勾配である。
【０２８１】
標本の勾配ΔＴｚｍ，ｚＬは式（３４）で表される。標本の勾配ΔＴｚｍ，ｚＬ（ｒ）の算出には、例えば、勾配降下法を用いることができる。
【数３４】
2021024420.tif 000036
ここで、
ｆ＊は、ｆの複素共役
δは、ゼロ除算を防ぐための正規化定数、
である。
【０２８２】
上述のように、推定標本３０は、複数の薄い層に置き換えられている。よって、薄い層の各々について、標本の勾配ΔＴｚｍ，ｚＬ（ｒ）を算出する必要がある。
【０２８３】
図１９は、標本の勾配と波面の伝搬を示す図である。推定標本３０が４つの薄い層に置き換えられた場合、つまりＮＩＭ＝４の場合について説明する。また、ｚＬ＝１とする。図１９（ａ）は、標本の勾配を示す図である。図１９（ｂ）は、波面の伝搬を示す図である。
【０２８４】
波面ｇｏｕｔｍ（ｒ）の算出には、振幅透過率Ｔ４（ｒ）が用いられている。振幅透過率Ｔ４（ｒ）は、推定された振幅透過率である。よって、この振幅透過率Ｔ４（ｒ）は、標本２１の振幅透過率と異なる。
【０２８５】
振幅透過率Ｔ４（ｒ）と標本２１の振幅透過率との差が大きくなるほど、波面ｇｏｕｔｍ（ｒ）と波面ｇ’ｏｕｔｍ，１（ｒ）との差も大きくなる。よって、波面ｇｏｕｔｍ（ｒ）と波面ｇ’ｏｕｔｍ，１（ｒ）との差は、振幅透過率Ｔ４（ｒ）と標本２１の振幅透過率との差を反映していると見なすことができる。
【０２８６】
波面ｆ４ｍ（ｒ）、振幅透過率Ｔ４（ｒ）、波面ｇｏｕｔｍ（ｒ）、及び波面ｇ’ｏｕｔｍ，１（ｒ）は、既知である。そこで、式（３４）でｚ＝４、ｚＬ＝１とすることで、図１９（ａ）に示すように、標本の勾配ΔＴ４ｍ，１（ｒ）を算出することができる。
【０２８７】
ｇ４ｍ（ｒ）と波面ｇｏｕｔｍ（ｒ）は同じなので、ｇ４ｍ（ｒ）の代わりに、波面ｇｏｕｔｍ（ｒ）を用いれば良い。また、ｇ’４ｍ，１（ｒ）はｇ’ｏｕｔｍ，１（ｒ）と同じなので、ｇ’４ｍ，１（ｒ）の代わりに、ｇ’ｏｕｔｍ，１（ｒ）を用いれば良い。
【０２８８】
次に、標本の勾配ΔＴ３ｍ，１（ｒ）を算出する。標本の勾配ΔＴ３ｍ，１（ｒ）の算出には、波面ｇ３ｍ（ｒ）の位置における波面が必要である。この波面を算出するためには、図１９（ｂ）に示すように、波面ｆ’４ｍ，１（ｒ）が必要である。
【０２８９】
波面ｆ’４ｍ，１（ｒ）は、式（３５）でｚ＝４、ｚＬ＝１とすることで算出できる。
【数３５】
2021024420.tif 000037
【０２９０】
次に、算出された波面ｆ’４ｍ，１（ｒ）を用いて、波面ｇ３ｍ，１（ｒ）の位置における波面を算出する。
【０２９１】
図２０は、標本の勾配と波面の伝搬を示す図である。図２０（ａ）は、波面の伝搬と標本の勾配を示す図である。図２０（ｂ）は、波面の伝搬を示す図である。
【０２９２】
図２０（ａ）に示すように、波面ｇ’３ｍ，１（ｒ）は、波面ｆ’４ｍ，１（ｒ）がΔｚだけ伝搬したときの波面である。これは、第４層から第３層への波面の伝搬である。
【０２９３】
上述のように、第３層から第４層への波面の伝搬は、式（２５）で表わされる。よって、式（２５）において、以下のようにすることで、波面ｇ’３ｍ，１（ｒ）を算出できる。
波面ｆ４ｍ（ｒ）を、波面ｇ’３ｍ，１（ｒ）に置き換える。
波面ｇ３ｍ（ｒ）を、波面ｆ’４ｍ，１（ｒ）に置き換える。
ΔＤ＝−Δｚとする。
【０２９４】
波面ｆ３ｍ（ｒ）、振幅透過率Ｔ３（ｒ）、波面ｇ３ｍ（ｒ）、及び波面ｇ’３ｍ，１（ｒ）は、既知である。そこで、式（３４）でｚ＝３、ｚＬ＝１とすることで、図２０（ｂ）に示すように、標本の勾配ΔＴ３ｍ，１（ｒ）を算出することができる。
【０２９５】
波面ｆ’３ｍ，１（ｒ）は、式（３５）でｚ＝３、ｚＬ＝１とすることで算出できる。
【０２９６】
第２層と第１層についても、第３層と同様にして、標本の勾配の算出を行えばよい。
【０２９７】
図２１は、標本の勾配を示す図である。図２１では、第１層における標本の勾配ΔＴ１ｍ，１（ｒ）、第２層における標本の勾配ΔＴ２ｍ，１（ｒ）、第３層における標本の勾配ΔＴ３ｍ，１（ｒ）、及び第４層における標本の勾配ΔＴ４ｍ，１（ｒ）、が算出されている。
【０２９８】
ステップＳ３１３で得られる標本の勾配ΔＴｚｍ，１（ｒ）は、ｍ番目の光源で照明し、且つ、第１層の位置における撮影画像と第１層の位置における推定標本の像で補正した時の標本の勾配である。標本の勾配ΔＴｚｍ，１（ｒ）は、全ての光源から射出された照明光によって決まる。よって、全ての光源について、標本の勾配ΔＴｚｍ，１（ｒ）を求めなくてはならない。
【０２９９】
ステップＳ３１４では、変数ｍの値が光源の数ＮＬＳと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ３１５が実行される。判断結果がＹＥＳの場合は、ステップＳ３０４が実行される。
【０３００】
（判断結果がＮＯの場合：ｍ≠ＮＬＳ）
判断結果がＮＯの場合、ステップＳ３１５で、変数ｍの値に１が加算される。ステップＳ３１５が終ると、ステップＳ３１１に戻る。
【０３０１】
ステップＳ３１５で、変数ｍの値が１つ増えている。そのため、別の光源について、ステップＳ３１１で波面ｕｚ’ｍ，１（ｒ）が算出され、ステップＳ３１２で波面ｇｏｕｔ’ｍ，１（ｒ）が算出され、ステップＳ３１３で標本の勾配ΔＴｚｍ，１（ｒ）が算出される。
【０３０２】
ステップＳ３１１、ステップＳ３１２、及びステップＳ３１３は、全ての光源について標本の勾配ΔＴｚｍ，１（ｒ）が求まるまで、繰り返し行われる。
【０３０３】
図２２は、標本の勾配を示す図である。図２２では、全ての光源について、標本の勾配ΔＴｚｍ．１（ｒ）が求められている。
【０３０４】
（判断結果がＹＥＳの場合：ｍ＝ＮＬＳ）
判断結果がＹＥＳの場合、ステップＳ３０４で、振幅透過率Ｔｚ（ｒ）が更新される。ステップＳ３０４は、推定標本を更新するステップするステップである。
【０３０５】
更新された振幅透過率Ｔｚ（ｒ）は、式（３６）で表される。
【数３６】
2021024420.tif 000038
ここで、
αは、標本の勾配の補正係数、
である。
【０３０６】
ステップＳ３０５では、変数Ｌの値が補正の回数ＮＣＲと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ３０６が実行される。判断結果がＹＥＳの場合は、ステップＳ３０が実行される。
【０３０７】
（判断結果がＮＯの場合：ｍ≠ＮＣＲ）
判断結果がＮＯの場合、ステップＳ３０６で、変数Ｌの値に１が加算される。ステップＳ３０６が終ると、ステップＳ３０１に戻る。
【０３０８】
ステップＳ３０１で、１からＮＩＭまでのなかから１つがランダムに選択される。選択された数字に基づいて、補正に用いる推定標本の像と取得位置が決まる。
【０３０９】
そして、ステップＳ３１１で波面ｕｚ’ｍ，１（ｒ）が算出され、ステップＳ３１２で波面ｇｏｕｔ’ｍ，１（ｒ）が算出され、ステップＳ３１３で標本の勾配ΔＴｚｍ，１（ｒ）が算出され、ステップＳ３０４で振幅透過率Ｔｚ（ｒ）が更新される。
【０３１０】
ステップＳ３０１、ステップＳ３０２、ステップＳ３０３、ステップＳ３０４、及びステップＳ３１０は、設定された回数の補正が終るまで、繰り返し行われる。
【０３１１】
（判断結果がＹＥＳの場合：ｍ＝ＮＣＲ）
判断結果がＹＥＳの場合、ステップＳ３０に戻る。更新された振幅透過率Ｔｚ（ｒ）で、ステップＳ３０からステップＳ３００までが実行される。
【０３１２】
ステップＳ３０からステップＳ３００までが繰り返し実行されることで、更新された振幅透過率Ｔｓ（ｒ）は、徐々に、標本２１の振幅透過率に近づく。すなわち、残差が小さくなる。やがて、残差は閾値よりも小さくなる。
【０３１３】
（判断結果がＹＥＳの場合：残差＜閾値）
ステップＳ１１０では、推定標本の屈折率分布が算出される。得られた振幅透過率Ｔｚ（ｒ）は、標本２１の振幅透過率と同一か、又は、略同一である。得られた振幅透過率Ｔｚ（ｒ）と式（１）から、屈折率分布ｎｚ（ｒ）が求まる。
【０３１４】
ステップＳ１１０で得られた屈折率分布ｎｚ（ｒ）を用いることで、推定標本の構造を再構成することができる。再構成された推定標本の構造は、例えば、表示装置に出力することができる。推定標本３０は、厚い標本である。第２のシミュレーションでは、厚い標本の構造は、推定標本の３次元構成を再構成することができる。
【０３１５】
上述のように、ステップＳ１１０で得られた振幅透過率Ｔｚ（ｒ）は、標本２１の振幅透過率と同一か、又は、略同一である。この場合、屈折率分布ｎｚ（ｒ）も、標本２１の屈折率分布と同一か、又は、略同一と見なすことができる。よって、再構成された推定標本３０の構造は、標本６の構造と同一か、又は、略同一と見なすことができる。
【０３１６】
第２のシミュレーションでは、ステップＳ２５０、ステップＳ２７０、及びステップＳ３１０が、繰り返し実行される。その結果、振幅透過率Ｔｚ（ｒ）が更新される。上述のように、ステップＳ２５０とステップＳ２７０は、推定標本の像を算出するステップである。ステップＳ３１０は、推定標本の屈折率分布を最適化するステップである。
【０３１７】
振幅透過率Ｔｚ（ｒ）は、推定標本を表している。よって、推定標本の像を算出するステップと推定標本の屈折率分布を最適化するステップが繰り返し実行されることで、推定標本が更新される。
【０３１８】
本実施形態の屈折率分布推定システムは、光源を有し、照明光学系は、コンデンサレンズと、第１開口部材と、を有し、結像光学系は、対物レンズと、結像レンズと、を有し、対物レンズの瞳位置に、第１開口部材の像が結像されることが好ましい。
【０３１９】
図２３は、本実施形態の屈折率分布推定システムを示す図である。図１（ａ）と同じ構成については同じ番号を付し、説明は省略する。
【０３２０】
屈折率分布推定システム４０は、光源４１を有する。照明光学系２は、コンデンサレンズ４２と、第１開口部材４３と、レンズ４４と、を有する。結像光学系３は、対物レンズ４５と、結像レンズ４６と、を有する。
【０３２１】
光源４１から射出された光は、レンズ４４で集光される。集光位置に、第１開口部材４３が配置されている。第１開口部材４３は、コンデンサレンズ４２の前側焦点位置に配置されている。よって、コンデンサレンズ４２から、平行光束が標本６に射出される。
【０３２２】
標本６は、対物レンズ４５の前側焦点位置に配置されている。標本６から対物レンズ４５に入射した光は、対物レンズ４５から平行光束となって射出される。対物レンズ４５から射出された平行光束は、結像レンズ４６で集光される。
【０３２３】
対物レンズ４５の前側焦点位置は、コンデンサレンズ４２の後側焦点位置と一致している。よって、対物レンズ４５の瞳位置Ｐｕに、第１開口部材４３の像が結像成される。照明光学系２と結像光学系３には、例えば、顕微鏡の光学系を用いることができる。
【０３２４】
本実施形態の屈折率分布推定システムでは、第１開口部材は、輪帯形状の透過部又は減光部を有することが好ましい。
【０３２５】
図２４は、第１開口部材を示す図である。図２４（ａ）は、第１例を示す図である。図２４（ｂ）は、第２例を示す図である。
【０３２６】
第１例では、第１開口部材５０は、輪帯形状の透過部５１と、円形の遮光部５２と、輪帯形状の遮光部５３と、を有する。透過部５１は、減光部にしても良い。
【０３２７】
第２例では、第２開口部材６０は、輪帯形状の透過部６１と、円形の遮光部６２と、輪帯形状の遮光部６３と、を有する。透過部６１は、減光部にしても良い。
【０３２８】
第１開口部材５０と第２開口部材６０では、透過部の位置と、透過部の幅が異なる。透過部５１は、透過部６１よりも外側に位置している。透過部５１の幅は、透過部６１の幅よりも狭い。
【０３２９】
本実施形態の屈折率分布推定システムは、第１開口部材とは異なる第２開口部材と、第１開口部材と第２開口部材とを切り替える移動機構と、を備えることが好ましい。
【０３３０】
図２５は、本実施形態の屈折率分布推定システムを示す図である。図９（ａ）と同じ構成については同じ番号を付し、説明は省略する。
【０３３１】
屈折率分布推定システム７０は、複数の開口部材を有する。例えば、第１開口部材５０と、第２開口部材６０と、を使用することができる。第１開口部材５０と第２開口部材６０は、移動機構７１に配置することができる。
【０３３２】
移動機構７１として、例えば、スライダーを用いることができる。この場合、スライダーを移動することで、第１開口部材５０と第２開口部材６０のどちらか一方を、光軸ＡＸ上に位置させることができる。
【０３３３】
移動機構７１として、例えば、ターレットを用いることができる。この場合、ターレットを回転することで、第１開口部材５０と第２開口部材６０のどちらか一方を、光軸ＡＸ上に位置させることができる。
【０３３４】
図２６は、第３のシミュレーションのフローチャートである。第２のフローチャートにおける処置と同じ処理につては、同じ番号を付し、説明は省略する。
【０３３５】
ステップＳ４００では、開口部材の数ＮＡＰが設定される。上述のように、屈折率分布推定システム７０では、複数の開口部材が使用される。
【０３３６】
ステップＳ４１０では、変数ｎの値が初期化される。後述のステップＳ２３１、ステップＳ５０１、ステップＳ５０２、ステップＳ７０、ステップＳ２９０、及びステップＳ６００は、少なくとも１つの開口部材に対して実行される。変数ｎは、これらのステップが実行された回数を表している。
【０３３７】
ステップＳ５００は、ステップＳ２３１、ステップＳ５０１、ステップＳ５０２、ステップＳ７０、ステップＳ２９０、ステップＳ５０３、ステップＳ５０４、ステップＳ５０５、及びステップＳ６００を有する。
【０３３８】
ステップＳ２３１では、推定標本３０における振幅透過率Ｔｚ（ｒ）に初期値が設定される。振幅透過率Ｔｚ（ｒ）は、ｚ番目の層における振幅透過率である。
【０３３９】
図２６では、ステップＳ２３１しか図示されていない。しかしながら、振幅透過率Ｔｚ（ｒ）を算出する具体的な処理は、第２のシミュレーションにおける振幅透過率Ｔｚ（ｒ）の算出処理と同じである。よって、振幅透過率Ｔｚ（ｒ）の数は、層の数と同じである。
【０３４０】
ステップＳ５０１では、推定標本の像Ｉｅｓｔｚｎ（ｒ）が算出される。推定標本の像Ｉｅｓｔｚｎ（ｒ）は、ｎ番目の開口部材を用いたときのｚ番目の層における推定標本の像である。
【０３４１】
図２６では、ステップＳ５０１しか図示されていない。しかしながら、推定標本の像Ｉｅｓｔｚｎ（ｒ）を算出する具体的な処理は、第２のシミュレーションにおける推定標本の像Ｉｅｓｔｚ（ｒ）の算出処理と同じである。よって、推定標本の像Ｉｅｓｔｚｎ（ｒ）の数は、層の数と同じである。
【０３４２】
推定標本の像Ｉｅｓｔｚｎ（ｒ）の算出では、以下の式（３７）〜式（４４）が用いられる。
【数３７】
2021024420.tif 000039
【数３８】
2021024420.tif 000040
【数３９】
2021024420.tif 000041
【数４０】
2021024420.tif 000042
ここで、
ΔＤは、隣り合う２つの層の間隔、
である。
【数４１】
2021024420.tif 000043
【数４２】
2021024420.tif 000044
ここで、
△Ｄは、推定標本の表面から撮影画像の取得位置までの距離、
である。
【数４３】
2021024420.tif 000045
【数４４】
2021024420.tif 000046
【０３４３】
ステップＳ７０では、残差と閾値との比較が行われる。判断結果がＮＯの場合は、ステップＳ２９０が実行される。判断結果がＹＥＳの場合は、ステップＳ１１０が実行される。
【０３４４】
（判断結果がＮＯの場合：残差≧閾値）
ステップＳ２９０で変数Ｌの値が初期化された後、ステップＳ６００が実行される。
【０３４５】
ステップＳ６０１で、標本の勾配ΔＴｚｎ，ｍ，ｚＬ（ｒ）が算出される。標本の勾配ΔＴｚｎ，ｍ，ｚＬ（ｒ）は、ｎ番目の開口部材を用いたときのｚ番目の層における標本の勾配である。
【０３４６】
標本の勾配ΔＴｚｎ，ｍ，ｚＬ（ｒ）では、ｍ番目の光源で標本が照明されている。また、ｚＬ番目の取得位置における撮影画像と、ｚＬ番目の層における推定標本の像を用いて、補正が行われている。
【０３４７】
標本の勾配ΔＴｚｎ，ｍ，ｚＬ（ｒ）の算出では、以下の式（４５）〜式（４９）が用いられる。
【数４５】
2021024420.tif 000047
【数４６】
2021024420.tif 000048
ここで、
△Ｄは、ＺＬ番目の層から推定標本の表面までの距離、
である。
【数４７】
2021024420.tif 000049
【数４８】
2021024420.tif 000050
【数４９】
2021024420.tif 000051
【０３４８】
ステップＳ６０１で得られる標本の勾配ΔＴｚｎ，ｍ，ＺＬ（ｒ）は、ｎ番目の開口部材を用い、ｍ番目の光源で照明し、且つ、ｚＬ番目の位置における撮影画像とｚＬ番目の位置における推定標本の像で補正した時の標本の勾配である。標本の勾配ΔＴｚｎ，ｍ，ＺＬ（ｒ）は、全ての光源から射出された照明光によって決まる。よって、全ての光源について、標本の勾配ΔＴｚｎ，ｍ，ｚＬ（ｒ）を求めなくてはならない。
【０３４９】
ステップＳ６０２では、変数ｍの値が光源の数ＮＬＳと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ６０３が実行される。判断結果がＹＥＳの場合は、ステップＳ６０４が実行される。
【０３５０】
（判断結果がＮＯの場合：ｍ≠ＮＬＳ）
判断結果がＮＯの場合、ステップＳ６０３で、変数ｍの値に１が加算される。ステップＳ６０３が終ると、ステップＳ３０１に戻る。
【０３５１】
ステップＳ６０３で、変数ｍの値が１つ増えている。そのため、別の光源について、ステップＳ６０１で標本の勾配Ｔｚｎ，ｍ，ｚＬ（ｒ）が算出される。
【０３５２】
（判断結果がＹＥＳの場合：ｍ＝ＮＬＳ）
判断結果がＹＥＳの場合、ステップＳ６０４で、振幅透過率Ｔｚ（ｒ）が更新される。ステップＳ６０４は、推定標本を更新するステップするステップである。
【０３５３】
更新された振幅透過率Ｔｚ（ｒ）は、式（５０）で表される。
【数５０】
2021024420.tif 000052
ここで、
αは、標本の勾配の補正係数、
である。
【０３５４】
ステップＳ５０３では、変数ｎの値が開口部材の数ＮＡＰと一致したか否かが判断される。判断結果がＮＯの場合は、ステップＳ５０４が実行される。判断結果がＹＥＳの場合は、ステップＳ５０５が実行される。
【０３５５】
（判断結果がＮＯの場合：ｎ≠ＮＡＰ）
判断結果がＮＯの場合、ステップＳ５０４で、変数ｎの値に１が加算される。ステップＳ５０４が終ると、ステップＳ２３１に戻る。
【０３５６】
ステップＳ５０４で、変数ｎの値が１つ増えている。そのため、別の開口部材について、ステップＳ５００が実行される。
【０３５７】
（判断結果がＹＥＳの場合：ｎ＝ＮＡＰ）
判断結果がＹＥＳの場合、ステップＳ５０５で、変数ｎの値に１が設定される。ステップＳ５０５が終ると、ステップＳ２３１に戻る。
【０３５８】
全ての開口部材についてステップＳ５００を実行しても、残差が閾値よりも大きい場合がある。このような場合は、再び、最初の開口部材を用いて、ステップＳ５００が実行される。
【０３５９】
（判断結果がＹＥＳの場合：残差＜閾値）
ステップＳ１１０では、推定標本の屈折率分布が算出される。得られた振幅透過率Ｔｚ（ｒ）は、標本２１の振幅透過率と同一か、又は、略同一である。得られた振幅透過率Ｔｚ（ｒ）と式（１５）から、屈折率分布ｎｚ（ｒ）が求まる。
【０３６０】
図２７は、開口部材と撮影画像を示す図である。図２７（ａ）は、第１の開口部材と撮影画像を示す図である。図２７（ｂ）は、第２の開口部材と撮影画像を示す図である。図２７（ｃ）は、第３の開口部材と撮影画像を示す図である。
【０３６１】
第１の開口部材、第２の開口部材及び第３の開口部材は、輪帯形状の透過部を有する。透過部の位置、透過部の幅は、各々の開口部材で異なる。
【０３６２】
撮影画像は、ＸＺ面における光学像の画像である。Ｘ方向は、光軸と直交する方向である。Ｚ方向は光軸方向である。よって、撮影画像では、線状の画像が光軸方向に積み重なっている。
【０３６３】
上述のシミュレーションでは、推定標本における振幅透過率Ｔｚ（ｒ）に、初期値が設定されている。第２のシミュレーションで説明したように、初期値の設定に、強度輸送方程式を用いることができる。
【０３６４】
図２８は、標本、撮影画像、及び初期値の画像を示す図である。図２８（ａ）は、標本を示す図である。図２８（ｂ）は、撮影画像を示す図である。図２８（ｃ）は、初期値の画像を示す図である。
【０３６５】
標本は、フォトニッククリスタルファイバー（以下、「ＰＣＦ」という）である。撮影画像は、ＸＺ面における光学像の画像である。初期値を示す画像は、式（１５）で算出した値に基づく画像である。
【０３６６】
図２９は、伝達可能な空間周波数帯域を示す図である。図２９（ａ）は、透過部の形状が円形の場合の空間周波数帯域を示す図である。図２９（ｂ）は、透過部の形状が輪帯形状の場合の空間周波数帯域を示す図である。
【０３６７】
２つの図において、横軸はＺ方向、縦軸はＸ方向、又はＹ方向である。Ｙ方向は、Ｚ方向とＸ方向の両方と直交する方向である。また、空間周波数帯域の算出では、開口数が１．４の対物レンズを用い、照明光学系の開口数を変化させている。
【０３６８】
標本の像が再生されるためには、標本の持つ空間周波数が結像面まで伝達されなくてはならない。伝達可能な空間周波数が決まる要因として、照明光学系の開口数と、結像光学系の開口数がある。
【０３６９】
透過部の形状が円形の場合、図２９（ａ）に示すように、破線で示す範囲の内側には、空間周波数帯域が存在しない。これに対して、透過部の形状が輪帯形状の場合、図２９（ｂ）に示すように、破線で示す範囲の内側には、空間周波数帯域が存在する。
【０３７０】
破線で示す範囲は、空間周波数が低い範囲である。よって、透過部の形状が輪帯形状の場合、低い空間周波数から高い空間周波数までを、伝達することができる。その結果、標本の像を、より正確に再生することができる。
【０３７１】
輪帯形状の位置、輪帯形状の幅によって、伝達可能な空間周波数は変化する。
【０３７２】
図３０は、標本を示す図である。図３０（ａ）は、標本の全体を示す断面図である。図３０（ｂ）は、標本の一部を示す断面図である。
【０３７３】
標本は、ＰＣＦである。ＰＣＦでは、クラッドの直径は１８０．０μｍで、穴の直径は６．０μｍで、穴の間隔は１２．９μｍである。クラッドの屈折率は１．４５６である。
【０３７４】
ＰＣＦは液体に浸されている。よって、穴は液体で満たされている。また、ＰＣＦの外側も、同じ液体で覆われている。液体の屈折率は１．４３６である。クラッドの屈折率は、クラッド以外の部分の屈折率よりも０．０２だけ高い。
【０３７５】
照明光の波長は１５００ｎｍである。また、対物レンズの開口数は１．４である。
【０３７６】
図３１は、撮影画像、伝達可能な空間周波数を示す図、再生像の全体の画像、及び再生像の一部の画像である。図３１（ａ）は、第１輪帯における撮影画像である。図３１（ｂ）は、第２輪帯における撮影画像である。図３１（ｃ）は、第３輪帯における撮影画像である。図３１（ｄ）は、第４輪帯における撮影画像である。
【０３７７】
輪帯の位置と輪帯の幅は、照明光学系の開口数で表すことができる。上述のように、４つの輪帯では、以下のようになっている。数値は、照明光学系の開口数である。左側の数値は、輪帯の内側を表している。右側の数値は、輪帯の外側を表している。
【０３７８】
第１輪帯：０．１４−０．５４。
第２輪帯：０．７−０．９。
第３輪帯：１．０−１．１。
第４輪帯：１．１５−１．２５。
【０３７９】
照明光学系の開口数の数値が大きくなるほど、輪帯は外側に位置する。よって、第１輪帯が最も内側に位置し、第４輪帯が最も外側に位置している。
【０３８０】
図３１（ａ）、図３１（ｂ）、図３１（ｃ）、図３１（ｄ）は、ＸＺ面における光学像の画像である。
【０３８１】
図３１（ａ）、図３１（ｂ）、図３１（ｃ）、及び図３１（ｄ）を比較すると、第１輪帯、第２輪帯、第３輪帯、第４輪帯の順で、画像が鮮明になっている。すなわち、輪帯が外側に位置するほど、画像が鮮明になっている。
【０３８２】
画像は、伝達可能な空間周波数帯域に含まれる空間周波数が高くなるほど鮮明になる。上述のように、輪帯が外側に位置するほど、画像が鮮明になっている。よって、輪帯が外側に位置するほど、伝達可能な空間周波数帯域に含まれる空間周波数が高くなる。
【０３８３】
このように、輪帯の位置と伝達可能な空間周波数には、相関がある。そこで、複数の輪帯を組み合わせることで、伝達可能な空間周波数帯域に含まれる空間周波数を変えることができる。
【０３８４】
図３１（ｅ）、図３１（ｆ）、図３１（ｇ）、及び図３１（ｈ）では、伝達可能な空間周波数を示し、輪帯の数と組み合わせは、以下のようになっている。
図３１（ｅ）：第１輪帯、第２輪帯、第３輪帯、第４輪帯。
図３１（ｆ）：第１輪帯、第４輪帯。
図３１（ｇ）：第４輪帯。
図３１（ｈ）：第１輪帯。
【０３８５】
図３１（ｅ）、図３１（ｆ）、図３１（ｇ）、及び図３１（ｈ）を比較すると、輪帯が外側に位置するほど、伝達可能な空間周波数に含まれる空間周波数は高くなる。また、輪帯の組み合わせ数が多くなるほど、伝達可能な空間周波数に含まれる空間周波数が多くなる。
【０３８６】
標本の像の再生では、撮影画像が用いられる。標本の像をより正確に再生するためには、撮影画像に、低い空間周波数から高い空間周波数までが含まれていると良い。
【０３８７】
図３１（ｉ）、図３１（ｊ）、図３１（ｋ）、図３１（ｌ）、図３１（ｍ）、図３１（ｎ）、図３１（ｏ）、及び図３１（ｐ）では、再生像を示し、輪帯の数と組み合わせは、以下のようになっている。
図３１（ｉ）、（ｍ）：第１輪帯、第２輪帯、第３輪帯、第４輪帯。
図３１（ｊ）、（ｎ）：第１輪帯、第４輪帯。
図３１（ｋ）、（ｏ）：第４輪帯。
図３１（ｌ）、（ｐ）：第１輪帯。
【０３８８】
図３１（ｉ）、図３１（ｊ）、図３１（ｋ）、及び図３１（ｌ）を比較すると、輪帯が外側に位置するほど、再生像が鮮明になっている。また、輪帯の組み合わせ数が多くなるほど、再生像が鮮明になっている。
【０３８９】
４つの輪帯を用いた場合、２つの輪帯を用いた場合に比べて、ＰＣＦの外側のアーチファクトが少ない。そのため、図３１（ｍ）と図３１（ｎ）から分かるように、４つの輪帯を用いた方が、実際の穴の形状に近くなっている。
【０３９０】
穴の径は、ＰＣＦの外径よりも小さい、穴の方が、ＰＣＦの外形（輪郭）よりも、高い空間周波数を有する。よって、第４輪帯を用いると、再生像において、高い空間周波数の部分を、鮮明に再生することができる。そのため、図３１（ｋ）、（ｏ）に示すように、穴の再構成はできている。
【０３９１】
しかしながら、第４輪帯を用いると、再生像において、低い空間周波数の部分については情報が欠落している。そのため、図３１（ｋ）に示すように、ＰＣＦの外形が上下に伸びていている。
【０３９２】
また、第１輪帯を用いると、図３１（ｌ）に示すようにＰＣＦの外形が、概ね再生できている。しかしながら、第１輪帯を用いると、再生像において、高い低い空間周波数の部分については情報が欠落している。そのため、図３１（ｐ）に示すように、穴の再構成ができていない。
【０３９３】
このように、複数の輪帯を組み合わせることで、像の再構成に利用できる空間周波数帯域を広げることができる。その結果、像の再構成の性能を向上させることができる。
【０３９４】
図３２は、標本を示す図と再生像の画像である。図３２（ａ）は、第１標本を示す図である。図３２（ｂ）は、第２標本を示す図である。図３２（ｃ）は、第１標本の再生像の画像である。図３２（ｄ）は、第２標本の再生像の画像である。
【０３９５】
第１標本は、５０μｍの細胞の塊である。第２標本は、１００μｍの細胞の塊である。照明光の波長は、１．５μｍである。１つの細胞の大きさは、１０μｍである。核の大きさは３μｍである。細胞質の屈折率は、１．３６である。 核の屈折率は、１．４である。溶液の屈折率１．３３である。
【０３９６】
図３２（ｃ）と図３２（ｄ）に示すように、第１標本と第２標本のいずれについても、鮮明な再生像が得られている。
【０３９７】
図３３は、標本を示す図、初期値の画像、及び再生像の画像である。図３３（ａ）は、標本を示す図である。図３３（ｂ）は、初期値の画像である。図３３（ｃ）は、標本の再生像の画像である。
【０３９８】
標本は、５０μｍの細胞塊である。図３３（ｃ）に示すように、３次元の標本でも鮮明で正確な再生像をえることができている。
【０３９９】
本実施形態の屈折率分布推定システムでは、推定標本像を算出するステップは、照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を計算し、複数の第１波面が推定標本を通過した後の複数の第２波面を計算し、複数の第２波面を用いて結像光学系の結像位置における複数の第１強度分布を計算し、複数の第１強度分布を足し合わせることにより、推定標本の像を算出し、標本の屈折率分布を最適化するステップは、複数の光源のそれぞれについて、複数の第２波面を撮影画像及び推定標本の像を用いて補正した複数の第２補正波面を算出し、複数の第２波面と複数の第２補正波面との誤差から推定標本の屈折率分布の勾配を算出し、屈折率分布の勾配を用いて推定標本の屈折率分布を最適化することが好ましい。
【０４００】
第１波面の計算は、ステップＳ４１で行われる。第２波面の計算は、ステップＳ４２で行われる。第１強度分布の計算は、ステップＳ４５、又はステップＳ２６３で行われる。推定標本の像の算出は、ステップＳ５０、又はステップＳ２７０で行われる。
【０４０１】
標本の屈折率分布を最適化するステップは、ステップＳ９０、又はステップＳ３１０である。第２補正波面の算出は、ステップＳ９１、又はステップＳ３１１で行われる。推定標本の屈折率分布の勾配の算出は、ステップＳ９２、又はステップＳ３１２で行われる。
【０４０２】
本実施形態の屈折率分布推定システムでは、推定標本像を算出するステップは、照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を計算するステップと、複数の第１波面が推定標本を通過した後の複数の第２波面を計算するステップと、複数の第２波面から結像光学系の標本側の焦点位置における複数の第３波面を計算するステップと、複数の第３波面と結像光学系の瞳関数とを用いて結像光学系の結像位置における複数の第４波面を計算し、複数の第４波面をそれぞれ２乗して複数の第１強度分布を計算するステップと、複数の第１強度分布を足し合わせることにより、推定標本像を算出するステップと、を有することが好ましい。
【０４０３】
第１波面を計算するステップは、ステップＳ４１である。第２波面を計算するステップは、ステップＳ４２である。第３波面を計算するステップは、ステップＳ４３、又はステップＳ２６２である。
【０４０４】
第１強度分布を計算するステップは、ステップＳ４４とステップＳ４５、又はステップＳ２６２とステップＳ２６３である。推定標本像を算出するステップは、ステップＳ５０、又はステップＳ２７０である。
【０４０５】
本実施形態の屈折率分布推定システムでは、標本の屈折率分布を最適化するステップは、複数の光源のそれぞれについて、複数の第２波面を撮影画像及び推定標本の像を用いて補正した複数の第２補正波面を算出し、複数の第２波面と複数の第２補正波面との誤差から推定標本の屈折率分布の勾配を算出し、屈折率分布の勾配を用いて推定標本の屈折率分布を最適化することが好ましい。
【０４０６】
本実施形態の屈折率分布推定システムでは、ＴＶ正則化を用いて振幅透過率を更新することが好ましい。
【０４０７】
ＴＶ正則化は、エッジを保持しつつ、振動成分を抑制できる。そのため、ノイズ除去やぼけ画像の修正などの画像処理で用いられる。
【０４０８】
振幅透過率Ｔｓ（ｒ）又は振幅透過率Ｔｚ（ｒ）を、式（５１）に示される最小化問題によって、Ｔ’ｓ（ｒ）又は振幅透過率Ｔ’ｚ（ｒ）に更新する。
【数５１】
2021024420.tif 000053
【０４０９】
右辺の第１項は、推定残差のＬ２ノルムを示す。右辺の第２項は、正則化項と呼ばれ、振幅透過率Ｔｓ（ｒ）又は振幅透過率Ｔｚ（ｒ）の局所変化が少なければ少ないほど、小さい値をとる性質の関数が一般的に用いられる。τは正則化パラメータと呼ばれる定数である。
【０４１０】
この正則化項として、式（５２）に示す隣接する画素間の推定値の差分の絶対値和を意味するＴＶ正則化項を加えると、エッジを保持したまま平滑化することができる。
【数５２】
2021024420.tif 000054
【０４１１】
図３４は、標本を示す図と再生像の画像である。図３４（ａ）は、標本を示す図である。図３４（ｂ）、（ｃ）は撮影画像である。図３４（ｄ）、（ｅ）、（ｆ）、（ｇ）、（ｈ）は、標本の再生像の画像である。
【０４１２】
標本は円筒である。円筒の直径は１０μｍである。円筒の屈折率は１．３６である。円筒は液体に浸されている。液体の屈折率は１．３３である。つまり、円筒の屈折率は、円筒以外の部分の屈折率よりも０．０３だけ高い。
【０４１３】
照明光の波長は１５００ｎｍである。また、対物レンズの開口数は１．４である。
【０４１４】
図３４（ｂ）、（ｃ）は、ＸＺ面における光学像の画像である。図３４（ｂ）の撮影画像では、照明光学系の開口数は１．０である。図３４（ｃ）の撮影画像では、照明光学系の開口数は０．５−１．０である。
【０４１５】
照明条件、初期値の有無、ＴＶ正則化の有無と、各図の関係を以下に示す。
【表１】
2021024420.tif 000055
【０４１６】
図３４（ｄ）と図３４（ｆ）の比較から、初期値を設定すると、より正確に像を再生することができることが分かる。
【０４１７】
図３４（ｄ）と図３４（ｅ）の比較から、ＴＶ正則化を行うと、エッジを保持しつつ振動成分を抑制できるため、ノイズ除去され自然な画像に再構成できることが分かる。
【０４１８】
図３４（ｇ）と図３４（ｈ）の比較から、輪帯を用いた照明だと、より正確に像を再生することができることが分かる。
【０４１９】
本発明には、以下の発明が含まれる。
（付記項１）
複数の方向から同時に入射する光線により標本を照明する照明光学系と、
標本の光学像を形成する結像光学系と、
結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、
撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、
プロセッサは、
標本の屈折率分布を含む推定標本を推定するステップと、
照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を用いて、結像光学系の結像位置における複数の第１強度分布を計算し、複数の第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、
複数の第１波面が推定標本を通過した後の複数の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、
推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、
更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする屈折率分布推定システム。
【０４２０】
（付記項２）
照明光学系は、空間的パーシャルコヒーレントな光で標本を照明することを特徴とする付記項１に記載の屈折率分布推定システム。
【産業上の利用可能性】
【０４２１】
以上のように、本発明は、光軸と直交する面内における分解能高く、標本以外の情報の影響が少ない屈折率分布推定システムに適している。
【符号の説明】
【０４２２】
１  屈折率分布推定システム
２  照明光学系
３  結像光学系
４  撮像素子
５  プロセッサ
６  標本
６ａ  表面
６’  光学像
１０  推定標本
１０ａ  表面
２０  屈折率分布推定システム
２１  標本
２１ａ  表面
２１’  光学像
２２  移動ステージ
３０  推定標本
３０ａ、３０ｂ  表面
４０  屈折率分布推定システム
４１  光源  
４２  コンデンサレンズ
４３  第１開口部材
４４  レンズ
４５  対物レンズ
４６  結像レンズ
５０  第１開口部材
６０  第２開口部材
５１、６１  透過部
５２、５３、６２、６３  遮光部
７０  屈折率分布推定システム
７１  移動機構
ＡＸ  光軸
ＬＯＮ、ＬＯＦＦ  光線
Ｆｏ  結像光学系の焦点位置
Ｐｕ  結像光学系の瞳位置（対物レンズの瞳位置）
ＩＬＬ  照明領域
ＰＯＢ  結像光学系の瞳
ＩＰ  結像面
【図１】
2021024420.tif 000056
【図２】
2021024420.tif 000057
【図３】
2021024420.tif 000058
【図４】
2021024420.tif 000059
【図５】
2021024420.tif 000060
【図６】
2021024420.tif 000061
【図７】
2021024420.tif 000062
【図８】
2021024420.tif 000063
【図９】
2021024420.tif 000064
【図１０】
2021024420.tif 000065
【図１１】
2021024420.tif 000066
【図１２】
2021024420.tif 000067
【図１３】
2021024420.tif 000068
【図１４】
2021024420.tif 000069
【図１５】
2021024420.tif 000070
【図１６】
2021024420.tif 000071
【図１７】
2021024420.tif 000072
【図１８】
2021024420.tif 000073
【図１９】
2021024420.tif 000074
【図２０】
2021024420.tif 000075
【図２１】
2021024420.tif 000076
【図２２】
2021024420.tif 000077
【図２３】
2021024420.tif 000078
【図２４】
2021024420.tif 000079
【図２５】
2021024420.tif 000080
【図２６】
2021024420.tif 000081
【図２７】
2021024420.tif 000082
【図２８】
2021024420.tif 000083
【図２９】
2021024420.tif 000084
【図３０】
2021024420.tif 000085
【図３１】
2021024420.tif 000086
【図３２】
2021024420.tif 000087
【図３３】
2021024420.tif 000088
【図３４】
2021024420.tif 000089
【手続補正書】
【提出日】令和3年8月2日(2021.8.2)
【手続補正１】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項１】
パーシャルコヒーレント照明により標本を照明する照明光学系と、
前記標本の光学像を形成する結像光学系と、
前記結像光学系により形成された前記標本の光学像から撮影画像を取得する撮像素子と、
前記撮影画像から前記標本の屈折率分布を再構成するプロセッサと、を備え、
前記プロセッサは、
前記標本の屈折率分布を含む推定標本を推定するステップと、
前記照明光学系の瞳の強度分布をモデル化した複数の光源のそれぞれに対応する第１波面を用いて、前記複数の光源のそれぞれに対応する前記結像光学系の結像位置における第１強度分布を計算し、前記複数の光源のそれぞれに対応する前記第１強度分布を足し合わせることにより、前記推定標本の像を算出するステップと、
前記複数の光源のそれぞれに対応する前記第１波面のそれぞれが前記推定標本を通過した後の第２波面、前記撮影画像及び前記推定標本の像を用いて前記推定標本の屈折率分布を最適化するステップと、
前記推定標本の像の算出と前記推定標本の屈折率分布の最適化とを繰り返して、前記推定標本を更新するステップと、
更新された前記推定標本の屈折率分布を用いて前記推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする屈折率分布推定システム。
【請求項２】
前記結像光学系の焦点位置と前記標本の位置との間隔を、前記結像光学系の光軸方向に変化させる駆動機構を更に備え、前記駆動機構により前記間隔を変化させて、複数の間隔のそれぞれに対応する複数の標本の撮影画像を取得し、
前記複数の間隔のそれぞれに対応する複数の推定標本の像を算出し、
前記複数の間隔のそれぞれにおいて、前記推定標本の屈折率分布を最適化することを特徴とする請求項１に記載の屈折率分布推定システム。
【請求項３】
光源を有し、
前記照明光学系は、コンデンサレンズと、第１開口部材と、を有し、
前記結像光学系は、対物レンズと、結像レンズと、を有し、
前記対物レンズの瞳位置に、前記第１開口部材の像が結像されることを特徴とする請求項１に記載の屈折率分布推定システム。
【請求項４】
前記第１開口部材は、輪帯形状の透過部又は減光部を有することを特徴とする請求項３に記載の屈折率分布推定システム。
【請求項５】
前記第１開口部材とは異なる第２開口部材と、
前記第１開口部材と前記第２開口部材とを切り替える移動機構と、を備えることを特徴とする請求項３に記載の屈折率分布推定システム。
【請求項６】
前記推定標本像を算出するステップは、
前記照明光学系の瞳の強度分布をモデル化した複数の光源から射出された複数の第１波面を計算し、
前記複数の第１波面が推定標本を通過した後の複数の第２波面を計算し、
前記複数の第２波面を用いて前記結像光学系の結像位置における複数の第１強度分布を計算し、前記複数の第１強度分布を足し合わせることにより、前記推定標本の像を算出し、
前記標本の屈折率分布を最適化するステップは、
前記複数の光源のそれぞれに対応する前記第２波面のそれぞれを処理対象として、前記第２波面と同じ前記光源に対応する前記撮影画像及び前記推定標本の像を用いて補正した第２補正波面を算出し、
前記複数の光源のそれぞれに対応する前記第２波面のそれぞれを処理対象として、同じ前記光源に対応する前記第２波面と前記第２補正波面との誤差から前記推定標本の屈折率分布の勾配を算出し、前記屈折率分布の勾配を用いて前記推定標本の屈折率分布を最適化することを特徴とする請求項１に記載の屈折率分布推定システム。
【請求項７】
前記推定標本像を算出するステップは、
前記複数の光源のそれぞれに対応する前記第１波面を計算するステップと、
前記第１波面のそれぞれに対応する、前記推定標本を通過した後の前記第２波面を計算するステップと、
前記第２波面のそれぞれに対応する、前記結像光学系の前記標本側の焦点位置における第３波面を計算するステップと、
前記第３波面のそれぞれと前記結像光学系の瞳関数とを用いて、前記第３波面のそれぞれに対応する、前記結像光学系の結像位置における第４波面を計算し、前記第４波面のそれぞれを２乗して前記第４波面のそれぞれに対応する第１強度分布を計算するステップと、
前記第４波面のそれぞれに対応する前記第１強度分布を足し合わせることにより、前記推定標本像を算出するステップと、を有することを特徴とする請求項１に記載の屈折率分布推定システム。
【請求項８】
前記標本の屈折率分布を最適化するステップは、
前記複数の光源のそれぞれに対応する前記第２波面のそれぞれを処理対象として、前記第２波面と同じ前記光源に対応する前記撮影画像及び前記推定標本の像を用いて補正した第２補正波面を算出し、
前記複数の光源のそれぞれに対応する前記第２波面のそれぞれを処理対象として、同じ前記光源に対応する前記第２波面と前記第２補正波面との誤差から前記推定標本の屈折率分布の勾配を算出し、前記屈折率分布の勾配を用いて前記推定標本の屈折率分布を最適化することを特徴とする請求項７に記載の屈折率分布推定システム。
【請求項９】
複数の方向から同時に入射する光線により標本を照明する照明光学系と、
前記標本の光学像を形成する結像光学系と、
前記結像光学系により形成された前記標本の光学像から撮影画像を取得する撮像素子と、
前記撮影画像から前記標本の屈折率分布を再構成するプロセッサと、を備え、
前記プロセッサは、
前記標本の屈折率分布を含む推定標本を推定するステップと、
前記照明光学系の瞳の強度分布をモデル化した複数の光源のそれぞれに対応する第１波面を用いて、前記複数の光源のそれぞれに対応する前記結像光学系の結像位置における第１強度分布を計算し、前記複数の光源のそれぞれに対応する前記第１強度分布を足し合わせることにより、前記推定標本の像を算出するステップと、
前記複数の光源のそれぞれに対応する前記第１波面のそれぞれが前記推定標本を通過した後の第２波面、前記撮影画像及び前記推定標本の像を用いて前記推定標本の屈折率分布を最適化するステップと、
前記推定標本の像の算出と前記推定標本の屈折率分布の最適化とを繰り返して、前記推定標本を更新するステップと、
更新された前記推定標本の屈折率分布を用いて前記推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする屈折率分布推定システム。
【請求項１０】
前記照明光学系は、空間的パーシャルコヒーレントな光で前記標本を照明することを特徴とする請求項９に記載の屈折率分布推定システム。
【手続補正２】
【補正対象書類名】明細書
【補正対象項目名】００２４
【補正方法】変更
【補正の内容】
【００２４】
上述した課題を解決し、目的を達成するために、本発明の少なくとも幾つかの実施形態に係る屈折率分布推定システムは、
パーシャルコヒーレント照明により標本を照明する照明光学系と、
標本の光学像を形成する結像光学系と、
結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、
撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、
プロセッサは、
標本の屈折率分布を含む推定標本を推定するステップと、
照明光学系の瞳の強度分布をモデル化した複数の光源のそれぞれに対応する第１波面を用いて、複数の光源のそれぞれに対応する結像光学系の結像位置における第１強度分布を計算し、複数の光源のそれぞれに対応する第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、
複数の光源のそれぞれに対応する第１波面のそれぞれが推定標本を通過した後の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、
推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、
更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする。
【手続補正３】
【補正対象書類名】明細書
【補正対象項目名】００２５
【補正方法】変更
【補正の内容】
【００２５】
本発明の少なくとも幾つかの実施形態に係る屈折率分布推定システムは、
複数の方向から同時に入射する光線により標本を照明する照明光学系と、
標本の光学像を形成する結像光学系と、
結像光学系により形成された標本の光学像から撮影画像を取得する撮像素子と、
撮影画像から標本の屈折率分布を再構成するプロセッサと、を備え、
プロセッサは、
標本の屈折率分布を含む推定標本を推定するステップと、
照明光学系の瞳の強度分布をモデル化した複数の光源のそれぞれに対応する第１波面を用いて、複数の光源のそれぞれに対応する結像光学系の結像位置における第１強度分布を計算し、複数の光源のそれぞれに対応する第１強度分布を足し合わせることにより、推定標本の像を算出するステップと、
複数の光源のそれぞれに対応する第１波面のそれぞれが推定標本を通過した後の第２波面、撮影画像及び推定標本の像を用いて推定標本の屈折率分布を最適化するステップと、
推定標本の像の算出と推定標本の屈折率分布の最適化とを繰り返して、推定標本を更新するステップと、
更新された推定標本の屈折率分布を用いて推定標本の構造を再構成して出力するステップと、を含む処理を行うことを特徴とする。
【手続補正４】
【補正対象書類名】明細書
【補正対象項目名】００４４
【補正方法】変更
【補正の内容】
【００４４】
結像面ＩＰには、撮像素子４の撮像面が位置している。撮像素子４によって、光学像６’の画像の取得が行われる。その結果、図１（ｂ）に示す撮影画像Ｉｍｅａ（ｒ）が得られる。ｒは、（ｘ，ｙ）の２次元座標を示す。
【手続補正５】
【補正対象書類名】明細書
【補正対象項目名】００５４
【補正方法】変更
【補正の内容】
【００５４】
図５では、１番目の光源からＮＬＳ番目までの光源が図示されている。光源は、照明光学系２の瞳位置に配置することができる。このようにすることで、照明光学系２の瞳の強度分布を、複数の光源でモデル化することができる。
【手続補正６】
【補正対象書類名】明細書
【補正対象項目名】０２７７
【補正方法】変更
【補正の内容】
【０２７７】
上述のように、波面ｕ’１ｍ（ｒ）を用いて、振幅透過率Ｔｚ（ｒ）を算出することができる。ただし、図１８（ａ）に示すように、振幅透過率Ｔｚ（ｒ）の算出には、波面ｇｏｕｔｍ（ｒ）の位置における波面が必要である。
【国際調査報告】
2021024420.tif 000090
2021024420.tif 000091
2021024420.tif 000092
2021024420.tif 000093