JP5965638B2 - 画像処理システムを備えた顕微鏡システム - Google Patents

Description

本発明は、画像処理システムを備えた顕微鏡システムに関する。

一般に、顕微鏡を用いて高倍率で取得される画像は被写界深度が浅い。そこで、例えば１つの標本について互いに異なる焦点面を有する複数の顕微鏡画像を取得し、各顕微鏡画像のうち合焦している画像を組み合わせて全焦点画像や３次元再構築画像を作成する技術が知られている。ここで、合焦している画像の選択、すなわち合焦しているか否かの判断には、画像のコントラストに基づく評価が用いられ得る。

全焦点画像や３次元再構成画像の作成において、画像中の各位置について焦点面が異なる複数の画像から最も合焦している画像を選択して標本の３次元形状を推定した後に、この３次元形状の推定値に対して最適化処理が必要となる。この最適化処理には、画素間の相関に基づき孤立点の推定誤差を低減することが含まれ得る。また、上記の選択ができなかった位置について、標本の形状を推定することが含まれ得る。

全焦点画像の作成について、例えば特許文献１には、全焦点画像を作成した後に回復フィルタによる処理を行うことに係る技術が開示されている。一般に、顕微鏡の光学系は、デジタルカメラの光学系に比べて像の拡大率が大きい。このため、顕微鏡撮影においては、カメラの撮像素子のサンプリング帯域に対して、顕微鏡の光学系の帯域がそれほど高くない場合がある。光学系の帯域は、光学系の開口数や倍率等によって変化し得る。例えば顕微鏡がズーム光学系を有する場合にも、その光学系の帯域は変化する。特許文献１に開示されている技術では、回復フィルタの係数は、光学系の帯域変化を考慮して対物レンズの倍率や開口数といった光学系の設定に応じて決められている。

特開平９−２９８６８２号公報

前記特許文献１に係る技術では、回復フィルタによる処理を行うことで、構造の回復及びエッジの強調が行われる。これは、上記したような標本の３次元形状の推定値に対する最適化処理とは異なる。上記したような最適化処理においても、光学系の帯域を考慮しないと、最適化処理に誤差が含まれ得る。

そこで本発明は、光学系が備える帯域を考慮した簡易な方法で標本の３次元形状の推定値に対する最適化処理を行える画像処理システムを備えた顕微鏡システムを提供することを目的とする。

前記目的を果たすため、本発明の一態様によれば、顕微鏡システムは、顕微鏡と、画像処理システムとを備えた顕微鏡システムであって、前記顕微鏡は、光源と、光学系と、画像取得部とを有し、前記光学系は、対物レンズと、変倍光学系を含む観察光学系とを含み、前記光源は、所定の波長の照明光を射出し、前記画像処理システムは、前記光学系を介して同一の標本における互いに異なる焦点位置を有する複数の画像を取得する取得部と、前記複数の画像のそれぞれに対し、これら画像内の画素毎のコントラスト評価を行い、各々の前記画素のコントラスト評価値を算出するコントラスト評価部と、前記コントラスト評価値に基づいて、夫々の前記画素のうちから合焦している前記画素を選択し、前記標本の３次元形状の候補値を推定する候補値推定部と、前記光学系の情報に基づいて、前記光学系のカットオフ周波数を取得するカットオフ周波数取得部と、前記カットオフ周波数に基づいて、前記候補値に対してデータ補正とデータ補間とのうち少なくとも一方を行い、修正候補値を算出する候補値修正部と、を有し、前記光学系の情報は、前記対物レンズの開口数と、前記観察光学系の総合倍率と、前記照明光の前記波長とを含むことを特徴とする。

本発明によれば、光学系の情報に基づいてカットオフ周波数を取得し、このカットオフ周波数に応じて、３次元形状の候補値に対してデータ補正とデータ補間とのうち少なくとも一方を行うので、顕微鏡の光学系の帯域を考慮した簡易な方法で標本の３次元形状の推定値に対する最適化処理を行える画像処理システムを備えた顕微鏡システムを提供できる。

本発明の各実施形態に係る顕微鏡システムの構成例を示すブロック図。 光学系に応じた周波数とコントラストとの関係を説明するための概略図。 第１の実施形態に係るデータ補正部における処理例を示すフローチャート。 コアリング処理を説明するための図であり、形状候補値に相当する原信号の例を表す図。 コアリング処理を説明するための図であり、移動平均と閾値との例を表す図。 コアリング処理を説明するための図であり、コアリング処理の結果例を表す図。 第１の実施形態に係るデータ補間部における処理例を示すフローチャート。

［第１の実施形態］
本発明の第１の実施形態について図面を参照して説明する。図１は、本実施形態に係る顕微鏡システム１００の構成例の概略を示す。この図に示すように、顕微鏡システム１００は、顕微鏡１１０を有する。顕微鏡１１０は、例えばデジタル顕微鏡である。顕微鏡１１０は、ＬＥＤ光源１１１と、照明光学系１１２と、光路制御素子１１３と、対物レンズ１１４と、標本面１１５と、観察光学系１１８と、撮像面１１９と、画像取得部１２０とを有する。観察光学系１１８は、ズーム光学系１１６と、結像光学系１１７とを有する。観察光路には、標本面１１５から撮像面１１９に向かって、対物レンズ１１４と、光路制御素子１１３と、ズーム光学系１１６と、結像光学系１１７とが、この順に配置されている。

ＬＥＤ光源１１１から射出された照明光は、照明光学系１１２を介して、光路制御素子１１３に入射する。光路制御素子１１３は、照明光を観察光路上の対物レンズ１１４に向けて反射する。照明光は、対物レンズ１１４を介して標本面１１５に配置された標本に照射される。

