Gene Cheung
Associate Professor, York University
26th September, 2018

Graph Spectral Image Processing.
1

Acknowledgement
Collaborators:• Y. Nakatsukasa (NII, Japan)• S. Muramatsu (Niigata, Japan)• A. Ortega (USC, USA)• D. Florencio (MSR, USA)
• P. Frossard (EPFL, Switzerland)• J. Liang, I. Bajic (SFU, Canada)• X. Wu (McMaster U, Canada)• V. Stankovic (U of Strathclyde, UK)• P. Le Callet (U of Nantes, France)• X. Liu (HIT, China)• W. Hu, J. Liu, Z. Guo, W. Gao (Peking U., China)• X. Ji, L. Fang (Tsinghua, China)• Y. Zhao (BJTU, China)• C.-W. Lin (National Tsing Hua University, Taiwan)• E. Peixoto, B. Macchiavello, E. M. Hung (U. Brasilia,
Brazil)2

VISION: SCIENCE TO
APPLICATIONS (VISTA)
www.yorku.ca/vista

3

Vision: Science to Applications Program
Funded by: Canada First Research Excellence Fund

• VISTA’s Mission: to support transdisciplinary vision research,bridging fundamental and applied vision research, bothcomputational and biological.
• VISTA involves CVR and 5 other York research groups, 5
Faculties, 49 investigators (so far), and over 50 partners(hospitals, industry, international research groups, etc.)• VISTA funds new faculty positions, research chairs, staff, postdoctoral fellows, graduate students, research grants, travel
scholarships and commercialization.

• Federal funding of $33.34M; partners bring total to $120M.
• For more info, see:
http://vista.info.yorku.ca/
4

VISTA Research Concept
 Interdisciplinary Axis 5

 Science to Application Axis 

VISTA @ Innovation York

Opportunities for International Collaboration

6

Opportunity

Description

AmountResearch Grants

Partner-based research projects$50,000 per project

Travel Awards

Travel/accommodations for faculty andtrainees to undertake collaborative
research with partner institutions$7,500-$10,000 per
visit

Prototyping Funds

Development of prototypes for Yorkresearchers in the vision sciences$50,000 per projectGraduate
Scholarships

Salary stipends for Masters and PhD$10,000 per year
students above other external or internal
funds

Postdoctoral
Fellowships

Salary plus research allowance andtravel funds over a 2-year period

For more information: yorku.ca/vista$55,000-$70,000
per year + $8,500
per yearVISTA Contacts• Doug Crawford, Scientific Director− jdc@yorku.ca• Richard Wildes, Associate Director− wildes@cse.yorku.ca• Jennifer Teague, Program Manager− jteague@yorku.ca

7Outline• GSP Fundamentals• GSP for Image Compression• Optimality of GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV• Deep GLR• Ongoing & Future Work8

Outline• GSP Fundamentals• GSP for Image Compression• Optimality of GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV• Deep GLR• Ongoing & Future Work
9

Digital Signal Processing• Discrete signals on regular data kernels.
• Ex.1: audio on regularly sampled timeline.
• Ex.2: image on 2D grid.

• Harmonic analysis tools (transforms, wavelets)
for diff. tasks:• Compression.
• Restoration.
• Segmentation, classification.

10

Smoothness of Signals•

Signals are often smooth.

• Notion of frequency, band-limited.
N 1• Ex.: DCT:
 
1 
X k   xn cos  n  k 
2 
n 0
N

2D DCT basisa xdesired signaltransform

 a0 
a 
 1
a   0  Typical pixel blocks have  almost no high frequency  components.
Compact signal
representation 0transform coeff.

11Graph Signal Processing• Signals on irregular data kernels described by graphs.
• Graph: nodes and edges.
• Edges reveals node-to-node relationships.

1. Data domain is naturally a graph.
• Ex: ages of users on social networks.

example graph-signal

2. Underlying data structure unknown.
•

Ex: images: 2D grid → structured graph.

Graph Signal Processing (GSP) addresses the problem ofprocessing signals that live on graphs.
[1] D. I. Shuman et al.,”The Emerging Field of Signal Processing on Graphs: Extending High-dimensional Data Analysis to Networks
and other Irregular Domains,” IEEE Signal Processing Magazine, vol.30, no.3, pp.83-98, 2013.

12

