Gene Cheung
National Institute of Informatics
8th June, 2018

Recent Advances in
Graph Spectral Image Processing.
GSP Workshop 6/8/2018

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)
GSP Workshop 6/8/2018

2

Why GSP for Image Processing?
• GSP: signals on irregular data kernels described bygraphs.
• Graph: nodes and edges.
• Edges reveals node-to-node relationships.

1. Graph captures underlying image statistics.

example graph-signal• Ex: GFT decorrelates signal for compression.

2. Graph captures (dis)similarities of pixels.
•

Ex: Bilateral filter weights based on Euclidean / photometricdistance.

GSP provides design space of spectrum forimage filtering.

3[1] G. Cheung et al. "Graph Spectral Image Processing," Proceedings of IEEE, vol.106, no.4, pp. 907-930, May 2018.

Outline• GSP for Image Compression• Optimality of GFT• Generalized GFT• Signed GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV• Deep GLR• Ongoing & Future Work
GSP Workshop 6/8/2018

4

Outline• GSP for Image Compression• Optimality of GFT• Generalized GFT• Signed GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV• Deep GLR• Ongoing & Future Work
GSP Workshop 6/8/2018

5

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.

6

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
GSP Workshop 6/8/2018

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

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 
8
 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

9
0
0

5

10

15

20

25

30

35

40

9
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.

10

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:

11

6.8dB
5.9dB
2.5dB
1.2dB

Subjective Results

HR-DCT

HR-SGFTGSP Workshop 6/8/2018

MR-GFT

12

Generalized GFT

Intra-prediction in H.264• Intra-prediction
𝑥0 𝑥1 𝑥2 𝑥3 𝑥4
Boundary pixel(predictor)

Predicted pixels 𝑥𝑖: prediction residuals• Statistical model for prediction residuals?

[1] W. Hu et al., “Intra-Prediction and Generalized Graph Fourier Transform for Image Coding,” IEEE Signal
Processing Letters, vol.22, no.11, pp. 1913-1917, November, 2015.

13

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

where,
predictor

predictor from decoded block x0  e1 k  1
xk  
 xk 1  ek 1  k  N
x1  e1   x0
x2  x1  e2
x N  x N 1  eN

x0 x1 x2
Prediction
residual

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

Fx  bc

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

x F bF c-10 e1 
 x0 
e 
 0 
0
 2



  
0
, b   , c  

0
 


 


0
 



e
1
0


 N 

GSP Workshop 6/8/2018

Optimalpredictor-1

 x0 F 1c     x0 

14

GGFT: Derivation of Optimal Edge Weights• Covariance matrix

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

T



0

  
 
 F E b b F 
 E r r T  E F 1bbT F
1

T

1

1

1



 22



 N2

2
1

1

2

2
3

3

N

Boundary
condition

T 1

T 1

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

1



1

1

1

2

 32









0

 22


2

diagonal matrix D

1

 22


1


1



1

 N2 1

 N2 1  N2

0



1

 N2


0 
Generalized graph  1Laplacian
0

0


 2

  1

0
0


0



1

 
 2  N   0 
0


1 
Defaults to ADST if σi2=1
 N2 

[1] J. Han et al., “Jointly Optimized Spatial Prediction and Block Transform for video and Image Coding,” IEEE TIP, April 2012.

[2] G. Strang, “The discrete cosine transform,” SIAM Review, vol. 41, pp.135–147, 1999.

15

Experimental Results• Test images: PWS images and natural images• Compare proposed intra-prediction (pIntra) + GGFT against:- edge-aware intra-prediction (eIntra) + DCT- eIntra + ADST
- eIntra + GFT

York University 2/26/2018

16

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

where,

negative correlation

k 1


xk   xl 1  el k  lx  e k 1 k 1  k  N
x1  
x2  x1  e2
xl  xl 1  elx 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  
l-th row 0 0 0 0  1
GSP Workshop 6/8/2018

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

SGFT: Derivation of Optimal Edge Weights
2

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


1



2
2



1



0

 22


1

1



2
3



1

 22


1

1

1

1

 l21 l21  l2

 l2





1

1

1



l



2
l


2

0

2
l 1





1



2
l


1

 N2

 22


0 

 









1 

 N2 



l-1

2

2
l



 l2

1

 l2

l

… N

Graph Laplacian withnegative edges &self-loops
l-th row(l+1)-th row

18[1] Weng-Tai Su, Gene Cheung, Chia-Wen Lin, "Graph Fourier Transform with Negative Edges for Depth
Image Coding," IEEE International Conference on Image Processing, Beijing, China, September, 2017.

SGFT: Derivation of Optimal Edge Weights
2