標本に照明光が照射されると、標本から観察光が生じる。ここで、観察光は反射光や蛍光等である。観察光は、光路制御素子１１３に入射する。光路制御素子１１３は、照明光の場合と異なり、観察光を透過させて、観察光をズーム光学系１１６及び結像光学系１１７を有する観察光学系１１８に入射させる。このように光路制御素子１１３は、入射光の特性に応じて、入射光を反射又は透過させる光学素子である。光路制御素子１１３には、例えば、入射光の偏光方向に応じて入射光を反射又は透過させるワイヤーグリッドや偏光ビームスプリッタ（ＰＢＳ）等の偏光素子が用いられる。また光路制御素子１１３には、例えば、入射光の周波数に応じて入射光を反射又は透過させるダイクロイックミラー等が用いられる。

観察光学系１１８は、観察光を撮像面１１９に集光させ、標本の像を撮像面１１９に結像させる。画像取得部１２０は、撮像面１１９に結像した画像に基づく画像信号を生成し、この画像信号を顕微鏡画像として出力する。

顕微鏡システム１００は、入力部２１０と、コントローラ２２０と、画像入力部２３０と、記憶部２４０と、コントラスト評価部２５０と、形状候補推定部２６０と、光学帯域情報取得部２７０と、パラメータ決定部２８０と、データ修正部２９０と、３Ｄ形状推定部３００と、画像合成部３１０とを有する。データ修正部２９０は、データ補正部２９２と、データ補間部２９４と、で構成された画像処理システム２００を含む。

入力部２１０は、例えば顕微鏡システム１００の操作に係るユーザの指示を受け取る。入力部２１０は、受け取ったユーザの指示をコントローラ２２０に出力する。コントローラ２２０は、顕微鏡１１０の各部の動作を制御する。コントローラ２２０は、顕微鏡１１０の制御において、入力部２１０から入力されたユーザの指示を用いる。本実施形態では、顕微鏡システム１００は、同一標本について、異なる焦点面において撮像した複数の顕微鏡画像を取得する。このため、コントローラ２２０は、顕微鏡１１０の光学系を制御して焦点面を徐々に変化させながら、各焦点面における標本の画像を画像取得部１２０に取得させる。具体的には、例えばコントローラ２２０は、顕微鏡１１０のステージの高さ、対物レンズの高さ又はフォーカスレンズの位置を変化させながら、画像取得部１２０に画像を取得させる。

コントローラ２２０は、光学帯域情報取得部２７０からの要求に応じて顕微鏡１１０の各種光学系のパラメータを光学帯域情報取得部２７０に出力する。光学系のパラメータには、例えば、対物光学系、照明光学系、ズーム光学系、結像光学系のパラメータが含まれる。なお、コントローラ２２０は、顕微鏡１１０の光学系を制御しているので、これら光学系のパラメータに係る情報を有している。

画像入力部２３０は、画像取得部１２０から顕微鏡画像を取得する。画像入力部２３０は、この顕微鏡画像を記憶部２４０に出力する。記憶部２４０は、入力された顕微鏡画像を一時記憶する。本実施形態では、画像取得部１２０はコントローラ２２０の制御下で同一の標本について焦点面が異なる複数の顕微鏡画像を取得するので、画像入力部２３０は焦点面が異なる複数の顕微鏡画像を取得し、記憶部２４０にはこれら複数の顕微鏡画像が記憶される。記憶部２４０は、コントラスト評価部２５０からの要求に応じて、記憶した顕微鏡画像をコントラスト評価部２５０に出力する。また、記憶部２４０は、画像合成部３１０からの要求に応じて、記憶した顕微鏡画像を画像合成部３１０に出力する。

コントラスト評価部２５０は、記憶部２４０から焦点面が異なる複数の顕微鏡画像を取得する。コントラスト評価部２５０は、複数の顕微鏡画像の各画素について高周波成分の強度を評価し、コントラスト評価値を算出する。コントラスト評価部２５０は、算出した各顕微鏡画像の各画素に係るコントラスト評価値を形状候補推定部２６０に出力する。

形状候補推定部２６０は、コントラスト評価部２５０から入力されたコントラスト評価値に基づいて、複数の顕微鏡画像のそれぞれについて画像中の各画素についての合焦を評価する。さらに、形状候補推定部２６０は、画像中の各画素について焦点面が異なる複数の画像のうち最も合焦している画像を選択する。形状候補推定部２６０は、最も合焦している画像が撮像された際の焦点面の情報をコントローラ２２０から取得し、その情報に基づいて画像中の各画素に対応する標本の高さを推定し、標本の形状に係る情報である形状候補値を算出する。コントラスト評価値に基づいて標本の高さを推定できなかった画素については、形状候補推定部２６０は、当該画素に対応する形状候補値を、推定不能を表す値とする。形状候補推定部２６０は、算出した形状候補値をデータ修正部２９０内のデータ補正部２９２に出力する。

光学帯域情報取得部２７０は、コントローラ２２０から各種光学系のパラメータを取得する。光学帯域情報取得部２７０は、取得した光学系のパラメータをパラメータ決定部２８０に出力する。パラメータ決定部２８０は、光学帯域情報取得部２７０から入力された光学系のパラメータに基づいて、後述する補正パラメータを算出する。パラメータ決定部２８０は、算出した補正パラメータをデータ修正部２９０内のデータ補正部２９２に出力する。また、パラメータ決定部２８０は、光学系のパラメータに基づいて、後述する補間パラメータを算出する。パラメータ決定部２８０は、算出した補間パラメータをデータ修正部２９０内のデータ補間部２９４に出力する。

