Proceedings of APSIPA Annual Summit and Conference 201712 - 15 December 2017, Malaysia

Multiscale Directional Transforms based on
Cosine-Sine Modulated Filter Banksfor Sparse Directional Image Representation
Yusuke Nomura∗ , Ryutaro Ogawa∗, Seisuke Kyochi∗, and Taizo Suzuki†∗ The University of Kitakyushu, Fukuoka, Japan

E-mail: s-kyochi@kitakyu-u.ac.jp† University of Tsukuba, Ibaraki, Japan

E-mail: taizo@cs.tsukuba.ac.jp

Abstract—This paper proposes multiscale directional transforms (MDTs) based on cosine-sine modulated ﬁlter banks(CSMFBs). Sparse image representation by directional transforms is necessary for image analysis and processing tasks
and has been extensively studied. Conventionally, cosine-sinemodulated ﬁlter banks (CSMFBs) have been proposed as one ofseparable directional transforms (SepDTs). Their computational
cost is much lower than non-SepDTs, and they can work better
than other SepDTs, e.g., dual-tree complex wavelet transforms(DTCWTs) in image processing applications. One drawbackof CSMFBs is a lack of multiscale directional selectivity, i.e.,
it cannot provide multiple scale directional atoms as in the
DTCWT frame, and thus ﬂexible image representation cannot be
achieved. In this work, we show a design method of multiscale
CSMFBs by extending modulated lapped transforms, which area subclass of CSMFBs. We conﬁrm its effectiveness in nonlinear
approximation and image denoising as a practical application.

I. I NTRODUCTION
Directional transforms (DTs) such as dual-tree complex wavelet transforms (DTCWTs) [?], [?], Contourlet [?],
Curvelet [?], and so on, have been widely studied for their efﬁciency of sparse image representation. Among DTs, DTCWTs
have gathered much attention due to their low computational
complexity [?], [?]. While two dimensional non-separable(non-Sep) DTs usually require high computational complexity,
DTCWTs, which is a subclass of 2D SepDTs, can save computational complexity by a parallel 2D separable implementationof two ﬁlter banks (FBs), i.e., ﬁltering along vertical andhorizontal directions. Thanks to DS, it can be successfully
applied to various kinds of practical image processing, suchas image denoising, image analysis, and image compression[?].
As well as the DTCWT, cosine-sine modulated ﬁlter banks(CSMFBs) have been shown as another SepDTs [?], [?], [?].
One problem of DTCWTs is the difﬁculty in designing M channel DTCWT with good ﬁlter performances, e.g., stopband attenuation or coding gain, since the fractional delay
requirement between two FBs has to be approximated. On
the other hand, CSMFBs require NO fractional delay (just theThis work was supported in part by JSPS Grants-in-Aid (16K18100 and
17K14683).

978-1-5386-1542-3@2017 APSIPA

modulation of a prototype ﬁlter), and thus it is easier to design
M -channel CSMFBs with better performance than M -channel
DTCWTs.
One problem on the conventional CSMFBs is a lack of multiscale directional selectivity, i.e., it cannot provide multiscaledirectional atoms in their frame. DTCWTs, which is known
as the class of multiscale directional transforms (MDTs), cansparsely express an input image with multiscale directional
atoms by cascading the ﬁrst and the second FBs to eachlowpass subbands. On the other hand, as it will be explained,
if CSMFBs are applied to images in the same way of theDTCWT, multiscale directional image representation cannot
be achieved.
This paper attempts to design MDTs based on CSMFBs. For
that, modulated lapped transforms (MLTs), which is a subclassof the CSMFB, are used. Although MDTs cannot be realized
by using MLTs directly, this paper clariﬁes that a shifted MLT
can realize the MDTs. The resulting transform is termed as“multiscale directional MLT (MDMLT)”. Experimental resultsof nonlinear approximation and image denoising show that the
MDMLT can reduce the ringing artifact compared with theconventional CSMFBs.
The rest of this paper is organized as follows. Sec. II brieﬂyreviews about the CSMFB. The MDMLTs are presented inSec. III. Sec. IV evaluates the MDMLT in nonlinear approximation and image denoising. Finally, this paper is concluded
in Sec. V.
−nNotations: H(z) is deﬁned as H(z) :=.
n h(n)z
−jωn
H(ω)
:=
h(n)e
,
H(ω)
:=
H(ω
,
ω
)
=
1
2
n

