WO2006046493A1

WO2006046493A1 - 超解像処理の高速化方法

Info

Publication number: WO2006046493A1
Application number: PCT/JP2005/019437
Authority: WO
Inventors: Masayuki Tanaka; Masatoshi Okutomi
Original assignee: Tokyo Institute Of Technology
Priority date: 2004-10-29
Filing date: 2005-10-17
Publication date: 2006-05-04
Also published as: US20090080805A1; US8009933B2; EP1811453A1; JP3837575B2; EP1811453A4; JP2006127241A

Abstract

推定回数である畳み込み演算の回数を減少させることによって超解像処理の高速化を実現した超解像処理の高速化方法を提供する。　位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法であって、高解像度画像空間に前記複数の低解像度画像の位置合わせを行い、位置合わせされた前記複数の低解像度画像の全ての画素を、前記高解像度画像空間において不等間隔にサンプリングされた画素とする第１のステップと、前記高解像度画像空間を所定の大きさを有する複数の小領域に分割する第２のステップと、第２のステップで分割された各小領域に対して、当該小領域内の所定の代表位置での推定値を当該小領域に存在する全ての画素の推定値とする第３のステップとを有する。

Description

超解像処理の高速化方法

技術分野

本発明は、複数の低解像度画明像から一つの高解像度画像を推定する超解像処理を高速化する超解像処理の高速化方法に関し、特に、再構成型超解像処理において、推定回数を低減することによって、超解像処理の書

高速化を実現させる超解像処理の高速化方法に関する。

背景技術

位置ずれのある複数の低解像度画像より一つの高解像度画像を推定する超解像処理に関しては、近年数多くの研究で報告されている（非特許文献 1参照）。例えば、非特許文献 2に開示されている M L (Maxinrnm - likelihood)法、非特許文献 3 に開示されている M A P (Maximum A Posterior)法や、非特許文献 4に開示されている P O C S (Projection Onto Convex Sets)法など様々な超解像処理方法が提案されている。

ML法とは、高解像度画像から推定された低解像度画像の画素値と実際に観測された画素値の二乗誤差を評価関数とし、評価関数を最小化するような高解像度画像を推定画像とする方法で、つまり、最尤推定の原理に基づく超解像処理方法である。

また、 MA P法とは、二乗誤差に高解像度画像の確率情報を付加した評価関数を最小化するような高解像度画像を推定する方法で、つまり、高解像度画像に対するある先見情報を利用して、事後確率を最大化する最適化問題として高解像度画像を推定する超解像処理方法である。そして、 P O C S法とは、高解像度画像と低解像度画像の画素値に関して連立方程式を作成し、その方程式を逐次的に解くことにより、高解像度画像を得る超解像処理方法である。

上述したいずれの超解像処理方法では、まず、高解像度画像を仮定し、そして仮定した高解像度画像から力メラモデルから得られる点広がり関数（ P S F関数）に基づき、全ての低解像度画像の画素毎に、その画素値を推定し、その推定値と観測された画素値（観測値）の差が小さくなるような高解像度画像を探索するという共通の特徴を有しており、これらの超解像処理方法は、再構成型超解像処理方法と呼ばれている。

再構成型超解像処理方法の共通の性質として、まず、高解像度画像の画素数の未知数があり非常に高次元の問題であり、そして、複数の低解像度画像の全ての画素について高解像度画像から低解像度画像を推定する必要がある、ということが挙げられる。

再構成型超解像、処理方法においては、高解像度画像の画素数の未知数の次元が非常に大きいため、解析的に高解像度画像を求めることは現実的ではなく、繰返し計算により高解像度画像が推定されている。また、繰返し計算では一度の繰返しに対して、全ての低解像度画像の画素に対して推定を行う必要があり、その計算コストが大きいという問題がよく知られている。つまり、再構成型超解像処理では、計算コストが大きいので、計算コストの低減が既存の超解像処理方法の主な課題であった。

また、超解像処理では、推定値と観測値の二乗誤差を推定誤差の評価関数として、その最適化計算の結果として高解像度画像を推定するようにしている。よって、その最適化計算のために、二乗誤差の評価関数や評価関数の微分値を計算する必要がある。

上述したように、既存の超解像処理方法では、二乗誤差の評価関数やその微分値を計算するために、複数の低解像度画像の全ての画素の観測値に対応する推定値を算出する必要がある。従って、複数の低解像度画像の総画素数回の推定を行う必要がある。推定計算は、カメラモデルから得られる伝達関数に相当する点広がり関数（ P S F関数）との畳み込み演算として定式化されているが、その畳み込み演算を複数の低解像度画像の全ての画素に対して行う必要があり、一般的に、超解像処理に用いられる複数の低解像度画像の総画素数は 8 0 0 0〜 8 0万にもなり、推定に必要な計算コストが非常に大きいという問題が有った。

本発明は、上述のような事情よりなされたものであり、本発明の目的は、推定回数である畳み込み演算の回数を減少させることによって超解像処理の高速化を実現した超解像処理の高速化方法を提供することにある。発明の開示

