サンプル値制御理論による画像の解像度変換

Resolution Conversion of Digital Images
via Sampled-Data Control Theory○永原 正章 (京都大学), 山本 裕 (京都大学)

M. Nagahara∗ , Y. Yamamoto∗ ,Kyoto University∗
Abstract For handling digital images with diﬀerent kinds of resolution, we have to convert
their resolution. Conventionally, interpolation methods (e.g., bicubic) are used for this purpose.
In applying these methods, it is generally assumed that the original analog images are fully bandlimited up to the Nyquist frequency. As a result, the converted images may be blurred. To
solve this problem, we introduce the sampled-data control theory. By this theory, we can obtain
interpolation taking account of the characteristic of original analog images.

1

はじめに解像度の異なるデバイス間で画像データを取り扱った帰着できることを示す．

2

問題設定

り，画像のサイズを変換するためには解像度の変換が必要となる．画像の解像度変換には， 最近傍 (nearest

図 1 のサンプルレート変換器を考える．ここで自然数neighbour) 法や双線形 (bilinear)，双 3 次 (bicubic) 法

M1 と M2 は既約とする．まずサンプル周期 M1 h の離などの補間法がよく用いられている [2]．最近傍法は，補間する点の輝度値を最も近傍の輝度値とする方法であvd-φdψdud↑ M1 K(z) ↓ M2 -り，0 次ホールドと同じ働きをする．一方，双線形法と双 3 次法はそれぞれ，1 次および 3 次補間の一種であり，双線形法では近傍の 4 点，双 3 次法では近傍の 16 点を用いて補間を行う．最近傍法では単純に近傍の値をホールドするだけなので，処理は単純であるが画質は低い．
一方，双線形法や双 3 次法は計算は複雑であるが，画質は良くなる．
これらの手法は補間すべき点が 1 次や 3 次の多項式で単純に近似できるとしている．これは，アナログの原画像がナイキスト周波数以下に完全に帯域制限されていると仮定した上での議論である．しかし実際には，アナログ画像がそのような性質を持っているとは言いがたく，

図 1: Sample-rate converter
散時間信号 vd がアップサンプラ
yd [l],↑ M1 : yd → xd : xd [k] =0,k = M1 l, l = 0, 1, . . .
otherwise

によってサンプル周期 h の離散時間信号 φd に変換される．次に信号 φd はディジタルフィルタ K(z) によって処理され，その出力 ψd がダウンサンプラ↓ M2 : xd → ud : ud [k] = xd [M2 k]

もとのアナログ画像の特性を活かした補間処理が望まれる．実際，上記の双線形法や双 3 次法による解像度変換
では，原画像の高周波成分が除去され，ぼけた画像になりやすい．
そこで本研究では，画像の解像度変換にサンプル値制御理論 [3] を導入する．すなわち，サンプル値制御理論によって最終的にサンプル周期 M2 h/M1 の信号 ud に変換される．このサンプルレート変換器を用いれば，画像の解像度を M1 /M2 倍することができ，有理数倍の解像度変換が達成される．
一方，サンプルレート変換器は図 2 のようにも構成でにより，もとのアナログ画像の周波数特性を考慮し，アきる．ここでシステム K1 (z)(↑ M1 ) はインターポレータ

ナログ原画像と処理画像の誤差系の H ∞ ノルムが最小であり，また (↓ M2 )K2 (z) はデシメータである．この構となるような補間フィルタを設計する．サンプル値系で成の利点はフィルタ K1 (z) と K2 (z) が前述のフィルタある誤差系に対して，ノルムが等しい離散時間の誤差系K(z) に比べ低い次数で構成が可能である点である．本
研究では，サンプルレート変換器の設計において，フィを導くことによって，設計問題は離散時間 H ∞ 問題にルタ K(z) を直接設計するのではなく，インターポレー

- ↑ M1- K1(z)- K2(z)- ↓ M2-Hhw+eyK2

図 2: Sample-rate converter
ションフィルタ K1 (z) およびデシメーションフィルタ

インターポレータの設計問題

する．図 3 のブロック線図を考える．入力信号 w ∈ L2wF1↑ M1φdζ−K1ζd

れ，ダウンサンプラ ↓ M2 によってサンプル周期 M2 hの離散時間信号 ud に変換される．次に ud は 0 次ホールド HM2 h により連続時間信号 ζ となる．一方，離散号に変換し，時間 m2 h (m2 ∈ Z+ ) だけ遅れたものを z+ edyd

HM2 h

時間信号 wd を 0 次ホールド Hh によって連続時間信zdShySM1 hud

ef号 wd はまずディジタルフィルタ K2 (z) によって処理さまずインターポレータの (準) 最適設計問題を定式化

e−m1 hs↓ M2

F1

