Polygon approximation with a circle

Question

There is an article describing how to detect a pupil from an eye photo.

• A.-H. Javadi, Z. Hakimi, M. Barati, V. Walsh, and L. Tcheang, "SET: a pupil detection method using sinusoidal approximation," Front. Neuroeng., vol. 8, pp. 1–10, Apr. 2015

On the 4th step of pupil detection they use Decomposition of each segment into sinusoidal components method that I do not understand.

They say that

We fit Equation 1 to the extracted points along the border to extract the free parameters $\omega$, $(h,k)$, $a$ and $b$.

But I do not understand how exactly they came from polygon segments to the ellipse free parameters.

Should I use some obvious computational algorithms for these step to get them?

Answer 1

They just curve fit the two equations $$ x(\omega) = h - b \sin(\omega) ~,~~ y(\omega) = k + a \cos(\omega) $$ using a nonlinear curve fitting algorithm. However, the accuracy of the fit depends on the choice of algorithm and you have to be careful to make sure that $a$ and $b$ are always positive (by adding constraints).

One way of approaching the problem is shown below (using Octave).

1) Set up the points

fig1 = figure()
gca
[x, y, buttons] = ginput(10)
plot(x, y, 'x-');
hold on;

2) Compute convex hull

k = convhull(x,y)
plot(x(k), y(k), "r-")

3) Remove points inside hull

xx = zeros(length(k)-1)
yy = zeros(length(k)-1)
for i=1:length(k)-1
  xx(i) = x(k(i));
  yy(i) = y(k(i));
end

4) Compute centroid of hull

x_cen = mean(xx);
y_cen = mean(yy);
plot(x_cen, y_cen, 'mo');

5) Compute angles and plot x,y vs. angle

for i=1:length(xx)
  theta(i) = atan2(yy(i)-y_cen, xx(i)-x_cen);
end
angle = theta*180/pi
fig2 = figure();
plot(angle, yy, 'r-'); hold on;
plot(angle, xx, 'b-');

6) Load optimization package

pkg load optim

7) Set up curve fit

indep_omega = theta(1:length(theta))
obs_x = xx(1:length(theta)) - x_cen
obs_y = yy(1:length(theta)) - y_cen

% Model functions
fx_omega = @ (p, omega) p(1) - p(2) * sin(omega)
fy_omega = @ (p, omega) p(1) + p(2) * cos(omega)

% Initial values
init_p_x = [0; 0.5]
init_p_y = [0; 0.5]

% Constraints (A.' * p + B >= 0)
A = [0; 1] ; B = 0;
settings = optimset ("inequc", {A, B});

8) Do curve fit

[p_x, model_values, cvg, outp] = nonlin_curvefit(fx_omega, init_p_x, indep_omega, obs_x', settings)
[p_y, model_values, cvg, outp] = nonlin_curvefit(fy_omega, init_p_y, indep_omega, obs_y', settings)