Graph Signal Processing
Research questions*:• Sampling: how to efficiently acquire / sense
a graph-signal?

node
edge

• Graph sampling theorems.

• Representation: Given graph-signal, how tocompactly represent it?
• Transforms, wavelets, dictionaries.

• Signal restoration: Given noisy and/or
partial graph-signal, how to recover it?
• Graph-signal priors.

*Graph Signal Processing Workshop, Philadelphia, US, May, 2016. https://alliance.seas.upenn.edu/~gsp16/wiki/index.php?n=Main.Program
13*Graph Signal Processing Workshop, Pittsburgh, US, May, 2017. https://gsp17.ece.cmu.edu/*Graph Signal Processing Workshop, Lausanne, Switzerland, June, 2018. https://gsp18.epfl.ch/undirected graph

Graph Fourier Transform (GFT)• Degree Matrix D: diagonal matrix w/ entry Di,ibeing sum of column entries in row i of A.
Di ,i   Ai , jj• Combinatorial Graph Laplacian L: L = D-A• L is symmetric (graph undirected).
• L is a high-pass filter.
• L is related to 2nd derivative.

1
3

4

 0
w
A   1, 2
 0

 0

w1, 2

0 0
1 0
0 1

1 0

 w1, 2
 0
D
 0

 0

0 0
w1, 2  1 0 0
0
2 0

0
0 1

 w1, 2 w
L   1, 2
 0

 0

0
1
0
0

 w1, 20
w1, 2  1  1 0 
1
2  1
0
1 1 
0

L3,: x   x2  2 x3  x4
f x  h   2 f x   f x  h f  x   lim
h 0*https://en.wikipedia.org/wiki/Second_derivative

1
2

1

Graph Laplacian:• Adjacency Matrix A: entry Ai,j has non-negativeedge weight wi,j connecting nodes i and j.

w1, 2

h2

14

Graph Spectrum from GFT• Graph Fourier Transform (GFT) is eigen-matrix of graph Laplacian L.
eigenvalues along diagonal

L  V  VT

GFT~
x  VT x

1

w1, 2

1

2

1
3

4

…

8

1st AC eigenvector

eigenvectors in columns
GFT coefficients

1. Edge weights affect shapes of eigenvectors.
2. Eigenvalues (≥ 0) as graph frequencies.
• Constant eigenvector is DC.
• # zero-crossings increases as λ increases.

• GFT defaults to DCT for un-weighted connected line.
• GFT defaults to DFT for un-weighted connected circle.
15

Variants of Graph Laplacians• Graph Fourier Transform (GFT) is eigen-matrix of graph Laplacian L.
eigenvalues along diagonal

L  V  VT

GFT

eigenvectors in columns• Other definitions of graph Laplacians:• Normalized graph Laplacian:Ln  D 1 / 2 LD 1 / 2  I  D 1 / 2 AD 1 / 2• Random walk graph Laplacian:1

1

Lrw  D L  I  D A• Generalized graph Laplacian [1]:Lg  L  D*

Characteristics:• Normalized.
• Symmetric.
• No DC component.
• Normalized.
• Asymmetric.
• Eigenvectors not orthog.

• Symmetric.
• L plus self loops.
• Defaults to DST, ADST.

[1] Wei Hu, Gene Cheung, Antonio Ortega, "Intra-Prediction and Generalized Graph Fourier Transform for Image Coding," IEEE
Signal Processing Letters, vol.22, no.11, pp. 1913-1917, November 2015.

16

GSP and Graph-related Research
GSP: SP framework that unifies concepts from multiple fields.

Partial Differential
Eq’ns

MachineLearning

Laplaceequationgraphical model,
manifold learning,classifier learning

LaplaceBeltrami
operator

Computer Graphics

Graph Signal
Processing* (GSP)
DSP

spectralclustering

Computer
VisionMax cut, graphtransformation

eigen-analysis ofgraph Laplacian,adjacency matrices

Combinatorial
Graph Theory

SpectralGraph Theory17

Outline• GSP Fundamentals• GSP for Image Compression• Optimality of GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV• Deep GLR• Ongoing & Future Work18

GFT for Image Compression• DCT are fixed basis. Can we do better?
• Idea: use adaptive GFT to improve sparsity [1].
1. Assign edge weight 1 to adjacent pixel pairs.
2. Assign edge weight 0 to sharp signal discontinuity.
3. Compute GFT for transform coding, transmit coeff.

