SphericalHarmonics.h

gehe zur Dokumentation dieser Datei

#ifndef _ZFXMATH_INCLUDE_SH_H_
#define _ZFXMATH_INCLUDE_SH_H_

namespace ZFXMath
{

namespace SphericalHarmonics
{

template<class PrecisionType>
struct TSample
{
    // Default constructor
    TSample() { };

    // Construct from spherical co-ordinates
    TSample( const PrecisionType u, const PrecisionType v ) :
        theta(u),
        phi(v)
    {
        x = PrecisionType( sin(theta) * cos(phi) );
        y = PrecisionType( sin(theta) * sin(phi) );
        z = PrecisionType( cos(theta) );
    }

    PrecisionType   theta, phi;

    PrecisionType   x, y, z;
};

template<class PrecisionType, class FuncValueType>
struct TSphericalFunction
{
    virtual FuncValueType   EvalFunction(const TSample<PrecisionType>& s) const = 0;
};

template<class T>
class TCoeffs
{
    public:
        template<class TAlloc>
        class TAllocator
        {
            public:
                TAllocator( int size )
                {
                    m_Size = size;
                    m_pMemory = new TAlloc[size];
                    m_pNextMem = m_pMemory;
                }

                ~TAllocator()
                {
                    delete[] m_pMemory;
                }

                TAlloc* GetMem( int amount )
                {
                    m_Size -= amount;
                    assert( m_Size >= 0 );

                    TAlloc* pMem = m_pNextMem;
                    m_pNextMem += amount;

                    return pMem;
                }

            private:
                TAlloc*             m_pMemory;
                TAlloc*             m_pNextMem;
                int                 m_Size;
        };
        typedef TAllocator<T>   Allocator;

        TCoeffs() : m_Coeffs(0), m_NbBands(0), m_NbCoeffs(0) { }

        TCoeffs(const int nb_bands) :
            m_Coeffs(0),
            m_NbBands(nb_bands),
            m_NbCoeffs(m_NbBands * m_NbBands),
            m_FromAllocator(false)
        {
            // Allocate the coeff list
            m_Coeffs = new T[m_NbCoeffs];
        }

        TCoeffs(const int nb_bands, Allocator* pAlloc) :
            m_Coeffs(0),
            m_NbCoeffs(nb_bands * nb_bands),
            m_FromAllocator(true)
        {
            // Allocate the coeff list
            m_Coeffs = pAlloc->GetMem( m_NbCoeffs );
        }

        ~TCoeffs()
        {
            // Stack allocator cleans itself up
            if (!m_FromAllocator)
                delete [] m_Coeffs;
        }

        TCoeffs(const TCoeffs& other) :
            m_Coeffs(0),
            m_NbBands(other.m_NbBands),
            m_NbCoeffs(other.m_NbCoeffs),
            m_FromAllocator(other.m_FromAllocator)
        {
            // Just copy the pointer in this case
            if (m_FromAllocator)
                m_Coeffs = other.m_Coeffs;

            else
            {
                // Make a copy when not using the stack allocator
                m_Coeffs = new T[m_NbCoeffs];
                for (int i = 0; i < m_NbCoeffs; i++)
                    m_Coeffs[i] = other.m_Coeffs[i];
            }
        }

        // Const/non-const accessors for band/arg
        T operator () (const int l, const int m) const
        {
            return (m_Coeffs[Check(l * (l + 1) + m)]);
        }
        T& operator () (const int l, const int m)
        {
            return (m_Coeffs[Check(l * (l + 1) + m)]);
        }

        // Const/non-const accessors by index
        T operator () (const int i) const
        {
            return (m_Coeffs[Check(i)]);
        }
        T& operator () (const int i)
        {
            return (m_Coeffs[Check(i)]);
        }

        int GetNbBands(void) const
        {
            return (m_NbBands);
        }

        int GetSize(void) const
        {
            return (m_NbCoeffs);
        }

    private:
        inline int Check(const int index) const
        {
            // Check bounds in debug build
            assert(index >= 0 && index < m_NbCoeffs);

            return (index);
        }

        // List of SH coefficients
        T*      m_Coeffs;

        // Number of bands used
        int     m_NbBands;

        // Number of coefficients (=bands^2)
        int     m_NbCoeffs;

        // Allocated from a special-purpose allocator?
        bool    m_FromAllocator;

};


// http://www.research.scea.com/gdc2003/spherical-harmonic-lighting.html
// Project a spherical function onto the spherical harmonic bases
//
// a. Spherical Function to project
// b. List of sample locations on the sphere
// c. List of y(l,m,theta,phi) for each sample
// d. Number of samples on the sphere
// e. Destination projection
template<class PrecisionType, class FuncValueType>
void Project(const TSphericalFunction<PrecisionType,FuncValueType>& func,
                const TSample<PrecisionType>* samples,
                const TCoeffs<PrecisionType>* coeffs,
                const int nb_samples,
                TCoeffs<FuncValueType>& dest )
{
    int nb_coeffs = dest.GetSize();

    // Check input
    assert(samples && coeffs);
    assert(nb_samples >= 0);
    assert(coeffs[0].GetSize() == nb_coeffs);

    // Clear out sums
    for (int i = 0; i < nb_coeffs; i++)
        dest(i) = 0;

    for (int i = 0; i < nb_samples; i++)
    {
        // Take the sample at this point on the sphere
        const TSample<PrecisionType>&   s = samples[i];
        const TCoeffs<PrecisionType>&   c = coeffs[i];
        FuncValueType sample = func.EvalFunction(s);

        // Sum the projection of this sample onto each SH basis
        for (int j = 0; j < nb_coeffs; j++)
            dest(j) += sample * FuncValueType( c(j) );
    }

    // Divide each coefficient by the number of samples and multiply by weights
    for (int i = 0; i < nb_coeffs; i++)
        dest(i) = dest(i) * FuncValueType( PrecisionType( (4 * PI) / nb_samples ) );
}

//---------------------------------------------------------------------------------------
// Evaluate Associated Legendre Polynomial
template<class PrecisionType>
PrecisionType P(const int l, const int m, const PrecisionType x)
{
    // Start with P(0,0) at 1
    PrecisionType pmm = 1;

    // First calculate P(m,m) since that is the only rule that requires results
    // from previous bands
    // No need to check for m>0 since SH function always gives positive m

    // Precalculate (1 - x^2)^0.5
    PrecisionType somx2 = (PrecisionType)sqrt(1 - x * x);

    // This calculates P(m,m). There are three terms in rule 2 that are being iteratively multiplied:
    //
    // 0: -1^m
    // 1: (2m-1)!!
    // 2: (1-x^2)^(m/2)
    //
    // Term 2 has been partly precalculated and the iterative multiplication by itself m times
    // completes the term.
    // The result of 2m-1 is always odd so the double factorial calculation multiplies every odd
    // number below 2m-1 together. So, term 3 is calculated using the 'fact' variable.
    PrecisionType fact = 1;
    for (int i = 1; i <= m; i++)
    {
        pmm *= -1 * fact * somx2;
        fact += 2;
    }

    // No need to go any further, rule 2 is satisfied
    if (l == m)
        return (pmm);

    // Since m<l in all remaining cases, all that is left is to raise the band until the required
    // l is found

    // Rule 3, use result of P(m,m) to calculate P(m,m+1)
    PrecisionType pmmp1 = x * (2 * m + 1) * pmm;

    // Is rule 3 satisfied?
    if (l == m + 1)
        return (pmmp1);

    // Finally, use rule 1 to calculate any remaining cases
    PrecisionType pll = 0;
    for (int ll = m + 2; ll <= l; ll++)
    {
        // Use result of two previous bands
        pll = PrecisionType( (x * (2.0f * ll - 1.0f) * pmmp1 - (ll + m - 1.0f) * pmm) / (ll - m) );

        // Shift the previous two bands up
        pmm = pmmp1;
        pmmp1 = pll;
    }

    return (pll);
}

// Evaluate Real Spherical Harmonic
template<class PrecisionType>
PrecisionType y(const int l, const int m, const PrecisionType theta, const PrecisionType phi)
{
    if (0 == m)
        return (PrecisionType(K<PrecisionType>(l, 0) * P(l, 0, (PrecisionType)cos(theta))));

    if (m > 0)
        return PrecisionType(sqrt(2.0) * K<PrecisionType>(l, m) * cos(m * phi) * P(l, m, (PrecisionType)cos(theta)));

    // m < 0, m is negated in call to K
    return PrecisionType(sqrt(2.0) * K<PrecisionType>(l, -m) * sin(-m * phi) * P(l, -m, (PrecisionType)cos(theta)));
}


// Re-normalisation constant for SH function
template<class PrecisionType>
PrecisionType K(const int l, const int m)
{
    struct FactorialGen
    {
        FactorialGen(const int count) :
        factorials(0),
            nb_factorials(count)
        {
            factorials = new int[nb_factorials];

            // Build a factorial LUT
            factorials[0] = 1;
            for (int i = 1; i < nb_factorials; i++)
                factorials[i] = i * factorials[i - 1];
        }

        ~FactorialGen(void)
        {
            delete [] factorials;
        }

        int operator () (const int i) const
        {
            // Check input for debug builds
            assert(i >= 0 && i < nb_factorials);
            return (factorials[i]);
        }

        int*    factorials;
        int     nb_factorials;
    };

    // Generate lut the first time this function is called
    // Note that this method doesn't interact well with explicit init/shutdown memory managers
    // Additionally, 32-bits is not enough to store any higher bands (7 max) !!!
    static FactorialGen factorial(14);

    // Note that |m| is not used here as the SH function always passes positive m
    return PrecisionType(sqrt(((2 * l + 1) * factorial(l - m)) / (4 * PI * factorial(l + m))));
}


}; // namespace SphericalHarmonics

}; // namespace ZFXMath

#endif //_ZFXMATH_INCLUDE_SH_H_

---