IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

1

Multi-resolution Graph Fourier Transform for
Compression of Piecewise Smooth Images
Wei Hu, Student Member, IEEE, Gene Cheung, Senior Member, IEEE, Antonio Ortega, Fellow, IEEE,
and Oscar C. Au, Fellow, IEEE

Abstract—Piecewise smooth (PWS) images (e.g., depth maps
or animation images) contain unique signal characteristics suchas sharp object boundaries and slowly-varying interior surfaces.
Leveraging on recent advances in graph signal processing, inthis paper we propose to compress PWS images using suitable
Graph Fourier Transforms (GFT) to minimize the total signal
representation cost of each pixel block, considering both thesparsity of the signal’s transform coefficients and the compactnessof transform description. Unlike fixed transforms such as the
Discrete Cosine Transform (DCT), we can adapt GFT to aparticular class of pixel blocks. In particular, we select one amonga defined search space of GFTs to minimize total representation
cost via our proposed algorithms, leveraging on graph optimization techniques such as spectral clustering and minimum graphcuts. Further, for practical implementation of GFT we introduce
two techniques to reduce computation complexity. First, at theencoder we low-pass filter and down-sample a high-resolution(HR) pixel block to obtain a low-resolution (LR) one, so thata LR-GFT can be employed. At the decoder, up-sampling andinterpolation are performed adaptively along HR boundaries
coded using arithmetic edge coding (AEC), so that sharp object
boundaries can be well preserved. Second, instead of computing
GFT from a graph in real-time via eigen-decomposition, themost popular LR-GFTs are pre-computed and stored in a table
for lookup during encoding and decoding. Using depth maps
and computer-graphics images as examples of PWS images,
experimental results show that our proposed multi-resolution(MR)-GFT scheme outperforms H.264 intra by 6.8 dB on averagein PSNR at the same bit rate.
Index Terms—Image Compression, Graph Fourier Transform,
Piecewise Smooth Images.

I. I NTRODUCTION
A popular approach to image compression is transform coding [1]: an image is first divided into non-overlapping blocksof pixels, with each block projected to a chosen transformdomain, and the resulting transform coefficients are quantized,
entropy-coded and transmitted to the decoder. While one canapply a transform to an entire image, most compression codecs(e.g., JPEG1 , H.264 [2]) employ transform coding block-byblock. This is because block-based transforms can adapt to
Copyright (c) 2013 IEEE. Personal use of this material is permitted.
However, permission to use this material for any other purposes must be
obtained from the IEEE by sending a request to pubs-permissions@ieee.org.
W. Hu and O. C. Au are with Department of Electronic and Computer
Engineering, The Hong Kong University of Science and Technology, Clear
Water Bay, Hong Kong, China (e-mail: huwei@ust.hk; eeau@ust.hk).
G. Cheung is with National Institute of Informatics (NII) and The Graduate
University for Advanced Studies (SOKENDAI), 2-1-2, Hitotsubashi, Chiyodaku, Tokyo, Japan 101–8430 (e-mail: cheung@nii.ac.jp).
A. Ortega is with Department of Electrical Engineering, Universityof Southern California, Los Angeles, CA 90089 USA (e-mail: antonio.ortega@sipi.usc.edu).
1 http://www.jpeg.org/jpeg/

the non-stationary statistics of natural images and are moreamenable to hardware implementation. We focus on blockbased transform coding of images in this paper.
A key to good compression performance is to select atransform that promotes sparsity in representation; i.e., uponsignal transformation there remain few non-zero transform
coefficients, resulting in coding gain. For a vector following a known correlation model (statistically described by acovariance matrix), the Karhunen-Loève Transform (KLT) is
the optimal transform in terms of decorrelating componentsof the vector2 [7]. In order to use the KLT for coding of
non-overlapping pixel blocks in an image, one would first
adaptively classify each block into one of several statisticalclasses, then apply the KLT corresponding to the identified
statistical class for decorrelation. However, description of theidentified statistical class for the coded block (and hence
the transformation employed) must be conveyed to the decoder for correct transform decoding. Given this transform
description overhead, how to suitably select a good set of
statistical classes for KLT transform coding to achieve goodoverall performance—in terms of both sparsity of signal
representation and compactness in transform description—isthe challenge.
The Discrete Cosine Transform (DCT), on the other hand,is a fixed transform with no signal adaptation, thus it requires
no description overhead. It is equivalent to the KLT for the
single statistical class where the correlation between each pairof adjacent pixels is assumed to be 1. The DCT is widelyused in image compression systems, such as the JPEG and
H.26x codecs [8]. However, the assumption that the pixel
correlation tends to 1 is not always true. In particular, adjacentpixels across sharp boundaries clearly have correlation muchsmaller than one. This means that for blocks that straddle
sharp boundaries, the DCT will lead to non-sparse signal
representation, resulting in poor coding performance.
In order to achieve a balance between signal representationand transform description, we propose to use graph Fouriertransforms (GFTs) [9] for transform coding of images. GFT
is a matrix of eigenvectors of the graph Laplacian matrix ofa constructed graph. For transform coding, the graph nodes
represent pixels in a coding block, and each edge weight represents correlation between two connected (adjacent) pixels.
2 Optimality in decorrelating a vector statistically, however, does not necessarily imply optimality in transform coding. It is shown that the KLT is
optimal if the input vector is jointly Gaussian or a mixture of Gaussians, forvariable-rate high-resolution coding [3], [4], [5]. For other sources, the KLT
may yield sub-optimal performance [6].

IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

Given a defined search space of graphs that connect pixels ina given block, an “optimal” graph (and hence corresponding“optimal” GFT) minimizes the total signal representation
cost, evaluating both the sparsity of the signal’s transform
coefficients and the compactness of transform description.

(a) Teddy(b) DudeFig. 1. Examples of PWS images: (a) a depth map of an image sequence
Teddy; (b) a computer-graphics image Dude.

In particular, we propose to use GFTs for transform coding of a class of images—piecewise smooth (PWS) images.
PWS images consist of sharp boundaries and smooth interior
surfaces. They have recently received renewed attention in theimage processing community, due to advances in depth sensing
technologies, such as Microsoft Kinect R , that capture fairly
accurate depth maps in real-time. A depth image is a per-pixel
map that measures the physical distance between objects in the
3D scene and the capturing camera. Because textural contents
of objects are not captured, depth images are smooth within
an object’s surface, and hence PWS. Depth maps can be used
for advanced image processing tasks, such as foreground /background segmentation, 3D motion estimation [10], depthimage-based rendering (DIBR) [11], etc. Other examples of
PWS images include animation images, magnetic resonanceimaging (MRI) images, dense motion fields [12], computergraphics images, etc. See Fig. 1 for examples of PWS images.
Unique signal characteristics of sharp boundaries and smoothinterior surfaces in PWS images mean that conventional coding schemes designed for natural images such as JPEG do
not offer good compression performance. We focus on thecompression of PWS images in this paper.
Our proposed GFT-based coding scheme for PWS images, called multi-resolution GFT (MR-GFT), has three main
contributions. First, we consider a large space of possible
graphs, including unweighted graphs and weighted graphs,
to exploit the generality afforded by GFT. Further, we define
an optimality criterion of GFT for compression that aims tominimize the total representation cost, including transformdomain representation and transform description. To enable
low GFT description cost, we restrict our feasible space of
GFTs to graphs with simple connectivity and a small discrete
set of edge weights. Discrete edge weights are derived froma statistical analysis under a model specifically designed tocharacterize PWS signals, for which we demonstrate that thederived GFT approximates the KLT.
Second, given a defined search space of graphs (or equivalently GFTs), we search for the optimal GFT for each pixel
block with our proposed efficient algorithms. In particular,
for ease of optimization we divide the search space into
two subspaces—GFTs on weighted graphs and GFTs on un-

2weighted graphs. For GFTs on weighted graphs, we formulate
a separation-deviation (SD) problem and solve it using aknown efficient graph cut algorithm [13]. Meanwhile, we
search for good GFTs on unweighted graphs via a greedy
algorithm, leveraging on a graph partitioning technique basedon spectral clustering [14].
Third, for practical implementation, we design two techniques to reduce computation complexity. In the first technique, we propose a multi-resolution (MR) approach, wheredetected object boundaries are encoded in the original highresolution (HR), and smooth surfaces are low-pass-filtered anddown-sampled to a low-resolution (LR) one, before performing
LR GFT for a sparse transform domain representation. At
the decoder, after recovering the LR block via inverse GFT,
we perform up-sampling and interpolation adaptively along
the encoded HR boundaries, so that sharp object boundaries
are well preserved. The key insight is that on average PWS
signals suffer very little energy loss during edge-adaptive lowpass filtering, which enables the low-pass filtering and downsampling of PWS images. This MR technique also enables us
to perform GFT on large blocks, resulting in large coding gain.
In the second technique, instead of computing GFT from agraph in real-time via eigen-decomposition of the graph Laplacian matrix, we pre-compute and store the most popular LRGFTs in a table for simple lookup during actual encoding anddecoding. Further, we exploit graph isomorphism to reduce thenumber of GFTs required for storage to a manageable size.
Using depth maps and graphics images as examples of PWS
images, experimental results show that our proposed MRGFT scheme outperforms H.264 intra (with intra predictionused) by 6.8 dB on average in PSNR. Further, when applying
compressed depth maps to DIBR, we achieve 2.2 dB gain on
average in synthesized view PSNR compared to H.264 intra.
The outline of the paper is as follows. We first discuss related work in Section II. We then overview our proposed MRGFT coding system in Section III and review the definition of
GFT in Section IV. Next, we present a mathematical analysis
of optimal GFTs in Section V. The problem formulations andalgorithms for optimal GFTs on two types of graph supports
are discussed in Section VI. We then outline the detailed
implementation in Section VII. Finally, experimental results
and conclusions are presented in Section VIII and Section IX,
respectively.
II. R ELATED W ORK
The field of PWS image compression encompasses diverse
source coding schemes, ranging from segmentation-based
techniques to popular transform approaches. In general, they
feature efficient representations of the geometric structures in
PWS images. We discuss them in order as follows.
Segmentation-based compression schemes segment the image into homogeneous partitions followed by coding of eachpartition. One of the most popular is quadtree-based compression [15], [16], [17], which recursively divides the image
into simple geometric regions. [15] designs the quadtree
segmentation in an optimal Rate-Distortion (RD) frameworkand then approximates each segment using a polynomial

IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

model separated by a linear boundary. The direction of theboundary, however, is chosen from a given limited discreteset, which may not be sufficiently descriptive to describe
arbitrarily shaped boundaries. In contrast, our proposed MRGFT can represent any boundary accurately and efficiently via
arithmetic edge coding (AEC) of HR edges [18], [19]. In [16]depth maps are modeled by piecewise-linear functions (called
platelets) separated by straight lines, which are adaptive toeach subdivision of the quadtree with variable sizes in a global
RD tradeoff. However, this representation has non-vanishingapproximation error, since depth maps are not exactly piecewise linear. In contrast, our proposed MR-GFT can accurately
represent PWS images in GFT domain.
Transform approaches are also designed for PWS image
compression. An ensemble of transforms exploit the geometryof PWS images, such as the wavelet-domain compression [20],[21], curvelet [22], and contourlet [23]. However, they are
all deployed over the entire image, and so are not easily
amenable to block-based processing for hardware-friendly
implementation. In contrast, our proposed MR-GFT is a blockbased transform coding approach.
GFT is first used for depth map coding in [24], whichempirically demonstrated that using GFT for transform codingof a depth pixel block with sharp boundaries, one can achievea sparser transform domain representation than with DCT.
(Sparsity in representation of depth block in GFT domain isalso used as a signal prior in depth image denoising in [25].)
However, [24] has the following shortcomings. First, given an
input pixel block, [24] considers only a single variant of GFT
based on an unweighted graph, while we consider a muchlarger space of possible graphs, including both unweighted
and weighted graphs. Second, graphs are deduced directlyfrom detected boundaries in [24], thus it is not clear if theconstruction of graphs is optimal in any sense. In contrast,
we formally define an optimality criterion that reflects representation cost and propose efficient algorithms to √search√forgraphs that lead to optimal GFTs. Finally, for each N × N
pixel block, [24] requires real-time eigen-decomposition of an
N × N graph Laplacian matrix to derive basis vectors at boththe encoder and decoder. This large computation burden alsomeans that GFT cannot be practically performed for larger
blocks in PWS images, which will otherwise result in larger
coding gain. To address this issue, we design two techniques
to reduce computation complexity, namely, the MR schemeand table lookup. Experimental results show that we achieve
5.9 dB gain on average in PSNR over [24].
Drawing a connection between GFT and a graphical statistical model called Gaussian Markov Random Field (GMRF),[26] provides a theoretical analysis on the optimality of GFT
in terms of decorrelating the input vector under the GMRF
model. Unlike [26], we select a GFT for each coding blockby choosing from a set of candidate GFTs using an optimality
criterion that considers both the sparsity of the GFT domain
signal representation and the compactness of the transform
description.
Image compression specifically tailored for depth maps, a
representative class of PWS images, has been studied extensively in recent years. Beyond the class of techniques that3

exploit the piecewise smooth characteristics of depth signals[16], [17], [21] as previously discussed, another class of
depth coding techniques are tailored specifically for virtual
view rendering via DIBR [27], [28], [29], [30], [31], [32].
Specifically, they optimize depth coding by considering theimpact of depth errors due to lossy image compression on the
DIBR-synthesized views. [27], [28] utilize synthesized view
distortion as metric for depth coding optimization, such asmode selection. [29], [30] use synthesized view distortion as
objective for transform domain sparsification (TDS). In [31],view synthesis is incorporated for improved depth prediction
in a RD optimized framework. [32] proposes down-sampling/ up-sampling techniques in an optimized fashion by incorporating the effect of depth re-sampling on view synthesis
quality. In contrast, since we focus on compression of general
PWS images, we do not consider synthesized view distortion
as an optimization metric. Nevertheless, our proposed method
implicitly leads to good synthesized views via edge-preservingof depth maps, as shown in Section VIII.
III. M ULTI - RESOLUTION G RAPH F OURIER T RANSFORMC ODING S YSTEM
Encoder
DCT
Original
PWS Signal

Intra
prediction
LP

K

Edge
Detection

LR-GFT(Lookup)Mode coeff.
DecisionQ
index

Entropy
Enoding

HR Edge
Map

HR Edge Map
Intraprediction
Interpolation
Reconstructed
PWS Signal

K+

index Entropy

InverseLR-GFT

DeodingNindex=0?

InverseDCT-1 coeff.

Q

Y

DecoderFig. 2.

MR-GFT coding system for PWS images.

We first provide an overview of our proposed MR-GFT
coding system for compression of PWS images, shown in
Fig. 2. Given a PWS image, we discuss the encoding anddecoding procedures as follows.
A. Encoder
At the encoder, we first detect prominent boundaries (large
inter-pixel intensity difference) in the HR image via hardthresholding of image gradients. The threshold is set basedon the mean and variance of the image, so that the boundary
detection is adaptive to the image statistics. We encode HR
boundaries losslessly for the adaptive intra prediction andinterpolation at the decoder (discussed later), using AEC [18],[19]. We apply AEC for the entire image, which avoids initialization for each block and efficiently encodes√√ long continuous
boundaries in the image. Then for each K N × K N target
pixel block considering a down-sampling factor K = 2, we
execute the following three steps.

IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

First, we perform edge-aware intra prediction as proposedin [33]. Different from the intra prediction in H.264 [2], edgeaware intra prediction efficiently reduces the energy of theprediction error by predicting within the confine of detected
HR boundaries, thus reducing bits required for coding of theresidual signal.
Second, we try two types of transforms for transform codingof the residual block: i) fixed DCT on the original HR residual
block (HR-DCT), ii) a pre-computed set of LR GFT (LRGFT) on the down-sampled LR residual block (including
LR weighted GFT and LR unweighted GFT, as discussed
in Section V). We then choose the one transform with thebest RD performance. Before transform coding using LR-GFT,
however,√ adaptively low-pass-filter√and down-sample√√ we first
the K N × K N block uniformly to a N × N block.
Low-pass filtering is first used to avoid aliasing caused bydown-sampling. We propose an edge-adaptive low-pass filterin the pixel domain for the preservation of sharp boundaries.
Specifically, a pixel is low-pass-filtered by taking average of
its neighbors on the same side of HR boundaries within a(2K − 1) × (2K − 1) window centering at the to-be-filteredpixel. The advantage of this edge-adaptive low-pass filtering isthat filtering across arbitrary-shape boundaries will not occur,
so pixels across boundaries will not contaminate each otherthrough filtering.
For the implementation of the HR-DCT and LR-GFT, we
pre-compute the optimal transforms (discussed in Section VI)
and store them in a lookup table (discussed in Section VII)a priori. During coding, we try each one and choose the one
with the best RD performance. The two types of transforms,
HR-DCT and LR-GFT, are employed to adapt to different
block characteristics. HR-DCT is suitable for blocks whereedge-adaptive low-pass filtering would result in non-negligible
energy loss. If very little energy is lost during low-passfiltering, LR-GFT would result in a larger coding gain. Notethat if a given block is smooth, the LR-GFT will default to
the DCT in LR, and would generally result in a larger gain
than HR-DCT due to down-sampling (the rates of transform
indices for both, i.e., the transform description overhead, arethe same in this case).
Third, after the RD-optimal transform is chosen from thetwo transform candidates, we quantize and entropy-encode theresulting transform coefficients for transmission to the decoder.
The transform index identifying the chosen transform is also
encoded, so that proper inverse transform can be performed at
the decoder.

4

the same side of boundaries within a (2K − 1) × (2K − 1)window centering√ at pixel√x.
Finally, the K N ×K N block is reconstructed by adding
the intra predictor to the residual block.
IV. G RAPH F OURIER T RANSFORM F OR I MAGES
Before we proceed to problem formulation and algorithms
for optimal GFTs, we first review the basics of GFT.
A graph G = {V, E, W} consists of a finite set of vertices Vwith cardinality |V| = N , a set of edges E connecting vertices,
and a weighted adjacency matrix W. W is a real N × N
matrix, where Wi,j is the weight assigned to the edge (i, j)connecting vertices i and j. We consider here only undirected
graphs, which correspond to symmetric weighted adjacency
matrices, i.e., Wi,j = Wj,i . We also assume weights are nonnegative, i.e., Wi,j ≥ 0.
While there exist different variants of Laplacian matrices,
we are interested in the unnormalized combinatorial graphLaplacian in this work, which is defined as L := D − W,
where D is the degree matrix—a diagonal matrix whose ith
diagonal element is the sum of all elements in the ith row
NP
of W, i.e., Di,i =Wi,j . Since the Laplacian matrix is
j=1

