# Why tensors?

### Contents

Images in *DIPlib 3* have pixel values generalized to tensors. Currently only
tensors of order up to 2 are supported. Thus, pixel values can be scalar (0^{th}-order
tensor), vector (1^{st}-order) or matrix (2^{nd} order). This covers all current
applications of multi-valued pixels, as far as we’re aware. On this page we show
some examples of these applications that show the notational simplicity that comes
with the tensor representation. In each of the examples below, please think about
how an equivalent program would look if one did not have access to a tensor
representation.

## Stain unmixing

In brightfield microscopy it is well-known that multiple stains on a slide generate a linear mixture of those stains’ color vectors (see Ruifrok & Johnston, Anal Quant Cytol Histol 23(4):291-9, 2001, who unfortunately called it “color deconvolution”, a name that has stuck even though this has nothing to do with deconvolution). Given, for example, two stains imaged in an RGB image, then the absorption in each of the R, G and B channels is:

where are the optical densities of the stains at a particular point, is the absorption of stain in the red channel, etc. The actual intensity measured by the camera is the transmittance , given by

(as a fraction of the incident light intensity); see Beer-Lambert law.

Thus, given `img`

, an RGB image recorded in a brightfield microscope, the intensity
of the background , and the mixing matrix , you would follow the steps above
in reverse:

dip::Image I{ I_R, I_G, I_B }; dip::Image t = dip::Log10( -img / I ); dip::Image S{ S_1R, S_1G, S_1B, S_2R, S_2G, S_2B }; S.ReshapeTensor( 3, 2 ); dip::Image U = dip::PseudoInverse( S ); dip::Image d = U * t;

Note that the most verbose part of the code above is to create the 0D (single-pixel)
images `S`

and `I`

, corresponding to the mixing matrix and the background intensity
vector . These are a single pixel because the same properties are valid across
the image. In the arithmetic operations, singleton-expansion causes these 0D images
to be expanded to the size of the image `img`

, such that each pixel is multiplied
by the same matrix (though this expansion happens implicitly, there’s not actually such an
image created in memory). Next, note that the operation `img / I`

is a per-element
division, such that `img[0]`

(red channel) is divided by `I[0]`

, `img[1]`

by `I[1]`

,
and `img[2]`

by `I[2]`

. In contrast, `U * t`

is a matrix multiplication.

The above is implemented in the functions `dip::BeerLambertMapping`

and `dip::UnmixStains`

.

## The structure tensor

The structure tensor (see Wikipedia) is defined by the gradient of the image:

where is the partial derivative of the image in the direction, and the overline indicates local averaging. The equation above shows the structure tensor for a 2D image, but the expression on the left is valid for any number of dimensions. The eigenvalues of describe the local neighborhood of a pixel. If we sort the eigenvalues such that , then it is possible to define a measure of anisotropy as or simply as , and an energy measure as . Anisotropy is high on lines and edges, and small on uniform areas. Furthermore, the direction of the eigenvector corresponding to the largest eigenvalue gives the orientation perpendicular to the line or edge.

The anisotropy measure generalizes to higher-dimensional images. For example, in a 3D image it is possible to distinguish lines (two large eigenvalues), planes (one large eigenvalue), and isotropic areas (all eigenvalues approximately equal).

To compute the structure tensor using *DIPlib*, we start with a scalar image `img`

:

dip::Image g = dip::Gradient( img ); dip::Image S = dip::Gauss( g * dip::Transpose( g ), 5 );

The image `S`

is a 2-by-2 tensor image if `img`

is a 2D image, or a 5-by-5 tensor
image if `img`

is a 5D image. Note that, since `g * Transpose(g)`

yields a
symmetric matrix, `S`

is stored in such a way that no duplicate matrix elements
are stored (nor computed). That is, for the 2D case, `S`

has 3 tensor elements,
corresponding to the two diagonal elements and the one unique off-diagonal element.
The constant 5 in the computation of `S`

is the sigma of the Gaussian smoothing.
This value needs to be adjusted depending on the scale of the structures we are
interested in.

dip::Image e = dip::Eigenvalues( S ); dip::Image anisoptropy = ( e[0] - e[1] ) / ( e[0] + e[1] );

Note that it is possible to compute the `anisotropy`

image more efficiently
(running through the image `e`

once instead of three times, and avoiding
the storage of two intermediate images). See Iterators.

The above is implemented in the functions `dip::StructureTensor`

and
`dip::StructureTensorAnalysis`

.

## The Harris corner detector

The well-known Harris corner detector (see Wikipedia) is based on the structure tensor. Pixels where the two eigenvalues are large are corners. It is common to use the following equivalence to avoid computation of the eigenvalues:

We can compute this in *DIPlib* as follows, again starting with a scalar
image `img`

:

dip::Image g = dip::Gradient( img ); dip::Image S = dip::Gauss( g * dip::Transpose( g ), 5 ); dip::Image trace = dip::Trace( S ); dip::Image corners = dip::Determinant( S ) - k * trace * trace;

Note again that is is possible to compute `corners`

more efficiently
by Iterators. That way, one can run through the image `S`

only once,
and avoid the temporary intermediate images.

The above is implemented in the function `dip::HarrisCornerDetector`

.

## Lucas-Kanade optical flow

The Lucas-Kanade solution to the optical flow problem (see
Wikipedia)
also involves the structure tensor. The problem to solve is ,
where is a matrix with and partial derivatives at each
pixel in a neighborhood, is a vector with derivatives at each
of those pixels, and is the velocity vector for the neighborhood.
This is rewritten as , where is the
structure tensor. Using *DIPlib*, and assuming a 3D image `img`

where
the third dimension is time, we can write:

dip::Image A = dip::Gradient( img, { 1.0 }, "best", {}, { true, true, false } ); dip::Image b = -dip::Derivative( img, { 0, 0, 1 } ); dip::Image ATA = dip::Gauss( A * dip::Transpose( A ), 5 ); dip::Image ATb = dip::Gauss( A * b, 5 ); dip::Image v = dip::Inverse( ATA ) * ATb;

The most complicated function call is on the first line. Up to now we had
been using all the default parameters to `dip::Gradient`

, but now we need
to set the `process`

parameter: a boolean array indicating along which
dimensions to compute the derivative. Since this parameter is towards the
end of the parameter list, we must fill out the other default values,
which we simply copied from the function declaration.

## Anisotropic diffusion

There are many more examples where per-pixel matrix algebra is useful and
*DIPlib* allows simple, efficient implementation. The last example we’ll
give here is from the function `dip::CoherenceEnhancingDiffusion`

.

The diffusion equation can be discretized along the time axis to yield an iterative update process described by

If is constant (spatially and in time), the above iterative process leads to Gaussian smoothing. By adjusting to be small at edges, an anisotropic, edge-enhancing diffusion is obtained. Coherence enhancing diffusion uses a tensor for . The above-mentioned function allows to construct this tensor in two ways, starting from (yet again!) the structure tensor . Here we apply an eigen decomposition, leading to a full matrix image (the eigenvectors), and a diagonal matrix image (the eigenvalues):

Next we compute and , and re-compose a tensor using these new eigenvalues:

Using *DIPlib* we can write:

dip::Image S = dip::StructureTensor( img ); // see above for how this is computed dip::Image E, V; dip::EigenDecomposition( S, E, V ); E = 1 / E; E /= dip::Trace( E ); dip::Image D = V * E * dip::Transpose( V ); img += lambda * dip::Divergence( D * dip::Gradient( img ));

Iterating this bit of code leads to a coherence enhancing diffusion simulation
on the image `img`

.