fitting an ellipse inside of the boundary of a polygon

I am struggling with fitting the biggest possible ellipse (biggest possible are of the ellipse) inside of a polygon. So the ellipse should not go outside of the polygon border, which is not the case for the common "best fitting ellipse" algorithms. The polygon/desired form is saved as binary mask. (i will enclose one so you can better understand the task).

I have already looked for algorithms for this task but was not able to find one that works here.
Do you know any working algorithm (preferrably for python) or have an idea what would be the best approach here or maybe even have a code for this?

thank you very much!

This problem is called the maximum volume inscribed ellipsoid, and can be solved in any dimension. The problem and its solution are presented in Convex optimization book by Boyd & Vandenberghe.

Below is a Python implementation. It uses pycddlib (cdd) to get the H-representation of the polygon, that is to say to write it as the solution of a set of inequalities Ax <= b. Then it uses cvxpy to solve the problem.

The implementation returns the shape matrix B of the elliipse and its center d. That means the ellipse is the set of points of the form Bu + d for u in the unit ball. You can write u = [cos(theta), sin(theta)] and you'll get a usual equation for an ellipse, useful for plotting.

import cvxpy as cp
import numpy as np
import cdd as pcdd

# vertices of the polyhedron
pts = np.array([
    [0, 0],
    [0, 1],
    [1, 0]    
])
npoints = pts.shape[0]

# make the V-representation of the points; you have to prepend the points with a column of ones
V = np.column_stack((np.ones(npoints), pts))
mat = pcdd.Matrix(V, number_type='fraction') # use fractions if possible
mat.rep_type = pcdd.RepType.GENERATOR
poly = pcdd.Polyhedron(mat)

# H-representation of the points: obtain matrix A and vector b such that 
# the polyhedron is defined by Ax <= b
H = poly.get_inequalities()
# get A and b
G = np.asarray(H[0:])
b = G[:,0]
A = -G[:,1:3]

# now we define the convex optimization problem
n = 2 # dimension
# variables
Bvar = cp.Variable((n, n), symmetric=True) # n x n symmetric matrix
dvar = cp.Variable(n)      # vector
# constraints
constraints = [] # to store the constraints
for i in range(npoints):
    constraints.append(
        cp.norm(Bvar @ A[i,:]) + A[i,:] @ dvar <= b[i] 
    )
# objective
objective = cp.Minimize(-cp.log_det(Bvar))
# problem
problem = cp.Problem(objective, constraints)
# solve problem
x = problem.solve(solver=cp.SCS)
# solutions
variables = problem.solution.primal_vars
print("The shape matrix B is")
print(Bvar.value)
print("The center d is")
print(dvar.value)

The shape matrix B is
[[ 0.32199634 -0.08628191]
 [-0.08628191  0.32195628]]
The center d is
[0.3333573  0.33331871]

from matplotlib import pyplot as plt
from math import pi

polygon = np.append(pts, [pts[0,:]], axis = 0)

B = Bvar.value
d = dvar.value
d = np.transpose(np.repeat([d], 100, axis = 0))
t = np.linspace(0, 2*pi, 100)
u = np.array([np.cos(t) , np.sin(t)])
coords = B @ u + d

f, ax = plt.subplots(1)
ax.plot(polygon[:,0], polygon[:,1])
ax.plot(coords[0,:], coords[1,:] )
ax.grid(color='lightgray', linestyle='--')
ax.show()