Gene Cheung, Xianming Liu
National Institute of Informatics
11th July, 2016

Graph Signal Processing for Image
Compression & Restoration(Part II).
ICME'16 Tutorial 07/11/2016

1

Outline (Part II)• Image Restoration using GSP Tools• Image Denoising• Soft Decoding of JPEG Compressed Images

ICME'16 Tutorial 07/11/2016

2

Outline (Part II)• Image Restoration using GSP Tools• Image Denoising• Sparsity Prior• Smoothness Prior• Soft Decoding of JPEG Compressed Images

ICME'16 Tutorial 07/11/2016

3

Introduction to PWS Image Denoising•

Limitations of current sensing technologies- acquired PWS images are often corrupted by non-negligible acquisition noise.

•

Denoising is an inverse imaging problem.
observation•y  xv

desired signal
noise

Signal prior is key to inverse imaging problems!
•

Depth images are PWS, self-similar.
ICME'16 Tutorial 07/11/2016

4

Existing Image Denoising Methods•

Local methods (e.g., bilateral filtering)•

Nonlocal image denoising
Buades et al, "A non-local algorithm for image denoising," CVPR 2005
- Assumption: nonlocal self-similarity•

Dictionary learning based
Elad et al, "Image denoising via sparse and redundant representation overlearneddictionaries," TIP 2006.

- represent a signal by the linear combination of a few atoms out of a dictionary
Other related works- Huhle et al, “Robust non-local denoising of colored depth data,” CVPR Workshop
2008- Tallon et al, “Upsampling and denoising of depth maps via joint segmentation,”
EUSIPCO 2012
ICME'16 Tutorial 07/11/2016

5Key Idea in Non-local GFT

Nonlocal self-similarity Local Piecewise Smoothness

unify in GFT domain
Challenges
Our method
1. Adapt to nonlocal statistics --- adapt to nonlocal statistics via nonlocal self-similarity
2. Characterize PWS--- characterize PWS via GFT representation+ learn GFT dictionary efficiently

ICME'16 Tutorial 07/11/2016

6

NL-GFT Algorithm
common GFT from avg. patch
observation iNNmin  yi  U i 2     i 0
U ,

i 1

Algorithm:

2

i 1

code vector for observation i

1.

Identify similar patches, compute avg patch.
（self-similarity）

2.

Given avg patch, use Gaussian kernel tocompute weights between adjacent pixels.

3.

Compute graph Fourier transform (GFT).

4.

Given GFT, soft thresholding on transform coeff.
for sparse representation.

[1] W. Hu, X. Li, G. Cheung, O. Au, "Depth Map Denoising using Graph-based Transform and Group Sparsity," IEEE International
Workshop on Multimedia Signal Processing, Pula (Sardinia), Italy, October, 2013. (Top 10% paper award.)

7

Justification of Sparsity Prior• GFT domain sparsity prior in objective function:
KKmin  yi  xi 2    xi 0
 , xi

i 1

2

i 1

• ”Argument”:• GFT approximates KLT if statistical model is GMRF and each graphweight captures correlation of 2 connected pixels [2, 3].
• Underlying “causes” of PWS signals are few; PWS signal can besparsely represented in GFT domain [4, 5].
[2] 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.
[3] W. Hu, G. Cheung, A. Ortega, O. Au, “Multi-resolution Graph Fourier Transform for Compression of Piecewise Smooth Images,”
IEEE Transactions on Image Processing, January 2015.
[4] G. Shen, W.-S. Kim, S.K. Narang, A. Ortega, J. Lee, and H. Wey, “Edge-adaptive transforms for efﬁcient depth map coding,” in IEEE
Picture Coding Symposium, Nagoya, Japan, December 2010.
[5] W. Hu, G. Cheung, X. Li, O. Au, “Depth Map Compression using Multi-resolution Graph-based Transform for Depth-image-based
Rendering,” IEEE International Conference on Image Processing, Orlando, FL, September 2012.

8

Experimental Results (1)• Setup:- Test Middleburry depth maps: Cones, Teddy, Sawtooth- Add Additive White Gaussian Noise- Compare against Bilateral Filtering (BF), Non-Local Means Denoising (NLM)
and Block-Matching 3D (BM3D)• Results– Up to 2.28dB improvement over BM3D.

ICME'16 Tutorial 07/11/2016

NLGFT

BM3D

NLM

BF
9

