Scalespace or Differential Geometry

Introduction

Differential geometry is geometry done using differential calculus or shape description through derivatives. It applies techniques from multivariable calculus and linear algebra to study features of curves and surfaces of dimension 1, 2, and higher. Although this article is limited to 2D, the differential geometry approach can be used for 1D, 2D, 3D and 4D signal and image processing and the digital elevation models (DEM) of geomorphometry. The basic principle is "what you describe is what you get".

This webpage gives an overview grayscale scale space or differential geometry as it can be used in digital image processing, from digital microscopy and High Content Screening (HCS) to medical imaging, geology, and astronomy or any other spatial and temporal structure in a dynamic space-time continuum.
The geometric properties of differential gray-value invariants allow for the selection of structural image features for analysis and quantification. Although this article is limited to 2D, the scale space approach can be used for 3D and even 4D image processing. Scalespace stacks (dimension, scale, order) can also be used in convolutional neural networks, such as in the M⁵ framework for exploring the cytome. The basic principle is "what you describe is what you get".
For an introduction to the spatial color model see color scale space, which is another application of scale space applied to color image processing in microscopy and High Content Screening.

This webpage was inspired by the pioneering work on High Content Screening and microscopy at the now defunct Biological Imaging laboratory (BIL) at Janssen Pharmaceutica in Beerse (Belgium).

Derivatives of a 1D function

A parabole y=x²+x and its first (f_x=2*x+1) and second order (f_xx=2) derivatives.

The principal interpretations of derivatives of a function f concern instantaneous rate of change and tangency. In mathematics the derivatives of a one-dimensional function f are denoted as:
f_x=dy/dx = 1^st order derivative (the slope of the function)
f_xx=d²y/dx² = 2^nd order derivative (the curvature of the function, or the rate of change of the slope)
f_xxx=d³y/dx³ = 3^rd order derivative (the rate of change of the curvature of the function)

At each point x at which the function f(x) is decreasing the slope of the tangent is negative -- that is, f_x is negative.
At each point x at which the function f(x) is increasing the slope of the tangent is positive -- that is, f_x is positive.
At each point x which is the top of a mountain or the bottom of a valley the tangent is horizontal, so its slope is zero -- that is, f_x = 0.
- If in addition f_xx > 0, then f(x) is at a minimum.
- If in addition f_xx < 0, then f(x) is at a maximum.

An image as a 2D function

3D plot
A gray scale 2D image and the gray values (intensities) plotted in a 3D view. The x,y axis is the position of the pixel in the 2D plane of the image, while the z-axis is the gray value or intensity of the pixel.

A 2D image can be visualized as a 3D landscape by plotting the surface height proportional to the pixel gray value. The principles of differential geometry can be used on 2D images, which becomes clear if we plot the intensity profile of the image in 3D view.
The intensity profile of the image constitutes a landscape of (dark) valleys and (bright) hills and ridges, which form the curves of a 2D function. The concepts of image processing become similar to the concepts being used in terrain analysis. Terrain analysis employs elevation data to describe the landscape, for basic visualization, modeling, or to support decision making. Geomorphometry is 'the science of quantitative land surface analysis or topographic quantification; its operational focus is the extraction of land-surface parameters and objects from digital elevation models (DEM)' (see Geomorphometry: Concepts, Software, Applications, Tomislav Hengl, Hannes I. Reuter, Elsevier, 2009, p. 4). Digital elevation models resemble a 3D plot of the intensity profile of a 2D image, which makes studying the intensity profile of an image similar to studying the topography of a terrain. The basic notions of geomorphometry and scale space consist of concepts such as:

Uphill/downhill
Contour lines (curves of constant elevation)
Steepness of slope
Peaks/Valleys (local extrema)

These translate into mathematical concepts us as:

Tangent plane
Normal vectors
Curvature

and gradient vectors (vectors of partial derivatives) which define/compute these features.

The 3D luminance profile of a 2D image (L) can be seen as a luminance landscape. L(x, y) is the luminance function (pixel values) of an image and the first order derivative with respect to x is Lx, and the derivative with respect to y is Ly. The geometric interpretation of the first order derivative or gradient at x₀ is the slope of the tangent line to the curve at point x₀: slope = Δy/Δx. The first derivative tells us how whether a function is increasing or decreasing, and by how much it is increasing or decreasing. The second derivative (curvature or concavity) with respect to x is Lxx, the second derivative with respect to y is Lyy.

