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
@chris en.wikipedia.org/wiki/Marching_cubes is 3D, based on en.wikipedia.org/wiki/Marching_squares mentioned there. Also Azim's answer there is interesting to read. PS: Under wikipedia matlab or java codes are also available. –  s.s.o Apr 8 at 18:49

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 '13 at 14:10
as in ImageMultiply[WatershedComponents[Blur[img, 20]] // Image, WatershedComponents[ColorNegate@ Blur[img, 20]] // Image] –  chris Jun 9 '13 at 14:17
@chris -- thanks, see updates –  bill s Jun 9 '13 at 15:12

Another possible path is to extend Vitaliy's function FindCrossings2D to 3D

In its current form, it is inefficient and somewhat buggy.

Identify all 3D extrema

u = GaussianRandomField[16, 3, Function[k, k^-1]] // Chop // GaussianFilter[#, 6] &;


fu = u // ListInterpolation[#, Method -> "Spline", InterpolationOrder -> 4] &;
fux = Function[{x, y, z}, D[fu[x, y, z], x] // Evaluate];
fuy = Function[{x, y, z}, D[fu[x, y, z], y] // Evaluate];
fuz = Function[{x, y, z}, D[fu[x, y, z], z] // Evaluate];


Identify the zero contours of the 3 components of the gradient

ContourPlot3D[{fux[x, y, z] == 0, fuy[x, y, z] == 0,
fuz[x, y, z] == 0}, {x, 1, 16}, {y, 1, 16}, {z, 1, 16}];


We can hack Vitaly's function to 3D

FindCrossings3D[{f_, g_, h_}, {x_, xmin_, xmax_}, {y_, ymin_, ymax_},
{z_, zmin_, zmax_}] := {x, y, z} /. (FindRoot[{f[x, y, z] == 0, g[x, y, z] == 0,
h[x, y, z] == 0},
{{x, #[[1]]}, {y, #[[2]]}, {z, #[[3]]}},PrecisionGoal->2] & /@
(ContourPlot3D[{f[x, y,z] == 0, g[x, y, z] == 0, h[x, y, z] == 0},
{x, xmin,xmax}, {y, ymin, ymax}, {z, zmin, zmax}, PlotPoints -> 8][[1, 1,;;;;10]]))


so that all the extrema are given by

 pts = FindCrossings3D[{fux, fuy, fuz}, {x, 1, 16}, {y, 1, 16}, {z, 1, 16}]//Quiet;


Note two minor modifications: I Plot fewer points than the default (PlotPoints -> 8) and only start from one point out of ten (;;;;10).

Then I can do

 Show[ListContourPlot3D[u], Graphics3D[{AbsolutePointSize[6], Red, Point /@ pts}]]


The remaining issue is the FindRoot does not converge well to unique extrema but rather to a bundle of them.

Sort out extrema of various type

In order to pick only extrema of a given type (maxima, minima, wall and tube saddles) we need to check their signature

fuxx = Function[{x, y, z}, D[fu[x, y, z], x, x] // Evaluate];
fuyy = Function[{x, y, z}, D[fu[x, y, z], y, y] // Evaluate];
fuzz = Function[{x, y, z}, D[fu[x, y, z], z, z] // Evaluate];
fuxy = Function[{x, y, z}, D[fu[x, y, z], x, y] // Evaluate];
fuxz = Function[{x, y, z}, D[fu[x, y, z], x, z] // Evaluate];
fuyz = Function[{x, y, z}, D[fu[x, y, z], y, z] // Evaluate];


Now we can compute the determinant, trace and minor of the hessian matrix.

pxx=  Map[fuxx @@ # &, pts]; pyy=  Map[fuxx @@ # &, pts]; pzz=  Map[fuxx @@ # &, pts];
pxy=  Map[fuxy @@ # &, pts];  pxz=  Map[fuxz @@ # &, pts]; pyz=  Map[fuyz @@ # &, pts];
dets =  Map[Det,{{pxx,pxy,pxz},{pxy,pyy,pyz},{pxz,pyz,pzz}}// Transpose[#, {3, 2, 1}] &];
trs = pxx+pyy+pzz;
minors= pxx pyy + pxx pzz + pyy pzz - pxy^2 - pxz^2 - pyz^2;


Selecting the maxima (-,-,-) and minima (+,+,+)

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

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


and the saddles as critical points which have a signature (-,-,+) and (-,+,+)

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


and

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


For instance the green points are minima:

 Show[ListContourPlot3D[u],
Graphics3D[{AbsolutePointSize[6], Red, Point /@ Union[pts]}],
Graphics3D[{AbsolutePointSize[6], Green, Point /@ Union[min]}]]


-

Let me start and answer my own question (with what I have so far) since it might trigger some interest before the bounty expires. The motivation is to find rapidly 3D maxima of the field.

Posible idea

The idea is to define 3D maxima, as the intersection of the union of 2D maxima (using 'MaxDetect') sliced in two different directions.

Starting with a 3D cube:

 u = GaussianRandomField[128, 3, Function[k, k^-1]] // Chop // GaussianFilter[#, 12] &;


I then build the set of 2D maxima across two directions:

 m1 = Map[ImageData@MaxDetect[Image[#]] &, u];
m2 = Transpose@Map[ImageData@MaxDetect[Image[#]] &, Transpose[u]]];


It looks like this

 GraphicsRow[{
m1[[2 ;; -1, 2 ;; -1, 2 ;; -1]] // Image3D // ImageAdjust,
m2[[2 ;; -1, 2 ;; -1, 2 ;; -1]] // Image3D // ImageAdjust}]


Note that I remove the edges for which MaxDetect seems to do something strange.

The 3D maxima are then approximately given by

(m1*m2)[[2 ;; -1, 2 ;; -1, 2 ;; -1]] // ListContourPlot3D


Note that it should also work in 4D!

This technique could clearly be applied to detect the minima of the field as well.

In fact this allows us to check that maxima and minima are indeed on either sides of iso-contours.

u = GaussianRandomField[32, 3, Function[k, k^-1]] // Chop // GaussianFilter[#, 6] &;
field = u // ListContourPlot3D[#, Mesh -> None] &;
m1 = Map[ImageData@MaxDetect[Image[#]] &, u];
m2 = Transpose@Map[ImageData@MaxDetect[Image[#]] &, Transpose[u]];
max = (m1*m2)[[2 ;; -1, 2 ;; -1, 2 ;; -1]] // ListContourPlot3D;
m11 = Map[ImageData@MaxDetect[Image[#]] &, -u];
m22 = Transpose@Map[ImageData@MaxDetect[Image[#]] &, Transpose[-u]];
min = (m11*m22)[[2 ;; -1, 2 ;; -1, 2 ;; -1]] // ListContourPlot3D;
{Show[field, min], Show[field, max]} // GraphicsRow


-

Here is an other image processing alternative. The main idea is as a pre-processing step using Gradient Filter. (I think most of them may work in 3D as well.)

img = Import["http://i.stack.imgur.com/g7TFl.png"];
imgG = ColorConvert[img, "Grayscale"];
imgB= ImageMultiply[MeanShiftFilter[imgG, 2, 0.3, MaxIterations -> 200], imgA]


imgC = ImageAdd[img, Thinning@Binarize[imgB, .70]]


saddle= ImageCorners[Thinning@Binarize[imgB, .70], MaxFeatures -> 150];


-

A third attempt would be to try the WatershedComponent Approach presented by bill s above, as in

u = GaussianRandomField[128, 3, Function[k, k^-1]] // Chop //GaussianFilter[#, 15] &;


so that the watershed is applied on each slice:

WS = Map[WatershedComponents[Image[#],  Method -> {"MinimumSaliency", .3}] &, u, 1];

Colorize[#] & /@ WS[[Range[16]]]


What remains to be sorted out it how to identify ids of each slice.

This might be possible via a transposition:

 WS2 = Map[WatershedComponents[Image[#],Method -> {"MinimumSaliency", .3}] &, Transpose[u], 1];

-