本発明は、位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法に関し、本発明の上記目的は、高解像度画像空間に前記複数の低解像度画像の位置合わせを行い、位置合わせされた前記複数の低解像度画像の全ての画素を、前記高解像度画像空間において不等間隔にサンプリングされた画素とする第 1 のステップと、前記高解像度画像空間を所定の大きさを有する複数の小領域に分割する第 2のステップと、第 2のステップで分割された各小領域に対して、当該小領域内の所定の代表位置での推定値を当該小領域に存在する全ての画素の推定値とする第 3のステップとを有することにより、或いは、前記小領域は正方形の形状を有することにより、或いは、前記正方形の一辺の大きさは、前記高解像度画像のピクセルの大きさの整数分の 1であることにより、或いは、第 3 のステップにおける前記推定値は、カメラモデルから得られる点広がり関数と前記高解像度画像の畳み込みとして推定して得られた値であることによつて効果的に達成される。

また、本発明の上記目的は、前記超解像処理に必要な評価関数は次の式で表され、

ただし、 I は前記小領域に対する評価関数であり、また、 Mは前記小領域に属する全ての画素の数であり、 f iは前記小領域内の位置座標（X ; ， y i )を有する i番目の画素の観測値を表し、（X _c ， y Jは前記小領域の代表位置を表し、 (x_e， _e)は前記代表位置の推定値を表すことにより、或いは、前記超解像処理に必要な評価関数は次の式で表され、

ここで、下記数式が成立し、

ただし、 I は前記小領域に対する評価関数であり、また、 Mは前記小領域内に属する全ての画素の数であり、 f ；は前記小領域内の位置座標（ X i ， y ； )を有する i 番目の画素の観測値を表し、（ X _c , y 。）は前記小領域の代表位置を表し、 ( , ）は前記代表位置の推定値を表し、 /は前記小領域内の画素の観測値の平均を表すことにより、或いは、前記超解像処理に必要な評価関数は次の式で表され、

ただし、 I は前記小領域に対する評価関数であり、また、 Mは前記小領域に属する全ての画素の数であり、 f iは前記小領域内の位置座標（X i ， _{y i} )を有する i 番目の画素の観測値を表し、（x ^ y jは前記小領域の代表位置を表し、は前記代表位置の推定値を表し、 —

X _c， y ；— y Jは重み関数で、前記小領域内の位置座標（X い y i )を有する i 番目の画素に対応する重みを表し、その値は位置座標（X i， y i ) を有する i番目の画素が前記代表位置（X _c， y c )から離れるに従って小さくなるようにすることにより、或いは、前記超解像処理に必要な評価関数は次の式で表され、

ここで、下記数式が成立し、

一 I ^M

} 30 ただし、 I は前記小領域に対する評価関数であり、また、 Mは前記小領域内に属する全ての画素の数であり、 f _;は前記小領域内の位置座標（

X ；， y i )を有する i 番目の画素の観測値を表し、（X c , y _c.)は前記小領域の代表位置を表し、）は前記代表位置の推定値を表し、 w (x i — X _c， y i — y Jは重み関数で、前記小領域内の位置座標（x い y を有する i 番目の画素に対応する重みを表し、その値は位置座標（X i， y Jを有する i番目の画素が前記代表位置（x _c, y _c)から離れるに従つて小さくなるようになつており、 Wは前記小領域内の各画素に対応する重みの合計であり、ァは前記小領域の重み付け平均を表すことによって一層効果的に達成される。図面の簡単な説明

第 1図は、従来の超解像処理方法による超解像処理の概念を説明するための模式図である。

第 2図は、本発明に係る超解像処理の高速化方法を適用した超解像処理の概念を説明するための模式図である。

第 3図（A) は真値とした高解像度画像を示す図である。第 3図（B ) は低解像度画像の一つを高解像度画像の大きさに拡大した図である。第 3図（C ) は従来の超解像処理方法による超解像処理の結果を示す図である。第 3図（D) は本発明に係る超解像処理の高速化方法を適用した超解像処理の結果を示す図である。

第 4図は、本発明に係る超解像処理の高速化方法を適用した超解像処理において、分割された小領域の大きさと、計算時間 * RM S誤差との関係を示す図である。

第 5図は、実画像（ 1 ) を用いて超解像処理結果を示す図である。第 5図（A) は低解像度画像を高解像度画像の大きさに拡大した図である。第 5図（ B ) は従来の超解像処理方法によって推定された高解像度画像を示す図である。第 5図（C ) は本発明の超解像処理方法によって推定された高解像度画像を示す図である。

第 6図は、実画像（ 2 ) を用いて超解像処理結果を示す図である。第 6図（A) は低解像度画像を高解像度画像の大きさに拡大した図である。第 6図（B ) は従来の超解像処理方法によって推定された高解像度画像を示す図である。第 6図（C ) は本発明の超解像処理方法によって推定された高解像度画像を示す図である。第 6図（D ) は P S Fの変形を考慮しながら従来の超解像処理方法によって推定された高解像度画像を示す図である。発明を実施するための最良の形態

以下、本発明を実施するための最良の形態を図面を参照して説明する先ず、本発明の詳細を説明する前に、従来の超解像処理の概要について説明する。

超解像処理とは、観測画像である位置ずれを含む複数の低解像度画像から、一つの高解像度画像を推定する処理である。

なお、ここでは、『復元したい高解像度画像の空間』を『高解像度画像空間』と称する。

位置ずれを含む複数の低解像度画像の位置合わせを高解像度画像空間に行った後に、全ての観測画像のデータ（つまり、全ての観測画像の画素）は、第 1 図（C ) に示されるように、高解像度画像空間における不規則間隔でサンプリングされたデータとして扱うことができる。つまり、高解像度画像空間に位置合わせされた複数の低解像度画像は、不等間隔にサンプリングされた画素と考えることができる。

