# Solid Ellipse Fitting on 2D Image Using Gradient Descent

I am attempting to fit an ellipse for a specific color, $\mu$ at grayscale, on an image that will cover as much of the region as possible, with the targeted color inside the ellipse. This is not contour fitting, or simply fitting edges, my goal is to get the largest possible ellipse that has, in it, the color I am looking for - $\mu$. I tried to construct an error function so I can take the partial derivatives, form a gradient and move in the direction of that gradient.

The formula for an ellipse, ignoring the the tilt for now, is

$$R(x,y;h,k,r_x,r_y) = \frac{(x-h)^2}{r_x^2} + \frac{(y-k)^2}{r_y^2} = 1$$

I know that for $x,y$ points outside $R$ will result in $>1$ and for the ones inside $<1$. I thought I could use this to create an error / objective function which I can minimize.

$$E = \sum_{i=1}^n R(x_i,y_i;h,k,r_x,r_y) (I(x_i,y_i)-\mu)^2$$

where $I(x,y)$ is the intensity function for the grayscale value at location $x,y$. Let's call $C_i = (I(x_i,y_i)-\mu)^2$

$$\frac{\partial }{\partial h} \bigg( \sum_{i=1}^n \frac{C_ix_i^2}{r_x^2} - \frac{2C_i x_i h}{r_x^2} + \frac{C_ih^2}{r_x^2} + ... \bigg) = \sum_{i=1}^n \frac{- 2C_i (x_i - h)}{r_x^2}$$

$$\frac{\partial }{\partial k} = \sum_{i=1}^n \frac{- 2C_i (y_i - k)}{r_y^2}$$

$$\frac{\partial }{\partial r_x} = \sum_{i=1}^n \frac{-2C_i x_i^2 + 4C_i x_i h - 2C_i h^2}{r_x^3} = \sum_{i=1}^n \frac{-2C_i(x_i^2 - 2x_i h + h^2)}{r_x^3} = \sum_{i=1}^n \frac{-2C_i(x_i-h)^2}{r_x^3}$$

$$\frac{\partial }{\partial r_y} = \sum_{i=1}^n \frac{-2C_i(y_i-k)^2}{r_y^3}$$

I implemented this method with;

from PIL import Image
img = Image.open('ellipsedented.png').convert('L')
A = np.array(img)

x = []; y=[]; I=[]
for i in  range(A.shape[1]):
for j in range(A.shape[0]):
x.append(float(i))
y.append(A.shape[0]-float(j)) # images have (0,0) on upper left
I.append(A[j,i])
x = np.array(x); y=np.array(y); I=np.array(I)

h = 50; k = 20; rx = 20; ry = 20; mu=150; eta = 0.000001
iter = 1000
C = ((I-mu)**2) + 0.01
for i in range(iter):
tmp = -2*C*(x-h) / rx**2
h_step = tmp.sum() * eta
tmp = -2*C*(y-k) / ry**2
k_step = tmp.sum() * eta
tmp = -2*C*(x-h)**2 / np.power(rx,3)
rx_step = tmp.sum() * eta
tmp = -2*C*(y-k)**2 / np.power(ry,3)
ry_step = tmp.sum() * eta
h = h - h_step
k = k - k_step
rx = rx - rx_step
ry = ry - ry_step
#break
print h,k,rx,ry

# for plots
from matplotlib.patches import Ellipse
ell = Ellipse(xy=[h,A.shape[0]-k], width=rx, height=ry)
fig, ax = plt.subplots()
plt.hold(True)
plt.imshow(img)
plt.show()


The results I am seeing:

116.348736677 69.0497488464 173.708913479 115.869883583


which is not great. More iterations moves ellipse more on top of the real one, but the model will also keep growing. Any help to get this approach to an acceptable working condition would be greatly appreciated. Are the partial derivatives correct? Can I change the formulation for the minimizer to make it more robust? I would like to get this approach to work as best as possible - it will be used as an example on the uses of gradient descent.

Thanks,

• your objective function seems ... not correct. in your device, the contour has obviously a great effect – reuns Sep 7 '16 at 9:39
• What do you want to do ? A good solution in general is given the ellipse estimated position, compute a mean value for its color inside AND color outside, and compute the (weighted) squared error between the image and your model of ellipse (with the color you just computed) – reuns Sep 7 '16 at 10:49
• And when I said weighted, I meant an high weight on the contour, but a low weight far outside the ellipse. – reuns Sep 7 '16 at 10:51
• Estimated position and size is all I need: color is specified from outside, chosen by the user. Say given grayscale 155, and the goal is to find an ellipse that covers as much of the ellipse with the 155'ish pixels in it (picture above) as possible - dents and all. – BB_ML Sep 7 '16 at 10:56
• Note: I made some correction to formulas and code - still not great, but results do not look too insane either. – BB_ML Sep 7 '16 at 10:58

An effective approach is by performing the Principal Component Analysis, which will give you the center and main axis of the ellipse. In 2D all computation can be done analytically.

