Commit 27ebbb6d authored by Laurent Belcour's avatar Laurent Belcour

Moments seem to work. I get strange behaviour for glossy materials

parent 04c2252b
...@@ -51,6 +51,7 @@ static const std::map<params::input, const param_info> input_map = { ...@@ -51,6 +51,7 @@ static const std::map<params::input, const param_info> input_map = {
/* 4D Params */ /* 4D Params */
{params::RUSIN_TH_PH_TD_PD, {"RUSIN_TH_PH_TD_PD", 4, "Complete Half angle parametrization"}}, {params::RUSIN_TH_PH_TD_PD, {"RUSIN_TH_PH_TD_PD", 4, "Complete Half angle parametrization"}},
{params::SPHERICAL_TL_PL_TV_PV, {"SPHERICAL_TL_PL_TV_PV", 4, "Complete classical parametrization"}}, {params::SPHERICAL_TL_PL_TV_PV, {"SPHERICAL_TL_PL_TV_PV", 4, "Complete classical parametrization"}},
{params::STEREOGRAPHIC, {"STEREOGRAPHIC", 4, "Light/View vector in stereographic projection"}},
/* 6D Params */ /* 6D Params */
{params::CARTESIAN, {"CARTESIAN", 6, "Complete vector parametrization"}} {params::CARTESIAN, {"CARTESIAN", 6, "Complete vector parametrization"}}
...@@ -211,6 +212,23 @@ void params::from_cartesian(const double* invec, params::input outtype, ...@@ -211,6 +212,23 @@ void params::from_cartesian(const double* invec, params::input outtype,
#endif #endif
break; break;
case params::STEREOGRAPHIC:
{
// Project the View vector invec[0,1,2] on a 2D direction on the
// surface outvec[0,1]
double dotVN = invec[2];
outvec[0] = invec[0] / (1.0+dotVN);
outvec[1] = invec[1] / (1.0+dotVN);
// Project the Light vector invec[0,1,2] on a 2D direction on the
// surface outvec[2,3]
double dotLN = invec[5];
outvec[2] = invec[3] / (1.0+dotLN);
outvec[3] = invec[4] / (1.0+dotLN);
break;
}
// 6D Parametrization // 6D Parametrization
case params::CARTESIAN: case params::CARTESIAN:
memcpy(outvec, invec, 6*sizeof(double)); memcpy(outvec, invec, 6*sizeof(double));
......
...@@ -8,13 +8,13 @@ ...@@ -8,13 +8,13 @@
#include <cassert> #include <cassert>
#include <iostream> #include <iostream>
/*! \brief a static class allowing to change from one parametrization /*! \class params
* to another.
* \ingroup core * \ingroup core
* \brief a static class allowing to change from one parametrization
* to another.
* *
* \details * Any function object or data object should have an associated
* Any \typedef function object or \typedef data object should have * parametrization.
* an associated parametrization.
* *
* We use the following convention to defined the tangent, normal and * We use the following convention to defined the tangent, normal and
* bi-normal of the surface: * bi-normal of the surface:
...@@ -31,17 +31,22 @@ class params ...@@ -31,17 +31,22 @@ class params
//! <em>unknown</em>. //! <em>unknown</em>.
enum input enum input
{ {
ROMEIRO_TH_TD, RUSIN_TH_PH_TD_PD, /*!< Half-angle parametrization as described in [Rusinkiewicz'98] */
RUSIN_TH_TD, RUSIN_TH_PH_TD,
RUSIN_TH_PH_TD, RUSIN_TH_TD_PD,
RUSIN_TH_TD_PD, RUSIN_TH_TD, /*!< Half-angle parametrization with no azimutal information */
RUSIN_TH_PH_TD_PD, COS_TH_TD,
COS_TH, COS_TH,
COS_TH_TD,
ISOTROPIC_TV_TL_DPHI, STEREOGRAPHIC, /*!< Stereographic projection of the Light and View vectors */
ISOTROPIC_TD_PD, // Difference between two directions such as R and H
CARTESIAN, SPHERICAL_TL_PL_TV_PV, /*!< Light and View vectors represented in spherical coordinates */
SPHERICAL_TL_PL_TV_PV, ISOTROPIC_TV_TL_DPHI, /*!< Light and View vectors represented in spherical coordinates,
with the difference of azimutal coordinates in the last component */
ISOTROPIC_TD_PD, /*!< Difference between two directions such as R and H */
CARTESIAN, /*!< Light and View vectors represented in cartesian coordinates */
UNKNOWN_INPUT UNKNOWN_INPUT
}; };
...@@ -107,14 +112,19 @@ class params ...@@ -107,14 +112,19 @@ class params
static void convert(const double* invec, params::input intype, static void convert(const double* invec, params::input intype,
params::input outtype, double* outvec) params::input outtype, double* outvec)
{ {
double temvec[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // Temp CARTESIAN vectors if(intype != outtype)
to_cartesian(invec, intype, temvec); {
#ifdef DEBUG // temporary CARTESIAN vector
std::cout << "<<DEBUG>> temp vec = [" double temvec[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
<< temvec[0] << ", " << temvec[1] << ", " << temvec[2] << "] => ["
<< temvec[3] << ", " << temvec[4] << ", " << temvec[5] << "]" << std::endl; to_cartesian(invec, intype, temvec);
#endif from_cartesian(temvec, outtype, outvec);
from_cartesian(temvec, outtype, outvec); }
else
{
int dim = dimension(outtype);
for(int i=0; i<dim; ++i) { outvec[i] = invec[i]; }
}
} }
//! \brief convert a input vector in a given parametrization to an //! \brief convert a input vector in a given parametrization to an
...@@ -162,13 +172,13 @@ class params ...@@ -162,13 +172,13 @@ class params
{ {
// Calculate the half vector // Calculate the half vector
double half[3]; double half[3];
half[0] = sin(theta_h)*sin(phi_h); half[0] = sin(theta_h)*cos(phi_h);
half[1] = sin(theta_h)*cos(phi_h); half[1] = sin(theta_h)*sin(phi_h);
half[2] = cos(theta_h); half[2] = cos(theta_h);
// Compute the light vector using the rotation formula. // Compute the light vector using the rotation formula.
out[0] = sin(theta_d)*sin(phi_d); out[0] = sin(theta_d)*cos(phi_d);
out[1] = sin(theta_d)*cos(phi_d); out[1] = sin(theta_d)*sin(phi_d);
out[2] = cos(theta_d); out[2] = cos(theta_d);
// Rotate the diff vector to get the output vector // Rotate the diff vector to get the output vector
......
...@@ -219,6 +219,8 @@ void lookup_brdf_val(double* brdf, double theta_in, double fi_in, ...@@ -219,6 +219,8 @@ void lookup_brdf_val(double* brdf, double theta_in, double fi_in,
theta_diff < 0.0 || theta_diff > 0.5*M_PI || theta_diff < 0.0 || theta_diff > 0.5*M_PI ||
fi_diff > M_PI) fi_diff > M_PI)
{ {
std::cerr << "<<ERROR>> the input vec is incorrect: TL = " << theta_in << ", PL = " << fi_in << ", TV = " << theta_out << ", PV = " << fi_out << std::endl;
std::cerr << "<<ERROR>> the input vec is incorrect: TH = " << theta_half << ", TD = " << theta_diff << ", PD = " << fi_diff << std::endl;
throw; //! \todo Add exception list throw; //! \todo Add exception list
} }
// Find index. // Find index.
......
...@@ -19,29 +19,29 @@ ...@@ -19,29 +19,29 @@
vec coord(vec V, vec L, vec X, vec Y, vec N) vec coord(vec V, vec L, vec X, vec Y, vec N)
{ {
vec pV = V-dot(V,N)*N; vec pV = V-dot(V,N)*N;
vec vCoord(2); vec vCoord(2);
vCoord[0] = dot(pV,X); vCoord[0] = dot(pV,X);
vCoord[1] = dot(pV,Y); vCoord[1] = dot(pV,Y);
vCoord /= (1.0+dot(V,N)); vCoord /= (1.0+dot(V,N));
vec pL = L-dot(L,N)*N; vec pL = L-dot(L,N)*N;
vec lCoord(2); vec lCoord(2);
lCoord[0] = dot(pL,X); lCoord[0] = dot(pL,X);
lCoord[1] = dot(pL,Y); lCoord[1] = dot(pL,Y);
lCoord /= (1.0+dot(L,N)); lCoord /= (1.0+dot(L,N));
if (norm(lCoord)>EPSILON) if (norm(lCoord)>EPSILON)
{ {
vec lDir = normalize(lCoord); vec lDir = normalize(lCoord);
vec temp(2); vec temp(2);
temp[0] = lDir[0]*vCoord[0] + lDir[1]*vCoord[1]; temp[0] = lDir[0]*vCoord[0] + lDir[1]*vCoord[1];
temp[1] = lDir[0]*vCoord[1] - lDir[1]*vCoord[0]; temp[1] = lDir[0]*vCoord[1] - lDir[1]*vCoord[0];
vCoord = temp; vCoord = temp;
} }
return vCoord; return vCoord;
} }
int main(int argc, char** argv) int main(int argc, char** argv)
...@@ -85,69 +85,71 @@ int main(int argc, char** argv) ...@@ -85,69 +85,71 @@ int main(int argc, char** argv)
const int nX = d->dimX(); const int nX = d->dimX();
const int nY = d->dimY(); const int nY = d->dimY();
// Normalisation factor for the integration // Normalisation factor for the integration
const int nb_theta_int = 90; const int nb_theta_int = 90;
const int nb_phi_int = 180; const int nb_phi_int = 360;
const double normalization = M_PI*M_PI / (double)(nb_phi_int*nb_theta_int); const double normalization = M_PI*M_PI / (double)(nb_phi_int*nb_theta_int);
double in_angle[4] = {0.0, 0.0, 0.0, 0.0} ; double in_angle[4] = {0.0, 0.0, 0.0, 0.0} ;
// Static data vec m_0(nY);
vec X(3); X[0] = 1; X[1] = 0; X[2] = 0; vec m_x(nY);
vec Y(3); Y[0] = 0; Y[1] = 1; Y[2] = 0; vec m_y(nY);
vec N(3); N[0] = 0; N[1] = 0; N[2] = 1; vec m_xx(nY);
vec m_xy(nY);
vec m_yy(nY);
// compute cumulants
vec k_x(nY);
vec k_y(nY);
vec k_xx(nY);
vec k_xy(nY);
vec k_yy(nY);
for(int theta_in=0; theta_in<90; theta_in++) for(int theta_in=0; theta_in<90; theta_in++)
{ {
in_angle[0] = theta_in * 0.5*M_PI / 90.0; in_angle[0] = theta_in * 0.5*M_PI / 90.0;
vec m_0(nY); m_0 = vec::Zero(nY); m_0 = vec::Zero(nY);
vec m_x(nY); m_0 = vec::Zero(nY); m_x = vec::Zero(nY);
vec m_y(nY); m_0 = vec::Zero(nY); m_y = vec::Zero(nY);
vec m_xx(nY); m_0 = vec::Zero(nY); m_xx = vec::Zero(nY);
vec m_xy(nY); m_0 = vec::Zero(nY); m_xy = vec::Zero(nY);
vec m_yy(nY); m_0 = vec::Zero(nY); m_yy = vec::Zero(nY);
// Integrate over the light hemisphere // Theta angle in [0 .. PI / 2]
for(int theta_out=0; theta_out<nb_theta_int; theta_out++) for(int theta_out=0; theta_out<nb_theta_int; theta_out++)
{ {
in_angle[2] = theta_out * 0.5*M_PI / (double)nb_theta_int; in_angle[2] = theta_out * 0.5*M_PI / (double)nb_theta_int;
// Azimutal angle in [-PI .. PI]
for(int phi_out=0; phi_out<nb_phi_int; phi_out++) for(int phi_out=0; phi_out<nb_phi_int; phi_out++)
{ {
in_angle[3] = phi_out * 2.0*M_PI / nb_phi_int; in_angle[3] = phi_out * 2.0*M_PI / (double)nb_phi_int - M_PI;
vec in(nX), cart(6), L(3), V(3); vec in(nX), stereographics(4);
params::convert(in_angle, params::SPHERICAL_TL_PL_TV_PV, data_param, &in[0]); 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]); params::convert(in_angle, params::SPHERICAL_TL_PL_TV_PV, params::STEREOGRAPHIC, &stereographics[0]);
L[0] = cart[0];
L[1] = cart[1];
L[2] = cart[2];
V[0] = cart[3];
V[1] = cart[4];
V[2] = cart[5];
// Copy the input vector // Copy the input vector
vec x = d->value(in); vec x = d->value(in);
// Get the projected 2D coordinate
vec xy = coord(V, L, X, Y, N);
for(int i=0; i<nY; ++i) for(int i=0; i<nY; ++i)
{ {
double val = x[i] /* cos(in_angle[2])*/ * normalization; double val = x[i] * /* cos(in_angle[2]) */ normalization;
m_0[i] += val; m_0[i] += val ;
m_x[i] += val * xy[0]; m_x[i] += val * stereographics[2];
m_y[i] += val * xy[1]; m_y[i] += val * stereographics[3];
m_xx[i] += val * xy[0] * xy[0]; m_xx[i] += val * stereographics[2] * stereographics[2];
m_xy[i] += val * xy[0] * xy[1]; m_xy[i] += val * stereographics[2] * stereographics[3];
m_yy[i] += val * xy[1] * xy[1]; m_yy[i] += val * stereographics[3] * stereographics[3];
} }
} }
} }
/*
for(int i=0; i<nY; ++i) for(int i=0; i<nY; ++i)
{ {
m_x[i] /= m_0[i]; m_x[i] /= m_0[i];
...@@ -156,13 +158,13 @@ int main(int argc, char** argv) ...@@ -156,13 +158,13 @@ int main(int argc, char** argv)
m_xy[i] /= m_0[i]; m_xy[i] /= m_0[i];
m_yy[i] /= m_0[i]; m_yy[i] /= m_0[i];
} }
*/
// compute cumulants // compute cumulants
vec k_x = m_x; k_x = m_x;
vec k_y = m_y; k_y = m_y;
vec k_xx = m_xx;// - product(m_x,m_x); k_xx = m_xx - product(m_x,m_x);
vec k_xy = m_xy;// - product(m_x,m_y); k_xy = m_xy - product(m_x,m_y);
vec k_yy = m_yy;// - product(m_y,m_y); k_yy = m_yy - product(m_y,m_y);
// Output the value into the file // Output the value into the file
file << in_angle[0] << "\t"; file << in_angle[0] << "\t";
......
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