# Directional Analytic Discrete Cosine Frames

Seisuke Kyochi , Taizo Suzuki, and Yuichi Tanaka S. Kyochi is with the
Department of Computer Science, Kogakuin University, Tokyo 163-8677, Japan
(e-mail: kyochi@cc.kogakuin.ac.jp).T. Suzuki is with the Faculty of
Engineering, Information and Systems, University of Tsukuba, Tsukuba 305-8573
Japan (e-mail:taizo@cs.tsukuba.ac.jp).Y. Tanaka is with the Department of
Computer and Information Sciences, Tokyo University of Agriculture and
Technology, Koganei, Tokyo, 184-8588 Japan (e-mail: ytnk@cc.tuat.ac.jp).

###### Abstract

Block frames called directional analytic discrete cosine frames (DADCFs) are
proposed for sparse image representation. In contrast to conventional
overlapped frames, the proposed DADCFs require a reduced amount of 1)
computational complexity, 2) memory usage, and 3) global memory access. These
characteristics are strongly required for current high-resolution image
processing. Specifically, we propose two DADCFs based on discrete cosine
transform (DCT) and discrete sine transform (DST). The first DADCF is
constructed from parallel separable transforms of DCT and DST, where the DST
is permuted by row. The second DADCF is also designed based on DCT and DST,
while the DST is customized to have no DC leakage property which is a
desirable property for image processing. Both DADCFs have rich directional
selectivity with slightly different characteristics each other and they can be
implemented as non-overlapping block-based transforms, especially they form
Parseval block frames with low transform redundancy. We perform experiments to
evaluate our DADCFs and compare them with conventional directional block
transforms in image recovery.

###### Index Terms:

Block frame, discrete cosine transform, directional selectivity, sparse image
representation

##  I Introduction

Sparse representation by frames has been an essential technique for image
analysis and processing [1, 2]. Various kinds of signal recovery tasks, e.g.,
denoising, deblurring, and restoration from compressive samples, can also be
realized by incorporating sparse image representation in convex and non-convex
optimization algorithms[3, 4].

Significant efforts have been made to construct efficient frames for sparse
image representation. Of particular focus has been directional frames, such as
curvelet [5], contourlet [6], directional filter banks (FBs) [7, 8, 9, 10,
11], and dual-tree complex wavelet transforms (DTCWTs) [12, 13, 14, 15, 16]
for 2D signals. Directional atoms111In this paper, atom is referred to as an
element $\mathbf{d}_{n}$ $(n=0,\ldots,N-1)$ of a frame
$\mathbf{D}=\begin{bmatrix}\mathbf{d}_{0}&\cdots&\mathbf{d}_{N-1}\end{bmatrix}\in\mathbb{R}^{M\times
N}$ [17]. are crucial for sparse image representation since images usually
contain edges and textures lying along oblique directions. Extended versions
of these frames for higher-dimensional signals, such as videos, have also been
proposed [18, 19, 20, 21]. In addition to directional frames, more general
systems such as dictionary [22, 23] and graph WTs/FBs [24] that capture highly
complex structures and non-local similarity have been proposed. Those
transforms can provide sparse image representation for various fine
components.

Although these existing frames and dictionaries have been successfully applied
to image processing, problems and limitations have been recognized in
practical situations. First, computational complexity for calculating sparse
coefficients is typically high due to the complicated algorithms involved,
such as 2D filtering [6], sparse coding with various iterative schemes [22,
23], and eigenvalue decomposition of a large-scale graph Laplacian [24].
Second, they typically require high transform redundancy which leads to a
large amount of memory usage to store the coefficients. Third, since the
supports of their atoms in those frames are overlapped with each other, they
require global memory access, which disrupts parallel computation. Although
recent digital devices have been increasing their computational power, the
resolution of captured images has also been increasing and sometimes multiple
images will also be taken at once for producing visually pleasant images like
those having low-noise and/or high-dynamic range. Hence, the computational
cost for image processing has to be kept as small as possible for avoiding
installing extra hard/software modules in such devices.

Block-based bases and frames, whose supports are identical or disjoint, are
thus highly desired due to their computational efficiency. They are still a
key for many image processing applications like video coding. In addition,
patch-based techniques based on block-based transforms, such as BM3D and
redundant (type-II) discrete cosine transform (DCT) [25, 26, 27], show their
effectiveness in image recovery. Nevertheless, directional block frames have
received less attention compared with overlapped frames and, unfortunately,
they have been believed that they cannot provide rich directional
selectivity222In this paper, “directional selectivity” is measured by the
number of distinguishable directional subbands for an $M$-channel 2D
transform.. However, we can realize such block-based transforms by carefully
choosing their building blocks.

In this paper, we focus on directional block frames and propose directional
analytic discrete cosine frames (DADCFs) based on DCT [28] and (type-II)
discrete sine transform (DST) [29]. They have the following advantages against
alternative directional block frames:

  * •

Directional selectivity of the DADCFs is much richer than that of existing
directional block frames.

  * •

DADCFs can be designed by appending the DST (or a DST-like transform) and
simple extra operations to the DCT, and thus are compatible with the DCT.

We introduce two types of DADCF, both forming Parseval frames. The first DADCF
contains the DCT and a row-wise permuted version of the DST, and the second
DADCF contains the DCT and a DST without DC leakage. The second one is called
regularity-constrained DADCF (RDADCF). In order to realize RDADCF, we propose
the DST without DC leakage, regularity-constrained DST (RDST), for the first
time. The DADCF and the RDADCF have different advantages. The DADCF provides
richer directional selectivity while the DC energy will be distributed over
several subbands. As it will be described later, the DC leakage can be avoided
by integrating the DADCF with Laplacian pyramid at the expense of redundancy.
In contrast, the RDADCF can structurally avoid the DC leakage problem by
incorporating the proposed RDST as its building block instead of the row-wise
permuted DST. We numerically compare two DADCFs with some existing approaches
in compressive sensing reconstruction.

The rest of this paper is organized as follows. Section II summarizes related
works. Section III reviews the conventional directional block bases and the
analyticity for images. Section IV explains the definition and a customization
for preventing DC leakage of the DADCF. Section V introduces the RDADCF.
Section VI evaluates the proposed DADCFs in compressive sensing
reconstruction. Section VII concludes with a brief summary.

###  I-A Notations

Bold-faced lower-case letters and upper-case letters are vectors and matrices,
respectively. The subscripts $h$ and $v$ are used to indicate variables
corresponding to horizontal and vertical directions, respectively. The other
mathematical notations are summarized in Table I.

##  II Related Works

Directional block bases and frames can be classified into two categories. One
is the fixed class, i.e., transforms equipped with directionally oriented
bases. This class of transforms includes discrete Fourier transform (DFT)
[30], discrete Hartley transform (DHT) [31], and real-valued conjugate-
symmetric Hadamard transform [32]. The other is the adaptive class, i.e., the
application of a non-directional block transform (e.g., the DCT) along
suitable oblique directions provided by preprocessing (e.g., edge analysis)
for each block [33, 34]. Applications of the latter class are relatively
limited because transform directions have to be determined from an input
signal in advance. For example, in signal recovery, degraded observations make
it difficult to find suitable directions. Our directional block frames
correspond to the fixed class.

The main problem with DFT and its variants is that they contain duplicated
atoms along the same direction in their basis and hence cannot provide rich
directional selectivity (i.e., the number of directional orientations in a
basis or a frame). This degrades the efficiency of signal analysis and
processing. In order to achieve richer directional selectivity, in this paper,
we extend the DCT to the DADCF. Definitely, the DCT is one of the most
effective block transforms for image processing tasks and is already
integrated into many digital devices. For example, video coding standards,
e.g. HEVC [35] and VVC [36], employ the various sizes of the (integer) DCT.
However, since it does not contain obliquely oriented atoms in its basis, it
cannot achieve rich directional selectivity.

In this paper, we reveal that by appending some extra modules, i.e., DST and
scaling/addition (and subtraction)/permutation (SAP) operations, to the DCT,
the resulting transform provides directionally oriented atoms and thus leads
to rich directional selectivity. Furthermore, since the DST can be designed by
the (row-wise) flipped and sign-altered version of the DCT, the implementation
cost for the proposed transform can be kept low, i.e., the total procedure can
be fully implemented by using the DCT and a few SAP operations.

A preliminary version of this work was presented in [37], which provides a
basic structure of the DADCF. In this paper, we newly introduce theory and
design algorithm of the RDADCF, and comprehensive experiments.

TABLE I: Notations

---
Notation | Terminology
$\mathbb{R}$ | Real numbers
$j$ | $\sqrt{-1}$
$\Omega_{N}$, $\Omega_{N_{1},N_{2}}$ |  $\\{0,\ldots,N\\}$, $\\{N_{1},\ldots,N_{2}\\}$
$\mathbb{R}^{N}$ |  $N$-dimensional real-valued vector space
$\mathbb{R}^{N_{v}\times N_{h}}$ |  $N_{v}\times N_{h}$ real-valued matrices
$\mathbf{I}$, $\mathbf{O}$ | Identity matrix, zero matrix
$\mathbf{A}_{N}$ |  $N\times N$ square matrix
$\mathbf{A}^{\top}$ | Transpose of $\mathbf{A}$
$x_{n}$ and $[\mathbf{x}]_{n}$ |  $n$-th element of a vector $\mathbf{x}$
$X_{m,n}$ and $[\mathbf{X}]_{m,n}$ |  $(m,n)$-th element of a matrix $\mathbf{X}$
$\mathbf{X}^{(m,n)}\in\mathbb{R}^{M\times M}$ |  $(m,n)$-th $M\times M$ subblock of $\mathbf{X}\in\mathbb{R}^{ML_{v}\times ML_{h}}$
$\mathrm{vec}(\mathbf{X})\in\mathbb{R}^{N_{v}N_{h}}$ |  | Vectorization of $\mathbf{X}\in\mathbb{R}^{N_{v}\times N_{h}}$
---
$x_{N_{v}n+m}=X_{m,n}$
$\mathrm{bvec}(\mathbf{X})$ |  | Block-wise vectorization of $\mathbf{X}\in\mathbb{R}^{ML_{v}\times ML_{h}}$
---
$[\mathrm{vec}(\mathbf{X}^{(0,0)})^{\top}\ \ldots\
\mathrm{vec}(\mathbf{X}^{(L_{v}-1,L_{h}-1)})^{\top}]^{\top}$
$\otimes$ | Kronecker product
$\mathcal{N}(\mathbf{A})$ | Null space of a matrix: $\mathbf{A}$
| $\mathrm{diag}(a_{0},\ldots,a_{N-1})$
---
$\mathrm{diag}(\mathbf{A}^{(\mathrm{0})},\ldots,\mathbf{A}^{({N-1})})$
Diagonal/block-diagonal matrices.
$H(z)$ |  $z$-transform: $\sum_{n}h(n)z^{-n}$
$H(\omega)$, $H(e^{j\omega})$, $\mathcal{F}[h]$ | Discrete-time Fourier transform: $\sum_{n}h(n)e^{-j\omega n}$
$H_{k}(z)$ ($k=0,\ldots,M-1$) |  $M$-channel subband filters
$H({\bm{\omega}})$ | $H(e^{j\omega_{v}},e^{j\omega_{h}})$
$H_{k_{v},k_{h}}({\bm{\omega}})$ |  | $H_{k_{v}}(\omega_{v})H_{k_{h}}(\omega_{h})$
---
$\|\mathbf{x}\|_{p}$ ($\mathbf{x}\in\mathbb{R}^{N}$, $p\in[1,\infty)$) |  |  $\ell_{p}$-norm, $\|\mathbf{x}\|_{p}=\left(\sum^{N-1}_{n=0}|x_{n}|^{p}\right)^{\tfrac{1}{p}}$
---
|  $\|\mathbf{x}\|_{1,2}$ ($\mathbf{x}\in\mathbb{R}^{KN}$)
---
|  $\ell_{1,2}$-mixed norm,
---
$\|\mathbf{x}\|_{1,2}=\sum^{K-1}_{k=0}\left(\sum^{N-1}_{n=0}|{x}_{Nk+n}|^{2}\right)^{\tfrac{1}{2}}$

##  III Preliminaries

###  III-A Conventional Block Bases

The DCT [28] is one of the most popular time-frequency transforms. Its
transform matrix $\mathbf{F}^{(\mathrm{C})}\in\mathbb{R}^{M\times M}$
($M=2^{m}$, $m\geq 1$)333For simplicity, we restrict the sizes of all the
block transforms to $M=2^{m}$ throughout this paper, but it is easily extended
to the general $M$. is defined as

$\displaystyle[\mathbf{F}^{(\mathrm{C})}]_{k,n}=$ $\displaystyle\
\alpha_{k}\sqrt{\frac{2}{M}}\cos(\theta_{k,n})$ (1)

where $k$ and $n$ are subband and time indices ($k,\ n\in\Omega_{M-1}$),
$\theta_{k,n}=\frac{\pi}{M}k\left(n+\frac{1}{2}\right)$,
$\alpha_{k}=\tfrac{1}{\sqrt{2}}$ for $k=0$ and $\alpha_{k}=1$ for otherwise.
For $\mathbf{x}=\mathrm{vec}(\mathbf{X})$ ($\mathbf{X}\in\mathbb{R}^{M\times
M}$), the 2D DCT is given by
$\mathbf{F}^{(\mathrm{C})}\otimes\mathbf{F}^{(\mathrm{C})}\
\in\mathbb{R}^{M^{2}\times M^{2}}$. Since the DCT is the approximation of the
Karhunen-Loève transform for a first-order Markov process with a correlation
coeffcient $\rho$ when $\rho\rightarrow 1$, the 2D DCT coefficients of natural
images tend to be sparse (i.e., its $\ell_{1}$ norm
$\|(\mathbf{F}^{(\mathrm{C})}\otimes\mathbf{F}^{(\mathrm{C})})\mathbf{x}\|_{1}$
is small). Thus, the DCT is widely applied to many applications, especially
for source coding. However, it is a separable transform and hence it lacks
directional selectivity. Formally, its 2D atom
$\mathbf{B}^{(k_{v},k_{h})}\in\mathbb{R}^{M\times M}$ in the DCT basis forms

$\displaystyle{B}^{(k_{v},k_{h})}_{n_{v},n_{h}}=\alpha_{k_{v}}\alpha_{k_{h}}\frac{2}{M}\cos\left(\theta_{k_{v},n_{v}}\right)\cos\left(\theta_{k_{h},n_{h}}\right),$
(2)

where $k_{d}$ and $n_{d}$ ($d\in\\{h,v\\}$) denote subband and spatial
indices, respectively ($k_{d},n_{d}\in\Omega_{M-1}$). Fig. 1(a) shows an
example of the 2D DCT atoms444In Figs. 1, 4, and 8, each atom is enlarged for
visualization.. Clearly, they “mix” two diagonal components along $45^{\circ}$
and $-45^{\circ}$ which reduce directional selectivity.

The 2D DFT can be regarded as a block transform with directional selectivity
because it is a complex-valued transform. Its 2D atoms
$\mathbf{B}^{(k_{v},k_{h})}\in\mathbb{R}^{M\times M}$ are represented as

$\displaystyle{B}^{(k_{v},k_{h})}_{n_{v},n_{h}}=\frac{1}{M}e^{j\left(\varphi_{k_{v},n_{v}}+\varphi_{k_{h},n_{h}}\right)},$
(3)

where $\varphi_{k,n}=\frac{2\pi}{M}kn$. As shown in Fig. 1(b), the DFT bases
can decompose diagonal components into different subbands.

There are some real-valued variants of the DFT [31, 32] that provide
directionally oriented atoms. For example, the DHT [31] can form a
directionally oriented basis by modifying some of the original DFT atoms
${B}^{(k_{v},k_{h})}_{n_{v},n_{h}}=\frac{1}{M}\mathrm{cas}(\varphi_{k_{v},n_{v}})\mathrm{cas}(\varphi_{k_{h},n_{h}})$,
$(\mathrm{cas}(\varphi_{k,n})=\cos(\varphi_{k,n})+\sin(\varphi_{k,n}))$ to
${B}^{(k_{v},k_{h},\pm 1)}_{n_{v},n_{h}}$ as