a real symmetric matrix, it admits a set of real eigenvalues{λl }l=0,1,...,N −1 with a complete set of orthonormal eigenvectors {ψl }l=0,1,...,N −1 , i.e., Lψl = λl ψl , for l = 0, 1, ..., N − 1.
We employ this Laplacian matrix for two reasons.
First, because elements in each row of L sum to zeroby construction, 0 is guaranteed to be an eigenvalue with
[1 . . . 1]T as the corresponding eigenvector. This means afrequency interpretation of GFT, where the eigenvalues λl ’s
are the graph frequencies, will always have a DC component,
which is beneficial for the compression of PWS images wheremost regions are smooth.
Second, GFT defaults to the well known DCT when definedfor a line graph (corresponding to the 1D DCT) or a 4connectivity graph (2D DCT) with all edge weights equal to1 [26]. That means GFT is at least as good as the DCT insparse signal representation if the weights are chosen in thisway. Due to the above two desirable properties, we use theunnormalized Laplacian matrix in our definition of GFT.
We note that the graph Laplacian can be used to describe
the total variation of the signal with respect to the graph; i.e.,
for any signal x ∈ RN residing on the vertices of a graph withthe graph Laplacian L, we can write [35]N

xT Lx =
B. Decoder
At the decoder, we first perform inverse quantization andinverse transform for the reconstruction of the residual block.
The transform index is used to identify the transform chosen
at the encoder for transform coding.
Secondly,
LR-GFT is employed, we up-sample the recon√ if √N×N LR residual block to the original resolution
structed√√
K N ×K N , and then fill in missing pixels via our proposed
image-based edge-adaptive interpolation [34], where a pixel x
is interpolated by taking average of its neighboring pixels on

N

1 XX
Wi,j (xi − xj )2 .
2 i=1 j=1

(1)

xT Lx is small when x has similar values at each pair ofvertices i and j connected by an edge, or when the weight Wi,jis small for an edge connecting i and j with dissimilar values.
Thus, a signal is smooth (mostly low-frequency components)
with respect to a graph if the edge weights capture thesimilarity of connected pixels in the signal. Since xT Lx issmall if a signal is smooth on the graph (thus a sparse
representation in the GFT domain), we will use xT Lx to selectgraphs that lead to good compression results in later sections.

IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

5

trade-off between the coding rate and total induced distortion.
Specifically, a typical system design seeks to minimize aweighted sum of rate and distortion for chosen weighting
parameters. However, assuming high bit rate, [36] shows thata uniform quantizer yields the following expected distortion:
D=
Fig. 3. An example of constructing GFT from a 2×2 pixel block. The verticalcontour separates pixel 1 and 2 from pixel 3 and 4, and a graph is constructedby connecting pixels on each side of the contour. The corresponding adjacency
matrix W, degree matrix D, Laplacian matrix L as well as the computed
GFT U are shown on the right.

In particular, we will discuss how we use (1) as the rate proxyof transform coefficients during the search for optimal GFTs
in Section VI.
The eigenvectors {ψl }l=0,1,...,N −1 of the Laplacian matrix
are then used to define the GFT. Formally, for any signal x ∈RN residing on the vertices of G, its GFT x̂ ∈ RN is definedin [9] asx̂(l) = < ψl , x > =NXψl∗ (n)x(n), l = 0, 1, ..., N − 1. (2)n=1

The inverse GFT follows asx(n) =

N
−1
X

(3)

(4)

where q is the quantization step size employed for eachcoefficient, and N is the toal number of coefficients, whichis the same for different orthogonal transforms. This indicatesthat the expected distortion does not change when considering
different transforms under the same assumptions.
Hence, we only need to consider the total coding rate. Givena pixel block x ∈ RN , the GFT representation derived from thegraph G = {V, E, W} constructed on x has two representation
costs: i) the cost of transform coefficient vector α denotedby Rα (x, W), and ii) the cost of transform description T
denoted by RT (W). We thus arrive at the following definitionof optimality for GFTs:
Definition For a given image block x ∈ RN under fixed
uniform quantization at high bit rate, an optimal GFT is theone that minimizes the total rate, i.e.:minRα (x, W) + RT (W)s.t.

Wi,j ∈ CW

x̂(l)ψl (n), n = 1, 2, ..., N.

N q2
,
12

∀ i, j ∈ V

(5)

l=0

Having defined the GFT, we give an example of how toconstruct a GFT given
graph defined for a√ an unweighted√
pixel block. Given a N × N pixel block, we first treat
each pixel in the block as a vertex in a graph G and connect it
to its four immediate neighbors, resulting in a 4-connectivity
graph. See Fig. 3 for an illustration. Second, if there is a large
discrepancy in values between two neighboring pixels, we
eliminate their connection. Given the connectivity graph, we
can define the adjacency matrix W, where Wi,j = Wj,i = 1
if pixel positions i and j are connected, and 0 otherwise. The
degree matrix D can then be computed. In the third step, using
computed W and D, we compute the graph√Laplacian√ matrix
L = D − W. We then stack pixels in the N × N patch
into a length-N vector x and compute the GFT according to(2).
V. O PTIMAL G RAPH F OURIER T RANSFORMS FOR
P IECEWISE S MOOTH I MAGES
We now define a notion of optimality of GFT for compression of PWS images. Towards a formal definition of an
optimization problem, we then define the search space for
GFT as a discrete set of graph supports with edge weights
drawn from a small discrete set. The weights are later derived
assuming a specific statistical model for PWS images, where
we demonstrate that the derived GFT approximates the KLT
under this model.
A. Optimality Definition for Graph Fourier Transforms
In lossy image compression, different coding systems arecompared based on their RD performance, which describes thewhere C is the feasible set of edge weights.
Note that while an edge weight could in general take on any
non-negative real value, we restrict weights to a small discrete
set C in order to enable low description cost RT for GFT. Thisis further discussed next.
B. Definition of Search Space for GFTs
To lower the description cost RT for GFT in optimization (5), instead of a full graph we assume a simpler 4connectivity graph, where each vertex is only connected toits four adjacent neighbors. Further, we view pixels that areneighbors in the 4-connectivity graph as random variables, andconsider a discrete set of possible weights to be chosen foreach edge connecting two pixels. These weights correspondto three representative classes of the correlation between two
pixels: 1) strong correlation between the pixels, which models
pixel pairs occurring in smooth regions of the foreground/ background; 2) zero correlation between the two pixels,
occurring when they straddle sharp boundaries between theforeground and background; and 3) weak correlation between
the pixels that correspond to distinctly different parts of thesame foreground / background object. See Fig. 4 for an
illustration. Correspondingly, the weights are assigned to be
C = {1, 0, c}, where c is a carefully chosen real value between
0 and 1 (to be formally derived next).
Having defined the edge weight set, for ease of computation
we further divide the optimization problem (5) into two subproblems with two corresponding non-overlapping GFT search
spaces:IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

6

general 2D setting where the locations of smooth and weaktransition regions S and P are known deterministically, andthe statistics of the random variables involved in the model
are known probabilistically.
Given the unique characteristics of PWS images, we assume a piecewise first-order autoregressive process x =[x1 , ..., xN ]T with independent and identically distributed(i.i.d.) standard normal noise ek ∼ N (0, 1), k = 1, ..., N [37]:

k=1
 η,
xk−1 + ek ,1 < k ≤ N, [k − 1, k] ∈ Sxk =xk−1 + g + ek , 1 < k ≤ N, [k − 1, k] ∈ P

Fig. 4. An intuitive illustration of how different transforms adapt to pixel
blocks with different pixel correlations in the depth map of Ballet.

Unweighted GFT (UGFT), with C = {1, 0} for blocks
that can be well described by pixel pairs with strong andzero correlations only; and• Weighted GFT (WGFT), with C = {1, c} for blocks thatcan be well described by pixel pairs with strong and weakcorrelations only.
If a block can be well described by pixel pairs with strong
correlation only (i.e., a smooth block), then the GFT defaults
to the DCT. See Fig. 4 for an illustration. In a nutshell, the
WGFT considers only graphs that describe a single connectedcomponent (i.e., only one DC coefficient in any WGFT
representation of the signal). The UGFT considers graphsof multiple connected components, where each component isconnected by edges with weight 1. (A connected component
with a pair of disconnected neighboring pixels appears rarely
in practice, and thus is not considered in the UGFT optimization.) Both search spaces are searched for a given input
coding block for the best possible signal representation. The
rationale for this division is twofold. First, if the coding block
is sufficiently small, then the likelihood of a block containing
pixel pairs with all the three classes of correlations is verysmall. Second, division into two independent searches for the
UGFT and WGFT leads to efficient search algorithms, which
will be presented in Section VI.
The key question now is: what is the most suitable value ofc, one that correctly models weak correlation in the third classof pixel pairs, so that an optimal GFT in the WGFT search
space has the best possible performance? We next derive theoptimal c from statistical analysis under a model designed
for PWS signals, where we demonstrate that the derived GFT
approximates the KLT for the corresponding class of signals.
•(6)
where we assume the first variable x1 to be η ∼ N (0, σ12 ).
In the given smooth region S, xk is assumed to be the sumof xk−1 and a standard normal noise, while across the known
weak transition region P (e.g., from one smooth region of theforeground to another, as shown in Fig. 4), xk is modeled
as the sum of xk−1 , a random gap g ∼ N (mg , σg2 ), and astandard normal noise.
Further, we assume there exists only one weak transition
in a code block. Since experimentally we use a sufficiently
small block size, it is very likely to have at most one weaktransition in a block in each dimension. Assuming the only
weak transition exists between xk−1 and xk , then (6) can be
expanded as:x1 = η
x2 − x1 = e2···(7)xk − xk−1 = g + ek···
xN − xN −1 = eN
and further written into the matrix form:
Fx = b,

