#include <Geometry/ProjectionParametersSpherical.hh>
#include <Geometry/Projection.hh>
#include "Base/Math/Random.hh"
#include "MathAlgo/PolynomialSolve.hh"
using namespace BIAS;
using namespace std;
int main(int argc, char *argv[]) {
bool testDistortion = true;
bool skipLinearMonom = !testDistortion;
if (argc<2) {
cout<<"error, no xml file given."<<endl;
return -2;
}
if (P.
Load(argv[1])!=0) {
cout<<"error, could not load xml."<<endl;
return -2;
}
if (ppS == NULL) {
cout<<"error, no spherical parameters."<<endl;
return -2;
}
unsigned int width, height;
cout<<
"maxcamangle of loaded projection is "<< ppS->
GetMaxCamAngle()<<endl;
int degree = 4;
std::vector<double> coefficients;
coefficients.resize(degree+1);
cout<<endl<<endl<<endl<<"now estimating polynomial"<<endl;
if (testDistortion) {
cout<<"distortion function is used ..."<<endl;
} else {
cout<<"undistortion function is used ..."<<endl;
}
cout<<"coefficients are ";
for (unsigned int i=0; i<coefficients.size(); i++) {
cout<<coefficients[i]<<" ";
}
cout<<endl;
double px=0, py=0;
cout<<endl<<endl<<endl<<"now initing with polynomial"<<endl;
coefficients,
!testDistortion);
cout<<endl<<endl<<endl<<"now estimating polynomial again "<<endl;
if (testDistortion) {
cout<<"distortion function is used ..."<<endl;
} else {
cout<<"undistortion function is used ..."<<endl;
}
cout<<"coefficients are ";
for (unsigned int i=0; i<coefficients.size(); i++) {
cout<<coefficients[i]<<" ";
}
cout<<endl;
cout<<"generating some 3D points with original projection .."<<endl;
vector<HomgPoint3D> points3D;
vector<HomgPoint2D> points2D;
for (unsigned int i=0; i<20; i++) {
if (X.NormL2()>0.1) {
points2D.push_back(x);
points3D.push_back(X);
}
}
double averageError = 0;
for (unsigned int i=0; i<points3D.size(); i++) {
ProjectionRay.Normalize();
if (testDistortion) {
if(
Equal(rPhiTheta[0], 0.0)) {
BIASERR("sphere coordinate transform failed for ray "
<< ProjectionRay);
}
double radius = PS.EvaluatePolynomial(rPhiTheta[2], coefficients);
result[0] = radius * cos(rPhiTheta[1]) + px;
result[1] = radius * sin(rPhiTheta[1]) + py;
cout<<"projects to "<<result<<", measured at "<<points2D[i]
<<" diff is "<<(result-points2D[i]).NormL2()<<endl;
averageError += (result-points2D[i]).NormL2();
cout<<"Reconstructed projects to "<<result<<", measured at "<<points2D[i]
<<" diff is "<<(result-points2D[i]).NormL2()<<endl;
} else {
points2D[i][1]-py,
0);
double radius = sqrt(ImageRay[0]*ImageRay[0] + ImageRay[1]*ImageRay[1]);
ImageRay[2] = -1.0*PS.EvaluatePolynomial(radius, coefficients);
ImageRay.Normalize();
if (ImageRay.ScalarProduct(ProjectionRay)<0.0)
ImageRay *= -1.0;
cout<<"ray error for manual polynomial projection in degree: ";
averageError += acos(ImageRay.ScalarProduct(ProjectionRay))*180.0/M_PI;
cout<<acos(ImageRay.ScalarProduct(ProjectionRay))*180.0/M_PI
<<" degree"<<endl;
cout<<"ray error for reconstructed interpolator projection in degree: ";
cout<<acos(ImageRay.ScalarProduct(pRay))*180.0/M_PI<<" degree"<<endl;
}
testpoint2 = points2D[i];
cout<<"ray error (inconsistency) for original project/unproject: ";
<<" degree"<<endl;
cout<<endl;
}
averageError /= points3D.size();
cout<<"Average angular error (manual) is "<<averageError <<endl;
return 0;
}