$\displaystyle{B}^{(k_{v},k_{h},\pm 1)}_{n_{v},n_{h}}=$
$\displaystyle\frac{1}{2}\left({B}^{(k_{v},k_{h})}_{n_{v},n_{h}}\pm{B}^{(M-k_{v},M-k_{h})}_{n_{v},n_{h}}\right)$
$\displaystyle=$ $\displaystyle\
\frac{1}{M}\begin{cases}\cos(\varphi_{k_{v},n_{v}}-\varphi_{k_{h},n_{h}})\\\
\sin(\varphi_{k_{v},n_{v}}+\varphi_{k_{h},n_{h}})\end{cases},$ (4)

where we assume $M\geq 4$ and $k_{v},\ k_{h}\neq 0,M/2$.

One problem shared by these conventional directional block transforms is that
they contain multiple atoms along the same direction in their basis. For the
$M\times M$ DFT and DHT, the number of distinguishable directions is
$2\left(\tfrac{M-2}{2}\right)^{2}$ compared to the number of atoms $M^{2}$.
They cannot provide rich directional selectivity, as shown in Fig. 1(b).

(a) DCT

(b) DFT (real part)

Figure 1: Atoms ${B}_{n_{v},n_{h}}^{(k_{v},k_{h})}$ in basis ($M=4$).

###  III-B Analyticity for Directional Selectivity

As explained in Section III-A, the 2D DCT cannot provide a directional image
representation. We explain this phenomenon in the 2D frequency domain. Let
$H_{k}(\omega)$ be a frequency spectrum of the $k$-th row of the DCT, i.e.,
$H_{k}(\omega)=\mathcal{F}[[\mathbf{F}_{C}]_{k,\cdot}]$. Since $H_{k}(\omega)$
is the frequency response of a real-valued filter, its spectrum symmetrically
distributes in both positive and negative $\omega$ (Fig. 2(a)). Thus, the
frequency spectrum of the 2D separable DCT $H_{k_{v},k_{h}}({\bm{\omega}})$
always has nonzero frequency responses in four quadrants, as in Fig. 2(c), and
it mixes $\pm 45^{\circ}$ frequency spectra for example.

In contrast, any spectrum of the DFT
$U_{k}(\omega)=\frac{1}{\sqrt{M}}\mathcal{F}[e^{-j\varphi_{k,\cdot}}]$
(complex-valued filter) has a frequency response in only positive (or
negative) $\omega$, as in Fig. 2(b). This property is referred to as
analyticity [12], i.e., $|U_{k}(\omega)|\approx 0$ for $\omega<0$ (or
$\omega>0$). Thus, frequency spectra of the 2D separable DFT
$U_{k_{v},k_{h}}({\bm{\omega}})$ are localized in one quadrant (Fig. 2(d)),
which indicates the directional subband.

Conventional separable directional WTs/FBs utilize analyticity. For example,
DTCWTs consist of two $M$-channel filter banks
$\\{H_{k}(\omega)\\}_{k=0}^{M-1}$ and $\\{G_{k}(\omega)\\}_{k=0}^{M-1}$, where
those complex combination satisfies analyticity as follows:

$\displaystyle H_{k}({\omega})=$ $\displaystyle\
\frac{1}{2}\left(U_{k}({\omega})+\overline{U_{k}({\omega})}\right),$
$\displaystyle G_{k}({\omega})=$ $\displaystyle\
\frac{1}{2j}\left(U_{k}({\omega})-\overline{U_{k}({\omega})}\right),$
$\displaystyle U_{k}({\omega})=$ $\displaystyle\
H_{k}({\omega})+jG_{k}({\omega}),\ |U_{k}({\omega})|\approx 0\ (\omega<0)$ (5)

Here, we assume that the frequency spectrum $U_{k}({\omega})$ distributes in
the positive frequency domain (Fig. 2(b)). Then, by using the 2D frequency
spectra of the complex-valued filters
$U_{k_{v},k_{h}}({\bm{\omega}})=U_{k_{v}}(\omega_{v})U_{k_{h}}(\omega_{h})$
and
$U_{k_{v},\overline{k_{h}}}({\bm{\omega}}):=U_{k_{v}}({\omega_{v}})\overline{U_{k_{h}}({\omega}_{h})}$,
the 2D directional frequency spectrum of the real-valued filter can be
designed as follows:

$\displaystyle\
\frac{1}{2}\left(U_{k_{v},k_{h}}({\bm{\omega}})+\overline{U_{k_{v},k_{h}}({\bm{\omega}})}\right)$
$\displaystyle=\
H_{k_{v}}({\omega}_{v})H_{k_{h}}({\omega}_{h})-G_{k_{v}}({\omega}_{v})G_{k_{h}}({\omega}_{h}),$
$\displaystyle\
\frac{1}{2}\left(U_{k_{v},\overline{k_{h}}}({\bm{\omega}})+\overline{U_{k_{v},\overline{k_{h}}}({\bm{\omega}})}\right)$
$\displaystyle=\
H_{k_{v}}({\omega}_{v})H_{k_{h}}({\omega}_{h})+G_{k_{v}}({\omega}_{v})G_{k_{h}}({\omega}_{h}).$
(6)

Considering (III-B), a directional frequency decomposition can be realized by
two 2D separable FBs followed by addition/subtraction, as in Fig. 3.
$M$-channel DTCWTs can distinguish $2M^{2}$ directional subbands.

(a) Real-valued filter

(b) Complex-valued filter

(c) $H_{k_{v},k_{h}}({\bm{\omega}})$

(d) $U_{k_{v},k_{h}}({\bm{\omega}})$

(e) $U_{k_{v},k_{h}}+\overline{U_{k_{v},k_{h}}}$

Figure 2: Example of frequency spectra (analytic and non-analytic filters).

Figure 3: Configurations for 2D DTCWTs. For $M=4$, 32 directional subbands can
be distinguished.

##  IV Directional Analytic Discrete Cosine Frames

This section introduces the DADCF. The definition of the DADCF is given in
Section IV-A. Directional selectivity of the DADCF is then discussed by
analyzing its atoms in Section IV-B. As it will be shown in Section IV-D, the
DADCF suffers from the DC leakage problem. One solution is given by
constructing the DADCF pyramid (the DADCF with Laplacian pyramid) in Section
IV-D.

###  IV-A Definition of Directional Analytic Discrete Cosine Frame

This section introduces DADCFs for 2D signals by extending the conventional
DCT.

###### Definition 1

The analysis operator of the DADCF $\mathbf{F}^{(\mathrm{D})}\ \in\
\mathbb{R}^{2M^{2}\times M^{2}}$ is defined as

$\displaystyle\mathbf{F}^{(\mathrm{D})}:=$ $\displaystyle\
\mathbf{P}^{(\mathrm{I})\top}\mathbf{W}^{(\mathrm{I})}\mathbf{P}^{(\mathrm{I})}\begin{bmatrix}\mathbf{F}^{(\mathrm{C})}\otimes\mathbf{F}^{(\mathrm{C})}\\\
\mathbf{F}^{(\mathrm{S})}\otimes\mathbf{F}^{(\mathrm{S})}\end{bmatrix},$
$\displaystyle\mathbf{W}^{(\mathrm{I})}=$ $\displaystyle\
\mathrm{diag}\left(\frac{1}{\sqrt{2}}\mathbf{I}_{2M-1},\frac{1}{2}\begin{bmatrix}\mathbf{I}_{(M-1)^{2}}&-\mathbf{I}_{(M-1)^{2}}\\\
\mathbf{I}_{(M-1)^{2}}&\mathbf{I}_{(M-1)^{2}}\end{bmatrix}\right),$ (7)

where $\mathbf{F}^{(\mathrm{C})}$ is defined in (1) and
$\mathbf{P}^{(\mathrm{I})}\in\mathbb{R}^{M^{2}\times M^{2}}$ is a permutation
matrix that places the $2M-1$ DCT and $2M-1$ DST coefficients associated with
the subband indices $k_{v}=0$ or $k_{h}=0$ to the first part, and the other
$2(M-1)^{2}$ coefficients associated with the subband indices $k_{v}\neq 0$
and $k_{h}\neq 0$ to the last (see Fig. 4(a)). $\mathbf{F}^{(\mathrm{S})}\
\in\mathbb{R}^{M\times M}$ is defined as

$\displaystyle[\mathbf{F}^{(\mathrm{S})}]_{k,n}=$ $\displaystyle\
\begin{cases}\sqrt{\frac{1}{M}}\sin\left(\pi\left(n+\frac{1}{2}\right)\right)&(k=0)\\\
\sqrt{\frac{2}{M}}\sin\left(\frac{\pi}{M}k\left(n+\frac{1}{2}\right)\right)&(k\neq
0)\end{cases}.$ (8)

$\mathbf{F}^{(\mathrm{S})}$ is nothing but the row-wise permuted version of
the DST. In this paper, we simply denote the row-wise permuted DST as the DST.
Because the DCT ($\mathbf{F}^{(\mathrm{C})}$) and the DST
($\mathbf{F}^{(\mathrm{S})}$) are orthogonal matrices, the DADCF is a Parseval
block frame:
$\mathbf{F}^{(\mathrm{D})\top}\mathbf{F}^{(\mathrm{D})}=\mathbf{I}_{M^{2}}$.

The construction flow of the DADCF is illustrated in Fig. 4(a). The DADCF
requires two block transforms, additions and subtractions between two
transforms, and scaling operations. Its computational cost is slightly higher
than conventional block transforms due to the SAP operations but much lower
than other overlapped frames and dictionaries, as mentioned in Section I. Its
redundancy ratio is 2: It is the same as the DFT and the DTCWTs [12, 32, 14,
15], and thus it can reduce the amount of memory usage compared with highly
redundant frames and dictionaries like those in [22, 23].

###### Remark 1

According to the basic knowledge on the DCT/DST, the DST
$\mathbf{F}^{(\mathrm{S})}\in\mathbb{R}^{M\times M}$ can be implemented as the
permuted and sign-altered version of the DCT
$\mathbf{F}^{(\mathrm{C})}\in\mathbb{R}^{M\times M}$, i.e.,
$\mathbf{F}^{(\mathrm{S})}=\mathbf{P}^{(\mathrm{II})}\mathbf{F}^{(\mathrm{C})}\mathrm{diag}(1,-1,\ldots,1,-1)$,
where $\mathbf{P}^{(\mathrm{II})}\in\mathbb{R}^{M\times M}$ denotes the
permutation matrix that arranges the rows of matrices in reverse order. Thus,
the DADCF can be implemented by the DCT with a few trivial SAP operations.

(a)

(b) ${C}_{n_{v},n_{h}}^{(k_{v},k_{h},1)}$,
${B}_{n_{v},n_{h}}^{(k_{v},k_{h},-1)}$

(c) ${S}_{n_{v},n_{h}}^{(k_{v},k_{h},1)}$,
${B}_{n_{v},n_{h}}^{(k_{v},k_{h},1)}$

Figure 4: (a) Procedure of the DADCF ($M=4$). (b) and (c): Atoms
${C}_{n_{v},n_{h}}^{(k_{v},k_{h},1)}$, ${S}_{n_{v},n_{h}}^{(k_{v},k_{h},1)}$
(red), and ${B}_{n_{v},n_{h}}^{(k_{v},k_{h},\pm 1)}$ (blue and green) in the
DADCF. The numbers indicate the rightmost subband indices in (a).

###  IV-B Directional Atoms in DADCF

Here, we examine the directional selectivity of the DADCF defined in (1). The
frequency spectra of the $k$-th rows of the DCT (1) and the DST (8) are given
by
$H_{k}({\omega}):=\mathcal{F}[[\mathbf{F}^{(\mathrm{C})}]_{k,\cdot}],G_{k}({\omega}):=\mathcal{F}[[\mathbf{F}^{(\mathrm{S})}]_{k,\cdot}]$,
where $k\geq 1$. Their complex combination

$\displaystyle
H_{k}({\omega})+jG_{k}({\omega})=\sqrt{\frac{2}{M}}\sum_{n=0}^{M-1}e^{j\theta_{k,n}}e^{-j\omega
n},$ (9)

which is the spectra of (9), approximately satisfies the analyticity, as shown
in Fig. 5(c). As a result, the DADCF is a directional transform with real
coefficients from the 2D DCT and DST followed by addition/subtraction
operations. Note that the frequency spectrum of
$H_{0}({\omega})+jG_{0}({\omega})$, i.e., low-pass spectrum, does not satisfy
the analyticity. As a result, the DADCF can distinguish $2(M-1)^{2}$
directional subbands.

(a) $H_{k}({\omega})$

(b) $G_{k}({\omega})$

(c) $H_{k}({\omega})+jG_{k}({\omega})$

Figure 5: Frequency spectra (frequency: $[0,2\pi]$, $M=8$): (a) DCT, (b) DST,
(c) the complex combination ($1\leq k\leq 7$).

Next, we show the atoms of the DADCF. Because the DADCF forms a Parseval block
frame, it is enough to examine the synthesis transform
$\begin{bmatrix}\mathbf{f}_{0}&\ldots&\mathbf{f}_{2M^{2}-1}\end{bmatrix}^{\top}:=\mathbf{F}^{(\mathrm{D})\top}$.
From (1), $\mathbf{F}^{(\mathrm{D})\top}$ is composed of 1) an atom in the 2D
DCT basis, 2) an atom in the 2D DST basis, or 3) directional atoms arising
from the addition/subtraction of 2D DCT/DST atoms. Let
$\mathbf{B}^{(k_{v},k_{h},1)},\
\mathbf{B}^{(k_{v},k_{h},-1)}\in\mathbb{R}^{M\times M}$ be two directional
atoms of the DADCF that correspond to the subband
$(k_{v},k_{h})\in\Omega_{1,M-1}\times\Omega_{1,M-1}$. These atoms can be
represented as

$\displaystyle{B}^{(k_{v},k_{h},\pm 1)}_{n_{v},n_{h}}=$ $\displaystyle\
{C}^{(k_{v},k_{h})}_{n_{v},n_{h}}\pm{S}^{(k_{v},k_{h})}_{n_{v},n_{h}}$
$\displaystyle=$ $\displaystyle\
\frac{2}{M}\cos\left(\theta_{k_{v},n_{v}}\mp\theta_{k_{h},n_{h}}\right),$ (10)

where
${C}^{(k_{v},k_{h})}_{n_{v},n_{h}}=[\mathbf{F}^{(\mathrm{C})}]_{k_{v},n_{v}}[\mathbf{F}^{(\mathrm{C})}]_{k_{h},n_{h}}$
and
${S}^{(k_{v},k_{h})}_{n_{v},n_{h}}=[\mathbf{F}^{(\mathrm{S})}]_{k_{v},n_{v}}[\mathbf{F}^{(\mathrm{S})}]_{k_{h},n_{h}}$.
In contrast to the DFT and the DHT bases (3) and (III-A), these 2D atoms lie
along various oblique directions, as illustrated in Figs. 4(b) and (c).

###  IV-C Lack of Regularity of DADCF

As previously shown, some 2D frequency responses
$U_{k_{v},k_{h}}({\bm{\omega}})+\overline{U_{k_{v},k_{h}}({\bm{\omega}})}$ and
$U_{k_{v},\overline{k_{h}}}({\bm{\omega}})+\overline{U_{k_{v},\overline{k_{h}}}({\bm{\omega}})}$
obtained from the DADCF do not decay at ${\bm{\omega}}=(0,0)$ which leads to
DC leakage.

Figs. 6(a) and (b) show an image decomposition example. The image used is
Zoneplate $\\{\mathbf{X}^{(i,j)}\\}_{i,j\in\Omega_{31}}$
($\mathbf{X}\in\mathbb{R}^{256\times 256}$) and its (half of the arranged)
DADCF coefficients $\\{\mathbf{x}_{2}\\}_{i,j\in\Omega_{31}}$ with $M=8$,
where $[\mathbf{x}_{1}^{\top}\
\mathbf{x}_{2}^{\top}]^{\top}=\mathbf{F}^{(\mathrm{D})}\
\mathrm{vec}(\mathbf{X}^{(i,j)})$ are shown in Fig. 6(b). We observe that the
DC leakage has been appeared and it leads to the reduction of the sparsity.

The DC leakage is due to the fact that the DST $\mathbf{F}^{(\mathrm{S})}$
loses regularity, as mathematically explained in the following. For a block
transform $\mathbf{F}\in\mathbb{R}^{M\times M}$, regularity condition [38] is
formulated as

