Anisotropic Diffusion with Structure Tensor

This numerical tour explores the structure tensor to represent the geometry of images and textures. It applies it to perform anisotropic image 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/');

Computing the Structure Tensor

The structure tensor is a field of symetric positive matrices that encodes the local orientation and anisotropy of an image.

It proposed for instance in the paper:

Michael Kass, Andrew P. Witkin, Analyzing Oriented Patterns. IJCAI, 1985: 944-952

Given an image f,
its structure tensor with scale σ>0
is defined as

T(x)=Gσ⋆( g(x)g(x)∗),whereg(x)=f(x)

where f(x)R2
is the gradient at x
of the image,

Gσ(x)=e−x22σ2

is the Gaussian kernel, ⋆
is the convolution operator, and v
is the rotation by π/2
of the vector vR2.

T(x)
is thus a positive definite matrix.

Load an image f.

n = 256;
name = 'fingerprint';
f = rescale(load_image(name,n));

Display it.

clf;
imageplot(f);

Up-sample the image to avoid aliazing before computing the structure tensor. This is important for textures containing high frequencies.

[X,Y] = meshgrid(1:n,1:n);
[XI,YI] = meshgrid(1:1/2:n+1/2,1:1/2:n+1/2); XI(:,end) = n; YI(end,:) = n;
f1 = interp2(X,Y,f,XI,YI);

Compute the gradient field f1(x).
We use a finite difference approximation.

G = grad(f1);

Compute the rank-1 tensor field associated to the gradient rotated by
π/2.

T0 = zeros(n*2,n*2,2,2);
T0(:,:,1,1) = G(:,:,2).^2;
T0(:,:,2,2) = G(:,:,1).^2;
T0(:,:,1,2) = -G(:,:,1).*G(:,:,2);
T0(:,:,2,1) = -G(:,:,1).*G(:,:,2);

Note: because of the square in the definition of the tensor, this causes its entry to be aliazed by a factor 2. In region of high frequency texture, the orientation would be aliazed without the up-sampling.

Blur the tensor field using a scale σ2.

sigma = 20;
T = perform_blurring(T0, sigma);

Sub-sample it.

T = T(1:2:end,1:2:end,:,:);

Display the tensor field.

options.sub = 8;
clf;
plot_tensor_field(T, [], options);

Since the tensor field a symmetric matrix field, it is possible to diagonalize it

T(x)=λ1(x)e1(x)e1(x)∗+λ2(x)e2(x)e2(x)∗,

where (e1,e2)
are the orthogonal eigenvector fields, 0λ1λ2
are the eigenvalues.

Diagonalize the tensor field.

[e1,e2,lambda1,lambda2] = perform_tensor_decomp(T);

Compute the energy and anisotropy

E(x)=λ1(x)+λ2(x)−−−−−−−−−−−andA(x)=λ1(x)λ2(x)−−−−−.

E = sqrt(lambda1+lambda2);
A = sqrt(lambda1./lambda2);

Display it.

clf;
imageplot({E A}, {'Energy', 'Anisotropy'});

Tensor Driven Anisotropic Diffusion

A tensor can be used as anisotropic metric to drive a diffusion PDE flow.

The best reference for such a flow is

J. Weickert, Anisotropic Diffusion in Image Processing. Teubner, Stuttgart, 1998.

We first derive the diffusion tensor field from the original tensor field. The idea is to remap non-linearly the eigenvalue to achieve various effects.

Here we make something very simple, we simply fix the anisotropy so that the diffusion in the main direction is proportional to the diffusion in the orthogonal direction.

Anisotropy.

rho = .1;

Tensor field.

T = perform_tensor_decomp(e1,e2,ones(n),ones(n)*rho);

Shortcut for the multiplication of tensor by vector field.

TensorMult = @(T,u)cat(3,  T(:,:,1,1).*u(:,:,1) + T(:,:,1,2).*u(:,:,2), ...
                            T(:,:,2,1).*u(:,:,1) + T(:,:,2,2).*u(:,:,2) );

We consider an anisotropic diffusion: