# Edge Detection¶


This numerical tour explores local differential operators (grad, div, laplacian) and their use to perform edge detection.

In [2]:
addpath('toolbox_signal')


## Diffusion and Convolution¶

To obtain robust edge detection method, it is required to first remove the noise and small scale features in the image. This can be achieved using a linear blurring kernel.

Size of the image.

In [3]:
n = 256*2;


Load an image $f_0$ of $N=n \times n$ pixels.

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


Display it.

In [5]:
clf;
imageplot(f0);


Blurring is achieved using convolution: $$f \star h(x) = \sum_y f(y-x) h(x)$$ where we assume periodic boundary condition.

This can be computed in $O(N\log(N))$ operations using the FFT, since $$g = f \star h \qarrq \forall \om, \quad \hat g(\om) = \hat f(\om) \hat h(\om).$$

In [6]:
cconv = @(f,h)real(ifft2(fft2(f).*fft2(h)));


Define a Gaussian blurring kernel of width $\si$: $$h_\si(x) = \frac{1}{Z} e^{ -\frac{x_1^2+x_2^2}{2\si^2} }$$ where $Z$ ensure that $\hat h(0)=1$.

In [7]:
t = [0:n/2 -n/2+1:-1];
[X2,X1] = meshgrid(t,t);
normalize = @(h)h/sum(h(:));
h = @(sigma)normalize( exp( -(X1.^2+X2.^2)/(2*sigma^2) ) );


Define blurring operator.

In [8]:
blur = @(f,sigma)cconv(f,h(sigma));


Exercise 1

Test blurring with several blurring size $\si$.

In [9]:
exo1()

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


The simplest edge detectors only make use of the first order derivatives.

For continuous functions, the gradient reads $$\nabla f(x) = \pa{ \pd{f(x)}{x_1}, \pd{f(x)}{x_2} } \in \RR^2.$$

We discretize this differential operator using first order finite differences. $$(\nabla f)_i = ( f_{i_1,i_2}-f_{i_1-1,i_2}, f_{i_1,i_2}-f_{i_1,i_2-1} ) \in \RR^2.$$ Note that for simplity we use periodic boundary conditions.

Compute its gradient, using (here decentered) finite differences.

In [11]:
s = [n 1:n-1];
nabla = @(f)cat(3, f-f(s,:), f-f(:,s));


One thus has $\nabla : \RR^N \mapsto \RR^{N \times 2}.$

In [12]:
v = nabla(f0);


One can display each of its components.

In [13]:
clf;
imageplot(v(:,:,1), 'd/dx', 1,2,1);
imageplot(v(:,:,2), 'd/dy', 1,2,2);


A simple edge detector is simply obtained by obtained the gradient magnitude of a smoothed image.

A very simple edge detector is obtained by simply thresholding the gradient magnitude above some $t>0$. The set $\Ee$ of edges is then $$\Ee = \enscond{x}{ d_\si(x) \geq t }$$ where we have defined $$d_\si(x) = \norm{\nabla f_\si(x)}, \qwhereq f_\si = f_0 \star h_\si.$$

Compute $d_\si$ for $\si=1$.

In [14]:
sigma = 1;
d = sqrt( sum(nabla(  blur(f0,sigma)  ).^2,3) );


Display it.

In [15]:
clf;
imageplot(d);


Exercise 2

For $\si=1$, study the influence of the threshold value $t$.

In [16]:
exo2()

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


Exercise 3

Study the influence of $\si$.

In [18]:
exo3()

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


## Zero-crossing of the Laplacian¶

Defining a Laplacian requires to define a divergence operator. The divergence operator maps vector field to images. For continuous vector fields $v(x) \in \RR^2$, it is defined as $$\text{div}(v)(x) = \pd{v_1(x)}{x_1} + \pd{v_2(x)}{x_2} \in \RR.$$ It is minus the adjoint of the gadient, i.e. $\text{div} = - \nabla^*$.

It is discretized, for $v=(v^1,v^2)$ as $$\text{div}(v)_i = v^1_{i_1+1,i_2} + v^2_{i_1,i_2+1}.$$

In [20]:
t = [2:n 1];
div = @(v)v(t,:,1)-v(:,:,1) + v(:,t,2)-v(:,:,2);


The Laplacian operatore is defined as $\Delta=\text{div} \circ \nabla = -\nabla^* \circ \nabla$. It is thus a negative symmetric operator.

In [21]:
delta = @(f)div(nabla(f));


Display $\Delta f_0$.

In [22]:
clf;
imageplot(delta(f0));


Check that the relation $\norm{\nabla f} = - \dotp{\Delta f}{f}.$

In [23]:
dotp = @(a,b)sum(a(:).*b(:));
fprintf('Should be 0: %.3i\n', dotp(nabla(f0), nabla(f0)) + dotp(delta(f0),f0) );

Should be 0: -8.879e-11


The zero crossing of the Laplacian is a well known edge detector. This requires first blurring the image (which is equivalent to blurring the laplacian). The set $\Ee$ of edges is defined as: $$\Ee = \enscond{x}{ \Delta f_\si(x) = 0 } \qwhereq f_\si = f_0 \star h_\si .$$

It was proposed by Marr and Hildreth:

Marr, D. and Hildreth, E., Theory of edge detection, In Proc. of the Royal Society London B, 207:187-217, 1980.

Display the zero crossing.

In [24]:
sigma = 4;
clf;
plot_levelset( delta(blur(f0,sigma)) ,0,f0);


Exercise 4

Study the influence of $\si$.

In [25]:
exo4()