$\displaystyle\begin{bmatrix}c&0&\cdots&0\end{bmatrix}^{\top}=\mathbf{F}\mathbf{1},$
(11)

where $c$ is some constant and
$\mathbf{1}=\begin{bmatrix}1&1&\ldots&1\end{bmatrix}^{\top}$. As shown Fig.
5(b), the DST $\mathbf{F}^{(\mathrm{S})}$ leads to the DC leakage. It can be
theoretically verified as in the following proposition.

###### Proposition 1

Let vectors $\\{\mathbf{s}_{k}\\}_{k=0}^{M-1}$ be the basis of $M\times M$
DST, i.e.,
$\begin{bmatrix}\mathbf{s}_{0}&\ldots&\mathbf{s}_{M-1}\end{bmatrix}=\mathbf{F}^{(\mathrm{S})\top}$.
Then,

$\displaystyle\langle\mathbf{s}_{k},\mathbf{1}\rangle=\begin{cases}\frac{\sqrt{2}}{\sqrt{M}\sin\left(\frac{\pi}{2M}k\right)}&(k=2\ell+1)\\\
0&(\mathrm{otherwise})\end{cases},$

where $\ell\in\Omega_{\frac{M}{2}-1}$.

###### Proof:

It is clear that $\langle\mathbf{s}_{0},\mathbf{1}\rangle=0$. For the other
cases,

$\displaystyle\langle\mathbf{s}_{k},\mathbf{1}\rangle=$ $\displaystyle\
\sqrt{\frac{2}{M}}\sum_{n=0}^{M-1}\sin\left(\frac{\pi}{M}k\left(n+\frac{1}{2}\right)\right)$
$\displaystyle=$ $\displaystyle\
\sqrt{\frac{2}{M}}\mathcal{I}\left[\sum_{n=0}^{M-1}e^{j\frac{\pi}{M}k\left(n+\frac{1}{2}\right)}\right]=\frac{(1-(-1)^{k})}{\sqrt{2M}\sin\left(\frac{\pi}{2M}k\right)},$

where $\mathcal{I}$ takes the imaginary part of a complex number. ∎

From the above proposition, the odd rows $(k=2\ell+1)$ of the DST produce
nonzero responses for a constant-valued signal, i.e., DC leakage.

###  IV-D DADCF Pyramid

To obtain sparser coefficients, we introduce the DADCF pyramid inspired by
[6]. The analysis operator of the DADCF pyramid $\mathbf{F}^{(\mathrm{DP})}$
is defined as:

$\displaystyle\mathbf{F}^{(\mathrm{DP})}=\begin{bmatrix}(\mathbf{D}\mathbf{M})^{\top}&(\mathbf{F}^{(\mathrm{D})}(\mathbf{I}-\widetilde{\mathbf{M}}\mathbf{D}^{\top}\mathbf{D}\mathbf{M}))^{\top}\end{bmatrix}^{\top},$
(12)

where $\mathbf{D}=\begin{bmatrix}1&0&\cdots&0\end{bmatrix}$ is the
downsampling operator (and thus $\mathbf{D}^{\top}$ corresponds to the
upsampling operator),
$\mathbf{M}=\mathbf{M}^{(\mathrm{0})}\otimes\mathbf{M}^{(\mathrm{0})}\
\in\mathbb{R}^{M^{2}\times M^{2}}$ is the averaging operator, where
$[\mathbf{M}^{(\mathrm{0})}]_{k,n}=\frac{1}{M}$,
$\widetilde{\mathbf{M}}=M{\mathbf{M}}$. By applying the DADCF pyramid to the
input block $\mathrm{vec}(\mathbf{X}^{(i,j)})$, we can obtain its average
value (denoted as ${x}_{L}$) and the DADCF coefficients of the DC-subtracted
input block (denoted as $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$) as
$\begin{bmatrix}{{x}}_{L}&{\mathbf{x}}_{1}^{\top}&{\mathbf{x}}_{2}^{\top}\end{bmatrix}^{\top}=\mathbf{F}^{(\mathrm{DP})}\mathrm{vec}(\mathbf{X}^{(i,j)})$.

For example, Fig. 6(c) shows the (half of) the transformed coefficients
$\\{\mathbf{x}_{2}\\}_{i,j\in\Omega_{31}}$. It is clear that sparser
coefficients can be obtained and the DADCF pyramid
$\mathbf{F}^{(\mathrm{DP})}$ is still invertible. By this operation, however,
the number of transformed coefficients is slightly increased from $2N^{2}$ to
$2N^{2}+(N/M)^{2}$ for $N\times N$ input images.

(a)

(b)

(c)

Figure 6: (a): Zoneplate, (b) (Half of) DADCF coefficients ($M=8$), (c) (Half
of) DADCF coefficients in the DADCF pyramid ones ($M=8$).

##  V Regularity-Constrained DADCF

In this section, we introduce another DADCF, called RDADCF. We introduce a
RDST in Section V-A and V-B. Then, in Section V-C, we propose the RDADCF,
which overcomes the problem of the DADCF, i.e., DC leakage, and saves the
number of the transformed coefficients fewer than the DADCF pyramid.

###  V-A Design of RDST

This section introduces a modified DST without DC leakage for constructing
RDADCF. For notation simplicity, we present steps for constructing the RDST
matrix $\mathbf{F}^{(\mathrm{RS})}$.

Step 1:

First, we define a modified DST $\mathbf{S}\in\mathbb{R}^{M\times M}$.

$\displaystyle[\mathbf{S}]_{k,n}=$ $\displaystyle\
\begin{cases}\sqrt{\frac{1}{M}}&(k=0)\\\
\sqrt{\frac{2}{M}}\sin\left(\frac{\pi}{M}k\left(n+\frac{1}{2}\right)\right)&(k\neq
0)\end{cases}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}\mathbf{s}_{0}&\mathbf{s}_{1}&\cdots&\mathbf{s}_{M-1}\end{bmatrix}^{\top}.$
(13)

In short, it is constructed by replacing the 0-th row of the DST with that of
the DCT. The modified DST satisfies the following property (see Appendix A for
its proof).

###### Proposition 2

$\mathrm{rank}(\mathbf{S})=M-1$.

Then, we further modify $\mathbf{S}$ in (Step 1:). From (11), in order to
impose the regularity condition on $\mathbf{S}$,
$\\{\mathbf{s}_{k}\\}_{k=1}^{M-1}$ should be orthogonal to $\mathbf{s}_{0}$.
Now, we orthogonalize the odd rows of $\mathbf{S}$ in the following way.

Step 2:

Set
$\widetilde{\mathbf{S}}^{(0)}=\begin{bmatrix}\mathbf{s}_{0}&\mathbf{0}&\mathbf{s}_{2}&\ldots&\mathbf{s}_{M-1}\end{bmatrix}^{\top}$.

Here, $\widetilde{\mathbf{S}}^{(0)}$ satisfies the following proposition (see
Appendix B for its proof).

###### Proposition 3

$\mathrm{rank}(\widetilde{\mathbf{S}}^{(0)})=M-1$.

From Proposition 3, there is only one zero singular value and its
corresponding right-singular vector (denoted as $\mathbf{v}^{(0)}$) belongs to
the null space of $\widetilde{\mathbf{S}}^{(0)}$. It implies that
$\widetilde{\mathbf{S}}^{(0)}\mathbf{v}^{(0)}=\mathbf{0}$, i.e.,
$\mathbf{v}^{(0)}$ satisfies the regularity condition.
$\widetilde{\mathbf{S}}^{(0)}$ is updated by replacing $\mathbf{0}$ to
$\mathbf{v}^{(0)}$.

Step 3:

Set
$\mathbf{S}^{(0)}=\begin{bmatrix}\mathbf{s}_{0}&\mathbf{v}^{(0)}&\mathbf{s}_{2}&\ldots&\mathbf{s}_{M-1}\end{bmatrix}^{\top}$.

Note that $\mathbf{v}^{(0)}$ can be explicitly represented as
$[\mathbf{v}^{(0)}]_{n}=\sqrt{\frac{1}{M}}(-1)^{n}=[\mathbf{F}^{(\mathrm{S})}]_{0,n}\left(=\sqrt{\frac{1}{M}}\sin\left(\pi\left(n+\frac{1}{2}\right)\right)\right)$
because the row of the DST $[\mathbf{F}^{(\mathrm{S})}]_{0,n}$ corresponding
to the highest frequency subband is orthogonal to $\\{\mathbf{s}_{0},\
\mathbf{s}_{2},\ \ldots,\ \mathbf{s}_{M-1}\\}$. It clearly follows that
$\mathrm{rank}(\mathbf{S}^{(0)})=M$.

Consequently, by repeating Steps 1 and 2, we can obtain the orthogonal matrix
$\mathbf{S}^{(M/2-1)}$ whose odd rows are replaced by the different ones from
the initial $\mathbf{S}^{(0)}$. A summary of the algorithm is given in
Algorithm 1. Hereafter $\mathbf{F}^{(\mathrm{RS})}:=\mathbf{S}^{(M/2-1)}$ is
termed as a RDST.

Algorithm 1 The design procedure for RDST

1: Set $\mathbf{S}$ is as in (Step 1:).

2: for $k=0$ to $M/2-1$ do

3: Set
$\widetilde{\mathbf{S}}^{(k)}=\begin{bmatrix}\ldots&\mathbf{s}_{2k}&\mathbf{0}&\mathbf{s}_{2k+2}&\ldots\end{bmatrix}^{\top}$.

4: Find the right-singular vector $\mathbf{v}^{(k)}$ corresponding to zero
singular value .

5: Set
$\mathbf{S}^{(k)}=\begin{bmatrix}\ldots&\mathbf{s}_{2k}&\mathbf{v}^{(k)}&\mathbf{s}_{2k+2}&\ldots\end{bmatrix}^{\top}$.

6: end for

7: Output $\mathbf{S}^{(M/2-1)}$.

The RDST satisfies the following properties (see Appendix C for its proof).

###### Proposition 4

Let $\mathbf{F}^{(\mathrm{RS})}\in\mathbb{R}^{M\times M}$ be the RDST.

  1. 1.

This satisfies the regularity condition, i.e.,

$\begin{bmatrix}c&0&\cdots&0\end{bmatrix}^{\top}=\mathbf{F}^{(\mathrm{RS})}\mathbf{1}.$

  2. 2.

Some rows of $\mathbf{F}^{(\mathrm{RS})}\in\mathbb{R}^{M\times M}$ are
identical with those in the DST matrix $\mathbf{F}^{(\mathrm{S})}$:
$[\mathbf{F}^{(\mathrm{RS})}]_{0,n}=\sqrt{\frac{1}{M}}$,
$[\mathbf{F}^{(\mathrm{RS})}]_{1,n}=[\mathbf{F}^{(\mathrm{S})}]_{0,n}$,
$[\mathbf{F}^{(\mathrm{RS})}]_{2\ell,n}=[\mathbf{F}^{(\mathrm{S})}]_{2\ell,n},$
where $\mathbf{F}^{(\mathrm{S})}\in\mathbb{R}^{M\times M}$ is the DST matrix
and $\ell\in\Omega_{\frac{M}{2}-1}$.

  3. 3.

The passband of the spectrum
$\mathcal{F}[[\mathbf{F}^{(\mathrm{RS})}]_{2\ell+1,\cdot}]$ is the same as
that of the DST $\mathcal{F}[[\mathbf{F}^{(\mathrm{S})}]_{2\ell+1,\cdot}]$
$(\ell\geq 2)$.

In Fig. 7(a), the red lines show the frequency spectra of the newly updated
rows ($k=0,3,5,7$) in the RDST ($M=8$) and the dashed gray lines show those of
the corresponding rows in the DST (the rest frequency spectra of the RDST are
identical to those of the DST). The frequency spectra of the RDST approximate
those of the original DST, but decay at zero frequency.

(a)

(b)

Figure 7: Frequency spectra (frequency: $[0,2\pi]$, $M=8$): (a) red lines:
$\mathcal{F}[[\mathbf{F}^{(\mathrm{RS})}]_{k,\cdot}]$, dashed gray lines:
$\mathcal{F}[[\mathbf{F}^{(\mathrm{S})}]_{k,\cdot}]$ ($k=0,3,5,7$), (b) red
lines:
$\mathcal{F}[[\mathbf{F}^{(\mathrm{C})}]_{k,\cdot}]+j\mathcal{F}[[\mathbf{F}^{(\mathrm{RS})}]_{k,\cdot}]$,
dashed gray lines:
$\mathcal{F}[[\mathbf{F}^{(\mathrm{C})}]_{k,\cdot}]+j\mathcal{F}[[\mathbf{F}^{(\mathrm{S})}]_{k,\cdot}]$
($k=3,5,7$).

###  V-B Implementation of RDST

From Proposition 4, the $\frac{M}{2}$ rows of the RDST
$\mathbf{F}^{(\mathrm{RS})}\in\mathbb{R}^{M\times M}$ are the same as the rows
of the original DST $\mathbf{F}^{(\mathrm{S})}\in\mathbb{R}^{M\times M}$, and
both matrices are orthogonal. Thus, we can derive that the RDST can be
implemented by the cascade of the DST and an orthogonal matrix as in the
following.

Let $\mathbf{F}^{(\mathrm{S,e})},\
\mathbf{F}^{(\mathrm{S,o})}\in\mathbb{R}^{\frac{M}{2}\times M}$ be the even
and odd rows of the $\mathbf{F}^{(\mathrm{S})}\in\mathbb{R}^{M\times M}$,
respectively. Then, the RDST $\mathbf{F}^{(\mathrm{RS})}\in\mathbb{R}^{M\times
M}$ can be expressed as:

$\displaystyle\mathbf{F}^{(\mathrm{RS})}=$ $\displaystyle\
\mathbf{P}^{(\mathrm{III})}\begin{bmatrix}\mathbf{F}^{(\mathrm{S,e})}\\\
\widetilde{\mathbf{F}}^{(\mathrm{S,o})}\end{bmatrix}=\mathbf{P}^{(\mathrm{III})}\mathrm{diag}(\mathbf{I}_{\frac{M}{2}},\mathbf{\Gamma}_{\frac{M}{2}})\begin{bmatrix}\mathbf{F}^{(\mathrm{S,e})}\\\
\mathbf{F}^{(\mathrm{S,o})}\end{bmatrix}$ $\displaystyle=$ $\displaystyle\
\mathbf{P}^{(\mathrm{III})}\mathrm{diag}(\mathbf{I}_{\frac{M}{2}},\mathbf{\Gamma}_{\frac{M}{2}})\mathbf{P}^{(\mathrm{IV})}{\mathbf{F}}^{(\mathrm{S})},$
(14)

where $\mathbf{P}^{(\mathrm{III})},\
\mathbf{P}^{(\mathrm{IV})}\in\mathbb{R}^{M\times M}$ are the permutation
matrices, and the matrix $\mathbf{\Gamma}_{\frac{M}{2}}$ is guaranteed to be
an orthogonal matrix because of orthogonality of the RDST and the DST. Since
$\mathbf{\Gamma}_{\frac{M}{2}}$ is an orthogonal matrix, it can be factorized
into $\frac{M(M-2)}{8}$ rotation matrices. Thus, the RDST is still a hardware-
friendly transform that can be implemented by the $\mathbf{F}^{(\mathrm{C})}$
with some trivial operations.

###  V-C Design of RDADCF

Finally, a RDADCF $\mathbf{F}^{(\mathrm{RD})}$ is defined using the RDST as
follows.

###### Definition 2

Let $\mathbf{F}^{(\mathrm{RD})}\in\mathbb{R}^{2M^{2}\times M^{2}}$ be the
analysis operator of the RDADCF defined as:

