【原】分段样条曲线

1.8分段样条曲线(标准形式)

复杂的曲线不能适当地用三次样条曲线来建模。它们通常是S型曲线，而复杂的曲线可能包含多个扭曲和转弯。一个选择是使用高阶样条曲线来建模曲线；但是，它们需要解高阶方程，这增加了系统的计算开销和时间延迟。而且，高阶样条对CPs的微小变化过于敏感，这通常是不可取的，因为我们通常希望通过调整它们的CPs来进行样条的微小改变，而不是在形状上进行大的改变。这种曲线最好使用多个端对端连接的三次样条曲线来建模。这些被称为分段样条曲线。

考虑四个给定的点 和 ，需要找到通过它们的分段样条曲线的方程。从本质上讲，这意味着不是一个通过四个点的三次样条曲线，而是需要找到通过每对点的三个单独的样条曲线，如图所示。

假设给定点的坐标是 和 。让三个三次曲线段分别在点 和 、 和 、以及 和 之间被称为 和 。与以前一样，让起始的三次方程式为 的形式。由于现在有三个曲线段，因此需要有三个不同的系数集，如下所示：

总共有12个不同的未知数，至少需要12个不同的方程来解决它们。

为了制定这12个方程，使用各种约束条件来确保三个单独的样条曲线段连接在一起形成一个单一的平滑曲线。第一个约束条件被称为连续条件，它规定了为了形成平滑曲线，三个样条在它们的连接点上应该物理上相遇（Hearn和Baker，1996年）。换句话说，样条应该经过点和，样条应该经过点和，样条应该经过点和。将点坐标代入各自的起始方程，可以得到以下六个方程。如果表示通过点的段，我们可以写成：

第二个要遵守的约束条件被称为连续条件，它规定为了形成平滑曲线，各个样条曲线段的斜率在它们的交点处应该相等（Hearn和Baker，1996年）。对样条方程求导，得到以下结果：

在这种情况下：在处的斜率 在处的斜率。如果表示点处的段的斜率，我们有：

重新排列：

此外，在处的斜率 在处的斜率。重新排列方程得到：

第三个要遵守的约束条件是连续条件，它规定为了形成平滑曲线，各个样条曲线段在它们的交点处的曲率应该相等（Hearn和Baker，1996年）。对样条方程进行二次导数运算，得到以下结果：

在这种情况下，在处的曲率 在处的曲率。如果表示点处的段的曲率，我们有：

重新排列：

此外，在处的曲率 在处的曲率。重新排列方程得到：

重新排列：

最后要考虑的约束条件涉及端点条件。样条的起始斜率和样条的结束斜率也应该是已知的，以便明确指定样条。设起始和结束点的斜率分别为和。因此，在这种情况下，，。将斜率值代入导数方程，得到以下结果：

为了找到这个系统的解，将所有12个方程都代入矩阵形式中：

解这个方程的解为：

示例：

找到通过点, 和的分段三次曲线方程。第一个和最后一个点的斜率分别为 4 和 -2。将给定的值代入解矩阵中：

解出系数并代入起始方程，得到所需的解。

验证: .

clear all; format compact; clc;
syms x;
x1 = 0; y1 = 1;
x2 = 2; y2 = 2;
x3 = 5; y3 = 0;
x4 = 8; y4 = 0;
s1 = 4; s2 = -2;
X1 = [x1, x2, x3, x4];
Y1 = [y1, y2, y3, y4];
Y = [y1; y2; y2; y3; y3; y4; 0; 0; 0; 0; s1; s2];
C = [1, x1, x1^2, x1^3, 0, 0, 0, 0, 0, 0, 0, 0;
 1, x2, x2^2, x2^3, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 1, x2, x2^2, x2^3, 0, 0, 0, 0;
 0, 0, 0, 0, 1, x3, x3^2, x3^3, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 1, x3, x3^2, x3^3;
 0, 0, 0, 0, 0, 0, 0, 0, 1, x4, x4^2, x4^3;
 0, -1, -2*x2, -3*x2^2, 0, 1, 2*x2, 3*x2^2, 0, 0, 0, 0;
 0, 0, 0, 0, 0, -1, -2*x3, -3*x3^2, 0, 1, 2*x3, 3*x3^2;
 0, 0, -2, -6*x2, 0, 0, 2, 6*x2, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, -2, -6*x3, 0, 0, 2, 6*x3;
 0, 1, 2*x1, 3*x1^2, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2*x4, 3*x4^2];
