From 70bd11350373f6ffd2f18133f8da89acd5e5cb0a Mon Sep 17 00:00:00 2001 From: Laurent Belcour Date: Thu, 29 Aug 2013 18:28:07 +0200 Subject: [PATCH] Updating the parametrization code --- sources/core/params.cpp | 21 ++++--- sources/core/params.h | 52 ++++++++--------- sources/core/rational_function.cpp | 16 +----- sources/softs/data2moments/main.cpp | 85 +++++++++++++++++++--------- sources/softs/generate_data/main.cpp | 7 ++- sources/softs/tests/main.cpp | 53 ++++++++++++----- 6 files changed, 143 insertions(+), 91 deletions(-) diff --git a/sources/core/params.cpp b/sources/core/params.cpp index 8649c82..b061025 100644 --- a/sources/core/params.cpp +++ b/sources/core/params.cpp @@ -162,7 +162,7 @@ void params::from_cartesian(const double* invec, params::input outtype, // 3D Parametrization case params::RUSIN_TH_PH_TD: outvec[0] = acos(half[2]); - outvec[1] = atan2(half[0], half[1]); + outvec[1] = atan2(half[1], half[0]); outvec[2] = acos(half[0]*invec[0] + half[1]*invec[1] + half[2]*invec[2]); break; case params::RUSIN_TH_TD_PD: @@ -172,11 +172,18 @@ void params::from_cartesian(const double* invec, params::input outtype, diff[0] = invec[0]; diff[1] = invec[1]; diff[2] = invec[2]; - rotate_normal(diff, -outvec[1]); - rotate_binormal(diff, -outvec[0]); + rotate_normal(diff, -atan2(half[1], half[0])); +#ifdef DEBUG + std::cout << "diff* = [ " << diff[0] << ", " << diff[1] << ", " << diff[2] << "]" << std::endl; +#endif + rotate_binormal(diff, -outvec[0]); +#ifdef DEBUG + std::cout << "half = [ " << half[0] << ", " << half[1] << ", " << half[2] << "]" << std::endl; + std::cout << "diff = [ " << diff[0] << ", " << diff[1] << ", " << diff[2] << "]" << std::endl; +#endif outvec[1] = acos(diff[2]); - outvec[2] = atan2(diff[0], diff[1]); + outvec[2] = atan2(diff[1], diff[0]); break; case params::ISOTROPIC_TV_TL_DPHI: outvec[0] = acos(invec[2]); @@ -193,8 +200,8 @@ void params::from_cartesian(const double* invec, params::input outtype, diff[0] = invec[0]; diff[1] = invec[1]; diff[2] = invec[2]; - rotate_normal(diff, -outvec[1]); rotate_binormal(diff, -outvec[0]); + rotate_normal(diff, -outvec[1]); outvec[2] = acos(diff[2]); outvec[3] = atan2(diff[0], diff[1]); @@ -202,9 +209,9 @@ void params::from_cartesian(const double* invec, params::input outtype, case params::SPHERICAL_TL_PL_TV_PV: outvec[0] = acos(invec[2]); - outvec[1] = atan2(invec[0], invec[1]); + outvec[1] = atan2(invec[1], invec[0]); outvec[2] = acos(invec[5]); - outvec[3] = atan2(invec[3], invec[4]); + outvec[3] = atan2(invec[4], invec[3]); #ifdef DEBUG std::cout << invec[2] << " - acos -> " << outvec[0] << std::endl; #endif diff --git a/sources/core/params.h b/sources/core/params.h index 93d8c81..6a01dc0 100644 --- a/sources/core/params.h +++ b/sources/core/params.h @@ -15,6 +15,12 @@ * \details * Any \typedef function object or \typedef data object should have * an associated parametrization. + * + * We use the following convention to defined the tangent, normal and + * bi-normal of the surface: + * * The normal is the upper vector (0, 0, 1) + * * The tangent direction is along x direction (1, 0, 0) + * * The bi-normal is along the y direction (0, 1, 0) */ class params { @@ -194,37 +200,31 @@ class params out[5] = cos(theta_v); } - //! \brief rotate a cartesian vector with respect to the normal of - //! theta degrees. - static void rotate_normal(double* vec, double theta) - { - const double cost = cos(theta); - const double sint = sin(theta); + //! \brief rotate a cartesian vector with respect to the normal of + //! theta degrees. + static void rotate_normal(double* vec, double theta) + { + const double cost = cos(theta); + const double sint = sin(theta); - const double temp = cost * vec[0] + sint * vec[1]; + const double temp = cost * vec[0] + sint * vec[1]; - vec[1] = cost * vec[1] - sint * vec[0]; - vec[0] = temp; - } + vec[1] = cost * vec[1] - sint * vec[0]; + vec[0] = temp; + } - //! \brief rotate a cartesian vector with respect to the bi-normal of - //! theta degrees. - static void rotate_binormal(double* vec, double theta) - { - const double cost = cos(theta); - const double sint = sin(theta); + //! \brief rotate a cartesian vector with respect to the bi-normal of + //! theta degrees. + static void rotate_binormal(double* vec, double theta) + { + const double cost = cos(theta); + const double sint = sin(theta); - const double temp = cost * vec[1] + sint * vec[2]; + const double temp = cost * vec[0] + sint * vec[2]; -#ifdef DEBUG - std::cout << acos(vec[2]) << std::endl; -#endif - vec[2] = cost * vec[2] - sint * vec[1]; -#ifdef DEBUG - std::cout << acos(vec[2]) << std::endl; -#endif - vec[1] = temp; - } + vec[2] = cost * vec[2] - sint * vec[0]; + vec[0] = temp; + } static void print_input_params(); diff --git a/sources/core/rational_function.cpp b/sources/core/rational_function.cpp index b026402..9a089c6 100644 --- a/sources/core/rational_function.cpp +++ b/sources/core/rational_function.cpp @@ -219,7 +219,8 @@ double rational_function_1d::p(const vec& x, int i) const { #ifdef POLYNOMIALS res *= pow(x[k], deg[k]) ; -#else // LEGENDRE +// res *= pow(2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5), deg[k]) ; +#else res *= legendre(2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5), deg[k]); #endif } @@ -228,18 +229,7 @@ double rational_function_1d::p(const vec& x, int i) const } double rational_function_1d::q(const vec& x, int i) const { - std::vector deg = index2degree(i); - double res = 1.0; - for(int k=0; kparametrization(); - int d_size = params::dimension(data_param); + const int nX = d->dimX(); + const int nY = d->dimY(); double in_angle[4] = {0.0, 0.0, 0.0, 0.0} ; - // Sample every degree - double dtheta = 0.5*M_PI / 90.0; - // Static data vec X(3); X[0] = 1; X[1] = 0; X[2] = 0; vec Y(3); Y[0] = 0; Y[1] = 1; Y[2] = 0; vec N(3); N[0] = 0; N[1] = 0; N[2] = 1; - // Moments - vec rawm0(d->dimY()); - vec rawm1(d->dimY()); for(int theta_in=0; theta_in<90; theta_in++) { in_angle[0] = theta_in * 0.5*M_PI / 90.0; + vec m_0(nY); m_0 = vec::Zero(nY); + vec m_x(nY); m_0 = vec::Zero(nY); + vec m_y(nY); m_0 = vec::Zero(nY); + vec m_xx(nY); m_0 = vec::Zero(nY); + vec m_xy(nY); m_0 = vec::Zero(nY); + vec m_yy(nY); m_0 = vec::Zero(nY); + // Integrate over the light hemisphere for(int theta_out=0; theta_out<90; theta_out++) { in_angle[2] = theta_out * 0.5*M_PI / 90.0; - for(int phi_out=0; phi_out<180; phi_out++) + for(int phi_out=0; phi_out<360; phi_out++) { - in_angle[3] = phi_out * 0.5*M_PI / 180.0; + in_angle[3] = phi_out * 2.0*M_PI / 360.0; - vec in(d_size), cart(6), L(3), V(3); + vec in(nX), cart(6), L(3), V(3); params::convert(in_angle, params::SPHERICAL_TL_PL_TV_PV, data_param, &in[0]); params::convert(in_angle, params::SPHERICAL_TL_PL_TV_PV, params::CARTESIAN, &cart[0]); L[0] = cart[0]; @@ -123,33 +125,62 @@ int main(int argc, char** argv) // Copy the input vector vec x = d->value(in); - // Get the projected 2D coordinate - vec xy = coord(V, L, X, Y, N); + // Get the projected 2D coordinate + vec xy = coord(V, L, X, Y, N); - for(int i=0; idimY(); ++i) + for(int i=0; idimY(); ++i) +/* + for(int i=0; i> rusin 3d: " << rusi << std::endl; std::cout << "<> cartesian: " << spherical << std::endl; nb_tests_failed++; - } - std::cout << "<> rusin 3d: " << rusi << std::endl; - std::cout << "<> cartesian: " << spherical << std::endl; - std::cout << std::endl; + } // Convert issue #1 vec cart2(6); spherical[0] = 0; spherical[1] = 0; spherical[2] = 1.51844; spherical[3] = -2.96706; params::convert(&spherical[0], params::SPHERICAL_TL_PL_TV_PV, params::CARTESIAN, &cart2[0]); - std::cout << "<> spherical before: " << spherical << std::endl; try { params::convert(&spherical[0], params::SPHERICAL_TL_PL_TV_PV, params::RUSIN_TH_TD_PD, &rusi[0]); @@ -122,16 +118,43 @@ int parametrization_tests() std::cout << "<> rusin 3d: " << rusi << std::endl; std::cout << "<> spherical: " << spherical << std::endl; nb_tests_failed++; - } - std::cout << "<> rusin after: " << rusi << std::endl; - std::cout << "<> spherical after: " << spherical << std::endl; - std::cout << "<> " << cart << " / " << cart2 << std::endl; - { - double dot = cart[0]*cart[3] + cart[1]*cart[4] + cart[2]*cart[5]; - double dot2 = cart2[0]*cart2[3] + cart2[1]*cart2[4] + cart2[2]*cart2[5]; - std::cout << "<> dot = " << dot << ", " << "dot2 = " << dot2 << std::endl; - } - std::cout << std::endl; + } + + + + /// Test on a known couple of directions + /// in = [0, 0, 1] out = [1, 0, 0] + + // Convert from Cartesian to spherical + cart[0] = 0; cart[1] = 0; cart[2] = 1; + cart[3] = 1; cart[4] = 0; cart[5] = 0; + params::convert(&cart[0], params::CARTESIAN, params::SPHERICAL_TL_PL_TV_PV, &spherical[0]); + std::cout << "From cartesian to spherical" << std::endl; + std::cout << spherical << std::endl << std::endl; + + params::convert(&cart[0], params::CARTESIAN, params::RUSIN_TH_TD_PD, &rusi[0]); + std::cout << "From cartesian to rusi" << std::endl; + std::cout << rusi << std::endl << std::endl; + + params::convert(&rusi[0], params::RUSIN_TH_TD_PD, params::CARTESIAN, &cart[0]); + std::cout << "From rusi to cartesian" << std::endl; + std::cout << cart << std::endl << std::endl; + + params::convert(&spherical[0], params::SPHERICAL_TL_PL_TV_PV, params::RUSIN_TH_TD_PD, &rusi[0]); + std::cout << "From spherical to rusi" << std::endl; + std::cout << rusi << std::endl << std::endl; + + params::convert(&rusi[0], params::RUSIN_TH_TD_PD, params::SPHERICAL_TL_PL_TV_PV, &spherical[0]); + std::cout << "From rusi to spherical" << std::endl; + std::cout << spherical << std::endl << std::endl; + + params::convert(&rusi[0], params::RUSIN_TH_TD_PD, params::CARTESIAN, &cart[0]); + std::cout << "From rusi to cartesian" << std::endl; + std::cout << cart << std::endl << std::endl; + + params::convert(&spherical[0], params::SPHERICAL_TL_PL_TV_PV, params::RUSIN_TH_TD_PD, &rusi[0]); + std::cout << "From spherical to rusi" << std::endl; + std::cout << rusi << std::endl << std::endl; return nb_tests_failed; } -- GitLab