双直交指向性離散コサイン変換の設計
Design of Biorthogonal Directional Discrete Cosine Transforms
市田智大
†京地清介
†鈴木大三
††田中雄一
†††† 北九州市立大学大学院国際環境工学研究科
†† 筑波大学システム情報系
††† 東京農工大学大学院工学研究院
Tomohiro ICHITA
†Seisuke KYOCHI
†Taizo SUZUKI
††Yuichi TANAKA
††††The University of Kitakyushu
††University of Tsukuba
†††Tokyo University of Agriculture and Technology
アブストラクト 本論文では，指向性離散コサイン変換
（DDCT: Directional Discrete Cosine Transform）を拡張 した双直交指向性離散コサイン変換（BDDCT: Biorthog-onal DDCT）を提案する．DDCT は，直交行列である離 散コサイン変換（DCT: Discrete Cosine Transform）と 離散サイン変換（DST: Discrete Sine Transform）の並列 可分型二次元変換によって構成され，画像のテクスチャが 持つ多方向の指向性を解析できる利点を持つ．本論文で は，DCT と DST が変調フィルタバンクの特殊例である とみなし，そのプロトタイプフィルタのパラメータを与 える正則行列を DCT/DST に乗じることで，指向性を維 持したまま，DDCT よりも高い性能（符号化利得等）を 有する双直交変換を実現する．BDDCT の性能評価とし て画像復元に適用し，DDCT に比べて高い復元性能を有 することを示す． 1 はじめに 近年，画像の解析や加工，復元技術のための画像疎表現 [1], [2] を与える画像変換として，詳細なテクスチャが有す る多様な指向性を捉える Curvelet [3]，Contourlet [4] や， 指向性フィルタバンク [5]，複素ウェーブレット・フィル タバンク [6], [7] 等の指向性変換が提案され，有効性が示 されている． しかし，従来の指向性変換は計算負荷を伴うため，計算 資源の十分でないデバイスで扱う際に問題が生じる．例 えば，Contourlet における二次元フィルタは一次元フィ ルタに比べて大きな演算量を要する．また複素ウェーブ レットのように，一次元フィルタによって実現できる手法 も存在するが，atom（基底・フレームの要素）が互いに 重複しているため（重複変換），並列処理化が困難である こと，あるいは画像全体のメモリアクセスが必要になる 点が問題として挙げられる． 前述の重複変換とは異なり，離散フーリエ変換（DFT: Discrete Fourier Transform）[8]，離 散 コ サ イ ン 変 換 （DCT: Discrete Cosine Transform）[9] といったブロッ ク変換は，atom が互いに重複しておらず，各々のブロッ クで変換が完結するため，並列処理化やメモリ使用効率に おいて，重複変換よりも有利であり，画像符号化をはじめ， 様々なアプリケーションに適用されている．しかし，指向 性を有するブロック変換に関しては重複変換と比較して 疎表現効率が劣るため，これまで十分な検討がされていな かった．そこで DFT や DCT に比べ，より豊富な指向性 （atom の方向数）を持つブロック変換として，指向性離 散コサイン変換（DDCT: Directional DCT）[10] が提案 されている．この変換は，DCT と DST（DST: Discrete Sine Transform）の並列可分型二次元変換により構成さ れ，ブロック変換の利点である演算量，並列演算，メモリ 使用効率面での優位性を保持しつつ，高い画像解析・処 理能力を実現している． 本研究では DDCT を拡張した双直交指向性離散コサ イン変換（BDDCT: Biorthogonal DDCT）を設計し， DDCT の性能を向上させる．DDCT を構成する DCT， DST は変調フィルタバンクの特殊例（プロトタイプフィ ルタの係数が全て 1）であるとみなすことができるため， DCT/DST にプロトタイプフィルタのパラメータを与え る正則行列を乗じて双直交化を行うことで，指向性を保 持したまま DDCT よりも高い性能（阻止域減衰量，符号 化利得等 [11]）を実現する． 以降，本論文を以下のように構成する．2 節では，従来の 代表的なブロック変換，指向性ブロック変換として DCT， DFT，DDCT を述べる．3 節で提案法である BDDCT を 導入する．性能評価として BDDCT を画像復元（欠損画 像補間）の実験に適用する．最後に 5 節で本論文を結ぶ． 第31回信号処理シンポジウム 2016年11月8日∼10日（関西大学）
(a) DCT (b) DFT （real part） 図 1 DCT, DFT の atom B(k1,k2) ∈ RM×M （M = 8）． 1.1 表記 R は実数の集合を表す．N 次元の実数係数ベクトル空 間をRN と表記する．N r行 Nc 列の実数係数行列の全 体をRNr×Nc とする．小文字・大文字の太字はそれぞれ ベクトル・行列を表す．j := √−1. ΩN :={0, . . . , N}， ΩN1,N2 :={N1, . . . , N2}．I は単位行列，A⊤は A の転置 行列，xi，[x]iはベクトル x の第 i 成分，Xi,j，[X]i,jは 行列 X の第 (i, j) 番目の要素とする．X(i,j) ∈ RM×M は行列 X ∈ RM L1×ML2 (i ∈ Ω L1−1, j ∈ ΩL2−1) の (i, j) 番目の小ブロック（M × M）とする．行列 X ∈ RNr×Nc に対して，vec(X)∈ RNrNcは X のベクトル化
（xNrj+i= Xi,j）とする．bvec(X)∈ R
(M L1)(M L2)は行列
X∈ RM L1×ML2のブロック毎のベクトル化（bvec(X) =
[vec(X(0,0))⊤ vec(X(1,0))⊤ . . . vec(X(L1−1,L2−1))⊤]⊤） とする．⊗ はクロネッカー積を表す． 2 前提事項 2.1 DCT，DFT M × M の DCT を FCとするとき，FCは次式で定義 される． [FC]k,n= αk √ 2 M cos(θk,n), θk,n= π Mk ( n +1 2 ) . (1) ここで，k は周波数変数，n は時間変数 (k, n∈ ΩM−1) と し，α0 = √12(k = 0)，αk = 1(k ̸= 0) とする．DCT は 信号の疎表現に優れており，圧縮等に適用されているが， 画像の詳細なテクスチャが有する多様な指向性を効率よ く表現することができないことが難点である．これはす なわち，二次元 DCT の基底に含まれる atom を B(k1,k2) n1,n2 とするとき， B(k1,k2) n1,n2 = αk1αk2 2 M cos (θk1,n1) cos (θk2,n2) , (2) のように指向性が無いことに起因する（図 1(a)）．一方， 指向性を有する代表的なブロック変換として DFT が挙げ られる．DFT を FFとするとき，FFは [FF]k,n= 1 √ Me jφk,n, φ k,n= 2π Mkn (3) のように表現され，二次元 DFT の atom は， B(k1,k2) n1,n2 = 1 Me j(φk1,n1+φk2,n2) (4) となることから，指向性をもつことがわかる（図 1(b)）． DFT は，指向性を有する atom を形成できるが，問題 点がある．例えば，図 1(b) の k1, k2= 1 と k1, k2= 7 を 見ると，同一帯域を示す atom が重複して含まれているこ とにより豊富な指向性が得られないため，画像解析の効 率が低下すると考えられる． 2.2 DDCT 2M2× M2の DDCT を F Dとし，定義を以下に示す． FD:=P⊤1 1 √ 2I1 1 2I2 − 1 2I2 1 √ 2I1 1 2I2 1 2I2 P1 [ FC⊗ FC FS⊗ FS ] , (5) ここで，FCは (1) の DCT であり，FS は次式で定義さ れる． [FS]k,n= √ 1 Msin (π M(M ) ( n +12)) (k = 0) √ 2 Msin (π Mk ( n +12)) (k̸= 0) . (6) これはすなわち，（type-II）DST を行に関して 1 行周期 シフトした行列である．I1 ∈ R(2M−1)×(2M−1)，I2 ∈ R(M−1)2×(M−1)2 は単位行列，P1= diag( bP1, bP1)，bP1 ∈ RM2×M2は置換行列とする（ bP⊤ 1Pb1= I）．具体的に bP1 は，二次元のサブバンドインデックスが k1 = 0 または k2= 0 である 2M− 1 個の DCT（または DST）係数を先 頭に，その他の k1̸= 0 かつ k2̸= 0 であるような (M −1)2 個の係数をその後ろに並べ替える行列である． DDCT の変換のフローを図 2 に示す．DCT と DST の 和/差である二つの可分型ブロック変換を並列処理するた め，DCT や DFT のブロック変換と比較すると，演算量 は多くなる．しかし，Curvelet や Contourlet 等の重複変 換と比較すると，演算量は少ない． 最後に，DDCT の atom の指向性を確認する．DDCT は定義より，FDを用いて F⊤DFD= I が成り立つ．これ は，DDCT がパーセヴァルフレームを形成することを意味 し，パーセヴァルフレームの atom は F⊤Dを調べればよい． DDCT の atom には，二次元 DCT と二次元 DST の atom の他に，二次元 DCT と DST の atom の和/差が含まれて
2D DCT 2D DST 'LUHFWLRQDO 6XEEDQG6LJQDOV 18 17 3 4 4 19 17 24 23 11 12 25 26 27 28 29 30 31 32 26 27 28 29 30 31 32 10 11 23 12 9 21 20 7 8 22 6 6 7 20 8 5 5 1 2 3 2 1 19 18 9 10 25 24 22 21 13 13 14 14 15 15 16 16 18 17 3 4 4 19 17 24 23 11 12 25 26 27 28 29 30 31 32 26 27 28 29 30 31 32 10 11 23 12 9 21 20 7 8 22 6 6 7 20 8 5 5 1 2 3 2 1 19 18 9 10 25 24 22 21 13 13 14 14 15 15 16 16 図 2 DDCT の変換フロー (M = 4) と二次元周波数分 割形状．灰色は上側と下側から出力されるサブ バンドを表す． (a) B(k1,k2,1) (b) B(k1,k2,−1) 図 3 DDCT の atom（M = 8）． おり，atom の和/差をそれぞれ B(k1,k2,1), B(k1,k2,−1) ∈ RM×M (k 1, k2)∈ Ω1,M−1× Ω1,M−1とすると，それらは B(k1,k2,±1) n1,n2 = C (k1,k2) n1,n2 ± S (k1,k2) n1,n2 = 2 M cos (θk1,n1∓ θk2,n2) , (7) のように表すことができる．図 3 に DDCT の atom を示 す．(7) と図 3 より，DDCT の atom は指向性を有してい ることがわかる．サブバンドインデックスが k1 = 0 ま たは k2 = 0 であるサブバンドは指向性を持たないため， DDCT の行列のサイズが 2M2× M2であるときの方向選 択数は 2(M− 1)2となる． 3 双直交指向性離散コサイン変換 前節では DDCT の構成について述べた．本節では， DDCT の指向性を保持しながら，更なる性能向上を目的 とした BDDCT の構成を説明する．DDCT は，DCT と DST の並列可分型二次元変換であり，使用する DCT と DST は変調フィルタバンクの特殊例（プロトタイプフィル タの係数が全て 1）であるとみなすことができる．DCT， DST のプロトタイプフィルタの係数をカスタマイズするこ とで，阻止域減衰量や符号化利得 [11] 等の指標に対して性 能の高いフィルタバンクを用いた並列可分型二次元変換を 設計する．以下で，DCT/DST のプロトタイプフィルタの 係数を変更して得られる双直交 DCT（BDCT: Biorthog-onal DCT），および双直交 DST（BDST: BiorthogBiorthog-onal DST）を定義する．パラメータ p = [ p0 · · · p2M−1 ] に対して行列 P を P = pM +(M/2)−1 p(M/2)−1 . .. . .. pM p0 p2M−1 pM−1 . .. . .. pM +(M/2) pM/2 (8) と設定する．(8) を使用することにより，BDCT を FBC FBC = FCP (9) と定義する．次に行列 Q を Q = pM +(M/2)−1 −p(M/2)−1 . .. . .. pM −p0 −p2M−1 pM−1 . .. . .. −pM +(M/2) pM/2 (10) と設定する．(10) を使用することにより，BDST を FBS FBS= FSQ (11) のように定義する．(9)，(11) の行列 P，Q によって，以 下のようにプロトタイプフィルタが与えられることが分 かる． [FBC]k,n= fk,n[FC]k,n (12) [FBS]k,n= fk,n[FS]k,n (13) ただし， fk,n= pM +(M/2)−1−n+ pM +(M/2)+n (k : even) (14) fk,n= pM +(M/2)−1−n− pM +(M/2)+n (k : odd) (15) である．以上の (9)，(11) を使用し，DDCT の FCを FBC に，FSを FBSに置き換えることで，BDDCT が構成され る．DDCT はパーセヴァルフレームであり，F⊤DFD = I であったが，BDDCT はプロトタイプフィルタ係数行列 を乗じたため，F−1BDFBD= I となるような双直交変換と
(a) B(k1,k2,1) (b) B(k1,k2,−1) 図 4 BDDCT の atom（M = 8）． なる．ここで，F−1BDは，以下のように計算できる． F−1BD:= [ (FBC⊗ FBC)−1(FBS⊗ FBS)−1 ] × P⊤ 1 1 √ 2I1 1 2I2 1 2I2 1 √ 2I1 1 2I2 − 1 2I2 P1. (16) M=8 の場合の BDDCT の atom を図 4 に示す．図 4 より 明らかに，DDCT の指向性を保持したまま，双直交への 一般化が実現できていることが分かる． 4 実験 本節では，提案法である BDDCT の性能を欠損画素補間 によって評価する．比較に使用した画像として Barbara，
Zoneplate，M andrill，Lena の 4 枚（図 5(a)–(d) 参照）
を用いた．欠損画素補間の性能比較のため，観測画像と して 4 枚の画像に 35，45，55%のそれぞれの確率でラ ンダムに欠損させたものを使用する（図 5(e)–(h) は原画 像の 45%の画素で構成された画像）．観測画像から原画 像を変換係数の L1ノルムの正則化を用いて復元を行う （詳細は付録参照）．比較する対象として，拡張元である DDCT の他に DCT，DFT，Discrete Hartley Transform （DHT）[12] を使用した．ブロックサイズは M = 8 に設 定した．BDDCT のプロトタイプフィルタの設計に関して は，画像の疎表現の効率を評価する符号化利得を，MAT-LAB の関数 “fminunc” を用いて最大化することで決定 した．DDCT，および最適化結果の BDDCT の符号化利 得を表 1 に，DCT/DST/BDCT/BDST の周波数応答を 図 6 に示す．また，文献 [10] で指摘されているように， DDCT/BDDCT では DC 漏れが生じるため，画像の各ブ ロックの平均を計算する行列 M = M0⊗ M0（ただし， [M0]k,n= 1/M ）を用いて，FBD(I− M) とした変換を 画像復元に適用する． 表 1 と図 6 に実験結果を示す．ほとんどの画像と欠損率 において，提案法である BDDCT は DCT，DFT，DHT， DDCT よりも優れた復元精度を示した．BDDCT の符号
(a) Barbara (b) Zoneplate (c) Mandrill (d) Lena
(e) (f) (g) (h) 図 5 (a)–(d): 原画像 （256× 256），(e)–(h): 観測画 像 （原画像の 45%の画素数）． 表 1 符号化利得 DDCT BDDCT 11.93 12.52 化利得が DDCT よりも高い，これはすなわち疎表現効率 が改善されたことに起因している． 5 結論 本論文では，画像解析のための指向性変換である DDCT の性能向上のために拡張を施した BDDCT を提案した． 拡張前の DDCT 同様に，ブロック処理が可能であり，従 来の重複フレーム処理よりも演算面で効率が良い．更に， DDCT の持つ方向選択数を保持したまま，更に良い性能 で指向性変換を実現できることが分かった． 実験では，BDDCT の性能を欠損画素補間を用いて評 価を行った．実験結果では，拡張前の DDCT よりも拡張 後の BDDCT の方が優れた復元精度が得られた． 付録 欠損画素補間の詳細なアルゴリズム 実験に使用する画像復元問題を以下のように定式化する． x∗= arg min x∈RNr Nc ρ∥Fx∥1+ ιC[0,1](x) + Fv(Φx), (17) ただし，ρ > 0，x = bvec(X)，X ∈ RNr，Φ は劣化過 程を示し，ιA(x) は指数関数1とする．F はブロック毎 の BDDCT によるため，文献 [10] と同様に，画像の各ブ ロックの平均を計算する行列 M = M0⊗ M0（ただし， [M0]k,n= 1/M ）を用いて，C[0,1]は各要素が [0, 1] に含 まれるベクトルの集合とする． Fv ∈ Γ(RN) は観測画像 v に対するデータ忠実項とする． 本論文では欠損画素補 1集合 A における指示関数は ι A(x) = 0, (x ∈ A), ιA(x) = ∞, (x /∈ A) と定義される．
(a) (b) (c) (d) 図 6 周波数応答（周波数:[0, 2π]）:(a)DCT，(b)DST， (c)BDCT，(d)BDST．（M = 8）． 間に適用するため，データ忠実項は ι{v}(x) とする．ただ し Φ (17) は，単位行列の対角成分において，欠損した画 素に対応する場所の 1 を 0 に変えることで定義する．本実 験では，パラメータ γ1，γ2(19) ，ρ (17) は，0.01，12γ11， 0.5 とし，終了条件は∥x(n+1)− x(n)∥2≤ 0.01 とした．
(17) を解くために， primdual splitting （PDS） al-gorithm [13] を用いる．ここで，以下の凸適化問題 x⋆∈ arg min x∈Rp f (x) + g(Lx), (18) を考える．ただし f ∈ Γ0(Rp)， g ∈ Γ0(Rq) (Γ0(RN) はRN 上の下半連続な真凸関数の集合とする [14])， L∈ Rp×q．このとき，最適解 x⋆は以下のように導出できる： xk+1:= proxγ1f[xk− γ1L ∗zk] zk+1:= proxγ2g∗[zk+ γ2L(2xk+1− xk)] , (19) ただし，prox は近接写像 [14]，h∗ は h の共役関数 [14]， L∗は L の随伴作要素とする．この実験の場合，関数 g， h，行列 L (18) は以下のように設定する： g(x) = ιC[0,1](x), h([z ⊤ 1 z⊤2]⊤) = ρ∥z1∥1+ ι{v}(z2), z1= Fx, z2= Φx, L = [ F Φ ] . (20) 最終的なアルゴリズムは Algorithm 1 にまとめる2． 2x ∈ RN に対して，[prox γ∥·∥1(x)]i = sign(xi) max{|xi| − γ, 0} （soft-thresholding），proxι C[0,1](x) は [0, 1] へのクリッピン グ， proxι {v}(x) = v，ただし v∈ R Nは観測画像． (a) DCT (b) DFT (c) DDCT (d) BDDCT 図 7 結果画像 （(a)，(b)，(c)：従来法，(d)：提案法）． Algorithm 1 (17) の最適化アルゴリズム 1: n = 0，x(0), z(0) 1 , z (0) 2 , γ1, γ2の設定. 2: while 終了条件が偽な場合 do 3: x(n+1)= prox γ1ιC[0,1](x (n)− γ 1(F⊤z(n)1 + Φ⊤z (n) 2 )) 4: t(n)1 = z(n)1 + γ2F(2x(n+1)− x(n)). 5: t(n)2 = z(n)2 + γ2Φ(2x(n+1)− x(n)). 6: ˆt(n)1 = prox1 γ2ρ∥·∥1 ( 1 γ2t (n) 1 ) . 7: ˆt(n)2 = prox1 γ2ι{v} ( 1 γ2t (n) 2 ) . 8: z(n+1)k = t(n)k − γ2ˆt(n)k (k = 1, 2). 9: n = n + 1. 10: end while 11: u(n)を出力. 参考文献
[1] J.-L. Starck, F. Murtagh, and J. Fadili, Sparse
Im-age and Signal Processing: Wavelets, Curvelets, Morphological Diversity, New York, NY, USA:
Cambridge University Press, 2010.
[2] M. V. Afonso, J. M. B.-Dias, and M. A. T. Figueiredo,“An augmented lagrangian approach to the constrained optimization formulation of imag-ing inverse problems,”IEEE Trans. Image Process., vol. 20, no. 3, pp. 681–695, Mar. 2011.
[3] J. L. Starck, E. J. Cand‘es, and D. L. Donoho, “ The curvelet transform for image denoising, ”IEEE
Trans. Image Process., vol. 11, no. 6, pp. 670–684,
表 2 復元誤差 Image: Barbara PSNR [dB] DCT DFT DHT DDCT BDDCT 35% : 10.77 29.06 28.29 27.56 29.06 29.69 45% : 9.689 26.62 26.35 25.64 27.27 27.68 55% : 8.812 24.18 24.26 23.74 25.35 25.65 Image: Mandrill PSNR [dB] DCT DFT DHT DDCT BDDCT 35% : 10.10 27.69 27.94 27.29 28.24 28.46 45% : 8.981 25.94 26.16 25.60 26.78 26.98 55% : 8.134 24.36 24.60 24.05 25.43 25.51 Image: Lena PSNR [dB] DCT DFT DHT DDCT BDDCT 35% : 11.73 30.27 28.58 28.04 29.37 29.88 45% : 10.65 27.49 26.40 25.90 27.15 27.46 55% : 9.777 25.50 24.63 24.09 25.54 25.80 Image: Zoneplate PSNR [dB] DCT DFT DHT DDCT BDDCT 35% : 8.803 18.51 19.38 18.51 21.01 22.38 45% : 7.727 16.05 17.36 16.48 18.53 19.59 55% : 6.866 13.53 15.41 14.51 16.08 16.91
[4] M. N. Do and M. Vetterli, “ The contourlet trans-form: an efficient directional multiresolution image representation, ”IEEE Trans. Image Process., vol. 14, no. 12, pp. 2091–2106, Dec. 2005.
[5] W. Q. Lim, “ Nonseparable shearlet transform, ”
IEEE Trans. Image Process., vol. 22, no. 5, pp.
2056–2065, May 2013.
[6] S. Kyochi and M. Ikehara, “ A class of near shift-invariant and orientation-selective transform based on delay-less oversampled evenstacked cosine-modulated filter banks, ”IEICE Trans. Fundam., vol. 93, no. 4, pp. 724–733, Apr. 2010.
[7] S. Kyochi, T. Uto, and M. Ikehara, “ Dual-tree complex wavelet transform arising from cosine-sine modulated filter banks, ”in Proc. IEEE Int. Symp.
Circuits and Syst. (ISCAS), May 2009.
[8] J. W. Cooley and J. W. Tukey,“ An algorithm for the machine calculation of complex Fourier series, ”
Mathematics of Computation, vol. 19, pp. 297–301,
Apr. 1965.
[9] Z. Wang, “ On computing the discrete fourier and
cosine transforms,”IEEE Trans. Acoust. Speech
Sig-nal Process., Vol. 33, No. 4, pp. 1341–1344, 1985.
[10] 京地 清介，鈴木 大三，田中 雄一，“ 指向性離散コサ
インフレームの設計 ”，電子情報通信学会技術研究 報告，July 2016．
[11] G. Strang and T. Q. Nguyen, Wavelets and Filter
Banks, Wellesley, MA: Wellesley-Cambridge, 1997.
[12] R. N. Bracewell,“The fast hartley transform,”Proc.
of the IEEE, vol. 72, no. 8, pp. 1010–1018, Aug.
1984.
[13] L. Condat,“A primal-dual splitting method for con-vex optimization involving lipschitzian, proximable and linear composite terms, ”J. Optimization
The-ory and Applications, Dec. 2013.
[14] H. H. Bauschke and P. L. Combettes, Convex
Analysis and Monotone Operator Theory in Hilbert Spaces, New York, NY, USA: Springer- Verlag,
![表 2 復元誤差 Image: Barbara PSNR [dB] DCT DFT DHT DDCT BDDCT 35% : 10.77 29.06 28.29 27.56 29.06 29.69 45% : 9.689 26.62 26.35 25.64 27.27 27.68 55% : 8.812 24.18 24.26 23.74 25.35 25.65 Image: Mandrill PSNR [dB] DCT DFT DHT DDCT BDDCT 35% : 10.10 27.69 27.94](https://thumb-ap.123doks.com/thumbv2/123deta/6492381.658474/6.892.80.430.124.546/表復元誤差276855881224182426237425Image.webp)
