davideberly / GeometricTools

A collection of source code for computing in the fields of mathematics, geometry, graphics, image analysis and physics.
Boost Software License 1.0
1.08k stars 202 forks source link

Does a continuous mapping from the matrices to the eigenvectors exist? #81

Closed melnikovsky closed 5 months ago

melnikovsky commented 6 months ago

It seems to me that the ridge detection algorithm may produce parasitic solutions when the sign of either P = Dot(U,DF) or Q = Dot(V,DF) changes without actually hitting a genuine root. This may happen because the eigenvectors are defined up to a factor and the vectors produced by the SymmetricEigensolver will jump at certain points. For example, the discontinuity will happen around the matrix: $\pmatrix{ 1 & 1 \\ 1 & 1}.$

Am I wrong?

davideberly commented 6 months ago

I do not know what you mean by "parasitic solutions" or by "genuine root" or by "eigenvectors are defined up to a factor." If I understand correctly, you appear to be questioning what happens when the eigenspace is 2-dimensional (the eigenvalues are equal).

The definition I provide in my ridges book and also in an early draft PDF of part of the book is the following where f(x,y) is the function of interest:

Df // gradient of f (vector of first-order partial derivatives)
D2f // Hessian of f (matrix of second-order partial derivatives)
a, b // eigenvalues of D2f with a < b
U, V // eigenvectors for a and b, respectively
D2f*U = a*U
D2f*V = b*V
P = Dot(U,Df)
Q = Dot(V,Df)
P = 0 and a < 0 and a < b // points (x,y) for which these are true are ridge points

In the book and PDF, I mention that a < b is required and that P and Q are (usually) discontinuous when a = b (the height field has an umbilic point when a = b).

melnikovsky commented 6 months ago

Sorry for the confusion, my concerns are not related to the degenerate case. Consider the following two Hessians at two neighbor points:

$$ \pmatrix{1 & 0.999 \\ 0.999 & 1} \quad\text{and}\quad \pmatrix{1.001 & 1 \\ 1 & 1} $$

Their eigenvalues are well separated ($\approx 0$ and $\approx 1$). The eigenvectors returned by SymmetricEigensolver2x2 are:

$$ \pmatrix{0.707107 \\ -0.707107}, \pmatrix{0.707107 \\ 0.707107} \quad\text{and}\quad \pmatrix{-0.70693 \\ 0.707284}, \pmatrix{-0.707284 \\ -0.70693} $$

Both P and V will change signs, but this does not mean they passed zero. Does it make sense?

davideberly commented 6 months ago

Thank you for the details.

Your first matrix has eigenvalues 1.999 and 0.001, so neither eigenvalue is negative. According to the ridge definition, the point at which that matrix was computed is not a ridge point, so there is no reason to test P = 0.

Your second matrix has eivenvalues approximately 2.0005 and 0.00499875, so neither eigenvalue is negative. Once again, the point at which the matrix was computed is not a ridge point.

The idea is that as you "walk across the ridge", the graph is "concave down". At the ridge point, the graph has a maximum in the direction you walk, and a < 0 indicates the graph has a maximum at that point in the direction you walk. If you were to reach that point and walk towards a peak, the graph does not have a maximum at that point in the tangential direction towards the peak.

davideberly commented 6 months ago

Ah yes, I believe I now understand what your concern is. I was focused more on the ridge definition that a < 0 and a < b. Similar ideas apply to defining a "valley" with b > 0 and b > a. Your example matrices involve valley point classification. Multiplying them by -1, now you have a ridge point classification.

Your focus is on the numerical eigensolver returning U0 at a point and U1 at a nearby point with U1 approximately -U0. Assuming real-number arithmetic, the sign should not matter. If P0 = Dot(U0,Df) = 0 at a point, then -P0 = Dot(-U0,Df) = 0. At the nearby point, P1 = Dot(U1,Df) = 0 and -P1 = Dot(-U1,Df) = 0. Using floating-point arithmetic, say P0 is a small positive number and -P1 is a small negative number. Now you have a sign change because of rounding errors, but it is not clear why you would be comparing nearby P-values when attempting to detect a ridge at a single point.

Because you mentioned SymmetricEigensolver, you are probably looking at the ridge sample application (which is applied to discrete images, not on a continuous domain). I will revisit the sample to understand whether there is a problem (implicitly using P-value sign changes?).

Thank you for posting an issue about this.

melnikovsky commented 6 months ago

Right, I was referring to ExtractRidgesConsole.cpp, is there some other place to look for the reference implementation? A "natural" approach to seek for the sign changes to locate the roots may fail when when the eigenvector changes from $\mathbf{U}$ to $-\mathbf{U}$. I'll probably submit a fix for that. And sorry for being unclear, my example really corresponds to a valley (course?) rather than to a ridge.

davideberly commented 6 months ago

That code was written originally about 30 years ago. I have not thought about the topic for many years.

Mathematically, the sign change problem can be eliminated by converting to a minimization problem: Minimize (Dot(U,Df))^2, which effectively transforms the problem so that a ridge curve for the height field becomes a 2D curve in the plane. The minimization approach makes it difficult to derive a robust algorithm.

If I remember correctly, I had tried using polynomial interpolation of images so that you have a piecewise-defined height field that is globally C^2 continuous. The best ridge detector was a mixture of symbolic arithmetic and floating-point arithmetic.

I used the term "valley", avoiding the word "course." The latter word has the connotation of "watersheds", but in the preface to the ridges book, I point out that watershed ridges are non-local, which is not desirable in an image processing setting. The German word "Thalweg" has also been used, but this is also a geological concept related to watersheds.

davideberly commented 6 months ago

I merged your changes with the master branch. Thank you for fixing the problems. I posted a version of the file with a new "Version" number and modified some code to match my formatting preferences.

davideberly commented 5 months ago

Closing this as it apparently addresses the concerns of the issue poster.

melnikovsky commented 5 months ago

Hi again, sorry for reopening this, is there a better place to post comments? Could you please take a look at the short note I wrote. I'll post a proof of concept implementation later. Your opinion would be gratefully appreciated.

davideberly commented 5 months ago

I have looked at your PDF. I will send you private email with my detailed feedback. I suspect you do not have a copy of the "Ridges" book. You or anyone else should not purchase one because the price is excessively large. I have sworn off book publishers.

The book has a discussion about equivalent ridge definitions, one of them effectively what you proposed. I can extract a section of the book as a PDF (which is allowed by copyright, even though I assigned the copyright to the publisher).

davideberly commented 5 months ago

Your "short note" is very well written! Closing this based on private email discussion.