通过给定一组数据点并反求控制点的NURBS曲线插值生成Matlab编程实例

不想不努力的菜菜

已于 2024-01-17 13:38:36 修改

阅读量2.2w

点赞数 92

分类专栏： Matlab 文章标签： nurbs matlab

于 2020-08-02 23:10:12 首次发布

本文链接：https://blog.csdn.net/shs1992shs/article/details/107748375

版权

Matlab 专栏收录该内容

1 篇文章 5 订阅

订阅专栏

0.引言

非均匀有理B样条(NURBS)是几何造型的一个十分有用的工具，以其良好的特性越来越受到工程界的重视。它可以用统一的数学形式表示规则曲面与自由曲面；有操纵控制顶点及权因子，为各种形状设计提供了充分的灵活性，有强有力的几何配套技术，能用于设计、分析与处理等各个环节，是非有理B样条形式以及有理、非有理贝塞尔形式良好的推广。NURBS方法已成为IGES标准、STEP标准中描述产品几何形状的唯一方法。
NURBS曲线曲面在实际应用中可以分为两种形式：一种是已知控制点求解曲线曲面上的点，称为正算问题。另一种情况是已知曲线曲面上一组数据点，即型值点，求解曲线曲面的控制点，称为反算问题。在实际应用中，常常是给出一组离散的型值点，要求构造通过该型值点的曲线曲面，即所谓的曲线曲面插值。

1.NURBS曲线方程

NURBS曲线方程的表达式，有理分式表达：

式中：
在这里插入图片描述

其中， $w_i$ ( $i$ =0,1,…,n)称为权或者权因子，分别与控制顶点 $d_i(i=0,1,.....,n)$ 相对应。首末权因子 $w_0$ ， $w_n > 0$ ，其余 $w_i \ge 0$ ，并且连续 $k$ 个权因子不能同时为零。顺序连接 $d_i$ 形成控制多边形。 ${N}_{i,k}(u)$ 是由节点矢量 $U=[u_0,u_1,\ldots,u_{n+k+1}]$ 计算确定的 $k$ 次规范 B 样条基函数， $i$ 表示序号， $k$ 是样条阶数， $k = 3$ ，即是三次曲线，。对于NURBS 曲线首尾数据点不相连的，也就是开曲线，通常取两端节点的重复度为 $r = k + 1$ ，即 $u_0=u_{1}=\ldots =u_k$ ， $u_{n+1}=u_{n+2}=\ldots =u_{n+k+1}$ 。在很多实际应用中，端节点值常取0、1，故NURBS 曲线的定义域可表示为 $u\in[u_k,u_{n+1}]\subset [0,1]$ 。由此可得，当 $n = k$ 时， $k$ 次NURBS 曲线即退化为 $k$ 次有理Bézier 曲线，其节点矢量与同次有理Bézier 曲线的端点具有相同的几何性质。若首末权因子均不为零，即子 $w_1，w_n > 0$ ，则曲线首末端点分别与控制多边形首末顶点重合，同时曲线在该点处与控制多边形首末边相切。

2.NURBS曲线插值生成

给定一组数据点，即型值点，生成通过这些型值点的的NURBS曲线，称为曲线的反算，一般包括三个步骤：（1）计算节点矢量；（2）计算边界条件；（3）反算控制顶点。

2.1 计算节点矢量