ここで、複数の観測画像（つまり、複数の低解像度画像）の位置合わせ後の i番目のデータ（つまり、位置合わせ後の i番目の画素）の位置を座標（X i , y i )とし、その画素値を f i とする。つまり、位置合わせ後の i番目の画素（以下、単に i番目の画素と称する）の高解像度画像空間における観測された画素値（以下、単に観測値と称する）は f ；である。複数の観測画像の位置合わせ後のデータの個数は、高解像度画像空間における不規則間隔でサンプリングされたデータの個数となり、つまり、全ての観測画像のピクセル数（画素数）に一致する。

一方、 '高解像度画像を h (x， y )と仮定した場合に、位置座標（X i， y _; )を有する i 番目の画素の画素値は、カメラモデルから得られる点広がり関数（ P S F関数）と高解像度画像 h ( X， y )の畳み込みとして推定することができる。この畳み込み演算によって得られる i番目の画素の推定画素値（以下、単に推定値と称する）を (

とする。

従って、再構成型超解像処理方法とは、複数の観測画像の位置合わせ後の全ての画素（総画素数は Nである）に対して、画素毎の観測値とその推定値の二乗誤差の合計が小さくなるように、高解像度画像 h ( x , y )を調整する方法である。

複数の観測画像の位置合わせ後の全ての画素の観測値とその推定値の二乗誤差の合計を数式で表せば、下記数 1 のようになる。

【数 1 】

J-ム —/(½ ）

/'=1 ただし、 I は、複数の観測画像の位置合わせ後の全ての画素の観測値とその推定値との二乗誤差の総和で、つまり、再構成型超解像処理における推定誤差を評価するための評価関数である。また、 Nは、高解像度画像空間における不規則間隔でサンプリングされたデータの個数であり、つまり、全ての観測画像（低解像度画像）の総画素数である。 f ；は、位置座標（ X _; , y i )を有する i 番目の画素の観測値（観測データ）である。また、は、位置座標（X i， y i )を有する i番目の画素の推定値である。

上記数 1 から分かるように、従来の超解像処理方法では、観測値（観測データ）と推定値を比較するためには、複数の観測画像の位置合わせ後の全ての画素毎に対し、つまり、高解像度画像空間における不規則間隔でサンプリングされたデータの個数 N回だけ、高解像度画像 h ( X， y )からを推定する必要がある。

しかしながら、高解像度画像 h ( X， y )から .を推定する、すなわち

、畳み込み演算の計算コストが大きいことがよく知られている。

要するに、上記数 1 に表される評価関数を用いる従来の超解像処理方法では、例えば、第 2図（A ) に示されたように、高解像度画像空間に複数の低解像度画像の位置合わせを行った後に、高解像度画像空間に観測された不規則間隔でサンプリングされたデータがあった場合に、前述の通り、観測データの個数回（つまり、位置合わせされた全ての低解像度画像の総画素数回）の推定計算（つまり、畳み込み演算）が必要となり、第 2図（A ) の場合では、高解像度画像空間に 8個の観測データがあるので、評価関数を算出するために、計 8回の畳み込み演算が必要である。

そこで、本発明では、推定回数である畳み込み演算の回数を減少させることにより、推定に必要な計算コストの減少を実現し、よつて超解像処理の高速化を可能にした。

次に、本発明に係る超解像処理の高速化方法について詳細に説明する本発明に係る超解像処理の高速化方法では、高解像度画像空間に位置ずれを含む複数の低解像度画像の位置合わせを行った後の全ての画素を、高解像度画像空間における不規則間隔で分布される画素とし、その高解像度画像空間を所定の大きさを有する複数の小領域に分割し、分割される各小領域毎に、その小領域内に属する全ての画素の推定値を、その小領域内の所定の位置を有する代表点（以下、単に代表位置とも称する ) の代表推定値（以下、単に代表値とも称する）で近似することにより、超解像処理に必要な推定回数を低減させることを最大の特徴としている。

ここで、例えば、第 2図（A ) に示されたように、高解像度画像空間に 8個の観測データがある場合に、本発明に係る超解像処理の高速化方法を適用すると、例えば、第 2図（B ) に示されるように、高解像度画像空間を 4分割することができる。つまり、高解像度画像空間を所定の大きさの小領域に分割する。

高解像度画像空間に属する画素が異なれば、つまり、画素の座標（_x， y )が異なれば、当然、画素の推定値も異なるが、本発明では、分割された各小領域の中ではその小領域内に属する全ての画素の推定値が変化レないと仮定し、つまり、その小領域内に属する全ての画素の推定値をその小領域内の代表点の代表値で近似する。

このようにすることにより、各小領域に存在する複数の観測データに対しては、つまり、各小領域に属する複数の画素に対しては、 1 回の推定演算を行うだけで、観測データと（近似）推定値との比較を行うことができる。要するに、一回の畳み込み演算による推定により、ある小領域内に存在する複数個の観測データの誤差計算、すなわち評価関数計算が可能となる。

分割されたある小領域、例えば、第 2図（B ) の網掛けの小領域を注目領域とすれば、注目領域内に 3点の観測データがあるが、本発明に係る超解像処理の高速化方法を適用すれば、 1回の推定演算により、この 1