Besides using a Cartesian coordinate system we can also use gauge coordinates such as the gradient or curvature gauge. Gauge coordinates are coordinates with respect to axes that are attached to the function geometry. A gauge coordinate system therefore is independent of the natural coordinate axes.
The differentiation of the 2D image is done by convolving the original image with the derivatives of a 2D Gaussian kernel or, equivalently, by multiplying the Fourier transform of the original image with the Fourier transform of the Gaussian and doing the reverse Fourier transform of the result.

The same principles can be extended to 3D images, but this is beyond the current scope of this webpage. The principles of differential geometry can also be used for the analysis of spacetime phenomena, but this is also beyond the scope of this article.

The Gaussian kernel

1D Gauss kernel
1D Gaussian, probability density function (PDF) of the normal distribution with standard deviation σ and variance σ²

2D Gauss kernel
2D Gaussian kernels, up to 2^nd order derivatives, bivariate normal distribution and equal standard deviation σ =σ_x=σ_y

The Gaussian kernel is the unique kernel for differentiation for an uncommitted front-end (eg. human eye). Scale-space uses the Gaussian function and its derivatives as a physical differential operator. In computational vision we replace the mathematical derivatives by fuzzy derivatives or Gaussian Derivatives. Differentiating an image is done by convolving the image in the spatial domain with a Gaussian kernel with a standard deviation σ. An alternative way is using the Fourier transform of a convolution, which is the pointwise product of Fourier transforms of the image and the Gaussian kernel. In other words, convolution in one domain (e.g., spatial domain) equals point-wise multiplication in the other domain (e.g., frequency domain) (see convolution theorem). The Gaussian kernel we will use here is isotropic, which means that the behavior of the function is in any direction the same. For 2D this means the Gaussian function is circular, for 3D it would look like a fuzzy sphere. When we would use unequal standard deviations in the different dimensions we would apply an anisotropic Gaussian kernel. The σ (sigma) parameter of the Gaussian kernel determines the size (scale) of the features that are resolved, ie. it is a measure for the aperture of the filter. Local image observations with Gaussian Derivatives are always done at a finite scale (σ) to be selected by the observer. The larger you make the standard deviation σ of the Gaussian kernel, the more the image gets blurred. In order to find the optimal size of the Gaussian filter we may apply many filters of different sizes and embed them in a structure called a (Gaussian) scale space pyramid. An exponential increase in scale was proposed by Koenderink in 1984 and a method for automatic scale detection was proposed by Lindeberg in 1998.

Notation and conventions

Evaluation of differential expressions

L is the image
L(x, y) is the luminance function (pixel values) of a 2D image
Lx = 1^st order x derivative
Ly = 1^st order y derivative
Lxx = 2^nd order x derivative, etc.

Li = 1st order Partial Differential Equation (PDE) according to summation convention
"i" sums over all dimensions, i.e. Li = Lx+Ly
Li is a vector:{Lx,Ly,Lz}. Summation convention on paired indices

Lii = 2^nd order, Lii = Lxx+Lyy
Lij = 2^nd order, Lij = Hessian matrix

Gradient and curvature gauge

The classical Euclidean image coordinate frame is an arbitrary convenient choice without much relevance. Attaching the coordinate frame to the local image structure makes the coordinate frame independent of arbitrary Euclidean transformations (rotations, scalings and translations).

The most important gauge coordinate systems are:

Gradient gauge (v, w) and derivatives: Lv, Lw, Lvw, etc.
Curvature gauge (p, q) and derivatives: Lp, Lq, Lpq, etc.

More information on gauge coordinates.

Derivatives of the Gaussian kernel

Zero order derivative

Zero order Gaussian
Zero order derivative of a 1D Gaussian kernel, sigma (σ) is 10.0.

First order derivative

First order Gaussian
First order derivative of a 1D Gaussian kernel, sigma (σ) is 10.0.

Second order derivative

Second order Gaussian
Second order derivative of a Gaussian kernel, sigma (σ) is 10.0.

Third order derivative

Third order Gaussian
Third order derivative of a 1D Gaussian kernel, sigma (σ) is 10.0

Differentiating with the Gaussian kernel