1• SGFT approx. KLT for anti-correlation.
• Self-loops ensure PSD (Gershgorin).
• 1st eigenvector is PWC.
• 2nd eigenvector is not constant.

1

 22



l-1

2

2
l



 l2

1

 l2

l

… N

• Negative Edges:• Improve classifier robustness [1].
• Enhance image contrast.

[1] Gene Cheung, Weng-Tai Su, Yu Mao, Chia-Wen Lin,
"Robust Semi-Supervised Graph Classifier Learning with
Negative Edge Weights," accepted to IEEE Transactions on
Signal and Information Processing over Networks, March 2018.

19

Summary of GFT for Image Coding• Optimality of GFT, Generalized GFT, Signed GFT for variants of
AR models.
• Selection of models trades off encoding cost of SI.
• 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.

20

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.
GSP Workshop 6/8/2018

21

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.

22

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
RGRG
BGR
BBB
RGR
G GRG
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
23

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)

24

Outline• GSP for Image Compression• Optimality of GFT• Generalized GFT• Signed GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV• Deep GLR• Ongoing & Future Work
GSP Workshop 6/8/2018

25

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.

26

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.

27

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

28

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 Work
GSP Workshop 6/8/2018

29

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.

30

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

GSP Workshop 6/8/2018

31

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.

32

Outline• GSP for Image Compression• Optimality of GFT• Generalized GFT• Signed GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer: Soft Decoding of JPEG Images• Reweighted Graph TV• Deep GLR• Ongoing & Future Work
GSP Workshop 6/8/2018

33

Soft Decoding of JPEG Images• Setting: JPEG compresses natural images:
1. Divide image into 8x8 blocks, DCT.
2. Perform DCT transform per block and quantize:
DCTqi  round (Yi / Qi ), Y =Ty
quantization parameter

8x8 pixel block

DCT Coefficients

3. Quantized DCT coeff entropy coded.

• Decoder: uncertainty in signal reconstruction:qi Qi  Yi  (qi  1) Qi , i  1, 2,, 64.

[1] A. Zakhor, “Iterative procedures for reduction of blocking effects in transform image coding,” IEEE Transactions on Circuits andSystems for Video Technology,, vol. 2, no. 1, pp. 91–95, Mar 1992.
[2] K. Bredies and M. Holler, “A total variation-based JPEG decompression model,” SIAM J. Img. Sci., vol. 5, no. 1, pp. 366–393, Mar. 2012.
[3] H. Chang, M. Ng, and T. Zeng, “Reducing artifacts in jpeg decompression via a learned dictionary,” IEEE Transactions on Signal
Processing,,vol. 62, no. 3, pp. 718–728, Feb 2014.

34

LERaG for Soft Decoding of JPEG Images• Problem: reconstruct image given indexed quant. bin in 8x8 DCT.

• Procedure:
1. Initialize per-block
MMSE sol’n via
Laplacian prior.

1 T
d min
   ~k2  k2  0k

