JP2020057163A - 点像分布関数推定方法及び点像分布関数推定装置 - Google Patents

Abstract

【課題】試料の観察時に迅速に点像分布関数を用いたボケの除去を行うこと。【解決手段】点像分布関数推定方法は、光学系を用いて試料Ｓを撮像した結果得られた試料画像データと、光学系を用いて円形の開口部１９の像を試料Ｓと同時に撮像した結果得られた参照画像データとを取得する取得ステップＳ１０２と、参照画像データと開口部１９の形状に対応した理想画像データとを基に光学系の点像分布関数を推定する推定ステップＳ１０４，Ｓ１０５，Ｓ１０６と、点像分布関数を用いて試料画像データをデコンボリューションすることにより出力用試料画像データを生成する生成ステップＳ１０９と、点像分布関数と出力用試料画像データとを出力する出力ステップＳ１１０と、を備える。【選択図】図５

Description

本発明の一側面は、光学系における点像分布関数を推定する点像分布関数推定方法及び点像分布関数推定装置に関する。

従来から、顕微鏡等の光学系における分解能の指標となる点像分布関数（PSF：Point Spread Function）を推定することが行われている。例えば、下記特許文献１には、点像分布関数を推定する装置が開示され、この装置は、環状の形状を有する図形が行列状に配置された参照チャートを撮像装置の光学系を構成するレンズを通して撮像し、撮像画像と１個の図形に対応した参照画像とを用いて点像分布関数を推定する。

特開２０１６−１６３３１２号公報

上述したような従来の装置では、仮想的に設定された点像分布関数を用いて参照画像を畳み込み演算することにより得られる中間画像と撮像画像との差分が小さくなるように点像分布関数が演算されている。この装置では、参照チャートを撮像した後に繰り返し演算により点像分布関数が推定されているので、試料の観察の際に迅速な点像分布関数の推定を実現することは困難である。

そこで、本発明の一側面は、かかる課題に鑑みてなされたものであり、試料の観察時に迅速な点像分布関数を用いたボケの除去を可能にする点像分布関数推定方法及び点像分布関数推定装置を提供することを課題とする。

本発明の一形態は、光学系を用いて試料を撮像した結果得られた試料画像データと、光学系を用いて既知の回転対称な形状の像を試料と同時に撮像した結果得られた参照画像データとを取得する取得ステップと、参照画像データと、形状の情報とを基に光学系の点像分布関数を推定する推定ステップと、点像分布関数を用いて試料画像データをデコンボリューションすることにより出力用試料画像データを生成する生成ステップと、点像分布関数と出力用試料画像データとを出力する出力ステップと、を備える点像分布関数推定方法である。

あるいは、本発明の他の形態は、試料の像及び既知の回転対称な形状の像を同時に導く光学系と、試料の像を光学系を介して撮像して試料画像データを取得すると同時に、既知の回転対称な形状の像を光学系を介して撮像して参照画像データを取得する撮像装置と、参照画像データと、形状の情報とを基に光学系の点像分布関数を推定し、点像分布関数を用いて試料画像データをデコンボリューションすることにより出力用試料画像データを生成し、点像分布関数と出力用試料画像データとを出力する演算装置と、を備える点像分布関数推定装置である。

上記いずれかの形態によれば、光学系を用いて試料を撮像して得られた試料画像データと回転対称な形状の像を撮像して得られた参照画像データとが同時に得られ、参照画像データと形状の情報とを基に推定された点像分布関数を用いて、試料画像データの画像からボケが除去された出力用試料画像データが生成され、点像分布関数及び出力用試料画像データが出力される。これにより、試料の観察時に迅速にボケが除去された試料の画像を得ることができる。

上記一形態あるいは上記他の形態においては、回転対称な形状は円形である、ことが好適である。この場合、光学系の角度調整等の煩雑な調整を必要とせずに、正確な点像分布関数を得ることができる。

また、上記一形態においては、取得ステップで、試料を保持する保持部及び光を通過させる円形の開口部を有する試料ホルダと、試料ホルダに向けて光を照射する光源とを用いて、試料画像データと参照画像データとを取得する、ことが好適である。また、上記他の形態においては、試料を保持する保持部及び光を通過させる円形の開口部を有する試料ホルダと、試料ホルダに向けて光を照射する光源と、をさらに備え、撮像装置は、光源から光が照射された試料の像、及び、開口部を通過する光によって生ずる開口部の像を撮像する、ことが好適である。かかる構成を備えれば、光が照射された試料の撮像と同時に、円形の開口部を通過する光によって形成される像を撮像することによって、簡易に試料画像データ及び参照画像データを得ることができる。

さらに、上記一形態においては、推定ステップで、点像分布関数がゼロ以上であることを拘束条件として点像分布関数を推定する、ことも好適である。さらに、上記他の形態においては、演算装置は、点像分布関数がゼロ以上であることを拘束条件として点像分布関数を推定する、ことも好適である。こうすれば、ノイズにロバストな点像分布関数の推定を実現することができる。

