The concept and application of spherical harmonic function: spherical harmonic function

The concept and application of spherical harmonic function

 

Properties of spherical harmonic function

Similar to the properties of Fourier series, spherical harmonic function is also based on orthogonal function, and the orthogonal basis of Fourier series is sin(nx) and cos(nx). The spherical harmonic function is the orthogonal base on the sphere. Its main properties are:

  • Standard orthogonality
  • Rotation invariance
  • The integral of the product of a function is equal to the dot product of its spherical harmonic coefficient vector

Standard Orthogonality:

amongRepresents two sets of spherical harmonic function bases.

 

Rotation invariance means:

 

The integral of function product is equal to the dot product of its spherical harmonic coefficient vector:

 

We list the expression of spherical harmonic function:

We can understand it as the Fourier series orthogonal basis of the processed plate. amongIs:

We can write the solution according to the formulaCode:

double Harmonics::K(const int l, const int m)
{
    //When m is less than 0, multiply - 1 to pass in
    return (sqrt(((2 * l + 1) * Factorial(l - m)) / (4 * PI * Factorial(l + m))));
}

//Solving the first order factorials
int Harmonics::Factorial(int v)
{
    if (v == 0)
        return (1);

    int result = v;
    while (--v > 0)
        result *= v;
    return (result);
}

 

Is the adjoint Legendre polynomial (ALP). First, Legendre function P(x) is the solution of Legendre differential equation

Legendre polynomial expression is:

Two parameters l, m are introduced with Legendre polynomials (ALP) to define based on Legendre polynomials. The expression is as follows:

When we solve ALP, we can use three recursions:

Through these three recursions, we can solve arbitrary order random Legendre polynomials. We write the following code according to the formula:

double Harmonics::P(const int l, const int m, const double x)
{
    // Formula 2 does not need recursive calculation
    if (l == m)
        return (pow(-1.0f, m) * DoubleFactorial(2 * m - 1) * pow(sqrt(1 - x * x), m));

    // Formula 3
    if (l == m + 1)
        return (x * (2 * m + 1) * P(m, m, x));

    // Formula 1
    return ((x * (2 * l - 1) * P(l - 1, m, x) - (l + m - 1) * P(l - 2, m, x)) / (l - m));
}

//Quadratic factorial
int Harmonics::DoubleFactorial(int x)
{
    if (x == 0 || x == -1)
        return (1);

    int result = x;
    while ((x -= 2) > 0)
        result *= x;
    return (result);
}

completeandAfter the calculation of, we can calculate the spherical harmonic function according to the expression. The code is as follows:

double Harmonics::y(const int l, const int m, const double theta, const double phi)
{
    if (m == 0)
        return (K(l, 0) * P(l, 0, cos(theta)));

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

    // When m is less than 0, multiply - 1 in advance and send it to K
    return (sqrt(2.0f) * K(l, -m) * sin(-m * phi) * P(l, -m, cos(theta)));
}

The expansion of the sphere is:

amongIs the spherical harmonic coefficient,Is a spherical harmonic function. This formula is very similar to Fourier series. Again, we solve the coefficientsThe process is very similar:

amongIs the original function. When we solve the coefficients of the corresponding order, we can reconstruct them, our reconstructed function is:

The higher the order l is. The reconstruction results are closer to the original function.

The solution code is as follows:

//Calculation coefficient
void Harmonics::Evaluate(const std::vector<Vertex>& vertices)
{
    int n = (degree_ + 1)*(degree_ + 1);
    coefs = vector<glm::vec3>(n,glm::vec3(0.0));
    
    //Integrate all sampling points on the sphere
    for (const Vertex& v : vertices)
    {
        vector<float> Y = Basis(v.theta,v.phi);
        for (int i = 0; i < n; i++)
        {
            //v.color is the color we sampled from the original function f(s)
            coefs[i] = coefs[i] + Y[i] * v.color;
        }
    }
    for (glm::vec3& coef : coefs)
    {
        coef = 4.0f*PI*coef / (float)vertices.size();
    }
}


//Reconstruction of f(x) with coefficients
glm::vec3 Harmonics::Render(const double theta,const double phi)
{
    int n = (degree_ + 1)*(degree_ + 1);

    vector<float> Y = Basis(theta,phi);
    glm::vec3 color=glm::vec3(0.0f);
    for (int i = 0; i < n; i++)
    {
        color = color + Y[i] * coefs[i];
    }
    return color;
}

vector<float> Harmonics::Basis(const double theta,const double phi)
{
    int n = (degree_ + 1)*(degree_ + 1);
    vector<float> Y(n);

    for(int l=0 ; l<=degree_; l++){
        for(int m = -1 * l; m <= l; m++){
            //Using index: n=l(l+1)+m
            Y[l*(l+1)+m] = y(l, m, theta, phi);
        }
    }
    return Y;
}

As the order increases, the effect is as follows:

When l=3, the difference between the effect and the real environment light is less than 1%. Low order is very suitable for ambient light map, ao map and so on. The results of our experiment are as follows:

The lower ball is l=0, 1, 2, 3. Above is the original function.

Tags: less

Posted on Thu, 18 Jun 2020 05:07:18 -0400 by Loathor__Healer