Color Image Denoising with Median Filtering¶

Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This numerical tour explores denoising of color images using a local multi-dimensional median. This is the sequel to the numerical tour <../tv_median/ Outliers and Median Denoiser>.

In [2]:

addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/multidim_7_median')

Multidimensional Median¶

|\sum_i abs(x(i)-y)|

This should be contrasted with the mean that minimizes the sum of squares.

|\sum_i (x(i)-y)^2|

This allows one to define a mutidimensional median for set of points |x(i)| in dimension |d| by replacing |abs| by the d-dimensional norm.

We define a Gaussian point cloud in 2D.

In [3]:

d = 2; % dimension
n = 1000; % number of points
X = randn(d,n);

We modify some points as positive outliers (to shift the mean).

In [4]:

p = 100; % number of outliers
sel = randperm(n); sel = sel(1:p); % index of outliers
X(:,sel) = rand(d,p)*50;

We can compute the mean point.

In [5]:

m = mean(X,2);

To compute the median in 2D, one needs to minimize the sum of norms. This is not as straightforward as the sum of squares, since there is no close form solution. One needs to use an iterative algorithm, for instance the re-weighted least squares, that computes weighted means.

number of iterations of the method

In [6]:

niter = 30;

initialize the median using the mean

In [7]:

med = m;
energy = [];
for i=1:niter
    % comute the distance from med to the points
    dist = sqrt( sum( (X-repmat(med,[1 n])).^2 ) );
    % compute the weight, take care of not dividing by 0
    weight = 1./max( dist, 1e-10 ); weight = weight/sum(weight);
    % compute the weighted mean
    med = sum( repmat(weight,[d 1]).*X, 2 );
    energy(end+1) = sum( dist );
end

We can display the decay of the L1 energy through the iterations.

In [8]:

clf;
plot(energy, '.-'); axis('tight')
set_label('Iteration', 'L1 energy');

We can display the points, the mean and the median.

In [9]:

clf;
hold('on');
plot(X(1,:), X(2,:), '.');
plot(m(1,:), m(2,:), 'k*');
plot(med(1,:), med(2,:), 'ro');
axis('tight');

Color Image Denoising using 1D Median¶

A median filter can be used to denoise a color image, by applying it to each channel of the image.

We load a color image, which is an array of size |[n,n,3]|.

In [10]:

name = 'flowers';
options.nbdims = 3;
n = 256;
M0 = load_image(name, n, options);
M0 = rescale(M0);

We create a colored impulse noise by taking two Gaussians of different standard deviations.

percent of strong Gaussian

In [11]:

rho = .4;

mask of pixel corrupted by strong gaussian

In [12]:

mask = repmat(rand(n,n)<rho, [1 1 3]);

deviation of the two Gaussian

In [13]:

sigma1 = .03; sigma2 = 1;

noise with two different Gaussians

In [14]:

noise = sigma1*randn(n,n,3).*(1-mask) + sigma2*rand(n,n,3).*mask;

Add the noise to the image.

In [15]:

M = M0+noise;
pnoisy = snr(M0,M);

Display the clean and noisy images.

In [16]:

clf;
imageplot(M0, 'Clean image', 1,2,1);
imageplot(clamp(M), strcat(['Noisy, SNR=' num2str(pnoisy)]), 1,2,2 );

In the following, we use a fixed window width.

In [17]:

k = 4;
w = 2*k+1;

Exercise 1

In [18]:

exo1()