Linear Diffusion
Linear Diffusion Flows
This tours studies linear diffusion PDEs, a.k.a. the heat equation. A good reference for diffusion flows in image processing is [Weickert].
Contents
Installing toolboxes and setting up the path. Heat Diffusion Explicit Solution using Convolution Bibliography Installing toolboxes and setting up the path.
You need to download the following files: signal toolbox and general toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general in your directory.
For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.
Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.
Execute this line only if you are using Matlab.
getd = @(p)path(p,path); % scilab users must *not* execute this Then you can add the toolboxes to the path.
getd('toolbox_signal/'); getd('toolbox_general/'); Heat Diffusion
The heat equation reads for a function and where (the solution at initial time ) is given.
The Laplacian operator reads Δ The flow is discretized in space by considering a discrete image of pixels.
n = 256; Load an image , that will be used to initialize the flow at time .
name = 'hibiscus'; f0 = load_image(name,n); f0 = rescale( sum(f0,3) ); Display it.
clf; imageplot(f0);
The flow is discretized in time using an explicit time-stepping Δ We use finite difference Laplacian Δ where we assume periodic boundary conditions, and where is the spacial step size.
h = 1/n; delta = @(f)1/h^2 * div(grad(f)); The step size should satisfy for the discretized flow to be stable.
The discrete solution converges to the continuous solution at time if both and under the condition .
Select a small enough step size.
tau = .5 * h^2/4; Final time.
T = 1e-3; Number of iterations.
niter = ceil(T/tau); Initialize the diffusion at time .
f = f0; One step of discrete diffusion.
f = f + tau * delta(f); Exercice 1: (the solution is exo1.m) Compute the solution to the heat equation.
exo1;
Explicit Solution using Convolution
The solution to the heat equation can be computed using a convolution where denotes the convolution of continuous functions and is a Gaussian kernel of width
One can thus approximate the solution using a discrete convolution. Convolutions can be computed in operations using the FFT, since
cconv = @(f,h)real(ifft2(fft2(f).*fft2(h))); Define a discrete Gaussian blurring kernel of width .
t = [0:n/2 -n/2+1:-1]; [X2,X1] = meshgrid(t,t); normalize = @(h)h/sum(h(:)); h = @(t)normalize( exp( -(X1.^2+X2.^2)/(4*t) ) ); Define blurring operator.
heat = @(f, t)cconv(f,h(t)); Example of blurring.
clf; imageplot(heat(f0,2));
Exercice 2: (the solution is exo2.m) Display the heat convolution for increasing values of .
exo2;
Bibliography
[Weickert98] Joachim Weickert, Anisotropic Diffusion in Image Processing, ECMI Series, Teubner-Verlag, Stuttgart, Germany, 1998.
Copyright (c) 2010 Gabriel Peyre