incomplete ellipse image fitting by python

Question

I am fitting an incomplete ellipse-an ellipse cap (it is an image already processed). the problem is that, as you can see, at the bottom of fitted ellipse, two small parts on the left and right are not included in the fitted ellipse. These two segments are quite important and I should include them in the fitted ellipse. I tried RANSAC in scikit-image and fitEllipse in opencv2, but none of them can do this, and as I know RANSAC cannot take weight, which I want to put on the two outlined parts. The ideal result would be that an ellipse will include the two segments, while the upper part will be slightly larger than the raw ellipse. The code I am using is available upon request.

raw ellipse

no ideal ellipse fitted

ideal ellipse fitted

Answer 1

I tried hacking this together and I got something close to what you want, but it's finicky with the initial conditions.

import skimage as sk
import numpy as np
import matplotlib.pyplot as plt
import lmfit
# Loading your image and removing what I believe are some artifacts from saving.
img = sk.io.imread("mLSEYuHD.png")
img = (img > 0).astype(int)
img = sk.morphology.binary_opening(
    sk.morphology.binary_opening(sk.morphology.binary_opening(img))
)
lbl = sk.measure.label(img)
lbl = (lbl == 1).astype(int)
# Use regionprops to get some initial estimates of the ellipse
props = sk.measure.regionprops(lbl)
props = props[0]
props.centroid, props.axis_major_length, props.axis_minor_length, props.orientation
# ((215.51666590357584, 240.36841261846985),
# 237.78103980008083,
# 236.05236540309176,
# 1.0263988084463667)

# Define the formula for an ellipse. If it's <=1, it's inside the ellipse.
def ellipse_formula(x, y, center_x, center_y, major, minor, angle_rad):
    cos_angle = np.cos(angle_rad)
    sin_angle = np.sin(angle_rad)

    term1 = ((x - center_x) * cos_angle + (y - center_y) * sin_angle) ** 2 / major**2
    term2 = ((x - center_x) * sin_angle - (y - center_y) * cos_angle) ** 2 / minor**2

    return term1 + term2

# Since i'm using lmfit to get the shape parameters, I'm initializing the parameter object with the guesses from the distribution. I found that these have a great influence in the fit, but perhaps, with some work, you'll get better initial guesses.

params = lmfit.Parameters()
params.add("cx", value=215, min=0)
params.add("cy", value=240, min=0)
params.add("major", value=300 / 2, min=0)
params.add("minor", value=300 / 2, min=0)
params.add("ang", value=0)

# Defining a residual. Essentially, this sums up the pixels outside the ellipse and tries to minimize that. But if it finds pixels that should be inside, it'll return a huge number.

def residual(params, img):
    cx = params["cx"].value
    cy = params["cy"].value
    maj = params["major"].value
    min = params["minor"].value
    ang = params["ang"].value

    x = np.arange(0, img.shape[1])
    y = np.arange(0, img.shape[0])
    xx, yy = np.meshgrid(x, y, indexing="ij")
    ell_mask = ellipse_formula(xx, yy, cx, cy, maj, min, ang) <= 1
    if np.sum(img[~ell_mask]) > 0:
        return 1000000000
    else:
        return np.sum(~img[ell_mask]) # Sum of background pixels inside the ellipse, which should be as low as possible

# I try to modify this to get the ellipse roughly around the shape I want. If it's too fa

plt.imshow(lbl)
x = np.arange(0, lbl.shape[1])
y = np.arange(0, lbl.shape[0])
xx, yy = np.meshgrid(x, y, indexing="ij")
ell_mask = (
    ellipse_formula(
        xx,
        yy,
        params["cx"],
        params["cy"],
        params["major"],
        params["minor"],
        params["ang"],
    )
    <= 1
)
plt.imshow(ell_mask, cmap="gray", alpha=0.5)
residual(params, lbl.astype(bool))

# Fitting and getting the optimized parameters

fit = lmfit.minimize(residual, params, args=[lbl.astype(bool)], method="nelder")

# Visualizing the output

plt.imshow(lbl)
x = np.arange(0, lbl.shape[1])
y = np.arange(0, lbl.shape[0])
xx, yy = np.meshgrid(x, y, indexing="ij")
ell_mask = (
    ellipse_formula(
        xx,
        yy,
        fit.params["cx"],
        fit.params["cy"],
        fit.params["major"],
        fit.params["minor"],
        fit.params["ang"],
    )
    <= 1
)
plt.imshow(ell_mask, cmap="gray", alpha=0.5)
residual(params, lbl.astype(bool))

This was my output. It's pretty close to what you wanted I think.