"焦点位置の異なる、L個L:2以上の自然数の被写体画像を撮像する撮像ステップと、前記撮像ステップにおいて撮像された像に基づいて、前記焦点位置ごとに前記被写体のグレースケール画像を取得するグレースケール画像取得ステップと、前記グレースケール画像取得ステップにおいて取得された前記グレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成ステップと、前記多重解像度画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、前記撮像ステップにおいて撮像された前記焦点位置ごとの像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成ステップと、を備えていることを特徴とする全焦点画像生成方法。前記最適焦点位置取得ステップにおいて、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置としてxiを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出することを特徴とする請求項1に記載の全焦点画像生成方法。前記焦点位置確率算出ステップが、前記多重解像度画像生成ステップにおいて生成された、前記焦点位置が異なる複数の前記解像度分解画像の変換係数同士を前記解像度分解画像中の画素ごとに比較し、比較した変換係数をその絶対値が最も大きいものから大きい順に所定の数だけ抽出すると共に、抽出した各変換係数に対応する前記解像度分解画像における前記焦点位置を抽出する抽出ステップと、前記複数の周波数成分のうち、低周波成分を除いた、複数の高周波成分のそれぞれにおける、画素位置iに対応する前記解像度分解画像中の位置である対応位置と、当該対応位置の近傍の位置である近傍位置とに関して、前記抽出ステップにおいて抽出された各変換係数に応じた重みを、その変換係数に対応する前記焦点位置における画素位置iに関する確率値として加算する確率算出ステップとを含んでいることを特徴とする請求項2に記載の全焦点画像生成方法。前記近傍位置をv,wとし、前記対応位置をvc,wcとし、νvwを前記近傍位置v,wに関して前記抽出ステップで抽出された変換係数とし、λ及びθを所定の定数とするとき、前記変換係数に応じた重みが、下記の数式1で表される値に比例することを特徴とする請求項3に記載の全焦点画像生成方法。前記最適焦点位置取得ステップが、下記の数式2で表されるメッセージmti→jを、全ての隣り合うi,j及びxjに関してt=0,1,2,...,Tの順に、メッセージが収束するt=Tまで算出するメッセージ算出ステップと、下記の数式3により前記最適な焦点位置xiを算出する焦点位置算出ステップとを含んでいることを特徴とする請求項2~4のいずれか1項に記載の全焦点画像生成方法。前記最適焦点位置取得ステップが、元の解像度でDixiにより算出されたコスト値から1回又は複数回のダウンサンプリングにより低解像度のコスト値を導出する低解像度化ステップと、前記低解像度化ステップにおいて導出された低解像度のコスト値とペナルティ関数に基づいて前記メッセージ算出ステップを実行する低解像度メッセージ算出ステップと、前記低解像度メッセージ算出ステップにおいて収束するまで算出されたメッセージを初期値として、より高解像度のコスト値に基づいて前記メッセージ算出ステップを実行する高解像度メッセージ算出ステップとを含んでいることを特徴とする請求項5に記載の全焦点画像生成方法。前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、前記画素位置ごとの信頼度を、当該画素位置における前記焦点位置確率分布の最大値が所定の閾値以下である場合に、前記最大値が所定の閾値を超える場合よりも低くなるように設定する信頼度設定ステップをさらに備え、Dixiが、前記信頼度設定ステップで設定された画素位置iにおける前記信頼度が低いほど大きくなるような関数であることを特徴とする請求項5又は6に記載の全焦点画像生成方法。Pixiを前記焦点位置確率分布の画素位置i、前記焦点位置xiにおける確率値とし、Riを前記信頼度設定ステップで設定された画素位置iにおける前記信頼度とするとき、Dixi及びVxi,xjが、下記の数式4及び5によって表されることを特徴とする請求項7に記載の全焦点画像生成方法。焦点位置の異なる、L個L:2以上の自然数の被写体画像から、被写体の合焦点画像を生成する全焦点画像生成装置であって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度変換画像生成手段と、前記多重解像度変換画像生成手段が生成した前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出手段と、前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得手段と、前記最適焦点位置取得手段が取得した前記最適な焦点位置に基づいて、前記焦点位置ごとの被写体の像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成手段と、を備えていることを特徴とする全焦点画像生成装置。前記最適焦点位置取得手段が、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置xiを割り当てた場合に、前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出することを特徴とする請求項9に記載の全焦点画像生成装置。前記焦点位置確率算出手段が、前記多重解像度変換画像生成手段が生成した、前記焦点位置が異なる複数の前記解像度分解画像の変換係数同士を前記解像度分解画像中の画素ごとに比較し、比較した変換係数をその絶対値が最も大きいものから大きい順に所定の数だけ抽出すると共に、抽出した各変換係数に対応する前記解像度分解画像における前記焦点位置を抽出する抽出手段と、前記複数の周波数成分のうち、低周波成分を除いた、複数の高周波成分のそれぞれにおける、画素位置iに対応する前記解像度分解画像中の位置である対応位置と、当該対応位置の近傍の位置である近傍位置とに関して、前記抽出ステップにおいて抽出された各変換係数に応じた重みを、その変換係数に対応する前記焦点位置における画素位置iに関する確率値として加算する確率算出手段とを含んでいることを特徴とする請求項10に記載の全焦点画像生成装置。前記最適焦点位置取得手段が、下記の数式6で表されるメッセージmti→jを、全ての隣り合うi,j及びxjに関してt=0,1,2,...,Tの順に、メッセージが収束するt=Tまで算出するメッセージ算出手段と、下記の数式7により前記最適な焦点位置xiを算出する焦点位置算出手段とを含んでいることを特徴とする請求項10又は11に記載の全焦点画像生成装置。焦点位置の異なる、L個L:2以上の自然数の被写体画像から、被写体の合焦点画像を生成するようにコンピュータを機能させるプログラムであって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度変換画像生成ステップと、前記多重解像度変換画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、前記撮像ステップにおいて撮像された前記焦点位置ごとの像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成ステップと、をコンピュータに実行させることを特徴とする全焦点画像生成プログラム。前記最適焦点位置取得ステップにおいて、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置としてxiを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出するようにコンピュータを機能させることを特徴とする請求項13に記載の全焦点画像生成プログラム。前記焦点位置確率算出ステップが、前記多重解像度変換画像生成ステップにおいて生成された、前記焦点位置が異なる複数の前記解像度分解画像の変換係数同士を前記解像度分解画像中の画素ごとに比較し、比較した変換係数をその絶対値が最も大きいものから大きい順に所定の数だけ抽出すると共に、抽出した各変換係数に対応する前記解像度分解画像における前記焦点位置を抽出する抽出ステップと、前記複数の周波数成分のうち、低周波成分を除いた、複数の高周波成分のそれぞれにおける、画素位置iに対応する前記解像度分解画像中の位置である対応位置と、当該対応位置の近傍の位置である近傍位置とに関して、前記抽出ステップにおいて抽出された各変換係数に応じた重みを、その変換係数に対応する前記焦点位置における画素位置iに関する確率値として加算する確率算出ステップとを含んでいることを特徴とする請求項14に記載の全焦点画像生成プログラム。前記最適焦点位置取得ステップが、下記の数式8で表されるメッセージmti→jを、全ての隣り合うi,j及びxjに関してt=0,1,2,...,Tの順に、メッセージが収束するt=Tまで算出するメッセージ算出ステップと、下記の数式9により前記最適な焦点位置xiを算出する焦点位置算出ステップとを含んでいることを特徴とする請求項14又は15に記載の全焦点画像生成プログラム。焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるL個L:2以上の自然数の被写体画像を撮像する撮像ステップと、前記撮像ステップにおいて撮像された像に基づいて、前記焦点位置ごとに前記被写体のグレースケール画像を取得するグレースケール画像取得ステップと、前記グレースケール画像取得ステップにおいて取得された前記グレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度変換画像生成ステップと、前記多重解像度変換画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得ステップと、を備えていることを特徴とする被写体高さ情報取得方法前記最適焦点位置取得ステップにおいて、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置としてxiを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出することを特徴とする請求項17に記載の被写体高さ情報取得方法。焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるL個L:2以上の自然数の被写体画像に基づいて、当該被写体の位置を取得する被写体高さ情報取得装置であって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度変換画像生成手段と、前記多重解像度変換画像生成手段が生成した前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出手段と、前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得手段と、前記最適焦点位置取得手段が取得した前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得手段と、を備えていることを特徴とする被写体高さ情報取得装置。前記最適焦点位置取得手段が、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置xiを割り当てた場合に、前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出することを特徴とする請求項19に記載の被写体高さ情報取得装置。焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるL個L:2以上の自然数の被写体画像に基づいて、当該被写体の高さ情報を取得するようにコンピュータを機能させるプログラムであって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度変換画像生成ステップと、前記多重解像度変換画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得ステップと、をコンピュータに実行させることを特徴とする被写体高さ情報取得プログラム。前記最適焦点位置取得ステップにおいて、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置としてxiを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出するようにコンピュータを機能させることを特徴とする請求項21に記載の被写体高さ情報取得プログラム。特許文献1~3のいずれも、上記のとおり、ウェーブレット値を利用したり、輝度値の変化量が大きいものを取り出したりして合焦位置の画素を抽出している。しかしながら、例えば、撮像時に画像ノイズ、飽和、ぼけなどの問題が生じた場合、撮像画像において、焦点が合っている領域ではコントラスト値が小さくなり、焦点が合っていない領域ではコントラスト値が大きくなることがある。このような問題が生じると、特許文献1や特許文献2を適用し、ウェーブレット値や輝度値の変化量の大小関係を比較して合焦位置を取得する場合、焦点が合っていない位置を合焦位置と誤って判定してしまうおそれがある。また、特許文献3は、画像を3つの領域に分け、領域ごとに異なる方法で焦点位置を取得することにより、上記の問題の解決を図っている。しかしながら、特許文献3によると、エッジ部分やその近傍から離隔した部分では、エッジ部分の画素とエッジ部分の近傍の画素とに基づく補間により合焦画素を取得する。したがって、補間対象となる画素がエッジ部分やその近傍から遠く離れていたり、これらの部分と深度方向に異なる位置にあったりすると、補間の際に生じる誤差が大きくなる。このため、仮にエッジ部分とその近傍からなる局所的な領域においては適切に焦点が合った画像を生成できたとしても、画像中の大域的な領域に亘って適切に焦点が合った画像を生成できないおそれがある。本発明の第1の目的は、焦点位置が異なる複数の撮像画像から画素位置ごとに焦点が合った画像を生成する方法であって、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って適切に焦点が合った画像を生成できる全焦点画像生成方法、並びに、当該方法に対応した全焦点画像生成装置及び全焦点画像生成プログラムを提供することにある。また、焦点位置が異なる複数の撮像画像から画素位置ごとに焦点が合っているものを抽出することは、焦点位置から換算される高さ情報を画素位置ごとに取得することと等しい。そして、上記のような画像ノイズ等による問題は、被写体の高さ情報を画素位置ごとに取得する場合にも全く同様の問題として表れる。本発明の第2の目的は、焦点位置が異なる複数の撮像画像に基づいて被写体の高さ情報を画素位置ごとに取得する方法であって、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って適切に被写体の高さ情報を取得できる被写体高さ情報取得方法、並びに、当該方法に対応した被写体高さ情報取得装置及び被写体高さ情報取得プログラムを提供することにある。第1の目的を達成するためになされた本発明の第1の観点によると、請求項1に記載の全焦点画像生成方法に係る発明は、焦点位置の異なる、L個L:2以上の自然数の被写体画像を撮像する撮像ステップと、前記撮像ステップにおいて撮像された像に基づいて、前記焦点位置ごとに前記被写体のグレースケール画像を取得するグレースケール画像取得ステップと、前記グレースケール画像取得ステップにおいて取得された前記グレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度画像生成ステップと、前記多重解像度画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、前記撮像ステップにおいて撮像された前記焦点位置ごとの像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成ステップと、を備えている。請求項1に記載の全焦点画像生成方法は、最適な焦点位置を、大域最適化を用いて近似的に導出する。そして、導出した最適な焦点位置を焦点が合った位置とし、当該焦点位置に対応する撮像画像から画素を取得して全焦点画像を生成する。最適な焦点位置を取得することは、画像中の大域的な領域に係る問題である。確率伝播法、平均場近似、グラフカットなどの大域最適化は、このような大域的な問題を最適に解決するための確率論的手法である。したがって、本発明によると、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って適切に焦点が合った画像を取得することができる。また、請求項1に記載の全焦点画像生成方法は、請求項2に記載の発明のように、前記最適焦点位置取得ステップにおいて、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置としてxiを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出してもよい。これによると、確率伝播法を用いて具体的に最適な焦点位置を導出できる。また、請求項2に記載の全焦点画像生成方法は、請求項3に記載の発明のように、前記焦点位置確率算出ステップが、前記多重解像度画像生成ステップにおいて生成された、前記焦点位置が異なる複数の前記解像度分解画像の変換係数同士を前記解像度分解画像中の画素ごとに比較し、比較した変換係数をその絶対値が最も大きいものから大きい順に所定の数だけ抽出すると共に、抽出した各変換係数に対応する前記解像度分解画像における前記焦点位置を抽出する抽出ステップと、前記複数の周波数成分のうち、低周波成分を除いた、複数の高周波成分のそれぞれにおける、画素位置iに対応する前記解像度分解画像中の位置である対応位置と、当該対応位置の近傍の位置である近傍位置とに関して、前記抽出ステップにおいて抽出された各変換係数に応じた重みを、その変換係数に対応する前記焦点位置における画素位置iに関する確率値として加算する確率算出ステップとを含んでいることが好ましい。これによると、解像度分解画像において各画素位置に対応する対応位置の近傍におけるウェーブレット係数に応じた重みを当該画素位置における確率値として加算するので、確率値に基づいて全焦点画像を生成する際、各画素位置とその近傍画素との局所的な整合性を考慮した画像を生成できる。また、請求項3に記載の全焦点画像生成方法は、請求項4に記載の発明のように、前記近傍位置をv,wとし、前記対応位置をvc,wcとし、νvwを前記近傍位置v,wに関して前記抽出ステップで抽出された変換係数とし、λ及びθを所定の定数とするとき、前記変換係数に応じた重みが、下記の数式1で表される値に比例するものとすることが好ましい。これによると、対応位置と近傍位置との距離に応じた適切な重みを算出し、各画素位置の確率値に加算できる。また、請求項2~4のいずれかに記載の全焦点画像生成方法は、請求項5に記載の発明のように、前記最適焦点位置取得ステップが、下記の数式2で表されるメッセージmti→jを、全ての隣り合うi,j及びxjに関してt=0,1,2,...,Tの順に、メッセージが収束するt=Tまで算出するメッセージ算出ステップと、下記の数式3により前記最適な焦点位置xiを算出する焦点位置算出ステップとを含んでいることが好ましい。これによると、メッセージを繰り返し算出することにより、確率伝播法に基づいて最適な焦点位置を適切に取得できる。また、請求項5に記載の全焦点画像生成方法は、請求項6に記載の発明のように、前記最適焦点位置取得ステップが、元の解像度でDixiにより算出されたコスト値から1回又は複数回のダウンサンプリングにより低解像度のコスト値を導出する低解像度化ステップと、前記低解像度化ステップにおいて導出された低解像度のコスト値とペナルティ関数に基づいて前記メッセージ算出ステップを実行する低解像度メッセージ算出ステップと、前記低解像度メッセージ算出ステップにおいて収束するまで算出されたメッセージを初期値として、より高解像度のコスト値に基づいて前記メッセージ算出ステップを実行する高解像度メッセージ算出ステップとを含んでいることが好ましい。これによると、高解像度メッセージ算出ステップでは、低解像度メッセージ算出ステップでのメッセージ算出結果を初期値としてメッセージを算出する。このため、高解像度では、より最適解に近い値を初期値としてメッセージが更新されていくため、局所的な解に陥りにくい。また、請求項5又は6に記載の全焦点画像生成方法は、請求項7に記載の発明のように、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、前記画素位置ごとの信頼度を、当該画素位置における前記焦点位置確率分布の最大値が所定の閾値以下である場合に、前記最大値が所定の閾値を超える場合よりも低くなるように設定する信頼度設定ステップをさらに備え、Dixiが、前記信頼度設定ステップで設定された画素位置iにおける前記信頼度が低いほど大きくなるような関数であることが好ましい。発明者の知見によると、撮像画像において、ノイズ等の撮像時の問題が発生している領域では、深度確率分布の最大値が比較的小さくなることが確認されている。これに基づき、上記のとおり、焦点位置確率分布の最大値と閾値とを比較し、最大値が閾値を超えない場合には閾値を超える場合と比べて信頼度を低く設定することにより、問題領域が合焦位置と誤って判定されないようにしている。また、請求項7に記載の全焦点画像生成方法は、請求項8に記載の発明のように、Pixiを前記焦点位置確率分布の画素位置i、前記焦点位置xiにおける確率値とし、Riを前記信頼度設定ステップで設定された画素位置iにおける前記信頼度とするとき、Dixi及びVxi,xjが、下記の数式4及び5によって表されることが好ましい。これによると、コスト関数が信頼度と焦点位置確率分布とに基づいて具体的に算出される。また、ペナルティ関数が、隣接する画素同士の焦点位置の差に基づいて具体的に算出される。また、第1の観点による請求項9に記載の全焦点画像生成装置は、焦点位置の異なる、L個L:2以上の自然数の被写体画像から、被写体の合焦点画像を生成する全焦点画像生成装置であって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度変換画像生成手段と、前記多重解像度変換画像生成手段が生成した前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出手段と、前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得手段と、前記最適焦点位置取得手段が取得した前記最適な焦点位置に基づいて、前記焦点位置ごとの被写体の像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成手段と、を備えている。また、請求項9に記載の全焦点画像生成装置は、請求項10に記載の発明のように、前記最適焦点位置取得手段が、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置xiを割り当てた場合に、前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出してもよい。また、請求項10に記載の全焦点画像生成装置は、請求項11に記載の発明のように、前記焦点位置確率算出手段が、前記多重解像度変換画像生成手段が生成した、前記焦点位置が異なる複数の前記解像度分解画像の変換係数同士を前記解像度分解画像中の画素ごとに比較し、比較した変換係数をその絶対値が最も大きいものから大きい順に所定の数だけ抽出すると共に、抽出した各変換係数に対応する前記解像度分解画像における前記焦点位置を抽出する抽出手段と、前記複数の周波数成分のうち、低周波成分を除いた、複数の高周波成分のそれぞれにおける、画素位置iに対応する前記解像度分解画像中の位置である対応位置と、当該対応位置の近傍の位置である近傍位置とに関して、前記抽出ステップにおいて抽出された各変換係数に応じた重みを、その変換係数に対応する前記焦点位置における画素位置iに関する確率値として加算する確率算出手段とを含んでいることが好ましい。また、請求項10又は11に記載の全焦点画像生成装置は、請求項12に記載の発明のように、前記最適焦点位置取得手段が、下記の数式2で表されるメッセージmti→jを、全ての隣り合うi,j及びxjに関してt=0,1,2,...,Tの順に、メッセージが収束するt=Tまで算出するメッセージ算出手段と、上記の数式3により前記最適な焦点位置xiを算出する焦点位置算出手段とを含んでいることが好ましい。また、第1の観点による請求項13に記載の全焦点画像生成プログラムは、焦点位置の異なる、L個L:2以上の自然数の被写体画像から、被写体の合焦点画像を生成するようにコンピュータを機能させるプログラムであって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度変換画像生成ステップと、前記多重解像度変換画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、前記撮像ステップにおいて撮像された前記焦点位置ごとの像から前記最適な焦点位置に対応する画素値を前記画素位置ごとに取得すると共に、取得した画素値に対応する画素を当該画素位置に配置して合焦点画像を生成する合焦画像合成ステップと、をコンピュータに実行させる。また、請求項13に記載の全焦点画像生成プログラムは、請求項14に記載の発明のように、前記最適焦点位置取得ステップにおいて、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置としてxiを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出するようにコンピュータを機能させてもよい。また、請求項14に記載の全焦点画像生成プログラムは、請求項15に記載の発明のように、前記焦点位置確率算出ステップが、前記多重解像度変換画像生成ステップにおいて生成された、前記焦点位置が異なる複数の前記解像度分解画像の変換係数同士を前記解像度分解画像中の画素ごとに比較し、比較した変換係数をその絶対値が最も大きいものから大きい順に所定の数だけ抽出すると共に、抽出した各変換係数に対応する前記解像度分解画像における前記焦点位置を抽出する抽出ステップと、前記複数の周波数成分のうち、低周波成分を除いた、複数の高周波成分のそれぞれにおける、画素位置iに対応する前記解像度分解画像中の位置である対応位置と、当該対応位置の近傍の位置である近傍位置とに関して、前記抽出ステップにおいて抽出された各変換係数に応じた重みを、その変換係数に対応する前記焦点位置における画素位置iに関する確率値として加算する確率算出ステップとを含んでいることが好ましい。また、請求項14又は15に記載の全焦点画像生成プログラムは、請求項16に記載の発明のように、前記最適焦点位置取得ステップが、上記の数式2で表されるメッセージmti→jを、全ての隣り合うi,j及びxjに関してt=0,1,2,...,Tの順に、メッセージが収束するt=Tまで算出するメッセージ算出ステップと、下記の数式3により前記最適な焦点位置xiを算出する焦点位置算出ステップとを含んでいることが好ましい。第2の目的を達成するためになされた本発明の第2の観点によると、請求項17に記載の被写体高さ情報取得方法に係る発明は、焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるL個L:2以上の自然数の被写体画像を撮像する撮像ステップと、前記撮像ステップにおいて撮像された像に基づいて、前記焦点位置ごとに前記被写体のグレースケール画像を取得するグレースケール画像取得ステップと、前記グレースケール画像取得ステップにおいて取得された前記グレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度変換画像生成ステップと、前記多重解像度変換画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得ステップと、前記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得ステップと、を備えている。請求項17に記載の被写体高さ情報取得方法は、被写体の撮像時に、各焦点位置と所定の深度方向に関するその焦点位置の位置情報とを関連付けておく。一方、第1の観点と同様に、最適な焦点位置を、大域最適化を用いて近似的に導出する。そして、導出した最適な焦点位置を焦点が合った位置とし、当該焦点位置から換算された高さ情報を画素位置ごとに取得する。最適な焦点位置を取得することは、画像中の大域的な領域に係る問題である。確率伝播法、平均場近似、グラフカットなどの大域最適化は、このような大域的な問題を最適に解決するための確率論的手法である。したがって、本発明によると、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って最適に取得された焦点位置に基づき、被写体の高さ情報を取得することができる。また、請求項17に記載の被写体高さ情報取得方法は、請求項18に記載の発明のように、前記最適焦点位置取得ステップにおいて、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置としてxiを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出してもよい。また、第2の観点による請求項19に記載の被写体高さ情報取得装置は、焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるL個L:2以上の自然数の被写体画像に基づいて、当該被写体の位置を取得する被写体高さ情報取得装置であって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度変換画像生成手段と、前記多重解像度変換画像生成手段が生成した前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出手段と、前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得手段と、前記最適焦点位置取得手段が取得した前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得手段と、を備えている。また、請求項19に記載の被写体高さ情報取得装置は、請求項20に記載の発明のように、前記最適焦点位置取得手段が、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置xiを割り当てた場合に、前記焦点位置確率算出手段が算出した前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出してもよい。また、第2の観点による請求項21に記載の被写体高さ情報取得プログラムは、焦点位置から換算された高さ情報と関連付けられた、焦点位置の異なるL個L:2以上の自然数の被写体画像に基づいて、当該被写体の高さ情報を取得するようにコンピュータを機能させるプログラムであって、被写体の前記焦点位置ごとのグレースケール画像に多重解像度変換を施すことにより、前記グレースケール画像が複数の周波数帯域の成分に分解された画像である解像度分解画像を、前記焦点位置ごとに生成する多重解像度変換画像生成ステップと、前記多重解像度変換画像生成ステップにおいて生成された前記解像度分解画像に基づいて、前記焦点位置に関する確率分布である焦点位置確率分布を、前記グレースケール画像の画素位置ごとに算出する焦点位置確率算出ステップと、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づいて、最適な前記焦点位置を、大域最適化を用いて前記画素位置ごとに取得する最適焦点位置取得ステップと、記最適焦点位置取得ステップにおいて取得された前記最適な焦点位置に基づいて、焦点位置から換算された高さ情報を、前記画素位置ごとに取得する高さ情報取得ステップと、をコンピュータに実行させる。また、請求項21に記載の被写体高さ情報取得プログラムは、請求項22に記載の発明のように、前記最適焦点位置取得ステップにおいて、i,jを前記画素位置を示す添え字とし、xi,xjをL以下の自然数を取る変数として、複数の画素位置iに関して前記焦点位置としてxiを割り当てた場合に、前記焦点位置確率算出ステップにおいて算出された前記焦点位置確率分布に基づく画素位置iごとのコスト関数Dixiと、xiとxjの差の絶対値が大きいほど大きい値を取るペナルティ関数Vxi,xjとを、前記複数の画素位置i及び隣り合うi,jに関して足し合わせた関数を評価関数Eとするとき、Eが最小になるような前記複数の画素位置iに関するxiの組み合わせを、確率伝播法を用いて、前記最適な焦点位置として近似的に導出するようにコンピュータを機能させてもよい。第1の観点による請求項1に記載の全焦点画像生成方法、請求項9に記載の全焦点画像生成装置、及び、請求項13に記載の全焦点画像生成プログラムによると、最適な焦点位置を、大域最適化を用いて近似的に導出する。そして、導出した最適な焦点位置を焦点が合った位置とし、当該焦点位置に対応する撮像画像から画素を取得して全焦点画像を生成する。最適な焦点位置を取得することは、画像中の大域的な領域に係る問題である。確率伝播法、平均場近似、グラフカットなどの大域最適化は、このような大域的な問題を最適に解決するための確率論的手法である。したがって、本発明によると、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って適切に焦点が合った画像を取得することができる。第2の観点による請求項17に記載の被写体高さ情報取得方法、請求項19に記載の被写体高さ情報取得装置、及び、請求項21に記載の被写体高さ情報取得プログラムによると、被写体の撮像時に、光軸方向の距離情報を、各焦点位置から換算し、高さ情報として、各焦点位置と関連付けておく。一方、第1の観点と同様に、最適な焦点位置を、大域最適化を用いて近似的に導出する。そして、導出した最適な焦点位置を焦点が合った位置とし、当該焦点位置に対応する高さ情報から、被写体の高さ情報を画素位置ごとに取得する。最適な焦点位置を取得することは、画像中の大域的な領域に係る問題である。確率伝播法、平均場近似、グラフカットなどの大域最適化は、このような大域的な問題を最適に解決するための確率論的手法である。したがって、本発明によると、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って最適に取得された焦点位置に基づき、被写体の高さ情報を適切に取得することができる。第1の実施形態本発明の全焦点画像生成方法、全焦点画像生成装置及び全焦点画像生成プログラムに係る一実施の形態である第1の実施形態について、図面に基づいて説明する。本実施形態の全焦点画像生成装置は、図1に示すように、サンプルSの拡大画像を撮像する顕微鏡1と顕微鏡1によって撮像された画像を処理する画像処理部100とを有している。顕微鏡1は、複数の対物レンズ2が固定されたレボルバ3と、サンプルSが設置されるステージ4と、ステージ4が上端に固定されたスライダ5と、スライダ5を鉛直方向に移動可能に支持するスライダ支持体7と、スライダ5をスライダ支持体7に対して鉛直方向に移動させるスライダ移動機構6とを有している。レボルバ3には、互いに異なる倍率の対物レンズ2が固定されている。レボルバ3を回転させることで、下方のステージ4側に向いた対物レンズ2を切り替えることができる。ステージ4の上面は水平であり、対物レンズ2の直下に水平に位置するようにサンプルSをステージ4上に置くことができる。スライダ移動機構6は、スライダ5を鉛直方向に微小に移動させる圧電アクチュエータと、鉛直方向に関するスライダ5の位置を高精度に検出する位置検出センサとを有している。位置検出センサは、スライダ5の位置を検出した結果を後述の制御部11へと出力する。また、顕微鏡1は、サンプルSへの照明光を射出する照明光学系8と、サンプルSからの光を撮像素子10に結像させる撮像光学系9と、結像した被写体像を光電変換して出力する撮像素子10とを有している。照明光学系8は、照明用の光源や光源からの照明光をサンプルSへと導く反射鏡等を有している。撮像光学系9は、対物レンズ2からの被写体像を撮像素子10へと導き、結像させるフォーカスレンズ9a図2参照などのレンズや反射鏡等を有している。被写体像は、対物レンズ2と撮像光学系9に設けられたレンズとの協働により、高倍率で拡大されて撮像素子10に結像される。撮像素子10は、撮像した被写体画像に対応する画像データを生成し、画像処理部100へと出力する。この画像データは、x方向、及び、x方向に直交するy方向に並んだ複数のカラー画素を示すデータからなる。水平方向に平行な所定の方向がx方向に対応し、水平方向に平行で上記所定の方向に直交する方向がy方向に対応する。この画像データにおいては、各画素の色が、光の三原色であるR赤、G緑及びB青のうち、一部の成分のみによって表されている。あるいは、各画素の色が、補色であるCyシアン,Mgマゼンタ、Yeイエロー及びGのうち、一部の成分のみによって表されていてもよい。顕微鏡1には、各部の構成を制御する制御部11が設けられている。制御部11は、スライダ移動機構6の駆動を制御し、ステージ4を鉛直方向に関してさまざまな位置に移動させる。ステージ4が移動すると、ステージ4上のサンプルSと対物レンズ2との距離が変動するため、ステージ4上に配置されたサンプルSからフォーカスレンズ9aまでの光の経路に沿った距離が変動する。一方、フォーカスレンズ9aの焦点距離が一定であり、フォーカスレンズ9aと撮像素子10との間の距離も一定であるとすると、被写体像が撮像素子10上にちょうど焦点が合うように結像するためには、被写体からフォーカスレンズ9aまでの光の経路に沿った距離が所定の大きさになる必要がある。以下、被写体において、フォーカスレンズ9aに対してちょうど焦点が合った位置を焦点位置と呼称する。サンプルSがフォーカスレンズ9aに対して焦点位置にあるとき、サンプルSの像は焦点が合った状態で撮像素子10に結像する。しかしながら、顕微鏡1で用いられる高倍率のレンズは、焦点深度が浅い。このため、サンプルSの表面に凹凸がある場合には、フォーカスレンズ9aまでの距離がサンプルSの表面の位置ごとに異なることになる。この距離の違いが焦点深度に対して無視できない大きさになると、サンプルSの表面上の特定の領域が焦点位置に配置されたとしても、特定の領域以外の領域は焦点位置からずれてしまい、画像がぼけてしまうことが生じ得る。例えば、図2に示すように、サンプルSの表面における領域s1と領域s2とでは、光の経路に沿った方向に関してフォーカスレンズ9aまでの距離D1及びD2が互いに異なる。したがって、この状態で領域s1は焦点が合っていたとしても、領域s2では焦点が合わないことがあり得る。被写体のあらゆる領域で焦点が合った画像を取得するためには、焦点位置が少しずつ異なった複数の画像から、焦点が合った領域の画素を集め、合成することが考えられる。そこで、制御部11は、スライダ移動機構6を制御し、焦点位置が少しずつ異なるようにサンプルSを配置させると共に、撮像素子10を制御し、各焦点位置でのサンプルSの被写体像を撮像させる。撮像された画像同士は、焦点位置が少しずつ異なるため、被写体において焦点が合っている領域が少しずつ異なったものとなる。制御部11は、撮像素子10が、焦点位置の異なる画像データを1つ、画像処理部100に出力するごとに、その焦点位置に対応するステージ4の位置を、スライダ移動機構6に設けられた位置検出センサの検出結果から取得し、画像処理部100へと出力する。これらの処理が繰り返されることにより、撮像素子10からは、焦点位置が異なるL個L:2以上の自然数の画像が出力されると共に、制御部11からは、L個の画像に関連付けられたL個のステージ位置情報が出力される。画像処理部100は、撮像素子10から入力されたL個の被写体画像に対応する画像データに基づいて、各画素に関して焦点が合った全焦点画像に対応する画像データを生成する。この処理のため、画像処理部100は、前段処理部110、グレースケール変換部101、多重解像度変換部102多重解像度変換画像生成手段、焦点位置確率算出部120焦点位置確率算出手段、信頼度設定部103、最適焦点位置取得部130最適焦点位置取得手段及び後段処理部140の各機能部を有している。画像処理部100は、CPUCentralProcessingUnit、ROMReadOnlyMemoryその他の各種のハードウェアと、ROMに格納されたプログラムその他のソフトウェアとから構成されている。このソフトウェアがハードウェアを、前段処理部110等の各機能部を実現するように機能させることにより、画像処理部100が実現されている。あるいは、各機能部における演算処理に特化した電子回路等により前段処理部110等が構成されていてもよいし、各機能部における演算処理に特化した電子回路等と、ハードウェアとソフトウェアの組み合わせによって各機能部を実現する構成とが互いに協働するように構成されていてもよい。画像処理部100に用いられるソフトウェアは、磁気的記録媒体や光学的記録媒体、USBメモリなどの各種記録媒体に格納されてもよい。そして、これらの記録媒体から読み出されたソフトウェアに従ってコンピュータが機能することにより、画像処理部100が実現されてもよい。画像処理部100が有する各機能部は、画像データを処理することにより、その画像データに対応する画像に対して各種の画像処理を施す。このように、各機能部が直接処理するのは画像データであるが、以下においては、冗長な表現を避けるため、各機能部が画像自身に対して直接、各種の画像処理を施すように表現する場合がある。このような表現は、その画像処理に対応する処理をその画像に対応する画像データに対して施すことを示すものとする。前段処理部110は、位置ずれ補正部111、撮像倍率補正部112及びRGB画像取得部113を有している。位置ずれ補正部111は、撮像素子10からのL個の被写体画像同士を比較し、水平方向に関する被写体の位置ずれを検出する。この位置ずれは、画像中のx方向及びy方向に関する位置ずれとなって表れる図3の画像IM1参照。位置ずれ補正部111は、被写体画像中の被写体の位置をずらすことにより、焦点位置が異なる被写体画像同士で被写体の位置が一致するように被写体画像を補正する。撮像倍率補正部112は、制御部11から入力されるステージ4の位置情報に基づいて、被写体画像の倍率を補正する。撮像光学系9の構成によっては、被写体の焦点位置が変動すると、撮像倍率が変動することがある。撮像倍率補正部112は、ステージ4の位置に関連付けて、焦点位置の変動に伴う撮像倍率の変動分に対応する拡大縮小率を保持している。これにより、撮像倍率補正部112は、制御部11から入力されるステージ4の位置に対応する拡大縮小率を導出する。そして、撮像倍率補正部112は、導出した拡大縮小率に従って、被写体画像を拡大又は縮小する。撮像倍率補正部112は、これにより、焦点位置が異なるL個の被写体画像同士で倍率が一致するように各被写体画像の倍率を補正する。RGB画像取得部113は、図3に示すように、位置ずれ補正部111及び撮像倍率補正部112による補正が施されたL個の被写体画像IM1に対応する画像データから、各画素においてR,G及びBの各色の成分値が全て決定されたL個の画像IM2に対応するRGB画像データIM2を生成する。画像IM1に対応する画像データには、上述のとおり、画素ごとにR,G及びBのいずれか一部の値のみが含まれている。したがって、RGB画像取得部113は、各画素において、足りない成分の値を、その画素の周囲に位置する画素の値に基づいて補間する。例えば、RGB画像取得部113は、補間の対象となる画素の周囲の画素値に重みを乗算したものを足し合わせることにより、補間した画素値を算出する。これにより、RGB画像取得部113は、焦点位置が異なるL個の被写体画像のそれぞれに関して、画像IM2に対応するRGB画像データを生成する。グレースケール変換部101は、図3に示すように、RGB画像取得部113が生成したL個の画像IM2に対応するRGB画像データを、L個のグレースケール画像IM3に対応するグレースケール画像データに変換する。変換方法はさまざまなものを採用できる。一例として、NTSCNationalTelevisionSystemCommittee方式によるRGB信号から輝度信号への変換方法を採用する場合、グレースケール変換部101は、以下の式を使用して、グレースケール画像の各画素値を算出する:Y=0.298912R0.586611G0.114478B。なお、Yは、グレースケール画像の各画素値であり、R,G及びBは、RGB画像の各画素におけるRの成分値,Gの成分値及びBの成分値である。多重解像度変換部102は、グレースケール変換部101が生成したグレースケール画像データに対応する画像に多重解像度変換を施す。この多重解像度変換には、ガウシアンピラミッド変換、ラプラシアンピラミッド変換などのさまざまな方法を採用できる。本実施形態では、離散ウェーブレット変換を採用する。離散ウェーブレット変換としては、HaarWabelet、DaubechiesWavelet,Cohen-Daubechies-FeauveauWaveletなどのうち、いずれかを用いた変換方法を使用する。離散ウェーブレット変換は、画像を高周波成分と低周波成分とに分解すると共に、分解した低周波成分をさらに高周波成分と低周波成分とに分解するという処理を繰り返すことで行われる。これにより、高周波成分及び低周波成分に分解された画像である複数の解像度分解画像が取得される。高周波成分及び低周波成分への分解には、各種のフィルタ処理が用いられる。このフィルタ処理に用いるフィルタは、画像のx方向に関するハイパスフィルタ及びローパスフィルタと、画像のy方向に関するハイパスフィルタ及びローパスフィルタとを含んでいる。以下、これらの各方向に関するフィルタを、xハイパスフィルタ、xローパスフィルタ、yハイパスフィルタ及びyローパスフィルタと呼称する。フィルタ処理を画像に施す具体的な方法は以下のとおりである。多重解像度変換部102は、複数のフィルタ行列を示すデータをフィルタごとに有している。例えば、多重解像度変換部102は、xハイパスフィルタに対応するフィルタ行列を示すデータや、yローパスフィルタに対応するフィルタ行列を示すデータを有している。フィルタ行列は複数のフィルタ係数から構成されており、各フィルタ係数は、フィルタが施される領域に含まれる画素ごとに設定されている。多重解像度変換部102は、フィルタ行列に含まれる各フィルタ係数を、その係数に対応する位置の画素値に乗算する。そして、多重解像度変換部102は、乗算した値を足し合わせ、その足し合わせた和を、当該フィルタ処理後の画素値とする。また、多重解像度変換部102は、フィルタ処理後の画像全体における画素数が、x方向及びy方向のそれぞれに関して、フィルタ処理前の画像全体における画素数の半分になるように各フィルタ処理を施す。多重解像度変換部102は、グレースケール画像変換部が生成した画像データに対して各種のフィルタ処理を施すことで、以下のように1回目のウェーブレット変換を施す。まず、多重解像度変換部102は、xハイパスフィルタを用いてフィルタ処理を施した後に、yローパスフィルタを用いてフィルタ処理を施す。これにより、グレースケール画像において、x方向に関しては高周波であり、y方向に関しては低周波である成分が取得される。以下、この成分をHL成分と呼称する。同様に、多重解像度変換部102は、xハイパスフィルタを用いてフィルタ処理を施した後に、yハイパスフィルタを用いてフィルタ処理を施す。これにより、x方向及びy方向のいずれに関しても高周波である成分が取得される。以下、この成分をHH成分と呼称する。また、多重解像度変換部102は、xローパスフィルタを用いてフィルタ処理を施した後に、yハイパスフィルタを用いてフィルタ処理を施す。これにより、x方向に関しては低周波であり、y方向に関しては高周波である成分が取得される。以下、この成分をLH成分と呼称する。さらに、多重解像度変換部102は、xローパスフィルタを用いてフィルタ処理を施した後に、yローパスフィルタを用いてフィルタ処理を施す。これにより、x方向及びy方向のいずれに関しても低周波である成分が取得される。以下、この成分をLL成分と呼称する。このように、ウェーブレット変換によって、グレースケール画像から、HH成分、HL成分、LH成分及びLL成分からなる、複数の周波数帯域の成分を有する解像度分解画像が生成される。解像度分解画像の各成分に含まれる、ウェーブレット変換によって導出された1つ1つの値を、ウェーブレット係数ウェーブレット変換によって生成された変換係数という。図4aは、1回目のウェーブレット変換により生成されたHH成分、HL成分、LH成分及びLL成分を、x方向及びy方向に対応するv方向及びw方向に、2行2列に並べたものである。これらの成分は、HH成分が右下隅に、HL成分が右上隅に、LH成分が左下隅に、LL成分が左上隅になるように配列されている。多重解像度変換部102は、1回目のウェーブレット変換によって生成されたLL成分に対し、2回目のウェーブレット変換を同様に施す。これにより、再びHH成分、HL成分、LH成分及びLL成分が得られる。このように、前回のウェーブレット変換によって得られたLL成分に対してさらにウェーブレット変換を施すことの繰り返しにより、HH成分、HL成分、LH成分及びLL成分の各成分がウェーブレット変換を施すたびに得られる。以下、ウェーブレット変換n:自然数をn回施すことをnレベルと呼称する。これらの各画像を、ウェーブレット変換ごとに、前回のウェーブレット変換におけるLL成分の位置に図4aに示す配列で配置したものが、図4bの画像IM4である。解像度分解画像における各成分をこのように配置したものを多重解像度表現という。画像IM4は、一例として、3レベルの多重解像度表現を示している。この多重解像度表現は、図4bに示すように、各ウェーブレット変換により得られるHH成分、HL成分及びLH成分と、最後のウェーブレット変換により得られるLL成分とからなる。各レベルのHH成分、HL成分及びLH成分を高周波成分といい、最終的に生成されたLL成分を低周波成分という。図4bのL1は1レベルの、L2は2レベルの、L3は3レベルの高周波成分を示している。多重解像度表現は、v方向に関する幅がグレースケール画像のx方向に関する幅と等しく、w方向に関する幅がグレースケール画像のy方向に関する幅と等しい。グレースケール画像の各画素値は、多重解像度表現において周波数帯域が異なる複数のウェーブレット係数へと分解される。このため、グレースケール画像における画素位置x,yと、多重解像度表現における位置v,wとは、一対一には対応せず、一体多の対応関係となる。例えば、図5に示すように、グレースケール画像IM3における画素P1は、多重解像度表現IM4における異なる周波数帯域のそれぞれにおいてP1に対応する位置であるQ11~Q34のウェーブレット係数と対応することとなる。このようにして、多重解像度変換部102は、焦点位置が異なるL個のグレースケール画像データのそれぞれから多重解像度表現に対応するデータを生成する。こうして得られたL個の多重解像度表現を焦点位置ごとに並べたものが、図6aの多重解像度表現IM4_1~IM4_Lである。多重解像度変換部102は、v方向及びw方向に関する位置ごとに、多重解像度表現IM4_1~IM4_Lの全てのウェーブレット係数を順に算出する。つまり、v方向及びw方向に関するある位置について、多重解像度表現IM4_1~IM4_Lの全てのウェーブレット係数を順に算出した後に、v方向及びw方向に関する次の位置について、多重解像度表現IM4_1~IM4_Lの全てのウェーブレット係数を順に算出する。例えば、図6aの二点鎖線R2で囲まれた各点は、v方向及びw方向に関する位置が同じである。また、二点鎖線R3で囲まれた各点は、v方向及びw方向に関する位置が同じであり、R2で囲まれた点とは位置が異なる。このとき、多重解像度変換部102は、例えば、R2で囲まれた点の位置に関して多重解像度表現IM4_1~IM4_Lの全てのウェーブレット係数を算出してから、R3で囲まれた点の位置に関して多重解像度表現IM4_1~IM4_Lの全てのウェーブレット係数を算出する。多重解像度変換部102がこのような順序でウェーブレット係数を算出する理由については、後述する。焦点位置確率算出部120は、多重解像度変換部102が生成した解像度分解画像から、グレースケール画像における画素位置ごとに、焦点位置に関する確率分布焦点位置確率分布を算出する。ここで、焦点位置に関する確率値とは、次のようなものとする。以下、ii:グレースケール画像の総画素数以下の自然数をグレースケール画像の画素位置を区別する添え字とし、xiをL以下の自然数を取る変数とする。このとき、焦点位置確率算出部120が算出する確率分布をPixiと表す。Pixiは、画素位置iにおいて、焦点が合っているものが、L個の画像のうち、xi番目の画像である確からしさを示す。言い換えれば、Pixiは、画素位置iにおいて被写体の表面が焦点位置に配置されているものが、L個の画像のうち、xi番目の画像であることの確からしさに相当する。このような確率分布を導出するため、焦点位置確率算出部120は、画像候補決定部121抽出手段及び確率分布算出部122確率算出手段を有している。画像候補決定部121は、多重解像度表現におけるv方向及びw方向に関する位置ごとに、L個の多重解像度表現から、確率分布を算出するための候補となるM個M:L未満の自然数の多重解像度表現を抽出する。候補の数は、パラメータチューニングによって適宜設定される。このM個の多重解像度表現は、v方向及びw方向に関する当該位置に関して、L個の画像から、ウェーブレット係数の絶対値が大きい上位M個を抽出することにより取得される。具体的には、画像候補決定部121は、当該位置におけるウェーブレット係数の絶対値を、焦点位置が異なる画像同士で比較する。そして、画像候補決定部121は、ウェーブレット係数の絶対値を、最も大きいものからM番目に大きいものまで抽出すると共に、抽出した値に対応するM個の多重解像度表現を抽出する。例えば、図6aの多重解像度表現IM4_1~IM_4において、二点鎖線R1、R2及びR3で囲まれた各位置のウェーブレット係数を、R1~R3のそれぞれにおいて比較する。これにより、ウェーブレット係数の絶対値が大きいものからM個を抽出した結果と、M個のウェーブレット係数に対応する多重解像度表現とが、一例として、図6bに表されている。本実施形態は、多重解像度変換部102と画像候補決定部121とが以下のように各自の処理を協働して実行することにより、ウェーブレット係数を格納するためのメモリとして、候補となる値を格納するのに必要な容量だけを用意すればよいように構成されている。多重解像度変換部102は、上述のとおり、v方向及びw方向に関して同じ位置のウェーブレット係数を、全ての多重解像度表現IM4_1~IM4_Lに関して順に算出する。一方、画像候補決定部121は、多重解像度変換部102が算出した1個目~M個目のウェーブレット係数を、候補格納用のメモリに順に格納していく。次に、画像候補決定部121は、多重解像度変換部102がM1個目のウェーブレット係数を算出した段階で、このM1個目のウェーブレット係数を既にメモリに格納したM個のウェーブレット係数と比較する。そして、画像候補決定部121は、メモリ内のM個の値のうち、M1個目の値より小さいものがあれば、そのメモリの値をM1個目の値に変更する。これを、多重解像度変換部102がL番目の多重解像度表現に関するウェーブレット係数を算出するまで繰り返すことで、画像候補決定部121は、v方向及びw方向に関する当該位置について、M個のウェーブレット係数の候補を決定する。このように処理することにより、画像候補を決定するために必要なメモリ容量が、v方向及びw方向に関する位置ごとにM個の値を格納可能な大きさに限定されている。確率分布算出部122は、画像候補決定部121が抽出したウェーブレット係数に基づいて、確率分布Pixiを算出する。確率分布算出部122は、画素位置iごとに、多重解像度表現における当該画素位置iに対応する位置である対応位置と、その対応位置の近傍の位置である近傍位置とにおけるM個の候補値を使用する。なお、M個の候補値は、各多重解像度表現に含まれる高周波成分に関するもののみが使用される。例えば、図5のグレースケール画像IM3における点P1が注目画素であるとする。多重解像度表現IM4において点P1に対応する高周波成分における対応位置は、Q11~Q33である。また、対応位置Q11~Q33の近傍位置としては、一例として、各対応位置にv方向及びw方向のそれぞれに関して隣接する4つの位置と、斜め方向に隣接する4つの位置とからなる8つの位置が用いられる。例えば、Q13に関しては、図5に示すように、対応位置であるQ13自身とその近傍位置とからなる縦3つ横3つからなる合計9個の位置における候補値が用いられる。したがって、Q11~Q33の全体としては、合計で81個の対応位置及び近傍位置における候補値が用いられる。確率分布算出部122は、対応位置又は近傍位置に対応するグレースケール画像の画素位置をiとし、その候補値に対応する多重解像度表現をIM_xiとしたとき、対応位置又は近傍位置の候補値を、確率値Pixiに投票する。例えば、図7に示すように、ある対応位置の第1の候補値、第2の候補値、...、第Mの候補値が、多重解像度IM4_2、IM_4、...、IM_nn:L以下の自然数に対応する候補値であったとする。このとき、この対応位置の第1の候補値、第2の候補値、...、第Mの候補値は、Pi2、Pi4、...、Pinに対して投票される。投票は、具体的には、以下のような演算によって行われる。まず、Pi1、Pi2、...、PiLは、全てのiに関して、初期値0が与えられる。そして、投票する側の位置である対応位置又は近傍位置を投票位置v,wとし、投票される側の位置である対応位置を投票中心vc,wcとし、対応位置又は近傍位置における候補値をνvwとするとき、候補値νvwに応じた投票値は、下記の数式6のように算出される。パラメータλ、θは、状況に応じて設定されるが、一例として、λ=0.1、θ=1.0などと設定される。minα,βは、αとβのうち、小さいものを表す。確率分布算出部122は、各候補値に関して数式10によって算出された投票値を、当該画素位置iにおける、候補値に対応する多重解像度表現IM_xiに関する確率値Pixiに加算する。なお、図5には、Q13が対応位置であるときの投票中心と、投票位置の一例とが示されている。確率分布算出部122は、画素位置iに関する全ての候補値について投票を完了すると、Pi1、Pi2、...、PiL全体のノルムが1になるようにPixiを正規化する。確率分布算出部122は、全てのiに関して投票を実行する。このように、対応位置のみならず、その近傍の近傍位置からも投票がなされると共に、投票中心との距離に応じた重みをつけて投票がなされる。このため、確率分布の算出の際、局所的な整合性を考慮することができる。以上のように、焦点位置確率算出部120が、グレースケール画像の各画素位置iにおける確率分布Pixiを算出する。Pixiは、上記の投票により求められることから、L個の画像のうち、画素位置iの周辺の投票値が比較的大きい画像において、大きい値を有することとなる。投票値が大きいことは、ウェーブレット係数の高周波成分において振幅が大きいことに対応する。したがって、通常、Pixiが最大値を取るxiが、図8aに示すように、焦点が合っている最も確からしい画像、つまり、焦点位置に対応する画像を示すと考えられる。ところが、撮像時に画像ノイズ、飽和、ぼけなどの問題が生じた場合、上述のように算出したPixiにおいて最大値を取るxiが、図8bに示すように、実際の焦点位置から大きく離れた画像を示す場合がある。この場合、Pixiの値をそのまま使用して焦点位置を見積もると、焦点が合っていない位置を焦点が合っているものと誤って見積もってしまうおそれがある。そこで、信頼度設定部103は、Pi1、Pi2、...、PiLのうちの最大値と所定の閾値との大小を比較する。そして、最大値が所定の閾値以下となる場合に、最大値が所定の閾値より大きい場合と比べて、その画素位置の信頼度を低く設定する。具体的には、信頼度設定部103は、Pixiの画素位置iにおける最大値PiMAXが所定の閾値εを超えた場合図8aの場合、画素位置iの信頼度RiをRi=1PiMAX-εと設定する。また、図8aに示すように、Pixiの値は変更しない。一方、Pixiの画素位置iにおける最大値PiMAXが所定の閾値ε以下であった場合図8bの場合、信頼度設定部103は、画素位置iが問題領域内であると判定し、Ri=1と設定する。さらに、Pi1、Pi2、...、PiLの値は、図8bに示すように、全て1/Lになるような均一分布に変更して設定する。最適焦点位置取得部130は、焦点位置確率算出部120が算出したPixi及び信頼度設定部103が設定したRiに基づいて、画素位置iごとに、焦点位置が異なるL個の画像のうち、いずれの画像が焦点の合っている最適な焦点位置に対応する画像かを取得する。最適な焦点位置に対応する画像を取得する際の問題は、画像における大域的な領域に亘る整合性をいかに保証するかにある。そこで、本実施形態では、かかる問題解決の手段として、以下のように、マルコフ確率場による定式化を採用する。最適な焦点位置に対応する画像を決めることは、各画素位置iに関して、1番目~L番目のいずれの画像が焦点位置に対応する画像かを決めることに相当する。つまり、1~Lの任意の自然数を取るxiに対して、最適な焦点位置に対応する画像を与えるxi=xiをどのように割り当てるかということに相当する。マルコフ確率場によると、ある割り当て{xi}を与える全体の確率が、各画素位置iに隣接する画素位置jにおける割り当てと画素位置iにおける割り当てとの関係のみから決定される。本実施形態においては、図9に示すように、x方向及びy方向のそれぞれに関して各画素位置iに隣接する画素位置j1~j4の4つにおける割り当てと画素位置iにおける割り当てとの関係のみから決定されるものとする。このようなマルコフ確率場において、最適な割り当て{xi}を与える問題は、下記の数式7で表される評価関数EXを最小化するXを決める問題に帰着することが知られている。EXは、割り当てX=x1,x2,...,xNに関する関数であることを示している。Nは、グレースケール画像の総画素数である。ここで、Aは、互いに隣接する一対の画素からなる集合を示している。つまり、jは、L以下の自然数であって、iに隣接する画素位置を示す添え字である。また、Dixiはコスト関数であり、画素位置iに関してxiを割り当てるときのコストを示す。本実施形態において、Dixiは、以下の数式8によって与えられる。このように、Dixiは、信頼度設定部103が設定した信頼度Riと、焦点位置確率算出部120が算出した確率値Pixiとの積で表される。Vxi,xjはペナルティ関数であり、互いに隣接する一対の画素位置i、jに、それぞれxi、xjが割り当てられたときのペナルティを示す。なお、Vxi,xjは、通常、互いに隣接する画素の割り当てが異なるほど大きいペナルティを与える関数であることから、平滑化項とも呼ばれる。本実施形態において、Vxi,xjは、以下の数式9によって与えられる。このように、Vxi,xjは、割り当てられたxi,xj同士の差の絶対値に比例する値を有している。EXを最小化する問題の解法には、確率伝播法BeliefPropagationが用いられる。確率伝播法によって、多くの場合に大域的な最適解に近い解が近似的に得られることが知られており、例えば、ステレオ視のための視差計算や、いわゆる画像領域分割の解法として用いられている。確率伝播法に基づいてEXの最小化問題を解くため、最適焦点位置取得部130は、階層化演算部131、メッセージ算出部132メッセージ算出手段及び最適焦点位置算出部133焦点位置算出手段を有している。まず、階層化演算部131は、上記の数式8に基づき、コスト関数Dixiを全ての画素位置iに関して算出する。次に、コスト関数Dixiをダウンサンプリングする。例えば、x方向に関して1/2、y方向に関して1/2の画素数のコスト関数Di1xi1を算出する。ダウンサンプリングは、x方向及びy方向に関して1つおきにコスト関数Dixiから値を取る間引き処理により実行してもよい。あるいは、注目画素とその周辺画素とにおける複数のDixiの値の和や加重平均を取ることを、x方向及びy方向に関して1つおきに実施してもよい。さらに、Di1xi1を同様にダウンサンプリングすることにより、Di2xi2を取得する。このように、階層化演算部131は、ダウンサンプリングをn回n:自然数実行することにより、複数階層の低解像度のコスト関数Di1xi1、Di2xi2、...、Dinxinを取得する。以下、ダウンサンプリングしていない元のコスト関数を、低解像度のコスト関数の表記と合わせて、Di0xi0と表す。メッセージ算出部132は、階層化演算部131が算出したコスト関数Di0xi0、Di1xi1、Di2xi2、...、Dinxinに基づいて、下記の数式10で表されるメッセージmti→jkを繰り返し計算する。mti→jk[xjk]は、t回目の繰り返しにおいて、画素位置iから隣接する画素位置jへの伝播メッセージ値を示す。xjkは1~Lのいずれかを取る変数であるため、mti→jkもL個の要素を持つ。Nei\jは、画素位置iに隣接する画素位置であって画素位置j以外のものからなる集合を意味する。kは、階層を示す0以上且つn以下の整数である。メッセージ算出部132は、まず、k=nのときのメッセージを、t=0,1,2,...の順に計算する。t=0のときの初期値m0i→jnは0ゼロに設定される。ここで、minαの下にxikが記載された演算子は、xikに関するαの最小値を与える。k=nのときのメッセージは、以下のように計算される。なお、tt≧1のメッセージを算出するときには、t-1以下のメッセージが全て計算されているとする。まず、下記の数式11に示すように、あるxjnに関して、xin=1、2、...、Lを順に代入することで、数式10におけるminのカッコ内を算出する。次に、数式11に基づいて算出した値のうち、最小値を取るものを抽出する。これが、当該xjnに関するmti→jn[xjn]の値となる。これを、xjn=1、2、...、Lに関して順に繰り返すことで、全てのxjnに関してmti→jn[xjn]が算出されることとなる。このような作業を、全ての隣接するi,jについて繰り返すことで、あるtにおける全てのメッセージmti→jn[xjn]を算出することができる。メッセージ算出部132は、このようなメッセージmti→jnの算出を、t=0,1,2,...の順に、mti→jnが収束するt=Tまで繰り返す。メッセージ算出部132は、次に、k=nより1階層高解像度のk=n―1のときのメッセージmti→jn-1を上記と同様に算出する。ここで、メッセージの初期値には、低解像度k=nにおいて算出されたメッセージの収束値mTi→jnが用いられる。収束値mTi→jnは、低解像度のときの画素位置i→画素位置jの経路に対応する高解像度k=n-1の全ての経路におけるメッセージの初期値として用いられる。例えば、図10に示すように、低解像度における1画素に対応する範囲が、高解像度における縦2×横2の4つの画素の範囲に対応したとする。このとき、低解像度における破線の矢印に示す経路に対応するメッセージの収束値mTi→jnが、この経路に対応する高解像度における4本の実線の矢印に示す経路に対応するメッセージの初期値として共通に使用される。そして、メッセージ算出部132は、k=n―1としたときのメッセージmti→jn-1を、t=0,1、2、...の順に、収束するまで繰り返す。この収束値は、上記と同様、さらに1階層高解像度k=n-2のメッセージの初期値として用いられる。以上のように、メッセージ算出部132は、k=n,n-1,...,2,1,0の順にメッセージmti→jkの算出を繰り返す。そして、最適焦点位置算出部133は、最終的に得られたメッセージの収束値mTi→j0を用いて、最適な焦点位置に対応する画像を示すxiを、下記の数式12に基づいて算出する。後段処理部140は、最適焦点位置取得部130が取得したxiに基づいて、全画素において焦点が合った全焦点画像を取得する。後段処理部140は、ノイズ除去部141及び全焦点画像取得部142合焦画像合成手段を有している。ノイズ除去部141は、最適焦点位置取得部130が取得したxiに対して、さらに、ガウシアンフィルタやメディアンフィルタを用いてノイズ除去処理を施す。xiは上記のとおり、大域的な最適化が可能な確率伝播法に基づいて取得されるため、低周波数の推定誤差が生じにくい。したがって、ここで使用されるフィルタとしては、比較的少ないフィルタカーネルのもので十分である。全焦点画像取得部142は、ノイズ除去部141がノイズを除去したxiに基づき、全焦点画像を取得する。xiは、画素位置iにおいて、L個の画像のうち、いずれの画像が焦点の合っているものに相当するかを示す。これに従い、全焦点画像取得部142は、図3の画像IM2のうち、各画素位置iにおいて、xi番目の画像に対応する画像の画素値を取得する。そして、全焦点画像取得部142は、画素位置iごとに取得した画素を対応する位置に配置し直すことにより、全焦点画像を取得する。以下、顕微鏡1及び画像処理部100が実行する処理の全体の流れについて、図11を参照しつつ説明する。図11の各ステップに従って実行される処理が、本発明の全焦点画像生成方法に対応する。また、これらの処理は、上記のとおり、画像処理部100を構成するハードウェアとソフトウェアとが協働して実行してもよい。この場合のソフトウェアが、本発明の全焦点画像生成プログラムに対応する。まず、顕微鏡1において、ステージ4の位置を移動させることにより、焦点位置が異なるL個の被写体画像の撮影が行われるステップS1。これらの被写体画像は、各焦点位置に対応したステージ4の位置情報と共に、画像処理部100へと送られる。次に、画像処理部100において、位置ずれ補正部111が、被写体の水平方向の位置ずれに対応する画像中のx方向及びy方向に関する位置ずれを補正するステップS2。次に、撮像倍率補正部112が、被写体画像ごとに顕微鏡1から入力されるステージ4の位置に基づいて被写体画像の倍率を補正するステップS3。次に、RGB画像取得部113が、位置ずれ補正及び倍率補正が施されたL個の画像IM1に基づいて、各画素においてR,G及びBの各色の成分値が全て決定されたL個の画像IM2を生成するステップS4。次に、グレースケール変換部101が、L個の画像IM2から、L個のグレースケール画像IM3を生成するステップS5グレースケール画像取得ステップ図3参照。次に、多重解像度変換部102が、グレースケール画像IM3にウェーブレット変換を施すステップS6。このとき、多重解像度変換部102は、v方向及びw方向に関する位置ごとに、多重解像度表現IM4_1~IM4_LのL個のウェーブレット係数を順に算出する図6参照。一方、画像候補決定部121は、多重解像度変換部102が、最初のM個を算出するまで、算出されたウェーブレット係数を順にメモリに格納するステップS7。次に、多重解像度変換部102が、M1個目以降のウェーブレット係数を算出するステップS8ごとに、画像候補決定部121が、既にメモリに格納したM個の値のうち、新たに算出されたウェーブレット係数より小さいものがないか判定するステップS9。そして、小さいものがあればステップS9、Yes、画像候補決定部121は、そのメモリの値を新たに算出されたウェーブレット係数に更新するステップS10。一方、小さいものがなければステップS9、No、画像候補決定部121は、当該位置に関して既にL個のウェーブレット係数が算出されたか否かを判定するステップS11。まだL個に達していなければステップS11、No、ステップS8の処理に戻る。一方、L個に達していればステップS11、Yes、画像候補決定部121は、当該位置でのウェーブレット係数の候補を、この時点でメモリに格納されているM個の値に決定するステップS12。なお、ステップS6及びS8が、本発明における多重解像度画像生成ステップに対応する。また、ステップS7及びS9~S12が抽出ステップに対応する。次に、画像候補決定部121は、ウェーブレット係数を算出するべき位置が残っているか否かを判定するステップS13。まだウェーブレット係数を算出するべき位置が残っていると判定された場合にはステップS13、Yes、ステップS6の処理に戻る。v方向及びw方向に関する全ての位置についてウェーブレット係数が算出されたと判定された場合にはステップS13、No、ステップS14の処理に移る。次に、確率分布算出部122が、画像候補決定部121が決定した候補値に基づいて、各画素位置iにおける確率分布Pixiを算出するステップS14確率算出ステップ。なお、ステップS7、S9~S12及びS14が、本発明における焦点位置確率算出ステップに対応する。次に、信頼度設定部103が、Pixiに基づいて画素位置iごとに信頼度Riを設定するステップS15信頼度設定ステップ。次に、階層化演算部131が、確率分布算出部122が算出した確率分布Pixiと信頼度設定部103が設定した信頼度Riとに基づいて、コスト関数Dixiを算出すると共に、低解像度のコスト関数Dikxikk=1,2,...,nを取得するステップS16低解像度化ステップ。次に、メッセージ算出部132が、k=n、n-1、...、1、0の順にメッセージmti→jk[xj]を算出するステップS17。メッセージ算出部132は、メッセージが収束するt=Tまで繰り返しメッセージを算出していくステップS18、No→ステップS17。メッセージが収束するとステップS18、Yes、メッセージ算出部132は、次の階層があればステップS19、Yes、次の階層に関してメッセージを算出するステップS17。次の階層がなければステップS19、No、最適焦点位置算出部133が、最終的に収束したメッセージmTi→j0を用いて、最適な焦点位置に対応する画像を示すxiを算出するステップS20焦点位置算出ステップ。なお、ステップS16~S20が、本発明における最適焦点位置取得ステップに対応する。また、ステップS17及びS18が、本発明におけるメッセージ算出ステップに対応する。次に、ノイズ除去部141が、xiからノイズを除去しステップS21、その後、全焦点画像取得部142が、ノイズ除去部141がノイズを除去したxiに基づいて、画素位置iごとに焦点が合った全焦点画像を取得するステップS22合焦点画像合成ステップ。以上説明した本実施形態によると、評価関数EXを最小化する問題を確率伝播法に基づいて解くことにより、最適な焦点位置に対応する画像を示すxiを画素位置iごとに取得する。評価関数EXを最小化することは画像中の大域的な領域に係る問題である。また、確率伝播法は、このような大域的な問題を最適に解決するための確率論的手法である。したがって、本実施形態によると、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って適切に焦点が合った画像を取得することができる。また、本実施形態においては、確率伝播法におけるメッセージの算出の際、低解像度でメッセージを算出した結果を初期値として、より高解像度でメッセージを算出する。このため、高解像度での算出の際、より最適解に近い値を初期値としてメッセージが更新されるので、さらに、局所的な解に陥りにくくなる。第2の実施形態以下、本発明の被写体高さ情報取得方法、被写体高さ情報取得装置及び被写体高さ情報取得プログラムに係る第2の実施形態について説明する。第2の実施形態の被写体高さ情報取得装置において、第1の実施形態の全焦点画像生成装置と異なるのは、後段処理部100の代わりに図12の後段処理部240が設けられることのみである。また、後段処理部240において、後段処理部100と異なるのは、焦点画像取得部142の代わりに被写体高さ情報取得部242が設けられることのみである。第2の実施形態の被写体高さ情報取得装置は、図11のステップS1~S21までの処理を実行した後、被写体高さ情報取得部242が、以下のとおり、xiに基づいて、被写体の位置を取得する。xiが表す被写体画像は、フォーカスレンズ9aに対し、被写体の表面において画素位置iに対応する領域がちょうど焦点が合うときに撮像された画像に対応する。そこで、被写体高さ情報取得部242は、顕微鏡1の制御部11から送られたステージ4の位置情報に基づき、xiが表す画像が撮影されたときのステージ4の位置を取得する。このステージ4の位置から、被写体の各領域における鉛直方向に関する相対的な位置が取得できる。例えば、図2において、フォーカスレンズ9aに対し、光の経路に沿った方向に関してちょうど距離D1のときに焦点が合うとする。つまり、図2の状態のとき、サンプルSの表面において、領域s1に焦点が合っている。一方で、領域s2に焦点が合うためには、領域s2がフォーカスレンズ9aに対し距離D1になるところまでサンプルSが移動しなければならない。そこで、領域s1に焦点が合っている状態におけるステージ4の位置と、領域s2に焦点が合っている状態におけるステージ4の位置との差から、領域s1の位置と領域s2の位置との差D1-D2を取得することができる。被写体高さ情報取得部242は、このようにして、各画素位置iに対応する被写体の表面上の各領域の相対的な高さ情報を取得することにより、被写体における表面全体の相対的な高さ情報を取得する。第2の実施形態においても、第1の実施形態と同様にxiを取得することから、画像ノイズ、飽和、ぼけなどのさまざまな撮像時の問題がある場合にも、画像中の大域的な領域に亘って適切に被写体の位置を取得することができる。なお、第2の実施形態に係る被写体高さ情報取得装置が図11のステップS1~S21を実行すると共に、被写体高さ情報取得部242がxiに基づいて被写体の位置を取得する処理が、本発明の被写体高さ情報取得方法に対応する。また、このような処理を実行するためのプログラムが、本発明の被写体高さ情報取得プログラムに対応する。なお、第2の実施形態に使用される各パラメータは、第1の実施形態に使用されるパラメータと異なっていてもよい。変形例以上は、本発明の一実施の形態についての説明であるが、本発明は上述の実施の形態に限られるものではなく、課題を解決するための手段に記載された範囲の限りにおいて、以下に例を挙げるとおり、様々な変更が可能なものである。[1]上述の実施の形態では、グレースケール変換部101が、RGB画像データからグレースケール画像データを生成する。しかし、撮像素子10が直接グレースケール画像を生成してもよい。また、撮像素子10が生成したG画素のみからなる画像をグレースケール画像として使用してもよい。[2]上述の実施の形態では、画像候補決定部121が、確率分布を算出するための候補となるM個M<Lの多重解像度表現を抽出する。しかし、L個全ての多重解像度表現から確率分布を算出してもよい。ただし、本発明者の知見によると、全ての多重解像度表現を使用する場合、画質が悪化するおそれがあるため、M個の候補を抽出したほうが好ましい。[3]上述の実施の形態においては、確率分布を算出するため、対応位置及び近傍位置のそれぞれにおけるM個の候補値から所定の重みを用いて算出した投票値を投票している。しかし、近傍位置から投票せず、対応位置のみから投票してもよい。[4]上述の実施の形態においては、コスト関数Dixiをダウンサンプリングした低解像度のコスト関数に基づいてメッセージを算出する。しかし、コスト関数をダウンサンプリングせず、そのまま用いてメッセージ関数を算出してもよい。また、コスト関数として、Pixiに信頼度Riを乗算したものでなく、Pixiそのものを使用してもよい。[5]上述の実施の形態においては、顕微鏡1において撮影された被写体画像を用いる場合に本発明が適用されている。しかし、デジタルカメラや、デジタルビデオカメラ、マイクロスコープ等、顕微鏡以外の光学系にも適用可能である。[6]上述の実施の形態においては、全焦点画像生成装置や被写体高さ情報取得装置が、顕微鏡1光学系と画像処理部100とを有している。しかし、本発明が、光学系とは独立に構成された画像処理装置として実現されてもよい。この場合、例えば、画像処理装置とは独立に構成された撮像装置において、焦点位置が異なる複数の画像が撮像された後、これらの画像に対応するデータが光学的記録媒体などの記録媒体に記録される。画像処理装置は、この記録媒体からデータを読み出す読み出し部を有しており、記録媒体から画像データを読み出すと共に、読み出した画像データに対して上述の画像処理部と同様に画像処理を施す。[7]上述の実施形態においては、大域最適化の方法として確率伝播法が採用されているが、平均場近似、グラフカットなど、その他の大域最適化が用いられてもよい。"