Differential geometry applies techniques from multivariable calculus and linear algebra to study features of curves and surfaces of dimension 1, 2, and higher.

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 and astronomy.

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. 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.

A parabole y=x

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^{2}y/dx^{2} = 2^{nd} order derivative
(the curvature of the function, or the rate of change of the slope)

f_{xxx}=d^{3}y/dx^{3} = 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**is negative._{x}

- At each point
**x**at which the function**f(x)**is**increasing**the slope of the tangent is positive -- that is,**f**is positive._{x} - 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**, then_{xx}> 0**f(x)**is at a**minimum**. - If in addition
**f**, then_{xx}< 0**f(x)**is at a**maximum**.

- If in addition

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.

The principles of derivatives 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 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 principle can be extended to 3D images, but this is beyond the current scope of this webpage.

The Gaussian kernel is the unique kernel for differentiation for an
uncommitted front-end (eg. human eye).

Differentiating an image is done by convolving the image
in the spatial domain with a Gaussian kernel with a standard deviation sigma.

The sigma parameter of the Gaussian kernel determines the size of the
features that are resolved, ie. it is a measure for the aperture of the filter.

Zero order derivative of a Gaussian kernel in 1D, sigma is 10.0.

First order derivative of a Gaussian kernel in 1D, sigma is 10.0.

Second order derivative of a Gaussian kernel in 1D, sigma is 10.0.

Third order derivative of a Gaussian kernel in 1D, sigma is 10.0

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 images.

The original 256x256 pixels image from which the following differentiated images are derived.

Zero order Gaussian derivative in x and y direction, just blurs the image. Sigma (scale) is 2.

Original resolution is 256^{2}

First order Gaussian derivative in x direction. Sigma (scale) is 2.

Second order Gaussian derivative in x direction. Sigma (scale) is 2.

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 = || L|| = Lw

Gradient^{2} = L_{i}L_{i} = L_{x}^{2} + L_{y}^{2}

Laplacian = Trace of Hessian matrix = trH_{f} = L_{ii} = L_{xx} + L_{yy}

Determinant of Hessian matrix = |H_{f}| = L_{xx}L_{yy}-L_{xy}^{2}

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^{2}) = Sqrt(L_{ij}L_{ji}) =
Sqrt(L_{xx}^{2} + 2L_{xy}^{2} + L_{yy}^{2})

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}^{2} + L_{y}^{2})
*K*_{p2} = L_{yy}/Sqrt(L_{x}^{2} + L_{y}^{2})

Inflexion strength = LvvLw^{5}-3LvvLvwLw^{4}

Umbilicity = U = (trH)^{2}/TrH^{2} - 1 = (L_{xx} + L_{yy})^{2} / D^{2} =
L_{ll}^{2} / L_{ij}L_{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)^{2}(2trH^{2}-(trH)^{2}))/(2trH^{2})^{2} =
(L_{xx}^{2} + L_{yy}^{2})/D^{4}

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

All resultants and discriminants are Euclidian and affine invariant

R_{12} = Corner strength = L_{x}^{2}L_{yy}-2L_{y}L_{x}L_{xy}+
L_{y}^{2}L_{xx} = LvvLw^{2}

D_{1} = 1

D_{2} = det*H* = -L_{xy}^{2}+L_{yy} + L_{xx} or
L_{xx}+L_{yy}-L_{xy}^{2}

Lx = 0

Ly = 0

LxxLyy-Lxy^{2} >= 0

|Lxx| > threshold

threshold determines salient features

Lx = 0

Ly = 0

LxxLyy-Lxy^{2} >= 0

Lxx < -threshold

threshold determines salient features

Lx = 0

Ly = 0

LxxLyy-Lxy^{2} >= 0

Lxx > threshold

threshold determines salient features

Lx = 0

Ly = 0

LxxLyy-Lxy^{2} < 0

Lx = 0

Ly = 0

LxxLyy-Lxy^{2} = 0

affine corner detection: LwwLw^{2}

with high curvature: LwwLw^{2} > threshold

Bright ridges: Lqq < -4.0 * threshold / sigma^{2}

Dark ridges: Lpp > 4.0 * threshold / sigma^{2}

threshold determines salient features

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 that are of importance as it has geometrical meaning
indepedent 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 H. Weyl (1885-1955) in 1918 in Sitzungsber. d. Preuss. Akad. d. Wissensch. 30 May 475 (OED2).

A special notation is used in image processing for this gauge system: the
*w,v*-coordinates. *w* is the coordinate in the gradient direction 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.

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 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_{xx}f_{yy}-f_{xy}^{2}=
LxxLyy-Lxy^{2}=LwwLvv-Lvw^{2}

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), named for Otto Hesse, who had used the term functional determinants.

Evaluation of differential expressions:

L = the 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

A way to approximate an image is to perform a Taylor expansion of the local
image structure. Locally, i.e. in a small neighborhoud 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.

Consider a 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:

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