またさらに、上記一形態においては、推定ステップで、参照画像データと形状の情報とをフーリエ変換した値と、フーリエ変換後の点像分布関数とに関する最小化問題を解くことにより、点像分布関数を推定する、ことも好適である。またさらに、上記他の形態においては、演算装置は、参照画像データと形状の情報とをフーリエ変換した値と、フーリエ変換後の点像分布関数とに関する最小化問題を解くことにより、点像分布関数を推定する、ことも好適である。かかる構成を採れば、高速な演算によって短時間で点像分布関数を推定することができる。

さらにまた、上記一形態においては、推定ステップで、フーリエ変換した値と、フーリエ変換後の点像分布関数と、点像分布関数のフーリエ空間における滑らかさを示す値とに関する最小化問題を解くことにより、点像分布関数を推定する、ことが好適である。さらにまた、上記他の形態においては、演算装置は、フーリエ変換した値と、フーリエ変換後の点像分布関数と、点像分布関数のフーリエ空間における滑らかさを示す値とに関する最小化問題を解くことにより、点像分布関数を推定する、ことが好適である。かかる構成を採れば、短時間で真値に近い点像分布関数を推定することができる。

また、上記一形態においては、推定ステップで、推定した点像分布関数を基に線像分布関数を計算し、線像分布関数を任意の関数でフィッティングし、フィッティングした線像分布関数を用いて点像分布関数を再構成する、ことが好適である。また、上記他の形態においては、演算装置は、推定した点像分布関数を基に線像分布関数を計算し、線像分布関数を任意の関数でフィッティングし、フィッティングした線像分布関数を用いて点像分布関数を再構成する、ことが好適である。こうすれば、実空間で滑らかな点像分布関数を得ることができ、推定精度をより向上できる。

さらに、上記一形態においては、推定ステップでは、線像分布関数を撮像面における複数の方向に沿って積分することによって複数の線像分布関数を計算し、複数の線像分布関数を任意の関数でフィッティングし、フィッティングした複数の線像分布関数を用いて点像分布関数を再構成する、ことも好適である。さらに、上記他の形態においては、演算装置は、線像分布関数を撮像面における複数の方向に沿って積分することによって複数の線像分布関数を計算し、複数の線像分布関数を任意の関数でフィッティングし、フィッティングした複数の線像分布関数を用いて点像分布関数を再構成する、ことも好適である。こうすれば、実空間で滑らかな点像分布関数を得ることができ、推定精度をより向上できる。

本発明の一側面によれば、試料の観察時に迅速に点像分布関数を用いたボケの除去が可能となる。

実施形態にかかる点像分布関数推定装置である観察システム１の概略構成図である。 図１の試料ホルダ７を結像用レンズ９側から見た平面図である。 図１の演算装置１３のハードウェア構成を示す図である。 図１の演算装置１３の点像分布関数の推定時に用いられるカーネルの構成を示す図である。 実施形態にかかる観察システム１による点像分布関数の推定から観察画像の出力までの手順を示すフローチャートである。 シミュレーションによって推定された点像分布関数の画像イメージを示す図である。 シミュレーションによって推定された点像分布関数の一次元プロファイルを示す図である。 シミュレーションによって推定された点像分布関数の画像イメージを示す図である。 シミュレーションによって推定された点像分布関数の一次元プロファイルを示す図である。 図９の一次元プロファイルを拡大して示す図である。 実験光学系の概略構成図である。 図１１に示す実験光学系を用いて推定された点像分布関数の画像イメージを示す図である。 図１１に示す実験光学系を用いて推定された点像分布関数の一次元プロファイルを示す図である。 図１１に示す実験光学系を用いて推定された点像分布関数の一次元プロファイルを示す図である。

以下、添付図面を参照して、本発明の実施形態について詳細に説明する。なお、説明において、同一要素又は同一機能を有する要素には、同一符号を用いることとし、重複する説明は省略する。

図１は、実施形態にかかる点像分布関数推定装置である観察システム１の概略構成図である。図１に示す観察システム１は、観察システム１に含まれる光学系に関する点像分布関数を推定すると同時に、推定した点像分布関数を用いて試料を撮像して得られた画像からボケを除去してその画像を出力する装置である。点像分布関数は、光学系の点光源に対する応答を表す関数であり、一般に、この点像分布関数が推定できれば光学系を介して得られた画像のボケを除去することができる。この観察システム１は、光源３、集光用レンズ（集光用光学系）５、試料ホルダ７、結像用レンズ（結像用光学系）９、カメラ（撮像装置）１１、及び演算装置１３を含んでいる。これらの光源３、集光用レンズ５、及び結像用レンズ９は、試料を観察するための顕微鏡（光学系）を構成する。