3点の観測データの比較が行える。つまり、 1回の畳み込み演算により 3点の観測データ分の推定を行ったことになる。

一方、従来の超解像処理方法では、この 3点の観測データとの比較を行うためには、観測データ数と同じ 3回の推定演算が必要である。つまり、この注目領域に関しては、本発明に係る超解像処理の高速化方法を適用することにより、従来の超解像処理方法と比べて、' 3倍の高速化を期待することができる。

注目領域に関する二乗誤差の合計、つまり、注目領域に対する評価関数を数式で表すと、従来の超解像処理方法が下記数 2、本発明に係る超解像処理の高速化方法が下記数 3 となる。

【数 2】

【数 3 】

ただし、数 2及ぴ数 3において、 I は、注目領域内の全ての画素の観測値とその推定値との二乗誤差の総和で、つまり、注目領域に対する評価関数である。また、 Mは、注目領域内における不等間隔でサンプリングされたデータの個数であり、つまり、注目領域に属する全ての画素の数である。 f ；は、注目領域内の位置座標（X i， y i )を有する i番目の画素の観測値（観測データ）を表す。 (X,., は、注目領域内の位置座標（X i， y i )を有する i番目の画素の推定値を表す。また、（x _e， y _c ) は、注目領域の代表位置を表す。 ( , ）は、その代表位置の推定値を表す。数 3からも分かるように、本発明では、注目領域内に観測データの個数 Mが変化しても、注目領域の代表位置の推定値は、一定値であるため、推定演算は 1回で済む。一方、数 2からも分かるように、従来の超解像処理方法では、注目領域内に観測データの個数回の推定演算、 (例えば、第 2図（ B ) の網掛けの注目領域の場合には、 3回の推定演算）を行わなければならない。

本発明に係る超解像処理の高速化方法では、超解像処理に必要な評価関数 I を算出する際の手順（後述するステップ 1〜ステップ 5 ) が以下のようになっている。なお、本発明では、超解像処理に必要な評価関数の微分値も、評価関数の算出手順と同じような手順で算出される。

ステップ 1 ：

高解像度画像空間を複数の所定の大きさを有する小領域に分割する。ステップ 2 ：

評価関数 I の初期値を設定する。つまり、 1 = 0 とする。

ステップ 3 ：

ある小領域に関して、その小領域に属する画素がある場合に、つまり、その小領域内に観測データがある場合に、その小領域の代表点（代表位置）だけに対して畳み込み演算による推定を行う。そして、小領域に属する全ての画素に対して、その画素の観測値（観測データ）と得られた代表点の推定値との差の二乗を評価関数 I に加えるようにする。一方、その小領域内に観測データがない場合に、評価関数 I に対して何もしないことにする。

ステップ 4 ：

次の小領域に移動する。

ステップ 5 ：分割された全ての小領域に対して、ステップ 3 とステップ 4を操り返す。

次に、本発明において、小領域の分割方法について説明する。

本発明に係る超解像処理の高速化方法では、高解像度画像空間を所定の形状を有する小領域に複数分割するようにしている。例えば、その一例として、高解像度画像空間を複数の正方形の小領域に分割することができる。

超解像処理における推定計算は、カメラモデルから得られる P S F (Po int Spread Funct i on：点広がり関数）との畳み込み演算として表現される。例えば、分割される小領域が正方形の形状を有する場合に、小領域（正方形）の一辺の大きさが、高解像度画像のピクセルの大きさの整数分の 1 (つまり、 1 Z L、ここで Lは整数である）であれば、 P S Fカーネルの種類は、ピクセル単位での位置合わせを考慮することにより、 L ²であるこ'とが分かる。

従って、本発明では、正方形の形状を有する小領域を、高解像度画像のピクセルの大きさの整数分の 1 とすることが好ましい。

以上は、正方形の形状を有する小領域について説明したが、本発明では、小領域が正方形という形状に限定されることがなく、例えば、高解像度画像空間を長方形などの他の幾何的な形状を有する小領域に分割することも可能である。

また、本発明では、下記数 4で表す評価関数 I を用いることも可能でめ O

【数 4】

I = M\f -f(x_c,y 4

ここで、下記数 5が成立する

【数 5 】

ただし、数 4及び数 5において、 I は注目領域に対する評価関数である。また、 Mは、注目領域内における不等間隔でサンプリングされたデータの個数であり、つまり、注目領域に属する全ての画素の数である。 f iは、注目領域内の位置座標（X i， y i )を有する i 番目の画素の観測値（観測データ）を表す。また、（X y Jは注目領域の代表位置を表す。はその代表位置の推定値を表す。 _ は注目領域内の画素の観測値の平均を表す。

ところで、本発明では、注目領域内（つまり、所定のある小領域内）で、その注目領域に属する全ての画素の推定値は一定であると近似しているが、その注目領域に属する画素が注目領域の代表位置から離れるに従い、その代表推定値の信頼性は小さくなる。このことを考慮して、本発明では、注目領域内の各画素の観測値とその推定値との二乗誤差に重みを付加することによって、下記数 6で表す評価関数 I を用いることができる。

【数 6】

M _Λ

= ―^, — Λ ) — f(x_c， y_c )

=1 ただし、数 6において、 I は注目領域に対する評価関数である。また、 Mは、注目領域内における不等間隔でサンプリングされたデータの個数であり、つまり、注目領域に属する全ての画素の数である。 f ；は、注目領域内の位置座標（X ；， y i )を有する i 番目の画素の観測値（観測データ）を表す。また、（x ^ y jは、注目領域の代表位置を表す。 (x_c，y_c)は、その代表位置の推定値を表す。 w i — x _c ' y i — y jは、重み関数で、注目領域内の位置座標（X i , y i )を有する i 番目の画素に対応する重みを表し、その値は位置座標（X i i y Jを有する i 番目の画素が代表位置（X _c， y Jから離れるに従って小さくなる。

