Structure tensor

The 2D structure tensor

Continuous version

For a function ${\displaystyle I}$ of two variables p = (x, y), the structure tensor is the 2×2 matrix

$${\displaystyle S_w(p)=\left[\begin{array}{cc}\int w(r)(I_x(p-r){)}^{2}dr& \int w(r)I_x(p-r)I_y(p-r)dr\\ \int w(r)I_x(p-r)I_y(p-r)dr& \int w(r)(I_y(p-r){)}^{2}dr\end{array}\right]}$$ where ${\displaystyle I_x}$ and ${\displaystyle I_y}$ are the partial derivatives of ${\displaystyle I}$ with respect to x and y; the integrals range over the plane ${\displaystyle {\mathbb{R}}^2}$; and w is some fixed "window function" (such as a Gaussian blur), a distribution on two variables. Note that the matrix ${\displaystyle S_w}$ is itself a function of p = (x, y).

The formula above can be written also as $S_w(p)=\int w(r)S_0(p-r)dr$, where ${\displaystyle S_0}$ is the matrix-valued function defined by $${\displaystyle S_0(p)=\left[\begin{array}{cc}(I_x(p){)}^2& I_x(p)I_y(p)\\ I_x(p)I_y(p)& (I_y(p){)}^2\end{array}\right]}$$

If the gradient ${\displaystyle \mathrm{\nabla }I=(I_x,I_y{)}^{\text{T}}}$ of ${\displaystyle I}$ is viewed as a 2×1 (single-column) matrix, where ${\displaystyle (\cdot {)}^{\text{T}}}$ denotes transpose operation, turning a row vector to a column vector, the matrix ${\displaystyle S_0}$ can be written as the matrix product ${\displaystyle (\mathrm{\nabla }I)(\mathrm{\nabla }I{)}^{\text{T}}}$ or tensor or outer product ${\displaystyle \mathrm{\nabla }I\otimes \mathrm{\nabla }I}$. Note however that the structure tensor ${\displaystyle S_w(p)}$ cannot be factored in this way in general except if ${\displaystyle w}$ is a Dirac delta function.

Discrete version

In image processing and other similar applications, the function ${\displaystyle I}$ is usually given as a discrete array of samples ${\displaystyle I[p]}$, where p is a pair of integer indices. The 2D structure tensor at a given pixel is usually taken to be the discrete sum

$${\displaystyle S_w[p]=\left[\begin{array}{cc}\sum _{r}w[r](I_x[p-r]{)}^2& \sum _{r}w[r]I_x[p-r]I_y[p-r]\\ \sum _{r}w[r]I_x[p-r]I_y[p-r]& \sum _{r}w[r](I_y[p-r]{)}^2\end{array}\right]}$$

Here the summation index r ranges over a finite set of index pairs (the "window", typically ${\displaystyle \{-m\dots +m\}\times \{-m\dots +m\}}$ for some m), and w[r] is a fixed "window weight" that depends on r, such that the sum of all weights is 1. The values ${\displaystyle I_x[p],I_y[p]}$ are the partial derivatives sampled at pixel p; which, for instance, may be estimated from by ${\displaystyle I}$ by finite difference formulas.

The formula of the structure tensor can be written also as $S_w[p]=\sum _{r}w[r]S_0[p-r]$, where ${\displaystyle S_0}$ is the matrix-valued array such that $${\displaystyle S_0[p]=\left[\begin{array}{cc}(I_x[p]{)}^2& I_x[p]I_y[p]\\ I_x[p]I_y[p]& (I_y[p]{)}^2\end{array}\right]}$$

Interpretation

The importance of the 2D structure tensor ${\displaystyle S_w}$ stems from the fact eigenvalues ${\displaystyle {\lambda }_1,{\lambda }_2}$ (which can be ordered so that ${\displaystyle {\lambda }_1\ge {\lambda }_2\ge 0}$) and the corresponding eigenvectors ${\displaystyle e_1,e_2}$ summarize the distribution of the gradient ${\displaystyle \mathrm{\nabla }I=(I_x,I_y)}$ of ${\displaystyle I}$ within the window defined by ${\displaystyle w}$ centered at ${\displaystyle p}$.

Namely, if ${\displaystyle {\lambda }_1>{\lambda }_2}$, then ${\displaystyle e_1}$ (or ${\displaystyle -e_1}$) is the direction that is maximally aligned with the gradient within the window.