図 4: Error system T2

K2 (z) を別々に設計する．
2.1z

e−m2 hsとおく．誤差を e := z − ζ とし ef := F2 e とする．入−力 wd から 誤差 ef までのシステムを T2 とおくと，設計問題は次のように定式化される:問題 2 安定かつ厳密にプロパーな F2 (s) および，アッ

図 3: Error system T1

プサンプラ定数 M2 ， 遅れステップ m2 ，正数 γ2 が与えられたとき，がアナログローパスフィルタ F1 (s) を通り y を得る．こT2  := sup

こで F1 (s) はアナログ信号 y の周波数特性を示していwdる．次にこの信号はサンプル周期 M1 h で動作するサンプラ SM1 h によって離散時間信号 yd に変換される．この離散時間信号 yd はアップサンプラ ↑ M1 によってサ∈2

T2 wd 2< γ2wd 2(2)を満たすディジタルフィルタ K2 (z) を求めよ．

ンプル周期 h の離散時間信号 φd に変換される．φd は3

サンプル値制御理論によるフィルタ設計

ディジタルフィルタ K1 (z) を通り，離散時間信号 ζd と3.1

離散時間 H ∞ 問題への帰着なる．一方，連続時間信号 y が時間 m1 h (m1 ∈ Z+ )だけ遅れたものをサンプル周期 h で離散化した信号をzd とおく．そして zd と ζd との差を ed := zd − ζd とし，w から ed までのシステムを T1 とおくと，フィルタ K1 (z) の設計問題は以下のように定式化される．
問題 1 安定かつ厳密にプロパーな F1 (s) および，アップサンプラ定数 M1 ， サンプリング周期 h，遅れステップ m1 ，正数 γ1 が与えられたとき，

T1 w2
T1  := sup< γ1w∈L2 w2(1)を満たすディジタルフィルタ K1 (z) を求めよ．

2.2デシメータの設計問題

ここではデシメータの (準) 最適設計問題を定式化する．図 4 のブロック線図を考える．ここでは，入力信号はサンプル周期 h の離散時間信号 wd ∈ 2 である．信

まず，離散時間リフティングを定義する．