$\displaystyle\mathbf{F}^{(\mathrm{RD})}:=$ $\displaystyle\
\mathbf{P}^{(\mathrm{V})\top}\mathbf{W}^{(\mathrm{II})}\mathbf{P}^{(\mathrm{V})}\begin{bmatrix}\mathbf{F}^{(\mathrm{C})}\otimes\mathbf{F}^{(\mathrm{C})}\\\
\mathbf{F}^{(\mathrm{RS})}\otimes\mathbf{F}^{(\mathrm{RS})}\end{bmatrix},$
$\displaystyle\mathbf{W}^{(\mathrm{II})}=$ $\displaystyle\
\mathrm{diag}\left(\frac{1}{\sqrt{2}}\mathbf{I}_{4M-4},\frac{1}{2}\begin{bmatrix}\mathbf{I}_{(M-2)^{2}}&-\mathbf{I}_{(M-2)^{2}}\\\
\mathbf{I}_{(M-2)^{2}}&\mathbf{I}_{(M-2)^{2}}\end{bmatrix}\right),$ (15)

where $\mathbf{P}^{(\mathrm{V})}\in\mathbb{R}^{2M^{2}\times 2M^{2}}$ is a
permutation matrix. $\mathbf{P}^{(\mathrm{V})}$ places the $4M-4$ DCT and DST
coefficients associated with the subband indices $k_{v}\in\\{0,1\\}$ or
$k_{h}\in\\{0,1\\}$ to the first part, and the other $2(M-2)^{2}$ coefficients
to the last (see Fig. 8(a)). Due to the orthogonality of the RDST, the RDADCF
clearly forms a Parseval block frame, i.e.,
$\mathbf{F}^{(\mathrm{RD})\top}\mathbf{F}^{(\mathrm{RD})}=\mathbf{I}_{M}$.

Now, we discuss the capability of the directional subband decomposition based
on the DCT and the RDST. Let $\mathbf{F}^{(\mathrm{C})},\
\mathbf{F}^{(\mathrm{RS})}\in\mathbb{R}^{M\times M}$ be the DCT and the RDST
matrices. As discussed in Section III-B, the complex combination
$[\mathbf{F}^{(\mathrm{C})}]_{k,\cdot}\pm
j[\mathbf{F}^{(\mathrm{RS})}]_{k,\cdot}$ should have a one-sided frequency
spectrum for directional subband decomposition. In the case of even $k\ (\geq
2)$, the rows of $[\mathbf{F}^{(\mathrm{RS})}]_{k,\cdot}$ are identical to
those of the DST. Therefore, the frequency spectrum
$[\mathbf{F}^{(\mathrm{C})}]_{k,\cdot}\pm
j[\mathbf{F}^{(\mathrm{RS})}]_{k,\cdot}$ is one-sided. In the case of odd $k\
(\geq 3)$, where the rows $[\mathbf{F}^{(\mathrm{RS})}]_{k,\cdot}$ are newly
designed in Algorithm 1, the frequency spectrum
$[\mathbf{F}^{(\mathrm{C})}]_{k,\cdot}\pm
j[\mathbf{F}^{(\mathrm{RS})}]_{k,\cdot}$ can be one-sided (Fig. 7(b)).

Analyticity of the RDADCF can be explained as follows. Let
$\\{\mathbf{s}_{\ell}^{(s)}\\}$ and $\\{\mathbf{s}_{\ell}^{(r)}\\}$ be the
rows of the DST and the RDST, respectively. For any odd $k$ ($\geq 3$),
$\mathbf{s}_{k}^{(r)}$ can be obtained by applying orthogonal projection to
$\mathbf{s}_{k}^{(s)}$ onto the orthogonal complement of
$\\{\mathbf{s}_{\ell}^{(s)}\\}_{\Omega_{M-1}\setminus\\{k\\}}$, as

$\displaystyle\mathbf{s}_{k}^{(r)}=\frac{\pm
1}{\eta_{k}}\left(\mathbf{s}_{k}^{(s)}-\sum_{\ell\in\Omega_{M-1}\setminus\\{k\\}}\langle\mathbf{s}_{\ell}^{(r)},\mathbf{s}_{k}^{(s)}\rangle\mathbf{s}_{\ell}^{(r)}\right),$
(16)

where $\eta_{k}$ is the normalization factor for $\mathbf{s}_{k}^{(r)}$ having
unit norm. Let
$[\mathbf{F}^{(\mathrm{W})}]_{k,n}=\frac{1}{\sqrt{M}}e^{-j\frac{2\pi}{M}kn}$
denote the DFT. Because $\mathbf{F}^{(\mathrm{W})}\mathbf{s}_{m}^{(r)}$ and
$\mathbf{F}^{(\mathrm{W})}\mathbf{s}_{k}^{(r)}$ have different passbands,
$|\langle\mathbf{s}_{m}^{(r)},\mathbf{s}_{k}^{(r)}\rangle|=|\langle\mathbf{F}^{(\mathrm{W})}\mathbf{s}_{m}^{(r)},\mathbf{F}^{(\mathrm{W})}\mathbf{s}_{k}^{(r)}\rangle|$
is small. Therefore, the spectrum of the $\mathbf{s}_{k}^{(r)}$ can
approximate $\mathbf{s}_{k}^{(s)}$ over the passband of
$\mathbf{s}_{k}^{(s)}$.

The atoms of the RDADCF lie along the $2(M-2)^{2}$ frequency directions, as
shown in Figs. 8(b) and (c), where
${C}^{(k_{v},k_{h})}_{n_{v},n_{h}}=[\mathbf{F}^{(\mathrm{C})}]_{k_{v},n_{v}}[\mathbf{F}^{(\mathrm{C})}]_{k_{h},n_{h}}$,
${S}^{(k_{v},k_{h})}_{n_{v},n_{h}}=[\mathbf{F}^{(\mathrm{RS})}]_{k_{v},n_{v}}[\mathbf{F}^{(\mathrm{RS})}]_{k_{h},n_{h}}$,
and ${B}^{(k_{v},k_{h},\pm
1)}_{n_{v},n_{h}}={C}^{(k_{v},k_{h})}_{n_{v},n_{h}}\pm{S}^{(k_{v},k_{h})}_{n_{v},n_{h}}$.
The number of directional selectivities of the RDADCF is slightly less than
the original DADCF. Since the RDADCF with $M=2$ cannot ensure directional
selectivity, we recommend $M=2^{m}$ where $m\geq 2$.

(a)

(b) ${C}_{n_{v},n_{h}}^{(k_{v},k_{h},1)}$,
${B}_{n_{v},n_{h}}^{(k_{v},k_{h},-1)}$

(c) ${S}_{n_{v},n_{h}}^{(k_{v},k_{h},1)}$,
${B}_{n_{v},n_{h}}^{(k_{v},k_{h},1)}$

Figure 8: (a) Procedure of the RDADCF ($M=4$). (b) and (c): Atoms
${C}_{n_{v},n_{h}}^{(k_{v},k_{h})}$, ${S}_{n_{v},n_{h}}^{(k_{v},k_{h})}$
(red), and ${B}_{n_{v},n_{h}}^{(k_{v},k_{h},\pm 1)}$ (blue and green) in the
RDADCF. The numbers indicate the rightmost subband indices in (a).

##  VI Experimental results

We evaluated the performance of the proposed DADCF pyramid (Section IV-D) and
RDADCF in compressive image sensing reconstruction [4], as an example of image
processing applications. $512\times 512$ pixel images in Fig. 9 were used as
the test set. Each incomplete observation
($\mathbf{y}=\mathrm{vec}(\mathbf{Y})\in\mathbb{R}^{512^{2}}$) is obtained by
Noiselet transform [39] ($\mathbf{\Phi}\in\mathbb{R}^{512^{2}\times 512^{2}}$)
followed by random sampling of 30%, 40%, 50%, and 60% pixels
($\mathbf{R}_{\mathrm{samp}}\in\mathbb{R}^{\mathcal{R}(512^{2}p)\times
512^{2}}$ where $\mathcal{R}$ is the rounding operator and $p=0.3,\ 0.4,\
0.5,\ 0.6$) in the presence of additive white Gaussian noise
($\mathbf{n}\in\mathbb{R}^{512^{2}}$) with the standard derivation
$\sigma=0.1$ as
$\mathbf{y}=\mathbf{R}_{\mathrm{samp}}\mathbf{\Phi}\mathbf{x}+\mathbf{n}$,
($\mathbf{x}=\mathrm{vec}(\mathbf{X})\in\mathbb{R}^{512^{2}}$). Figs. 9(a)–(d)
indicate the estimated latent images by using the Moore-Penrose pseudo inverse
of
$\widetilde{\mathbf{\Phi}}^{\dagger}=\mathbf{\Phi}^{\top}\mathbf{R}_{\mathrm{samp}}^{\top}$
($\widetilde{\mathbf{\Phi}}=\mathbf{R}_{\mathrm{samp}}\mathbf{\Phi}$) in the
case of $p=0.3$.

Up to now, many block transform-based methods for image recovery have been
proposed, such as BM3D, patch-based redundant DCT approaches, and so on [25,
26, 27]. For fair comparison, we simply evaluate directional block transforms
in two image recovery problems:

  * •

Problem 1: image recovery based on sparsity of block-wise transformed
coefficients.

  * •

Problem 2: image recovery based on sparsity of block-wise transformed
coefficients and weighted total variation (WTV) for block boundaries as
presented in [40].

The cost function for these two problems is described as follows:

$\displaystyle\mathbf{x}^{\star}=\mathop{\rm
arg~{}min}\limits_{\mathbf{x}\in\mathbb{R}^{512^{2}}}$
$\displaystyle\|\mathbf{F}\mathbf{P}_{\mathrm{v2bv}}{\mathbf{x}}\|_{1}+\rho\|\widetilde{\mathbf{W}}_{\mathrm{b}}\mathbf{D}_{\mathrm{hv}}\mathbf{x}\|_{1,2}$
$\displaystyle+\iota_{C_{[0,1]}}(\mathbf{x})+F_{\mathbf{y}}(\widetilde{\mathbf{\Phi}}\mathbf{x}),$
(17)

where $\mathbf{P}_{\mathrm{v2bv}}$ is the permutation matrix permuting a
vectorized version of a matrix to a block-wise-vectorized one
$\mathbf{P}_{\mathrm{v2bv}}{\mathbf{x}}=\mathrm{bvec}(\mathbf{X})$,
$\mathbf{F}=\mathbf{I}^{(\mathrm{V})}\otimes\mathbf{F}^{(\mathrm{DP})}$ or
$\mathbf{F}=\mathbf{I}^{(\mathrm{V})}\otimes\mathbf{F}^{(\mathrm{RD})}$
($\mathbf{I}^{(\mathrm{V})}=\mathbf{I}_{512^{2}/M^{2}}$), and
$\iota_{A}(\mathbf{x})$ is the indicator function555Indicator function of set
$A$ is defined as $\iota_{A}(\mathbf{x})=0,\ (\mathbf{x}\in
A),\quad\iota_{A}(\mathbf{x})=\infty,\ (\mathbf{x}\notin A)$. of a set $A$.
$C_{[0,1]}$ is the set of vectors whose entries are within $[0,1]$. The data-
fidelity function was set as
$F_{\mathbf{y}}=\iota_{\mathcal{B}(\mathbf{y},\epsilon)}$
($\mathcal{B}(\mathbf{y},\epsilon):=\\{\mathbf{x}\in\mathbb{R}^{M}|\|\mathbf{x}-\mathbf{y}\|_{2}\leq\epsilon\\}$)
is the indicator function defined by the $\ell_{2}$-norm ball. The radius was
set as $\epsilon=\|\mathbf{x}_{o}-\mathbf{y}\|_{2}$, where $\mathbf{x}_{o}$ is
an original image.
$\mathbf{D}_{\mathrm{hv}}=\begin{bmatrix}\mathbf{D}_{\mathrm{v}}^{\top}&\mathbf{D}_{\mathrm{h}}^{\top}\end{bmatrix}^{\top}\in\mathbb{R}^{(2\cdot
512^{2})\times 512^{2}}$ denotes the vertical and horizontal difference
operator.
$\widetilde{\mathbf{W}}_{b}=\mathbf{I}^{(\mathrm{V})}\otimes\mathbf{W}_{b}$,
where $\mathbf{W}_{b}\in\mathbb{R}^{M\times M}$ is the weighting matrix for
block boundary as $[\mathbf{W}_{b}]_{m,n}=0$ $(1\leq m,n\leq M-2)$,
$[\mathbf{W}_{b}]_{m,n}=1$ (otherwise) . The cost functions with $\rho=0$ and
$\rho=1$ correspond to Problem 1 and 2, respectively. The detailed algorithm
used in the experiments is given in Appendix E.

For comparison, we also used the DCT, the DFT, and the DHT in (VI). The block
size is set to $M=8,\ 16,\ 32$.

Fig. 10 shows the resulting images of the proposed and conventional transforms
obtained in the case of sampling rate $p=0.4$. As these figures show
(particularly in the dashed red boxes), the DCT cannot recover directional
textures precisely. Table II shows the numerical results. In most cases, the
DADCF pyramid or the RDADCF outperformed the DCT, the DFT, and the DHT in
terms of the reconstruction error (PSNR). The RDADCF recovers the images
better than the DADCF pyramid, especially for Monarch and Parrot (smooth
images), due to its regularity property. In fact, the DCT is superior to the
DADCF pyramid and the RDADCF in some cases. However, since the DADCF pyramid
and the RDADCF are compatible with the DCT, we can select the DCT, the DADCF
pyramid, or the RDADCF by using or bypassing the DST/RDST and the SAP
operations, depending on the input image.

(a) Barbara

(b) Mandrill

(c) Monarch

(d) Parrot

Figure 9: (a)–(d): Original images (256 $\times$ 256) and recovered images by
$\widetilde{\Phi}^{\dagger}$ (sampling rate $p=0.3$).

(a) DCT

(b) DFT

(c) DADCF

(d) RDADCF

(e) DCT+WTV

(f) DFT+WTV

(g) DADCF+WTV

(h) RDADCF+WTV

Figure 10: Zoomed resulting images reconstructed from 60% noiselet
coefficients ((a)–(d)) and 30% noiselet coefficients ((e)–(h)). The size and
the decomposition level of the transforms is $M=8$ and $J=2$, respectively.
TABLE II: Numerical results of compressive sensing reconstruction.

---
Problem 1: Image recovery based on sparsity of transformed coefficients |
| Block size: $M=8$ | Block size: $M=16$ | Block size: $M=32$
PSNR [dB] | DCT | DFT | DHT | DADCFP | RDADCF | DCT | DFT | DHT | DADCFP | RDADCF | DCT | DFT | DHT | DADCFP | RDADCF
Image: Barbara |
60%: 10.77 | 32.76 | 31.76 | 31.00 | 32.51 | 32.97 | 33.81 | 32.54 | 31.71 | 33.14 | 33.89 | 33.81 | 33.06 | 32.17 | 33.61 | 34.19
50%: 10.22 | 30.10 | 29.63 | 28.84 | 30.40 | 30.69 | 31.43 | 30.47 | 29.61 | 31.16 | 31.85 | 31.45 | 30.96 | 30.14 | 31.57 | 32.14
40%: 9.689 | 27.59 | 27.58 | 26.78 | 28.47 | 28.62 | 29.10 | 28.44 | 27.62 | 29.25 | 29.91 | 29.26 | 28.87 | 28.04 | 29.75 | 30.25
30%: 9.255 | 25.11 | 25.46 | 24.66 | 26.51 | 26.46 | 26.70 | 26.38 | 25.58 | 27.25 | 27.78 | 27.02 | 26.79 | 25.95 | 27.70 | 28.24
Image: Mandrill |
60%: 10.77 | 25.68 | 25.95 | 25.29 | 26.64 | 26.54 | 25.99 | 26.47 | 25.76 | 27.04 | 27.05 | 25.99 | 26.66 | 25.94 | 27.25 | 27.32
50%: 10.22 | 23.89 | 24.21 | 23.55 | 25.05 | 24.88 | 24.27 | 24.78 | 24.08 | 25.44 | 25.43 | 24.32 | 24.91 | 24.25 | 25.66 | 25.70
40%: 9.689 | 22.18 | 22.57 | 21.97 | 23.59 | 23.34 | 22.68 | 23.16 | 22.55 | 23.94 | 23.91 | 22.77 | 23.30 | 22.68 | 24.15 | 24.20
30%: 9.255 | 20.54 | 21.01 | 20.40 | 22.21 | 21.87 | 21.20 | 21.68 | 21.15 | 22.56 | 22.51 | 21.39 | 21.82 | 21.25 | 22.76 | 22.80
Image: Monach |
60%: 10.77 | 37.43 | 35.40 | 34.37 | 36.31 | 37.22 | 37.23 | 35.14 | 34.09 | 36.09 | 36.77 | 36.52 | 35.39 | 34.36 | 36.19 | 36.59
50%: 10.22 | 34.70 | 32.72 | 31.66 | 33.99 | 34.73 | 34.61 | 32.62 | 31.52 | 33.78 | 34.41 | 34.10 | 32.93 | 31.91 | 33.94 | 34.31
40%: 9.689 | 31.65 | 29.88 | 28.77 | 31.61 | 32.06 | 31.93 | 30.10 | 29.03 | 31.53 | 32.04 | 31.66 | 30.60 | 29.57 | 31.88 | 32.21
30%: 9.255 | 28.23 | 26.92 | 25.79 | 29.05 | 29.19 | 29.04 | 27.44 | 26.44 | 29.07 | 29.44 | 28.92 | 28.09 | 27.15 | 29.49 | 29.77
Image: Parrot |
60%: 10.77 | 39.25 | 38.03 | 37.28 | 38.46 | 39.40 | 39.46 | 37.94 | 37.15 | 38.05 | 39.35 | 39.10 | 37.96 | 37.11 | 37.98 | 39.11
50%: 10.22 | 37.13 | 35.91 | 34.95 | 36.67 | 37.58 | 37.47 | 35.84 | 34.91 | 36.27 | 37.58 | 37.07 | 35.85 | 34.90 | 36.13 | 37.30
40%: 9.689 | 34.57 | 33.48 | 32.49 | 34.54 | 35.35 | 35.19 | 33.57 | 32.57 | 34.23 | 35.56 | 34.78 | 33.63 | 32.62 | 34.11 | 35.28
30%: 9.255 | 31.60 | 30.79 | 29.69 | 32.30 | 32.85 | 32.47 | 31.10 | 30.11 | 32.04 | 33.26 | 32.21 | 31.09 | 30.07 | 31.88 | 33.09