データ修正部２９０が備えるデータ補正部２９２は、形状候補推定部２６０から入力された形状候補値に対してノイズコアリングを行い、形状候補値のノイズを除去する。データ補正部２９２は、ノイズコアリングを行う際に、後に詳述するように、パラメータ決定部２８０から入力された補正パラメータを用いる。データ補正部２９２は、ノイズが除去された形状候補値であるノイズ除去形状候補値を、データ補間部２９４に出力する。

データ修正部２９０が備えるデータ補間部２９４は、データ補正部２９２から入力されたノイズ除去形状候補値のうち、推定不能を表す値となっている画素について、データの補間を行う。データ補間部２９４は、データを補間する際に、後に詳述するようにパラメータ決定部２８０から入力された補間パラメータを用いる。データ補間部２９４は、ノイズが除去され推定不能な画素の値が補間された形状候補値である補間形状候補値を、３Ｄ形状推定部３００に出力する。

３Ｄ形状推定部３００は、データ補間部２９４から入力された補間形状候補値に基づいて、高さ情報の最適化を行い、標本の３次元形状の推定値を確定する。３Ｄ形状推定部３００は、確定した標本の３次元形状を画像合成部３１０に出力する。画像合成部３１０は、３Ｄ形状推定部３００から入力された標本の３次元形状と、記憶部２４０から取得した複数の顕微鏡画像とに基づいて、焦点面が異なる複数の顕微鏡画像を合成して合成画像を作成する。この合成画像は、例えば３次元再構成画像や、全焦点画像等である。画像合成部３１０は、作成した合成画像を、例えば表示させるために表示部に出力したり、例えば記憶させるために記憶装置に出力したりする。

本実施形態に係る顕微鏡システム１００の動作を説明する。標本は、標本面１１５に配置される。コントローラ２２０は、顕微鏡１１０を制御する。コントローラ２２０は、例えば標本面１１５の位置を光軸方向に徐々に変化させ、標本に対する光学系の焦点面を徐々に変化させる。具体的には、例えばコントローラ２２０は、顕微鏡１１０のステージの高さ、対物レンズの高さ又はフォーカスレンズの位置を変化させる。このとき、コントローラ２２０は、各焦点面における標本の顕微鏡画像を画像取得部１２０に順次取得させる。画像入力部２３０は、画像取得部１２０から各焦点面における標本の顕微鏡画像を取得する。画像入力部２３０は、取得した顕微鏡画像を記憶部２４０に出力する。記憶部２４０は、画像入力部２３０から入力された複数の顕微鏡画像を記憶する。

コントラスト評価部２５０は、記憶部２４０から焦点面が異なる複数の顕微鏡画像を取得する。コントラスト評価部２５０は、複数の顕微鏡画像の各画素について高周波成分の強度を評価し、コントラスト評価値を算出する。コントラスト評価部２５０は、算出した各顕微鏡画像の各画素に係るコントラスト評価値を形状候補推定部２６０に出力する。

形状候補推定部２６０は、コントラスト評価部２５０から入力されたコントラスト評価値に基づいて、複数の顕微鏡画像のそれぞれについて画像中の各画素についての合焦を評価する。形状候補推定部２６０は、例えばコントラストが高いほど合焦度が高いと評価する。さらに、形状候補推定部２６０は、画像中の各画素について焦点面が異なる複数の画像のうち最も合焦している画像を選択する。形状候補推定部２６０は、最も合焦している画像が撮像された際の焦点面の情報をコントローラ２２０から取得する。形状候補推定部２６０は、コントローラ２２０から取得した情報に基づいて画像中の各画素に対応する標本の高さを推定し、標本の形状に係る情報である形状候補値Ｐ（ｉ，ｊ）を算出する。ここで（ｉ，ｊ）は座標を表す。形状候補値Ｐ（ｉ，ｊ）は、例えば座標（ｉ，ｊ）における標本の高さを表す。コントラスト評価値に基づいて標本の高さを推定できなかった場合は、形状候補推定部２６０は、当該画素に対応する形状候補値Ｐ（ｉ，ｊ）を、推定不能を表す値とする。形状候補推定部２６０は、算出した形状候補値Ｐ（ｉ，ｊ）をデータ補正部２９２に出力する。

コントローラ２２０は、光学帯域情報取得部２７０の要求に応じて顕微鏡画像取得時の顕微鏡１１０の各種光学系に係る情報を光学帯域情報取得部２７０に出力する。ここで、光学帯域情報取得部２７０は、光学系のパラメータである対物レンズ１１４の開口数ＮＡ_ｏｂｊと観察光学系１１８の総合倍率である像倍率βと、ＬＥＤ光源１１１が射出する照明光の波長λとを取得する。光学帯域情報取得部２７０は、これらパラメータ決定部２８０に出力する。

パラメータ決定部２８０は、取得した光学系パラメータを用いて、例えば下記式（１）に基づいてカットオフ周波数Ｆｃを算出する。
ここで、ＮＡ_ｏｂｊは開口数、βは倍率、λは波長を表している。なお、開口数ＮＡ_ｏｂｊについては、厳密には対物レンズの開口数のみならず照明側のコンデンサレンズの開口数も関係する。しかしながらここでは、照明がインコヒーレントであり、標本が散乱体であり、照明側の開口数を考慮しなくてもよい場合を想定している。

なお、像倍率βが低い場合、すなわち低倍率のときは、高周波成分が多く含まれ、像倍率βが高い場合、すなわち高倍率のときは、高周波成分は減少する。上記式（１）は、このこともよく表している。なお、波長λは、照明光の波長によらず、例えば５５０ｎｍといった所定の値としてもよい。

