Bezier曲线及其de casteljau算法 matlab实现

最新推荐文章于 2024-01-18 11:30:09 发布

mw_1422102031

最新推荐文章于 2024-01-18 11:30:09 发布

阅读量3.7k

点赞数 11

分类专栏：数学matlab程序文章标签：线性代数几何学算法

本文链接：https://blog.csdn.net/mw_1422102031/article/details/107498738

版权

数学matlab程序专栏收录该内容

14 篇文章 7 订阅

订阅专栏

一、Bezier曲线的定义

给定n+1个空间向量 $P_i\in R^3 (i= 0,1,...,n)$ ，称n次参数曲线段
$\sum^n_{i=0}P_iB^n_i(t) ,0\leq t \leq1 .（1）$

为一条n次Bezier曲线. $P_i$ 称为控制顶点，依次用直线段连接相邻的两个 $P_i$ 所得的n边折线多边形称为Bezier多边形或者控制多边形.
下图为3次Bezier曲线，其中折线为相应的控制多边形.
n=3
下面给出matlab代码：

function bezier
%bezier曲线画图
px=[0,1,2,3];
py=[0,4,1,5];                              %控制顶点(0,0),(1,4),(2,1),(3,5)
M = 40;
hx = 1/M;                                    %将[0,1]区间M等分
x = (0:hx:1)';
n=length(px)-1;
Px = 0;
Py = 0;
for i = 1:n+1
    Px = Px + px(i)*B(x,n,i-1)        %将控制顶点与Bernstein基函数相乘得到bezier曲线
    Py = Py + py(i)*B(x,n,i-1)
end
plot(Px,Py,'r',px,py,'b')

end

function y = B(x,n,i)
y = k(n,i).*(x.^i).*((1-x).^(n-i));
end

function y = k(n,i)
y1 = factorial(n);            %n的阶乘
y2 = factorial(i)*factorial(n-i);
y = y1/y2;
end

二、Bezier曲线的性质

几何不变性和仿射不变性. 因为Bernstein基函数满足单位分解性，对Bezier曲线（1）式进行仿射变换，即用线性变换M和平移c作用，得到新的曲线
$\ P^*(t) = MP(t)+c = M\sum^{n}_{i = 0}P_iB^{n}_i(t) + c\sum^{n}_{i = 0}B^{n}_i(t) \\= \sum^{n}_{i = 0}MP_iB^{n}_i(t) + \sum^{n}_{i = 0}cB^{n}_i(t) = \sum^n_{i=0}P^*_iB^n_i(t)$
为对原Bezier曲线（1）式的控制顶点 $P_i$ 进行仿射变换后得到的新的控制顶点 $P^*_i(i=0,1,...,n)$ 对应的Bezier曲线.这说明，Bezier曲线不依赖于坐标系的选取，是几何不变的.
凸包性质. 从几何上看，意味着bezier曲线落在了控制多边形的凸包之中。
端点性质.
$P(0) = P_0,P(1) = P_n$
$P'(0)=n(P_1-P_0),P'(1) = n(P_n-P_{n-1}).$
移动n次Bezier曲线（1）式的第i个控制顶点 $P_i$ ,将对整条曲线产生影响，其中参数为 $t=\frac {i}{n}$ 的那个点 $P(\frac {i}{n})$ 处的影响最大.这是因为相应的基函数 $B^n_i(t)$ 在 $t=\frac {i}{n}$ 达到最大值。

三、Bezier曲线的de Casteljau算法和几何作图法

在构造出一条Bezier曲线之后，常常需要计算曲线上对应某个参数值 $t^*$ 的点 $P(t^*)$ .数值计算的方法由直接代入参数值即可得到所求点的坐标.借助于Bezier基函数，我们还可以采用几何作图的发放计算Bezier曲线上的点.
由Bernstein基函数的递推性质
$B^n_i(t) = (1-t)B^{n-1}_i(t)+tB^{n-1}_{i-1}(t),\\ \\B^{n-1}_{-1}(t)=B^{n-1}_{n}(t)=0, i=0,1,...,n$
可知
$\sum^n_{i=0}P_iB^n_i(t)$
$=P_0(1-t)B^{n-1}_0(t)+ \sum^{n-1}_{i=0}[P_iB^{n-1}_i(t)+t B^{n-1}_{i-1}(t)] +P_n(t) B^{n-1}_{n-1}(t) \\=(1-t) \sum^{n-1}_{i=0}P_iB^{n-1}_i(t)+t \sum^{n-1}_{i=0}P_{i+1}B^{n-1}_i(t) \\= \sum^{n-1}_{i=0}[(1-t)P_i+tP_{i+1}] B^{n-1}_i(t)$

这说明，一条n次Bezier曲线可以表示为分别由前后n个控制顶点决定的两条n-1次Bezier曲线的线性组合.由此可得Bezier曲线上某一点递归求值的 de Casteljau算法