2. Solve per-patch signal
restoration problem w/ 2priors:1. Sparsity prior
2. Graph-signal
smoothness prior[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.

35

Soft Decoding Algorithm w/ Prior Mixture• Objective: fidelity termsparsity priorgraph-signal smoothness
prior

graph-signal,
code vector

quantization bin constraint• Optimization:
1. Laplacian prior provides an initial estimation;
2. Fix x and solve for α;
3. Fix α and solve for x.

GSP Workshop 6/8/201836

Evolution of 2nd Eigenvector• 2nd eigenvector becomes more PWS:

• PWS means:
1. better pixel clusters,
2. smaller Fiedler number (2nd eigenvalue),
3. Smaller smoothness penalty term.

GSP Workshop 6/8/2018

37

Experimental Setup• Compared methods• BM3D: well-known denoising algorithm• KSVD: with a large enough over-complete dictionary (100x4000); ourmethod uses a much smaller one (100x400).
• ANCE: non-local self similarity [Zhang et al. TIP14]• DicTV: Sparsity + TV [Chang et al, TSP15]• SSRQC: Low rank + Quantization constraint [Zhao et al. TCSVT16]York University 2/26/2018

38

PSNR / SSIM Comparison

GSP Workshop 6/8/2018

39

Subjective Quality Evaluation40

Subjective Quality Evaluation

41

Outline• GSP for Image Compression• Optimality of GFT• Generalized GFT• Signed GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV• Deep GLR• Ongoing & Future Work
GSP Workshop 6/8/2018

42

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.

43

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.

44

Outline• GSP for Image Compression• Optimality of GFT• Generalized GFT• Signed GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV: Image Deblurring• Deep GLR• Ongoing & Future Work
GSP Workshop 6/8/2018

45

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.
GSP Workshop 6/8/2018

46

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

47

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.
GSP Workshop 6/8/201848

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

Conventional Graph TV:

Reweighted Graph TV:

GSP Workshop 6/8/2018

49

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.
GSP Workshop 6/8/201850

Spectral domain• Deduction for spectrum of GTV

New weightfunction
GSP Workshop 6/8/2018

51

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

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

GSP Workshop 6/8/2018

52

Spectral domain
Secondeigenvector

Relative
λeigenvalues 𝑘
𝜆2• We do 1D illustrative experiments to compute Graph
Spectrum with conventional Graph Laplacian and GTV
Laplacian.
• The parameter are the same for both. 4 neighbors and¥sigma=0.3
GSP Workshop 6/8/2018

53

Spectral domain
Secondeigenvector

Relative
λeigenvalues 𝑘
𝜆2After addingnoise

GSP Workshop 6/8/2018

54

Spectral domain
Secondeigenvector

Relative
λeigenvalues 𝑘
𝜆2

After
Blurring

GSP Workshop 6/8/2018

55

Spectral Domain• Weights in RGTV is functions of graph signal. Therefore,
there is no fixed graph spectrum for RGTV.
• Initialize weights for RGTV like GTV, and then updateweights and spectrum. We analyze the gradualtransformation of spectrum iteratively.

GSP Workshop 6/8/2018

56

Spectral domain

GSP Workshop 6/8/2018

57

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.

GSP Workshop 6/8/2018

58

Extensions• Accelerated algorithm for Gaussian blur:1𝐱ො = argminx(𝐈 + 𝑎 ⋅ 𝐋Γ )𝐱 − 𝐛 22 + 𝛽 𝐱 𝑅𝐺𝑇𝑉
2𝐱ො =𝑔(𝐋Γ )𝐛2
𝑔 𝐋Γ + 2𝛽 ⋅ 𝐋Γ= 𝐔Γ𝑔(𝚲Γ )
𝐔Γ𝑇 𝐛
2
𝑔 𝚲Γ + 2𝛽 ⋅ 𝚲ΓSolve it via accelerated graph filter
via Lanczos method [1][1] A. Susnjara, N. Perraudin, D. Kressner, and P. Vandergheynst, “Accelerated filtering on graphs using lanczos method,” arXiv preprint
arXiv:1509.04537, 2015.

59

Workflow

Blurry Image

Skeleton
Image
Reconstruction

Reconstruction

Kernel Estimation

GSP Workshop 6/8/2018

60

Experimental Results

GSP Workshop 6/8/2018

61

Experimental Results

GSP Workshop 6/8/2018

62Outline• GSP for Image Compression• Optimality of GFT• Generalized GFT• Signed GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV: 3D Point Cloud Denoising• Deep GLR• Ongoing & Future Work
GSP Workshop 6/8/2018

63

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

64

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    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
65& Total Variation,” accepted to IEEE 20th International Workshop on Multimedia Signal Processing, Vancouver, Canada, August 2018.

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.

66

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.
GSP Workshop 6/8/2018

67

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.
GSP Workshop 6/8/2018

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

GSP Workshop 6/8/2018

69

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

GSP Workshop 6/8/2018

70

Experimental Results – Numerical Comparison

GSP Workshop 6/8/201871Outline• GSP for Image Compression• Optimality of GFT• Generalized GFT• Signed GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV• Deep GLR• Ongoing & Future Work
GSP Workshop 6/8/201872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
Bilateral weights:

Q: what is the “most appropriate” graph?

 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.

73

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..

74

Unrolling Graph Laplacian Regularizer

75[1] J. Pang et al., “Deep Graph Laplacian Regularization,“ submitted to ECCV, Feb 2018.

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.
76[1] J. Pang et al., “Deep Graph Laplacian Regularization,“ submitted to ECCV, Feb 2018.

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 GCNNaussian denoiser: Residual learning of deep for image denoising,” TIP 2017.
[2] Marc Lebrun et al, “The noise clinic: a blind image denoising algorithm,” IPOL 2015.

77

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.

78

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.

DeepGLR79

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.

DeepGLR80

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.

DeepGLR81

Outline• GSP for Image Compression• Optimality of GFT• Generalized GFT• Signed GFT• GSP for Inverse Imaging• Graph Laplacian Regularizer• Reweighted Graph TV: 3D Point Cloud Denoising• Deep GLR• Ongoing & Future Work
GSP Workshop 6/8/201882

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.

GSP Workshop 6/8/2018

83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.

• Metric learning for edge weight computation for graphbased binary classifiers.
• Model-guided learning safeguard against worst-case /adversary noise?

GSP Workshop 6/8/2018

84

Q&A
• Email: gene.cheung@ieee.org• Homepage: http://research.nii.ac.jp/~cheung/

GSP Workshop 6/8/201885