---
Problem 2: Image recovery based on sparsity of transformed coefficients and weighted total variation for block boundaries |
| Block size: $M=8$ | Block size: $M=16$ | Block size: $M=32$
PSNR [dB] | DCT | DFT | DHT | DADCFP | RDADCF | DCT | DFT | DHT | DADCFP | RDADCF | DCT | DFT | DHT | DADCFP | RDADCF
Image: Barbara |
60%: 10.77 | 32.27 | 31.98 | 31.48 | 32.32 | 32.62 | 33.15 | 32.60 | 31.93 | 32.96 | 33.47 | 33.38 | 33.15 | 32.30 | 33.49 | 33.97
50%: 10.22 | 29.95 | 30.00 | 29.52 | 30.29 | 30.53 | 30.89 | 30.69 | 29.98 | 31.07 | 31.54 | 31.09 | 31.14 | 30.32 | 31.53 | 31.99
40%: 9.689 | 27.91 | 28.17 | 27.70 | 28.49 | 28.64 | 28.76 | 28.89 | 28.22 | 29.27 | 29.69 | 29.05 | 29.24 | 28.45 | 29.78 | 30.16
30%: 9.255 | 26.06 | 26.42 | 26.01 | 26.74 | 26.82 | 26.88 | 27.11 | 26.49 | 27.42 | 27.82 | 27.11 | 27.30 | 26.53 | 27.86 | 28.30
Image: Mandrill |
60%: 10.77 | 26.85 | 26.95 | 26.59 | 27.23 | 27.21 | 26.46 | 26.91 | 26.37 | 27.26 | 27.25 | 26.20 | 26.86 | 26.20 | 27.29 | 27.34
50%: 10.22 | 25.24 | 25.42 | 25.07 | 25.68 | 25.65 | 24.89 | 25.35 | 24.82 | 25.67 | 25.67 | 24.61 | 25.22 | 24.63 | 25.76 | 25.79
40%: 9.689 | 23.84 | 24.03 | 23.71 | 24.31 | 24.25 | 23.47 | 23.90 | 23.43 | 24.25 | 24.26 | 23.21 | 23.73 | 23.19 | 24.29 | 24.33
30%: 9.255 | 22.46 | 22.67 | 22.41 | 22.96 | 22.86 | 22.16 | 22.58 | 22.18 | 22.92 | 22.90 | 21.91 | 22.38 | 21.92 | 22.95 | 22.99
Image: Monach |
60%: 10.77 | 40.60 | 39.55 | 39.18 | 39.49 | 40.10 | 39.26 | 38.01 | 37.29 | 38.08 | 38.67 | 37.63 | 36.96 | 36.09 | 37.20 | 37.59
50%: 10.22 | 38.70 | 37.49 | 37.04 | 37.50 | 38.13 | 37.01 | 35.68 | 34.89 | 35.87 | 36.48 | 35.34 | 34.58 | 33.66 | 34.95 | 35.34
40%: 9.689 | 36.49 | 35.12 | 34.62 | 35.21 | 35.82 | 34.69 | 33.27 | 32.46 | 33.56 | 34.12 | 33.04 | 32.23 | 31.37 | 32.85 | 33.17
30%: 9.255 | 33.76 | 32.33 | 31.84 | 32.58 | 33.13 | 31.97 | 30.40 | 29.64 | 30.89 | 31.44 | 30.39 | 29.59 | 28.71 | 30.34 | 30.64
Image: Parrot |
60%: 10.77 | 40.83 | 40.24 | 39.91 | 40.20 | 40.71 | 40.34 | 39.56 | 39.00 | 39.35 | 40.19 | 39.66 | 39.07 | 38.37 | 38.84 | 39.69
50%: 10.22 | 39.33 | 38.73 | 38.34 | 38.68 | 39.27 | 38.72 | 37.85 | 37.19 | 37.72 | 38.62 | 37.89 | 37.21 | 36.40 | 37.09 | 38.02
40%: 9.689 | 37.43 | 36.83 | 36.41 | 36.84 | 37.45 | 36.77 | 35.91 | 35.21 | 35.77 | 36.77 | 35.77 | 35.17 | 34.32 | 35.09 | 36.07
30%: 9.255 | 35.32 | 34.69 | 34.25 | 34.70 | 35.32 | 34.55 | 33.71 | 33.00 | 33.64 | 34.64 | 33.47 | 32.77 | 31.91 | 32.85 | 33.95

##  VII Concluding Remarks

In this paper, we proposed the DADCF and the RDADCF for directional image
representation by extending the DCT. Since they are Parseval block frames with
low redundancy, they can deliver computational efficiency for practical
situations. Furthermore, unlike the conventional directional block transforms,
they can finely decompose the frequency plane and provide richer directional
atoms. Comparing both the DADCF and the RDADCF, the DADCF provides richer
directional selectivity than the RDADCF. However, in practice, the slightly
redundant DADCF pyramid should be used instead of the DADCF to avoid the DC
leakage and perform good image processing, i.e., the RDADCF can save more
amount of memory usage than the DADCF (pyramid). Also, they can be easily
implemented by appending trivial operations (the DST or the RDST, and the SAP
operations) to the DCT. Moreover, the DST can be implemented based on the DCT
with the permutation and sign-alternation operations and the RDST based on the
DST and one additional orthogonal matrix with the size of the half. Since the
DCT is integrated into many existing digital devices, the system modification
for the proposed method is minimal. Note that the DADCF and the RDADCF are
compatible with the DCT. Depending on applications, we can switch the
DCT/DADCF/RDADCF by using the DST/RDST and SAP operations.

We evaluated the DADCF pyramid and the RDADCF in compressive image sensing
reconstruction as a practical application. The experimental results showed
that, for both fine textures and smooth images, they could achieve higher
numerical qualities than the DCT, the DFT, and the DHT. Furthermore, it was
shown that the RDADCF could recover smooth regions better than the DADCF
pyramid due to its regularity property.

##  Appendix A Proof for Proposition 2

###### Proof:

We first introduce the following lemma.

###### Lemma 1

  1. 1.

The elements in the upper-right triangle
$[\mathbf{S}\mathbf{S}^{\top}]_{k_{v},k_{h}}=[\langle\mathbf{s}_{k_{v}},\mathbf{s}_{k_{h}}\rangle]_{k_{v},k_{h}}$
are expressed as

$\displaystyle\langle\mathbf{s}_{k_{v}},\mathbf{s}_{k_{h}}\rangle=\begin{cases}\frac{\sqrt{2}}{M\sin\left(\frac{\pi}{2M}k_{h}\right)}&(k_{v}=0\
\mathrm{and}\ k_{h}=2\ell+1)\\\ 1&(k_{v}=k_{h})\\\
0&(\mathrm{otherwise})\end{cases},$

where $\ell\in\Omega_{\frac{M}{2}-1}$. For example, for $M=4$,

$\displaystyle\mathbf{S}\mathbf{S}^{\top}\approx\begin{bmatrix}1&0.9239&0&0.3827\\\
0.9239&1&0&0\\\ 0&0&1&0\\\ 0.3827&0&0&1\end{bmatrix}.$ (18)

  2. 2.

$\sum^{M/2-1}_{\ell=0}\langle\mathbf{s}_{0},\mathbf{s}_{2\ell+1}\rangle^{2}=1$.

###### Proof:

  1. 1.

It is clear that $\langle\mathbf{s}_{0},\mathbf{s}_{0}\rangle=1$ and
$\langle\mathbf{s}_{k_{v}},\mathbf{s}_{k_{h}}\rangle=\delta(k_{v}-k_{h})$ for
$k_{v},k_{h}\in\Omega_{1,M-1}$ because $\\{\mathbf{s}_{k}\\}_{k=1}^{M-1}$ are
the rows of the DST $\mathbf{F}^{(\mathrm{S})}$. In the other cases, it is
clear from Proposition 1.

  2. 2.

$\sum^{M/2-1}_{\ell=0}\langle\mathbf{s}_{0},\mathbf{s}_{2\ell+1}\rangle^{2}=\left\|\mathbf{F}^{(\mathrm{S})}\frac{1}{\sqrt{M}}\mathbf{1}\right\|_{2}^{2}=1.$

∎

Let $\widehat{\mathbf{S}}$ be
$\widehat{\mathbf{S}}=\mathbf{S}\mathbf{S}^{\top}=\begin{bmatrix}\widehat{\mathbf{s}}_{0}&\widehat{\mathbf{s}}_{1}&\cdots&\widehat{\mathbf{s}}_{M-1}\end{bmatrix}^{\top}$.
From Lemma 1 2),
$\widehat{\mathbf{s}}_{0}-\sum^{M/2-1}_{\ell=0}\langle\mathbf{s}_{0},\mathbf{s}_{2\ell+1}\rangle\widehat{\mathbf{s}}_{2\ell+1}=\mathbf{0}.$
This implies that $\mathrm{rank}(\widehat{\mathbf{S}})=M-1$, and so
$\mathrm{rank}(\mathbf{S})=M-1$. ∎

##  Appendix B Proof for Proposition 3

###### Proof:

From Proposition 1, the elements in the upper-right triangle of
$\widehat{\mathbf{S}}^{(0)}=\mathbf{S}^{(0)}\mathbf{S}^{(0)\top}=\begin{bmatrix}\widehat{\mathbf{s}}_{0}^{(0)}&\cdots&\widehat{\mathbf{s}}_{M-1}^{(0)}\end{bmatrix}^{\top}$
are as follows:

$\displaystyle[\widehat{\mathbf{S}}^{(0)}]_{k_{v},k_{h}}=\begin{cases}\frac{\sqrt{2}}{M\sin\left(\frac{\pi}{2M}k_{h}\right)}&(k_{v}=0\
\mathrm{and}\ k_{h}=2\ell+1)\\\ 1&(k_{v}=k_{h}\ \mathrm{and}\ k_{v}\neq 1)\\\
0&(\mathrm{otherwise})\end{cases},$

where $1\leq\ell\leq\frac{M}{2}-1$. For example, for $M=4$,

$\displaystyle\widehat{\mathbf{S}}^{(0)}\approx\begin{bmatrix}1&0&0&0.3827\\\
0&0&0&0\\\ 0&0&1&0\\\ 0.3827&0&0&1\end{bmatrix}.$ (19)

Then, we can derive that

$\displaystyle\widehat{\mathbf{s}}_{0}^{(0)}-\sum^{M/2-1}_{\ell=1}\langle\mathbf{s}_{0},\mathbf{s}_{2\ell+1}\rangle\widehat{\mathbf{s}}_{2\ell+1}^{(0)}\neq\mathbf{0}.$
(20)

Thus, it can be concluded that $\mathrm{rank}(\mathbf{S}^{(0)})=M-1$. ∎

##  Appendix C Proof for Proposition 4

###### Proof:

The statements 1) and 2) are clearly true. We only show the proof for 3).

For any even row $2k\geq 2$,
$\mathcal{F}[[\mathbf{F}^{(\mathrm{RS})}]_{2k,\cdot}]$ is exactly the same as
$\mathcal{F}[[\mathbf{F}^{(\mathrm{S})}]_{2k,\cdot}]$. Thus, it is enough to
show the case of odd rows $[\mathbf{F}^{(\mathrm{RS})}]_{2k+1,\cdot}$
$(2k+1\geq 3)$.

$[\mathbf{F}^{(\mathrm{RS})}]_{2k+1,n}$ is the same as
$[\mathbf{s}^{(k+1)}_{2k+1}]_{n}$ of
$\mathbf{S}^{(k+1)}=\begin{bmatrix}\mathbf{s}_{0}^{(k+1)}&\ldots&\mathbf{s}_{M-1}^{(k+1)}\end{bmatrix}^{\top}$
in Step 5 of Algorithm 1. Let
$\mathbf{T}^{(k)}=\begin{bmatrix}\mathbf{t}_{0}^{(k)}&\ldots&\mathbf{t}_{M-1}^{(k)}\end{bmatrix}$
be the inverse matrix of $\mathbf{S}^{(k)}$. Since $\mathbf{s}^{(k+1)}_{2k+1}$
is designed to be orthogonal to
$\\{\mathbf{s}^{(k)}_{n}\\}_{\Omega_{M-1}\setminus\\{2k+1\\}}$ in Algorithm 1,
$\mathbf{s}^{(k+1)}_{2k+1}$ can be expressed with a linear combination of
$\\{\mathbf{t}^{(k)}_{n}\\}_{\Omega_{M-1}}$ as:

$\displaystyle\mathbf{s}^{(k+1)}_{2k+1}=$
$\displaystyle\mathbf{T}^{(k)}\mathbf{S}^{(k)}\mathbf{s}^{(k+1)}_{2k+1}=\sum_{m=0}^{M-1}\langle\mathbf{s}^{(k)}_{m},\mathbf{s}^{(k+1)}_{2k+1}\rangle\mathbf{t}^{(k)}_{m}$
$\displaystyle=$
$\displaystyle\langle\mathbf{s}^{(k)}_{2k+1},\mathbf{s}^{(k+1)}_{2k+1}\rangle\mathbf{t}^{(k)}_{2k+1}.$
(21)

Here, we use the following lemma (see its proof in Appendix D).

###### Lemma 2

The passband of the frequency response of $\mathbf{t}^{(k)}_{\ell}$ of
$\mathbf{T}^{(k)}=\begin{bmatrix}\mathbf{t}_{0}^{(k)}&\ldots&\mathbf{t}_{M-1}^{(k)}\end{bmatrix}$
is the same as that of $\mathbf{s}^{(k)}_{\ell}$ of
$\mathbf{S}^{(k)}=\begin{bmatrix}\mathbf{s}_{0}^{(k)}&\ldots&\mathbf{s}_{M-1}^{(k)}\end{bmatrix}^{\top}$.

From Lemma 2, the passband of the frequency response of
$\mathbf{s}^{(k+1)}_{2k+1}$ is located at the same position as that of
$\mathbf{t}^{(k)}_{2k+1}$ and
$\mathbf{s}^{(k)}_{2k+1}=([\mathbf{F}^{(\mathrm{S})}]_{2k+1,\cdot})^{\top}$.
Consequently, we conclude that statement 3) is true. ∎

##  Appendix D Proof for Lemma 2

###### Proof:

