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