bugfix: prevent integer overflow in EllipseModel

[mark-boer]

Description

This is a bugfix to issue #5042. It makes sure that the EllipseModel and CircleModel do not cause integer overflows

Checklist

I have a question, while trying to implement a unit test I ran into another bug. The test is still implemented, but is marked with a pytest.mark.skip. How would you like to handle this?

Also, does the code need to pass any kind of static analysis? Like black or flake8?

For reviewers

• 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]

Hello @mark-boer! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2021-01-21 02:12:32 UTC

[mark-boer]

Found the reason the test I wrote is failing, the input data chosen, leads to a ridiculously small eigenvalue, which in turn led to significant floating point errors in the matrix inverse. It is probably a solution to first subtract the centroid of the input data and then add this centroid to the center of the Circle. Also I found that the algebraic method used, is not actually a proper least square fit. I'll create a new issue.

[mark-boer]

The behaviour could be improved by using something like this:

A = np.append(xy * 2, np.ones((xy.shape[0], 1)), axis=1)
f = np.sum(xy ** 2, axis=1)
C, _, _, _ = np.linalg.lstsq(A, f, rcond=None)
center = C[0:2]
distances = scipy.spatial.minkowski_distance(center,xy)
r = np.sqrt(np.mean(distances ** 2))

See e.g. sphere fit i wrote a year ago, this contains 2 fits, one fast algebraic one and one iterative one:
https://github.com/mark-boer/geomfitty/blob/4f24071f524302a95144dbd91bc84fa33c5dc26b/geomfitty/fit3d.py#L42

[stefanv]

Thanks for the fix, @mark-boer. I did a quick check and it looks like vstack does type promotion already. Could you explain the rationale behind the fix (casting the data to float)?

Before merging, we should remove the failing test and address that as part of your other issue.

[mark-boer]

The vstack on which line do you mean? The EllipseModel overflows in the following line:

S1 = D1.T @ D1

For more info on what goes wrong see issue #5042.

[jni]

@stefanv I think the point of the cast might be to keep things at the same float bit depth for both D1 and D2? (ie both float32 or both float64, depending on input dtype).

[jni]

I agree with @stefanv that ideally the skipped test should instead get added in a separate PR that fixes the issue, but I'm not that fussed if we instead just merge this as-is and remove the skip marker in a follow-up PR.

@scikit-image/core any others up for review?

Thank you @mark-boer!

[mark-boer]

@jni, I changed the input data so it cast to a float of at least 32 bits. But if the user chooses to call the function with a float64, it is not cast down to a float32. If that doesn't make sense, we can also just forcibly cast it to a float64.

I will remove the failing test this evening and create a new PR with the fix to this failing test.

[jni]

I changed the input data so it cast to a float of at least 32 bits. But if the user chooses to call the function with a float64, it is not cast down to a float32. If that doesn't make sense, we can also just forcibly cast it to a float64.

Makes perfect sense, actually! I think the fact that dtype was set for the second vstack and not the first confused @stefanv upon reading the code. Took me a while as well after reading his comment. Perhaps adding a dtype= to the vstack calls, even if redundant, will make the code clearer? Or maybe just add a comment to the effect of what you just wrote?

[mark-boer]

The dtype, is actually an argument of the np.ones, even though this is not really necessary , I thought that saved an additional cast.

edit: I will add a comment to the code 😉 that explains the why

[rfezzani]

What about

Suggested change

	D2 = np.vstack([x, y, np.ones(len(x), dtype=float_type)]).T
	D2 = np.vstack([x, y, np.ones_like(x)]).T

[mark-boer]

That's even better :-)

[stefanv]

Thanks, this is very readable now @mark-boer. I left some minor suggestions, but they're mostly stylistic.

[jni]

Just gonna approve this some more for good measure. =P

[stefanv]

Thank you very much, @mark-boer