Using inbuilt RANSAC function for circle fit in 2D data

66 views (last 30 days)
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?
The function is more or less described on Slide 21 but i do not know how to implement it?
Thanks for help!

Accepted Answer

Image Analyst
Image Analyst on 9 May 2019
  2 Comments
Matlab User
Matlab User on 9 May 2019
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
Image Analyst on 3 Oct 2020
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.

Sign in to comment.

More Answers (5)

Image Analyst
Image Analyst on 12 Feb 2021
@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
  1 Comment
Said
Said on 16 Feb 2021
Hello ;
did you manage to modify the code and making working for noisy circles as well ?
Thanks
Said

Sign in to comment.


Said
Said on 11 Feb 2021
Edited: Said on 11 Feb 2021
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
  2 Comments
Said
Said on 11 Feb 2021
Edited: Said on 11 Feb 2021
for the fitting just apply any fitting algo; there are a lot
you can use this one for example
Regards
Said
Image Analyst
Image Analyst on 12 Feb 2021
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.

Sign in to comment.


Said
Said on 11 Feb 2021
Edited: Said on 11 Feb 2021
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
  1 Comment
Image Analyst
Image Analyst on 11 Feb 2021
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

Sign in to comment.


Said
Said on 11 Feb 2021
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
  1 Comment
Image Analyst
Image Analyst on 11 Feb 2021
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?

Sign in to comment.


Said
Said on 11 Feb 2021
sorry here is the png

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!