计算机与数学 —— 使用高斯勒让得求积来获得样条的长度

这篇博客介绍了如何使用高斯勒让得求积（Legendre-Gauss quadrature）来获得样条函数的长度的方法。

---

背景

在我构建海飞丝插件（Hair Strand Plugin）时，由于每一根发束的构建都是使用了Hermit spline，在某些约束运算的时候需要获得样条的长度。因此查找了一些资料，最终发现在数值方法中，高斯勒让得求积的方法的表现比较好。

总体来讲，只需要5次采样就能够获得一个合理的结果。

高斯勒让得求积

当我们想要求解函数∫1−1f(x)dx的积分值时，往往可以使用∑ni=1wif(xi)来进行近似运算，其中wi,i=1...n为权重值。高斯求积仅仅当函数f(x)可以由在区间[-1, 1]的多项式近似时才能获得准确的近似解，这种方法并不适合函数具有奇异点的情况。

于是乎，我们可以把函数f(x)写作f(x)=W(x)g(x)，其中g(x)是近似多项式，W(x)是已知的权重函数，这样我们就有

∫1−1f(x)dx=∫1−1W(x)g(x)dx≈∑i=1nw′ig(xi)

而针对于高斯勒让得求积来说，针对于权重函数W(x)=1时，关联多项式为Pn(x)。

wi=2(1−x2i)[P′n(xi)2]，其中xi是Pn(x)的第i個根。

而Pn(x)函数如下：

Pn(x)=∏1≤i≤n(x−xi)

针对于五个采样点的计算，其wi和Pi的值如下：

Pi wi 0 128255 ±245−1470√√21 322+1370√900 ±245+1470√√21 322−1370√900

这些值可以提前计算出来，到时候直接使用：

    struct FLegendreGaussCoefficient    {        float Abscissa;        float Weight;    };    static const FLegendreGaussCoefficient LegendreGaussCoefficients[] =    {        { 0.0f, 0.5688889f },        { -0.5384693f, 0.47862867f },        { 0.5384693f, 0.47862867f },        { -0.90617985f, 0.23692688f },        { 0.90617985f, 0.23692688f }    };

代码实现

了解了算法之后其实代码实现就很简单了，输入四个向量P0,P1,T0,T1分别代表开始与结束点的位置和切线，返回这段样条的长度。

代码如下：

float LengthOfSegment(FVector P0, FVector P1, FVector T0, FVector T1){    FVector Coeff1 = ((P0 - P1) * 2.0f + T0 + T1) * 3.0f;    FVector Coeff2 = (P1 - P0) * 6.0f - T0 * 4.0f - T1 * 2.0f;    FVector Coeff3 = T0;    float HalfParam = 0.5f;    float Length = 0.0f;    for (const auto& LegendreGaussCoefficient : LegendreGaussCoefficients)    {        // Calculate derivative at each Legendre-Gauss sample, and perform a weighted sum        const float Alpha = HalfParam * (1.0f + LegendreGaussCoefficient.Abscissa);        const FVector Derivative = ((Coeff1 * Alpha + Coeff2) * Alpha + Coeff3);        Length += Derivative.Size() * LegendreGaussCoefficient.Weight;    }    Length *= HalfParam;    return Length;}

<全文完>