In your particular case, the ellipse is perturbed by dents. For the most accurate results, you should extract the outline pixels and perform robust ellipse fitting using one of the published methods. Anyway, the above method will provide a pretty good starting approximation.

• PCA can be used on this problem, yes; I would like to use gradient descent to minimize an objective function however, computing partial derivatives, starting from an initial point and moving in the direction of gradient.. Is the objective convex? Maybe $\log$ is better for $R$? .. Etc.. – BB_ML Sep 7 '16 at 5:41
• Least-squares won't perform better than PCA (and be much slower) because of the dent. You need a robust method. – Yves Daoust Sep 7 '16 at 6:32
• Hi - thanks for the suggestion on PCA. I might look at that after this, first I'd like to get this approach to an acceptable working condition. Is it impossible for this GD to converge to roughly the correct answer? One would think the method, while trying to minimize the objective function, make the ellipse as large as possible to fit most of the $\mu$ colored points in it. I also made some additions to the question - for the derivatives and code. Thanks, – BB_ML Sep 7 '16 at 8:33
• @BB_ML: well, your objective function will lead you nowhere (more precisely the solution that minimizes $R$ degenerates to a single point, the center of the ellipse). You can minimize the sum of (squared) distance from the outline points to the ellipse instead. Alternatively, $(R-1)^2$. If you want to use all pixels (no just the outline), I guess that you must switch to a constrained minimization problem. – Yves Daoust Sep 7 '16 at 8:48
• @mathreadler: actually the tool that you need is the ellipse of inertia, but as this term is infrequent I preferred to speak of PCA. – Yves Daoust Sep 7 '16 at 14:17

A nother approach if you want to create a kind of envelope - to focus on fitting in the outer points into the ellise is this approach. If we first for each point build a matrix as the outer product of (column) vectors:

$${\bf M = (v-v_c)(v-v_c)^T}$$

where $\bf v_c$ is an estimated center-point of the object. We can then be sure that all $\bf M$ in the image will be symmetric positive semi-definite. We can create a sequence of matrix exponent pairs, say for example $\left(k+2,\frac{1}{k}\right), k\in\mathbb N$. We then apply the first to the individual terms and the second on the sum. If one wants to one can verify that for large $k$ this will asymptotically give the same effect as the max-estimator.

The low k:s focus more on fitting co-variances, and the high focus more on fitting the max-eigensystem.

I believe I have something that could give me the global minimum when ellipse is right on the dented ellipse image. I make it so that $C_i = -\log (|I - \mu|) / M$ where $M$ is a constant, like 50. So I look for 150, all 155 pixels become $-\log(|155-150|/50)$, this gives me negatives for unwanted colors, positives for wanted ones. At the same time in $E$ I use $\log R_i C_i$. Taking the log of $R$ gives me negatives inside the ellipse, positive outside it. The minimum of this multiplication is where the ellipse is on wanted points, because that's the only place where negative multiplies positive (inside) and positive multiplies negative (outside).

Code

import pandas as pd
from matplotlib.patches import Ellipse
from PIL import Image

img = Image.open('ellipsedented.png').convert('L')
A = np.array(img)

x = []; y=[]; I=[]
for i in  range(A.shape[1]):
for j in range(A.shape[0]):
x.append(float(i))
y.append(A.shape[0]-float(j)) # imalarin (0,0) noktasi ust solda
I.append(A[j,i])
x = np.array(x); y=np.array(y); I=np.array(I)

h = 20; k = 20; rx = 30; ry = 30; mu=150; eta = 0.01
M = 50; iter = 1000

C =  -1*np.log((np.abs(I-mu) / M) + 0.001)

def R(xx,yy,hh,kk,rxx,ryy):
return ((xx-hh)**2 / rxx**2) + ((yy-kk)**2 / ryy**2)

def plot_ellipse_image(h,k,rx,ry,i):
ell = Ellipse(xy=[h,A.shape[0]-k], width=rx, height=ry, alpha=0.2)
fig, ax = plt.subplots()
plt.hold(True)
plt.imshow(img)
plt.savefig('ellipse_fitting_%d.png' % i)

for i in range(iter):
old = np.array([h,k,rx,ry])
r = R(x,y,h,k,rx,ry)
tmp = (-2*C*(x-h) / rx**2) / r
h_step = pd.Series(tmp).fillna(0).sum() * eta
tmp = (-2*C*(y-k) / ry**2) / r
k_step = pd.Series(tmp).fillna(0).sum() * eta
tmp = (-2*C*(x-h)**2 / np.power(rx,3)) / r
rx_step = pd.Series(tmp).fillna(0).sum() * eta / delay
tmp = (-2*C*(y-k)**2 / np.power(ry,3)) / r
ry_step = pd.Series(tmp).fillna(0).sum() * eta / delay
h = h - h_step
k = k - k_step
rx = rx - rx_step
ry = ry - ry_step
new = np.array([h,k,rx,ry])
if i % 100 == 0:
plot_ellipse_image(h,k,rx,ry,i)
print h,k,rx,ry,np.abs((new-old).sum()), i
if np.abs((new-old).sum()) < 0.40: break