−j(ω1 n1 +ω2 n2 ).
Hk1 ,k2 (ω):=
n1 ,n2 h(n1 , n2 )ej 2π
N
Hk1 (ω1 )Hk2 (ω2 ). WN = e . The M × L FB means
the M -channel ﬁlter bank with the ﬁlter length of L. IN ,JN ∈ RN ×N are the identity and the reversal identity,
respectively. diag(A0 , . . . , AN −1 ) is a diagonal matrix. ⊗denotes the Kronecker product.
II. P RELIMINARIESA CSMFB is constructed two maximally decimated perfectreconstruction FBs, as illustrated in Fig. 1. The ﬁlter coef-

1830

APSIPA ASC 2017

Proceedings of APSIPA Annual Summit and Conference 201712 - 15 December 2017, Malaysia/HYHO/HYHO
//

)LUVW
)%
6HSDUDEOH
'
7UDQVIRUP

/+
+/
++

//
)LUVW
)%
6HSDUDEOH
'
7UDQVIRUP





//
6HFRQG)%
6HSDUDEOH
'
7UDQVIRUP

/+
+/
++

/+
+/
++

//





6HFRQG)%
6HSDUDEOH
'
7UDQVIRUP

/+
+/
++

'LUHFWLRQDO
6XEEDQGV(2)ﬁcients of {Hk (z), Hk (z)} of CSMFBs are expressed asfollows [?]:(1)hk (n) := 2p(n)ck,n = 2p(n) cos(θ(k, n))(2)hk (n) := 2p(n)sk,n = 2p(n) sin(θ(k, n))πθ(k, n) =
M




1
N −1
k+
n−
+ ϕk
2
2

(1)
(2)
(3)

where k = 0, . . . , M − 1, N = 2mM is the ﬁlter length,
n = 0, . . . , N − 1, ϕk = (−1)k π4 and p(n) is the prototype
ﬁlter. Equivalently, these ﬁlters can be represented in the zdomain as follows:
(1)
Hk (z) := ck Uk (z) + ck Vk (z)(2)
Hk (z) := −j(ck Uk (z) − ck Vk (z))




k+ 1
−(k+ 1 )Uk (z) = P zW2M 2 , Vk (z) = P zW2M 2
(1)

(1)

(2)

'LUHFWLRQDO
6XEEDQGV

as illustrated in Fig. 2. Finally, directional subband signals
with twice number of samples of the input image can beobtained. For example, in the case of the number of channel[1][1]
M = 2, in order to obtain the Y1,1,1 (ω) and Y1,1,−1 (ω), thefollowing operation is performed to input image;


[1]
(1)
(2)
Yk1 ,k2 ,±1 (ω) = Hk1,k2 (ω) ± Hk1,k2 (ω) X(ω)

(Uk1 (ω1 )Uk2 (ω2 ) + Vk1 (ω1 )Vk2 (ω2 )) X(ω)=(Uk1 (ω1 )Vk2 (ω2 ) + Vk1 (ω1 )Uk2 (ω2 )) X(ω)
In the case of the DTCWT, directionality can be achieved atany scale as described in the conﬁguration in Fig. 2. However,this procedure does not work in the CSMFB case at higher
decomposition levels J > 1. For example, in the second
scale J = 2, from the deﬁnition of the subband ﬁlters in the
CSMFB, the following equation can be veriﬁed.
(1)

