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.