また、本発明では、注目領域内の各画素の二乗誤差を変形させると、つまり、上記数 6の定数成分を除いて下記数 7を得ることができる。つまり、本発明では、下記数 7で表す評価関数 I を用いることも可能である。

【数 7】

i = w\f-f(x_c,y_c)\²

ここで、下記数 8及ぴ数 9が成立する。

【数 8】

【数 9】

一 1 ^Μ

ゾ ^:=^ 2^^ (^— ， ― c) .

/=1

ただし、数 7、数 8及ぴ数 9において、 I は注目領域に対する評価関数である。また、 Μは、注目領域内における不等間隔でサンプリングされたデータの個数であり、つまり、注目領域に属する全ての画素の数である。 f iは、注目領域内の位置座標（_{X i}， y i )を有する i 番目の画素の観測値（観測データ）を表す。また、（X c , y Jは注目領域の代表位置を表す。， _c)は、その代表位置の推定値を表す。 w i x i — x ^ y i 一 y _c )は、重み関数で、注目領域内の位置座標（x ^ y i )を有する i 番目の画素に対応する重みを表し、その値は位置座標（X い y を有する i番目の画素が代表位置（X y _c )から離れるに従って小さくなる。 Wは注目領域内の各画素に対応する重みの合計である。 /は注目領域の重み付け平均を表す。

また、観測画像が与えられて、そして分割される小領域を設定した場合に、 Wおよび/は一意に決まるので、超解像処理の間は変化しない。そして、数 6 と数 7を比較した場合に、明らかに数 7の計算量が小さいことが分かる。

上述したように、本発明に係る超解像処理の高速化方法では、繰返しごとに必要になる計算の効率化を目指していることが特徴である。

要するに、本発明に係る超解像処理の高速化方法では、まず、高解像度画像空間に位置ずれを含む複数の低解像度画像の位置合わせを行い、位置合わせ後の観測された複数の低解像度画像を、高解像度画像空間において不等間隔にサンプリングされた画素として扱い、次に、高解像度画像空間を複数の小領域に分割し、分割された各小領域内において高解像度画像から推定される値は一定であると近似することにより、超解像処理の高速化を図るようにしている。

なお、本発明に係る超解像処理の高速化方法においては、位置ずれを含む複数の低解像度画像の位置合わせのためのマッチングに関して、例えば非特許文献 5に開示されている多くの方法があるので、これら既存のマッチング方法を利用することができる。

次に、本発明による超解像処理の顕著な高速化効果を確認するために、まず、合成画像に対して、本発明に係る超解像処理の高速化方法を適用した超解像処理の結果と、従来の超解像処理方法による超解像処理の結果について説明する。 7

従来の超解像処理方法による超解像処理の結果と、本発明に係る超解像処理の高速化方法を適用した超解像処理の結果とを比較するために、 I S OZD I S 1 2 6 4 0 ( I S O 4 0 0 )の基準画像の一部を真の高解像度画像と仮定し、ガウシアン P S F ( _CT = 0. 3 )を仮定して 6 4枚の合成画像（8 0 X 6 0 )を低解像度画像とした。

第 3図（A) に真値とした高解像度画像を、第 3図（B ) に低解像度画像の一つを高解像度画像の大きさに拡大した図を、それぞれ示す。高解像度画像の倍率は 3. 2 とし、 2 5 6 X 1 9 2の大きさの高解像度画像を MA P法により高解像度画像を再構成した。 MA P法の拘束項には 4近傍の MR Fを仮定し（非特許文献 6参照）、拘束パラメータ αは 0 . 0 5 とした。また、最適化計算にはフレッチャー · リーべス（ Fletcher- Reevess) の共役勾配法を利用し（非特許文献 7参照）、初期画像は基準となる低解像度画像をバイキュービック補間により拡大した画像を利用している。主な計算の条件は表 1 にまとめる。

【表 1 】計算条件

第 3図（C ) に従来の超解像処理方法による超解像処理の結果を、第 3図（D) に本発明に係る超解像処理の高速化方法を適用した超解像処理の結果を、それぞれ示す。 8

まず、主観評価として、第 3図（C) 、第 3図（D ) ともに、第 3図 ( B ) と比較して解像感が上がっていることが確認される。また、第 3 図（C) と第 3図（D) における違いは見られない。

次に、十分なメモリが利用できる従来の超解像処理方法による超解像処理（以下、単に従来の超解像処理方法（ a ) とも称する）、十分なメモリが利用できない従来の超解像処理方法による超解像処理（以下、単に従来の超解像処理方法（ b ) とも称する）、及び本発明に係る超解像処理の高速化方法を適用した超解像処理（以下、単に本発明の超解像処理方法（ c ) とも称する）の 3種類超解像処理方法に関して、 P S F力一ネルの個数、畳み込み演算の回数、計算時間を比較するために、 P S Fカーネルの生成に要する前処理時間、繰返し最適化に要する時間、おょぴ高解像度画像の画質を定量的に評価するための真値との RM S誤差を、表 2にそれぞれ示す。

