# Bilateral Filtering¶


This tour explores edge preserving filtering and in particular the bilateral filter, with applications to denoising and detail enhancement.

In [67]:
warning off
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/denoisingadv_8_bilateral')
warning on


## Gaussian Linear Filtering¶

The most basic filtering operation is the Gaussian filtering, that tends to blur edges.

Image size.

In [3]:
n = 256*2;


Load an image.

In [4]:
name = 'hibiscus';
f0 = load_image(name, n);
f0 = rescale(crop( sum(f0,3) ,n));


A Gaussian filter of variance $\si$ reads $$\Gg_\si(f)(x) = \frac{1}{Z} \sum_y G_\si(x-y) f(y) \qwhereq Z = \sum_y G_\si(y),$$ and where the Gaussian kernel is defined as: $$G_\si(x) = e^{-\frac{\norm{x}^2}{2\si^2}}.$$

A convolution can be computed either over the spacial domain in $O(N\si_x^2)$ operations or over the Fourier domain in $O(N \log(N))$ operations. Depending on the value of $\si_x$, one should prefer either Fourier or spacial domain. For simplicity, we consider here the Fourier domain (and hence periodic boundary conditions).

Define a Gaussian function, centered at the top left corner (because it corresponds to the 0 point for the FFT).

In [5]:
x = [0:n/2-1, -n/2:-1];
[Y,X] = meshgrid(x,x);
GaussianFilt = @(s)exp( (-X.^2-Y.^2)/(2*s^2) );


Display a recentered example of Gaussian.

In [6]:
clf;
imageplot(fftshift(GaussianFilt(40)));


Define a shortcut to perform linear Gaussian filtering over the Fourier domain. This function is able to process in parallel a 3D block |F| by filtering each |F(:,:,i)|.

In [7]:
Filter = @(F,s)real( ifft2( fft2(F).*repmat( fft2(GaussianFilt(s)), [1 1 size(F,3)] ) ) );


Example of filtering $\Gg_\si(f_0)$.

In [8]:
clf;
imageplot( Filter(f0,5) );


Exercise 1

Display a filtering $\Gg_\si(f_0)$ with increasing size $\si$.

In [9]:
exo1()

In [10]:
%% Insert your code here.


## Bilateral Filter¶

The bilateral filter is a spacially varying filter that better preserves edges than the Gaussian filter.

It was first introduced in:

Carlo Tomasi and Roberto Manduchi, Bilateral Filtering for Gray and Color Images, Proceedings of the ICCV 1998

A very good set of ressources concerning the bilateral filter can be found on <http://people.csail.mit.edu/sparis/bf/ Sylvain Paris web page>. It includes a fast implementation, research and review papers.

Given a spacial width $\si_x$ and a value width $\si_v$, the filter opterates as: $$\Bb_{\si_x,\si_v}(f)(x) = \frac{1}{Z_x} \sum_y G_{\si_x}( x-y ) G_{\si_v}(f(x)-f(y)) f(y)$$ where the normalizing constant is $$Z_x = \sum_y G_{\si_x}( x-y ) G_{\si_v}(f(x)-f(y)).$$

At a given pixel $x$, it corresponds to an averaring with the data-dependant kernel $G_{\si_x}( x-y ) G_{\si_v}(f(x)-f(y))$.

## Bilateral Filter by Stacking¶

Implementing the bilateral filter directly over the spacial domain requires $O( N \si_x^2 )$ operations where $N$ is the number of pixels.

A fast approximate implementation is proposed in:

Fast Bilateral Filtering for the Display of High-Dynamic-Range Images, Fredo Durand and Julie Dorsey, SIGGRAPH 2002.

It exploits the fact that for all pixels $x$ with the same value $v=f(x)$, the bilateral filter can be written as a ratio of convolution $$\Bb_{\si_x,\si_v}(f)(x) = F_v(x) = \frac{ [G_{\si_x} \star ( f \cdot W_v )](x) }{ [G_{\si_x} \star W_v](x) }$$ where the weight map reads $$W_v(x) = G_{\si_v}(v-f(x)).$$

Instead of computing all possible weight maps $W_v$ for all possible pixel values $v$, one considers a subset $\{v_i\}_{i=1}^p$ of $p$ values and computes the weights $\{ W_{v_i} \}_i$.

Using these convolutions, one thus optains the maps $\{ F_{v_i} \}_i$, that are combined to obtain an approximation of $\Bb_{\si_x,\si_v}(f)$.

The computation time of the method is $O(p N\log(N))$ over the Fourier domain (the method implemented in this tour) and $O(p N \si_x^2)$ over the spacial domain.

Value of $\sigma_x$:

In [11]:
sx = 5;


Value of $\sigma_v$:

