Solid Ellipse Fitting on 2D Image Using Gradient Descent

Question

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()
ax.add_patch(ell)
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,

Answer 1

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()
     ax.add_patch(ell)
     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

Answer 2

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.

Answer 3

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.