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

UPDATE 2014

WRI has addressed this in Mathematica 10

2D case

In 2D, the saddle points can be identified as the intersections the two watersheds for the field and minus the field.

Using this post, we can generate a 2D Gaussian random field as

 u = GaussianRandomField[n = 256, 2, Function[k, 1/k Exp[-1/25 k^2]]] //Chop;
 u /= Max[Flatten@u];

Now the intersection of the two watersheds are simply given by

ttp = WatershedComponents@Image[u] /. x_?NumberQ :> If[x == 0, 1, 0];
ttm = WatershedComponents@Image[-u] /. x_?NumberQ :> If[x == 0, -1, 0];
tts = -ttp*ttm;

We can check that we identify correctly their intersection via the HighlightImage function:

HighlightImage[ImageAdd[ImageAdjust@Image@u, ImageAdd[Image@ ttp, Image@ttm]] // 
ImageAdjust, Image@tts,Method -> {"DiskMarkers", 2},"HighlightColor" -> Purple]

Mathematica graphics

The minima and maxima are found using MinDetect and MaxDetect

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

This can be shown to work as follows:

{HighlightImage[ImageAdjust@Image@u, um, Method -> {"DiskMarkers", 2.5}],
HighlightImage[ImageAdjust@Image@u, uM, Method -> {"DiskMarkers", 2.5},
"HighlightColor" -> Purple]} // Row

Mathematica graphics

Note that as expected, they sit on the critical lines:

u = GaussianRandomField[n = 64*16, 2, Function[k, k Exp[-1/15 k^2]]] //Chop;
tt = Map[If[# == 0, 1, 0] &, u // Image // WatershedComponents, {2}] //
SparseArray; im = tt // Image;
tt = Map[If[# == 0, 1, 0] &, -u // Image // 
WatershedComponents, {2}] // SparseArray; im2 = tt // Image;
HighlightImage[HighlightImage[HighlightImage[HighlightImage[
ReliefImage[u, ColorFunction -> ColorData["ThermometerColors"]], im],
 u // Image // MaxDetect, Method -> {"DiskMarkers", 2}, "HighlightColor" -> Red], 
 im2, "HighlightColor" -> Blue],
 -u // Image // MaxDetect, Method -> {"DiskMarkers", 2}, "HighlightColor" -> Purple]

Mathematica graphics

The branching points of the two sets of critical lines can identified using MorphologicalTransform on the watersheds:

{ImageAdd[
HighlightImage[ttp // Image, 
MorphologicalTransform[ttp, "SkeletonBranchPoints"] // Image,
Method -> {"DiskMarkers", 1}, "HighlightColor" -> Red],
ImageAdjust@Image@(u)],
ImageAdd[
HighlightImage[-ttm // Image, 
MorphologicalTransform[-ttm, "SkeletonBranchPoints"] // Image, 
Method -> {"DiskMarkers", 1}, "HighlightColor" -> Blue], 
ImageAdjust@Image@(-u)]} // Row

Mathematica graphics

3D case

In version 10 Mathematica MaxDetect and MinDetect are now working in 3D as follows

u = GaussianRandomField[n = 128, 3, Function[k, Exp[- k^2]]] // Chop;    
u // Image3D // ImageAdjust

Mathematica graphics

{MaxDetect[Image3D[u]], MinDetect[Image3D[u]]} // Row

Mathematica graphics

Let us check the projection on a more rough 3D field

u1 = GaussianRandomField[n = 32*4, 3, Function[k, Exp[-1/15 k^2]]] // Chop;
pl1 = HighlightImage @@ {Plus @@ u1 // Image // ImageAdjust, 
Plus @@ ImageData[MaxDetect[Image3D[u1]]] // Image} 

Mathematica graphics

We can also access directly the coordinates of the extrema

 u = GaussianRandomField[n = 32, 3, Function[k, Exp[-k^2]]] // Chop;
 MaxDetect[u] // SparseArray // ArrayRules

(*
{{1,1,2}->1,{1,1,20}->1,{1,1,32}->1,{1,12,32}->1,{1,23,4}->1,{1,27,32}->1,{1,30,20}->1,{19,11,1}->1,{19,11,32}->1,{20,32,13}->1,{22,32,1}->1,{22,32,30}->1,{27,27,32}->1,{29,22,3}->1,{,,_}->0} *)

With the new functionality Watershed of mathematica 10.0.1 It is now possible to get directly the 1D manifold connecting saddle points together as follows:

Now the corresponding watersheds obeys

ttp = WatershedComponents@Image3D[u] /. x_?NumberQ :> If[x == 0, 1, 0];
ImageAdjust@Image3D@tts    

Mathematica graphics

ttm = WatershedComponents@Image3D[-u] /. x_?NumberQ :> If[x == 0, 1, 0];

Mathematica graphics

tts = ttp*ttm;
ImageAdjust@Image3D@tts    

Mathematica graphics

Following this great answer we can define a function to trace the skeleton as well

intersections[ws_] := Module[{unique, dims},
dims = Dimensions[ws];
Reap[
Do[
  If[ws[[i, j, k]] == 0,
    unique = (Flatten[
      ws[[i - 1 ;; i + 1, j - 1 ;; j + 1, k - 1 ;; k + 1]]] // 
      Union)~Drop~1;
    If[Length[unique] >= 3, Sow[{i, j, k, unique}]]
  ], 
  {i, 2, dims[[1]] - 1}, 
  {j, 2, dims[[2]] - 1}, 
  {k, 2, dims[[3]] - 1}
  ]
 ][[2, 1]]
]

Then the skeleton of a cube is simply given by

skl[cube_] := Module[{list = intersections[cube]},
         SparseArray[#[[1 ;; 3]] -> 1 & /@ list] ]

So that

u = GaussianRandomField[n = 32*6, 3, Function[k, 1/k Exp[-1/2 k^2]]] //Chop;
dat = WatershedComponents@Image3D[u];
dat2 = WatershedComponents@Image3D[-u];

produces the set of critical lines connecting peaks and minima to saddles:

Row[Image3D[Normal[skl[#]]] & /@ {dat, dat2}] // Rasterize

Mathematica graphics

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

3 Answers 3

up vote 11 down vote accepted
+500

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

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

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

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.