WO2006118352A1

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

Info

Publication number: WO2006118352A1
Application number: PCT/JP2006/309437
Authority: WO
Inventors: Masayuki Tanaka; Masatoshi Okutomi
Original assignee: Tokyo Institute Of Technology
Priority date: 2005-05-02
Filing date: 2006-05-01
Publication date: 2006-11-09
Also published as: CN101164082B; US8041151B2; CN101164082A; US20090080806A1; EP1879146A1; JP4126378B2; JP2006309649A

Abstract

再構成型超解像処理において、評価関数と評価関数の高解像度画像に対する微分の計算を高速化することによって超解像処理の高速化を実現した超解像処理の高速化方法を提供する。　位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法であって、評価関数と前記評価関数の前記高解像度画像に対する微分の計算を、基本的な画像演算の組合せにより計算するようにする。

Description

明細書超解像処理の高速化方法技術分野

本発明は、複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化する超解像処理の高速化方法に関し、特に、再構成型超解像処理において、超解像処理の計算コストを大幅に削減することにより、超解像処理の高速化を実現させる超解像処理の高速化方法に関する。背景技術

位置ずれのある複数の低解像度画像より一つの髙解像度画像を推定する超解像処理に関しては、近年数多くの研究で報告されている（非特許文献 1参照）。例えば、非特許文献 2に開示されている M L (Maximum - 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関数）に基づき、全ての低解像度画像の画素毎に、その画素値を推定し、その推定値と観測された画素値（観測値）の差が小さくなるような高解像度画像を探索するという共通の特徴を有しており、これらの超解像処理方法は、再構成型超解像処理方法と呼ばれている。

上述のように、再構成型超解像処理は、高解像度面像に対して定義される評価関数の最適化問題として定式化される。つまり、再構成型超解像処理は、推定された低解像度画像と観測された観測画像の二乗誤差に基づく評価関数の最適化問題に帰着される。

再構成型超解像処理では、未知数の非常に大きな最適化計算であるため、その評価関数の最適化には、例えば最急降下法などの繰返し計算法がよく利用される。このとき、繰返しごとに、評価関数と評価関数の高解像度画像に対する微分を計算する必要がある。しかし、この繰返し計算のコストが非常に大きい。また、繰返し計算では一度の繰返しに対して、全ての低解像度画像の画素に対して推定を行う必要があり、その計算コストも非常に大きい。

つまり、再構成型超解像処理では、計算コストが大きいので、計算コストの低減が大きな課題となっている。

この課題を解決するために、本発明の発明者らは、評価関数を計算するための画素値推定計算回数を低減させるという観点から、日本国特許出願 2 0 0 4— 3 1 6 1 5 4に開示されている『超解像処理の高速化方法』を提案している。つまり、複数の低解像度画像を位置合わせし、高解像度画像空間に離散化点を設定し、離散化点近傍に対応する画素の平均画素値を利用することを特徴としている。このとき、離散化点近傍に対応する画素の個数を重みとして考慮することにより、推定精度を落とさず高速に計算可能であることを示した。

つまり、本発明の発明者らは、日本国特許出願 2 0 0 4— 3 1 6 1 5 4において、少ない低解像度画素の画素値推定計算で、精度を落とさずに髙解像度画像を再構成できる評価関数を利用した超解像処理の高速化手法を提案した。

より詳細に説明すると、日本国特許出願 2 0 0 4— 3 1 6 1 5 4に開示されている『超解像処理の高速化方法』では、レジストレーシヨン後 (位置合わせ後）の複数の低解像度画像の画素が高解像度画像空間において不等間隔にサンプリングされた画素と考えられ、解像度画像空間を小領域に分割し、その小領域に含まれる画素の平均画素値を利用することを特徴としている。この小領域に対する評価関数は、下記数 1 のように表される。

【数 1 】

ここで、下記数 2が成立する。

【数 2】

ただし、 I は小領域に対する評価関数を、 Mは小領域に含まれる画素の個数を、 f ；は小領域内の i 番目の画素の画素値を、（ X c , y。）は小領域の代表位置を、 Λ ,Λ)は小領域の代表位置に封する画素値の推定値を、それぞれ表す。実際の計算では、全ての小領域に対する評価関数の和が超解像処理全体の評価関数となる。

そして、本発明の発明者らは、さらに、日本国特許出願 2 0 0 4— 3 1 6 1 5 4に提案されている評価関数を利用し、その評価関数と評価関数の高解像度画像に対する微分の計算を高速化させる'という観点から、本発明に係る超解像処理の高速化方法を提案する。

要するに、本発明は、上述のような事情よりなされたものであり、本発明の目的は、再構成型超解像処理において、評価関数と評価関数の高解像度画像に対する微分の計算を高速化することによって超解像処理の高速化を実現した超解像処理の高速化方法を提供することにある。発明の開示

本発明は、位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法に関し、本発明の上記目的は、評価関数と前記評価関数の前記高解像度画像に対する微分の計算を、基本的な画像演算の組合せにより計算することによって効果的に'達成される。

また、本発明の上記目的は、前記評価関数は、空間領域の前記高解像度画像に対して、次のように定式化され、

7(h)=|b * h -f + c * h 前記評価関数の空間領域の前記高解像度画像に対する微分は、次のように定式化され、 = 2b'*(w®(b*h-f)) + 2Qfc'*c*h ただし、 hは空間領域の前記高解像度画像のベクトル表現を、 bは p s F画像のベクトル表現を、 fは平均観測画像のベクトル表現を、 wは重み画像のべクトル表現を、 cは前記高解像度画像の事前情報を表すカーネル画像のべクトル表現を、 αは拘束の強さを表す拘束パラメータを、 *は畳込み積分を、 ®は要素毎の乗算を、それぞれ表すことによって一層効果的に達成される。 . で、二つのノノレム

₂は、下記のように定義され、

XI = (w®x)

2 一 T

X„ =X X ただし、は複素共役を、 τは転置を、それぞれ表す。

また、 b'，c'は、下記のように定義され、

ただし、 Fはフーリエ変換を、 F— ¹は逆フーリエ変換を、それぞれ表す。 ■ '

前記評価関数の空間領域の前記高解像度画像に対する微分は、次のように定式化され、 w ® [b (8) [h]] -f)]]+2or F~^l [c'(g)c ® [h]]

ただし、 hは空間領域の前記高解像度画像のベクトル表現を、 bは p s F画像のべクトル表現を、 bは bのフーリエ変換を、 cは前記高解像度画像の事前情報を表すカーネル画像のべクトル表現を、 £は cのフーリェ変換を、 b'は b'のフーリエ変換を、 έ'は c'のフーリエ変換を、 f は平均観測画像のベクトル表現を、 wは重み画像のベクトル表現を、 _α は拘束の強さを表す拘束パラメータを、 *は畳込み積分を、 ®は要素毎の乗算を、 Fはフーリエ変換を、 F — ¹は逆フーリエ変換を、それぞれ表すことによって一層効果的に達成される。

.また、本発明の上記目的は、前記評価関数は、周波数領域の前記高解像度画像に対して、次のように定式化され、

前記評価関数の周波数領域の前記高解像度画像に対する微分は、うに定式化され、 w ® IF- ¹ [£ ® hl-f)l+2-c'®c ® h

ただし、 hは周波数領域の前記高解像度画像のベクトル表現で、 bは p

S F画像のベクトル表現を、 bは bのフーリエ変換を、 b 'は b'のフーリエ変換を、 f は平均観測画像のべクトル表現を、 wは重み画像のベタトル表現を、 cは前記高解像度画像の事前情報を表すカーネル画像のベタトル表現を、 cは cのフーリエ変換を、 έ'は c'のフーリエ換を、 _α は拘束の強さを表す拘束パラメータを、 ρはフーリエ変換による定数倍の違いを規格化する定数を、 *は畳込み積分を、 ®は要素毎の乗算を、 Fはフーリエ変換を、 F ¹は逆フーリエ変換を、それぞれ表すことによって一層効果的に達成される。図面の簡単な説明

第 1 図は、本発明実施例 1 の超解像処理の高速化方法（方法 1 ) に使用される評価関数の誤差項及ぴその誤差項の高解像度画像に対する微分の計算手順を示すブロック図である。