Convolving an image with a Gaussian kernel and derivatives of the Gaussian kernel, highlights features of the image with a size larger than the chosen scale sigma (σ).
The examples shown here are for 2D images, but can be extended to 3D (spatial), 4D (spatiotemporal) and spatiospectral dimensions.

Original image
The original 256x256 pixels image L(x, y) from which the following differentiated images are derived.

Zero order derivative Zero order x derivative, 3D Plot
Zero order Gaussian derivative in x and y direction (Lxy), just blurs the image. Sigma (σ, scale) is 2.
Original resolution is 256²

First order x derivative First order x derivative, 3D Plot
First order Gaussian derivative in x direction (Lx). Sigma (σ, scale) is 2.

Second order x derivative Second order x derivative, 3D Plot
Second order Gaussian derivative in x direction (Lxx). Sigma (σ, scale) is 2.

Invariants

Invariant differential feature detectors are polynomial combinations of image derivatives, which exhibit invariance under some chosen group of transformations.

2D grayscale

First Order

Intensity gradient 3D plot of intensity gradient
The image on the left shows the response intensity of the Canny filter (bright means high gradient strength), which is also shown in 3D on the right, the scale, sigma (σ), is 3.

Intensity gradient = Canny edge detector = ||∇L|| = Lw
For an explanation of Lv, Lw, etc. see the section on gauge coordinates.

Second Order

Second order derivatives yield isophote and flowline curvature and cornerness of an image. Isophotes are lines drawn through areas of constant brightness while flowlines are lines drawn through areas of changing brightness.

Gradient² = L_iL_i = L_x² + L_y²
Marr edge detector = ∇²L = L_xx + L_yy
Canny edge detector = L_ww ≈ 0
Laplacian of Gaussian = Trace of Hessian matrix = trH_f = L_ii = L_xx + L_yy
Determinant of the Hessian matrix (blob detection) = |H_f| = L_xxL_yy-L_xy²
Euclidian shortening flow = Lvv = Ls
Ridge strength = Lvv
Derivative of gradient in gradient direction = Lww (if Lww=0, Lw maximal)
Deviation from flatness = D = Sqrt(tr H²) = Sqrt(L_ijL_ji) = Sqrt(L_xx² + 2L_xy² + L_yy²)
Isophote curvature = K = -Lvv / Lw
Flow-line curvature = -Lvw / Lw
Isophote density per unit luminance = -Lww / Lw

Principal intensity surface curvature in principal coordinates in 2D:
K_p1 = L_xx/Sqrt(L_x² + L_y²)
K_p2 = L_yy/Sqrt(L_x² + L_y²)

Inflexion strength = LvvLw⁵-3LvvLvwLw⁴

Umbilicity = U = (trH)²/TrH² - 1 = (L_xx + L_yy)² / D² = L_ll² / L_ijL_ji

Pure hyperbolic: U = -1
Hyperbolic region: -1 < U < 0
Parabolic region: U = 0
Elliptic region: 0 < U < 1
Spherical point: U = 1

Elongatedness = E = ((trH)²(2trH²-(trH)²))/(2trH²)² = (L_xx² + L_yy²)/D⁴

Higher Order

Third order derivatives provide T-junction detection. T-junctions in an image occur at occlusion points which occur at a point where a contour ends or emerges because there is another object in front of the contour.

T-junction = ((LvvLvvLvwLvw+LvvLvvLwwLww)/Lw²-2LvvLvvvLvw/Lw-2LvvLvvwLww/Lw+(LvvvLvvv+LvvwLvvw))/Lw²

Resultants and Discriminants

All resultants and discriminants are Euclidian and affine invariant

R₁₂ = Corner strength = L_x²L_yy-2L_yL_xL_xy+ L_y²L_xx = LvvLw²

D₁ = 1
D₂ = detH = -L_xy²+L_yy + L_xx or L_xx+L_yy-L_xy²

Partial Derivative Equations

A partial differential equation (PDE) is a differential equation that contains unknown multivariable functions and their partial derivatives. PDEs can be used to describe a wide variety of phenomena in multidimensional systems such as images.

Extrema

Lx = 0
Ly = 0
LxxLyy-Lxy² >= 0
|Lxx| > threshold
threshold determines salient features

Maxima

Lx = 0
Ly = 0
LxxLyy-Lxy² >= 0
Lxx < -threshold
threshold determines salient features

Minima

