Heat Diffusion
This numerical tour explores local differential operators (grad, div, laplacian) and their use to perform heat diffusion.
Contents
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/');
Gradient, Laplacian
Local differential operators like gradient, divergence and laplacian are the building blocks for variational image processing using PDEs.
Load an image \(f_0\) of \(N=n \times n\) pixels.
n = 256;
f0 = load_image('boat',n);
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 finite differences.
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}. \)
v = nabla(f0);
One can display each of its components.
clf; imageplot(v(:,:,1), 'd/dx', 1,2,1); imageplot(v(:,:,2), 'd/dy', 1,2,2);
One can also display it using a color image.
clf; imageplot(v);
One can display its magnitude \(\norm{\nabla f(x)}\), which is large near edges.
clf; imageplot( sqrt( sum3(v.^2,3) ) );
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}. \]
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.
delta = @(f)div(nabla(f));
Display \(\Delta f_0\).
clf; imageplot(delta(f0));
Check that the relation \( \norm{\nabla f} = - \dotp{\Delta f}{f}. \)
dotp = @(a,b)sum(a(:).*b(:));
fprintf('Should be 0: %.3i\n', dotp(nabla(f0), nabla(f0)) + dotp(delta(f0),f0) );
Should be 0: -2.012e-07
Heat Diffusion
The heat equation is an elliptic linear PDE \[ \pdd{f}{t} = \Delta f \qwhereq f_{t=0}=f_0. \] It corresponds to the gradient descent flow of the energy \(E(f)=\norm{\nabla f}\).
It can be discretized in space using a finite difference Laplacian \(\Delta\) and in time using an explicit Euler time stepping, which leads to the following iterative scheme \[ f^{(0)} = f_0 \qandq f^{(\ell+1)} = f^{(\ell)} + \tau \Delta f^{(\ell)}\] where \(\tau>0\) is a small time step, and \(f^{(\ell)}\) approximate the solution at time \(t=\ell\tau\).
Final smoothing time.
T = 10;
Step size \(\tau\) for the discretization. It should be smaller than 1/4 for the discretization to be stable. A smaller time step means more iterations but a higher precision.
tau = .05;
Number of iteration.
niter = ceil(T/tau);
Initialize the diffusion process at time \(t=0\).
f = f0;
Perform the descent step.
f = f + tau*delta(f);
Exercice 1: (the solution is exo1.m) Compute a Heat diffusion flow by performing a gradient descent of the square norme of the gradient. Record in a vector E the decay of the Sobolev energy \(E(f^{(\ell)})\).
exo1;
Display the decay of the Sobolev energy through the iterations. It should converges toward 0.
clf; plot(log10(E/E(1))); title('log(E)'); axis('tight');