光源３は、試料ホルダ７に保持された試料に光を照射するＬＥＤ（Light EmittingDiode）、ランプ光源等の光源装置である。この光源３としては、発光波長が特定の波長のものには限定されず、可視光領域、Ｘ線領域等の任意の発光波長のものが選択されうる。集光用レンズ５は、光源３から照射された光を試料ホルダ７に向けて集光する光学系である。結像用レンズ９は、試料ホルダ７に保持された試料を透過した光をカメラ１１に結像させるための光学系である。

試料ホルダ７は、観察対象の試料を保持するための部材であり、集光用レンズ５と結像用レンズ９との間に配置されている。図２は、試料ホルダ７を結像用レンズ９側から見た平面図である。このように、試料ホルダ７は、平板状部材１５によって構成され、その平板状部材１５の中央に設けられた、試料Ｓが保持される光透過部（保持部）１７と、平板状部材１５において光透過部１７の近傍で光透過部１７を囲むように点対称な位置に設けられた４つの円形の開口部１９とを有する。この開口部１９は、既知の回転対称な形状を有していれば、円形の形状には限定されずリング状等であってもよいが、より正確な点像分布関数を得るという観点からは円形であることが好ましい。また、開口部１９の位置および数は任意に変更されてもよく、開口部１９の代わりに、例えば、光透過部１７内に円形状の不透明の（あるいは半透明の）物体が配置されていてもよい。光透過部１７は、光源３からの光を透過する材料によって構成され、例えば、光源３としてＸ線源が用いられる場合にはシリコン窒化膜（ＳｉＮ膜）によって構成される。このような構成の試料ホルダ７は、集光用レンズ５と結像用レンズ９との間において、集光用レンズ５によって集光された光が光透過部１７及び４つの開口部１９の全体に向けて照射されるように、平板状部材１５が光の照射方向にほぼ垂直になるように配置される。

図１に戻って、カメラ１１は、試料ホルダ７によって保持されて光源３からの光が透過した試料Ｓの像と、試料ホルダ７の開口部１９を光が通過することによって生ずる開口部１９の像とを、結像用レンズ９を介して同時に撮像し、それらの像を表す２次元画像データを出力する撮像装置である。このカメラ１１は、例えば、ＣＣＤ（Charge Coupled Device）イメージセンサ、ＣＭＯＳ（ComplementaryMetal Oxide Semiconductor）イメージセンサ等である。演算装置１３は、カメラ１１に対して有線接続あるいは無線接続によってデータ通信可能に接続され、カメラ１１から出力された２次元画像データを用いて、観察システム１の光学系に関する点像分布関数を推定するとともに、推定した点像分布関数を用いてボケが除去された試料Ｓの画像データを出力する装置である。演算装置１３の機能の詳細については後述する。

図３は、演算装置１３のハードウェア構成を示している。図３に示すように、演算装置１３は、物理的には、プロセッサであるＣＰＵ（Central Processing Unit）３０１、記録媒体であるＲＡＭ（Random Access Memory）３０２又はＲＯＭ（Read Only Memory）３０３、通信モジュール３０４、及び入出力モジュール３０６等を含んだコンピュータ等であり、各々は電気的に接続されている。演算装置１３の各機能は、ＣＰＵ３０１及びＲＡＭ３０２等のハードウェア上にプログラム等を読み込ませることにより、ＣＰＵ３０１の制御のもとで、通信モジュール３０４、及び入出力モジュール３０６等を動作させるとともに、ＲＡＭ３０２におけるデータの読み出し及び書き込みを行うことで実現される。なお、演算装置１３は、入出力デバイスとして、ディスプレイ、キーボード、マウス、タッチパネルディスプレイ等を含んでいてもよいし、ハードディスクドライブ、半導体メモリ等のデータ記録装置を含んでいてもよい。また、演算装置１３は、単一の装置であってもよいし、相互にデータ通信可能に構成された複数の装置であってもよい。

以下、演算装置１３による点像分布関数の推定機能の詳細を説明する。演算装置１３は、フーリエ空間において最小化問題を解くことにより点像分布関数を推定する第１の計算手法による推定機能と、実空間において最小化問題を解くことにより点像分布関数を推定する第２の計算手法による推定機能とを有している。演算装置１３で実行される計算手法は、カメラ１１によって得られる試料の像におけるＳ／Ｎの高低に応じて、第１の計算手法と第２の計算手法のうちから選択される。この選択は、演算装置１３に対するユーザからの入力に応じて行われてもよいし、光源３の発光強度等を基にＳ／Ｎの高低を判断することにより演算装置１３により自動的に行われてもよい。演算装置１３は、点像分布関数の推定処理を行う際には、カメラ１１から出力された２次元画像データの中から試料Ｓの像の部分を切り出して試料画像データとして予め取得し、２次元画像データの中から１つの開口部１９の像の部分を切り出して参照画像データとして予め取得する。加えて、演算装置１３は、開口部１９の既知の円形の形状を基に予め設定された開口部１９の理想像に対応した理想画像データ（形状の情報）を取得する。例えば、開口部１９の場合は、予め計測した開口部１９の半径に対応する半径を有する円の内部が白であり、その円の外部が黒である画像データを、理想画像データとして取得する。開口部１９の代わりに円形の不透明物体を用いる場合は、その物体の半径に対応する半径を有する円の内部が黒であり、その円の外部が白である画像データを、理想画像データとして取得する。