(8)

where


1
−1


0
F=
0


0
0

0
1
..
.
0

0
0
..
.
−1

0
0

0
0

0
0
..
.
−1


   00η e2  00   
 ..   ..    0
,b =  .  + .
 ek  g 0    .  .
 ..   .. 
0
1
eN
0

(9)

0
1
1
..
.
1

···
···
···
..
.
···

(10)

0
0

0
0

0
1
..
.
0

Since F is invertible with

1
1
1
F−1 = 
.
.
.
1

0
0
1
..
.
1


0
0
0
,
.. 

.
1

we have

C. Derivation of Optimal Edge Weights for Weak Correlation

x = F−1 b.

For simplicity, we consider the derivation of the optimaledge weight c in one dimension. We note, however, that the
optimality of our derived edge weight c—by extension the
optimality of the corresponding GFT—carries over to the more

In order to approximate the optimal decorrelation propertyof the KLT for a particular ensemble of 1D signals as x under
our specific model, we first compute the covariance matrix of
x. Assuming µ = E[x], we have(11)IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

7

C = E[(x − µ)(x − µ)T ]= E[xxT ] − µµT(12)= E[F−1 bbT (FT )−1 ] − µµT= F−1 E[bbT ](FT )−1 − µµT ,where µ takes the form

0,
1≤i<k
µi =
mg , k ≤ i ≤ N(13)Further computation gives σ2 0
1

0 1

···
···

0
0

···
···

0
0





..




.
 0 0 ··· 1
T
0
···
0 
E[bb ] =  0 0 ··· 0 σ2 +m2 +1 0 ··· 0 gg 0 0 ··· 001 ··· 0 ..
.
0 0 ···

0

0(14)
Fig. 5.

A one-dimensional graph model.

0 ··· 1pairs with similar weak correlation based on the discontinuity
magnitude, and then compute the sample variance of each classσ1 σ1σ1 ···
σ1
σ12
···
σ12
of pixel pairs as the estimation of σg2 .
..
..
..
..

 .. ..
To summarize, an optimal graph—one that leads to a GFT
 . .

.
.
.
 σ12 σ12 +1 σ12.+2 ··· σ12 +k−2
σ12 +k−2 ···σ12 +k−2 thatapproximates the KLT—is one with edge weights 1 except
 2 2
σ σ +1 σ12 +2 ··· σ12 +σg2 +k−1 σ12 +σg2 +k−1 ··· σ12 +σg2 +k−1 C=
theweight(17) of the edge that connects two pixels in the 12 12 σ1 σ1 +1 σ12 +2 ··· σ12 +σg2 +k−1 σ12 +σg2 +k ··· σ12 +σg2 +k  transition region P.
 . .

..
..
..
..
 . .

.2 2.
.
.
.2
.
2
2
2
2
2
2
σ1 σ1 +1 σ1 +2 ··· σ1 +σg +k−1 σ1 +σg +k ··· σ1 +σg +N −1
VI. A DAPTIVE S ELECTION OF G RAPH F OURIER(15)
T RANSFORMS
In order to relate the KLT to GFT, we consider the corHaving derived the optimal edge weight c for weakly
responding precision matrix Q, i.e., the inverse of C. Q
correlated
pixel pair in a coding block—and by extension thecomputes to a tridiagonal matrix:
optimalWGFT
for decorrelation, ideally one can simply use 1+ 12 −1
σ1that optimal WGFT for transform coding. However, this does
 not account for the cost of transform description RT in (5).
 −1 2 −1
.. .. ..
 To account for both the signal representation cost in the GFT

.
.
.


−1 2
−1
 domain and the desription cost of the chosen GFT, as discussed

1
1
−1
+1
−


2 +1
2 +1
σg
σg
Q = C−1 =  in Section V, we search for the optimal GFT that minimizes11
− σ2 +1 σ2 +1 +1 −1
 (5), and further divide the search space into two subspaces:

g
g


−1
2 −1
 1) WGFT on weighted and connected graphs, and 2) UGFT
.. .. ..

. . .
on unweighted and disconnected graphs. In this section we−1 2 −1−1 1
present the respective problem formulations and algorithms in(16)
the two subspaces in detail.
1
The first element 1 + σ2 is close to 1, since the variance1
σ12 of the first pixel (in the absence of any other information)
tends to be very large. Then for the 2-connectivity graph in A. Weighted Graph Fourier Transform
The WGFT is derived from a weighted and connected graphFig. 5, if we assign Wk−1,k = σ21+1 and all the other edgeg
G={V, E, W}, where the feasible weight set is C = {1, c}.
weights 1, Q is approximately equivalent to the combinatorialFor
easeof optimization, we first propose plausible proxies
Laplacian matrix L from its definition in Section IV. Since
for
the
two
rate terms in (5), i.e., the rates for transform
Q shares the same set of eigenvectors with C [26], i.e., thecoefficients
andtransform description. Then we cast the opbasis vectors of the KLT, the derived GFT is approximately
timization
problemof WGFT as a separation-deviation (SD)the KLT for the class of signals with the covariance matrix C
problem,
in
order
to
leverage on the well known and efficientdefined as in (15).
algorithm
developed
for the SD problem in [13].
Hence, the optimal edge weight for weak correlation that
1)
Rate
Proxy
of
Transform
Coefficients: Transform coefleads to the decorrelation GFT isficients
consist
of
DC
and
AC
coefficients.
Since the WGFT
1.
(17) is constructed on a connected graph (i.e., only one connected
c = Wk−1,k = 2
σg + 1component in the graph), the zero eigenvalue is of multiplicity
2
In practice, we estimate σg from collected pixel pairs of weak one, and each WGFT produces only one DC coefficient forcorrelation from training images. Specifically, we classify pixel a given block x. The cost of DC coefficients of the WGFT
Finally C is computed to be
 2 2
2
2

IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 20158

is thus the same. We then approximate the cost of quantized
AC coefficients (corresponding to non-zero eigenvalues), fora given quantization step size q and without consideration ofrounding for simplicity, as follows:xT Lx/q 2 = xT (N
−1X

λl ψl ψl T )x/q 2

l=0

=

N
−1
X

λl (xT ψl )(ψl T x)/q 2(18)

We = Wv1 (e),v2 (e) . A link is drawn between two nodes inthe dual graph if the two corresponding edges in the originalgraph share the same vertex as an endpoint, or traverse thesame contour in the coding block. See Fig. 6 for an exampleof the conversion.
Given the node set V d in the dual graph, (19) can now be
rewritten as
XRα (x, W) = ρWe (xv1 (e) − xv2 (e) )2 ,
(20)
e∈V d

l=0

=

N
−1
X

λl (αl /q)2 ,l=0

where αl , l = 1, ..., N − 1, is the l-th transform coefficient.
In words, xT Lx/q 2 is an eigenvalue-weighted sum of squared
quantized transform coefficients. This also means that the DC
coefficient is not reflected in (18). By minimizing (18), we
suppress high-frequency coefficients. Recall that xT Lx canalso be written in the form in (1) in terms of edge weights
and adjacent pixel differences. We thus propose the following
proxy for the cost of quantized AC coefficients:NRα (x, W) =

=

N

1 XX
Wi,j (xi − xj )2 /q 2
2 i=1 j=1
NN
1 XX
Wi,j (xi − xj )2 ,
ρ
2 i=1 j=1

(19)

where ρ = 1/q 2 . We note that [38] proposed to use (19) as thecost function for all GFT coefficients, even for disconnected
graphs. As discussed, (19) excludes the cost of DC coefficients, which means that (19) fails to capture the variable costof DC coefficients in cases where the graph contains a variable
number of connected components. In our WGFT search space,
we avoid this problem since our search space contains only
graphs that have a single connected component, because 0 is
excluded as a possible choice for an edge weight. Note thatwithout this restriction, the all zero connection graph would
be chosen as the optimal choice, which obviously does notlead to an optimal GFT in general.

Fig. 6. An illustration for converting edges in original graph to vertices inthe dual graph for SD formulation.

Further, for ease of later algorithm development, we convert
edges in the graph G to nodes in a dual graph3 G d , and define
the AC coefficient cost (19) on G d instead. More specifically,
we map each edge e that connects vertices v1 (e) and v2 (e)in G to a node4 ved in G d , and assign node ved with the value
3 The dual graph is also termed the line graph in the mathematical disciplineof graph theory.
4 To avoid confusion, we use the terminology of vertex and edge in theoriginal graph, and node and link in the dual graph.

