Nonlocal image processing techniques:overview and
advance developments
Vladimir Katkovnik
Department of Signal Processing,Tampere University of Technology,Tampere,Finland
November 7,2009
ICIP 2009,Cairo,Egypt
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 1
/184
Part I:Basics,Modeling and Algorithms
1
From local to nonlocal approximations;
2
Nonlocal means,block matching ltering;
3
Highorder local and nonlocal models;
4
Block matching ltering with collaborative ltering (BM3D algorithm);
5
Redundancy and multiple model nonlocal approximations;
6
Applications:denoising,color image denoising,deblurring,demosaicking;
7
Development of blockmatching:shapeadaptive patches and adaptive PCA
transforms.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 2
/184
Introduction
In local image reconstruction algorithms usually use observations in a
neighborhood of a pixel of the interest;
In the nonlocal techniques,algorithms analyze data in largeand collects
the observations from the whole image looking for similar features;
The evolution of the nonlocal techniques from the sample nonlocal means
(NL) to the transform domain processing is an evolution to higherorder
models;
The latter algorithms are redundant,the data are processed by overlapping
blocks,and multiple estimates obtained for each pixel are fused (aggregated)
into the nal image estimates.
Katkovnik,V.,A.Foi,K.Egiazarian,and J.Astola,From local kernel to
nonlocal multiplemodel image denoising,to appear Int.J.Computer
Vision,2009,http://www.cs.tut./~lasip/#ref_papers.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 3
/184
Introduction
In local image reconstruction algorithms usually use observations in a
neighborhood of a pixel of the interest;
In the nonlocal techniques,algorithms analyze data in largeand collects
the observations from the whole image looking for similar features;
The evolution of the nonlocal techniques from the sample nonlocal means
(NL) to the transform domain processing is an evolution to higherorder
models;
The latter algorithms are redundant,the data are processed by overlapping
blocks,and multiple estimates obtained for each pixel are fused (aggregated)
into the nal image estimates.
Katkovnik,V.,A.Foi,K.Egiazarian,and J.Astola,From local kernel to
nonlocal multiplemodel image denoising,to appear Int.J.Computer
Vision,2009,http://www.cs.tut./~lasip/#ref_papers.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 3
/184
Introduction
In local image reconstruction algorithms usually use observations in a
neighborhood of a pixel of the interest;
In the nonlocal techniques,algorithms analyze data in largeand collects
the observations from the whole image looking for similar features;
The evolution of the nonlocal techniques from the sample nonlocal means
(NL) to the transform domain processing is an evolution to higherorder
models;
The latter algorithms are redundant,the data are processed by overlapping
blocks,and multiple estimates obtained for each pixel are fused (aggregated)
into the nal image estimates.
Katkovnik,V.,A.Foi,K.Egiazarian,and J.Astola,From local kernel to
nonlocal multiplemodel image denoising,to appear Int.J.Computer
Vision,2009,http://www.cs.tut./~lasip/#ref_papers.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 3
/184
Introduction
In local image reconstruction algorithms usually use observations in a
neighborhood of a pixel of the interest;
In the nonlocal techniques,algorithms analyze data in largeand collects
the observations from the whole image looking for similar features;
The evolution of the nonlocal techniques from the sample nonlocal means
(NL) to the transform domain processing is an evolution to higherorder
models;
The latter algorithms are redundant,the data are processed by overlapping
blocks,and multiple estimates obtained for each pixel are fused (aggregated)
into the nal image estimates.
Katkovnik,V.,A.Foi,K.Egiazarian,and J.Astola,From local kernel to
nonlocal multiplemodel image denoising,to appear Int.J.Computer
Vision,2009,http://www.cs.tut./~lasip/#ref_papers.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 3
/184
Table of Classication
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 4
/184
Image denoising:observation model
Let we have independent random observation pairs fz
i
,x
i
g
n
i =1
given for
simplicity in additive form
z
i
= y
i
+ε
i
,
where y
i
= y(x
i
) is a signal of interest,x
i
2 R
2
and ε
i
= ε(x
i
) is an additive
noise.
The denoising problem is to reconstruct y(x
i
) from fz
i
g
n
i =1
.
Variational and heuristic approaches.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 5
/184
Image denoising:variational approach
The variational approach formalizes an image reconstruction as an optimization
problem
ˆy = argmin
y
1
σ
2
jjz yjj
2

{z
}
delity
+µ pen(y)

