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)