Using inbuilt RANSAC function for circle fit in 2D data

I would like to fit a circle with a predifend radius r to a 2D dataset using the inbuild RANSAC function.

My Dataset is for example:

data=[
    0.0135    0.2333
    0.0133    0.2336
    0.0036    0.2345
    0.0034    0.2343
    0.0230    0.2347
    0.0231    0.2348
    0.0130    0.2337
    0.0134    0.2336
    0.0032    0.2342
    0.0229    0.2342
   -0.0050    0.2375
    0.0035    0.2344
    0.0233    0.2352
   -0.0058    0.2382
    0.0111    0.2361
    0.0321    0.2382
    0.0325    0.2384
    0.0014    0.2367
   -0.0051    0.2382
    0.0235    0.2365
    0.0060    0.2402
    0.0130    0.2413
   -0.0096    0.2410
   -0.0054    0.2380
   -0.0100    0.2410
    0.0321    0.2379
    0.0360    0.2410
    0.0132    0.2349
    0.0320    0.2381
    0.0363    0.2411
   -0.0088    0.2410
   -0.0056    0.2381
    0.0037    0.2362
    0.0232    0.2368
   -0.0101    0.2409
    0.0272    0.2424
    0.0353    0.2407
    0.0368    0.2407
    0.0332    0.2381
   -0.0046    0.2388
   -0.0163    0.2445
   -0.0091    0.2416
    0.0347    0.2410
   -0.0166    0.2446
   -0.0124    0.2379
    0.0424    0.2444
    0.0229    0.2405
    0.0428    0.2446
    0.0391    0.2391
    0.0306    0.2386
   -0.0168    0.2445
    0.0396    0.2383
    0.0423    0.2433
   -0.0084    0.2414
   -0.0166    0.2449
    0.0344    0.2414
    0.0129    0.2360
    0.0425    0.2446
    0.0033    0.2372
   -0.0205    0.2506
    0.0227    0.2374
   -0.0221    0.2490
   -0.0171    0.2453
   -0.0220    0.2489
   -0.0221    0.2503
    0.0294    0.2387
   -0.0221    0.2485
    0.0468    0.2509
   -0.0223    0.2509
    0.0275    0.2401
   -0.0039    0.2386
   -0.0009    0.2400
   -0.0223    0.2499
   -0.0165    0.2457
   -0.0220    0.2488
   -0.0213    0.2501
    0.0424    0.2453
    0.0450    0.2458
   -0.0085    0.2422
    0.0333    0.2420
   -0.0227    0.2504
   -0.0224    0.2494
   -0.0227    0.2509
   -0.0226    0.2495
   -0.0052    0.2391
   -0.0194    0.2538
    0.0467    0.2504
    0.0293    0.2383
   -0.0207    0.2503
   -0.0162    0.2466
    0.0418    0.2466
   -0.0243    0.2526
   -0.0201    0.2503
   -0.0228    0.2492
    0.0463    0.2507
   -0.0248    0.2525];

%Plot it
plot(data(:,1),data(:,2),'o')
axis([-0.05 0.1, 0.2 0.35])
grid on

Now i would like to use the inbuild RANSAC function

[model,inlierIdx] = ransac(data,fitFcn,distFcn,sampleSize,maxDistance)

OK data is clear, sampleSize = 3, as a circle requires minimum 3 points to be defined, maxDistance can be changed depending on the noise level (how noise the data are).

How do I need to define the fitFcn and the distFcn?

In addition i would like do predefine the radius of the circel for example with 0.05 for my given data above?

A hint could by the lectur note from Thomas Opsahl available via: https://www.uio.no/studier/emner/matnat/its/UNIK4690/v16/forelesninger/lecture_3_3-robust-estimation-with-ransac.pdf

The function is more or less described on Slide 21 but i do not know how to implement it?

Thanks for help!

[Szymon Bogdani]

Have you found anything regarding this problem? Im in the same situation.

Answer 1

[Image Analyst]

Why RANSAC? Why not use the FAQ: https://matlab.fandom.com/wiki/FAQ#How_can_I_fit_a_circle_to_a_set_of_XY_data.3F

[Matlab User]

Because a least-squares circle fitting treats all points equally, so a random noise-point shifts the estimated circle away from the desired solution: compare slide 17 with 23 from Thomas Opsahl lecture notes and I often have such noise outliers.

[Image Analyst]

If you want to fit a circle to minimize the residuals perpendicular to the fitted curve, then it becomes a lot more complicated. Here's the paper on that (attached). I don't have MATLAB code for it but if anyone makes some, it would be nice if you uploaded it to the File Exchange, or back here.