{z
}
penalty
,µ > 0.
The delity term follows from a statistical noise model and the penalty is a prior
for y.
Parametric and nonparametric formulations.
Typical parametric models are of the form
y(x) =
∑
k
c
k
φ
k
(x).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 6
/184
Image denoising:heuristic approach (local means)
The weighted least square method gives the following criterion
J(x
0
,C) =
∑
s
w
h
(x
0
X
s
)e
2
s
,e
s
= z
s
C
0
,
where w is a signalindependent window function,
w
h
(x) = w(x/h),
and C
0
is an approximation of y(x) at x = x
0
.
Minimizing J on C
0
:
ˆ
C
0
(x) = arg min
C
0
J(x
0
,C
0
) =) ˆy
h
(x
0
) =
ˆ
C
0
(x) =
∑
s
w
h
(x
0
X
s
)z
s
∑
s
w
h
(x
0
X
s
)
.
Example of signalindependent weight
w
h
(x
0
X
s
) = e
jjx
0
X
s
jj
2
h
2
,h > 0.
Example of signaldependent weight
w
h
(x
0
X
s
,y
0
Y
s
) = e
jjx
0
Xs jj
2
h
2
jy
0
Ys j
2
γ
,γ,h > 0.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 7
/184
Image denoising:heuristic approach (local means)
The weighted least square method gives the following criterion
J(x
0
,C) =
∑
s
w
h
(x
0
X
s
)e
2
s
,e
s
= z
s
C
0
,
where w is a signalindependent window function,
w
h
(x) = w(x/h),
and C
0
is an approximation of y(x) at x = x
0
.
Minimizing J on C
0
:
ˆ
C
0
(x) = arg min
C
0
J(x
0
,C
0
) =) ˆy
h
(x
0
) =
ˆ
C
0
(x) =
∑
s
w
h
(x
0
X
s
)z
s
∑
s
w
h
(x
0
X
s
)
.
Example of signalindependent weight
w
h
(x
0
X
s
) = e
jjx
0
X
s
jj
2
h
2
,h > 0.
Example of signaldependent weight
w
h
(x
0
X
s
,y
0
Y
s
) = e
jjx
0
Xs jj
2
h
2
jy
0
Ys j
2
γ
,γ,h > 0.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 7
/184
Image denoising:heuristic approach (local means)
The weighted least square method gives the following criterion
J(x
0
,C) =
∑
s
w
h
(x
0
X
s
)e
2
s
,e
s
= z
s
C
0
,
where w is a signalindependent window function,
w
h
(x) = w(x/h),
and C
0
is an approximation of y(x) at x = x
0
.
Minimizing J on C
0
:
ˆ
C
0
(x) = arg min
C
0
J(x
0
,C
0
) =) ˆy
h
(x
0
) =
ˆ
C
0
(x) =
∑
s
w
h
(x
0
X
s
)z
s
∑
s
w
h
(x
0
X
s
)
.
Example of signalindependent weight
w
h
(x
0
X
s
) = e
jjx
0
X
s
jj
2
h
2
,h > 0.
Example of signaldependent weight
w
h
(x
0
X
s
,y
0
Y
s
) = e
jjx
0
Xs jj
2
h
2
jy
0
Ys j
2
γ
,γ,h > 0.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 7
/184
Image denoising:heuristic approach (local means)
The weighted least square method gives the following criterion
J(x
0
,C) =
∑
s
w
h
(x
0
X
s
)e
2
s
,e
s
= z
s
C
0
,
where w is a signalindependent window function,
w
h
(x) = w(x/h),
and C
0
is an approximation of y(x) at x = x
0
.
Minimizing J on C
0
:
ˆ
C
0
(x) = arg min
C
0
J(x
0
,C
0
) =) ˆy
h
(x
0
) =
ˆ
C
0
(x) =
∑
s
w
h
(x
0
X
s
)z
s
∑
s
w
h
(x
0
X
s
)
.
Example of signalindependent weight
w
h
(x
0
X
s
) = e
jjx
0
X
s
jj
2
h
2
,h > 0.
Example of signaldependent weight
w
h
(x
0
X
s
,y
0
Y
s
) = e
jjx
0
Xs jj
2
h
2
jy
0
Ys j
2
γ
,γ,h > 0.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 7
/184
Image denoising:heuristic approach (nonlocal means)
Criterion
J
h,x
0 (C) =
∑
s
w
h
(y
0
Y
s
)[z
s
C
0
]
2
,y
0
= y(x
0
),
where the weights w
h
depend on the distance between the signal values at
the observation points y
s
and the desirable point y
0
= y(x
0
).
Minimization,min
C
J
h,x
0 (C),gives the nonlocal means estimate
ˆy
h
(x
0
) =
ˆ
C
0
(x
0
) =
1
∑
s
w
h
(y
0
Y
s
)
∑
s
w
h
(y
0
Y
s
)z
s
,
where y
0
= y(x
0
).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 8
/184
Image denoising:heuristic approach (nonlocal means)
Criterion
J
h,x
0 (C) =
∑
s
w
h
(y
0
Y
s
)[z
s
C
0
]
2
,y
0
= y(x
0
),
where the weights w
h
depend on the distance between the signal values at
the observation points y
s
and the desirable point y
0
= y(x
0
).
Minimization,min
C
J
h,x
0 (C),gives the nonlocal means estimate
ˆy
h
(x
0
) =
ˆ
C
0
(x
0
) =
1
∑
s
w
h
(y
0
Y
s
)
∑
s
w
h
(y
0
Y
s
)z
s
,
where y
0
= y(x
0
).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 8
/184
Local versus Nonlocal Supports
Local versus nonlocal supports for zero and rst order polynomial tting:local
III;nonlocal zeroorder model I [II;nonlocal rstorder model I.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 9
/184
Local Pointwise Approximations
Three key slogans associated with these techniques:
Locality;
Anisotropy;
Adaptivity.
The locality means that there is a neighborhood where the image intensity
is well approximated some continuous basis functions.
The anisotropy means that a good local approximation can be achieved
only in a nonsymmetric neighborhood.
The adaptivity means that both the size and the shape should be data
adaptive.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 10
/184
Local Polynomial Approximation
In a neighborhood of x the Taylor series gives for y(X
s
):
y(X
s
)'y(x) y
(1)
x
1
(x)(x
1
X
1,s
) y
(1)
x
2
(x)(x
2
X
2,s
)+
y
(2)
x
1
,x
2
(x)(x
1
X
1,s
)(x
2
X
2,s
) +
1
2
y
(2)
x
1
(x)(x
1
X
1,s
)
2
+
1
2
y
(2)
x
2
(x)(x
2
X
2,s
)
2
....
Since the function y and the derivatives are unknown we look for
tting data y(X
s
) in the form
y(x,X
s
)'C
0
C
1
(x
1
X
1,s
) C
2
(x
2
X
2,s
) +
C
12
(x
1
X
1,s
)(x
2
X
2,s
) +
1
2
C
11
(x
1
X
1,s
)
2
+
1
2
C
22
(x
2
X
2,s
)
2
...
where the coe¢ cients C
0
,C
1
and C
2
give estimates for y(x),y
(1)
x
1
(x)
and y
(1)
x
2
(x),and C
12
,C
11
,C
22
give estimates for y
(2)
x
1
,x
2
(x),y
(2)
x
1
(x) and
y
(2)
x
2
(x).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 11
/184
How to formalize a local t?
The weighted least square method gives the following criterion
J(x,C) =
∑
s
w
h
(x X
s
)e
2
s
,
e
s
= y
s
y(x,X
s
),
where w is a window function,w
h
(x) = w(x/h),and
C = (C
0
,C
1
,....)
T
.
Minimizing J on C:
ˆ
C(x) = arg min
C
J(x,C).
The notation
ˆ
C(x) emphasizes that the estimate of C depends on x.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 12
/184
Estimates
The vector parameter
ˆ
C(x) immediately gives the estimates of the
function y and the derivatives y
(r )
:
ˆy(x) =
ˆ
C
0
(x),
ˆy
(1)
x
1
(x) =
ˆ
C
1
(x),
ˆy
(1)
x
2
(x) =
ˆ
C
2
(x),
.......
The conventional windows can be used:Kaiser,Hamming,Bartlett,
Blackman,Chebyshev,etc.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 13
/184
Convolutional estimates
Using the standard notation for multidimensional convolution,the
estimates can be represented in the following compact form
ˆy
h
(x) = (g
h
~z)(x) =
∑
X
1,s
,X
2,s
g
h
(x
1
X
1,s
,x
2
X
2,s
)z(X
1,s
,X
2,s
),
ˆy
(r )
h
(x) = (g
(r )
h
~z)(x),x 2 X
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 14
/184
2D examples (demo_CreateLPAKernels.m)
The smoothing kernel g
h
and its amplitude frequency characteristic jG
h
j:
Gaussian window,m = [2,2].
The lowpass lter with a peak of the frequency characteristic at ¯ω = 0.
Figure:
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 15
/184
2D examples (cont)
The di¤erentiation kernel g
(1,0)
h
and its amplitude frequency characteristic
jG
(1,0)
h
j:Gaussian window,m = [2,2].
The bandpass lter jG
(1,0)
h
j = 0 at ¯ω = 0.
Figure:
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 16
/184
2D examples (cont)
The di¤erentiation kernel g
(0,1)
h
and its amplitude frequency characteristic
jG
(0,1)
h
j:Gaussian window,m = [2,2].
The bandpass lter jG
(0,1)
h
j = 0 at ¯ω = 0.
Figure:
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 17
/184
2D examples (cont)
The di¤erentiation kernel g
(2,0)
h
and its amplitude frequency characteristic
jG
(2,0)
h
j:Gaussian window,m = [2,2].
The bandpass lter jG
(2,0)
h
j = 0 at ¯ω = 0.
Figure:
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 18
/184
2D examples (cont)
The di¤erentiation kernel g
(0,2)
h
and its amplitude frequency characteristic
jG
(0,2)
h
j:Gaussian window,m = [2,2].
The bandpass lter jG
(0,2)
h
j = 0 at ¯ω = 0.
Figure:
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 19
/184
2D examples
The di¤erentiation kernel g
(1,1)
h
and its amplitude frequency characteristic
jG
(1,1)
h
j:Gaussian window,m = [2,2].
The bandpass lter jG
(1,1)
h
j = 0 at ¯ω = 0.
Figure:
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 20
/184
ICI Adaptive Window Size
We use pointwise nonparametric regression methods known under a
generic name Lepski
0
s approach and developed by Lepski O.,Nemirovski
A.,Goldenshluger A.,Spokoiny V.
Overall,the algorithm searches for the largest local vicinity of the
point of estimation x where the LPA assumptions t well to the data.
The estimates are calculated for a few scales and compared.
The adaptive scale is dened as the largest one of those for which the
estimate does not di¤er signicantly from the estimates corresponding
to the smaller scales.
We use this methods in the form known as the intersection of
condence intervals (ICI ) rule.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 21
/184
ICI Adaptive Window Size
We use pointwise nonparametric regression methods known under a
generic name Lepski
0
s approach and developed by Lepski O.,Nemirovski
A.,Goldenshluger A.,Spokoiny V.
Overall,the algorithm searches for the largest local vicinity of the
point of estimation x where the LPA assumptions t well to the data.
The estimates are calculated for a few scales and compared.
The adaptive scale is dened as the largest one of those for which the
estimate does not di¤er signicantly from the estimates corresponding
to the smaller scales.
We use this methods in the form known as the intersection of
condence intervals (ICI ) rule.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 21
/184
ICI Adaptive Window Size
We use pointwise nonparametric regression methods known under a
generic name Lepski
0
s approach and developed by Lepski O.,Nemirovski
A.,Goldenshluger A.,Spokoiny V.
Overall,the algorithm searches for the largest local vicinity of the
point of estimation x where the LPA assumptions t well to the data.
The estimates are calculated for a few scales and compared.
The adaptive scale is dened as the largest one of those for which the
estimate does not di¤er signicantly from the estimates corresponding
to the smaller scales.
We use this methods in the form known as the intersection of
condence intervals (ICI ) rule.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 21
/184
ICI Adaptive Window Size
We use pointwise nonparametric regression methods known under a
generic name Lepski
0
s approach and developed by Lepski O.,Nemirovski
A.,Goldenshluger A.,Spokoiny V.
Overall,the algorithm searches for the largest local vicinity of the
point of estimation x where the LPA assumptions t well to the data.
The estimates are calculated for a few scales and compared.
The adaptive scale is dened as the largest one of those for which the
estimate does not di¤er signicantly from the estimates corresponding
to the smaller scales.
We use this methods in the form known as the intersection of
condence intervals (ICI ) rule.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 21
/184
LPA accuracy
e
ˆy
h
(x)
= y(x) ˆy
h
(x) = y(x) (g
h
~(y +ε))(x) =
y(x) (g
h
~y)(x)

{z
}
bias
(g
h
~ε)(x)

{z
}
random error
.
The variance
σ
2
ˆy
h
= σ
2
∑
s
g
2
h
(X
s
).
The ICI rule gives h minimizing MSE
l
ˆy
h
= Efe
2
ˆy
h
(x)
g = bias
2
+σ
2
ˆy
h
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 22
/184
ICI rule (minimum MSE window size)
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 23
/184
ICI rule (minimum MSE window size)
H = fh
1
< h
2
<...< h
J
g.
D
j
= [ ˆy
h
j
(x) Γ σ
ˆy
h
(x),ˆy
h
j
(x) +Γ σ
ˆy
h
j
(x)],j = 1,...,J,
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 24
/184
A layout of the adaptive scale LPAICI algorithm
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 25
/184
Anisotropy:starshaped neighborhood
Consider a ball (disk) of the radius h dening a spherical neighborhood of
x,
B
h
= fu:jjx ujj hg.
Introduce a sectorial partition of this ball with K nonoverlapping sectors
having x as a common vertex.
The adaptivity ICI technique is used to nd the varying scale for each
sector independently.We obtain K estimates ˆy
h,θ
i
(x),i = 1,...,K,with K
nonoverlapping supports S
h
θ
i
covering the ball B
h
.
The union of the supports of these sectors is a starshaped set.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 26
/184
Adaptive size sectorial neighborhoods are obtained by the
LPAICI algorithm
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 27
/184
These are ICI adaptive directional window sizes calculated for 8 directions
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 28
/184
These ICI adaptive directional estimates are fused in the nal one using
the multiwindow method
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 29
/184
Fusing (aggregation) of directional estimates
Using K sectors we obtain K independently derived estimates for each x.
The multi window (fused) estimates are exploited to obtain the unique
nal estimate ˆy(x) from the partial directional ones.
With the inverse variances as the weights for the linear fusing it gives
ˆy(x) =
∑
K
j =1
λ
j
ˆy
θ
j
(x),ˆy
θ
j
(x) = ˆy
h,θ
j
(x)j
h=h
+
(x,θ
j
)
,
λ
j
= σ
2
j
(x)/
∑
K
i =1
σ
2
i
(x),
where
ˆy
h,θ
j
(x) = (g
h,θ
j
~z)(x),σ
2
j
(x) = σ
2
ˆy
,θ
j
(x,h)j
h=h
+
(x,θ
j
)
,
σ
2
ˆy,
θ
j
(x,h) = σ
2
∑
x
g
2
h,θ
j
(x).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 30
/184
Fusing (aggregation) of directional estimates
The weights λ
j
are datadriven adaptive as σ
2
j
(x) depend on the
adaptive pointwise h
+
(x,θ
j
).
Assuming that the supports of the kernels g
θ
j
,h
are not overlapping and
neglecting that the kernels have the point x as a common one we obtain
for the variance of the fused estimate
σ
2
ˆy
(x) =
∑
K
j =1
λ
2
j
σ
2
j
(x) =
1
∑
K
j =1
σ
2
j
(x)
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 31
/184
ShapeAdaptive DCT Transform (Local Multipoint)
This approach to estimation for a point x can be roughly described as the
following four stage procedure.
Stage I (spatial adaptation):For every x 2 X,dene a neighborhood
˜
U
+
x
of x where a simple loworder polynomial model ts the data;
Stage II (order selection):apply some localized transform (parametric
series model) to the data on the set
˜
U
+
x
,use thresholding operator
(model selection procedure) in order to identify the signicant (i.e.
nonzero) elements of the transform (and thus the order of the
parametric model).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 32
/184
ShapeAdaptive DCT Transform (Local Multipoint)
This approach to estimation for a point x can be roughly described as the
following four stage procedure.
Stage I (spatial adaptation):For every x 2 X,dene a neighborhood
˜
U
+
x
of x where a simple loworder polynomial model ts the data;
Stage II (order selection):apply some localized transform (parametric
series model) to the data on the set
˜
U
+
x
,use thresholding operator
(model selection procedure) in order to identify the signicant (i.e.
nonzero) elements of the transform (and thus the order of the
parametric model).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 32
/184
ShapeAdaptive DCT Transform
Stage III (approximation):Calculate,by inversetransformation of the
signicant elements only,the corresponding estimates ˆy
˜
U
+
x
(v) of the
signal for all v 2
˜
U
+
x
.These ˆy
˜
U
x
are calculated for all x 2 X.
Stage IV (aggregation):Let x 2 X and I
x
=
x 2 X:x 2
˜
U
+
x
be
the set of the centers of the neighborhoods which have x as a
common point.The nal estimate ˆy(x) is calculated as an aggregate
of
ˆy
˜
U
+
x
(x)
x2I
x
.
One key aspect in this procedure is that by demanding the local t of
a loworder polynomial model,we are able to avoid the presence of
singularities,discontinuities,or sharp transitions within the transform
support
˜
U
+
x
.In this way,we increase the sparsity in the transform
domain,improving the e¤ectiveness of thresholding.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 33
/184
ShapeAdaptive DCT Transform
Stage III (approximation):Calculate,by inversetransformation of the
signicant elements only,the corresponding estimates ˆy
˜
U
+
x
(v) of the
signal for all v 2
˜
U
+
x
.These ˆy
˜
U
x
are calculated for all x 2 X.
Stage IV (aggregation):Let x 2 X and I
x
=
x 2 X:x 2
˜
U
+
x
be
the set of the centers of the neighborhoods which have x as a
common point.The nal estimate ˆy(x) is calculated as an aggregate
of
ˆy
˜
U
+
x
(x)
x2I
x
.
One key aspect in this procedure is that by demanding the local t of
a loworder polynomial model,we are able to avoid the presence of
singularities,discontinuities,or sharp transitions within the transform
support
˜
U
+
x
.In this way,we increase the sparsity in the transform
domain,improving the e¤ectiveness of thresholding.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 33
/184
ShapeAdaptive DCT Transform
Stage III (approximation):Calculate,by inversetransformation of the
signicant elements only,the corresponding estimates ˆy
˜
U
+
x
(v) of the
signal for all v 2
˜
U
+
x
.These ˆy
˜
U
x
are calculated for all x 2 X.
Stage IV (aggregation):Let x 2 X and I
x
=
x 2 X:x 2
˜
U
+
x
be
the set of the centers of the neighborhoods which have x as a
common point.The nal estimate ˆy(x) is calculated as an aggregate
of
ˆy
˜
U
+
x
(x)
x2I
x
.
One key aspect in this procedure is that by demanding the local t of
a loworder polynomial model,we are able to avoid the presence of
singularities,discontinuities,or sharp transitions within the transform
support
˜
U
+
x
.In this way,we increase the sparsity in the transform
domain,improving the e¤ectiveness of thresholding.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 33
/184
Illustrations
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 34
/184
ShapeAdaptive DCT Transform
Illustration of the shapeadaptive DCT transform and its inverse.
Transformation is computed by cascaded application of onedimensional
varyinglength DCT transforms,along the columns and along the rows.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 35
/184
Multipoint Estimate Fusing (fusingaggregation)
In order to obtain a single global estimate ˆy:X!R dened on the whole
image domain,all the local patch estimates are averaged together using
adaptive weights w
x
2 R in the following convex combination:
ˆy =
∑
x2X
w
x
ˆy
˜
U
+
x
jX
∑
x2X
w
x
χ
˜
U
+
x
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 36
/184
Implementation of SADCT
http://www.cs.tut./~foi/SADCT/
demo_SADCT_denoising.m
demo_SADCT_color_denoising.m
demo_SADCT_deblurring.m
demo_SADCT_deblocking.m
demo_SADCT_inverse_halftoning.m
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 37
/184
Aggregation (multipoint) in action
A crosssection of length 31 pixels from Peppers testimage ( σ = 25).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 38
/184
NLmeans lters
A.Buades,et al.,A review of image denoising algorithms,with a new one,
SIAM Multiscale Modeling and Simulation,vol.4,2005.
The nonlocal means (NLmeans) as they are introduced by Buades et al.
(2005) have been given in a di¤erent form where these weights calculated
over spatial neighborhoods of the points x
0
and x
s
.
These neighborhoodwise di¤erences can be interpreted as more reliable way
to estimate y
0
y
s
from the noise samples alone.
Then,the nonlocal mean estimate is calculated in a pointwise manner as the
weighted mean with the weights dened by the proximity measure between
the image patches used in the estimate.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 39
/184
NLmeans lters
A.Buades,et al.,A review of image denoising algorithms,with a new one,
SIAM Multiscale Modeling and Simulation,vol.4,2005.
The nonlocal means (NLmeans) as they are introduced by Buades et al.
(2005) have been given in a di¤erent form where these weights calculated
over spatial neighborhoods of the points x
0
and x
s
.
These neighborhoodwise di¤erences can be interpreted as more reliable way
to estimate y
0
y
s
from the noise samples alone.
Then,the nonlocal mean estimate is calculated in a pointwise manner as the
weighted mean with the weights dened by the proximity measure between
the image patches used in the estimate.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 39
/184
NLmeans lters
A.Buades,et al.,A review of image denoising algorithms,with a new one,
SIAM Multiscale Modeling and Simulation,vol.4,2005.
The nonlocal means (NLmeans) as they are introduced by Buades et al.
(2005) have been given in a di¤erent form where these weights calculated
over spatial neighborhoods of the points x
0
and x
s
.
These neighborhoodwise di¤erences can be interpreted as more reliable way
to estimate y
0
y
s
from the noise samples alone.
Then,the nonlocal mean estimate is calculated in a pointwise manner as the
weighted mean with the weights dened by the proximity measure between
the image patches used in the estimate.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 39
/184
NLmeans lters (cont.)
This estimation can be formalized as minimization of the weighted squared
residual criterion
I
h,x
(C) =
∑
s
w
h,s
(x
0
,x
s
)[z
s
C]
2
,
with,say,Gaussian weights
w
h,s
(x
0
,x
s
) = e
∑
v2V
(z(x
0
+v )z(x
s
+v))
2
h
2
dened by the Euclidean distance between the observations z in
Vneighborhoods of the points x
0
and x
s
,V being a xed neighborhood of
x
0
.
The nonlocal means estimate is calculated as the weighted mean
ˆy
h
(x
0
) =
∑
s
g
h,s
(x
0
)z
s
,g
h,s
(x
0
) =
w
h,s
(x
0
,x
s
)
∑
s
w
h,s
(x
0
,x
s
)
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 40
/184
NLmeans lters (cont.)
This estimation can be formalized as minimization of the weighted squared
residual criterion
I
h,x
(C) =
∑
s
w
h,s
(x
0
,x
s
)[z
s
C]
2
,
with,say,Gaussian weights
w
h,s
(x
0
,x
s
) = e
∑
v2V
(z(x
0
+v )z(x
s
+v))
2
h
2
dened by the Euclidean distance between the observations z in
Vneighborhoods of the points x
0
and x
s
,V being a xed neighborhood of
x
0
.
The nonlocal means estimate is calculated as the weighted mean
ˆy
h
(x
0
) =
∑
s
g
h,s
(x
0
)z
s
,g
h,s
(x
0
) =
w
h,s
(x
0
,x
s
)
∑
s
w
h,s
(x
0
,x
s
)
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 40
/184
Transform domain NLmeans:multipoint estimates
We consider nonlocal estimates by use of the transforms enabling the adaptive
highorder approximations of the windowed data.
Let the image be dened on a regular 2D grid X.Consider a windowing
C = fX
r
,r = 1,...,N
s
g of X with N
s
blocks (uniform windows)
X
r
X of size n
r
n
r
such that [
N
s
r =1
X
r
= X.
The noisefree data y (x) and the noisy data z(x) windowed on X
r
are
arranged in n
r
n
r
blocks denoted as Y
r
and Z
r
,respectively.
Typically,the blocks may overlap and we use transforms in conjunction with
the concept of the redundancy of natural signals.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 41
/184
Transform domain NLmeans:multipoint estimates
We consider nonlocal estimates by use of the transforms enabling the adaptive
highorder approximations of the windowed data.
Let the image be dened on a regular 2D grid X.Consider a windowing
C = fX
r
,r = 1,...,N
s
g of X with N
s
blocks (uniform windows)
X
r
X of size n
r
n
r
such that [
N
s
r =1
X
r
= X.
The noisefree data y (x) and the noisy data z(x) windowed on X
r
are
arranged in n
r
n
r
blocks denoted as Y
r
and Z
r
,respectively.
Typically,the blocks may overlap and we use transforms in conjunction with
the concept of the redundancy of natural signals.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 41
/184
Transform domain NLmeans (cont.)
Mainly these are the 2D discrete Fourier and cosine transforms (DFT and
DCT),orthogonal polynomials,and wavelet transforms.The transform,
denoted as T
2D
r
,is applied for each window X
r
independently as
θ
r
= T
2D
r
(Y
r
),
h
θ
r
= D
r
Y
r
D
T
r
i
r = 1,...,N
s
,
where θ
r
is the spectrum of Y
r
,and D
r
are orthonormal matrices
D
r
D
T
r
= D
T
r
D
r
= I.
The equality enclosed in square brackets holds when the transform T
2D
r
is
realized as a separable composition of 1D transforms,each computed by
matrix multiplication against an n
r
n
r
orthogonal matrix D
r
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 42
/184
Transform domain NLmeans (cont.)
Mainly these are the 2D discrete Fourier and cosine transforms (DFT and
DCT),orthogonal polynomials,and wavelet transforms.The transform,
denoted as T
2D
r
,is applied for each window X
r
independently as
θ
r
= T
2D
r
(Y
r
),
h
θ
r
= D
r
Y
r
D
T
r
i
r = 1,...,N
s
,
where θ
r
is the spectrum of Y
r
,and D
r
are orthonormal matrices
D
r
D
T
r
= D
T
r
D
r
= I.
The equality enclosed in square brackets holds when the transform T
2D
r
is
realized as a separable composition of 1D transforms,each computed by
matrix multiplication against an n
r
n
r
orthogonal matrix D
r
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 42
/184
Transform domain NLmeans (cont.)
The inverse T
2D
1
r
of T
2D
r
denes the signal from the spectrum as
Y
r
= T
2D
1
r
(θ
r
),
h
Y
r
= D
T
r
θ
r
D
r
i
r = 1,...,N
s
.
The noisy spectrum of the noisy signal is dened as
˜
θ
r
= T
2D
r
(
Z
r
)
,
h
˜
θ
r
= D
r
Z
r
D
T
r
i
r = 1,...,N
s
.
The signal y is sparse if it can be well approximated by a small number of
nonzero elements of the spectrum θ
r
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 43
/184
Transform domain NLmeans (cont.)
The inverse T
2D
1
r
of T
2D
r
denes the signal from the spectrum as
Y
r
= T
2D
1
r
(θ
r
),
h
Y
r
= D
T
r
θ
r
D
r
i
r = 1,...,N
s
.
The noisy spectrum of the noisy signal is dened as
˜
θ
r
= T
2D
r
(
Z
r
)
,
h
˜
θ
r
= D
r
Z
r
D
T
r
i
r = 1,...,N
s
.
The signal y is sparse if it can be well approximated by a small number of
nonzero elements of the spectrum θ
r
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 43
/184
Transform domain NLmeans (cont.)
The inverse T
2D
1
r
of T
2D
r
denes the signal from the spectrum as
Y
r
= T
2D
1
r
(θ
r
),
h
Y
r
= D
T
r
θ
r
D
r
i
r = 1,...,N
s
.
The noisy spectrum of the noisy signal is dened as
˜
θ
r
= T
2D
r
(
Z
r
)
,
h
˜
θ
r
= D
r
Z
r
D
T
r
i
r = 1,...,N
s
.
The signal y is sparse if it can be well approximated by a small number of
nonzero elements of the spectrum θ
r
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 43
/184
Transform domain NLmeans calculation:single model
Let Z
r
be a reference window and Z
s
be all others.
The nonlocal estimate for rth widow is formalized as minimization of the
weighted criterion
I
h,Z
r
(θ) =
∑
s
w
h
(Z
r
,Z
s
)jjZ
s
D
T
θDjj
2
2
,
with,say,Gaussian weights
w
h
(Z
r
,Z
s
) = exp(jjZ
r
Z
s
jj
2
2
/h
2
).
The estimate for r th reference window is
ˆ
θ
r
= arg min
θ
I
h,Z
r
(θ).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 44
/184
Transform domain NLmeans calculation:single model
(cont)
The minimum conditions
∂
∂θ
I
h,Z
r
(θ) = 0 =)
∑
s
w
h
(Z
r
,Z
s
)(DZ
s
D
T
θ) = 0 =)
ˆ
θ
r
=
∑
s
w
h
(Z
r
,Z
s
)
˜
θ
s
∑
s
w
h
(Z
r
,Z
s
)
,
˜
θ
r
= DZ
r
D
T
.
The solution is the weighted mean of the noisy spectrums
˜
θ
s
versus the weighted
mean of the means in the standard NLmeans.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 45
/184
Transform domain NLmeans:single model with prior
(penalty)
All denoising is because of this weighted mean.It can be essentially improved
imposing a prior on the spectrum θ.
The nonlocal estimate for rth widow is formalized as minimization of the
weighted criterion
I
h,Z
r
(θ) =
1
σ
2
∑
s
w
h
(Z
r
,Z
s
)jjZ
s
D
T
θDjj
2
2
+µ
2
pen(θ).
pen(θ) = jjθjj
0
is l
0
norm is equal to a number of nonzero items in the
matrix θ.jjθjj
0
is a measure of the model complexity.
Another norms also can be used jjθjj
1
=
∑
jθ
ij
j,jjθjj
2
=
q
∑
θ
2
ij
.
A minimal complexity estimate is calculated as
ˆ
θ
r
= arg min
θ
1
σ
2
∑
s
w
h
(Z
r
,Z
s
)jj
˜
θ
s
θjj
2
2
+µ
2
jjθjj
0
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 46
/184
Transform domain NLmeans:single model with prior
(penalty)
All denoising is because of this weighted mean.It can be essentially improved
imposing a prior on the spectrum θ.
The nonlocal estimate for rth widow is formalized as minimization of the
weighted criterion
I
h,Z
r
(θ) =
1
σ
2
∑
s
w
h
(Z
r
,Z
s
)jjZ
s
D
T
θDjj
2
2
+µ
2
pen(θ).
pen(θ) = jjθjj
0
is l
0
norm is equal to a number of nonzero items in the
matrix θ.jjθjj
0
is a measure of the model complexity.
Another norms also can be used jjθjj
1
=
∑
jθ
ij
j,jjθjj
2
=
q
∑
θ
2
ij
.
A minimal complexity estimate is calculated as
ˆ
θ
r
= arg min
θ
1
σ
2
∑
s
w
h
(Z
r
,Z
s
)jj
˜
θ
s
θjj
2
2
+µ
2
jjθjj
0
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 46
/184
Transform domain NLmeans:single model with prior
(penalty)
All denoising is because of this weighted mean.It can be essentially improved
imposing a prior on the spectrum θ.
The nonlocal estimate for rth widow is formalized as minimization of the
weighted criterion
I
h,Z
r
(θ) =
1
σ
2
∑
s
w
h
(Z
r
,Z
s
)jjZ
s
D
T
θDjj
2
2
+µ
2
pen(θ).
pen(θ) = jjθjj
0
is l
0
norm is equal to a number of nonzero items in the
matrix θ.jjθjj
0
is a measure of the model complexity.
Another norms also can be used jjθjj
1
=
∑
jθ
ij
j,jjθjj
2
=
q
∑
θ
2
ij
.
A minimal complexity estimate is calculated as
ˆ
θ
r
= arg min
θ
1
σ
2
∑
s
w
h
(Z
r
,Z
s
)jj
˜
θ
s
θjj
2
2
+µ
2
jjθjj
0
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 46
/184
Transform domain NLmeans:single model with prior
(penalty)
All denoising is because of this weighted mean.It can be essentially improved
imposing a prior on the spectrum θ.
The nonlocal estimate for rth widow is formalized as minimization of the
weighted criterion
I
h,Z
r
(θ) =
1
σ
2
∑
s
w
h
(Z
r
,Z
s
)jjZ
s
D
T
θDjj
2
2
+µ
2
pen(θ).
pen(θ) = jjθjj
0
is l
0
norm is equal to a number of nonzero items in the
matrix θ.jjθjj
0
is a measure of the model complexity.
Another norms also can be used jjθjj
1
=
∑
jθ
ij
j,jjθjj
2
=
q
∑
θ
2
ij
.
A minimal complexity estimate is calculated as
ˆ
θ
r
= arg min
θ
1
σ
2
∑
s
w
h
(Z
r
,Z
s
)jj
˜
θ
s
θjj
2
2
+µ
2
jjθjj
0
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 46
/184
Transform domain NLmeans:single model with prior,
single window
For a single window the minimal complexity model is calculated as
ˆ
θ
r
= arg min
θ
1
σ
2
jj
˜
θ
r
θjj
2
2
+µ
2
jjθjj
0
.
It can be shown that this solution gives the socalled hardthresholding (ltering)
ˆ
θ
r
(i,j ) =
˜
θ
r
(i,j ),if j
˜
θ
r
(i,j )j σµ,
0,if j
˜
θ
r
(i,j )j < σµ.
We use the symbol Υ() for the thresholding operation
ˆ
θ
r
(i,j ) = Υ(
˜
θ
r
(i,j )).
When
ˆ
θ
r
found the windowed estimate is calculated as
ˆ
Y
r
= D
T
ˆ
θ
r
D.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 47
/184
Transform domain NLmeans biasedness and variance
Let us evaluate the e¢ ciency of the thresholding ltering.
For the noisy data we have (before ltering)
Ef
˜
θ
r
g = EfD
r
Z
r
D
T
r
g = D
r
EfZ
r
gD
T
r
= D
r
Y
r
D
T
r
= θ
r
,
thus θ
r
Ef
˜
θ
r
g = 0.
For the hardthresholded data
Ef
ˆ
θ
r
(i,j )g =
θ
r
(i,j ),if j
˜
θ
r
(i,j )j σµ,
0,if j
˜
θ
r
(i,j )j < σµ,
=)
θ
r
(i,j ) Ef
ˆ
θ
r
(i,j )g =
0,if j
˜
θ
r
(i,j )j σµ,
θ
r
,if j
˜
θ
r
(i,j )j < σµ.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 48
/184
Transform domain NLmeans biasedness and variance
(cont)
The variance of the thresholded data
Δ
r
θ =
˜
θ
r
θ
r
= D
r
Z
r
D
T
r
D
r
Y
r
D
T
r
= D
r
E
r
D
T
r
;
S
var
=
∑
i,j
var f
˜
θ
r
(i,j )g =
∑
i,j
E(Δ
r
θ(i,j ))
2
= σ
2
n
2
r
.
S
var,thr
=
∑
i,j
E(
ˆ
θ
r
(i,j ))
2
= σ
2
(n
2
r
n
2
r,0
) = σ
2
N
x
r
har
,
where n
2
r,0
is a number of j
ˆ
θ
r
j = 0,N
x
r
har
is the number of retained (nonzero)
coe¢ cients after hardthresholding.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 49
/184
Transform domain NLmeans with prior:single model,
many windows
It can be veried that
I
h,Z
r
(θ) =
1
σ
2
∑
s
w
h
(Z
r
,Z
s
)jjZ
s
D
T
θDjj
2
2
+µ
2
pen(θ) can be rewritten
as
I
h,Z
r
(θ) =
1
σ
2
∑
s
w
h
(Z
r
,Z
s
)jj
¯
θ
r
θjj
2
2
+µ
2
pen(θ) +const,
where
¯
θ
r
=
∑
s
w
h
(Z
r
,Z
s
)
˜
θ
s
/
∑
s
w
h
(Z
r
,Z
s
).
Then the minimization of I
h,Z
r
(θ) and estimation for the reference block becomes
ˆ
θ = arg min
θ
jj
¯
θ
r
θjj
2
2
+µ
2
r
pen(θ),
ˆ
Y
r
= T
2D
1
ˆ
θ
.
where µ
2
r
= σ
2
µ
2
/
∑
s
w
h
(Z
r
,Z
s
).
If pen(θ) = jjθjj
0
the estimate
ˆ
θ is a thresholded version of
¯
θ
r
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 50
/184
Block matching and 3D ltering (BM3D):nonlocal
ltering with multiple models
K.Dabov,et al.,Image denoising by sparse 3D transformdomain
collaborative ltering,IEEE Tran.IP,vol.16,no.8,2007.

Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 51
/184
Stage I (block matching,grouping)
The grouping is a concept of collecting similar ddimensional fragments of a
given signal into a d +1dimensional data structure that we term group.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 52
/184
Distance between the blocks
d (Z
r
,Z
s
) =
Υ
0
˜
θ
r
Υ
0
˜
θ
s
2
2
,
where Υ
0
is the hardthresholding operator.
The grouping rule
S
ht
r
=
f
x 2 X:d
(
Z
r
,Z
s
)
τ
g
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 53
/184
Stage II (collaborative 3D ltering)
A collaborative ltering is realized as shrinkage in a transform domain.
Assuming 2 +1dimensional groups of similar signal fragments are already
formed,this approach comprises of the following steps.
Apply a 2 +1dimensional linear transform on a group.
Shrink (e.g.by soft and hardthresholding,Wiener ltering) the transform
coe¢ cients to attenuate the noise.
Invert the linear transform to produce estimates of all grouped fragments.
Return the ltered fragments to their original locations.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 54
/184
Stage II (collaborative 3D ltering)
A collaborative ltering is realized as shrinkage in a transform domain.
Assuming 2 +1dimensional groups of similar signal fragments are already
formed,this approach comprises of the following steps.
Apply a 2 +1dimensional linear transform on a group.
Shrink (e.g.by soft and hardthresholding,Wiener ltering) the transform
coe¢ cients to attenuate the noise.
Invert the linear transform to produce estimates of all grouped fragments.
Return the ltered fragments to their original locations.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 54
/184
Stage II (collaborative 3D ltering)
A collaborative ltering is realized as shrinkage in a transform domain.
Assuming 2 +1dimensional groups of similar signal fragments are already
formed,this approach comprises of the following steps.
Apply a 2 +1dimensional linear transform on a group.
Shrink (e.g.by soft and hardthresholding,Wiener ltering) the transform
coe¢ cients to attenuate the noise.
Invert the linear transform to produce estimates of all grouped fragments.
Return the ltered fragments to their original locations.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 54
/184
Stage II (collaborative 3D ltering)
A collaborative ltering is realized as shrinkage in a transform domain.
Assuming 2 +1dimensional groups of similar signal fragments are already
formed,this approach comprises of the following steps.
Apply a 2 +1dimensional linear transform on a group.
Shrink (e.g.by soft and hardthresholding,Wiener ltering) the transform
coe¢ cients to attenuate the noise.
Invert the linear transform to produce estimates of all grouped fragments.
Return the ltered fragments to their original locations.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 54
/184
Stage II (collaborative 3D ltering)
A collaborative ltering is realized as shrinkage in a transform domain.
Assuming 2 +1dimensional groups of similar signal fragments are already
formed,this approach comprises of the following steps.
Apply a 2 +1dimensional linear transform on a group.
Shrink (e.g.by soft and hardthresholding,Wiener ltering) the transform
coe¢ cients to attenuate the noise.
Invert the linear transform to produce estimates of all grouped fragments.
Return the ltered fragments to their original locations.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 54
/184
Collaborative ltering implementation
The collaborative ltering of Z
S
ht
xr
is realized by hardthresholding in 3D transform
domain.
The 3D linear transform,denoted T
ht
3D
,is expected to take advantage of the
sparsity for the true signal group Y
S
ht
xr
.
This allows for e¤ective noise attenuation by hardthresholding,followed by
inverse transform that yields a 3D array of blockwise estimates
b
Y
ht
S
ht
x
r
= T
ht
1
3D
Υ
T
ht
3D
Z
S
ht
x
r
,
where Υ is a hardthreshold operator with threshold µ
3D
σ.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 55
/184
Stage III (fusion/aggregation of ltered fragments)
The estimates of Stage II are heavily overlapping.
Calculate the pointwise estimates using the variances of the estimates in the
fragments.
The nal estimate by is computed by a weighted average of the blockwise
estimates
b
Y
ht,x
r
x2S
ht
x
r
using the special weights w
ht
x
R
by
basic
(x) =
∑
x
r
2X
∑
x
m
2S
ht
x
r
w
ht
x
r
b
Y
ht,x
r
x
m
(x)
∑
x
r
2X
∑
x
m
2S
ht
x
r
w
ht
x
r
χ
x
m
(x)
,8x 2 X,
w
ht
x
r
=
(
1
σ
2
N
xr
har
,if N
x
r
har
1
1,otherwise
where χ
x
m
:X!f0,1g is the characteristic function of the square support of a
block located at x
m
2 X,and the blockwise estimates
b
Y
ht,x
r
x
m
are zeropadded
outside of their support,where N
x
r
har
is the number of retained (nonzero)
coe¢ cients after hardthresholding.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 56
/184
BM3D Algorithm Flowchart
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 57
/184
Features of the Wieners Part of BM3D
(1) Hardthresholding is replaced by the Wiener ltering.We dene the empirical
Wiener shrinkage coe¢ cients as
W
S
wie
xr
=
T
wie
3D
b
Y
basic
S
wie
x
r
2
T
wie
3D
b
Y
basic
S
wie
xr
2
+σ
2
,
b
Y
wie
S
wie
xr
= T
wie
1
3D
W
S
wie
xr
T
wie
3D
Z
S
wie
xr
.
(2) The weights in the aggregation are calculated as
by
nal
(
x
)
=
∑
x
r
2X
∑
x
m
2S
wie
x
r
w
wie
x
r
b
Y
wie
S
wie
x
m
∑
x
r
2X
∑
x
m
2S
wie
x
r
w
wie
x
r
χ
x
m
(x)
,8x 2 X,w
wie
x
r
= σ
2
W
S
wie
x
r
2
2
.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 58
/184
Examples of performance
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 59
/184
Sparsity of Collaborative 3D HardThresholding
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 60
/184
Nonlocal BM3D vs SADCT:value of nonlocality
10
15
20
25
29
30
31
32
33
34
35
PSNR (dB)
10
15
20
25
32
33
34
35
36
PSNR (dB)
10
15
20
25
32
33
34
35
36
PSNR (dB)
σ σ σ
(a) Barbara (b) Lena (c) House
PSNR as a function of σ.The proposed method ("squares"),BLSGSM
(circles),exemplarbased denoising ("x"),KSVD denoising ("diamonds"),and
pointwise SADCT denoising ("stars").
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 61
/184
Dependency of the output PSNR (dB) on the used
transforms
Transform
Boats
Lena
σ = 25
T
ht
2D
T
wie
2D
T
1D
T
ht
2D
T
wie
2D
T
1D
Haar
29.31
29.84
29.91
31.24
31.93
32.08
Db2
29.22
29.83
29.90
31.19
31.97
32.06
Db4
29.34
29.88
29.89
31.31
32.01
32.06
Db6
29.30
29.86
29.89
31.28
31.98
32.06
Bior1.3
29.42
29.87
29.90
31.35
31.96
32.06
Bior1.5
29.43
29.88
29.90
31.37
31.97
32.06
WHT
29.22
29.84
29.88
31.24
32.00
32.07
DCT
29.35
29.91
29.88
31.42
32.08
32.07
DST
29.33
29.91
29.79
31.36
31.97
31.92
DC+rand
29.07
29.75
29.88
31.06
31.88
32.06
DConly


28.03


30.65
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 62
/184
Experimental comparison
The algorithms are applied on a set of 10 di¤erent test images corrupted by
additive white Gaussian noise with standard deviations σ = 5,15,20,25,35.
The comparison is made in terms of both PSNR and mean structural similarity
index map (MSSIM).
Z.Wang,et al."Image quality assessment:From error visibility to structural
similarity,"IEEE Tran.IP.13,no.4,2004.

We can see that three algorithms based on the collaborative ltering paradigm
occupy the topthree places also in this comparison.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 63
/184
Experimental comparison
PSNR results obtained for 10 di¤erent images with 5 di¤erent levels of noise.The
numbers are the results obtained by the BM3DSAPCA algorithm.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 64
/184
Experimental comparison
MSSIM results obtained for 10 di¤erent images with 5 di¤erent levels of noise.
The numbers are the results obtained by the BM3DSAPCA algorithm.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 65
/184
Image Denoising with BM3D
http://www.cs.tut./~foi/GCFBM3D/
(1) BM3D.m
[PSNR,y_est] = BM3D(y,z,sigma,prole,print_to_screen)
(2) CBM3D.m ( color version)
[PSNR,yRGB_est] = CBM3D(yRGB,zRGB,sigma,prole,print_to_screen,
colorspace).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 66
/184
Color denoising
The key idea is the following:the structures (e.g.objects,edges,details,etc.)
which determine the adaptive shapes are the same across all three color channels.
The method is implemented after transformation in the"opponent"(luminance
and chrominances) colorspace
0
@
y
1
y
2
y
3
1
A
=
0
B
@
1
3
1
3
1
3
1
p
6
0
1
p
6
1
3
p
2
2
3
p
2
1
3
p
2
1
C
A
0
@
y
r
y
g
y
b
1
A
.
For all three color channels the same grouping is used dened for the luminance
channel.
K.Dabov,et al.,Color image denoising via sparse 3D collaborative
ltering with grouping constraint in luminancechrominance space,Proc.
IEEE Int.Conf.Image Process.,ICIP 2007,San Antonio,TX,USA,2007.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 67
/184
Color Image Denoising with CBM3D.m
Flowchart of the colorimage denoising method.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 68
/184
Color Image Denoising Example
Top row contains noisefree Y,Cb,and Cr channels and the bottom row contains
corresponding noisy ones (σ = 22).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 69
/184
Color Image Denoising Example
On the left:noisy (σ =35) House and a fragment of it;on the right:the
corresponding denoised image (PSNR=31.58 dB) and fragment.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 70
/184
Image Deblurring:setting of the problem
Deblurring problem:reconstruct y
i
from fz
i
,x
i
g
n
i =1
,
z
i
= (y ~v)(x
i
) +σε
i
,
where ε
i
= ε(x
i
) i.i.d.N(0,1).
Standard approach:regularized inverse
ˆy = arg min
y
jjz (y ~v)jj
2
2
+µ pen(y),λ > 0.
An invariant µ is a principal limitation.
In what follows we use V(f ) = Ffvg and Z(f ) = Ffzg,Y(f ) = Ffyg.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 71
/184
BM3D based deblurring
Our approach is based on two inverse estimates:
Regularized inverse with small regularization parameter complemented by
3DBM ltering:in the frequency domain it is calculated as
ˆ
Y
RI
(f ) =
V
(f )
jV(f )j
2
+α
2
Z(f ),
Wiener inverse complemented by BM3D ltering:in the frequency domain
it is calculated as
ˆ
Y
RWI
(f ) =
V
(f )
jV(f )j
2
+n
1
n
2
σ
2
/j
ˆ
Y
RI
(f )j
2
Z(f ).
This ltering implements (imitates) e¤ects of the varying adaptive
regularization parameter µ.

K.Dabov,et al.,Image restoration by sparse 3D transformdomain
collaborative ltering,Proc.SPIE Electronic Imaging 08,no.681207,San
Jose,2008.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 72
/184
BM3D based deblurring
Our approach is based on two inverse estimates:
Regularized inverse with small regularization parameter complemented by
3DBM ltering:in the frequency domain it is calculated as
ˆ
Y
RI
(f ) =
V
(f )
jV(f )j
2
+α
2
Z(f ),
Wiener inverse complemented by BM3D ltering:in the frequency domain
it is calculated as
ˆ
Y
RWI
(f ) =
V
(f )
jV(f )j
2
+n
1
n
2
σ
2
/j
ˆ
Y
RI
(f )j
2
Z(f ).
This ltering implements (imitates) e¤ects of the varying adaptive
regularization parameter µ.

K.Dabov,et al.,Image restoration by sparse 3D transformdomain
collaborative ltering,Proc.SPIE Electronic Imaging 08,no.681207,San
Jose,2008.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 72
/184
Image Deblurring with BM3D (BM3DDEB.m)
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 73
/184
Color Filter Array Demosaicking
Bayer color lter array.
A.Danielyan,et al.,Proc.Int.Workshop on Local and NonLocal Approx.in
Image Process.,LNLA 2009,Tuusula,Finland,August 2009.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 74
/184
NoiseFree Demosaicking
High correlation between channels.It is likely that color channels is going to have
similar texture and edge locations.

L.Zhang and X.Wu,Color demosaicking via directional linear minimum mean
squareerror estimation,IEEE Trans.IP.,vol.14,no.12,2005.
Paliy,et al.Denoising and Interpolation of Noisy Bayer Data with Adaptive
CrossColor Filters,SPIEIS&T Electronic Imaging,Visual Communications and
Image Processing 2008,vol.6822,San Jose,2008.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 75
/184
NoiseFree CFA Demosaicking.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 76
/184
Noisy CFA Demosaicking
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 77
/184
How to deal with noise?
Denoising AFTER demosaicking.
No adequate noise model!Interpolation changes statistical model of the noise in
a complex and hardly computable form;
Filter should rely only on constraints reecting the general a priori knowledge
about the image structure.
JOINT demosaickingdenoising.(Paliy et al.,Zhang et.al)
Modify demosaicking method to be robust against noise;
Design is di¢ cult!Combining antagonistic procedures.Denoising s smoothing,
while demosaicking/interpolation s reconstruction of high frequency details.
Can we do Denoising BEFORE demosaicking?
Mosaic structure violates assumptions about local smoothness of the natural
images on which lters were relying;
Split(R,G
1
,G
2
,B) > Denoise > Combine  this leads to smoothing ne
details.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 78
/184
Modern CrossColor Denoising Approaches
What has changed in novel algorithms?
Local smoothing is not critical anymore.Filters exploiting nonlocal similarity of
small image patches.

Zhang,et all.,PCAbased Spatial Adaptive Denoising of CFA Images for
SingleSensor Digital Cameras,IEEE Trans.IP,vol.18,no.4,2009.
A.Danielyan,et al.,Crosscolor BM3D ltering of noisy raw data,Proc.Int.
Workshop on Local and NonLocal Approx.in Image Process.,LNLA 2009,
Tuusula,Finland,pp.125129,2009.
Foi,A.,Clipped noisy images:heteroskedastic modeling and practical
denoising,Signal Processing,vol.89,no.12,pp.26092629,Dec.2009.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 79
/184
BM3D with colorconstrained grouping
Block grouping in BM3D modeling.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 80
/184
Experiment 1.Denoising and demosaicking (GN).
Proposed
Zhang,2009
Proposed
Zhang,2009
σ
5/255
12/255
R
G
B
37.8
39.1
37.5
36.8
38.0
36.6
33.9
34.6
33.8
32.6
33.2
32.6
R
G
B
32.7
34.5
32.8
31.7
33.2
31.9
29.6
30.5
29.8
28.5
29.3
28.7
R
G
B
35.7
36.8
36.3
34.8
35.9
35.4
31.8
32.5
32.6
30.9
31.6
31.6
R
G
B
37.7
39.3
38.2
37.3
38.7
37.6
34.4
35.5
34.6
33.9
34.8
33.9
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 81
/184
Experiment 2.GaussianPoissonian Noise and clipping
Proposed
Zhang,2009
σ
2
(y) = ay +b
a = 0.004,b = 0.02
2
R
G
B
34.1
34.9
34.3
32.7
33.5
33.1
R
G
B
29.5
30.4
29.7
28.3
29.2
28.5
R
G
B
32.0
32.6
32.7
31.0
31.7
31.7
R
G
B
34.2
35.4
34.8
33.7
34.7
34.2
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 82
/184
Gaussian noise
From left to right:ground truth,proposed denoising + interpolation
(Zhang,2005),denoising (Zhang,2009) + interpolation (Zhang,2005),
Gaussian noise,σ = 12/255.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 83
/184
Signaldependent noise
From left to right:ground truth,proposed denoising + interpolation
(Zhang,2005),denoising (Zhang,2009) + interpolation (Zhang,2005),
signal dependent noise (a = 0.004,b = 0.02 ).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 84
/184
BM3D algorithm with adaptiveshape neighborhoods
The algorithm uses grouping of adaptiveshape neighborhoods whose
surrounding square supersets have been found similar by a blockmatching
procedure.
The data dened on these grouped neighborhoods is stacked together,
resulting in 3D structures which are generalized cylinders with
adaptiveshape cross sections.
These 3D groups are characterized by a high correlation along all the three
dimensions.
A 3D decorrelating transform is computed as a separable composition of the
ShapeAdaptive DCT (SADCT) and a 1D orthonormal transform.
K.Dabov,A.Foi,V.Katkovnik,and K.Egiazarian,A Nonlocal and ShapeAdaptive TransformDomain Collaborative
Filtering,Proc.2008 Int.Workshop on Local and NonLocal Approximation in Image Processing,LNLA 2008,
Lausanne,Switzerland,2008.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 85
/184
Illustration of applying 3D transform on a group of shapeadaptive
neighborhoods.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 86
/184
BM3D SADCT BM3DSADCT
27.93 27.51 27.95
BM3Ds good reconstruction of textures and regular image structures and
SADCTs good reconstruction of sharp edges of varying curvature and image
singularities.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 87
/184
BM3D algorithm with shapeadaptive PCA
A proper selection of the transform is crucial element for ensuring the
success of the transformbased methods.
This problem,known under di¤erent names as best basis,dictionary,or prior
selection,has been a subject of intensive study from the very beginning of
the development and application of estimation/approximation methods.In
particular,the use of bases adaptive to the data at hand is of special interest.
This latest version of BM3D algorithm incorporating a shapeadaptive
Principal Component Analysis (PCA) as part of the 3D transform.For a
3D group of adaptiveshape image patches,a shapeadaptive PCA basis is
obtained by eigenvalue decomposition of an empirical secondmoment
matrix computed from these patches.
Overall 3D transform is a separable composition of the PCA (applied on
each image patch) and a xed orthogonal 1D transform in the third
dimension.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 88
/184
Flowchart of BM3DSAPCA
K.Dabov,et al.,BM3D Image Denoising with ShapeAdaptive Principal Component Analysis,Proc.Workshop on Signal
Processing with Adaptive Sparse Structured Representations (SPARS09),SaintMalo,France,April 2009
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 89
/184
Principal Components of PCA in BM3DSAPCA
Illustration of the PCs (shown on the right side).The green overlay is used to
show the found similar neighborhoods used to form a 3D group.
The PCs are listed in decreasing magnitude of their corresponding eigenvalues.
The rst few PCs have the strongest similarity with the noisefree signal in the
neighborhood.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 90
/184
PCA calculations
.Each of the 2D neighborhoods in the group is vectorized with vectors a
s
of the length N
e
.The sample estimate of the covariance matrix N
e
N
e
is calculated as
ˆ
C = [a
1
,a
2
,...,a
N
g
] [a
1
,a
2
,...,a
N
g
]
T
,
where N
g
is a number of the patches in the group.
The PCA eigenvalue decomposition yields
U
T
ˆ
CU = S = diagfs
1
,...,s
N
e
g.
The largest eigenvalues such that s
j
> th σ
2
are selected and the
corresponding orthonormal columns of U are used as a set of orthonormal
bases.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 91
/184
Illustrations
The use of a datadriven adaptive transform for the collaborative ltering results
in a further improvement of the denoising performance,especially in preserving
image details and introducing very few artifacts SAPCA.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 92
/184
Numerical comparison
15
20
25
35
2.5
2
1.5
1
0.5
0
Noise standard deviation
Difference in PSNR [dB]
34.43
33.20
32.23
30.72
Lena testimage
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 93
/184
Numerical comparison
15
20
25
35
2.5
2
1.5
1
0.5
0
Noise standard deviation
Difference in PSNR [dB]
32.37
30.88
29.81
28.17
Cameraman testimage
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 94
/184
Part II:Some Applications
NonGaussian Image Processing;
Video Processing with BM3D;
Compressive Sensing;
Image Resizing;
Video SuperResolution;
BM3D Denoising and Root Sharpening
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 95
/184
NonGaussian image processing:observation model
Consider a signal dependent observation model in the form
z(x) = y (x) +σ(y (x)) ξ(x),x 2 X,
where the errors (noise) η(x) = z(x) y (x) are random independent and the
standard deviation of these observations is modeled by stdfη(x)g = σ(y (x)).
It is assumed that Efξ(x)g = 0,Efξ
2
(x)g = 1.
For Poissonian data:z(x) P(y(x)),Efz(x)g = y(x),varfz(x)g = y(x),
σ (y (x)) =
q
y(x).
If the noise is an additive mix of Gaussian and Poissonian random variables,then
varfz(x)g = σ
2
+γy(x).
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 96
/184
NonGaussian image processing:observation model
For imaging sensors (chargecoupled device (CCD) or complementary
metaloxidesemiconductor (CMOS) imager)
varfz(x)g = ay(x) +b,y(x) b/a,a,b > 0.
For clipped data
¯z(x) = maxf0,minf1,z(x)gg,¯y(x) = Ef¯z(x)g 6= y(x),
with standard deviation curves
It is important to compute the functions ¯y and
¯
σ given and y,and vice versa.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 97
/184
Variance stabilization:as a general approach
The Anscombe transform for Poissonian signals
˜y(x) = 2
r
1 +
3
8
y(x),y(x) =
8
3
(˜y
2
(x)/4 1),
gives var f˜y(x)g'1.˜y(x) is approximately Gaussian provided y(x) > 20.
In general a stabilizing transform ˜y = f (¯z) is derived as follows.Using the
Taylor series
˜y = f (¯z) = f (¯y +(¯z ¯y))'f (¯y) +f
0
(¯y)(¯z ¯y) )
std(˜y) = f
0
(¯y)std(¯y)
Assume that std(˜y) = c then
f (t) =
Z
t
0
c
std(¯y)
d ¯y.
Foi,A.,Clipped noisy images:heteroskedastic modeling and practical
denoising,Signal Processing,vol.89,no.12,pp.26092629,December 2009.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 98
/184
Variance stabilization denoising
Algorithm consists from the following four steps:
Variance stabilization transform;
Denoising clipped data using homoskedastic Gaussian noise algorithm;
Inverse of the variance stabilizing transform;
Declipping.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 99
/184
Examples
Clipped noisy data.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 100
/184
Examples
Figure:
Denoising clipped data.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 101
/184
Examples
Denoising and declipping.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 102
/184
Examples
Crosssections of observations and estimates.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 103
/184
PSNR (dB) values for the denoised D(¯z),denoised and declipped ˆy,and
rangeconstrained estimates.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 104
/184
Examples
noisy
denoisedclipped
denoiseddeclipped
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 105
/184
Examples
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 106
/184
Video Processing with BM3D (VBM3D.m)
It implements a video denoising method that enables highly e¤ective noise
attenuation at near realtime execution times.
It is based on enhanced sparse representation in local 3D transform domain.
As the noisy video is processed in blockwise manner,the sparsity enhancement is
achieved by grouping 2D fragments similar to the current one into a 3D data
array that we call group.
K.Dabov,et al.,Video denoising by sparse 3D transformdomain collaborative
ltering,Proc.15th European Signal Processing Conference,EUSIPCO 2007,
Poznan,Poland,September 2007.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 107
/184
Video Processing with BM3D (cont.)
Illustration of the proposed predictive search blockmatching.
Each pixel represents a block located at it;grey pixels denote blocks that are part
of the search neighborhood;red pixels denote blocks matched as similar.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 108
/184
Video Processing with BM3D (cont.)
The grouping is realized as a predictivesearch blockmatching,similar to
techniques used for motion estimation;
For each formed group,we apply collaborative ltering in order to take
advantage of the correlation between grouped blocks.We realize this
ltering by a 3D transformdomain shrinkage (hardthresholding and Wiener
ltering).The collaborative ltering produces estimates of all grouped
blocks.
Since these estimates are overlapping in general,we aggregate them by a
weighted average in order to form a nonredundant video estimate.
A signicant improvement of this approach is the use of 2step algorithm
where an intermediate estimate is produced by grouping and collaborative
hardthresholding and then used for improving the grouping and for applying
collaborative empirical Wiener ltering.
The experimental results demonstrate the stateoftheart denoising
performance and subjective visual quality.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 109
/184
Video Processing with BM3D (cont.)
For frame t,which contains the reference block,the search neighborhood is
(nonadaptive) xed around the reference block.
For frame t+1,the blocks matched as similar to in the frame t are used to
determine the centers of the small 3 3 neighborhoods,whose union forms
the overall adaptive neighborhood for this frame.
For frame t+2,the motion compensated predictive search uses the motion
vectors (shown as white arrows) between the matched blocks of the previous
two frames to determine the centers of the 3 3 neighborhoods.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 110
/184
Flowchart of the algorithm
Flowchart of the proposed BM3D video denoising method.
The operation enclosed by dashed lines are repeated for each reference block.
Grouping is illustrated by showing a reference block marked with Rand the
matched ones in a temporal window of 5 frames.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 111
/184
Videodemo,Bus
BUS!
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 112
/184
Videodemo,Flowers
Flowers!
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 113
/184
Videodemo,Man
Man!
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 114
/184
Compressive sensing (CS):conventional approach
The basic setting:an unknown signal of interest is observed (sensed)
through a limited number linear functionals;
These observations can be considered as an incomplete portion of the
spectrum of the signal with respect to a given linear transform;
It is assumed that the signal can be represented sparsely with respect to a
di¤erent relevant basis (e.g.,wavelets);
The algorithms rely on convex optimization with a penalty expressed by the
`
0
or`
1
norms which is exploited to enable the assumed sparsity;
It results in parametric modeling of the solution and in problems that are
then solved by mathematical programming algorithms.
Candes,E.,J.Romberg,and T.Tao,Robust uncertainty principles:exact signal
reconstruction from highly incomplete frequency information,IEEE Trans.Inf.
Theory,vol.52,no.2,2006.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 115
/184
Compressive sensing (CS):our approach
We replace the traditional parametric modeling used in CS by a
nonparametric one.
The nonparametric modeling is implemented by the use of spatially adaptive
lters.
The logic behind of this approach is as follows.
The regularization imposed by the`
0
or`
1
norms (or by more general
criteria) is essentially only as a tool for design of some nonlinear ltering.
Let us replace this implicit regularization by explicit ltering,exploiting
spatially adaptive lters sensitive to image features and details.
Then,CS signal reconstruction is realized by a recursive algorithm based on
spatially adaptive image denoising.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 116
/184
Examples:observations
Sample domain Ω for the FFT spectrum::22 radial lines,11 radial lines,
90 degrees limitedangle with 61 radial lines.
Vladimir Katkovnik (TUT)
Nonlocal imaging
November 7,2009 ICIP 2009,Cairo,Egypt 117
/184
Examples:reconstructions
Clockwise from topleft:backprojection estimates for 22 radial lines,11 radial
Comments 0
Log in to post a comment