(57)【要約】 【課題】  強度輸送方程式をフーリエ級数展開して伝播光の位相を求める際にコンピュータメモリの増大を抑制する方法を提供する。【解決手段】  伝播距離を変えて伝播光の強度を複数種類取得し、前記複数種類の強度と伝播距離の情報から強度輸送方程式を算出した後、強度輸送方程式をフーリエ級数展開し、伝播光の位相の空間周波数成分を求める際に、特定の範囲だけで位相の空間周波数成分を求め、こうして求めた位相の空間周波数成分を逆フーリエ変換して伝播光の位相を算出する。【選択図】  図８
2019200080.tif 000002
【特許請求の範囲】
【請求項１】
伝播する光に対して、伝播距離を変えて強度をディテクタで複数計測し、前記複数の強度データから強度輸送方程式を導出したのち、前記強度輸送方程式のフーリエ変換から求めた伝播光の空間周波数成分を逆フーリエ変換することで伝播光の位相を算出する方法において、
前記位相に対して、初期位相を設定する初期位相設定ステップと、
前記位相の空間周波数成分を求める範囲を設定する空間周波数範囲設定ステップと、
前記空間周波数範囲内で伝播光の位相の空間周波数成分を求める空間周波数成分算出ステップと、
前記伝播光の位相の空間周波数成分を逆フーリエ変換して位相を算出する位相算出ステップと、
前記算出位相と前記初期位相を比較し、両者の位相がほぼ等しいと判断されれば、前記算出位相を伝播光の位相とし、そうでなければ、前記算出位相を前記初期位相に設定し直し、前記空間周波数範囲を包含する空間周波数範囲を設定したのち空間周波数成分算出ステップに戻るステップと、
を有することを特徴とする位相算出方法。
【請求項２】
伝播する光に対して、伝播距離を変えて強度をディテクタで複数計測し、前記複数の強度データから導出した強度輸送方程式をフーリエ級数展開法で解くことで伝播光の位相を算出する方法において、
前記位相に対して、初期位相を設定する初期位相設定ステップと、
前記位相の空間周波数成分を求める範囲を設定する空間周波数範囲設定ステップと、
前記空間周波数範囲内で伝播光の位相の空間周波数成分を求める空間周波数成分算出ステップと、
前記伝播光の位相の空間周波数成分を逆フーリエ変換して位相を算出する位相算出ステップと、
前記算出位相と前記初期位相を比較し、両者の位相がほぼ等しいと判断されれば、前記算出位相を伝播光の位相とし、そうでなければ、前記算出位相を前記初期位相に設定し直し、前記空間周波数範囲を包含する空間周波数範囲を設定したのち空間周波数成分算出ステップに戻るステップと、
を有することを特徴とする位相算出方法。
【請求項３】
請求項１または２に記載の位相算出方法において、第一、及び第二の空間周波数範囲は矩形であることを特徴とする位相算出方法。
【請求項４】
請求項１または２に記載の位相算出方法において、第一、及び第二の空間周波数範囲は円形であることを特徴とする位相算出方法。
【請求項５】
請求項１または２に記載の位相算出方法において、第一、及び第二の空間周波数範囲の最大値は、強度データの最大空間周波数の半分の値であることを特徴とする位相算出方法。
【請求項６】
結像光学系の光学像をフォーカスを変えてディテクタで複数計測し、前記複数の光学像データから強度輸送方程式を導出したのち、前記強度輸送方程式のフーリエ変換から求めた光学像の空間周波数成分を逆フーリエ変換することで光学像の位相を算出する方法において、
結像光学系の瞳を基準にして光学像の空間周波数成分を求め、前記光学像の空間周波数成分を逆フーリエ変換して光学像の位相を算出することを特徴とする位相算出方法。
【請求項７】
結像光学系の光学像をフォーカスを変えてディテクタで複数計測し、前記複数の光学像データから導出した強度輸送方程式をフーリエ級数展開法で解くことで光学像の位相を算出する方法において、
結像光学系の瞳を基準にして光学像の空間周波数成分を求め、前記光学像の空間周波数成分を逆フーリエ変換して光学像の位相を算出することを特徴とする位相算出方法。
【請求項８】
請求項６または７に記載の位相算出方法において、瞳の２倍の大きさの円形領域内で光学像の空間周波数成分を求めることを特徴とする位相算出方法。
【請求項９】
請求項６または７に記載の位相算出方法において、瞳の２倍の大きさの円に外接する矩形領域内で光学像の空間周波数成分を求めることを特徴とする位相算出方法。
【請求項１０】
請求項１乃至４のいずれか１項に記載の位相計測方法を実行する計算機に接続された位相測定器。
【請求項１１】
請求項６乃至９のいずれか１項に記載の位相計測方法を実行する計算機に接続された顕微鏡。
【請求項１２】
請求項１１に記載の顕微鏡において、照明のコヒーレンスファクターが０．３以下であることを特徴とした顕微鏡。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、伝播光の位相を算出する方法、及び、位相を算出する装置に関する。
【背景技術】
【０００２】
光は振幅と位相で記述される。位相は非常に高速に振動するので、通常のディテクタでは振動を追いきれず、我々が観測できるのはある時間内での積算強度である。そこで、光の位相を測定するには、干渉計を用いるのが一般的である。
【０００３】
一方、強度輸送方程式（Ｔｒａｎｓｐｏｒｔ  ｏｆ  Ｉｎｔｅｎｓｉｔｙ  Ｅｑｕａｔｉｏｎ；  ＴＩＥ）は伝播光の強度と位相を関係付けている。強度輸送方程式に関しては特許文献１に詳しく記載されているが、強度輸送方程式を解けば伝播光の位相を算出することができる。強度輸送方程式に入力すべきデータは、伝播距離を変えて取得した複数の強度データから得ることができる。そのため、比較的簡素な測定で強度輸送方程式を得ることができる。
【０００４】
しかし、強度輸送方程式は数学的に高度な式で、強度輸送方程式を解くのは容易ではない。特に、特許文献１記載の強度輸送方程式の解法には二種類の問題がある。一番目の問題は、強度データの大部分において強度が０、もしくは、強度が非常に弱い領域で占められている場合、強度輸送方程式の解が不安定になるという問題である。強度が０だと位相が不定になるからである。二番目の問題は、取得した強度データ境界に位相分布があるとき、強度輸送方程式の解の精度が低下するという問題である。
【０００５】
第一の問題を避けるためには背景が明るいことが好ましく、第二の問題を避けるためには位相分布がデータ中央に孤立していることが好ましい。それでもなお、第一の問題を避けきれず、観測強度に強度がほぼ０の領域が含まれてしまう場合がある。強度がある閾値（例えば、最大強度の１％）以下の領域が非常に小さい（例えば面積で０．１％以下）ならば、特許文献１に記述されているように、前記領域の値を適当な値（例えば、前記閾値と同じ値）に変えればよい。
【０００６】
もし、強度がある閾値以下の領域がデータの０．１％以上を占めると、特許文献１記載の方法では強度輸送方程式の解の精度が低下する恐れがある。このような場合は、非特許文献１に記載されている強度輸送方程式のフーリエ級数展開法を用いればよいことがわかっている。強度輸送方程式のフーリエ級数展開法では、伝播光の位相の空間周波数成分を求めることができる。位相の空間周波数成分を逆フーリエ変換すれば伝播光の位相を求めることができる。
【０００７】
しかし、非特許文献１記載の強度輸送方程式のフーリエ級数展開法では膨大なコンピュータメモリを必要とする。具体的には、強度データがＭ×Ｎピクセルであるときは、ＭＮ行ＭＮ列の行列の逆行列を求める必要がある。例えば、５１２×５１２ピクセルの強度データが与えられたときは５１２２行５１２２列の行列の逆行列を求めなければならない。これは倍精度の変数を用いれば約２７５ＧＢに相当するデータ量であり、通常のコンピュータでは計算が困難である。よって、非特許文献１記載の強度輸送方程式のフーリエ級数展開法は理論的には有効な方法であるが、実用的には実装が困難であった。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】米国特許第６，８８５，４４２号
【非特許文献】
【０００９】
【非特許文献１】Ｔ．  Ｅ．  Ｇｕｒｅｙｅｖ  ａｎｄ  Ｋ．  Ａ．  Ｎｕｇｅｎｔ著，  「Ｐｈａｓｅ  ｒｅｔｒｉｅｖａｌ  ｗｉｔｈ  ｔｈｅ  ｔｒａｎｓｐｏｒｔ−ｏｆ−ｉｎｔｅｎｓｉｔｙ  ｅｑｕａｔｉｏｎ．  ＩＩ．  Ｏｒｔｈｏｇｏｎａｌ  ｓｅｒｉｅｓ  ｓｏｌｕｔｉｏｎ  ｆｏｒ  ｎｏｎｕｎｉｆｏｒｍ  ｉｌｌｕｍｉｎａｔｉｏｎ」  Ｊｏｕｒｎａｌ  ｏｆ  Ｏｐｔｉｃａｌ  Ｓｏｃｉｅｔｙ  ｏｆ  Ａｍｅｒｉｃａ  Ａ．Ｖｏｌ．  １３，  アメリカ合衆国，１９９６年，ｐｐ．１６７０−１６８２
【発明の概要】
【発明が解決しようとする課題】
【００１０】
強度輸送方程式をフーリエ級数展開して伝播光の位相を求める際、コンピュータメモリの増大が実用上の問題になっており、データが非常に小さいときしか解が得られなかった。
【課題を解決するための手段】
【００１１】
強度輸送方程式をフーリエ級数展開し、伝播光の位相の空間周波数成分を求める際に、特定の範囲だけで位相の空間周波数成分を求める。こうして求めた位相の空間周波数成分を逆フーリエ変換して伝播光の位相を算出する。
【００１２】
本発明の更なる目的又はその他の特徴は、以下、添付図面を参照して説明される好ましい実施例によって明らかにされるであろう。
【発明の効果】
【００１３】
本発明によれば、強度輸送方程式をフーリエ級数展開して伝播光の位相を求める際に問題となっていたコンピュータメモリの増大を抑制することができる。
【図面の簡単な説明】
【００１４】
【図１】本発明の一側面としてのコンピュータ構成を示す図である。
【図２】本発明の構成を示したフローチャートである。
【図３】本発明の別の側面を示したフローチャートである。
【図４】本発明の実施の形態で強度輸送方程式を解く際に使用する強度データである。（ａ）と（ｂ）では伝播距離に１ｍｍの差がある。
【図５】（ａ）光伝播シミュレーションから得られた伝播光の位相。（ｂ）本発明の解法で求めた伝播光の位相。（ｃ）（ａ）と（ｂ）の差分。
【図６】本発明の実施の形態で強度輸送方程式を解く際に使用する強度データで、顕微鏡の光学像である。（ａ）と（ｂ）では伝播距離に１ｍｍの差がある。
【図７】（ａ）シミュレーションで得られた光学像の位相。（ｂ）数２を解くことで得られた光学像の位相。（ｃ）本発明の解法で得られた光学像の位相。
【図８】本発明の一側面としての位相計測装置の概略ブロック図である。
【発明を実施するための形態】
【００１５】
本発明は、強度輸送方程式を用いた伝播光の位相算出方法、及び、位相算出装置に関する。特に、光学系を介して得られる光学像の位相を強度輸送方程式から算出する方法、及び、装置に関する。
【００１６】
本発明で開示される概念は、数学的にモデル化することができる。従って、本発明は、コンピュータ・システムのソフトウエア機能として実装することができる。コンピュータ・システムのソフトウエア機能は、実行可能なソフトウエア・コードを有するプログラミングを含み、本実施形態では、強度輸送方程式を解く。ソフトウエア・コードは、コンピュータ・システムのプロセッサによって実行される。ソフトウエア・コード動作中において、コード又は関連データ記録は、コンピュータ・プラットフォームに格納される。
【００１７】
但し、ソフトウエア・コードは、他の場所に格納される、或いは、適切なコンピュータ・システムにロードされることもある。従って、ソフトウエア・コードは、１つ又は複数のモジュールとして、少なくとも１つのコンピュータ読み取り可能な記録媒体で保持することができる。本発明は、上述したコードという形式で記述することが可能であり、１つ又は複数のソフトウエア製品として機能させることができる。
【００１８】
本明細書で用いる光学系について説明する。座標系は二種類ある。第一の座標系は、光強度が観測される面の座標であり、本明細書では（ｘ，ｙ）で記述する。第二の座標系は、座標系（ｘ，ｙ）とは互いにフーリエ変換の関係にある座標系で、本明細書では（ｆ，ｇ）で記述する。（ｆ，ｇ）は空間周波数とも呼ばれる。光軸をｚ軸とする。
強度輸送方程式について説明する。光の強度をＩ、光の位相をφとする。光の波長をλ、光の波数ｋ＝λ／（２π）とする。便宜上、∇⊥＝（∂／∂ｘ，∂／∂ｙ）という演算子を使うと、強度輸送方程式は以下の数１で与えられる。
【００１９】
【数１】
2019200080.tif 000003
【００２０】
なお、数１中の点（ドット）はベクトルの内積を表す。ｚ＝ｚ１で取得した強度データをＩ１（ｘ，ｙ）、ｚ＝ｚ２で取得した強度データをＩ２（ｘ，ｙ）とすれば、∂Ｉ／∂ｚ＝（Ｉ１−Ｉ２）／（ｚ１−ｚ２）、Ｉ＝（Ｉ１＋Ｉ２）／２と近似できる。よって、伝播距離が異なる位置で強度データを最低２種類取得すれば強度輸送方程式を得ることができる。
【００２１】
以下、数１の解法を説明する。まず、強度データはＭ行Ｎ列のピクセルからなり、行方向のピクセル位置をｍ、列方向のピクセル位置をｎで表す。フーリエ変換と逆フーリエ変換をそれぞれＦＴとＦＴ−１とする。特許文献１記載の解法では、特別な条件下でしか成立しない近似を使っている。近年では、特許文献１記載の解法を発展させ、数２のように繰り返しループで強度輸送方程式を解く方法が提案されている。
【００２２】
【数２】
2019200080.tif 000004
【００２３】
ただし、∇―２ＦＦＴは、以下の数３で与えられる。
【００２４】
【数３】
2019200080.tif 000005
【００２５】
数２では、ｊ＝１からスタートし、通常初期位相φ０は初期位相であるが通常は０とする。繰り返しループ回数は通常１０回以下である。
【００２６】
数２を見ると、１／Ｉという演算が必要となる。そのため、強度Ｉが０に近い値を持つと、１／Ｉが大きくなる。その後、∇−２ＦＦＴ∇⊥という演算を行うことによって、１／Ｉの影響が増幅され解が低下する。しかも、その低下した解に対して繰り返しループを行うのでますます解が低下する。
【００２７】
一方、非特許文献１では強度輸送方程式の異なる解法を提案している。強度Ｉのフーリエ変換を
2019200080.tif 000006
【００２８】
、位相φのフーリエ変換を
2019200080.tif 000007
【００２９】
、数１の右辺のフーリエ変換を
2019200080.tif 000008
【００３０】
とする。数値計算ではデータは離散化されている。ｐを０以上Ｍ未満の整数、ｑを０以上Ｎ未満の整数とすれば、ｐ行ｑ列目のデータを下付き添え字ｐ，ｑで表すとする。さて、数１の両辺をフーリエ変換すれば数４を得ることができる。
【００３１】
【数４】
2019200080.tif 000009
【００３２】
数４は数１をフーリエ級数展開することと等価で、強度輸送方程式のフーリエ級数展開法と呼ばれる。数４を解けば位相φのフーリエ変換
2019200080.tif 000010
【００３３】
を求めることができる。位相φを求めるには
2019200080.tif 000011
【００３４】
を逆フーリエ変換すればよい。
【００３５】
数４には１／Ｉという演算が必要ない。よって、強度データに強度が０に近い領域が数％含まれていても精度が低下することなく強度輸送方程式を解くことができる。しかし、
2019200080.tif 000012
【００３６】
のピクセル数はＭＮであるので、ＭＮ行ＭＮ列の行列に対して逆行列を求める必要があり、膨大なコンピュータメモリが必要となる。
【００３７】
そこで、発明者は、コンピュータメモリを低減して数４を解く方法を考えることにした。もっとも単純かつ効果的な方法として、
2019200080.tif 000013
【００３８】
を算出する空間周波数範囲を限定する方法を思いついた。計算しうる全ての空間周波数において
2019200080.tif 000014
【００３９】
を算出するから、ＭＮ行ＭＮ列の行列に対して逆行列を求める必要が生じる。しかし、空間周波数を限定して
2019200080.tif 000015
【００４０】
を算出すれば、行列サイズを低減することができる。問題は、どの範囲まで空間周波数を限定して
2019200080.tif 000016
【００４１】
を算出するかであるが、この問題を議論する前に、本発明に係る強度輸送方程式を解くためのコンピュータの構成とフローチャートを説明する。
本発明の実施形態に係る強度輸送方程式を解くためのコンピュータの構成を、図１を用いて説明する。
【００４２】
コンピュータ１は、バス配線１０、制御部２０、表示部３０、記憶部４０、信号入出力部５０、入力部６０及び媒体インターフェース７０を備える。
【００４３】
制御部２０、表示部３０、記憶部４０、信号入出力部５０、入力部６０及び媒体インターフェース７０は、バス配線１０を介して相互に接続されている。媒体インターフェース７０は、記録媒体８０を接続可能に構成されている。
【００４４】
記憶部４０には、波長情報４０ａ、伝播距離情報４０ｂ、強度情報４０ｃ、ピクセルサイズ情報４０ｄ、ＮＡ情報４０ｅ、及び強度輸送方程式ソルバー４０ｆが記憶されている。波長情報４０ａは、伝播光の波長である。伝播距離情報４０ｂは、複数の強度データを取得した相対的な伝播距離である。強度情報４０ｃは、伝播距離に対応した伝播光の強度である。ピクセルサイズ情報４０ｄは、ディテクタのピクセルサイズである。ＮＡ情報４０ｅは、光学系の像側開口数である。強度輸送方程式ソルバー４０ｆは、強度輸送方程式を解くためのプログラムである。
【００４５】
制御部２０は、例えば、ＣＰＵ，ＧＰＵ，ＤＳＰ又はマイコンなどであり、一時記憶のためのキャッシュメモリをさらに含む。表示部３０は、例えば、ＣＲＴディスプレイや液晶ディスプレイなどの表示デバイスである。記憶部４０は、例えば、メモリやハードディスクなどである。入力部６０は、例えば、キーボードやマウスなどである。媒体インターフェース７０は、例えば、ＣＤ−ＲＯＭドライブやＵＳＢインターフェースなどである。記録媒体８０は、ＣＤ−ＲＯＭやＵＳＢメモリなどである。  次に、強度輸送方程式ソルバーを実行して伝播光の位相を算出する処理の流れを、図２のフローチャートを用いて説明する。
【００４６】
ステップＳ１０１では、制御部２０は、波長情報４０ａ、伝播距離４０ｂ、強度情報４０ｃ、ピクセルサイズ情報４０ｄ、そして、ＮＡ情報４０ｅをロードし、計算条件を設定する。
【００４７】
また、強度輸送方程式ソルバー４０ｅが記録された記録媒体８０が、媒体インターフェース７０に接続される。そして、強度輸送方程式ソルバー４０ｆは、インストールされ、制御部２０を介して記憶部４０に記憶される。
【００４８】
入力部６０には、利用者により、強度輸送方程式ソルバー４０ｆの起動命令が入力される。制御部２０は、強度輸送方程式ソルバー４０ｆの起動命令を受け取り、その起動命令に基づいて、記憶部４０を参照し、強度輸送方程式ソルバー４０ｆを起動する。
【００４９】
ステップＳ１０２では、制御部２０は、
2019200080.tif 000017
【００５０】
を算出する空間周波数範囲を設定する。
【００５１】
すなわち、強度輸送方程式ソルバー４０ｆがステップ１で設定された計算条件から
2019200080.tif 000018
【００５２】
を算出する空間周波数範囲を設定する。利用者が
2019200080.tif 000019
【００５３】
を算出する空間周波数範囲が設定することも可能である。
【００５４】
ステップＳ１０３では、コンピュータ１の制御部２０は、
2019200080.tif 000020
【００５５】
を算出する空間周波数範囲を限定して数４を解く。
【００５６】
すなわち、強度輸送方程式ソルバー４０ｆは、波長情報４０ａ、伝播距離４０ｂ、強度情報４０ｃ、そして、ピクセルサイズ情報４０ｄを参照し、かつ、ステップＳ１０２で設定された空間周波数範囲に限定して数４を解く。結果として得られた
2019200080.tif 000021
【００５７】
を記憶部４０に記憶される。
【００５８】
ステップＳ１０４では、制御部２０は、ステップＳ１０３で得られた
2019200080.tif 000022
【００５９】
をフーリエ変換する。
【００６０】
すなわち、制御部２０は、ステップＳ１０３で得られた
2019200080.tif 000023
【００６１】
をフーリエ変換し、最終的なアウトプットとして位相φを算出する。位相φは表示部３０に表示され、利用者が結果を確認することができる。
【００６２】
以上が本実施形態の強度輸送方程式ソルバー４０ｆを用いて伝播光の強度データから伝播光の位相を回復する処理である。本発明の特徴は図２記載のフローチャートのステップＳ１０２と、ステップＳ１０２の結果を受けて実行されるステップＳ１０３である。特に重要なのはステップＳ１０２である。Ｓ１０１で設定された計算条件から、
2019200080.tif 000024
【００６３】
が取りうる最大空間周波数は一意に決まる。ステップＳ１０２では、
2019200080.tif 000025
【００６４】
をどの範囲まで、もしくは、どの領域内で算出するかを決定する。すなわち、求めるべき
2019200080.tif 000026
【００６５】
の空間周波数領域を設定する。ステップＳ１０２では、利用者が
2019200080.tif 000027
【００６６】
の空間周波数領域を設定しても良いが、制御部２０がステップＳ１０１で設定された計算条件から
2019200080.tif 000028
【００６７】
の空間周波数領域を限定することができる。以下で、
2019200080.tif 000029
【００６８】
の空間周波数領域を限定する特徴的な３種類の方法を説明する。
【００６９】
第一の方法は、トライアンドエラーである。自然画像には低空間周波数成分が支配的であることは一般的に良く知られている。そこで、（ｆ２＋ｇ２）１／２≦ｒの範囲だけで
2019200080.tif 000030
【００７０】
を算出し、そのフーリエ変換としてφを求める。ｒを最初小さい値に設定し、徐々に大きくしていきφを求める。ｒが大きくなるにしたがってφは一定の値に近づいていくので、φが一定の値に近づくまでｒを大きくしていけばよい。この手続きは実際には図２のフローチャートというよりは、図３のフローチャートで実行されるべきである。
【００７１】
図３記載のフローチャートは図２記載のフローチャートから内容が理解できるので、詳細説明は省略する。なお、ステップＳ２０８では、φｊとφｊ−１の差分のＬ２ノルムがある閾値以下かどうかを判定している。当該ステップＳ２０８は、φｊとφｊ−１が十分類似しているかどうかを判定するステップで、その目的が達成できるのであればＬ２ノルム以外の評価量を使ってよい。ｊｍａｘはＭの半分、もしくは、Ｎの半分程度の値でよい。また、δｒはトライアンドエラーの繰り返しの中で、空間周波数の範囲をどれくらい大きくすべきかを表す量で、１ピクセルに相当する大きさであることが好ましい。ステップＳ２０５では、（ｆ２＋ｇ２）１／２≦ｒｊの範囲で
2019200080.tif 000031
【００７２】
を算出しているが、実際には、｜ｆ｜＜ｒｊ、｜ｇ｜＜ｒｊという正方領域で
2019200080.tif 000032
【００７３】
を算出してもよい。
【００７４】
第二の方法は、数値計算の原理に基づく方法である。数４には、
2019200080.tif 000033
【００７５】
という値がでてくるが、０以上ｐ−ｍはＭ未満である必要があり、ｑ−ｎも０以上Ｎ未満である必要がある。そのためには、ｐはＭ／２以下であるのが好ましく、同様にｑもＮ／２以下であることが好ましい。その結果、最大空間周波数の半分の範囲まで
2019200080.tif 000034
【００７６】
を算出すればよい。
【００７７】
第三の方法は、ハードウエアに限定を加える方法である。強度輸送方程式は伝播光の位相を算出することができる。しかし、単純な伝播光ではなく、結像光学系を介して得られた光学像の位相を強度輸送方程式で算出する状況だけを考える。実用的な例は、明視野顕微鏡でほぼ透明で孤立した位相物体の位相を観察する場合である。
光学像は光学系の瞳を通った光によって形成される。光学系の射出瞳（単に瞳ともいう）は像面とフーリエ変換の関係にあることに着目する。瞳の半径を１に規格化すれば、瞳とは（ｆ２＋ｇ２）１／２≦１の領域と置き換えることができる。瞳内を通った電場が像面に到達し、その絶対値の２乗として光学像（強度）Ｉが得られるわけである。よって、光学像Ｉのフーリエ変換
2019200080.tif 000035
【００７８】
は、瞳関数の自己相関で表され、瞳の２倍の領域内、すなわち（ｆ２＋ｇ２）１／２≦２でしか値を持たない。同じ原理で、
2019200080.tif 000036
【００７９】
も（ｆ２＋ｇ２）１／２≦２でしか値を持たない。逆に言えば、（ｆ２＋ｇ２）１／２＞２ならば
2019200080.tif 000037
【００８０】
である。よって、（ｆ２＋ｇ２）１／２≦２の範囲内だけで
2019200080.tif 000038
【００８１】
を求め、それ以外では
2019200080.tif 000039
【００８２】
とするのは良い近似であることが予想できる。もし、（ｆ２＋ｇ２）１／２≦２の範囲内だけでしか
2019200080.tif 000040
【００８３】
を求めないのであれば、瞳面座標全体で
2019200080.tif 000041
【００８４】
を求める必要がなくなるので、数値計算負荷を低減することが可能となる。なお、数値計算的には、｜ｆ｜＜２、｜ｇ｜＜２の正方領域内で
2019200080.tif 000042
【００８５】
を求めるほうが圧倒的に楽である。  以下の実施例では、本発明の撮像装置、及び、撮像方法、本発明の効果、更なる応用方法について、図を用いて詳細に説明する。
［実施例１］
シミュレーションにより本発明の効果を確認する。光源の波長を６３３ｎｍとする。ピクセルサイズが１９．５μｍのＣＣＤをある位置に配置し取得した強度分布は図４（a）のようになった。さらにＣＣＤを光の進行方向に１ｍｍずらして取得した強度分布は図４（ｂ）のようになった。以上の情報をもとに、制御部２０は、波長情報４０ａを６３３ｎｍ、伝播距離４０ｂを１ｍｍ、強度情報４０ｃを図４（ａ）と（ｂ）、そして、ピクセルサイズ情報４０ｄを１９．５μｍに設定する。光学系を介して光学像を得たわけではないので、ＮＡ情報４０ｅは設定されない。ここで、強度輸送方程式ソルバーを起動し、伝播光の位相回復を行う。
【００８６】
図４（ａ）と（ｂ）の強度データには５１２×５１２ピクセルあるので、非特許文献１に従って数４を解くには５１２２行５１２２列の行列の逆行列が必要となるが、発明者のコンピュータではこの計算はメモリ不足のため無理であった。
【００８７】
そこで、図３記載のフローチャートに従い、位相φを算出することにした。δｒは１ピクセルに相当するようにセットした。また、数値計算の容易さから、｜ｆ｜＜ｒｊ、｜ｇ｜＜ｒｊという正方領域で
2019200080.tif 000043
【００８８】
を算出した。実際に計算を行うと、ｊ＝１５で解がほとんど変化しなくなったので、繰り返しループを終了した。
【００８９】
シミュレーションにより伝播光の位相を求めると、図５（ａ）のようになった。一方、図３記載のフローチャートに従い強度輸送方程式を解いた結果が、図５（ｂ）である。図５（ａ）と（ｂ）の差分が図５（ｃ）で、両者にはほとんど差がないことがわかる。すなわち、十分な精度で解が得られていることがわかる。
【００９０】
本実施例では、ｊ＝１５で繰り返し計算を終えたので、（２ｊ＋１）２行（２ｊ＋１）２列、すなわち、９６１行９６１列の行列の逆行列を解けばよい。理論的には、５１２２行５１２２列の行列の逆行列を求める必要があったことを考えると、行列の要素数を３１４／５１２４＝１．３×１０−５まで減らしていることになる。つまり、要素数を約０．００１％まで低減している。それにもかかわらず位相回復誤差は約０．１％であり、十分な精度で位相が回復できていることがわかる。
【００９１】
［実施例２］
デジタル明視野顕微鏡を考える。光源の波長は６５０ｎｍで、コヒーレンス度の高い光源を用いていると仮定する。顕微鏡の像側開口数は０．７５／８０である。ピクセルサイズは４．２μｍである。略ベストフォーカスの位置で光学像を取得し（図６（ａ））、その位置から１μｍデフォーカスした位置で別の光学像を取得した（図６（ｂ））。
【００９２】
以上の情報をもとに、制御部２０は、波長情報４０ａを６５０ｎｍ、伝播距離４０ｂを１μｍ、強度情報４０ｃを図６（ａ）と（ｂ）、ピクセルサイズ情報４０ｄを４．２μｍ、そして、ＮＡ情報４０ｅを０．７５／８０と設定する。ここで、強度輸送方程式ソルバーを起動し、光学像の位相を算出する。
【００９３】
まず、ピクセルサイズの逆数をピクセル数で割り、その値に波長をかけさらにＮＡで割ると瞳面でピクセルサイズを求めることができる。こうして求めた瞳座標は瞳半径が１に規格化されている。本実施例では、瞳面でのピクセルサイズは０．０６４５となった。瞳の２倍の円の直径は６２ピクセルに相当する。そこで、ｆ方向に６２ピクセル分、ｇ方向に６２ピクセル分の正方領域内だけで
2019200080.tif 000044
【００９４】
を求め、その逆フーリエ変換として光学像の位相を算出した。参考までに、数２で光学像の位相を算出した。
【００９５】
図７（ａ）はシミュレーションで算出した正解位相である。図７（ｂ）は従来方法、すなわち、数２を用いて１０回ループを回して算出した位相で、図７（ｃ）は本発明の方法を用いて算出した位相である。明らかに従来方式では位相が正しく算出されていないことがわかる。
【００９６】
図６の光学像は、吸収がある物体の光学像である。そのため、強度が非常に弱い領域があり、数２を用いると解の精度が低下する。特に、繰り返しループで解の精度が低下していく。ところが、数４では１／Ｉの演算が入っていないため、解が極端に低下するということはない。さらに、本発明によれば数４を解く際に必要となる行列の要素数を大幅に減らすことができる。
【００９７】
顕微鏡の被検体が位相物体で、前記位相物体に吸収があるときや、位相が略１８０度反転している領域があると、観測強度に強度がほぼ０の領域が含まれてしまうリスクがある。本発明によれば、そのような場合でも解の精度の劣化を最小限に抑えて光学像の位相を算出することができる。
【００９８】
なお、強度輸送方程式は理論的にはコヒーレント光に対して成立する式である。しかし、実際にはコヒーレンス度が高ければコヒーレント光でなくても高い精度で成立する式である。顕微鏡では部分コヒーレント照明が用いられる。部分コヒーレント照明のコヒーレンス度を定量化するためにしばしばコヒーレンスファクターという量が使われる。コヒーレンスファクターとは、照明光の開口数を結像光学系の物体側開口数で割った値である。コヒーレンスファクターが０．３以下ならば、強度輸送方程式を解くことで物体の位相が高精度で算出できることを確かめた。
【００９９】
［実施例３］
本発明の強度輸送方程式ソルバーを用いて伝播光の位相を計測する装置を図８を用いて説明する。ここで、図８は、計測装置１００の概略ブロック図である。
【０１００】
計測装置１００は、ディテクタ１１０とステージ１２０からなり、それらは、コンピュータ１により制御されている。
【０１０１】
伝播光１３０の光軸をｚ軸とする。伝播光１３０の強度はディテクタ１１０で計測される。ディテクタ１１０は、ｘｙｚステージ１２０で保持されている。
【０１０２】
伝播光１３０がディテクタ１１０に入射するようにステージ１２０を駆動する。伝播光１３０がディテクタ１１０に入射するようアライメントされれば、コンピュータ１の制御部２０から強度取得の命令を信号入出力部５０を介してディテクタ１１０に送る。ディテクタ１１０で得られた光強度信号はコンピュータ１の信号入出力部５０に転送され、強度情報４０ｃとして記憶部４０に記憶される。
【０１０３】
次に、コンピュータ１の制御部２０は、ステージ１２０をｚ方向に所定の量駆動する命令を送り、ステージ１２０を駆動する。ステージ駆動量は記憶部４０に送られ伝播距離情報４０ｂとして保存される。その後、コンピュータ１の制御部２０から強度取得の命令を信号入出力部５０を介してディテクタ１１０に送る。ディテクタ１１０で得られた光強度信号はコンピュータ１の信号入出力部５０に転送され、強度情報４０ｃとして記憶部４０に記憶される。
【０１０４】
波長情報４０ａは使用する光源の波長に応じて使用者が入力部６０を介して入力する。また、ピクセルサイズ情報４０ｄはディテクタ１１０のカタログスペックに応じて使用者が入力部６０を介して入力する。
【０１０５】
以上、強度輸送方程式を解くために必要な情報を得ることができたので、コンピュータ１は図２、もしくは、図３記載のフローチャートに従い強度輸送方程式を解き、伝播光の位相を算出する。より詳しくは、前述の実施例１、もしくは、実施例２を参照していただきたい。
【０１０６】
以上、本発明の好ましい実施例について説明したが、本発明はこれらの実施例に限定されないことは言うまでもなく、その要旨の範囲内で種々の変形及び変更が可能である。例えば、本発明では２種類の強度分布を用いて強度輸送方程式を解いたが、３種類以上の強度分布を用いて強度輸送方程式を解くことも可能である。光の波長も可視領域に限定されず、例えば波長がＸ線であってもかまわない。
本発明は、伝播光の位相を強度輸送方程式から算出する方法及び装置に適用できる。
【符号の説明】
【０１０７】
１  コンピュータ
１００  計測装置
【図１】
2019200080.tif 000045
【図２】
2019200080.tif 000046
【図３】
2019200080.tif 000047
【図４】
2019200080.tif 000048
【図５】
2019200080.tif 000049
【図６】
2019200080.tif 000050
【図７】
2019200080.tif 000051
【図８】
2019200080.tif 000052