This forum is disabled, please visit https://forum.opencv.org
 | 1 | initial version |
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;
}
cout << "[a b c] ="<<EVec.col(index)<<endl;
cout << "[d e f] ="<<T*EVec.col(index)<<endl;
return abs(d);
}
| | 2 | No.2 Revision |
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] ="<<EVec.col(index)<<endl;
="<<tmp13 abs(tmp46.at<double="">(2))<<endl; cout="" <<="" "[d="" e="" f]="<<tmp46/abs(tmp46.at<double>(2))<<endl;
}
else
{
double a=tmp13.at<double>(0),b=tmp13.at<double>(1),c=tmp13.at<double>(2);
double d=tmp46.at<double>(0),e=tmp46.at<double>(1),f=tmp46.at<double>(2);</p>
<pre><code>double ac=abs(a+b+c-d-e+f);
cout << "[d e f] ="<<T*EVec.col(index)<<endl;
<< " [a="" b="" c]="<<a/mx/mx/ac<<" \t"<<b="" mx="" my="" ac<<"\t"<<c="" my="" my="" ac="" <<endl;="" cout="" <<="" "[d="" e="" f]="<<(-b/mx+d/mx-2*a/mx)/ac<<" \t"<<(-b="" my+d="" my-2*c="" my)="" ac<<"\t"<<-1<<endl;="" <="" code="">
}
return abs(d);
}
}
| | 3 | No.3 Revision |
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(tmp46.at<double="">(2))<<endl; cout="" <<="" "[d="" e="" f]="<<tmp46/abs(tmp46.at<double>(2))<<endl;
="<<tmp13/abs(tmp46.at<double>(2))<<endl;
cout << "[d e f] ="<<tmp46/abs(tmp46.at<double>(2))<<endl;
}
else
{
double a=tmp13.at<double>(0),b=tmp13.at<double>(1),c=tmp13.at<double>(2);
a=tmp13.at<double>(0),b=tmp13.at<double>(1),c=tmp13.at<double>(2);
double d=tmp46.at<double>(0),e=tmp46.at<double>(1),f=tmp46.at<double>(2);</p>
<pre><code>double d=tmp46.at<double>(0),e=tmp46.at<double>(1),f=tmp46.at<double>(2);
double ac=abs(a+b+c-d-e+f);
cout << " [a="" b="" c]="<<a/mx/mx/ac<<" \t"<<b="" mx="" my="" ac<<"\t"<<c="" my="" my="" ac="" <<endl;="" cout="" <<="" "[d="" e="" f]="<<(-b/mx+d/mx-2*a/mx)/ac<<" \t"<<(-b="" my+d="" my-2*c="" my)="" ac<<"\t"<<-1<<endl;="" <="" code="">
<< "[a b c] ="<<a/mx/mx/ac<<"\t"<<b/mx/my/ac<<"\t"<<c/my/my/ac <<endl;
cout << "[d e f] ="<<(-b/mx+d/mx-2*a/mx)/ac<<"\t"<<(-b/my+d/my-2*c/my)/ac<<"\t"<<-1<<endl;
}
return abs(d);
}