In particular, if ${\displaystyle {\lambda }_1>0,{\lambda }_2=0}$ then the gradient is always a multiple of ${\displaystyle e_1}$ (positive, negative or zero); this is the case if and only if ${\displaystyle I}$ within the window varies along the direction ${\displaystyle e_1}$ but is constant along ${\displaystyle e_2}$. This condition of eigenvalues is also called linear symmetry condition because then the iso-curves of ${\displaystyle I}$ consist in parallel lines, i.e there exists a one dimensional function ${\displaystyle g}$ which can generate the two dimensional function ${\displaystyle I}$ as ${\displaystyle I(x,y)=g(d^{\text{T}}p)}$ for some constant vector ${\displaystyle d=(d_x,d_y{)}^T}$ and the coordinates ${\displaystyle p=(x,y{)}^T}$.

If ${\displaystyle {\lambda }_1={\lambda }_2}$, on the other hand, the gradient in the window has no predominant direction; which happens, for instance, when the image has rotational symmetry within that window. This condition of eigenvalues is also called balanced body, or directional equilibrium condition because it holds when all gradient directions in the window are equally frequent/probable.

Furthermore, the condition ${\displaystyle {\lambda }_1={\lambda }_2=0}$ happens if and only if the function ${\displaystyle I}$ is constant (${\displaystyle \mathrm{\nabla }I=(0,0)}$) within ${\displaystyle W}$.

More generally, the value of ${\displaystyle {\lambda }_k}$, for k=1 or k=2, is the ${\displaystyle w}$-weighted average, in the neighborhood of p, of the square of the directional derivative of ${\displaystyle I}$ along ${\displaystyle e_k}$. The relative discrepancy between the two eigenvalues of ${\displaystyle S_w}$ is an indicator of the degree of anisotropy of the gradient in the window, namely how strongly is it biased towards a particular direction (and its opposite). This attribute can be quantified by the coherence, defined as

$${\displaystyle c_w={\left(\frac{{\lambda }_1-{\lambda }_2}{{\lambda }_1+{\lambda }_2}\right)}^2}$$

if ${\displaystyle {\lambda }_2>0}$. This quantity is 1 when the gradient is totally aligned, and 0 when it has no preferred direction. The formula is undefined, even in the limit, when the image is constant in the window (${\displaystyle {\lambda }_1={\lambda }_2=0}$). Some authors define it as 0 in that case.

Note that the average of the gradient ${\displaystyle \mathrm{\nabla }I}$ inside the window is not a good indicator of anisotropy. Aligned but oppositely oriented gradient vectors would cancel out in this average, whereas in the structure tensor they are properly added together. This is a reason for why ${\displaystyle (\mathrm{\nabla }I)(\mathrm{\nabla }I{)}^{\text{T}}}$ is used in the averaging of the structure tensor to optimize the direction instead of ${\displaystyle \mathrm{\nabla }I}$.

By expanding the effective radius of the window function ${\displaystyle w}$ (that is, increasing its variance), one can make the structure tensor more robust in the face of noise, at the cost of diminished spatial resolution. The formal basis for this property is described in more detail below, where it is shown that a multi-scale formulation of the structure tensor, referred to as the multi-scale structure tensor, constitutes a true multi-scale representation of directional data under variations of the spatial extent of the window function.

Complex version

The interpretation and implementation of the 2D structure tensor becomes particularly accessible using complex numbers. The structure tensor consists in 3 real numbers

$${\displaystyle S_w(p)=\left[\begin{array}{cc}{\mu }_{20}& {\mu }_{11}\\ {\mu }_{11}& {\mu }_{02}\end{array}\right]}$$

