贝塞尔曲线与de Casteljau算法
贝塞尔曲线与de Casteljau算法
一、简介
前言
在贝塞尔曲线原理、推导及Matlab实现这篇文章中，详细地介绍了贝塞尔曲线的原理、推导过程以及Matlab实现。文章中计算贝塞尔曲线所采用的方法是定义法，该方法简洁易懂，不过其中的二项式系数\(\left( \begin{matrix} n \\ i \end{matrix} \right) = \frac{n !}{i ! \cdot (n - i) !}\)看似简单，实则在高阶时会带来非常大的计算量。对于程序设计来说，这样的开销显然是不太能够接受的。所以，本文将会引用一个新的算法——de Casteljau算法。
趣闻
de Casteljau算法是雪铁龙公司一位名为Paul de Faget de Casteljau的工程师在1963年的发明，理论上来说Paul de Faget de Casteljau才是贝塞尔曲线的“第一位发明者”。他发现对于汽车的几何外形设计，与其定义一条以长度为参数的曲线，不如定义一系列控制点作为控制曲线的参数，通过控制点的位置变化来影响曲线形态的变化。这一举动对于汽车设计产业而言无疑是创新性的，所以在初期遭到了公司内部的质疑，但最终这位年轻的工程师还是说服公司采纳了他的算法。该算法经采用后，最终取得了成功，雪铁龙公司要求Casteljau对de Casteljau算法严格保密，直到8年之后，这位工程师的成果才为世人所知。
同期，尽管雪铁龙公司对de Casteljau算法进行了严密的保护，但其竞争对手雷诺公司还是了解到Casteljau的工作，并且也开始对该方向进行研究。1966年，雷诺公司的工程师Pierre Bézier和他的团队独立重现了Casteljau的研究成果，并将这条曲线正式命名为Bézier curve，即我们熟知的贝塞尔曲线。与雪铁龙公司的态度相反，雷诺公司允许Bézier发表他的研究成果，所以更为世人所熟知的是贝塞尔曲线而不是de Casteljau算法。
二、原理
相比定义法在求取高阶贝塞尔曲线时带来的大量计算和数值不稳定性，de Casteljau算法又快又稳定，并且更加贴近贝塞尔曲线的特性。
de Casteljau算法可以看作是重复线性插值（repeated linear interpolation）的过程。
重复线性插值
重复线性插值是指在数据点之间使用线性插值来估算其他数据点的值，通常是在时间序列或连续数据中使用的一种插值方法。
而线性插值是一种简单的插值方法，它假定在两个已知数据点之间的值变化是线性的。在重复的线性插值中，这个过程会在数据点之间不断重复进行，以填补整个数据序列中的缺失值或生成更密集的数据点。
重复的线性插值的过程可以通过下述的三张过程图来直观地感受到：
拥有上述知识后，开始着手构建贝塞尔曲线。
一次曲线（线性）
一次曲线贝塞尔曲线由两个端点\(P_0\)和\(P_1\)组成，曲线由线段直接连接这两个点，是一条最简单的线性曲线。
\(P(t)\)描述一条由\(P_0\)至\(P_1\)的直线，即使用了标准的线性插值：
\[P(t) = (1 - t) P_0 + t P_1, t \in [0, 1] \tag1 \]
二次曲线
二次贝塞尔曲线由定点\(P_0, P_1, P_2\)确定，通过三次线性插值来构建：
\[\begin{matrix} P_0^1 = (1 - t) P_0 + t P_1 \\ P_1^1 = (1 - t) P_1 + t P_2 \\ P(t) = P_0^2 = (1 - t) P_0^1 + t P_1^1 \end{matrix} \tag2 \]式\((2)\)也可以表述为：
- \(P_0, P_1\)之间的连续点\(P_0^1\)，描述一条线性贝塞尔曲线；
- \(P_1, P_2\)之间的连续点\(P_1^1\)，描述一条线性贝塞尔曲线；
- \(P_0^1, P_1^1\)之间的连续点\(P_0^2\)，描述一条二次贝塞尔曲线。
![二次贝塞尔曲线演示动画，t在[0,1]区间.gif](/i/l/?n=24&i=blog/3078251/202407/3078251-20240704170630149-667424969.gif)
三次曲线
三次曲线由定点\(P_0, P_1, P_2, P_3\)确定，通过六次线性插值来构建：
\[\begin{matrix} P_0^1 = (1 - t) P_0 + t P_1 \\ P_1^1 = (1 - t) P_1 + t P_2 \\ P_2^1 = (1 - t) P_2 + t P_3 \\ P_0^2 = (1 - t) P_0^1 + t P_1^1 \\ P_1^2 = (1 - t) P_1^1 + t P_2^1 \\ P(t) = P_0^3 = (1 - t) P_0^2 + t P_1^2 \\ \end{matrix} \tag3 \]式\((2)\)也可以表述为：
- \(P_0, P_1\)之间的连续点\(P_0^1\)，描述一条线性贝塞尔曲线；
- \(P_1, P_2\)之间的连续点\(P_1^1\)，描述一条线性贝塞尔曲线；
- \(P_2, P_3\)之间的连续点\(P_2^1\)，描述一条线性贝塞尔曲线；
- \(P_0^1, P_1^1\)之间的连续点\(P_0^2\)，描述一条二次贝塞尔曲线；
- \(P_1^1, P_2^1\)之间的连续点\(P_1^2\)，描述一条二次贝塞尔曲线；
- \(P_0^2, P_1^2\)之间的连续点\(P_0^3\)，描述一条三次贝塞尔曲线。
高次曲线
为构建更高次的曲线，则通过更多次的线性插值来构建，这就是重复线性插值。
如，四次曲线通过十次线性插值：
![四次贝塞尔曲线演示动画，t在[0,1]区间.gif](/i/l/?n=24&i=blog/3078251/202407/3078251-20240704170713632-721194624.gif)
五次曲线通过十五次线性插值：
![五次贝塞尔曲线演示动画，t在[0,1]区间.gif](/i/l/?n=24&i=blog/3078251/202407/3078251-20240704170740830-603385477.gif)
推广
对上文总结得出，只要对控制多边形上的各点进行线性插值，就可以得到更少的点和线段。
比如有四个控制点的三次曲线，对四个点进行线性插值，得到三个点和对应的三条线段；对三个点进行线性插值，得到两个点和对应的两条线段；最终对两个点进行线性插值，得到一个可以描述三次曲线的目标点。
该过程可以用下图自下而上的金字塔构建方法形象地给出：
故对于任意一点，它都可以由下式表示：
\[P_i^j(t)= (1 - t) P_i^{j - 1} + t P_{i + 1}^{j - 1}, i = 0, 1, \cdots , n; j = 0, 1, \cdots , i \tag4 \]特性
观察贝塞尔曲线，从几何意义上来看，当曲线参数\(t = 0\)时，对应曲线的第\(0\)个（首个）控制点，即\(P(0) = P_0\)；当\(t = 1\)时，对应曲线的第\(n\)个（末尾）控制点，即\(P(1) = P_n\)。
这便是贝塞尔曲线的端点插值特性。
三、应用
Matlab实现
main.m
%% 初始化
clc
clear
%% 输入
n = input("请输入曲线阶数：\n");
P = pos_input(n); % 获取控制点坐标
%% 调用贝塞尔曲线计算函数
t = linspace(0, 1, 1000); % 生成一系列t参数，范围在[0, 1]
curve = BC_deCasteljau(P, t); % 计算贝塞尔曲线
%% 绘制贝塞尔曲线
plot(curve(:, 1), curve(:, 2), '-b'); % 绘制贝塞尔曲线
hold on;
plot(P(:, 1), P(:, 2), '-ro'); % 绘制控制点
title(sprintf('贝塞尔曲线 (阶数： %d)', n));
xlabel('X');
ylabel('Y');
grid on;
legend('贝塞尔曲线', '控制点');
hold off;
BC_deCasteljau.m（计算贝塞尔曲线）
function B = BC_deCasteljau(control_points, t)
% 通过de Casteljau算法计算贝塞尔曲线
n = size(control_points, 1) - 1; % 通过控制点个数获取阶数
B = zeros(length(t), size(control_points, 2)); % 初始化曲线矩阵
for k = 1 : length(t) % 每个参数t都对应一个曲线坐标
curve = control_points; % 复制控制点坐标
for i = 1 : n % 通过线性插值从1次到n次，如同金字塔第1层至第n层
for j = 1 : n - i + 1 % i次曲线对应的点数，如同金字塔第n层长度
% i次曲线第j点由（i - 1）次曲线的第j点和第j + 1点线性插值得到
curve(j, :) = (1 - t(k)) * curve(j, :) + t(k) * curve(j + 1, :);
end
end
% 取n次曲线的第1个（也是唯一一个）坐标点作为第k个t参数对应的曲线坐标
% 如同金字塔塔尖的点
B(k, :) = curve(1, :);
end
end
pos_input.m（选择输入控制点 OR 随机生成）
function P = pos_input(I)
% 输入贝塞尔曲线控制点
P = ones(I, 2); % 初始化坐标矩阵
switch input("是否使用随机点？（范围[-100, 100], y = Yes, n = No）：\n", 's')
case 'y'
P = generate_random_points(I + 1);
case 'n'
for i = 0 : I
P(i + 1, :) = input("请输入输入坐标点P" + num2str(i) + "（以一维矩阵形式输入）: \n");
end
otherwise
disp("请选择 y 或者 n噢~\n");
return
end
end
generate_random_points.m
function random_points = generate_random_points(n)
% 生成n个随机点的坐标，范围为[-100, 100]
random_points = zeros(n, 2); % 初始化存储点坐标的数组
% 生成随机坐标
for i = 1 : n
% 生成x坐标
random_points(i, 1) = rand * 200 - 100; % [-100, 100]
% 生成y坐标
random_points(i, 2) = rand * 200 - 100; % [-100, 100]
end
end
曲线绘制