Lx = 0
Ly = 0
LxxLyy-Lxy² >= 0
Lxx > threshold
threshold determines salient features

Saddle points

Lx = 0
Ly = 0
LxxLyy-Lxy² < 0

Parabolic points

Lx = 0
Ly = 0
LxxLyy-Lxy² = 0

Corners

affine corner detection: LwwLw²
with high curvature: LwwLw² > threshold

Ridges

Bright ridges: Lqq < -4.0 * threshold / sigma²
Dark ridges: Lpp > 4.0 * threshold / sigma²
threshold determines salient features

Gauge coordinates

A convenient way to calculate features invariant under Euclidean coordinate transformations is done by using gauge coordinates such as the gradient or curvature gauge. Any combination of derivatives with respect to gauge coordinates is invariant under Euclidean transformations.

Coordinate axes in space can be chosen arbitrarily. The space and its properties do not change when we change the coordinate axes. But if coordinate systems are arbitrary, what then is the value of derivatives in the axis direction? The individual derivatives (f_x, f_xy, ...) are of no use. It is the collection of all derivatives in the N-jet (collection of all spatial image derivatives up to order N) that are of importance as it has geometrical meaning independent of the arbitrary coordinate system.
Gauge coordinates are coordinates with respect to axes that are attached to the function geometry. A gauge coordinate system therefore is independent of the natural coordinate axes.
Examples of gauge systems are the gradient gauge and the principal curvature gauge, in which all points in the image have their own gauge coordinate system attached to it.

The word gauge (in gauge theory) was introduced as the German word maßstab by Hermann Weyl (1885-1955 CE). In 1918, he introduced the notion of gauge, and gave the first example of what is now known as a gauge theory.

The classical Euclidean image coordinate frame is an arbitrary convenient choice. Attaching the coordinate frame to the local image structure makes the coordinate frame independent of arbitrary Euclidean transformations (rotations, scalings and translations).
The most important gauge coordinate systems are:

Gradient gauge (v, w) where one of the coordinate axes (w) is aligned with the image gradient vector (flowline).
Curvature gauge (p, q) where the coordinate axes are aligned are aligned with the direction of minimal and maximal directional second order derivative.

The Gradient Gauge

Gradient gauge coordinates, where w is the coordinate in the gradient direction (flowline, maximal change of the intensity)
and v is the coordinate perpendicular to the gradient direction (isophote).

A gauge is a "coordinate system" that varies depending on one's "location" with respect to some "base space" or "parameter space", a gauge transform is a change of coordinates applied to each such location, and a gauge theory is a model for some physical or mathematical system to which gauge transforms can be applied (and is typically gauge invariant, in that all physically meaningful quantities are left unchanged (or transform naturally) under gauge transformations).
Consider intrinsic geometry: every point is described in such a way, that if we have the same structure, or local landscape form, no matter the rotation, the description is always the same. This can be accomplished by setting up in each point a dedicated coordinate frame which is determined by some special local directions given by the landscape locally itself.
In each point separately the local coordinate frame is fixed in such a way that one frame vector points to the direction of maximal change of the intensity, and the other perpendicular to it (90 degrees clockwise). We can now fix locally the direction for our new intrinsic local coordinate frame (v and w in the image) and this set of local directions is called a gauge, the new frame forms the gauge coordinates.

A special notation is used in image processing for this gauge system: the w,v-coordinates. w is the coordinate in the gradient direction (maximal change of the intensity) and v is the coordinate perpendicular to the gradient direction.

{v,w} are gradient gauge coordinates, w along and v perpendicular on the gradient; w and v are unit vectors. Lw is the derivative of the luminance field in the w direction.
Lv = 0, Lw = |Grad L|
Lvv, Lwv, Lww = 2^nd order, etc.
Lww is also called the SDGD (Second Derivative in Gradient Direction).
Lw gives us information about the local change in intensity of an image, which is a first order derivative phenomenon.

The Curvature Gauge

Curvature Gauge

In situations where the gradient gauge vanishes (Lx=0, Ly=0), the second order differentials may indicate that something interesting is happening, for instance the points along the axis of a line. A gauge system suitable for this phenomenon is the curvature gauge. The eigenvectors of the Hessian matrix H_L then provide a powerful way of defining a gauge system.

The local second order differentials around a point in a 2D image can be grouped into a Hessian matrix in which all second order derivatives are gathered, it describes the local second-order profile of intensity variations in a 2D image.

