Revision history [back]

Hello,

i am trying to implement an ellipse fitting algorithm according to the following paper: Link Paper

I went trough every matrix step by step. I could observe the first major differences at this part of the code:

M = C.inv() * (S1 - S2(S3.inv()S2.t()));

Acutally, I don't know why. From my point of view the code is similar to the Matlab script. Going on, by computing the eigenvalue and the eigenvector i observed completly different results (see below).

Eigenvalue in MatLab:

Eigenvector in MatLab:

This are the results for my C++ implementation:

Could someone tell me why this is happening? Am i missing something or did I do something wrong?

In the following my implementation and the MatLab implementation:

C++ Code:

float LeastSquaresEllipseFit(unsigned int *xCors, unsigned int *yCors,
unsigned int startID, unsigned int stopID){
    unsigned int length = (stopID - startID) + 1;
//Design matrix split into quadratic and linear part
cv::Mat D1(length, 3, CV_32F); cv::Mat D2(length, 3, CV_32F);
//Scatter matrix split into quadratic, combined and linear part
cv::Mat S1(3, 3, CV_32F); cv::Mat S2(3, 3, CV_32F); cv::Mat S3(3, 3, CV_32F);

vector<unsigned int> X;
vector<unsigned int> Y;

for (size_t i = 0; i < length; i++)
{
    X.push_back(xCors[i]);
    Y.push_back(yCors[i]);
}

for (unsigned int i = 0; i < length; i++)
{
    D1.at<float>(i, 0) = xCors[i] * xCors[i];
    D1.at<float>(i, 1) = xCors[i] * yCors[i];
    D1.at<float>(i, 2) = yCors[i] * yCors[i];
    D2.at<float>(i, 0) = xCors[i];
    D2.at<float>(i, 1) = yCors[i];
    D2.at<float>(i, 2) = 1;
}

cv::Mat T(3, 3, CV_32F);
cv::Mat M(3, 3, CV_32F);
cv::Mat C(3, 3, CV_32F); C.setTo(0);

C.at<float>(0, 2) = 2;
C.at<float>(1, 1) = -1;
C.at<float>(2, 0) = 2;

S1 = D1.t() * D1;
S2 = D1.t() * D2;
S3 = D2.t() * D2;

//cout << C << endl;

//Solution 1:
M = C.inv() * (S1 - S2*(S3.inv()*S2.t()));

//cout << "M " << M << endl;

cv::Mat EVal, EVec;

cv::eigen(M, EVal, EVec);

cout << "Eigenvalue: " << EVal << endl;
cout << "Eigenvector: " << EVec << endl;

float d = 0;
return abs(d);

MatLab Code:

function a = fit_ellipse(x, y)
for i = 1:1:length(x)
    D1(i,1:3) = [x(i).^2, x(i)*y(i), y(i).^2];
    D2(i,1:3) = [x(i), y(i), 1];        
end
C = [0,0,2;0 -1, 0; 2,0,0];
S1 = D1' * D1;                  
S2 = D1' * D2;                  
S3 = D2' * D2;
M = inv(C)*(S1-S2*inv(S3)*S2');
[evec, eval] = eig(M)
cond = 4 * evec(1, :) .* evec(3, :) - evec(2, :) .^2; % evaluate a’Ca
a1 = evec(:, find(cond > 0));   % eigenvector for min. pos. eigenvalue
a = [a1; T * a1];               % ellipse coefficients

Here are the datapoints: http://cloud.directupload.net/hqaW

Thank you very much for any helpfull suggestions.