レンズの変調伝達関数（Ｍｏｄｕｌａｔｉｏｎ Ｔｒａｎｓｆｅｒ Ｆｕｎｃｔｉｏｎ；ＭＴＦ）を考慮しない場合における、像倍率βに応じた周波数とコントラストとの関係の概略を図２に示す。この図では、βが１の場合と５の場合とを示している。ここで、コントラストＣは、下記式（２）で表される。
ここで、Ｌ０は明暗のチャートにおける周波数０における明の輝度レベルを、Ｄ０は周波数０における暗の輝度レベルを、Ｌｆは周波数ｆにおける明領域の輝度レベルを、Ｄｆは周波数ｆにおける暗領域の輝度レベルをそれぞれ表す。なお、コントラストが０となる周波数が、カットオフ周波数Ｆｃと定義されている。

パラメータ決定部２８０は、カットオフ周波数Ｆｃに基づいて、データ補正部２９２で用いる補正パラメータｍ，ｎ，ｗ（ｋ，ｌ）を決定する。パラメータ決定部２８０は、例えばカットオフ周波数Ｆｃと補正パラメータｍ，ｎ，ｗ（ｋ，ｌ）との関係を示すルックアップテーブルを記憶している。パラメータ決定部２８０は、このルックアップテーブルを参照して、カットオフ周波数Ｆｃに基づいて補正パラメータｍ，ｎ，ｗ（ｋ，ｌ）を決定する。ここで、補正パラメータｍ，ｎは、カットオフ周波数Ｆｃが低いほど大きな値となる。また、ｗ（ｋ，ｌ）としては、後述の式（３）においてｍ，ｎの値が大きいときに重みが減少しないような関数が与えられる。パラメータ決定部２８０は、決定した補正パラメータｍ，ｎ，ｗ（ｋ，ｌ）をデータ補正部２９２に出力する。

また、パラメータ決定部２８０は、カットオフ周波数Ｆｃに基づいて、データ補間部２９４で用いる補間パラメータσ_ｋ，σ_ｌを決定する。パラメータ決定部２８０は、例えばカットオフ周波数Ｆｃと補間パラメータσ_ｋ，σ_ｌとの関係を示すルックアップテーブルを記憶している。パラメータ決定部２８０は、このルックアップテーブルを参照して、カットオフ周波数Ｆｃに基づいて補間パラメータσ_ｋ，σ_ｌを決定する。ここで、補間パラメータσ_ｋ，σ_ｌは、カットオフ周波数Ｆｃが低いほど大きな値となる。パラメータ決定部２８０は、決定した補間パラメータσ_ｋ，σ_ｌをデータ補間部２９４に出力する。

データ補正部２９２は、形状候補値Ｐ（ｉ，ｊ）に対してノイズ又は孤立点の除去を行うノイズ・孤立点除去処理を行う。本実施形態では、ノイズ・孤立点除去処理は、コアリング処理によって行われる。ノイズ・孤立点除去処理を表すフローチャートを図３に示す。

ステップＳ１１０において、データ補正部２９２は、形状候補値Ｐ（ｉ，ｊ）を読み込む。本実施形態では、顕微鏡画像のサイズは、横が０からｐまでのｐ＋１画素、縦が０からｑまでのｑ＋１画素であるものとする。ステップＳ１２０において、データ補正部２９２は、補正パラメータｍ，ｎ，ｗ（ｋ，ｌ）を読み込む。

本実施形態では、図３に示すように、ステップＳ１３１乃至ステップＳ１３４において、画像の全画素に対応する形状候補値Ｐ（ｉ，ｊ）について順に以下の処理を行う。ステップＳ１３１において、データ補正部２９２は、（ｉ，ｊ）を含む領域の基準値Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）を下記式（３）に基づいて算出する。
上記式（３）のとおり、基準値Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）は、この領域における平均値を表す。また、上記式（３）において、パラメータ決定部２８０が決定した補正パラメータｍ，ｎ，ｗ（ｋ，ｌ）が用いられている。すなわち、上記式（３）は、カットオフ周波数Ｆｃに応じて変化する。

ステップＳ１３２において、データ補正部２９２は、形状候補値Ｐ（ｉ，ｊ）と基準値Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）との差が所定の閾値よりも小さいか否かを判定する。形状候補値Ｐ（ｉ，ｊ）と基準値Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）との差が所定の閾値Ｔｈ_ｒ−１よりも小さいとき、つまり、（｜Ｐ（ｉ，ｊ）−Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）｜＜Ｔｈ_ｒ−１）の判断により処理はステップＳ１３４に進められる。
なお、ここで閾値Ｔｈ_ｒ−１は、基準値の誤差範囲内に入るか否かの基準等、経験則に基づいて決定される。

一方、形状候補値Ｐ（ｉ，ｊ）と基準値Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）との差が所定の閾値よりも小さくないとき、ステップＳ１３３において、データ補正部２９２は、形状候補値Ｐ（ｉ，ｊ）は孤立点であるか否かを判定する。形状候補値Ｐ（ｉ，ｊ）は孤立点であるとき、処理はステップＳ１３４に進められる。
ここで、孤立点の判定の具体的な方法は、（｜Ｐ（ｉ，ｊ）−Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）｜＞Ｔｈ_ｒ−２）により判断され、Ｔｈ_ｒ−２は所定の複数画素領域における分散値から設定されるものである。より具体的には、例えば、分散値がσであった場合、Ｔｈ_ｒ−２は±２σとして設定し、孤立点を判断する。

ステップＳ１３４において、データ補正部２９２は、形状候補値Ｐ（ｉ，ｊ）の値を基準値Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）に設定する。これらステップＳ１３１乃至ステップＳ１３４の処理が、全ての画素について行われる。すなわち、本処理は、上記の所定の閾値をΔＴとし、処理後の形状候補値をノイズ除去形状候補値Ｐ´（ｉ，ｊ）としたときに、下記式（４）で表される。