まず、演算装置１３の第１の計算手法による推定機能について説明する。

この機能では、参照画像データのフーリエ変換後の値Ｇと、理想画像データのフーリエ変換後の値Ｆとを用いて、点像分布関数のフーリエ変換後の値Ｐを、下記式（１）及び下記式（２）に示すフーリエ空間での最小化問題を解くことにより推定される。




上記式（１）、（２）中、添え字「＊」は推定解であることを示し、添え字「Ｒ」、「Ｉ」はそれぞれ実成分及び虚成分を示し、添え字「ｉ」、「ｊ」はフーリエ空間での周波数位置を示し、「ε」は任意の実定数、「λ」は予め適切な値に設定された実定数である。なお、上記式（１）、（２）の右辺第２項は、点像分布関数の値Ｐのフーリエ空間における滑らかさを示す値であり、具体的には、値ＰのグラディエントのＬ２ノルムの二乗によって計算され、上記式（１）、（２）を解くことは、値Ｐ，Ｇ，Ｆと滑らかさを示す値とに関する最小化問題を解くことを意味する。

そして、演算装置１３は、上記のようにして推定した点像分布関数のフーリエ変換値の実成分Ｐ_Ｒ及び虚成分Ｐ_Ｉから点像分布関数のフーリエ変換値Ｐを下記式（３）；


を用いて算出する。上記式（３）中、「ｉ」は虚数単位を示す。さらに、演算装置１３は、点像分布関数のフーリエ変換値Ｐを逆フーリエ変換することにより、下記式（４）；


を用いて、点像分布関数の推定値ｐ^＊を算出する。

上記第１の計算手法による推定の原理について説明すると、一般に、光学系を介して実際に観測される像ｇは、その光学系の点像分布関数ｐによって理想的な像ｆがぼけた像となってしまい、下記式（５）に示すように、像ｆと点像分布関数ｐとのコンボリューション（畳み込み演算）によって表される。


また、畳み込み定理により、フーリエ変換後の値Ｇ，Ｐ，Ｆは、下記式（６）；


の関係を有することが知られている。この関係から、従来の点像分布関数の推定手法では、下記式（７）；


を用いて、点像分布関数を推定していた（逆フィルタ法）。

上記逆フィルタ法では、ノイズに対して脆弱であることが知られており、逆フィルタ法を改善したウィーナフィルタ法が、ノイズに対してロバストな推定手法として用いられている。このウィーナフィルタ法では、下記式（８）；


を用いて点像分布関数が推定される。上記式（８）中、「ε」は任意の実定数である。ウィーナフィルタ法では「ε」を大きくすることによりノイズ耐性を改善することができるが、その一方で「ε」が大きすぎると点像分布関数の推定精度が下がってしまう。そこで、本実施形態での第１の計算手法では、フーリエ空間において点像分布関数の実成分及び虚成分がそれぞれフーリエ空間において滑らかであることも拘束条件として点像分布関数を推定する。すなわち、式（１）、（２）において、右辺の第１項はウィーナフィルタ法での推定のための式（８）に対応する項であり、右辺の第２項は滑らかさを拘束条件として要求する項であり、実定数λが大きいほど滑らかな解が得られる。本実施形態では、実定数εを小さく保ちながら、実定数λを適切な値に調整することにより、点像分布関数の解の真値からの乖離を小さく保ちつつ推定処理のノイズ耐性を改善することができる。

特に、演算装置１３は、式（１）、（２）をフーリエ変換を用いて高速に解くことができる。すなわち、式（１）、（２）を解くことは、下記式（９）の形式の最小化問題を解くことと同義である。


上記式（９）を変数ｘ_ｉ，ｊについて微分してゼロに等しいとすれば、下記式（１０）；


が導かれる。この式（１０）は、変数ｘと図４に示すようなカーネルとの畳み込み演算をした値が変数ｙに一致することを意味している。よって、演算装置１３は、点像分布関数のフーリエ変換値Ｐを、式（１０）の形式に変形してから上述した逆フィルタ法を用いて解くことにより高速に求めることができる。

次に、演算装置１３の第２の計算手法による推定機能について説明する。

この機能では、点像分布関数の実空間の値ｐが、参照画像データの実空間の値ｇと、理想画像データの実空間の値ｆとを用いて、下記式（１１）に示す実空間での最小化問題を解くことにより推定される。


