(57)【要約】【課題】光学系のさまざまな収差の影響にロバストな画像ブレ推定法を提供する。【解決手段】ブレを含む観測画像の周波数領域において、ブレの方向に平行な方向には周波数成分が減衰し、ブレの方向に垂直な方向には周波数成分が減衰しない、というブレの周波数特性に基づいて、ブレ角度候補を一つあるいは複数選択するステップと、選択されたブレ角度候補によって限定された探索範囲において、ブレの長さ、あるいは、ブレの長さ及び角度を推定するステップと、からなる。【選択図】図１
2009187378.tif 000002
【特許請求の範囲】
【請求項１】
ブレを含む観測画像の周波数領域において、ブレの方向に平行な方向には周波数成分が減衰し、ブレの方向に垂直な方向には周波数成分が減衰しない、というブレの周波数特性に基づいて、ブレ角度候補を一つあるいは複数選択するステップと、
選択されたブレ角度候補によって限定された探索範囲において、ブレの長さ、あるいは、ブレの長さ及び角度を推定するステップと、
からなるブレ情報推定方法。
【請求項２】
前記ブレ角度候補を選択するステップは、ブレ方向に対応する角度を変数に持ち、真のブレ角度において特徴値をとる評価関数を用いるものである、請求項１に記載のブレ情報推定方法。
【請求項３】
前記ブレ角度候補を選択するステップは、関数Ｒ（θ）
2009187378.tif 000080
を用いるものである、請求項１，２いずれかに記載のブレ情報推定方法。
ここで、Ｔθは、周波数領域で基準軸と角度θをなして延出するテンプレート、|G（ξ，η）|は観測画像のスペクトル振幅、ｆ１（ｘ）はｘに関する単調非減少、かつ定数でない関数であり、角度θは、テンプレートを周波数領域における直流成分に対応する位置のＤＣ成分画素を中心として回転させることで変化させる。
【請求項４】
前記ブレ角度候補を選択するステップは、ブレ方向に対応する角度を変数に持ち、真のブレ角度において特徴値をとる評価関数を用いるものであり、
前記評価関数は、前記関数Ｒ（θ）、および／あるいは、関数Ｒ（θ＋９０deg）の関数である、請求項３に記載のブレ情報推定方法。
【請求項５】
前記評価関数Ｊ（θ）は、
2009187378.tif 000081
である、請求項４に記載のブレ情報推定方法。
【請求項６】
前記ブレ角度候補を選択するステップにおいて、ＤＣ成分画素を含む少なくとも一つのテンプレートと、ＤＣ成分画素を含まない少なくとも一つのテンプレートを用いてそれぞれＲ（θ）を求める、請求項３乃至５いずれかに記載のブレ情報推定方法。
【請求項７】
前記テンプレートは、一定の幅で延出する帯状テンプレートである、請求項３乃至６いずれかに記載のブレ情報推定方法。
【請求項８】
前記テンプレートは、ＤＣ成分画素を通る直線に対して線対称である、請求項３乃至７いずれかに記載のブレ情報推定方法。
【請求項９】
前記テンプレートは、幅方向に重み付けすることで設計される、請求項３乃至８いずれかに記載のブレ情報推定方法。
【請求項１０】
前記テンプレートは、前記重みが一定である、請求項９に記載のブレ情報推定方法。
【請求項１１】
前記テンプレートは、一定の幅寸法を備え、当該幅寸法をパラメータとして設計される、請求項１０に記載のブレ情報推定方法。
【請求項１２】
前記テンプレートは、幅方向に重みが変化する、請求項９に記載のブレ情報推定方法。
【請求項１３】
前記テンプレートは、幅方向の中間部位の重みを小さくしてなる、請求項１２に記載のブレ情報推定方法。
【請求項１４】
前記重みを軽くした中間部位にＤＣ成分画素が含まれる、請求項１３に記載のブレ情報推定方法。
【請求項１５】
前記重みを軽くした中間部位の幅方向両側に重みが一定の領域を有する、請求項１４，１５いずれかに記載のブレ情報推定方法。
【請求項１６】
前記中間部位の重みが０である、請求項１３乃至１５いずれかに記載のブレ情報推定方法。
【請求項１７】
前記テンプレートは、一定の全体の幅寸法と、一定の中間部位の幅寸法を備え、全体の幅寸法と中間部位の幅寸法とをパラメータとして設計される、請求項１６に記載のブレ情報推定方法。
【請求項１８】
前記ブレ角度候補を選択するステップは、前記評価関数Ｊ（θ）の出力として得られた特徴値を平滑化するステップを備えており、平滑化された特徴値に基づいてブレ角度候補を選択する、請求項１乃至１７いずれかに記載のブレ情報推定方法。
【請求項１９】
前記ブレの長さ、あるいはブレの長さ及び角度を推定するステップは、スペクトル振幅がゼロとなる周期あるいは、周期及び方向を検出するものである、請求項１乃至１８いずれかに記載のブレ情報推定方法。
【請求項２０】
ブレを含む観測画像の周波数領域において、ブレの方向に平行な方向には周波数成分が減衰し、ブレの方向に垂直な方向には周波数成分が減衰しない、というブレの周波数特性に基づいて、ブレ角度候補を一つあるいは複数選択する手段と、
選択されたブレ角度候補によって限定された探索範囲において、ブレの長さ、あるいは、ブレの長さ及び角度を推定する手段と、
からなるブレ情報推定装置。
【請求項２１】
ブレを含む観測画像におけるブレ情報を推定するためにコンピュータを、
ブレを含む観測画像の周波数領域において、ブレの方向に平行な方向には周波数成分が減衰し、ブレの方向に垂直な方向には周波数成分が減衰しない、というブレの周波数特性に基づいて、ブレ角度候補を一つあるいは複数選択する手段と、
選択されたブレ角度候補によって限定された探索範囲において、ブレの長さ、あるいは、ブレの長さ及び角度を推定する手段と、
として機能させるコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ブレ情報推定方法に関するものである。
【背景技術】
【０００２】
「ブレ」とは、光学系の視線と被写体の相対運動によって露光時間中に生じる画像移動量（motion blur）を言う。ブレを等速直線運動と仮定した場合は、推定すべきパラメータは、ブレの長さと方向（角度）である。ブレによって画像は劣化するが、その寄与は、周波数領域において周期的な暗線パターン（以下、縞模様と呼ぶ）を持ち、この縞模様が観測画像のフーリエスペクトルにも現れることが良く知られている。この縞模様は、縞模様の方向と垂直方向がブレの方向に対応し、縞模様間隔とブレ長さが対応する（具体的には縞模様間隔はブレ長さに反比例する）、という特徴を備えている。
【０００３】
Liらは画像のフーリエスペクトルに含まれるノイズ統計量について考察し、耐ノイズ性を得るため、ブレ長さを推定する前処理として、２次元フーリエスペクトルをブレ方向に沿って（縞模様に直交する方向に）平滑化することを提案した(非特許文献１)。具体的には角度を変えてスペクトルを積算し、極値が最も顕著に現れる方向が、縞模様を垂直に切っていると想定できることを利用して、これをブレ角度としている。これにより、SNRの低い画像からでもブレ情報が推定できるようになった。しかし、程度の激しい収差やデフォーカスがある場合、その各々がフーリエスペクトルにおいて極値を与えることから、上記の方法では推定が正しく行なえない、という問題点があった。
【０００４】
米司らは観測画像のスペクトル振幅の周期的に現れる極小値を、「くし型」のテンプレートを多数用意して、各々と相関を取ることで、被写体の影響やノイズに対するロバスト性を向上させた（特許文献１、非特許文献２）。縞模様検出には本来不要である、中央を通るテンプレートを、正しく長さを推定するための正則化項として用いており、パラメータチューニングが不要であるという利点がある。
【０００５】
Sakanoらはスペクトル振幅の局所勾配が、この暗直線に対して、理想的には直交することに着目し、Hough変換によって直線を検出し、ノイズの多い画像に対しても安定してブレパラメータを推定できることを示した（非特許文献３）。
【０００６】
このように、従来のブレ検出技術は、あらゆるテクスチャ(模様) の被写体に対して、また、ノイズが多い画像に対して、この縞模様を安定して検出する技術に力点が置かれてきた。しかし、従来のブレ検出技術が画像の劣化要因として、ブレ以外に本格的に考慮してきたのは、白色雑音のみであった。光学系が歪んでいる場合、収差と呼ばれる劣化要因が別途生じる。収差の中には、フーリエ領域において縞模様を持つもの（例えば、非点収差）があり、そのような収差を強く含んだ画像に対しては、縞模様を安定して検出する従来技術は適用できない。
【０００７】
ブレ角度推定方法に関しては、縞模様を利用しない手法も提案されている。橋本らは画像のフーリエスペクトルにおいて、直交する２方向について強度比をとり、値が最大となる位置をブレ方向と推定した（非特許文献４）。この方法は、物体の空間周波数情報が、ブレに沿った方向には強く減衰し、ブレに垂直な方向にはほとんど失われないという直線的なブレの性質を利用したものであり、縞模様は利用していない。非特許文献４では、これに続いて、違和感のない復元画像を与えるブレ長さを求めている。具体的には、さまざまなブレ長さを仮定した画像復元フィルタを設計し、これを元に物体の復元処理を行なう。この画像のうち、平均的に最も空間周波数が高く、かつリンギングの起こっていない復元画像を与えるブレ長さを、推定値とした。この手法はきわめて計算量が多い。また、方向性の強いテクスチャ、収差があると推定性能が劣化する。
【０００８】
関連する他の従来技術として、ブラインドデコンボリューションがある。これは、画像のみを基に、劣化要素と被写体を推定するアルゴリズムである。しかしながら、ブラインドデコンボリューションは、画像から、劣化要素を「収差、ブレ情報を含んだ総体」として推定するものであって、劣化要素を収差とブレ要素に分離する、という観点では報告例がなく、ブラインドデコンボリューションによって、収差を含んだ画像からブレを推定することはできない。
【特許文献１】特開２００６−２２１３４７号
【非特許文献１】QiangLI and Yasuo YOSHIDA. Parameter estimation and restoration for motion blurredimages (special section on digital signal processing). IEICE transactions onfundamentals of electronics, communications and computer sciences, Vol. 80, No.8, pp. 1430.1437, 19970825.
【非特許文献２】米司健一, 田中正行, 奥富正敏. 直線的手ブレ画像復元のためのPSF パラメータ推定手法(卒論セッション). 情報処理学会研究報告. CVIM,[コンピュータビジョンとイメージメディア], Vol. 2005, No. 38, pp. 47.52, 2005.
【非特許文献３】M.SAKANO, N. SUETAKE, and E. UCHINO. A PSF Estimation Based on Hough Transform ConcerningGradient Vector for Noisy and Motion Blurred Images. IEICE TRANSACTIONS onInformation and Systems, Vol. 90, No. 1, pp. 182.190, 2007.
【非特許文献４】橋本正一, 斎藤英雄. Psf のパラメータ分布を推定するシフトバリアントなぼけ画像の復元法. 電子情報通信学会論文誌. D-II, 情報・システム, II-情報処理, Vol. 77, No. 4, pp. 719.728,19940425.
【発明の開示】
【発明が解決しようとする課題】
【０００９】
本発明は、周波数領域において縞模様を伴う劣化要素（例えば収差やデフォーカス）にロバストな画像ブレ推定法を提供することを目的とするものである。
【課題を解決するための手段】
【００１０】
かかる課題を解決するために本発明が採用した技術手段は、ブレを含む観測画像の周波数領域において、ブレの方向に平行な方向には周波数成分が減衰し、ブレの方向に垂直な方向には周波数成分が減衰しない、というブレの周波数特性に基づいて、ブレ角度候補を一つあるいは複数選択するステップと、選択されたブレ角度候補によって限定された探索範囲において、ブレの長さ、あるいは、ブレの長さ及び角度を推定するステップと、からなるブレ情報推定方法、である。本発明は、縞模様に注目せずに、まずブレ角度を推定し、推定されたブレ角度によって限定された探索範囲でブレ長さを推定する点に特徴を備えている。ブレ画像の周波数領域は、典型的には、ＦＦＴによって得られるフーリエ領域であるが、ＤＣＴ等によるものでもよい。
【００１１】
一つの態様では、前記ブレ角度候補を選択するステップは、ブレ方向に対応する角度を変数に持ち、真のブレ角度において特徴値をとる評価関数を用いる。特徴値が出力された時の入力である角度がブレ角度候補となる。典型的な態様では、前記特徴値は最大値あるいは最小値である。例えば、ブレ角度の真値において、評価関数の出力は最大値となり、出力値が大きい順に入力である角度を１つ以上選択してブレ角度候補とする。
【００１２】
一つの態様では、前記ブレ角度候補を選択するステップは、関数Ｒ（θ）
【数１】
2009187378.tif 000003
を用いるものである。
ここで、Ｔθは、周波数領域（ξ，η）で基準軸と角度θをなして延出するテンプレート、|G（ξ，η）|は観測画像のスペクトル振幅、ｆ１（ｘ）はｘに関する単調非減少かつ、定数でない関数であり、角度θは、テンプレートを周波数領域における直流成分に対応する位置を中心として回転させることで変化させる。
角度θを設定するための基準軸は、後述する態様では、周波数領域の鉛直軸（η軸）であるが、他の方向の軸（例えば、水平軸（ξ軸））でもよい。
一つの態様では、単調非減少関数ｆ１は、単調増加関数であり、典型的な一つの態様では、ｆ１は対数関数（ｆ１（ｘ）＝ｌｏｇｘ）である。単調増加関数ｆ１は二乗関数（ｆ１（ｘ）＝ｘ２）でもよい。また、一つの態様では、単調非減少関数ｆ１は、ｆ１（ｘ）＝ｘである。
ある特徴抽出のための型をテンプレートと呼び、本明細書においては、特に、画像のスペクトル振幅から、スペクトル振幅の方向性を抽出するための型としてテンプレートＴθを用いている。具体的には、関数Ｒ（θ）によって、周波数領域における２次元関数であるテンプレートＴθによって設定される領域のスペクトル振幅を積算する。後述するように、テンプレートＴθは重みを備えてもよく、スペクトル振幅に重みが掛けられた値を用いて積算される。
前記角度θは、前記テンプレートを、周波数領域における直流成分に対応する位置（具体的には、周波数の直流成分に対応する画素であるＤＣ成分画素）を中心として回転させることで変化させる。通常、ＤＣ成分画素は周波数領域の中心に位置する。周波数領域へ変換した画像において、直流成分が画像の中心に来るように、画像の象限を入れ替える関数を用いることは当業者に周知である。
【００１３】
一つの態様では、前記評価関数は、前記関数Ｒ（θ）、および／あるいは、関数Ｒ（θ＋９０deg）の関数である。
一つの態様では、評価関数Ｊ(θ)は、角度θにおいて、Ｒ（θ）について単調非減少し、Ｒ（θ＋９０deg）について単調非増加し、定数でない関数である。この場合、評価関数Ｊ(θ)は、ブレ角度の真値において最大値を取る。より具体的には、評価関数Ｊ（θ）は、Ｒ（θ）について単調増加し、Ｒ（θ＋９０deg）について単調減少する関数である。このような関数としては、
【数２】
2009187378.tif 000004
が例示される。尚、上記関数の分母と分子を逆にしたものや、−（マイナス）を付したものが、同様に評価関数として用いられ得ることは当業者に容易に理解される。また、評価関数Ｊ（θ）は、Ｒ（θ）とＲ（θ＋９０deg）の比率によるものに限定されず、例えば、評価関数Ｊ（θ）＝Ｒ（θ）−Ｒ（θ＋９０deg）であってもよい。
一つの態様では、Ｊ（θ）は、ブレ角度真値の近傍において最大値を取る。この場合、角度θについて、Ｊ（θ）を平滑化した結果、ブレ角度の真値においてＪ（θ）は最大値を取る。平滑化の手法については公知の技術が多く利用でき、移動トリム平均フィルタや、移動メディアン平均フィルタなどが例示される。
一つの態様では、評価関数Ｊ(θ)は、Ｒ（θ）について単調非増加し、Ｒ（θ＋９０deg）について単調非減少し、定数でない関数である。この場合、評価関数Ｊ(θ)は、ブレ角度の真値において最小値を取る。より具体的な態様では、評価値Ｊ（θ）は、角度θにおいて、Ｒ（θ）について単調減少し、Ｒ（θ＋９０deg）について単調増加する関数である。
一つの態様では、Ｊ（θ）は、ブレ角度真値の近傍において最小値を取る。この場合、角度θについて、Ｊ（θ）を平滑化した結果、ブレ角度の真値においてＪ（θ）は最小値を取る。平滑化の手法については公知の技術が多く利用でき、移動トリム平均フィルタや、移動メディアン平均フィルタなどが例示される。
一つの態様では、評価関数Ｊ(θ)は、Ｒ（θ）、Ｒ（θ＋９０deg）のいずれか一方を用いるものでもよい。すなわち、ブレの方向に平行な方向には周波数成分が減衰し、ブレの方向に垂直な方向には周波数成分が減衰しない、というブレの周波数特性において、「ブレの方向に平行な方向には周波数成分が減衰する」あるいは「ブレの方向に垂直な方向には周波数成分が減衰しない」の一方を用いるものでもよい。
【００１４】
テンプレートは、一つの態様では、ブレ画像の周波数領域の画像全体に亘って延びているが、テンプレートは、画像全体幅よりも短尺であってもよい。一つの態様では、前記テンプレートは、一定の幅で延出する帯状テンプレートである。一つの態様では、前記テンプレートは、ＤＣ成分画素を通る直線に対して線対称である。換言すると、テンプレートは、その中央が周波数領域におけるＤＣ成分画素を通る。Ｒ（θ）の計算は１つのテンプレートによって行われ、１つのテンプレートがＮ本の帯からなる場合に、Ｎ本帯テンプレートと言う。１つのテンプレートの領域内に、重みが０の帯領域がＭ本ある場合、（Ｍ＋１）本帯テンプレートと言う。「帯の本数」は、重みの設計を通して設定される。テンプレートの領域とは、重みが０でない領域全体を包含する領域、すなわち包絡域のことを言う。
一つの態様では、前記テンプレートは、ＤＣ成分画素を含む少なくとも一つのテンプレートを含む。一つの態様では、前記テンプレートは、ＤＣ成分画素を含まない少なくとも一つのテンプレートを含む。ここで、ＤＣ成分画素を含むとは、ＤＣ成分画素における重みが０でないことを意味し、ＤＣ成分画素を含まないとはＤＣ成分画素における重みが０であることを意味する。
【００１５】
一つの態様では、前記特徴値を、複数のテンプレートを用いて求める。すな
わち、複数種類の帯状テンプレートを用いて前記特徴値Ｊ（θ）を求め、各帯状
テンプレートを用いて取得された特徴値Ｊ（θ）からブレ角度候補を各々１つあるいは複数選択し、その集合和を、前記複数のテンプレートによるブレ角度候補とする。より具体的な一つの態様では、前記複数のテンプレートは、ＤＣ成分画素を含む少なくとも一つのテンプレートと、ＤＣ成分画素を含まない少なくとも一つのテンプレートを含み、前記ブレ角度候補を選択するステップにおいて、ＤＣ成分画素を含む少なくとも一つのテンプレートと、ＤＣ成分画素を含まない少なくとも一つのテンプレートを用いてそれぞれＲ（θ）を求める。
【００１６】
一つの態様では、前記テンプレートは、幅方向（延出方向に直交する方向）に重み付けすることで設計される。一つの態様では、前記テンプレートは、前記重みが一定である。さらに、一つの態様では、前記テンプレートは、一定の幅寸法を備え、当該幅寸法をパラメータとして設計される。一つの態様では、前記幅寸法は、設計評価値ＪＴの確率的期待値を最大化することで最適化される。ここで、ＪＴは、marginJ/δＪである。一つの態様では、前記幅寸法は、２５６×２５６画素の画像に対して約１１画素、例えば、９〜１３画素の範囲の比率で設計される。一つの態様では、前記幅寸法は１画素である。
【００１７】
一つの態様では、前記テンプレートは、幅方向に重みが変化する。
幅方向に重みが変化するテンプレートにおいて、一つの態様では、前記テンプレートは、幅方向の中間部位の重みを小さくしてなる。さらなる態様では、前記重みを軽くした中間部位にＤＣ成分画素が含まれる。さらなる態様では、前記重みを軽くした中間部位の幅方向両側に重みが一定の領域を有する。一つの典型的な態様では、前記中間部位の重みが０である。さらなる態様では、帯状テンプレートは、全体の幅寸法と中間部位の幅寸法とをパラメータとして設計される。一つの態様では、この当該パラメータは、marginJによって決定される。
幅方向に重みが変化するテンプレートにおいて、一つの態様では、前記テンプレートは、幅方向の中間部位の重みを大きくしてなる。さらなる態様では、前記重みを大きくした中間部位にＤＣ成分画素が含まれる。例えば、３本帯テンプレートにおいて、中央の帯の部分の重みを両側の帯の重みよりも大きく設定することもできる。
【００１８】
一つの態様では、前記ブレ角度候補を選択するステップは、前記評価関数Ｊ（θ）の出力として得られた特徴値を平滑化するステップを備えており、平滑化された特徴値に基づいてブレ角度候補を選択する。
【００１９】
一つの態様では、前記ブレの長さ、あるいはブレの長さ及び角度を推定するステップは、スペクトル振幅がゼロとなる周期及び方向を検出するものである。
【００２０】
本発明に係るブレ情報推定方法を実行するためのハードウェア構成としては、入力手段、記憶手段、演算手段、出力手段等を備えたコンピュータが用いられ、本発明は、ブレ情報推定装置、ブレ情報推定のためのコンピュータプログラムとしても提供される。
すなわち、
ブレを含む観測画像の周波数領域において、ブレの方向に平行な方向には周波数成分が減衰し、ブレの方向に垂直な方向には周波数成分が減衰しない、というブレの周波数特性に基づいて、ブレ角度候補を一つあるいは複数選択する手段と、
選択されたブレ角度候補によって限定された探索範囲において、ブレの長さ、あるいは、ブレの長さ及び角度を推定する手段と、
からなるブレ情報推定装置、
ブレを含む観測画像におけるブレ情報を推定するためにコンピュータを、
ブレを含む観測画像の周波数領域において、ブレの方向に平行な方向には周波数成分が減衰し、ブレの方向に垂直な方向には周波数成分が減衰しない、というブレの周波数特性に基づいて、ブレ角度候補を一つあるいは複数選択する手段と、
選択されたブレ角度候補によって限定された探索範囲において、ブレの長さ、あるいは、ブレの長さ及び角度を推定する手段と、
として機能させるコンピュータプログラム、
として提供される。
【発明の効果】
【００２１】
本発明は、縞模様に注目せずに、まずブレ角度を推定し、推定されたブレ角度によって限定された探索範囲でブレ長さを推定するので、縞模様を伴う劣化要素、例えば収差やデフォーカス、に対してロバストである。角度推定テンプレートの設計を最適化することで、方向性の強い劣化要素、例えば収差や、方向性の強いテクスチャを持った被写体、に対してロバストとなる。ＤＣ成分画素を含むテンプレートと、含まないテンプレートを併用してＲ（θ）を取得するものでは、両者の異なる性質を互いに補完することで、ブレ角度を限定する際に、正解（すなわち、ブレ角度の真値）がブレ角度候補から漏れるおそれを減じることができる。
【発明を実施するための最良の形態】
【００２２】
１．本発明の概要
１．１  基本的な流れ
本発明の目標は、ブレを等速直線運動であると仮定し、ブレの長さρと、角度θを推定することにある。
ステップ１：撮影画像のフーリエスペクトル振幅を計算する。
まず撮影画像を適切に窓処理する。窓処理としては、直交型２次元ハニング窓が例示されるが、窓処理はこれに限定されるものではない。窓処理した画像をｇとする（図２６（Ａ）参照）。画像ｇのフーリエスペクトル振幅｜Ｇ｜を次式で求める（図２６（Ｂ）参照）
【数３】
2009187378.tif 000005
Ｆ［ｘ］で、ｘの２次元フーリエ変換を表わす。ここでは、振幅は、画像スペクトルの絶対値として求められる。
【００２３】
ステップ２：フーリエスペクトル振幅の角度依存性評価値Ｒ（θ）を計算する。Ｒ（θ）は、ブレ画像のスペクトル振幅と、ブレ画像の周波数領域で鉛直軸方向と角度θをなして延出するテンプレートと、の相関値であると考えることができる。図２に示すように、帯状のテンプレートに沿ってスペクトル振幅を積算する。幅１８０degの連結した区間内部で、θを変えて、ＤＣ成分画素を中心としてテンプレートを回転させながら、R（θ）、R（θ＋９０deg）を計算する。θは、例えば、５deg毎に計算されるが、θの変化量は当業者において適宜設定される。以上により、１次元変数θに対して関数値R（θ）が求まる。また、これにより、R（θ＋ ９０deg）が求まる。
【００２４】
ステップ３：角度推定評価値J（θ）を計算する。ステップ２で得られたR（θ）を元に、関数値J（θ）を計算する。
【数４】
2009187378.tif 000006
ブレ方向θにおいて、R（θ）は、ブレと垂直方向であり、周波数成分は保存され、R（θ＋９０deg）はブレに平行方向であり、周波数成分が減衰する。したがって、ブレ方向θ＝θtrueでJ（θ）は最大となる。
【００２５】
ステップ４：評価値J（θ）を平滑化する。J（θ）を平滑化する。平滑化は、異常に大きなノイズを低減するのが目的である。平滑化処理としては、移動メディアンフィルタ、移動トリム平均フィルタなどが例示される。smooth{f}は、関数fを平滑化する関数を表わすものとする。
【数５】
2009187378.tif 000007
【００２６】
ステップ５：ブレ角度候補を絞る。J（θ ）を大きい順にN 個選び、そのJを与える角度θの集合を｛θ（ハット）｝とおく。（ハット）は推定結果、あるいは予測値であることを表わす。括弧｛x｝は、複数の要素変数xからなる変数配列であることを表わす。Sort｛f｝は、１つあるいは複数の要素を含む配列fを昇順に並びかえる関数を表わすものとする。endは配列要素の識別子の最後であることを表わす。
【数６】
2009187378.tif 000008
【００２７】
ステップ６：ブレの長さ、角度を推定する。Ｎ個のθによって限定された探索範囲で、ブレ長さ、あるいは、ブレ長さと角度を推定する。すなわち、長さρについては、全範囲、角度θ については、ステップ５ で得られたブレ角度候補｛θ（ハット）｝のみを探索領域とする。これを概念的に図６に示す。一つの態様では、Ｎ＝１であり、この場合は、長さρのみの推定が行われる。ステップ６には、縞模様を検出することでブレ情報を推定する従来技術（特許文献１、非特許文献１，２、３に開示された手法を用いることができる。後述する実施形態では、非特許文献２（特許文献１）に記載された手法を採用している。
【００２８】
本発明によって角度の候補が限定されている場合に、非特許文献３に開示された手法によれば、次のようにブレ情報を推定できる。画像はNxM画素とし、M=Nとする。M＝Nでない場合においても、非特許文献３に開示された手法を用いて同様に推定が可能である。混乱を避けるため、記号は本明細で再定義した。
準備ステップ：矩形部分空間の定義
準備ステップは、逐次実施する必要はない。多数回の画像処理に対しても最初の1回のみ実施すればよい
座標（ξ,η）を通り、η軸と角度θをなす直線の、原点からの距離qは、次式で与えられることが知られている。
【数７】
2009187378.tif 000009
このqは、ξとηを定数と見れば、パラメータ空間(θ,q)において、曲線を描き、これはハフ曲線として知られている。ここで、空間(θ,q)を小さな矩形領域に分割し、各々を矩形部分空間とよぶことにする。すなわち、θを△θ間隔で、qを△q間隔で分割する。△θは、本発明における角度推定ステップにおいて利用した角度刻み幅を用いることが例示される。具体的には、5度が例示される。また、△qは、1画素が例示される。
ステップ１：観測画像のフーリエスペクトル振幅｜G（ξ,η）｜を計算する。
ステップ２：｜G｜の局所勾配θg(ξ,η)を計算する。｜G｜の各画素において、2ｌ-１画素四方の正方形領域近傍の画素値を用いて、次式により与えられる。ただし、左辺値は0以上360度未満の角度とする。
【数８】
2009187378.tif 000010
ステップ３：暗線を検出する。
準備ステップにおいて定義した部分矩形領域に対し、投票を行なう。投票とは、曲線qが空間(θ,q)の部分矩形領域を通過したとき、次式で与える値uを得点（スコア）として、その部分矩形領域に与える作業を言う。
【数９】
2009187378.tif 000011
｜G（ξ,η）｜の全定義域、すなわち全画素に対して投票が完了したら、投票結果、すなわちこのときの部分矩形領域の得点をuf(θ,q)とする。このufを用いて、次式で与える関数hを求める。
【数１０】
2009187378.tif 000012
wは例えば1と設定する。
h（θ,q）を（θ,q）空間内で平滑化し、hs（θ,q）を得る。平滑化は、（θ,q）で定義される2次元正規分布との畳み込み演算によって行なわれる。平均を（0,0）とし、分散は（△θ2, △q2）とした正規分布が例示される。
本発明により限定された角度候補{θ（ハット）}のみを探索範囲として、hs（θ,q）が最大となるパラメータ（θ,q）を求める。これをθ（ハット）、およびq（ハット）とする。
ステップ４：ブレ情報を推定する。
ρ（ハット）＝N /q（ハット）を計算する。
ブレ長さの推定結果を、ρ（ハット）とする。ブレ角度の推定結果はθ（ハット）とする。
以上で、ブレ情報の推定が完了する。
【００２９】
また、長さ推定は、フーリエスペクトルの縞模様の間隔を利用するものに限定されず、非特許文献４に開示された手法を用いても良い。
本発明によって角度の候補が限定されている場合に、非特許文献４に開示された手法によれば、次のようにブレ情報を推定できる。混乱を避けるため、記号は本明細で再定義した。
ステップ１：本発明により限定された角度候補{θ（ハット）}の要素をθ（ハット）iとおき、そのおのおのについて、PSFのスペクトルを計算する。具体的には、長さρと、θ（ハット）iと、に対応するブレのPSFのスペクトルは式(2.5)で計算されるH(ξ,η)を、DC成分画素を中心として角度θ（ハット）iだけ回転させたものとなる。これを、H（ξ,η;ρ,i）と書くことにする。
ステップ２：既知の観測画像g、あるいはgの一部（部分領域）と、既知のPSFスペクトルH（ξ,η;ρ,i）と、から公知技術であるウィーナーフィルタによって復元画像を求める。復元画像を、f（チルダ）（ρ,i）とする。具体的には、次式で与えられる。*は複素共役を表わす。ΓはパワースペクトルのSNRを表わす定数である。Γを画像から推定する公知技術もあるが、典型的には、0.01など小さな値に設定すればよい。IFFT{f}は、関数fの２次元 逆フーリエ変換を表わす。
【数１１】
2009187378.tif 000013
ステップ３：f（チルダ）（ρ,i）と、f（チルダ）（ρ+Δρ,i）と、の相関値C（ρ,i）を、ρおよびiを変えながら求める。Δρは、例えば1画素とするが、より大きな値でも良い。
【数１２】
2009187378.tif 000014
ステップ４：C（ρ,i）の最大値を与えるパラメータの組ρ（ハット）と、i（ハット）を求める。
ステップ５：ブレ長さの推定結果を、ρ（ハット）とする。ブレ角度の推定結果はθ（ハット）i（ハット）とする。
以上で、ブレ情報の推定が完了する。
【００３０】
１．２  帯状テンプレート
フーリエ面を規定する直交座標をξ、ηとし、η軸に平行なテンプレートを定義する。このテンプレートを角度θ回転することによって、角度推定を行なう。
帯状テンプレートを、「テンプレートがＤＣ成分画素を含むか否か」で分類する。「帯の本数」、「ＤＣ成分画素を含むか否か」の2つの観点にから、ブレ角度推定テンプレートを分類したものを図５に示す。ＤＣ成分画素を含む線分について線対称であるという条件を付加し、（設計パラメータが少ないという意味で）単純なテンプレートを設計しようとすると、ＤＣ成分画素を含むテンプレートの代表が  1本帯テンプレート、ＤＣ成分画素を含まないテンプレートの代表が２本帯テンプレートとなる。図１１に示すように、帯状テンプレートは、外幅と内幅の２つのパラメータによって特定され、内幅＝０の場合は１本帯であり、内幅が＞０の場合は２本帯を構成する。一定でない重みを持ったテンプレートや、３本帯、４本帯テンプレートも設計可能である。
【００３１】
［ＤＣ成分画素を含むテンプレートの設計］
ＤＣ成分画素を含むテンプレートは、重みwaはna個のパラメータpa_1,pa_2,,,pa_naによって形状を規定する（図６（Ａ））。waはこれらのパラメータの関数として与えられ、wa(ξ,η; pa_1,pa_2,..,pa_na)と書くことができる。waはDC成分画素に対応する座標において非ゼロの値をとり、典型的には正値を取る。waに関する制限は、DC成分画素に対応する座標における値以外にはない。パラメータ個数naについては、1以上の整数であること以外に制限はない。ただし、典型的には、ごく少数、例えば唯一つのパラメータで規定されるwaであっても、テンプレート設計評価値JTを用いた設計の効果が得られる。
設計方法としては、後述する2.4節に記載の方法で設計することができる。より汎用的には、評価値JTを最大化するパラメータの組み合わせ(pa_1,pa_2,,,pa_na)optimalを設計解とする。
テンプレートの設計例を説明する。簡単のため、ただ1つのパラメータTmaxによって形状を規定できるテンプレートを設計する（図６（Ｃ））。このテンプレートは次の特徴を持つ：重みwaは全幅Tmaxにわたって一律；重みwaはフーリエ面においてDC成分画素を通る直線に対して線対称である。設計すべきパラメータはTmaxのみであり、評価値JTを最大化するTmaxを設計解とする。
【００３２】
［DC成分画素を含まないテンプレートの設計］
DC成分画素を含まないテンプレートは、重みwbはnb個のパラメータpb_1,pb_2,,,pb_nbによって形状を規定する（図６（Ｂ））。wbはこれらのパラメータの関数として与えられ、wb(ξ,η; pb_1,pb_2,..,pb_nb)と書くことができる。wbはDC成分画素に対応する座標においてゼロをとる。wbに関する制限は、DC成分画素に対応する座標における値以外にはない。パラメータ個数nbについては、1以上の整数であること以外に制限はない。ただし、典型的には、ごく少数、例えば2つのパラメータで規定されるwbであっても、marginJを用いた設計の効果が得られる。
後述する2.6節記載の方法で設計することができる。より汎用的には、次の条件１を満たすパラメータの組(pb_1,pb_2,,,pb_nb)optimal
を設計解とする。
＊ 条件１）予測されるブレ長さ範囲内で、常にmarginJ＞０となるパラメータ
の組。
条件１を満たすパラメータの組み合わせが複数ある場合には、これらを全て選んでも良いし、任意の一つ、あるいは複数を選んでも良い。あるいは、条件１を満たすパラメータの組み合わせの中で、重みwbを変数とする関数値が特徴値をとるようなパラメータの組み合わせを一つ、あるいは複数選んでも良い。あるいは、典型的な例として、条件１を満たすパラメータの組み合わせの中で、重みｗｂの積算値が小さい方から1つ、あるいは複数選んでもよい。
テンプレートの設計例を説明する。簡単のため、2つのパラメータTmax、Tminのみによって形状を規定できるテンプレートを設計する（図６（Ｄ））。設計すべきパラメータはTmax、Tminの計２つである。後述するように予測されるブレ長さ範囲内で常にmarginJ＞０となるような、最大のTmin を求める。これは、条件１を満たすパラメータの組み合わせの中で、重みｗｂの積算値が最小となるパラメータの組を１つ選んでいることに対応する。
【００３３】
一つの態様では、ＤＣ成分画素を含むテンプレートと、含まないテンプレートを併用することで、両者の異なる性質を互いに補完する。
ＤＣ成分画素を含むテンプレートは、信号を無駄なく拾いやすいという利点がある一方、ノイズを拾いやすいという欠点がある。すなわち、信号も雑音も大きい。これに対して、ＤＣ成分画素を含まないテンプレートは、ノイズを拾いにくい可能性がある（ノイズの種類に依存する。具体的には、ノイズの空間周波数依存性次第である）一方、信号も拾いにくいという欠点がある。すなわち、信号も雑音も小さい。
ブレ情報は、空間周波数ゼロ（すなわち直流成分）で最大値をとる。したがって、テンプレートがＤＣ成分画素に感度をもつと、信号を拾いやすく、テンプレートがＤＣ成分画素に感度をもたないと、信号を拾いにくい。ブレ推定問題におけるノイズ（収差、被写体の影響、デフォーカス＝ピンボケなど）は、空間周波数ゼロで最大値をとる。したがって、テンプレートがＤＣ成分画素に感度をもつと、ノイズを拾いやすく、テンプレートがＤＣ成分画素に感度をもたないと、ノイズを拾いにくい。
ＤＣ成分画素を含むテンプレートと、含まないテンプレートを併用することで、ブレ角度を限定する際に、正解（すなわち、真のブレ角度）がブレ角度候補から漏れる危険性が減ることを期待できる。
【００３４】
２．本発明に係るブレ推定方法の実施形態
本発明に係るブレ推定方法の実施形態について詳述する。2.1 節にてブレを画像から推定する基本原理に触れた後、従来手法による解法とその問題点を2.2 節にまとめる。本実施形態では、限定された探索範囲でのブレ長さ推定、あるいは、ブレ長さ及び角度推定に、2.2に述べる従来手法を用いる。続いて、2.3 節で、この問題点を解決する枠組みを提供する。この枠組みにもとづき、2.4 節では問題の定式化を行なう。2.5節は、２種類のテンプレートを併用する意義と方法論を与える。以上をまとめ、提案するブレ推定手法全体の手順を2.6節に整理する。2.7節で従来手法との比較シミュレーションを説明する。
【００３５】
2.1 ブレ推定の原理
物体ｏが劣化過程を経て画像ｇとして得られるモデルは一般に、＊を２次元畳みこみ演算子として、
【数１３】
2009187378.tif 000015
と表わされる。劣化過程のうち、hはPSFと呼ばれ、物体上の1 点の、画像平面上における強度分布を表わす。nは加法ノイズである。式(2.1) をフーリエ変換すると
【数１４】
2009187378.tif 000016
とかける。G, H, O, N はそれぞれ、g, h, o, n の２次元フーリエ変換である。視線が露光時間中に等速直線運動した場合のPSF は線分としてモデル化できる。ブレの長さをρ、角度θ =0 すなわちx 方向とすると、PSF は
【数１５】
2009187378.tif 000017
と書かれる。δ(y)はディラックのデルタ関数である。画像の画素数がN×N
pixelであるとすれば、式(2.3)をフーリエ変換すると
【数１６】
2009187378.tif 000018
となる。明らかに、この関数は周期
【数１７】
2009187378.tif 000019
ごとにゼロ点を持つ。また、η方向には一定であることから、このゼロ点は実際には暗線となる。このゼロ点により、フーリエ領域において縞模様が形成される。以後、ブレに由来するこの周期的なゼロ点を単に「縞模様」と呼ぶことにする。観測画像のスペクトル振幅は式(2.2)より容易に
【数１８】
2009187378.tif 000020
となるが、縞模様は観測画像のスペクトル振幅｜G｜にも保存されることが分かる。周波数領域におけるブレ推定が利用している原理は、観測画像のスペクトル振幅に現れる縞模様の周期と方向を調べることによって、ブレのPSFの長さρと角度θを検出しようというものである。
【００３６】
2.2 従来手法と問題点
本節では、被写体のテクスチャによらず安定してブレ角度・長さ推定が可能な非特許文献２（特許文献１）に開示された手法の原理と手順を紹介する。続いて、この手法が、デフォーカスや特定の収差に弱いことを指摘する。
非特許文献２（特許文献１）の手法は、図７に示す２次元テンプレートTρ，θにより、観測画像の２次元フーリエスペクトル振幅|G(ξ，η)| における暗縞模様を検出するものである。具体的には、両者の相関を計算し、式(2.8)に示す評価値を最小化するパラメータ(ｖp，θp) を推定する。ｖpはテンプレートに最も適合した縞模様の間隔であり、式(2.6)に従ってブレ長さρp=N/ｖpに対応する。すなわち、ブレ情報は、(ρp，θp) として得られる。
【数１９】
2009187378.tif 000021
なお、対数をとっているのは、高周波帯域での縞模様を推定に利用するためである。
この手法は、方向性の強いテクスチャや収差には耐性がある。なぜならば、縞模様のみに着目しており、スペクトルの方向性には全く注目していないためである。しかし一方で、縞模様を伴う劣化要素には耐性が低いことが予想される。縞模様を伴う劣化要素にはデフォーカス、非点収差、球面収差などがあるが、このうち非点収差を与えた画像のフーリエスペクトル振幅を図８に示す。0degと90deg方向に縞模様が見え、これはブレに由来する縞模様とは無関係の方向に走っていることが確認できる。この収差に対し、評価値Jρ，θがよく適合する可能性があり、このとき、ブレ推定値は真値とは大幅に異なるものとなる。これが、縞模様のみに着目したブレ推定手法に共通する問題点である。
【００３７】
2.3 収差・デフォーカス耐性向上のためのブレ推定の枠組み
2.3.1 課題解決の方針
本発明が解決しようとする課題は、画像からのブレ情報推定法において、これまで検討されてこなかった劣化要素、特に収差に対する推定性能を改善することである。本実施形態では、これを解決する枠組みを与える。枠組みは、本節で述べる２つの方針からなり、その概要を図１に示す。
方針１）ブレ角度を限定してから長さを推定する
式(2.6)
から明らかなように、ブレの長さはフーリエスペクトル振幅に現れる縞模様の間隔ｖと対応付けられるため、その推定については縞模様を検出する方針を本実施形態でも採用する。一方、ブレの角度を推定するには、縞模様を検出する必要はないことに注目する。収差・デフォーカスの影響を受けると、観測画像のスペクトルは縞模様を持ち、方向依存性を持つ。縞模様の検出を行なうことなく、角度を推定することがもし可能であれば、その推定角度に探索範囲を限定して縞模様の間隔を調べ、ブレの長さを推定することができる。角度候補が正しく限定されていれば、縞模様を伴う劣化要素に関する耐性向上が期待できる。実装上は、図３に示すように、ブレに由来する縞模様に平行
なテンプレートに限定して、観測画像スペクトル振幅との相関を計算すればよいことになる。
【００３８】
方針２）ブレ角度推定に用いるテンプレートを最適化する
方針１は、スペクトル振幅の縞模様に着目せずにブレ角度を推定することを要求する。その方法として、スペクトルの方向依存性を調べることが挙げられる。具体的には、観測画像のスペクトル振幅と、図２に示すような帯状のテンプレートとの相関を調べ、最大となる角度を求めればよい。このテンプレートを、角度θを推定するテンプレートという意味でTθと呼ぶ。スペクトルの方向依存性を調べる方法は、方向依存性を持った劣化要素にできるだけロバストでなければならない。このようなテンプレートを用いたブレ角度推定において、特に、方向依存性の強い劣化要素に対しての耐性を向上させるには、テンプレートを最適化することが望ましい。この課題達成を、本実施形態ではテンプレートTθを、ブレ推定問題のSNR を高めるという意味で最適化する問題に帰着する。
【００３９】
提案する手法は次の手続きをとることによって、縞模様を検出することなく角度を推定し、最終的に異方性ノイズ（方向依存性を持つ劣化要素）の影響を受けずにブレ情報を推定する。
(1) 角度依存性のあるノイズの影響を受けにくい角度検出テンプレートTθを設計する。
(2) Tθを用いてブレ角度評価値を計算し、これを角度方向に平滑化する。
(3) 平滑化済みの評価値の良い候補を選択する。
(4) 限定された解空間に対し、Tρ，θを用いてブレの角度と長さを推定する。
【００４０】
2.3.2 手順の詳細
［角度評価値の計算］
ブレ角度を推定するための評価値を調べる方法を述べる。
観測画像のフーリエスペクトル振幅|G|の方向性は
【数２０】
2009187378.tif 000022
によって調べることができる。画像の空間スペクトル振幅は指数関数的に減少することが知られており、このままでは高周波成分が推定に使えず、都合が悪い。そこで、振幅の対数をとり、
【数２１】
2009187378.tif 000023
によって方向依存性を調べることにする。画像は画素によりサンプリングされているから、この関数R(θ)は、角度検出テンプレートTθとの相関を取った出力と考えることができる。
【数２２】
2009187378.tif 000024
直線的なブレの周波数特性は、ブレの方向θtrueに平行な方向には高周波成分が減衰し、直交する方向には全く減衰しない。これは、式(2.7) より、観測画像にも反映される。図１０で説明すれば、R(θ)は、角度θ=θtrueで最大値をとり、θ=θtrue+90degで最小値を取ることが期待される。したがって、θ方向とθ+90deg方向のスペクトルの強さの比が最大となる方向が、ブレ方向となると考えられる。そこで、ブレ角度推定評価値として
【数２３】
2009187378.tif 000025
を定義する。
【００４１】
［角度推定テンプレートTθの設計パラメータ］
テンプレートTθを具体的に設計する。パラメータとして、図１１に示す(Tmax，Tmin) を採用する。これはそれぞれバンドパスフィルタの上側遮断周波数、下側遮断周波数を意味する。パラメータと幅W との関係は、次式で表わされる。
【数２４】
2009187378.tif 000026
Tmaxの設計方法は、2.4 節に述べる。Tminについては2.6 節に述べる。なお、本実施形態においてはテンプレートの重みを一様とする。これは、問題を簡単化するためである。重みの設計には、異方性ノイズのスペクトル分布の情報が必要となる。
【００４２】
［角度推定評価値の平滑化］
ブレを直線運動に限定すると、ブレに由来するスペクトルの角度依存性は単峰性となる。ここで重要なのは最大値を与える角度のみであり、些細な成分は意味を持たない。従って、平滑化を行ない、大局的な最大値を与える角度を求める。
【数２５】
2009187378.tif 000027
平滑化に用いるフィルタは移動メディアンフィルタ、移動トリム平均フィルタなどが適する。図１２に、平滑化する前の角度推定評価値と、移動トリム平均フィルタによって平滑化したものを示す。平滑化前の評価値は多峰性であり、囲みで示したように、大角度を隔てた位置に複数の最大値が見られる場合がある。平滑化によって、大局的に見た最大値を与える角度を推定できる。
【００４３】
［角度候補の選択］
平滑化した角度推定評価値の、大きい方から（Nθ≧1) 個を選び、これを与える角度を選択する。この角度を、推定値の候補とする。
【数２６】
2009187378.tif 000028
Nθの増加は、必ずしも推定性能を改善しない。なぜならば、角度候補を限定せずに縞模様のみ利用してブレ情報を推定する手法に漸近し、2.2 節にて指摘した問題点が生じるためである。設計されたテンプレートTθのブレ角度推定性能が十分高ければ、候補の数Nθ=1とすべきと考えられる。
【００４４】
［ブレ情報の推定］
ブレ情報(ρ，θ)は、式(2.8) に従って推定する。すなわち、画像のフーリエスペクトル振幅|G|と、テンプレートTρ，θとの相関を計算し、最小評価値を与える角度と長さとして求める。角度候補が1つであった場合、ここで求めるのは長さのみとなる。提案する手法が従来手法と異なるのは、このときの探索範囲であり、角度推定結果に基づき次のとおり限定する。
＊ 長さρについては、全範囲
＊ 角度θについては、既に限定した推定角度候補{θ（ハット）} 
これを、模式的に図４に示す。図４に示すように、角度θの候補が絞り込まれる（網掛けした扇形領域）。本発明ではこの内部についてのみ探索し、最も適合度が高いと判定された(ρ,θ)を最終的な解と見なす。従来技術では角度θを絞り込まずに全ての(ρ,θ)を探索するため、例えば収差に由来する強い縞模様を適合度が高いと判定する場合がある（灰色の正方形領域）。このように、縞模様を利用した長さ推定を開始する段階で、探索範囲を限定しておくことで、縞模様を伴う劣化要素に対する耐性向上を期待する。
【００４５】
2.4 モデルにもとづくブレ角度推定テンプレート設計法
前節で利用を提案したブレ角度推定テンプレートTθはパラメータTmax，Tminをもつ。本節では、その設計方法を提示する。まず、ブレ角度推定過程をモデル化し、これに現れるノイズを、方向依存性の有無に関して分類する。続いて、これを定式化し、仮想的な観測画像のスペクトルを求める。これに基づいてテンプレートを設計する。
【００４６】
2.4.1 ノイズ特性の整理
画像からブレを推定する問題においては、ブレが検出すべき信号であり、そのほかの要素は全てノイズと考えられる（図１３）。このノイズを、方向性を持つ（異方性）か、持たない（等方性）かに従って分類する。本実施形態におけるモデル化の方法を併せて表2.1 に示す。
【表１】
2009187378.tif 000029
【００４７】
［等方性ノイズ］
（ア）被写体
画像信号をフラクタルと考えると、画像に含まれる長い空間周期（低周波）の特徴を少ないパラメータでよく表現することが指摘されている。自然画像のフーリエスペクトル振幅は、空間周波数fの逆数に概ね反比例する（q=1）。より具体的には、多くの自然画像に対して、
【数２７】
2009187378.tif 000030
に比例する。
【００４８】
（イ）回折
光学系は有限の開口（口径）をもつため、これを通過した光波は回折（diffraction）を受ける。たとえ光学系が理想的（無収差）であっても、回折効果は避けることができない。影響はデフォーカスに比べるとはるかに小さいが、考慮することにする。口径D、焦点距離fl、波長λ なる円形開口の無収差光学系によるPSFは次式で与えられる。
【数２８】
2009187378.tif 000031
J1は第一種ベッセル関数である。
【００４９】
（ウ）デフォーカス
理想的な結像点からの光軸方向の変位であり、俗にピンボケと呼ばれる。図１４に示すように、口径D、焦点距離fl、デフォーカスδf のとき、検出器位置における光束の半径rlは全く幾何光学的に、
【数２９】
2009187378.tif 000032
と求まる。このrlは錯乱円半径と呼ばれる。これを用いてデフォーカスによるPSFは次式のpillbox型の関数で表現できる。
【数３０】
2009187378.tif 000033
【００５０】
［異方性ノイズ］
異方性ノイズは強度の角度分布が多岐にわたり、モデル化に必要な情報が得られない場合が多い。例えば、被写体のテクスチャや、コマ収差の有無、などである。本実施形態では、これらを詳細にモデル化することをせず、きわめて単純化したモデルを提示し、利用する。具体的には2.4.7 節に示す。
【００５１】
［ブレ推定におけるノイズ対策上の問題点］
定性的には、等方性ノイズと異方性ノイズとで、その影響を低減するアプローチが正反対であり、最適設計点が不明である、という問題がある。
等方性ノイズは高空間周波数帯域で、ブレ信号を減衰する働きがある。従って、このノイズの影響を低減するには、より低空間周波数帯域でブレ角度推定を行なえばよい。具体的には、
−Tmin（テンプレート間空隙）を小さく設計すること、
−テンプレート幅W をある程度大きく設計すること、
によって達成されると考えられる。
異方性ノイズの影響を低減するには、低空間周波数帯域でのブレ角度推定を避けたり、ノイズの影響を選択的に除去したりすればよい。具体的には、
−Tmin（テンプレート間空隙）を大きく設計すること、
−角度推定評価値を強い異常値の存在を考慮して平滑化すること、
−角度推定値候補を複数選ぶこと、
によって達成されると考えられる。
【００５２】
2.4.2 角度推定テンプレートTθ設計方針
テンプレートを設計する目的で、図１５Ａ、１５Ｂ、１５Ｃに示すモデルを用いる。まず、比較的モデル化しやすい等方性ノイズの各要素を定量化し、これらの寄与をHiso なるスペクトルで代表させる。一方、ブレについてはモデル化できるから、信号が等方性ノイズによって劣化された仮想的な画像のスペクトルGisoを考えることができる。その振幅|Giso|と、さまざまなパラメータ(Tmax，Tmin) をもったテンプレートTθとの相関値を計算すると、ブレ推定における余裕を、等方性ノイズの影響の考慮を済ませたかたちで見積もることができる。続いて、簡単化した異方性ノイズモデルHanisoに対する感度を定義し、余裕と感度の比率を求める。この比率は等方性/異方性ノイズに対するロバスト性指標と解釈できるため、これをTθの設計評価値として採用する。
【００５３】
以上の設計を実現するための方法論を次節以降に展開する。まず準備として、ノイズを等方性/異方性に分類し、ブレ角度推定モデルを構築する。具体的には、信号（ブレ）と等方性ノイズをモデル化し、評価値のマージンを求める。続いて、異方性ノイズが発生したときの評価値変分を求める。
【００５４】
2.4.3 ノイズ特性のモデル化
2.4.1 節で列挙したノイズ要素を、周波数領域でモデル化する。各要素はそれぞれ、パラメータを持つが、ここでは既知としておく。パラメータの設定方法は2.4.9 節にまとめて記す。fで空間周波数を表わすことにする。すなわち、
【数３１】
2009187378.tif 000034
とする。2次元画像をフーリエ変換すると、2次元値となるが、これを別の画像とみなすこともできる。すなわち、フーリエ面もまた画素によって構成されているとみなすことができる。そこで、一般性を失わずにNxM画素の画像に対してξの定義域を[ −N/2,
N/2   ）、ηの定義域は[
−M/2, M/2 ）とする。このとき原点（ξ,η）=（０,０）がDC成分画素に対応する。ξ方向のナイキスト周波数はN/2となる。η方向のナイキスト周波数はM/2となる。なお、典型的にはM=Nである。このとき、fnyquist = N/2となる。
【００５５】
［被写体］
フラクタルブラウンモデルにもとづく物体のスペクトルは
【数３２】
2009187378.tif 000035
とかける。qはフラクタル次元であり、2.4.9
節で設定される。
ここで、比例定数c は被写体に依存する。c を定量化する1 つの方法は、被写体に関する典型的な、適当に大きな値、例えば105 と設定することである。別の方法は、Oiso の最大値が1 となるよう正規化することである。これは被写体に依存しない設計であり、本明細ではこの方法を採用する。Oiso の最大値は
【数３３】
2009187378.tif 000036
で与えられるが、周波数領域における周波数の最小値は0（DC
成分）であり、このときフラクタルモデルが定義できない。これは、fmin = 1 とすることで解決する。1 は画素によるサンプリング周波数である。f < fmin = 1 なる非常に小さな周波数はf = fmin = 1と仮定する。以上によりOiso を次のように再定義する。
【数３４】
2009187378.tif 000037
図１５Ｂに、q=1 の場合のスペクトルを示す。
【００５６】
［回折］
回折によるスペクトルは口径D、焦点距離fl を用いて次式で表わされる。
【数３５】
2009187378.tif 000038
ここで、
【数３６】
2009187378.tif 000039
は、コヒーレント照明におけるカットオフ周波数である。光学系仕様が既知である場合を除き、f0は求まらない。この場合の設定方法は2.4.9 節に示すが、ここで示す方法に従ってf0を設定した場合のスペクトルを図１５Ｂに示す。空間周波数の増加に伴い緩やかにゲインが落ちてゆくが、ゼロ点はない。
【００５７】
［デフォーカス］
式(2.21) をフーリエ変換し、
【数３７】
2009187378.tif 000040
はデフォーカスによるスペクトルを表わす。ここで、J1は１次の第一種ベッセル関数である。256×256[pixel]の画像に対してrl=5pixel としたときのスペクトルを図１５Ｂに示す。ゲインの低下は被写体や回折の影響よりも大きいことがわかる。また、空間周波数の増加に伴い急激にゲインが落ち、周期的なゼロ点が同心円状に現れていることがわかる。第一暗環の半径r1stDarkは次式で与えられる。
【数３８】
2009187378.tif 000041
実際に生じているデフォーカスδf が未知の場合の設定方法は2.4.9
節に示す。
【００５８】
［球面収差］
本実施形態では、球面収差の効果は、デフォーカスに含める。これは、定量的に程度が把握できないためである。しかしこれは実用上問題ない。収差の分類に用いられるZernike多項式においては両者は独立であるものの、実際のデフォーカスは、Zernike多項式の球面収差項とデフォーカス項の線形和で表わされるためである。
【００５９】
［検出器による等方性ノイズ］
μspec noiseはスペクトルノイズ振幅の平均値である。画像取得モデルにおける加法ノイズw(x，y) が平均０、標準偏差σωの正規白色雑音であると仮定すれば、そのスペクトルノイズの振幅は
【数３９】
2009187378.tif 000042
とかける。図１５Ｂにこれを示す。
【００６０】
2.4.4   Jiso：ブレ角度推定評価値の導出
ブレのPSF のスペクトルHblurは式(2.5)で計算されるH(ξ,η)を、DC 成分画素を中心として角度θ だけ回転させたものとする。ブレおよび等方性ノイズHiso を考慮した仮想的な画像のスペクトルGiso は、等方性ノイズの各パラメータ（上付き文字で示した）が与えられれば、ブレ長さρの関数として次式で与えられる。
ただし右辺に現れる全ての項を定量化することは必須ではない。必須なのはHblurだけであり、すなわちGiso
= Hblurが、最小限の構成となる。ただし、より適切な評価関数を設計する目的で、Hδfの定量化が望ましい。また、数値計算の条件を向上する目的でμspec noiseの定量化が望ましい。また、更に適切な評価関数を設計する目的で、被写体の影響Oiso、および／あるいは回折の影響Hdiffの定量化が望ましい。また、右辺に現れる項以外の項目を追加することは妨げない。
【数４０】
2009187378.tif 000043
このスペクトルに対し、テンプレートとの相関が直ちに計算でき
【数４１】
2009187378.tif 000044
となる。ブレ角度θ∈(−90deg，+90deg]と、ブレ長さρ> 0に対して、Riso(θ, ρ)を数値的に求めれば、Gisoの角度推定評価値として
【数４２】
2009187378.tif 000045
が得られる。
ただし、式(2.29) における各要素のパラメータの設定によっては、Giso(f, ρ) < 1
となることがあり、このとき、式(2.30) によって定義されるRiso(θ,ρ) は負値をとる。このような場合には式(2.40) によって定義される、仮想的なブレ角度推定評価値Jisoは、ブレ角度真値において最小値を取ってしまい都合が悪い。このような場合は、Jiso を次のように、一般性を失うことなく再定義する。
例えば、負号をつけて再定義した、
【数４３】
2009187378.tif 000046
や、逆数を取って再定義した、
【数４４】
2009187378.tif 000047
が採用できる。
以降の議論は、特に断りがない限り、式(2.31) と(2.31-1)と(2.31-2) のいずれに対しても、成立する。最終的に得られるテンプレート評価値JT は、これらいずれについても同一となる。
【００６１】
2.4.5   JT：異方性ノイズに対するロバスト指標の定義
前節で導入したJisoは、ブレと等方性ノイズのみが存在する仮想的な観測モデルにもとづく角度推定評価値であった。本節では、テンプレートパラメータ設計評価値JT を、異方性ノイズに対するロバスト指標として定義する。
【数４５】
2009187378.tif 000048
ここで、funcは匿名関数であり、具体的な関数形を2.4.6 節で与える。このマージンを、異方性ノイズに対する感度で除したものを評価値とする。感度とは、異方性ノイズδnが存在するときの角度推定評価値Jの摂動δJと定義する。つまり、
【数４６】
2009187378.tif 000049
とし、この感度δJを用いてテンプレートの評価値を
【数４７】
2009187378.tif 000050
と定義する。パラメータは上付き添え字で表わす。例えば
【数４８】
2009187378.tif 000051
は、デフォーカスパラメータδf =0、検出器ノイズパラメータσw
=1 である条件のもとで構築したモデルにもとづく評価値JTであることを表わす。
【００６２】
2.4.6 marginJ ：等方性モデルにおける余裕の導出
Jisoの最大値と最小値の差が大きいほど、余裕が大きいと考えられる。そこで、大局的振幅を利用する方法が考えられる。すなわち、marginJを、その最大値と最小値の差として定義する方法である。
【数４９】
2009187378.tif 000052
この方法は直感的で分かりやすいが、ブレ角度真値近傍での挙動が考慮されていない点で、適切な定義とはいえない。例えば、図１６Ａに示す2つのJisoは、式(2.36)により定義されたmarginJは同一である。ここで、異方性ノイズが混入した場合を考えると、ある角度位置において局所的なピークが現れる。同じノイズに対して、真値近傍がなだらかな場合（同図左）はθ（ハット）= θtrueとならず、推定に失敗している。一方、真値近傍での勾配が急峻な場合（同図右）にはθ（ハット）= θtrue であり、正しく推定できている。これは、後者の方が異方性ノイズにロバストであることを示しており、marginJ を導出する目的で注目する角度範囲を、真値近傍に限定すべきことを示唆している。
【００６３】
そこで、これを具現化したmarginJの定義を与える。Jisoのブレ角度真値近傍区間Sθにおける最大値と、区間端における値の差をmarginJと定義する。これを図１６Ｂに示す。正しく推定できている同図右ではmarginJが大きく、推定に失敗している同図左では、marginJが小さいことが確認でき、異方性ノイズに対する余裕として適切な定義であると考えられる。
以上を数式で表わせば、
【数５０】
2009187378.tif 000053
第一項はブレ角度真値における角度評価値であり、多くの場合Jiso(ρ, θ = θtrue) に一致する。また、第二項は真値からδθs だけ離れた角度における評価値であり、多くの場合minθ∈Sθ Jiso(ρ, θ) に一致する。しかし、デフォーカスが大きい場合や、Tmax
が大きい場合には、これらが一致しない場合もある。δθs の決定方法は2.4.9
節に記す。
さらに、近傍Sθの端部における、Jisoの勾配を表わす指標として再定義したmarginJは、有力な典型例として利用できる。具体的には、式(2.38)において、Jisoの最大値を調べる範囲を、真値近傍に、より限定して再定義する。すなわち、
【数５１】
2009187378.tif 000054
近傍Sθmは、近傍Sθに比べてより狭い範囲として定義される。典型的には、
【数５２】
2009187378.tif 000055
が利用できる。
【００６４】
2.4.7 δJ：異方性ノイズに対する感度δJ の導出
本節では、テンプレートの異方性ノイズに対する感度δJ を見積もる。marginJ
が同一の場合、δJ が小さいほど異方性ノイズの影響を受けにくいと期待できる。
【００６５】
まず、モデル化を進めるために異方性ノイズに関して次の仮定をおく。
仮定１：異方性ノイズは特定の角度においてのみ正値をとる。そのフーリエスペクトルは、DC 成分画素を通る、幅無限小の直線である。
この仮定にもとづき、角度θ = θ0 に生じる異方性ノイズの角度依存性を図１７左に示す。また、フーリエスペクトル振幅を同図右に示す。
【００６６】
ブレによって空間周波数が減衰する方向、すなわち周波数領域においてθ = θtrue
+90deg 方向に異方性ノイズが加わる場合を想定する。テンプレートと、このノイズの相関値をδR(> 0) とおく。
まず、評価値Jiso(θ, ρ) が式(2.42) で定義される場合を考える。このとき、ノイズを想定しない場合の評価値Jiso(θtrue,ρ)noiseless からの変化量をΔJ(ρ) とすれば、
【数５３】
2009187378.tif 000056
とかける。２行目から３行目の式変形でδR << |Riso(θtrue +90deg, ρ)noiseless|
を仮定した。３行目から４行目の式変形でRiso(θtrue
+90deg, ρ)noiseless . Riso(θtrue,ρ) を仮定した。
式を整理すると、
【数５４】
2009187378.tif 000057
と表現できる。
一方、Riso(θ, ρ) < 0 なる場合には、評価値Jiso(θ, ρ) が式(2.32) に従い再定義された。この場合には、式(2.40) においてRiso(θtrue +90deg, ρ)noiseless
をRiso(θtrue
+90deg, ρ)noiseless と読み替えれば、式(2.41-1) の代わりに
【数５５】
2009187378.tif 000058
が成立する。
式(2.41-1),(2.41-2) のいずれにおいても、摂動ΔJ(ρ) は負であることがわかっているから、減少量をδJ(ρ) とおけば、
【数５６】
2009187378.tif 000059
となる。これは式(2.43),(2.44) のいずれの定義においても成立する。
δJ(ρ) は、ブレ真値角度における評価値の、異方性ノイズδR による減少量であるから、異方性ノイズに対する感度を表わす指標として採用できる。
【００６７】
［１本帯テンプレート(Tmax，0) の場合］
先の仮定１の下では、δRはTmaxに依らない定数となるため、定数をc´を用いて
【数５７】
2009187378.tif 000060
とかける。テンプレート間でδJの比較することが目的であるため、定数c´の絶対値は任意の正の実数としてよい。
【００６８】
［２本帯テンプレート(Tmax，Tmin) の場合］
２本帯テンプレート(Tmax，Tmin)の場合のδJの導出は困難である。なぜならば、δJが、異方性ノイズのフーリエスペクトル分布に強く依存するためである。異方性ノイズのスペクトル分布を定量化することは、一般に容易ではない。本明細書2.5 節では代替手法として、δJを用いずに２本帯テンプレートを設計する手法を提案する。
【００６９】
2.4.8   JTの導出
式(2.43)を式(2.34)に代入すれば、１本帯テンプレート(Tmax，0) の評価値JT は、式(2.38) で導入したmarginJを用いて
【数５８】
2009187378.tif 000061
と計算できる（c =1/c´とした）。この式にもとづいて数値的に求めたJTを図１８に示す。式(2.43)
における比例定数c =1 とし、パラメータδf,
σw は適当に与えた。
【００７０】
状態量推定においては、求めたJTを最大化するTmaxを求める際、ブレ長さに不確定性があることを考慮しなければならない。以下の２つの基準が採用できる。
【００７１】
best average基準評価値の期待値を最大化するテンプレート
【数５９】
2009187378.tif 000062
ただし、E{f}で、関数fの確率的期待値を表わすとする。
minmax基準信頼区間S 内の評価値最小値を最大化するテンプレート
【数６０】
2009187378.tif 000063
実際の設計結果は6.2.3 節に示す。
【００７２】
2.4.9 実際のテンプレート設計のためのパラメータ設定方法
2.4.3 節では、JTを導出するためのモデルにおけるノイズ要素を定式化したが、各々のパラメータは既知としていた。提案するJTにもとづいてブレ角度推定テンプレートを設計するには、これらのパラメータを実際に設定しなければならない。同様に2.4.6 節では、marginJを定量化するために必要なパラメータを既知としていた。本節では、これら未確定パラメータの設定方法を述べる。表2.2 に、これらパラメータを整理して示す。なお、上付き添え字dgで、本発明による設計値を表わすことにする。例えば、パラメータδfdgで、本発明により設計された、デフォーカスパラメータδfであることを表わす。
【表２】
2009187378.tif 000064
【００７３】
［marginJを定義する近傍角度δθsの設定方法］
ブレ角度真値からδθだけずれたテンプレートTρ,θで、長さ推定する場合を考える。縞模様と平行でないテンプレートで縞模様を検出することになるため、δθを0から大きくしていくと、やがて長さ推定が正常に行なえなくなる。この限界値をδθs とする。このような角度は、図１９に示す角度δθとして見積もることができる。すなわち、N×Npixelの画像に対し
【数６１】
2009187378.tif 000065
と計算できる。
【００７４】
さて、式(2.47) は明らかに、δθがブレ長さρに依存することを示している。しかし、設計パラメータθsは、ブレ長さによらず一定値としたい。その理由は、異方性ノイズがこの近傍角度内部に生じる影響をブレ長さによらず一定とするためである。異方性ノイズに対する感度δJは、その導出の過程で異方性ノイズの分布を仮定しないものであったため、近傍角度についてはブレ長さによらず一定であることを前提としていた。以上から、設計点を特定のブレ長さに固定する。最も短いブレρminを設計点に定めた。すなわち、角度δθsを、
【数６２】
2009187378.tif 000066
と定めた。こうして決められた角度は、テンプレートTρ,θのブレ長さに関する全定義域を包含する。
【００７５】
［フラクタル次元q設定方法］
物体の影響は周波数の高域側での落ち込みが少なく、加法ノイズに埋もれるため、デフォーカス量の影響に比べるとJTへの寄与が小さい。フラクタル次元は画像に依存するが、これを予測することは困難であるため、よく利用される値としてq =1を用いる。
【００７６】
［回折に関するパラメータf0の設定方法］
式(2.24) によって回折によるスペクトルを求めるには、光学系の口径D、焦点距離fl、観測波長λが既知でなければならない。通常、このようなパラメータは利用できる。柔軟構造光学系においても、設計仕様から直ちに求めることができる。
しかし、もし利用できない場合、例えば画像のみ与えられ、光学系に関する先験知識がない場合はどうすればよいだろうか。この場合は、式(2.27)で与えられるPSFの第一暗環がちょうど1画素に収まると仮定する。宇宙機搭載の光学系や多くのデジタルカメラにおいては、実際にこのように設計される場合が多い。この設計は、ピクセルサイズ（１画素の物理的な寸法）をpとすれば、
【数６３】
2009187378.tif 000067
とかかれる。一方、長さp間隔でサンプリングされたデータのナイキスト周波数はfnyquist = 1/(2p)であるから、pを消去して
【数６４】
2009187378.tif 000068
が得られる。
このような設計になっていない光学系も存在するが、後述するデフォーカスの影響に比べると軽微であるため、回折の効果を厳密に設定する必要はない。
【００７７】
［デフォーカスパラメータrlの設定方法］
光学系がたわむので、デフォーカスはゼロとはならないと考えるのが妥当である。しかし、定量的に把握するのは現実的ではない。
ひとつの方針として、確率的期待値を利用した方法が考えられる。熱・光学・構造設計データからデフォーカスの確率分布がモデル化できる場合、これにもとづいてJTの確率的期待値をとる。しかし、そのような確率分布も、利用できる場合は限定的と思われる。
デフォーカス量に関する確率分布を必要としない方法として、推定するブレの帯域に着目した方法を提案する。デフォーカスによるパワースペクトルは、その第一暗環の内側に83.8%が集中していることが知られており、これより外側、すなわち高周波帯域での情報をブレ推定に積極的に利用するのは、SNRの観点からは適切ではない。そこで、第一暗環より内側のスペクトル振幅だけでブレ推定が行なえるような、最大のデフォーカス量を設定値とする。具体的には、次の手続きに従えばよい。
【００７８】
Step1) 推定しようとする最短ブレ長さρminを決める。
Step2) ρminに対応する縞模様間隔を計算する。具体的には
【数６５】
2009187378.tif 000069
となる。νmaxは、テンプレートTρ,θの縞模様間隔の最大値となる。
Step3) νmaxを第一暗環の半径とするデフォーカスを求め、これをJTへの設定値とする式(2.27)と式(2.51)を等しいとおけば、直ちに
【数６６】
2009187378.tif 000070
が得られる。これが、設定すべきデフォーカスパラメータである。
以上の手続きは、ブレ角度推定が条件よく行なえる範囲で、最大のデフォーカスを想定していることが特徴である。
rdglが、実際にどの程度の焦点はずれに対応するかは、次のように確認できる。式(2.52) と式(2(l)
.20) とからrlを消去すれば、式(2.52) に対応する焦点はずれは
【数６７】
2009187378.tif 000071
であることがわかる。
実際に生じうるデフォーカスについては注目していないため、デフォーカスの真値が式(2.53) とは異なる値になる場合がある。6.2.2 節ではこの影響を、数値計算により検証する。
【００７９】
［検出器ノイズσw の設定方法］
検出器ノイズは、検出器の仕様や、想定される光量から適当に求める。例えば、検出器が8bit であればSNR が256 = 28を超えることはなく、これが上界を与える。あまり厳密に設定する必要はないが、ゼロと設定することは許容されない。その定量的な理由は、σw =0 のとき、式(2.30) より明らかなように、デフォーカスやブレに由来する縞模様のゼロ点（暗線）位置において、対数の中身がゼロとなるためである。画素値の最大値で正規化した検出器ノイズσwの設定例として、量子化ノイズのみを考える。画像が単色n bit階調であれば、量子化bit幅の半幅を1σとすればσw = 0.5 / (2n)、あるいは、3σとれば
σw
= 0.5 / (2n)/3となる。
【００８０】
2.5 角度推定の有効性をさらに高めるための２本帯テンプレート設計手法
2.4.7 節において、Tmin>0なるテンプレートについては、その評価値JTの定量化が困難であることを示した。これは、異方性ノイズのスペクトル分布の普遍的な定量化が不可能であること、分布の様子（例えばカットオフ周波数）に依存して、推定性能が大幅に異なることに起因した。しかし一方で、特定の収差に対してはTmin > 0なるテンプレートを利用することで、良好な推定結果を残すことができることが分かっている。
以上から、ブレ角度真値の見逃しが減ることを期待して、１本帯テンプレートと２本帯テンプレートを1つずつ用いてそれぞれブレ角度推定候補を求める手法を提案する。２本帯テンプレートに求められる性質は次の２つである。
(1) JTにもとづいて設計した１本帯テンプレート(Tdgmax,0)と、なるべく帯域が異なり、定性能が異なることを期待できる。
(2) 単独での推定性能が何らかの意味で保証されている。
(1) は、Tminを最大化することで達成できる。すなわち、
【数６８】
2009187378.tif 000072
を設計値とする場合である（このとき、幅W =1 となる）。しかしこのテンプレートはmarginJ < 0となる可能性がある。つまり、ブレ角度真値近傍区間で評価値が最大とならないことを意味し、正常な推定結果をもたらさず、用件(2) を満たさない。そこで、ブレ長さ予測値の信頼区間において、marginJ > 0 となるような、最大のTminを設計値とすることを提案する。数式的には、
【数６９】
2009187378.tif 000073
すなわち、Tmaxを、JTによる設計値Tmaxdgに固定したまま、TminをTmaxdgから減らしていったとき、ブレ長さに関する信頼区間Sρ全域においてmarginJ＞０となるようなTminを、設計値Tmindgとする。信頼区間Sρは、与えられる確率分布に対し、適当な有意水準を設定すれば定量化できる。
この手法で設計した２本帯テンプレートは、先の要件(1)(2) をともに満たすため、1本帯テンプレートと併用することで推定性能が高まると期待できる。
【００８１】
2.6 収差にロバストな画像ブレ推定手法
本章では、ブレ推定の原理と従来手法の原理を概説し、縞模様のみに着目してブレの角度と長さを同時に推定する手法が、周波数領域において縞模様を伴う収差があると利用できなくなるという問題点を指摘した。そこで、まずスペクトルの角度依存性に着目する手法でこの問題点の解決を試みた。この手法は、角度依存性を調べるために使うテンプレートTθの幅に性能が強く依存することを考慮し、ブレ角度推定モデルを用いた設計法を展開した。設計パラメータは、図１１に示した(Tmax,Tmin)であった。
【００８２】
以上、「方向性の強い収差」および「周波数領域において縞模様を伴う収差」の双方に対してロバストな画像ブレ推定手法は、以下のようにまとめられる。全体のブロック図を図２４Ａ、処理フローにおいて、テンプレートTθを設計するフローを図２４Ｂに、テンプレートTθを用いてブレ推定を行なうフローを図２４Ｃに示す。
【００８３】
［テンプレートTθを設計するフロー］
DC成分画素を含むテンプレートTθを設計する作業は、具体的には、テンプレートの重みwaを規定するna個のパラメータの組pa_1, pa_2, ..,pa_naを決定する作業である。
DC成分画素を含まないテンプレートTθを設計する作業は、具体的には、テンプレートの重みwbを規定するnb個のパラメータの組pb_1, pb_2, ..,pb_nbを決定する作業である。
【００８４】
step1：Gisoの計算
式(2.29)に従って、Gisoを、空間周波数fと、ブレの長さρの関数値として数値的に求める。Gisoの定量化に必要なパラメータは、2.4.9節に記載した方法に従って設定する。具体的には、[表2]
に示すフラクタル次元q、回折のパラメータf0、デフォーカスのパラメータδfdg、検出器ノイズのパラメータσwである。なお、ρの定義域をρmin以上ρmaxとし、この範囲で適当な間隔毎にGisoを求める。ここで、間隔は典型的には1画素とする。ρminは典型的には1画素とする。ρmaxは典型的には、NxN画素の画像に対してN/4程度とする。それぞれ、ブレ長さを推定するテンプレートの定義域と一致させる必要はない。
【００８５】
step2：Risoの計算
式(2.30)に従ってRisoを、ブレ角度θと、ブレ長さρの関数値として数値的に求める。式中に現れるテンプレートTθは、パラメータpa_1,pa_2,..,pa_naが未確定のままである。これらna個のパラメータを、適当な定義域内部を適当な間隔で変化させ、そのたびに、対応するRisoを数値的に求める。すなわち、ここでRiso（θ,ρ;pa_1,pa_2,..,pa_na）が得られる。パラメータの定義域は、パラメータが唯一つの場合の典型的な例として、最小値を1画素、最大値をNxN画素の画像に対してN/2程度とする。パラメータを変化させる間隔は、パラメータが唯一つの場合の典型的な例として、1画素とする。パラメータが複数ある場合は、まず定義域と間隔を適当に定め、設計結果が多峰性を示すなど、詳細を把握する必要がある場合に、間隔をより狭めて再設計する方法や、設計結果が定義域周辺部において、一定値に収束してない場合に、定義域を広げて再設計する方法が例示される。
【００８６】
step3A：δJの計算
δJの計算は、DC成分画素を含むテンプレートを設計する場合にのみ実施する。
概要を図１６Aに、具体的には[0067][0068]に記載した方法に従い、δJを計算する。 既に数値的に得られているRiso（θ,ρ; pa_1,pa_2,..,pa_na）を式(2.43)に代入することで、δJ（ρ; pa_1, pa_2,..,pa_na）を数値的に得る。式中、c’は1として一般性を失わない。
【００８７】
step3B：Jisoの計算
Jisoの計算は、DC成分画素を含むテンプレートを設計する場合にのみ実施する。[0060]に記載の方法でJisoを計算する。Riso（θ,ρ;pa_1,pa_2,..,pa_na）を式(2.31)に代入し、Jiso（θ,ρ;pa_1,pa_2,..,pa_na）を数値的に求める。
【００８８】
step4：marginJの計算
marginJの計算は、DC成分画素を含むテンプレートを設計する場合にのみ実施する。概要を図１５Cに、具体的には[0062][0063]に記載した方法でmarginJを計算する。まず、2.4.9節に記載した方法に従い、パラメータを設定する。具体的には、近傍角度δθsを 式(2.48)に従って設定する。これと、既に得られたJiso（θ,ρ; pa_1,pa_2,..,pa_na）を式(2.38)に代入して、marginJ(ρ; pa_1,pa_2,..,pa_na)が数値的に得られる。ここで、θtrueはmarginJを、異なるテンプレート同士で比較するために導入される基準であり、例えば0として一般性を失わない。
【００８９】
step5：設計評価値JTの計算
JTの計算は、DC成分画素を含むテンプレートを設計する場合にのみ実施する。
数値的に得られているmarginJ(ρ; pa_1,pa_2,..,pa_na)と、Riso（θ,ρ;
pa_1,pa_2,..,pa_na）とを式(2.45)に代入し、JT(ρ; pa_1,pa_2,..,pa_na)を数値的に求める。
【００９０】
step6A：DC成分画素を含むテンプレートを規定するパラメータの決定
[0031]に記載したように、JT(ρ; pa_1,pa_2,..,pa_na)の、ブレ長さに関する確率的期待値を最大化するパラメータの組を(pa_1,pa_2,..,pa_na)optimalとし、これを以って、テンプレートの重みwaを規定するna個のパラメータの組とする。
パラメータの組は、確率的期待値の大きいものから順に複数選んでも良いが、典型的には最大のもの1組選ぶ。
式(2.45)は、waが唯一つのパラメータTmaxを持つ場合に、このTmaxを設計する数式を例示したものである。式(2.46)に示したように、ある信頼区間におけるJTの最悪値を最大化するようなパラメータを選んでも良い。
なお、期待値を求める際、あるいは信頼区間を定義する際、ブレの長さに関する確率分布が利用できる場合は、その既知の確率分布を、用いればよい。
一方、ブレの長さに関する確率分布が利用できない場合は、例えば正規分布を仮定し、大きな誤差分散σρ2=σρb2を設定すれば、ブレ長さの予測値に依存しないようなテンプレートを設計できる。ただ一つの設計パラメータTmaxをもつテンプレートを、256x256画素の画像に適用する場合、σρb = 45pixelとすることで、Tmaxは約11pixelとなり、ブレ長さの予測値に依存しなくなる（σρを0から大きくしていき、45pixelとなったときに依存性が無くなった）。一般にNxN画素の画像に適用する場合、
【数７０】
2009187378.tif 000074
とすれば、設計結果は、
【数７１】
2009187378.tif 000075
となる。ここで、roundは丸め関数を表わし、典型的には小数点以下を四捨五入する関数とする。
【００９１】
step6B：DC成分画素を含まないテンプレートを規定するパラメータの決定
ブレの長さに関する確率分布において、適当な信頼区間を設計し、この区間内部でmarginJ>0となるような、テンプレートのパラメータを決定する。信頼区間は、よく知られた分布、例えば正規分布においては、有意水準を与えることで一意に設計できる。有意水準は典型的には5%や10%が利用できる。典型的には、[0032]に記載した条件２）に従い、重みwbの積算値を最小とするパラメータの組を1つ選ぶ。ただし、重みwbの積算値の小さいものから順に複数選んでも良い。
【００９２】
ブレの長さに関する確率分布が利用できる場合は、その既知の確率分布を、用いればよい。
ブレの長さに関する確率分布が利用できない場合は、例えば正規分布を仮定し、大きな誤差分散σρ2=σρb2を設定すれば、ブレ長さの予測値に依存しないようなテンプレートを設計できる。２つの設計パラメータ（Tmax,Tmin）をもつテンプレートを、256x256画素の画像に適用する場合、σρb = 45pixelとすることで、Tmax=11pixel,Tmin=9pixelなる設計値が得られる。一般にNxN画素の画像に適用する場合、σρbとTmaxはstep6Aにおける場合と同じ値を用いればよい。Tminについては、
【数７２】
2009187378.tif 000076
となる。
【００９３】
［テンプレートTθを用いてブレ推定を行なうフロー］
観測画像g(x,y)と、n(≧1)個のテンプレートTθが利用できる状態を初期状態とする。テンプレートTθは一部、また全部が設計したものでも良いし、一部、また全部が別の手段によって得られたものでも良い。
【００９４】
Step1：スペクトル振幅|G(ξ,η)|の計算
観測画像g(x,y)を適切に窓処理し、2次元フーリエ変換し、結果の絶対値を計算することで、スペクトル振幅|G(ξ,η)|を得る。[0022] に記載した通りである。
Step2、Step3、Step4においては、利用するn個のテンプレートTθの各々に対し、同様の処理を実施する。
【００９５】
Step2：R（θ）の計算
利用するn個のテンプレートTθの各々に対し、同様の処理を実施する。
スペクトル振幅|G|とテンプレートTθとから、式(2.11)に従って、R(θ)を計算する。
【００９６】
Step3：J（θ）の計算と平滑化
利用するn個のテンプレートTθの各々に対し、同様の処理を実施する。
R（θ）から、式(2.12)に従って、ブレ角度推定評価値J（θ）を計算する。
このJ（θ）を[0042]と式(2.14)に従って平滑化する。典型的には、移動メディアンフィルタや移動トリム平均フィルタを利用する。フィルタサイズは、典型的には角度に換算して、半幅5degから15deg程度に対応させる。
【００９７】
Step4：ブレ角度の推定
利用するn個のテンプレートTθの各々に対し、同様の処理を実施する。
J（θ）から、[0043]と式(2.15)に従って、Nθ（≧１）個の角度候補を選定する。典型的には、Nθ＝１とすればよい。ただし、平滑化を実施した結果、複数の異なる角度θが、J（θ）に同一の、最適値を与える場合がある。このような場合には、典型的には、これらを全て選択する。ただし、任意の1個を選んでもよい。最終的に選んだ角度候補の個数を同じくNθと再定義する。得られた候補は{θ（ハット）}iとかく。ただし、i=1,2,..,Nθとする。
【００９８】
Step5：ブレ角度候補の融合
{θ（ハット）}iの和集合を、n個のテンプレートTθによる、角度候補とする。これを、添え字を省略して{θ（ハット）}とする。
【数７３】
2009187378.tif 000077
【００９９】
Step6：ブレ角度と長さ推定
{θ （ハット）}をブレ角度に関する探索範囲として、ブレ角度とブレ長さを推
定する。{θ（ハット）}がただ1 つの要素からなる場合には、ブレ長さのみ推定
する。最終的な角度の推定値は{θ（ハット）}そのものとなる。
ブレ長さは、例えば、[0028]、あるいは[0029]、あるいは[0036]に記載した方法によって推定することができる。[0036]に記載した方法を用いる場合には、スペクトル振幅|G|の縞模様を検出することによってブレの角度と長さを推定するテンプレートTρ ,θ を用いて、式(2.8)で定義される評価値Jρ ,θ を最小化するようなテンプレートTρ ,θ のパラメータに対応するブレ角度と長さを、推定値とすることができる。
なお、テンプレートTρ ,θの設計に当たっては、角度と長さの定義域を決めておく必要がある。角度については、典型的には、本発明に係るブレ角度推定テンプレートTθと同一とすればよい。長さについて、定義域を[ρmin,ρmax]とする。ρminは、６画素が例示される。これより小さくすることは可能だが、余り利点がない。ρmaxは、NxN画素の画像に対して、N/2程度以下、例えばN/3程度に設定することが例示される。
以上により、ブレ角度θ（ハット）と、ブレ長さρ（ハット）が推定値として得られる。
【０１００】
ここで提案したブレ推定手法の実装上の特徴を定性的に述べる。テンプレート評価値JTのパラメータは2.4.9 節の方法で事前に計算しておくことができる。状態量推定において逐次、演算が必要なのは、この評価値のほかには、状態量推定フィルタから得られるブレの長さ予測値と、その誤差分散のみでありその他の状態量や観測量が不要である。また、もしブレの長さ予測値と、その誤差分散情報も利用できない場合には、大きな誤差分散を設定することで、固定のテンプレートが得られるため、これを用いることができる。
【０１０１】
2.7 本発明に係る推定方法と非特許文献２の推定方法との比較シミュレーション
［実験条件］
実験は、全て数値計算によった。被写体画像を100枚生成し、これに、次に示す推定すべきブレと、推定を乱す劣化要素として次に示す収差、次に示す検出器ノイズの影響を、理論式にもとづいて与えた。デフォーカス、および回折の効果は与えていない。以上は全て、空間領域で実施した。以上により、観測画像gを100枚生成した。各観測画像について、ブレ情報の推定を行なった。推定手法は、次に示す本発明に係る推定方法と次に示す非特許文献2の推定方法とを比較した。
比較は、100枚の画像に対するブレ情報の推定結果を、角度推定値と長さ推定値とに分けて度数分布表示した。
［被写体］
図２５（Ａ）に示す山岳地帯の画像（高解像度衛星画像：4845x5333画素）をランダムな位置で100 枚切り出し、部分画像（256x256 画素、白黒8bit）を100枚生成。これを被写体とした。
［ブレ］
次のパラメータで等速直線運動によるブレを与えた
方向：各観測画像の水平軸から反時計回りに20度をなす方向（θ=20deg）
長さ：24画素（ρ=24pixel）
［収差］
規格化FRINGE Zernike多項式の、5次の係数を2πとし、その他の次数の係数を全て0とするような、非点収差を与えた。これは、かなり強い非点収差である。具体的な印加方法は、まず上述のとおり一意に定義できる波面収差をフーリエ変換することで、この収差に対応するPSFを求め、これを観測画像と畳み込んだ。
［検出器ノイズ］
SNR100とし、白色正規雑音を加えた。
【０１０２】
図２５（Ａ）から切り出した被写体画像の一例（256x256画素）を図２５（Ｂ）に示す。図２５（Ｂ）に対し、ブレ、収差、検出器ノイズを印加した観測画像（256x256画素）を図２６（Ａ）に示す。図２６（Ａ）の観測画像のスペクトル振幅を図２６（Ｂ）に示す。参考として、観測画像のスペクトル振幅（図２６（Ｂ））に含まれる、推定すべきブレのスペクトル振幅を図２７（Ａ）に示す。参考として、観測画像のスペクトル振幅（図２６（Ｂ））に含まれる、非点収差のスペクトル振幅を図２７（Ｂ）に示す。
【０１０３】
［本発明に係る推定方法］
角度推定テンプレートTθを表３に示す。
【表３】
2009187378.tif 000078
Tθの関数形を決定するのに用いたパラメータを表４に示す。
【表４】
2009187378.tif 000079
「長さあるいは角度と長さ推定テンプレートTρθ」については、非特許文献2の推定方法と同一とした。
【０１０４】
［非特許文献2の推定方法］
非特許文献2において数値計算に用いられているテンプレートと同一のものを採用した。具体的には、非特許文献2においてρp(ξ,η;ν,θ)と呼ばれている関数をTρθとして採用した。ただし、角度推定は−85deg〜＋９０degの範囲で5度間隔とした。図２８に、数値実験で用いた、「長さあるいは角度と長さ推定テンプレートTρθ」の、縞模様間隔νに関するパラメータを示す。
【０１０５】
［シミュレーション結果］
図２９に結果を示す。従来手法では収差による影響が大きく出たが、本発明に係る手法では、収差の影響が限定的である。本発明では、ブレ長さ推定に先立ってブレ角度を限定したことで縞模様を伴う収差にロバストであることが示された。
【産業上の利用可能性】
【０１０６】
本発明は、画像からブレ情報を推定することに用いられる。本発明は、通常のコンピュータで実現可能であり、特別なハードウェアや重い計算処理は必要としないので、実用化は容易である。本発明の計算処理の大半は並列処理が可能であるため、実時間処理を達成することも可能である。本発明は例えば、撮影に使う光学系と被写体との相対運動に関する情報が乏しい場合に、撮影した画像のみを情報源として、ブレ情報を推定することを目的として利用され、特に、光学系の性能が不確定な環境下における、有用性が高い。ここで、光学系の性能が不確定となる環境としては、十分な校正、及びメンテナンスが行えない環境、光学系が変質したり、力学的、熱的、化学的に鏡筒が変形するような過酷な環境、より具体的には、宇宙空間、深海、災害地など人間が立ち入ることが困難な場所、が例示される。推定したブレ情報を用いた応用例として、画像復元と運動推定が例示される。ブレだけでなく収差やデフォーカスなど、様々な劣化要素を含んだ画像から復元画像を得ることは一般に困難であるが、収差にロバストなブレ推定手法を与える本発明によれば、まずブレ情報のみを推定し、これによって復元に用いる処理の自由度を制限することによって、より条件の良い劣化要素の推定と、画像復元が可能になる。また、本発明は、運動推定における利用価値が特に高い。被写体やカメラ自身の運動推定を行う場合、実時間で計算を行うことが必須であり、ブレ情報推定誤差がゼロを平均とする正規分布であると都合がよい。本発明は、並列計算による高速化が容易であり、実時間処理が可能である。また、ブレ情報推定誤差の偏りや異常推定値が少ないため、ブレ情報推定誤差がゼロを平均とする正規分布であると仮定しても、運動推定に失敗する恐れが少ない。
【図面の簡単な説明】
【０１０７】
【図１】さまざまな収差にロバストな画像ブレ推定手法の概略を示すブロック図である。
【図２】ブレ角度を推定するテンプレートを示す。図２下図に示すように、テンプレートの幅寸法としては様々な寸法を採用し得る。
【図３】ブレに由来する縞模様に平行な縞模様を持ったテンプレートを例示する。
【図４】ぶれ情報パラメータ(ρ,θ)を極座標表示した説明図である。
【図５】「帯の本数」、「ＤＣ成分画素を含むか否か」の2つの観点によるテンプレートの分類を示す。
【図６】（Ａ）DC成分画素を含むテンプレートの分布（縦軸waが重みであり、wa(ξ=0)は０でない）である。（Ｂ）DC成分画素を含まないテンプレートの分布（縦軸wbが重みであり、wb(ξ=0)=0）である。（Ｃ）DC成分画素を含むテンプレートの分布の例（ただ1つのパラメータTmaxで重み分布waが規定される１本帯テンプレート）である。（Ｄ）DC成分画素を含まないテンプレートの分布の例（2つのパラメータTmax,Tminで重み分布wbが規定される２本帯テンプレート）である。
【図７】非特許文献２に開示されたブレ角度・長さ推定テンプレートTρ，θ(ξ，η；ρ，θ)（長さρ, 角度θのブレに対応したもの）の例（左）と、その断面（右）。
【図８】ブレと非点収差を与えた画像のフーリエスペクトル振幅（左）と、その模式的説明図（右）である。
【図９】本発明に係るブレ推定方法のフローチャートである。
【図１０】フーリエスペクトル振幅|G|と角度推定テンプレートTθの相関を計算する様子。
【図１１】R（θ）を計算する角度推定テンプレートTθの断面図である。テンプレートはフーリエ面におけるＤＣ成分画素に対して対称とする。
【図１２】平滑化前の角度推定評価値J（上）と、移動トリム平均によって平滑化した評価値Jsmoothed（下）の例。異なる３種類の画像について示した。
【図１３】一般的なブレ推定モデル。ブレに対し各種劣化要素（破線囲み）が作用した結果、観測画像が得られる。この観測画像からブレを推定するのがブレ推定過程である。
【図１４】デフォーカスによるPSF の説明図。
【図１５Ａ】角度推定テンプレートTθ設計のためのブレ角度推定モデル。等方性ノイズブロックHisoの影響を受けて観測画像スペクトルGisoが生成される。続いて、異方性ノイズブロックHanisoの影響を受け、最終的な観測画像が得られる。
【図１５Ｂ】図１５Ａにおいて、marginJの導出を説明する図である。
【図１５Ｃ】図１５Ａにおいて、δＪの導出を説明する図である。
【図１６Ａ】振幅の等しい2 種類のJiso。真値近傍での勾配変化が緩やかな場合（左）と、同変化が急峻な場合（右）に、同じ異方性ノイズを与えると、推定結果が異なる。
【図１６Ｂ】真値のδθs 近傍で定義した評価値差を比較すると、真値近傍での勾配変化が緩やかな場合に小さく（左）、同変化が急峻な場合に大きい（右）
【図１７】仮定１にもとづく、異方性ノイズの角度依存性（左）と、フーリエスペクトル（右）
【図１８】σw=6.5e.4,df=0 のブレ長さ依存性（図１１に示すテンプレートのパラメータを変えて４種類(J) 示(T) した）。
【図１９】近傍角度δθs の決定方法（空間領域における解釈（左）と、周波数領域における解釈（右）は、等価である）。
【図２０】デフォーカスパラメータrl の設定方法概念図。推定しようとする最短のブレに対応する縞模様間隔νmaxと等しいような、第一暗環半径r1stDark をもつデフォーカスとする。
【図２１】平滑化した角度推定評価値Jsmoothedの単峰性/多峰性にもとづく推定結果の信頼度推定。
【図２２】テンプレートの設計手順。まず１本帯テンプレート(Tmax, Tmin)＝（Tdgmax,０)を設計する。その後(Tmax,Tmin)=(Tdgmax ,Tdgmin)を設計する。
【図２３】ブレ角度推定テンプレートの設計方針を示すブロック図である。設計方針は、次のように拡張して考えることができる。DC成分画素を含むテンプレートはmarginJおよびδJによって決定される。DC成分画素を含まないテンプレートはmarginJによって決定される。
【図２４Ａ】収差にロバストな画像ブレ推定手法の全体を示すブロック図である。
【図２４Ｂ】収差にロバストな画像ブレ推定手法の処理フローにおいて、テンプレートTθを設計するフローチャートである。
【図２４Ｃ】収差にロバストな画像ブレ推定手法の処理フローにおいて、テンプレートTθを用いてブレ推定を行なうフローチャートである。
【図２５】（Ａ）被写体として用いた高解像度衛星画像（4845x5333画素）である。（Ｂ）図２５（Ａ）から切り出した被写体画像の一例（256x256画素）である。
【図２６】（Ａ）図２５（Ｂ）に対し、ブレ、収差、検出器ノイズを印加した観測画像（256x256画素）である。（Ｂ）図２６（Ａ）の観測画像のスペクトル振幅を示す。
【図２７】（Ａ）観測画像のスペクトル振幅（図２６（Ｂ））に含まれる、推定すべきブレのスペクトル振幅を示す。（Ｂ）観測画像のスペクトル振幅（図２６（Ｂ））に含まれる、非点収差のスペクトル振幅を示す。
【図２８】ブレ長さ推定テンプレートの明線間隔と検出しようとするブレの長さとの関係を示す図である。
【図２９】本発明に係る推定方法と非特許文献２の推定方法との比較シミュレーション結果を示す。
【図１】
2009187378.tif 000082
【図３】
2009187378.tif 000083
【図６】
2009187378.tif 000084
【図１３】
2009187378.tif 000085
【図１４】
2009187378.tif 000086
【図１７】
2009187378.tif 000087
【図２３】
2009187378.tif 000088
【図２４Ｂ】
2009187378.tif 000089
【図２４Ｃ】
2009187378.tif 000090
【図２８】
2009187378.tif 000091
【図２】
2009187378.tif 000092
【図４】
2009187378.tif 000093
【図５】
2009187378.tif 000094
【図７】
2009187378.tif 000095
【図８】
2009187378.tif 000096
【図９】
2009187378.tif 000097
【図１０】
2009187378.tif 000098
【図１１】
2009187378.tif 000099
【図１２】
2009187378.tif 000100
【図１５Ａ】
2009187378.tif 000101
【図１５Ｂ】
2009187378.tif 000102
【図１５Ｃ】
2009187378.tif 000103
【図１６Ａ】
2009187378.tif 000104
【図１６Ｂ】
2009187378.tif 000105
【図１８】
2009187378.tif 000106
【図１９】
2009187378.tif 000107
【図２０】
2009187378.tif 000108
【図２１】
2009187378.tif 000109
【図２２】
2009187378.tif 000110
【図２４Ａ】
2009187378.tif 000111
【図２５】
2009187378.tif 000112
【図２６】
2009187378.tif 000113
【図２７】
2009187378.tif 000114
【図２９】
2009187378.tif 000115