Experimental Results (2)• Setup:- Test Middleburry depth maps: Cones, Teddy, Sawtooth- Add Additive White Gaussian Noise- Compare agaist Bilateral Filtering (BF), Non-Local Means Denoising(NLM) and Block-Matching 3D (BM3D)• Results– Up to 2.28dB improvement over BM3D.

ICME'16 Tutorial 07/11/2016

NLGFT

BM3D

NLM

BF

10

10

Outline (Part II)• Image Restoration using GSP Tools• Image Denoising• Sparsity Prior• Smoothness Prior• Soft Decoding of JPEG Compressed Images• Joint Denoising / Contrast Enhancement

ICME'16 Tutorial 07/11/2016

11Motivation (I)• Image denoising—a basic restoration problem:
observationnoisez  uedesired signal• It is under-determined, needs image priors for regularization:
fidelity term

min z  u 2   prior(u)2prior termu

• Graph Laplacian regularizer: should be small for target patch u

SG (u)  u T Lu

L  DA
graph Laplacian matrix• Many works use Gaussian kernel to compute graph weights [1, 6]:  dist i, j 2 
wij  exp
2





dist(i, j ) is some distance metric between pixels i and j[6] D. Shuman et al., “The emerging field of signal processing on graphs: extending high-dimensional data analysis to networks andother irregular domains,” IEEE Signal Processing Magazine, vol. 30, no. 3, pp. 83–98, 2013.

12

Motivation (II)• However…
a.

Why is SG ( x)  u T Lu a good prior?

Tub. How to design the optimal Lu for restoration?

Tc. u Lu performs particularly well on PWS images, why?
• We answer these basic questions by viewing:•

discrete graph as samples of high-dimensional manifold.

approximate

discrete graphcontinuous manifold[7] Jiahao Pang, Gene Cheung, Antonio Ortega, Oscar C. Au, "Optimal Graph Laplacian Regularization for Natural Image Denoising," IEEE
International Conference on Acoustics, Speech and Signal Processing, Brisbane, Australia, April, 2015.
[8] Jiahao Pang and Gene Cheung, “Graph Laplacian Regularization for Inverse Imaging: Analysis in the Continuous Domain,” arXiv 13
preprint, arXiv:1604.07948, 2016.

Our Contributions
1.

We show SG (u)  u T Lu converges to a continuous functional S ,T
analysis of S explains the mechanism of u Lu for inverse imaging
Graph Laplacian
regularizer S

convergeG

2.

We derive the optimal graph Laplacian regularizer for denoising, which is
discriminant for small noise and robust when very noisy.
Non-local self-similarity and
MMSE formulation

3.

A continuous functionalS for regularization

obtainOptimal regularizer SG

We interpret graph Laplacian regularization as anisotropic diffusion, showthat it not only smooths but may also sharpens the image, promote
piecewise smooth images

ICME'16 Tutorial 07/11/2016

14