【表 2】従来の超解像処理方法と本発明の超解像処理方法との比較

(合成画像を用いた場合）

また、表 2には、計算速度を比較するために、本発明に係る超解像処理の高速化方法を適用した超解像処理の計算時間を 1 としたときの比率も併せて示してある。なお、超解像処理に計算機の C P Uは、 2. 8 [G H z ] (Pentium4 ) (登録商標）を使用した。

表 2から計算時間の比較を行うと、従来の超解像処理方法（ a ) と比較して 4 4. 1倍の高速化、また、従来の超解像処理方法（ b ) と比較して 2. 8倍の高速化が図られていることが分かる。また、出力画像（ S構成された高解像度画像）に関しても、本発明に係る超解像処理の高速化方法と従来の超解像処理方法との比較では、 RM S誤差が若干大きくなつているものの、その差は僅か 0. 1 6 であり、ほぼ変化していないことが確認された。

次に、本発明に係る超解像処理の高速化方法で近似を行う小領域の大きさ（つまり、 1 ZL、ここで Lは整数である）と、再構成された高解像度画像の精度 ' 高解像度画像の算出に必要な計算時間との'比較を行うため、小領域の大きさ（ 1 / L) をパラメータとして、同様の超解像処理を行った。そのときの出力画像（再構成された高解像度画像）の真値との RM S誤差 · 高解像度画像の算出に必要な計算時間と分割される小領域の大きさとの関係を第 4図に示す。

第 4図から分かるように、 R M S誤差と計算時間には小領域の大きさをパラメータとしてトレードオフの関係があることが確認された。また、計算時間が L = 6 (要するに、分割される小領域の大きさ 1 ZL == 1 / 6 ) あたりで収束しているのは、本発明に係る超解像処理の高速化方法での畳み込み演算回数が必ず従来の超解像処理方法での畳み込み演算回数以下であるため、畳み込み演算回数の上限に達しているためである。そして、 RM S誤差に関しても、 L = 6以降での変化は小さレ、ことも分かる。

次に、実画像を用いて、本発明に係る超解像処理の高速化方法を適用した超解像処理の結果と、従来の超解像処理方法による超解像処理の結果について説明する。

S o n y (登録商標） V X 2 0 0 0 という撮像装置（カメラ）を利用して動画像を撮影し、カラー画像を 8 ビット階調のグレースケールに変換した後に、表 1 の条件によって超解像処理を行った。なお、画像間の位置ずれが平行移動となるようにカメラを固定し、撮影対象を平行に移動させた。また、位置合わせには、非特許文献 8に開示されている二次元同時マッチングを利用した。

第 5図（A) に低解像度画像を高解像度画像の大きさに拡大した図を、第 5図（B ) に従来の超解像処理方法によって推定された高解像度画像を、第 5図（C ) に本発明の超解像処理方法によって推定された高解像度画像を、それぞれ示す。第 5図（A) 、第 5図（B ) 及び第 5図（ C ) では、それぞれ左側に全体図を、右側に一部を拡大した拡大図を示してある。

第 5図（A) 、第 5図（B ) 及び第 5図（ C) を比較すると、第 5図 ( B ) 、第 5図（ C) ともに、第 5図（A) の低解像度画像に比べて解像感が上がっている。主観評価としては、第 5図（B ) 'と第 5図（C ) では、拡大図で比較しても同程度である。真値との比較はできないが、従来の超解像処理方法による出力画像と本発明の超解像処理方法による出力画像との違いを表す画像全体の RM S誤差は 2. 5であった。

また、合成画像を用いて従来の超解像処理方法及ぴ本発明の超解像処理方法による超解像処理の比較と同様に、実画像（ 1 ) を用いて従来の超解像処理方法及び本発明の超解像処理方法による超解像処理に必要な計算時間の比較を表 3にまとめる。

【表 3】 2

従来の超解像処理方法と本発明の超角^ [象処理方法との比較

f実画像（1)を用いた場合）

表 3から分かるように、実画像（ 1 ) の場合でも、本発明の超解像処理方法は、従来の超解像処理方法（ a ) と比較して 3 4 . 1倍、また、従来の超解像処理方法（ b ) と比較して 3 . 2倍、高速になっている。

次に、手持ちのカメラで撮影した動画像に対して、同様に超解像処理を行った。観測画像（低解像度画像）を撮影した撮像装置は手持ちカメラであるので、観測画像間の位置合わせは、射影変形を仮定して位置合わせを行った。

位置合わせの推定には、非特許文献 9に開示されている同時推定法を利用した。射影的な変形を考えた場合に、高解像度画像から観測画像への伝達関数である P S Fも、射影変形に合わせて厳密には変形させる必要がある（非特許文献 6参照）。しかしながら、 P S Fの変形を無視できる程度に、射影的な変形は小さいと仮定して、本発明の超解像処理方法により超解像処理を行った。

第 6図（A ) に低解像度画像を高解像度画像の大きさに拡大した図を、第 6図（B ) に従来の超解像処理方法によって推定された高解像度画像を、第 6図（C ) に本発明の超解像処理方法によって推定された高解像度画像を、それぞれ示す。また、第 6図（D ) に P S Fの変形を考慮しながら、従来の超解像処理方法によって推定された高解像度画像を示す。なお、第 6図（A) 、第 6図（B ) 、第 6図（C ) 及び第 6図（D ) では、それぞれ左側に全体図を、右側に一部を拡大した拡大図を示してある。