本実施形態で用いられるコアリング処理の概念を図４Ａ，４Ｂ，４Ｃを参照して説明する。図４Ａは、形状候補値Ｐ（ｉ，ｊ）に相当する原信号を表す。この原信号に対して、上記式（３）で算出される平均値に相当する移動平均が、図４Ｂに一点鎖線で示されている。また、移動平均に対して上記の所定の閾値ΔＴに相当する閾値を加算又は減算した値が図４Ｂに破線で示されている。この場合、上記式（４）で表されるように、図４Ｂにおいて原信号が２本の破線の間に位置する場合には、原信号は、一点鎖線で示された移動平均に置き換えられる。その結果、図４Ｃに示すような結果が得られる。なお、この図において丸印は移動平均に置き換えられた値を示す。このように、コアリング処理は、小振幅信号と判断される変動成分を抑制し、エラーを削除する効果を有する。

データ補正部２９２は、形状候補値Ｐ（ｉ，ｊ）に対して図３を参照して説明したノイズ・孤立点除去処理を施したもの、すなわち、ノイズ除去形状候補値Ｐ´（ｉ，ｊ）を、データ補間部２９４に出力する。

データ補間部２９４は、データ補正部２９２から入力されたノイズ除去形状候補値Ｐ´（ｉ，ｊ）について、推定不能を表す値となっているデータについて補間を行う。推定不能とは、コントラスト評価部２５０が算出したコントラスト評価値に基づいて形状候補推定部２６０が形状候補値Ｐ（ｉ，ｊ）を算出する際に、形状候補推定部２６０が画像の合焦を特定できなかったことを意味する。すなわち推定不能とは、当該画素について複数の顕微鏡画像のうち何れの顕微鏡画像のコントラスト評価値も、所定の合焦を示す条件に適合しなかったことを表す。

ノイズ除去形状候補値Ｐ´（ｉ，ｊ）について推定不能を表す値の周囲が推定不能でない場合、すなわち、例えば５画素×５画素の領域のうち１画素のみが推定不能である場合には、データ補間部２９４は、当該推定不能のデータについて近隣のデータを用いて補間を行う。この際、データ補間部２９４は、データの補間に例えば双１次補間（ｂｉｌｉｎｅａｒ ｉｎｔｅｒｐｏｌａｔｉｏｎ）や、双３次補間（ｂｉｃｕｂｉｃ ｉｎｔｅｒｐｏｌａｔｉｏｎ）を用いることができる。

一方、ノイズ除去形状候補値Ｐ´（ｉ，ｊ）について推定不能を表す値が連続して存在する場合には、データ補間部２９４は、近隣のデータとの相関を表す関数に基づいて、推定不能となっているデータの補間を行う。すなわち、推定不能となっている箇所の周辺の分布を仮定することで、当該箇所の値を推定する。本実施形態では、補間にはカーネル回帰法が用いられる。この際、データ補間部２９４は、パラメータ決定部２８０から入力された補間パラメータσ_ｋ，σ_ｌを用いる。この補間の処理例を図５に示すフローチャートを参照して説明する。

ステップＳ２１０において、データ補間部２９４は、ノイズ除去形状候補値Ｐ´（ｉ，ｊ）を読み込む。ステップＳ２２０において、データ補間部２９４は、補間パラメータσ_ｋ，σ_ｌを読み込む。続いてデータ補間部２９４は、補間データＱ（ｉ，ｊ）の算出を行う。ここで、補間データＱ（ｉ，ｊ）は、下記式（５）で表される。
ここで、Ｎはサンプリング数を表し、下記式（６）で表される。
また、Ｃ（ｋ，ｌ）は、下記式（７）で表される。
式（７）に示すようにＣ（ｋ，ｌ）は、補間パラメータσ_ｋ，σ_ｌに応じて決定される。また、Ａは、所定の係数である。

ステップＳ２３１において、データ補間部２９４は、変数Ａを更新する。ステップＳ２３２において、データ補間部２９４は、ノイズ除去形状候補値Ｐ´（ｉ，ｊ）に、上記式（５）乃至（７）に基づいてガウスカーネルを重畳する。ステップＳ２３３において、データ補間部２９４は、ステップＳ２３２で得られた値が所定の収束条件を満たしているか否かを判定する。収束条件を満たしていれば、処理はステップＳ２４０に進められる。一方、満たしていなければ、所定の回数であるＤ回まで、ステップＳ２３１乃至ステップＳ２３３の処理を繰り返す。すなわち、収束条件を満たすまで、ステップＳ２３１で変数Ａの値を変更しながら、ステップＳ２３２で各変数Ａに対する補間データＱ（ｉ，ｊ）を算出し、ステップＳ２３３で算出した補間データＱ（ｉ，ｊ）が収束条件を満たすか否か判定する。

ステップＳ２３３における判定で、補間データＱ（ｉ，ｊ）が収束条件を満たしていると判定されたら、ステップＳ２４０において、データ補間部２９４は、条件を満たした補間データＱ（ｉ，ｊ）に基づいて拡張データを生成する。ステップＳ２５０において、データ補間部２９４は、ノイズ除去形状候補値Ｐ´（ｉ，ｊ）のうち推定不能であったデータに生成した拡張データを割り当てて、補間形状候補値Ｐ´´（ｉ，ｊ）を生成する。データ補間部２９４は、生成した補間形状候補値Ｐ´´（ｉ，ｊ）を３Ｄ形状推定部３００に出力する。