$\begin{cases}P^0_i(t)\equiv P^0_i =P_i,i = 0.1....,n, \\ P^r_i(t) = (1-t)P^{r-1}_i(t)+tP^{r-1}_{i+1}(t) , \\ r=0,1,...,n ,i = 0.1....,n-r. \end{cases}$
于是
$\sum^{n-1}_{i=0}P^1_iB^{n-1}_i(t)=... = \sum^{n-r}_{i=0}P^r_iB^{n-r}_i(t)=... =P^n_0(t),$
$P^n_0(t)$ 即为所求的 $P (t)$ .若记 $P^0=(P_0,P_1,\cdots,P_n)^T,P^r(t)= (P_0^r,P_1^r,\cdots,P_{n-r}^r)^T$ ,则可以把de Castaljau算法表示为
$P^r(t)= M_r(t) \cdots M_2(t)M_1(t)P^0,$
其中
$M_r(t)=\begin{pmatrix} 1-t & t & 0 & \cdots & 0 & 0 \\ 0 & 1-t & t & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1-t & t & 0 \\ 0 & 0 & \cdots & 0 &1-t & 0 \\ \end{pmatrix}$
是一个 $\times (n-r+2)$ 阶矩阵.

四、de Casteljau算法的matlab程序实现：


px=[0,1,2,3];
py=[0,4,1,5];                             %控制顶点(0,0),(1,4),(2,1),(3,5)
M=100;                                     %将[0,1]区间M等分
hx = 1/M;                                  %每个小区间长度
n=length(px)-1;
xx=[];yy=[];
for x = 0:hx:1
    xx1 = px';
    yy1 = py';
    for i = 1:n
        M(i,i) = 1-x;
        M(i,i+1) = x;                             %生成迭代矩阵
    end
    for i = 1:n
        M = M(1:n+1-i,1:n+2-i);
        xx1 = M*xx1;
        yy1 = M*yy1;
    end
     xx(end+1)=xx1(1);                     %xx存储曲线的横坐标
     yy(end+1)=yy1(1);                     %yy存储曲线的纵坐标
end
plot(xx,yy,px,py,px,py,'b*');
hold on;
%以下是当x = 0.3时的具体实现
xx1 = px';
yy1 = py';
x = 0.3;
for i = 1:n
    M(i,i) = 1-x;
    M(i,i+1) = x;
end
for i = 1:n
    M = M(1:n+1-i,1:n+2-i);
    xx1 = M*xx1;
    yy1 = M*yy1;
    plot(xx1,yy1,xx1,yy1,'o');hold on  %将每一步得到的结果画出来
end

程序结果如下：
在这里插入图片描述
接下来附上我第一次上课时交的编程作业（傻乎乎的程序）：

px=[0,1,2,3];py=[0,4,1,5];
A=[px;py];N=size(A,2)-1;
x=zeros(1,N+1);y=zeros(1,N+1);
xx=[];yy=[];
for t=0:0.001:1;
      for i=1:N
        for j=0:N-i
          if i==1
            x(j+1) = (1 - t)*A(1,j+1) + t*A(1,j+2);
            y(j+1) = (1 - t)*A(2,j+1) + t*A(2,j+2);
            continue;
          end
         x(j+1) = x(j+1) * (1 - t) + x(j + 2) * t;
         y(j+1) = y(j+1) * (1 - t) + y(j + 2) * t;
    end
  end
  xx(end+1)=x(1);
  yy(end+1)=y(1);
end
plot(A(1,:),A(2,:),'mo',A(1,:),A(2,:),'b','linewidth',2);
hold on;grid on;
axis tight;  
plot(xx,yy,'r','linewidth',2);
hold on;
t=0.3;
for j=0:N-1
    x(j+1) = (1 - t)*A(1,j+1) + t*A(1,j+2);
    y(j+1) = (1 - t)*A(2,j+1) + t*A(2,j+2);
    plot(x(j+1),y(j+1),'-go','linewidth',2);hold on;
end
plot(x(1:N),y(1:N),'k-','linewidth',2);hold on;
for j=0:N-2
   x(j+1) = x(j+1) * (1 - t) + x(j + 2) * t;
   y(j+1) = y(j+1) * (1 - t) + y(j + 2) * t;
   plot(x(j+1),y(j+1),'-co','linewidth',2);hold on;
end
plot(x(1:N-1),y(1:N-1),'g-','linewidth',2);hold on;
for j=0:N-3
   x(j+1) = x(j+1) * (1 - t) + x(j + 2) * t;
   y(j+1) = y(j+1) * (1 - t) + y(j + 2) * t;
   plot(x(j+1),y(j+1),'-ko','linewidth',2);hold on;
end

程序结果为：
在这里插入图片描述

五、参考文献：计算几何教程；王仁宏，李崇军，朱春钢编著

mw_1422102031

关注

11
点赞
踩
43

收藏

觉得还不错? 一键收藏
0
评论
Bezier曲线及其de casteljau算法 matlab实现

一、Bezier曲线的定义      给定n+1个空间向量Pi∈R3(i=0,1,...,n)P_i\in R^3 (i= 0,1,...,n)Pi∈R3(i=0,1,...,n)，称n次参数曲线段P(t)=∑i=0nPiBin(t),0≤t≤1P(t) = \sum^n_{i=0}P_iB^n_i(t) ,0\leq t \leq1P(t)=i=0∑nPiBin(t),0≤t≤1为一条n次Bezier曲线.PiP_iPi
复制链接

扫一扫