第 2図は、本発明実施例 1 の超解像処理の高速化方法（方法 1 ) に使用される評価関数の拘束項及びその拘束項の高解像度画像に対する微分の計算手順を示すプロック図である。

第 3図は、本発明実施例 2の超解像処理の高速化方法（方法 2 ) に使用される評価関数の誤差項及びその誤差項の高解像度画像に対する微分の計算手順を示すプロック図である。

第 4図は、本発明実施例 2の超解像処理の高速化方法（方法 2 ) に使用される評価関数の拘束項及ぴその拘束項の高解像度画像に対する微分の'計算手順を示すプロック図である。

第 5図は、本発明実施例 3の超解像処理の高速化方法（方法 3 ) に使用される評価関数の誤差項及びその誤差項の周波数領域での高解像度画像に対する微分の計算手順を示すプロック図である。

第 6図は、本発明実施例 3の超解像処理の高速化方法（方法 3 ) に使用される評価関数の拘束項及びその拘束項の周波数領域での高解像度画像に対する微分の計算手順を示すプロック図である。

第 7図は、本発明に係る超解像処理の高速化方法を適用した超解像処理の結果を示す図である。第 7図（A ) は、数 3で表す『平均画素値による評価関数』に基づく超解像処理の結果を、第 7図（B ) は、本発明実施例 1 の超解像処理の高速化方法を適用した超解像処理の結果を、第 7図（C ) は、本発明実施例 2の超解像処理の高速化方法を適用した超解像処理の結果を、第 7図（D ) は、本発明実施例 3 の超解像処理の高速化方法を適用した超解像処理の結果を、それぞれ示している。発明を実施するための最良の形態

