Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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"] &

Mathematica graphics

Indeed

um = Image[u] // MinDetect;
uM = Image[u] // MaxDetect;
ImageAdd[ImageAdd[iu, uM], um]

Mathematica graphics

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}}]]

Mathematica graphics

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}],
 Graphics[{AbsolutePointSize[6], Green, Point /@ saddles}]}]

Mathematica graphics

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

share|improve this question
    
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
1  
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
1  
@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
show 2 more comments

5 Answers

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

Start with a Gaussian random field

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

Build its gradient

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}];

Imgur

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}]]

Imgur

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]}]]

Imgur

share|improve this answer
add comment

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"];
imgA= ColorNegate@GradientFilter[imgG, 3] // ImageAdjust
imgB= ImageMultiply[MeanShiftFilter[imgG, 2, 0.3, MaxIterations -> 200], imgA]

image A and B

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

enter image description here

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

enter image description here

share|improve this answer
add comment

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

Mathematica graphics

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

share|improve this answer
add comment

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}]

Imgur

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

Imgur

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

Imgur

share|improve this answer
add comment

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"]

Gaussian random field

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

WatershedComponents[img] // Colorize

watershed components for random field

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

watershed components for random field with blurring

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

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

enter image description here

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

ImageMultiply[img, outline]

enter image description here

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.

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.