First, consider the case of $k=1$. Since
$\mathbf{T}^{(1)}=\begin{bmatrix}\mathbf{t}_{0}^{(1)}&\ldots&\mathbf{t}_{M-1}^{(1)}\end{bmatrix}$
is the inverse of
$\mathbf{S}^{(1)}=\begin{bmatrix}\mathbf{s}_{0}^{(1)}&\ldots&\mathbf{s}_{M-1}^{(1)}\end{bmatrix}^{\top}$,
each $\mathbf{t}_{m}^{(1)}$ can be represented as
$\mathbf{t}_{m}^{(1)}=\sum_{n=0}^{M-1}\langle\mathbf{t}_{n}^{(1)},\mathbf{t}_{m}^{(1)}\rangle\mathbf{s}_{n}^{(1)}$.
Therefore, it is enough to show that
$|[(\mathbf{T}^{(1)})^{\top}\mathbf{T}^{(1)}]_{m,n}|=|\langle\mathbf{t}_{m}^{(1)},\mathbf{t}_{n}^{(1)}\rangle|\ll|\langle\mathbf{t}_{m}^{(1)},\mathbf{t}_{m}^{(1)}\rangle|$.

For that, we consider the eigenvalue decomposition of
$\mathbf{S}^{(1)}\mathbf{S}^{(1)\top}=\mathbf{U}^{(1)}\mathbf{D}^{(1)}\mathbf{U}^{(1)\top}$,
then calculate
$\mathbf{T}^{(1)\top}\mathbf{T}^{(1)}=\mathbf{U}^{(1)}(\mathbf{D}^{(1)})^{-1}\mathbf{U}^{(1)\top}$,
where
$\mathbf{U}^{(1)}=\begin{bmatrix}\mathbf{u}_{0}&\ldots&\mathbf{u}_{M-1}\end{bmatrix}$
and $\mathbf{D}^{(1)}=\mathrm{diag}(\lambda_{0},\ldots,\lambda_{M-1})$ are
some orthogonal and diagonal matrices consisting of eigenvectors and
eigenvalues, respectively. Similar to Lemma 1, we can derive
$\widehat{\mathbf{S}}^{(1)}=\mathbf{S}^{(1)}\mathbf{S}^{(1)\top}$ forms

$\displaystyle\widehat{\mathbf{S}}^{(1)}=$ $\displaystyle\
\begin{bmatrix}1&0&0&\widehat{s}_{3}&0&\widehat{s}_{5}&\cdots\\\
0&1&0&0&0&0&\cdots\\\ 0&0&1&0&0&0&\cdots\\\
\widehat{s}_{3}&0&0&1&0&0&\cdots\\\ 0&0&0&0&1&0&\cdots\\\
\widehat{s}_{5}&0&0&0&0&1&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\\ \end{bmatrix},$ (22)

where
$\widehat{s}_{2\ell+1}=\frac{\sqrt{2}}{M\sin\left(\frac{\pi}{2M}(2\ell+1)\right)}\
(\ell\geq 1)$. Let us consider some eigenvalue $\lambda$ and its corresponding
eigenvectors $\mathbf{u}$ of $\widehat{\mathbf{S}}^{(1)}$. Note that all the
eigenvalues are positive $\lambda_{n}>0$, since
$\mathrm{rank}(\widehat{\mathbf{S}}^{(1)})=M$. Suppose an eigenvalue
$\lambda=1$, then its eigenvector
$\mathbf{u}=\begin{bmatrix}{u}_{0}&\ldots&{u}_{M-1}\end{bmatrix}^{\top}$
should satisfy

$\displaystyle\widehat{\mathbf{S}}^{(1)}\begin{bmatrix}{u}_{0}&\ldots&{u}_{M-1}\end{bmatrix}^{\top}=\begin{bmatrix}{u}_{0}&\ldots&{u}_{M-1}\end{bmatrix}^{\top}$
$\displaystyle\Rightarrow$
$\displaystyle\begin{cases}{u}_{3}=-\sum_{\ell=2}^{M/2-1}\frac{\widehat{s}_{2\ell+1}}{\widehat{s}_{3}}{u}_{2\ell+1}\\\
{u}_{0}=0\end{cases}$ $\displaystyle\Rightarrow$
$\displaystyle\mathbf{u}=\begin{bmatrix}0&u_{1}&u_{2}&-\sum_{\ell=2}^{M/2-1}\frac{\widehat{s}_{2\ell+1}}{\widehat{s}_{3}}{u}_{2\ell+1}&{u}_{4}&\cdots\end{bmatrix}^{\top}$
$\displaystyle\Rightarrow$
$\displaystyle\mathbf{u}\in\mathrm{span}\left\\{\mathbf{u}_{1},\ldots,\mathbf{u}_{M-2}\right\\},$
(23)

where

$\displaystyle\mathbf{u}_{m}=$ $\displaystyle\
\begin{cases}{\bm{\delta}_{m}}&(m=1,2)\\\ {\bm{\delta}_{m+1}}&(m=2\ell+1,m\geq
3)\\\
{\bm{\delta}_{m+1}}-\frac{\widehat{s}_{m+1}}{\widehat{s}_{3}}{\bm{\delta}_{3}}&(m=2\ell,m\geq
3)\\\ \end{cases},$ (24)

where ${\bm{\delta}_{m}}\in\mathbb{R}^{M}$ $(m\in\Omega_{M-1})$ consists of
$[{\bm{\delta}}_{m}]_{m}=1$ and $[{\bm{\delta}}_{m}]_{n}=0$ $(m\neq n)$. Since
$\mathbf{U}^{(1)}$ should be an orthogonal matrix, but the vectors
$\\{\mathbf{u}_{m}\\}_{m=2\ell,m\geq 3}$ are not orthogonal yet, Gram-Schmidt
orthonormalization is applied to them.

Next, we consider the case of $\lambda\neq 1$.

$\displaystyle\widehat{\mathbf{S}}^{(1)}\begin{bmatrix}{u}_{0}&\ldots&{u}_{M-1}\end{bmatrix}^{\top}=\lambda\begin{bmatrix}{u}_{0}&\ldots&{u}_{M-1}\end{bmatrix}^{\top}$
$\displaystyle\Rightarrow$
$\displaystyle\begin{cases}\sum_{\ell=1}^{M/2-1}\widehat{s}_{2\ell+1}{u}_{2\ell+1}=(\lambda-1){u}_{0}\
(u_{0}\neq 0)\\\ \widehat{s}_{2\ell+1}{u}_{0}=(\lambda-1){u}_{2\ell+1}\\\
(\lambda-1){u}_{\ell}=0\end{cases}$ $\displaystyle\Rightarrow$
$\displaystyle\lambda^{2}-2\lambda+\left(1-\sum_{\ell=1}^{M/2-1}\widehat{s}_{2\ell+1}^{2}\right)=0$
(25)

Thus, $\lambda_{0}=1+\sqrt{\sum_{\ell=1}^{M/2-1}\widehat{s}_{2\ell+1}^{2}},\
\lambda_{M-1}=1-\sqrt{\sum_{\ell=1}^{M/2-1}\widehat{s}_{2\ell+1}^{2}}$. For
$\lambda_{0}$ and $\lambda_{M-1}$, the eigenvectors $\mathbf{u}_{0}$ and
$\mathbf{u}_{M-1}$ can be found as

$\displaystyle\mathbf{u}_{0}^{(1)}=$ $\displaystyle\
\frac{1}{\sqrt{2}}{\bm{\delta}}_{0}+\sum_{\ell=1}^{M/2-1}\frac{\widehat{s}_{2\ell+1}}{\sqrt{2}\sqrt{\sum_{\ell=1}^{M/2-1}\widehat{s}_{2\ell+1}^{2}}}{\bm{\delta}}_{2\ell+1},$
$\displaystyle\mathbf{u}_{M-1}^{(1)}=$ $\displaystyle\
\frac{1}{\sqrt{2}}{\bm{\delta}}_{0}-\sum_{\ell=1}^{M/2-1}\frac{\widehat{s}_{2\ell+1}}{\sqrt{2}\sqrt{\sum_{\ell=1}^{M/2-1}\widehat{s}_{2\ell+1}^{2}}}{\bm{\delta}}_{2\ell+1}.$
(26)

$\\{\lambda_{n}\\}$ and $\\{\mathbf{u}^{(1)}_{n}\\}$ give us the eigenvalue
decomposition of $\widehat{\mathbf{S}}^{(1)}$. For example, when $M=8$,
${\mathbf{D}}^{(1)}=\mathrm{diag}(\lambda_{0},1,1,1,1,1,1,\lambda_{M-1})$ and

$\displaystyle\mathbf{U}^{(1)}=$
$\displaystyle\begin{bmatrix}\frac{1}{\sqrt{2}}&0&0&0&0&0&0&\frac{1}{\sqrt{2}}\\\
0&1&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0\\\
\frac{A_{3}}{\sqrt{2}}&0&0&0&u_{3,4}&0&u_{3,6}&-\frac{A_{3}}{\sqrt{2}}\\\
0&0&0&1&0&0&0&0\\\
\frac{A_{5}}{\sqrt{2}}&0&0&0&u_{5,4}&0&0&-\frac{A_{5}}{\sqrt{2}}\\\
0&0&0&0&0&1&0&0\\\
\frac{A_{7}}{\sqrt{2}}&0&0&0&0&0&u_{7,6}&-\frac{A_{7}}{\sqrt{2}}\end{bmatrix},$
(27)

where
$A_{2\ell+1}=\frac{\widehat{s}_{2\ell+1}}{\sqrt{\widehat{s}_{3}^{2}+\widehat{s}_{5}^{2}+\widehat{s}_{7}^{2}}}\
(<1)$ and $u_{3,4},\ u_{3,6},\ u_{5,4},\ u_{7,6}$ are the elements after
orthogonalization. Then, the elements in the upper-right triangle of
$\widehat{\mathbf{S}}^{(1)}$ are

$\displaystyle[\widehat{\mathbf{S}}^{(1)}]_{m,n}=[{\mathbf{U}}^{(1)}{\mathbf{D}}^{(1)}{\mathbf{U}}^{(1)\top}]_{m,n}$
$\displaystyle=$ $\displaystyle\
\begin{cases}\frac{1}{2}(\lambda_{0}+\lambda_{M-1})\quad(m=n=1)\\\
\frac{A_{n}}{2}(\lambda_{0}-\lambda_{M-1})\quad(m=0,n=2\ell+1\geq 3)\\\
\frac{A_{m}^{2}}{2}(\lambda_{0}+\lambda_{M-1})+\Delta_{m,m}\quad(m=n=2\ell+1\geq
3)\\\ \frac{A_{m}A_{n}}{2}(\lambda_{0}+\lambda_{M-1})+\Delta_{m,n}\\\
(m=2\ell_{m}+1,\ n=2\ell_{n}+1,\ \ell_{m}\neq\ell_{n},\ m,\ n\geq 3)\\\
1\quad(m=n=1,\ 2,\ 2\ell\ (\ell\geq 2))\\\
0\quad\mathrm{otherwise}\end{cases},$ (28)

where $\Delta_{m,n}$ contains the result of multiplication. The upper-right
triangle elements of
$\widehat{\mathbf{T}}^{(1)}={\mathbf{U}}^{(1)}({\mathbf{D}}^{(1)})^{-1}{\mathbf{U}}^{(1)\top}$
are

$\displaystyle[\widehat{\mathbf{T}}^{(1)}]_{m,n}=[{\mathbf{U}}^{(1)}({\mathbf{D}}^{(1)})^{-1}{\mathbf{U}}^{(1)\top}]_{m,n}$
$\displaystyle=$ $\displaystyle\
\begin{cases}\frac{1}{2}(\frac{1}{\lambda_{0}}+\frac{1}{\lambda_{M-1}})\quad(m=n=1)\\\
\frac{A_{n}}{2}(\frac{1}{\lambda_{0}}-\frac{1}{\lambda_{M-1}})\quad(m=0,n=2\ell+1\geq
3)\\\
\frac{A_{m}^{2}}{2}(\frac{1}{\lambda_{0}}+\frac{1}{\lambda_{M-1}})+\Delta_{m,m}\quad(m=n=2\ell+1\geq
3)\\\
\frac{A_{m}A_{n}}{2}(\frac{1}{\lambda_{0}}+\frac{1}{\lambda_{M-1}})+\Delta_{m,n}\\\
(m=2\ell_{m}+1,\ n=2\ell_{n}+1,\ \ell_{m}\neq\ell_{n},\ m,\ n\geq 3)\\\
1\quad(m=n=1,\ 2,\ 2\ell\ (\ell\geq 2))\\\
0\quad\mathrm{otherwise}\end{cases}.$ (29)

From Lemma 1, it follows that

$\displaystyle\frac{1}{\lambda_{0}\lambda_{M-1}}=\frac{1}{1-(\sum_{\ell=1}^{M/2-1}\widehat{s}_{2\ell+1}^{2})}=\frac{1}{\widehat{s}_{1}^{2}}=\frac{(M\sin(\frac{\pi}{2M}))^{2}}{2}.$
(30)

Since
$M\sin(\frac{\pi}{2M})=\frac{\pi}{2}\frac{2M}{\pi}\sin(\frac{1}{\frac{2M}{\pi}})$,
$x\sin(\frac{1}{x})$ monotonically increases over $[\frac{4}{\pi},\infty)$ and
$x\sin(\frac{1}{x})\xrightarrow{x\rightarrow\infty}1$, then

$\displaystyle 1.1716\approx\frac{(4\sin(\frac{\pi}{8}))^{2}}{2}<$
$\displaystyle\frac{(M\sin(\frac{\pi}{2M}))^{2}}{2}\xrightarrow{M\rightarrow\infty}\frac{\pi^{2}}{8}\approx
1.2337,$ $\displaystyle\frac{1}{\lambda_{0}\lambda_{M-1}}=$ $\displaystyle\
1+\epsilon\quad(\epsilon<\frac{1}{4}),$ (31)

and
$\frac{1}{\lambda_{0}}+\frac{1}{\lambda_{M-1}}=\frac{\lambda_{0}+\lambda_{M-1}}{\lambda_{0}\lambda_{M-1}}=(\lambda_{0}+\lambda_{M-1})(1+\epsilon)$,
$\frac{1}{\lambda_{0}}-\frac{1}{\lambda_{M-1}}=-\frac{\lambda_{0}-\lambda_{M-1}}{\lambda_{0}\lambda_{M-1}}=-(\lambda_{0}-\lambda_{M-1})(1+\epsilon)$.
Thus, by substituting (22), (D), and (D) into (D), we can derive

$\displaystyle[\widehat{\mathbf{T}}^{(1)}]_{m,n}=[{\mathbf{U}}^{(1)}({\mathbf{D}}^{(1)})^{-1}{\mathbf{U}}^{(1)\top}]_{m,n}$
$\displaystyle=$ $\displaystyle\ \begin{cases}1+\epsilon>1\quad(m=n=1)\\\
-\widehat{s}_{n}(1+\epsilon)\quad(m=0,n=2\ell+1\geq 3)\\\
1+\frac{A_{m}^{2}}{2}\epsilon>1\quad(m=n=2\ell+1\geq 3)\\\
\frac{A_{m}A_{n}}{2}\epsilon<\frac{1}{8}\\\ (m=2\ell_{m}+1,\ n=2\ell_{n}+1,\
\ell_{m}\neq\ell_{n},\ m,\ n\geq 3)\\\ 1\quad(m=n=1,\ 2,\ 2\ell\ (\ell\geq
2))\\\ 0\quad\mathrm{otherwise}\end{cases}.$ (32)

Here, let
$\rho(M,n)=\widehat{s}_{n}=\frac{\sqrt{2}}{(M\sin(\frac{\pi}{2M}n))}$. Since
we assume that the size $M$ for the RDADCF is $M\geq 4$,

$\displaystyle\frac{2}{5}>0.3827\approx\rho(4,3)>\rho(M,n)>\rho(M+1,n),$
$\displaystyle\frac{2}{5}>0.3827\approx\rho(4,3)>\rho(M,n)>\rho(M,n+1),$ (33)

thus $|-\widehat{s}_{n}(1+\epsilon)|<\frac{2}{5}\frac{5}{4}=\frac{1}{2}$.
Finally, we conclude that
$|[\mathbf{T}^{(1)\top}\mathbf{T}^{(1)}]_{m,n}|=|\langle\mathbf{t}_{m}^{(1)},\mathbf{t}_{n}^{(1)}\rangle|\ll|\langle\mathbf{t}_{m}^{(1)},\mathbf{t}_{m}^{(1)}\rangle|$,
which implies the passband of each $\mathbf{t}_{m}^{(1)}$ is the same as
$\mathbf{s}_{m}^{(1)}$.

