# Identifying critical points of 2/3D image/cubes

Upshot

I am interested in identifying critical points of a 3D field/cubes (maxima, minima, tube-like and wall-like saddle points) and 2D field/image (maxima, minima, saddle points). I.e. the generalisation of How to find all the local minima/maxima in a range.

Attempt(s)

The mathematica functions MinDetect and MaxDetect address partially my request in 2D, as illustrated (using the GaussianRandomField function defined here)

u = GaussianRandomField[192] // Chop // GaussianFilter[#, 8] &;
iu = Image[u] // ImageAdjust //Colorize[#, ColorFunction -> "ThermometerColors"] &


Indeed

um = Image[u] // MinDetect;
uM = Image[u] // MaxDetect;


Shows that maxima and minima have been detected. On the other hand, saddle points are still to be found.

A related query was made there, which is a possible starting point since critical points are the simultaneous zeros of the (2 or 3 components of the) gradient. Let us try and follow this method:

fu = u // ListInterpolation;
fux = Function[{x, y}, D[fu[x, y], x] // Evaluate];
fuy = Function[{x, y}, D[fu[x, y], y] // Evaluate];

fuxx = Function[{x, y}, D[fu[x, y], {x, 2}] // Evaluate];
fuyy = Function[{x, y}, D[fu[x, y], {y, 2}] // Evaluate];
fuxy = Function[{x, y}, D[fu[x, y], x, y] // Evaluate];

pts = FindAllCrossings2D[{fux[x, y], fuy[x, y]}, {x, 1, Dimensions[u][[1]]},
{y, 1, Dimensions[u][[2]]},
Method -> {"Newton", "StepControl" -> "LineSearch"},
PlotPoints -> 256, WorkingPrecision -> 20] // Chop;

dets = Map[fuxx @@ # &, pts] Map[fuyy @@ # &, pts] -
Map[fuxy @@ # &, pts]^2;
trs = Map[fuxx @@ # &, pts] + Map[fuyy @@ # &, pts];


Selecting saddles as critical points which have a Hessian determinant negative.

w = Map[If[# < 0, 1, 0] &, dets]*Map[If[# > 0, 1, 0] &, trs];
saddles = Most /@ Select[Transpose[Join[Transpose[pts], {w}]], #[[3]] == 1 &];


and the maxima and minima

w = Map[If[# > 0, 1, 0] &, dets]*Map[If[# < 0, 1, 0] &, trs];
max = Most /@
Select[Transpose[Join[Transpose[pts], {w}]], #[[3]] == 1 &];
w = Map[If[# > 0, 1, 0] &, dets]*Map[If[# > 0, 1, 0] &, trs];
min = Most /@
Select[Transpose[Join[Transpose[pts], {w}]], #[[3]] == 1 &];


Check that they are indeed at the intersection of the zero gradient contours

ContourPlot[{fux[x, y], fuy[x, y]}, {x, 1, 256}, {y, 1, 256},
Contours -> {0}, ContourShading -> False,
ContourStyle -> {Red, AbsoluteThickness[1.5]},
Epilog -> Join[{{AbsolutePointSize[6], Green, Point /@ saddles},
{AbsolutePointSize[6], Red, Point /@ max},
{AbsolutePointSize[6], Blue, Point /@ min}}]]


and that they correspond to critical points of the underlying field

 Show[{ContourPlot[fu[x, y], {x, 1, 256}, {y, 1, 256}, Contours -> 35,
ColorFunction -> "ThermometerColors"],
Graphics[{AbsolutePointSize[6], Purple, Point /@ max}],
Graphics[{AbsolutePointSize[6], Gold, Point /@ min}],


Questions

i) would you know of a more efficient way of finding the saddle points in 2D?

ii) could this algorithm be generalized to 3D without significant loss of efficiency? (involves possibly writing a 3D version of FindAllCrossings2D)

iii) how robust can it be in terms of the smoothness of the field?

Here I am after an algorithm which would be efficient (I would like to identify thousands of critical points).

-
Note that the above method missed at least one extremum @ coordinates $\approx$ {0,55} as it does not properly account for the periodicity of the field. –  chris Aug 27 '12 at 21:38
Do you mean GaussianRandomField in the question (which was slow) or the improved addendum in your answer? I think you mean the latter, in which case you can link to it directly... –  rm -rf Aug 27 '12 at 21:45
@R.M I guess for this example speed does not matter. I don't know how to link to a subsection of a question. –  chris Aug 27 '12 at 21:50
In using ListInterpolate you introduced certain smoothness to your data which comes with the interpolation scheme. While you used this to get a function which you can derivate this maybe covers some critical points. Why don't you stay on your discrete data and check the huge amount of literature about extreme/sattle point finding algorithms? As an example maybe see this paper here. –  halirutan Sep 3 '12 at 3:34

If you are willing to accept an approximate answer, then this can be attacked in a simple way using the image processing functionality. Methods like the Watershed algorithm are intended to partition images into pieces, and can be used to locate (approximate) boundaries between the basins of critical points. For example, start with the image above

img = Import["http://i.stack.imgur.com/g7TFl.png"]


The segmentation algorithm separates the image into pieces. The Colorize command displays the segments in false color.

WatershedComponents[img] // Colorize


If too many regions are detected (the image is reasonably noisy) then it would make sense to filter before applying the watershed.

WatershedComponents[Blur[img, 10]] // Colorize


In order to look at the regions (without the fake color)

outline = Binarize[WatershedComponents[Blur[img, 10]] // Colorize, 0.1]


Superimposing this over the original image shows how the segmentation works:

ImageMultiply[img, outline]


Advantages of this kind of method are that it is quick, generalizes nicely, and is flexible. Disadvantages are that it is not really doing exactly the same thing as requested. To be clear: it is not finding saddle points per se, it is finding places where the image segments nicely. These are similar, but not identical goals.

If more control is desired over the process, it is possible to combine this approach with the analysis. After making the list of maxima and minima, you can enter these as optional parameters into the Watershed algorithm: it will then take each of these "markers" as a point in a separate region and segment accordingly.

-
In 2D you can define the saddle points as the intersection of watersheds of the field and minus the field. –  chris Jun 9 at 14:10
as in ImageMultiply[WatershedComponents[Blur[img, 20]] // Image, WatershedComponents[ColorNegate@ Blur[img, 20]] // Image] –  chris Jun 9 at 14:17
@chris -- thanks, see updates –  bill s Jun 9 at 15:12