実画像（ 2 ) の場合も同様に、第 6図（A) と比較して、第 6図（B ) 、第 6図（C ) 及び第 6図（D) は解像感が上がっている。また、第 6図（B ) 、第 6図（C) 及び第 6図（D) のそれぞれの拡大図で比較しても、その違いは殆どない。

P S Fの変形を無視した場合に、実画像（ 2 ) を用いて従来の超解像処理方法及び本発明の超解像処理方法による超解像処理に必要な計算時間の比較を表 4にまとめる。

【表 4】従来の超解像処理方法と本発明の超解像処理方法との比較

(実画像（2)を用いた場合）

表 4から分かるように、実画像（ 1 ) を用いた場合とほぼ同様に、本発明の超解像処理方法により、従来の超解像処理方法（ a ) と比較して 3 2. 5倍、また、従来の超解像処理方法（ b ) と比較して 3. 1倍、計算時間が短縮されていることが確認される。

上述したように、本発明では、複数の低解像度画像の全ての画素に対して推定を行うのではなく、高解像度画像空間を複数の小領域に分割し、分割された各小領域内で推定値が一定であるという近似を行うことで、つまり、その小領域内に含まれている画素座標位置に対しては、その小領域の代表点の推定値を近似的に小領域内に存在する全ての画素に対する推定値とみなすことで、推定回数を減少させ、よって計算コストを減少させることを最大の特徴としている。

要するに、本発明では、一つの小領域に対して一回の推定さえ行えば (その小領域の代表点の推定）、小領域内のピクセルの画素値と推定値を比較することが可能である。よって、推定回数の減少分が直接に超解像処理の高速化につながる。また、上記のように、合成画像及び 2種類の実画像に対して、本発明に係る超解像処理の高速化方法を適用することによって、超解像処理の高速化が実現されたことを確認した。

なお、本発明に係る超解像処理の高速化方法を適用して、 M A P法に関しての高解像度画像からの低解像度画像の画素の推定の回数を減少させることが好ましいが、本発明に係る超解像処理の高速化方法は、 M A P法への適用に限らず、他の超解像処理方法、例えば M L法、 P O C S 法に関しても適用可能である。産業上の利用可能性

以上のように、本発明に係る超解像処理の高速化方法を用いれば、まず、高解像度画像空間に位置ずれを含む複数の低解像度画像の位置合わせを行い、位置合わせ後の観測された複数の低解像度画像を、高解像度画像空間において不等間隔にサンプリングされた画素として扱い、そして、高解像度画像空間を複数の小領域に分割し、分割された各小領域内において高解像度画像から推定される値は一定であると近似することにより、超解像処理の高速化を実現できたという優れた効果を奏する。また、合成画像及び 2種類の実画像に対して、本発明に係る超解像処理の高速化方法を適用することによって、超解像処理の高速化が実現されたことを確認した。

ぐ参考文献一覧 >

非特許文献 1 ：

サングシー . ピー. (Sung C. P. ) ' ミンケィ . ピー . （Min K.P.

) 共著、「スーパーレゾルーシヨンイメージリコンストラクション

：ァーテク二カノレーノ一ビュー Super-Resolution Image Reconstruction: A Technical Overviewリ」， I E E E シグナノレプロタ . マガジン (IEEE Signal Proc. Magazine) ，第 26 卷，第 3 号， p.21-36, 2003年

非特許文献 2 ：

ビー . シー . トム（ B. C. Tom ) ' エイ . ケィ . 力トサゲロス（

A. K. Katsaggelos) 共著、「リコンストラクションオファーハイレゾルーションイメージバイシマルティニァスレジストレイシヨンレストレイシヨンアンドインターポレイシヨンオフローレゾルーションィメ ¹ ^ンス ( Reconstruct 1 on of a high-resolution image by simultaneous registration, restoration, and interpolation of low-resolution images ) J , プロク . I E E E イント . コンフ . イメージプロセシング（Proc. IEEE Int. Conf. Image Processing) ，第 2巻， p.539 - 542, 1995年非特許文献 3 ：

アール. アール. シュルツ (R. R. Schulz) · アーノレ. エル. スティブンスン (R. L. Stevenson) 共著、「ェクストラクシヨンオフハイレゾル一シヨンフレームズフロムビデオシーケンス ( Extraction of high-resolution frames from video sequences) 」， I ji E £ トフンス . イメージプロセシング (IEEE Trans. Image Processing) ，第 5 卷， p.996-1011, 1996年非特許文献 4 ：

ェイチ.スターク（H. Stark) · ピー.ォスコウイ（f. Oskoui) 共著、「ノ、ィレゾノレーシヨンイメージリカノくリーフロムイメージプレーンアレーズユーシングコンベックスプロジェクシヨンズ（ High resolution image recovery from image-plane arrays, using convex projections) 」，ジエー. ォプ卜 . ソク . ェム. エイ ' (J. Opt. Soc. Am. A) ，第 6巻， p.1715 - 1726, 1989年非特許文献 5 ：

ヒ、、一.ジトーノ (B. Zitova) ' ジエー . フノレッサノレ（J.Flusser) 共著、「イメージレジストレーシヨンメソッドズ：ァサーベイ（ Image registration methods: a survey; 」 , メ ¹ ~、ノゾン卜ヒンヨンコンピューティング (Image and Vision Computing) ，第 21 卷，第 11号， p.977 1000, 2003年非特許文献 6 ：

