Extrema and saddle point of 3D field at different scales

Question

I have a scalar 3D field $f(x, y, z)$ with $x,y,z$ on a regular grid. I would like to know the location of the maxima, minima, saddle points and their relation as a function of a smoothing scale.

For that, I'm convolving my field $f$ with a filtering function that is a Gaussian with a variance $\sigma$. I vary $\sigma$ from 0 to the size of my grid. I want to be able to:

1. track the location of extrema and saddle point as a function of $\sigma$
2. locate the merging events between two saddle points / a saddle point and an extremum
3. be able to tell which saddle point / extremum merged into which

An extremum / saddle point is defined as

$$ \nabla f = 0 $$

Extrema are distinguished using the Hessian:

$$ H_{i,j} = \nabla_i\nabla_jf = \mathrm{diag}(e_1, e_2, e_3)$$

where $e_1, e_2, e_3$ are the eigenvalues of the Hessian. If all of them are negative (resp. positive), the point is a maximum (resp. minimum). If one and only one is positive (resp. negative), the point is a saddle point which is a maximum (resp. minimum) in 2 dimensions.

Ridges are defined as curves on which

$$ \nabla f = \mathrm{diag}(0,e_1,e_2)$$ where $e_1, e_2 \neq 0$.

For now, I am finding the location of the extrema and saddle point by doing a quadratic interpolation $P(x,y,z)$ on each point of my grid and by computing its first and second derivatives. However I have two issues:

1. each time I smooth at a different scale, I don't use the information about the location of the extrema and saddle point at the previous scale
2. I don't know how to robustly find which extrema is linked to which saddle point, and hence, I can't tell robustly which point is merging into which.

Do you have any advice about this problem?

NB: it has to be fast enough to be able to run in less than a month on a 100x100x100 grid, at ~100 smoothing scale for 100,000 grids.

• Raphael: Welcome to Computer Science! The title you have chosen is not well suited to representing your question. Please take some time to improve it; we have collected some advice here. Thank you!
• Raphael: Do you have any ideas for helpful tags?
• Nick Alger: @Raphael Perhaps digital-morse-theory ? Probably this problem is equivalent to constructing Reeb graph for the 4D function $g(x,y,z,\sigma) := f_\sigma(x,y,z)$. A few years ago I saw a talk where someone was doing something like this computationally, but I don't remember the details.
• cphyc: Thanks for the tip, I'm trying to create an algorithm based on that. I'll update the feed once I've got some results.
• cphyc: For the record, there is the Disperse project (github.com/thierry-sousbie/DisPerSE) as well as diamorse (github.com/AppliedMathematicsANU/diamorse). Unfortunately none of them are 4D enable yet.