上記式（１１）に示すように、値ｐがゼロ以上であることを拘束条件として、値ｆと値ｐとの畳み込み演算の値が値ｇに近くなるように、点像分布関数の値ｐが推定される。つまり、点像分布関数は非負でなければならないので、観測された像を再現し、かつ、非負であるような点像分布関数が推定される。この式の解は、ラグランジュ（Lagrange）の未定定数法を用いて拘束条件無しの最小化問題に帰着することで、最急降下法などの勾配法を用いて求めることができる。この第２の計算手法は、第１の計算手法に比較してノイズにロバストである一方、第２の計算手法を用いた際の演算時間は逐次近似法を用いて繰り返し演算により解を推定する必要があるために第１の計算手法に比較して長くなる傾向にある。なお、観測された画像の画素値のデータの統計性を考慮し、例えば、データがポアソン分布に従うときは、演算装置１３は、第２の計算手法として、下記式（１２）；


を用いて値ｐの解を求めることができ、この場合はより正確な解を求めることができる。

上記式（１２）を用いた第２の計算方法は、以下の考え方に基づいている。まず、単一画素について考える。画素値がポアソン分布に従うとする。ポアソン分布は平均値λをパラメータとして含む。つまり、計測を多数回繰り返し、注目している単一画素の値を平均するとλに漸近する。画素値としてｇが得られる確率は、下記式（１３）；


によって計算される値となる。次に複数の画素を考える。画素の位置を表す添え字をｉとすると、ｉ番目の画素で画素値ｇ_ｉが得られる確率は、下記式（１４）；


によって計算される値となる。ただしλ_ｉは画素ｉにおける画素値の平均値である。画素ｉ＝０〜Ｎ−１において画素値ｇ_０〜ｇ_Ｎ−１が得られる確率は、下記式（１５）；


によって計算される。これは、総乗記号Πを用いて、下記式（１６）；


のように表される。

今、ぼけのない像をｆ、点像分布関数をｐとすれば、観測される画像はｆとｐとのコンボリューションとなる。これは統計的揺らぎを考慮しない場合である。ポアソン統計に従うとすれば、画素ｉの値は平均がコンボリューションのポアソン分布に従う。つまり、上記式（１６）のλとしてｆとｐとのコンボリューションを代入すれば上記式（１２）を得ることができる。

さらに、演算装置１３は、第１の計算手法による推定機能及び第２の計算手法による推定機能の他、両推定機能によって求められた点像分布関数の値ｐがノイズが多く不十分な場合には次のようなフィルタ処理を行う機能（フィルタ機能）も有する。このフィルタ機能では、顕微鏡の典型的な線像分布関数がガウス関数の重ね合わせであるという性質が利用される。すなわち、演算装置１３は、推定した点像分布関数の値ｐを撮像面上のある角度θの方向に沿って順投影（線積分）することにより線像分布関数の値ｈを計算する。そして、演算装置１３は、得られた線像分布関数の値ｈをガウス関数によってフィッティングし、フィッティングにより線像分布関数の値ｈを再取得する。このとき、演算装置１３は、複数の線像分布関数の値ｈを、予め適切に設定された角度ステップΔθで０度から１８０度の範囲で変更された複数の角度θで順投影することにより計算し、計算した複数の線像分布関数の値ｈを対象にフィッティングを実行する。さらに、演算装置１３は、フィッティングにより再取得した複数の線像分布関数ｈを再構成することによって最終的に推定された点像分布関数の値ｐを取得する。この点像分布関数への再構成は、複数の角度θの線像分布関数を投影像としてフィルタドバックプロジェクション（Filtered Back Projection）法等のＣＴ（Computed Tomography）画像再構成のアルゴリズムを用いて実行される。これにより、点像分布関数の先鋭度を損ねることなく推定結果を改善することができる。なお、フィッティングに用いる関数はガウスライク関数には限定されず、事前に系の線像分布関数形が予測できればその関数形の関数を使ってフィッティングすることが好適である。

以下、上述した観察システム１を用いた点像分布関数の推定方法の手順、及び、試料Ｓの観察画像の取得方法の手順について説明する。図５は、観察システム１による点像分布関数の推定から観察画像の出力までの手順を示すフローチャートである。

まず、ユーザの所定の操作により試料Ｓの観察が開始されると、観察システム１の光学系を介して試料Ｓの像と試料ホルダ７の開口部１９の像とが同時に撮像される（ステップＳ１０１）。これにより、カメラ１１によって生成された２次元画像データを基に、演算装置１３において試料画像データと参照画像データとが同時に取得される（ステップＳ１０２）。