where a label We assigned to node e in the dual graph G d isthe edge weight Wv1 (e),v2 (e) in the original graph G. A label
assignment W to a node e in a graph inducing difference costis typically called a deviation cost in an SD problem [13].
2) Rate Proxy of Transform Description: Weighted graphs
for a pixel block are defined by edges assigned weak correlation c; the rest of the edges are assigned weight 1. Further,edges of weak correlation tend to collectively outline an
object’s contour, i.e., fractional weight edges tend to appear
consecutively. Hence, a rare WGFT—one selected with low
probability p or high entropy cost −p log p—will have weakcorrelation edges scattered throughout the block. We thus
propose a rate proxy to penalize these rare GFTs. Specifically,
we use the label differential cost of nodes connected by linksddEP in the dual graph G as the rate proxy of WGFT, i.e.,(e,s)∈E d |We − Ws |. Assigning different labels to nodes
connected by a link will induce a separation cost in an SD
problem [13].
Further, since pixels in a PWS image are more likely to bestrongly correlated, we penalize a block with larger number
of weak correlation edges by adding another term to the rate
proxy:X
X
RT (W) =|We − Ws | +γρ(1 − We ), (21)(e,s)∈E de∈V dwhere ρ is the same one as in (19), and γ is used to assign theimportance of the second term relative to the first term. While
the first term specifies the cost of encoding a given distributionof edges (consecutive fractional weight edges consume fewer
bits as discussed earlier), the second term essentially counts
the number of fractional weight edges in the block to encode.
The actual relative importance between these two terms, i.e.,
the assignment of γ, depends on the particular edge encoding
scheme.
3) Problem Formulation for WGFT: Collecting the two
proposed rate proxies together, we can now formulate theoptimal WGFT problem in the dual graph as an SD problem
as follows:
min ρ
W

s.t.

X

[We (xv1 (e) − xv2 (e) )2 + γ(1 − We )] +e∈V dWe ∈ {1, c}

X|We − Ws |(e,s)∈E d∀ e ∈ V d.
(22)

The problem is an SD problem, because the first and second
term can be interpreted as a deviation cost, while the third term
can be interpreted as a separation cost. We thus employ thealgorithm in [13] to efficiently solve the SD problem in (22).

IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

9

Unlike WGFT, in the UGFT case we do not have closed
form expressions for Rα (x, W) and RT (W) in terms ofW. The reason is that it is difficult to approximate the firstterm in (5)—the rate of UGFT coefficients—using a simple
proxy as done in WGFT. As discussed earlier, the rate proxyin (19) captures only the cost of AC coefficients, ignoring
the cost of DC coefficients, which could be variable for
UGFTs corresponding to graphs with several disconnectedcomponents. Due to this combinatorial nature, we develop a
greedy algorithm based on spectral clustering for the UGFT
search sub-problem, which computes the actual rate via a
divide-and-conquer strategy. We first formulate the UGFT
search problem, and then elaborate on the proposed greedy
algorithm.
1) Problem Reformulation: Since UGFTs are defined on
disconnected and unweighted graphs with the feasible weight
set C = {1, 0}, we modify the constraint in (5) for the UGFT
search problem as follows:
min
W

s.t.

Rα (x, W) + RT (W)Wi,j ∈ {1, 0}∀ i, j ∈ V(23)2) A Greedy Algorithm for UGFT: Instead of exhaustive
search, we develop a greedy algorithm combined with spectralclustering in order to efficiently solve (23). Spectral clustering[35], identification of clusters of similar data, takes advantageof the eigenvectors of the graph Laplacian derived from asimilarity graph of the data. Among the family of spectralclustering algorithms, normalized cuts [14] is a very popularone. We hence employ normalized cuts to identify clusters ofsimilar pixels for UGFT.
For a given pixel block x, the normalized cuts algorithm
consists of four steps: i) construct a similarity graph G s ={V s , E s , Ws } on x, which measures the similarity betweenspixels i and j using Wi,j; ii) acquire the generalized eigenvectors ψ of the graph Laplacian by solving (Ds − Ws )ψ =λDs ψ; iii) bipartition x using the eigenvector with the smallest
non-zero eigenvalue; iv) recursively repartition the segmented
clusters if necessary.
Leveraging on normalized cuts, we develop a greedy algorithm to search for a locally optimal UGFT. The key idea is torecursively partition the pixel block using normalized cuts until
the resulting representation cost increases. For each iteration,
we perform the following three steps on each pixel block x.
5 ”Strongly polynomial time” means that the computation complexity of analgorithm in question is polynomial in the size of the input bits that describe
the instance of the optimization problem. See [39] for details.

which describes the similarity in pixel intensities6 . σw controls
the sensitivity of the similarity measure to the range ofintensity differences. From this graph, we use normalized cuts
to partition the block into two clusters.
Second, we build a graph G u where pixels in the samecluster are connected with weight 1 while those in different
clusters are disconnected. The UGFT is then derived from G u .
Thirdly, we compute the total representation cost as definedin the objective of (23). We encode transform coefficients
via entropy coding to compute the coding rate, and encode
disconnected edges in the graph via AEC to compute thetransform description cost. Note that, instead of using a proxy
as in the search for WGFT, we compute the actual cost oftransform description, since it is more accurate while only
inducing moderate complexity to the greedy algorithm. If thecurrent representation cost is smaller than that of the previous
iteration, then further repartition is performed.
In practice, few iterations are required due to the small
size of the coding block, and hence this greedy algorithm iscomputation-efficient.
C. Illustrative Examples: Optimal GFTs for Different Blocks& at Different Rates
Having elaborated on the problem formulations and algorithms of WGFT and UGFT, we provide examples to illustrate
optimal GFTs for different classes of blocks and at different
rates.

75

60

60

60

75

60

60

60

65

61

61

60

65

61

61

60

75

60

60

75

60

60

60

65

61

61

60

60

65

61

61

60

W=1
W=0.13

250

300

200

200

150
100



B. Unweighted Graph Fourier Transform

First, we construct a similarity graph G s = {V s , E s , Ws }son x. Each weight Wi,jis defined as in [35]:−|xi − xj |2
s
,
(24)
Wi,j
= exp
2
σw



Since the defined deviation term is linear and the separation
term is bi-linear in (22), it is solvable in strongly polynomial
time5 using the algorithm in [13] based on minimum graphcuts. Specifically, the running time of this algorithm for a 4×4
block (adopted in our coding system) is on the order of 10−2
seconds on an Intel Core(TM) i5-4570 CPU 3.20GHz machine(around 2.7 minutes for a 512 × 512 image).

100

0

50
-100
-200
-2

0
0

2



（a）

4

6

8

-50
-2

0

2



4

6

8

（b）

Fig. 7. An illustration of UGFT and WGFT. (a) An example block containing
pixel pairs with strong and zero correlations, and its signal representation in
the UGFT domain. (b) An example block containing pixel pairs with strong
and weak correlations, and its signal representation in the WGFT domain.

We first provide an example in Fig. 7 to show how weemploy the two flavors of GFTs to capture various correlations
6 Note that the similarity graph G s used for normalized cuts is a different
graph than the ones used to define GFTs, and thus the edge weights are not
restricted to be in a small discrete set C.

IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 201510in pixel blocks. Fig. 7(a) shows an example block containing pixel pairs with strong and zero correlations. UGFT ischosen during mode decision based on RD costs, where theblock is divided into two separate partitions and a connected
graph is constructed in each partition. The resulting transform
representation is one DC coefficient for the description ofeach partition, as shown at the bottom of Fig. 7(a). Fig. 7(b)shows an example block containing pixel pairs with strong
and weak correlations. Accordingly WGFT is chosen frommode decision, where a graph containing edge weights c (herec = 0.13) is constructed on the block. The resulting transform
coefficients consist of one large DC term and one small AC
term, shown at the bottom of Fig. 7(b).

40

11

11

11

40

11

11

11

40

11

11

11

40

28

11

11

40

28

11

11

40

28

11

11

40

28

28

11

40

28

28

11

40

28

28

11

40

40

40

40

40

40

40

40

40

40

40

40

storing most popular GFTs in a table for simple lookups.
Graph isomorphism is further exploited to optimize the table.
We first detail the construction of GFT lookup table, and thencompare the complexity of table lookup against that of online
eigen-decomposition.
The table size could be very large if we simply store all theused GFTs. A large table would lead to high requirement instorage and expensive lookup. Hence, we propose to constructa GFT lookup table of relatively small size. Specifically, we
have three strategies:
1) We perform GFT on a LR block (4 × 4), which admits
a smaller space of GFT variants to begin with.
2) Only the most frequently used LR-GFTs are stored.
3) Exploiting graph isomorphism, only one LR-GFT isstored for all isomorphic graphs.
Due to self-similarity in images, the same or similar structures are likely to recur throughout. Hence, the underlying LRGFTs with respect to those structures are frequently used. We
thus store only the most popular LR-GFTs in our lookup table,
while covering a large fraction of the total used LR-GFTs.
We store one GFT for all isomorphic graphs. Two graphs G
and H are isomorphic if there exists such a function mapping
vertices of G to vertices of H:
f : V(G) → V(H),(a) high bitrate(b) medium bitrate(c) low bitrate(25)

that any two vertices i and j of G are adjacent in G if and only
if f (i) and f (j) are adjacent in H. Intuitively, the structuresof isomorphic graphs are the same.

Fig. 8. An illustration of optimal UGFTs at different rates. For a given
pixel block, the figure shows the optimal graph construction (all the edges
connecting pixels are assigned weight 1) and resulting transform coefficientsat (a) high bit rate, (b) medium bit rate and (c) low bit rate.