(2)

Fk (z) = z N −1 Hk (z −1 ), Fk (z) = z N −1 Hk (z −1 ).
(4)where P (z) is the z-transform of the prototype ﬁlter and ck =(k+ 1 ) N −1
ejϕk W2M 2 2 . The frequency spectra Uk (ω) and Vk (ω),
which are modulations of the prototype ﬁlter, are lying only inthe positive and the negative domain, respectively. According
to the above setting, if the ﬁrst system satisﬁes the perfectreconstruction (PR) property, the whole system also satisﬁes
the PR property [?].
DTCWTs also consists of two maximally decimated FBs
as shown in Fig. 1 [?]. However, the pairs of subband ﬁlters
should satisfy the fractional sample delay condition. Since itis difﬁcult to approximate the fractional delay with FIR ﬁlters,
the resulting ﬁlter performance and efﬁciency in practicalapplication are worse than CSMFBs.

(2)

(2)

[2]

[2]

=Yk1 ,k2 ,±1 (ω) + Δ = Yk1 ,k2 ,±1 (ω),(5)

where Δ includes the frequency components except the(k1 , k2 ) subband. Fig. 3(b) shows (the half of) the second level
directional subband coefﬁcients of the test image Zoneplate(Fig. 3(a)). The CSMFB cannot distinguish ±45◦ directions
at the second level. Therefore, CSMFBs cannot be extended
to multiscale decomposition straightforwardly.
B. Subband Filters of MDMLT
To design a MDT based on the CSMFB, we use modulated
lapped transforms the prototype ﬁlter p(n) and the phase termθ̂(k, n) of which are deﬁned as follows:(1)ĥk (n) := p(n)ĉk,n = p(n) cos(θ̂(k, n))(2)ĥk (n) := p(n)ŝk,n = p(n) sin(θ̂(k, n))




1
π
π
sin
n+
p(n) = −
M
2M
2



π
1
M +1
θ̂(k, n) =
k+
n+
M
2
2

A. Lack of directionality
Before introducing the proposed transform, this section
mentions the shortage of the CSMFB, i.e., lack of multiscaledirectionality. To obtain the directional subband signal with
CSMFBs at the ﬁrst scale J = 1, the ﬁrst and the second FB
are separably performed along vertical and horizontal directions. Then, output signals from each FBs are added/subtracted(1)

Hk1 ,k2 (2ω)H0,0 (ω)X(ω) ± Hk1 ,k2 (2ω)H0,0 (ω)X(ω)

III. M ULTISCALE M ODULATED L APPED T RANSFORMS

978-1-5386-1542-3@2017 APSIPA

Fig. 2. The conﬁguration to obtain directional subband signals in the 2D
frequency plane (the number of channel: M = 2).

Fig. 1. M -channel CSMFBs
(1)





(6)

For the number of channel M , the ﬁlter length of the MLT is
2M . Since it can be veriﬁed that ck,n = (−1)σ(k) ĉk,N −1−n(σ(k) = 0 or 1, ck,n is in (1)), the MLT can be considered

1831

APSIPA ASC 2017.

Proceedings of APSIPA Annual Summit and Conference 201712 - 15 December 2017, Malaysia

as a subclass of the CSMFB. Based on the prototype ﬁlter,
we derive the following relationship between the ﬁrst and thesecond lowpass ﬁlters:
Theorem 1:




M
M
(2)
= p n−
ŝ0,n− M
n−
ĥ0
2
2
2
(1)

= − p(n)ĉ0,n = −ĥ0 (n)(7)

Proof: From the following equations:



1
π 1
n+
= −p(n)ŝ0,n− M = sin
2
M2
2





M
π
M +1
−M
= ĉ0,n ,p n−= − sinn+22M2
we can derive the theorem.
This indicates that the 2D frequency response of the 2D
separable ﬁlter of the ﬁrst and the second FBs are identical:
(1)
(1)
(1)
Ĥ0,0 (ω) = Ĥ0 (ω1 )Ĥ0 (ω2 )
(2)
(2)
(2)
= Ĥ0 (ω1 )Ĥ0 (ω2 ) = Ĥ0,0 (ω).

(b)

(c)

expression from the viewpoint of a matrix multiplication as:
J −1

max

H = WT[Jmax ]T[J] T[1] ,
1 IN 2W= √2 IN 2

J=2

 



[1]
C[1] ⊗ C[1]
IN 2
T1
[1]
, T =
[1] =
−IN 2
S[1] ⊗ S[1]
T2

[J]

[J]

[J]

T[J] = diag(T1 , T2 ), T1 = diag(C[J] ⊗ C[J] , IN 2 −( NJ )2 )M[J]

T2 = diag(S[J] ⊗ S[J] , IN 2 −( NJ )2 ),

(11)

M

(8)

With this in mind, we deﬁne the (analysis) subband ﬁlters of
MDMLT as:(1)ĥk (n) = p(n)ĉk,n⎧M⎪⎨p n − 2 ŝ0,n− M2(2)ĥk (n) = p(n)ŝ0,n⎪⎩p(n)ŝk,n(a)Fig. 3. (a) Zone plate, (b) (half of) second level transformed coeff. of the
CSMFB, (c) (half of) second level transformed coeff. of the MDMLT.

(k = 0, J < Jmax )(k = 0, J = Jmax ) (9)(k = 0)According to the above setting, we can realize the MDMLT,
which can be veriﬁed by analyzing the 2D frequency responses
at the higher scale. For example, at the second scale J = 2,
we can derive as

[J]

J

[J]

J

where the matrices C , S ∈ RN/M ×N/M denote 2D(1)(2)
separable transform of {ĥk (n)}, {ĥk (n)} with downsampling at J-th scale. Then, the inverse transform F can be givenby the Moore-Penrose pseudo-inverse matrix as:
−1 
F = H H
H .
(12)
Since the MLT is paraunitary FB, i.e., (C[J] ⊗ C[J] ) (C[J] ⊗C[J] ) = I, (S[J] ⊗ S[J] ) (S[J] ⊗ S[J] ) = I,⎛ J −1
⎞−1
J −1
max
max


[J]
[J]
T2
T2 ⎠ HF = ⎝I +
J=1−1= V (I + D)

J=0




V H ,

(13)

where the unitary matrix V and the diagonal matrixDthe eigenvalue decomposition as:
 are obtained byJmax −1 [J] 
(1)
(2)
(1)
Jmax −1 [J]= (Ĥk1 ,k2 (2ω) ± Ĥk1 ,k2 (2ω))Ĥ0,0 (ω)X(ω)
T2T2 = VDV .
J=0
J=0

(1)
1) Two-level Decomposition Case: As shown in (13), the(Uk1 (2ω1 )Uk2 (2ω2 ) + Vk1 (2ω1 )Vk2 (2ω2 )) Ĥ0,0 (ω)X(ω)=
MDMLT requires the eigenvalue
decomposition of the huge(1)
(Uk1 (2ω1 )Vk2 (2ω2 ) + Vk1 (2ω1 )Uk2 (2ω2 )) Ĥ0,0 (ω)X(ω)
Jmax −1 [J]  Jmax −1 [J] matrixto ﬁnd the maT2
T2
[2]
J=0
J=0
(10)
= Yk1 ,k2 ,±1 (ω).
trices V and D, which consumes signiﬁcant computational
It indicates that the proposed MDMLT, the redundancy ratio is cost. For this problem, if we restrict that the number of the2, has multiscale directionality. Fig. 3(c) shows (the half of) decomposition level for MDMLT is 2, more computationallythe second level directional subband coefﬁcients of the test efﬁcient MDMLT can be designed. The 2-level MDMLT isimage Zoneplate (Fig. 3(a)). In this case, the MDMLT can expressed as:
distinguish ±45◦ directions at the second level.
H = WT[2] T[1]−1
C. Inverse Transform of MDMLTF = I + (S[1] ⊗ S[1] ) (S[1] ⊗ S[1] )
H
(14)
(1)

