Polygon approximation with a circle

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

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?

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)