where ${\mu }_{20}=\int (w(r)(I_x(p-r){)}^{2}dr$ , ${\mu }_{02}=\int (w(r)(I_y(p-r){)}^{2}dr$ and ${\mu }_{11}=\int w(r)I_x(p-r)I_y(p-r)dr$ in which integrals can be replaced by summations for discrete representation. Using Parseval's identity it is clear that the three real numbers are the second order moments of the power spectrum of ${\displaystyle I}$. The following second order complex moment of the power spectrum of ${\displaystyle I}$ can then be written as

$${\displaystyle {\kappa }_{20}={\mu }_{20}-{\mu }_{02}+i2{\mu }_{11}=\int w(r)(I_x(p-r)+i{I}_y(p-r){)}^{2}dr=({\lambda }_1-{\lambda }_2)\exp (i2\varphi )}$$

where ${\displaystyle i=\sqrt{-1}}$ and ${\displaystyle \varphi }$ is the direction angle of the most significant eigenvector of the structure tensor ${\displaystyle \varphi =\mathrm{\angle }{e}_1}$ whereas ${\displaystyle {\lambda }_1}$ and ${\displaystyle {\lambda }_2}$ are the most and the least significant eigenvalues. From, this it follows that ${\displaystyle {\kappa }_{20}}$ contains both a certainty ${\displaystyle |{\kappa }_{20}|={\lambda }_1-{\lambda }_2}$ and the optimal direction in double angle representation since it is a complex number consisting of two real numbers. It follows also that if the gradient is represented as a complex number, and is remapped by squaring (i.e. the argument angles of the complex gradient is doubled), then averaging acts as an optimizer in the mapped domain, since it directly delivers both the optimal direction (in double angle representation) and the associated certainty. The complex number represents thus how much linear structure (linear symmetry) there is in image ${\displaystyle I}$, and the complex number is obtained directly by averaging the gradient in its (complex) double angle representation without computing the eigenvalues and the eigenvectors explicitly.

Likewise the following second order complex moment of the power spectrum of ${\displaystyle I}$, which happens to be always real because ${\displaystyle I}$ is real,

$${\displaystyle {\kappa }_{11}={\mu }_{20}+{\mu }_{02}=\int w(r)|I_x(p-r)+i{I}_y(p-r){|}^{2}dr={\lambda }_1+{\lambda }_2}$$

can be obtained, with ${\displaystyle {\lambda }_1}$ and ${\displaystyle {\lambda }_2}$ being the eigenvalues as before. Notice that this time the magnitude of the complex gradient is squared (which is always real).

However, decomposing the structure tensor in its eigenvectors yields its tensor components as

$${\displaystyle S_w(p)={\lambda }_1{e}_1{e}_1^{\text{T}}+{\lambda }_2{e}_2{e}_2^{\text{T}}=({\lambda }_1-{\lambda }_2)e_1{e}_1^{\text{T}}+{\lambda }_2(e_1{e}_1^{\text{T}}+e_2{e}_2^{\text{T}})=({\lambda }_1-{\lambda }_2)e_1{e}_1^{\text{T}}+{\lambda }_{2}E}$$

where ${\displaystyle E}$ is the identity matrix in 2D because the two eigenvectors are always orthogonal (and sum to unity). The first term in the last expression of the decomposition, ${\displaystyle ({\lambda }_1-{\lambda }_2)e_1{e}_1^{\text{T}}}$, represents the linear symmetry component of the structure tensor containing all directional information (as a rank-1 matrix), whereas the second term represents the balanced body component of the tensor, which lacks any directional information (containing an identity matrix ${\displaystyle E}$). To know how much directional information there is in ${\displaystyle I}$ is then the same as checking how large ${\displaystyle {\lambda }_1-{\lambda }_2}$ is compared to ${\displaystyle {\lambda }_2}$.

Evidently, ${\displaystyle {\kappa }_{20}}$ is the complex equivalent of the first term in the tensor decomposition, whereas $${\displaystyle \frac{1}{2}(|{\kappa }_{20}|-{\kappa }_{11})={\lambda }_2}$$is the equivalent of the second term. Thus the two scalars, comprising three real numbers,

$${\displaystyle \begin{array}{cccc}{\kappa }_{20}& =& ({\lambda }_1-{\lambda }_2)\exp (i2\varphi )& =w\ast (h\ast I{)}^2\\ {\kappa }_{11}& =& {\lambda }_1+{\lambda }_2& =w\ast |h\ast I{|}^2\end{array}}$$ where ${\displaystyle h(x,y)=(x+iy)\exp (-(x^2+y^2)/(2{\sigma }^2))}$ is the (complex) gradient filter, and ${\displaystyle \ast }$ is convolution, constitute a complex representation of the 2D Structure Tensor. As discussed here and elsewhere ${\displaystyle w}$ defines the local image which is usually a Gaussian (with a certain variance defining the outer scale), and ${\displaystyle \sigma }$ is the (inner scale) parameter determining the effective frequency range in which the orientation ${\displaystyle 2\varphi }$ is to be estimated.

The elegance of the complex representation stems from that the two components of the structure tensor can be obtained as averages and independently. In turn, this means that ${\displaystyle {\kappa }_{20}}$ and ${\displaystyle {\kappa }_{11}}$ can be used in a scale space representation to describe the evidence for presence of unique orientation and the evidence for the alternative hypothesis, the presence of multiple balanced orientations, without computing the eigenvectors and eigenvalues. A functional, such as squaring the complex numbers have to this date not been shown to exist for structure tensors with dimensions higher than two. In Bigun 91, it has been put forward with due argument that this is because complex numbers are commutative algebras whereas quaternions, the possible candidate to construct such a functional by, constitute a non-commutative algebra.

The complex representation of the structure tensor is frequently used in fingerprint analysis to obtain direction maps containing certainties which in turn are used to enhance them, to find the locations of the global (cores and deltas) and local (minutia) singularities, as well as automatically evaluate the quality of the fingerprints.