若要让一条 $K$ 次NURBS 曲线通过一组给定的型值点 $P_i(i=0,1,.....,n)$ ，不仅要保证曲线的首末端点与型值点重合，还应保证 $P_i$ 依次与构造曲线定义域内的节点 $u_{i+k}(i=0,1,.....,n)$ 一一对应。通常需要对型值点进行参数化处理，以确定型值点 $P_i$ 的参数值 $u_{i+k}(i=0,1,.....,n)$ 。具有 $n + 1$ 个型值点 $P_i$ 的 $k$ 次NURBS曲线将由 $n + 3$ 个控制顶点 $D_i(i=0,1,.....,n,n+1,n+2)$ 及其权因子 $w_i$ 以及节点矢量 $U=[u_0,u_1,\ldots,u_{n+k+3}]$ 定义。节点矢量最后的下标是 $n + k + 3$ ，是因为节点的数量要比控制顶点的数量多重复度 $r = k + 1$ 个，而控制顶点数量比型值点数量多2个，所以节点矢量最后下标是 $= n + 2 + k + 1 = n + k + 3$ （ $n + 1$ 个型值点，下标从0开始， $n$ 结束）。首末节点的重复度 $r = k + 1$ 。前 $k + 1$ 个节点取值为0，后 $k + 1$ 个节点取值为1。
引入符号 $\bigtriangleup$ ，将每个节点区间长度表示为： $\bigtriangleup_i=u_{i+1}-u_i$ 。常用的参数化方法有以下几种，以三次曲线（ $k = 3$ ）为例：
（1）均匀参数化法
使节点在参数轴上呈现等距离分布，也就是应用向前差分表示每个节点区间的长度，即 $\bigtriangleup_i=u_{i+1}-u_i=$ 正常数，为保证节点矢量在[0,1]的范围内，对节点规范化处理，保证 $u_0,u_n]=[0,1]$ ，具体的参数化方法为：
在这里插入图片描述
这种均匀参数化方法适用的类型仅为给定数据点的多边形各边（或弦长）接近相等的情形，具有一定的局限性。
（2）积累弦长参数化法
该方法克服了均匀参数化法的弊端，如若遇到数据点多边形各弦长分布不均匀时，应用该方法可以避免均匀参数化法出现的问题，并如实地表现出数据点依据多边形各弦长的分布情形，且能得到光顺性较好的插值曲线，是最常采用的参数化方法，该参数化方法满足：
在这里插入图片描述
（3）向心参数化法
这种参数化方法如实反映了数据点相邻弦线的折拐情况。

（4）修正弦长参数化法
这种参数化方法对绝对曲率较大，与实际弧长相比弦长又偏短的曲线段起到了修正的作用，且生成的插值曲线具有较好的光顺性。
在这里插入图片描述

2.2 边界条件的确定

下式表示的线性方程组由 $n + 1$ 个矢量方程组成，具有 $n + 3 个$ 未知控制顶点，上述方程组中的节点矢量 $U=[u_0,u_1,\ldots,u_{n+k+3}]$ 已确定，但未知顶点数大于方程数，因此求解该方程组，需要增加两个适当的边界条件，而常用的边界条件有：（1）切矢条件；（2）自由端点条件；（3）闭曲线条件。常采用切矢条件。以三次曲线（ $k = 3$ ）为例进行说明。在这里插入图片描述
（1）切矢条件
要求首末端点的切线方向固定:

式中： $\bigtriangleup_i=u_{i+1}-u_i$ ， $P^{'}_0$ ， $P^{'}_n$ 为首末端点切矢。
（2）自由端点条件
自由端点条件要求首末端点具有零曲率：

（3）闭曲线条件
要求曲线首末端点重合且二阶连续：在这里插入图片描述

2.3 控制顶点的反算

1、对于 $c^2$ 连续的NURBS 三次闭曲线，由于其首末数据点相重， $n - 1$ 个方程足以求出 $n + 1$ 个控制顶点，所以NURBS 三次闭曲线控制顶点反算的矩阵表达式为:
在这里插入图片描述
式中： $d_i$ 为控制顶点， $\bigtriangleup_i=u_{i+1}-u_i (i=0,1,\ldots,n-2)$
2、对于NURBS 三次开曲线，则还需要添加两个边界条件满足的附加方程。NURBS 三次开曲线两端点的重复度取4，所以NURBS 三次曲线的首末控制顶点就是其首末端的型值点，即： $d_0=P_0$ ， $d_{n+2}=P_n$ 。边界条件以切矢条件为例，应用矩阵的形式给出NURBS 三次开曲线控制顶点反算的表达式：
在这里插入图片描述
式中： $d_i$ 为控制顶点， $\bigtriangleup_i=u_{i+1}-u_i (i=0,1,\ldots,n)$ ，则：

在这里插入图片描述
求解上述线性方程组，即可得到所有未知的控制顶点。根据要求选取合适的权因子，即可最终确定NURBS 曲线的表达式。

3 Matlab编程实例

本程序既适用于二维平面曲线，也适用于三维空间曲线，且有良好的鲁棒性。
data1.txt

data2.txt