Answer 2

[Image Analyst]

@Said, thanks. That's outstanding. I'm attaching a full demo (using noisy but "good" data that did not need RANSAC) in case anyone is interested:

% Demo to fit an ellipse through a set of noisy (x,y) coordinates.
% Initialization / cleanup steps:
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 22;
fprintf('Beginning to run %s.m ...\n', mfilename);

%----------------------------------------------------------------------------------------------------
% Create sample noisy data.
numPoints = 1000;
angles = linspace(0, 360, numPoints);
% Make ellipse centered at origin with x radius of 15 and y radius of 6.
x = 15 * cosd(angles) + rand(1, numPoints);
y = 6 * sind(angles) + rand(1, numPoints);
% Rotate the ellipse that is centered at the origin by 30 degrees.  Ref: https://en.wikipedia.org/wiki/Rotation_matrix
thetaDegrees = 30;
rotationMatrix = [cosd(thetaDegrees), -sind(thetaDegrees); sind(thetaDegrees), cosd(thetaDegrees)];
xy = [x(:), y(:)];
% Now do the rotation by doing a matrix multiply of the data by the rotation matrix.
xyRotated = xy * rotationMatrix;
x = xyRotated(:, 1);
y = xyRotated(:, 2);
% Now offset to put the center at x = 35, y = 75.
Cx = 35;
Cy = 75;
x = x + Cx;
y = y + Cy;
% Plot input data:
plot(x, y, 'b.');
xline(Cx, 'Color', 'g', 'LineWidth', 2); % Draw line through center.
yline(Cy, 'Color', 'g', 'LineWidth', 2); % Draw line through center.
grid on;
axis equal;
title('Input Data', 'FontSize', fontSize);
% Make into NumPoints-by-2 tall array so we can call the fitellipse() function.
data = [x(:), y(:)];

%----------------------------------------------------------------------------------------------------
% Do the fit:
% FITELLIPSE  Least-squares fit of ellipse to 2D points.
%        A = FITELLIPSE(X,Y) returns the parameters of the best-fit
%        ellipse to 2D points (X,Y).
%        The returned vector A contains the center, radii, and orientation
%        of the ellipse, stored as (Cx, Cy, Rx, Ry, theta_radians)
%
% Authors: Andrew Fitzgibbon, Maurizio Pilu, Bob Fisher
% Reference: "Direct Least Squares Fitting of Ellipses", IEEE T-PAMI, 1999
fittedEllipseParameters = fitellipse(x, y)
%----------------------------------------------------------------------------------------------------

% Prepare for construction of the fitted ellipse by extracting the parameters from the returned argument.
CxFit = fittedEllipseParameters(1)
CyFit = fittedEllipseParameters(2)
Rx = fittedEllipseParameters(3)
Ry = fittedEllipseParameters(4)
theta_radians = fittedEllipseParameters(5)
% Negate the angle because the function gives you the negative of the actual angle for some reason.
thetaDegrees = -rad2deg(theta_radians) % Negate and convert to degrees.
% First make the ellipse that is NOT rotated.
xFit = Rx * cosd(angles);
yFit = Ry * sind(angles);
% Compute the rotation matrix.  Ref: https://en.wikipedia.org/wiki/Rotation_matrix
rotationMatrix = [cosd(thetaDegrees), -sind(thetaDegrees); sind(thetaDegrees), cosd(thetaDegrees)];
xy = [xFit(:), yFit(:)]; % Load x and y into N-by-2 matrix.
% Do the rotation by doing a matrix multiply of the data by the rotation matrix.
xyRotated = xy * rotationMatrix;
% Extract back out the rotated x and y.
xFit = xyRotated(:, 1);
yFit = xyRotated(:, 2);
% Now offset to put the center at where it should be.
xFit = xFit + CxFit;
yFit = yFit + CyFit;

% Now plot the fitted results:
hold on;
plot(xFit, yFit, 'r-', 'LineWidth', 3);
darkGreen = [0, 0.5, 0];
xline(Cx, 'Color', darkGreen, 'LineWidth', 2); % Draw line through center.
yline(Cy, 'Color', darkGreen, 'LineWidth', 2); % Draw line through center.
grid on;
title('Ellipse Fitted Through Noisy Data', 'FontSize', fontSize);
legend('Data', 'Cx', 'Cy', 'Fit', 'Fit Cx', 'Fit Cy', 'Location', 'southwest');
fprintf('Done running %s.m\n', mfilename);
% END OF DEMO SCRIPT.


%=========================================================================================================================================================
% Reference: http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/FITZGIBBON/ELLIPSE/
%  Copyright (c) 1999, 2005, Andrew Fitzgibbon, Maurizio Pilu, Bob Fisher
%
% Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
function a = fitellipse(X,Y)

% FITELLIPSE  Least-squares fit of ellipse to 2D points.
%        A = FITELLIPSE(X,Y) returns the parameters of the best-fit
%        ellipse to 2D points (X,Y).
%        The returned vector A contains the center, radii, and orientation
%        of the ellipse, stored as (Cx, Cy, Rx, Ry, theta_radians)
%
% Authors: Andrew Fitzgibbon, Maurizio Pilu, Bob Fisher
% Reference: "Direct Least Squares Fitting of Ellipses", IEEE T-PAMI, 1999
%
% This is a more bulletproof version than that in the paper, incorporating
% scaling to reduce roundoff error, correction of behaviour when the input
% data are on a perfect hyperbola, and returns the geometric parameters
% of the ellipse, rather than the coefficients of the quadratic form.
%
%  Example:  Run fitellipse without any arguments to get a demo
if nargin == 0
	% Create an ellipse
	t = linspace(0,2);
	
	Rx = 300
	Ry = 200
	Cx = 250
	Cy = 150
	Rotation = .4 % Radians
	
	x = Rx * cos(t);
	y = Ry * sin(t);
	nx = x*cos(Rotation)-y*sin(Rotation) + Cx;
	ny = x*sin(Rotation)+y*cos(Rotation) + Cy;
	% Draw it
	plot(nx,ny,'o');
	% Fit it
	fitellipse(nx,ny)
	% Note it returns (Rotation - pi/2) and swapped radii, this is fine.
	return
end

% normalize data
mx = mean(X);
my = mean(Y);
sx = (max(X)-min(X))/2;
sy = (max(Y)-min(Y))/2;

x = (X-mx)/sx;
y = (Y-my)/sy;

% Force to column vectors
x = x(:);
y = y(:);

% Build design matrix
D = [ x.*x  x.*y  y.*y  x  y  ones(size(x)) ];

% Build scatter matrix
S = D'*D;

% Build 6x6 constraint matrix
C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;

% Solve eigensystem
[gevec, geval] = eig(S,C);

% Find the negative eigenvalue
I = find(real(diag(geval)) < 1e-8 & ~isinf(diag(geval)));

% Extract eigenvector corresponding to negative eigenvalue
A = real(gevec(:,I));

% unnormalize
par = [
	A(1)*sy*sy,   ...
	A(2)*sx*sy,   ...
	A(3)*sx*sx,   ...
	-2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy,   ...
	-A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy,   ...
	A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my   ...
	- A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my   ...
	+ A(6)*sx*sx*sy*sy   ...
	]';

% Convert to geometric radii, and centers

thetarad = 0.5*atan2(par(2),par(1) - par(3));
cost = cos(thetarad);
sint = sin(thetarad);
sin_squared = sint.*sint;
cos_squared = cost.*cost;
cos_sin = sint .* cost;

Ao = par(6);
Au =   par(4) .* cost + par(5) .* sint;
Av = - par(4) .* sint + par(5) .* cost;
Auu = par(1) .* cos_squared + par(3) .* sin_squared + par(2) .* cos_sin;
Avv = par(1) .* sin_squared + par(3) .* cos_squared - par(2) .* cos_sin;

% ROTATED = [Ao Au Av Auu Avv]

tuCentre = - Au./(2.*Auu);
tvCentre = - Av./(2.*Avv);
wCentre = Ao - Auu.*tuCentre.*tuCentre - Avv.*tvCentre.*tvCentre;

uCentre = tuCentre .* cost - tvCentre .* sint;
vCentre = tuCentre .* sint + tvCentre .* cost;

Ru = -wCentre./Auu;
Rv = -wCentre./Avv;

Ru = sqrt(abs(Ru)).*sign(Ru);
Rv = sqrt(abs(Rv)).*sign(Rv);

a = [uCentre, vCentre, Ru, Rv, thetarad];
end

[Said]

Hello ;

did you manage to modify the code and making working for noisy circles as well ?

Thanks

Said

Answer 3

[Said]

Ah no; RANSAC is not supposed to fit your data; Ransac gives you only the best candidate points for fitting the data. That why when you have outliers you have to apply RANSAC first before any ellipse fitting algorithm.

here is the fitted data

[Said]

for the fitting just apply any fitting algo; there are a lot

you can use this one for example

http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/FITZGIBBON/ELLIPSE/

Regards

Said

[Image Analyst]