Next, we provide an illustrative example of optimal UGFTs
for the same pixel block but at different target bit rates. Givena pixel block with three smooth regions as shown in Fig. 8, theminimal representation cost at high bit rates is achieved when
it is divided into three connected components corresponding
to the three smooth regions. This is because this derived
UGFT results in only three DC coefficients, leading to theminimal and significantly smaller cost of quantized transform
coefficients (compared to other UGFTs resulting in manyhigh frequency coefficients) at high bit rates, which together
with the description cost results in the minimal representation
cost. At medium rates, the larger quantization parameter (QP)
quantizes more coefficients to zero, and a simpler UGFT
with similar quantized transform coefficient cost but smallertransform description cost would be more attractive, resulting
in the graph in Fig. 8(b). Finally, at low rates, the large QP
quantizes most coefficients to zero, and the simplest UGFT isbest, as shown in Fig. 8(c).
VII. FAST I MPLEMENTATION OF G RAPH F OURIER
T RANSFORMS
In this section we propose a fast implementation of GFT
for practical deployment. As stated earlier, the online eigendecomposition for the construction of GFT is a hurdle to realtime implementation. We avoid this by pre-computing and

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Fig. 9. An illustration of the 8 underlying isomorphic graphs with respect
to the same graph structure, i.e., the same GFT after appropriate permutationof vertex labels.

For 4-connected graphs constructed on a pixel block, we canconstruct 8 isomorphic graphs with respect to the same graphstructure, thus resulting in the same GFT after appropriate
permutation of vertex labels. See Fig. 9 for an illustration.
The corresponding mapping functions f include reflection androtation of graph vertices. Specifically, a graph can be mapped
to four of its isomorphic graphs by reflection with respect tothe central horizontal, vertical, diagonal and anti-diagonal axes
respectively, and three of its isomorphic graphs by different
degrees of rotation. Hence, we can reduce the table size by up
to a factor of 8 via graph isomorphism. During the encoding
and decoding, the actual GFT, if not directly stored, is retrieved
by mapping from the stored GFT.
With all the above factors considered, the lookup table canbe reduced to a manageable size. Experimentally, we collect
LR-GFTs from ten PWS images with diverse structures toIEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

B. Selection of Transform Modes
We first investigate the selection of transform modes for test
images. Fig. 10 shows an example of the mode decision in
7 http://vision.middlebury.edu/stereo/data/scenes2003/
8 http://www.cvlab.cs.tsukuba.ac.jp/dataset/tsukubastereo.php

Further, we observe that GFTs corresponding to contiguous boundaries are used more frequently, which verifies ourassumption for the proxy of transform description in (21).
C. Objective Comparison in RD Performance
We now compare the performance of the proposed MR-GFT
scheme against SAW, HR-UGFT, HR-DCT and MR-UGFT for
PWS images. Fig. 11 presents the RD performance of theseschemes for four test images with a typical PSNR range. The
proposed MR-GFT achieves significant gain in PSNR over awide range of bit rate. On average we have 6.8 dB gain over
HR-DCT, 5.9 dB gain over HR-UGFT, 2.5 dB gain over SAW,
and 1.2 dB gain over MR-UGFT.
Cones

Teddy

50

45

45
PSNR(dB)

50

40

35
MR−GFT
MR−UGFTSAW
HR−UGFT
HR−DCT

30

25
0.06

0.08

0.1
0.12
0.14
Bit per pixel (bpp)

0.16

40
35

MR−GFT
MR−UGFTSAW
HR−UGFT
HR−DCT

30
25
0.08

0.18

0.1

(a)

0.12
0.14
Bit per pixel (bpp)

0.16

0.18

(b)

Dude

Tsukuba

46

50

44
45

42
40
38
MR−GFT
MR−UGFTSAW
HR−UGFT
HR−DCT

36
34
32
30
0.08

0.1

0.12
0.14
Bit per pixel (bpp)

0.16

0.18

PSNR(dB)To evaluate our proposed MR-GFT coding system for PWS
signals, we use four test image sequences: 448 × 368 depth
maps of Teddy and Cones7 , and graphics images of Dude
(800 × 480) and Tsukuba8 (640 × 480).
For the three transform candidates in the proposed coding
system, the block size of HR-DCT is 8 × 8, and that of LRWGFT and LR-UGFT is 4 × 4. The weighting parameter γ
in (22) is empirically assigned 3 in all the experiments. We
note that our experimental results are not very sensitive to thespecific choice of γ. For the calculation of the fractional edge
weight c, we collect pixel pairs of weak correlation from thetraining images, and compute c to be 0.13 via (17).
We compare coding performance of our proposed scheme
against four compression schemes: H.264 intra (HR-DCT) [2],the GFT coding (HR-UGFT) in [24], the shape-adaptive
wavelet (SAW) in [21], and our previous work MR-UGFT in[34]. Note that intra prediction is used for all schemes.

(a) The original(b) Transform Modes
Fig. 10. The selection of transform modes for a portion of Teddy. Blocks
in red choose WGFT and blocks in blue choose UGFT.

PSNR(dB)

VIII. E XPERIMENTAL R ESULTS
A. Experimental Setup

Teddy. It can be observed that blocks containing pixel pairs
with weak correlation (weak boundaries) in red, generallychoose LR-WGFT, while blocks containing pixel pairs withzero correlation (strong boundaries) in blue, generally choose
LR-UGFT. This verifies our design of WGFT and UGFT
for the representation of blocks containing pixels pairs withdifferent correlations. Note that some blocks with strong /weak boundaries do not choose LR-WGFT / LR-UGFT, e.g.,
the blocks containing black holes. This is because thoseblocks lose much energy during the low-pass filtering if LRWGFT or LR-UGFT is selected, which greatly degrades thereconstruction quality. In these cases, HR-DCT will be chosen
for better preservation of details.

PSNR(dB)build the lookup table. Note that, the training images aredifferent from the testing images in the final experimentation.
Statistical results show that 126 most popular LR-GFTs out ofa total 3749 LR-GFTs cover 88% of actual computed optimal
LR-GFTs. We thus set the table size L to 126 to store the mostpopular LR-GFTs. Further, we encode the table indices by
Huffman Coding [40] based on the used frequency estimatedfrom test images, which is known at both the encoder anddecoder.
Having constructed the lookup table, during encoding wesearch for the GFT (including UGFT and WGFT) for agiven block by finding out the one that gives the best RD
performance. Different from previous optimization (5), we
consider the resulting distortion for a given GFT during tablelookup. Recall that the expected distortion remains the same ina statistical sense when we design GFTs for particular classesof statistical signals. In contrast, during run-time when thedesigned GFTs are fixed in the table, we can consider the
distortion of each block deterministically. We then transmit thetable index losslessly to indicate which LR-GFT is employed
for the given block, so that the decoder is able to identify thecorrect inverse transform.
We now compare the complexity of using eigendecomposition for the derivation of GFT and that of tablelookup. The computational complexity √of eigen-decomposition√at both the encoder and decoder for a N × N (N = 16 in
our case) coding block is O(N 3 ). In contrast, the complexityof table lookup is O(L) (L = 126 in our setting) at theencoder and O(1) at the decoder. Hence, table lookup leads tosignificant reduction in complexity compared against on-line
eigen-decomposition at both the encoder and decoder.

11

40
35

MR−GFT
MR−UGFTSAW
HR−UGFT
HR−DCT

30
25
0.04

0.06

0.08
0.1
0.12
Bit per pixel (bpp)

0.14

0.16

(c)
(d)
Fig. 11. RD performance comparison among different compression schemes
for PWS images.

The gain comes from three improvements: 1) the added edge
weight c, which leads to sparser GFT-domain representation

IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

via proper characterization of weak correlation; 2) minimization of the total representation cost, which results in minimalrate; and 3) the MR scheme, which down-samples each HR
block to LR and thus reducing coding bits.
D. Subjective Comparison(a) HR-DCT

12boundaries among all methods, which validates the effectiveness of the edge-preserving MR-GFT approach. In the images
reconstructed by HR-DCT, boundaries are severely corrupted
since the assumption of the near unity inter-pixel correlation in
the DCT is not valid along strong boundaries in PWS images.
Though GFT is employed in HR-UGFT, some boundaries in
the restored images by HR-UGFT are still corrupted. This isbecause the DCT is chosen in those regions as GFT consumes
more bits due to the boundary coding. Images compressed by
SAW remain sharp along strong boundaries. However, weakboundaries are often blurred. In contrast, MR-GFT preserves
both strong and weak boundaries, since we properly model
weak correlation using the fractional weight c.

(b) HR-UGFT

E. Application to Depth-Image-Based Rendering(c) SAW(d) MR-GFT(e) HR-DCT(f) HR-UGFT(g) SAW(h) MR-GFT

We conduct further experiments on depth maps, which isan important class of PWS images. Instead of being observed
directly, depth maps generally facilitate various end applications, such as virtual view synthesis via depth-image-based
rendering (DIBR). Note that in this case, one can replace thedepth map distortion metric in (4) with a synthesized view
distortion metric. The synthesized view distortion still does notchange using different transforms, because the position error
in the synthesized view is a linear function of the distortion in
the depth map under some assumptions [41]. Hence, the GFT
training remains the same for depth map coding for DIBR.
During run-time, one can replace the depth map distortion
with a synthesized view distortion metric during table lookup.
We use a simple implementation of 3D warping [11] toperform DIBR. Fig. 13 presents the RD curves for Teddy andCones, where the PSNR of synthesized views is evaluated atvarious total rates of stereo depth maps. On average we achieve
2.2dB gain over HR-DCT, 1.8dB gain over HR-UGFT, 1.2dB
gain over SAW, and 1.0dB gain over MR-UGFT. Further, inFig. 14 we show the virtual views of Teddy and Cones synthesized from stereo depth maps compressed using HR-DCT
and MR-GFT at the same bit rate. MR-GFT is observed to
produce more pleasant synthesized images, with fewer ringing
artifacts and corrupted boundaries. The good performance ismostly due to the well-preserved depth map boundaries by
MR-GFT, which plays a critical role in DIBR.
Teddy

