'''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
