结构张量扩散

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.
• Computing the Structure Tensor
• Tensor Driven Anisotropic Diffusion
• Curve and Anisotropy Detection

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:

f(x,t)t=div(T(x)f(x,t)).

The PDE is discretized in time using a explicit time stepping. The time step. should be small enough for the diffusion to be stable.

dt = .2;

First initialize the image to diffuse at time t=0 .

f1 = f;

Perform the update of the solution.

f1 = f1 + dt * div( TensorMult(T, grad(f1) ) );

Exercice 1: (the solution is exo1.m) Perform the full diffusion up to a large enough time.

exo1;

Another useful application of diffusion PDEs is to perform the visualisation of vector and tensor fields.

Exercice 2: (the solution is exo2.m) Do the same thing, but start from a noise image rand(n) at time t=0 . At each iteration, maintain the gray value histograms by performing an histogram equalization using perform_hist_eq.

exo2;

Exercice 3: (the solution is exo3.m) Perform the same diffusion, but for different anisotropies.

exo3;

Exercice 4: (the solution is exo4.m) Design a simple vector field, and use this method to vizualize it using a filtered noise.

exo4;

Curve and Anisotropy Detection

The structure tensor can be used to enhance the anisotropy of point clouds.

n = 160;

Compute a noisy image representing partial curves.

f = zeros(n);
r = [.3 .6 .4]*n;
c = [[.5;.5] [0;0] [.8;.8]]*n;
for i=1:length(r)
    q = round(r(i)*5);
    t = rand(q,1)*2*pi;
    x = round( c(1,i) + r(i)*cos(t) ); y = round( c(2,i) + r(i)*sin(t) );
    I = find(x>0 & x<=n & y>0 & y<=n);
    x = x(I); y = y(I);
    f(x+(y-1)*n) = 1;
end

Add noise.

f = f + rand(n)*.4;

Display.

clf;
imageplot(f);

Exercice 5: (the solution is exo5.m) Compute the structure tensor.

exo5;

Display energy and anisotropy.

[e1,e2,lambda1,lambda2] = perform_tensor_decomp(T);
E = sqrt(lambda1+lambda2);
A = sqrt(lambda1./lambda2);
clf;
imageplot({E A}, {'Energy', 'Anisotropy'});

Display only the tensor in the region with high energy.

I = find( E<=mean(E(:))*1.1 );
T1 = T;
T1([I; I+n^2; I+2*n^2; I+3*n^2]) = 0;
options.sub = 4;
clf;
plot_tensor_field(T1, [], options);