For $k=2$, $\mathbf{S}^{(2)}$ forms as in (22) with $\widehat{s}_{3}=0$. With
the same discussion when $k=1$, it can be derived that lower bounds of the
diagonal elements of $|[\widehat{\mathbf{T}}^{(2)}]_{m,m}|$ are 1 and upper
bounds of the elements $|[\widehat{\mathbf{T}}^{(2)}]_{m,n}|$ ($m\neq n$) are
$\frac{1}{2}$ or $\frac{1}{8}$, as in (D). Thus,
$|\langle\mathbf{t}_{m}^{(2)},\mathbf{t}_{n}^{(2)}\rangle|\ll|\langle\mathbf{t}_{m}^{(2)},\mathbf{t}_{m}^{(2)}\rangle|$.
This is the end of proof for Lemma 2. ∎

##  Appendix E Detailed Algorithm of Image Recovery Used in Experiments

To solve (VI), the primal-dual splitting (PDS) algorithm [41, 42] is used.
Consider the following convex optimization problem to find

$\displaystyle\mathbf{x}^{\star}\in\mathop{\rm
arg~{}min}\limits_{\mathbf{x}\in\mathbb{R}^{N_{1}}}f(\mathbf{x})+g(\mathbf{L}\mathbf{x}),$
(34)

where $f\in\Gamma_{0}(\mathbb{R}^{N_{1}})$,
$g\in\Gamma_{0}(\mathbb{R}^{N_{2}})$ ($\Gamma_{0}(\mathbb{R}^{N_{2}})$ is the
set of proper lower semicontinuous convex functions [43] on $\mathbb{R}^{N}$),
and $\mathbf{L}\in\mathbb{R}^{{N_{2}}\times{N_{1}}}$. Then, the optimal
solution $\mathbf{x}^{\star}$, can be obtained as

$\displaystyle\begin{cases}\mathbf{x}^{(n+1)}:=\mathrm{prox}_{\gamma_{1}f}[\mathbf{x}^{(n)}-\gamma_{1}\mathbf{L}^{\top}\mathbf{z}^{(n)}]\\\
\mathbf{z}^{(n+1)}:=\mathrm{prox}_{\gamma_{2}g^{\ast}}[\mathbf{z}^{(n)}+\gamma_{2}\mathbf{L}(2\mathbf{x}^{(n+1)}-\mathbf{x}^{(n)})]\end{cases},$
(35)

where $\mathrm{prox}$ denotes the proximal operator [43], $g^{\ast}$ is the
conjugate function [43] of $g$. In the experiments, the parameters
$\gamma_{1}$ and $\gamma_{2}$ in (35), are chosen as 0.01 and
$\frac{1}{12\gamma_{1}}$. For Problem 2 ($\rho=1$ in (VI)), the functions $f$
and $g$, and the matrix $\mathbf{L}$ in (34), are set as

$\displaystyle f(\mathbf{x})=\iota_{C_{[0,1]}}(\mathbf{x}),$ $\displaystyle
g([\mathbf{z}_{1}^{\top}\ \mathbf{z}_{2}^{\top}\
\mathbf{z}_{3}^{\top}]^{\top})=\|\mathbf{z}_{1}\|_{1}+\|\mathbf{z}_{2}\|_{1,2}+\iota_{\\{\mathbf{y}\\}}(\mathbf{z}_{3}),$
$\displaystyle\mathbf{z}_{1}=\mathbf{F}\mathbf{P}_{\mathrm{v2bv}}\mathbf{x},\
\mathbf{z}_{2}=\widetilde{\mathbf{W}}_{\mathrm{b}}\mathbf{D}_{\mathrm{hv}}\mathbf{x},\
\mathbf{z}_{3}=\widetilde{\mathbf{\Phi}}{\mathbf{x}},$
$\displaystyle\mathbf{L}=\begin{bmatrix}(\mathbf{F}\mathbf{P}_{\mathrm{v2bv}})^{\top}&(\widetilde{\mathbf{W}}_{\mathrm{b}}\mathbf{D}_{\mathrm{hv}})^{\top}&\widetilde{\mathbf{\Phi}}^{\top}\end{bmatrix}^{\top}.$
(36)

The resulting solver for (VI) is described in Algorithm 2666For
$\mathbf{x}\in\mathbb{R}^{N}$,
$[\mathrm{prox}_{\gamma\|\cdot\|_{1}}(\mathbf{x})]_{i}=\mathrm{sign}(x_{i})\max\\{|x_{i}|-\gamma,0\\}$
(soft-thresholding), $\mathrm{prox}_{\gamma\|\cdot\|_{1,2}}(\mathbf{x})$ is
the group soft-thresholding [44],
$\mathrm{prox}_{\iota_{C_{[0,1]}}}(\mathbf{x})$ is the clipping operation to
$[0,1]$, and
$\mathrm{prox}_{\iota_{\\{\mathbf{v}\\}}}(\mathbf{x})=\mathbf{y}$, where
$\mathbf{y}\in\mathbb{R}^{N}$ is an observation. . The stopping criterion is
$\|\mathbf{x}^{(n+1)}-\mathbf{x}^{(n)}\|_{2}\leq 0.01$. The algorithm fof
Problem 1 ($\rho=0$ in (VI)) can be designed by removing the terms and steps
(Step 5, 8, and 10) relating to $\mathbf{z}_{2}$, $\mathbf{t}_{2}$,
${\mathbf{t}}_{2}^{(n)}$, and $\hat{\mathbf{t}}_{2}^{(n)}$ from Algorithm 2.

Algorithm 2 Solver for (VI)

1: set $n=0$ and choose $\mathbf{x}^{(0)}$, $\mathbf{z}_{1}^{(0)}$,
$\mathbf{z}_{2}^{(0)}$, $\gamma_{1}$, $\gamma_{2}$.

2: while stop criterion is not satisfied do

3:
$\mathbf{x}^{(n+1)}=\mathrm{prox}_{\gamma_{1}\iota_{C_{[0,1]}}}(\mathbf{x}^{(n)}-\gamma_{1}((\mathbf{F}\mathbf{P}_{\mathrm{v2bv}})^{\top}\mathbf{z}_{1}^{(n)}+(\widetilde{\mathbf{W}}_{\mathrm{b}}\mathbf{D}_{\mathrm{hv}})^{\top}\mathbf{z}_{2}^{(n)}+\widetilde{\mathbf{\Phi}}^{\top}\mathbf{z}_{3}^{(n)}))$

4:
$\mathbf{t}_{1}^{(n)}=\mathbf{z}_{1}^{(n)}+\gamma_{2}\mathbf{F}\mathbf{P}_{\mathrm{v2bv}}(2\mathbf{x}^{(n+1)}-\mathbf{x}^{(n)})$.

5:
$\mathbf{t}_{2}^{(n)}=\mathbf{z}_{2}^{(n)}+\gamma_{2}\widetilde{\mathbf{W}}_{\mathrm{b}}\mathbf{D}_{\mathrm{hv}}(2\mathbf{x}^{(n+1)}-\mathbf{x}^{(n)})$.

6:
$\mathbf{t}_{3}^{(n)}=\mathbf{z}_{3}^{(n)}+\gamma_{2}\widetilde{\mathbf{\Phi}}(2\mathbf{x}^{(n+1)}-\mathbf{x}^{(n)})$.

7:
$\hat{\mathbf{t}}_{1}^{(n)}=\mathrm{prox}_{\frac{1}{\gamma_{2}}\|\cdot\|_{1}}\left(\frac{1}{\gamma_{2}}\mathbf{t}_{1}^{(n)}\right)$.

8:
$\hat{\mathbf{t}}_{2}^{(n)}=\mathrm{prox}_{\frac{1}{\gamma_{2}}\|\cdot\|_{1,2}}\left(\frac{1}{\gamma_{2}}\mathbf{t}_{2}^{(n)}\right)$.

9:
$\hat{\mathbf{t}}_{3}^{(n)}=\mathrm{prox}_{\frac{1}{\gamma_{2}}\iota_{\\{\mathbf{y}\\}}}\left(\frac{1}{\gamma_{2}}\mathbf{t}_{3}^{(n)}\right)$.

10:
$\mathbf{z}_{k}^{(n+1)}=\mathbf{t}_{k}^{(n)}-\gamma_{2}\hat{\mathbf{t}}_{k}^{(n)}$
($k=1,2,3$).

11: $n=n+1$.

12: end while

13: Output $\mathbf{u}^{(n)}$.

## References

  * [1] M. Elad, _Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing_ , 1st ed. Springer Publishing Company, Incorporated, 2010.
  * [2] J.-L. Starck, F. Murtagh, and J. Fadili, _Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity_. New York, NY, USA: Cambridge University Press, 2010.
  * [3] P. L. Combettes and J.-C. Pesquet, “A proximal decomposition method for solving convex variational inverse problems,” _Inverse Problems_ , vol. 24, no. 6, p. 065014, Nov. 2008.
  * [4] M. V. Afonso, J. M. B.-Dias, and M. A. T. Figueiredo, “An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems,” _IEEE Trans. Image Process._ , vol. 20, no. 3, pp. 681–695, Mar. 2011.
  * [5] J. L. Starck, E. J. Candès, and D. L. Donoho, “The curvelet transform for image denoising,” _IEEE Trans. Image Process._ , vol. 11, no. 6, pp. 670–684, Jun. 2002.
  * [6] M. N. Do and M. Vetterli, “The contourlet transform: an efficient directional multiresolution image representation,” _IEEE Trans. Image Process._ , vol. 14, no. 12, pp. 2091–2106, Dec. 2005.
  * [7] G. Easley, D. Labate, and W.-Q. Lim, “Sparse directional image representations using the discrete shearlet transform,” _Applied and Computational Harmonic Analysis_ , vol. 25, no. 1, pp. 25 – 46, 2008.
  * [8] Y. Tanaka, M. Ikehara, and T. Q. Nguyen, “Multiresolution image representation using combined 2-D and 1-D directional filter banks,” _IEEE Trans. Image Process._ , vol. 18, no. 2, pp. 269–280, Feb. 2009.
  * [9] W. Q. Lim, “Nonseparable shearlet transform,” _IEEE Trans. Image Process._ , vol. 22, no. 5, pp. 2056–2065, May 2013.
  * [10] C. Zhiyu and S. Muramatsu, “Poisson denoising with multiple directional LOTs,” in _Proc. IEEE Int. Conf. Acoust., Speech and Signal Process. (ICASSP)_ , May 2014.
  * [11] Y. Shi, X. Yang, and Y. Guo, “Translation invariant directional framelet transform combined with gabor filters for image denoising,” _IEEE Trans. Image Process._ , vol. 23, no. 1, pp. 44–55, Jan. 2014.
  * [12] I. W. Selesnick, R. G. Baraniuk, and N. C. Kingsbury, “The dual-tree complex wavelet transform,” _IEEE Signal Process. Mag._ , vol. 22, no. 6, pp. 123–151, Nov. 2005.
  * [13] L. Liang and H. Liu, “Dual-tree cosine-modulated filter bank with linear-phase individual filters: An alternative shift-invariant and directional-selective transform,” _IEEE Trans. Image Process._ , vol. 22, no. 12, pp. 5168–5180, Dec. 2013.
  * [14] S. Kyochi and M. Ikehara, “A class of near shift-invariant and orientation-selective transform based on delay-less oversampled even-stacked cosine-modulated filter banks,” _IEICE Trans. Fundam._ , vol. 93, no. 4, pp. 724–733, Apr. 2010.
  * [15] 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.
  * [16] R. Ishibashi, T. Suzuki, S. Kyochi, and H. Kudo, “Image boundary extension with mean value for cosine-sine modulated lapped/block transforms,” _IEEE Trans. Circuits and Syst. for Video Technol._ , vol. 29, no. 1, pp. 1–11, Jan. 2019.
  * [17] S. Mallat, _A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way_ , 3rd ed. Academic Press, 2008.
  * [18] Y. M. Lu and M. N. Do, “Multidimensional directional filter banks and surfacelets,” _IEEE Trans. Image Process._ , vol. 16, no. 4, pp. 918–931, Apr. 2007.
  * [19] J. Yang, Y. Wang, W. Xu, and Q. Dai, “Image and video denoising using adaptive dual-tree discrete wavelet packets,” _IEEE Trans. Circuits and Syst. for Video Technol._ , vol. 19, no. 5, pp. 642–655, May 2009.
  * [20] T. T. Nguyen and H. Chauris, “Uniform discrete curvelet transform,” _IEEE Trans. Signal Process._ , vol. 58, no. 7, pp. 3618–3634, Jul. 2010\.
  * [21] S. Held, M. Storath, P. Massopust, and B. Forster, “Steerable wavelet frames based on the riesz transform,” _IEEE Trans. Image Process._ , vol. 19, no. 3, pp. 653–667, Mar. 2010.
  * [22] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation,” _IEEE Trans. Signal Process._ , vol. 54, no. 11, pp. 4311–4322, Nov. 2006.
  * [23] R. Rubinstein, T. Peleg, and M. Elad, “Analysis K-SVD: A dictionary-learning algorithm for the analysis sparse model,” _IEEE Trans. Signal Process._ , vol. 61, no. 3, pp. 661–677, Feb. 2013.
  * [24] Y. Tanaka and A. Sakiyama, “$M$-channel oversampled graph filter banks,” _IEEE Trans. Signal Process._ , vol. 62, no. 14, pp. 3578–3590, Jul. 2014\.
  * [25] A. Danielyan, V. Katkovnik, and K. Egiazarian, “BM3D frames and variational image deblurring,” _IEEE Trans. Image Process._ , vol. 21, no. 4, pp. 1715–1728, Apr. 2012.
  * [26] S. Fujita, N. Fukushima, M. Kimura, and Y. Ishibashi, “Randomized redundant DCT: Efficient denoising by using random subsampling of dct patches,” in _Proc. SIGGRAPH Asia 2015 Technical Briefs_ , 2015.
  * [27] Y. Wang, C. Xu, S. You, C. Xu, and D. Tao, “DCT regularized extreme visual recovery,” _IEEE Trans. Image Process._ , vol. 26, no. 7, pp. 3360–3371, July 2017.
  * [28] K. R. Rao and P. Yip, _Discrete Cosine Transform: Algorithms, Advantages, Applications_. San Diego, CA, USA: Academic Press Professional, Inc., 1990.
  * [29] S. Martucci, “Symmetric convolution and the discrete sine and cosine transforms,” _IEEE Trans. Signal Process._ , vol. 42, no. 5, pp. 1038–1051, 1994.
  * [30] 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.
  * [31] R. N. Bracewell, “The fast Hartley transform,” _Proc. of the IEEE_ , vol. 72, no. 8, pp. 1010–1018, Aug. 1984.
  * [32] S. Kyochi and Y. Tanaka, “General factorization of conjugate-symmetric hadamard transforms,” _IEEE Trans. Signal Process._ , vol. 62, no. 13, pp. 3379–3392, Jul. 2014.
  * [33] H. Xu, J. Xu, and F. Wu, “Lifting-based directional DCT-like transform for image coding,” _IEEE Trans. Circuits and Syst. for Video Technol._ , vol. 17, no. 10, pp. 1325–1335, Oct. 2007.
  * [34] B. Zeng and J. Fu, “Directional discrete cosine transforms; a new framework for image coding,” _IEEE Trans. Circuits and Syst. for Video Technol._ , vol. 18, no. 3, pp. 305–313, Mar. 2008.
  * [35] G. J. Sullivan, J. R. Ohm, W. J. Han, and T. Wiegand, “Overview of the high efficiency video coding (HEVC) standard,” _IEEE Trans. Circuits and Syst. for Video Technol._ , vol. 22, no. 12, pp. 1649–1668, Dec. 2012.
  * [36] B. Bross, Y.-K. Wang, Y. Ye, S. Liu, J. Chen, G. J. Sullivan, and J.-R. Ohm, “Overview of the versatile video coding (VVC) standard and its applications,” _IEEE Trans. Circuits and Syst. for Video Technol._ , vol. 31, no. 10, pp. 3736–3764, 2021.
  * [37] T. Ichita, S. Kyochi, T. Suzuki and Y. Tanaka, “Directional discrete cosine transforms arising from discrete cosine and sine transform for directional block-wise image representation,” in _Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)_ , 2017.
  * [38] G. Strang and T. Nguyen, _Wavelets and Filter Banks_. Wellesley-Cambridge Press, 1996.
  * [39] R. Coifman, F. Geshwind, and Y. Meyer, “Noiselets,” _Applied and Computational Harmonic Analysis_ , vol. 10, no. 1, pp. 27 – 44, 2001.
  * [40] F. Alter, S. Y. Durand, and J. Froment, “Deblocking DCT-based compressed images with weighted total variation,” in _Proc. IEEE Int. Conf. Acoust., Speech, and Signal Process._ , 2004.
  * [41] L. Condat, “A primal-dual splitting method for convex optimization involving lipschitzian, proximable and linear composite terms,” _J. Optimization Theory and Applications_ , Dec. 2013.
  * [42] B. C. Vu, “A splitting algorithm for dual monotone inclusions involving cocoercive operators,” _Adv. Comput. Math._ , vol. 38, no. 3, pp. 667–681, Nov. 2011.
  * [43] H. H. Bauschke and P. L. Combettes, _Convex Analysis and Monotone Operator Theory in Hilbert Spaces_. New York, NY, USA: Springer-Verlag, 2011.
  * [44] F. Bach, R. Jenatton, J. Mairal, and G. Obozinski, “Optimization with sparsity-inducing penalties,” _Found. Trends Mach. Learn._ , vol. 4, no. 1, pp. 1–106, Jan. 2012.


