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:
- track the location of extrema and saddle point as a function of $\sigma$
- locate the merging events between two saddle points / a saddle point and an extremum
- 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:
- 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
- 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.