以下、本発明を実施するための最良の形態を図面を参照して説明する本発明では、日本国特許出願 2 0 0 4— 3 1 6 1 5 4に開示されている評価関数を利用し、その評価関数と評価関数の高解像度画像に対する微分の計算を基本的な画像演算の組合せにより計算することにより、また、更にフーリエ変換を応用することにより、超解像処理の高速化を実現するようにしている。

ところで、超解像処理において、位置ずれのある複数枚の観測画像（低解像度画像）の各画素は、レジストレーシヨン（位置合わせ）により、高解像度画像空間のある位置に対応付けられる。つまり、レジストレーシヨンの後では、複数枚の観測画像は、高解像度空間内で不等間隔にサンプリングされた画素であると考えることができる。

この不等間隔にサンプリングされた画素位置（以下、単に観測画素位置とも称する）を高解像度画像の画素位置（以下、単に、高解像度画素位置とも称する）で近似することを考える。このとき、ある高解像度画素位置に近似される観測画素（観測画素位置）は、複数ある場合が考えられる。逆に、近似される観測画素（観測画素位置）が存在しない高解像度画素位置も存在する。ここで、各高解像度画素位置に近似された複数の観測画素の平均画素値をそれぞれ計算することにより、 1つの画像を生成することができる。この画像を平均観測画像と呼ぶことにする。平均観測画像は、その画素間隔（画素数）が高解像度画像と等しく、高解像度画像をカメラモデルから得られる点広がり関数（P S F関数）でぼかした画像に相当する。ただし、近似される観測画素が存在しない画素位置に対しては、その画素値は定義されない。各高解像度画素位置に近似された観測画素の個数も、同様に 1つの画像となる。その画像を重み画像と呼ぶことにするまた、 P S Fも 1つの画像と考えることにする。そうすると、評価関数と評価関数の高解像度画像に対する微分を計算するに際して、考慮しなければならない画像は、高解像度画像 h [ i ， j ]、 P S F画像 b [ i , j ]、平均観測画像 f [ i ， j ]、重み画像 w [ i ，； j ]の四つである。本発明では、この 4種類の画像間の画像演算の組合せにより、評価関数と評価関数の高解像度画像に対する分の計算を行うようにしている。

具体的には、本発明では、日本国特許出願 2 0 0 4 — 3 1 6 1 5 4に開示されている評価関数において、小領域の大きさを高解像度画像の画素と等しく設定し、小領域の代表位置を画素中心とすることを前提とする。

このとき、超解像処理全体の評価関数は、下記数 3で表すように、誤差項（数 3の右辺の第 1項）と拘束項（数 3の右辺の第 2項）とからなる。なお、下記数 3で表す評価関数を『平均画素値による評価関数』とも称する。

【数 3】 0

ここで、 hは高解像度画像のべクトル表現を、 Nは小領域の数を、は i 番目の小領域の代表位置に対応する P S Fのべクトル表現を、 w iは i 番目の小領域に含まれる画素の数を、は i 番目の小領域に含まれる画素の平均画素値を、 Cは高解像度画像の事前情報を表すカーネルを、 αは拘束の強さを表す拘束パラメータを、 *は畳込み演算を、それぞれ表す。

なお、上述したように、小領域として高解像度画像の面素の大きさを設定しているため、小領域の数 Νは、高解像度画像の画素の数と等しく、小領域の代表位置は高解像度画像の画素中心と等しい。また、拘束パラメータ αが 0の場合（ α = 0 ) も考えられ、このとき、超解像処理全体の評価関数は、謓差項（数 3の右辺の第 1項）のみからなるものとする。

本発明に係る超解像処理の高速化方法は、数 3に示すような、誤差項と拘束項とからなる評価関数に基づく再構成型超解像処理に適用される本発明では、上述した 4種類の画像、即ち、高解像度画像、 P S F画像、平均観測画像、重み画像を用いて、数 3で表す『平均画素値による評価関数』を再定義して得られたものを超解像処理全体の評価関数とし、再定義された評価関数及びその微分の計算方法（詳細は実施例で説明する）によって、評価関数及びその微分の計算を行えば、超解像処理の高速化を実現することができる。また、本発明では、 4種類の画像で再定義された評価関数は、共通して誤差項と拘束項とから成る。 1

