CircleModel does not minimize squared distances

[mark-boer]

Description

While working on the PR #5179 for a integer overflow in CircleModel I recognized that the CircleModel does not properly minimize the square of the residuals. There is an excellent explanation of this in the scipy-cookbook. It currently is extra sensitive to outliers, because it minimizes: sum( (Ri**2 - Rc**2)**2)

I didn't have a look at the EllipseModel, there might be a similar problem there as well.

Way to reproduce

x = np.r_[36, 36, 19, 18, 33, 26]
y = np.r_[14, 10, 28, 31, 18, 26]
xy = np.stack([x,y], axis=1).astype(float)

model = CircleModel()
model.estimate(xy)

print(np.sum(model.residuals(xy) ** 2))

def func(params):
    model = CircleModel()
    model.params = tuple(params)
    return model.residuals(xy)

res = scipy.optimize.least_squares(func, x0=model.params)

model.params = tuple(res.x)

print(np.sum(model.residuals(xy) ** 2))

This results in:

5.1976653458467945
4.042521898322268

Version information

# Paste the output of the following python commands
from __future__ import print_function
import sys; print(sys.version)
import platform; print(platform.platform())
import skimage; print("scikit-image version: {}".format(skimage.__version__))
import numpy; print("numpy version: {}".format(numpy.__version__))

3.7.6 (tags/v3.7.6:43364a7ae0, Dec 19 2019, 00:42:30) [MSC v.1916 64 bit (AMD64)]
Windows-10-10.0.19041-SP0
scikit-image version: 0.16.2
numpy version: 1.17.4

[sciunto]

thanks for reporting. As a side note, it would be interesting to compare with scikit learn. See #2698
PS: SL does not provide a circle model afaik.

[sciunto]

I'm on it.

[jni]

That's a great reference page @mark-boer, thanks! And thanks @sciunto for tackling this! 💪

[mark-boer]

If you want you can also have look here, it is a sphere fit, which is almost equivalent to the circle fit and the code is a little bit cleaner than the scipy-cookbook code:

https://github.com/mark-boer/geomfitty/blob/4f24071f524302a95144dbd91bc84fa33c5dc26b/geomfitty/fit3d.py#L42

It contains a fast_fit, which is algebraic and a iterative fit which calculates the Least squared distances. (tho is see it contains a small error, when calculating radius; it uses the changed error function to calculate the center and then the least squares to calculate the radius). I added a code snippit to my PR that solves this:

#5179 (comment)

[sciunto]

I implemented a ransac circlemodel v2 that I'm currently benchmarking. My implementation is based on optimize.leastsq (mentioned to me the most efficient) in the scipy cookbook.

So far, my version is 3 times slower than the current scikit-image version, but I notice a significant improvement regarding both the radius and the center coordinates; it is way more reproducible and closer to the expected result, especially for the radius (10 times more accurate).