Skip to content

GitLab

A
alta
Commit 70bd1135 authored by Laurent Belcour's avatar Laurent Belcour
Browse files
Options

Updating the parametrization code

parent b7b79d8f
Show whitespace changes
Inline Side-by-side
Showing with 143 additions and 91 deletions
sources/core/params.cpp
View file @ 70bd1135
......@@ -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
......
sources/core/params.h
View file @ 70bd1135
......@@ -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
{
......@@ -214,16 +220,10 @@ class params
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();
......
sources/core/rational_function.cpp
View file @ 70bd1135
......@@ -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<int> deg = index2degree(i);
double res = 1.0;
for(int k=0; k<dimX(); ++k)
{
#ifdef POLYNOMIALS
res *= pow(x[k], deg[k]) ;
#else // LEGENDRE
res *= legendre(2.0*((x[k] - _min[k]) / (_max[k]-_min[k]) - 0.5), deg[k]);
#endif
}
return res ;
return p(x, i);
}
// Overload the function operator
......
sources/softs/data2moments/main.cpp
View file @ 70bd1135
......@@ -82,35 +82,37 @@ int main(int argc, char** argv)
{
// Data parametrization
params::input data_param = d->parametrization();
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];
......@@ -126,30 +128,59 @@ int main(int argc, char** argv)
// Get the projected 2D coordinate
vec xy = coord(V, L, X, Y, N);
for(int i=0; i<d->dimY(); ++i)
for(int i=0; i<nY; ++i)
{
double val = x[i] * cos(in_angle[2]);
rawm0[i] += val;
rawm1[i] += val * xy[0];
const double normalization = 360.0*90.0;
double val = x[i] * cos(in_angle[2]) / normalization;
m_0[i] += val;
m_x[i] += val * xy[0];
m_y[i] += val * xy[1];
m_xx[i] += val * xy[0] * xy[0];
m_xy[i] += val * xy[0] * xy[1];
m_yy[i] += val * xy[1] * xy[1];
}
}
}
for(int i=0; i<d->dimY(); ++i)
/*
for(int i=0; i<nY; ++i)
{
rawm0[i] /= 180.0*90.0;
rawm1[i] /= 180.0*90.0 * rawm0[i];
m_x[i] /= m_0[i];
m_y[i] /= m_0[i];
m_xx[i] /= m_0[i];
m_xy[i] /= m_0[i];
m_yy[i] /= m_0[i];
}
*/
// compute cumulants
vec k_x = m_x;
vec k_y = m_y;
vec k_xx = m_xx;// - product(m_x,m_x);
vec k_xy = m_xy;// - product(m_x,m_y);
vec k_yy = m_yy;// - product(m_y,m_y);
// Output the value into the file
file << in_angle[0] << "\t";
for(int i=0; i<rawm0.size(); ++i)
file << rawm0[i] << "\t";
for(int i=0; i<nY; ++i)
file << m_0[i] << "\t";
for(int i=0; i<nY; ++i)
file << k_x[i] << "\t";
for(int i=0; i<nY; ++i)
file << k_y[i] << "\t";
for(int i=0; i<nY; ++i)
file << k_xx[i] << "\t";
for(int i=0; i<nY; ++i)
file << k_xy[i] << "\t";
for(int i=0; i<nY; ++i)
file << k_yy[i] << "\t";
for(int i=0; i<rawm1.size(); ++i)
file << rawm1[i] << "\t";
file << std::endl;
}
......
sources/softs/generate_data/main.cpp
View file @ 70bd1135
......@@ -30,9 +30,10 @@ int main(int argc, char** argv)
//f << "#VS 2" << std::endl;
for(int i=0; i<nbx; ++i)
{
const float x = i / (float)nbx ;
//const float y = 100.0f * exp(-10.0 * x*x) * x*x - 0.01 *x*x*x + 0.1 ;
const float y = (1.0) / (1.0E-10 + x);
const float x = 10.0 * i / (float)nbx ;
const float xp = (x - 9);
const float y = 1000.0f * exp(- x*x) * x*x + 00.1 * exp(-100.0 * xp*xp) *x*x*x + 0.1 ;
//const float y = (1.0) / (1.0E-10 + x);
f << x << "\t" << y << /*"\t" << y*0.9f << "\t" << y*1.1f <<*/ std::endl ;
}
......
sources/softs/tests/main.cpp
View file @ 70bd1135
......@@ -100,16 +100,12 @@ int parametrization_tests()
std::cout << "<<DEBUG>> cartesian: " << spherical << std::endl;
nb_tests_failed++;
}
std::cout << "<<DEBUG>> rusin 3d: " << rusi << std::endl;
std::cout << "<<DEBUG>> 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 << "<<DEBUG>> spherical before: " << spherical << std::endl;
try
{
params::convert(&spherical[0], params::SPHERICAL_TL_PL_TV_PV, params::RUSIN_TH_TD_PD, &rusi[0]);
......@@ -123,15 +119,42 @@ int parametrization_tests()
std::cout << "<<DEBUG>> spherical: " << spherical << std::endl;
nb_tests_failed++;
}
std::cout << "<<DEBUG>> rusin after: " << rusi << std::endl;
std::cout << "<<DEBUG>> spherical after: " << spherical << std::endl;
std::cout << "<<DEBUG>> " << 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 << "<<DEBUG>> 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;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment