Circle Fit Demo generated using MATLAB MCP Server with VSCode & MATLAB 2025b

README.md

circle_fit_demo

Compact MATLAB demo for fitting circles and ellipses to noisy pupil-boundary samples (useful for pupil localization / gaze tracking).

This repository contains a single script, circle_fit_demo.m, which generates synthetic examples of boundary samples (full noisy circle, a projected circle that becomes an ellipse, and a partial arc with outliers) and demonstrates three fitting methods:

• Algebraic least-squares circle fit (closed form)
• Geometric nonlinear ciarcle fit (minimizes radial residuals with fminsearch)
• Direct least-squares ellipse fit (Fitzgibbon et al.)

Why this is useful

• When tracking pupils from eye-camera data, edge detectors often return boundary points at the transition between sclera and pupil. Depending on viewing angle and occlusion, those points form noisy circles, ellipses, or partial arcs. This demo helps compare fit quality and decide when to use a circle vs an ellipse.

Requirements

• MATLAB R2016b or newer
• No special toolboxes required (uses fminsearch available in base MATLAB)

How to run

1. Open MATLAB and change to the folder containing the demo:

cd('~username//Documents//MATLAB')
circle_fit_demo

What the script does

• Generates three examples: noisy circle, projected circle (appears as an ellipse), and a partial arc with outliers.
• Fits with the three methods above, plots the samples and fitted curves, and prints numeric summaries to the Command Window.

Example output (sample run):

Circle: true center [5.00 -2.00], r=3.50
 Algebraic fit center [4.990 -2.000], r=3.497
 Geometric fit center [4.990 -2.000], r=3.496

Proj example: geometric circle center [0.001 0.003], r=1.465
 Ellipse center [-0.001 -0.004], a=1.803 b=1.079 angle=22.22 deg

Partial: alg center [-0.274 2.306], r=1.566; geom center [-0.244 2.315], r=1.533

Notes and next steps

• For real video/frame processing, extract boundary points per frame (edge detection, thresholding) and feed them into the fitting functions.
• Add a RANSAC wrapper to be robust to outliers, and temporal filtering / Kalman smoothing for stable pupil center across frames.
• To recover 3D gaze direction from an observed ellipse, fit the ellipse per frame and apply an inverse projection model (requires camera calibration and an eye model).

circle_fit_demo.m

% circle_fit_demo.m
% Demo: fit circles and ellipses to noisy boundary samples
% Run: open and run this script in MATLAB (R2016b or newer)

rng(1);
% Examples to run
examples = {@example_circle, @example_ellipse_projected, @example_partial_and_noisy};

figure('Name','Circle/Ellipse Fit Demo','Position',[100 100 1100 400]);
for k=1:numel(examples)
    subplot(1,3,k);
    examples{k}();
end

%% Example implementations
function example_circle()
    % Perfect circle with Gaussian noise
    cx = 5; cy = -2; r = 3.5;
    t = linspace(0,2*pi,120)';
    pts = [cx + r*cos(t), cy + r*sin(t)];
    sigma = 0.05; % small gaussian noise
    pts = pts + sigma*randn(size(pts));

    [cA, rA] = fit_circle_algebraic(pts);
    [cG, rG] = fit_circle_geometric(pts);

    plot_data_and_fits(pts, cA, rA, cG, rG);
    title('Noisy Circle'); axis equal;
    fprintf('Circle: true center [%.2f %.2f], r=%.2f\n',cx,cy,r);
    fprintf(' Algebraic fit center [%.3f %.3f], r=%.3f\n',cA(1),cA(2),rA);
    fprintf(' Geometric fit center [%.3f %.3f], r=%.3f\n\n',cG(1),cG(2),rG);
end

function example_ellipse_projected()
    % Simulate a 3D circular pupil viewed at an angle -> ellipse
    cx = 0; cy = 0; r = 1.8;
    t = linspace(0,2*pi,140)';
    circle = [r*cos(t), r*sin(t)];
    % Simulate tilt by scaling y axis and rotating
    tilt = 0.6; % foreshortening factor
    A = [1 0; 0 tilt];
    theta = pi/8; R = [cos(theta) -sin(theta); sin(theta) cos(theta)];
    M = R*A; pts = (M*circle')';
    pts = pts + 0.03*randn(size(pts));

    % Fit circle (will be biased) and ellipse
    [cA, rA] = fit_circle_algebraic(pts);
    [cG, rG] = fit_circle_geometric(pts);
    el = fit_ellipse_direct(pts(:,1), pts(:,2));

    plot(pts(:,1), pts(:,2), '.k'); hold on;
    plot_circle(cG, rG, 'b--',1.5);
    plot_ellipse(el, 'r',1.5);
    legend('samples','geom circle','ellipse','Location','best');
    axis equal; title('Projected Circle -> Ellipse');
    fprintf('Proj example: geometric circle center [%.3f %.3f], r=%.3f\n',cG(1),cG(2),rG);
    fprintf(' Ellipse center [%.3f %.3f], a=%.3f b=%.3f angle=%.2f deg\n\n',el.X0_in,el.Y0_in,el.long_axis,el.short_axis,rad2deg(el.phi));
end

function example_partial_and_noisy()
    % Partial arc with outliers and radial jitter (harder case)
    cx = -1; cy = 2; r = 2.2;
    t = linspace(-0.6,1.4,80)';
    pts = [cx + r*cos(t), cy + r*sin(t)];
    % radial jitter and random outliers
    pts = pts + 0.06*randn(size(pts));
    % add some spurious points (outliers)
    pts = [pts; cx+3*(rand(10,1)-0.5), cy+3*(rand(10,1)-0.5)];

    [cA, rA] = fit_circle_algebraic(pts);
    [cG, rG] = fit_circle_geometric(pts);

    plot_data_and_fits(pts, cA, rA, cG, rG);
    title('Partial Arc + Outliers'); axis equal;
    fprintf('Partial: alg center [%.3f %.3f], r=%.3f; geom center [%.3f %.3f], r=%.3f\n\n',cA(1),cA(2),rA,cG(1),cG(2),rG);
end

%% Helper: plotting and fits
function plot_data_and_fits(pts, cA, rA, cG, rG)
    plot(pts(:,1), pts(:,2), '.k'); hold on;
    plot_circle(cA, rA, 'g',1.2);
    plot_circle(cG, rG, 'b--',1.5);
    legend('samples','algebraic','geometric','Location','best'); hold off;
end

function plot_circle(c, r, style, lw)
    th = linspace(0,2*pi,200);
    plot(c(1)+r*cos(th), c(2)+r*sin(th), style,'LineWidth',lw);
end

function plot_ellipse(el, style, lw)
    th = linspace(0,2*pi,300);
    x = el.long_axis*cos(th); y = el.short_axis*sin(th);
    % rotation
    R = [cos(el.phi) -sin(el.phi); sin(el.phi) cos(el.phi)];
    p = R*[x; y];
    plot(p(1,:)+el.X0_in, p(2,:)+el.Y0_in, style,'LineWidth',lw);
end

%% Circle fits
function [center, R] = fit_circle_algebraic(pts)
    x = pts(:,1); y = pts(:,2);
    A = [x, y, ones(size(x))];
    b = -(x.^2 + y.^2);
    params = A\b; % solves in least-squares sense
    a = params(1); bpar = params(2); c = params(3);
    center = -[a/2, bpar/2];
    R = sqrt((a^2 + bpar^2)/4 - c);
end

function [center, R] = fit_circle_geometric(pts)
    % uses fminsearch to minimize radial residuals
    x = pts(:,1); y = pts(:,2);
    % init from algebraic
    [c0, r0] = fit_circle_algebraic(pts);
    p0 = [c0(:)' r0];
    opt = optimset('Display','off');
    p = fminsearch(@(p) sum((sqrt((x-p(1)).^2+(y-p(2)).^2)-p(3)).^2), p0, opt);
    center = p(1:2); R = p(3);
end

%% Ellipse fit (direct least squares by Fitzgibbon et al.)
function ellipse_t = fit_ellipse_direct(x,y)
    % Inputs: column vectors x,y
    x = x(:); y = y(:);
    D = [x.^2, x.*y, y.^2, x, y, ones(size(x))];
    S = D'*D;
    % Constraint matrix C (for conic to be ellipse)
    C = zeros(6); C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;
    % Solve generalized eigenproblem
    [eigvec, eigval] = eig(inv(S)*C);
    eigval = diag(eigval);
    % find positive definite solution (ellipse constraint 4AC-B^2>0)
    cond = 4*eigvec(1,:).*eigvec(3,:) - eigvec(2,:).^2;
    idx = find(cond>0);
    if isempty(idx)
        % fallback: use smallest eigenvalue
        [~,idx] = min(abs(eigval));
    else
        idx = idx(1);
    end
    a = eigvec(:,idx);
    % unpack conic
    A=a(1); B=a(2); Cc=a(3); D=a(4); E=a(5); F=a(6);
    % compute center
    den = B^2 - 4*A*Cc;
    X0 = (2*Cc*D - B*E)/den;
    Y0 = (2*A*E - B*D)/den;
    % major/minor axes and angle
    up = 2*(A*E^2 + Cc*D^2 + F*B^2 - 2*B*D*E - 4*A*Cc*F);
    down1 = (B^2 - 4*A*Cc)*( (Cc - A)*sqrt(1 + (B/(A-Cc)).^2) - (Cc + A) );
    down2 = (B^2 - 4*A*Cc)*( (A - Cc)*sqrt(1 + (B/(A-Cc)).^2) - (Cc + A) );
    a_len = sqrt(up/down1);
    b_len = sqrt(up/down2);
    % angle
    if B==0 && A < Cc
        phi = 0;
    else
        phi = 0.5 * atan2(B, A - Cc);
    end
    % Prepare output
    ellipse_t.X0_in = X0; ellipse_t.Y0_in = Y0;
    ellipse_t.long_axis = max(a_len,b_len); ellipse_t.short_axis = min(a_len,b_len);
    ellipse_t.phi = phi;
    ellipse_t.a = A; ellipse_t.b = B; ellipse_t.c = Cc; ellipse_t.d = D; ellipse_t.e = E; ellipse_t.f = F;
end

prompts.md

Prompt #1

Hi - a problem that is common to solve is automated tracking the pupil of the eye given the a sample of points that a computer vision system will pick out from a video system. The video task will pick out points on the boundary between the white of the eye and the pupil. Ultimately the objective is to determine where the eye is pointed in space.

Can you demonstrate a circular regression model from a set samples of points around a circle with a variety of noise. Some of the circles would be ellipses because of the viewing angle.

Prompt #2

I'd like to see this run in MATLAB - please start it using the MCP server

Prompt #3

Can you give me a good short description of this code and the problem it tries to solve. I'm going to put this in LinkedIn. Something in MarkDown would be great