分段三次Hermite插值

上一节中我们提到，在工程实际中，出于数值稳定性以及龙格现象的考虑， 人们常用一些分段的低次多项式来完成多个节点(采样点)的插值。为了使连接处尽量平滑，人们找到了 Hermite 插值多项式来保证 \(C^1\) 的连续性。 对于前后两个采样点 \(x_{i}, x_{i+1}\) 根据它们的函数值 \(f(x_i), f(x_{i+1})\) 和一阶导数值 \(f'(x_i), f'(x_{i+1})\)，我们正好可以得到一个三次的多项式。 那么对于 \(x_0, \cdots, x_n\) 个采样点，我们可以将整个 \([x_0, x_n]\) 分成 \(n\) 段三次多项式。这就是所谓的分段三次 Hermite 插值 (PCHIP, Piecewise Cubic Hermite Interpolation)。

1. 朴素分段三次 Hermite 插值

我们用 \(x_0, \cdots, x_n\) 共 \(n+1\) 个采样点，将闭区间 \([x_0, x_n]\) 分成 \(n\) 段。对于其中某个特定的子区间 \([x_i, x_{i+1}]\)， 我们称 \(x_i\) 为它的左端点，函数值为 \(f(x_i)\)，导数为 \(f'(x_i)\)。\(x_{i+1}\) 为右端点，函数值为 \(f(x_{i+1})\)，导数为 \(f'(x_{i+1})\)。 直接套用上一节中的式(1) 和式(2)就可以写出一个三次多项式。

为了保证算法的数值稳定性，避免过大的数相乘导致浮点数精度丢失，工业界引入了一个局部的归一化变量 \(t\)，将区间 \([x_i, x_{i+1}]\) 映射到 \([0, 1]\) 上。

在 \([0, 1]\) 的标准区间上，套用 Hermite 插值多项式公式，可以用四个基函数的线性组合，写出区间 \([x_i, x_{i+1}]\) 上的局部三次多项式 \(H_3(t)\):

这四个基函数 \(h_{00}(t), h_{10}(t), h_{01}(t), h_{11}(t)\) 分别控制着左右两个端点的函数值和斜率，在计算机图形学中也被称为 Hermite 样条基。

我们定义了 类 PieceCubicHermite， 分别用三个数组 mXs, mYs, mDYs 记录各个采样点及其函数值、导数值。并在其成员函数 Evaluate 中， 先扫描 mXs 确认输入点 \(x\) 所在的子区间，再套用上式 (\(\ref{f1},\ref{f2}\)) 完成求解。 如果点 \(x\) 恰好与某个采样点 \(x_i\) 重合，我们就直接返回 mYs[i] 作为插值。

        Scalar Evaluate(Scalar x) const {
            if (x <= mXs.front()) return mYs.front();
            if (x >= mXs.back())  return mYs.back();
            // 定位 x 所述区间 [x_i, x_{i+1}]
            int i = 1;
            for (; i < mXs.size(); ++i) {
                if (x == mXs[i]) return mYs[i];
                if (x < mXs[i])  break;
            }
            i--;
            Scalar x0 = mXs[i];     Scalar y0 = mYs[i];     Scalar d0 = mDYs[i];
            Scalar x1 = mXs[i + 1]; Scalar y1 = mYs[i + 1]; Scalar d1 = mDYs[i + 1];

            // 归一化映射：将局部区间映射到标准区间 t \in [0, 1]
            Scalar h = x1 - x0;
            Scalar t = (x - x0) / h;
            Scalar t2 = t * t;
            Scalar t3 = t2 * t;

            Scalar h00 =  2.0 * t3 - 3.0 * t2 + 1.0; // 控制左端点值
            Scalar h10 =        t3 - 2.0 * t2 + t;   // 控制左端点斜率
            Scalar h01 = -2.0 * t3 + 3.0 * t2;       // 控制右端点值
            Scalar h11 =        t3 -       t2;       // 控制右端点斜率
            return h00 * y0 + h10 * h * d0 + h01 * y1 + h11 * h * d1;
        }

2. Fritsch & Carlson导数估计

