Extremas and saddle point of 3D field at different scales
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 maximas, minimas, saddle point 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 extremas 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 $$ Extremas 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. 281
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 extremas 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 extremas 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. 835