From 399277570bc9fbe58d64dee077853698b4c2c01f Mon Sep 17 00:00:00 2001 From: cnr-isti-vclab Date: Fri, 16 May 2008 10:36:35 +0000 Subject: [PATCH] First release --- vcg/math/legendre.h | 175 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 175 insertions(+) create mode 100644 vcg/math/legendre.h diff --git a/vcg/math/legendre.h b/vcg/math/legendre.h new file mode 100644 index 00000000..984c3e03 --- /dev/null +++ b/vcg/math/legendre.h @@ -0,0 +1,175 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2006 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ + +/**************************************************************************** + +$Log: not supported by cvs2svn $ + +****************************************************************************/ + +#ifndef __VCGLIB_LEGENDRE_H +#define __VCGLIB_LEGENDRE_H + +#include "vcg/math/base.h" + +namespace vcg { +namespace math { + +/* + * Contrary to their definition, the Associated Legendre Polynomials presented here are + * only computed for positive m values. + * + */ +template +class Legendre { + +private : + + /** + * Legendre Polynomial three term Recurrence Relation + */ + static inline ScalarType legendre_next(unsigned l, ScalarType P_lm1, ScalarType P_lm2, ScalarType x) + { + return ((2 * l + 1) * x * P_lm1 - l * P_lm2) / (l + 1); + } + + /** + * Associated Legendre Polynomial three term Recurrence Relation. + * Raises the band index. + */ + static inline double legendre_next(unsigned l, unsigned m, ScalarType P_lm1, ScalarType P_lm2, ScalarType x) + { + return ((2 * l + 1) * x * P_lm1 - (l + m) * P_lm2) / (l + 1 - m); + } + + /** + * Recurrence relation to compute P_m_(m+1) given P_m_m at the same x + */ + static inline double legendre_P_m_mplusone(unsigned m, ScalarType p_m_m, ScalarType x) + { + return x * (2.0 * m + 1.0) * p_m_m; + } + + /** + * Starting relation to compute P_m_m according to the formula: + * + * pow(-1, m) * double_factorial(2 * m - 1) * pow(1 - x*x, abs(m)/2) + * + * which becomes, if x = cos(theta) : + * + * pow(-1, m) * double_factorial(2 * m - 1) * pow(sin(theta), abs(m)/2) + */ + static inline double legendre_P_m_m(unsigned m, ScalarType sin_theta) + { + ScalarType p_m_m = 1.0; + + if (m > 0) + { + ScalarType fact = 1.0; //Double factorial here + for (unsigned i = 1; i <= m; ++i) + { + p_m_m *= fact * sin_theta; //raising sin_theta to the power of m/2 + fact += 2.0; + } + + if (m&1) //odd m + { + // Condon-Shortley Phase term + p_m_m *= -1; + } + } + + return p_m_m; + } + + static inline double legendre_P_l(unsigned l, ScalarType x) + { + ScalarType p0 = 1; + ScalarType p1 = x; + + if (l == 0) return p0; + + for (unsigned n = 1; n < l; ++n) + { + Swap(p0, p1); + p1 = legendre_next(n, p0, p1, x); + } + + return p1; + } + + /** + * Computes the Associated Legendre Polynomial for any given + * positive m and l, with m <= l and -1 <= x <= 1. + */ + static inline double legendre_P_l_m(unsigned l, unsigned m, ScalarType cos_theta, ScalarType sin_theta) + { + if(m > l) return 0; + if(m == 0) return legendre_P_l(l, cos_theta); //OK + else + { + ScalarType p_m_m = legendre_P_m_m(m, sin_theta); //OK + + if (l == m) return p_m_m; + + ScalarType p_m_mplusone = legendre_P_m_mplusone(m, p_m_m, cos_theta); //OK + + if (l == m + 1) return p_m_mplusone; + + unsigned n = m + 1; + + while(n < l) + { + Swap(p_m_m, p_m_mplusone); + p_m_mplusone = legendre_next(n, m, p_m_m, p_m_mplusone, cos_theta); + ++n; + } + + return p_m_mplusone; + } + } + +public : + + static double Polynomial(unsigned l, ScalarType x) + { + assert (x <= 1 && x >= -1); + return legendre_P_l(l, x); + } + + static double AssociatedPolynomial(unsigned l, unsigned m, ScalarType x) + { + assert (m <= l && x <= 1 && x >= -1); + return legendre_P_l_m(l, m, x, Sqrt(1.0 - x * x) ); + } + + static double AssociatedPolynomial(unsigned l, unsigned m, ScalarType cos_theta, ScalarType sin_theta) + { + assert (m <= l && cos_theta <= 1 && cos_theta >= -1 && sin_theta <= 1 && sin_theta >= -1); + return legendre_P_l_m(l, m, cos_theta, Abs(sin_theta)); + } +}; + +}} //vcg::math namespace + +#endif \ No newline at end of file