(57)【要約】【課題】  再構成型超解像処理において、評価関数と評価関数の高解像度画像に対する微分の計算を高速化することによって超解像処理の高速化を実現した超解像処理の高速化方法を提供する。【解決手段】  位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法であって、評価関数と前記評価関数の前記高解像度画像に対する微分の計算を、基本的な画像演算の組合せにより計算するようにする。【選択図】  図１
2006309649.tif 000002
【特許請求の範囲】
【請求項１】
位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法であって、
評価関数と前記評価関数の前記高解像度画像に対する微分の計算を、基本的な画像演算の組合せにより計算することを特徴とする超解像処理の高速化方法。
【請求項２】
前記評価関数は、空間領域の前記高解像度画像に対して、次のように定式化され、
2006309649.tif 000116
前記評価関数の空間領域の前記高解像度画像に対する微分は、次のように定式化され、
2006309649.tif 000117
ただし、
2006309649.tif 000118
は空間領域の前記高解像度画像のベクトル表現を、
2006309649.tif 000119
はＰＳＦ画像のベクトル表現を、
2006309649.tif 000120
は平均観測画像のベクトル表現を、
2006309649.tif 000121
は重み画像のベクトル表現を、
2006309649.tif 000122
は前記高解像度画像の事前情報を表すカーネル画像のベクトル表現を、αは拘束の強さを表す拘束パラメータを、＊は畳込み積分を、
2006309649.tif 000123
は要素毎の乗算を、それぞれ表し、
二つのノルム
2006309649.tif 000124
は、下記のように定義され、
2006309649.tif 000125
ただし、￣は複素共役を、Ｔは転置を、それぞれ表し、
また、
2006309649.tif 000126
は、下記のように定義され、
2006309649.tif 000127
ただし、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換を、それぞれ表す請求項１に記載の超解像処理の高速化方法。
【請求項３】
前記評価関数は、空間領域の前記高解像度画像に対して、次のように定式化され、
2006309649.tif 000128
前記評価関数の空間領域の前記高解像度画像に対する微分は、次のように定式化され、
2006309649.tif 000129
ただし、
2006309649.tif 000130
は空間領域の前記高解像度画像のベクトル表現を、
2006309649.tif 000131
はＰＳＦ画像のベクトル表現を、
2006309649.tif 000132
のフーリエ変換を、
2006309649.tif 000133
は前記高解像度画像の事前情報を表すカーネル画像のベクトル表現を、
2006309649.tif 000134
のフーリエ変換を、
2006309649.tif 000135
のフーリエ変換を、
2006309649.tif 000136
のフーリエ変換を、
2006309649.tif 000137
は平均観測画像のベクトル表現を、
2006309649.tif 000138
は重み画像のベクトル表現を、αは拘束の強さを表す拘束パラメータを、＊は畳込み積分を、
2006309649.tif 000139
は要素毎の乗算を、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換を、それぞれ表し、
二つのノルム
2006309649.tif 000140
は、下記のように定義され、
2006309649.tif 000141
ただし、￣は複素共役を、Ｔは転置を、それぞれ表し、
また、
2006309649.tif 000142
は、下記のように定義され、
2006309649.tif 000143
ただし、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換を、それぞれ表す請求項１に記載の超解像処理の高速化方法。
【請求項４】
前記評価関数は、周波数領域の前記高解像度画像に対して、次のように定式化され、
2006309649.tif 000144
前記評価関数の周波数領域の前記高解像度画像に対する微分は、次のように定式化され、
2006309649.tif 000145
ただし、
2006309649.tif 000146
は周波数領域の前記高解像度画像のベクトル表現で、
2006309649.tif 000147
はＰＳＦ画像のベクトル表現を、
2006309649.tif 000148
のフーリエ変換を、
2006309649.tif 000149
のフーリエ変換を、
2006309649.tif 000150
は平均観測画像のベクトル表現を、
2006309649.tif 000151
は重み画像のベクトル表現を、
2006309649.tif 000152
は前記高解像度画像の事前情報を表すカーネル画像のベクトル表現を、
2006309649.tif 000153
のフーリエ変換を、
2006309649.tif 000154
のフーリエ変換を、αは拘束の強さを表す拘束パラメータを、ρはフーリエ変換による定数倍の違いを規格化する定数を、＊は畳込み積分を、
2006309649.tif 000155
は要素毎の乗算を、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換を、それぞれ表し、
二つのノルム
2006309649.tif 000156
は、下記のように定義され、
2006309649.tif 000157
ただし、￣は複素共役を、Ｔは転置を、それぞれ表し、
また、
2006309649.tif 000158
は、下記のように定義され、
2006309649.tif 000159
ただし、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換を、それぞれ表す請求項１に記載の超解像処理の高速化方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化する超解像処理の高速化方法に関し、特に、再構成型超解像処理において、超解像処理の計算コストを大幅に削減することにより、超解像処理の高速化を実現させる超解像処理の高速化方法に関する。
【背景技術】
【０００２】
位置ずれのある複数の低解像度画像より一つの高解像度画像を推定する超解像処理に関しては、近年数多くの研究で報告されている（非特許文献１参照）。例えば、非特許文献２に開示されているＭＬ(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）共著、「リコンストラクション  オフ  アー  ハイレゾルーション  イメージ  バイ  シマルテイニアス  レジストレイション  レストレイション  アンド  インターポレイション  オフ  ローレゾルーション  イメージズ（Reconstructionof a high-resolution image by simultaneous registration, restoration, andinterpolation 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年
【非特許文献５】エス.ベイカ（S.Baker）・アイ.マシューズ（I.Matthews）共著、「ルーカスーカナデ  ２０イヤーズ  オン：  ア  ユニファイング  フレームワーク（Lucas-Kanade 20 Years On: AUnifying Framework）」,インターナショナル  ジャーナル  オフ  コンピュータ  ビジョン（International Journal of ComputerVision）,第56巻,第3号,p.221-255,2004年
【発明の開示】
【発明が解決しようとする課題】
【０００７】
上述のように、再構成型超解像処理は、高解像度画像に対して定義される評価関数の最適化問題として定式化される。つまり、再構成型超解像処理は、推定された低解像度画像と観測された観測画像の二乗誤差に基づく評価関数の最適化問題に帰着される。
【０００８】
再構成型超解像処理では、未知数の非常に大きな最適化計算であるため、その評価関数の最適化には、例えば最急降下法などの繰返し計算法がよく利用される。このとき、繰返しごとに、評価関数と評価関数の高解像度画像に対する微分を計算する必要がある。しかし、この繰返し計算のコストが非常に大きい。また、繰返し計算では一度の繰返しに対して、全ての低解像度画像の画素に対して推定を行う必要があり、その計算コストも非常に大きい。
【０００９】
つまり、再構成型超解像処理では、計算コストが大きいので、計算コストの低減が大きな課題となっている。
【００１０】
この課題を解決するために、本発明の発明者らは、評価関数を計算するための画素値推定計算回数を低減させるという観点から、特願２００４−３１６１５４に開示されている『超解像処理の高速化方法』を提案している。つまり、複数の低解像度画像を位置合わせし、高解像度画像空間に離散化点を設定し、離散化点近傍に対応する画素の平均画素値を利用することを特徴としている。このとき、離散化点近傍に対応する画素の個数を重みとして考慮することにより、推定精度を落とさず高速に計算可能であることを示した。
【００１１】
つまり、本発明の発明者らは、特願２００４−３１６１５４において、少ない低解像度画素の画素値推定計算で、精度を落とさずに高解像度画像を再構成できる評価関数を利用した超解像処理の高速化手法を提案した。
【００１２】
より詳細に説明すると、特願２００４−３１６１５４に開示されている『超解像処理の高速化方法』では、レジストレーション後（位置合わせ後）の複数の低解像度画像の画素が高解像度画像空間において不等間隔にサンプリングされた画素と考えられ、高解像度画像空間を小領域に分割し、その小領域に含まれる画素の平均画素値を利用することを特徴としている。この小領域に対する評価関数は、下記数１のように表される。
【００１３】
【数１】
2006309649.tif 000003
ここで、下記数２が成立する。
【００１４】
【数２】
2006309649.tif 000004
ただし、Ｉは小領域に対する評価関数を、Ｍは小領域に含まれる画素の個数を、ｆｉは小領域内のｉ番目の画素の画素値を、(ｘｃ,ｙｃ)は小領域の代表位置を、
2006309649.tif 000005
は小領域の代表位置に対する画素値の推定値を、それぞれ表す。実際の計算では、全ての小領域に対する評価関数の和が超解像処理全体の評価関数となる。
【００１５】
そして、本発明の発明者らは、さらに、特願２００４−３１６１５４に提案されている評価関数を利用し、その評価関数と評価関数の高解像度画像に対する微分の計算を高速化させるという観点から、本発明に係る超解像処理の高速化方法を提案する。
【００１６】
要するに、本発明は、上述のような事情よりなされたものであり、本発明の目的は、再構成型超解像処理において、評価関数と評価関数の高解像度画像に対する微分の計算を高速化することによって超解像処理の高速化を実現した超解像処理の高速化方法を提供することにある。
【課題を解決するための手段】
【００１７】
本発明は、位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法に関し、本発明の上記目的は、評価関数と前記評価関数の前記高解像度画像に対する微分の計算を、基本的な画像演算の組合せにより計算することによって効果的に達成される。
【００１８】
また、本発明の上記目的は、前記評価関数は、空間領域の前記高解像度画像に対して、次のように定式化され、
2006309649.tif 000006
前記評価関数の空間領域の前記高解像度画像に対する微分は、次のように定式化され、
2006309649.tif 000007
ただし、
2006309649.tif 000008
は空間領域の前記高解像度画像のベクトル表現を、
2006309649.tif 000009
はＰＳＦ画像のベクトル表現を、
2006309649.tif 000010
は平均観測画像のベクトル表現を、
2006309649.tif 000011
は重み画像のベクトル表現を、
2006309649.tif 000012
は前記高解像度画像の事前情報を表すカーネル画像のベクトル表現を、αは拘束の強さを表す拘束パラメータを、＊は畳込み積分を、
2006309649.tif 000013
は要素毎の乗算を、それぞれ表すことによって一層効果的に達成される。
【００１９】
ここで、二つのノルム
2006309649.tif 000014
は、下記のように定義され、
2006309649.tif 000015
ただし、￣は複素共役を、Ｔは転置を、それぞれ表す。
【００２０】
また、
2006309649.tif 000016
は、下記のように定義され、
2006309649.tif 000017
ただし、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換を、それぞれ表す。
【００２１】
また、本発明の上記目的は、前記評価関数は、空間領域の前記高解像度画像に対して、次のように定式化され、
2006309649.tif 000018
前記評価関数の空間領域の前記高解像度画像に対する微分は、次のように定式化され、
2006309649.tif 000019
ただし、
2006309649.tif 000020
は空間領域の前記高解像度画像のベクトル表現を、
2006309649.tif 000021
はＰＳＦ画像のベクトル表現を、
2006309649.tif 000022
のフーリエ変換を、
2006309649.tif 000023
は前記高解像度画像の事前情報を表すカーネル画像のベクトル表現を、
2006309649.tif 000024
のフーリエ変換を、
2006309649.tif 000025
のフーリエ変換を、
2006309649.tif 000026
のフーリエ変換を、
2006309649.tif 000027
は平均観測画像のベクトル表現を、
2006309649.tif 000028
は重み画像のベクトル表現を、αは拘束の強さを表す拘束パラメータを、＊は畳込み積分を、
2006309649.tif 000029
は要素毎の乗算を、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換を、それぞれ表すことによって一層効果的に達成される。
【００２２】
また、本発明の上記目的は、前記評価関数は、周波数領域の前記高解像度画像に対して、次のように定式化され、
2006309649.tif 000030
前記評価関数の周波数領域の前記高解像度画像に対する微分は、次のように定式化され、
2006309649.tif 000031
ただし、
2006309649.tif 000032
は周波数領域の前記高解像度画像のベクトル表現で、
2006309649.tif 000033
はＰＳＦ画像のベクトル表現を、
2006309649.tif 000034
のフーリエ変換を、
2006309649.tif 000035
のフーリエ変換を、
2006309649.tif 000036
は平均観測画像のベクトル表現を、
2006309649.tif 000037
は重み画像のベクトル表現を、
2006309649.tif 000038
は前記高解像度画像の事前情報を表すカーネル画像のベクトル表現を、
2006309649.tif 000039
のフーリエ変換を、
2006309649.tif 000040
のフーリエ変換を、αは拘束の強さを表す拘束パラメータを、ρはフーリエ変換による定数倍の違いを規格化する定数を、＊は畳込み積分を、
2006309649.tif 000041
は要素毎の乗算を、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換を、それぞれ表すことによって一層効果的に達成される。
【発明の効果】
【００２３】
本発明に係る超解像処理の高速化方法を用いれば、評価関数と評価関数の高解像画像に対する微分の計算を、基本的な画像の演算の組合せにより計算することにより、超解像処理に必要な計算コストを大幅に削減し、よって超解像処理の高速化を実現することができたという優れた効果を奏する。
【００２４】
また、実画像に対して、本発明に係る超解像処理の高速化方法を適用することによって、超解像処理の高速化が実現されたことを確認した。
【発明を実施するための最良の形態】
【００２５】
以下、本発明を実施するための最良の形態を図面を参照して説明する。
【００２６】
本発明では、特願２００４−３１６１５４に開示されている評価関数を利用し、その評価関数と評価関数の高解像度画像に対する微分の計算を基本的な画像演算の組合せにより計算することにより、また、更にフーリエ変換を応用することにより、超解像処理の高速化を実現するようにしている。
【００２７】
ところで、超解像処理において、位置ずれのある複数枚の観測画像（低解像度画像）の各画素は、レジストレーション（位置合わせ）により、高解像度画像空間のある位置に対応付けられる。つまり、レジストレーションの後では、複数枚の観測画像は、高解像度空間内で不等間隔にサンプリングされた画素であると考えることができる。
【００２８】
この不等間隔にサンプリングされた画素位置（以下、単に観測画素位置とも称する）を高解像度画像の画素位置（以下、単に、高解像度画素位置とも称する）で近似することを考える。このとき、ある高解像度画素位置に近似される観測画素（観測画素位置）は、複数ある場合が考えられる。逆に、近似される観測画素（観測画素位置）が存在しない高解像度画素位置も存在する。
【００２９】
ここで、各高解像度画素位置に近似された複数の観測画素の平均画素値をそれぞれ計算することにより、１つの画像を生成することができる。この画像を平均観測画像と呼ぶことにする。平均観測画像は、その画素間隔（画素数）が高解像度画像と等しく、高解像度画像をカメラモデルから得られる点広がり関数（ＰＳＦ関数）でぼかした画像に相当する。ただし、近似される観測画素が存在しない画素位置に対しては、その画素値は定義されない。各高解像度画素位置に近似された観測画素の個数も、同様に１つの画像となる。その画像を重み画像と呼ぶことにする。
【００３０】
また、ＰＳＦも１つの画像と考えることにする。そうすると、評価関数と評価関数の高解像度画像に対する微分を計算するに際して、考慮しなければならない画像は、高解像度画像ｈ[ｉ,ｊ]、ＰＳＦ画像ｂ[ｉ,ｊ]、平均観測画像ｆ[ｉ,ｊ]、重み画像ｗ[ｉ,ｊ]の四つである。本発明では、この４種類の画像間の画像演算の組合せにより、評価関数と評価関数の高解像度画像に対する微分の計算を行うようにしている。
【００３１】
具体的には、本発明では、特願２００４−３１６１５４に開示されている評価関数において、小領域の大きさを高解像度画像の画素と等しく設定し、小領域の代表位置を画素中心とすることを前提とする。
【００３２】
このとき、超解像処理全体の評価関数は、下記数３で表すように、誤差項（数３の右辺の第１項）と拘束項（数３の右辺の第２項）とからなる。なお、下記数３で表す評価関数を『平均画素値による評価関数』とも称する。
【００３３】
【数３】
2006309649.tif 000042
ここで、
2006309649.tif 000043
は高解像度画像のベクトル表現を、Ｎは小領域の数を、
2006309649.tif 000044
はｉ番目の小領域の代表位置に対応するＰＳＦのベクトル表現を、ｗｉはｉ番目の小領域に含まれる画素の数を、
2006309649.tif 000045
はｉ番目の小領域に含まれる画素の平均画素値を、
2006309649.tif 000046
は高解像度画像の事前情報を表すカーネルを、αは拘束の強さを表す拘束パラメータを、＊は畳込み演算を、それぞれ表す。
【００３４】
なお、上述したように、小領域として高解像度画像の画素の大きさを設定しているため、小領域の数Ｎは、高解像度画像の画素の数と等しく、小領域の代表位置は高解像度画像の画素中心と等しい。また、拘束パラメータαが０の場合（α＝０）も考えられ、このとき、超解像処理全体の評価関数は、誤差項（数３の右辺の第１項）のみからなるものとする。
【００３５】
本発明に係る超解像処理の高速化方法は、数３に示すような、誤差項と拘束項とからなる評価関数に基づく再構成型超解像処理に適用される。
【００３６】
本発明では、上述した４種類の画像、即ち、高解像度画像、ＰＳＦ画像、平均観測画像、重み画像を用いて、数３で表す『平均画素値による評価関数』を再定義して得られたものを超解像処理全体の評価関数とし、再定義された評価関数及びその微分の計算方法（詳細は実施例で説明する）によって、評価関数及びその微分の計算を行えば、超解像処理の高速化を実現することができる。また、本発明では、４種類の画像で再定義された評価関数は、共通して誤差項と拘束項とから成る。
【実施例１】
【００３７】
実施例１の超解像処理の高速化方法（以下、単に、方法１とも称する）では、上述した４種類の画像を利用し、高解像度画像の空間表現（以下、単に、空間領域の高解像度画像
2006309649.tif 000047
、又は、高解像度画像
2006309649.tif 000048
とも称する）に対して、数３で表す『平均画素値による評価関数』を再定義する。また、再定義された評価関数と評価関数の高解像度画像に対する微分の計算も、全て空間領域においてなされる。
【００３８】
まず、実施例１の超解像処理の高速化方法では、上述した４種類の画像を用いて、超解像処理全体の評価関数が下記数４で、また、評価関数の高解像度画像に対する微分が下記数５で、それぞれ定式化される。
【００３９】
なお、数４の右辺の第１項を方法１の誤差項と、数４の右辺の第２項を方法１の拘束項と、数５の右辺の第１項を方法１の誤差項の高解像度画像に対する微分と、数５の右辺の第２項を方法１の拘束項の高解像度画像に対する微分と、それぞれ称する。
【００４０】
【数４】
2006309649.tif 000049
【００４１】
【数５】
2006309649.tif 000050
ここで、
2006309649.tif 000051
は高解像度画像のベクトル表現を、
2006309649.tif 000052
はＰＳＦ画像のベクトル表現を、
2006309649.tif 000053
は平均観測画像のベクトル表現を、
2006309649.tif 000054
は重み画像のベクトル表現を、
2006309649.tif 000055
は高解像度画像の事前情報を表すカーネル画像のベクトル表現を、αは拘束の強さを表す拘束パラメータを、それぞれ表す。また、＊は畳込み積分を表し、
2006309649.tif 000056
は要素毎の乗算を表す。
【００４２】
さらに、二つのノルム
2006309649.tif 000057
は、下記数６、数７のように定義される。
【００４３】
【数６】
2006309649.tif 000058
【００４４】
【数７】
2006309649.tif 000059
ここで、￣は複素共役を、Ｔは転置を、それぞれ表す。
【００４５】
また、
2006309649.tif 000060
は、下記数８、数９のように定義される。
【００４６】
【数８】
2006309649.tif 000061
【００４７】
【数９】
2006309649.tif 000062
ここで、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換を、それぞれ表す。
【００４８】
次に、方法１の誤差項及びその誤差項の高解像度画像に対する微分の計算手順をブロック図の図１に示す。また、方法１の拘束項及びその拘束項の高解像度画像に対する微分の計算手順をブロック図の図２に示す。
【００４９】
なお、図１及び図２において、＊は畳込み積分演算で、＋は要素毎の加算で、
2006309649.tif 000063
は要素毎の乗算で、Σは要素の総和を計算する演算である。
【００５０】
図１のブロック図に示されるように、方法１の誤差項及びその誤差項の高解像度画像に対する微分の計算手順は次のようになっている。
ステップ１Ａ：
高解像度画像
2006309649.tif 000064
とＰＳＦ画像
2006309649.tif 000065
との畳込み積分演算を行う。
ステップ１Ｂ：
ステップ１Ａでの畳込み積分演算結果から平均観測画像
2006309649.tif 000066
を引くように演算する。
ステップ１Ｃ：
ステップ１Ｂでの演算結果と重み画像
2006309649.tif 000067
との乗算を行う。
ステップ１Ｄ：
ステップ１Ｃでの乗算結果とＰＳＦ画像
2006309649.tif 000068
との畳込み積分演算を行うことにより、方法１の評価関数の誤差項の微分が得られる。
ステップ１Ｅ：
ステップ１Ｂでの演算結果とステップ１Ｃでの乗算結果との乗算を行う。
ステップ１Ｆ：
ステップ１Ｅでの乗算結果に基づいて要素の総和を計算する演算により、方法１の評価関数の誤差項の値が得られる。
【００５１】
また、図２のブロック図に示されるように、方法１の拘束項及びその拘束項の高解像度画像に対する微分の計算手順は次のようになっている。
ステップ１ａ：
高解像度画像
2006309649.tif 000069
と高解像度画像の事前情報を表すカーネル画像
2006309649.tif 000070
との畳込み積分演算を行う。
ステップ１ｂ：
ステップ１ａでの畳込み積分演算結果とカーネル画像
2006309649.tif 000071
との畳込み積分演算を行うことにより、方法１の評価関数の拘束項の微分が得られる。
ステップ１ｃ：
ステップ１ａでの畳込み積分演算結果とステップ１ａでの畳込み積分演算結果との乗算を行う。
ステップ１ｄ：
ステップ１ｃでの乗算結果に基づいて要素の総和を計算する演算により、方法１の評価関数の拘束項の値が得られる。
【実施例２】
【００５２】
実施例２の超解像処理の高速化方法（以下、単に、方法２とも称する）では、上述した４種類の画像を利用し、高解像度画像の空間表現に対して、数３で表す『平均画素値による評価関数』を再定義する。方法１との違いとして、方法２は、畳込み演算の計算をフーリエ変換を利用して計算する点である。
【００５３】
実施例２の超解像処理の高速化方法では、上述した４種類の画像を用いて、超解像処理全体の評価関数が下記数１０で、また、評価関数の高解像度画像に対する微分が下記数１１で、それぞれ定式化される。
【００５４】
なお、数１０の右辺の第１項を方法２の誤差項と、数１０の右辺の第２項を方法２の拘束項と、数１１の右辺の第１項を方法２の誤差項の高解像度画像に対する微分と、数１１の右辺の第２項を方法２の拘束項の高解像度画像に対する微分と、それぞれ称する。
【００５５】
【数１０】
2006309649.tif 000072
【００５６】
【数１１】
2006309649.tif 000073
ここで、
2006309649.tif 000074
のフーリエ変換を、
2006309649.tif 000075
のフーリエ変換を、
2006309649.tif 000076
のフーリエ変換を、
2006309649.tif 000077
のフーリエ変換を、それぞれ表す。また、数１０及び数１１において、他のそれぞれの記号は、方法１と同じであるため、その説明を省略する。
【００５７】
次に、方法２の誤差項及びその誤差項の高解像度画像に対する微分の計算手順をブロック図の図３に示す。また、方法２の拘束項及びその拘束項の高解像度画像に対する微分の計算手順をブロック図の図４に示す。
【００５８】
なお、図３及び図４において、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換をそれぞれ表す。また、＋は要素毎の加算で、
2006309649.tif 000078
は要素毎の乗算で、Σは要素の総和を計算する演算である。
【００５９】
図３のブロック図に示されるように、方法２の誤差項及びその誤差項の高解像度画像に対する微分の計算手順は次のようになっている。
ステップ２Ａ：
高解像度画像
2006309649.tif 000079
のフーリエ変換を行う。
ステップ２Ｂ：
ステップ２Ａでのフーリエ変換結果とＰＳＦ画像
2006309649.tif 000080
との乗算を行う。
ステップ２Ｃ：
ステップ２Ｂでの乗算結果の逆フーリエ変換を行う。
ステップ２Ｄ：
ステップ２Ｃでの逆フーリエ変換結果から平均観測画像
2006309649.tif 000081
を引くように演算する。
ステップ２Ｅ：
ステップ２Ｄでの演算結果と重み画像
2006309649.tif 000082
との乗算を行う。
ステップ２Ｆ：
ステップ２Ｅでの乗算結果のフーリエ変換を行う。
ステップ２Ｇ：
ステップ２Ｆでのフーリエ変換結果とＰＳＦ画像
2006309649.tif 000083
との乗算を行う。
ステップ２Ｈ：
ステップ２Ｇでの乗算結果の逆フーリエ変換を行うことにより、方法２の評価関数の誤差項の微分が得られる。
ステップ２Ｉ：
ステップ２Ｄでの演算結果とステップ２Ｅでの乗算結果との乗算を行う。
ステップ２Ｊ：
ステップ２Ｉでの乗算結果に基づいて要素の総和を計算する演算により、方法２の評価関数の誤差項の値が得られる。
【００６０】
また、図４のブロック図に示されるように、方法２の拘束項及びその拘束項の高解像度画像に対する微分の計算手順は次のようになっている。
ステップ２ａ：
高解像度画像
2006309649.tif 000084
のフーリエ変換を行う。
ステップ２ｂ：
ステップ２ａでのフーリエ変換結果とカーネル画像
2006309649.tif 000085
との乗算を行う。
ステップ２ｃ：
ステップ２ｂでの乗算結果の逆フーリエ変換を行う。
ステップ２ｄ：
ステップ２ｃでの逆フーリエ変換結果に対して、フーリエ変換を行う。
ステップ２ｅ：
ステップ２ｄでのフーリエ変換結果とカーネル画像
2006309649.tif 000086
との乗算を行う。
ステップ２ｆ：
ステップ２ｅでの乗算結果の逆フーリエ変換を行うことにより、方法２の評価関数の拘束項の微分が得られる。
ステップ２ｇ：
ステップ２ｃでの逆フーリエ変換結果とステップ２ｃでの逆フーリエ変換結果との乗算を行う。
ステップ２ｈ：
ステップ２ｇでの乗算結果に基づいて要素の総和を計算する演算により、方法２の評価関数の拘束項の値が得られる。
【実施例３】
【００６１】
実施例１及び実施例２では、空間領域の高解像度画像
2006309649.tif 000087
に対して、数３で表す『平均画素値による評価関数』が再定義され、空間領域の高解像度画像に対する最適化を行う方法を述べた。
【００６２】
ところで、よく知られているように、空間領域の高解像度画像
2006309649.tif 000088
とそのフーリエ変換に対応する周波数領域の高解像度画像
2006309649.tif 000089
（以下、単に、周波数領域の高解像度画像
2006309649.tif 000090
とも称する）には、一対一の関係がある。よって、周波数領域での高解像度画像を最適化することと、空間領域での高解像度画像を最適化することは、等価な処理であると言える。
【００６３】
そのため、空間領域の高解像度画像
2006309649.tif 000091
に対してではなく、周波数領域の高解像度画像
2006309649.tif 000092
に対して最適化計算を行い、高解像度画像を再構成することも可能である。このように、周波数領域での高解像度画像
2006309649.tif 000093
に対して最適化を行うことを考えると、その評価関数
2006309649.tif 000094
は下記数１２のように表される。
【００６４】
つまり、実施例３の超解像処理の高速化方法（以下、単に、方法３とも称する）では、数３で表す『平均画素値による評価関数』が、高解像度画像の周波数表現、即ち、周波数領域での高解像度画像
2006309649.tif 000095
に対して再定義される。方法１や方法２と比べて、方法３は評価関数が高解像度画像の周波数表現に対して再定義されていることが大きな特徴である。
【００６５】
実施例３の超解像処理の高速化方法では、超解像処理全体の評価関数が下記数１２で、また、評価関数の周波数領域での高解像度画像
2006309649.tif 000096
に対する微分
2006309649.tif 000097
が下記数１３で、それぞれ定式化される。
【００６６】
なお、数１２の右辺の第１項を方法３の誤差項と、数１２の右辺の第２項を方法３の拘束項と、数１３の右辺の第１項を方法３の誤差項の周波数領域での高解像度画像
2006309649.tif 000098
に対する微分と、数１３の右辺の第２項を方法３の拘束項の周波数領域での高解像度画像
2006309649.tif 000099
に対する微分と、それぞれ称する。
【００６７】
【数１２】
2006309649.tif 000100
【００６８】
【数１３】
2006309649.tif 000101
ここで、
2006309649.tif 000102
は周波数領域での高解像度画像のベクトル表現で、ρはフーリエ変換による定数倍の違いを規格化する定数である。また、数１２及び数１３において、他のそれぞれの記号は、方法１、方法２と同じであるため、その説明を省略する。
【００６９】
次に、方法３の誤差項及びその誤差項の周波数領域での高解像度画像に対する微分の計算手順をブロック図の図５に示す。また、方法３の拘束項及びその拘束項の周波数領域での高解像度画像に対する微分の計算手順をブロック図の図６に示す。
【００７０】
なお、図５及び図６において、Ｆはフーリエ変換を、Ｆ−１は逆フーリエ変換をそれぞれ表す。また、＋は要素毎の加算で、
2006309649.tif 000103
は要素毎の乗算で、Σは要素の総和を計算する演算で、Conj.は各要素の複素共役をとる演算である。
【００７１】
図５のブロック図に示されるように、方法３の誤差項及びその誤差項の周波数領域での高解像度画像
2006309649.tif 000104
に対する微分の計算手順は次のようになっている。
ステップ３Ａ：
周波数領域での高解像度画像
2006309649.tif 000105
とＰＳＦ画像
2006309649.tif 000106
との乗算を行う。
ステップ３Ｂ：
ステップ３Ａでの乗算結果の逆フーリエ変換を行う。
ステップ３Ｃ：
ステップ３Ｂでの逆フーリエ変換結果から平均観測画像
2006309649.tif 000107
を引くように演算する。
ステップ３Ｄ：
ステップ３Ｃでの演算結果と重み画像
2006309649.tif 000108
との乗算を行う。
ステップ３Ｅ：
ステップ３Ｄでの乗算結果のフーリエ変換を行う。
ステップ３Ｆ：
ステップ３Ｅでのフーリエ変換結果とＰＳＦ画像
2006309649.tif 000109
との乗算を行うことにより、方法３の評価関数の誤差項の周波数領域での高解像度画像
2006309649.tif 000110
に対する微分が得られる。
ステップ３Ｇ：
ステップ３Ｃでの演算結果とステップ３Ｄでの乗算結果との乗算を行う。
ステップ３Ｈ：
ステップ３Ｇでの乗算結果に基づいて要素の総和を計算する演算により、方法３の評価関数の誤差項の値が得られる。
【００７２】
また、図６のブロック図に示されるように、方法３の拘束項及びその拘束項の周波数領域での高解像度画像
2006309649.tif 000111
に対する微分の計算手順は次のようになっている。
ステップ３ａ：
周波数領域での高解像度画像
2006309649.tif 000112
とカーネル画像
2006309649.tif 000113
との乗算を行う。
ステップ３ｂ：
ステップ３ａでの乗算結果とカーネル画像
2006309649.tif 000114
との乗算を行うことにより、方法３の評価関数の拘束項の周波数領域での高解像度画像
2006309649.tif 000115
に対する微分が得られる。
ステップ３ｃ：
ステップ３ａでの乗算結果に対して、各要素の複素共役をとる演算を行う。
ステップ３ｄ：
ステップ３ｃでの演算結果とステップ３ａでの乗算結果との乗算を行う。
ステップ３ｅ：
ステップ３ｄでの乗算結果に基づいて要素の総和を計算する演算により、方法３の評価関数の拘束項の値が得られる。
【００７３】
以上のように、本発明の具体的な実施例について説明したが、次に、手持ちカメラにより撮影された実画像に対して、本発明に係る超解像処理の高速化方法を適用し、本発明の有効性を確認する。超解像処理には、１６枚（ベイヤー配列）の観測画像を低解像度画像として利用した。低解像度画像の６０×６０の領域を、倍率４×４に超解像処理し、２４０×２４０の高解像度画像を再構成した。ＰＳＦはガウシアンＰＳＦを利用した。拘束カーネルとして、ラプラシアン拘束カーネルを利用した。最適化計算は最急降下法を利用し、繰返し計算回数は２０回とした。レジストレーションは、手持ちカメラであることを考慮し、非特許文献５に開示されている「射影変換を仮定した勾配法」を利用した。なお、計算には、２.８[ＧＨｚ]（Pentium４）（登録商標）のＣＰＵを使用した。
【００７４】
図７には、本発明の発明者らが提案した特願２００４−３１６１５４に開示されている『超解像処理の高速化方法』を適用した超解像処理の結果と、本発明に係る超解像処理の高速化方法を適用した超解像処理の結果を示している。
【００７５】
具体的に、図７（Ａ）は、数３で表す『平均画素値による評価関数』に基づく超解像処理の結果を、図７（Ｂ）は、本発明実施例１の超解像処理の高速化方法を適用した超解像処理の結果を、図７（Ｃ）は、本発明実施例２の超解像処理の高速化方法を適用した超解像処理の結果を、図７（Ｄ）は、本発明実施例３の超解像処理の高速化方法を適用した超解像処理の結果を、それぞれ示している。
【００７６】
図７（Ａ）、（Ｂ）、（Ｃ）、（Ｄ）に示す４種類の超解像処理結果は、主観評価としてほぼ同一である。次に、超解像処理の計算時間を見てみる。図７（Ａ）の超解像処理結果を得るための計算時間は９.６２[sec]で、図７（Ｂ）の超解像処理結果を得るための計算時間は１.３８[sec]で、図７（Ｃ）の超解像処理結果を得るための計算時間は１.０１[sec]で、
図７（Ｄ）の超解像処理結果を得るための計算時間は０.５９[sec]である。
【００７７】
ところで、特願２００４−３１６１５４に開示されているように、特願２００４−３１６１５４に係る超解像処理の高速化方法を適用した超解像処理、つまり、数３で表す『平均画素値による評価関数』に基づく超解像処理は、従来の超解像処理より、高速になっていることが明らかである。上記の超解像処理の計算時間を比較すると、特願２００４−３１６１５４に開示されている方法による超解像処理より、本発明に係る超解像処理の高速化方法による超解像処理は、更なる高速化を実現したことが良く分かる。
【００７８】
つまり、『平均画素値による評価関数』に基づく超解像処理と比較して、本発明の方法１は計算時間で６.９７倍、方法２は計算時間で９.５２倍、方法３は計算時間で１６.３倍、計算が高速になっていることが確認された。
【図面の簡単な説明】
【００７９】
【図１】本発明実施例１の超解像処理の高速化方法（方法１）に使用される評価関数の誤差項及びその誤差項の高解像度画像に対する微分の計算手順を示すブロック図である。
【図２】本発明実施例１の超解像処理の高速化方法（方法１）に使用される評価関数の拘束項及びその拘束項の高解像度画像に対する微分の計算手順を示すブロック図である。
【図３】本発明実施例２の超解像処理の高速化方法（方法２）に使用される評価関数の誤差項及びその誤差項の高解像度画像に対する微分の計算手順を示すブロック図である。
【図４】本発明実施例２の超解像処理の高速化方法（方法２）に使用される評価関数の拘束項及びその拘束項の高解像度画像に対する微分の計算手順を示すブロック図である。
【図５】本発明実施例３の超解像処理の高速化方法（方法３）に使用される評価関数の誤差項及びその誤差項の周波数領域での高解像度画像に対する微分の計算手順を示すブロック図である。
【図６】本発明実施例３の超解像処理の高速化方法（方法３）に使用される評価関数の拘束項及びその拘束項の周波数領域での高解像度画像に対する微分の計算手順を示すブロック図である。
【図７】本発明に係る超解像処理の高速化方法を適用した超解像処理の結果を示す図である。図７（Ａ）は、数３で表す『平均画素値による評価関数』に基づく超解像処理の結果を、図７（Ｂ）は、本発明実施例１の超解像処理の高速化方法を適用した超解像処理の結果を、図７（Ｃ）は、本発明実施例２の超解像処理の高速化方法を適用した超解像処理の結果を、図７（Ｄ）は、本発明実施例３の超解像処理の高速化方法を適用した超解像処理の結果を、それぞれ示している。
【図１】
2006309649.tif 000160
【図２】
2006309649.tif 000161
【図３】
2006309649.tif 000162
【図４】
2006309649.tif 000163
【図５】
2006309649.tif 000164
【図６】
2006309649.tif 000165
【図７】
2006309649.tif 000166