【実施例 1 】

実施例 1 の超解像処理の高速化方法（以下、単に、方法 1 とも称する ) では、上述した 4種類の画像を利用し、高解像度画像の空間表現（以下、単に、空間領域の高解像度画像 h、又は、高解像度画像 hとも称する）に対して、数 3で表す『平均画素値による評価関数』を再定義する。また、再定義された評価関数と評価関数の高解像度画像に対する微分の計算も、全て空間領域においてなされる。

まず、実施例 1 の超解像処理の高速化方法では、 '上述した 4種類の画像を用いて、超解像処理全体の評価関数が下記数 4で、また、評価関数の高解像度画像に対する微分が下記数 5 で、それぞれ定式化される。

なお、数 4の右辺の第 1項を方法 1 の誤差項と、数 4の右辺の第 2項を方法 1 の拘束項と、数 5の右辺の第 1項を方法 1 の誤差項の高解像度画像に対する微分と、数 5の右辺の第 2項を方法 1 の拘束項の高解像度画像に対する微分と、それぞれ称する。

【数 4】

1(h) =||b * h - f||² +a ||c * h|

【数 5】

7

— =2b ^,*(w ® (b *h -f))+2o c^,* c * h ここで、 hは高解像度画像のベクトル表現を、 bは P S F画像のベタトル表現を、 f は平均観測画像のべクトル表現を、 wは重み画像のベタトル表現を、 cは高解像度画像の事前情報を表すカーネル画像のべクトル表現を、 αは拘束の強さを表す拘束パラメータを、それぞれ表す。また、 *は畳込み積分を表し、 ®は要素毎の乗算を表す。 2 さらに、二つのノルム

は、下記数 6、数 7のように定義される

【数 6】

【数 7】二 X X ここで、は複素共役を、 Tは転置を、それぞれ表す。

また、 b',C'は、下記数 8、数 9のように定義される。

【数 8】

【数 9】

c' = [ ]

ここで、 Fはフーリエ変換を、 F ¹は逆フーリエ変換を、それぞれ表す。

次に、方法 1の誤差項及ぴその誤差項の高解像度画像に対する微分の計算手順をブロック図の第 1 図に示す。また、方法 1 の拘束項及ぴその拘束項の高解像度画像に対する微分の計算手順をプロック図の第 2図に示す。

なお、第 1 図及び第 2図において、 *は畳込み積分演算で、 +は要素毎の加算で、 ®は要素毎の乗算で、 ∑は要素の総和を計算する演算である。

第 1図のプロック図に示さ.れるように、方法 1 の誤差項及びその誤差項の高解像度画像に対する微分の計算手順は次のようになっている。ステップ 1 A ：

高解像度画像 hと P S F画像 bとの畳込み積分演算を行う。

ステップ 1 B ：

ステップ 1 Aでの畳込み積分演算結果から平均観測画像 f を引くように演算する。

ステップ 1 C ：

ステップ 1 Bでの'演算結果と重み画像 Wとの乗算を行う。

ステップ 1 D ：

ステップ 1 Cでの乗算結果と P S F画像 b 'との畳'込み積分演算を行うことにより、方法 1 の評価関数の誤差項の微分が得られる。

ステップ 1 E ：

ステップ 1 Bでの演算結果とステップ 1 Cでの乗算結果との乗算を行

5 。

ステップ 1 F ：

ステップ.1 Eでの乗算結果に基づいて要素の総和を計算する演算により、方法 1 の評価関数の誤差項の値が得られる。

また、第 2図のブロック図に示されるように、方法 1 の拘束項及ぴその拘束項の高解像度画像に対する微分の計算手順は次のようになつている。

ステップ 1 a ：

高解像度画像 hと高解像度画像の事前情報を表すカーネル画像 cとの畳込み積分演算を行う。

ステップ 1 b ：

ステップ 1 a での畳込み積分演算結果とカーネル画像 C 'との畳込み 4

積分演算を行うことにより、方法 1 の評価関数の拘束項の微分が得られる。

ステップ 1 c ：

ステップ 1 a での畳込み積分演算結果とステップ 1 a での畳込み積分演算結果との乗算を行う。

ステップ 1 d ：

ステップ 1 c での乗算結果に基づいて要素の総和を計算する演算により、方法 1 の評価関数の拘束項の値が得られる。

【実施例 2】

実施例 2の超解像処理の高速化方法（以下、単に、方法 2 とも称する ) では、上述した 4種類の画像を利用し、高解像度画像の空間表現に対して、数 3で表す『平均画素値による評価関数』を再定義する。方法 1 との違いとして、方法 2は、畳込み演算の計算をフーリエ変換を利用して計算する点である。