~x  VT x

GFT

4. Transmit bits (contour) to identify chosen GFT to

decoder (overhead of GFT).
[1] G. Shen et al., “Edge-adaptive Transforms for Efficient Depth Map Coding,”
IEEE Picture Coding Symposium, Nagoya, Japan, December 2010.
[2] W. Hu, G. Cheung, X. Li, O. Au, “Depth Map Compression using Multi-resolution Graph-based Transformfor Depth-image-based Rendering,” IEEE International Conference on Image Processing, Orlando, FL, September 2012.

19

GFT: Derivation of Optimal Edge Weights• Assume a 1D 1st-order autoregressive (AR) process

where,

0-mean r.v. with large var. σ2

k 1
xk  
 xk 1  ek 1  k  N

0-mean r.v. with var. σk2

x1  
x2  x1  e2
x N  x N 1  eN

1 0 0 0 0
 1 1 0 0 0
0   0 0
F
 0 0 1 1 0
0 0 0   0 0 0 0  1

0
 
e 
0
 2

0
, b   
0
 
 
0
 

1
eN 
20

GFT: Derivation of Optimal Edge Weights• Covariance matrix

C  E x   x   



T





1



  

 F 1 E b bT F 1

T

 

 

 E x x  E F bb F
T 2 0  0 


2
0


2

E bb T  
 
 0 

2
 0  0  N 

1 T

T

• Precision matrix (tri-diagonal)1
 1

 2  22
  1
  22Q  C 1  
0





0


large



1

 22

1

1

2

 32


2









0
1

 22


1


1



1

 N2 1

 N2 1  N2

0



1

 N2


0 
Graph Laplacian matrix!
 
0  L
1 
1
1
1
 2



N 
N
3 …
2
1
1 
21
 N2 
2
2

2
3

2
N

GFT for PWS Image Coding• Graph Laplacian ≈ Precision Matrix → GFT approx.
Karhunen-Loeve Transform (KLT).
• Encode blocks with signal-decorrelation GFT.
Rate of transform description T

Rate of transform coefficient vector• To limit the description costa small discrete set, restrict weights to
3500

Histogram of inter-pixel difference

strong

3000

2500- "1": strong correlation in smooth regions

2000- "0": zero correlation in sharp boundaries

1500- "c": weak correlation in slowly-varying parts

1000

zero

weak

500

22
0
0

5

10

15

20

25

30

35

40

22
45

50

Transform Comparison

Transform Representation

Transform Description

Karhunen-Loeve
Transform (KLT)“Sparsest” signal representation given
available statistical modelCan be expensive (if poorly
structured)

Discrete Cosine
Transform (DCT)

non-sparse signal representation
across sharp boundaries

little (fixed transform)

Graph Fourier
Transform (GFT)

minimizes the total rate of signal’s transform representation &transform description[1] Wei Hu, Gene Cheung, Antonio Ortega, Oscar Au, "Multiresolution Graph Fourier Transform for Compression of
Piecewise Smooth Images," IEEE Transactions on Image Processing, vol.24, no.1, pp.419-433, January 2015.

23

Experimentation• Setup- Test images: depth maps of Teddy and Cones, and graphics images of Dude and Tsukuba.

- Compare against: HR-DCT, HR-SGFT, SAW, MR-SGFT in H.264.

• Results

HR-DCT:
HR-SGFT:SAW:
MR-SGFT:

24

6.8dB
5.9dB
2.5dB
1.2dB

Subjective Results

HR-DCT

HR-SGFTMR-GFT

25

Summary of GFT for Image Coding• Optimality of GFT for AR model.
• Variants of GFT for prediction residuals, anti-correlated pixels.
• Fast implementation (w/o eigen-decomposition) via Graph
Lifting Transform (GLT) [1] or Fast Graph Fourier
Transform (FGFT) [2].

[1] Y.-H. Chao et al., "Edge-Adaptive Depth Map Coding with Lifting Transform on Graphs," 31st Picture Coding
Symposium, Cairns, Australia, May, 2015.
[2] L. Le Magoarou et al., "Approximate Fast Graph Fourier Transforms via Multilayer Sparse Approximations,"
IEEE TSIPN, May, 2018.

26

Graph-Signal Sampling / Encoding for 3D
Point Cloud• Problem: Point clouds require encoding specific 3D coordinates.
• Assumption: smooth 2D manifold in 3D space.
• Proposal: progressive 3D geometry rep. as series of graph-signals.

MIT 𝑑𝑎𝑡𝑎𝑠𝑒𝑡 ∗1. adaptively identifies new samples on the manifold surface, and
2. encodes them efficiently as graph-signals.

• Example:
1.

Interpolate 𝑖𝑡ℎ iteration samples (black circles) to a continuous kernel (mesh),an approximation of the target surface S.

2.

New sample locations, knots (squares), are located on the kernel surface.

3.

Signed distances between knots and S are recorded as sample values.

4.

Sample values (green circles) are encoded as a graph-signal via GFT.
27

Graph-Signal Sampling / Encoding for 3D
Point Cloud• Experimental Results:(a) Dataset1(b) Dataset2[1] M. Zhao, G. Cheung, D. Florencio, X. Ji, "Progressive Graph-Signal Sampling and Encoding for Static 3D Geometry
Representation," IEEE International Conference on Image Processing, Beijing, China, September, 2017.

28

Pre-Demosiac Light Field Image Compression
Using Graph Lifting Transform• Problem: Sub-aperture images in Light field data are huge.
Raw Lenselet
Image

Demosaicked
Image

Calibrated Color
Image

Sub-apertureB B B B Images

GBGBGBGB
RGR
RGR
BGRBGBRGRBGRGG
RRRR

BBBBBBBBBBBBBB
RB
RB
RB
RB
RB
RB
RB
RB
RG
RBG
RBG
RBG
RRG
RB
B
B
B
B
B
G
G
G
G
G
G
G
G
G
R
GB
GB
GB
GB
GB
GB B BGRGRGR
GRGG
GBG
GBG
RGBG
R
R
RG
RB
G
B
B
B
R
B
B
GB
GB
G
G
G
G
G
G
RB
G
G
G
G
G
G
G
BGBGBGBBBBB
GGRGGRGGRGGR
GRGG
R
B
B
B
B
B
B
G
G
G
G
G
G
G
G
GB
RB
GB
GB
GB
GB
GB
GB
GRGRGRGRGRGRGRGRGR
G
G
G
G
G
G
G
G
R
BGBGBGG G G G G G Calibration
GGGGGGGG(Scaling,
Demosaicing
Transition,
Rotation)BBBB
GBGBGBGBRGRGRGRG
BR
BBB
RGR
G GRG
RRRR

BBBB
GBGBGBGBRGRGRGRG
BR
BBB
RGR
G GRG
RRRR

×3
× ≈ 1.5• Proposal: postpone demosiacking to decoder.
Raw Lenselet
Image

Calibrated LenseletImageG

G R G R G R
B G B G B G
G R G R G R
B G B G B G

GR
BG

GRG
BGB

R
G

GR
Re-arranged on B G

GRG
BGB

R Re-arrangedG
G

Calibrated Image

B

in 4D space

Sub-aperture
space

Image
Codinggraph-based
lifting transform

G

R
G

RG
B R
RGB

G

RR
GB
R B GG
29Pre-Demosiac Light Field Image Compression Using
Graph Lifting Transform• Experimental Results:
Dataset: EPFL light field image dataset
Baseline: All-intra HEVC coding in YUV4:2:0 and RGB 4:4:4[1] Y.-H. Chao, G. Cheung, A. Ortega, "Pre-Demosiac Light Field Image Compression Using Graph Lifting

Transform," IEEE Int’l Conf. on Image Processing, Beijing, China, September, 2017. (Best student paper award)

30

Outline• GSP Fundamentals• GSP for Image Compression• Optimality of GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV• Deep GLR• Ongoing & Future Work31

Graph Laplacian Regularizer• x T L x (graph Laplacian regularizer) [1]) is one smoothness measure.
1
2
x T Lx   wi , j xi  x j    k ~xk2
signal contains
2 i, jk• Signal Denoising:observation• MAP Formulation:
fidelity term

signal smooth innodal domainy  x v

desired signal
noise
pixel intensity diff.

min y  x 2   x T L x
2x

update edge
weights

mostly low graph freq.

I L  x *  ypixel location diff.

 x x 2   l l 2 ij i j 22 wi , j  exp
