CloudRunner

... share your algorithms as a web service

Algorithm: Ellipse Fitting

Description:
Simulation of different Ellipse Fitting algorithms:
* Least-Squares with algebraic distance
* Least-Squares with orthogonal distance
* Gradient-Weighted Least-Squares
* M-Estimator (Cauchy)
* Least Median of Squares
Paper: Zhengyou Zhang. Parameter estimation techniques: a tutorial with application to conic fitting. Image and Vision Computing, 15(1):59 – 76, 1997.
Tags: Ellipse Fitting Least-Squares
Usage: Algorithm is public.
Viewed 4487 times, called 222 times
Upload:
Empty star Empty star Empty star Empty star Empty star
0 votes
Sebastian Dingler
02/17/2015 1:18 a.m. (version #2)

Select Version

Version #2: 02/17/2015 1:18 a.m.
Changelog:
Cleaned the code a bit.

Version #1: 02/04/2015 12:35 p.m.
No changelog available

Run Algorithm

: x center of ellipse
: y center of ellipse
: Major axis
: Minor axis
: Angle between major axis a and x axis
: Variance of measurement noise
: Interval of measurement sampling (must be between [0,2*pi])
: Number of measurements
: Number of outlier
: Bool to switch between simulation with outlier or without
: Allow cached result?
Close

Please Wait

Computation is running...

Plots

showing plot # of

Result Value(s)


		

Resource Usage

Execution of the algorithm took seconds.
Peak memory usage for running the algorithm was kilobytes (total usage: kilobytes).

Run ID / Link

ID of this run:
Link:

Result from Cache

Result has been delivered from cache and computed on (UTC).

Matlab log

(show)

		

Error

Using this algorithm in your local MATLAB environment is easy: Click here for instructions!

Usage Instructions for CloudRunner Client

  1. Download the CloudRunner Client by clicking here and place the downloaded file in your MATLAB working directory.

  2. Inside MATLAB, initialize the CloudRunner Client by calling CloudRunner:
    >> CloudRunner

    A login dialog will prompt for your CloudRunner mail address and password. For a start, you can leave the dialog empty and just click "Connect".

    Alternatively, you can provide the login credentials (or empty strings to skip login) as a parameter and hence skip the login dialog. This is useful when using CloudRunner in non-interactive scripts.
    >> CloudRunner('mail@example.com', 'password')

  3. Select this algorithm by its URL. Selecting an algorithm creates the lambda function that proxies calls to the algorithm to the server for execution:
    >> CloudRunnerSelect('http://www.cloudrunner.eu/algorithm/172/ellipse-fitting/version/2/')

    For the sake of convenience, you can also use the algorithm ID instead of its URL for public algorithms.

  4. Call functions from the algorithm like any regular local function.

Note: You can find further information on the help page.

Source Code

File:

 1 function r = XY2xy(X, alpha, Xc)
 2 %XY2XY Transforms ellipse that it is symmetric to x- and y-axis
 3 %
 4 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 5 %        Karlsruhe Institute of Technology (KIT), Germany
 6 %
 7 % DATE   22.12.2014
 8 R = [cos(alpha) sin(alpha);-sin(alpha) cos(alpha)];
 9 r = R*(X-Xc);
10 end
 1 function r = xy2XY2(x, alpha, Xc)
 2 %XY2XY2 Rotates ellipse back (inverse function of XY2xy.m)
 3 %
 4 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 5 %        Karlsruhe Institute of Technology (KIT), Germany
 6 %
 7 % DATE   22.12.2014
 8 R = [cos(alpha) sin(alpha);-sin(alpha) cos(alpha)];
 9 r = pinv(R)*x+Xc;
10 end
 1 function [ output_args ] = simFitEllipseOutlier(  X_c,Y_c,a,b,alpha,var,sampleInterval,nMeasurements,nOutlier)
 2 %SIMFITELLIPSEOUTLINER Simulates the fitting of an ellipse with outlier
 3 %
 4 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 5 %        Karlsruhe Institute of Technology (KIT), Germany
 6 %
 7 % DATE   22.12.2014
 8 r=plotEllipse(X_c, Y_c, a, b, rad2deg(alpha), 100,'black',3);
 9 hold on;
10 noise=sampleEllipseWithOutlier(X_c,Y_c,a,b,alpha,var,sampleInterval,nMeasurements,nOutlier);
11 plot(noise(1,:),noise(2,:),'O','Color','black');
12 hold on;
13 % Ellipse fitting with Least-Squares based on orthogonal distance
14 rORTHO=fitWithLSO(noise(1,:),noise(2,:));
15 r=plotEllipse(rORTHO(1), rORTHO(2), rORTHO(3), rORTHO(4), rad2deg(rORTHO(5)), 100,'red',1);
16 hold on;
17 % Ellipse fitting with Least-Squares based on algebraic distance
18 p=fitWithLSA(noise(1,:),noise(2,:));
19 [x0,y0,a,b,angle]=Matrix2AngleForm(p(1),p(2),1-p(1),p(3),p(4),p(5));
20 r=plotEllipse(x0, y0, a, b, rad2deg(angle), 100,'green',1);
21 hold on;
22 % Ellipse fitting with Weighted Least-Squares based on algebraic distance
23 p=fitWithGWLS(noise(1,:),noise(2,:));
24 [x0,y0,a,b,angle]=Matrix2AngleForm(p(1),p(2),1-p(1),p(3),p(4),p(5));
25 r=plotEllipse(x0, y0, a, b, rad2deg(angle), 100,'blue',1);
26 hold on;
27 % Ellipse fitting with M-Estimator using Cauchy loss function
28 p=fitWithMEstimator(noise(1,:),noise(2,:));
29 [x0,y0,a,b,angle]=Matrix2AngleForm(p(1),p(2),1-p(1),p(3),p(4),p(5));
30 r=plotEllipse(x0, y0, a, b, rad2deg(angle), 100,'yellow',1);
31 hold on;
32 % Ellipse fitting with Least Median of Squares
33 p=fitWithLSoS(noise,15);
34 r=plotEllipse(p(1), p(2), p(3), p(4), rad2deg(p(5)), 100,'black',1);
35 legend('GroundTruth','Noise','LSO','LSA','LSWA','M-Estimator','LSoM')
36 ylim([-30 30]);
37 end
 1 function [ output_args ] = simFitEllipse( X_c,Y_c,a,b,alpha,var,sampleInterval,nMeasurements)
 2 %SAMPLEELLIPSE Simulates the fitting of an ellipse without outlier
 3 %
 4 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 5 %        Karlsruhe Institute of Technology (KIT), Germany
 6 %
 7 % DATE   22.12.2014
 8 r = plotEllipse(X_c, Y_c, a, b, rad2deg(alpha), 100,'black',4);
 9 hold on;
10 noise = sampleEllipse(X_c,Y_c,a,b,alpha,var,sampleInterval,nMeasurements);
11 plot(noise(1,:),noise(2,:),'O','Color','black');
12 hold on;
13 % Ellipse fitting with Least-Squares based on orthogonal distance
14 rORTHO = fitWithLSO(noise(1,:),noise(2,:));
15 r = plotEllipse(rORTHO(1), rORTHO(2), rORTHO(3), rORTHO(4), rad2deg(rORTHO(5)), 100,'red',1);
16 hold on;
17 % Ellipse fitting with Least-Squares based on algebraic distance
18 p = fitWithLSA(noise(1,:),noise(2,:));
19 [x0,y0,a,b,alpha] = Matrix2AngleForm(p(1),p(2),1-p(1),p(3),p(4),p(5));
20 r = plotEllipse(x0, y0, a, b, rad2deg(alpha), 100,'green',1);
21 hold on;
22 % Ellipse fitting with Weighted Least-Squares based on algebraic distance
23 p = fitWithGWLS(noise(1,:),noise(2,:));
24 [x0,y0,a,b,alpha] = Matrix2AngleForm(p(1),p(2),1-p(1),p(3),p(4),p(5));
25 r = plotEllipse(x0, y0, a, b, rad2deg(alpha), 100,'yellow',1);
26 
27 legend('GroundTruth','Noise','L-S with orthogonal distance','L-S with algebraic distance','Gradient Weighted L-S')
28 ylim([-30 30]);
29 end
 1 function [ r ] = sampleEllipseWithOutlier(X_c, Y_c, a, b, alpha, var, sampleInterval, nSample, nOutlier)
 2 %SAMPLEELLIPSEWITHOUTLINER Samples a ellipse with outliers
 3 % Samples a ellipse of the form:
 4 % x(t)=x_c + a * cos(t) * cos(alpha) - b * sin(t) * sin(alpha) + var*randn;
 5 % y(t)=y_c + a * cos(t) * sin(alpha) + b * sin(t) * cos(alpha) + var*randn;
 6 % 
 7 % Plus adds uniformly distributed outliers 
 8 %
 9 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
10 %        Karlsruhe Institute of Technology (KIT), Germany
11 %
12 % DATE   22.12.2014
13 [XY] = sampleEllipse(X_c,Y_c,a,b,alpha,var,sampleInterval,nSample);
14 % add uniformly distributed outliers 
15 X = [XY(1,:) (a*2) * (rand(1,nOutlier)-0.5)+X_c];
16 Y = [XY(2,:) (a*2) * (rand(1,nOutlier)-0.5)+Y_c];
17 r = [X;Y];
18 end
 1 function [r] = sampleEllipse(x_c, y_c, a, b, phi, var, sampleInterval, nSamples)
 2 %SAMPLEELLIPSE Samples a ellipse
 3 % Samples a ellipse of the form:
 4 % x(t)=x_c + a * cos(t) * cos(alpha) - b * sin(t) * sin(alpha) + var*randn;
 5 % y(t)=y_c + a * cos(t) * sin(alpha) + b * sin(t) * cos(alpha) + var*randn;
 6 % 
 7 %
 8 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 9 %        Karlsruhe Institute of Technology (KIT), Germany
10 %
11 % DATE   22.12.2014
12 t=linspace(sampleInterval(1),sampleInterval(2),nSamples);
13 for i=1:length(t)
14     x(i)=x_c+a*cos(t(i))*cos(phi)-b*sin(t(i))*sin(phi)+sqrt(var)*randn;
15     y(i)=y_c+a*cos(t(i))*sin(phi)+b*sin(t(i))*cos(phi)+sqrt(var)*randn;
16 end
17 r=[x;y];
18 end
The file is not a text file or too large to be displayed here. You can still download the file using the link below.
 1 function p = plotEllipse( x, y, a, b, angle, steps, lColor, lWidth)
 2 %SAMPLEELLIPSE Plots a ellipse
 3 %
 4 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 5 %        Karlsruhe Institute of Technology (KIT), Germany
 6 %
 7 % DATE   22.12.2014
 8 p = getPointsOnEllipse(x, y, a, b, angle, steps);
 9 plot(p(:,1), p(:,2),'Color',lColor,'LineWidth',lWidth);
10 axis equal;
11 end
 1 function [x0, y0, a, b, angle] = Matrix2AngleForm(A, B, C, D, E, F)
 2 %MATRIX2ANGLEFORM Computes the angle form of matrix form of an ellipse
 3 %  See: http://en.wikipedia.org/wiki/Matrix_representation_of_conic_sections#Reduced_equation
 4 %
 5 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 6 %        Karlsruhe Institute of Technology (KIT), Germany
 7 %
 8 % DATE   22.12.2014
 9 
10 B = 2*B;
11 D = 2*D;
12 E = 2*E;
13 AQ = [A B/2 D/2;B/2 C E/2;D/2 E/2 F];
14 A33 = [A B/2;B/2 C];
15 % Center
16 center = pinv(A33)*[-D/2;-E/2];
17 x0 = center(1);
18 y0 = center(2);
19 % Major and minor axis
20 e = eig(A33);
21 [V,D] = eig(A33);
22 a = sqrt(-det(AQ)/(det(A33)*e(1)));
23 b = sqrt(-det(AQ)/(det(A33)*e(2)));
24 % angle
25 angle = acos(-dot(V(:,1),[1;0])/norm(V(:,1)));
26 end
 1 % Simulation of different ellipse fitting algorithms
 2 %
 3 % Input:
 4 % x_c            - x center of ellipse
 5 % y_c            - y center of ellipse
 6 % a              - major axis a
 7 % b              - minor axis b
 8 % alpha          - angle of ellipse
 9 % var            - variance of measurement noise
10 % sampleInterval - interval of measurements must be between [0;2*pi]
11 % nMeasurements  - # of measurements
12 % nOutlier       - # of outlier
13 % withOutlier    - bool switch for simulation with or withour outlier
14 %
15 % Example:
16 % >> main(0,0,24,12,0,0.25,[0.5*pi 1.5*pi],200,30,0)
17 %
18 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
19 %        Karlsruhe Institute of Technology (KIT), Germany
20 %
21 % DATE   22.12.2014
22 function main(x_c, y_c, a, b, alpha, var, sampleInterval, nMeasurements, nOutlier, withOutlier)
23 if withOutlier
24     simFitEllipseOutlier(x_c,y_c,a,b,alpha,var,sampleInterval, nMeasurements, nOutlier);
25 else
26     simFitEllipse(x_c,y_c,a,b,alpha,var,sampleInterval, nMeasurements);
27 end
28 end
 1 function r = getPointsOnEllipse(x, y, a, b, angle, steps)
 2 %GETPOINTSONELLIPSE Draws points from an ellipse
 3 %
 4 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 5 %        Karlsruhe Institute of Technology (KIT), Germany
 6 %
 7 % DATE   22.12.2014
 8 
 9 beta = angle * (pi / 180);
10 sinbeta = sin(beta);
11 cosbeta = cos(beta);
12 
13 alpha = linspace(0, 360, steps)' .* (pi / 180);
14 sinalpha = sin(alpha);
15 cosalpha = cos(alpha);
16 
17 X = x + (a * cosalpha * cosbeta - b * sinalpha * sinbeta);
18 Y = y + (a * cosalpha * sinbeta + b * sinalpha * cosbeta);
19 
20 r = [X Y];
21 end
 1 function r = getOrthoPoint(xi, yi, a, b)
 2 %GETORTHOPOINT Compute orthogonal Point on ellipse
 3 %
 4 % [1] Sung Joon Ahn, W. Rauh, and M. Recknagel, "Ellipse fitting and
 5 % parameter assessment of circular object targets for robot vision",
 6 % Intelligent Robots and Systems, 1999.
 7 %
 8 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 9 %        Karlsruhe Institute of Technology (KIT), Germany
10 %
11 % DATE   22.12.2014
12 
13 %Orthogonal contacting point on ellipse
14 xk1 = ([xi;yi]*a*b)/sqrt(b^2*xi^2+a^2*yi^2);
15 if abs(xi)<a
16     xk2=[xi;sign(yi)*(b/a)*sqrt(a^2-xi^2)];
17 else
18     xk2=[sign(xi)*a;0];
19 end
20 
21 x_e = 0.5*(xk1+xk2); % x_0
22 for i=1:4
23     x = x_e(1);
24     y = x_e(2);
25     Qk = [b^2*x a^2*y;(a^2-b^2)*y+b^2*yi (a^2-b^2)*x-a^2*xi];
26     fk = [0.5*(a^2*y^2+b^2*x^2-a^2*b^2);b^2*x*(yi-y)-a^2*y*(xi-x)];
27     x_e = x_e-pinv(Qk)*fk;
28 end
29 r = x_e;
30 end
 1 function p = fitWithMEstimator(x, y)
 2 %FITWITHMESTIMATOR Ellipse fitting with a M-Estimor
 3 %
 4 % [1] Zhengyou Zhang, "Parameter estimation techniques: a tutorial with
 5 % application to conic fitting",
 6 % Image and Vision Computing, 15(1):59 - 76, 1997.
 7 %
 8 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 9 %        Karlsruhe Institute of Technology (KIT), Germany
10 %
11 % DATE   22.12.2014
12 X = [];
13 y2 = [];
14 for i = 1 : length(x)
15     X = [X;x(i)^2-y(i)^2 2*x(i)*y(i) 2*x(i) 2*y(i) 1];
16     y2 = [y2;-y(i)^2];
17 end
18 % M-Estimator
19 p=robustfit(X,y2,'cauchy','','off');
20 end
 1 function r = fitWithLSoS(measurements, p2)
 2 %FITWITHLSOS Ellipse fitting with Least Median of Square
 3 %
 4 % [1] Zhengyou Zhang, "Parameter estimation techniques: a tutorial with
 5 % application to conic fitting",
 6 % Image and Vision Computing, 15(1):59 - 76, 1997.
 7 %
 8 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 9 %        Karlsruhe Institute of Technology (KIT), Germany
10 %
11 % DATE   22.12.2014
12 
13 % 1. Draw m random subsamples of p different points
14 m=round(length(measurements)/p2);
15 p=[];
16 
17 for i = 1 : m
18     J=datasample(measurements',p2);
19     % 2. For each subsample use LSA for fitting a conic
20     p(i,:) = fitWithLSA(J(:,1),J(:,2));
21     [x0,y0,a,b,alpha]=Matrix2AngleForm(p(i,1),p(i,2),1-p(i,1),p(i,3),p(i,4),p(i,5));
22     % 3. Compute the residual for all measurements
23     for j = 1 : length(measurements)
24         xy = XY2xy([measurements(1,j);measurements(2,j)], alpha, [x0;y0]);
25         x_new = getOrthoPoint(xy(1), xy(2), a, b);
26         X_new = xy2XY2(x_new, alpha, [x0;y0]);
27         r(j) = norm([measurements(1,j);measurements(2,j)]-X_new);
28     end
29     M(i) = median(r.^2);
30 end
31 % Take median
32 [v,idx] = min(M);
33 [x0,y0,a,b,alpha] = Matrix2AngleForm(p(idx,1),p(idx,2),1-p(idx,1),p(idx,3),p(idx,4),p(idx,5));
34 r = [x0;y0;a;b;alpha];
35 end
 1 function fit = fitWithLSO( Xi, Yi)
 2 %FITWITHLSO Ellipse fitting with Least-Squares based on orthogonal distance
 3 %
 4 % [1] Sung Joon Ahn, W. Rauh, and M. Recknagel, "Ellipse fitting and
 5 % parameter assessment of circular object targets for robot vision",
 6 % Intelligent Robots and Systems, 1999.
 7 %
 8 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 9 %        Karlsruhe Institute of Technology (KIT), Germany
10 %
11 % DATE   22.12.2014
12 
13 % initial guess
14 p = fitWithLSA(Xi,Yi);
15 [Xc,Yc,a,b,alpha] = Matrix2AngleForm(p(1),p(2),1-p(1),p(3),p(4),p(5));
16 
17 % step size
18 lambda=0.1;
19 
20 for k=1:50
21     J = [];
22     X_new2 = [];
23     for i=1:length(Xi)
24         r = XY2xy([Xi(i);Yi(i)],alpha,[Xc;Yc]);
25         x_new = getOrthoPoint(r(1),r(2),a,b);
26         X_new = xy2XY2(x_new, alpha, [Xc;Yc]);
27         xi = XY2xy([Xi(i);Yi(i)],alpha,[Xc;Yc]);
28         J = [J;calcJacobianMatrix(a,b,x_new(1),x_new(2),alpha,xi(1),xi(2))];
29         X_new2 = [X_new2; [Xi(i);Yi(i)]-X_new];
30     end
31     r=-pinv(J) * X_new2;
32     % update 
33     Xc=Xc-lambda*r(1);
34     Yc=Yc-lambda*r(2);
35     a=a-lambda*r(3);
36     b=b-lambda*r(4);
37     alpha=alpha-lambda*r(5);
38 end
39 
40 fit = [Xc,Yc,a,b,alpha];
41 end
42 
43 % Jacobian matrix at the orthogonal contacting point on ellipse
44 function r = calcJacobianMatrix(a, b, x, y, alpha, xi, yi)
45 C = cos(alpha);
46 S = sin(alpha);
47 R = [C S;-S C];
48 B1 = [b^2 * x * C-a^2 * y * S; b^2 * (yi-y) * C+a^2 * (xi-x) * S];
49 B2 = [b^2 * x * S+a^2 * y * C;b^2 * (yi-y) * S-a^2 * (xi-x) * C];
50 B3 = [a * (b^2-y^2); 2 * a * y * (xi-x)];
51 B4 = [b * (a^2-x^2);-2 * b * x * (yi-y)];
52 B5 = [(a^2-b^2) * x * y; (a^2-b^2) * (x^2 - y^2 - x * xi + y * yi)];
53 B = [B1 B2 B3 B4 B5];
54 Qk = [b^2 * x a^2*y; (a^2-b^2)*y+b^2*yi (a^2-b^2)*x-a^2*xi];
55 r = pinv(R)*pinv(Qk)*B;
56 end
 1 function [p, X, y2] = fitWithLSA(x, y)
 2 %FITWITHLSA Ellipse Fitting with Least-Squares based on algebraic distance
 3 % 
 4 % [1] Zhengyou Zhang, "Parameter estimation techniques: a tutorial with 
 5 % application to conic fitting", 
 6 % Image and Vision Computing, 15(1):59 - 76, 1997.
 7 %
 8 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 9 %        Karlsruhe Institute of Technology (KIT), Germany
10 %
11 % DATE   22.12.2014
12 
13 X = [];
14 y2 = [];
15 for i = 1 : length(x)
16     X = [X;x(i)^2-y(i)^2 2*x(i)*y(i) 2*x(i) 2*y(i) 1];
17     y2 = [y2;-y(i)^2];
18 end
19 % Least-Square
20 p=pinv(X' * X) * X' * y2;
21 end
 1 function p = fitWithGWLS(x, y)
 2 %FITWITHGWLS Ellipse fitting with Gradient Weighted Least-Squares
 3 % 
 4 % [1] Zhengyou Zhang, "Parameter estimation techniques: a tutorial with 
 5 % application to conic fitting", 
 6 % Image and Vision Computing, 15(1):59 - 76, 1997.
 7 %
 8 % AUTHOR Sebastian Dingler <s.dingler@gmail.com>
 9 %        Karlsruhe Institute of Technology (KIT), Germany
10 %
11 % DATE   22.12.2014
12 
13 % initial guess
14 [p0,X,y2] = fitWithLSA(x,y);
15 
16 
17 W = zeros(length(x));
18 p = p0;
19 for k = 1 : 10
20     for i = 1 : length(x)
21         % Gradient
22         partialX = 2 * p(1) * x(i) + 2 * p(2) * y(i) + 2 * p(3);
23         partialY = -2 * p(1) * y(i) + 2 * p(2) * x(i) + 2 * p(4) + 2 * y(i);
24         W(i,i)=partialX^2+partialY^2;
25     end
26     % Weighted Least-Squares
27     p_new=pinv(X' * pinv(W) * X) * X' * pinv(W) * y2;
28     
29     % stop criterium
30     if norm(p-p_new)<0.1
31         break
32     else
33         p=p_new;
34     end
35 end
Download algorithm (17 files) as ZIP

Comments

Please login to post a comment.