Gene Cheung
National Institute of Informatics
6th July, 2015

Depth Image Coding & Processing
Part 3: Depth Image Processing.
COST Training School 7/06/2015

1

Outline• Depth Image Denoising• Graph Sparsity Prior• Graph-signal Smoothness Prior• Bit-depth Enhancement

COST Training School 7/06/2015

2

Outline• Depth Image Denoising• Graph Sparsity Prior• Graph-signal Smoothness Prior• Bit-depth Enhancement

COST Training School 7/06/2015

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.
COST Training School 7/06/2015

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
COST Training School 7/06/2015

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

COST Training School 7/06/2015

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.

COST Training School 7/06/2015

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.

COST Training School 7/06/2015

NLGFT

BM3D

NLM

BF

10

10

Outline• Depth Image Denoising• Graph Sparsity Prior• Graph-signal Smoothness Prior• Bit-depth Enhancement

COST Training School 7/06/2015

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

min y  x 2   prior(x)2prior termx

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

SG (x)  xT Lx

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 dist  i, j 2 
wij  exp 

2






Motivation (II)• However…
a.

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

b.

Why using Gaussian kernel for edge weights?

c.

How to design a discriminant x Lx for restoration?

T• 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, Gene Cheung, Wei Hu, Oscar C. Au, "Redefining Self-Similarity in Natural Images for Denoising Using Graph Signal 13
Gradient," APSIPA ASC, Siem Reap, Cambodia, December, 2014.

 dist  i, j 2 
wij  exp 

2






Our Contributions
1.

Using Gaussian kernel to compute graph weights, SG (x)  x T Lx
converges to a continuous functional S .
Graph Laplacian
regularizer S

convergeG

2.

Analysis of functional S provides understanding of how signals arebeing discriminated and to what extent; careful graph construction leads
to discriminant signal prior.
Select edge weights via“feature functions”

3.

A continuous functionalS for regularization

obtainEffective regularizer SG

We derive the optimal graph Laplacian regularizer for denoising, which is
discriminant for small noise and robust when very noisy.
COST Training School 7/06/2015

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 continuousNfeature functions{ f n }n 1Sample { f n }nN1 to obtainN
the discrete {f nD }n1

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

Obtain continuousfunctional S (h)Compute the weights
and Laplacian L  R M M

CONVERGE

Graph Laplacian
Regularizer SG (h D )• Different features { f n }nN1 lead to different regularization behavior!
COST Training School 7/06/2015

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

Features{ f n }nN1

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

G  R 22

converge

Functional

Regularizer
SG (h D )S  ( h)• (Freely) choose N continuous functions
f n ( x, y ) :   R , 1  n  N𝑦𝑓𝑛 (𝑥, 𝑦)called feature functions, for example• intensity for gray-scale image ( N  1 )• R, G, B channels for color image ( N  3 )COST Training School 7/06/2015

Ω
𝑂
𝑥
17

Graph Construction (II)• Sampling f n at positions in  gives N
discretized feature functions
f nD  [ f n ( x1 , y1 ) f n ( x2 , y2 )  f n ( xM , yM )]T
𝑓𝑛 (𝑥, 𝑦)𝑦𝑓𝑛 (𝑥𝑖 , 𝑦𝑖 )Roadmap
Features{ f n }nN1

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

GR

22

converge

Functional
𝒔𝑖

𝑂

Ω

S  ( h)Regularizer
SG (h D )𝑥

• For each sample si , define a length N vectorv i  [f1D (i ) f 2D (i )  f ND (i)]T• Build a graph G with M vertices; each sample si  has avertex Vi
COST Training School 7/06/201518

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

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

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

M

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

converge

FunctionalS  ( h)Regularizer
SG (h D )“Distance” between two features
dij2  v i  v j2
2

where r   Cr and Cr is a constant• G is an r-neighborhood graph, i.e., no edge connecting two
vertices with distance greater than r

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

Features{ f n }nN1

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

G  R 22

converge

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