exp

 2

 12




2





smoothness prior

Bilateral filter weights

linear system of eqn’s w/ sparse, symmetric PD matrix[1] P. Milanfar, “A Tour of Modern Image Filtering: New Insights and Methods, Both Practical and Theoretical,” IEEE Signal
Processing Magazine, vol.30, no.1, pp.106-128, January 2013.

32

Graph Laplacian Regularizer

 x x 2   l l 2 
i
j 2 

 i j 2wi , j  exp
exp
2

 2

1




2





• x T Lx x promotes piecewise smooth (PWS) signal behavior [1].
2




x

x
1
1

2
2
i
j
T


x Lx x   wi , j xi  x j    ui , j exp
x

x
 i
j
2
2 i, j
2 i, j• Spectral Clustering [2]:
T
v
Ln v
*
v  arg min Tv
v v

Rayleigh quotient

s.t. v T v1  0• v1 minimizes obj → Sol’n is 2nd eigenvector of Ln.
• 2nd eigenvalue—Fiedler number—measures “connectedness”.
• PWS signal = 2 clusters of similar nodes → xTLx  0[1] X. Liu, G. Cheung, X. Wu, D. Zhao, "Random Walk Graph Laplacian based Smoothness Prior for Soft Decoding of JPEG Images,"
IEEE Transactions on Image Processing, vol.26, no.2, pp.509-524, February 2017.
[2] J. Shi and J. Malik, “Normalized cuts and image segmentation,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 22, no. 8, pp. 888–905,
Aug. 2000.

33

Analysis of Graph Laplacian Regularizer• [1] shows SG (u)  u T Lu converges to continuous functional S
and objective becomes:• Solution can be implemented as anisotropic diffusion:feature function vectorv i  [f1 (i ) f 2 (i )  f N (i)]T
distance

diffusivity

dij2  v i  v j2
2edge weightwij  ( i  j )   (dij )• it not only smooths but may also sharpens the image,• promote piecewise smooth images, like total variation (TV).

metricN
T
space G   f n ·f n[1] J. Pang, G. Cheung, "Graph Laplacian Regularization for Inverse Imaging: Analysis in the Continuous Domain," IEEE
Transactions on Image Processing, vol. 26, no.4, pp.1770-1785, April 2017.

n 1

34

Outline• GSP for Image Compression• Optimality of GFT• Generalized GFT• Signed GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer: Image Denoising• Reweighted Graph TV:• Deep GLR• Ongoing & Future Work35

Optimal Graph Laplacian Regularization for Denoising• 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.

Spatialf1D (i )   2   ·xi

2.

Compute feature functions, leading to edge weights and Laplacian.

3.

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

T

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

Intensity

1q

f3 D 

1
K  
2e2g

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).
[1] J. Pang, G. Cheung, "Graph Laplacian Regularization for Inverse Imaging: Analysis in the Continuous Domain," IEEE
Transactions on Image Processing, vol. 26, no.4, pp.1770-1785, April 2017.

36

K 1
k 0

pk

Denoising Experiments (natural images)• 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
37Denoising Experiments (depth images)• Subjective comparisons (  I  30 )Original

Noisy, 18.66 dB

BM3D, 33.26 dB NLGBT, 33.41dB OGLR, 34.32 dB[1] W. Hu et al., "Depth Map Denoising using Graph-based Transform and Group Sparsity," IEEE International Workshop on
Multimedia Signal Processing, Pula (Sardinia), Italy, October, 2013.

38

Outline• GSP Fundamentals• GSP for Image Compression• Optimality of GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV• Deep GLR• Ongoing & Future Work
39

Reweighted Graph Total Variation• TV on graphs.
Gradient of nodes on the graph:

Conventional Graph TV:

pixel intensity difference x x 2 
i
j 2 

wi , j  exp

 12





Reweighted Graph TV:[1] M. Hidane, O. Lezoray, and A. Elmoataz, “Nonlinear multilayered representation of graph-signals,” in Journal of Mathematical Imaging
and Vision, February 2013, vol. 45, no.2, pp. 114–137.
[2] P. Berger, G. Hannak, and G. Matz, “Graph signal recovery via primal-dual algorithms for total variation minimization,” in IEEE Journal
on Selected Topics in Signal Processing, September 2017, vol. 11, no.6, pp. 842–855.