次に、演算装置１３により、試料画像データのＳ／Ｎの高低が判断される（ステップＳ１０３）。この判断は、例えば、光源３の発光強度等に基づいて判断される。判断の結果、Ｓ／Ｎが比較的高い場合には（ステップＳ１０３；Ｙｅｓ）、第１の計算手法が選択され、演算装置１３により、第１の計算手法を用いてフーリエ空間での点像分布関数の値Ｐが計算される（ステップＳ１０４）。その後、演算装置１３により、値Ｐが逆フーリエ変換されることにより、実空間の点像分布関数の推定値ｐが計算される（ステップＳ１０５）。一方で、Ｓ／Ｎが比較的低いと判断された場合には（ステップＳ１０３；Ｎｏ）、第２の計算手法が選択され、演算装置１３により、第２の計算手法を用いて実空間での点像分布関数の推定値ｐが計算される（ステップＳ１０６）。

さらに、演算装置１３により、計算した点像分布関数の推定値ｐが良好であるかが判断される（ステップＳ１０７）。この判断は、推定値ｐにおけるノイズの高低に基づいて行われる。より具体的には、コントラスト、あるいは点像分布周辺部の分散等を指標として判断される。判断の結果、推定値ｐが良好でない場合には（ステップＳ１０７；Ｎｏ）、演算装置１３によって、フィルタ処理が施されて点像分布関数の値ｐが改善される（ステップＳ１０８）。一方、推定値ｐが良好であると判断された場合には（ステップＳ１０７；Ｙｅｓ）、演算装置１３によって、ステップＳ１０５またはステップＳ１０６で推定された値ｐがそのまま利用される。

次に、演算装置１３により、推定した点像分布関数の値ｐを用いて、試料画像データがデコンボリューション（畳み込み演算の逆演算）されることにより、出力用試料画像データが生成される（ステップＳ１０９）。そして、演算装置１３により、推定した点像分布関数の値ｐと、出力用試料画像データとが、ディスプレイ等の出力デバイスに画像イメージとして同時に出力される（ステップＳ１１０）。これにより、ユーザに点像分布関数の妥当性をリアルタイムで確認させると同時に、ボケの除去された試料Ｓの観察画像を出力することができる。

以上説明した観察システム１及びそれを用いた点像分布関数推定方法によれば、光学系を用いて試料Ｓを撮像して得られた試料画像データと開口部１９の像を撮像して得られた参照画像データとが同時に得られ、参照画像データと開口部１９の形状に対応した理想画像データとを基に推定された点像分布関数を用いて、試料画像データの画像からボケが除去された出力用試料画像データが生成され、点像分布関数及び出力用試料画像データが画像イメージとして出力される。これにより、試料Ｓの観察時に迅速にボケが除去された試料Ｓの画像を得ることができる。

また、本実施形態においては、試料Ｓを保持する光透過部１７及び円形の開口部１９を有する試料ホルダ７と、試料ホルダ７に向けて光を照射する光源３とを用いて、試料画像データと参照画像データとが取得されている。このような構成により、光が照射された試料Ｓの撮像と同時に、円形の開口部１９を通過する光によって形成される像を撮像することによって、簡易に試料画像データ及び参照画像データを得ることができる。

さらに、本実施形態における第２の計算手法を用いた点像分布関数の推定においては、点像分布関数がゼロ以上であることを拘束条件として推定されている。こうすれば、ノイズにロバストな点像分布関数の推定を実現することができる。

また、本実施形態における第１の計算手法を用いた点像分布関数の推定おいては、フーリエ変換した値Ｐ，Ｇ，Ｆと点像分布関数のフーリエ空間における滑らかさを示す値とに関する最小化問題を解くことにより、点像分布関数が推定されている。かかる手法によれば、高速な演算によって短時間で真値に近い点像分布関数を推定することができる。

また、本実施形態では、推定した点像分布関数の値が不十分な場合には、その値にさらにフィルタ処理が施されている。こうすれば、実空間で滑らかな点像分布関数を得ることができ、点像分布関数の推定精度をより向上できる。

次に、本実施形態の効果をシミュレーションによって確認した結果を示す。シミュレーションでは、真の点像分布関数を標準偏差２．５ピクセルのガウス関数として、半径が４０ピクセルの円形開口像を理想画像データとして用いて第１の計算手法及び第２の計算手法を用いて点像分布関数を推定した。このシミュレーションでは円形開口像の最大画素値の１０％を標準偏差とするガウスノイズを加算したものを参照画像データに設定した。

図６において、（ａ）部には真の点像分布関数の画像イメージを示し、（ｂ）部には、従来のウィーナフィルタ法によって推定した比較例１の点像分布関数の画像イメージを示し、（ｃ）部には、本実施形態の第２の計算手法によって推定した点像分布関数の画像イメージを示している。また、図７には、真の点像分布関数、比較例１における点像分布関数、及び本実施形態の第２の計算手法によって推定した点像分布関数の一次元方向のプロファイルを示している。ここでは、比較例１では、実定数ε＝１０００と設定した。この結果により、ウィーナフィルタ法に比較して第２の計算手法のほうが真値に近い良好な関数が得られていることが分かる。この第２の計算手法は、演算時間は比較的長いが、ノイズにロバストな手法であるため、光源強度の弱い光学系に好適である。