The partial second order derivatives of a 2D image L(x,y) are represented like L_xx, L_xy, etc. .
The Hessian matrix of second order derivatives is a symmetric matrix and can be transformed to a diagonal matrix, using a rotation matrix. The coordinate frame that diagonalizes the Hessian matrix is often denoted as the (p,q)-frame, it is built from the eigenvectors of the Hessian matrix.
The eigenvectors of the Hessian, p and q, point in the directions of the minimal and maximal directional second order derivatives.
The eigenvalues of the Hessian are the second order differentials f_pp and f_qq or Lpp and Lqq.
{p,q} are the Hessian gauge coordinates, this diagonalizes the Hessian (eigenvalues of the Hessian matrix), by definition Lpq = 0

The trace (trH_f) of a matrix equals the sum of its diagonal elements, it is a sort of Laplacian, it might be called the mean curvature and it describes the curvature at patch boundaries.
The sign of the trace indicates whether creases are ridges or valleys and its magnitude reflects the angle at which two linear patches meet.
The trace of the Hessian=trH_f=f_xx+f_yy= Lxx+Lyy=Lww+Lvv

The determinant (|H_f|) of the Hessian matrix is:
|H_f|=f_xxf_yy-f_xy²= LxxLyy-Lxy²=LwwLvv-Lvw²

The determinant is sometimes called deviation of flatness, and it indicates the overall change in first-order information; magnitude and direction.

The term Hessian was coined by James Joseph Sylvester (1814-1897 CE), named for Otto Hesse (1811-1874 CE), who had used the term functional determinants.

Taylor series of a function

Locally an image is characterized with its Taylor expansion. A way to approximate an image is to perform a Taylor expansion of the local image structure. Locally, i.e. in a small neighborhood of a point, an image (as any function) can be represented by its Taylor expansion. A Taylor series is a series expansion of a function, it is the most basic tool for approximation of functions. It is a representation of a function as an infinite sum of terms that are calculated from the values of the function's derivatives at a single point.

A Taylor expansion is useful for calculating image values between pixels and extrapolation. Truncating the Taylor expansion at some finite number N, we obtain an approximated representation of the function (image). Since the complete expansion (N=infinity) is a unique description of the function f, the approximated expansion can be the same for different images.
The local N-jet is definded as the complete set of images which all give the same Taylor expansion at that point up to order N. This is the set of images which are equal in the first N derivatives at that point.
The scale has be to taken into account in the Taylor-expansion, The scale is the level of resolution at which we look at an image, only resolving features at a certain level of resolution (the aperture of the filter in image analysis). All derivatives are scaled derivatives.

Taylor series of a 1D-function

Consider a one-dimensional function f(x) that is smooth - when we say 'smooth', what we mean is that its derivatives exist and are bounded (for the following discussion, we need f to have (n+1)derivatives). We would like to approximate f(x) near the point x=a, and we can do it as follows with a 1-D Taylor series:

If we look at the infinite sum (let n-> infinity), then the resulting infinite sum is called the Taylor series of f(x) about x=a.
This result is also know as Taylor's Theorem.

If a=0, the expansion is known as a Maclaurin series, which is a Taylor series expansion of a function about 0. The third order Taylor expansion (and beyond) becomes:

Maclaurin series

An alternative form of the 1-D Taylor series may be obtained by letting

so that

Substitute this result into () to give

Taylor series of a 2D-function

We want to generalize Taylor expansions to functions of more than one variable so we can use it in image processing and analysis. Our approach will be limited to 2D images, but can be expanded to 3D images also. Let us consider a two-dimensional function f(x, y). We want to approximate f(x₀+Δx, y₀+Δy) for Δx and Δy near zero. We rewrite f(x₀+Δx, y₀+Δy) = F(1) as F(t) = f(x₀+tΔx, y₀+tΔy) where x₀, y₀, Δx and Δy are constants, so that F is now a function of a single variable t. Now we can use out single variable formula with t₀=0 and Δt=1 (e.g. Maclaurin series). For this we need to compute various derivatives of F(t) at t = 0, by applying the chain rule to our formula.

A Taylor series of a function in two variables f(x,y), such as a 2D-image, can be given by expanding f(x+Δx, y+Δy):

