第 5 章超解像技術 超解像 (Super Resolution) 技術とは 図 5.1(a) に示すような画像を構成するドット ( 画素 ) の数が少ない低画質な画像 ( 低解像度画像 ) から 図 5.2(b) に示すような画像を構成するドットの数が多い高画質な画像 ( 高解像度画像 ) を生成

第 5 章 超解像技術

超解像（Super Resolution）技術とは、図 5.1（a）に示すような画像を構成するドット（画素）の数が 少ない低画質な画像（低解像度画像）から、図 5.2（b）に示すような画像を構成するドットの数が多 い高画質な画像（高解像度画像）を生成する技術のことである 1,2）_{。ここで、注意してほしいのは、} 超解像技術で生成された高解像度画像は元の低解像度画像の拡大画像であるが、一般の補間 技術で生成した拡大画像（図 5.1（c））は高解像度画像ではない。超解像度技術には、複数枚の 低解像度画像から 1 枚の高解像度画像を生成するマルチフレーム超解像技術 14）_{と、類似画像} の学習による 1 枚の低解像度画像から１枚の高解像度画像を生成するシングルフレーム超解像 技術15, 16）_{などがある。本章ではそれらの技術について述べる。} 図5.1 （a）低解像度画像、（b）高解像度画像、（c）拡大画像

5.1 マルチフレーム超解像技術

画像を拡大する場合、画素と画素間に新たな画素を埋め込む必要がある。従来の補間技術では、 周辺の既存の画素値を用いて線形、または非線形関数で新しい画素値を推定している。スムーズ な拡大画像が得られるが、高解像度画像ではない。一方、マルチフレーム超解像技術の場合、周 辺の既存の画素値を用いて画素間の隙間を埋め込むのではなく、フレーム間の微小な位置ズレを subpixelの精度で推定（レジストレーション）し、他のフレームの画素で画素間の隙間を埋め込む。 そのため、高解像度画像が得られる。その様子を図5.2に示す。 マルチフレーム超解像技術を用いた例について説明する。実験に、図 5.3 に示す 16 枚の低解 像度画像を入力画像として用いた。第 4 章で述べた Optical Flow 法で推定した各フレームの位 置ズレを、それぞれの図の下に示す。本実験では、幾何学的変形は平行移動のみである。また、 第 1 番目の画像を基準画像とした。

図5.2 マルチフレーム超解像技術の基本概念図 0 0 0.023486 0.20736 -0.0077691 0.51619 -0.024578 0.7646 0.21182 -0.002746 0.22784 0.22399 0.24828 0.511 0.23901 0.77166 0.479 -0.03149 0.5056 0.23181 0.50438 0.51625 0.49163 0.77207 0.73852 -0.075242 0.72692 0.26881 0.73636 0.51597 0.73397 0.76685 図 5.3 マルチフレーム超解像実験に用いた低解像度画像（16 枚） 低解像度画像を横と縦 2 倍ずつ拡大した画像（0 を挿入したアップコンバータ）を、図 5.4（中 段）に示す。4 画素のうち、1 画素のみが濃度値が入っているが、他の 3 画素は濃度値が入って いない、 隙間だらけの画像となっているため、その隙間を埋め込むために、従来補間技術が使

われている。バイキュービック補間法で補間した結果を、図 5.4（下段）に示す。補間だけでは、た だ画像が拡大しただけで、中の文字が読めないことがわかる。即ち、解像度は向上していない。

図 5.4 補間による拡大画像

第 3 章 画像復元技術

3.1 劣化画像と点広がり関数