ディー.ケープル（D. Capel) 著、「イメージモザィキングアンドスーノ一レゾリューション (Image Mosaicing and Super-resolut i on) J ，スプリンガー (Springer) ， 2003年非特許文献 7 ：

ダプリュ . ェイチ. プレス（W. H.Pi-ess) ' エス .エイ .テゥコルスカイ（ S. A. Teukolsky) 他共著、「ニューメリカルレシピズイン C + + (Numerical Recipes in C++) 」 .，ケンブリッジュニノく一シティ一プレス (Cambridge University Press) , 2002年非特許文献 8 ：

清水 · 奥富共著、「領域ベースマッチングのための 2次元同時サブピクセル推定法」、信学論（D- 11)，第 J87- D- II 巻，第 2 号， p.554- 564, 2004 年非特許文献 9 ：

清水 · 矢野 · 奥富共著、「 2次元サブピクセル同時推定を拡張した画像変形 Nパラメータ同時推定」、情報処理学会研究報告， 2004- CVIM- 143，第 2004卷，第 26号， p.81-88, 2004年

Claims

請求の範囲

1 . 位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法であって、高解像度画像空間に前記複数の低解像度画像の位置合わせを行い、位置合わせされた前記複数の低解像度画像の全ての画素を、前記高解像度画像空間において不等間隔にサンプリングされた画素とする第 1 のステップと、

前記高解像度画像空間を所定の大きさを有する複数の小領域に分割する第 2のステップと、

第 2のステップで分割された各小領域に対して、当該小領域内の所定の代表位置での推定値を当該小領域に存在する全ての画素の推定値とする第 3 のステップと、

を有することを特徴とする超解像処理の高速化方法。

2 . 前記小領域は正方形の形状を有する請求の範囲第 1項に記載の超解像処理の高速化方法。

3 . 前記正方形の一辺の大きさは、前記高解像度画像のピクセルの大きさの整数分の 1である請求の範囲第 2項に記載の超解像処理の高速化方法。

4 . 第 3のステップにおける前記推定値は、カメラモデルから得られる点広がり関数と前記 '高解像度画像の畳み込みとして推定して得られた値である請求の範固第 1項乃至請求の範囲第 3項のいずれかに記載の超解像処理の高速化方法。

5 . 前記超解像処理に必要な評価関数は次の式で表され、

ただし、 I は前記小領域に対する評価関数であり、また、 Μは前記小領域に属する全ての画素の数であり、 f _;は前記小領域内の位置座標（3C i , y i )を有する i 番目の画素の観測値を表し、（ X _c， y c )は前記小領域の代表位置を表し、 /( ，）は前記代表位置の推定値を表す請求の範囲第 1項乃至請求の範囲第 4項のいずれかに記載の超解像処理の高速化方法。

6 . 前記超解像処理に必要な評価関数は次の式で表され、

I = M f - f(x_c,y_c) ここで、下記数式が成立し、

U

f=—Y f.

^J M ' ただし、 I は前記小領域に対する評価関数であり、また、 Mは前記小領域内に属する全ての画素の数であり、 f iは前記小領域内の位置座標（ X i， y ； )を有する i 番目の画素の観測値を表し、（ X _c， y _e )は前記小領域の代表位置を表し、）は前記代表位置の推定値を表し、ァは前記小領域内の画素の観測値の平均を表す請求の範囲第 1項乃至請求の範囲第 4項のいずれかに記載の超解像処理の高速化方法。

7 . 前記超解像処理に必要な評価関数は次の式で表され

ただし、 I は前記小領域に対する評価関数であり、また、 Mは前記小領域に属する全ての画素の数であり、 f ；は前記小領域内の位置座標（χ ; , y ； )を有する i 番目の画素の観測値を表し、（ X _c， y c )は前記小領域の代表位置を表し、 ( , は前記代表位置の推定値を表し、 w ( _{X i} — X _c , y ； - y Jは重み関数で、前記小領域内の位置座標（X ；， y i )を有する i 番目の画素に対応する重みを表し、その値は位置座標（X い y _; ) を有する i番目の画素が前記代表位置（ χ _c， y c )から離れるに従って小さくなるようにする請求の範囲第 1項乃至請求の範囲第 4項のいずれかに記載の超解像処理の高速化方法。

8 . 前記超解像処理に必要な評価関数は次の式で表され、

ここで、下記数式が成立し、

ただし、 I は前記小領域に対する評価関数であり、また、 Mは前記小領域内に属する全ての画素の数であり、 f _;は前記小領域内の位置座標（ x i， y i )を有する i 番目の画素の観測値を表し、（X _c， y 。）は前記小領域の代表位置を表し、）は前記代表位置の推定値を表し、 w ( x i 一 _X c， y i — _y Jは重み関数で、前記小領域内の位置座標（X _;， y i ) を有する i番目の画素に対応する重みを表し、その値は位置座標（X i , y を有する i番目の画素が前記代表位置（X _c， y _e )から離れるに従つて小さくなるようになっており、 Wは前記小領域内の各画素に対応する重みの合計であり、 /は前記小領域の重み付け平均を表す請求の範囲第 1項乃至請求の範囲第 4項のいずれかに記載の超解像処理の高速化方法