Let us now take this step by step. Given $f(x,y)$ we want the formula for the Taylor series expansion around $x=x_0$ and $y=y_0$ .
First, we rewrite Δx as α₀ and Δy as α₁ and parameterise:
$x(t) = x_0 + \alpha_0 t$ and $y(t) = y_0 + \alpha_1 t$ .
Then, because Δx, Δy, α₀ and α₁ are constants, we now have a function of a single variable t:
$f(x,y) = f(x(t), y(t)) = f(t)$

The Taylor series expansion of $f(t)$ at $t = 0$ (Maclaurin series):
$f(t) = f(0) + t f'(0) + \dfrac{t^2}{2} f''(0) + ...$

Applying the chain rule we get:
$\dfrac{df}{dt} = \dfrac{\partial f}{\partial x} \dfrac{dx}{dt} + \dfrac{\partial f}{\partial y} \dfrac{dy}{dx} = \alpha_0 \dfrac{\partial f}{\partial x} + \alpha_1 \dfrac{\partial f}{\partial y}$

$\dfrac{d^2 f}{dt^2} = \dfrac{d}{dt} \left ( \dfrac{df}{dt} \right ) = \dfrac{\partial}{\partial x} \left ( \alpha_0 \dfrac{\partial f}{\partial x} + \alpha_1 \dfrac{\partial f}{\partial y} \right ) \dfrac{dx}{dt} + \dfrac{\partial}{\partial y} \left (\alpha_0 \dfrac{\partial f}{\partial x} + \alpha_1 \dfrac{\partial f}{\partial y} \right )\dfrac{dy}{dt} = \alpha_0^2 \dfrac{\partial^2 f}{\partial x^2} + 2 \alpha_0 \alpha_1 \dfrac{\partial^2 f}{\partial x \partial y} + \alpha_1^2 \dfrac{\partial^2 f}{\partial y^2}$

Substituting this in the Taylor expansion gives:
$f(t) = f(0) + t \alpha_0 \dfrac{\partial f}{\partial x} + t \alpha_1 \dfrac{\partial f}{\partial y} + \dfrac{1}{2} \left [ \alpha_0^2 t^2 \dfrac{\partial^2 f}{\partial x^2} + 2 \alpha_0 \alpha_1 t^2 \dfrac{\partial^2 t}{\partial x \partial y} + \alpha_1^2 t^2 \dfrac{\partial^2 f}{\partial y^2} \right ] + ...$

Since $\alpha_0 t = x - x_0$ and $\alpha_1 t = y - y_0$ the Taylor expansion now becomes:
$f(t) = f(0) + (x- x_0) \dfrac{\partial f}{\partial x} + (y - y_0) \dfrac{\partial f}{\partial y} + \dfrac{1}{2} \left [ (x - x_0)^2 \dfrac{\partial^2 f}{\partial x^2} + 2 (x - x_0)(y - y_0) \dfrac{\partial^2 f}{\partial x \partial y} + (y - y_0)^2 \dfrac{\partial^2 f}{\partial y^2} \right ] + ...$

and finally:
$f(x, y) = f(x_0, y_0) + (x- x_0) \dfrac{\partial f}{\partial x} + (y - y_0) \dfrac{\partial f}{\partial y} + \dfrac{1}{2} \left [ (x - x_0)^2 \dfrac{\partial^2 f}{\partial x^2} + 2 (x - x_0)(y - y_0) \dfrac{\partial^2 f}{\partial x \partial y} + (y - y_0)^2 \dfrac{\partial^2 f}{\partial y^2} \right ] + ...$

Practical applications in image analysis

A Taylor series approximates continuous functions by using derivative information of different orders and in the same way the image intensity function of an image could be approximated. By using diferent orders of derivatives all possible shapes of the image could be classified and detected. The Taylor approximation is being done by fitting a polynomial to the pixel intensities in a patch or region of the image. By using a "sliding window" algorithm for patch-based operations we try to determine the fit with the Taylor expansion up to the n^th order. At each postion of the sliding window we assign the polynomial's derivatives to the pixel at the window's center. The window slides over the entire image, thereby creating a new intensity plot proportional to the response of the polynomial to each patch in the image. The aperture of the Gaussian filter defines the size of the patch and determines the size of the features which will be detected. Keep in mind, this approach can also be used for higher dimensional spatial and spatiotemporal data.

Minima and Maxima