% 给定一组数据点，即型值点，反求其控制顶点，并构造通过这些型值点的2D或3D nurbs曲线，首末型值点与控制顶点重合，称为nurbs曲线插值
% 具有n个型值点的k次nurbs插值曲线，将有n+2个未知控制顶点，首末节点取重复度r=k+1,从而具有n+2+k+1个节点
% 三次Nurbs曲线的首末节点取重复度r=k+1=4，u1=u2=u3=u4=0,Un+3=Un+4=Un+5=Un+6=1
clear
pt=load('data2.txt'); 
[row,column]=size(pt);  %每行代表一个型值点，采用列向量表示x,y,z坐标，便于线性方程组求解
n=length(pt);  %型值点数量
k=3;  %样条阶数
U=zeros(1,n+k+3);  %节点矢量
%第一步：计算节点矢量********************
if(column == 2)   % 2D curve
    x=pt(:,1);
    y=pt(:,2);
else              % 3D curve
    x=pt(:,1);
    y=pt(:,2);
    z=pt(:,3);
end

%n个数据点，利用规范积累弦长法进行型值点参数化处理
temp=zeros(1,n-1);
for i=1:n-1  %n个数据点，n-1段弦长
   if(column == 2)    % 2D curve
   temp(i)=sqrt((x(i+1)-x(i))^2+(y(i+1)-y(i))^2);
   else               % 3D curve
   temp(i)=sqrt((x(i+1)-x(i))^2+(y(i+1)-y(i))^2+(z(i+1)-z(i))^2);
   end
end
sumtemp=sum(temp);%顺序相邻两型值点之间距离的和,即总弦长
for i=1:k+1      %前k+1个节点为0
    U(i)=0;
end
for i=n+k:n+k+3  %后k+1个节点为1
    U(i)=1;
end
 %n个数据点，内节点为n-2个，U(k+1)作为初始值
for i=k+1:n+k-2
    U(i+1)=U(i)+temp(i-k)/sumtemp;  %n-1段弦长，U（k+1）到U（n+k）共n个节点
end

%第二步：反算n+2个控制点，采用切矢边界*****************************
%控制顶点的首末端点和给定型值点的首末端点重合
if(column == 2)
    dpt1=[0 1];%给定首数据点切矢
    dptn=[-1 0];%给定末数据点切矢
else
    dpt1=[0 0 1];%给定首数据点切矢
    dptn=[-1 0 0];%给定末数据点切矢
end
dU=zeros(1,n+k+3);  %节点增量，△U=Ui+1-Ui
for i=k+1:n+k-1
    dU(i)=U(i+1)-U(i);
end

%求解线性方程组获得控制顶点向量，A*D=E,A为系数矩阵，元素为B样条基函数的值；D是控制顶点列向量；E是列向量
A=zeros(n);
if(column == 2)
    E=zeros(n,2);  %2D curve
else
    E=zeros(n,3);  %3D curve
end
A(1,1)=1;  %切矢条件a1=1,b1=c1=0
A(n,n)=1;  %切矢条件an=bn=0,cn=1
E(1,:)=pt(1,:)+(dU(4)/3)*dpt1;    %首端点条件
E(n,:)=pt(n,:)-(dU(n+2)/3)*dptn;  %末端点条件
%计算系数矩阵A的元素a,b,c以及列向量E的元素e的值
for i=2:n-1  
    A(i,i-1)=dU(i+3).^2/(dU(i+1)+dU(i+2)+dU(i+3));                  %a的值
    A(i,i)=dU(i+3)*(dU(i+1)+dU(i+2))/(dU(i+1)+dU(i+2)+dU(i+3))+...  %b的值
        dU(i+2)*(dU(i+3)+dU(i+4))/(dU(i+2)+dU(i+3)+dU(i+4));
    A(i,i+1)=dU(i+2).^2/(dU(i+2)+dU(i+3)+dU(i+4));                  %c的值
    E(i,:)=(dU(i+2)+dU(i+3))*pt(i,:);                               %e的值
end
D=A\E;  %解方程组，获得去除首末端点的控制顶点向量
D=[pt(1,:);D;pt(n,:)]; %控制顶点向量，控制顶点比数据点多两个，加上首末端点
[s,t]=size(D);  %s：控制点数量