図８において、（ａ）部には真の点像分布関数の画像イメージを示し、（ｂ）部には、従来のウィーナフィルタ法（実定数ε＝５００）によって推定した比較例２の点像分布関数の画像イメージを示し、（ｃ）部には、従来のウィーナフィルタ法（実定数ε＝２０００）によって推定した比較例３の点像分布関数の画像イメージを示し、（ｄ）部には、本実施形態の第１の計算手法によって推定した点像分布関数の画像イメージを示している。また、図９には、真の点像分布関数、比較例２，３における点像分布関数、及び本実施形態の第１の計算手法によって推定した点像分布関数の一次元方向のプロファイルを示し、図１０には、図９の一次元プロファイルを拡大して示している。ここでは、第１の計算手法では、実定数ε＝５００、実定数λ＝１０と設定した。この結果により、ウィーナフィルタ法に比較して第１の計算手法のほうが真値に近い良好な関数が得られていることが分かる。特に、実定数εを小さく保ったまま実定数λとして適切な値を与えることで、点像分布関数の先鋭度を損なうことなく良好な推定が行われていることが確認できた。この第１の計算手法は、演算時間がウィーナフィルタ法の高々定数倍程度であり高速である。また、第１の計算手法は光源強度の比較的強い光学系に好適である。

次に、図１１に示すような実験光学系１０１を用いて本実施形態の効果を評価した結果を示す。図１１に示すように、実験光学系１０１としては、本実施形態を想定して、ＬＥＤである光源３とカメラ１１との間に、Ｎｉｋｏｎ社製の４倍対物レンズ９ａと焦点距離８０ｍｍの平凸レンズ９ｂとを組み合わせた結像用レンズを配置したものを用いた。そして、光源３と対物レンズ９ａとの間に直径１ｍｍの円形開口部を有する遮蔽板２１を配置し、光源３から照射された光によって生じる円形開口像をカメラ１１を用いて撮像し、その結果得られた画像データを基に演算装置１３を用いて点像分布関数を計算した。さらには、従来手法であるエッジ法も評価するため、実験光学系１０１を用いて、カッターナイフの刃を撮像した結果得られる画像データを基にエッジ法によって線像分布関数を計算した。エッジ法とは、ナイフエッジあるいは金属ワイヤのエッジ等を撮像した結果得られる画像データの一次元分布を計算することにより線像分布関数を得る手法である。

図１２には、実験光学系１０１を用いて推定された線像分布関数の画像イメージを示し、（ａ）部には本実施形態の第２の計算手法による推定結果を示し、（ａ）部には従来のウィーナフィルタ法（実定数ε＝１）による推定結果を示し、（ｃ）部には本実施形態の第１の計算手法（実定数ε＝１、実定数λ＝３）による推定結果を示している。この結果より、従来法による推定結果においては輝度の高い中心部の周囲に値のばらつきが見られるが、本実施形態の推定結果においては値のばらつきが少なく良好な結果が得られていることが分かる。

図１３には、実験光学系１０１を用いて第２の計算手法によって推定された点像分布関数の一次元プロファイルを、エッジ法によって取得された線像分布関数と比較して示している。この結果から、第２の計算手法による推定結果はエッジ法による計算結果とよく一致していることが明らかになった。

図１４には、実験光学系１０１を用いて第１の計算手法によって推定された点像分布関数の一次元プロファイルを、従来のウィーナフィルタ法（実定数ε＝１）による推定結果、及びエッジ法によって取得された線像分布関数と比較して示している。従来のウィーナフィルタ法ではノイズの影響によりピークの周囲に値のばらつきが見られるが、第１の計算手法では、フーリエ空間における滑らかさが要求されることにより、このような値のばらつきが改善されていることが明らかとなった。また、第１の計算手法による推定結果もエッジ法による計算結果とよく一致していることが分かった。

以上、本発明の種々の実施形態について説明したが、本発明は上記実施形態に限定されるものではなく、各請求項に記載した要旨を変更しない範囲で変形し、又は他のものに適用したものであってもよい。

上記実施形態では、試料Ｓの像とともに試料ホルダ７の円形開口部１９の像を撮像していたが、回転対称な形状の像を撮像できればそれに限定されるものではなく、例えば、試料Ｓの像とともに蛍光ビーズ等の球状の物体の像を撮像してもよいし、円筒形の物体の像を撮像してもよい。シャープなエッジを撮像できるという点では円筒形の物体を撮像することも好ましい。

１…観察システム、３…光源（光学系）、５…集光用レンズ（光学系）、７…試料ホルダ、９…結像用レンズ（光学系）、１１…カメラ（撮像装置）、１３…演算装置、１７…光透過部（保持部）、１９…開口部、Ｓ…試料。

Claims (16)

