Ask Your Question
0

C++ cv::eigen and matlab eig() return different results

asked 2016-08-12 10:43:56 -0500

Ice_T02 gravatar image

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:

image description

Eigenvector in MatLab:

image description

This are the results for my C++ implementation:

image description

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.

edit retag flag offensive close merge delete

Comments

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?

LorenaGdL gravatar imageLorenaGdL ( 2016-08-12 13:15:48 -0500 )edit

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.

Ice_T02 gravatar imageIce_T02 ( 2016-08-12 15:53:25 -0500 )edit

2 answers

Sort by » oldest newest most voted
2

answered 2016-08-13 10:05:05 -0500

LBerger gravatar image

updated 2016-08-16 12:47:02 -0500

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)
edit flag offensive delete link more

Comments

Thank you very much for your answer.

Ice_T02 gravatar imageIce_T02 ( 2016-08-16 03:01:09 -0500 )edit

can you pleas post values for S1 S2 S3 T (inv(S3)S2') and S2T and M using matlab?

LBerger gravatar imageLBerger ( 2016-08-16 04:44:50 -0500 )edit

Find all the matrices in my second post.

Ice_T02 gravatar imageIce_T02 ( 2016-08-16 08:13:04 -0500 )edit
0

answered 2016-08-16 07:52:48 -0500

Ice_T02 gravatar image

updated 2016-08-16 07:53:27 -0500

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: image description

S2: image description

S3: image description

T: image description

S2*T: image description

M: image description

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)
edit flag offensive delete link more

Comments

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.

Ice_T02 gravatar imageIce_T02 ( 2016-08-16 11:51:32 -0500 )edit

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

LBerger gravatar imageLBerger ( 2016-08-16 12:59:37 -0500 )edit

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]
Ice_T02 gravatar imageIce_T02 ( 2016-08-17 02:32:46 -0500 )edit

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]
LBerger gravatar imageLBerger ( 2016-08-17 02:48:48 -0500 )edit

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
LBerger gravatar imageLBerger ( 2016-08-17 02:49:18 -0500 )edit

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.

Ice_T02 gravatar imageIce_T02 ( 2016-08-17 06:24:37 -0500 )edit

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);
Ice_T02 gravatar imageIce_T02 ( 2016-08-17 06:27:16 -0500 )edit

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)

LBerger gravatar imageLBerger ( 2016-08-17 12:05:21 -0500 )edit

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.

Ice_T02 gravatar imageIce_T02 ( 2016-08-17 13:26:38 -0500 )edit
Login/Signup to Answer

Question Tools

1 follower

Stats

Asked: 2016-08-12 10:43:56 -0500

Seen: 1,130 times

Last updated: Aug 16 '16