A = inv(C)*Y ;
a1=A(1); b1=A(2); c1=A(3); d1=A(4);
a2=A(5); b2=A(6); c2=A(7); d2=A(8);
a3=A(9); b3=A(10); c3=A(11); d3=A(12);
fprintf('Equations for segments : \n');
yA = a1 + b1*x + c1*x^2 + d1*x^3; yA = vpa(yA, 3)
yB = a2 + b2*x + c2*x^2 + d2*x^3; yB = vpa(yB, 3)
yC = a3 + b3*x + c3*x^2 + d3*x^3; yC = vpa(yC, 3)
%plotting
xx = 0:0.1:9;
yp1 = subs(yA, x, xx);
yp2 = subs(yB, x, xx);
yp3 = subs(yC, x, xx);
plot(xx, yp1, 'k-', xx, yp2, 'k--', xx, yp3, 'k-.');
axis([0, 8, -5, 10]); grid on; hold on;
plot(X1, Y1, 'ko');
scatter(X1, Y1, 20, 'r', 'filled');
text(0.5,1,'P_1');
text(2.5,2,'P_2');
text(5.5,-0.5,'P_3');
text(7.5,0,'P_4');
legend('A', 'B', 'C'); xlabel('x'); ylabel('y');
hold off;
%verification
vrf1 = eval(subs(yA, x, x1)) %should return y1
vrf2 = eval(subs(yB, x, x2)) %should return y2
vrf3 = eval(subs(yC, x, x3)) %should return y3

Equations for segments : 
yA =
0.447*x^3 - 2.64*x^2 + 4.0*x + 1.0
yB =
0.0473*x^3 - 0.244*x^2 - 0.801*x + 4.2
yC =
- 0.189*x^3 + 3.29*x^2 - 18.5*x + 33.7
vrf1 =
     1
vrf2 =
    2.0000
vrf3 =
   9.4587e-11

从图可以观察到， 和 是三个不同的样条，具有不同的形状，这是预期的，因为它们具有不同的方程。然而，连续性条件已经限制它们以一种方式形成了一个单一的平滑曲线，仅在区间到内。超出它们各自的区间后，它们分别发散成不同的轨迹。

• eval: 评估一个表达式
• legend: 使用文本字符串在图表中指定不同的颜色或线型

1.9 分段样条曲线（参数形式）

为了结束本章并解释一个非常重要的领域转换概念，我们将讨论更复杂的分段样条曲线形式，涉及参数方程。由于现在每个样条段都需要用两个方程（ vs. 和 vs. ）来表示，未知数的数量实际上增加到了24个。然而，使用约束条件形成方程的基本思想仍然保持不变，可以将其视为上一节中解释的思想的扩展版本。为了降低情况的复杂性，在这里使用了一个简化的假设： vs. 的关系是线性的，而不是三次的。这将减少未知数的总数，使情况更容易理解。然而，讨论这一部分的主要原因是要让读者意识到，当在空间域中指定条件并且需要在参数域中获得所需的方程时（或反之亦然），不能简单地将值替代到约束方程中（如在前一种情况下所做的那样），而需要将其从一个域转换到另一个域（如下所示将进行说明）。

与之前一样，假设给定点的坐标为和。让三个三次曲线段分别指定为和。这次我们假设关系是线性的，所以起始方程的形式为：。

因此，各自分段样条的方程变为：

对于样条，在起始点时，在结束点时。这提供了解：

对于样条，在起始点时，在结束点时。 这提供了解：

对于样条，在起始点时，在结束点时。这提供了解：

这些解得出得相对较容易，因此现在将注意力转向更困难的问题，即处理三次 vs. 关系。在这里，预期的起始关系形式是：

分段样条曲线的方程如下：

这些方程的一阶导数如下：

这些方程的二阶导数如下：

对于第一个与 连续性条件相关的约束，需要注意在样条曲线 中，在起始点 时，，在结束点 时，。这提供了以下解：

对于样条曲线 ，在起始点 时，，在结束点 时，。这提供了以下解：

对于样条曲线 ，在起始点 时，，在结束点 时，。这提供了以下解：

第二个要遵守的约束是连续性条件， 该条件要求在形成平滑曲线时，各个样条段的斜率在它们的交点处相等。在这种情况下，点 处的样条 的斜率应该等于点 处的样条 的斜率，即 。但斜率在空间域（空间中的物理斜率）相等，而上面所示的导数方程是在参数域（ 对 ）中计算的，因此它们不能简单地相等，而需要一种从一个域到另一个域的转换。为了得到转换因子，需要利用微分的链式法则。

根据微分的链式法则：

这明确了所需的转换因子：空间域中的导数等于参数域中的导数乘以特定曲线段的缩放因子 。类似地，空间域中的二阶导数等于参数域中的二阶导数乘以缩放因子 。

将这些乘数代入 约束方程中： 在 处的斜率 = 在 处的斜率 [在空间域中]

重新排列：

