(57)【要約】【課題】焦点位置が異なる複数の撮像画像から画素位置ごとに焦点が合った画像及び高さ情報を、大域的な領域に亘って適切に焦点が合うように生成する。【解決手段】焦点位置が異なる複数の画像を取得（ステップＳ１）後、グレースケール画像の取得、多重解像度画像の生成、確率分布の生成を行う（ステップＳ２〜Ｓ１４）。確率分布に基づくコスト関数とペナルティ関数とを足し合わせた評価関数を焦点位置の割り当てＸに対してＥ（Ｘ）とするとき、別の割り当てＸ＊（ステップＳ１０５）に対してＥ（Ｘ＊）＜Ｅ（Ｘ）となるか否かを判定する（ステップＳ１０７）。これを満たす場合に、Ｘ＊をＸに代入して再びステップＳ１０５及びＳ１０７を実行する。最適な焦点位置から、画素位置ごとに焦点が合った画像及び高さ情報を生成する（ステップＳ２２）。【選択図】図１４
2014010803.tif 000002
【特許請求の範囲】
【請求項１】
焦点位置の異なる、Ｌ個（Ｌ：２以上の自然数）の被写体画像を撮像する撮像ステップと、
前記撮像ステップにおいて撮像された像に基づいて、前記焦点位置ごとに前記被写体のグレースケール画像を取得するグレースケール画像取得ステップと、
前記グレースケール画像取得ステップにおいて取得された前記グレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成ステップと、
前記多重解像度画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、
前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得ステップと、
前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、前記撮像ステップにおいて撮像された前記焦点位置ごとの像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成ステップと、
を備えており、
前記最適焦点位置取得ステップにおいて、
Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、
あるＸ＊を割り当てる割り当てステップと、Ｅ（Ｘ＊）を算出する算出ステップと、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定ステップとを実行すると共に、前記判定ステップにおいてＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定された場合に、Ｘ＊を新たなＸとして前記割り当てステップ、前記算出ステップ及び前記判定ステップを再び実行することを特徴とする全焦点画像生成方法。
ただし、
Ｎ：前記最適な焦点位置を取得するべき全画素数
ｉ，ｊ：Ｎ以下の自然数
ｘｉ，ｘｉ＊：Ｌ以下の自然数
Ｘ＝（ｘ１，ｘ２，…，ｘＮ）
Ｘ＊＝（ｘ１＊，ｘ２＊，…，ｘＮ＊）≠Ｘ
とする。
【請求項２】
前記割り当てステップが、
α及びβ（α≠β）をＬ以下の自然数からそれぞれ選択する選択ステップと、
ｘｋ∈｛α,β｝を満たすｋ番目の画素位置（ｋ：Ｎ以下の自然数）に対応するノード、αに対応するターミナルノードα、βに対応するターミナルノードβ及びこれらのノード同士を結ぶエッジからなるグラフの最小カットを導出する導出ステップと、
前記導出ステップにおいて導出された最小カットに基づいて、ターミナルノードαから切断されたノードに対応する割り当てをβとし、ターミナルノードβから切断されたノードに対応する割り当てをαとすることにより、ＸからＸ＊を取得する取得ステップとを含んでおり、
前記割り当てステップにおいて、一のα及びβの組み合わせに関して前記選択ステップ、前記導出ステップ及び前記取得ステップを実行した後、これらのステップを別の前記組み合わせに関して再び実行することを特徴とする請求項１に記載の全焦点画像生成方法。
【請求項３】
前記導出ステップにおいて、
ｉ及びｊ番目の画素位置に対応するノード同士を結ぶリンクのコストがＶ（ｘｉ，ｘｊ）で表され、画素位置に対応するノードと前記ターミナルノードαを結ぶリンクのコスト、並びに、画素位置に対応するノードと前記ターミナルノードβを結ぶリンクのコストが数式１のｔliで表されるとき、前記グラフの最小カットとして、切断されるノードのコストの総和が最小となるようなカットを導出することを特徴とする請求項２に記載の全焦点画像生成方法。
【数１】
2014010803.tif 000012
ただし、
Ｎｉ：ｉ番目の画素位置に隣接する画素位置の番号の集合
ｌ∈｛α，β｝
とする。
【請求項４】
焦点位置の異なる、Ｌ個（Ｌ：２以上の自然数）の被写体画像から、被写体の合焦点画像を生成する全焦点画像生成装置であって、
被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成手段と、
前記多重解像度画像生成手段が生成した前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出手段と、
前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得手段と、
前記最適焦点位置取得手段が取得した前記最適な焦点位置に基づいて、前記焦点位置ごとの被写体の像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成手段と、
を備えており、
前記最適焦点位置取得手段が、
Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、あるＸ＊を割り当てる割り当て手段、Ｅ（Ｘ＊）を算出する算出手段、及び、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定手段を有すると共に、前記判定手段がＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定した場合に、Ｘ＊を新たなＸとして前記割り当て手段、前記算出手段及び前記判定手段による処理を再び実行させることを特徴とする全焦点画像生成装置。
ただし、
Ｎ：前記最適な焦点位置を取得するべき全画素数
ｉ，ｊ：Ｎ以下の自然数
ｘｉ，ｘｉ＊：Ｌ以下の自然数
Ｘ＝（ｘ１，ｘ２，…，ｘＮ）
Ｘ＊＝（ｘ１＊，ｘ２＊，…，ｘＮ＊）≠Ｘ
とする。
【請求項５】
焦点位置の異なる、Ｌ個（Ｌ：２以上の自然数）の被写体画像から、被写体の合焦点画像を生成するようにコンピュータを機能させるプログラムであって、
被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成ステップと、
前記多重解像度画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、
前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得ステップと、
前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、前記撮像ステップにおいて撮像された前記焦点位置ごとの像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成ステップと、
をコンピュータに実行させると共に、
前記最適焦点位置取得ステップにおいて、
Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、
あるＸ＊を割り当てる割り当てステップと、Ｅ（Ｘ＊）を算出する算出ステップと、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定ステップとをコンピュータに実行させると共に、前記判定ステップにおいてＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定された場合に、Ｘ＊を新たなＸとして前記割り当てステップ、前記算出ステップ及び前記判定ステップを再びコンピュータに実行させることを特徴とする全焦点画像生成プログラム。
ただし、
Ｎ：前記最適な焦点位置を取得するべき全画素数
ｉ，ｊ：Ｎ以下の自然数
ｘｉ，ｘｉ＊：Ｌ以下の自然数
Ｘ＝（ｘ１，ｘ２，…，ｘＮ）
Ｘ＊＝（ｘ１＊，ｘ２＊，…，ｘＮ＊）≠Ｘ
とする。
【請求項６】
焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるＬ個（Ｌ：２以上の自然数）の被写体画像を撮像する撮像ステップと、
前記撮像ステップにおいて撮像された像に基づいて、前記焦点位置ごとに前記被写体のグレースケール画像を取得するグレースケール画像取得ステップと、
前記グレースケール画像取得ステップにおいて取得された前記グレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成ステップと、
前記多重解像度画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、
前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得ステップと、
前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得ステップと、
を備えており、
前記最適焦点位置取得ステップにおいて、
Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、
あるＸ＊を割り当てる割り当てステップと、Ｅ（Ｘ＊）を算出する算出ステップと、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定ステップとを実行すると共に、前記判定ステップにおいてＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定された場合に、Ｘ＊を新たなＸとして前記割り当てステップ、前記算出ステップ及び前記判定ステップを再び実行することを特徴とする被写体高さ情報取得方法。
ただし、
Ｎ：前記最適な焦点位置を取得するべき全画素数
ｉ，ｊ：Ｎ以下の自然数
ｘｉ，ｘｉ＊：Ｌ以下の自然数
Ｘ＝（ｘ１，ｘ２，…，ｘＮ）
Ｘ＊＝（ｘ１＊，ｘ２＊，…，ｘＮ＊）≠Ｘ
とする。
【請求項７】
焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるＬ個（Ｌ：２以上の自然数）の被写体画像に基づいて、当該被写体の位置を取得する被写体高さ情報取得装置であって、
被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成手段と、
前記多重解像度画像生成手段が生成した前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出手段と、
前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得手段と、
前記最適焦点位置取得手段が取得した前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得手段と、
を備えており、
前記最適焦点位置取得手段が、
Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、あるＸ＊を割り当てる割り当て手段、Ｅ（Ｘ＊）を算出する算出手段、及び、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定手段を有すると共に、前記判定手段がＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定した場合に、Ｘ＊を新たなＸとして前記割り当て手段、前記算出手段及び前記判定手段による処理を再び実行させることを特徴とする被写体高さ情報取得装置。
ただし、
Ｎ：前記最適な焦点位置を取得するべき全画素数
ｉ，ｊ：Ｎ以下の自然数
ｘｉ，ｘｉ＊：Ｌ以下の自然数
Ｘ＝（ｘ１，ｘ２，…，ｘＮ）
Ｘ＊＝（ｘ１＊，ｘ２＊，…，ｘＮ＊）≠Ｘ
とする。
【請求項８】
焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるＬ個（Ｌ：２以上の自然数）の被写体画像に基づいて、当該被写体の高さ情報を取得するようにコンピュータを機能させるプログラムであって、
被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成ステップと、
前記多重解像度画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、
前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得ステップと、
前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得ステップと、
をコンピュータに実行させると共に、
前記最適焦点位置取得ステップにおいて、
Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、
あるＸ＊を割り当てる割り当てステップと、Ｅ（Ｘ＊）を算出する算出ステップと、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定ステップとをコンピュータに実行させると共に、前記判定ステップにおいてＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定された場合に、Ｘ＊を新たなＸとして前記割り当てステップ、前記算出ステップ及び前記判定ステップを再びコンピュータに実行させることを特徴とする被写体高さ情報取得プログラム。
ただし、
Ｎ：前記最適な焦点位置を取得するべき全画素数
ｉ，ｊ：Ｎ以下の自然数
ｘｉ，ｘｉ＊：Ｌ以下の自然数
Ｘ＝（ｘ１，ｘ２，…，ｘＮ）
Ｘ＊＝（ｘ１＊，ｘ２＊，…，ｘＮ＊）≠Ｘ
とする。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、全焦点画像生成方法、特に、複数の焦点位置において撮像された画像に基づく全焦点画像生成方法、全焦点画像生成装置、全焦点画像生成プログラム、被写体高さ情報取得方法、被写体高さ情報取得装置及び被写体高さ情報取得プログラムに関する。
【背景技術】
【０００２】
レンズを介して被写体を撮像する際、焦点の合っている部分と合っていない部分が同じ視野内に含まれることにより、撮像された画像において特定の領域以外がぼけた画像となることが知られている。特に、高倍率のレンズほど焦点深度が浅くなり、このような現象が生じやすい。例えば、高倍率の対物レンズ２を使用した顕微鏡１で標本等の被写体を撮像すると、標本の凹凸等によって、画像中のある領域では焦点が合っている一方、別の領域では焦点が合っておらず、画像から被写体の全体像を正確に捉えることができないという事態が生じ得る。この問題を解決するため、レンズ又は被写体を光軸方向に移動させて異なる焦点位置で被写体を撮像し、これにより得られた複数の画像に基づいて、画像内の全ての領域で焦点が合った全焦点画像を生成するという処理がなされている。
【０００３】
上記処理のために、従来、さまざまな画像処理技術が提案されている。特許文献１では、異なる焦点位置において撮像された複数の画像をそれぞれウェーブレット変換し、変換した多重解像表現の各画素において、焦点位置が異なる画像間で最大の振幅値を取るウェーブレット値を抽出する。そして、抽出したウェーブレット値を各画素に配置した１つの合成画像を生成し、その合成画像に逆ウェーブレット変換を施すことで、全焦点画像を形成している。
【０００４】
特許文献２では、特許文献１と同様に、異なる焦点位置において撮像された複数の画像をウェーブレット変換する。その後、注目画素の近傍画素同士でウェーブレット係数を足し合わせた値を、焦点位置が異なる画像同士で比較し、その値が最大となる画像から画素を抽出して全焦点画像を合成している。なお、特許文献２では、ウェーブレット変換後の各サブバンドに対し、処理対象となる各画素を中心とした縦５画素、横５画素の画素範囲にモルフォロジー処理を施す。これにより、ノイズの低減等を実現し、高周波領域及び低周波領域の境界が強調された全焦点画像を生成している。
【０００５】
特許文献３では、異なる焦点位置において撮像された複数の画像のそれぞれに関し、エッジ部分、エッジ部分の近傍、及び、エッジ部分でもその近傍でもない部分の３つの領域に分け、領域ごとに異なる方法で合焦位置を取得している。まず、注目画素とその周囲の画素との輝度値の変化量を、焦点位置が異なる画像同士で比較し、最大値を取得する。そして、取得した最大値が基準となる変化量より大きい場合に、その注目画素をエッジ部分とすると共に、その最大値を取る画像の位置を合焦位置とする。一方、エッジ部分の近傍に位置する画素に関しては、深度方向に関する輝度値の変曲点を合焦位置とする。さらに、エッジ部分でもエッジ部分の近傍でもない画素に関しては、エッジ部分とエッジ部分の近傍とにおける合焦位置の画素として取得した画素からの補間により、焦点が合った画素を導出している。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開平６−２４３２５０号公報
【特許文献２】特開２０１０−１２９０７７号公報
【特許文献３】特開２０１０−１６６２４７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００７】
特許文献１〜３のいずれも、上記のとおり、ウェーブレット値を利用したり、輝度値の変化量が大きいものを取り出したりして合焦位置の画素を抽出している。しかしながら、例えば、撮像時に画像ノイズ、飽和、ぼけなどの問題が生じた場合、撮像画像において、焦点が合っている領域ではコントラスト値が小さくなり、焦点が合っていない領域ではコントラスト値が大きくなることがある。このような問題が生じると、特許文献１や特許文献２を適用し、ウェーブレット値や輝度値の変化量の大小関係を比較して合焦位置を取得する場合、焦点が合っていない位置を合焦位置と誤って判定してしまうおそれがある。
【０００８】
また、特許文献３は、画像を３つの領域に分け、領域ごとに異なる方法で焦点位置を取得することにより、上記の問題の解決を図っている。しかしながら、特許文献３によると、エッジ部分やその近傍から離隔した部分では、エッジ部分の画素とエッジ部分の近傍の画素とに基づく補間により合焦画素を取得する。したがって、補間対象となる画素がエッジ部分やその近傍から遠く離れていたり、これらの部分と深度方向に異なる位置にあったりすると、補間の際に生じる誤差が大きくなる。このため、仮にエッジ部分とその近傍からなる局所的な領域においては適切に焦点が合った画像を生成できたとしても、画像中の大域的な領域に亘って適切に焦点が合った画像を生成できないおそれがある。
【０００９】
本発明の第１の目的は、焦点位置が異なる複数の撮像画像から画素位置ごとに焦点が合った画像を生成する方法であって、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って適切に焦点が合った画像を生成できる全焦点画像生成方法、並びに、当該方法に対応した全焦点画像生成装置及び全焦点画像生成プログラムを提供することにある。
【００１０】
また、焦点位置が異なる複数の撮像画像から画素位置ごとに焦点が合っているものを抽出することは、焦点位置から換算される高さ情報を画素位置ごとに取得することと等しい。そして、上記のような画像ノイズ等による問題は、被写体の高さ情報を画素位置ごとに取得する場合にも全く同様の問題として表れる。
【００１１】
本発明の第２の目的は、焦点位置が異なる複数の撮像画像に基づいて被写体の高さ情報を画素位置ごとに取得する方法であって、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って適切に被写体の高さ情報を取得できる被写体高さ情報取得方法、並びに、当該方法に対応した被写体高さ情報取得装置及び被写体高さ情報取得プログラムを提供することにある。
【課題を解決するための手段】
【００１２】
第１の目的を達成するためになされた本発明の第１の観点によると、本発明の全焦点画像生成方法は、焦点位置の異なる、Ｌ個（Ｌ：２以上の自然数）の被写体画像を撮像する撮像ステップと、前記撮像ステップにおいて撮像された像に基づいて、前記焦点位置ごとに前記被写体のグレースケール画像を取得するグレースケール画像取得ステップと、前記グレースケール画像取得ステップにおいて取得された前記グレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成ステップと、前記多重解像度画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、前記撮像ステップにおいて撮像された前記焦点位置ごとの像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成ステップと、を備えており、前記最適焦点位置取得ステップにおいて、Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、あるＸ＊を割り当てる割り当てステップと、Ｅ（Ｘ＊）を算出する算出ステップと、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定ステップとを実行すると共に、前記判定ステップにおいてＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定された場合に、Ｘ＊を新たなＸとして前記割り当てステップ、前記算出ステップ及び前記判定ステップを再び実行する。ただし、Ｎ：前記最適な焦点位置を取得するべき全画素数、ｉ，ｊ：Ｎ以下の自然数、ｘｉ，ｘｉ＊：Ｌ以下の自然数、Ｘ＝（ｘ１，ｘ２，…，ｘＮ）、Ｘ＊＝（ｘ１＊，ｘ２＊，…，ｘＮ＊）≠Ｘとする。
【００１３】
本発明の全焦点画像生成装置は、焦点位置の異なる、Ｌ個（Ｌ：２以上の自然数）の被写体画像から、被写体の合焦点画像を生成する全焦点画像生成装置であって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成手段と、前記多重解像度画像生成手段が生成した前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出手段と、前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得手段と、前記最適焦点位置取得手段が取得した前記最適な焦点位置に基づいて、前記焦点位置ごとの被写体の像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成手段と、を備えており、前記最適焦点位置取得手段が、Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、あるＸ＊を割り当てる割り当て手段、Ｅ（Ｘ＊）を算出する算出手段、及び、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定手段を有すると共に、前記判定手段がＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定した場合に、Ｘ＊を新たなＸとして前記割り当て手段、前記算出手段及び前記判定手段による処理を再び実行させる。
【００１４】
本発明の全焦点画像生成プログラムは、焦点位置の異なる、Ｌ個（Ｌ：２以上の自然数）の被写体画像から、被写体の合焦点画像を生成するようにコンピュータを機能させるプログラムであって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成ステップと、前記多重解像度画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、前記撮像ステップにおいて撮像された前記焦点位置ごとの像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成ステップと、をコンピュータに実行させると共に、前記最適焦点位置取得ステップにおいて、Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、あるＸ＊を割り当てる割り当てステップと、Ｅ（Ｘ＊）を算出する算出ステップと、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定ステップとをコンピュータに実行させると共に、前記判定ステップにおいてＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定された場合に、Ｘ＊を新たなＸとして前記割り当てステップ、前記算出ステップ及び前記判定ステップを再びコンピュータに実行させる。
【００１５】
第１の観点に係る本発明の全焦点画像生成方法、全焦点画像生成装置及びプログラムによると、Ｘ＊を割り当て、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定すると共に、これを満たす場合には、Ｘ＊を新たなＸとして再び上記の割り当て及び判定を実行する。これにより、現状のＸより評価関数が小さくなるようなＸ＊が見つかれば、これをあらためてＸとした上で上記の処理を繰り返すので、評価関数が最小（極小）となる最適な焦点位置に到達できる。評価関数の大きさは、隣接する画素位置同士の焦点位置の割り当て方に依存するため、ひいては、画素位置同士を繋いでいくことで構成される広範囲の割り当て方に依存する。このため、評価関数が最小（極小）となるような焦点位置を取得することにより、大域的な領域に亘って適切に焦点が合った画像を取得できる。なお、Ｅ（Ｘ＊）＜Ｅ（Ｘ）となることがある程度期待できる、Ｘ＊を割り当てる方法としては、グラフカット、平均場近似を利用する方法など、さまざまなものが考えられる。また、最初のＸ＊を割り当てる際には確率伝播法を用い、次のＸ＊を割り当てる段階からグラフカットなどその他の方法を用いてもよい。
【００１６】
また、本発明においては、前記割り当てステップが、α及びβ（α≠β）をＬ以下の自然数からそれぞれ選択する選択ステップと、ｘｋ∈｛α,β｝を満たすｋ番目の画素位置（ｋ：Ｎ以下の自然数）に対応するノード、αに対応するターミナルノードα、βに対応するターミナルノードβ及びこれらのノード同士を結ぶエッジからなるグラフの最小カットを導出する導出ステップと、前記導出ステップにおいて導出された最小カットに基づいて、ターミナルノードαから切断されたノードに対応する割り当てをβとし、ターミナルノードβから切断されたノードに対応する割り当てをαとすることにより、ＸからＸ＊を取得する取得ステップとを含んでおり、前記割り当てステップにおいて、一のα及びβの組み合わせに関して前記選択ステップ、前記導出ステップ及び前記取得ステップを実行した後、これらのステップを別の前記組み合わせに関して再び実行することが好ましい。これによると、グラフの最小カットに対応するＸ＊を求めることにより、評価関数がより小さくなることが期待できる。したがって、評価関数が最小となるような割り当てに到達しやすくなる。また、種々のα及びβの組み合わせに関して割り当てＸ＊を求めることで、評価関数がより小さくなるような割り当てを導出できる。
【００１７】
また、本発明においては、前記導出ステップにおいて、画素位置に対応するノード同士を結ぶリンクのコストがＶ（ｘｉ，ｘｊ）で表され、画素位置に対応するノードと前記ターミナルノードαを結ぶリンクのコスト、並びに、画素位置に対応するノードと前記ターミナルノードβを結ぶリンクのコストが数式１のｔliで表されるとき、前記グラフの最小カットとして、切断されるノードのコストの総和が最小となるようなカットを導出することが好ましい。これによると、最小カットを導出することで、評価関数が小さくなるようなＸ＊を割り当てることができる。ただし、Ｎｉ：ｉ番目の画素位置に隣接する画素位置の番号の集合、ｌ∈｛α，β｝とする。
【数１】
2014010803.tif 000003
【００１８】
第２の目的を達成するためになされた本発明の第２の観点によると、本発明の被写体高さ情報取得方法は、焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるＬ個（Ｌ：２以上の自然数）の被写体画像を撮像する撮像ステップと、前記撮像ステップにおいて撮像された像に基づいて、前記焦点位置ごとに前記被写体のグレースケール画像を取得するグレースケール画像取得ステップと、前記グレースケール画像取得ステップにおいて取得された前記グレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成ステップと、前記多重解像度画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得ステップと、を備えており、前記最適焦点位置取得ステップにおいて、Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、あるＸ＊を割り当てる割り当てステップと、Ｅ（Ｘ＊）を算出する算出ステップと、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定ステップとを実行すると共に、前記判定ステップにおいてＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定された場合に、Ｘ＊を新たなＸとして前記割り当てステップ、前記算出ステップ及び前記判定ステップを再び実行する。ただし、Ｎ：前記最適な焦点位置を取得するべき全画素数、ｉ，ｊ：Ｎ以下の自然数、ｘｉ，ｘｉ＊：Ｌ以下の自然数、Ｘ＝（ｘ１，ｘ２，…，ｘＮ）、Ｘ＊＝（ｘ１＊，ｘ２＊，…，ｘＮ＊）≠Ｘとする。
【００１９】
本発明の被写体高さ情報取得装置は、焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるＬ個（Ｌ：２以上の自然数）の被写体画像に基づいて、当該被写体の位置を取得する装置であって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成手段と、前記多重解像度画像生成手段が生成した前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出手段と、前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得手段と、前記最適焦点位置取得手段が取得した前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得手段と、を備えており、前記最適焦点位置取得手段が、Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、あるＸ＊を割り当てる割り当て手段、Ｅ（Ｘ＊）を算出する算出手段、及び、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定手段を有すると共に、前記判定手段がＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定した場合に、Ｘ＊を新たなＸとして前記割り当て手段、前記算出手段及び前記判定手段による処理を再び実行させる。
【００２０】
本発明の被写体高さ情報取得プログラムは、焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるＬ個（Ｌ：２以上の自然数）の被写体画像に基づいて、当該被写体の高さ情報を取得するようにコンピュータを機能させるプログラムであって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成ステップと、前記多重解像度画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得ステップと、をコンピュータに実行させると共に、前記最適焦点位置取得ステップにおいて、Ｎ個の画素位置に対して前記焦点位置としてＸを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づくｉ番目の画素位置のコスト関数Ｄｉ（ｘｉ）と、ｘｉとｘｊの差の絶対値が大きいほど大きい値を取るペナルティ関数Ｖ（ｘｉ，ｘｊ）とを、前記Ｎ個の画素位置並びに隣り合うｉ及びｊ番目の画素位置に関して足し合わせた関数を評価関数Ｅ（Ｘ）とするとき、あるＸ＊を割り当てる割り当てステップと、Ｅ（Ｘ＊）を算出する算出ステップと、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定する判定ステップとをコンピュータに実行させると共に、前記判定ステップにおいてＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすと判定された場合に、Ｘ＊を新たなＸとして前記割り当てステップ、前記算出ステップ及び前記判定ステップを再びコンピュータに実行させる。
【００２１】
第２の観点に係る本発明の被写体高さ情報取得方法、被写体高さ情報取得装置及びプログラムによると、第１の観点と同様に、評価関数が最小（極小）となる最適な焦点位置を取得することができる。そして、評価関数の大きさが最小（極小）となるような焦点位置を取得することにより、画像の大域的な領域に亘って適切に被写体の高さ情報を取得できる。
【発明の効果】
【００２２】
第１の観点に係る本発明の全焦点画像生成方法、全焦点画像生成装置及び全焦点画像生成プログラムによると、Ｘ＊を割り当て、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定すると共に、これを満たす場合には、Ｘ＊を新たなＸとして再び上記の割り当て及び判定を実行する。これにより、現状のＸより評価関数が小さくなるようなＸ＊が見つかれば、これをあらためてＸとした上で上記の処理を繰り返すので、評価関数が最小（極小）となる最適な焦点位置に到達できる。評価関数の大きさは、隣接する画素位置同士の焦点位置の割り当て方に依存するため、ひいては、画素位置同士を繋いでいくことで構成される広範囲の割り当て方に依存する。このため、評価関数が最小（極小）となるような焦点位置を取得することにより、大域的な領域に亘って適切に焦点が合った画像を取得できる。
【００２３】
第２の観点に係る本発明の被写体高さ情報取得方法、被写体高さ情報取得装置及びプログラムによると、第１の観点と同様に、評価関数が最小（極小）となる最適な焦点位置を取得することができる。そして、評価関数の大きさが最小（極小）となるような焦点位置を取得することにより、画像の大域的な領域に亘って適切に被写体の高さ情報を取得できる。
【図面の簡単な説明】
【００２４】
【図１】第１の参考形態に係る全焦点画像生成装置の構成を示すブロック図である。
【図２】サンプル、フォーカスレンズ、対物レンズ及び撮像素子の位置関係を光の経路に沿った方向に関して示す模式図である。
【図３】画像処理部が施す変換処理による変換前後の画像を示す模式図である。
【図４】図４（ａ）は、１回のウェーブレット変換により生成される解像度分解画像に係る多解像度表現を示す。図４（ｂ）は、３回のウェーブレット変換により生成される解像度分解画像に係る多解像度表現を示す。
【図５】ウェーブレット変換の変換前後の画像における画素の対応を示す。
【図６】図６（ａ）は、Ｌ個の多重解像度表現を焦点位置ごとに並べたものである。図６（ｂ）は、候補値として抽出されたウェーブレット係数と、その候補値に対応する多重解像度表現との一例を示している。
【図７】候補値に基づいて確率値に投票する方法を示す模式図である。
【図８】図８（ａ）は、確率分布の一例を示すグラフである。図８（ｂ）は、確率分布の別の一例を示すグラフである。
【図９】縦３×横３の９個の画素配列を示す模式図である。白丸が各画素を示している。
【図１０】高解像度での画素配置と低解像度の１画素との関係を示す模式図である。
【図１１】顕微鏡及び画像処理部が実行する処理の全体の流れを示すフローチャートである。
【図１２】第２の参考形態において第１の参考形態と異なる部分を示すブロック図である。
【図１３】本発明の一実施の形態に係る全焦点画像生成装置において、第１の参考形態と異なる部分を示すブロック図である。
【図１４】本発明の一実施の形態に係る全焦点画像生成装置が実行する処理の全体の流れを示すフローチャートである。
【図１５】図１４の処理において使用されるグラフの一例である。
【図１６】図１５のグラフに係る最小カットの一例である。
【発明を実施するための形態】
【００２５】
（第１の参考形態）
本発明に強く関連する第１及び第２の参考形態について、図面に基づいて説明する。第１及び第２の参考形態において、後述する本発明に係る一実施の形態と異なるのは、最適焦点位置取得部１３０の部分のみである。このため、以下の第１及び第２の参考形態に係る説明のうち、最適焦点位置取得部１３０以外の部分に関して、本発明に係る一実施の形態の構成に関する説明として参照されたい。第１の参考形態の全焦点画像生成装置は、図１に示すように、サンプルＳの拡大画像を撮像する顕微鏡１と顕微鏡１によって撮像された画像を処理する画像処理部１００とを有している。顕微鏡１は、複数の対物レンズ２が固定されたレボルバ３と、サンプルＳが設置されるステージ４と、ステージ４が上端に固定されたスライダ５と、スライダ５を鉛直方向に移動可能に支持するスライダ支持体７と、スライダ５をスライダ支持体７に対して鉛直方向に移動させるスライダ移動機構６とを有している。
【００２６】
レボルバ３には、互いに異なる倍率の対物レンズ２が固定されている。レボルバ３を回転させることで、下方のステージ４側に向いた対物レンズ２を切り替えることができる。ステージ４の上面は水平であり、対物レンズ２の直下に水平に位置するようにサンプルＳをステージ４上に置くことができる。スライダ移動機構６は、スライダ５を鉛直方向に微小に移動させる圧電アクチュエータと、鉛直方向に関するスライダ５の位置を高精度に検出する位置検出センサとを有している。位置検出センサは、スライダ５の位置を検出した結果を後述の制御部１１へと出力する。
【００２７】
また、顕微鏡１は、サンプルＳへの照明光を射出する照明光学系８と、サンプルＳからの光を撮像素子１０に結像させる撮像光学系９と、結像した被写体像を光電変換して出力する撮像素子１０とを有している。照明光学系８は、照明用の光源や光源からの照明光をサンプルＳへと導く反射鏡等を有している。撮像光学系９は、対物レンズ２からの被写体像を撮像素子１０へと導き、結像させるフォーカスレンズ９ａ（図２参照）などのレンズや反射鏡等を有している。被写体像は、対物レンズ２と撮像光学系９に設けられたレンズとの協働により、高倍率で拡大されて撮像素子１０に結像される。撮像素子１０は、撮像した被写体画像に対応する画像データを生成し、画像処理部１００へと出力する。この画像データは、ｘ方向、及び、ｘ方向に直交するｙ方向に並んだ複数のカラー画素を示すデータからなる。水平方向に平行な所定の方向がｘ方向に対応し、水平方向に平行で上記所定の方向に直交する方向がｙ方向に対応する。この画像データにおいては、各画素の色が、光の三原色であるＲ（赤）、Ｇ（緑）及びＢ（青）のうち、一部の成分のみによって表されている。あるいは、各画素の色が、補色であるＣｙ（シアン），Ｍｇ（マゼンタ）、Ｙｅ（イエロー）及びＧのうち、一部の成分のみによって表されていてもよい。
【００２８】
顕微鏡１には、各部の構成を制御する制御部１１が設けられている。制御部１１は、スライダ移動機構６の駆動を制御し、ステージ４を鉛直方向に関してさまざまな位置に移動させる。ステージ４が移動すると、ステージ４上のサンプルＳと対物レンズ２との距離が変動するため、ステージ４上に配置されたサンプルＳからフォーカスレンズ９ａまでの光の経路に沿った距離が変動する。一方、フォーカスレンズ９ａの焦点距離が一定であり、フォーカスレンズ９ａと撮像素子１０との間の距離も一定であるとすると、被写体像が撮像素子１０上にちょうど焦点が合うように結像するためには、被写体からフォーカスレンズ９ａまでの光の経路に沿った距離が所定の大きさになる必要がある。以下、被写体において、フォーカスレンズ９ａに対してちょうど焦点が合った位置を「焦点位置」と呼称する。
【００２９】
サンプルＳがフォーカスレンズ９ａに対して焦点位置にあるとき、サンプルＳの像は焦点が合った状態で撮像素子１０に結像する。しかしながら、顕微鏡１で用いられる高倍率のレンズは、焦点深度が浅い。このため、サンプルＳの表面に凹凸がある場合には、フォーカスレンズ９ａまでの距離がサンプルＳの表面の位置ごとに異なることになる。この距離の違いが焦点深度に対して無視できない大きさになると、サンプルＳの表面上の特定の領域が焦点位置に配置されたとしても、特定の領域以外の領域は焦点位置からずれてしまい、画像がぼけてしまうことが生じ得る。例えば、図２に示すように、サンプルＳの表面における領域ｓ１と領域ｓ２とでは、光の経路に沿った方向に関してフォーカスレンズ９ａまでの距離Ｄ１及びＤ２が互いに異なる。したがって、この状態で領域ｓ１は焦点が合っていたとしても、領域ｓ２では焦点が合わないことがあり得る。被写体のあらゆる領域で焦点が合った画像を取得するためには、焦点位置が少しずつ異なった複数の画像から、焦点が合った領域の画素を集め、合成することが考えられる。
【００３０】
そこで、制御部１１は、スライダ移動機構６を制御し、焦点位置が少しずつ異なるようにサンプルＳを配置させると共に、撮像素子１０を制御し、各焦点位置でのサンプルＳの被写体像を撮像させる。撮像された画像同士は、焦点位置が少しずつ異なるため、被写体において焦点が合っている領域が少しずつ異なったものとなる。制御部１１は、撮像素子１０が、焦点位置の異なる画像データを１つ、画像処理部１００に出力するごとに、その焦点位置に対応するステージ４の位置を、スライダ移動機構６に設けられた位置検出センサの検出結果から取得し、画像処理部１００へと出力する。これらの処理が繰り返されることにより、撮像素子１０からは、焦点位置が異なるＬ個（Ｌ：２以上の自然数）の画像が出力されると共に、制御部１１からは、Ｌ個の画像に関連付けられたＬ個のステージ位置情報が出力される。
【００３１】
画像処理部１００は、撮像素子１０から入力されたＬ個の被写体画像に対応する画像データに基づいて、各画素に関して焦点が合った全焦点画像に対応する画像データを生成する。この処理のため、画像処理部１００は、前段処理部１１０、グレースケール変換部１０１、多重解像度変換部１０２（多重解像度画像生成手段）、焦点位置確率算出部１２０（焦点位置確率算出手段）、信頼度設定部１０３、最適焦点位置取得部１３０及び後段処理部１４０の各機能部を有している。
【００３２】
画像処理部１００は、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）その他の各種のハードウェアと、ＲＯＭに格納されたプログラムその他のソフトウェアとから構成されている。このソフトウェアがハードウェアを、前段処理部１１０等の各機能部を実現するように機能させることにより、画像処理部１００が実現されている。あるいは、各機能部における演算処理に特化した電子回路等により前段処理部１１０等が構成されていてもよいし、各機能部における演算処理に特化した電子回路等と、ハードウェアとソフトウェアの組み合わせによって各機能部を実現する構成とが互いに協働するように構成されていてもよい。画像処理部１００に用いられるソフトウェアは、磁気的記録媒体や光学的記録媒体、ＵＳＢメモリなどの各種記録媒体に格納されてもよい。そして、これらの記録媒体から読み出されたソフトウェアに従ってコンピュータが機能することにより、画像処理部１００が実現されてもよい。
【００３３】
画像処理部１００が有する各機能部は、画像データを処理することにより、その画像データに対応する画像に対して各種の画像処理を施す。このように、各機能部が直接処理するのは画像データであるが、以下においては、冗長な表現を避けるため、各機能部が画像自身に対して直接、各種の画像処理を施すように表現する場合がある。このような表現は、その画像処理に対応する処理をその画像に対応する画像データに対して施すことを示すものとする。
【００３４】
前段処理部１１０は、位置ずれ補正部１１１、撮像倍率補正部１１２及びＲＧＢ画像取得部１１３を有している。位置ずれ補正部１１１は、撮像素子１０からのＬ個の被写体画像同士を比較し、水平方向に関する被写体の位置ずれを検出する。この位置ずれは、画像中のｘ方向及びｙ方向に関する位置ずれとなって表れる（図３の画像ＩＭ１参照）。位置ずれ補正部１１１は、被写体画像中の被写体の位置をずらすことにより、焦点位置が異なる被写体画像同士で被写体の位置が一致するように被写体画像を補正する。
【００３５】
撮像倍率補正部１１２は、制御部１１から入力されるステージ４の位置情報に基づいて、被写体画像の倍率を補正する。撮像光学系９の構成によっては、被写体の焦点位置が変動すると、撮像倍率が変動することがある。撮像倍率補正部１１２は、ステージ４の位置に関連付けて、焦点位置の変動に伴う撮像倍率の変動分に対応する拡大縮小率を保持している。これにより、撮像倍率補正部１１２は、制御部１１から入力されるステージ４の位置に対応する拡大縮小率を導出する。そして、撮像倍率補正部１１２は、導出した拡大縮小率に従って、被写体画像を拡大又は縮小する。撮像倍率補正部１１２は、これにより、焦点位置が異なるＬ個の被写体画像同士で倍率が一致するように各被写体画像の倍率を補正する。
【００３６】
ＲＧＢ画像取得部１１３は、図３に示すように、位置ずれ補正部１１１及び撮像倍率補正部１１２による補正が施されたＬ個の被写体画像ＩＭ１に対応する画像データから、各画素においてＲ，Ｇ及びＢの各色の成分値が全て決定されたＬ個の画像ＩＭ２に対応するＲＧＢ画像データＩＭ２を生成する。画像ＩＭ１に対応する画像データには、上述のとおり、画素ごとにＲ，Ｇ及びＢのいずれか一部の値のみが含まれている。したがって、ＲＧＢ画像取得部１１３は、各画素において、足りない成分の値を、その画素の周囲に位置する画素の値に基づいて補間する。例えば、ＲＧＢ画像取得部１１３は、補間の対象となる画素の周囲の画素値に重みを乗算したものを足し合わせることにより、補間した画素値を算出する。これにより、ＲＧＢ画像取得部１１３は、焦点位置が異なるＬ個の被写体画像のそれぞれに関して、画像ＩＭ２に対応するＲＧＢ画像データを生成する。
【００３７】
グレースケール変換部１０１は、図３に示すように、ＲＧＢ画像取得部１１３が生成したＬ個の画像ＩＭ２に対応するＲＧＢ画像データを、Ｌ個のグレースケール画像ＩＭ３に対応するグレースケール画像データに変換する。変換方法はさまざまなものを採用できる。一例として、ＮＴＳＣ（National Television System Committee）方式によるＲＧＢ信号から輝度信号への変換方法を採用する場合、グレースケール変換部１０１は、以下の式を使用して、グレースケール画像の各画素値を算出する：Ｙ＝０．２９８９１２＊Ｒ＋０．５８６６１１＊Ｇ＋０．１１４４７８＊Ｂ。なお、Ｙは、グレースケール画像の各画素値であり、Ｒ，Ｇ及びＢは、ＲＧＢ画像の各画素におけるＲの成分値，Ｇの成分値及びＢの成分値である。
【００３８】
多重解像度変換部１０２は、グレースケール変換部１０１が生成したグレースケール画像データに対応する画像に多重解像度変換を施す。この多重解像度変換には、ガウシアンピラミッド変換、ラプラシアンピラミッド変換などのさまざまな方法を採用できる。本参考形態では、離散ウェーブレット変換を採用する。離散ウェーブレット変換としては、Ｈａａｒ  Ｗａｂｅｌｅｔ、Ｄａｕｂｅｃｈｉｅｓ  Ｗａｖｅｌｅｔ，Ｃｏｈｅｎ−Ｄａｕｂｅｃｈｉｅｓ−Ｆｅａｕｖｅａｕ  Ｗａｖｅｌｅｔなどのうち、いずれかを用いた変換方法を使用する。
【００３９】
離散ウェーブレット変換は、画像を高周波成分と低周波成分とに分解すると共に、分解した低周波成分をさらに高周波成分と低周波成分とに分解するという処理を繰り返すことで行われる。これにより、高周波成分及び低周波成分に分解された画像である複数の解像度分解画像が取得される。高周波成分及び低周波成分への分解には、各種のフィルタ処理が用いられる。このフィルタ処理に用いるフィルタは、画像のｘ方向に関するハイパスフィルタ及びローパスフィルタと、画像のｙ方向に関するハイパスフィルタ及びローパスフィルタとを含んでいる。以下、これらの各方向に関するフィルタを、ｘハイパスフィルタ、ｘローパスフィルタ、ｙハイパスフィルタ及びｙローパスフィルタと呼称する。
【００４０】
フィルタ処理を画像に施す具体的な方法は以下のとおりである。多重解像度変換部１０２は、複数のフィルタ行列を示すデータをフィルタごとに有している。例えば、多重解像度変換部１０２は、ｘハイパスフィルタに対応するフィルタ行列を示すデータや、ｙローパスフィルタに対応するフィルタ行列を示すデータを有している。フィルタ行列は複数のフィルタ係数から構成されており、各フィルタ係数は、フィルタが施される領域に含まれる画素ごとに設定されている。多重解像度変換部１０２は、フィルタ行列に含まれる各フィルタ係数を、その係数に対応する位置の画素値に乗算する。そして、多重解像度変換部１０２は、乗算した値を足し合わせ、その足し合わせた和を、当該フィルタ処理後の画素値とする。また、多重解像度変換部１０２は、フィルタ処理後の画像全体における画素数が、ｘ方向及びｙ方向のそれぞれに関して、フィルタ処理前の画像全体における画素数の半分になるように各フィルタ処理を施す。
【００４１】
多重解像度変換部１０２は、グレースケール画像変換部が生成した画像データに対して各種のフィルタ処理を施すことで、以下のように１回目のウェーブレット変換を施す。まず、多重解像度変換部１０２は、ｘハイパスフィルタを用いてフィルタ処理を施した後に、ｙローパスフィルタを用いてフィルタ処理を施す。これにより、グレースケール画像において、ｘ方向に関しては高周波であり、ｙ方向に関しては低周波である成分が取得される。以下、この成分をＨＬ成分と呼称する。
【００４２】
同様に、多重解像度変換部１０２は、ｘハイパスフィルタを用いてフィルタ処理を施した後に、ｙハイパスフィルタを用いてフィルタ処理を施す。これにより、ｘ方向及びｙ方向のいずれに関しても高周波である成分が取得される。以下、この成分をＨＨ成分と呼称する。また、多重解像度変換部１０２は、ｘローパスフィルタを用いてフィルタ処理を施した後に、ｙハイパスフィルタを用いてフィルタ処理を施す。これにより、ｘ方向に関しては低周波であり、ｙ方向に関しては高周波である成分が取得される。以下、この成分をＬＨ成分と呼称する。さらに、多重解像度変換部１０２は、ｘローパスフィルタを用いてフィルタ処理を施した後に、ｙローパスフィルタを用いてフィルタ処理を施す。これにより、ｘ方向及びｙ方向のいずれに関しても低周波である成分が取得される。以下、この成分をＬＬ成分と呼称する。このように、ウェーブレット変換によって、グレースケール画像から、ＨＨ成分、ＨＬ成分、ＬＨ成分及びＬＬ成分からなる、複数の周波数帯域の成分を有する解像度分解画像が生成される。解像度分解画像の各成分に含まれる、ウェーブレット変換によって導出された１つ１つの値を、ウェーブレット係数（ウェーブレット変換によって生成された変換係数）という。
【００４３】
図４（ａ）は、１回目のウェーブレット変換により生成されたＨＨ成分、ＨＬ成分、ＬＨ成分及びＬＬ成分を、ｘ方向及びｙ方向に対応するｖ方向及びｗ方向に、２行２列に並べたものである。これらの成分は、ＨＨ成分が右下隅に、ＨＬ成分が右上隅に、ＬＨ成分が左下隅に、ＬＬ成分が左上隅になるように配列されている。
【００４４】
多重解像度変換部１０２は、１回目のウェーブレット変換によって生成されたＬＬ成分に対し、２回目のウェーブレット変換を同様に施す。これにより、再びＨＨ成分、ＨＬ成分、ＬＨ成分及びＬＬ成分が得られる。このように、前回のウェーブレット変換によって得られたＬＬ成分に対してさらにウェーブレット変換を施すことの繰り返しにより、ＨＨ成分、ＨＬ成分、ＬＨ成分及びＬＬ成分の各成分がウェーブレット変換を施すたびに得られる。以下、ウェーブレット変換（ｎ：自然数）をｎ回施すことをｎレベルと呼称する。これらの各画像を、ウェーブレット変換ごとに、前回のウェーブレット変換におけるＬＬ成分の位置に図４（ａ）に示す配列で配置したものが、図４（ｂ）の画像ＩＭ４である。解像度分解画像における各成分をこのように配置したものを多重解像度表現という。画像ＩＭ４は、一例として、３レベルの多重解像度表現を示している。
【００４５】
この多重解像度表現は、図４（ｂ）に示すように、各ウェーブレット変換により得られるＨＨ成分、ＨＬ成分及びＬＨ成分と、最後のウェーブレット変換により得られるＬＬ成分とからなる。各レベルのＨＨ成分、ＨＬ成分及びＬＨ成分を高周波成分といい、最終的に生成されたＬＬ成分を低周波成分という。図４（ｂ）のＬ１は１レベルの、Ｌ２は２レベルの、Ｌ３は３レベルの高周波成分を示している。多重解像度表現は、ｖ方向に関する幅がグレースケール画像のｘ方向に関する幅と等しく、ｗ方向に関する幅がグレースケール画像のｙ方向に関する幅と等しい。グレースケール画像の各画素値は、多重解像度表現において周波数帯域が異なる複数のウェーブレット係数へと分解される。このため、グレースケール画像における画素位置（ｘ，ｙ）と、多重解像度表現における位置（ｖ，ｗ）とは、一対一には対応せず、一体多の対応関係となる。例えば、図５に示すように、グレースケール画像ＩＭ３における画素Ｐ１は、多重解像度表現ＩＭ４における異なる周波数帯域のそれぞれにおいてＰ１に対応する位置であるＱ１１〜Ｑ３４のウェーブレット係数と対応することとなる。
【００４６】
このようにして、多重解像度変換部１０２は、焦点位置が異なるＬ個のグレースケール画像データのそれぞれから多重解像度表現に対応するデータを生成する。こうして得られたＬ個の多重解像度表現を焦点位置ごとに並べたものが、図６（ａ）の多重解像度表現ＩＭ４＿１〜ＩＭ４＿Ｌである。多重解像度変換部１０２は、ｖ方向及びｗ方向に関する位置ごとに、多重解像度表現ＩＭ４＿１〜ＩＭ４＿Ｌの全てのウェーブレット係数を順に算出する。つまり、ｖ方向及びｗ方向に関するある位置について、多重解像度表現ＩＭ４＿１〜ＩＭ４＿Ｌの全てのウェーブレット係数を順に算出した後に、ｖ方向及びｗ方向に関する次の位置について、多重解像度表現ＩＭ４＿１〜ＩＭ４＿Ｌの全てのウェーブレット係数を順に算出する。
【００４７】
例えば、図６（ａ）の二点鎖線Ｒ２で囲まれた各点は、ｖ方向及びｗ方向に関する位置が同じである。また、二点鎖線Ｒ３で囲まれた各点は、ｖ方向及びｗ方向に関する位置が同じであり、Ｒ２で囲まれた点とは位置が異なる。このとき、多重解像度変換部１０２は、例えば、Ｒ２で囲まれた点の位置に関して多重解像度表現ＩＭ４＿１〜ＩＭ４＿Ｌの全てのウェーブレット係数を算出してから、Ｒ３で囲まれた点の位置に関して多重解像度表現ＩＭ４＿１〜ＩＭ４＿Ｌの全てのウェーブレット係数を算出する。多重解像度変換部１０２がこのような順序でウェーブレット係数を算出する理由については、後述する。
【００４８】
焦点位置確率算出部１２０は、多重解像度変換部１０２が生成した解像度分解画像から、グレースケール画像における画素位置ごとに、焦点位置に関する確率分布（焦点位置確率分布）を算出する。ここで、焦点位置に関する確率値とは、次のようなものとする。以下、ｉ（ｉ：グレースケール画像の総画素数以下の自然数）をグレースケール画像の画素位置を区別する添え字とし、ｘｉをＬ以下の自然数を取る変数とする。このとき、焦点位置確率算出部１２０が算出する確率分布をＰｉ（ｘｉ）と表す。Ｐｉ（ｘｉ）は、画素位置ｉにおいて、焦点が合っているものが、Ｌ個の画像のうち、ｘｉ番目の画像である確からしさを示す。言い換えれば、Ｐｉ（ｘｉ）は、画素位置ｉにおいて被写体の表面が焦点位置に配置されているものが、Ｌ個の画像のうち、ｘｉ番目の画像であることの確からしさに相当する。
【００４９】
このような確率分布を導出するため、焦点位置確率算出部１２０は、画像候補決定部１２１（抽出手段）及び確率分布算出部１２２（確率算出手段）を有している。画像候補決定部１２１は、多重解像度表現におけるｖ方向及びｗ方向に関する位置ごとに、Ｌ個の多重解像度表現から、確率分布を算出するための候補となるＭ個（Ｍ：Ｌ未満の自然数）の多重解像度表現を抽出する。候補の数は、パラメータチューニングによって適宜設定される。このＭ個の多重解像度表現は、ｖ方向及びｗ方向に関する当該位置に関して、Ｌ個の画像から、ウェーブレット係数の絶対値が大きい上位Ｍ個を抽出することにより取得される。具体的には、画像候補決定部１２１は、当該位置におけるウェーブレット係数の絶対値を、焦点位置が異なる画像同士で比較する。そして、画像候補決定部１２１は、ウェーブレット係数の絶対値を、最も大きいものからＭ番目に大きいものまで抽出すると共に、抽出した値に対応するＭ個の多重解像度表現を抽出する。例えば、図６（ａ）の多重解像度表現ＩＭ４＿１〜ＩＭ＿４において、二点鎖線Ｒ１、Ｒ２及びＲ３で囲まれた各位置のウェーブレット係数を、Ｒ１〜Ｒ３のそれぞれにおいて比較する。これにより、ウェーブレット係数の絶対値が大きいものからＭ個を抽出した結果と、Ｍ個のウェーブレット係数に対応する多重解像度表現とが、一例として、図６（ｂ）に表されている。
【００５０】
本参考形態は、多重解像度変換部１０２と画像候補決定部１２１とが以下のように各自の処理を協働して実行することにより、ウェーブレット係数を格納するためのメモリとして、候補となる値を格納するのに必要な容量だけを用意すればよいように構成されている。
【００５１】
多重解像度変換部１０２は、上述のとおり、ｖ方向及びｗ方向に関して同じ位置のウェーブレット係数を、全ての多重解像度表現ＩＭ４＿１〜ＩＭ４＿Ｌに関して順に算出する。一方、画像候補決定部１２１は、多重解像度変換部１０２が算出した１個目〜Ｍ個目のウェーブレット係数を、候補格納用のメモリに順に格納していく。次に、画像候補決定部１２１は、多重解像度変換部１０２がＭ＋１個目のウェーブレット係数を算出した段階で、このＭ＋１個目のウェーブレット係数を既にメモリに格納したＭ個のウェーブレット係数と比較する。そして、画像候補決定部１２１は、メモリ内のＭ個の値のうち、Ｍ＋１個目の値より小さいものがあれば、そのメモリの値をＭ＋１個目の値に変更する。これを、多重解像度変換部１０２がＬ番目の多重解像度表現に関するウェーブレット係数を算出するまで繰り返すことで、画像候補決定部１２１は、ｖ方向及びｗ方向に関する当該位置について、Ｍ個のウェーブレット係数の候補を決定する。このように処理することにより、画像候補を決定するために必要なメモリ容量が、ｖ方向及びｗ方向に関する位置ごとにＭ個の値を格納可能な大きさに限定されている。
【００５２】
確率分布算出部１２２は、画像候補決定部１２１が抽出したウェーブレット係数に基づいて、確率分布Ｐｉ（ｘｉ）を算出する。確率分布算出部１２２は、画素位置ｉごとに、多重解像度表現における当該画素位置ｉに対応する位置である対応位置と、その対応位置の近傍の位置である近傍位置とにおけるＭ個の候補値を使用する。なお、Ｍ個の候補値は、各多重解像度表現に含まれる高周波成分に関するもののみが使用される。
【００５３】
例えば、図５のグレースケール画像ＩＭ３における点Ｐ１が注目画素であるとする。多重解像度表現ＩＭ４において点Ｐ１に対応する高周波成分における対応位置は、Ｑ１１〜Ｑ３３である。また、対応位置Ｑ１１〜Ｑ３３の近傍位置としては、一例として、各対応位置にｖ方向及びｗ方向のそれぞれに関して隣接する４つの位置と、斜め方向に隣接する４つの位置とからなる８つの位置が用いられる。例えば、Ｑ１３に関しては、図５に示すように、対応位置であるＱ１３自身とその近傍位置とからなる縦３つ横３つからなる合計９個の位置における候補値が用いられる。したがって、Ｑ１１〜Ｑ３３の全体としては、合計で８１個の対応位置及び近傍位置における候補値が用いられる。
【００５４】
確率分布算出部１２２は、対応位置又は近傍位置に対応するグレースケール画像の画素位置をｉとし、その候補値に対応する多重解像度表現をＩＭ＿ｘｉとしたとき、対応位置又は近傍位置の候補値を、確率値Ｐｉ（ｘｉ）に投票する。例えば、図７に示すように、ある対応位置の第１の候補値、第２の候補値、…、第Ｍの候補値が、多重解像度ＩＭ４＿２、ＩＭ＿４、…、ＩＭ＿ｎ（ｎ：Ｌ以下の自然数）に対応する候補値であったとする。このとき、この対応位置の第１の候補値、第２の候補値、…、第Ｍの候補値は、Ｐｉ（２）、Ｐｉ（４）、…、Ｐｉ（ｎ）に対して投票される。
【００５５】
投票は、具体的には、以下のような演算によって行われる。まず、Ｐｉ（１）、Ｐｉ（２）、…、Ｐｉ（Ｌ）は、全てのｉに関して、初期値０が与えられる。そして、投票する側の位置である対応位置又は近傍位置を投票位置（ｖ，ｗ）とし、投票される側の位置である対応位置を投票中心（ｖｃ，ｗｃ）とし、対応位置又は近傍位置における候補値をνｖｗとするとき、候補値νｖｗに応じた投票値は、下記の数式２のように算出される。パラメータλ、θは、状況に応じて設定されるが、一例として、λ＝０．１、θ＝１．０などと設定される。ｍｉｎ（α，β）は、αとβのうち、小さいものを表す。
【数２】
2014010803.tif 000004
【００５６】
確率分布算出部１２２は、各候補値に関して数式２によって算出された投票値を、当該画素位置ｉにおける、候補値に対応する多重解像度表現ＩＭ＿ｘｉに関する確率値Ｐｉ（ｘｉ）に加算する。なお、図５には、Ｑ１３が対応位置であるときの投票中心と、投票位置の一例とが示されている。確率分布算出部１２２は、画素位置ｉに関する全ての候補値について投票を完了すると、Ｐｉ（１）、Ｐｉ（２）、…、Ｐｉ（Ｌ）全体のノルムが１になるようにＰｉ（ｘｉ）を正規化する。確率分布算出部１２２は、全てのｉに関して投票を実行する。このように、対応位置のみならず、その近傍の近傍位置からも投票がなされると共に、投票中心との距離に応じた重みをつけて投票がなされる。このため、確率分布の算出の際、局所的な整合性を考慮することができる。
【００５７】
以上のように、焦点位置確率算出部１２０が、グレースケール画像の各画素位置ｉにおける確率分布Ｐｉ（ｘｉ）を算出する。Ｐｉ（ｘｉ）は、上記の投票により求められることから、Ｌ個の画像のうち、画素位置ｉの周辺の投票値が比較的大きい画像において、大きい値を有することとなる。投票値が大きいことは、ウェーブレット係数の高周波成分において振幅が大きいことに対応する。したがって、通常、Ｐｉ（ｘｉ）が最大値を取るｘｉが、図８（ａ）に示すように、焦点が合っている最も確からしい画像、つまり、焦点位置に対応する画像を示すと考えられる。ところが、撮像時に画像ノイズ、飽和、ぼけなどの問題が生じた場合、上述のように算出したＰｉ（ｘｉ）において最大値を取るｘｉが、図８（ｂ）に示すように、実際の焦点位置から大きく離れた画像を示す場合がある。この場合、Ｐｉ（ｘｉ）の値をそのまま使用して焦点位置を見積もると、焦点が合っていない位置を焦点が合っているものと誤って見積もってしまうおそれがある。
【００５８】
そこで、信頼度設定部１０３は、Ｐｉ（１）、Ｐｉ（２）、…、Ｐｉ（Ｌ）のうちの最大値と所定の閾値との大小を比較する。そして、最大値が所定の閾値以下となる場合に、最大値が所定の閾値より大きい場合と比べて、その画素位置の信頼度を低く設定する。具体的には、信頼度設定部１０３は、Ｐｉ（ｘｉ）の画素位置ｉにおける最大値ＰｉＭＡＸが所定の閾値εを超えた場合（図８（ａ）の場合）、画素位置ｉの信頼度ＲｉをＲｉ＝１＋（ＰｉＭＡＸ−ε）と設定する。また、図８（ａ）に示すように、Ｐｉ（ｘｉ）の値は変更しない。一方、Ｐｉ（ｘｉ）の画素位置ｉにおける最大値ＰｉＭＡＸが所定の閾値ε以下であった場合（図８（ｂ）の場合）、信頼度設定部１０３は、画素位置ｉが問題領域内であると判定し、Ｒｉ＝１と設定する。さらに、Ｐｉ（１）、Ｐｉ（２）、…、Ｐｉ（Ｌ）の値は、図８（ｂ）に示すように、全て１／Ｌになるような均一分布に変更して設定する。
【００５９】
最適焦点位置取得部１３０は、焦点位置確率算出部１２０が算出したＰｉ（ｘｉ）及び信頼度設定部１０３が設定したＲｉに基づいて、画素位置ｉごとに、焦点位置が異なるＬ個の画像のうち、いずれの画像が焦点の合っている最適な焦点位置に対応する画像かを取得する。最適な焦点位置に対応する画像を取得する際の問題は、画像における大域的な領域に亘る整合性をいかに保証するかにある。そこで、本参考形態では、かかる問題解決の手段として、以下のように、マルコフ確率場による定式化を採用する。最適な焦点位置に対応する画像を決めることは、各画素位置ｉに関して、１番目〜Ｌ番目のいずれの画像が焦点位置に対応する画像かを決めることに相当する。つまり、１〜Ｌの任意の自然数を取るｘｉに対して、最適な焦点位置に対応する画像を与えるｘｉ＝ｘｉ＊をどのように割り当てるかということに相当する。
【００６０】
マルコフ確率場によると、ある割り当て｛ｘｉ｝を与える全体の確率が、各画素位置ｉに隣接する画素位置ｊにおける割り当てと画素位置ｉにおける割り当てとの関係のみから決定される。本参考形態においては、図９に示すように、ｘ方向及びｙ方向のそれぞれに関して各画素位置ｉに隣接する画素位置ｊ１〜ｊ４の４つにおける割り当てと画素位置ｉにおける割り当てとの関係のみから決定されるものとする。
【００６１】
このようなマルコフ確率場において、最適な割り当て｛ｘｉ＊｝を与える問題は、下記の数式３で表される評価関数Ｅ（Ｘ）を最小化するＸを決める問題に帰着することが知られている。Ｅ（Ｘ）は、割り当てＸ＝（ｘ１，ｘ２，…，ｘＮ）に関する関数であることを示している。Ｎは、グレースケール画像の総画素数である。
【数３】
2014010803.tif 000005
【００６２】
ここで、Ａは、互いに隣接する一対の画素からなる集合を示している。つまり、ｊは、Ｌ以下の自然数であって、ｉに隣接する画素位置を示す添え字である。また、Ｄｉ（ｘｉ）はコスト関数であり、画素位置ｉに関してｘｉを割り当てるときのコストを示す。本参考形態において、Ｄｉ（ｘｉ）は、以下の数式４によって与えられる。このように、Ｄｉ（ｘｉ）は、信頼度設定部１０３が設定した信頼度Ｒｉと、焦点位置確率算出部１２０が算出した確率値Ｐｉ（ｘｉ）との積で表される。
【数４】
2014010803.tif 000006
【００６３】
Ｖ（ｘｉ，ｘｊ）はペナルティ関数であり、互いに隣接する一対の画素位置ｉ、ｊに、それぞれｘｉ、ｘｊが割り当てられたときのペナルティを示す。なお、Ｖ（ｘｉ，ｘｊ）は、通常、互いに隣接する画素の割り当てが異なるほど大きいペナルティを与える関数であることから、平滑化項とも呼ばれる。本参考形態において、Ｖ（ｘｉ，ｘｊ）は、以下の数式５によって与えられる。このように、Ｖ（ｘｉ，ｘｊ）は、割り当てられたｘｉ，ｘｊ同士の差の絶対値に比例する値を有している。
【数５】
2014010803.tif 000007
【００６４】
Ｅ（Ｘ）を最小化する問題の解法には、確率伝播法（Belief Propagation）が用いられる。確率伝播法によって、多くの場合に大域的な最適解に近い解が近似的に得られることが知られており、例えば、ステレオ視のための視差計算や、いわゆる画像領域分割の解法として用いられている。確率伝播法に基づいてＥ（Ｘ）の最小化問題を解くため、最適焦点位置取得部１３０は、階層化演算部１３１、メッセージ算出部１３２（メッセージ算出手段）及び最適焦点位置算出部１３３（焦点位置算出手段）を有している。
【００６５】
まず、階層化演算部１３１は、上記の数式４に基づき、コスト関数Ｄｉ（ｘｉ）を全ての画素位置ｉに関して算出する。次に、コスト関数Ｄｉ（ｘｉ）をダウンサンプリングする。例えば、ｘ方向に関して１／２、ｙ方向に関して１／２の画素数のコスト関数Ｄｉ１（ｘｉ１）を算出する。ダウンサンプリングは、ｘ方向及びｙ方向に関して１つおきにコスト関数Ｄｉ（ｘｉ）から値を取る間引き処理により実行してもよい。あるいは、注目画素とその周辺画素とにおける複数のＤｉ（ｘｉ）の値の和や加重平均を取ることを、ｘ方向及びｙ方向に関して１つおきに実施してもよい。さらに、Ｄｉ１（ｘｉ１）を同様にダウンサンプリングすることにより、Ｄｉ２（ｘｉ２）を取得する。このように、階層化演算部１３１は、ダウンサンプリングをｎ回（ｎ：自然数）実行することにより、複数階層の低解像度のコスト関数Ｄｉ１（ｘｉ１）、Ｄｉ２（ｘｉ２）、…、Ｄｉｎ（ｘｉｎ）を取得する。以下、ダウンサンプリングしていない元のコスト関数を、低解像度のコスト関数の表記と合わせて、Ｄｉ０（ｘｉ０）と表す。
【００６６】
メッセージ算出部１３２は、階層化演算部１３１が算出したコスト関数Ｄｉ０（ｘｉ０）、Ｄｉ１（ｘｉ１）、Ｄｉ２（ｘｉ２）、…、Ｄｉｎ（ｘｉｎ）に基づいて、下記の数式６で表されるメッセージｍｔｉ→ｊｋを繰り返し計算する。ｍｔｉ→ｊｋ［ｘｊｋ］は、ｔ回目の繰り返しにおいて、画素位置ｉから隣接する画素位置ｊへの伝播メッセージ値を示す。ｘｊｋは１〜Ｌのいずれかを取る変数であるため、ｍｔｉ→ｊｋもＬ個の要素を持つ。Ｎｅ（ｉ）＼ｊは、画素位置ｉに隣接する画素位置であって画素位置ｊ以外のものからなる集合を意味する。ｋは、階層を示す０以上且つｎ以下の整数である。
【数６】
2014010803.tif 000008
【００６７】
メッセージ算出部１３２は、まず、ｋ＝ｎのときのメッセージを、ｔ＝０，１，２，…の順に計算する。ｔ＝０のときの初期値ｍ０ｉ→ｊｎは０（ゼロ）に設定される。ここで、ｍｉｎ（α）の下にｘｉｋが記載された演算子は、ｘｉｋに関するαの最小値を与える。ｋ＝ｎのときのメッセージは、以下のように計算される。なお、ｔ（ｔ≧１）のメッセージを算出するときには、ｔ−１以下のメッセージが全て計算されているとする。まず、下記の数式７に示すように、あるｘｊｎに関して、ｘｉｎ＝１、２、…、Ｌを順に代入することで、数式６におけるｍｉｎのカッコ内を算出する。次に、数式７に基づいて算出した値のうち、最小値を取るものを抽出する。これが、当該ｘｊｎに関するｍｔｉ→ｊｎ［ｘｊｎ］の値となる。これを、ｘｊｎ＝１、２、…、Ｌに関して順に繰り返すことで、全てのｘｊｎに関してｍｔｉ→ｊｎ［ｘｊｎ］が算出されることとなる。このような作業を、全ての隣接するｉ，ｊについて繰り返すことで、あるｔにおける全てのメッセージｍｔｉ→ｊｎ［ｘｊｎ］を算出することができる。メッセージ算出部１３２は、このようなメッセージｍｔｉ→ｊｎの算出を、ｔ＝０，１，２，…の順に、ｍｔｉ→ｊｎが収束するｔ＝Ｔまで繰り返す。
【数７】
2014010803.tif 000009
【００６８】
メッセージ算出部１３２は、次に、ｋ＝ｎより１階層高解像度のｋ＝ｎ―１のときのメッセージｍｔｉ→ｊｎ−１を上記と同様に算出する。ここで、メッセージの初期値には、低解像度（ｋ＝ｎ）において算出されたメッセージの収束値ｍＴｉ→ｊｎが用いられる。収束値ｍＴｉ→ｊｎは、低解像度のときの画素位置ｉ→画素位置ｊの経路に対応する高解像度ｋ＝ｎ−１の全ての経路におけるメッセージの初期値として用いられる。例えば、図１０に示すように、低解像度における１画素に対応する範囲が、高解像度における縦２×横２の４つの画素の範囲に対応したとする。このとき、低解像度における破線の矢印に示す経路に対応するメッセージの収束値ｍＴｉ→ｊｎが、この経路に対応する高解像度における４本の実線の矢印に示す経路に対応するメッセージの初期値として共通に使用される。そして、メッセージ算出部１３２は、ｋ＝ｎ―１としたときのメッセージｍｔｉ→ｊｎ−１を、ｔ＝０，１、２、…の順に、収束するまで繰り返す。この収束値は、上記と同様、さらに１階層高解像度（ｋ＝ｎ−２）のメッセージの初期値として用いられる。以上のように、メッセージ算出部１３２は、ｋ＝ｎ，ｎ−１，…，２，１，０の順にメッセージｍｔｉ→ｊｋの算出を繰り返す。
【００６９】
そして、最適焦点位置算出部１３３は、最終的に得られたメッセージの収束値ｍＴｉ→ｊ０を用いて、最適な焦点位置に対応する画像を示すｘｉ＊を、下記の数式８に基づいて算出する。
【数８】
2014010803.tif 000010
【００７０】
後段処理部１４０は、最適焦点位置取得部１３０が取得したｘｉ＊に基づいて、全画素において焦点が合った全焦点画像を取得する。後段処理部１４０は、ノイズ除去部１４１及び全焦点画像取得部１４２（合焦画像合成手段）を有している。ノイズ除去部１４１は、最適焦点位置取得部１３０が取得したｘｉ＊に対して、さらに、ガウシアンフィルタやメディアンフィルタを用いてノイズ除去処理を施す。ｘｉ＊は上記のとおり、大域的な最適化が可能な確率伝播法に基づいて取得されるため、低周波数の推定誤差が生じにくい。したがって、ここで使用されるフィルタとしては、比較的少ないフィルタカーネルのもので十分である。
【００７１】
全焦点画像取得部１４２は、ノイズ除去部１４１がノイズを除去したｘｉ＊に基づき、全焦点画像を取得する。ｘｉ＊は、画素位置ｉにおいて、Ｌ個の画像のうち、いずれの画像が焦点の合っているものに相当するかを示す。これに従い、全焦点画像取得部１４２は、図３の画像ＩＭ２のうち、各画素位置ｉにおいて、ｘｉ＊番目の画像に対応する画像の画素値を取得する。そして、全焦点画像取得部１４２は、画素位置ｉごとに取得した画素を対応する位置に配置し直すことにより、全焦点画像を取得する。
【００７２】
以下、顕微鏡１及び画像処理部１００が実行する処理の全体の流れについて、図１１を参照しつつ説明する。図１１の各ステップに従って実行される処理が、本発明の全焦点画像生成方法に対応する。また、これらの処理は、上記のとおり、画像処理部１００を構成するハードウェアとソフトウェアとが協働して実行してもよい。この場合のソフトウェアが、本発明の全焦点画像生成プログラムに対応する。
【００７３】
まず、顕微鏡１において、ステージ４の位置を移動させることにより、焦点位置が異なるＬ個の被写体画像の撮影が行われる（ステップＳ１）。これらの被写体画像は、各焦点位置に対応したステージ４の位置情報と共に、画像処理部１００へと送られる。次に、画像処理部１００において、位置ずれ補正部１１１が、被写体の水平方向の位置ずれに対応する画像中のｘ方向及びｙ方向に関する位置ずれを補正する（ステップＳ２）。次に、撮像倍率補正部１１２が、被写体画像ごとに顕微鏡１から入力されるステージ４の位置に基づいて被写体画像の倍率を補正する（ステップＳ３）。次に、ＲＧＢ画像取得部１１３が、位置ずれ補正及び倍率補正が施されたＬ個の画像ＩＭ１に基づいて、各画素においてＲ，Ｇ及びＢの各色の成分値が全て決定されたＬ個の画像ＩＭ２を生成する（ステップＳ４）。
【００７４】
次に、グレースケール変換部１０１が、Ｌ個の画像ＩＭ２から、Ｌ個のグレースケール画像ＩＭ３を生成する（ステップＳ５（グレースケール画像取得ステップ）;図３参照）。次に、多重解像度変換部１０２が、グレースケール画像ＩＭ３にウェーブレット変換を施す（ステップＳ６）。このとき、多重解像度変換部１０２は、ｖ方向及びｗ方向に関する位置ごとに、多重解像度表現ＩＭ４＿１〜ＩＭ４＿ＬのＬ個のウェーブレット係数を順に算出する（図６参照）。一方、画像候補決定部１２１は、多重解像度変換部１０２が、最初のＭ個を算出するまで、算出されたウェーブレット係数を順にメモリに格納する（ステップＳ７）。
【００７５】
次に、多重解像度変換部１０２が、Ｍ＋１個目以降のウェーブレット係数を算出する（ステップＳ８）ごとに、画像候補決定部１２１が、既にメモリに格納したＭ個の値のうち、新たに算出されたウェーブレット係数より小さいものがないか判定する（ステップＳ９）。そして、小さいものがあれば（ステップＳ９、Ｙｅｓ）、画像候補決定部１２１は、そのメモリの値を新たに算出されたウェーブレット係数に更新する（ステップＳ１０）。一方、小さいものがなければ（ステップＳ９、Ｎｏ）、画像候補決定部１２１は、当該位置に関して既にＬ個のウェーブレット係数が算出されたか否かを判定する（ステップＳ１１）。まだＬ個に達していなければ（ステップＳ１１、Ｎｏ）、ステップＳ８の処理に戻る。一方、Ｌ個に達していれば（ステップＳ１１、Ｙｅｓ）、画像候補決定部１２１は、当該位置でのウェーブレット係数の候補を、この時点でメモリに格納されているＭ個の値に決定する（ステップＳ１２）。なお、ステップＳ６及びＳ８が、本発明における多重解像度画像生成ステップに対応する。また、ステップＳ７及びＳ９〜Ｓ１２が抽出ステップに対応する。
【００７６】
次に、画像候補決定部１２１は、ウェーブレット係数を算出するべき位置が残っているか否かを判定する（ステップＳ１３）。まだウェーブレット係数を算出するべき位置が残っていると判定された場合には（ステップＳ１３、Ｙｅｓ）、ステップＳ６の処理に戻る。ｖ方向及びｗ方向に関する全ての位置についてウェーブレット係数が算出されたと判定された場合には（ステップＳ１３、Ｎｏ）、ステップＳ１４の処理に移る。
【００７７】
次に、確率分布算出部１２２が、画像候補決定部１２１が決定した候補値に基づいて、各画素位置ｉにおける確率分布Ｐｉ（ｘｉ）を算出する（ステップＳ１４（確率算出ステップ））。なお、ステップＳ７、Ｓ９〜Ｓ１２及びＳ１４が、本発明における焦点位置確率算出ステップに対応する。次に、信頼度設定部１０３が、Ｐｉ（ｘｉ）に基づいて画素位置ｉごとに信頼度Ｒｉを設定する（ステップＳ１５（信頼度設定ステップ））。次に、階層化演算部１３１が、確率分布算出部１２２が算出した確率分布Ｐｉ（ｘｉ）と信頼度設定部１０３が設定した信頼度Ｒｉとに基づいて、コスト関数Ｄｉ（ｘｉ）を算出すると共に、低解像度のコスト関数Ｄｉｋ（ｘｉｋ）（ｋ＝１，２，…，ｎ）を取得する（ステップＳ１６（低解像度化ステップ））。
【００７８】
次に、メッセージ算出部１３２が、ｋ＝ｎ、ｎ−１、…、１、０の順にメッセージｍｔｉ→ｊｋ［ｘｊ］を算出する（ステップＳ１７）。メッセージ算出部１３２は、メッセージが収束するｔ＝Ｔまで繰り返しメッセージを算出していく（ステップＳ１８、Ｎｏ→ステップＳ１７）。メッセージが収束すると（ステップＳ１８、Ｙｅｓ）、メッセージ算出部１３２は、次の階層があれば（ステップＳ１９、Ｙｅｓ）、次の階層に関してメッセージを算出する（ステップＳ１７）。次の階層がなければ（ステップＳ１９、Ｎｏ）、最適焦点位置算出部１３３が、最終的に収束したメッセージｍＴｉ→ｊ０を用いて、最適な焦点位置に対応する画像を示すｘｉ＊を算出する（ステップＳ２０（焦点位置算出ステップ））。なお、ステップＳ１６〜Ｓ２０が、本発明における最適焦点位置取得ステップに対応する。また、ステップＳ１７及びＳ１８が、本発明におけるメッセージ算出ステップに対応する。
【００７９】
次に、ノイズ除去部１４１が、ｘｉ＊からノイズを除去し（ステップＳ２１）、その後、全焦点画像取得部１４２が、ノイズ除去部１４１がノイズを除去したｘｉ＊に基づいて、画素位置ｉごとに焦点が合った全焦点画像を取得する（ステップＳ２２（合焦点画像合成ステップ））。
【００８０】
以上説明した本参考形態によると、評価関数Ｅ（Ｘ）を最小化する問題を確率伝播法に基づいて解くことにより、最適な焦点位置に対応する画像を示すｘｉ＊を画素位置ｉごとに取得する。評価関数Ｅ（Ｘ）を最小化することは画像中の大域的な領域に係る問題である。また、確率伝播法は、このような大域的な問題を最適に解決するための確率論的手法である。したがって、本参考形態によると、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って適切に焦点が合った画像を取得することができる。
【００８１】
また、本参考形態においては、確率伝播法におけるメッセージの算出の際、低解像度でメッセージを算出した結果を初期値として、より高解像度でメッセージを算出する。このため、高解像度での算出の際、より最適解に近い値を初期値としてメッセージが更新されるので、さらに、局所的な解に陥りにくくなる。
【００８２】
（第２の参考形態）
以下、本発明の被写体高さ情報取得方法、被写体高さ情報取得装置及び被写体高さ情報取得プログラムに係る第２の参考形態について説明する。第２の参考形態の被写体高さ情報取得装置において、第１の参考形態の全焦点画像生成装置と異なるのは、後段処理部１００の代わりに図１２の後段処理部２４０が設けられることのみである。また、後段処理部２４０において、後段処理部１００と異なるのは、焦点画像取得部１４２の代わりに被写体高さ情報取得部２４２が設けられることのみである。第２の参考形態の被写体高さ情報取得装置は、図１１のステップＳ１〜Ｓ２１までの処理を実行した後、被写体高さ情報取得部２４２が、以下のとおり、ｘｉ＊に基づいて、被写体の位置を取得する。
【００８３】
ｘｉ＊が表す被写体画像は、フォーカスレンズ９ａに対し、被写体の表面において画素位置ｉに対応する領域がちょうど焦点が合うときに撮像された画像に対応する。そこで、被写体高さ情報取得部２４２は、顕微鏡１の制御部１１から送られたステージ４の位置情報に基づき、ｘｉ＊が表す画像が撮影されたときのステージ４の位置を取得する。このステージ４の位置から、被写体の各領域における鉛直方向に関する相対的な位置が取得できる。
【００８４】
例えば、図２において、フォーカスレンズ９ａに対し、光の経路に沿った方向に関してちょうど距離Ｄ１のときに焦点が合うとする。つまり、図２の状態のとき、サンプルＳの表面において、領域ｓ１に焦点が合っている。一方で、領域ｓ２に焦点が合うためには、領域ｓ２がフォーカスレンズ９ａに対し距離Ｄ１になるところまでサンプルＳが移動しなければならない。そこで、領域ｓ１に焦点が合っている状態におけるステージ４の位置と、領域ｓ２に焦点が合っている状態におけるステージ４の位置との差から、領域ｓ１の位置と領域ｓ２の位置との差Ｄ１−Ｄ２を取得することができる。被写体高さ情報取得部２４２は、このようにして、各画素位置ｉに対応する被写体の表面上の各領域の相対的な高さ情報を取得することにより、被写体における表面全体の相対的な高さ情報を取得する。
【００８５】
第２の参考形態においても、第１の参考形態と同様にｘｉ＊を取得することから、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って適切に被写体の位置を取得することができる。なお、第２の参考形態に係る被写体高さ情報取得装置が図１１のステップＳ１〜Ｓ２１を実行すると共に、被写体高さ情報取得部２４２がｘｉ＊に基づいて被写体の位置を取得する処理が、本発明の被写体高さ情報取得方法に対応する。また、このような処理を実行するためのプログラムが、本発明の被写体高さ情報取得プログラムに対応する。なお、第２の参考形態に使用される各パラメータは、第１の参考形態に使用されるパラメータと異なっていてもよい。
【００８６】
（本発明に係る実施形態）
以下、本発明に係る一実施の形態について説明する。本実施形態が、第１又は第２の参考形態と異なるのは、最適焦点位置取得部１３０の代わりに図１３に示す最適焦点位置取得部３３０を設けたことであり、その他については第１又は第２の参考形態と同様である。第１の参考形態における最適焦点位置取得部１３０を最適焦点位置取得部３３０に代えたものが、本発明の第１の観点に係る実施形態に対応する。また、第２の参考形態における最適焦点位置取得部１３０を最適焦点位置取得部３３０に代えたものが、本発明の第２の観点に係る実施形態に対応する。以下では、主に最適焦点位置取得部３３０の構成及び機能について説明し、それ以外の構成や条件等については第１又は第２の参考形態と同様であるため、説明を適宜省略する。本実施形態に係る全焦点画像生成装置においては、図１１のステップＳ１〜ステップＳ１５が実行された後、最適焦点位置取得部３３０が、確率分布Ｐｉ（ｘｉ）及び信頼度Ｒｉに基づいて最適な焦点位置を導出する。この目的のため、最適焦点位置取得部３３０は、グラフ演算部３３１、評価関数算出部３３２及び判定部３３３を有している。これらの各部が実行する演算処理について、図１４〜図１６を参照しつつ説明する。
【００８７】
図１４のステップＳ１０１において、グラフ演算部３３１は、変数Ｘ＝（ｘ１，ｘ２，…，ｘＮ）に初期値を割り当てる。初期値は、各画素位置ｉに関し、その初期値において確率分布が最大となるような値とする。つまり、画素位置ｉに関して割り当てた初期値をａｉとするとき、各画素位置について、ａｉ以外のＬ以下のいずれの自然数ｂに対しても、Ｐｉ（ａｉ）がＰｉ（ｂ）以上となる。
【００８８】
次に、グラフ演算部３３１は、フラグ変数ｆを偽に初期化する（ステップＳ１０２）。例えば、ｆが０（ゼロ）であるときに偽であり、ｆが１であるときに真であるとする。次に、グラフ演算部３３１は、Ｌ以下の自然数から、互いに異なる変数α及びβをそれぞれ選択する（ステップＳ１０３）。後述の通り、ステップＳ１０３〜Ｓ１０９の処理は、ステップＳ１１０の条件分岐により、全ての｛α，β｝の組み合わせに関して繰り返される。そこで、ステップＳ１０３では、今までに選択されていない｛α，β｝の組み合わせが選択される。選択は、あらかじめ決められた順番でなされてもよい。例えば、｛１，２｝→｛１，３｝→…→｛１，Ｌ｝→｛２，３｝→｛２，４｝→…→｛２，Ｌ｝→｛３，４｝→｛３，５｝→…→｛３，Ｌ｝→｛４，５｝→…の順に、全ての組み合わせに関して選択していってもよい。または、ランダムな順番で選択してもよい。
【００８９】
次に、グラフ演算部３３１は、選択された｛α，β｝に基づいて構築されるグラフＧの最小カットを導出する（ステップＳ１０４）。最小カットを導出する手順は以下の通りである。グラフＧは、図１５に示すように、ノード及びエッジからなる。ノードには、ｘｉ∈｛α，β｝を満たす各画素位置ｉに一対一に対応するノードＮと、αに対応する１個のターミナルノードＴαと、βに対応する１個のターミナルノードＴβとが含まれている。図１５に示すように、ｘｋ≠α且つｘｋ≠βを満たす画素位置ｋに関しては、グラフＧにノードとして含まれない。エッジは、ノード同士を繋ぐリンクであり、ｖ方向又はｗ方向に隣接する２つのノードＮ同士を繋ぐｎ−リンクと、ノードＮとターミナルノードＴα又はＴβを結ぶｔ−リンクとが含まれている。ｎ−リンクは、ｖ方向又はｗ方向に隣り合う２つの画素位置に対応するノード同士をそれぞれ結ぶように形成される。図１５に示すように、グラフにノードとして含まれない画素位置ｋに関しては、ｎ−リンクが形成されない。ｔ−リンクは、全てのノードとターミナルノードαとを一対一に結ぶと共に、全てのノードとターミナルノードβとを一対一に結ぶように形成される。
【００９０】
ｎ−リンク及びｔ−リンクのそれぞれにはコストが設定される。ｎ−リンクのコストには、そのリンクが繋いでいる２つのノードに対応する画素位置ｉ，ｊに対して、数式５のペナルティ関数Ｖ（ｘｉ，ｘｊ）によって与えられる値が設定される。ｔ−リンクのコストは、そのリンクがターミナルノードと繋いでいるノードに対応する画素位置ｉに対して、数式９が示すｔliによって与えられる。ここで、Ｄｉ（l）及びＶ（ｌ，ｘｊ）は、数式４及び５にｘｉ＝ｌを代入した値にそれぞれ対応する。ターミナルノードαとのリンクのコストを算出する場合はｌ＝αが用いられ、ターミナルノードβとのリンクのコストを算出する場合はｌ＝βが用いられる。Ｎｉは、画素位置ｉに隣接する画素位置の集合である。したがって、数式９のΣは、画素位置ｉに隣接する画素位置ｊに関し、ｘｊがαでもβでもないものについて足し合わせることを示す。
【数９】
2014010803.tif 000011
【００９１】
以上のように設定されるグラフに関し、グラフ演算部３３１は、最小カットを導出する。最小カットは、以下のように定義される。グラフＧを、各リンクを切断することにより、部分グラフＧ１と部分グラフＧ２とに分割することを想定する。部分グラフＧ１は、Ｔα及びノードＮ並びにこれら同士を結ぶエッジのみからなる。部分グラフＧ２は、グラフＧから部分グラフＧ１を分割した際の残りのノードを全て含み、Ｔβ及びノードＮ並びにこれら同士を結ぶエッジのみからなる。図１６は、その一例である。図１５には、図１６のように分割する場合に切断すべきリンクが「×」で示されている。
【００９２】
このとき、切断されるリンクに設定されたコストの総和が最小となるような切断が、最小カットとされる。つまり、切断されるｎ−リンクについては数式５を用いてコストが算出されると共に、切断されるｔ−リンクについては数式９を用いてコストが算出される。そして、算出されたコストが全て足し合わされることにより、コストの総和が求められる。このコストの総和が最小となるような切断が、最小カットである。
【００９３】
このような最小カットを導出するには、グラフＧを各リンクに対応するコストを容量とするネットワークグラフとみなした場合における、ターミナルノードαからターミナルノードβへの最大フローを利用すればよいことが知られている。グラフ演算部３３１は、このような最大フローを求める方法のいずれかを使用して、最小カットを導出する（導出ステップ）。一例としては、フォード・ファルカーソンのアルゴリズム（FORD-Fulkerson algorithm）が使用される。最大フローを求めることにより最小カットを導出する手順の概略は以下の通りである。ターミナルノードα及びターミナルノードβ間を結ぶ１本の経路であって、経路の各リンクに容量の空きがあるものに関して流れを形成できる。このような流れをフローに加えていくことを、いずれの経路にも空きがなくなるまで繰り返す。これによって得られた全フローの総和を取得する。これが、最大フローとなる。最大フローは最小カットのコストの総和と等しい。このため、切断によるコストの総和が最大フローと一致するようなカットの方法を探索することにより、最小カットを導出できる。
【００９４】
次に、グラフ演算部３３１は、最小カットに基づき、Ｘ＊＝（ｘ１＊，ｘ２＊，…，ｘＮ＊）を以下のように割り当てる（ステップＳ１０５；割り当てステップ）。Ｘにおいて、ｘｉ∈｛α，β｝を満たす画素位置ｉに関しては、最小カットにおいて部分グラフＧ１に含まれるノード（ターミナルノードβから切断されるノード）についてｘｉ＊＝αと割り当てると共に、最小カットにおいて部分グラフＧ２に含まれるノード（ターミナルノードαから切断されるノード）についてｘｉ＊＝βと割り当てる。また、ｘｋ≠α且つｘｋ≠βとなる画素位置ｋに関しては、ｘｋ＊＝ｘｋと割り当てる。
【００９５】
次に、評価関数算出部３３２が、数式３を用いて、Ｅ（Ｘ）及びＥ（Ｘ＊）をそれぞれ算出する（ステップＳ１０６；算出ステップ）。さらに、判定部３３３が、Ｅ（Ｘ）及びＥ（Ｘ＊）の大小を比較する（ステップＳ１０７；判定ステップ）。そして、Ｅ（Ｘ＊）がＥ（Ｘ）より小さいと判定部３３３が判定した場合（ステップＳ１０７、Ｙｅｓ）、最適焦点位置取得部３３０は、Ｘ＊を変数Ｘに代入する（ステップＳ１０８）と共に、フラグ変数ｆを真とする。これにより、フラグ変数ｆが真であることは、少なくとも１回、ＸをＸ＊に更新したことを示すこととなる。また、フラグ変数ｆが偽であることは、１回もＸをＸ＊に更新していないことを示すこととなる。そして、ステップＳ１１０の処理に進む。一方、ステップＳ１０７において、Ｅ（Ｘ＊）がＥ（Ｘ）以上であると判定部３３３が判定した場合（ステップＳ１０７、Ｎｏ）、そのままステップＳ１１０の処理に進む。
【００９６】
次に、最適焦点位置取得部３３０は、｛α，β｝の全ての組み合わせについてステップＳ１０３〜Ｓ１１０を実行したか否かを判定する（ステップＳ１１０）。そして、まだ実行していない組み合わせがあると判定した場合（ステップＳ１１０、Ｎｏ）には、最適焦点位置取得部３３０は、ステップＳ１０３の処理に戻る。全ての組み合わせについて実行したと判定した場合（ステップＳ１１０、Ｙｅｓ）には、最適焦点位置取得部３３０は、ステップＳ１１１の処理に進む。ステップＳ１１１において、最適焦点位置取得部３３０は、フラグ変数ｆが真であるか否かを判定する。フラグ変数ｆが真、つまり、Ｘが少なくとも１回は更新されたと判定した場合（ステップＳ１１１、Ｙｅｓ）には、最適焦点位置取得部３３０は、ステップＳ１０２の処理に戻る。一方、フラグ変数ｆが偽、つまり、Ｘが一度も更新されなかったと判定した場合（ステップＳ１１１、Ｎｏ）には、最適焦点位置取得部３３０は、ステップＳ１１２の処理に移る。ステップＳ１１２において、最適焦点位置取得部３３０は、最新のＸを最適な焦点位置とすると共に、ステップＳ２１及びＳ２２を実行する。
【００９７】
以上説明した本実施形態によると、ステップＳ１０５においてＸ＊を割り当て、ステップＳ１０７においてＥ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定すると共に、これを満たす場合には、ステップＳ１０８において、Ｘ＊を新たなＸとする。そして再び、ステップＳ１０３〜Ｓ１０７の処理を実行する。これにより、現状のＸより評価関数が小さいＸ＊が見つかれば、これをあらためてＸとすることを繰り返すので、評価関数の大きさは、隣接する画素位置同士の焦点位置の割り当て方に依存するため、ひいては、画素位置同士を繋いでいくことで構成される広範囲の割り当て方に依存する。このため、評価関数が最小となるような焦点位置を取得することにより、第１の観点に応じた実施形態では、大域的な領域に亘って適切に焦点が合った画像を取得できる。また、第２の観点に応じた実施形態では、画像の大域的な領域に亘って適切に被写体の高さ情報を取得できる。
【００９８】
また、本実施形態においては、ステップＳ１０３〜Ｓ１０５において、グラフの最小カットに対応するＸ＊を求める。これにより、Ｅ（Ｘ＊）がＥ（Ｘ）より小さくなることが期待できる。したがって、評価関数が最小（極小）となるような割り当てに到達しやすくなる。
【００９９】
また、本実施形態においては、一のα及びβの組み合わせに関してステップＳ１０３〜１０９を実行した後、これらのステップを別のα及びβの組み合わせに関して再び実行する。したがって、種々のα及びβの組み合わせに関して割り当てＸ＊を求めることで、評価関数が小さくなるような割り当てを導出できる。
【０１００】
なお、本実施形態においては、ｎ−リンクのコストをＶ（ｘｉ，ｘｊ）とし、ｔ−リンクのコストを数式９としている。これにより、最小カットを導出することで、数式３で表される評価関数Ｅ（Ｘ＊）が小さくなるようなＸ＊を割り当てることができる。
【０１０１】
（変形例）
以上は、本発明の一実施の形態についての説明であるが、本発明は上述の実施の形態に限られるものではなく、課題を解決するための手段に記載された範囲の限りにおいて、以下に例を挙げるとおり、様々な変更が可能なものである。
【０１０２】
［１］上述の実施の形態では、グレースケール変換部１０１が、ＲＧＢ画像データからグレースケール画像データを生成する。しかし、撮像素子１０が直接グレースケール画像を生成してもよい。また、撮像素子１０が生成したＧ画素のみからなる画像をグレースケール画像として使用してもよい。
【０１０３】
［２］上述の実施の形態では、画像候補決定部１２１が、確率分布を算出するための候補となるＭ個（Ｍ＜Ｌ）の多重解像度表現を抽出する。しかし、Ｌ個全ての多重解像度表現から確率分布を算出してもよい。ただし、本発明者の知見によると、全ての多重解像度表現を使用する場合、画質が悪化するおそれがあるため、Ｍ個の候補を抽出したほうが好ましい。
【０１０４】
［３］上述の実施の形態においては、確率分布を算出するため、対応位置及び近傍位置のそれぞれにおけるＭ個の候補値から所定の重みを用いて算出した投票値を投票している。しかし、近傍位置から投票せず、対応位置のみから投票してもよい。
【０１０５】
［４］上述の参考形態においては、コスト関数Ｄｉ（ｘｉ）をダウンサンプリングした低解像度のコスト関数に基づいてメッセージを算出する。しかし、コスト関数をダウンサンプリングせず、そのまま用いてメッセージ関数を算出してもよい。また、コスト関数として、Ｐｉ（ｘｉ）に信頼度Ｒｉを乗算したものでなく、Ｐｉ（ｘｉ）そのものを使用してもよい。
【０１０６】
［５］上述の実施の形態においては、顕微鏡１において撮影された被写体画像を用いる場合に本発明が適用されている。しかし、デジタルカメラや、デジタルビデオカメラ、マイクロスコープ等、顕微鏡以外の光学系にも適用可能である。
【０１０７】
［６］上述の実施の形態においては、全焦点画像生成装置や被写体高さ情報取得装置が、顕微鏡１（光学系）と画像処理部１００とを有している。しかし、本発明が、光学系とは独立に構成された画像処理装置として実現されてもよい。この場合、例えば、画像処理装置とは独立に構成された撮像装置において、焦点位置が異なる複数の画像が撮像された後、これらの画像に対応するデータが光学的記録媒体などの記録媒体に記録される。画像処理装置は、この記録媒体からデータを読み出す読み出し部を有しており、記録媒体から画像データを読み出すと共に、読み出した画像データに対して上述の画像処理部と同様に画像処理を施す。
【０１０８】
［７］上述の実施形態においては、大域最適化の具体的な方法としてグラフカットを利用している。上述の通り、グラフカットを使用することにより、おおよそＥ（Ｘ＊）＜Ｅ（Ｘ）となるＸ＊を割り当てることが期待できる。しかし、Ｘ＊を割り当てる際、この他の方法を採用してもよい。例えば、第１の参考形態における最適焦点位置取得部１３０が実施する確率伝播法に基づき、Ｘ＊が割り当てられてもよい（数式８参照）。これによると、Ｅ（Ｘ＊）を当初から最小値に近接させることが期待できる。例えば、図１４のステップＳ１０１の後、確率伝播法（図１１のステップＳ１６〜Ｓ２０、数式８）に基づき、Ｘ＊を割り当てる。そして、Ｅ（Ｘ＊）＜Ｅ（Ｘ）を満たすか否かを判定し、これを満たす場合には、Ｘ＊をＸに代入する。その後、ステップＳ１０２以降の処理を実行する。これにより、Ｘの初期値として、評価関数が最小値に近いものを使用できる。または、平均場近似を用いてＸ＊を割り当てても、Ｅ（Ｘ＊）＜Ｅ（Ｘ）となることが期待できる。
【０１０９】
以上のように、Ｘ＊を割り当てる方法として、Ｅ（Ｘ＊）＜Ｅ（Ｘ）となることが期待できる種々の方法が採用され得る。これらの方法において、Ｅ（Ｘ＊）＜Ｅ（Ｘ）となるようなＸ＊を割り当てられるかどうかの確度はさまざまである。しかしながら、いずれの方法が用いられるにせよ、“Ｘ＊の割り当て→Ｅ（Ｘ＊）＜Ｅ（Ｘ）の判定→これを満たす場合にＸにＸ＊を代入”との処理を繰り返すことにより、Ｅ（Ｘ）が最小（極小）となるような最適なＸに確実に近づけることができる。
【符号の説明】
【０１１０】
１…顕微鏡、１００…画像処理部、１１０…前段処理部、１０１…グレースケール変換部、１０２…多重解像度変換部、１２０…焦点位置確率算出部、１０３…信頼度設定部、３３０…最適焦点位置取得部、１４０…後段処理部
【図１】
2014010803.tif 000013
【図２】
2014010803.tif 000014
【図３】
2014010803.tif 000015
【図４】
2014010803.tif 000016
【図５】
2014010803.tif 000017
【図６】
2014010803.tif 000018
【図７】
2014010803.tif 000019
【図８】
2014010803.tif 000020
【図９】
2014010803.tif 000021
【図１０】
2014010803.tif 000022
【図１１】
2014010803.tif 000023
【図１２】
2014010803.tif 000024
【図１３】
2014010803.tif 000025
【図１４】
2014010803.tif 000026
【図１５】
2014010803.tif 000027
【図１６】
2014010803.tif 000028