A necessary condition that a function f assumes a relative maximum or a relative minimum at (x₀, y₀) is that:
f_x = f_y = 0 at (x₀,y₀)
If this condition is true then in summary the following criteria can be used.
All subsequent images were analyzed with sigma (σ)= 5.0.

Relative Maximum

Relative maximum

f_x = 0, f_y = 0 and f_xxf_yy > f²_xy
or i.o.w.
f_x = 0, f_y = 0 and f_xxf_yy - f²_xy > 0
or i.o.w.
Lx = 0, Ly = 0 and LxxLyy - Lxy² > 0

Relative Minimum

Relative minimum

f_x = 0, f_y = 0 and f_xxf_yy < f²_xy
or i.o.w.
f_x = 0, f_y = 0 and f_xxf_yy - f²_xy < 0
or i.o.w.
Lx = 0, Ly = 0 and LxxLyy - Lxy² < 0

Saddle Point

Saddle point

f_x = 0, f_y = 0 and f_xxf_yy = f²_xy
or i.o.w.
f_x = 0, f_y = 0 and f_xxf_yy - f²_xy = 0
or i.o.w.
Lx = 0, Ly = 0 and LxxLyy - Lxy² = 0

Elliptic Patches

Relative Maximum Region

Relative maximum region

f_xx < 0 and f_xxf_yy > f²_xy
or i.o.w.
f_xx < 0 and f_xxf_yy - f²_xy > 0
or i.o.w.
Lxx < 0 and LxxLyy - Lxy² > 0
or i.o.w.
Lww < 0 and LvvLww - Lvw² > 0

Relative Minimum Region

Relative minimum region

f_xx > 0, f_xxf_yy > f²_xy
or i.o.w.
f_xx > 0, f_xxf_yy - f²_xy > 0
or i.o.w.
Lxx > 0 and LxxLyy - Lxy² > 0
or i.o.w.
Lww > 0, LvvLww - Lvw² > 0

Saddle Region

Saddle point region

f_xxf_yy < f²_xy
or i.o.w.
f_xxf_yy - f²_xy < 0
or i.o.w.
LxxLyy - Lxy² < 0
or i.o.w.
LvvLww - Lvw² < 0

Comment

By using for instance Lxx/L in the above mentioned equations, a normalisation relative to the image intensity is obtained.
By selecting for instance Lxx/L > Threshold or Lxx/L < Threshold, a noise reduction is performed, as only salient features are selected.

Tracing lines and ridges

Dark ridge differentiated
Tracing the neurites of a neuron (left) with differential geometry as an example of a line detector (result on the right).
The equation used for detecting a dark ridge, with sigma (σ)=2.0 and threshold 1.0:
Lpp > 4.0 * threshold / σ²

The purpose of a line filter is the discrimination of line structures from other image structures.
A line or a bar is a second order differential structure (see Hessian matrix). The points along the axis have a gradient near zero. A gauge system that is capable of capturing these local phenomena is the above mentioned principal curvature gauge.
As we have seen, the eigenvectors of the Hessian matrix, p and q, point in the directions of the minimal and maximal directional second order derivatives.
The difference between Lqq and Lpp determines the directionality of the local image structure.
If Lqq = Lpp, there is no preferred orientation; if Lqq >> Lpp this indicates a line structure.
In short, we have a line when Lqq - Lpp > Th and Lx = 0, Ly = 0. The parameter Th determines the saliency of the line (see also C. Steger (1998)).

The larger the difference between the trace of the Hessian (trH_f) and the determinant (|H_f|), the larger is the directionality of the local 2D image structure and as such the probability for a line.

With respect to the Cartesian coordinate frame (x,y):

p = (f_xx-f_yy- sqrt((trH_f)²-4*|H_f|)) / (2*f_xy)

q = (f_xx-f_yy+ sqrt((trH_f)²-4*|H_f|)) / (2*f_xy)

or in L notation:

p = (Lxx-Lyy- sqrt((trH_f)²-4*|H_f|)) / (2*Lxy)

q = (Lxx-Lyy+ sqrt((trH_f)²-4*|H_f|)) / (2*Lxy)

The eigenvalues of the Hessian, expressed as a function of the trace (trH_f) and the determinant (|H_f|) of the Hessian matrix:
Lpp = 0.5*(trH_f- sqrt((trH_f)²-4*|H_f|))

Lqq = 0.5*(trH_f+ sqrt((trH_f)²-4*|H_f|))