Graph-Based Image Processing• Graph for image restoration• Each pixel corresponds to a vertex in a graph (denote # of pixels as M ).

e.g., graph of a 5×5 patch,(not necessarily be a grid graph)• Regard the image as a signal defined on a weighted graph.
• With proper graph configuration, construct filter for image (graph signal)using prior knowledge (i.e., smooth on the graph).

15

Road Map
Continuous Domain

Discrete Domain

Choose the continuous SAMPLE
exemplar functions{ f n }nN1

Sample { f n }nN1 to obtainN
the discrete {f n }n1

Get metric space G  R 22
on point-by-point basis

Compute the weights
and Laplacian L  R M M

Obtain continuousfunctional S (u )CONVERGE

Graph Laplacian
Regularizer SG (u)• Different exemplars { f n }nN1 lead to different regularization behavior!
ICME'16 Tutorial 07/11/2016

16

Graph Construction (I)Roadmap• First, define:• 2D domain   R 2—shape of an image patch•   si  [ xi yi ] | si  ,1  i  M T— M uniformly distributed
random samples on  ,pixel locations in our work

Exemplars{ f n }nN1

SamplesN{f n }n1MatrixGraph weights,
and L  R M M

G  R 22

converge

FunctionalS (u )

Regularizer
SG (u)• (Freely) choose N continuous functions
f n ( x, y ) :   R , 1  n  Ncalled exemplar functions, for example• intensity for gray-scale image ( N  1 )• R, G, B channels for color image ( N  3 )ICME'16 Tutorial 07/11/2016

17

Graph Construction (II)• Sampling f n at positions in  gives N
discretized exemplar functions
f n  [ f n ( x1 , y1 ) f n ( x2 , y2 )  f n ( xM , yM )]T

Roadmap
Exemplars{ f n }nN1

SamplesN{f n }n1MatrixGraph weights,
and L  R M M

G  R 22

converge

FunctionalS (u )

Regularizer
SG (u)• For each sample si  , define a length N vectorv i  [f1 (i ) f 2 (i )  f N (i )]T• Build a graph G with M vertices; each sample si   has avertex Vi
ICME'16 Tutorial 07/11/201618

Graph Construction (III)• Weight between vertices Vi and V jdegree before normalization

i   j 1 (dij )Roadmap
Exemplars{ f n }nN1

SamplesN{f n }n1MatrixGraph weights,
and L  R M M

M

wij  ( i  j )   ( dij )normalization factor Thresholded Gaussian kernel d2 
exp   2  d  r , (d )   2 
0otherwiseG  R 22

converge

Functional

Regularizer
SG (u)S (u )“Distance” between two vertices
dij2  v i  v j• G is an r-neighborhood graph:no edge connecting two vertices with distance greater than r

2
2

19

Graph Construction (IV)
Roadmap• Our graph G is very general• e.g., one can derive that the popular
2D grid graph is a special case of ours

Exemplars{ f n }nN1

SamplesN{f n }n1MatrixGraph weights,
and L  R M M

GR

22

converge

FunctionalS (u )

• A — (i, j ) -th entry is wijm—diagonal
entryisD j 1 wij

Regularizer
SG (u)

unnormalized Graph LaplacianL  DA• u ( x, y ) :   R is a continuous imageu  [u ( x1 , y1 ) u ( x2 , y2 )  u ( xM , yM )]T — discrete version of u ( x, y )• SG (u)  u T Lu — graph Laplacian
regularizer
ICME'16 Tutorial 07/11/2016

20

Convergence of the Graph Laplacian Regularizer (I)Roadmap• The continuous counterpart of SG is afunctional S for image on domain S (u )   u G u
T

1



 det G 

Exemplars{ f n }nN1

SamplesN{f n }n1MatrixGraph weights,
and L  R M M

2 1

ds

G  R 22

u  [ xu  y u ]T is the gradient of uconverge

FunctionalS (u )

• G is a 2-by-2 matrix-valued function:  N   x f n 2
n 1
G N

  n 1  x f n · y f n  f · f 
  f ·f   f  Regularizer
SG (u)

N

n 1

x n

2

N

n 1

y n

y n

N

n 1

n

T
n

Structure tensor [9] of theNgradients {f n ( x, y )}n 1N• G is computed from {f n }n 1 on a point-by-point basis[9] H. Knutsson, C.-F. Westin, and M. Andersson, “Representing local structure using tensors ii,”
in Image Analysis. Springer, 2011, vol. 6688, pp. 545–556.

21

Convergence of the Graph Laplacian Regularizer (II)Roadmap• Theorem : convergence of SG to S

lim SG  u  ~ S  u 

M 
r 0

Exemplars{ f n }nN1

SamplesN{f n }n1MatrixGraph weights,
and L  R M M

G  R 221) number of samples M increases2) neighborhood r shrinks

converge

Functional“~” means there exist a constantsuch that equality holds.

S (u )

Regularizer
SG (u)• With results of [10], we proved it by viewing a graph as proxy of an
N -dimensional Riemannian manifold
Vertex

Coordinate on Ω

Coordinate on N-D manifoldVisi   xi , yi 

v i  [f1 (i ) f 2 (i )  f N (i )]TICME'16 Tutorial 07/11/2016[10] M. Hein, “Uniform convergence of adaptive graph-based regularization,” in Learning Theory. Springer, 2006, pp. 50–64.

22Interpretation of Graph Laplacian Regularizer (I)• SG converges to S , with S , any new insights we gain on SG ??
• Inspect the equations carefully…

S (u )   u G u
T



1

 det G 

2 1

ds

N

G   f n ·f n Tn 1

SG (u)  u T Lu

• 3 observations:•u T G 1u measures length of u in a metric space built by G !

• The eigen-space of G reflects dominant directions of {f n }nN1• S integrates the norm of gradient
ICME'16 Tutorial 07/11/2016

23

Justification of Graph Laplacian Regularizer (II)• Metric space defined by G ?

• At a certain location ( x, y ) on the image