In [12]:
sv = .2;


Number $p$ of stacks. The complexity of the algorithm is proportional to the number of stacks.

In [13]:
p = 10;


Function to build the weight stack $\{ W_{v_i} \}_i$ for $v_i$ uniformly distributed in $[0,1]$.

In [14]:
Gaussian = @(x,sigma)exp( -x.^2 / (2*sigma^2) );
WeightMap = @(f0,sv)Gaussian( repmat(f0, [1 1 p]) - repmat( reshape((0:p-1)/(p-1), [1 1 p]) , [n n 1]), sv );


Compute the weight stack map. Each |W(:,:,i)| is the map $W_{v_i}$ associated to the pixel value $v_i$.

In [15]:
W = WeightMap(f0,sv);


Exercise 2

Display several weights $W_{v_i}$.

In [16]:
exo2()

In [17]:
%% Insert your code here.


Shortcut to compute the bilateral stack $\{ F_{v_i} \}_i$.

In [18]:
bileteral_stack_tmp = @(f0,sx,W)Filter(W.*repmat(f0, [1 1 p]), sx) ./ Filter(W, sx);
bileteral_stack = @(f0,sx,sv)bileteral_stack_tmp(f0,sx,WeightMap(f0,sv));


Compute the bilateral stack $\{ F_{v_i} \}_i$.

In [19]:
F = bileteral_stack(f0,sx,sv);


Exercise 3

Display several stacks $F_{v_i}$.

In [20]:
exo3()

In [21]:
%% Insert your code here.


Destacking corresponds to selecting a layer $I(x) \in \{1,\ldots,p\}$ at each pixel. $$f_I(x) = F_{ v_{I(x)} }(x).$$

Shortcut for de-stacking using a set of indexes.

In [22]:
[y,x] = meshgrid(1:n,1:n);
indexing = @(F,I)F(I);
destacking = @(F,I)indexing(F,x + (y-1)*n + (I-1)*n^2);


The simplest reconstruction method performs the destacking using the nearest neighbor value: $$I(x) = \uargmin{ 1 \leq i \leq p } \abs{f(x)-v_i}.$$

Shortcut for performing de-stacking by nearest neighbor interpolation.

In [23]:
bilateral_nn = @(f0,sx,sv)destacking( bileteral_stack(f0,sx,sv), round( f0*(p-1) ) + 1 );


Display.

In [24]:
fNN = bilateral_nn(f0,sx,sv);
clf;
imageplot( fNN );


A better reconstruction is obtained by using a first order linear interpolation to perform the destacking. $$\Bb_{\si_x,\si_v}(f)(x) \approx (1-\la(x))f_{I(x)}(x) + \la(x) f_{I(x)+1}(x)$$ where $I(x)$ and $\la(x)$ are defined as $$f(x) = (1-\la(x)) v_{I(x)} + \la(x) v_{I(x)+1} \qwhereq \la(x) \in [0,1).$$

Shortcut for the linear interpolation reconstruction.

In [25]:
frac = @(x)x-floor(x);
lininterp = @(f1,f2,Fr)f1.*(1-Fr) + f2.*Fr;
bilateral_lin1 = @(F,f0)lininterp( destacking(F, clamp(floor(f0*(p-1)) + 1,1,p) ), ...
destacking(F, clamp(ceil(f0*(p-1)) + 1,1,p) ), ...
frac(f0*(p-1)) );
bilateral_lin = @(f0,sx,sv)bilateral_lin1(bileteral_stack(f0,sx,sv), f0);


Exercise 4

Compare nearest-neighbor and linear destacking.

In [26]:
exo4()

In [27]:
%% Insert your code here.


Exercise 5

Study the influence of $\sigma_x$ on the filter, for a fixed $\sigma_v=0.2$.

In [28]:
exo5()

In [29]:
%% Insert your code here.


Exercise 6

Study the influence of $\sigma_v$ on the filter, for a fixed $\sigma_x=8$.

In [30]:
exo6()

In [31]:
%% Insert your code here.


Note that this stack implementation of the bilateral filter can be quite inacurate, in particular if $p$ is small. A more precise fast implementation is described in:

A Fast Approximation of the Bilateral Filter using a Signal Processing Approach, Sylvain Paris and Fr do Durand, International Journal of Computer Vision (IJCV'09)

## Denoising using the Bilateral Filter¶

The first basic application of the bilateral filter is for denoising. It performs a filtering that respects edges.

Noise level.

In [32]:
mu = .05;


Create a noisy image $f = f_0 + w$ where $w$ is a Gaussian white noise of variance $\mu^2$.

In [33]:
f = f0 + randn(n,n)*mu;


Display the noisy image.

In [34]:
clf;
imageplot(clamp(f));