球谐函数的概念与应用：球谐函数

球谐函数的概念与应用：球谐函数简介

 

球谐函数性质

和傅里叶级数性质类似，球谐函数也是以正交函数作为基底，傅里叶级数的正交基底为sin(nx)和cos(nx)。球谐函数则是球面上的正交基底。其主要性质有：

• 标准正交性
• 旋转不变性
• 函数乘积的积分等于其球谐系数向量的点积

标准正交性表示：

其中，表示两组球谐函数基底。

 

旋转不变性表示：

 

函数乘积的积分等于其球谐系数向量的点积：

 

我们列出球谐函数的表达式：

我们可以把它理解为加工版的傅里叶级数正交基底。其中为：

我们可以根据公式写出求解的代码：

double Harmonics::K(const int l, const int m)
{
    //m小于0的情况在外部乘上-1传入
    return (sqrt(((2 * l + 1) * Factorial(l - m)) / (4 * PI * Factorial(l + m))));
}

//求解一次阶乘
int Harmonics::Factorial(int v)
{
    if (v == 0)
        return (1);

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

 

为伴随勒让德多项式(ALP)。首先勒让德函数P(x)是勒让德微分方程的解：

勒让德多项式表达式为：

随勒让德多项式(ALP)引入了两个参数l，m来基于勒让德多项式定义，其表达式为：

我们求解ALP时，可以使用三条递归式：

通过这三条递归式，可以求解任意阶的随勒让德多项式。我们根据公式写出如下代码：

double Harmonics::P(const int l, const int m, const double x)
{
    // 公式二不需要递归计算
    if (l == m)
        return (pow(-1.0f, m) * DoubleFactorial(2 * m - 1) * pow(sqrt(1 - x * x), m));

    // 公式三
    if (l == m + 1)
        return (x * (2 * m + 1) * P(m, m, x));

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

//二次阶乘
int Harmonics::DoubleFactorial(int x)
{
    if (x == 0 || x == -1)
        return (1);

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

完成和的计算后，我们可以根据表达式计算球谐函数了，下面给出代码：

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)));

    // 当m小于0时，预先乘上-1，再传给K
    return (sqrt(2.0f) * K(l, -m) * sin(-m * phi) * P(l, -m, cos(theta)));
}

在球面的展开式为：

其中是球谐系数，为球谐函数。这式子和傅里叶级数是非常相似的。同样，我们求解系数的过程也十分相似：

其中为原函数。当我们求解完成相应阶数的系数后，可以重构,我们重构后的函数为：

其中阶数l越高。重构的结果越逼近原函数。

我们求解代码如下：

//计算系数
void Harmonics::Evaluate(const std::vector<Vertex>& vertices)
{
    int n = (degree_ + 1)*(degree_ + 1);
    coefs = vector<glm::vec3>(n,glm::vec3(0.0));
    
    //对球面上的所有采样点进行积分
    for (const Vertex& v : vertices)
    {
        vector<float> Y = Basis(v.theta,v.phi);
        for (int i = 0; i < n; i++)
        {
            //v.color是我们从原函数f(s)中采样到的颜色
            coefs[i] = coefs[i] + Y[i] * v.color;
        }
    }
    for (glm::vec3& coef : coefs)
    {
        coef = 4.0f*PI*coef / (float)vertices.size();
    }
}


//用系数重构f(x)
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++){
            //利用索引:n=l(l+1)+m
            Y[l*(l+1)+m] = y(l, m, theta, phi);
        }
    }
    return Y;
}

随着阶数提高，效果如下：

其中l=3的情况下，其效果和环境光真实差距已经不到1%。低阶非常适合应用于环境光贴图，以及ao贴图等。我们实验尝试的效果为：

其中下方球分别为l=0，1，2，3，4。上方则为原始函数。