LN := (↓ N )[1 z . . . z N −1 ]T ,−1
L−1. . . z −N +1 ](↑ N )N := [1 z

離散時間 LTI システム P (z) = C(zI − A)−1 B + D に対して，システムの離散時間リフティングをLN (P ) := LN P L−1N
⎡NA
⎢
⎢
C
⎢
⎢
=⎢
⎢ CA
⎢
..
⎢.
⎣CAN −1

AN −1 B...

AB

D
CB..
.
N −2CAB..

.

..

.

..

...

CB.

B

⎤

⎥
0 ⎥⎥
⎥⎥⎥
⎥
⎥
⎦
D

とおく．また Fi (s) (i = 1, 2) の実現を (Ai , Bi , Ci , Di )定理の公式の詳細な計算は [5, 6] を参照されたい．

とおき (ただし D1 = 0)，次の作用素を定義する．

B1 : L2[0,h) → Rn ,B1 w :=

 h0eA1 (h−τ )定理 1 より，問題 1 および問題 2 は，それぞれ，離散時間 LTI フィルタ K1 および K2 を求める 1 ブロック

B1 w(τ )dτ,および K2 から K1 および K2 を得るには，次の公式

C2 : Rn → L2[0,h) , (C2 x)(t) := C2 eA2 t x,
D2 : L2[0,h) → L2[0,h) , t

(D2 w)(t) := D2 w(t) +C2 e0A2 (t−τ )の離散時間 H ∞ 最適化問題に帰着された．得られた K1を用いればよい [5]．

K1 (z) = [1 z −1 . . . z −M1 +1 ]K1 (z M1 ),B2 w(τ )dτK2 (z) = K2 (z M2 )[1 z . . . z M2 −1 ]T
上式より K2 の因果性が失われる可能性があるが，画像

このとき，次の定理が成り立つ．
定理 1

1. インターポレータ設計の誤差系 T1 に対して次の関係式が成り立つ．
T1  = (Dm1 − K1 S)Fd1 (3)

処理への応用では，フィルタは非因果的であってもかまわない．

3.2

サンプルレート変換器

前節で設計したフィルタ K1 (z)，K2 (z) を用いて，図 1
ここでのサンプルレート変換器を図 2 のように構成し，K(z) =Dm1 := LM1 (zFd1 := LM1
Ad1 := eA1 h ,−m1),

S := [1 0 . . . 0]

Ad1

Bd1C0

K1 (z)K2 (z) とおく．このとき，誤差系
T := F2 e−mhs − HM2 h (↓ M2 )K(z)(↑ M1 )SM1 h F1,を考える．ここで，m = m1 + m2 である．設計された

T
Bd1 Bd1= B1 B1∗ ,K1 (z)，K2 (z) に対して T1  = γ1 ，T2  = γ2 とおく．

K1 := LM1 (K1 )[1 0 . . . 0]T

誤差系 T のノルムに対して，次の命題が成り立つ．

2. デシメータ設計の誤差系 T2 に対して次の関係式

命題 1が成り立つ．

T  ≤ γ1 γ2 + α1 γ1 + α2 γ2 + α3T2  = Fd2 (Dm2 − H K2 )

(4)

ここで

α1 = F2 Hh ,

ここで

Dm2 := LM2 (z −m2 ),
Ad2

Bd2

Cd2

Ad2 := eA2 h ,Dd2 hBd2 :=0,および AB ≤ AB を用いれば (5) が成り立つことがわかる．

eA2 τ dτ B,

T1 T2 = F2 e−m2 hs Hh T1 + T2 Sh e−m1 hs F1 − T∗

K2 := [1 0 . . . 0]LM2 (K2 )証明

次式と誘導ノルムの性質 A + B ≤ A + B

証明[Cd2 Dd2 ] [Cd2 Dd2 ] = [C2 D2 ] [C2 D2 ],T+ F2 e−mhs (I − Hh Sh )F12

1. サンプル値系 T1 ∈ B(L2 , 2 ) に対して， この命題より，サンプルレート変換器の性能 T  が評
LTI 離散時関系 Td1 ∈ B(2 , 2 ) が存在して，T1 2 = 価できる．すなわち，T  < γ となるようなサンプル∗T1 T1∗  = Td1 Td1
 = Td1 2 が成り立つ．ここレート変換器を設計したい場合，まず T1 に対して最適で Td1 = (Dm1 − K1 S)Fd1 である．
化を行い，T  = γ が得られたとすると，この γ を12. 上と同様にして，Td2 = Fd2 (Dm2 − H K2 ) に対し∗て，T2 2 = T2∗ T2  = Td2
Td2  = Td2 2 が成り立つ．

2α2 = Sh F1 ,α3 = F2 (I − Hh Sh )F1 

H := [1 1 . . . 1]T

Fd2 := LM2(5)

11

用いて

γ1 γ2 + α1 γ1 + α2 γ2 + α3 < γを満たすような γ2 に対して T2  < γ2 となるようにT2 を最適化すればよい．

3.32 次元フィルタの構成

前節で設計したフィルタは 1 次元である．これを画像処理に適用するために，処理を 2 次元化する必要がある．そこで，サンプルレート変換の 2 次元処理を分離可能 (separable)[4] なものとする．すなわち 2 次元フィル
タ K2D をK2D (z1 , z2 ) = K(z1 )K(z2 )図 6: Bicubic interpolation

で構成する

4

設計例
ここでは，提案手法による解像度変換と従来法の Bicu-bic による解像度変換を比較する．設計のパラメータは以下のとおりである．
F1 (s) =1,(s + 1)2F2 (s) =h = 1,

M1 = 2,m1 = 2,m2 = 31,
s+1

M2 = 3,図 5 のオリジナル画像に対して，M1 /M2 = 2/3 の解像度変換を行う．図 6 は Bicubic 法による解像度変換により得られた画像である．この処理には Photoshop CS[1] を用いた．いっぽう，提案手法による解像度変換を
図 7 に示す．従来法のぼけが提案手法によって改善されていることがわかる．

図 7: Sampled-data H ∞ interpolation

参考文献[1] Adobe Systems, Photoshop CS, Adobe Systems(2004)[2] K. R. Castleman, Digital Image Processing, Prentice Hall, (1996)[3] T. Chen and B. A. Francis, Optimal SampledData Control Systems, Springer, New York (1995)[4] W.-S. Lu and A. Antoniou, Two-dimensional Digital Filters, Marcel Dekker (1992)[5] M. Nagahara and Y. Yamamoto, “A new design
for sample-rate converters,” Proc. of 39th IEEE
Conf. on Decision and Control, pp. 4296–4301(2000)

図 5: Original image[6] M. Nagahara and Y. Yamamoto, “Optimal designof fractional delay ﬁlters,” Proc. of 42nd IEEE

5

おわりに
本研究は画像処理における解像度変換のためのサンプ

Conf. on Decision and Control, pp. 6539–6544(2003)

ルレート変換器の設計をサンプル値 H ∞ 設計問題とし[7] P. P. Vidyanathan, Multirate Systems and Filterて定式化し，誤差系をノルムに関して等価変換すること

Banks, Prentice Hall, Englewood Cliﬀs, (1993)

により離散時間 H ∞ 最適化問題へと帰着した．また，数値例により，従来法よりもより良い処理が行えることを示した．