実施例 2 'の超解像処理の高速化方法では、上述した 4種類の画像を用いて、超解像処理全体の評価関数が下記数 1 0で、また、評価関数の高解像度画像に対する微分が下記数 1 1で、それぞれ定式化される。

なお、数 1 0の右辺の第 1項を方法 2の誤差項と、数 1 0の右辺の第 2項を方法 2の拘束項と、数 1 1 の右辺の第 1項を方法 2の誤差項の高解像度画像に対する微分と、数 1 1 の右辺の第 2項を方法 2の拘束項の高解像度画像に対する微分と、それぞれ称する。

【数 1 0 】

1(h) = F~^l\b 0 [hll— f I + ^_1 [c ® hl] 【数 1 b ® [h]]-f)]]+2o -¹ c'(8)c ® F[h]]

ここで、 bは bのフーリエ変換を、 cは cのフーリエ変換を、 b'は b' のフーリエ変換を、は c'のフーリエ変換を、それぞれ表す。また、数 1 0及び数 1 1 において、他のそれぞれの記号は、方法 1 と同じであるため、その説明を省略する。 ■

次に、方法 2の誤差項及びその誤差項の高解像度画像に対する微分の計算手順をブロック図の第 3図に示す。また、方法 2の拘束項及びその拘束項の高解像度画像に対する微分の計算手順をプロック図の第 4図に示す

なお、第 3 図及ぴ第 4図において、 Fはフーリエ変換を、 F ¹は逆フーリエ変換をそれぞれ表す。また、 +は要素毎の加算で、 ®は要素毎の乗算で、 ∑は要素の総和を計算する演算である。

第 3図のプロック図に示されるように、方法 2の誤差項及ぴその誤差項の高解像度画像に対する微分の計算手順は次のようになっている。ステップ 2 A ：

髙解像度画像 hのフーリエ変換を行う。

ステップ 2 B ：

ステップ 2 Aでのフーリェ変換結果と P S F画像 bとの乗算を行う。ステップ 2 C ：

ステップ 2 Bでの乗算結果の逆フーリエ変換を行う。

ステップ 2 D ：

ステップ 2 Cでの逆フーリェ変換結果から平'均観測画像 f を.引くように演算する。ステップ 2 E ：

ステップ 2 Dでの演算結果と重み画像 Wとの乗算を行う。

ステップ 2 F ：

ステップ 2 Eでの乗算結果のフーリエ変換を行う。

ステップ 2 G ：

ステップ 2 Fでのフーリエ変換結果と P S F画像 b 'との乗算を行うステップ 2 H ：

ステップ 2 Gでの乗算結果の逆フーリエ変換を行うことにより、方法 2の評価関数の誤差項の微分が得られる。

ステップ 2 I ：

ステップ 2 Dでの演算結果とステップ 2 Eでの乗算結果との乗算を行 Ό 。

ステップ 2 J ：

ステップ 2 I での乗算結果に基づいて要素の総和を計算する演算により、方法 2の評価関数の誤差項の値が得られる。

また、第 4図のプロック図に示されるように、方法 2の拘束項及びその拘束項の高解像度画像に対する微分の計算手順は次のようになつている。

ステップ 2 a ：

高解像度画像 hのフーリエ変換を行う。

ステップ 2 b ：

ステップ 2 a でのフーリエ変換結果とカーネル画像 £との乗算を行うステップ 2 c ステップ 2 bでの乗算結果の逆フーリエ変換を行う。

ステップ 2 d ：

ステップ 2 c での逆フーリエ変換結果に対して、フーリエ変換を行う

ステップ 2 e ：

ステップ 2 d でのフーリェ変換結果とカーネル画像との乗算を行

0。 ―

ステップ 2 f ：

ステップ 2 eでの乗算結果の逆フーリエ変換を行うことにより、方法 2の評価関数の拘束項の微分が得られる。

ステップ 2 g ：

ステ.ップ 2 c での逆フーリエ変換結果とステップ 2 c での逆フーリエ変換結果との乗算を行う。

ステップ 2 h ：

ステップ 2 gでの乗算結果に基づいて要素の総和を計算する演算により、方法 2の評価関数の拘束項の値が得られる。

【実施例 3】

実施例 1及び実施例 2では、空間領域の高解像度画像 hに対して、数 3で表す『平均画素値による評価関数』が再定義され、空間領域の高解像度画像に対する最適化を行う方法を述べた。

