# C++ cv::eigen and matlab eig() return different results 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


Thank you very much for any helpfull suggestions.

edit retag close merge delete

If your M matrices are different, then I guess the title of the question is incorrect. Have you checked that the inverse of C and S3 are the same in both programs?

Yes, C and S3 are acutally the same in both programs. I also splitted the line M=inv(C)(S1-S2inv(S3)*S2') into multiple lines. It's fine until i multiply it with the inverse of C ... but i have no clue why.

Sort by » oldest newest most voted

First you should work with double instead of float. Yes there is a problem. I think problem is in S1, S2, S3 using scilab(matlab clone) I have got this

--> S1
S1  =

1.666D+13   6.640D+12   4.037D+12
6.640D+12   4.037D+12   3.387D+12
4.037D+12   3.387D+12   3.961D+12

--> S2
S2  =

3.988D+10   1.686D+10   1.013D+08
1.686D+10   1.211D+10   48020669.
1.211D+10   1.358D+10   48544327.

--> S3
S3  =

1.013D+08   48020669.   288422.
48020669.   48544327.   192469.
288422.     192469.     1156.


Data are not good for numerical computation. You have to change origin and scale in your data : normalize between -1 +1 It will give you :

Eigenvalue: [-2.118390551250616;
-117.9218332186111;
-122.9426765681663]
Eigenvector: [0.7071067312014415, 0.0001447164181614355, 0.7071068163628255;
-0.007183290097389554, 0.9999498482088568, 0.006978639585240439;
0.707070343765731, 0.01001399641712809, -0.707072308072356]
cond =[1.999845198382039, -0.9998939021741439, -1.999951296207895]
[a b c] =[0.7071067312014415;
-0.007183290097389554;
0.707070343765731]
[d e f] =[-0.001143313356372834;
-5.692680060736142e-005;
-0.6524152428074719]


Parameters are given in new axis. If you want in original system you have to rescale it

with this code

float LeastSquaresEllipseFit(vector<double> xCors, vector<double> yCors){
//Design matrix split into quadratic and linear part
cv::Mat D1(xCors.size(), 3, CV_64FC1); cv::Mat D2(xCors.size(), 3, CV_64FC1);
//Scatter matrix split into quadratic, combined and linear part
cv::Mat S1(3, 3, CV_64FC1); cv::Mat S2(3, 3, CV_64FC1); cv::Mat S3(3, 3, CV_64FC1);

double mx=0,my=0;
for (size_t i = 0; i < xCors.size(); i++)
{
mx += xCors[i];
my += yCors[i];
}
mx /= xCors.size();
my /= xCors.size();
cout << mx << "\t"<<my<<endl;
if (mx!=0)
for (size_t i = 0; i < xCors.size(); i++)
{
xCors[i] = (xCors[i]-mx)/mx;
yCors[i] = (yCors[i]-my)/my;
}

for (size_t i = 0; i < xCors.size(); i++)
{
D1.at<double>(i, 0) = xCors[i] * xCors[i];
D1.at<double>(i, 1) = xCors[i] * yCors[i];
D1.at<double>(i, 2) = yCors[i] * yCors[i];
D2.at<double>(i, 0) = xCors[i];
D2.at<double>(i, 1) = yCors[i];
D2.at<double>(i, 2) = 1;
}

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

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

S1 = D1.t() * D1;
S2 = D1.t() * D2;
S3 = D2.t() * D2;
T=-S3.inv(DECOMP_SVD)*S2.t();
//Solution 1:
M = C.inv() * (S1 + S2*(T));

cv::Mat EVal, EVec;

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

cout << "Eigenvalue: " << EVal << endl;
cout << "Eigenvector: " << EVec << endl;
Mat cond = 4 * EVec.row(0).mul(EVec.row(2)) - EVec.row(1).mul(EVec.row(1));
cout<<"cond ="<<cond<<endl;
int index=-1;
for (int i = 0; i < cond.cols; i++)
{
if (cond.at<double>(i)>0)
index =i;
}
Mat tmp13=EVec.col(index);
Mat tmp46=T*EVec.col(index);
if (mx == 0)
{
cout << "[a b c] ="<<tmp13/abs ...
more

can you pleas post values for S1 S2 S3 T (inv(S3)S2') and S2T and M using matlab? EDIT: I think i acted too fast. I took a closer look now and compared my results with the one of Matlab once again. I know the matrizes of matlab are not normalized. So i think they are not compareable. However i dont get your results either.

Matlab Results:

S1: S2: S3: T: S2*T: M: New C++ Code:

float LeastSquaresEllipseFit(unsigned int *xCors, unsigned int *yCors,
unsigned int startID, unsigned int stopID){
unsigned int length = (stopID - startID) + 1;
double mx = 0, my = 0;
cv::Mat XY(length, 2, CV_64F);
cv::Mat D1(length, 3, CV_64F);
cv::Mat D2(length, 3, CV_64F);
cv::Mat S1(3, 3, CV_64F);
cv::Mat S2(3, 3, CV_64F);
cv::Mat S3(3, 3, CV_64F);

for (unsigned int i = 0; i < length; i++)
{
mx += xCors[i];
my += yCors[i];
}

mx /= length;
my /= length;

if (mx != 0 && my!= 0)
for (unsigned int i = 0; i < length; i++)
{
XY.at<double>(i, 0) = (xCors[i] - mx) / mx;
XY.at<double>(i, 1) = (yCors[i] - my) / my;
}

//Compute design matrix
for (unsigned int i = 0; i < length; i++)
{
D1.at<double>(i, 0) = XY.at<double>(i, 0) * XY.at<double>(i, 0);    //x²
D1.at<double>(i, 1) = XY.at<double>(i, 0) * XY.at<double>(i, 1);    //xy
D1.at<double>(i, 2) = XY.at<double>(i, 1) * XY.at<double>(i, 1);    //y²
D2.at<double>(i, 0) = XY.at<double>(i, 0);                          //x
D2.at<double>(i, 1) = XY.at<double>(i, 1);                          //y
D2.at<double>(i, 2) = 1;
}

cv::Mat T(3, 3, CV_64F);
cv::Mat M(3, 3, CV_64F);
cv::Mat C(3, 3, CV_64F);
cv::Mat EVal, EVec;

//Fill constraint matrix
C.setTo(0);
C.at<double>(0, 2) = 2;
C.at<double>(1, 1) = -1;
C.at<double>(2, 0) = 2;

//Solve Ax² + Bxy + Cy² + Dx + Ey + F = 0
S1 = D1.t() * D1;
S2 = D1.t() * D2;
S3 = D2.t() * D2;
T = -S3.inv(cv::DECOMP_SVD) * S2.t();
M = C.inv() * (S1 + S2*(T));
cv::eigen(M, EVal, EVec);

//Condition: 4ac - b² > 0
cv::Mat cond = 4 * EVec.row(0).mul(EVec.row(2)) - EVec.row(1).mul(EVec.row(1));

#ifdef DEBUG
cout << "S1: " << S1 << endl;
cout << "S2: " << S2 << endl;
cout << "S3: " << S2 << endl;
cout << "T: " << T << endl;
cout << "S2 * T: " << S2 * T << endl;
cout << "M" << M << endl;
cout << "Eigenvalue: " << EVal << endl;
cout << "Eigenvector: " << EVec << endl;
cout << "Condition: " << cond << endl;
#endif

//Find eigenvalue based on condition
unsigned int index;
for (unsigned int i = 0; i < cond.cols; i++)
{
if (cond.at<double>(i) > 0)
index = i;
}

cv::Mat abc = EVec.col(index);
cv::Mat def = T * EVec.col(index);
double a = abc.at<double>(0, 0);
double b = abc.at<double>(1, 0);
double c = abc.at<double>(2, 0);
double d = def.at<double>(0, 0);
double e = def.at<double>(1, 0);
double f = def.at<double>(2, 0);

#ifdef DEBUG
cout << "[a b c] =" << abc << endl;
cout << "[d e f] =" << def << endl;
#endif

//http://math ...
more

Acutally i dont know what A is ... . When i run the matlab code from above i get the following for the coefficients: a = 0.3957; 4.9221e-05; 0.9184; 197.4691; 305.8196; -2.6958e+04 When i run the code in C++ i get: 0.707106...; -0.00255...; 0.707102...; -1.266152...e-5;-1.2737...e-5; 0.65189... I rechcked all matrices and they are similar in Matlab and C++. Only the computed eigenvalue and eigenvector deviates.

I modify y answer to have matlab results. May be it's better fit than matlab. last coefficient (f in paper) is -1. Divide matlab answer by -2.6958e+04 and you can compare.

About matrix now with image value cannot be compared easily. if you want to investigate more write matlab value in file. condition matrice for S3 is bad

ok, i have used your solution. Do you have the following values?

abc [5.627e-05; -4.4906e-05; 5.110e-05]
def [-0.0124905; -0.0059619; -1]


no

249.5   166.496
S1 = [320.2876262431536, 0.003645860382940951, 117.7085706879299;
0.003645860382940951, 117.7085706879299, -0.0205616053943787;
117.7085706879299, -0.0205616053943787, 426.9798963942378]
S2 = [-0.0004317050533206857, -0.002324935496871142, 470.5604716446925;
-0.002324935496828849, -0.008002495716534254, -0.008341223021411146;
-0.008002495716534698, -0.008425705572179287, 595.1871952006982]
S3 = [470.5604716446925, -0.008341223021411146, 7.638334409421077e-014;
-0.008341223021411146, 595.1871952006982, 6.548928066507642e-013;
7.638334409421077e-014, 6.548928066507642e-013, 1156]


and

M= [-62.28393582258224, -0.008133565471408984, 60.26858244189097;
-0.007041198823250678, -117.7085705086586, 0.01626713094281796;
64.37082112960803, 0.003520599411625336, -62.28393582258224]
Eigenvalue: [-2.015353094784753;
-117.7085093378262;
-122.5525797212122]
Eigenvector: [0.7071067752811925, 4.971163862746497e-005, 0.7071067853444671;
-0.002553782180504198, 0.9999936552424071, 0.002483479717073073;
0.7071021754654753, 0.003561882042342094, -0.7071024158128075]
cond =[1.999980434547407, -0.9999866022570985, -1.999993832290309]
[a b c] =1.495147e-005  -8.091897e-008  3.357502e-005
[d e f] =-7.447244e-003 -1.115990e-002  -1


Okay i got nearly the same results in C++ and Matlab now. Here are my final results for C++:

[a b c] =1.495147e-005  -8.091897e-008  3.357502e-005
[d e f] =-7.447244e-003 -1.115990e-002  -1


... and the results from Matlab:

[a b c] = 1.4678e-005  1.8257e-009  3.4064e-005
[d e f] = -7.3e-003 -1.13e-002  0.9999


Actually i am not certain why b is negative in C++ and positive in Matlab. Same applies for f.

Unfortunatly i can't post it in one comment. The new matlab code is here:

Xmean = mean(x);   % the centroid of the data set
Ymean = mean(y);

D1 = [(x(:)-Xmean).^2, (x(:)-Xmean).*(y(:)-Ymean),...
(y(:)-Ymean).^2];
D2 = [x(:)-Xmean, y(:)-Ymean, ones(size(x,1),1)];
S1 = D1'*D1;
S2 = D1'*D2;
S3 = D2'*D2;
T = -inv(S3)*S2';
M = S1 + S2*T;
M = [M(3,:)./2; -M(2,:); M(1,:)./2];
[evec,eval] = eig(M);
cond = 4*evec(1,:).*evec(3,:)-evec(2,:).^2;
A1 = evec(:,find(cond>0));
A = [A1; T*A1];
A4 = A(4)-2*A(1)*Xmean-A(2)*Ymean;
A5 = A(5)-2*A(3)*Ymean-A(2)*Xmean;
A6 = A(6)+A(1)*Xmean^2+A(3)*Ymean^2+...
A(2)*Xmean*Ymean-A(4)*Xmean-A(5)*Ymean;
A(4) = A4;  A(5) = A5;  A(6) = A6;
A = A/norm(A);


You can edit your previous question and post matlab code in. About -8e--8 and 1.8257e-009 diference is very low relative to 1 (Did you forget a minus before 0.999? in matlab result)

I had too less characters available. I computed the ellipse center from c++ code and it seems to work. Thank you very much for all the effort.

Official site

GitHub

Wiki

Documentation