(57)【要約】【課題】  推定回数である畳み込み演算の回数を減少させることによって超解像処理の高速化を実現した超解像処理の高速化方法を提供する。【解決手段】  位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法であって、高解像度画像空間に前記複数の低解像度画像の位置合わせを行い、位置合わせされた前記複数の低解像度画像の全ての画素を、前記高解像度画像空間において不等間隔にサンプリングされた画素とする第１のステップと、前記高解像度画像空間を所定の大きさを有する複数の小領域に分割する第２のステップと、第２のステップで分割された各小領域に対して、当該小領域内の所定の代表位置での推定値を当該小領域に存在する全ての画素の推定値とする第３のステップとを有する。【選択図】  図２
2006127241.tif 000002
【特許請求の範囲】
【請求項１】
位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法であって、
高解像度画像空間に前記複数の低解像度画像の位置合わせを行い、位置合わせされた前記複数の低解像度画像の全ての画素を、前記高解像度画像空間において不等間隔にサンプリングされた画素とする第１のステップと、
前記高解像度画像空間を所定の大きさを有する複数の小領域に分割する第２のステップと、
第２のステップで分割された各小領域に対して、当該小領域内の所定の代表位置での推定値を当該小領域に存在する全ての画素の推定値とする第３のステップと、
を有することを特徴とする超解像処理の高速化方法。
【請求項２】
前記小領域は正方形の形状を有する請求項１に記載の超解像処理の高速化方法。
【請求項３】
前記正方形の一辺の大きさは、前記高解像度画像のピクセルの大きさの整数分の１である請求項２に記載の超解像処理の高速化方法。
【請求項４】
第３のステップにおける前記推定値は、カメラモデルから得られる点広がり関数と前記高解像度画像の畳み込みとして推定して得られた値である請求項１乃至請求項３のいずれかに記載の超解像処理の高速化方法。
【請求項５】
前記超解像処理に必要な評価関数は次の式で表され、
2006127241.tif 000040
ただし、Ｉは前記小領域に対する評価関数であり、また、Ｍは前記小領域に属する全ての画素の数であり、ｆｉは前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値を表し、(ｘｃ,ｙｃ)は前記小領域の代表位置を表し、
2006127241.tif 000041
は前記代表位置の推定値を表す請求項１乃至請求項４のいずれかに記載の超解像処理の高速化方法。
【請求項６】
前記超解像処理に必要な評価関数は次の式で表され、
2006127241.tif 000042
ここで、下記数式が成立し、
2006127241.tif 000043
ただし、Ｉは前記小領域に対する評価関数であり、また、Ｍは前記小領域内に属する全ての画素の数であり、ｆｉは前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値を表し、(ｘｃ,ｙｃ)は前記小領域の代表位置を表し、
2006127241.tif 000044
は前記代表位置の推定値を表し、
2006127241.tif 000045
は前記小領域内の画素の観測値の平均を表す請求項１乃至請求項４のいずれかに記載の超解像処理の高速化方法。
【請求項７】
前記超解像処理に必要な評価関数は次の式で表され、
2006127241.tif 000046
ただし、Ｉは前記小領域に対する評価関数であり、また、Ｍは前記小領域に属する全ての画素の数であり、ｆｉは前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値を表し、(ｘｃ,ｙｃ)は前記小領域の代表位置を表し、
2006127241.tif 000047
は前記代表位置の推定値を表し、ｗ(ｘｉ−ｘｃ,ｙｉ−ｙｃ)は重み関数で、前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素に対応する重みを表し、その値は位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素が前記代表位置(ｘｃ,ｙｃ)から離れるに従って小さくなるようにする請求項１乃至請求項４のいずれかに記載の超解像処理の高速化方法。
【請求項８】
前記超解像処理に必要な評価関数は次の式で表され、
2006127241.tif 000048
ここで、下記数式が成立し、
2006127241.tif 000049
ただし、Ｉは前記小領域に対する評価関数であり、また、Ｍは前記小領域内に属する全ての画素の数であり、ｆｉは前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値を表し、(ｘｃ,ｙｃ)は前記小領域の代表位置を表し、
2006127241.tif 000050
は前記代表位置の推定値を表し、ｗ(ｘｉ−ｘｃ,ｙｉ−ｙｃ)は重み関数で、前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素に対応する重みを表し、その値は位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素が前記代表位置(ｘｃ,ｙｃ)から離れるに従って小さくなるようになっており、Ｗは前記小領域内の各画素に対応する重みの合計であり、
2006127241.tif 000051
は前記小領域の重み付け平均を表す請求項１乃至請求項４のいずれかに記載の超解像処理の高速化方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化する超解像処理の高速化方法に関し、特に、再構成型超解像処理において、推定回数を低減することによって、超解像処理の高速化を実現させる超解像処理の高速化方法に関する。
【背景技術】
【０００２】
位置ずれのある複数の低解像度画像より一つの高解像度画像を推定する超解像処理に関しては、近年数多くの研究で報告されている（非特許文献１参照）。例えば、非特許文献２に開示されているＭＬ(Maximum-likelihood)法、非特許文献３に開示されているＭＡＰ(Maximum  A  Posterior)法や、非特許文献４に開示されているＰＯＣＳ(Projection  Onto  Convex  Sets)法など様々な超解像処理方法が提案されている。
【０００３】
ＭＬ法とは、高解像度画像から推定された低解像度画像の画素値と実際に観測された画素値の二乗誤差を評価関数とし、評価関数を最小化するような高解像度画像を推定画像とする方法で、つまり、最尤推定の原理に基づく超解像処理方法である。
【０００４】
また、ＭＡＰ法とは、二乗誤差に高解像度画像の確率情報を付加した評価関数を最小化するような高解像度画像を推定する方法で、つまり、高解像度画像に対するある先見情報を利用して、事後確率を最大化する最適化問題として高解像度画像を推定する超解像処理方法である。
【０００５】
そして、ＰＯＣＳ法とは、高解像度画像と低解像度画像の画素値に関して連立方程式を作成し、その方程式を逐次的に解くことにより、高解像度画像を得る超解像処理方法である。
【０００６】
上述したいずれの超解像処理方法では、まず、高解像度画像を仮定し、そして仮定した高解像度画像からカメラモデルから得られる点広がり関数（ＰＳＦ関数）に基づき、全ての低解像度画像の画素毎に、その画素値を推定し、その推定値と観測された画素値（観測値）の差が小さくなるような高解像度画像を探索するという共通の特徴を有しており、これらの超解像処理方法は、再構成型超解像処理方法と呼ばれている。
【非特許文献１】サング  シー．ピー．（Sung C.P.）・ミン  ケイ．ピー．（Min K.P.）共著、「スーパーレゾルーション  イメージ  リコンストラクション:アー  テクニカル  オーバービュー（Super-ResolutionImage Reconstruction: A Technical Overview）」,ＩＥＥＥ  シグナル  プロク．  マガジン（IEEE SignalProc. Magazine）,第26巻,第3号, p.21-36，2003年
【非特許文献２】ビー．シー．トム（B.C.Tom）・エイ．ケイ．カトサゲロス（A.K.Katsaggelos）共著、「リコンストラクション  オフ  アー  ハイレゾルーション  イメージ  バイ  シマルテイニアス  レジストレイション  レストレイション  アンド  インターポレイション  オフ  ローレゾルーション  イメージズ（Reconstruction of a high-resolution image by simultaneousregistration, restoration, and interpolation of low-resolution images）」,プロク．  ＩＥＥＥ  イント．  コンフ．  イメージ  プロセシング（Proc.IEEE Int. Conf. Image Processing）,第2巻,p.539-542,1995年
【非特許文献３】アール.アール.シュルツ（R.R.Schulz）・アール.エル.スティブンスン（R.L.Stevenson）共著、「エクストラクション  オフ  ハイレゾルーション  フレームズ  フロム  ビデオ  シーケンス（Extractionof high-resolution frames from video sequences）」,ＩＥＥＥ  トランス.  イメージ  プロセシング（IEEETrans. Image Processing）,第5巻,p.996-1011,1996年
【非特許文献４】エイチ.スターク（H.Stark）・ピー.オスコウイ（P.Oskoui）共著、「ハイ  レゾルーション  イメージ  リカバリー  フロム  イメージプレーン  アレーズ  ユーシング  コンベックス  プロジェクションズ（Highresolution image recovery from image-plane arrays, using convex projections）」,ジェー.  オプト.  ソク.  エム.  エイ（J.Opt. Soc. Am. A）,第6巻,p.1715-1726,1989年
【非特許文献５】ビー.ジトーバ（B.Zitova）・ジェー.フルッサル（J.Flusser）共著、「イメージ  レジストレーション  メソッドズ：  ア  サーベイ（Imageregistration methods: a survey）」,イメージ  アンド  ビジョン  コンピューティング（Image and VisionComputing）,第21巻,第11号,p.977-1000,2003年
【非特許文献６】ディー.ケープル（D.Capel）著、「イメージ  モザイキング  アンド  スーパーレゾリューション（ImageMosaicing and Super-resolution）」,スプリンガー（Springer）,2003年
【非特許文献７】ダブリュ.エイチ.プレス（W.H.Press）・エス.エイ.テウコルスカイ（S.A.Teukolsky）他共著、「ニューメリカル  レシピズ  イン  Ｃ＋＋（Numerical Recipes inC++）」, ケンブリッジ  ユニバーシティー  プレス（Cambridge University Press）,2002年
【非特許文献８】清水・奥富共著、「領域ベースマッチングのための２次元同時サブピクセル推定法」、信学論(D-II),第J87-D-II巻,第2号,p.554-564,2004年
【非特許文献９】清水・矢野・奥富共著、「２次元サブピクセル同時推定を拡張した画像変形Ｎパラメータ同時推定」、情報処理学会研究報告,2004-CVIM-143,第2004巻,第26号,p.81-88,2004年
【発明の開示】
【発明が解決しようとする課題】
【０００７】
再構成型超解像処理方法の共通の性質として、まず、高解像度画像の画素数の未知数があり非常に高次元の問題であり、そして、複数の低解像度画像の全ての画素について高解像度画像から低解像度画像を推定する必要がある、ということが挙げられる。
【０００８】
再構成型超解像処理方法においては、高解像度画像の画素数の未知数の次元が非常に大きいため、解析的に高解像度画像を求めることは現実的ではなく、繰返し計算により高解像度画像が推定されている。また、繰返し計算では一度の繰返しに対して、全ての低解像度画像の画素に対して推定を行う必要があり、その計算コストが大きいという問題がよく知られている。つまり、再構成型超解像処理では、計算コストが大きいので、計算コストの低減が既存の超解像処理方法の主な課題であった。
【０００９】
また、超解像処理では、推定値と観測値の二乗誤差を推定誤差の評価関数として、その最適化計算の結果として高解像度画像を推定するようにしている。よって、その最適化計算のために、二乗誤差の評価関数や評価関数の微分値を計算する必要がある。
【００１０】
上述したように、既存の超解像処理方法では、二乗誤差の評価関数やその微分値を計算するために、複数の低解像度画像の全ての画素の観測値に対応する推定値を算出する必要がある。従って、複数の低解像度画像の総画素数回の推定を行う必要がある。推定計算は、カメラモデルから得られる伝達関数に相当する点広がり関数（ＰＳＦ関数）との畳み込み演算として定式化されているが、その畳み込み演算を複数の低解像度画像の全ての画素に対して行う必要があり、一般的に、超解像処理に用いられる複数の低解像度画像の総画素数は８０００〜８０万にもなり、推定に必要な計算コストが非常に大きいという問題が有った。
【００１１】
本発明は、上述のような事情よりなされたものであり、本発明の目的は、推定回数である畳み込み演算の回数を減少させることによって超解像処理の高速化を実現した超解像処理の高速化方法を提供することにある。
【課題を解決するための手段】
【００１２】
本発明は、位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法に関し、本発明の上記目的は、高解像度画像空間に前記複数の低解像度画像の位置合わせを行い、位置合わせされた前記複数の低解像度画像の全ての画素を、前記高解像度画像空間において不等間隔にサンプリングされた画素とする第１のステップと、前記高解像度画像空間を所定の大きさを有する複数の小領域に分割する第２のステップと、第２のステップで分割された各小領域に対して、当該小領域内の所定の代表位置での推定値を当該小領域に存在する全ての画素の推定値とする第３のステップとを有することにより、或いは、前記小領域は正方形の形状を有することにより、或いは、前記正方形の一辺の大きさは、前記高解像度画像のピクセルの大きさの整数分の１であることにより、或いは、第３のステップにおける前記推定値は、カメラモデルから得られる点広がり関数と前記高解像度画像の畳み込みとして推定して得られた値であることによって効果的に達成される。
【００１３】
また、本発明の上記目的は、前記超解像処理に必要な評価関数は次の式で表され、
2006127241.tif 000003
ただし、Ｉは前記小領域に対する評価関数であり、また、Ｍは前記小領域に属する全ての画素の数であり、ｆｉは前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値を表し、(ｘｃ,ｙｃ)は前記小領域の代表位置を表し、
2006127241.tif 000004
は前記代表位置の推定値を表すことにより、或いは、前記超解像処理に必要な評価関数は次の式で表され、
2006127241.tif 000005
ここで、下記数式が成立し、
2006127241.tif 000006
ただし、Ｉは前記小領域に対する評価関数であり、また、Ｍは前記小領域内に属する全ての画素の数であり、ｆｉは前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値を表し、(ｘｃ,ｙｃ)は前記小領域の代表位置を表し、
2006127241.tif 000007
は前記代表位置の推定値を表し、
2006127241.tif 000008
は前記小領域内の画素の観測値の平均を表すことにより、或いは、前記超解像処理に必要な評価関数は次の式で表され、
2006127241.tif 000009
ただし、Ｉは前記小領域に対する評価関数であり、また、Ｍは前記小領域に属する全ての画素の数であり、ｆｉは前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値を表し、(ｘｃ,ｙｃ)は前記小領域の代表位置を表し、
2006127241.tif 000010
は前記代表位置の推定値を表し、ｗ(ｘｉ−ｘｃ,ｙｉ−ｙｃ)は重み関数で、前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素に対応する重みを表し、その値は位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素が前記代表位置(ｘｃ,ｙｃ)から離れるに従って小さくなるようにすることにより、或いは、前記超解像処理に必要な評価関数は次の式で表され、
2006127241.tif 000011
ここで、下記数式が成立し、
2006127241.tif 000012
ただし、Ｉは前記小領域に対する評価関数であり、また、Ｍは前記小領域内に属する全ての画素の数であり、ｆｉは前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値を表し、(ｘｃ,ｙｃ)は前記小領域の代表位置を表し、
2006127241.tif 000013
は前記代表位置の推定値を表し、ｗ(ｘｉ−ｘｃ,ｙｉ−ｙｃ)は重み関数で、前記小領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素に対応する重みを表し、その値は位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素が前記代表位置(ｘｃ,ｙｃ)から離れるに従って小さくなるようになっており、Ｗは前記小領域内の各画素に対応する重みの合計であり、
2006127241.tif 000014
は前記小領域の重み付け平均を表すことによって一層効果的に達成される。
【発明の効果】
【００１４】
本発明に係る超解像処理の高速化方法を用いれば、まず、高解像度画像空間に位置ずれを含む複数の低解像度画像の位置合わせを行い、位置合わせ後の観測された複数の低解像度画像を、高解像度画像空間において不等間隔にサンプリングされた画素として扱い、そして、高解像度画像空間を複数の小領域に分割し、分割された各小領域内において高解像度画像から推定される値は一定であると近似することにより、超解像処理の高速化を実現できたという優れた効果を奏する。
【００１５】
また、合成画像及び２種類の実画像に対して、本発明に係る超解像処理の高速化方法を適用することによって、超解像処理の高速化が実現されたことを確認した。
【発明を実施するための最良の形態】
【００１６】
以下、本発明を実施するための最良の形態を図面を参照して説明する。
【００１７】
先ず、本発明の詳細を説明する前に、従来の超解像処理の概要について説明する。
【００１８】
超解像処理とは、観測画像である位置ずれを含む複数の低解像度画像から、一つの高解像度画像を推定する処理である。
【００１９】
なお、ここでは、『復元したい高解像度画像の空間』を『高解像度画像空間』と称する。
【００２０】
位置ずれを含む複数の低解像度画像の位置合わせを高解像度画像空間に行った後に、全ての観測画像のデータ（つまり、全ての観測画像の画素）は、図１（Ｃ）に示されるように、高解像度画像空間における不規則間隔でサンプリングされたデータとして扱うことができる。つまり、高解像度画像空間に位置合わせされた複数の低解像度画像は、不等間隔にサンプリングされた画素と考えることができる。
【００２１】
ここで、複数の観測画像（つまり、複数の低解像度画像）の位置合わせ後のｉ番目のデータ（つまり、位置合わせ後のｉ番目の画素）の位置を座標(ｘｉ,ｙｉ)とし、その画素値をｆｉとする。つまり、位置合わせ後のｉ番目の画素（以下、単にｉ番目の画素と称する）の高解像度画像空間における観測された画素値（以下、単に観測値と称する）はｆｉである。
【００２２】
複数の観測画像の位置合わせ後のデータの個数は、高解像度画像空間における不規則間隔でサンプリングされたデータの個数となり、つまり、全ての観測画像のピクセル数（画素数）に一致する。
【００２３】
一方、高解像度画像をｈ(ｘ,ｙ)と仮定した場合に、位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の画素値は、カメラモデルから得られる点広がり関数（ＰＳＦ関数）と高解像度画像ｈ(ｘ,ｙ)の畳み込みとして推定することができる。この畳み込み演算によって得られるｉ番目の画素の推定画素値（以下、単に推定値と称する）を
2006127241.tif 000015
とする。
【００２４】
従って、再構成型超解像処理方法とは、複数の観測画像の位置合わせ後の全ての画素（総画素数はＮである）に対して、画素毎の観測値とその推定値の二乗誤差の合計が小さくなるように、高解像度画像ｈ(ｘ,ｙ)を調整する方法である。
【００２５】
複数の観測画像の位置合わせ後の全ての画素の観測値とその推定値の二乗誤差の合計を数式で表せば、下記数１のようになる。
【００２６】
【数１】
2006127241.tif 000016
ただし、Ｉは、複数の観測画像の位置合わせ後の全ての画素の観測値とその推定値との二乗誤差の総和で、つまり、再構成型超解像処理における推定誤差を評価するための評価関数である。また、Ｎは、高解像度画像空間における不規則間隔でサンプリングされたデータの個数であり、つまり、全ての観測画像（低解像度画像）の総画素数である。ｆｉは、位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値（観測データ）である。また、
2006127241.tif 000017
は、位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の推定値である。
【００２７】
上記数１から分かるように、従来の超解像処理方法では、観測値（観測データ）と推定値を比較するためには、複数の観測画像の位置合わせ後の全ての画素毎に対し、つまり、高解像度画像空間における不規則間隔でサンプリングされたデータの個数Ｎ回だけ、高解像度画像ｈ(ｘ,ｙ)から
2006127241.tif 000018
を推定する必要がある。
【００２８】
しかしながら、高解像度画像ｈ(ｘ,ｙ)から
2006127241.tif 000019
を推定する、すなわち、畳み込み演算の計算コストが大きいことがよく知られている。
【００２９】
要するに、上記数１に表される評価関数を用いる従来の超解像処理方法では、例えば、図２（Ａ）に示されたように、高解像度画像空間に複数の低解像度画像の位置合わせを行った後に、高解像度画像空間に観測された不規則間隔でサンプリングされたデータがあった場合に、前述の通り、観測データの個数回（つまり、位置合わせされた全ての低解像度画像の総画素数回）の推定計算（つまり、畳み込み演算）が必要となり、図２（Ａ）の場合では、高解像度画像空間に８個の観測データがあるので、評価関数を算出するために、計８回の畳み込み演算が必要である。
【００３０】
そこで、本発明では、推定回数である畳み込み演算の回数を減少させることにより、推定に必要な計算コストの減少を実現し、よって超解像処理の高速化を可能にした。
【００３１】
次に、本発明に係る超解像処理の高速化方法について詳細に説明する。
【００３２】
本発明に係る超解像処理の高速化方法では、高解像度画像空間に位置ずれを含む複数の低解像度画像の位置合わせを行った後の全ての画素を、高解像度画像空間における不規則間隔で分布される画素とし、その高解像度画像空間を所定の大きさを有する複数の小領域に分割し、分割される各小領域毎に、その小領域内に属する全ての画素の推定値を、その小領域内の所定の位置を有する代表点（以下、単に代表位置とも称する）の代表推定値（以下、単に代表値とも称する）で近似することにより、超解像処理に必要な推定回数を低減させることを最大の特徴としている。
【００３３】
ここで、例えば、図２（Ａ）に示されたように、高解像度画像空間に８個の観測データがある場合に、本発明に係る超解像処理の高速化方法を適用すると、例えば、図２（Ｂ）に示されるように、高解像度画像空間を４分割することができる。つまり、高解像度画像空間を所定の大きさの小領域に分割する。
【００３４】
高解像度画像空間に属する画素が異なれば、つまり、画素の座標(ｘ,ｙ)が異なれば、当然、画素の推定値も異なるが、本発明では、分割された各小領域の中ではその小領域内に属する全ての画素の推定値が変化しないと仮定し、つまり、その小領域内に属する全ての画素の推定値をその小領域内の代表点の代表値で近似する。
【００３５】
このようにすることにより、各小領域に存在する複数の観測データに対しては、つまり、各小領域に属する複数の画素に対しては、１回の推定演算を行うだけで、観測データと（近似）推定値との比較を行うことができる。要するに、一回の畳み込み演算による推定により、ある小領域内に存在する複数個の観測データの誤差計算、すなわち評価関数計算が可能となる。
【００３６】
分割されたある小領域、例えば、図２（Ｂ）の網掛けの小領域を注目領域とすれば、注目領域内に３点の観測データがあるが、本発明に係る超解像処理の高速化方法を適用すれば、１回の推定演算により、この３点の観測データの比較が行える。つまり、１回の畳み込み演算により３点の観測データ分の推定を行ったことになる。
【００３７】
一方、従来の超解像処理方法では、この３点の観測データとの比較を行うためには、観測データ数と同じ３回の推定演算が必要である。つまり、この注目領域に関しては、本発明に係る超解像処理の高速化方法を適用することにより、従来の超解像処理方法と比べて、３倍の高速化を期待することができる。
【００３８】
注目領域に関する二乗誤差の合計、つまり、注目領域に対する評価関数を数式で表すと、従来の超解像処理方法が下記数２、本発明に係る超解像処理の高速化方法が下記数３となる。
【００３９】
【数２】
2006127241.tif 000020
【００４０】
【数３】
2006127241.tif 000021
ただし、数２及び数３において、Ｉは、注目領域内の全ての画素の観測値とその推定値との二乗誤差の総和で、つまり、注目領域に対する評価関数である。また、Ｍは、注目領域内における不等間隔でサンプリングされたデータの個数であり、つまり、注目領域に属する全ての画素の数である。ｆｉは、注目領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値（観測データ）を表す。
2006127241.tif 000022
は、注目領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の推定値を表す。また、(ｘｃ,ｙｃ)は、注目領域の代表位置を表す。
2006127241.tif 000023
は、その代表位置の推定値を表す。
【００４１】
数３からも分かるように、本発明では、注目領域内に観測データの個数Ｍが変化しても、注目領域の代表位置の推定値は、一定値であるため、推定演算は１回で済む。一方、数２からも分かるように、従来の超解像処理方法では、注目領域内に観測データの個数回の推定演算、（例えば、図２（Ｂ）の網掛けの注目領域の場合には、３回の推定演算）を行わなければならない。
【００４２】
本発明に係る超解像処理の高速化方法では、超解像処理に必要な評価関数Ｉを算出する際の手順（後述するステップ１〜ステップ５）が以下のようになっている。なお、本発明では、超解像処理に必要な評価関数の微分値も、評価関数の算出手順と同じような手順で算出される。
ステップ１：
高解像度画像空間を複数の所定の大きさを有する小領域に分割する。
ステップ２：
評価関数Ｉの初期値を設定する。つまり、Ｉ＝０とする。
ステップ３：
ある小領域に関して、その小領域に属する画素がある場合に、つまり、その小領域内に観測データがある場合に、その小領域の代表点（代表位置）だけに対して畳み込み演算による推定を行う。そして、小領域に属する全ての画素に対して、その画素の観測値（観測データ）と得られた代表点の推定値との差の二乗を評価関数Ｉに加えるようにする。一方、その小領域内に観測データがない場合に、評価関数Ｉに対して何もしないことにする。
ステップ４：
次の小領域に移動する。
ステップ５：
分割された全ての小領域に対して、ステップ３とステップ４を繰り返す。
【００４３】
次に、本発明において、小領域の分割方法について説明する。
【００４４】
本発明に係る超解像処理の高速化方法では、高解像度画像空間を所定の形状を有する小領域に複数分割するようにしている。例えば、その一例として、高解像度画像空間を複数の正方形の小領域に分割することができる。
【００４５】
超解像処理における推定計算は、カメラモデルから得られるＰＳＦ(Point  Spread  Function：点広がり関数)との畳み込み演算として表現される。例えば、分割される小領域が正方形の形状を有する場合に、小領域（正方形）の一辺の大きさが、高解像度画像のピクセルの大きさの整数分の１(つまり、１／Ｌ、ここでＬは整数である)であれば、ＰＳＦカーネルの種類は、ピクセル単位での位置合わせを考慮することにより、Ｌ２であることが分かる。
【００４６】
従って、本発明では、正方形の形状を有する小領域を、高解像度画像のピクセルの大きさの整数分の１とすることが好ましい。
【００４７】
以上は、正方形の形状を有する小領域について説明したが、本発明では、小領域が正方形という形状に限定されることがなく、例えば、高解像度画像空間を長方形などの他の幾何的な形状を有する小領域に分割することも可能である。
【００４８】
また、本発明では、下記数４で表す評価関数Ｉを用いることも可能である。
【００４９】
【数４】
2006127241.tif 000024
ここで、下記数５が成立する。
【００５０】
【数５】
2006127241.tif 000025
ただし、数４及び数５において、Ｉは注目領域に対する評価関数である。また、Ｍは、注目領域内における不等間隔でサンプリングされたデータの個数であり、つまり、注目領域に属する全ての画素の数である。ｆｉは、注目領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値（観測データ）を表す。また、(ｘｃ,ｙｃ)は注目領域の代表位置を表す。
2006127241.tif 000026
はその代表位置の推定値を表す。
2006127241.tif 000027
は注目領域内の画素の観測値の平均を表す。
【００５１】
ところで、本発明では、注目領域内（つまり、所定のある小領域内）で、その注目領域に属する全ての画素の推定値は一定であると近似しているが、その注目領域に属する画素が注目領域の代表位置から離れるに従い、その代表推定値の信頼性は小さくなる。このことを考慮して、本発明では、注目領域内の各画素の観測値とその推定値との二乗誤差に重みを付加することによって、下記数６で表す評価関数Ｉを用いることができる。
【００５２】
【数６】
2006127241.tif 000028
ただし、数６において、Ｉは注目領域に対する評価関数である。また、Ｍは、注目領域内における不等間隔でサンプリングされたデータの個数であり、つまり、注目領域に属する全ての画素の数である。ｆｉは、注目領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値（観測データ）を表す。また、(ｘｃ,ｙｃ)は、注目領域の代表位置を表す。
2006127241.tif 000029
は、その代表位置の推定値を表す。ｗ(ｘｉ−ｘｃ,ｙｉ−ｙｃ)は、重み関数で、注目領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素に対応する重みを表し、その値は位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素が代表位置(ｘｃ,ｙｃ)から離れるに従って小さくなる。
【００５３】
また、本発明では、注目領域内の各画素の二乗誤差を変形させると、つまり、上記数６の定数成分を除いて下記数７を得ることができる。つまり、本発明では、下記数７で表す評価関数Ｉを用いることも可能である。
【００５４】
【数７】
2006127241.tif 000030
ここで、下記数８及び数９が成立する。
【００５５】
【数８】
2006127241.tif 000031
【００５６】
【数９】
2006127241.tif 000032
ただし、数７、数８及び数９において、Ｉは注目領域に対する評価関数である。また、Ｍは、注目領域内における不等間隔でサンプリングされたデータの個数であり、つまり、注目領域に属する全ての画素の数である。ｆｉは、注目領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素の観測値（観測データ）を表す。また、(ｘｃ,ｙｃ)は注目領域の代表位置を表す。
2006127241.tif 000033
は、その代表位置の推定値を表す。ｗ(ｘｉ−ｘｃ,ｙｉ−ｙｃ)は、重み関数で、注目領域内の位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素に対応する重みを表し、その値は位置座標(ｘｉ,ｙｉ)を有するｉ番目の画素が代表位置(ｘｃ,ｙｃ)から離れるに従って小さくなる。Ｗは注目領域内の各画素に対応する重みの合計である。
2006127241.tif 000034
は注目領域の重み付け平均を表す。
【００５７】
また、観測画像が与えられて、そして分割される小領域を設定した場合に、Ｗおよび
2006127241.tif 000035
は一意に決まるので、超解像処理の間は変化しない。そして、数６と数７を比較した場合に、明らかに数７の計算量が小さいことが分かる。
【００５８】
上述したように、本発明に係る超解像処理の高速化方法では、繰返しごとに必要になる計算の効率化を目指していることが特徴である。
【００５９】
要するに、本発明に係る超解像処理の高速化方法では、まず、高解像度画像空間に位置ずれを含む複数の低解像度画像の位置合わせを行い、位置合わせ後の観測された複数の低解像度画像を、高解像度画像空間において不等間隔にサンプリングされた画素として扱い、次に、高解像度画像空間を複数の小領域に分割し、分割された各小領域内において高解像度画像から推定される値は一定であると近似することにより、超解像処理の高速化を図るようにしている。
【００６０】
なお、本発明に係る超解像処理の高速化方法においては、位置ずれを含む複数の低解像度画像の位置合わせのためのマッチングに関して、例えば非特許文献５に開示されている多くの方法があるので、これら既存のマッチング方法を利用することができる。
【００６１】
次に、本発明による超解像処理の顕著な高速化効果を確認するために、まず、合成画像に対して、本発明に係る超解像処理の高速化方法を適用した超解像処理の結果と、従来の超解像処理方法による超解像処理の結果について説明する。
【００６２】
従来の超解像処理方法による超解像処理の結果と、本発明に係る超解像処理の高速化方法を適用した超解像処理の結果とを比較するために、ＩＳＯ／ＤＩＳ１２６４０(ＩＳＯ４００)の基準画像の一部を真の高解像度画像と仮定し、ガウシアンＰＳＦ(σ＝０.３)を仮定して６４枚の合成画像(８０×６０)を低解像度画像とした。
【００６３】
図３（Ａ）に真値とした高解像度画像を、図３（Ｂ）に低解像度画像の一つを高解像度画像の大きさに拡大した図を、それぞれ示す。高解像度画像の倍率は３.２とし、２５６×１９２の大きさの高解像度画像をＭＡＰ法により高解像度画像を再構成した。ＭＡＰ法の拘束項には４近傍のＭＲＦを仮定し（非特許文献６参照）、拘束パラメータαは０.０５とした。また、最適化計算にはフレッチャー・リーベス（Fletcher-Reevess）の共役勾配法を利用し（非特許文献７参照）、初期画像は基準となる低解像度画像をバイキュービック補間により拡大した画像を利用している。主な計算の条件は表１にまとめる。
【００６４】
【表１】
2006127241.tif 000036
図３（Ｃ）に従来の超解像処理方法による超解像処理の結果を、図３（Ｄ）に本発明に係る超解像処理の高速化方法を適用した超解像処理の結果を、それぞれ示す。
【００６５】
まず、主観評価として、図３（Ｃ）、図３（Ｄ）ともに、図３（Ｂ）と比較して解像感が上がっていることが確認される。また、図３（Ｃ）と図３（Ｄ）における違いは見られない。
【００６６】
次に、十分なメモリが利用できる従来の超解像処理方法による超解像処理（以下、単に従来の超解像処理方法（ａ）とも称する）、十分なメモリが利用できない従来の超解像処理方法による超解像処理（以下、単に従来の超解像処理方法（ｂ）とも称する）、及び本発明に係る超解像処理の高速化方法を適用した超解像処理（以下、単に本発明の超解像処理方法（ｃ）とも称する）の３種類超解像処理方法に関して、ＰＳＦカーネルの個数、畳み込み演算の回数、計算時間を比較するために、ＰＳＦカーネルの生成に要する前処理時間、繰返し最適化に要する時間、および高解像度画像の画質を定量的に評価するための真値とのＲＭＳ誤差を、表２にそれぞれ示す。
【００６７】
【表２】
2006127241.tif 000037
また、表２には、計算速度を比較するために、本発明に係る超解像処理の高速化方法を適用した超解像処理の計算時間を１としたときの比率も併せて示してある。なお、超解像処理に計算機のＣＰＵは、２.８[ＧＨｚ](Pentium４)（登録商標）を使用した。
【００６８】
表２から計算時間の比較を行うと、従来の超解像処理方法（ａ）と比較して４４.１倍の高速化、また、従来の超解像処理方法（ｂ）と比較して２.８倍の高速化が図られていることが分かる。また、出力画像（再構成された高解像度画像）に関しても、本発明に係る超解像処理の高速化方法と従来の超解像処理方法との比較では、ＲＭＳ誤差が若干大きくなっているものの、その差は僅か０.１６であり、ほぼ変化していないことが確認された。
【００６９】
次に、本発明に係る超解像処理の高速化方法で近似を行う小領域の大きさ(つまり、１／Ｌ、ここでＬは整数である)と、再構成された高解像度画像の精度・高解像度画像の算出に必要な計算時間との比較を行うため、小領域の大きさ（１／Ｌ）をパラメータとして、同様の超解像処理を行った。そのときの出力画像（再構成された高解像度画像）の真値とのＲＭＳ誤差・高解像度画像の算出に必要な計算時間と分割される小領域の大きさとの関係を図４に示す。
【００７０】
図４から分かるように、ＲＭＳ誤差と計算時間には小領域の大きさをパラメータとしてトレードオフの関係があることが確認された。また、計算時間がＬ＝６（要するに、分割される小領域の大きさ１／Ｌ＝１／６）あたりで収束しているのは、本発明に係る超解像処理の高速化方法での畳み込み演算回数が必ず従来の超解像処理方法での畳み込み演算回数以下であるため、畳み込み演算回数の上限に達しているためである。そして、ＲＭＳ誤差に関しても、Ｌ＝６以降での変化は小さいことも分かる。
【００７１】
次に、実画像を用いて、本発明に係る超解像処理の高速化方法を適用した超解像処理の結果と、従来の超解像処理方法による超解像処理の結果について説明する。
【００７２】
Ｓｏｎｙ（登録商標）ＶＸ２０００という撮像装置（カメラ）を利用して動画像を撮影し、カラー画像を８ビット階調のグレースケールに変換した後に、表１の条件によって超解像処理を行った。なお、画像間の位置ずれが平行移動となるようにカメラを固定し、撮影対象を平行に移動させた。また、位置合わせには、非特許文献８に開示されている二次元同時マッチングを利用した。
【００７３】
図５（Ａ）に低解像度画像を高解像度画像の大きさに拡大した図を、図５（Ｂ）に従来の超解像処理方法によって推定された高解像度画像を、図５（Ｃ）に本発明の超解像処理方法によって推定された高解像度画像を、それぞれ示す。図５（Ａ）、図５（Ｂ）及び図５（Ｃ）では、それぞれ左側に全体図を、右側に一部を拡大した拡大図を示してある。
【００７４】
図５（Ａ）、図５（Ｂ）及び図５（Ｃ）を比較すると、図５（Ｂ）、図５（Ｃ）ともに、図５（Ａ）の低解像度画像に比べて解像感が上がっている。主観評価としては、図５（Ｂ）と図５（Ｃ）では、拡大図で比較しても同程度である。真値との比較はできないが、従来の超解像処理方法による出力画像と本発明の超解像処理方法による出力画像との違いを表す画像全体のＲＭＳ誤差は２.５であった。
【００７５】
また、合成画像を用いて従来の超解像処理方法及び本発明の超解像処理方法による超解像処理の比較と同様に、実画像（１）を用いて従来の超解像処理方法及び本発明の超解像処理方法による超解像処理に必要な計算時間の比較を表３にまとめる。
【００７６】
【表３】
2006127241.tif 000038
表３から分かるように、実画像（１）の場合でも、本発明の超解像処理方法は、従来の超解像処理方法（ａ）と比較して３４.１倍、また、従来の超解像処理方法（ｂ）と比較して３.２倍、高速になっている。
【００７７】
次に、手持ちのカメラで撮影した動画像に対して、同様に超解像処理を行った。観測画像（低解像度画像）を撮影した撮像装置は手持ちカメラであるので、観測画像間の位置合わせは、射影変形を仮定して位置合わせを行った。
【００７８】
位置合わせの推定には、非特許文献９に開示されている同時推定法を利用した。射影的な変形を考えた場合に、高解像度画像から観測画像への伝達関数であるＰＳＦも、射影変形に合わせて厳密には変形させる必要がある（非特許文献６参照）。しかしながら、ＰＳＦの変形を無視できる程度に、射影的な変形は小さいと仮定して、本発明の超解像処理方法により超解像処理を行った。
【００７９】
図６（Ａ）に低解像度画像を高解像度画像の大きさに拡大した図を、図６（Ｂ）に従来の超解像処理方法によって推定された高解像度画像を、図６（Ｃ）に本発明の超解像処理方法によって推定された高解像度画像を、それぞれ示す。また、図６（Ｄ）にＰＳＦの変形を考慮しながら、従来の超解像処理方法によって推定された高解像度画像を示す。なお、図６（Ａ）、図６（Ｂ）、図６（Ｃ）及び図６（Ｄ）では、それぞれ左側に全体図を、右側に一部を拡大した拡大図を示してある。
【００８０】
実画像（２）の場合も同様に、図６（Ａ）と比較して、図６（Ｂ）、図６（Ｃ）及び図６（Ｄ）は解像感が上がっている。また、図６（Ｂ）、図６（Ｃ）及び図６（Ｄ）のそれぞれの拡大図で比較しても、その違いは殆どない。
【００８１】
ＰＳＦの変形を無視した場合に、実画像（２）を用いて従来の超解像処理方法及び本発明の超解像処理方法による超解像処理に必要な計算時間の比較を表４にまとめる。
【００８２】
【表４】
2006127241.tif 000039
表４から分かるように、実画像（１）を用いた場合とほぼ同様に、本発明の超解像処理方法により、従来の超解像処理方法（ａ）と比較して３２.５倍、また、従来の超解像処理方法（ｂ）と比較して３.１倍、計算時間が短縮されていることが確認される。
【００８３】
上述したように、本発明では、複数の低解像度画像の全ての画素に対して推定を行うのではなく、高解像度画像空間を複数の小領域に分割し、分割された各小領域内で推定値が一定であるという近似を行うことで、つまり、その小領域内に含まれている画素座標位置に対しては、その小領域の代表点の推定値を近似的に小領域内に存在する全ての画素に対する推定値とみなすことで、推定回数を減少させ、よって計算コストを減少させることを最大の特徴としている。
【００８４】
要するに、本発明では、一つの小領域に対して一回の推定さえ行えば（その小領域の代表点の推定）、小領域内のピクセルの画素値と推定値を比較することが可能である。よって、推定回数の減少分が直接に超解像処理の高速化につながる。また、上記のように、合成画像及び２種類の実画像に対して、本発明に係る超解像処理の高速化方法を適用することによって、超解像処理の高速化が実現されたことを確認した。
【００８５】
なお、本発明に係る超解像処理の高速化方法を適用して、ＭＡＰ法に関しての高解像度画像からの低解像度画像の画素の推定の回数を減少させることが好ましいが、本発明に係る超解像処理の高速化方法は、ＭＡＰ法への適用に限らず、他の超解像処理方法、例えばＭＬ法、ＰＯＣＳ法に関しても適用可能である。
【図面の簡単な説明】
【００８６】
【図１】従来の超解像処理方法による超解像処理の概念を説明するための模式図である。
【図２】本発明に係る超解像処理の高速化方法を適用した超解像処理の概念を説明するための模式図である。
【図３】図３（Ａ）は真値とした高解像度画像を示す図である。図３（Ｂ）は低解像度画像の一つを高解像度画像の大きさに拡大した図である。図３（Ｃ）は従来の超解像処理方法による超解像処理の結果を示す図である。図３（Ｄ）は本発明に係る超解像処理の高速化方法を適用した超解像処理の結果を示す図である。
【図４】本発明に係る超解像処理の高速化方法を適用した超解像処理において、分割された小領域の大きさと、計算時間・ＲＭＳ誤差との関係を示す図である。
【図５】実画像（１）を用いて超解像処理結果を示す図である。図５（Ａ）は低解像度画像を高解像度画像の大きさに拡大した図である。図５（Ｂ）は従来の超解像処理方法によって推定された高解像度画像を示す図である。図５（Ｃ）は本発明の超解像処理方法によって推定された高解像度画像を示す図である。
【図６】実画像（２）を用いて超解像処理結果を示す図である。図６（Ａ）は低解像度画像を高解像度画像の大きさに拡大した図である。図６（Ｂ）は従来の超解像処理方法によって推定された高解像度画像を示す図である。図６（Ｃ）は本発明の超解像処理方法によって推定された高解像度画像を示す図である。図６（Ｄ）はＰＳＦの変形を考慮しながら従来の超解像処理方法によって推定された高解像度画像を示す図である。
【図１】
2006127241.tif 000052
【図２】
2006127241.tif 000053
【図３】
2006127241.tif 000054
【図４】
2006127241.tif 000055
【図５】
2006127241.tif 000056
【図６】
2006127241.tif 000057