３Ｄ形状推定部３００は、データ補間部２９４から入力された補間形状候補値Ｐ´´（ｉ，ｊ）に基づいて、高さ情報の最適化を行い、標本の３次元形状を推定する。３Ｄ形状推定部３００は、推定した標本の３次元形状を画像合成部３１０に出力する。

画像合成部３１０は、３Ｄ形状推定部３００から入力された標本の３次元形状と、記憶部２４０から取得した複数の顕微鏡画像とに基づいて、焦点面が異なる複数の顕微鏡画像を合成して合成画像を作成する。例えば合成画像が３次元再構成画像であれば、３次元形状とその３次元形状の各部に係る合焦している画像を合成させた合成画像を作成する。また、例えば合成画像が全焦点画像であれば、各画素についてその高さに相当する焦点面を有する顕微鏡画像から抽出した画像を組み合わせて、全ての画素について合焦している画像を合成する。画像合成部３１０は、作成した合成画像を表示部や記憶装置に出力する。通常の顕微鏡画像は被写界深度が浅いため、被写界深度よりも大きな高さがある標本の像は、通常の顕微鏡画像によってはユーザに認識されにくい。これに対して３次元再構成画像や全焦点画像によれば、被写界深度よりも大きな高さがある標本の像が、容易にユーザに認識され得る。

このように、例えば画像入力部２３０は、光学系を介して同一の標本における互いに異なる焦点位置を有する複数の画像を取得する取得部として機能する。例えばコントラスト評価部２５０は、複数の画像のそれぞれに対し、画素毎のコントラスト評価を行い各画素のコントラスト評価値を算出するコントラスト評価部として機能する。例えば形状候補推定部２６０は、コントラスト評価値に基づいて、各々の画素のうちから合焦している画素を選択し、標本の３次元形状の候補値を推定する候補値推定部として機能する。例えば光学帯域情報取得部２７０は、光学系の情報に基づいて、光学系のカットオフ周波数を取得するカットオフ周波数取得部として機能する。例えばデータ修正部２９０は、カットオフ周波数に基づいて、候補値に対してデータ補正とデータ補間とのうち少なくとも一方を行い、修正候補値を算出する候補値修正部として機能する。例えばデータ補正部２９２は、候補値が表す局所領域の値の相関を用いて修正候補値を算出する修正候補値算出部として機能する。例えば画像合成部３１０は、修正候補値と複数の画像とに基づいて全焦点画像を作成する全焦点画像作成部として機能する。また、例えば画像合成部３１０は、修正候補値と複数の画像とに基づいて、３次元再構成画像を作成する３次元再構成画像作成部として機能する。例えば３Ｄ形状推定部３００は、修正候補値に最適化を施して、標本の３次元形状を推定する３次元形状推定部として機能する。

本実施形態によれば、データ補正部２９２における処理の結果、顕微鏡システム１００におけるノイズ及び推定処理に起因するエラーが効果的に低減させられる。この際、顕微鏡１１０の光学系のカットオフ周波数Ｆｃに基づいて補正パラメータｍ，ｎが決定されている。ここで、補正パラメータｍ，ｎはカットオフ周波数Ｆｃが低いほど大きな値となっているので、上記式（３）において、カットオフ周波数Ｆｃが低いほど広い領域の形状候補値Ｐ（ｉ，ｊ）に基づいて基準値Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）が算出され、カットオフ周波数Ｆｃが高いほど狭い領域の形状候補値Ｐ（ｉ，ｊ）に基づいて基準値Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）が算出される。すなわち、顕微鏡１１０の光学系のカットオフ周波数Ｆｃに応じた最適な基準値Ｐ_ａｖｅ（ｉ，ｊ，ｍ，ｎ）が算出されることになる。その結果、光学系のカットオフ周波数Ｆｃを考慮しない場合に比べて、精度よくノイズが低減させられ得る。すなわち、形状候補値Ｐ（ｉ，ｊ）が過度に平滑化されることも無く、ノイズが多い場合でも入力信号が過度に高周波信号であると評価されることも無い。

また、データ補間部２９４における処理では、近隣のデータの相関が仮定される際に画像信号が本来有している周波数帯域の情報が用いられる。すなわち、周波数帯域に応じた最適化されたガウスカーネルが生成され、コントラスト評価値に基づいて推定不能であった位置における標本の高さの値が推定され得る。この際、光学系のカットオフ周波数Ｆｃに基づいて補間パラメータσ_ｋ，σ_ｌが与えられるので、これらの値を変化させながら収束値が探索される場合に比べて、計算量が減少するため処理速度が上昇し、かつ、計算結果が誤った値に集束することが防止され得る。ここで、補間パラメータσ_ｋ，σ_ｌはカットオフ周波数Ｆｃが低いほど大きな値となるので、上記式（５）において、カットオフ周波数Ｆｃが低いほど広い領域のノイズ除去形状候補値Ｐ´（ｉ，ｊ）に基づいて補間データＱ（ｉ，ｊ）が算出され、カットオフ周波数Ｆｃが高いほど狭い領域のノイズ除去形状候補値Ｐ´（ｉ，ｊ）に基づいて補間データＱ（ｉ，ｊ）が算出される。すなわち、過度に平滑化されることなく、エッジ構造の評価も適正に行われ、ノイズが多い場合でも入力信号が過度に高周波信号と評価されることも無い。

なお、顕微鏡１１０の光学系にズーム光学系が含まれている場合、ズーム光学系の焦点距離に応じて開口数が変化し顕微鏡画像の帯域が変化するので、本実施形態は特に効果を奏する。なお、上記の各式は一例であり、上述の効果が得られれば、これらの式に限定されず、他の式が用いられてももちろんよい。また、上述の実施形態では、画素毎に処理が行われているが、複数の画素からなる領域ごとに処理が行われてもよい。

