I am cross-posting this from GitHub as I am not sure whether the issue is a problem with my data or a bug. https://github.com/scikit-image/scikit-image/issues/6699
I have thousands of elliptical features in my microscopy data that I want to fit models to using skimage. The model fails on some for no obvious reason. Here's code to reproduce:
import numpy as np
from skimage.measure import EllipseModel
import plotly.express as px
good_x1 = [779.026, 778.125, 776.953, 776.195, 775.617, 775.068, 774.127, 773.696, 773.305, 773.113, 773.088, 773.233, 773.449, 773.913, 774.344, 774.625, 775.179, 775.777, 776.254, 777.039, 777.926, 778.945, 780.023, 781.059, 781.973, 782.777, 783.244, 783.922, 784.995, 785.825, 786.196, 786.486, 786.65, 786.614, 786.482, 786.153, 785.749, 785.507, 784.901, 784.482, 783.879, 782.809, 781.965, 780.998, 780.001, 779.026]
good_y1 = [309.143, 309.432, 309.912, 310.35, 310.46, 311.087, 312.099, 312.879, 314.085, 315.012, 315.995, 316.948, 318.166, 319.044, 319.751, 320.283, 320.794, 321.34, 321.505, 321.908, 322.254, 322.478, 322.467, 322.243, 321.929, 321.561, 321.449, 320.891, 319.995, 318.905, 318.07, 316.872, 315.97, 315.037, 313.883, 312.943, 312.17, 311.623, 311.093, 310.477, 310.151, 309.54, 309.18, 309.027, 309.022, 309.143]
good_x2 = [434.959, 434.0, 433.012, 432.093, 430.938, 430.279, 429.847, 429.535, 429.257, 429.031, 428.843, 429.0, 429.348, 429.872, 430.313, 431.048, 432.189, 433.043, 434.003, 434.971, 435.769, 436.199, 436.743, 437.263, 437.824, 438.017, 438.018, 437.831, 437.449, 437.29, 436.807, 436.255, 435.776, 434.959]
good_y2 = [215.849, 216.001, 215.929, 215.684, 215.09, 214.615, 214.117, 213.631, 212.903, 211.992, 211.017, 210.0, 209.39, 208.857, 208.587, 208.087, 207.57, 207.247, 207.135, 207.2, 207.565, 207.73, 208.248, 208.819, 210.055, 210.998, 212.001, 212.952, 213.687, 214.168, 214.781, 215.333, 215.49, 215.849]
good_x3 = [1666.998, 1666.014, 1665.206, 1664.689, 1664.302, 1663.977, 1663.969, 1664.293, 1664.527, 1665.09, 1665.929, 1667.048, 1668.016, 1668.658, 1669.171, 1669.638, 1669.599, 1668.995, 1667.916, 1666.998]
good_y3 = [85.023, 85.07, 85.414, 85.685, 86.245, 86.994, 88.004, 88.835, 89.364, 89.862, 90.302, 90.338, 90.034, 89.491, 89.134, 87.917, 86.807, 86.004, 85.251, 85.023]
bad_x1 = [1541.221, 1541.848, 1543.009, 1544.15, 1544.962, 1545.777, 1545.943, 1545.786, 1545.103, 1543.986, 1543.14, 1541.968, 1541.094, 1540.765, 1540.799, 1541.221]
bad_y1 = [1254.78, 1255.29, 1255.535, 1255.395, 1254.945, 1253.922, 1253.0, 1252.063, 1250.892, 1250.374, 1250.401, 1250.959, 1252.049, 1252.968, 1254.069, 1254.78]
bad_x2 = [1739.079, 1738.567, 1738.392, 1738.118, 1738.17, 1738.782, 1739.302, 1740.179, 1741.013, 1741.999, 1742.997, 1743.423, 1744.178, 1743.811, 1743.735, 1743.595, 1743.308, 1742.834, 1742.342, 1741.813, 1740.998, 1739.995, 1739.079]
bad_y2 = [329.807, 329.316, 328.814, 327.989, 327.061, 325.853, 325.22, 324.478, 324.115, 324.078, 324.154, 324.49, 324.753, 325.994, 326.902, 327.679, 328.143, 328.836, 329.41, 329.628, 329.99, 330.067, 329.807]
bad_x3 = [992.001, 991.057, 989.879, 989.599, 989.252, 989.286, 989.894, 991.05, 991.983, 992.806, 993.286, 993.846, 994.32, 994.481, 994.088, 992.959, 992.001]
bad_y3 = [136.048, 136.19, 136.883, 137.551, 138.053, 138.929, 140.102, 140.767, 140.846, 140.551, 140.416, 139.851, 139.115, 137.938, 136.94, 136.168, 136.048]
xy = [(good_x1, good_y1, 'good_1'), (good_x2, good_y2, 'good_2'), (good_x3, good_y3, 'good_3'),
(bad_x1, bad_y1, 'bad_1'), (bad_x2, bad_y2, 'bad_2'), (bad_x3, bad_y3, 'bad_3')]
for ii in xy:
points = list(zip(ii[0], ii[1]))
a_points = np.array(points)
model = EllipseModel()
if model.estimate(a_points) == False:
fig = px.line(x= ii[0], y= ii[1], title='model fitting failed for ' + ii[2])
fig.show()
try:
xc, yc, a, b, theta = model.params
print(model.params)
ellipse_centre = (xc, yc)
residuals = model.residuals(a_points)
print(residuals)
except Exception as e: print(e)
else:
fig = px.line(x= ii[0], y= ii[1], title='model fitting successful for ' + ii[2])
fig.show()
xc, yc, a, b, theta = model.params
print(model.params)
ellipse_centre = (xc, yc)
residuals = model.residuals(a_points)
print(residuals)
Visually these features all seem elliptical and there is no difference between them in the length of the point lists or other properties. I think it's a bug but would appreciate a 2nd opinion, thanks.
Addition 2023-01-25: There is another implementation of the same algorithm underlying the skimage function at SciPy, here: https://scipython.com/blog/direct-linear-least-squares-fitting-of-an-ellipse/
I have tested this on the same data and it variously fits, misfits or throws errors on my sets of points. Interesting it does this in a slightly different pattern to the skimage function. I think this problem is a fundamental flaw in the algorithm and not a bug.
Your x, y data points differ from their means by only a small amount, so the algorithms you are using are running into numerical errors which lead to failures in the linear algebra routines used in the fitting process.
You can solve this by subtracting off the means (i.e. spatially shifting the ellipses) so that there is a greater relative difference between the values: e.g. the following works for me:
import numpy as np
from skimage.measure import EllipseModel
import plotly.express as px
good_x1 = [779.026, 778.125, 776.953, 776.195, 775.617, 775.068, 774.127, 773.696, 773.305, 773.113, 773.088, 773.233, 773.449, 773.913, 774.344, 774.625, 775.179, 775.777, 776.254, 777.039, 777.926, 778.945, 780.023, 781.059, 781.973, 782.777, 783.244, 783.922, 784.995, 785.825, 786.196, 786.486, 786.65, 786.614, 786.482, 786.153, 785.749, 785.507, 784.901, 784.482, 783.879, 782.809, 781.965, 780.998, 780.001, 779.026]
good_y1 = [309.143, 309.432, 309.912, 310.35, 310.46, 311.087, 312.099, 312.879, 314.085, 315.012, 315.995, 316.948, 318.166, 319.044, 319.751, 320.283, 320.794, 321.34, 321.505, 321.908, 322.254, 322.478, 322.467, 322.243, 321.929, 321.561, 321.449, 320.891, 319.995, 318.905, 318.07, 316.872, 315.97, 315.037, 313.883, 312.943, 312.17, 311.623, 311.093, 310.477, 310.151, 309.54, 309.18, 309.027, 309.022, 309.143]
good_x2 = [434.959, 434.0, 433.012, 432.093, 430.938, 430.279, 429.847, 429.535, 429.257, 429.031, 428.843, 429.0, 429.348, 429.872, 430.313, 431.048, 432.189, 433.043, 434.003, 434.971, 435.769, 436.199, 436.743, 437.263, 437.824, 438.017, 438.018, 437.831, 437.449, 437.29, 436.807, 436.255, 435.776, 434.959]
good_y2 = [215.849, 216.001, 215.929, 215.684, 215.09, 214.615, 214.117, 213.631, 212.903, 211.992, 211.017, 210.0, 209.39, 208.857, 208.587, 208.087, 207.57, 207.247, 207.135, 207.2, 207.565, 207.73, 208.248, 208.819, 210.055, 210.998, 212.001, 212.952, 213.687, 214.168, 214.781, 215.333, 215.49, 215.849]
good_x3 = [1666.998, 1666.014, 1665.206, 1664.689, 1664.302, 1663.977, 1663.969, 1664.293, 1664.527, 1665.09, 1665.929, 1667.048, 1668.016, 1668.658, 1669.171, 1669.638, 1669.599, 1668.995, 1667.916, 1666.998]
good_y3 = [85.023, 85.07, 85.414, 85.685, 86.245, 86.994, 88.004, 88.835, 89.364, 89.862, 90.302, 90.338, 90.034, 89.491, 89.134, 87.917, 86.807, 86.004, 85.251, 85.023]
bad_x1 = [1541.221, 1541.848, 1543.009, 1544.15, 1544.962, 1545.777, 1545.943, 1545.786, 1545.103, 1543.986, 1543.14, 1541.968, 1541.094, 1540.765, 1540.799, 1541.221]
bad_y1 = [1254.78, 1255.29, 1255.535, 1255.395, 1254.945, 1253.922, 1253.0, 1252.063, 1250.892, 1250.374, 1250.401, 1250.959, 1252.049, 1252.968, 1254.069, 1254.78]
bad_x2 = [1739.079, 1738.567, 1738.392, 1738.118, 1738.17, 1738.782, 1739.302, 1740.179, 1741.013, 1741.999, 1742.997, 1743.423, 1744.178, 1743.811, 1743.735, 1743.595, 1743.308, 1742.834, 1742.342, 1741.813, 1740.998, 1739.995, 1739.079]
bad_y2 = [329.807, 329.316, 328.814, 327.989, 327.061, 325.853, 325.22, 324.478, 324.115, 324.078, 324.154, 324.49, 324.753, 325.994, 326.902, 327.679, 328.143, 328.836, 329.41, 329.628, 329.99, 330.067, 329.807]
bad_x3 = [992.001, 991.057, 989.879, 989.599, 989.252, 989.286, 989.894, 991.05, 991.983, 992.806, 993.286, 993.846, 994.32, 994.481, 994.088, 992.959, 992.001]
bad_y3 = [136.048, 136.19, 136.883, 137.551, 138.053, 138.929, 140.102, 140.767, 140.846, 140.551, 140.416, 139.851, 139.115, 137.938, 136.94, 136.168, 136.048]
xy = [(good_x1, good_y1, 'good_1'), (good_x2, good_y2, 'good_2'), (good_x3, good_y3, 'good_3'),
(bad_x1, bad_y1, 'bad_1'), (bad_x2, bad_y2, 'bad_2'), (bad_x3, bad_y3, 'bad_3')]
for ii in xy:
x = np.array(ii[0])
y = np.array(ii[1])
x = x - np.mean(x)
y = y - np.mean(y)
points = list(zip(x, y))
a_points = np.array(points)
model = EllipseModel()
if model.estimate(a_points) == False:
fig = px.line(x=x, y=y, title='model fitting failed for ' + ii[2])
fig.show()
try:
xc, yc, a, b, theta = model.params
print(model.params)
ellipse_centre = (xc, yc)
residuals = model.residuals(a_points)
print(residuals)
except Exception as e: print(e)
else:
fig = px.line(x=x, y=y, title='model fitting successful for ' + ii[2])
fig.show()
xc, yc, a, b, theta = model.params
print(model.params)
ellipse_centre = (xc, yc)
residuals = model.residuals(a_points)
print(residuals)