Ellipsis fitting: Check the condition number before computiing eigen vectors #3103

hmaarrfk · 2018-05-24T20:49:03Z

Under certain conditions, the result of a matrix we have to invert becomes ill conditioneed. The 32 bit solver and 64 bit solver converge on different answers.
We first check if that matrix is ill conditioned before we move forward with the computation.

See #3091

Checklist

[It's fine to submit PRs which are a work in progress! But before they are merged, all PRs should provide:]

Clean style in the spirit of PEP8
Docstrings for all functions
Gallery example in ./doc/examples (new features only)
Unit tests

[For detailed information on these and other aspects see scikit-image contribution guidelines]

References

[If this is a bug-fix or enhancement, it closes issue # ]
[If this is a new feature, it implements the following paper: ]

For reviewers

(Don't remove the checklist below.)

Check that the PR title is short, concise, and will make sense 1 year
later.
Check that new functions are imported in corresponding __init__.py.
Check that new features, API changes, and deprecations are mentioned in
doc/release/release_dev.rst.

pep8speaks · 2018-05-24T20:49:05Z

Hello @hmaarrfk! Thanks for updating the PR.

There are no PEP8 issues in the file skimage/measure/fit.py !
In the file skimage/measure/tests/test_fit.py, following are the PEP8 issues :

Line 238:80: E501 line too long (83 > 79 characters)
Line 241:80: E501 line too long (82 > 79 characters)
Line 242:80: E501 line too long (81 > 79 characters)
Line 243:80: E501 line too long (80 > 79 characters)

Comment last updated on September 20, 2018 at 03:00 Hours UTC

codecov-io · 2018-05-24T21:49:58Z

Codecov Report

Merging #3103 into master will decrease coverage by 0.09%.
The diff coverage is 100%.

@@           Coverage Diff            @@
##           master   #3103     +/-   ##
========================================
- Coverage   86.79%   86.7%   -0.1%     
========================================
  Files         341     341             
  Lines       27444   27429     -15     
========================================
- Hits        23820   23782     -38     
- Misses       3624    3647     +23

Impacted Files	Coverage Δ
skimage/measure/tests/test_fit.py	`100% <100%> (ø)`	⬆️
skimage/viewer/qt.py	`25.71% <0%> (-62.86%)`	⬇️
skimage/viewer/__init__.py	`80% <0%> (-20%)`	⬇️
skimage/segmentation/random_walker_segmentation.py	`92.01% <0%> (-1.88%)`	⬇️
skimage/draw/_random_shapes.py	`95.78% <0%> (-1.06%)`	⬇️
skimage/feature/tests/test_register_translation.py	`100% <0%> (ø)`	⬆️
skimage/feature/register_translation.py	`100% <0%> (+6.41%)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 6af7e2e...6ce1e2e. Read the comment docs.

stefanv · 2018-05-24T22:51:01Z

skimage/measure/fit.py

@@ -434,6 +434,11 @@ def estimate(self, data):
        except np.linalg.LinAlgError:  # LinAlgError: Singular matrix
            return False

+        # For some reason the above doesn't catch the error on linux 32 bit
+        # https://stackoverflow.com/questions/13249108/efficient-pythonic-check-for-singular-matrix
+        if np.linalg.cond(M) > 1 / np.finfo(np.double).eps:


This check can be done before the inversion is attempted? It simplifies the logic a bit.

I'm not sure. in fact, when I computed cond(M) in the debugger with import pdb; pdb.set_trace(), the condition number was already Inf. I don't really want to cause other overflow errors or "old numpy specific errors"

From the docs

eps | (float) The smallest representable positive number such that 1.0 + eps != 1.0. Type of eps is an appropriate floating point type.

They don't seem to have a eps^-1. Thought I could do
eps**-1 if you like.

stefanv · 2018-05-24T22:51:33Z

It seems like the test suite already fails, but it may be worth inserting an explicit test for this, potentially one that also fails on 64-bit.

matthew-brett · 2018-05-24T23:12:26Z

Should this be a bug report to numpy too?

stefanv · 2018-05-24T23:37:10Z

@matthew-brett I've raised this with the NumPy folks; they think it's not NumPy's responsibility.

Overall, it is a bad idea to use inverse if you can use solve. Could we rewrite the solution here to do solves instead?

hmaarrfk · 2018-05-24T23:43:25Z

Just for you two, I actually stepped through to figure out what was happening.

S1 and C1 are well conditioned. It is the result of the matrix multiplications and addition that is illconditionned. Typically this is caught a little lower

        # M*|a b c >=l|a b c >. Find eigenvalues and eigenvectors
        # from this equation [eqn. 28]
        eig_vals, eig_vecs = np.linalg.eig(M)

        # eigenvector must meet constraint 4ac - b^2 to be valid.
        cond = 4 * np.multiply(eig_vecs[0, :], eig_vecs[2, :]) \
               - np.power(eig_vecs[1, :], 2)
        a1 = eig_vecs[:, (cond > 0)]
        # seeks for empty matrix
        if 0 in a1.shape or len(a1.ravel()) != 3:
            return False

On 64 bit:

(Pdb) cond
array([-0.28717554, -0.71271829, -0.00108459])

On 32 bit

(Pdb) cond
array([-1.38708836,  0.84396406, -0.03172746])

It just happened that on 64 bit the error didn't show up. But I don't think that the eig algorithm will complain if you have a ill-conditionned matrix.

Therefore, the tests in both 64 bit and 32 bit are caught by my added logic.

hmaarrfk · 2018-05-24T23:56:22Z

Technically, since the matrix is so small, I think we can simply check that if one eig_vals is small and also return False.

Anyway, I think checking the condition number is more "mathematical"

hmaarrfk · 2018-05-25T01:42:54Z

@stefanv. I probably agree that doing something like finding tge best eigen value is probably verbose.

From ‘numpy.linalg.solve’

a must be square and of full-rank

We would need to check the condition number of the matrix anyway otherwise risk having an other linalg error

hmaarrfk · 2018-05-25T03:08:41Z

Finally @stefanv. I quickly skimmed the paper and the code seems to be following the formalism the paper is using.

I suggest that we leave the reformulation to somebody that wants to derive it. (It May be a paper for itself)

emmanuelle · 2018-09-01T16:49:39Z

So, should we merge this PR?

hmaarrfk · 2018-09-01T22:36:48Z

I think so! I can rebase if you want.

jni · 2018-09-19T08:04:28Z

@hmaarrfk I can't quite tell from the conversation whether you have a test case that will fail on 64 bit? It would be nice for this to be tested.

hmaarrfk · 2018-09-19T14:52:35Z

Sorry, I guess I coudln't think of one. I'll try to bring out my pen and paper this weekend. Remind me.

hmaarrfk · 2018-09-20T00:51:54Z

@jni

=============================================== FAILURES ================================================
__________________________________ test_ellipse_model_estimate_failers __________________________________

    def test_ellipse_model_estimate_failers():
        # estimate parameters of real data
        model = EllipseModel()
        assert not model.estimate(np.ones((5, 2)))
        # Before PR https://github.com/scikit-image/scikit-image/pull/3103
        # This next line would return True on 32bit
        assert not model.estimate(np.array([[50, 80], [51, 81], [52, 80]]))
        # This is the same test that causes the model to
        # return true on 64 bit prooving that this was also an issue on
        # 64 bit arch
>       assert not model.estimate(np.array([[50, 80], [51, 81], [52, 80.00000000001]]))
E       assert not True
E        +  where True = <bound method EllipseModel.estimate of <skimage.measure.fit.EllipseModel object at 0x7f422d07e5c0>>(array([[ 50.,  80.],\n       [ 51.,  81.],\n       [ 52.,  80.]]))
E        +    where <bound method EllipseModel.estimate of <skimage.measure.fit.EllipseModel object at 0x7f422d07e5c0>> = <skimage.measure.fit.EllipseModel object at 0x7f422d07e5c0>.estimate
E        +    and   array([[ 50.,  80.],\n       [ 51.,  81.],\n       [ 52.,  80.]]) = <built-in function array>([[50, 80], [51, 81], [52, 80.00000000001]])
E        +      where <built-in function array> = np.array

skimage/measure/tests/test_fit.py:233: AssertionError
================================== 1 failed, 22 passed in 0.72 seconds ==================================

You can go ahead and force push to remove my "Reverting" of the fix once travis shows that the tests fail. I'm not really a fan of using eps. Apparently it is pretty slow the first time you call it.

Also, you can precompute the inv of C1, and paste it in instead. It is nice round numbers so it stays readable. I figured that adding that in would delay this PR getting merged.

hmaarrfk · 2018-09-20T01:06:47Z

Failed on Appveyor: https://ci.appveyor.com/project/scikit-image/scikit-image/build/1.0.1014/job/qqf67dhmccr8f0so#L8458
Failed on Travis python 3.6: https://travis-ci.org/scikit-image/scikit-image/jobs/430813979#L2803
But passed on python 3.7: https://travis-ci.org/scikit-image/scikit-image/jobs/430813978

hmaarrfk · 2018-09-20T01:20:01Z

Python 3.7 finally failed. I guess adding all those tests was worth it.
https://travis-ci.org/scikit-image/scikit-image/jobs/430817941#L2839

hmaarrfk · 2018-09-20T03:51:27Z

@jni Here is what recreates the issue:

from skimage.measure import EllipseModel
import numpy as np
model = EllipseModel()
print(model.estimate(np.array([[50, 80], [51, 81], [52, 80.00000000001]])))
print(model.estimate(np.array([[50, 80], [51, 81], [52, 80.0001000001]])))
print(model.estimate(np.array([[50, 80], [51, 81], [52, 80.000000001]])))
print(model.estimate(np.array([[50, 80], [51, 81], [52, 80.00000001]])))
print(model.estimate(np.array([[50, 80], [51, 81], [52, 80.0000001]])))
print(model.estimate(np.array([[50, 80], [51, 81], [52, 80.000001]])))
print(model.estimate(np.array([[50, 80], [51, 81], [52, 80.00001]])))

False
False
True
True
True
True
False

I thought you could simply check the condition number, but it fails the doctest, which should clearly pass....

hmaarrfk · 2018-09-20T04:17:38Z

part of me thinks we should chek the number of input points. The problem of fitting an ellipse is ilposed with less than 5 points

jni · 2018-09-24T02:44:49Z

skimage/measure/tests/test_fit.py

+    # 64 bit arch
+    assert not model.estimate(np.array([[50, 80], [51, 81], [52, 80.00000000001]]))
+    # Unfortunately, it seems to depend on the version of python that is used
+    # Therefore, test for many cases


@hmaarrfk based on your comment, shouldn't this test be:

assert not all(model.estimate(np.array([[50, 80], [51, 81], [52, 80 + 10**(-n)]])) for n in range(4, 12))

?

@jni, refactored as per your suggestion. problem remains :/

Which problem? Do they all pass for some Python versions?

Well my fix causes the doctest to fail. But the doctest seems to make sense. https://travis-ci.org/scikit-image/scikit-image/jobs/432304643#L2823

I'm glad to make the test xfail for code self-documentation.

OOOOoooooOOOOhhhh. Very interesting. Hmm, I think ellipse-fitting to a line should work, actually. =\

…ipsis fitting Explicitely check the condition number of the matrix before computing the eigen vectors

Borda · 2019-04-26T23:15:28Z

skimage/measure/tests/test_fit.py

-    assert not model.estimate(np.array([[50, 80], [51, 81], [52, 80]]))
+    assert not model.estimate(np.array([[50, 80],
+                                        [51, 81],
+                                        [52, 80 + epsilon]]))


according to the comment, what is the epsilon number (it is not clear from the review scope)

sorry, quite often, epsilon is a small number in math. I could add a comment, but this wasn't the right approach anyway, and gave the same results as before.

hmaarrfk force-pushed the bugfix_linalg_error_ellipsis_model branch 2 times, most recently from 9fd7a12 to 4d40a5d Compare May 24, 2018 21:32

stefanv reviewed May 24, 2018

View reviewed changes

hmaarrfk changed the title ~~Fixes one bug on a 32 bit machine not caught by np~~ Ellipsis fitting: Check the condition number before computiing eigen vectors May 24, 2018

soupault added the type: bug label May 25, 2018

hmaarrfk force-pushed the bugfix_linalg_error_ellipsis_model branch 5 times, most recently from e2b7248 to d8c098c Compare June 19, 2018 02:06

hmaarrfk force-pushed the bugfix_linalg_error_ellipsis_model branch from d8c098c to e81aadb Compare September 1, 2018 22:37

hmaarrfk force-pushed the bugfix_linalg_error_ellipsis_model branch 2 times, most recently from 5381ec2 to 3eee483 Compare September 20, 2018 00:49

hmaarrfk force-pushed the bugfix_linalg_error_ellipsis_model branch from 3eee483 to 34395c3 Compare September 20, 2018 01:08

hmaarrfk force-pushed the bugfix_linalg_error_ellipsis_model branch 3 times, most recently from 12666fc to 9b5aae0 Compare September 20, 2018 03:00

jni reviewed Sep 24, 2018

View reviewed changes

hmaarrfk force-pushed the bugfix_linalg_error_ellipsis_model branch from 9b5aae0 to 6ce1e2e Compare September 24, 2018 02:58

hmaarrfk added 2 commits September 24, 2018 02:58

TST: Add a 64 bit failer for ellipsis fitting

5888d4e

Fixes scikit-image#3091: check for condition number before during ell…

6ce1e2e

…ipsis fitting Explicitely check the condition number of the matrix before computing the eigen vectors

Borda reviewed Apr 26, 2019

View reviewed changes

rfezzani added 🩹 type: Bug fix Fixes unexpected or incorrect behavior and removed type: bug labels Feb 22, 2020

Base automatically changed from master to main February 18, 2021 18:23

mkcor mentioned this pull request Feb 22, 2021

2021's calendar of community management #5169

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Ellipsis fitting: Check the condition number before computiing eigen vectors #3103

Ellipsis fitting: Check the condition number before computiing eigen vectors #3103

hmaarrfk commented May 24, 2018 •

edited

pep8speaks commented May 24, 2018 •

edited

codecov-io commented May 24, 2018 •

edited

stefanv May 24, 2018

hmaarrfk May 24, 2018

stefanv commented May 24, 2018

matthew-brett commented May 24, 2018

stefanv commented May 24, 2018

hmaarrfk commented May 24, 2018 •

edited

hmaarrfk commented May 24, 2018

hmaarrfk commented May 25, 2018

hmaarrfk commented May 25, 2018

emmanuelle commented Sep 1, 2018

hmaarrfk commented Sep 1, 2018

jni commented Sep 19, 2018

hmaarrfk commented Sep 19, 2018

hmaarrfk commented Sep 20, 2018

hmaarrfk commented Sep 20, 2018

hmaarrfk commented Sep 20, 2018

hmaarrfk commented Sep 20, 2018

hmaarrfk commented Sep 20, 2018

jni Sep 24, 2018

hmaarrfk Sep 24, 2018

jni Sep 24, 2018

hmaarrfk Sep 24, 2018

jni Sep 24, 2018

Borda Apr 26, 2019

hmaarrfk Apr 26, 2019

Ellipsis fitting: Check the condition number before computiing eigen vectors #3103

Are you sure you want to change the base?

Ellipsis fitting: Check the condition number before computiing eigen vectors #3103

Conversation

hmaarrfk commented May 24, 2018 • edited

Checklist

References

For reviewers

pep8speaks commented May 24, 2018 • edited

Comment last updated on September 20, 2018 at 03:00 Hours UTC

codecov-io commented May 24, 2018 • edited

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

stefanv commented May 24, 2018

matthew-brett commented May 24, 2018

stefanv commented May 24, 2018

hmaarrfk commented May 24, 2018 • edited

hmaarrfk commented May 24, 2018

hmaarrfk commented May 25, 2018

hmaarrfk commented May 25, 2018

emmanuelle commented Sep 1, 2018

hmaarrfk commented Sep 1, 2018

jni commented Sep 19, 2018

hmaarrfk commented Sep 19, 2018

hmaarrfk commented Sep 20, 2018

hmaarrfk commented Sep 20, 2018

hmaarrfk commented Sep 20, 2018

hmaarrfk commented Sep 20, 2018

hmaarrfk commented Sep 20, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

hmaarrfk commented May 24, 2018 •

edited

pep8speaks commented May 24, 2018 •

edited

codecov-io commented May 24, 2018 •

edited

hmaarrfk commented May 24, 2018 •

edited