𝜕𝑦Ellipses are contours (isolines),N
reflects how concentrate {f n }n 1

Green dots are {f n }n 1
N

𝑙

𝑂

𝜕𝑥
l: dominant direction,
eigenvector corresponds tothe largest eigenvalue of GS (u )   u G u
T



1



det G



2 1

N

ds

ICME'16 Tutorial 07/11/2016G   f n ·f n T
n 1

24

Justification of Graph Laplacian Regularizer (III)• The 2D metric space provides a clear picture of what signals arebeing discriminated and to what extent, on a point-by-point basis in
the continuous domain.
𝜕𝑦

𝑂

𝜕𝑥
𝑙

𝑙

(a)

𝜕𝑦

𝜕𝑦

𝜕𝑥

ground-truth

𝜕𝑥

𝑂

𝑂
𝑙

(b)

(c)

• Both (a)(b) are correct, but (b) is more discriminant,(c) is discriminant but incorrect• Lesson: when ground-truth is unknown, one should design a
discriminant metric space only to the extent that estimates ofground-truth are reliable!
ICME'16 Tutorial 07/11/2016

25

Noise Modeling in Gradient Domain• For a M  M noisy patch z 0  R , identify K  1 similar patcheson the noisy image, the K patches {z k }kK01 form a cluster
M(i )• On patch z k , gradient at pixel i is g k .
K 1
{
g
}
• Drop superscript i , model the noisy gradients k k 0 as

g k  g  ek , 0  k  K  1
Unknown ground-truth

Noise term, follows 2D Gaussian2with zero-mean and covariance  g I• PDF of g k given ground-truth g (likelihood) is simply

 12
Pr (g k | g ) 
exp   2 g  g k 2 
2
 2

2 g
g


1

ICME'16 Tutorial 07/11/2016

26

Seeking for the Optimal Metric Space (I)g𝜕𝑦• We first establish an ideal metric space assuming weknow ground truth: g𝜕𝑥

G I (g )  gg T   I

𝑙𝑂

It is discriminant to g  0 , smaller  makes the space more skewed
K 1• With noisy gradients {g k }k 0 seek for the optimal metric spaceG  arg min  2 G  G I g  F  Pr g | g 
2

*

R

G



K 1
k k 0

G *   2 G I g   Pr g | g k k 0 dg
R

ICME'16 Tutorial 07/11/2016

K 1

dg

(1)
27

Seeking for the Optimal Metric Space (II)• Assume the prior Pr (g ) is a 2D Gaussian with covariance  p2 I , then derivePr g | g k k 0

K 1



2
 1

exp   2 g  g  2 
2
2
 2


1

where the “ensemble” mean g and variance  2 are

g

1


 
K    p2g2  k 0 k
K 1

2
g
2
g

2

K   g2  pnoise variance of gk• Carrying out the integral in (1) gives the optimal metric spaceTG  gg  ( 2   )I
・

(2)

・
2
• Intuition: If noise  is small, gg dominates and G is discriminant;2
if  is large, ( 2   )I dominates, G ・ defaults to Euclidean space!
T

ICME'16 Tutorial 07/11/2016

28

From Metric Space to Graph LaplacianT• The structure of G  gg  ( 2   )I allows us to assign N  3
exemplar functions, such that they lead to the optimal metric space:・

f 2・ (i )   2   ·yi

f1・ (i )   2   ·xi
f3 ・ 

1
K  
2
g

2
p

K 1
k 0

zk

— Spatial— Intensity• f1・ (i ) and f 2・ (i ) correspond to the term ( 2   )I in G ・T・3• f (i ) leads to the term gg in G ・• Withf ・ 3
i i1 , compute the optimal graph Laplacian

ICME'16 Tutorial 07/11/2016

29

Graph Laplacian Regularization as Anisotropic Diffusion• Our denoising problem is u  arg min z  u 2   ·u T Lu
2・u

• Continuous counterpart is u ・  arg min u  z0    · u T Duds
2

where

DG

1

 det G 

u



2 1

• Differentiate with respect to u and equate it to zerou*  z0   div D u* • Equivalent to marching z0 forward in time using tensor diffusion

 t u  div  Du  , u (s, t  0)  z0 (s)with step size  and tensor diffusivity D

ICME'16 Tutorial 07/11/2016

30

Graph Laplacian Regularization as P-M Diffusion (I)T• Assuming small noise, with G  gg  ( 2   )I , tensor diffusion is
simplified as・ t u  G2 div   ( u 2 )u 

