diff --git a/vcg/math/spherical_harmonics.h b/vcg/math/spherical_harmonics.h
new file mode 100644
index 00000000..f100f8d3
--- /dev/null
+++ b/vcg/math/spherical_harmonics.h
@@ -0,0 +1,155 @@
+/****************************************************************************
+* 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.                                                         *
+*                                                                           *
+****************************************************************************/
+
+#ifndef __VCGLIB_SPHERICAL_HARMONICS_H
+#define __VCGLIB_SPHERICAL_HARMONICS_H
+
+#include "vcg/math/base.h"
+#include "vcg/math/legendre.h"
+#include "vcg/math/factorial.h"
+
+namespace vcg{
+namespace math{
+
+/**
+ * Although the Real Spherical Harmonic Function is correctly defined over any 
+ * 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
+ * for positive m only.
+ */
+template <typename ScalarType>
+class SphericalHarmonics{
+
+private :
+	inline static ScalarType scaling_factor(unsigned l, unsigned 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)
+	{
+		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)
+	{
+		return scaling_factor(l, m) * Legendre<ScalarType>::AssociatedPolynomial(l, m, Cos(theta), Sin(theta)) * Sin(m * phi);
+	}
+	
+	ScalarType * coefficients;
+	int max_band;			
+	
+public :
+	
+	/**
+	 * Returns the Real Spherical Harmonic Function
+	 * 
+	 * l is any positive integer,
+	 * m is such that -l <= m <= l
+	 * theta is inside [0, PI]
+	 * phi is inside [0, 2*PI]
+	 */
+	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);
+		
+		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 return SQRT_TWO * complex_spherical_harmonic_im(l, -m, theta, phi);
+	}
+	
+	typedef ScalarType (*PolarFunction) (ScalarType theta, ScalarType phi);
+	
+	static SphericalHarmonics Project(PolarFunction fun, unsigned max_band, unsigned n_samples)
+	{
+		const ScalarType weight = 4 * M_PI;
+		
+		unsigned sqrt_n_samples = (unsigned int) Sqrt((int)n_samples);
+		unsigned actual_n_samples = sqrt_n_samples * sqrt_n_samples;
+		unsigned n_coeff = max_band * max_band;
+		
+		ScalarType one_over_n = 1.0/(ScalarType)sqrt_n_samples;
+		
+		RandomGenerator rand;
+		SphericalHarmonics sph;
+
+		sph.coefficients = new ScalarType[n_coeff];
+		sph.max_band = max_band;
+		
+		int i = 0;
+		
+		for (unsigned k = 0; k < n_coeff; k++ ) sph.coefficients[k] = 0;
+		
+		for (unsigned a = 0; a < sqrt_n_samples; ++a )
+		{
+			for (unsigned b = 0; b < sqrt_n_samples; ++b)
+			{
+				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 theta = 2.0 * Acos(Sqrt(1.0 - x));
+				ScalarType phi = 2.0 * M_PI * y;
+				
+				for (int l = 0; l < (int)max_band; ++l)
+				{
+					for (int m = -l; m <= l; ++m)
+					{
+						int index = l * (l+1) + m;
+						sph.coefficients[index] += fun(theta, phi) * Real(l, m, theta, phi);
+					}
+				}
+				i++;
+			}
+		}
+		
+		ScalarType factor = weight / actual_n_samples;
+		for(i = 0; i < (int)n_coeff; ++i)
+		{
+			sph.coefficients[i] *= factor;
+		}
+		
+		return sph;
+	}
+	
+	ScalarType operator()(ScalarType theta, ScalarType phi)
+	{
+		ScalarType f = 0;
+		
+		for (int l = 0; l < max_band; ++l)
+		{
+			for (int m = -l; m <= l; ++m)
+			{
+				int index = l * (l+1) + m;
+				f += (coefficients[index] * Real(l, m, theta, phi));
+			}
+		}
+		
+		return f;
+	}
+};
+
+}} //namespace vcg::math
+
+#endif