Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This numerical tour explores the structure tensor to represent the geometry of images and textures. It applies it to perform anisotropic image diffusion. A good reference for diffusion flows in image processing is Weickert98.
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/pde_3_diffusion_tensor')
We define here a few features (convolution, gradient, etc.) that will be used in the sequel.
Size of the image of $N=n \times n$ pixels.
n = 256;
Load an image $f$.
name = 'hibiscus';
f = load_image(name,n);
f = rescale( sum(f,3) );
Display it.
clf;
imageplot(f);
We define circular convolution $$ (f \star h)_i = \sum_j f_j h_{i-j}. $$ Note that here, $f$ can be multi-channel, in which case each channel is convolved with $h$. This will be useful to blur tensor fields.
cconv = @(f,h)real(ifft2(fft2(f).*repmat(fft2(h),[1 1 size(f,3)])));
Define a Gaussian blurring kernel of width $\si$: $$ h_\si(x) = \frac{1}{Z} e^{ -\frac{x_1^2+x_2^2}{2\si^2} }$$ where $Z$ ensures that $\hat h_\si(0)=1$.
t = [0:n/2 -n/2+1:-1];
[X2,X1] = meshgrid(t,t);
normalize = @(h)h/sum(h(:));
h = @(sigma)normalize( exp( -(X1.^2+X2.^2)/(2*sigma^2) ) );
Define the convolution with $h_\si$.
blur = @(f,sigma)cconv(f,h(sigma));
We use in the following a centered finite difference approximation of $\nabla f$, which is a vector field in $\RR^{n \times n \times 2}$.
options.order = 2;
nabla = @(f)grad(f,options);
We define the tensor product associated to a vector $u = (u_1,u_2), v=(u_1,u_2) \in \RR^{2}$ as the symetric matrix $$ u \otimes v = u v^* = \begin{pmatrix} u_1 v_1 & v_1 u_2 \\ u_1 v_2 & u_2 v_2 \end{pmatrix}$ \in \RR^{2 \times 2}. $$ It is extended to vector fields $ (u(x))_x \in \RR^{N \times 2} $ as $$ (u \otimes v)(x) = u(x) \otimes v(x) $$
A tensor field $T$ is a collection of symmetric positive definite matrices $T(x) \in \RR^{2 \times 2}$.
A simple way to build a tensor field is by auto-tensorization of a vector field $u(x)$, i.e. $T = u \otimes u$.
Define a shortcut for $u \otimes u$ (we make use of symmetry to only store 3 components).
tensorize = @(u)cat(3, u(:,:,1).^2, u(:,:,2).^2, u(:,:,1).*u(:,:,2));
Rotate a tensor field by $\pi/2$ (for display only).
rotate = @(T)cat(3, T(:,:,2), T(:,:,1), -T(:,:,3));
The structure tensor is a field of symetric positive matrices that encodes the local orientation and anisotropy of an image.
It was initially introduced for corner detection HarSteph88 Forstner86 and oriented texture analysis KassWit85.
Given an image $f$, its structure tensor with scale $ \sigma>0 $ is defined as $$ T_\si = h_\si \star T_0 \qwhereq T_0 = \nabla f \otimes \nabla f. $$ For each location $x$, $T_\si(x)$ is thus a positive definite matrix.
T = @(f,sigma)blur( tensorize( nabla(f) ), sigma);
The matrix $T_\si(x)$ can be understood as the local covariance matrix of the set of gradient vector around $x$.
Another way to get some insight about this tensor field is to consider a localized version $f_x$ of the image around point $x$, defined by $f_x(y) = h_\si(x-y)^{1/2} f(y)$, which is close to zero when $y$ is far away from $x$. One has the following Taylor expansion of the $L^2$ norm between two close enough localizations: $$ \norm{f_x - f_{x+\de}}^2 = \de^* T_\si(x) \de + O(\norm{\de}^3). $$
To better understand the behavior of $T_\si$ as a function of $\si$, one can computes its Taylor expansion for small $\si$ $$ T_\si(x) = T_0(x) + \si^2 Hf(x)^2 + O(\si^3), $$ where $Hf(x) \in \RR^{2 \times 2}$ is the Hessian matrix of $f$ at point $x$. This shows that when $\si$ increases, the intial rank-1 tensor $T_0(x)$ becomes full rank because it integrates energy from $Hf(x)^2$.
A convenient way to display a tensor field such as $T_\si$ is to draw an ellispe $\Ee_x$ at each pixel $x$ as the (scaled and translated) unit ball of the tensor $$ \Ee_x = \enscond{\de \in \RR^2}{ \de^* T_\si(x) \de \leq 1 }. $$ This allows one to visualize the anisotropy and orientation encoded in the tensor field.
Display $T_\si$ for $\si=0.1$ (the tensors are almost rank-1):
options.sub = 8;
clf; sigma = .1;
plot_tensor_field(rotate(T(f,sigma)), f, options);
title(['\sigma=' num2str(sigma)]);
For $\si=4$:
clf; sigma = 4;
plot_tensor_field(rotate(T(f,sigma)), f, options);
title(['\sigma=' num2str(sigma)]);