(57)【要約】【課題】任意のフォーカス位置での画像を高精度に推定する。【解決手段】画像推定方法は、被写体と撮像光学系の合焦位置との撮像光軸方向での距離が互いに異なるＮ枚の入力画像データを準備するステップと、該入力画像データを用いて特定距離での推定画像データを算出する推定ステップとを有する。推定ステップでは、Ｍ枚の入力画像データを逆フーリエ変換してＭ枚の逆フーリエ変換データを生成し、該逆フーリエ変換データのそれぞれのうち撮像光学系の瞳半径の２倍より大きい領域の空間周波数データを０にして空間周波数選択データを生成し、該空間周波数選択データを用いて、光学像の空間周波数データを表現する式として、ｚの（Ｍ−１）次以下の多項式を生成し、特定距離に対応する光学像の空間周波数データである特定空間周波数データを多項式を用いて算出し、特定空間周波数データをフーリエ変換して推定画像データを算出する。【選択図】図９
2014002636.tif 000002
【特許請求の範囲】
【請求項１】
撮像光学系により形成された被写体の光学像を撮像することで取得された画像データであって、前記被写体と前記撮像光学系の合焦位置との撮像光軸方向での距離が互いに異なるＮ回の撮像により取得されたＮ枚の入力画像データを準備する準備ステップと、
前記Ｎ枚の入力画像データを用いて、前記Ｎ回の撮像における前記距離とは異なる特定距離での推定画像データを算出する推定ステップとを有し、
該推定ステップは、
前記Ｎ枚の入力画像データのうち、Ｎ以下であるＭ枚の入力画像データを前記撮像光軸方向に直交する２方向に逆フーリエ変換してＭ枚の逆フーリエ変換データを生成するステップと、
前記Ｍ枚の逆フーリエ変換データのそれぞれのうち、前記撮像光学系の瞳半径の２倍より大きい領域の空間周波数データを０にすることで、Ｍ枚の空間周波数選択データを生成するステップと、
前記Ｍ枚の空間周波数選択データを用いて、前記光学像の空間周波数データを表現する式として、ｚの（Ｍ−１）次以下の多項式を生成するステップと、
前記特定距離に対応する前記光学像の空間周波数データである特定空間周波数データを前記多項式を用いて算出するステップと、
前記特定空間周波数データをフーリエ変換して前記推定画像データを算出するステップとを有することを特徴とする画像推定方法。
【請求項２】
Ｎが３であることを特徴とする請求項１に記載の画像推定方法。
【請求項３】
前記多項式を生成するステップにおいて、
前記多項式は（Ｍ−１）次の多項式であり、
前記Ｎ枚の入力画像データのうち、前記距離の前記特定距離に対する差が小さい順にＭ枚の入力画像データから得られた前記Ｍ枚の空間周波数選択データを用いて、Ｍ元１次方程式を解くことで、前記光学像の空間周波数データを前記（Ｍ−１）次の多項式で表すことを特徴とする請求項１または２に記載の画像推定方法。
【請求項４】
前記多項式を生成するステップにおいて、
前記多項式は（Ｍ−１）より小さいＭ′次の多項式であり、
前記Ｎ枚の入力画像データのうち、前記距離の前記特定距離に対する差が小さい順にＭ枚の入力画像データから得られた前記Ｍ枚の空間周波数選択データを用いて、前記光学像の空間周波数を前記Ｍ′次の多項式にフィッティングすることを特徴とする請求項１または２に記載の画像推定方法。
【請求項５】
前記Ｎ回の撮像において、前記被写体に部分コヒーレントな照明光が照射されることを特徴とする請求項１から４のいずれか１項に記載の画像推定方法。
【請求項６】
前記多項式は、２次の多項式であることを特徴とする請求項１から５のいずれか１項に記載の画像推定方法。
【請求項７】
コンピュータに、
撮像光学系により形成された被写体の光学像を撮像することで取得された画像データであって、前記被写体と前記撮像光学系の合焦位置との撮像光軸方向での距離が互いに異なるＮ回の撮像により取得されたＮ枚の入力画像データを取り込ませて、
前記Ｎ枚の入力画像データを用いて、前記Ｎ回の撮像における前記距離とは異なる特定距離での推定画像データを算出する画像推定処理を実行させるコンピュータプログラムであって、
前記画像推定処理は、
前記Ｎ枚の入力画像データのうち、Ｎ以下であるＭ枚の入力画像データを前記撮像光軸方向に直交する２方向に逆フーリエ変換してＭ枚の逆フーリエ変換データを生成するステップと、
前記Ｍ枚の逆フーリエ変換データのそれぞれのうち、前記撮像光学系の瞳半径の２倍より大きい領域の空間周波数データを０にすることで、Ｍ枚の空間周波数選択データを生成するステップと、
前記Ｍ枚の空間周波数選択データを用いて、前記光学像の空間周波数データを表現する式として、ｚの（Ｍ−１）次以下の多項式を生成するステップと、
前記特定距離に対応する前記光学像の空間周波数データである特定空間周波数データを前記多項式を用いて算出するステップと、
前記特定空間周波数データをフーリエ変換して前記推定画像データを算出するステップとを有することを特徴とする画像推定プログラム。
【請求項８】
撮像光学系により形成された被写体の光学像を撮像することで取得された画像データであって、前記被写体と前記撮像光学系の合焦位置との撮像光軸方向での距離が互いに異なるＮ回の撮像により取得されたＮ枚の入力画像データが入力される入力部と、
前記Ｎ枚の入力画像データを用いて、前記Ｎ回の撮像における前記距離とは異なる特定距離での推定画像データを算出する画像推定処理を行う処理部と有し、
前記処理部は、
前記Ｎ枚の入力画像データのうち、Ｎ以下であるＭ枚の入力画像データを前記撮像光軸方向に直交する２方向に逆フーリエ変換してＭ枚の逆フーリエ変換データを生成し、
前記Ｍ枚の逆フーリエ変換データのそれぞれのうち、前記撮像光学系の瞳半径の２倍より大きい領域の空間周波数データを０にすることで、Ｍ枚の空間周波数選択データを生成し、
前記Ｍ枚の空間周波数選択データを用いて、前記光学像の空間周波数データを表現する式として、ｚの（Ｍ−１）次以下の多項式を生成し、
前記特定距離に対応する前記光学像の空間周波数データである特定空間周波数データを前記多項式を用いて算出し、
前記特定空間周波数データをフーリエ変換して前記推定画像データを算出することを特徴とする画像推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、任意のフォーカス位置での画像を推定する画像推定技術に関する。
【背景技術】
【０００２】
病理診断支援システムの１つであるバーチャルスライドシステムでは、生体試料（被写体）を撮像して取得したデジタル画像を診断支援に用いる。撮像では、生体試料は、光学素子（カバーグラス）で覆われたプレパラートとして取り扱われる。バーチャルスライドシステムは、顕微光学系（撮像光学系）、撮像素子および処理装置によって構成され、プレパラートの光学像をデジタル画像に変換して保存する。
【０００３】
このようなバーチャルスライドシステムを用いる場合、医師等の診断者が参照できるのは、撮像を行ったフォーカス位置での画像のみである。診断者が画像から試料の立体的構造を推定して病理判断をするためには、互いに異なるフォーカス位置での撮像により取得された複数枚の画像が必要である。
【０００４】
ただし、多数枚の画像を取得するための撮像には長時間を要し、保存される画像データ量も膨大になるため、取得画像の枚数は最小限に抑えることが望ましい。一方、取得画像の枚数が少なすぎると、診断者が求めるフォーカス位置での画像を提供できない可能性がある。そこで、取得画像の枚数を削減しつつ任意のフォーカス位置での画像の提供を可能とするために、取得画像から任意のフォーカス位置での画像を推定する方法が提案されている。
【０００５】
特許文献１には、複数のフォーカス位置で取得した画像に対し、光学系の焦点ぼけフィルタを利用することで画像を推定する方法が開示されている。この方法では、フォーカス位置の変化によって生じるぼけを焦点ぼけフィルタを用いて調整することで、所望のフォーカス位置での画像を推定する。
【０００６】
また、非特許文献１には、被写体の光学像をフォーカス位置ｚの関数として表現し、光学像をｚの多項式で展開することにより近似的に光学像を推定する方法が開示されている。この方法では、複数のフォーカス位置で取得した画像から、各ピクセルにおける展開係数を求めることで任意のフォーカス位置での画像を推定する。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００１−２２３８７４号公報
【非特許文献】
【０００８】
【非特許文献１】Kenji Yamazoe and Andrew R. Neureuther 共著, ”Modeling of through-focus aerial image with aberration and imaginary mask edge effects in optical lithography simulation” Applied Optics, Vol.50, No.20, pp.3570-3578,10 July 2011, USA
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、特許文献１にて開示された方法では、光学系の焦点ぼけフィルタを予め知っておく必要があり、そのための予備計測等の事前情報の取得に手間がかかる。
【００１０】
また、非特許文献１にて開示された方法では、多項式の展開係数が、実際に撮像したデジタル画像に含まれるノイズの影響を受けるため、画像の推定精度が低下する可能性がある。
【００１１】
本発明は、フォーカス位置を異ならせての多数回の撮像や手間のかかる事前情報の取得が不要であり、かつ撮像により取得された画像がノイズを含んでいても任意のフォーカス位置での画像を高精度に推定できる画像推定方法、プログラムおよび装置を提供する。
【課題を解決するための手段】
【００１２】
本発明の一側面としての画像推定方法は、撮像光学系により形成された被写体の光学像を撮像することで取得された画像データであって、被写体と撮像光学系の合焦位置との撮像光軸方向での距離が互いに異なるＮ回の撮像により取得されたＮ枚の入力画像データを準備するステップと、Ｎ枚の入力画像データを用いて、Ｎ回の撮像における上記距離とは異なる特定距離での推定画像データを算出する推定ステップとを有する。該推定ステップは、Ｎ枚の入力画像データのうち、Ｎ以下であるＭ枚の入力画像データを撮像光軸方向に直交する２方向に逆フーリエ変換してＭ枚の逆フーリエ変換データを生成するステップと、Ｍ枚の逆フーリエ変換データのそれぞれのうち、撮像光学系の瞳半径の２倍より大きい領域の空間周波数データを０にすることで、Ｍ枚の空間周波数選択データを生成するステップと、Ｍ枚の空間周波数選択データを用いて、光学像の空間周波数データを表現する式として、ｚの（Ｍ−１）次以下の多項式を生成するステップと、特定距離に対応する光学像の空間周波数データである特定空間周波数データを多項式を用いて算出するステップと、特定空間周波数データをフーリエ変換して推定画像データを算出するステップとを有することを特徴とする。
【００１３】
また、本発明の一側面としての画像処理プログラムは、コンピュータに、撮像光学系により形成された被写体の光学像を撮像することで取得された画像データであって、前記被写体と前記撮像光学系の合焦位置との撮像光軸方向での距離が互いに異なるＮ回の撮像により取得されたＮ枚の入力画像データを取り込ませて、Ｎ枚の入力画像データを用いて、Ｎ回の撮像における上記距離とは異なる特定距離での推定画像データを算出する画像推定処理を実行させるコンピュータプログラムである。該画像推定処理は、Ｎ枚の入力画像データのうち、Ｎ以下であるＭ枚の入力画像データを撮像光軸方向に直交する２方向に逆フーリエ変換してＭ枚の逆フーリエ変換データを生成するステップと、Ｍ枚の逆フーリエ変換データのそれぞれのうち、撮像光学系の瞳半径の２倍より大きい領域の空間周波数データを０にすることで、Ｍ枚の空間周波数選択データを生成するステップと、Ｍ枚の空間周波数選択データを用いて、光学像の空間周波数データを表現する式として、ｚの（Ｍ−１）次以下の多項式を生成するステップと、特定距離に対応する光学像の空間周波数データである特定空間周波数データを多項式を用いて算出するステップと、特定空間周波数データをフーリエ変換して推定画像データを算出するステップとを有することを特徴とする。
【００１４】
さらに、本発明の他の一側面としての画像推定装置は、撮像光学系により形成された被写体の光学像を撮像することで取得された画像データであって、被写体と撮像光学系の合焦位置との撮像光軸方向での距離が互いに異なるＮ回の撮像により取得されたＮ枚の入力画像データが入力される入力部と、Ｎ枚の入力画像データを用いて、Ｎ回の撮像における上記距離とは異なる特定距離での推定画像データを算出する画像推定処理を行う処理部と有する。処理部は、Ｎ枚の入力画像データのうち、Ｎ以下であるＭ枚の入力画像データを撮像光軸方向に直交する２方向に逆フーリエ変換してＭ枚の逆フーリエ変換データを生成し、Ｍ枚の逆フーリエ変換データのそれぞれのうち、撮像光学系の瞳半径の２倍より大きい領域の空間周波数データを０にすることで、Ｍ枚の空間周波数選択データを生成し、Ｍ枚の空間周波数選択データを用いて、光学像の空間周波数データを表現する式として、ｚの（Ｍ−１）次以下の多項式を生成し、特定距離に対応する光学像の空間周波数データである特定空間周波数データを前記多項式を用いて算出し、特定空間周波数データをフーリエ変換して推定画像データを算出することを特徴とする。
【発明の効果】
【００１５】
本発明によれば、撮像時間の短縮、画像データ量の削減および事前情報の不要化を図りつつ、取得画像がノイズを含んでいる場合でも任意のフォーカス位置での画像を高精度に推定することができる。
【図面の簡単な説明】
【００１６】
【図１】本発明の実施例であるバーチャルスライドシステムの構成を示す図。
【図２】実施例のバーチャルスライドシステムで行ったシミュレーションの条件を示す図であり、（ａ）は有効光源、（ｂ）は被写体、（ｃ）撮像光学系の収差を示す。
【図３】上記シミュレーションで用いた光学像を示す図。
【図４】実施例において、画像の空間周波数をｚの２次式で表したときの係数を示す図。
【図５】（ａ）は実施例において推定したベストフォーカス画像を示す図であり、（ｂ）はシミュレーションにより求めたベストフォーカス画像を示す図である。
【図６】実施例において推定した画像をピーク信号対雑音比で評価した結果を示す図。
【図７】（ａ）は実施例において推定したベストフォーカス画像を示す図であり、（ｂ）シミュレーションにより求めたベストフォーカス画像を示す図である。
【図８】実施例において推定した画像をピーク信号対雑音比で評価した結果を示す図。
【図９】実施例における画像推定処理を示すフローチャート。
【発明を実施するための形態】
【００１７】
以下、本発明の実施例について図面を参照しながら説明する。
【００１８】
まず、本発明の実施例である画像推定方法を適用するバーチャルスライドシステム（顕微鏡システム）に設けられる撮像光学系に関する座標系について説明する。撮像光学系の座標系には、第１の座標系と第２の座標系とがある。
【００１９】
第１の座標系は、物体面と像面の座標を表すための座標系であり、本実施例では（ｘ，ｙ）で表す。物体面上の大きさと像面上の大きさは撮像光学系の倍率に応じて異なるが、以下では説明を簡単にするため、物体面上の大きさに撮像光学系の倍率を乗じて、像面上の大きさと１：１で対応させる。このため、物体面の座標系と像面の座標系も１：１で対応する。
【００２０】
第２の座標系は、光源面と撮像光学系の瞳面の座標を表すための座標系であり、本実施例では（ｆ，ｇ）で表す。（ｆ，ｇ）は、撮像光学系の瞳の半径が１になるように規格化した座標系である。
【００２１】
撮像光学系の光軸方向（ｚ方向：以下、撮像光軸方向という）における合焦位置（ベストフォーカス位置）、すなわち撮像光学系に対して撮像素子と共役な位置を、原点ｚ＝０とする。この原点に対する撮像光軸方向での位置ｚは、撮像光学系の合焦位置と被写体との撮像光軸方向での距離であり、フォーカスずれ量（デフォーカス量）に相当する。被写体の光学像はｘ，ｙ，ｚの関数として表され、本実施例では該関数をＩ（ｘ，ｙ，ｚ）と記す。
【００２２】
図１には、上記バーチャルスライドシステム（顕微鏡システム）の概略構成を示している。バーチャルスライドシステムは、撮像ユニット１００と、制御ユニット２００と、情報処理ユニット４００とから構成される。
【００２３】
制御ユニット２００は、搬送部２０１と制御部２０２とにより構成される。搬送部２０１は、制御部２０２の命令に基づいて、被写体１０３を、撮像ユニット１００に含まれる可動ステージ１０２上に移動させる。
【００２４】
撮像ユニット１００は、可動ステージ１０２上に搭載された被写体１０３を撮像してその画像を取得する撮像装置であり、照明系１０１、可動ステージ１０２、撮像光学系１０４および撮像素子１０５により構成されている。可動ステージ１０２は、制御部２０２からの命令に従って撮像光軸方向に移動することが可能である。この可動ステージ１０２の移動により、上述した撮像光学系の合焦位置（原点）と被写体１０３との撮像光軸方向での距離であるデフォーカス量を変更して複数回（Ｎ回）の撮像を行い、複数枚（Ｎ枚）の画像（入力画像データ）を取得することができる。つまり、被写体１０３が搭載された可動ステージ１０２を移動させる本実施例では、デフォーカス量は、可動ステージ１０２の位置に対応する。
【００２５】
なお、可動ステージ１０２は、撮像光軸方向に直交する方向に移動する機能を備えていてもよい。また、本実施例では、固定された撮像光学系１０４に対して可動ステージ１０２とともに被写体１０３が撮像光学系１０４の撮像光軸方向に移動する場合について説明するが、固定された被写体に対して撮像光学系がその撮像光軸方向に移動するように構成してもよい。この構成によっても、撮像光学系の合焦位置と被写体との撮像光軸方向での距離（デフォーカス量）を変えることができる。
【００２６】
また、可動ステージ１０２上の被写体１０３には、照明系１０１から部分コヒーレント光としての照明光が照射される。そして、撮像光学系１０４によって、被写体１０３の拡大光学像が撮像素子１０５上に形成される。撮像素子１０５は、ＣＣＤセンサまたはＣＭＯＳセンサにより構成され、被写体１０３の拡大光学像を光電変換する。撮像素子１０５から出力される電気信号（アナログ信号）は、情報処理ユニット４００に伝送される。
【００２７】
画像推定装置としての情報処理ユニット４００は、コンピュータ４０１、画像処理部４０２、データ保存部４０３および表示部４０４により構成される。画像処理部４０２は、撮像素子１０５から送られてきたアナログ電気信号をデジタル信号に変換して入力画像データを生成（取得）する。入力画像データは、画素ごとの輝度を示す輝度信号ということもできる。さらに、画像処理部４０２は、入力画像データに対して、ノイズ除去や圧縮等の各種画像処理を行うとともに、後述する画像推定処理を行い、その結果である推定画像データをコンピュータ４０１に伝送する。画像処理部４０２は、入力部および処理部として機能する。
【００２８】
コンピュータ４０１は、画像処理部４０２から送られてきた推定画像データを、データ保存部４０３に転送し、これを保存させる。病理診断時には、コンピュータ４０１が保存部４０３から推定画像データを読み出し、表示部４０４に表示させる。
【００２９】
コンピュータ４０１、画像処理部４０２、データ保存部４０３、表示部４０４および制御部２０２は、１つのコンピュータに含まれていてもよい。また、入力画像データや推定画像データ等の画像データを、外部のサーバー（例えば、病院に設置されたサーバー）に保存してもよい。この場合、複数の診断者が保存された画像データにアクセスし、これを見ることができる。
【００３０】
以下、バーチャルスライドシステムが行う一連の処理について説明する。この処理は、主としてコンピュータ４０１および画像推定処理用コンピュータとしての画像処理部４０２が、コンピュータプログラムに従って実行する。
【００３１】
まず、システムの使用者（検査者等）は、可動ステージ１０２を撮像光軸方向に移動させてデフォーカス量を異ならせた撮像を行い画像を取得する枚数Ｎ（つまりは撮像の回数Ｎ）と可動ステージ１０２の位置ｚｊ（１≦ｊ≦Ｎ）をコンピュータ４０１に入力する。ここでのｚｊの範囲は、最小値ｚｍｉｎから最大値ｚｍａｘまでとする。
【００３２】
コンピュータ４０１は、制御部２０２に可動ステージ１０２の駆動命令を送り、可動ステージ１０２を位置ｚｊまで移動させ、被写体の光学像Ｉ（ｘ，ｙ，ｚｊ）を撮像素子１０５により撮像する。以下の説明では、可動ステージ１０２の位置ｚｊを、それによってフォーカス状態が変わることから、フォーカス位置という。以上の動作を指定された全てのフォーカス位置ｚｊで行うことで、Ｎ回の撮像によるＮ枚の入力画像データ（以下、単に入力画像という）が取得される。Ｎ枚の入力画像は、データ保存部４０３に保存される。
【００３３】
以下、Ｎ枚の入力画像を用いて、任意のフォーカス位置（特定距離）ｚｓにおける画像を推定する画像推定方法（画像推定処理）について説明する。画像は被写体の光学像であるので、画像を推定することは光学像を推定することと等価である。逆に言えば、光学像が分かれば画像が分かる。任意の位置ｚｓにおける画像は、推定画像データに相当し、以下、単に推定画像ともいう。
【００３４】
先に説明した非特許文献１によれば、あるフォーカス位置ｚにおける光学像Ｉ（ｘ，ｙ，ｚ）は、
【００３５】
【数１】
2014002636.tif 000003
【００３６】
と展開可能である。式（１）は、ｚの２次式で光学像を表現している。非特許文献１では、さらに高次の多項式に展開できることも開示されているが、説明を簡単にするため、ここでは光学像をｚの２次式で表す。
【００３７】
非特許文献１では、式（１）のように光学像を展開するために、フォーカス位置を変えて取得した３つの光学像（入力画像）を必要としている。これは、式（１）には３つの未知数Ｑ（ｘ，ｙ）、Ｌ（ｘ，ｙ）およびＣ（ｘ，ｙ）があり、それらを求めるには方程式が３つ必要だからである。このため、異なるフォーカス位置で取得した３つの光学像Ｉ（ｘ，ｙ，ｚ１）、Ｉ（ｘ，ｙ，ｚ２），Ｉ（ｘ，ｙ，ｚ３）が必要である。
【００３８】
非特許文献１では、異なるフォーカス位置で取得した３つの光学像Ｉ（ｘ，ｙ，ｚｉ）（ただし１≦ｉ≦３）から、Ｑ（ｘ，ｙ）、Ｌ（ｘ，ｙ）およびＣ（ｘ，ｙ）を直接求めていた。このため、これら３つの光学像Ｉ（ｘ，ｙ，ｚｉ）にノイズが含まれていると、Ｑ（ｘ，ｙ）、Ｌ（ｘ，ｙ）およびＣ（ｘ，ｙ）もノイズの影響を受けてしまい、結果として任意のフォーカス位置における光学像の推定精度も低下する。
【００３９】
そこで本実施例では、ノイズの影響を低減しつつ、任意のフォーカス位置における光学像を推定する方法を提供する。
【００４０】
光学像Ｉ（ｘ，ｙ，ｚｉ）は、撮像光学系１０４によって形成された光学像である。撮像光学系１０４には瞳があるため、光学像の空間周波数成分には上限がある。具体的には、光学像をｘｙ方向（撮像光軸方向に直交する２方向）にフーリエ変換した場合、瞳半径の２倍以上に相当する空間周波数は存在しない。光学像Ｉ（ｘ，ｙ，ｚｉ）をｘｙ方向に逆フーリエ変換した関数を、
【００４１】
【数２】
2014002636.tif 000004
【００４２】
とする。もし（ｆ２＋ｇ２）１／２＞２ならば、
【００４３】
【数３】
2014002636.tif 000005
【００４４】
にならなければならない。つまり、撮像光学系１０４の瞳半径の２倍より大きい領域の空間周波数を０にする必要がある。このため、（ｆ２＋ｇ２）１／２＞２にて、
【００４５】
【数４】
2014002636.tif 000006
【００４６】
であるときは、それはノイズに起因したものであり、本来の光学像には存在しない空間周波数成分であるので無視してよい。
【００４７】
式（１）の両辺をｘｙ方向に逆フーリエ変換すると、
【００４８】
【数５】
2014002636.tif 000007
【００４９】
となる。上述したように（ｆ２＋ｇ２）１／２＞２ならば、
【００５０】
【数６】
2014002636.tif 000008
【００５１】
でなければならないので、（ｆ２＋ｇ２）１／２＞２ならば、
【００５２】
【数７】
2014002636.tif 000009
【００５３】
も０でなければならない。よって、（ｆ２＋ｇ２）１／２≦２の範囲内で、
【００５４】
【数８】
2014002636.tif 000010
【００５５】
を求め、その後、
【００５６】
【数９】
2014002636.tif 000011
【００５７】
を求める。こうして得られた
【００５８】
【数１０】
2014002636.tif 000012
【００５９】
をフーリエ変換すれば、ノイズの影響を低減した任意のフォーカス位置ｚでの光学像を推定することが可能となる。
【００６０】
３枚の画像から、
【００６１】
【数１１】
2014002636.tif 000013
【００６２】
を求めるには、（ｆ２＋ｇ２）１／２≦２の範囲内の各（ｆ，ｇ）で３元１次方程式を解けばよい。３枚より多い入力画像を用いる場合は、フィッティングを用いれば、
【００６３】
【数１２】
2014002636.tif 000014
【００６４】
を求めることができる。
【００６５】
具体的な処理を以下説明する。データ保存部４０３には、フォーカス位置を異ならせて取得したＮ枚の画像データＩ（ｘ，ｙ，ｚｊ）（１≦ｊ≦Ｎ）が保存されているとする。Ｎは３以上である（Ｎ≧３）。
使用者が推定画像を得たい任意のフォーカス位置ｚｓをコンピュータ４０１に入力すると、コンピュータ４０１は、３つ以上のｚｊの中からｚｓに最も近いｚｋを探す。そして、Ｉ（ｘ，ｙ，ｚｋ−１）、Ｉ（ｘ，ｙ，ｚｋ）およびＩ（ｘ，ｙ，ｚｋ＋１）の３枚（Ｍ枚：ＭはＮ以下）の入力画像データを画像処理部４０２にロードする。
画像処理部４０２は、入力画像データＩ（ｘ，ｙ，ｚｋ−１）、Ｉ（ｘ，ｙ，ｚｋ）およびＩ（ｘ，ｙ，ｚｋ＋１）の逆フーリエ変換データである
【００６６】
【数１３】
2014002636.tif 000015
【００６７】
を計算する。次に、画像処理部４０２は、（ｆ２＋ｇ２）１／２≦２の範囲内で、空間周波数選択データである
【００６８】
【数１４】
2014002636.tif 000016
【００６９】
を求める。その後、画像処理部４０２は、式（２）にｚｓの値を代入して、フォーカス位置ｚｓでの光学像の空間周波数データ（特定空間周波数データ）である
【００７０】
【数１５】
2014002636.tif 000017
【００７１】
を求める。最後に、画像処理部４０２は、
【００７２】
【数１６】
2014002636.tif 000018
【００７３】
をフーリエ変換し、フォーカス位置ｚｓでの推定画像データを取得（算出）する。
【００７４】
以上の説明では、フォーカス位置を変えて取得した３枚の入力画像から光学像の空間周波数をｚの２次関数（２次の多項式）として表した。より高精度な画像推定では、フォーカス位置を変えて取得したＭ枚の入力画像から、光学像の空間周波数をｚの（Ｍ−１）次関数（（Ｍ−１）次の多項式）として表してもよい。また、フォーカス位置を変えて取得したＭ枚の入力画像から、光学像の空間周波数をｚの（Ｍ−１）次より小さいＭ′次の関数にフィッティングしてもよい。例えば、フォーカス位置を変えて取得した５枚の画像から、光学像の空間周波数をｚの２次関数としてフィッティングして、
【００７５】
【数１７】
2014002636.tif 000019
【００７６】
を求めてもよい。フィッティングを用いれば、同じ多項式の次数でも誤差の影響を減らすことができる。つまり、光学像の空間周波数は、ｚの（Ｍ−１）次以下の多項式で表せばよい。
【００７７】
以下の具体的なシミュレーション例から明らかになるが、本実施例の画像推定方法は、従来の方法では達成が困難な画像推定が容易になる。特に、以下の３つの特徴がある。
【００７８】
第１に、本実施例の方法は位相物体にも適用が可能である。デフォーカス量を変えたとき、位相物体はベストフォーカス位置に対して非対称な光学像を形成するが、このような場合でも本実施例の方法によって画像推定が可能である。
【００７９】
第２に、本実施例の方法は、撮像光学系が収差を有していても適用が可能である。撮像光学系に収差があれば、やはりデフォーカス量を変えたときに光学像はベストフォーカス位置に対して非対称な変化をするが、このような場合でも本実施例の方法によって画像推定が可能である。
【００８０】
第３に、本実施例の方法は、ｚに関して不等間隔で取得した複数のフォーカス位置での入力画像からも任意のフォーカス位置での推定画像を算出することが可能である。  従来から、データのアップサンプリングにはフーリエ変換が用いられた。フォーカス位置を変えて取得した複数枚の入力画像を積層して３次元データとみなせば、フーリエ変換を用いてアップサンプリングが可能である。例えば、Δｚの間隔でフォーカス位置を変えて取得した入力画像から、フーリエ変換を用いてΔｚより小さい間隔でフォーカス位置が異なる光学像を算出可能である。
【００８１】
しかしながら、フーリエ変換を用いるには周期境界条件を満たす必要がある。物体が位相物体であったり撮像光学系が収差を有していたりすると、周期境界条件を満たさないため、フーリエ変換によるアップサンプリングの精度は著しく低下する。さらに、不等間隔で取得したデータに対しては、高速フーリエ変換（Fast Fourier Transform; FFT）が使えないため、処理スピードが著しく低下する。上記の観点から、本実施例の画像推定方法は、従来のフーリエ変換を用いたアップサンプリング方法よりも優位性があることが分かる。
【００８２】
非特許文献１を広く解釈すれば、ｚに関して不等間隔で取得されたデータをもとに任意のフォーカス位置における光学像を推定することは可能である。しかし、非特許文献１では計算の簡素化を図るため、ｚに関して等間隔でデータを取得することしか言及していない。
【００８３】
焦点ぼけフィルタに対しても優位性がある。焦点ぼけフィルタはコヒーレント結像やインコヒーレント結像時に定義されるのが一般的であり、顕微鏡のように部分コヒーレント結像で光学像が得られる撮像系で焦点ぼけフィルタを定義するのは非常に困難である。しかし、本発明は部分コヒーレント結像系を介して得られた光学像に関しても適用可能である。
【００８４】
以下、具体的なシミュレーション例により、上述した実施例の効果について説明する。
［シミュレーション例１］
シミュレーション例１では、撮像条件として、撮像光学系１０４の開口数（ＮＡ）を０．７とし、照明系１０１から被写体１０３に照射される照明光の波長λを５５０ｎｍとする。照明系１０１の開口数を撮像光学系１０４の開口数で割った値をコヒーレンスファクターσと定義する。
【００８５】
さらに、被写体１０３がないときに撮像光学系１０４の瞳面上に形成される照明の光強度分布を、有効光源と称する。本シミュレーション例において、部分コヒーレントな光源が形成する有効光源を、図２（ａ）に示すように設定した。図２（ａ）において、半径１の白い円は、σが１であることを示す。白い輪帯領域内の各ピクセルは、空間的にコヒーレントな点光源を表し、互いに異なるピクセルはインコヒーレントであるとする。ここでは、被写体１０３として、図２（ｂ）のような複数の矩形物体（矩形開口）を考える。背景は透光性を有し、すべての矩形物体の光強度の透過率は８０％である。各矩形物体を透過する光は背景に対して位相差を与え、さらにそれぞれの矩形物体は背景に対して互いに異なる位相差を与える。
【００８６】
左上の矩形物体は、背面に対して１７０度の位相差を与える。ある矩形物体の右隣の矩形物体は、その左隣の矩形物体よりも１０度少ない位相差を与える。ある矩形物体の下の矩形物体は、その上の矩形物体よりも６０度少ない位相差を与える。このルールに従えば、右下の矩形物体は−１８０度の位相差を与える。
【００８７】
撮像光学系１０４には、図２（ｃ）に示したような収差があるとする。図２（ｃ）において、単位はラジアンである。ここで、白い円は撮像光学系１０４の瞳を表す。
【００８８】
本シミュレーション例では、デフォーカス量（つまりはフォーカス位置）の単位ＲＵを０．５λ／ＮＡ２とする。つまり、５６１ｎｍを１ＲＵと記す。
【００８９】
本シミュレーション例において、−０．９ＲＵ、０．３ＲＵおよび０．９ＲＵでの画像（入力画像）を計算した結果、それぞれ図３（ａ）〜（ｃ）に示す画像が得られた。（ａ）は−０．９ＲＵのフォーカス位置の画像、（ｂ）は０．３ＲＵのフォーカス位置の画像、（ｃ）は０．９ＲＵのフォーカス位置の画像である。なお、図３（ａ）〜（ｃ）に示す画像はいずれも、５１２×５１２ピクセルの画像である。
【００９０】
次に、図３に示した画像を逆フーリエ変換し、さらに式（２）における
【００９１】
【数１８】
2014002636.tif 000020
【００９２】
を求めると、図４（ａ）〜（ｆ）に示す結果が得られた。（ａ）はｚ２項の係数である
【００９３】
【数１９】
2014002636.tif 000021
【００９４】
の実数部を示し、（ｂ）は該ｚ２項の係数の虚数部を示す。（ｃ）はｚの１乗項の係数である
【数２０】
2014002636.tif 000022
【００９５】
の実数部を示し、（ｄ）は該ｚの１乗項の係数の虚数部を示す。（ｅ）は定数項である
【００９６】
【数２１】
2014002636.tif 000023
【００９７】
の実数部を示し、（ｆ）は該定数項の虚数部を示す。また、図４（ａ）〜（ｆ）において、白い円は、
（ｆ２＋ｇ２）１／２＝２
を表す。すなわち、図４において白い円の外側はノイズとみなしてよく、０とした。
【００９８】
瞳の外側にある高い空間周波数は光学像形成に寄与しないため、瞳面で考慮すべきピクセル数は像面のピクセル数より少ない。例えば、図３（ａ）〜（ｃ）の画像と図４（ａ）〜（ｆ）の画像（逆フーリエ変換画像データ）とのピクセル数を比較する。図３の画像は５１２×５１２ピクセルであるのに対し、図４の画像は６４×６４ピクセルである。よって、非特許文献１の方法を用いれば、３元１次方程式を５１２×５１２回解かねばならない。これに対して、本実施例では、３元１次方程式（Ｍ元１次方程式）を６４×６４回解くだけでよい（実際には白い円の中は３１２５ピクセルしかないので、３１２５回方程式を解くだけでよい）。
このように、本実施例ではノイズの影響を低減するのみならず、解くべき方程式の数を大幅に減らすこともできる。
【００９９】
次に、一例として、ベストフォーカス位置（ｚ＝０）の画像を推定する。そのためには、
【０１００】
【数２２】
2014002636.tif 000024
【０１０１】
を求めればよいが、これは、
【０１０２】
【数２３】
2014002636.tif 000025
【０１０３】
そのものである。このため、
【０１０４】
【数２４】
2014002636.tif 000026
【０１０５】
をフーリエ変換すれば、ベストフォーカス位置の光学像、つまりは得たい推定画像としてのベストフォーカス画像を求めることができる。
【０１０６】
このようにして推定したベストフォーカス画像の例を、図５（ａ）に示す。また、参考画像として、別途求めたベストフォーカス画像を図５（ｂ）に示す。図５（ａ）の画像と図５（ｂ）の画像との差分の絶対値のうち最大値は、０．０１３であった。すなわち、最大で約１％程度の誤差しかないことが分かる。これは人間の目が認識するには困難な量の誤差である。実際、図５（ａ）の画像と図５（ｂ）の画像との違いを判別するのは難しい。
【０１０７】
図５（ａ）の画像のより定量的な評価として、図５（ｂ）の画像に対するピーク信号対雑音比（ＰＳＮＲ：Peak Signal-to-Noise Ratio）を求めると、値は５５．５デシベルとなった。一般的にピーク信号対雑音比が３５デシベル以上であるならば人間の目ではほとんど同じに見えると言われている。このため、本実施例の方法で十分な推定精度を達成できていると言える。
【０１０８】
さらに、フォーカス位置が−１．０ＲＵから１．０ＲＵまでの範囲にある各画像を推定し、また各フォーカス位置で推定画像とは別に参考画像を算出した。このとき、図６に示すように、全てのフォーカス位置で推定画像の参考画像に対するピーク信号対雑音比（ＰＳＮＲ）が５０デシベルを超え、高精度な画像推定ができた。なお、表示の都合上、ピーク信号対雑音比が１００を超えた場合には、それを１００とした。
［シミュレーション例２］
シミュレーション例１と同じ撮像条件で、５枚の入力画像から光学像の周波数成分をｚの４次関数として表し、任意のフォーカス位置における光学像を推定する。シミュレーションにより、−２．４ＲＵ、−０．６ＲＵ、０．６ＲＵ、１．２ＲＵおよび２．７ＲＵのフォーカス位置での入力画像を計算した。これら５枚の入力画像をフーリエ変換し、５元１次方程式を（ｆ２＋ｇ２）１／２≦２の範囲内の各（ｆ，ｇ）で解くことで、
【０１０９】
【数２５】
2014002636.tif 000027
【０１１０】
をｚの４次関数として表した。次に、
【０１１１】
【数２６】
2014002636.tif 000028
【０１１２】
にｚ＝０を代入して、ベストフォーカス画像の推定結果を図７（ａ）に示す。また、参考画像として、別途求めたベストフォーカス画像を図７（ｂ）に示す。図７（ａ）の画像と図７（ｂ）の画像との差分の絶対値のうち最大値は、０．００９７であった。本シミュレーション例では、入力画像を取得するフォーカス位置間の間隔（デフォーカス量差）がシミュレーション例１に比べて広いが、多項式の項数を増やすことにより、推定精度が向上した。図７（ａ）の画像の図７（ｂ）の画像に対するピーク信号対雑音比は５９．２デシベルとなった。これにより、本実施例の画像推定方法は、十分な推定精度を達成できていると言える。
【０１１３】
さらに、フォーカス位置が−３．０ＲＵから３．０ＲＵまでの範囲にある各画像を推定し、また各フォーカス位置で推定画像とは別に参考画像を算出した。このときの各フォーカス位置で推定画像の参考画像に対するピーク信号対雑音比（ＰＳＮＲ）を図８に示す。取得した入力画像間のデフォーカス差が大きいため、一部で精度の低下も見られるが、−３．０ＲＵと−２．７ＲＵのフォーカス位置での推定画像以外の推定画像はピーク信号対雑音比が４０デシベルを超えており、高精度な画像推定ができた。なお、表示の都合上、ピーク信号対雑音比が１００を超えた場合は、それを１００とした。
【０１１４】
上記５枚の画像より、最小二乗法を用いて（ｆ２＋ｇ２）１／２≦２の領域で、
【０１１５】
【数２７】
2014002636.tif 000029
【０１１６】
をｚの２次関数として表すことも可能である。結果の図示は省略するが、ベストフォーカス位置での画像を推定することができた。
【０１１７】
以上説明した画像推定処理を画像処理部４０２が実行するためのコンピュータプログラムである画像推定プログラムを図９のフローチャートに示す。ここでは、予め互いに異なるフォーカス位置での撮像により取得したＮ枚（Ｎ≧３）の入力画像Ｉ（ｘ，ｙ，ｚｊ）（ただし、１≦ｊ≦Ｎ）がデータ保存部４０３に保存（準備）されているとする。また、使用者により、予め推定画像を得たい任意のフォーカス位置（以下、特定フォーカス位置という）ｚｓと、光学像の空間周波数を多項式で表すときのその多項式の次数Ｍと、フィッティングで多項式を表すかどうかの指示とがコンピュータに入力されたものとする。
【０１１８】
ステップＳ１００では、画像処理部４０２は、データ保存部４０３から画像推定に使用するＭ枚の入力画像のデータをロードする。すなわち、フォーカス位置ｚｊの中からフォーカス位置ｚｓに最も近い値を持つｚｋを探す。フィッティングせずに光学像の空間周波数を多項式で表すときは、ｚｋの近傍で取得したＭ枚の入力画像のデータをデータ保存部４０３からロードする。フィッティングを用いて光学像の空間周波数を多項式で表すときは、ｚｋの近傍で取得した入力画像のデータを、フィッティングせずに光学像の空間周波数を多項式で表す場合のＭより大きい枚数だけデータ保存部４０３からロードする。
【０１１９】
次にステップＳ１０１では、画像処理部４０２は、ステップＳ１００で得られたＭ枚の入力画像データをｘｙ方向に逆フーリエ変換する。
【０１２０】
次にステップＳ１０２では、画像処理部４０２は、光学像の空間周波数をｚの多項式で表す。具体的には、コンピュータは、ステップＳ１０１で得られた各画像のフーリエ変換から、（ｆ２＋ｇ２）１／２≦２の範囲内だけを抽出する。言い換えれば、撮像光学系１０４の瞳半径の２倍より大きい（（ｆ２＋ｇ２）１／２＞２の）領域の空間周波数データを０にして、空間周波数選択データを生成する。そして、その空間周波数選択データを用いて、光学像の空間周波数をｚの（Ｍ−１）次の多項式で表す。多項式算出の際にフィッティングを用いるかどうかは、上述したように使用者が事前にコンピュータ４０１に入力した指示に従う。
【０１２１】
続いてステップＳ１０３では、画像処理部４０２は、特定フォーカス位置ｚｓにおける光学像の空間周波数である特定空間周波数データを算出する。すなわち、コンピュータは、ステップＳ１０２でｚの（Ｍ−１）次の多項式として表した光学像の空間周波数のｚにｚｓを代入し、特定フォーカス位置ｚｓにおける光学像の空間周波数である
【０１２２】
【数２８】
2014002636.tif 000030
【０１２３】
を算出する。このとき、前述したように、（ｆ２＋ｇ２）１／２＞２の領域では光学像の空間周波数は０である。
【０１２４】
最後にステップＳ１０４では、画像処理部４０２は、ステップＳ１０３で求めた特定フォーカス位置ｚｓにおける光学像の空間周波数をフーリエ変換して、特定フォーカス位置ｚｓの光学像に対応する推定画像を算出する。
【０１２５】
このようにして生成された推定画像は、表示部４０４に表示されてもよいし、データ保存部４０３に保存されてもよい。
【０１２６】
以上説明したように、本実施例の画像推定方法によれば、フォーカス位置を異ならせての多数回の撮像や手間のかかる事前情報の取得が不要であり、かつ撮像により取得された画像がノイズを含んでいても任意のフォーカス位置での画像を高精度に推定できる。つまり、撮像時間の短縮、画像データ量の削減および事前情報の不要化を図りつつ、取得画像がノイズを含んでいる場合でも任意のフォーカス位置での画像を高精度に推定することができる。
【０１２７】
以上説明した各実施例は代表的な例にすぎず、本発明の実施に際しては、各実施例に対して種々の変形や変更が可能である。
【産業上の利用可能性】
【０１２８】
顕微鏡システム等において任意のフォーカス位置での画像を高精度に推定できる画像推定技術を提供できる。
【符号の説明】
【０１２９】
１０３  被写体
１０４  撮像光学系
１０５  撮像素子
４０２  画像処理部
【図１】
2014002636.tif 000031
【図２】
2014002636.tif 000032
【図３】
2014002636.tif 000033
【図４】
2014002636.tif 000034
【図５】
2014002636.tif 000035
【図６】
2014002636.tif 000036
【図７】
2014002636.tif 000037
【図８】
2014002636.tif 000038
【図９】
2014002636.tif 000039