1. 光学系を用いて試料を撮像した結果得られた試料画像データと、前記光学系を用いて既知の回転対称な形状の像を前記試料と同時に撮像した結果得られた参照画像データとを取得する取得ステップと、
   前記参照画像データと、前記形状の情報とを基に前記光学系の点像分布関数を推定する推定ステップと、
   前記点像分布関数を用いて前記試料画像データをデコンボリューションすることにより出力用試料画像データを生成する生成ステップと、
   前記点像分布関数と前記出力用試料画像データとを出力する出力ステップと、
   を備える点像分布関数推定方法。
2. 前記回転対称な形状は円形である、
   請求項１に記載の点像分布関数推定方法。
3. 前記取得ステップでは、試料を保持する保持部及び光を通過させる円形の開口部を有する試料ホルダと、試料ホルダに向けて光を照射する光源とを用いて、前記試料画像データと前記参照画像データとを取得する、
   請求項２に記載の点像分布関数推定方法。
4. 前記推定ステップでは、前記点像分布関数がゼロ以上であることを拘束条件として前記点像分布関数を推定する、
   請求項１〜３のいずれか１項に記載の点像分布関数推定方法。
5. 前記推定ステップでは、前記参照画像データと前記形状の情報とをフーリエ変換した値と、フーリエ変換後の前記点像分布関数とに関する最小化問題を解くことにより、前記点像分布関数を推定する、
   請求項１〜３のいずれか１項に記載の点像分布関数推定方法。
6. 前記推定ステップでは、前記フーリエ変換した値と、前記フーリエ変換後の前記点像分布関数と、前記点像分布関数のフーリエ空間における滑らかさを示す値とに関する最小化問題を解くことにより、前記点像分布関数を推定する、
   請求項５に記載の点像分布関数推定方法。
7. 前記推定ステップでは、推定した前記点像分布関数を基に線像分布関数を計算し、前記線像分布関数を任意の関数でフィッティングし、フィッティングした前記線像分布関数を用いて点像分布関数を再構成する、
   請求項１〜６のいずれか１項に記載の点像分布関数推定方法。
8. 前記推定ステップでは、前記線像分布関数を撮像面における複数の方向に沿って積分することによって複数の線像分布関数を計算し、前記複数の線像分布関数を任意の関数でフィッティングし、フィッティングした前記複数の線像分布関数を用いて点像分布関数を再構成する、
   請求項７に記載の点像分布関数推定方法。
9. 試料の像及び既知の回転対称な形状の像を同時に導く光学系と、
   前記試料の像を前記光学系を介して撮像して試料画像データを取得すると同時に、前記既知の回転対称な形状の像を前記光学系を介して撮像して参照画像データを取得する撮像装置と、
   前記参照画像データと、前記形状の情報とを基に前記光学系の点像分布関数を推定し、前記点像分布関数を用いて前記試料画像データをデコンボリューションすることにより出力用試料画像データを生成し、前記点像分布関数と前記出力用試料画像データとを出力する演算装置と、
   を備える点像分布関数推定装置。
10. 前記回転対称な形状は円形である、
    請求項９に記載の点像分布関数推定装置。
11. 試料を保持する保持部及び光を通過させる円形の開口部を有する試料ホルダと、
    前記試料ホルダに向けて光を照射する光源と、
    をさらに備え、
    前記撮像装置は、前記光源から前記光が照射された前記試料の像、及び、前記開口部を通過する前記光によって生ずる前記開口部の像を撮像する、
    請求項１０に記載の点像分布関数推定装置。
12. 前記演算装置は、前記点像分布関数がゼロ以上であることを拘束条件として前記点像分布関数を推定する、
    請求項９〜１１のいずれか１項に記載の点像分布関数推定装置。
13. 前記演算装置は、前記参照画像データと前記形状の情報とをフーリエ変換した値と、フーリエ変換後の前記点像分布関数とに関する最小化問題を解くことにより、前記点像分布関数を推定する、
    請求項９〜１１のいずれか１項に記載の点像分布関数推定装置。
14. 前記演算装置は、前記フーリエ変換した値と、前記フーリエ変換後の前記点像分布関数と、前記点像分布関数のフーリエ空間における滑らかさを示す値とに関する最小化問題を解くことにより、前記点像分布関数を推定する、
    請求項１３に記載の点像分布関数推定装置。
15. 前記演算装置は、推定した前記点像分布関数を基に線像分布関数を計算し、前記線像分布関数を任意の関数でフィッティングし、フィッティングした前記線像分布関数を用いて点像分布関数を再構成する、
    請求項９〜１４のいずれか１項に記載の点像分布関数推定装置。
16. 前記演算装置は、前記線像分布関数を撮像面における複数の方向に沿って積分することによって複数の線像分布関数を計算し、前記複数の線像分布関数を任意の関数でフィッティングし、フィッティングした前記複数の線像分布関数を用いて点像分布関数を再構成する、
    請求項１５に記載の点像分布関数推定装置。