(1)

(2)

(2)

Ĥk1 ,k2 (2ω)Ĥ0,0 (ω)X(ω) ± Ĥk1 ,k2 (2ω)Ĥ0,0 (ω)X(ω)

Now, we consider the inverse transform of the proposed(1)
MDMLT. While the ﬁrst FB {ĥk (n)} is invertible by setting(1)(2)
synthesis ﬁlter {fˆk (n)} as in (4), the second FB {ĥk (n)}M
is not since its lowpass ﬁlter is shifted by 2 samples.
To derive the inverse transform, we ﬁrst describe the forward transform with Jmax level decomposition by using the978-1-5386-1542-3@2017 APSIPA

From the eigenvalue decomposition of S[1]U[1] D[1] (V[1] ) , we can derive−1
F = I + (S[1] ⊗ S[1] ) (S[1] ⊗ S[1] )H

== (V[1] ⊗ V[1] )(I + D[1] )−1 (V[1] ⊗ V[1] ) H .

(15)

1832

APSIPA ASC 2017

Proceedings of APSIPA Annual Summit and Conference 2017(b) Barbara(a) Lena

12 - 15 December 2017, Malaysia(c) Pepper(a) Original

Fig. 4. Original images(b) Noisy

TABLE IT HE NUMERICAL RESULTS OF NLA
Zoneplate (256 × 256)
Num. of coeff.
2000
3000
8 × 22 DTCWT
12.57
14.87
8 × 16 CSMFB
14.91
18.93
64 × 128 CSMFB
11.11
11.93
8 × 16 MDMLT (Prop.)
18.06
20.51
Barbara (512 × 512)
Num. of coeff.
5000
10000
8 × 22 DTCWT
23.58
26.96
8 × 16 CSMFB
24.18
28.53
64 × 128 CSMFB
25.83
27.93
8 × 16 MDMLT (Prop.)
27.85
28.92

4000
16.64
21.61
12.71
22.99
15000
28.84
30.87
29.38
30.87

5000
17.83
24.22
13.47
25.02

(c) DTCWT

20000
30.21
32.57
30.55
32.36

(e) MDMLTTABLE II
T HE NUMERICAL RESULTS OF IMAGE DENOISING (PSNR [ D B])

This indicates that, except the element-wise multiplication(I + D[1] )−1 , each procedure is implemented as separable
transform, and thus the computational complexity can besigniﬁcantly reduced.
IV. E XPERIMENTAL R ESULTS

Test image
Noisy image
8 × 22 DTCWT [?]64 × 128 CSMFB [?]8 × 16 MDMLT

Zoneplate
22.12
23.34
23.47
27.83

Lena
22.12
28.88
28.55
29.79

Barbara
22.12
27.10
26.97
28.43

Pepper
22.12
28.96
28.25
29.70

than the DTCWT due to ﬂexible multiscale directional image
representation.

A. Nonlinear Approximation