where diffusivity is a scalar2u 2  ( u 2 )   1  2 

G 


 1.5and constantG   2  

It is called the Perona-Malik diffusion• Denote J ( u 2 )   ( u 2 ) u 2 , diffusion scheme rewritten as

 t u  G2   ( u 2 ) u  J ( u 2 ) u 
decomposed as two independent diffusion processes.
ICME'16 Tutorial 07/11/2016

31

Graph Laplacian Regularization as P-M Diffusion (II)• A closer look..

 : direction parallel to u : direction perpendicular to u
Denote a diffusion
process along edges

Denote a diffusion
process across edges

 t u  G2   ( u 2 ) u  J ( u 2 ) u 
Diffusivity along edges

Diffusivity across edges• Along edges: forward diffusion, i.e., smoothing, as2
u 2  ( u 2 )   1  2 

G 


 1.5

0

ICME'16 Tutorial 07/11/2016

32

Graph Laplacian Regularization as P-M Diffusion (III)• Across edges: forward-backward diffusion, define T  G

2(1   )

u 2  T  J   0 , forward diffusion
If   1

u 2  T  J   0 , no diffusion
u 2  T  J   0 , backward diffusion (sharpening)If   1, becomes a discretization of TV

If   1, forward diffusion to smooth the image• Not only smooth but also sharpen images→PWS2u 2 
Diffusivity across edges J1( u )  1  2 
2

G 



  2.5

 2 u 2

2

·1 
(


1)


G2



 t u  G2   ( u 2 ) u  J ( u 2 ) u 
Denote a diffusion process across edges
ICME'16 Tutorial 07/11/2016

33

Experimentation (I)• We develop an iterative patch-based algorithm• Optimal Graph Laplacian Regularization (OGLR) for denoising• Step 1: Search for similar patches• Step 2: Compute the optimal graph Laplacian• Step 3: Patch-based denoising• Step 4: Denoised image aggregation• Corruption model: i.i.d. Additive White Gaussian Noise (AWGN)• Measurements: PSNR (in dB), SSIM
• Natural images..

• Test images: Lena, Barbara, Peppers, Mandrill• Compared with state-of-the-arts: K-SVD, BM3D, PLOW
ICME'16 Tutorial 07/11/201634

Experimentation (II)• Objective results

K-SVD: PSNR

BM3D: PSNR

K-SVD: SSIM

BM3D: SSIM

PLOW: PSNROGLR: PSNR

PLOW: SSIM

OGLR: SSIM

0.3 dB better than BM3D!
ICME'16 Tutorial 07/11/2016

35

Experimentation (III)• Subjective comparisons (  I  40 )Original

Noisy, 16.48 dB

K-SVD, 26.84 dB

BM3D, 27.99 dB

PLOW, 28.11 dB

OGLR, 28.35 dB

ICME'16 Tutorial 07/11/201636

Experimentation (IV)• Piecewise smooth images..
• Test images: Cones, Teddy, Art, Moebius, Aloe
• Compared with state-of-the-arts: BM3D, NLGBT• Objective results

BM3D: PSNRNLGBT: PSNROGLR: PSNR

BM3D: SSIM

NLGBT: SSIM

OGLR: SSIM

1.6 dB better than NLGBT!

37

Experimentation (V)• Subjective comparisons (  I  30 )Original

Noisy, 18.66 dB

BM3D, 33.26 dB NLGBT, 33.41dB OGLR, 34.32 dB
ICME'16 Tutorial 07/11/2016

38

Summary of Image Denoising via Graph Smoothness Prior• Inverse imaging problems are ill-posed; we use graph Laplacian regularizer
as image prior• Graph Laplacian regularizer converges to a continuous functional, analysisof the functional explains the mechanisms and implications of graphLaplacian regularizer• We describe a methodology to derive the optimal edge weights given
nonlocal noisy gradient observations• By interpreting graph Laplacian regularization as anisotropic diffusion, weshow that it not only smooth images but may also sharpen images,
promoting piecewise smooth results• Our algorithm performs competitively with state-of-the-art methods fornatural images, and out-perform them for piecewise smooth images

ICME'16 Tutorial 07/11/2016

39

Outline (Part II)• Image Restoration using GSP Tools• Image Denoising• Sparsity Prior• Smoothness Prior• Soft Decoding of JPEG Compressed Images

ICME'16 Tutorial 07/11/201640