%第三步：求出Nurbs曲线*****************************
dt=0.001;  %Nurbs曲线插值密度，dt越小，越光滑
P=[];      %Nurbs曲线
syms dx;   %定义节点间的递增变量为变量dx
for i=k+1:s
    u=U;
    d=sym(D);
    %每次迭代后将d恢复初值
    for m=1:k
        for j=i-k:i-m
            alpha(j)=(dx-u(j+m))/(u(j+k+1)-u(j+m));
            d(j,:)=(1-alpha(j))*d(j,:)+alpha(j)*d(j+1,:);
        end
    end
    %带入变量dx求出各节点区间[ui,ui+1]内m=3的控制点d
    M=subs(d(i-k,:),dx,(u(i):dt:(u(i+1)-dt))');    %节点区间内的插值点替换dx
    P=[P;double(M)];
end
M=subs(d(s-k,:),dx,1);
P=[P;double(M)];%补上最后一个点

%第四步、绘制给定型值点、控制点及其Nurbs曲线*****************************
if(column == 2)  %2D curve
    plot(pt(:,1),pt(:,2),'*r');  %绘制型值点
    hold on;
    plot(D(:,1),D(:,2),'b-o');   %绘制控制点
    hold on;
    plot(P(:,1),P(:,2),'r');     %绘制nurbs曲线
    hold on;  
else             %3D curve
    plot3(pt(:,1),pt(:,2),pt(:,3),'*r'); %绘制型值点
    hold on;
    plot3(D(:,1),D(:,2),D(:,3),'b-o');   %绘制控制点
    hold on;
    plot3(P(:,1),P(:,2),P(:,3),'r');     %绘制nurbs曲线
    hold on;
end
axis equal;
grid on;

二维NURBS曲线插值：
在这里插入图片描述
三维NURBS曲线插值:

4 优化效率

上述程序在反求出控制点和节点矢量后，进行NURBS曲线生成时的效率较低，如果数据点过多，会十分影响运行速度，为此可以调用外部库函数提高NURBS曲线生成效率。这里使用的是 D.M. Spink提供的NURBS Toolbox，里面有大量的关于NURBS的函数。

%nurbs_toolbox库调用实例
%注释掉程序第三步的求解nurbs曲线的代码，运行以下代码，并修改第四步代码，使用nurbs_toolbox库函数绘制nurbs曲线函数
%使用nurbs_toolbox库函数nrbmak(coefs,knots)生成nurbs曲线，coefs是控制点向量；knots是节点向量
pts=[D(:,1)';D(:,2)';D(:,3)'];%nrbmak函数采用的行向量存储x,y,z坐标，上述程序是采用列向量存储x,y,z坐标，在此进行转换
crv=nrbmak(pts,U);
nrbplot(crv,1000);  %第二个参数是插值点的数量，值越大，则nurbs曲线越光滑

%第四步、绘制给定型值点、控制点及其Nurbs曲线*****************************
if(column == 2)  %2D curve
    plot(pt(:,1),pt(:,2),'*r');  %绘制型值点
    hold on;
    plot(D(:,1),D(:,2),'b-o');   %绘制控制点
    hold on;
%     plot(P(:,1),P(:,2),'r');     %绘制nurbs曲线
    nrbplot(crv,1000);    %使用nurbs_toolbox库函数绘制nurbs曲线
    hold on;  
else             %3D curve
    plot3(pt(:,1),pt(:,2),pt(:,3),'*r'); %绘制型值点
    hold on;
    plot3(D(:,1),D(:,2),D(:,3),'b-o');   %绘制控制点
    hold on;
%     plot3(P(:,1),P(:,2),P(:,3),'r');     %绘制nurbs曲线
    nrbplot(crv,1000);    %使用nurbs_toolbox库函数绘制nurbs曲线
    hold on;
end

参考文献

1、Piegl L, Tiller W. The NURBS book[M]. 2nd ed. Berlin Heidelberg: Springer-Verlag,1997
2、叶丽, 谢明红. 采用用积累弦长法拟合3次NURBS曲线[J]. 华侨大学学报(自然科学版), 2010, 31(4): 383-387
3、彭小军. 曲线曲面的NURBS 造型技术与数控仿真[D]. 西安: 长安大学, 2013