40

Reweighted Graph Total Variation• RGTV is separable. analyze each node pair.
• Promotes bi-modal inter-pixel differences.

[1] Y. Bai, G. Cheung, X. Liu, W. Gao, "Blind Image Deblurring via Reweighted Graph Total Variation," IEEE International Conference on
Acoustics, Speech and Signal Processing, Calgary, Alberta, April, 2018.

41

Outline• GSP Fundamentals• GSP for Image Compression• Optimality of GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV: Image Deblurring• Deep GLR• Ongoing & Future Work
42

Background for Image Deblurring• Image blur is a common image degradation.
• Typically, blur process is modeled:𝑦 = 𝑘𝑥
where y is the blurry image, k is the blur kernel, x is theoriginal sharp image.

• Blind-image deblurring focuses on estimating blur kernel k.
• Given k, problem becomes de-convolution.
43

Observation• Skeleton image:• PWS image keeping only structural edges.
• Proxy to estimate blur kernel k.

44

Observation• Examine statistical properties of local patches:• Edge weight distribution of a fully connected graph.

• Skeleton Image enjoys both Sharpness and bi-modal Weight distribution,
thus useful to estimate blur kernel.
45

Key Idea• Propose a Reweighted Graph Total Variation (RGTV) to
promote a skeleton image patch.

Conventional Graph TV:

Reweighted Graph TV:

46

Our algorithm• The optimization function can be written as follows,መ = argmin 𝜑 𝐱𝐤 − 𝐛 + 𝜇1 ∙ 𝜃𝑥 𝐱 + 𝜇2 ∙ 𝜃𝑘 (𝐤)
𝐱ෝ, 𝐤
𝐱,𝐤

• Assume 𝐿2 norm for fidelity term 𝜑(∙).
• 𝜃𝑥 ∙ = 𝑅𝐺𝑇𝑉(∙).
• 𝜃𝑘 ∙ = | ∙ |2 , assuming zero mean Gaussian distribution of k.
• RGTV is non-differentiable and non-convex.
Solution:• Solve x and k alternatingly.
• For x, spectral interpretation of GTV, fast spectral filter.
47

Spectral domain• Deduction for spectrum of GTV

New weightfunction48

Spectral domain
New
Adjacency
matrix 𝚪• Explanation:x GTV    i , j xi  x j 

2( i , j )E xT L x

Graph L1 spectrum

49

Our algorithm
Alternating Iterative algorithm:1
෡ − 𝐛 2 + 𝛽 𝐱 𝑅𝐺𝑇𝑉
𝐱ො = argminx
𝐱⊗𝒌
2
൞
2
መ = argmink 𝛻𝐱ො ⊗ 𝐤 − 𝛻𝐛 22 + 𝜇 𝐤 22
𝐤2෡𝐱ො = argminx 𝐱 ⊗ 𝒌 − 𝛻𝐛 2 + 𝛽 ⋅ 𝐱 𝐓 𝐋Γ 𝐱

෡ 𝑇𝐊
෡ + 2β ⋅ 𝐋Γ 𝐱ො = 𝐊
෡ 𝐓𝐛𝐊System of linear equations.
Efficiently solved via conjugate gradient.

50

Workflow

Blurry Image

Skeleton
Image
Reconstruction

Reconstruction

Kernel Estimation

51

Experimental Results

52Experimental Results

53Outline• GSP Fundamentals• GSP for Image Compression• Optimality of GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV: 3D Point Cloud Denoising• Deep GLR• Ongoing & Future Work54

GTV for Point Cloud Denoising• Acquisition of point cloud introduces noise.
• Point cloud is irregularly sampled 2D manifold in 3D space.
• Not appropriate to apply GTV directly on 3D coordinates [1].
• only a singular 3D point has zero GTV value.

 y y f (x )  f (x )
i

i 1

i

i

i

i 1

 n nii 1 1iyf(x)geometrysmoothness

functionalsmoothnessxx• Proposal: Apply GTV is to the surface normals of 3D point
cloud—a generalization of TV to 3D geometry.

[1] Y. Schoenenberger, J. Paratte, and P. Vandergheynst, “Graph-based denoising for time-varying point clouds,” inIEEE 3DTV-Conference, 2015, pp. 1–4

55