在实际工程中，我们往往只能采样一堆离散的函数值 \(f(x_i)\)，因为成本、采样原理等诸多因素，可能得不到采样点的一阶导数。 那么我们还想用分段三次Hermite插值，就要想办法自动估计一个合理的 \(f'(x_{i})\)。据说 Matlab 的 pchip 方法采用的就是一种保形(shape-preserving)的估计方法。

假设我们有离散数据点 \(x_i, f(x_i)\)。PCHIP 算法首先计算子区间的长度 \(h_i\) 和割线斜率 \(s_{i}\): $$ h_{i}=x_{i+1}-x_{i},\quad s_{i}=\frac{f(x_{i+1})-f(x_i)}{h_i} $$

对于区间内部的一个采样点 \(x_i, i \neq 0,n\)，它左边的割线斜率是 \(s_{i-1}\)，右边的割线斜率是 \(s_{i}\)。PCHIP 算法会根据这两个斜率的关系，分三种情况来估计该点的导数 \(f'(x_{i})\):

• \(s_{i-1} \cdot s_i ≤ 0\)，即两边的割线斜率异号，或者其中一个为 0。这意味着当前的采样点 \(x_i\) 处于一个局部的波峰、波谷，或处于一段平台上。 此时，PCHIP 就会要求 \(f'(x_i) = 0\)。
• \(s_{i-1}\) 与 \(s_i\) 同号并且相邻区间长度相等(\(h_{i-1} = h_i\))。此时 PCHIP 会采用调和平均数来估计导数 $$ f'(x_i) = \frac{2}{\frac{1}{s_{i-1}} + \frac{1}{s_i}} = \frac{2s_{i-1}s_i}{s_{i-1} + s_i} $$
• 一般情况下，相邻区间的长度不相等，此时采用加权的调和平均数 $$ f'(x_i) = \frac{w_1 + w_2}{\frac{w_1}{s_{i-1}} + \frac{w_2}{s_i}}, \qquad w_1 = 2h_i + h_{i-1}, \qquad w_2 = h_i + 2h_{i-1} $$

对于区间两端的采样点 \(x_0, x_n\) 由于缺少一边的割线斜率，需要特别估计一个导数值。以左端点 \(x_0\) 为例，假设它右侧的两个区间割线斜率为 \(s_0, s_1\)， 那么 PCHIP 先通过单侧二次抛物线写出初始估计值:

如果 \(s_{-1}\) 与 \(s_0\) 异号，说明抛物线在 \(x_0\) 处掉头，会导致插值曲线出现鼓包现象。此时 PCHIP 会令 \(f'(x_0) = 0\)。 如果同号并且 \(|s_{-1}| > 3|s_0|\) 则做截断，令 \(f'(x_0) = 3 s_0\)。其它情况则令 \(f'(x_0) = s_{-1}\)。对称的区间的右端点 \(x_n\) 也可以通过类似的讨论写出 \(f'(x_n)\)。 这种导数估计方法是 Fritsch 和 Carlson 在 1980 年提出来的，所以我们将上述估计方法的实现命名为 FritschCarlsonDerivatives。 这种方法的优点是能够尽可能的抑制函数值阶跃变换时产生的超调震荡。

3. Catmull & Rom 导数估计

在控制领域，人们不希望看到测量曲线出现超调震荡，这很可能导致控制信号也跟着出现超调。在计算机图形里面，人们会追求曲线的圆滑度，也会使用抛物线拟合的方法， 经典的 Catmull-Rom 样条就是一种常用的方法。

假设我们要估计内部节点 \(x_{i}\) 处的未知导数 \(f'(x_{i})\)，抛物线拟合就是要用 \(x_{i}\) 及 \(x_{i-1}, x_{i+1}\)，在区间的局部构建的二次抛物线 \(P(x) = ax^2 + bx + c\)。 那么我们就用 \(P'(x_i) = 2a x_i + b\) 作为导数 \(f'(x_i)\) 的估计值。实际上人们并没有尝试解出系数 \(a,b,c\)，而是将这三个点代入拉格朗日多项式求导，整理之后得到一个加权算术平均数:

对于端点 \(x_0, x_n\)，我们按照 PCHIP 的策略，采用单侧二次抛物线外推，但不考虑外推割线斜率过大的情况，我们实现了 CatmullRomDerivatives。 右侧是两种导数估计下的对于一个阶跃信号的插值，两种方法的特点还是很明显的。

4. 完