Oh, okay. I thought it was something like fitPolynomialRANSAC() in the Computer Vision Toolbox which has the fitting part built in to it. Thanks for the link to the ellipse fitting code. It's useful. And if your data is good to start with you won't even need to run RANSAC on it first.

Answer 4

[Said]

Attached the RANSAC for Ellipse

you can adapt it to a circle data model as circle is just special case of ellipse where the semi-axis are equal

Regards

Said

[Image Analyst]

OK, so what did I do wrong?

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 22;
fprintf('Beginning to run %s.m ...\n', mfilename);

%----------------------------------------------------------------------------------------------------
% Create sample noisy data.
numPoints = 1000;
angles = linspace(0, 360, numPoints);
x = 15 * cosd(angles) + rand(1, numPoints);
y = 5 * sind(angles) + rand(1, numPoints);
% Plot input data:
subplot(2, 1, 1);
plot(x, y, 'b.');
grid on;
axis equal;
title('Input Data', 'FontSize', fontSize);
% Make into NumPoints-by-2 tall array.
data = [x(:), y(:)];

%----------------------------------------------------------------------------------------------------
% Do the RANSAC fit:
filterData = ellipseDataFilter_RANSAC(data)
% Plot results:
subplot(2, 1, 2);
plot(filterData(:, 1), filterData(:, 2), 'r-');
grid on;
title('Fitted Data', 'FontSize', fontSize);
fprintf('Done running %s.m\n', mfilename);

%===================================================================================================
function filterData=ellipseDataFilter_RANSAC(data)
% Do ellipse scatter data filtering for ellipse fitting by RANSAC method.
% Author: Zhenyu Yuan
% Date: 2016/7/26
% Ref:  http://www.cnblogs.com/yingying0907/archive/2012/10/13/2722149.html
%       Extract RANSAC filtering in ellipsefit.m and make some modification
% Parameter initialization
nSampLen = 3;               % Set the number of points on which the model is based
nDataLen = size(data, 1);   % Data length
nIter = 50;                 % Maximum number of iterations
dThreshold = 2;             % Threshold
nMaxInlyerCount=-1;         % Minimum points
A=zeros([2 1]);
B=zeros([2 1]);
P=zeros([2 1]);
%  Main loop
for i = 1:nIter
	ind = ceil(nDataLen .* rand(1, nSampLen)); %Sampling, select nSampLen different points
	%  Build the model, store the coordinate points, the focal point and a point across the ellipse needed for modeling
	%Ellipse definition equation: the sum of the distances to two fixed points is constant
	A(:,1)=data(ind(1),:);    % focus
	B(:,1)=data(ind(2),:);    % focus
	P(:,1)=data(ind(3),:);    % Point on ellipse
	DIST=sqrt((P(1,1)-A(1,1)).^2+(P(2,1)-A(2,1)).^2)+sqrt((P(1,1)-B(1,1)).^2+(P(2,1)-B(2,1)).^2);
	xx=[];
	nCurInlyerCount=0;        % The number of initial points is 0
	%  Does it fit the model?
	for k=1:nDataLen
		%         CurModel=[A(1,1)   A(2,1)  B(1,1)  B(2,1)  DIST ];
		pdist=sqrt((data(k,1)-A(1,1)).^2+(data(k,2)-A(2,1)).^2)+sqrt((data(k,1)-B(1,1)).^2+(data(k,2)-B(2,1)).^2);
		CurMask =(abs(DIST-pdist)< dThreshold);     % The point whose distance to the straight line is less than the threshold
		%conforms to the model and is marked as 1
		nCurInlyerCount =nCurInlyerCount+CurMask;   % Count the number of points conforming to the ellipse model
		if(CurMask==1)
			xx =[xx;data(k,:)];
		end
	end
	% Choose the best model
	if nCurInlyerCount > nMaxInlyerCount   % The model with the most points that fits the model is the best model
		nMaxInlyerCount = nCurInlyerCount;
		%             Ellipse_mask = CurMask;
		%              Ellipse_model = CurModel;
		%              Ellipse_points = [A B P];
		filterData =xx;
	end
end
end

Answer 5

[Said]

you have done nothing wrong

your data is an ellipse so Ransac filter found most of it well fitted with an ellipse.

just use equal axis to plot the data against the filtred one as I did here in the attached figure

Said

[Image Analyst]

I can't open your figure. It just launches a second instance of MATLAB but doesn't show the figure. Please just post a PNG screenshot of your data and fitted ellipse like I did.

But with my code, the output fitted data is every bit as noisy as the input training data. It's not a perfectly smooth ellipse so I must have done something wrong. Can you post your complete demo code like I did?

Answer 6

[Said]

sorry here is the png