［第１の実施形態の変形例］
第１の実施形態の変形例について説明する。ここでは、第１の実施形態との相違点について説明し、同一の部分については同一の符号を付してその説明を省略する。第１の実施形態におけるデータ補間部２９４による処理では、図５を参照して説明したステップＳ２３１乃至ステップＳ２３３のループ処理において、上記式（７）におけるσ_ｋ，σ_ｌには補間パラメータσ_ｋ，σ_ｌが設定され、これらの値は変化させない。

これに対して本変形例では、σ_ｋ，σ_ｌもステップＳ２３１において変化させられながら収束値が探索される。このため本変形例では、パラメータ決定部２８０は、補間パラメータσ_ｋ，σ_ｌとして、取り得る範囲又は確率密度関数をデータ補間部２９４に出力する。データ補間部２９４は、ステップＳ２３１乃至ステップＳ２３３のループ処理において、パラメータ決定部２８０から入力された補間パラメータσ_ｋ，σ_ｌとしての取り得る範囲又は確率密度関数に基づいて、σ_ｋ，σ_ｌも変化させながら収束値を探索する。その他の動作は第１の実施形態と同様である。

本変形例によれば、第１の実施形態よりも処理量は増加するが、補間データＱ（ｉ，ｊ）は、第１の実施形態よりも適した収束値に集束し得る。本変形例においてもパラメータ決定部２８０は、光学系のカットオフ周波数Ｆｃに基づいて補間パラメータσ_ｋ，σ_ｌとしての取り得る範囲又は確率密度関数を決定するので、第１の実施形態と同様の効果が得られる。

［第２の実施形態］
本発明の第２の実施形態について説明する。ここでは、第１の実施形態との相違点について説明し、同一の部分については、同一の符号を付してその説明を省略する。本実施形態では、データ補正部２９２においてノイズ除去のためにバイラテラルフィルタが用いられる。本実施形態で用いられるバイラテラルフィルタを表す式を、下記式（８）に示す。
ここで、Ｃ（ｋ，ｌ）は距離相関を規定する因子であり、Ｓ（Ｐ_１−Ｐ_２）は異なる画素間における画素レベルの差に起因する相関を規定する因子である。Ｃ（ｋ，ｌ）及びＳ（Ｐ_１−Ｐ_２）にどのような分布関数を用いるかによって、生成される画像の鮮鋭度、Ｓ／Ｎ比が変化する。

本実施形態では、例えばＣ（ｋ，ｌ）及びＳ（Ｐ_１−Ｐ_２）をガウス分布に基づく関数とする。すなわち、Ｃ（ｋ，ｌ）を例えば下記式（９）とする。
ここで、σ_ｋ，σ_ｌは補正パラメータであり、Ａは所定の定数である。補正パラメータσ_ｋ，σ_ｌは、第１の実施形態における補間パラメータσ_ｋ，σ_ｌと同じである。また、Ｓ（Ｐ_１−Ｐ_２）を例えば下記式（１０）とする。
ここで、σ_Ｐは補正パラメータであり、Ｂは所定の定数である。本実施形態では、パラメータ決定部２８０は、補正パラメータσ_Ｐをも光学系のカットオフ周波数Ｆｃに基づいてルックアップテーブルを参照して決定する。ここで、補正パラメータσ_Ｐは、カットオフ周波数Ｆｃが低いほど大きな値となる。

第１実施形態の場合と同様に、光学系のカットオフ周波数Ｆｃの情報が取得されれば、顕微鏡画像が本来持っている鮮鋭度が推定され得る。例えば、カットオフ周波数Ｆｃが低い場合には、長距離の相関が重視されるように、Ｃ（ｋ，ｌ）が設定され、さらに、近隣のデータで急激な段差が生じないという仮定に基づいてＳ（Ｐ_１−Ｐ_２）が設定される。このように、例えばＳ（Ｐ_１−Ｐ_２）は、離れた２点の値の相関である第１の相関として機能する。例えばＣ（ｋ，ｌ）は、距離による相関である第２の相関として機能する。

本実施形態では、近隣データの相関が仮定される際に、顕微鏡画像が本来持っている周波数帯域の情報が用いられる。この近隣データの相関に基づいて、バイラテラルフィルタが設定される。その結果、本実施形態によれば、形状候補値Ｐ（ｉ，ｊ）のノイズ及びエラーが効果的に低減させられたノイズ除去形状候補値Ｐ´（ｉ，ｊ）が取得され得る。

なお、本実施形態においても第１の実施形態の第１の変形例と同様に、補正パラメータσ_ｋ，σ_ｌ，σ_Pを確立密度変数としてもよい。この場合も本実施形態と同様の効果が得られる。

［第２の実施形態の第１の変形例］
第２の実施形態の第１の変形例について説明する。ここでは、第２の実施形態との相違点について説明し、同一の部分については、同一の符号を付してその説明を省略する。本変形例では、データ補正部２９２においてノイズ除去のためにトリラテラルフィルタが用いられる。本実施形態で用いられるトリラテラルフィルタを表す式を、下記式（１１）に示す。
ここで、ＰΔ（ｉ，ｊ，ｋ，ｌ）は、下記式（１２）で表される。
また、Ｎ（ｉ，ｊ，ｋ，ｌ）は、下記式（１３）で表される。
ここで、Ｑ（ｉ，ｊ）は勾配ベクトルを表し、下記式（１４）で表される。
ここで、Ｐ_ｆ（ｉ，ｊ，ｋ，ｌ）は、下記式（１５）で表される。
ここで、Ｑ（ｉ，ｊ）_ｉは勾配の水平成分を表し、Ｑ（ｉ，ｊ）_ｊは勾配の垂直成分を表している。