ところで、よく知られているように、空間領域の高解像度画像 hとそのフーリエ変換に対応する周波数領域の高解像度画像 h (以下、単に、周波数領域の高解像庳画像とも称する）には、一対一の関係がある。よって、周波数領域での高解像度画像を最適化することと、空間領域での高解像度画像を最適化することは、等価な処理であると言える。

そのため、空間領域の高解像度画像 hに対してではなく、周波数領域の髙解像度画像 hに対して最適化計算を行い、高解像度画像を再構成することも可能である。このように、周波数領域での高解像度画像 hに対して最適化を行うことを考えると、その評価関数 0) は下記数 1 2のように表される。

つまり、実施例 3の超解像処理の高速化方法（以下、単に、方法 3 とも称する）では、数 3で表す『平均画素値による評価関数』が、高解像度画像の周波数表現、即ち、周波数領域での高解像度画像 hに対して再定義される。方法 1や方法 2 と比べて、方法 3は評価関数が高解像度画像の周波数表現に対して再定義されていることが大きな特徴である。

. 実施例 3の超解像処理の高速化方法では、超解像処理全体の評価関数が下記数 1 2で、また、評価関数の周波数領域での高解像度画像 hに対

31

する微分^が下記数 1 3で、それぞれ定式化される。

oh

なお、数 1 2の右辺の第 1項を方法 3の誤差項と、数 1 2の右辺の第 2項を方法 3の拘束項と、数 1 3の右辺の第 1項を方法 3の誤差項の周波数領域での高解像度画像 hに対する微分と、数 1 3の右辺の第 2項を方法' 3の拘束項の周波数領域での高解像度画像 hに対する微分と、それぞれ称する。

【数 1 2】

【数 1 3】 c^,(2)c ( ) h

ここで、 hは周波数領域での高解像度画像のベクトル表現で、 _P はフ —リエ変換による'定数倍の違いを規格化する定数である。また、数 1 2 及ぴ数 1 3において、他のそれぞれの記号は、方法 1 、方法 2 と同じであるため、その説明を省略する。

次に、方法 3の誤差項及ぴその誤差項の周波数領域での高解像度画像に対する微分の計算手順をブロック図.の第 5図に示す。また、方法 3の拘束項及びその拘束項の周波数領域での高解像度画像に対する微分の計算手順をプロック図の第 6図に示す。

なお、第 5図及び第 6図において、 Fはフーリエ変換を、 F ¹は逆フーリエ変換をそれぞれ表す。また、 +は要素毎の加算で、 ®は要素毎の乗算で、 ∑は要素の総和を計算する演算で、 Con j .は各要素の複素共役をとる演算である。

第 5図のブロック図に示されるように、方法 3.の誤差項及ぴその誤差項の周波数領域での高解像度画像 hに対する微分の計算手順は次のようになっている。

ステップ 3 A ：周波数領域での高解像度画像 hと P S F画像 bとの乗算を行う。

ステップ 3 B ：

ステップ 3 Aでの乗算結果の逆フーリエ変換を行う。

ステップ 3 C ：

ステップ 3 Bでの逆フーリエ変換結果から平均観測画像 f を引くように演算する。

ステップ 3 D ：ステップ 3 Cでの演算結果と重み画像 Wとの乗算を行う。

ステップ 3 E ：

ステップ 3 Dでの乗算結果のフーリエ変換を行う。

ステップ 3 F ：

ステップ 3 Eでのフーリェ変換結果と P S F画像 D 'との乗算を行うことにより、方法 3の評価関数の誤差項の周波数領域での高解像度画像 hに対する微分が得られる。

ステップ 3 G ：

ステップ 3 Cでの演算結果とステップ 3 Dでの乗算結果との乗算を行う。

ステップ 3 H ：

ステップ 3 Gでの乗算結果に基づいて要素の総和を計算する演算により、方法 3の評価関数の誤差項の値が得られる。

また、第 6図のブロック図に示されるように、方法 3の拘束項及びその拘束項の周波数領域での高解像度画像 ϋに対する微分の計算手順は次のようになっている。

ステップ 3 a ：

周波数領域での高解像度画像 hとカーネル画像 έとの乗算を行う。ステップ 3 b ：

ステップ 3 a での乗算結果とカーネル画像との乗算を行うことにより、方法 3の評価関数の拘束項の周波数領域での高解像度画像 hに対する微分が得られる。

ステップ 3 c ：

ステップ 3 a での乗算結果に対して、各要素の複素共役をと.る演算を行う。 2

ステップ 3 d ：

ステップ 3 c での演算結果とステップ 3 a での乗算結果との乗算を行

5。 .

ステップ 3 e ：

ステップ 3 dでの乗算結果に基づいて要素の総和を計算する演算により、方法 3の評価関数の拘束項の値が得られる。

以上のように、本発明の具体的な実施例について説明したが、次に、手持ちカメラにより撮影された実画像に対して、本発明に係る超解像処理の高速化方法を適用し、本発明の有効性を確認する。超解像処理には、 1 6枚（ペイヤー配列）の観測画像を低解像度画像として利用した。低解像度画像の 6 0 X 6 0の領域を、倍率 4 X 4に超解像処理し、 2 4 0 X 2 4 0の高解像度画像を再構成した。 P S Fはガウシアン P S Fを利用した。拘束カーネルとして、ラプラシアン拘束カーネルを利用した。最適化計算は蕞急降下法を利用し、繰返し計算回数は 2 0回とした。レジストレーシヨンは、手持ちカメラである.ことを考慮し、非特許文献 5 に開示されている「射影変換を仮定した勾配法」を利用した。なお、計算には、 2. 8 [G H z ] (Pentium4 ) (登録商標）の C P Uを使用した。

第 7図には、本発明の発明者らが提案した日本国特許出願 2 0 0 4 — 3 1 6 1 5 4に開示されている『超解像処理の高速化方法』を適用した超解像処理の結果と、本発明に係る超解像処理の高速化方法を適用した超解像処理の結果を示している。

具体的に、第 7図（A) は、数 3で表す『平均画素値による評価関数』に基づく超解像処理の結果を、第 7図（B ) は、本発明実施.例 1 の超解像処理の高速化方法を適用した超解像処理の結果を、第 7図（ C ) は、本発明実施例 2の超解像処理の高速化方法を適用した超解像処理の結果を、第 7図（D) は、本発明実施例 3 の超解像処理の高速化方法を適用した超解像処理の結果を、それぞれ示している ο

第 7図（A) 、第 7図（B ) 、第 7図 ( C ) 、第 7図（ D ) に示す 4 種類の超解像処理結果は、主観評価としてほぼ |pj一である。次に、超解像処理の計算時間を見てみる。第 7図（ A ) の解像処理結果を得るための計算時間は 9. 6 2 [sec]で、第 7図 ( B ) の超解像処理結果を得るための計算時間は 1 . 3 8 [sec]で、第 7図 ( C ) の超解像処理結果を得るための計算時間は 1 . 0 1 [sec]で、図 7 (D ) の超解像処理結果を得るための計算時間は 0. 5 9 [sec]である ο

ところで、日本国特許出願 2 0 0 4— 3 1 6 1 5 4に開示されているように、日本国特許出願 2 0 0 4— 3 1 6 1 5 4に係る超解像処理の高速化方法を適用した超解像処理、つまりヽ数 3で表す『平均画素値による評価関数』に基づく超解像処理は、従来の超解像処理より、高速になつていることが明らかである。上記の超解像処理の計算時間を比較すると、日本国特許出願 2 0 0 4— - 3 1 6 1 5 4に開示されている方法による超解像処理より、本発明に係る超解像処理の高速化方法による超解像処理は、更なる高速化を実現したことが良く分かる 0

つまり、『平均画素値による評価関数』に基づ <超解像処理と比較して、本発明の方法 1 は計算時間で 6. 9 7倍、方法 2は計算時間で 9. 5 2倍、方法 3は計算時間で 1 6. 3倍、計算が高速になっていることが確認された。産業上の利用可能性

以上のように、本発明に係る超解像処理の高速化方法を用いれば、評価関数と評価関数の高解像画像に対する微分の計算を、基本的な画像の演算の組合せにより計算することにより、超解像処理に必要な計算コストを大幅に削滅し、よって超解像処理の高速化を実現することができたという優れた効果を奏する。 ·

また、実画像に対して、本発明に係る超解像処理の高速化方法を適用することによって、超解像処理の高速化が実現されたことを確認した。

<参考文献一覧 >

非特許文献 1 ：

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

) 共著、「スーパーレゾノレーシヨンイメージリコンストラクション：ァーテク二力 /レォー /ヽービュー ( Super-Resolution Image Reconstruction: A Technical Overview) J , I E E E シグナルフロク . マガジン (IEEE Signal Proc. Magazine) ，第 26 卷，第 3 号， p.21-36, 2003年非特許文献 2 ：

ビー . シー . トム（ B. C. Tom ) ' エイ . ケィ . 力トサゲロス（ A. K. Katsaggelos) 共著、「リコンストラクションオファーハイレゾルーシヨンイメージバイシマルティニァスレジストレイシヨンレストレイシヨンアンドインターポレイシヨンオフローレゾノレーシヨンィメ一シス (Reconstruction 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 E E E トランス. イメージプロセシング (IEEE Trans. Image Processing) ，第 5 巻， p.996 - 1011, 1996年非特許文献 4 ：

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

