Identifying critical points of 2/3D image/cubes

Question

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;
ImageAdd[ImageAdd[iu, uM], um]

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

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