V. C ONCLUSIONSThis section evaluates the efﬁciency of the MDMLT by
NLA [?]. NLA is often used for the evaluation of sparse
representation by reconstructing images from the K-largest
K−largest=(magnitude) transformed coefﬁcients y = Fx −−−−−−−→ ycoeﬃcientsFx. Here we examined the 2-level DTCWTs with the FB
size of 8 × 22, the 1-level CSMFB with the size of 8 × 16,
and the 2-level MDMLT with the size of 8 × 16. We also
compare the 1-level 64 × 128 CSMFB of which 2D lowpass
subband is the same as that of the 2-level 8 × 16 MDMLT.
All the redundancy ratios are 2. Zoneplate (256×256 pixels)
and Barbara (512×512 pixels), which have rich directional
textures, are used as test images (512×512 pixels).
Table II shows the reconstruction error (PSNR [dB]) fromthe number of coefﬁcients retained. It shows the reconstructed
images by the proposed method have consistently higher
PSNRs compared to those by the conventional ones.
B. Image Denoising
In this section, the performances of the 2-level 8 × 22
DTCWT, the 1-level 64 × 128 CSMFB, the 2-level 8 × 16
MDMLT in image denoising is veriﬁed. The test images
Zoneplate, Lena, Barbara, and Pepper, are added the Gaussian
noise with variance σ 2 (σ = 20). Directional transforms ﬁrst
decompose the noisy images, then hard-thresholding with 3σ
is performed. From Fig. 5, since the CSMFB requires a large
support of atoms to decompose the lowpass frequency region
ﬁnely, the resulting image of the CSMFB is signiﬁcantly
degraded with severe ringing artifacts. On the other hand, thedenoising result of the MDMLT provides better visual quality

978-1-5386-1542-3@2017 APSIPA(d) CSMFB

Fig. 5. Original, noisy, and denoise images of Barbara.

In this paper, we proposed the MDMLT as a multiscaledirectional transform. To construct the MDMLT, we ﬁrst
derived the relationship between the lowpass ﬁlter of the ﬁrst
FB and the shifted lowpass ﬁlter of the second FB. Then,
we designed the MDMLT based on the MLT and the shifted
MLT. Finally, in the NLA and image denoising simulations,
it was clariﬁed that the proposed MDMLT provides better
numerical and visual reconstruction quality compared with theconventional CSMFBs and the DTCWTs.
R EFERENCES[1] R. G. Baraniuk I. W. Selesnick and N. G. Kingsbury, “The dual-treecomplex wavelet transform,” IEEE Signal Process. Mag., pp. 123–151,
Nov. 2005.
[2] I. Bayram and I.W. Selesnick, “On the dual-tree complex wavelet packet
and M -band transforms,” IEEE Trans. Signal Process., vol. 56, no. 6,
pp. 2298–2310, June 2008.
[3] M. N. Do and M. Vetterli, “The contourlet transform : An efﬁcient
directional multiresolution image representation,” IEEE Trans. Image
Process., vol. 14, no. 12, pp. 2091–2106, 2005.
[4] E. J. Candès and D. L. Donoho, “New tight frames of curvelets
and optimal representations of objects with piecewise C 2 singularities,”
Comm. Pure Appl. Math, vol. 56, no. 2, pp. 219–266, Feb. 2004.
[5] S. Kyochi, T. Uto, and M. Ikehara, “Dual-tree complex wavelet transformarising from cosine-sine modulated ﬁlter banks,” in Proc. IEEE Int. Symp.
Circuits Syst., May 2009, pp. 2189–2192.
[6] S. Kyochi and M. Ikehara, “A class of near shift-invariant and orientationselective transform based on delay-less oversampled evenstacked cosinemodulated ﬁlter banks,” IEICE Trans. Fundam. Electron, Commun.
Comput,, vol. E93-A, pp. 724–733, Apr. 2010.
[7] L. Liang and H. Liu, “Dual-tree cosine-modulated ﬁlter bank with linearphase individual ﬁlters: An alternative shift-invariant and directionalselective transform,” IEEE Trans. Image Process., vol. 22, no. 12, pp.
5168–5180, Dec. 2013.

1833

APSIPA ASC 2017

Proceedings of APSIPA Annual Summit and Conference 201712 - 15 December 2017, Malaysia[8] Stphane Mallat, A Wavelet Tour of Signal Processing, Third Edition: The
Sparse Way, Academic Press, 3rd edition, 2008.

978-1-5386-1542-3@2017 APSIPA1834

APSIPA ASC 2017