Regularizer
SG (h D )unnormalized Graph LaplacianL  DA• h( x, y ) :   R is some continuous candidate function— discrete version of h( x, y )h D  [h( x1 , y1 ) h( x2 , y2 )  h( xM , yM )]T• SG (h D )  (h D )T Lh D — graph Laplacian regularizer, functional in R M

20

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

1



 det G  dxdy

Features{ f n }nN1

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

2 1

G  R 22h  [ x h  y h]T is the gradient of hconverge

Functional• G is a 2-by-2 matrix:  N   x f n 2
n 1
G N

  n 1  x f n · y f n

S  ( h)  f · f 
  f · f    f  Regularizer
SG (h D )

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 1• G is computed from {f n }n 1 on a point-by-point basis
N[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)• Theorem : convergence of SG to S

Roadmaplim SG  h D  ~ S  h 

Features{ f n }nN1

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

M 
 0

1) number of samples M increases2) neighborhood r   Cr shrinks

G  R 22

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

FunctionalS  ( h)Regularizer
SG (h D )• 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  [f1D (i ) f 2D (i )  f ND (i)]T

COST Training School 7/06/2015[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 h    h  G h 
T



1

 det G  dx dyG   n 1 f n · f n 
N

2 1

T

SG (h D )  (h D )T Lh D

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

• The eigen-space of G reflects dominant directions of {f n }nN1• S integrates the gradient norm
COST Training School 7/06/2015

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 (h)   (h)T G 1 (h)dxdy
G   n 1 f n · f n 

COST Training School 7/06/2015

N

T

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!
COST Training School 7/06/2015

25

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

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

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

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

COST Training School 7/06/2015

26

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

G 0 (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 spaceposterior prob. of ground truth∆ is the whole gradient domainG  arg min  G  G 0 (g ) F Pr g | g k k 0 dg
2

・



G



K 1



 G   G 0 (g )·Pr g | g k k 0 dg
・

COST Training School 7/06/2015

K 1

(1)
27

Seeking for the Optimal Metric Space (II)• Assume the prior Pr (g ) is a 2D Gaussian with covariance  g2 I
we derive





2
 1
Pr g | g k k 0 
exp   2 g  g  2 
2
2
 2

where the “ensemble” mean g  and variance  2 are

g 

K 1

1

1

g
2  k 0 k
K 1

K   e2  g 
K    g2
2

2
e
2
e

noise variance of gk• Carrying out the integral in (1) gives the optimal metric space

G ・  g  g T  ( 2   )I (2)
T
・
2
g
g

• Intuition: If noise
is small,   dominates and G is discriminant;2
if  is large, ( 2   )I dominates, G ・ defaults to Euclidean space!
COST Training School 7/06/2015

28

From Metric Space to Graph Laplacian• The structure of G ・  g  g T  ( 2   )I allows us to select N  3feature functions, such that they lead to the optimal metric space:

f1D (i )   2   ·xi
f3 D 

f 2D (i )   2   ·yi

1K  
2
e

2
g

K 1
k 0

pk

— Spatial— Intensity• f1D (i ) and f 2D (i ) correspond to the term ( 2   )I in G ・D• f3 (i ) leads to the term g  g T in G ・• Our work is closely-related to joint (or cross) bilateral filtering, with theaveraging of similar patches as guidance image.
• However, we adapt to noise, resulting in robust weight estimates.
COST Training School 7/06/2015

29Formulation and Algorithm• Adopt a patch-based recovery framework, for a noisy patch p 0
1. Find K  1 patches similar to p 0 in terms of Euclidean distance.
2.

Compute the feature functions, leading to edge weights and Laplacian.

3.

Solve the unconstrained quadratic optimization:q*  arg min p0  q 2  qT Lq  q  I  L  p0
1

2

q

to obtain the denoised patch.
• Aggregate denoised patches to form an updated image.

• Denoise the image iteratively to gradually enhance its quality.
• Optimal Graph Laplacian Regularization for Denoising (OGLRD).
COST Training School 7/06/201530

Experimentation (I)• Test images: Lena, Boats, Peppers and Airplane• i.i.d. Additive White Gaussian Noise (AWGN)• Compare OGLRD to NLM and BF

1.5 dB better than NLM!

31

Experimental Results (II)• Visual comparisons (  n  25 ) of fragments

BF

NLM

OGLRD
32

Experimental Results (III)• Some visual results when  n  30

BeforeAfter

BeforeAfter

COST Training School 7/06/2015

BeforeAfter
33

Summary• Image denoising is an ill-posed problem; we use graph Laplacian
regularizer as prior for regularization.
• Graph Laplacian regularizer with Gaussian kernel weights converges to acontinuous functional.
• Analysis of the continuous functional provides theoretical justification of
why and to what extent the graph Laplacian regularizer can bediscriminant.
• We describe a methodology to derive the optimal edge weights given
nonlocal noisy gradient observations.

• Our denoising algorithm with graph Laplacian regularizer and gradientbased similarity out-performs NLM by up to 1.5 dB.

COST Training School 7/06/2015

34

Outline• Depth Image Denoising• Graph Sparsity Prior• Graph-signal Smoothness Prior• Bit-depth Enhancement

COST Training School 7/06/2015

35

Image Bit-depth Enhancement
Problem:
low bit-depth (LBD) image y—a quantizedversion of underlying HBD image

COST Training School 7/06/2015[11] Rudin et al, “Nonlinear total variation based noise removal algorithms”, Elsevier, 1992

an estimate of theoriginal HBD image

36

Image Bit-depth Enhancement
Objective: find x̂ that minimizes mean-squared-error (MSE),

x̂

MMSE

 arg min  x̂  x
x̂

22

f x | y  dx

squared err
posterior prob of HBD
signal x given LBD signal yPosterior:
Likelihood: equals to 1 iff xi quantizes
to yi

Smoothness prior: HBD signal issmoothConventional smoothness (e.g., total Variation)
are signal-independent  over-smoothing

Question: what’s a good signal smoothness prior?
COST Training School 7/06/201537

Image Bit-depth
Enhancement
Graph-signal smoothness prior

L is the graph Laplacian matrix
describing inter-pixel similarities*
Reconstruct smooth signal without blurring edges• MMSE problem is now well posed, but difficult to solve.

[12] P. Wan, G. Cheung, D. Florencio, C. Zhang, O. Au, “Image Bit-depth Enhancement via Maximum-a-Posteriori Estimation of Graph
AC Component," IEEE International Conference on Image Processing, Paris, France, October, 2014. (Top 10% accepted paper recognition)

38

Image Bit-depth Enhancement• MAP finds smoothest solution in feasible space.
• Can have arbitrarily large MSE!

Smoothest feasible signal
is basically DC.

COST Training School 7/06/2015

39

Image Bit-depth Enhancementx x A  x D

Proposed ACDC Algorithm:

Why better?

• Compute edge weights from quantized signal.

• Posterior ≈ likelihood * prior• Compute MAP solution of AC signal.

• Likelihood of AC f ( y |x A ) is less skewed• Compute MMSE solution of DC signal given AC signal. (integrating over possible DC values).

COST Training School 7/06/201540

Experiments
Numerical comparison:

COST Training School 7/06/2015

41

Experiments
Visual comparison:

42

Summary• Inverse imaging requires good signal priors.
• Depth Image Denoising• Graph Sparsity Prior (probabilistic interpretation)• Graph-signal Smoothness Prior (deterministic interpretation)• Bit-depth Enhancement• Instead of fidelity term, restricted feasible space due toquantization bin constraints (as likelihood term).

COST Training School 7/06/2015

43

Conclusion
Depth Image Coding & Processing• Coding: graph Fourier Transform (GFT), generalized graphFourier Transform (GGFT)• Denoising: graph sparsity prior, graph-signal smoothness prior

Future Work• Natural image coding using graph-based transforms.
• Depth image denoising / interpolation for non-AWGN noise.
• Apps: Given depth images, foreground / backgroundsegmentation, tracking, face modeling, etc.

COST Training School 7/06/201544