31

28.5
28
PSNR(dB)(j) HR-UGFT

29

30
PSNR(dB)(i) HR-DCT

29
MR−GFT
MR−UGFTSAW
HR−UGFT
HR−DCT

28
27
26
0.1

(k) SAW(l) MR-GFT
Fig. 12. The subjective quality comparison among different compression
schemes. (a) ∼ (d) Teddy at 0.10 bpp; (e) ∼ (h) Cones at 0.13 bpp; (i) ∼(l) Dude at 0.13 bpp.

Fig. 12 demonstrates images reconstructed from different
schemes for Teddy, Cones and Dude. We observe that
our MR-GFT produces the cleanest images with the sharpest

Cones

32

0.15

0.2
0.25
Bit per pixel (bpp)

0.3

0.35

27.5
27
MR−GFT
MR−UGFTSAW
HR−UGFT
HR−DCT

26.5
26
25.5

0.2

0.25
0.3
Bit per pixel (bpp)

0.35

0.4

(a)
(b)
Fig. 13. RD performance comparison among different compression schemes
for depth maps tailored for DIBR.

IX. C ONCLUSION
We propose a multi-resolution (MR) graph Fourier transform (GFT) coding scheme for compression of piecewiseIEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015(a) HR-DCT(b) MR-GFT(c) HR-DCT(d) MR-GFT
Fig. 14. Subjective quality comparison of DIBR-synthesized images among
different compression schemes. (a)(b) Teddy; (c)(d) Cones.

smooth (PWS) images. Unlike fixed transforms such as the
DCT, the defined optimal GFT is adaptive to each local blockby minimizing the total representation cost, considering boththe sparsity of the signal’s transform coefficients and thecompactness of transform description. We develop efficientalgorithms to search for optimal GFTs in a defined search
space, based on graph optimization techniques such as spectralclustering and minimum graph cuts. Further, we introduce
two techniques for practical implementation of GFT. One
is the MR scheme where GFT is deployed over a lowpass filtered and down-sampled version of a high-resolution
block. The other is the pre-computation of the most popular
GFTs in a stored table for simple lookup instead of realtime eigen-decomposition. Experimental results show that the
proposed scheme outperforms H.264 intra by 6.8dB in PSNR
on average at the same bit rate. By extension, while we tailor
our proposed approach for PWS image compression, it ispossible to implement our proposal as a coding mode during
compression of general images, so that when a code block isdeemed PWS, our coding scheme can be deployed.

R EFERENCES[1] V. K. Goyal, “Theoretical foundations of transform coding,” in IEEE
Signal Processing Magazine, vol. 18, no.5, September 2001, pp. 9–21.
[2] T. Wiegand, G. Sullivan, G. Bjontegaard, and A. Luthra, “Overview
of the H.264/AVC video coding standard,” in IEEE Transactions on
Circuits and Systems for Video Technology, vol. 13, no.7, July 2003, pp.
560–576.
[3] J. Huang and P. Schultheiss, “Block quantization of correlated gaussian
random variables,” in IEEE Trans. on Comm., vol. 11, no. 3, 1963, p.
289296.
[4] V. Goyal, J. Zhuang, and M. Vetterli, “Transform coding with backward
adaptive updates,” in IEEE Trans. on Inf. Th., vol. 46, no.4, 2000, p.
16231633.
[5] S. Jana and P. Moulin, “Optimality of klt for high-rate transform
coding of gaussian vector-scale mixtures: Application to reconstruction,
estimation and classification,” in IEEE Trans. on Inf. Th., vol. 52, no.9,
2006, pp. 4049–4067.
[6] M. Effros, H. Feng, and K. Zeger, “Suboptimality of the karhunen-loeve
transform for transform coding,” in IEEE Trans. on Inf. Th., vol. 50, no.8,
2004, p. 16051619.
[7] A. K. Jain, Fundamentals of Digital Image Processing. Upper Saddle
River, NJ: Prentice-Hall, 1989.

13[8] H. S. Malvar, A. Hallapuro, M. Karczewicz, and L. Kerofsky, “Lowcomplexity transform and quantization in h.264/avc,” in IEEE Transactions on Circuits and Systems for Video Technology, vol. 13, no.7, July
2003, pp. 598–603.
[9] D. Hammond, P. Vandergheynst, and R. Gribonval, “Wavelets on graphs
via spectral graph theory,” in Elsevier: Appplied and Computational
Harmonic Analysis, vol. 30, April 2010, pp. 129–150.
[10] P. Wan, Y. Feng, G. Cheung, I. Bajic, O. Au, and Y. Ji, “3D motion
in visual saliency modeling,” in IEEE International Conference on
Acoustics, Speech and Signal Processing, Vancouver, Canada, May
2013.
[11] D. Tian, P.-L. Lai, P. Lopez, and C. Gomila, “View synthesis techniquesfor 3D video,” in Applications of Digial Image Processing XXXII,
Proceedings of the SPIE, vol. 7443, (2009), February 2009, pp. 74 430T–
74 430T–11.
[12] S. C. Han and C. I. Podilchuk, “Video compression with dense motion
fields,” in IEEE Transactions on Image Processing, vol. 10, no.11,
November 2001, pp. 1605–1612.
[13] D. S. Hochbaum, “An efficient and effective tool for image segmentation,
total variations and regularization,” in SSVM’11 Proceedings of the Thirdinternational conference on Scale Space and Variational Methods in
Computer Vision, 2011, pp. 338–349.
[14] J. Shi and J. Malik, “Normalized cuts and image segmentation,” in IEEE
Transactions on Pattern Analysis and Machine Intelligence, vol. 22,
no.8, August 2000.
[15] R. Shukla, P. Dragotti, M. N. Do, and M. Vetterli, “Rate-distortion
optimized tree-structured compression algorithms for piecewise polynomial images,” in IEEE Transactions on Image Processing, vol. 14,
no.3, March 2005, pp. 343–359.
[16] Y. Morvan, P. de With, and D. Farin, “Platelets-based coding of depth
maps for the transmission of multiview images,” in SPIE Stereoscopic
Displays and Applications, San Jose, CA, January 2006.
[17] P. Merkle, Y. Morvan, A. Smolic, D. Farin, K. Muller, P. de With,and T. Wiegand, “The effects of multiview depth video compression
on multiview rendering,” in Signal Processing: Image Communication,
vol. 24, 2009, pp. 73–88.
[18] I. Daribo, G. Cheung, and D. Florencio, “Arithmetic edge coding forarbitrarily shaped sub-block motion prediction in depth video compression,” in IEEE International Conference on Image Processing, Orlando,
FL, September 2012.
[19] I. Daribo, D. Florencio, and G. Cheung, “Arbitrarily shaped motion
prediction for depth video compression using arithmetic edge coding,”
in IEEE Transactions on Image Processing, vol. 23, no. 11, November
2014, pp. 4696–4708.
[20] M. B. Wakin, J. K. Romberg, H. Choi, and R. G. Baraniuk, “Waveletdomain approximation and compression of piecewise smooth images,”in IEEE Transactions on Image Processing, vol. 15, no.5, May 2006,
pp. 1071 – 1087.
[21] M. Maitre and M. N. Do, “Depth and depth-color coding using shapeadaptive wavelets,” in Journal of Visual Communication and Image
Representation, vol. 21, July 2010, pp. 513–522.
[22] E. J. Candes and D. L. Donoho, “New tight frames of curvelets andoptimal representations of objects with piecewise-c2 singularities,” inComm. Pure and Appl. Math., 2002, pp. 219–266.
[23] M. N. Do and M. Vetterli, “The contourlet transform: an efficient
directional multiresolution image representation,” in IEEE Transactions
on Image Processing, vol. 14, no.12, December 2005, pp. 2091–2106.
[24] G. Shen, W.-S. Kim, S. Narang, A. Ortega, J. Lee, and H. Wey, “Edgeadaptive transforms for efficient depth map coding,” in IEEE Picture
Coding Symposium, Nagoya, Japan, December 2010.
[25] W. Hu, X. Li, G. Cheung, and O. Au, “Depth map denoising using graphbased transform and group sparsity,” in IEEE International Workshop on
Multimedia Signal Processing, Pula, Italy, October 2013.
[26] C. Zhang and D. Florencio, “Anaylzing the optimality of predictive
transform coding using graph-based models,” in IEEE Signal Processing
Letters, vol. 20, no. 1, January 2013, pp. 106–109.
[27] W.-S. Kim, A. Ortega, P. Lai, D. Tian, and C. Gomila, “Depth map
coding with distortion estimation of rendered view,” in SPIE Visual
Information Processing and Communication, San Jose, CA, January
2010.
[28] N.-M. Cheung, D. Tian, A. Vetro, and H. Sun, “On modeling therendering error in 3d video,” in IEEE International Conference on Image
Processing, Orlando, FL, September 2012.
[29] G. Cheung, A. Kubota, and A. Ortega, “Sparse representation of
depth maps for efficient transform coding,” in IEEE Picture Coding
Symposium, Nagoya, Japan, December 2010.

IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 201514[30] G. Cheung, J. Ishida, A. Kubota, and A. Ortega, “Transform domain
sparsification of depth maps using iterative quadratic programming,” inIEEE International Conference on Image Processing, Brussels, Belgium,
September 2011.
[31] S. Yea and A. Vetro, “View synthesis prediction for multiview videocoding,” in Signal Process., Image Commun., vol. 24, no.1-2, January
2009, pp. 89–100.
[32] D. B. Graziosi, N. M. M. Rodrigues, S. M. M. de Faria, D. Tian, and
A. Vetro, “Analysis of depth map resampling filters for depth-based
3d video coding,” in International Conference on Telecommunications,Casablanca, May 2013.
[33] G. Shen, W.-S. Kim, A. Ortega, J. Lee, and H. Wey, “Edge-aware intra
prediction for depth-map coding,” in IEEE International Conference on
Image Processing, Hong Kong, September 2010.
[34] W. Hu, G. Cheung, X. Li, and O. Au, “Depth map compression using
multi-resolution graph-based transform for depth-image-based rendering,” in IEEE International Conference on Image Processing, Orlando,
FL, September 2012.
[35] U. Luxburg, “A tutorial on spectral clustering,” in Statistics and Computing, vol. 17, no.4, December 2007, pp. 395–416.
[36] R. M. Gray and D. L. Neuhoff, “Quantization,” in IEEE Transactions
on Information Theory, vol. 44, no.6, October 1998.
[37] D. Geman and C. Yang, “Nonlinear image recovery with half-quadratic
regularization,” in IEEE Transactions on Image Processing, vol. 4, no.7,
July 1995, pp. 923–945.
[38] W. S. Kim, S. K. Narang, and A. Ortega, “Graph based transforms for
depth video coding,” in IEEE International Conference on Acoustics,
Speech, and Signal Processing, Kyoto, Japan, March 2012.
[39] K. S. C. Papadimitriou, “Combinatorial optimization: Algorithm andcomplexity,” in Dover Publications, Inc., 1998, p. 165.
[40] D. A. Huffman, “A method for the construction of minimum-redundancy
codes,” in Proceedings of the IRE, vol. 40, no. 9, 1952, pp. 1098–1101.
[41] W.-S. Kim, A. Ortega, P. Lai, D. Tian, and C. Gomila, “Depth map
distortion analysis for view rendering and depth coding,” in IEEE
International Conf. on Image Processing, Cairo, Egypt, November 2009.

Gene Cheung (M’00—SM’07) received the B.S.
degree in electrical engineering from Cornell University in 1995, and the M.S. and Ph.D. degrees inelectrical engineering and computer science from the
University of California, Berkeley, in 1998 and 2000,
respectively.
He was a senior researcher in Hewlett-Packard
Laboratories Japan, Tokyo, from 2000 till 2009. He
is now an associate professor in National Instituteof Informatics in Tokyo, Japan.
His research interests include image & video
representation, immersive visual communication and graph signal processing.
He has published over 140 international conference and journal publications.
He has served as associate editor for IEEE Transactions on Multimedia(2007–2011) and currently serves as associate editor for DSP Applications
Column in IEEE Signal Processing Magazine, APSIPA journal on signal &information processing and SPIE Journal of Electronic Imaging, and as area
editor for EURASIP Signal Processing: Image Communication. He serves as
the lead guest editor of the special issue on “Interactive Media Processing
for Immersive Communication” in IEEE Journal on Special Topics on Signal
Processing, target published date March 2015. He currently serves as memberof the Multimedia Signal Processing Technical Committee (MMSP-TC) in
IEEE Signal Processing Society (2012-2014). He has also served as area chair
in IEEE International Conference on Image Processing (ICIP) 2010, 20122013, technical program co-chair of International Packet Video Workshop(PV) 2010, track co-chair for Multimedia Signal Processing track in IEEE
International Conference on Multimedia and Expo (ICME) 2011, symposium
co-chair for CSSMA Symposium in IEEE GLOBECOM 2012, and area chairfor ICME 2013. He was invited as plenary speaker for IEEE International
Workshop on Multimedia Signal Processing (MMSP) 2013 on the topic “3Dvisual communication: media representation, transport and rendering”. He isa co-author of best student paper award in IEEE Workshop on Streaming andMedia Communications 2011 (in conjunction with ICME 2011), best paper
finalists in ICME 2011 and ICIP 2011, best paper runner-up award in ICME
2012, and best student paper award in ICIP 2013.

Wei Hu received the B.S. degree in electrical engineering from the University of Science and Technology of China, Hefei, China, in 2010. She iscurrently pursuing the Ph.D. degree in electronic andcomputer engineering from The Hong Kong University of Science and Technology, Kowloon, HongKong. Her research interests include graph signal
processing, 2D / 3D image & video representation,coding and processing.

Antonio Ortega (S’91–M’95–SM’00–F’07) received the Telecommunications Engineering degreein 1989 from the Universidad Politecnica de Madrid,
Madrid, Spain, and the Ph.D. degree in electricalengineering in 1994 from Columbia University, New
York, NY, where he was supported by a Fulbrightscholarship.
Since 1994, he has been with the Department
of Electrical Engineering Systems, University ofSouthern California (USC), Los Angeles, where heis currently a Professor. He was a Director of the
Signal and Image Processing Institute and is currently the Associate Chair ofElectrical Engineering Systems with the USC. His work with the USC has
been or is being funded by agencies such as the National Science Foundation,
the National Aeronautics and Space Administration, and the DOE, as well
as a number of companies. Over 25 Ph.D. students have completed their
Ph.D. thesis work under his supervision at the USC, and his work has led toover 250 publications in international conferences and journals. His researchinterests are in the areas of multimedia compression, communications, andsignal analysis. His recent work is focusing on distributed compression,
multiview coding, error tolerant compression, wavelet-based signal analysis,
and information representation in wireless sensor networks.
Dr. Ortega is a member of the Association for Computing Machinery. He
has been the Chair of the Image and Multidimensional Signal Processing
Technical Committee and a member of the Board of Governors of the IEEE
Signal Processing Society (2002). He has been Technical Program Co-chairof the IEEE International Conference on Image Processing 2008, the IEEE
International Workshop on Multimedia Signal Processing 1998, and the IEEE
International Conference on Multimedia and Expo 2002. He has been the
Associate Editor for the IEEE Transactions On Image Processing, the IEEE
Signal Processing Letters, and the EURASIP Journal on Advances in Signal
Processing. He received the National Science Foundation CAREER award,
the 1997 IEEE Communications Society Leonard G. Abraham Prize Paper
Award, the IEEE Signal Processing Society 1999 Magazine Award, and the2006 EURASIP Journal of Advances in Signal Processing Best Paper Award.

IEEE TRANSACTIONS ON IMAGE PROCESSING, JANUARY 2015

Oscar C. Au received his B.A.Sc. from Univ. of
Toronto in 1986, his M.A. and Ph.D. from Princeton
Univ. in 1988 and 1991 respectively. After being
a postdoc in Princeton for 1 year, he joined Hong
Kong Univ. of Science and Technology (HKUST) asan Assistant Professor in 1992. He is/was a Professor of Dept. of Electronic and Computer Engineering, Director of Multimedia Technology ResearchCenter, and Director of Computer Engineering in
HKUST. He is a core member of the State Key Lab
on Advanced Displays and Optoelectronic Technology.
His main research contributions are on video/image coding and processing,watermarking/light weight encryption, speech/audio processing. Research topics include fast motion estimation for H.261/3/4/5, MPEG-1/2/4, and AVS, optimal and fast sub-optimal rate control, mode decision, transcoding, denoising,
deinterlacing, post-processing, multi-view coding, view interpolation, depth
estimation, 3DTV, scalable video coding, distributed video coding, subpixel
rendering, JPEG/JPEG2000, HDR imaging, compressive sensing, halftone
image data hiding, GPU-processing, software-hardware co-design, etc. He
has published 70 technical journal papers, 370+ conference papers, 3 bookchapters, and 70+ contributions to international standards. His fast motion
estimation algorithms were accepted into the ISO/IEC 14496-7 MPEG-4
international video coding standard and the China AVS-M standard. His lightweight encryption and error resilience algorithms are accepted into the China
AVS standard. He was Chair of Screen Content Coding AdHoc Group in
JCTVC for HEVC. He has 25+ granted US patents and is applying for
70+ more on his signal processing techniques. He has performed forensic
investigation and stood as an expert witness in Hong Kong courts many times.
Dr. Au is a Fellow of IEEE and HKIE. He is/was Associate/Senior Editorsof several IEEE journals (TCSVT, TIP, TCAS1, SPL, JSTSP, SigView) andnon-IEEE journals (JVCIR, JSPS, TSIP, JMM, JFI, and SWJ). He is guest
editor of some special issues in JSTSP and TCSVT. He is/was BoG member
and Vice President Technical Activity of APSIPA. He is/was Chair of 3
technical committees: IEEE CAS MSA TC, IEEE SPS MMSP TC, andAPSIPA IVM TC. He is a member of 5 other TCs: IEEE CAS VSPC TC,
DSP TC, IEEE SPS IVMSP TC, IFS TC, and IEEE ComSoc MMC TC.
He served on 2 steering committees: IEEE TMM, and IEEE ICME. He also
served on organizing committee of many conferences including ISCAS 1997,
ICASSP 2003, ISO/IEC 71st MPEG in Jan 2005, ICIP 2010, etc. He was/willbe General Chair of several conferences: PCM 2007, ICME 2010, PV 2010,
MMSP 2015, APSIPA ASC 2015, and ICME 2017. He won 5 best paper
awards: SiPS 2007, PCM 2007, MMSP 2012, ICIP 2013, and MMSP 2013.
He was IEEE Distinguished Lecturer (DL) in 2009 and 2010, APSIPA DL in2013 and 2014, and has been keynote speaker multiple times.

15