このトリラテラルフィルタは、第２の実施形態で用いられるバイラテラルフィルタを勾配∇Ｐ（ｉ，ｊ）に適用したものである。∇Ｐ（ｉ，ｊ）を導入することにより、インパルスノイズ、すなわち、孤立変動成分を強く抑制することができる。

本変形例でも、第２の実施形態の場合と同様に、光学系のカットオフ周波数Ｆｃに応じて決定されたＣ（ｋ，ｌ）及びＳ（Ｐ_１−Ｐ_２）が用いられる。その結果、第２の実施形態の場合と同様の効果が得られる。

なお、本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記実施形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除しても、発明が解決しようとする課題の欄で述べられた課題が解決でき、かつ、発明の効果が得られる場合には、この構成要素が削除された構成も発明として抽出され得る。さらに、異なる実施形態にわたる構成要素を適宜組み合わせてもよい。

１００…顕微鏡システム、１１０…顕微鏡、１１１…ＬＥＤ光源、１１２…照明光学系、１１３…光路制御素子、１１４…対物レンズ、１１５…標本面、１１６…ズーム光学系、１１７…結像光学系、１１８…観察光学系、１１９…撮像面、１２０…画像取得部、２００…画像処理システム、２１０…入力部、２２０…コントローラ、２３０…画像入力部、２４０…記憶部、２５０…コントラスト評価部、２６０…形状候補推定部、２７０…光学帯域情報取得部、２８０…パラメータ決定部、２９０…データ修正部、２９２…データ補正部、２９４…データ補間部、３００…３Ｄ形状推定部、３１０…画像合成部。

Claims (14)

1. 顕微鏡と、画像処理システムとを備えた顕微鏡システムであって、
   前記顕微鏡は、
   光源と、光学系と、画像取得部とを有し、
   前記光学系は、対物レンズと、変倍光学系を含む観察光学系とを含み、
   前記光源は、所定の波長の照明光を射出し、
   前記画像処理システムは、
   前記光学系を介して同一の標本における互いに異なる焦点位置を有する複数の画像を取得する取得部と、
   前記複数の画像のそれぞれに対し、これら画像内の画素毎のコントラスト評価を行い、各々の前記画素のコントラスト評価値を算出するコントラスト評価部と、
   前記コントラスト評価値に基づいて、夫々の前記画素のうちから合焦している前記画素を選択し、前記標本の３次元形状の候補値を推定する候補値推定部と、
   前記光学系の情報に基づいて、前記光学系のカットオフ周波数を取得するカットオフ周波数取得部と、
   前記カットオフ周波数に基づいて、前記候補値に対してデータ補正とデータ補間とのうち少なくとも一方を行い、修正候補値を算出する候補値修正部と、
   を有し、
   前記光学系の情報は、前記対物レンズの開口数と、前記観察光学系の総合倍率と、前記照明光の前記波長とを含む
   顕微鏡システム。
2. 前記候補値修正部は、前記候補値が表す局所領域の値の相関を用いて前記修正候補値を算出する修正候補値算出部を含み、
   前記局所領域の面積は、前記カットオフ周波数に応じて決定される、
   ことを特徴とする請求項１に記載の顕微鏡システム。
3. 前記データ補正は、前記候補値が表す局所領域の値の相関を用いて小振幅信号と判断される変動成分を抑制することを特徴とする請求項１に記載の顕微鏡システム。
4. 前記相関は、離れた２点の値の相関である第１の相関と、距離による相関である第２の相関とを含み、
   前記第１の相関と前記第２の相関とは、前記カットオフ周波数に応じて決定される、
   ことを特徴とする請求項２に記載の顕微鏡システム。
5. 前記修正候補値算出部は、前記第１の相関と前記第２の相関とを含むバイラテラルフィルタを用いて前記修正候補値を算出することを特徴とする請求項４に記載の顕微鏡システム。
6. 前記相関は、確率密度関数に基づくことを特徴とする請求項４に記載の顕微鏡システム。
7. 前記データ補正は、前記候補値に含まれる孤立変動成分を抑制することを含むことを特徴とする請求項２に記載の顕微鏡システム。
8. 前記データ補間は、距離による相関を用いて、前記候補値のうち決定できなかった値を補間することを含むことを特徴とする請求項１乃至７のうち何れか１項に記載の顕微鏡システム。
9. 前記修正候補値と前記複数の画像とに基づいて、全焦点画像を作成する全焦点画像作成部をさらに具備することを特徴とする請求項１乃至８のうち何れか１項に記載の顕微鏡システム。
10. 前記修正候補値と前記複数の画像とに基づいて、３次元再構成画像を作成する３次元再構成画像作成部をさらに具備することを特徴とする請求項１乃至８のうち何れか１項に記載の顕微鏡システム。
11. 前記光学系の情報の、対物レンズの開口数をＮＡ_ｏｂｊ 、前記観察光学系の総合倍率をβ、前記光源の所定の波長をλとすると、
    前記カットオフ周波数取得部は、前記カットオフ周波数をＦｃとしたときに、
    に基づいて、前記カットオフ周波数を取得する、
    ことを特徴とする請求項１乃至１０のうち何れか１項に記載の顕微鏡システム。
12. 前記修正候補値に最適化を施して前記標本の３次元形状を推定する３次元形状推定部をさらに具備する請求項１乃至１１のうち何れか１項に記載の顕微鏡システム。
13. 前記面積は、前記カットオフ周波数が低いほど広くなるように決定されることを特徴とする請求項２に記載の顕微鏡システム。
14. 前記距離による相関は、前記カットオフ周波数が低いほど広い領域の値の重みが高くなるように決定されることを特徴とする請求項８に記載の顕微鏡システム。