在 处的斜率 = 在 处的斜率 [在空间域中]

重新排列：

以类似的方式，使用域转换乘数来处理 约束方程： 在 处的曲率 = 在 处的曲率 [在空间域中]

重新排列并两边除以2：

在 处的曲率 = 在 处的曲率 [在空间域中]

重新排列并两边除以2：

还需要为端点条件考虑域转换因子。 设 为段 的起点斜率， 为段 的终点斜率。然后：

将所有九个方程都写成 的形式：

这个方程的解是：。

示例

找到通过 和 的分段参数化立方曲线方程。首尾点的斜率分别为2和1。假设 和 之间存在线性关系。

由于 关系是线性的，让起始方程的形式为：。分别为不同部分的方程：

将给定的点代入上述方程并解出未知系数：

所需的 关系是：

由于 关系是三次的，让起始方程的形式为：。分别为不同部分的方程：

将给定的值代入约束矩阵：

解出并代入起始方程中得到：

验证：对于 ： 产生 ， 产生 ；对于 ： 产生 ， 产生 ；对于 ： 产生 ， 产生

这些方程表示了通过给定点 和 的分段参数化立方曲线的 关系。这些方程通过解未知系数得到，并通过验证确保它们满足了给定的点和斜率条件。

clear all; clc;
syms t;
x1 = 0; y1 = 1;
x2 = 2; y2 = 2;
x3 = 5; y3 = 0;
x4 = 8; y4 = 0;
s1 = 2; s2 = 1;
C=[ 1, 1, 1, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 1, 1, 1, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 1, 1, 1;
 (x3-x2), 2*(x3-x2), 3*(x3-x2), -(x2-x1), 0, 0, 0, 0, 0;
 0, 0, 0, (x4-x3), 2*(x4-x3), 3*(x4-x3), -(x3-x2), 0, 0;
 0, (x3-x2)^2, 3*(x3-x2)^2, 0, -(x2-x1)^2, 0, 0, 0, 0;
 0, 0, 0, 0, (x4-x3)^2, 3*(x4-x3)^2, 0, -(x3-x2)^2, 0;
 1, 0, 0, 0, 0, 0, 0, 0, 0;
 0, 0, 0, 0, 0, 0, 1, 2, 3
];
G = [y2-y1, y3-y2, y4-y3, 0, 0, 0, 0, s1*(x2-x1), s2*(x4-x3)];
A = inv(C)*G';
aA = y1; bA = A(1); cA = A(2); dA = A(3);
aB = y2; bB = A(4); cB = A(5); dB = A(6);
aC = y3; bC = A(7); cC = A(8); dC = A(9);
fprintf('Equations of segments :\n')
xA = x1 + (x2-x1)*t; xA = vpa(xA)
yA = aA + bA*t + cA*t^2 + dA*t^3; yA = vpa(yA, 3)
xB = x2 + (x3-x2)*t; xB = vpa(xB)
yB = aB + bB*t + cB*t^2 + dB*t^3; yB = vpa(yB, 3)
xC = x3 + (x4-x3)*t; xC = vpa(xC)
yC = aC + bC*t + cC*t^2 + dC*t^3; yC = vpa(yC, 3)
%plotting
tt = linspace(0,1);
xa = subs(xA, t, tt);
ya = subs(yA, t, tt);
xb = subs(xB, t, tt);
yb = subs(yB, t, tt);
xc = subs(xC, t, tt);
yc = subs(yC, t, tt);
subplot(131); plot(tt,xa, 'k-', tt, xb, 'k--', tt, xc, 'k-.');
xlabel('t'); ylabel('x'); title('t - x'); axis square;
subplot(132); plot(tt,ya, 'k-', tt, yb, 'k--', tt, yc, 'k-.');
xlabel('t'); ylabel('y'); title('t - y'); axis square;
subplot(133); X = [x1 x2 x3 x4]; Y = [y1 y2 y3 y4];
plot(xa,ya,'k-', xb,yb, 'k--',xc,yc,'k-.', X, Y, 'ko'); hold on;
scatter(X, Y, 20, 'r', 'filled'); grid;
text(0.5,1,'P_1');
text(2.5,2,'P_2');
text(5.5,0,'P_3');
text(7,0,'P_4');
legend('A', 'B', 'C');
xlabel('x'); ylabel('y'); title('x - y'); axis square;
hold off;

Equations of segments :
xA =
2.0*t
yA =
1.16*t^3 - 4.16*t^2 + 4.0*t + 1.0
xB =
3.0*t + 2.0
yB =
0.803*t^3 - 1.54*t^2 - 1.26*t + 2.0
xC =
3.0*t + 5.0
yC =
1.07*t^3 + 0.868*t^2 - 1.93*t