manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Spherical Harmonics are templatized on the number of coefficients

This commit is contained in:
Paolo Cignoni 2008-10-08 14:01:34 +00:00
parent 2da37bd5f7
commit 39e2cf2b3e
1 changed files with 27 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

57
vcg/math/spherical_harmonics.h
View File

@ -35,12 +35,12 @@ namespace vcg{
namespace math{ namespace math{
/** /**
* Although the Real Spherical Harmonic Function is correctly defined over any * Although the Real Spherical Harmonic Function is correctly defined over any
* positive l and any -l <= m <= l, the two internal functions computing the * positive l and any -l <= m <= l, the two internal functions computing the
* imaginary and real parts of the Complex Spherical Harmonic Functions are defined * imaginary and real parts of the Complex Spherical Harmonic Functions are defined
* for positive m only. * for positive m only.
*/ */
template <typename ScalarType> template <typename ScalarType, int MAX_BAND>
class SphericalHarmonics{ class SphericalHarmonics{
private : private :
@ -48,25 +48,25 @@ private :
{ {
return Sqrt( ( (2.0*l + 1.0) * Factorial<ScalarType>(l-m) ) / (4.0 * M_PI * Factorial<ScalarType>(l + m)) );; return Sqrt( ( (2.0*l + 1.0) * Factorial<ScalarType>(l-m) ) / (4.0 * M_PI * Factorial<ScalarType>(l + m)) );;
} }
inline static ScalarType complex_spherical_harmonic_re(unsigned l, unsigned m, ScalarType theta, ScalarType phi) inline static ScalarType complex_spherical_harmonic_re(unsigned l, unsigned m, ScalarType theta, ScalarType phi)
{ {
return scaling_factor(l, m) * Legendre<ScalarType>::AssociatedPolynomial(l, m, Cos(theta), Sin(theta)) * Cos(m * phi); return scaling_factor(l, m) * Legendre<ScalarType>::AssociatedPolynomial(l, m, Cos(theta), Sin(theta)) * Cos(m * phi);
} }
inline static ScalarType complex_spherical_harmonic_im(unsigned l, unsigned m, ScalarType theta, ScalarType phi) inline static ScalarType complex_spherical_harmonic_im(unsigned l, unsigned m, ScalarType theta, ScalarType phi)
{ {
return scaling_factor(l, m) * Legendre<ScalarType>::AssociatedPolynomial(l, m, Cos(theta), Sin(theta)) * Sin(m * phi); return scaling_factor(l, m) * Legendre<ScalarType>::AssociatedPolynomial(l, m, Cos(theta), Sin(theta)) * Sin(m * phi);
} }
ScalarType * coefficients; ScalarType coefficients[MAX_BAND * MAX_BAND];
int max_band; static const int max_band = MAX_BAND;
public : public :
/** /**
* Returns the Real Spherical Harmonic Function * Returns the Real Spherical Harmonic Function
* *
* l is any positive integer, * l is any positive integer,
* m is such that -l <= m <= l * m is such that -l <= m <= l
* theta is inside [0, PI] * theta is inside [0, PI]
@ -75,46 +75,43 @@ public :
static ScalarType Real(unsigned l, int m, ScalarType theta, ScalarType phi) static ScalarType Real(unsigned l, int m, ScalarType theta, ScalarType phi)
{ {
assert((int)-l <= m && m <= (int)l && theta >= 0 && theta <= M_PI && phi >= 0 && phi <= 2 * M_PI); assert((int)-l <= m && m <= (int)l && theta >= 0 && theta <= M_PI && phi >= 0 && phi <= 2 * M_PI);
if (m > 0) return SQRT_TWO * complex_spherical_harmonic_re(l, m, theta, phi); if (m > 0) return SQRT_TWO * complex_spherical_harmonic_re(l, m, theta, phi);
else if (m == 0) return scaling_factor(l, 0) * Legendre<ScalarType>::Polynomial(l, Cos(theta)); else if (m == 0) return scaling_factor(l, 0) * Legendre<ScalarType>::Polynomial(l, Cos(theta));
else return SQRT_TWO * complex_spherical_harmonic_im(l, -m, theta, phi); else return SQRT_TWO * complex_spherical_harmonic_im(l, -m, theta, phi);
} }
typedef ScalarType (*PolarFunction) (ScalarType theta, ScalarType phi); typedef ScalarType (*PolarFunction) (ScalarType theta, ScalarType phi);
static SphericalHarmonics Project(PolarFunction fun, unsigned max_band, unsigned n_samples) static SphericalHarmonics Project(PolarFunction fun, unsigned n_samples)
{ {
const ScalarType weight = 4 * M_PI; const ScalarType weight = 4 * M_PI;
unsigned sqrt_n_samples = (unsigned int) Sqrt((int)n_samples); unsigned sqrt_n_samples = (unsigned int) Sqrt((int)n_samples);
unsigned actual_n_samples = sqrt_n_samples * sqrt_n_samples; unsigned actual_n_samples = sqrt_n_samples * sqrt_n_samples;
unsigned n_coeff = max_band * max_band; unsigned n_coeff = MAX_BAND * MAX_BAND;
ScalarType one_over_n = 1.0/(ScalarType)sqrt_n_samples; ScalarType one_over_n = 1.0/(ScalarType)sqrt_n_samples;
RandomGenerator rand; RandomGenerator rand;
SphericalHarmonics sph; SphericalHarmonics sph;
sph.coefficients = new ScalarType[n_coeff];
sph.max_band = max_band;
int i = 0; int i = 0;
for (unsigned k = 0; k < n_coeff; k++ ) sph.coefficients[k] = 0; for (unsigned k = 0; k < n_coeff; k++ ) sph.coefficients[k] = 0;
for (unsigned a = 0; a < sqrt_n_samples; ++a ) for (unsigned a = 0; a < sqrt_n_samples; ++a )
{ {
for (unsigned b = 0; b < sqrt_n_samples; ++b) for (unsigned b = 0; b < sqrt_n_samples; ++b)
{ {
ScalarType x = (a + rand(INT_MAX)/(ScalarType)INT_MAX) * one_over_n; ScalarType x = (a + rand(INT_MAX)/(ScalarType)INT_MAX) * one_over_n;
ScalarType y = (b + rand(INT_MAX)/(ScalarType)INT_MAX) * one_over_n; ScalarType y = (b + rand(INT_MAX)/(ScalarType)INT_MAX) * one_over_n;
ScalarType theta = 2.0 * Acos(Sqrt(1.0 - x)); ScalarType theta = 2.0 * Acos(Sqrt(1.0 - x));
ScalarType phi = 2.0 * M_PI * y; ScalarType phi = 2.0 * M_PI * y;
for (int l = 0; l < (int)max_band; ++l) for (int l = 0; l < (int)max_band; ++l)
{ {
for (int m = -l; m <= l; ++m) for (int m = -l; m <= l; ++m)
@ -126,20 +123,20 @@ public :
i++; i++;
} }
} }
ScalarType factor = weight / actual_n_samples; ScalarType factor = weight / actual_n_samples;
for(i = 0; i < (int)n_coeff; ++i) for(i = 0; i < (int)n_coeff; ++i)
{ {
sph.coefficients[i] *= factor; sph.coefficients[i] *= factor;
} }
return sph; return sph;
} }
ScalarType operator()(ScalarType theta, ScalarType phi) ScalarType operator()(ScalarType theta, ScalarType phi)
{ {
ScalarType f = 0; ScalarType f = 0;
for (int l = 0; l < max_band; ++l) for (int l = 0; l < max_band; ++l)
{ {
for (int m = -l; m <= l; ++m) for (int m = -l; m <= l; ++m)
@ -148,7 +145,7 @@ public :
f += (coefficients[index] * Real(l, m, theta, phi)); f += (coefficients[index] * Real(l, m, theta, phi));
} }
} }
return f; return f;
} }
}; };