Commit 6b21de4d authored by Pierre Kraemer's avatar Pierre Kraemer
Browse files

Merge branch 'develop' into 'develop'

code bivariate polynomials

See merge request !75
parents d7fa9652 b2421372
/* /*
* spericalHarmonics.h * bivariatePolynomials.h
* *
* Created on: Oct 2, 2013 * Created on: July, 2015
* Author: sauvage * Author: sauvage
*/ */
#ifndef __SPHERICALHARMONICS_H__ #ifndef __BIVARIATEPOLYNOMIALS_H__
#define __SPHERICALHARMONICS_H__ #define __BIVARIATEPOLYNOMIALS_H__
#include <cassert> #include <cassert>
#include <cmath> #include <cmath>
#include <iostream> #include <iostream>
#include <Eigen/Core> //#include <Eigen/Core>
#include <Eigen/Dense> //#include <Eigen/Dense>
#include <Eigen/Cholesky> //#include <Eigen/Cholesky>
namespace CGoGN namespace CGoGN
{ {
...@@ -23,84 +23,71 @@ namespace Utils ...@@ -23,84 +23,71 @@ namespace Utils
{ {
template <typename Tscalar, typename Tcoef> template <typename Tscalar, typename Tcoef>
class SphericalHarmonics class BivariatePolynomials
{ {
private : private :
static const int max_resolution = 10 ; // max possible resolution for any object of that class const int degree ; // (bi-)degree of the polynomial (=4)
static int resolution ; // actual resolution level const int nb_coefs ; // number of coefs (=15)
static int nb_coefs ; // number of coefs = (resolution + 1) * (resolution + 1) Tcoef* coefs; // table of coefficients
static unsigned long cpt_instances; // number of instances of the class
static Tscalar K_tab [(max_resolution+1)*(max_resolution+1)]; // table containing constants K
static Tscalar F_tab [(max_resolution+1)*(max_resolution+1)]; // table for computing the functions : P or Y or y
Tcoef* coefs; // table of coefficients
public : public :
// construction, destruction and initialization // construction, destruction and initialization
SphericalHarmonics(); BivariatePolynomials();
SphericalHarmonics(SphericalHarmonics const &); BivariatePolynomials(BivariatePolynomials const &);
~SphericalHarmonics(); ~BivariatePolynomials();
static void set_level(int res_level); int get_degree () { return degree; }
static void set_nb_coefs (int nb_coefs); int get_nb_coefs () { return nb_coefs; }
static int get_resolution () { return resolution; }
static int get_nb_coefs () { return nb_coefs; }
// evaluation // evaluation
static void set_eval_direction (Tscalar theta, Tscalar phi) ; // fix the direction in which the SH has to be evaluated Tcoef evaluate_at (Tscalar u, Tscalar v) const;
static void set_eval_direction (Tscalar x, Tscalar y, Tscalar z) ; // fix the direction in which the SH has to be evaluated Tcoef evaluate_at (const Geom::Vec3<Tscalar>& tu, const Geom::Vec3<Tscalar>& tv, const Geom::Vec3<Tscalar>& n, const Geom::Vec3<Tscalar>& eval_dir) const;
Tcoef evaluate () const; // evaluates at a fixed direction
Tcoef evaluate_at (Tscalar theta, Tscalar phi) const; // eval spherical coordinates
Tcoef evaluate_at (Tscalar x, Tscalar y, Tscalar z) const; // eval cartesian coordinates
// I/O // I/O
const Tcoef& get_coef (int l, int m) const {assert ((l>=0 && l <=resolution) || !" maybe you forgot to call set_level()"); assert (m >= (-l) && m <= l); return get_coef(index(l,m));} const Tcoef& get_coef (int du, int dv) const {assert ((du>=0 && dv <=degree && dv>=0 && dv <=degree) || !" bi-degree is incorrect"); return get_coef(index(l,m));}
Tcoef& get_coef (int l, int m) {assert ((l>=0 && l <=resolution) || !" maybe you forgot to call set_level()"); assert (m >= (-l) && m <= l); return get_coef(index(l,m));} Tcoef& get_coef (int du, int dv) {assert ((du>=0 && dv <=degree && dv>=0 && dv <=degree) || !" bi-degree is incorrect"); return get_coef(index(l,m));}
template <typename TS,typename TC> friend std::ostream & operator<< (std::ostream & os, const SphericalHarmonics<TS,TC> & sh);
// operators Tcoef* get_coef_tab () {return coefs;}
void operator= (const SphericalHarmonics<Tscalar,Tcoef>&);
void operator+= (const SphericalHarmonics<Tscalar,Tcoef>&);
SphericalHarmonics<Tscalar,Tcoef> operator+ (const SphericalHarmonics<Tscalar,Tcoef>&) const;
void operator-= (const SphericalHarmonics<Tscalar,Tcoef> &);
SphericalHarmonics<Tscalar,Tcoef> operator- (const SphericalHarmonics<Tscalar,Tcoef>&) const;
void operator*= (Tscalar);
SphericalHarmonics<Tscalar,Tcoef> operator* (Tscalar) const;
void operator/= (Tscalar);
SphericalHarmonics<Tscalar,Tcoef> operator/ (Tscalar) const;
std::string CGoGNnameOfType() const { return "SphericalHarmonics"; } template <typename TS,typename TC> friend std::ostream & operator<< (std::ostream & os, const BivariatePolynomials<TS,TC> & sh);
// static void copy_K_tab (Tscalar tab[]) ; // obsolete, was used for shaders // operators
// void operator= (const BivariatePolynomials<Tscalar,Tcoef>&);
// void operator+= (const BivariatePolynomials<Tscalar,Tcoef>&);
// BivariatePolynomials<Tscalar,Tcoef> operator+ (const BivariatePolynomials<Tscalar,Tcoef>&) const;
// void operator-= (const BivariatePolynomials<Tscalar,Tcoef> &);
// BivariatePolynomials<Tscalar,Tcoef> operator- (const BivariatePolynomials<Tscalar,Tcoef>&) const;
// void operator*= (Tscalar);
// BivariatePolynomials<Tscalar,Tcoef> operator* (Tscalar) const;
// void operator/= (Tscalar);
// BivariatePolynomials<Tscalar,Tcoef> operator/ (Tscalar) const;
std::string CGoGNnameOfType() const { return "BivariatePolynomials"; }
// fitting // fitting
template <typename Tdirection, typename Tchannel> // template <typename Tdirection, typename Tchannel>
void fit_to_data(int n, Tdirection* t_theta, Tdirection* t_phi, Tchannel* t_R, Tchannel* t_G, Tchannel* t_B, double lambda); // void fit_to_data(int n, Tdirection* t_theta, Tdirection* t_phi, Tchannel* t_R, Tchannel* t_G, Tchannel* t_B, double lambda);
template <typename Tdirection, typename Tchannel> // template <typename Tdirection, typename Tchannel>
void fit_to_data(int n, Tdirection* t_x, Tdirection* t_y, Tdirection* t_z, Tchannel* t_R, Tchannel* t_G, Tchannel* t_B, double lambda); // void fit_to_data(int n, Tdirection* t_x, Tdirection* t_y, Tdirection* t_z, Tchannel* t_R, Tchannel* t_G, Tchannel* t_B, double lambda);
private : private :
static int index (int l, int m) {return l*(l+1)+m;} static int index (int du, int dv) {int n=du+dv; return n*(n+1)/2 + dv;}
// order of coefs : 1,u,v,u2, uv, v2, u3, u2v, uv2, v3, u4, u3v, u2v2, uv3, v4
// evaluation : // Thus i = n*(n+1)/2 + dv where n =
static void init_K_tab (); // compute the normalization constants K_l^m and store them into K_tab
static void compute_P_tab (Tscalar t); // Compute Legendre Polynomials at parameter t and store them in F_tab (only for m>=0)
static void compute_y_tab (Tscalar phi); // Compute the real basis functions y_l^m at (theta, phi) and store them in F_tab (compute_P_tab must have been called before)
const Tcoef& get_coef (int i) const {assert ((i>=0 && i<nb_coefs ) || !" maybe you forgot to call set_level()"); return coefs[i];} const Tcoef& get_coef (int i) const {assert ((i>=0 && i<nb_coefs ) || !" incorrect index"); return coefs[i];}
Tcoef& get_coef (int i) {assert ((i>=0 && i<nb_coefs ) || !" maybe you forgot to call set_level()"); return coefs[i];} Tcoef& get_coef (int i) {assert ((i>=0 && i<nb_coefs ) || !" incorrect index"); return coefs[i];}
// fitting // fitting
template <typename Tchannel> // template <typename Tchannel>
void fit_to_data(int n, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>& mM, Tchannel* t_R, Tchannel* t_G, Tchannel* t_B, double lambda); // void fit_to_data(int n, Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>& mM, Tchannel* t_R, Tchannel* t_G, Tchannel* t_B, double lambda);
}; };
} // namespace Utils } // namespace Utils
} // namespace CGoGN } // namespace CGoGN
#include "Utils/sphericalHarmonics.hpp" #include "Utils/bivariatePolynomials.hpp"
#endif /* __SPHERICALHARMONICS_H__ */ #endif /* __BIVARIATEPOLYNOMIALS_H__ */
/* /*
* sphericalHarmonics.hpp * bivariatePolynomials.hpp
* *
* Created on: Oct 2, 2013 * Created on: July, 2015
* Author: sauvage * Author: sauvage
*/ */
...@@ -11,12 +11,6 @@ namespace CGoGN ...@@ -11,12 +11,6 @@ namespace CGoGN
namespace Utils namespace Utils
{ {
template <typename Tscalar,typename Tcoef> int SphericalHarmonics<Tscalar,Tcoef>::resolution = -1;
template <typename Tscalar,typename Tcoef> int SphericalHarmonics<Tscalar,Tcoef>::nb_coefs = -1;
template <typename Tscalar,typename Tcoef> unsigned long SphericalHarmonics<Tscalar,Tcoef>::cpt_instances = 0;
template <typename Tscalar,typename Tcoef> Tscalar SphericalHarmonics<Tscalar,Tcoef>::K_tab[(max_resolution+1)*(max_resolution+1)];
template <typename Tscalar,typename Tcoef> Tscalar SphericalHarmonics<Tscalar,Tcoef>::F_tab[(max_resolution+1)*(max_resolution+1)];
/************************************************************************* /*************************************************************************
...@@ -24,49 +18,27 @@ construction, destruction and initialization ...@@ -24,49 +18,27 @@ construction, destruction and initialization
**************************************************************************/ **************************************************************************/
template <typename Tscalar,typename Tcoef> template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef>::SphericalHarmonics() BivariatePolynomials<Tscalar,Tcoef>::BivariatePolynomials() : degree(4), nb_coefs(15)
{ {
assert ( (nb_coefs > 0) || !" maybe you forgot to call set_level()"); assert ( (nb_coefs > 0) || !" negative # of coefs");
coefs = new Tcoef[nb_coefs]; coefs = new Tcoef[nb_coefs];
++ cpt_instances;
for (int i = 0; i < nb_coefs; i++) for (int i = 0; i < nb_coefs; i++)
coefs[i] = Tcoef (0); coefs[i] = Tcoef (0);
} }
template <typename Tscalar,typename Tcoef> template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef>::SphericalHarmonics(SphericalHarmonics const & r) BivariatePolynomials<Tscalar,Tcoef>::BivariatePolynomials(BivariatePolynomials const & r) : degree(r.degree), nb_coefs (r.nb_coefs)
{ {
assert ( (nb_coefs > 0) || !" maybe you forgot to call set_level()"); assert ( (nb_coefs > 0) || !" negative # of coefs");
coefs = new Tcoef[nb_coefs]; coefs = new Tcoef[nb_coefs];
++cpt_instances;
for (int i = 0; i < nb_coefs; i++) for (int i = 0; i < nb_coefs; i++)
coefs[i] = r.coefs[i]; coefs[i] = r.coefs[i];
} }
template <typename Tscalar,typename Tcoef> template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef>::~SphericalHarmonics() BivariatePolynomials<Tscalar,Tcoef>::~BivariatePolynomials()
{ {
delete[] coefs; delete[] coefs;
--cpt_instances;
}
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::set_level(int res_level)
{
assert(res_level >= 0 && res_level < max_resolution);
assert(cpt_instances == 0);
resolution = res_level;
nb_coefs = (resolution + 1) * (resolution + 1);
init_K_tab();
}
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::set_nb_coefs(int nbc)
{
assert(nbc > 0);
int sq = ceil(sqrt(nbc)) ;
assert(sq*sq == nbc || !"Number of coefs does not fill the last level") ;
set_level(sq-1) ;
} }
/************************************************************************* /*************************************************************************
...@@ -74,117 +46,46 @@ evaluation ...@@ -74,117 +46,46 @@ evaluation
**************************************************************************/ **************************************************************************/
template <typename Tscalar,typename Tcoef> template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::set_eval_direction (Tscalar theta, Tscalar phi) Tcoef BivariatePolynomials<Tscalar,Tcoef>::evaluate_at (Tscalar u, Tscalar v) const
{ {
compute_P_tab(cos(theta)); Tcoef r (0);
compute_y_tab(phi); Tscalar u2 = u*u;
} Tscalar u3 = u2*u;
Tscalar u4 = u3*u;
template <typename Tscalar,typename Tcoef> Tscalar v2 = v*v;
void SphericalHarmonics<Tscalar,Tcoef>::set_eval_direction (Tscalar x, Tscalar y, Tscalar z) Tscalar v3 = v2*v;
{ Tscalar v4 = v3*v;
compute_P_tab(z);
r += coefs[0];
Tscalar phi (0); r += coefs[1] * u;
if ((x*x + y*y) > 0.0) r += coefs[2] * v;
phi = atan2(y, x); r += coefs[3] * u2;
r += coefs[4] * u * v;
compute_y_tab(phi); r += coefs[5] * v2;
} r += coefs[6] * u3;
r += coefs[7] * u2 * v;
template <typename Tscalar,typename Tcoef> r += coefs[8] * u * v2;
Tcoef SphericalHarmonics<Tscalar,Tcoef>::evaluate () const r += coefs[9] * v3;
{ r += coefs[10] * u4;
Tcoef r (0); // (0.0,0.0,0.0); // TODO : use Tcoef (0) r += coefs[11] * u3 * v;
for (int i = 0; i < nb_coefs; i++) r += coefs[12] * u2 * v2;
{ r += coefs[13] * u * v3;
r += coefs[i] * F_tab[i]; r += coefs[14] * v4;
}
return r; return r;
} }
template <typename Tscalar,typename Tcoef> template <typename Tscalar,typename Tcoef>
Tcoef SphericalHarmonics<Tscalar,Tcoef>::evaluate_at (Tscalar theta, Tscalar phi) const Tcoef BivariatePolynomials<Tscalar,Tcoef>::evaluate_at (const Geom::Vec3<Tscalar>& tu, const Geom::Vec3<Tscalar>& tv, const Geom::Vec3<Tscalar>& n, const Geom::Vec3<Tscalar>& eval_dir) const;
{
set_eval_direction(theta, phi);
return evaluate();
}
template <typename Tscalar,typename Tcoef>
Tcoef SphericalHarmonics<Tscalar,Tcoef>::evaluate_at (Tscalar x, Tscalar y, Tscalar z) const
{
set_eval_direction(x, y, z);
return evaluate();
}
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::init_K_tab ()
{ {
for(int l = 0; l <= resolution; l++) Tscalar u = tu*eval_dir;
Tscalar v = tv*eval_dir;
if (n*eval_dir < 0)
{ {
// recursive computation of the squares Tscalar d = sqrt(u*u + v*v);
K_tab[index(l,0)] = (2*l+1) / (4*M_PI); u/=d;
for (int m = 1; m <= l; m++) v/=d;
K_tab[index(l,m)] = K_tab[index(l,m-1)] / (l-m+1) / (l+m);
// square root + symmetry
K_tab[index(l,0)] = sqrt(K_tab[index(l,0)]);
for (int m = 1; m <= l; m++)
{
K_tab[index(l,m)] = sqrt(K_tab[index(l,m)]);
K_tab[index(l,-m)] = K_tab[index(l,m)];
}
} }
} return evaluate_at(u,v);
/* obsolete : was used for shaders
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::copy_K_tab (Tscalar tab[])
{
assert ( (nb_coefs>0) || !" maybe you forgot to call set_level()");
for (unsigned int i = 0 ; i < nb_coefs ; ++i)
tab[i] = K_tab[i] ;
}
*/
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::compute_P_tab (Tscalar t)
{
// if (t<0) {t=-t*t;} else {t=t*t;} // for plotting only : expand the param near equator
F_tab[index(0,0)] = 1;
for (int l = 1; l <= resolution; l++)
{
F_tab[index(l,l)] = (1-2*l) * sqrt(1-t*t) * F_tab[index(l-1,l-1)]; // first diago
F_tab[index(l,l-1)] = t * (2*l-1) * F_tab[index(l-1,l-1)];// second diago
for (int m = 0; m <= l-2; m++)
{// remaining of the line under the 2 diago
F_tab[index(l,m)] = t * (2*l-1) / (float) (l-m) * F_tab[index(l-1,m)] - (l+m-1) / (float) (l-m) * F_tab[index(l-2,m)];
}
}
}
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::compute_y_tab (Tscalar phi)
{
for (int l = 0; l <= resolution; l++)
{
F_tab[index(l,0)] *= K_tab[index(l,0)]; // remove for plotting
}
for (int m = 1; m <= resolution; m++)
{
Tscalar cos_m_phi = cos ( m * phi );
Tscalar sin_m_phi = sin ( m * phi );
for (int l = m; l <= resolution; l++)
{
F_tab[index(l,m)] *= M_SQRT2; // sqrt(2.0); // remove for plotting
F_tab[index(l,m)] *= K_tab[index(l,m)]; // remove for plotting
F_tab[index(l,-m)] = F_tab[index(l,m)] * sin_m_phi ; // store the values for -m<0 in the upper triangle
F_tab[index(l,m)] *= cos_m_phi;
}
}
} }
/************************************************************************* /*************************************************************************
...@@ -192,12 +93,12 @@ I/O ...@@ -192,12 +93,12 @@ I/O
**************************************************************************/ **************************************************************************/
template <typename Tscalar,typename Tcoef> template <typename Tscalar,typename Tcoef>
std::ostream & operator << (std::ostream & os, const SphericalHarmonics<Tscalar,Tcoef> & sh) std::ostream & operator << (std::ostream & os, const BivariatePolynomials<Tscalar,Tcoef> & p)
{ {
for (int l = 0; l <= sh.resolution; l++) for (int d = 0; d <= sh.degree; d++)
{ {
for (int m = -l; m <= l; m++) for (int du = d; du >=0 ; du--)
os << sh.get_coef(l,m) << "\t"; os << p.get_coef(du,d-du) << "\t";
os << std::endl; os << std::endl;
} }
return os; return os;
...@@ -207,196 +108,11 @@ std::ostream & operator << (std::ostream & os, const SphericalHarmonics<Tscalar, ...@@ -207,196 +108,11 @@ std::ostream & operator << (std::ostream & os, const SphericalHarmonics<Tscalar,
operators operators
**************************************************************************/ **************************************************************************/
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::operator = (const SphericalHarmonics<Tscalar,Tcoef>& sh)
{
for (int i = 0; i < nb_coefs; i++)
this->coefs[i] = sh.coefs[i];
}
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::operator += (const SphericalHarmonics<Tscalar,Tcoef>& sh)
{
for (int i = 0; i < nb_coefs; i++)
this->coefs[i] += sh.coefs[i];
}
template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator + (const SphericalHarmonics<Tscalar,Tcoef>& sh) const
{
SphericalHarmonics<Tscalar,Tcoef> res(*this);
res += sh;
return res;
}
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::operator -= (const SphericalHarmonics<Tscalar,Tcoef>& sh)
{
for (int i = 0; i < nb_coefs; i++)
this->coefs[i] -= sh.coefs[i];
}
template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator - (const SphericalHarmonics<Tscalar,Tcoef>& sh) const
{
SphericalHarmonics<Tscalar,Tcoef> res(*this);
res -= sh;
return res;
}
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::operator *= (Tscalar s)
{
for (int i = 0; i < nb_coefs; i++)
this->coefs[i] *= s;
}
template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator * (Tscalar s) const
{
SphericalHarmonics<Tscalar,Tcoef> res(*this);
res *= s;
return res;
}
template <typename Tscalar,typename Tcoef>
void SphericalHarmonics<Tscalar,Tcoef>::operator /= (Tscalar s)
{
for (int i = 0; i < nb_coefs; i++)
this->coefs[i] /= s;
}
template <typename Tscalar,typename Tcoef>
SphericalHarmonics<Tscalar,Tcoef> SphericalHarmonics<Tscalar,Tcoef>::operator / (Tscalar s) const
{
SphericalHarmonics<Tscalar,Tcoef> res(*this);
res /= s;
return res;
}
/************************************************************************* /*************************************************************************
fitting fitting
**************************************************************************/ **************************************************************************/
template <typename Tscalar,typename Tcoef>
template <typename Tdirection, typename Tchannel>
void SphericalHarmonics<Tscalar,Tcoef>::fit_to_data(
int n,
Tdirection* t_theta, Tdirection* t_phi,
Tchannel* t_R, Tchannel* t_G, Tchannel* t_B,
double lambda)
{
Eigen::MatrixXd mM (nb_coefs,n); // matrix with basis function values, evaluated for all directions
// compute mM
for (int p = 0; p < n; ++p)
{
set_eval_direction(t_theta[p], t_phi[p]);
for (int i = 0; i < nb_coefs; ++i)
{
mM(i,p) = F_tab[i];
}
}
fit_to_data(n, mM, t_R, t_G, t_B, lambda);
}
template <typename Tscalar,typename Tcoef>
template <typename Tdirection, typename Tchannel>
void SphericalHarmonics<Tscalar,Tcoef>::fit_to_data(
int n,
Tdirection* t_x, Tdirection* t_y, Tdirection* t_z,
Tchannel* t_R, Tchannel* t_G, Tchannel* t_B,
double lambda)
{
Eigen::MatrixXd mM (nb_coefs,n); // matrix with basis function values, evaluated for all directions
// compute mM
for (int p=0; p<n; ++p)
{
set_eval_direction(t_x[p],t_y[p],t_z[p]);
for (int i=0; i<nb_coefs; ++i)
{
mM(i,p) = F_tab[i];
}
}
fit_to_data(n, mM, t_R, t_G, t_B, lambda);
}
template <typename Tscalar,typename Tcoef>
template <typename Tchannel>
void SphericalHarmonics<Tscalar,Tcoef>::fit_to_data(
int n,