画像復元（restoration）とは、ボケやブレのある画像を元のきれいな画像に復元することである。 即ち、画像が劣化した過程を逆に辿る処理を行う。第 1 章の式（1.1）において、焦点ズレや手ブ レなどの画像劣化のみを考慮する場合、劣化画像 g は式（3.1）のように表すことができる。 ) , ( * ) , ( ) , ( ) , ( ) , ( y x h y x f d d h y x f y x g = − − =

_{∫ ∫}

₋∞_{∞ −}∞_∞

ξ

η

ξ

η

ξ

η

（3.1） ここで、f は劣化前の画像であり、h は劣化を表す劣化関数である。*は畳み込み演算を表す。ここ で、注意してほしいのは、式（3.1）の線形畳み込み演算で表せるのは、劣化関数 h が線形で位置 不変（space-invariant）という特性を持つ場合に限る（劣化関数 h が space-variant の場合、非線形 問題となるので、本書では取り扱わない）。劣化前の画像 f（x, y）が式（3.2）に示す点光源（2 次元 デルタ関数）の場合、劣化画像は式（3.3）のようになり、h（x, y）となる。したがって、h（x, y）は点広 がり関数（PSF：Point Spread function）とも呼ぶ。

∫ ∫

−∞∞ −∞∞ = ⎩ ⎨ ⎧∞ = = = 1 ) , ( 0 ) 0 , 0 ( ) , ( ) , ( ) , ( dxdy y x otherwise y x y x y x f

δ

δ

（3.2） ) , ( ) , ( * ) , ( ) , (x y x y h x y h x y g =

δ

= （3.3） 焦点ボケによる点広がり関数 h（x, y）は、一般に式（3.4）のようなガウシアン関数で表すことができ る22）_。 ) 2 exp( 2 1 ) , ( ₂ 2 2 2

_σ

πσ

y x y x h = − + （3.4） 図3.1 焦点ボケによる劣化画像例（（a）劣化前の原画像、（b）劣化画像（σ=2）、（b）劣化画像 （σ=4））

ここで、σ は焦点ボケの程度または点の広がり程度を表すパラメータである。σ は大きいほどより焦

点がボケている。図3.1 に焦点ボケによる画像の劣化例を示す。

MATLAB で以下のように焦点ボケの点広がり関数（PSF）を生成することができる。なお、関数 fspecial の使い方は help fspecial で調べてほしい。

PSF=fspecial（'gaussian', size, sigma）;

図3.1（a）は劣化前の原画像であり、図 3.1（b）と（c）はそれぞれ σ が 2 と 4 のときの焦点ボケ画 像である。 カメラのブレによる点広がり関数は式（3.5）のように近似することができる22）_。 ⎪ ⎩ ⎪ ⎨ ⎧ > + ≤ + = 2 sin cos , 0 2 sin cos , 1 ) , ( w y x w y x w y x h

θ

θ

θ

θ

θ （3.5） ここで、θはブレの方向（反時計回り）を表すパラメータで、w はブレによるカメラの移動距離を表 すパラメータである。MATLAB で、以下のようにカメラブレによる点広がり関数（PSF）を生成する ことができる。 PSF= fspecial（'motion', w,THETA）; 図3.2 にカメラのブレによる画像の劣化例を示す。図 3.2（a）は劣化前の原画像であり、図 3.2 （b）はブレによる点広がり関数であり、図 3.2（c）はブレによる劣化画像を示す。 図3.2 カメラブレによる劣化画像例、（a）劣化前の原画像、（b）カメラブレによる点広がり関数 （w=31, THITA=11）、（c）カメラブレによる劣化画像

3.2 逆フィルタと Wiener フィルタ

画像復元は、劣化画像 g（x, y）から原画像 f（x, y）を推定する逆問題である。式（3.1）をフーリエ 変換すると、畳み込み演算は式（3.6）のような掛け算になる。 ) , ( ) , ( ) , (u v F u v H u v G = ⋅ （3.6）

v はそれぞれ x 方向、y 方向の空間周波数を表す。

点広がり関数 h が既知の場合、G（u, v）に 1/H（u, v）をかけることにより、F（u, v）を求めることが できる。 ) , ( 1 ) , ( ) , ( ˆ v u H v u G v u F = ⋅ （3.7） また、Fˆ(u,v)を逆フーリエ変換することにより、原画像 fˆ(x,y)を推定することができる。1/H（u, v）は周波数領域における逆フィルタと呼ばれる。 一方、劣化画像にノイズが含まれている場合、 ) , ( ) , ( ) , ( ) , (x y f x y h d d n x y g =

∫ ∫

₋∞_{∞ −}∞_∞ −

ξ

−

η

ξ

η

ξ

η

+ （3.8） となり、そのフーリエ変換は式（3.9）のようになる。 ) , ( ) , ( ) , ( ) , (u v F u v H u v N u v G = ⋅ + （3.9） ここで、n（x, y）と N（u, v）はそれぞれノイズとそのフーリエ変換である。 式（3.8）と（3.9）のノイズ画像に逆フィルタをかけると、Fˆ(u,v)は、 ) , ( ) , ( ) , ( ) , ( ) , ( ) , ( ) , ( ) , ( ˆ v u H v u N v u F v u H v u N v u H v u G v u F = + = + （3.10） となる。H（u, v）が零に近い周波数領域において、わずかなノイズでも式（3.10）の第 2 項が増幅さ れてしまい、非常にノイズの多い復元画像になってしまう。 ノイズ問題を最小二乗範囲で解決しようとするのが、Wienerフィルタである。Wienerフィルタは式 （3.11）に示す原画像と復元画像との平均二乗誤差が、最も小さくなるような復元フィルタである。原 画像や付随する雑音成分が、弱定常場に属すると仮定できるときには、Wienerフィルタの空間周波 数特性は次式で与えられる。 ) , ( / ) , ( ) , ( ) , ( ) , ( ₂ * v u S v u S v u H v u H v u WF f n + = （3.11） ここで、H*(u,v)はH( vu, )の複素共役を表し、S_f ( vu, )とS_n( vu, )はそれぞれ原画像とノ イズのパワースペクトルである。式（3.11）の分母の第 2 項は劣化画像のノイズ対信号比であり、 一般に既知ではないので、適当な定数を代入して復元を行う。その項が0 の場合、単純な逆フィ ルタとなる。Wiener フィルタによる画像復元の流れ図を図 3.3 に示す。 逆フィルタと Wiener フィルタによる画像復元例を図 3.4 に示す。図 3.4（a）は劣化前の原画像 （未知）であり、図3.4（b）はノイズを含む劣化画像（焦点ボケ）である。逆フィルタと Wiener フィルタ で復元した画像は、それぞれ図 3.4（c）と 3.4（d）に示す。図 3.4（c）に示すように、逆フィルタの場 合、ノイズが増幅されてしまい、復元画像はノイズ画像になっている。一方、Wiener フィルタの場 合、比較的ノイズの影響をほとんど受けず、復元されていることがわかる。

第 4 章 幾何学的変換の基本、レジストレーション技術、

補間技術

第1章の画像モデルに述べているように、実際の画像は様々な要因で変形（幾何学的変換）され て観測される場合が多い。特に、マルチフレーム超解像技術において、各フレーム間の幾何学的 変換を推定し、補正を行う必要がある。また、計算機上で画像を幾何学的変換すると、抜け画素が 出るので、補間を行う必要がある。幾何学的変換のパラメータ推定（registration）と補間は、マルチ フレーム超解像技術の最も基本的要素技術である。本章では、幾何学的変換の基本、幾何学的変 換パラメータを推定するレジストレーション技術と補間技術について述べる。

4.1 幾何学的変換

超解像技術において、最も考慮される幾何学的変換は、平行移動（Translation）、回転 （Rotation）、拡大・縮小（Scale）である。 4.1.1 平行移動（Translation）

x軸方向の移動量をu、y軸方向の移動量をvとすると、点（x, y）から点（x’, y’）への移動は式（4.1）

のように表される。パラメータはuとvである。 ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ v u y x y x ' ' （4.1） 4.1.2 回転（Rotation） 原点を中心に、半時計周りに角度θだけ回転する変換は、式（4.2）のように表される。パラメータ はθである。 ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ y x y x

θ

θ

θ

θ

cos sin sin cos ' ' （4.2） 4.1.3 拡大・縮小（Scale） sx 、syは、それぞれx方向、y方向の拡大（または縮小）率（Scale）とすると、拡大・縮小変換は式 （4.3）のように表される。一般に、sx=sy =sであるので、パラメータはsである。 ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ y x s s y x y x 0 0 ' ' （4.3） 上記の3つの幾何学変換をまとめると、以下のように表すことができる。推定するパラメータはs、 θ、u、vの4つである。

⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ v u y x s y x

θ

θ

θ

θ

cos sin sin cos ' ' （4.4） なお、スキューや鏡映などを含めたアフィン変換の場合、推定するパラメータは6つになる。更に、 射影変換の場合、推定するパラメータは8つになる。マルチフレーム超解像技術において、サブピ クセルの精度で平行移動のパラメータを推定するのが一般的であるが、回転角度とスケールを入 れて推定する場合もある。

4.2 レジストレーション技術（registration）

複数枚（マルチフレーム）画像を用いて高解像度画像を生成するためには、フレーム間の微小 な位置ズレを推定し、レジストレーション（アライメント）を行う必要がある。本項では、Optical Flow、 最小二乗法、主成分分析法を用いた 3 つのレジストレーション技術について紹介する。 4.2.1 Optical Flow 法 マルチフレーム超解像技術における位置合わせでは、使用する画像は連続画像なので、画 像間の回転量は 0 と近似できることから、平行移動パラメータのみを求めるだけでも良いと考えら れる。物体が平行移動のみで濃度値変化はないものと仮定し、x 方向の移動量 u と y 方向の移動 量v は、以下の誤差関数の最小化により求めることができる11）_。 =

∑

+ + − y x y x I v y u x I v u E , 2 1 2( , ) ( , )) ( ) , ( （4.5） ここで、I1とI2はそれぞれ時刻t と時刻 t +Δt で観測された画像である。移動量 u と v が非常に小 さい場合、I2はテイラー級数展開で式（4.6）のように近似することができる。 y y x I v x y x I u y x I v y u x I ∂ ∂ + ∂ ∂ + ≈ + + , ) ( , ) ( , ) ( , ) ( ₂ 2 2 2 （4.6） また、式（4.5）の誤差関数は式（4.7）のように書き直すことができる。 2 , 1 2 2 2 ( , ) ) , ( ) , ( ) , ( ) , (

∑

_⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − ∂ ∂ + ∂ ∂ + = y x y x I y y x I v x y x I u y x I v u E （4.7） 式（4.7）をそれぞれ u と v について偏微分し、その偏微分が 0 となる u、v を求めると、以下のよう になる。 _⎟⎟=A−1b ⎠ ⎞ ⎜⎜ ⎝ ⎛ v u （4.8） ここで、 =

∑

∇ ⋅∇ y x T y x I y x I , 2 2 ) , ( ) , ( A （4.9） =−

∑

∇

(

−

)

y x y x I y x I y x I , 2 2 1 ) , ( ) , ( ) , ( b （4.10） 実際に移動量が十分に微小でない場合、上記の推定値に誤差が生じてしまう。Optical Flow を用いた微小変位の推定プログラムを List4_1.m に示す。入力画像（img2）を基準画像（img1）に 位置を合わせるプログラムである。その出力例を図 4.1 に示す。 なお、より高精度に移動量を推定するために、推定した移動量で画像 I2を移動させ、新しい画

像I2を用いて移動量を再度推定し、そのプロセスを反復的に行うことで、推定誤差を小さくするこ とができる。また、Gaussian Pyramid を用いると更に有効である。反復的に推定するプログラムは 第 5 章のマルチフレーム超解像技術に示す。 図 4.1 Optical Flow による位置合わせ List4_1.m %========================================================================% % オプティカルフローによって位置合わせパラメータを推定する List4_1.m % % img2 (入力画像) → img1(基準画像) % % function [RP] = optical_flow(img1,img2) % % 作成：立命館大学知的画像処理研究室 笹谷 聡 % % 作成日：2010.9.17 % % Copyright 立命館大学情報理工学部陳研究室 % %========================================================================% % clear all; clc; % 画像の読み込み img1 = double(imread('image\eia_lr_2_1.pgm','pgm')); % 基準画像 img2 = double(imread('image\eia_lr_2_2.pgm','pgm')); % 入力画像

[ysize xsize color_num]=size(img1); %画像サイズの取得 s = [1 2 1; 0 0 0; -1 -2 -1]; %ソーベルフィルタ XX = conv2(img1,s,'same'); %x軸方向への微分 YY = conv2(img1,s','same'); %y軸方向への微分 XX = reshape(XX,ysize*xsize,1); YY = reshape(YY,ysize*xsize,1); A = [XX YY]; IMG2 = reshape(img2,ysize*xsize,1); IMG1 = reshape(img1,ysize*xsize,1); b = A'*(IMG2-IMG1); f = -1*inv(A'*A)*b %オプティカルフローの計算 %img2をパラメータ分だけ変形した画像を作成