so that

Substitute this result into () to give

A Taylor series of a function in two variables f(x,y) is given by

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.

A necessary condition that a function f assumes a relative maximum or a
relative minimum at (x_{0}, y_{0}) is that:

f_{x} = f_{y} = 0 at (x_{0},y_{0})

If this condition is true then in summary the following criteria can be used.

All subsequent images were analyzed with sigma = 5.0 .

f_{x} = 0, f_{y} = 0 and f_{xx}f_{yy} > f^{2}_{xy}

or i.o.w.

f_{x} = 0, f_{y} = 0 and f_{xx}f_{yy} - f^{2}_{xy} > 0

or i.o.w.

Lx = 0, Ly = 0 and LxxLyy - Lxy^{2} > 0

f_{x} = 0, f_{y} = 0 and f_{xx}f_{yy} < f^{2}_{xy}

or i.o.w.

f_{x} = 0, f_{y} = 0 and f_{xx}f_{yy} - f^{2}_{xy} < 0

or i.o.w.

Lx = 0, Ly = 0 and LxxLyy - Lxy^{2} < 0

f_{x} = 0, f_{y} = 0 and f_{xx}f_{yy} = f^{2}_{xy}

or i.o.w.

f_{x} = 0, f_{y} = 0 and f_{xx}f_{yy} - f^{2}_{xy} = 0

or i.o.w.

Lx = 0, Ly = 0 and LxxLyy - Lxy^{2} = 0

f_{xx} < 0 and f_{xx}f_{yy} > f^{2}_{xy}

or i.o.w.

f_{xx} < 0 and f_{xx}f_{yy} - f^{2}_{xy} > 0

or i.o.w.

Lxx < 0 and LxxLyy - Lxy^{2} > 0

or i.o.w.

Lww < 0 and LvvLww - Lvw^{2} > 0

f_{xx} > 0, f_{xx}f_{yy} > f^{2}_{xy}

or i.o.w.

f_{xx} > 0, f_{xx}f_{yy} - f^{2}_{xy} > 0

or i.o.w.

Lxx > 0 and LxxLyy - Lxy^{2} > 0

or i.o.w.

Lww > 0, LvvLww - Lvw^{2} > 0

f_{xx}f_{yy} < f^{2}_{xy}

or i.o.w.

f_{xx}f_{yy} - f^{2}_{xy} < 0

or i.o.w.

LxxLyy - Lxy^{2} < 0

or i.o.w.

LvvLww - Lvw^{2} < 0

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 the neurites of a neuron with differential geometry as an example of a line detector.

The equation used for detecting a dark ridge, with sigma 2.0 and threshold 1.0:

Lpp > 4.0 * threshold / s^{2}

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.

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})^{2}-4*|H_{f}|)) / (2*f_{xy})

q = (f_{xx}-f_{yy}+
sqrt((trH_{f})^{2}-4*|H_{f}|)) / (2*f_{xy})

p = (Lxx-Lyy-
sqrt((trH_{f})^{2}-4*|H_{f}|)) / (2*Lxy)

q = (Lxx-Lyy+
sqrt((trH_{f})^{2}-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})^{2}-4*|H_{f}|))

Lqq = 0.5*(trH_{f}+
sqrt((trH_{f})^{2}-4*|H_{f}|))

The eigenvalues of the Hessian matrix in L notation, expressed in Cartesian coordinates:

Lpp = 0.5*(Lyy+Lxx-
sqrt((Lyy+Lxx)^{2}-4*(LyyLxx-Lyx^{2})))

Lqq = 0.5*(Lyy+Lxx+
sqrt((Lyy+Lxx)^{2}-4*(LyyLxx-Lyx^{2})))

The eigenvalues of the Hessian matrix in L notation expressed in gauge coordinates:

Lpp = 0.5*(Lww+Lvv-
sqrt((Lww+Lvv)^{2}-4*(LwwLvv-Lvw^{2})))

Lqq = 0.5*(Lww+Lvv+
sqrt((Lww+Lvv)^{2}-4*(LwwLvv-Lvw^{2})))

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

- B.M. ter Haar Romeny, "Front-End Vision and Multi-2Scale Image Analysis", Springer, 2003. ISBN 1-24020-21507-20 (paperback), 1-24020-21503-28 (Hardcover)
- 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.
- 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.

- Color Differential Geometry or the Spatial Color Model
- Nyquist Sampling in Digital Microscopy
- Image Analysis gives a clear view in research
- The Basics of Microscopy
- Digital Cameras and Microscopy
- Common Imaging Artefacts
- Application of linear scale space and the spatial color model in light microscopy
- Automated Tiled Multi-mode Image Acquisition and Processing Applied to Pharmaceutical Research
- The M
^{5}framework for exploring the cytome

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), such as Frans Cornelissen, Hugo Geerts, Jan-Mark Geusebroek, Roger Nuyens, Rony Nuydens, Luk Ver Donck 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