Algorithm Overview• Use graph total variation (GTV) of surface normals over the K-NN graph:
n𝒊

n𝒋

𝑗
𝑖
• Denoising problem as l2-norm fidelity plus GTV of surface normals:

min q  p 2    wi , j n i  n j
2p ,n

i , jE

1• Surface normal estimation of n𝒊 is a nonlinear function of p𝒊 and neighbors.

Proposal:
1. Partition point cloud into two independent classes (say red and blue).
2. When computing surface normal for a red node, use only neighboring blue points.

3. Solve convex optimization for red (blue) nodes alternately.
[1] C. Dinesh, G. Cheung, I. V. Bajic, C. Yang, “Fast 3D Point Cloud Denoising via Bipartite Graph Approximation& Total Variation,” IEEE 20th International Workshop on Multimedia Signal Processing, Vancouver, Canada, August 2018.

56

Bipartite Graph Approx. & Normal Def’n
Step 1: bipartite graph approx. of k-NN graph.
Normal vector estimation
at a red node

Step 2: define red nodes’ normals using blue nodes.

Ai is a constant matrix and bi is a constant vector with respect to pi[1] J. Zeng, G. Cheung, A. Ortega, "Bipartite Approximation for Graph Wavelet Signal Decomposition," IEEE
Transactions on Signal Processing, vol.65, no.20, pp.5466-5480, October 2017.

57

Convex Optimization Formulation• After computing normals for each red node, construct a new k-NN
graph for red nodes only.
• For a red node graph, objective is a l2 -l1 -norm minimization w/
linear constraints:

min q  p 2    m i , j

mi, j  n i - n j
n i  A i pi  bi

2

p ,m

Solution:2q p 2 • ADMM: minp ,mi , jErmi , jEri, j 11
2m  Bp  v

Bp  v  m  u 2  const
2• p-minimization:• m-minimization:
Proximal gradient descent• Alternately update red and blue graphs until convergence.
58

Experimental Setup• 4 competing local methods: APSS [1], RIMLS [2], AWLOP [3], MRPCA[4]• 7 pint cloud datasets used: Bunny, Gargoyle, DC, Daratrch, Anchor,
Lordquas, Fandisk, Laurana• Metrics: point to point error (C2C) and point to plane error (C2P)• Gaussian noise with zero mean, standard deviation 𝜎 of 0.1 and 0.3.
[1] G. Guennebaud and M. Gross, “Algebraic point set surfaces,” ACM Transactions on Graphics (TOG),
vol. 26, no. 3, p. 23, 2007.
[2] A. C. Oztireli, G. Guennebaud, and M. Gross, “Feature preserving point set surfaces based on nonlinear kernel regression,” in Computer Graphics Forum, vol. 28, no. 2, 2009, pp. 493–501.
[3]H. Huang, S. Wu, M. Gong, D. Cohen-Or, U. Ascher, and H. R. Zhang, “Edge-aware point set
resampling,” ACM Transactions on Graphics, vol. 32, no. 1, p. 9, 2013.
[4] E. Mattei and A. Castrodad, “Point cloud denoising via moving RPCA,” in Computer Graphics Forum,
vol. 36, no. 8, 2017, pp. 123–137.
59

Experimental Results – Visual Comparison
Anchor model (𝜎=0.3)

60

Experimental Results – Visual Comparison
Daratech model (𝜎=0.3)

61

Experimental Results – Numerical Comparison

62

Outline• GSP Fundamentals• GSP for Image Compression• Optimality of GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV: 3D Point Cloud Denoising• Deep GLR• Ongoing & Future Work63

Unrolling Graph Laplacian Regularizer• Recall MAP formulation of denoising problem with quadratic
graph Laplacian regularizer:2
min y  x 2   x T L xxfidelity termsmoothness prior• Solution is system of linear equations:I L  x *  ylinear system of eqn’s w/ sparse, symmetric PD matrix

Q: what is the “most appropriate” graph?

Bilateral weights:
 x x 2   l l 2 
i
j 2 

 i j 2wi , j  exp
exp

 2

 12




2





[1] J. Pang, G. Cheung, "Graph Laplacian Regularization for Inverse Imaging: Analysis in the Continuous Domain," IEEE
Transactions on Image Processing, vol. 26, no.4, pp.1770-1785, April 2017.

64

