平面中的样条线

---

这个示例展示了如何使用曲线拟合工具箱指令 spmak, spcrv, cscvn 构造平面中的样条曲线。这包括切线的绘制以及曲面包围区域的计算。

一个简单的样条曲线

曲线拟合工具箱可以处理向量值样条。 d 向量值单变量样条曲线提供了 d 空间中的曲线。 在这种模式下，d = 2 是最常见的，因为它给出了平面曲线。

这是一个示例，其中构造并绘制了具有二维系数的样条曲线。

using TyPlot
using TyCurveFitting
using TyMath
knots = [1; 1:9; 9]'
curve = spmak(knots, transpose(repeat([0 0; 1 0; 1 1; 0 1], 2, 1)))
t = LinRange(2, 8, 121)
value = fnval(curve, t')
plot(value[1, :], value[2, :], linewidth=2)
axis([-0.4 1.4 -0.2 1.2])
axis("equal")
title("A Spline Curve");
hold("on")

一个谨慎的词

您可能注意到这个示例并没有使用 fnplt 绘制曲线，而是使用 fnval 绘制了曲线上的一些点，这是代码：

t = LinRange(2,8,121) values = fnval(curve,t) plot(values[1,:],values[2,:],linewidth=2)

对这个样条曲线直接使用 fnplt 会得到下图中的红色部分。

fnplt(curve, "r", 0.5)
title("The Full Spline Curve, in Red")

解释呢？

样条是 4 阶的，节点序列中的末端节点

knots

knots = 1×11 adjoint(::Vector{Int64}) with eltype Int64:
 
1  1  2  3  4  5  6  7  8  9  9       

重数仅为 2。因此，这个节点序列的所有 4 阶 B 样条在基本区间的端点处为 0，这使得曲线在原点开始与结束。

一个补救措施

由于在这种情况下，我们只对与参数区间 [3 .. 7] 对应的曲线段感兴趣，因此我们可以使用 fnbrk 提取该部分，然后使用 fnplt 将其绘制成黄色就很容易了。

mycv = fnbrk(curve, [3, 7])
fnplt(mycv, "y", 2.5)
title("The Spline Curve of Interest, in Yellow")

这条曲线包围的区域

由于您现在有一个样条曲线，即 mycv，它描述了曲线（仅此而已），您可以轻松计算此闭合曲线所包围的面积，如下所示。

area = diff(fnval(fnint(fncmb(fncmb(mycv, [0 1]), '*', fnder(fncmb(mycv, [1 0])))), fnbrk(mycv, "interval")), dims=2)

area = 1×1 Matrix{Float64}:

 -0.8333333333333334

稍加努力，您就可以将其识别为积分的值

int y(t) d(x(t)) = int y(t) Dx(t) dt

在样条曲线 mycv 的基本区间上，其中 (x(t),y(t)) := fnval(mycv,t) 曲线上对应于参数值 t 的点。 这里，fncmb(mycv,[1,0]), fncmb(mycv,[0,1]) 描述了样条曲线的两个分量，即标量值样条 x 和 y。

此外，曲线大致是一个半径为 1/2 的圆。 因此，您会期望一个大致的区域，

pi/4

ans = 0.7853981633974483

添加一些切向量

我们重新绘制曲线，并在某些点绘制曲线的切向量。

hold("off")
fnplt(mycv, "y", 2.5)
hold("on")
t = (3:0.4:6.2)';
cv = fnval(curve, t);
cdv = fnval(fnder(curve), t);
quiver(cv[1, :], cv[2, :], cdv[1, :], cdv[2, :], length=0.5, color="r");
title("A Spline Curve With Some Tangents")
axis([-0.4 1.4 -0.2 1.2])
axis("equal")

曲线与直线的交点

如果您想确定这条样条曲线与直线 y = x 的交点，下面的代码会将它们提供给您，并在曲线内绘制该直线的线段：

cuts = fnval(mycv, mean(fnzeros(fncmb(fncmb(mycv, [0 1]), '-', fncmb(mycv, [1 0]))), dims=1))
plot(cuts[1, :], cuts[2, :], "y", linewidth=2.5)
hold("off")
title("A Spline Curve With Some Tangents and a Cut Across")

SPRV：控制多边形和相应的样条曲线

样条曲线广泛用于插图的生成，其中只需要某种粗略想象形状的平滑曲线。 为此，曲线拟合工具箱包含一个特殊命令 spcrv，它可以独立于工具箱的其余部分使用。

给定平面中的点序列和可选的 k 阶，spcrv 通过重复中点节点插入生成 k 阶样条曲线，其控制多边形由给定序列指定。

下图显示了这样一个控制多边形，以及对应的 3 阶样条曲线。

points = transpose([0 0; 1 0; 1 1; 0 2; -1 1; -1 0; 0 -1; 0 -2])
values = spcrv(points, 3)

plot(points[1, :], points[2, :], "k")
axis([-2 2.25 -2.1 2.2])
hold("on")
plot(values[1, :], values[2, :], "r", linewidth=1.5)
legend(["Control Polygon", "Quadratic Spline Curve"],loc="southeast")

请注意，曲线在其中点接触控制多边形的每个部分，并遵循控制多边形所勾勒的形状。

提升阶

提高阶数 k 将使曲线远离控制多边形并使其更平滑，但也更短。 在这里，我们添加了相应的 4 阶样条曲线。

value4 = spcrv(points, 4)
plot(value4[1, :], value4[2, :], "b", linewidth=2)
legend(["Control Polygon", "Quadratic Spline Curve", "Cubic Spline Curve"],loc="southeast")

CSCVN

另一方面，要获得插值曲线，可以使用 cscvn 命令，它提供参数化的“自然”三次样条曲线。

fnplt(cscvn(points), "g", 1.5)
legend(["Control Polygon", "Quadratic Spline Curve", "Cubic Spline Curve", "Interpolating Spline Curve"],loc="southeast")

通过在第二个控制点 (1,0) 附近添加点 (.95,-.05)，我们可以创建一个插值样条曲线，该曲线在该处转动得更快。

np = size(points, 2);
fnplt(cscvn([points[:, 1] [0.95; -0.05] points[:, 2:np]]), "m", 1.5);
plot(0.95, -0.05, "*");
legend(["Control Polygon", "Quadratic Spline Curve",
   "Cubic Spline Curve", "Interpolating Spline Curve",
   "Faster Turning Near (1,0)"],loc="southeast")
hold("off")

RSCVN

您还可以获得由圆弧组成的切线连续曲线，该圆弧穿过平面中给定的点序列，并且可以选择在这些点处与给定的法线方向正交。 命令 rscvn 提供了这样的曲线。

例如下面生成一个圆

c = rscvn([-1 1 -1; 0 0 0], [1 1; 0 0])
fnplt(c)
axis([-1.05 1.05 -1.05 1.05])
axis("equal")
axis("off")

c 是仅由两部分组成的二次有理样条，如以下命令所示。

form, order, breaks = fnbrk(c, "f", "o", "b")

form = "rBform"

order = 3

breaks = 1×3 Matrix{Int64}:
 0  2  4

这样的工具仅仅使用一些点就可以生成让人震惊的结果。比如，下图是一个 Belfast 的 Ulster 博物馆中 Bronze Triskele Medallion 展品的设计，这是很久之前使用圆弧片段完成的。

pp = [zeros(1, 7); 5.4 3 6.9 2.75 2.5 0.5 5]
alpha = 2 * pi / 3
ca = cos(alpha)
sa = sin(alpha)
c = [ca sa; -sa ca]
d = [0 0 0.05 -0.05; 1 -1 0.98 0.98]
d = [d c * d]
yin = rscvn([pp[:, [7; 1:3]] c * pp[:, 3:4] pp[:, 3]], d[:, [1, 2, 1, 4, 7, 5, 1]])
fnplt(yin)
hold("on")
fnplt(fncmb(yin, c))
fnplt(fncmb(yin, c'))
yang = rscvn([pp[:, 6] -pp[:, 6] pp[:, 5] c * pp[:, 4]], [d[:, [2, 1, 1]] c[:, 2]])
fnplt(yang), fnplt(fncmb(yang, c)), fnplt(fncmb(yang, c'))
axis([-7.2 7.2 -7.2 7.2])
axis("equal")
axis("off")
hold("off")