エス . ベイカ（S. Baker) ' アイ .マシューズ ( I. Matthews) 共著、「ノレ一カスーカナデ 2 0ィヤーズオン：ァユニファイングフレームワーク 1, Lucas-Kanade 20 Years On: A Unifying Framework) 」，ンターナショナルジャーナルオフコンピュータビジョン（ International Journal of Computer Vision) ，第 56巻，第 3 ， p.221 - 255, 2004年

Claims

請求の範囲

1 . 位置ずれを含む複数の低解像度画像から一つの高解像度画像を推定する超解像処理を高速化するための超解像処理の高速化方法であって、評価関数と前記評価関数の前記高解像度画像に対する微分の計算を、基本的な画像演算の組合せにより計算することを特徴とする超解像処理の高速化方法。

2 . 前記評価関数は、空間領域の前記高解像度画像に対して、次のように定式化され、

前記評価関数の空間領域の前記高解像度画像に対する微分は、次のように定式化され、

dl

= 2b'*(w ®ib *h— f)) + 2 c'*c*h ただし、 hは空間領域の前記高解像度画像のベクトル表現を、 Dは p s F画像の.ベクトル表現を、 f は平均観測画像のベクトル表現を、 wは重み画像のべクトル表現を、 cは前記高解像度画像の事前情報を表す力一ネル画像のべクトル表現を、 αは拘束の強さを表す拘束パラメータを、 *は畳込み積分を、 ®は要素毎の乗算を、それぞれ表し、つのノノレム

は、下記のように定義され、 X ( w

®x) I II² —^T

XL = X ただし、 —は複素共役を、 Tは転置を、それぞれ表し、

また、 b'，c'は、下記のように定義され、

ただし、 Fはフーリエ変換を、 F— は逆フーリエ変換を、それぞれ表す請求の範囲第 1項に記載の超解像処理の高速化方法。

3 . 前記評価関数は、空間領域の前記高解像度画像に対して、次のように定式化され、 =| 一¹ L^b® ] f +a ^_1[c(8) [h]1 前記評価関数の空間領域の前記高解像度画像に対する微分は、次のように定式化され、 w<S)[F — 1 b ® F[h]] - f )]]+ 2 a F一¹ [c'(8)c ®F\ ^ ただし、 hは空間領域の前記高解像度画像のベクトル表現を、 bは p S F画像のベクトル表現を、 bは bのフーリエ変換を、 cは前記高解像度画像の事前情報を表すカーネル画象のべクトル表現を、 έは cのフーリエ変換を、 b'は b'のフーリエ変換を、 ^は c'のフーリエ変換を、 f は平均観測画像のベタトル表現を、 wは重み画像のべクトル表現を、ひは拘束の強さを表す拘束パラメータを、 *は畳込み積分を、 ®は要素毎の乗算を、 Fはフーリエ変換を、 F'— ¹は逆フーリエ変換を.、それぞれ表し、 2

二つのノノレム X X

w' は、下記のように定義され

2 =^{X X} ただし、は複素共役を、 Tは転置を、それぞれ表し、

また、 b'，c'は、下記のように定義され、

b'= - ¹ [兩]

c' - j

ただし、 Fはフーリエ変換を、 F ¹は逆フーリエ変換を、それぞれ表す請求の範囲第 1項に記載の超解像処理の高速化方法。

4. 前記評価関数は、周波数領域の前記高解像度画像に対して、次のように'定式化され、 c®h

前記評価関数の周波数領域の前記高解像度画像に対する微分は、次のように定式化され、

:2b*® w ® I —¹ [b ® h]-f)]+2-c'(8)c ® h

ah

ただし、 hは周波数領域の前記高解像度画像のベクトル表現で、 bは p s F画像のベクトル表現を、 bは bのフーリエ変換を、 b'は b'のフ一リエ変換を、 fは平均観測画像のベクトル表現を、 wは重み面像のベクトル表現を、 cは前記高解像度画像の事前情報を表すカーネル画像のベクトル表現を、 £は 0のフーリエ変換を、 £'は c'のフーリエ変換を、 _αは拘束の強さを表す拘束パラメータを、 Ρ はフーリエ変換による定数倍の違いを規格化する定数を、 *は畳込み積分を、 ®は要素毎の乗算を

、 Fはフーリエ変換を、 F — ¹は逆フーリエ変換を、それぞれ表し、

2

二つのノノレム

は、下記のように定義され、

X

X!

2 ^{= Χ Χ}

ただし、は複素共役を、 Τは転置を、それぞれ表し、

また、 b',c'は、下記のように定義され、

ただし、 Fはフーリエ変換を、 F — ¹は逆フーリエ変換を、それぞれ表す請求の範囲第 1項に記載の超解像処理の高速化方法。