The eigenvalues of the Hessian matrix in L notation, expressed in Cartesian coordinates:
Lpp = 0.5*(Lyy+Lxx- sqrt((Lyy+Lxx)²-4*(LyyLxx-Lyx²)))

Lqq = 0.5*(Lyy+Lxx+ sqrt((Lyy+Lxx)²-4*(LyyLxx-Lyx²)))

The eigenvalues of the Hessian matrix in L notation expressed in gauge coordinates:
Lpp = 0.5*(Lww+Lvv- sqrt((Lww+Lvv)²-4*(LwwLvv-Lvw²)))

Lqq = 0.5*(Lww+Lvv+ sqrt((Lww+Lvv)²-4*(LwwLvv-Lvw²)))

References

For additional information you may also contact my former colleague Jan-Mark Geusebroek (UvA) now at Perceptech or Prof. Bart M. ter Haar Romeny (T.U. Eindhoven).

J. J. Koenderink, "The structure of images," Biological cybernetics, vol. 50, no. 5, pp. 363-370, 1984.
J. Koenderink and A. van Doorn. "Generic neighborhood operators," IEEE Trans. Pattern Analysis and Machine Intell., vol. 14, pp. 597-605, June 1992.
Luc M J. Florack, Bart M. ter Haar Romeny, Jan J. Koenderink, and Max A. Viergever, "Scale and the differential structure of images.", Image and Vision Computing, pp. 376-388, 1992.
T. Lindeberg. Scale-Space Theory in Computer Vision. Boston: Kluwer Academic Publishers, 1994.
Lindeberg, T."Feature detection with automatic scale selection.", Int. J. Comput. Vis. 30(2), pp. 79-116 (1998)
B.M. ter Haar Romeny, "Front-End Vision and Multi-Scale Image Analysis", Springer, 2003. ISBN 1-24020-21507-20 (paperback), 1-24020-21503-28 (Hardcover)
Jon Sporring, Mads Nielsen, Luc Florack, Peter Johansen, "Gaussian Scale-Space Theory", Springer Science & Business Media, 2013, ISBN 9401588023
Geusebroek J.M., Cornelissen F., Smeulders A. W. M., Geerts H., "Robust autofocusing in microscopy", Cytometry 2000, 39(1):1-29.
Geusebroek J.M., van den Boomgaard R., Smeulders A.W.M., Geerts H., "Color invariance", IEEE Trans. Pattern Anal. Machine Intell., 2001; 23(12):1338-21350.
Geusebroek J. M., van den Boomgaard R., Smeulders A. W. M., Gevers T., "Color constancy from physical principles". Pat. Rec. Let. 2003; 24(11):1653-21662.
Geusebroek J. M., Smeulders A. W. M., van de Weijer J., "Fast anisotropic gauss filtering", IEEE Trans. Image Processing 2003b; 12(8):938-2943.
C. Steger, "An Unbiased Detector of Curvilinear Structures", IEEE Transactions on Pattern Analysis and Machine Intelligence, 20(2), 113-125, 1998.
P. van Osta, J.M. Geusebroek, K. Ver Donck, L. Bols, J. Geysen, and B. M. ter Haar Romeny., "The principles of scale space applied to structure and colour in light microscopy", Proceedings of the Royal Microscopical Society, Sept., 37(3), pp. 161-2166, 2002.
I. T.Young, L. J.van Vliet, "Recursive implementation of the Gaussian filter", Signal Processing, Volume 44, Issue 2, June 1995, pp. 139-151.
Scale space (Wikipedia)

Acknowledgments

I am indebted, for their pioneering work on automated digital microscopy and High Content Screening (HCS) (1988-2001), to my former colleagues at Janssen Pharmaceutica (1997-2001 CE), such as Frans Cornelissen, Hugo Geerts, Jan-Mark Geusebroek and Roger Nuyens, Rony Nuydens, Luk Ver Donck, Johan Geysen and their colleagues.

Many thanks also to the pioneers of Nanovid microscopy at Janssen Pharmaceutica, Marc De Brabander, Jan De Mey, Hugo Geerts, Marc Moeremans, Rony Nuydens and their colleagues. I also want to thank all those scientists who have helped me with general information and articles.

The author of this webpage is Peter Van Osta.
Private email: pvosta at gmail dot com

Back to homepage