# An Automated Approach for the Recognition of Bengali License Plates

Md Abdullah Al Nasim1, Atiqul Islam Chowdhury2, Jannatun Naeem Muna3 and
Faisal Muhammad Shah4 Department of Computer Science and Engineering
Ahsanullah University of Science and Technology1,3,4, United International
University2
Dhaka, Bangladesh
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS>
###### Abstract

Automatic Number Plate Recognition (ALPR) is a system for automatically
identifying the license plates of any vehicle. This process is important for
tracking, ticketing, and any billing system, among other things. With the use
of information and communication technology (ICT), all systems are being
automated, including the vehicle tracking system. This study proposes a hybrid
method for detecting license plates using characters from them. Our captured
image information was used for the recognition procedure in Bangladeshi
vehicles, which is the topic of this study. Here, for license plate detection,
the YOLO model was used where 81% was correctly predicted. And then, for
license plate segmentation, Otsu’s Thresholding was used and eventually, for
character recognition, the CNN model was applied. This model will allow the
vehicle’s automated license plate detection system to avoid any misuse of it.

###### Index Terms:

License plate, YOLO, CNN, recognition, segmentation

##  I Introduction

Automatic Number Plate Recognition (ALPR) is a system by which the vehicles
can be detected using their license plate automatically through this system
without any direct human involvement. Basically, the manual process is so
redundant that, by testing them randomly, the police officer must keep a
record on the number of vehicles registered. ALPR represents that recognising
the characters from the license plate properly and also detect them accurately
[1]. ALPR system is a must one nowadays for vehicle speed measurement,
automated toll collection, prevention of traffic tules violation, security
controlling car parking in empty spaces, stolen vehicle search, control
traffic congestion cost and many more. ALPR system basically works on object
identification, preprocessing of that object image and correctly recognition
of the pattern of that image. Different kinds of filters, morphological
operations, transformation methods, pattern recognition are used in order to
pre-process the image and then we can apply the CNN to recognize the letters
correctly from the plates.

Figure 1: License plate of a vehicle

Bangladesh Road Transport Authority (BRTA) provides a registration number for
every vehicle. The vehicle category is indicated by the vehicle class letter
(KA, KHA, GHA, CHA) [2]. A lot of researches has taken place in the field of
ALPR. From 2013, for every vehicle, BRTA provides a common shape of licence
plate. That is called a digital license plate or Retro Reflective Licence
Plate. This kind of digital plate contains two rows. At lower row, there are
numbers and at the upper row, there are alphabets. In the lower row, there
remain two separate parts containing six digits. The last four digits from the
second part represent the vehicle number in Bengali and the first two digits
from the first part represent the vehicle class. On the other hand for the
upper row, a short form of the city name remains in the first word. Then for
metropolitan areas, the metro is the second word and for non-metropolitan
areas, there is no need to use the metro word in the license pates [3]. For
the last word of the upper row that represents the category of the vehicle
that means which type of vehicle is this.

##  II Related Works

Automatic License Plate Recognition (ALPR) is a type of technology that
enables computer systems to study automatically the registration number
(license number) of vehicles from digital pictures. Studying automatically the
registration number means transforming the pixels of the digital image into
the ASCII text of the number plate. Number plate detection is a popular work
and there were so many experiments regarding this. In this section, some
papers which are related and relevant to our work are being discussed.

For this research, many papers were followed which deal with Bengali Number
Plate Detection and segmentation. In [4], they applied chain code and neural
network. Sobel filter, morphological operations and connected component
analysis were used to extract the license plate, by scanning the image they
segmented the characters and finally, recognition was done by chain code. In
[3], they used edge detection, regional localization and character
segmentation for image processing and morphological operations in order to
determine Bengali characters of the licence plate. They used an algorithm
based on connected components.

In [5], they worked on real-time Bengali license plate detection by using
YOLOv3 algorithm foe localizing plate and recognizing the characters as well.
In [6], they used different image processing techniques to pre-process the
image like resizing, image binarization, enhancement, noise filtering etc.
Basically, they used a deep learning tool ‘The Convolutional Neural network’
to recognize the plate character and they also used thresholding to extract
the plate region. A paper used deep object detection convolutional neural
network [7] and another one used Local binary pattern (LBP), histogram
matching technique and matching bounding box technique [8].

In [9], they basically showed a comparison between Gabor filter bank and
fuzzified Gabor filter and the most interesting research they showed that
fuzzified Gabor filter is better to detect the complex license plate than the
Gabor filter bank. In [10] and [11], they proposed a new systematic method for
ALPR including a smart parking service also. In [10], they showed that their
proposed method has higher accuracy than filters Gaussian blur and filter 2D
and in [11], they preprocessed the image using sensor tools and convert the
whole image to a binary image and they also used template matching, Sobel
filter, dilation and median filter for ALPR.

After analysing those papers, it was found that all images need to be
preprocessed and then segmentation based different techniques need to be used
to achieve a good result. So, we developed a model covering all these steps in
this experiment.

##  III Methodology

For the complete detection process of characters from the license plate, three
stages were split. Initially, the You Only Look Once (YOLO) paradigm was used
to separate license plates from videos (mainly cars, buses on the road). Then
the image was cropped using the predicted coordination of the YOLO model. And
finally, using the CNN model, the digit of the segmented character was
predicted. We have used our own dataset in our research for the purpose of
training and testing. Our dataset comprises more than 37,700 images, which is
enormous for any experiment. The dataset was divided into 70% for training,
20% for validation and 10% for testing.

The complete experiment follows some measures that provide improved
performance for license plates and digits from them to be identified. The flow
diagram of the proposed system of license plate recognition is shown in Fig 2.

Figure 2: Flow diagram of our experiment

###  III-A Data Collection

The first step was to collect data of vehicles from roads. The images were
about the front or backside of vehicles. Also, data were collected for digit
recognition. The data were collected using our own mobile camera by recording
video which was converted into the image frame. And from the video, different
frames were separated which contain the image of the license plate.

###  III-B Data Preprocessing and Annotation

After collecting the data, the deletion of all unclear and unimportant images
was done. Then annotation was completed for the number plate of all vehicles
in the first dataset. To annotate the dataset, we utilize labelimg. The
dataset is manually annotated. After annotation, we use Additive Gaussian
Noise, Gaussian Blur, Salt-Pepper Noise, Coarse Dropout, Contrast
Normalization for license plate recognition.

###  III-C Apply YOLOv3 and Train Dataset

Objectness was used by YOLOv3 for bounding box prediction and cost function
measurement. For each bounding box using logistic regression, YOLOv3 predicts
an object score. The cost function is calculated differently in YOLOv3. And
this model was used for detecting the number plates which were trained on the
first dataset after the annotation.

###  III-D Licence Plate Detection and Localization

Following that, after training, the best-predicted model was used. The
experiment could save the best epoch of training. After 27 epochs, the
targeted detection model was found while training the data.

###  III-E Crop Licence Plate using Predicted Bounding Box

The saved model was used to detect the licence plate. And then, predicted
bounding box coordinate was used to crop licence plates from images.

###  III-F Apply Otsu’s Thresholding for Character Segmentation

The next step was to complete the segmentation. For segmentation, Otsu’s
thresholding algorithm was used. Otsu’s thresholding is an old but effective
one for the segmentation process.

###  III-G Apply CNN model for Character Recognition

A deep learning classification model, which is our CNN model, is used for
character recognition. For the segmented portion of images, the model iss
applied. With Keras library augmentation, Convolutional Neural Network(CNN) is
used where input conv1 size is 64*64*3(rows*cols*filters), kernel size is 5*5
where size-out is 60*60*64, and conv2 size-out is 56*56*64 where we use 5*5
kernel. We use a maxpool size of 2*2 after that. The layer of Conv3 is
28*28*64 which runs up to conv7 where kernel size is 3*3, turning to flatten
the 4036 layer. After flattening, the license plate is then associated with 32
dense neurons where we got the accuracy of character recognition is 99.43%.

The YOLOv3 model was trained as the best precision after 27 epochs and saved
after that epoch. In the training process, the model takes about 4.5 hours.
The method is for a license plate to be detected from an image. And the
updated CNN is trained at 100 epochs for character recognition. The training
stage takes approximately 8.5 hours.

##  IV Experimental Result

As we mentioned in our system, for this experiment, a total of 3 phases were
developed. For license plate detection, YOLOv3 was used in the first step.
Then, in the second level, segmentation was completed. And the last phase was
about the license plates’ identification of characters. Here, all the
observations of these three phases are illustrated.

Licence plate detection using YOLOv3 had an output accuracy of 0.81 mAP (Mean
Average Precision). That means 81% guess/prediction is correct. After training
with a lot of data, the test phase showed us a better result for license plate
detection. Images from different angles were used for testing purpose. And our
proposed model worked much better regarding those issues. Figure 3 and 4 show
the detection of a license plate from an image.

Figure 3: License plate recognition output (120 degree angle) Figure 4:
License plate recognition output (90 degree angle) Figure 5: Analysis of
YOLOv3 Model for License Plate Detection

These figures show that the YOLOv3 model provides a good output for the
detection of number plate at any angles. Besides this, from figure 5, there
are many different terms which shows a graphical view of the detection
process. GIoU measures the overlap between two boundaries of Bengali license
plate. Here, loss degrades in every step which helps us to lead a better
accuracy. In figure 5 of GIoU, we have predicted the bounding box with label
bounding box figure with shape and orientation of coverage area. It determines
the localization of Bengali license plate which is differentiable. In this
graph, objectness means how accurately we bounded the box which contains
license plate. Here our estimated object detection of binary cross-entropy
loss is decreasing which means our predicted license plate detection’s
bounding box works better. The classification works to classify how many
objects we detect here. We just detect Bengali license plates in this study,
thus no classification is required. In 3770 test images, there are 3580 images
where we find there are license plate and predict that there are 3270 images
were detected correctly. Here we accurately predict 3270 images and 170 images
which can not detect license plate accurately. Prediction of bounding license
plate detection and actual bounding box for license plate detection which is
called Mean Average Precision (mAP) is 0.81.

After that, the images of license plates were cropped and segmentation was
applied for the better look of every character on the plate. The last phase is
character recognition from those plates. Character recognition is not an easy
step as sometimes images are blurry, sometimes cropped etc. So these issues
need to be maintained for any experiment. Then we utilized our modified CNN to
recognize the characters, which had a validation accuracy of 99.43%. We have
experimented different ratio which are shown in Table I.

Figure 6: Segmentation of the license plate image

Figure 7: Recognition of the characters from the license plates TABLE I: Ratio of Training & Testing Data Ratio of Training & Testing Data | Accuracy (%)
---|---
70:30 | 91.72
80:20 | 99.43
85:15 | 97.14

Finally, it is obvious that the identification of license plates and
characters fits best with our model. It was a tough problem for the detection
system, as we used our own dataset from the Bangladesh vehicle. But despite
those problems, our hybrid model shows strong performance.

##  V Conclusion

Automatic Number Plate detection is a whole package of capturing the license
plate from the vehicle and to recognize it accurately. Basically, the purpose
of this system is to detect the Bengali licence plate to identify the vehicles
properly who violates the rules of the road. The whole procedure was divided
into three subprocedures. First YOLOv3 model was used for the detection of a
license plate from the selected images. Then predicted bounding box was used
for cropping the images of the license plate. For segmenting the character,
Otsu’s thresholding method was applied. Then our own modified deep learning
classification CNN model was used for character recognition. At the end, we
can say that this research provides a better result for the overall
recognition of license plate and character. As a large number of different
vehicles are now added and to handle them properly, providing this device is
vital for every corner of the city, the future work of this scheme will be
much brighter.

## References

  * [1] S. Abdullah, M. M. Hasan, and S. M. S. Islam, “Yolo-based three-stage network for bangla license plate recognition in dhaka metropolitan city,” in _2018 International Conference on Bangla Speech and Language Processing (ICBSLP)_. IEEE, 2018, pp. 1–6.
  * [2] G. Rabbani, M. A. Islam, M. A. Azim, M. K. Islam, and M. M. Rahman, “Bangladeshi license plate detection and recognition with morphological operation and convolution neural network,” in _2018 21st International Conference of Computer and Information Technology (ICCIT)_. IEEE, 2018, pp. 1–5.
  * [3] E. Ahmed, A. Islam, F. Sarker, M. N. Huda, and K. Abdullah-Al-Mamun, “A road to independent living with smart homes for people with disabilities,” in _2016 5th International Conference on Informatics, Electronics and Vision (ICIEV)_. IEEE, 2016, pp. 472–477.
  * [4] A. K. Ghosh, S. K. Sharma, M. N. Islam, S. Biswas, and S. Akter, “Automatic license plate recognition (alpr) for bangladeshi vehicles,” _Global Journal of Computer Science and Technology_ , 2011.
  * [5] S. Abdullah, M. M. Hasan, and S. M. S. Islam, “Yolo-based three-stage network for bangla license plate recognition in dhaka metropolitan city,” in _2018 International Conference on Bangla Speech and Language Processing (ICBSLP)_. IEEE, 2018, pp. 1–6.
  * [6] G. Rabbani, M. A. Islam, M. A. Azim, M. K. Islam, and M. M. Rahman, “Bangladeshi license plate detection and recognition with morphological operation and convolution neural network,” in _2018 21st International Conference of Computer and Information Technology (ICCIT)_. IEEE, 2018, pp. 1–5.
  * [7] A. Islam Chowdhury, M. Munem Shahriar, A. Islam, E. Ahmed, A. Karim, and M. Rezwanul Islam, “An automated system in atm booth using face encoding and emotion recognition process,” in _2020 2nd International Conference on Image Processing and Machine Vision_ , 2020, pp. 57–62.
  * [8] A. K. Bachchan, A. Gorai, and P. Gupta, “Automatic license plate recognition using local binary pattern and histogram matching,” in _International Conference on Intelligent Computing_. Springer, 2017, pp. 22–34.
  * [9] V. Tadic, Z. Kiraly, P. Odry, Z. Trpovski, and T. Loncar-Turukalo, “Comparison of gabor filter bank and fuzzified gabor filter for license plate detection,” _Acta Polytechnica Hungarica_ , vol. 17, no. 1, pp. 1–21, 2020\.
  * [10] N. d. V. Dalarmelina, M. A. Teixeira, and R. I. Meneguette, “A real-time automatic plate recognition system based on optical character recognition and wireless sensor networks for its,” _Sensors_ , vol. 20, no. 1, p. 55, 2020\.
  * [11] P. Dutta, A. Chakraborty, S. Bhadury, A. Sinha, S. Mitra, and S. Das, “Smart usage of open source license plate detection and using iot tools for private garage and parking solutions,” _Available at SSRN 3533565_ , 2020.