Unrolling Graph Laplacian Regularizer• Deep Graph Laplacian Regularization:
1. Learn features f’s using CNN.
2. Compute distance from features.
3. Compute edge weights using Gaussian kernel.
4. Construct graph, solve QP.
featurespre-filter

weight para

[1] M. McCann et al., “Convolutional Neural Networks for Inverse Problems in Imaging,“ IEEE SPM, Nov. 2017.
[2] K. Gregor and Y. LeCun, “Learning fast approximations of sparse coding,” in Proc. 27th Int. Conf. Machine Learning, 2010..

65

Unrolling Graph Laplacian Regularizer

66[1] J. Zeng et al., “Deep Graph Laplacian Regularization,“ submitted to arXiv, July 2018. (https://arxiv.org/abs/1807.11637 )Unrolling Graph Laplacian Regularizer• Graph Model guarantees numerical stability of solution:

I L  x  y*• Thm 1: condition number κ of matrix satisfies [1]:maximum node degree• Observation: By restricting search space of CNN to degree-bounded graphs,
we achieve robust learning.
67[1] J. Zeng et al., “Deep Graph Laplacian Regularization,“ submitted to arXiv, July 2018. (https://arxiv.org/abs/1807.11637 )Experimental Results – Numerical Comparison• Trained on AWGN on 5 images, patches of size 26-by-26.
• Batch size is 4, model is trained for 200 epochs.
• Trained for both known and blind noise variance.

[1] Kai Zhang et al, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” TIP 2017.
[2] Marc Lebrun et al, “The noise clinic: a blind image denoising algorithm,” IPOL 2015.

68Experimental Results – Numerical Comparison• DeepGLR has average PSNR of 0.34 dB higher than CDnCNN [1].
• Model-based provides robustness against overfitting.

[1] Kai Zhang et al, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” TIP 2017.
[2] Marc Lebrun et al, “The noise clinic: a blind image denoising algorithm,” IPOL 2015.

69

Experimental Results – Visual Comparison• trained on Gaussian noise, tested on low-light images in (RENOIR).
• Competing methods: DnCNN [1], noise clinic [2].
• outperformed DnCNN by 5.52 dB, and noise clinic by 1.87 dB.

DnCNN

clinic[1] Kai Zhang et al, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” TIP 2017.
[2] Marc Lebrun et al, “The noise clinic: a blind image denoising algorithm,” IPOL 2015.

DeepGLR70

Experimental Results – Visual Comparison• trained on Gaussian noise, tested on low-light images in (RENOIR).
• Competing methods: DnCNN [1], noise clinic [2].
• outperformed DnCNN by 5.52 dB, and noise clinic by 1.87 dB.

DnCNN

clinic[1] Kai Zhang et al, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” TIP 2017.
[2] Marc Lebrun et al, “The noise clinic: a blind image denoising algorithm,” IPOL 2015.

DeepGLR71

Experimental Results – Visual Comparison• trained on Gaussian noise, tested on low-light images in (RENOIR).
• Competing methods: DnCNN [1], noise clinic [2].
• outperformed DnCNN by 5.52 dB, and noise clinic by 1.87 dB.

DnCNN

clinic[1] Kai Zhang et al, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” TIP 2017.
[2] Marc Lebrun et al, “The noise clinic: a blind image denoising algorithm,” IPOL 2015.

DeepGLR72

Outline• GSP Fundamentals• GSP for Image Compression• Optimality of GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV: 3D Point Cloud Denoising• Deep GLR• Ongoing & Future Work73

Summary• Variants of GFTs for optimal decorrelation• GFT, GGFT, SGFT
• Selection of statistical model vs. encoding cost of side information• GSP for Inverse Imaging• PWS-promoting Graph Laplacian Regularizer, RGTV
• Spectral interpretation of GTV, RGTV

• Graph-based model restricts search space of DNN.
• Robustness against overfitting.

74

Ongoing & Future Work• Unrolling of graph-based convex optimization.
• Unrolling of ADMM, proximal gradient with GTV prior, convex set
constraints.
• Learn (sparse) connectivity, edge weights.
• Learn features from RGBD images for depth inpainting / denoising.

• Model-guided learning safeguard against worst-case /adversary noise?

75

Q&A
• Email: genec@yorku.ca• Homepage: https://www.eecs.yorku.ca/~genec/index.html

76