三维闭合B样条曲线拟合算法Matlab代码

在做相关项目需要解决B样条插值问题，记录如下

感谢以下参考资料的帮助

% ref: 闭合 B 样条曲线控制点的快速求解算法及应用 % http://www.doc88.com/p-5714423317458.html % https://blog.csdn.net/liumangmao1314/article/details/54588155

=============================

%计算B样条曲线控制点

function px = LU_B1(CPnum, V)        a = 1;     b = 4;     c = 1;     d = 1;     e = 1;          f = zeros(CPnum-1,1);     g = zeros(CPnum-2,1);     h = zeros(CPnum,1);     k = zeros(CPnum-1,1);

    % get h[] & f[]     h(1) = b;     for i=1:CPnum-2         f(i) = a/h(i);         h(i+1) = b - f(i)*c;     end

    % get g[] & f[n-1]     g(1) = d/h(1);     for i=1:CPnum-3         g(i+1) = -g(i)*c/h(i+1);     end     f(CPnum-1) = ( a-g(CPnum-2)*c )/h(CPnum-1);

    % get k[] & h[n]     k(1) = e;     for i=1:CPnum-3         k(i+1) = -f(i)*k(i);     end     k(CPnum-1) = c - f(CPnum-2)*k(CPnum-2);

    gksum = 0;     for i=1:CPnum-2         gksum = gksum + g(i)*k(i);     end     h(CPnum) = b - gksum - f(CPnum-1)*c;

    % 矩阵求解过程，追的过程     x = zeros(CPnum,1);        x(1) = 6*V(end);          for i=1:CPnum-2               x(i+1) = 6*V(i) - f(i)*x(i);            end

    gxsum = 0;    

    for i=1:CPnum-2         gxsum = gxsum + g(i)*x(i);             end          x(CPnum) = 6*V(CPnum-1) - gxsum - f(CPnum-1)*x(CPnum-1);              % 赶的过程     px = zeros(CPnum+2,1);              px(CPnum) = x(CPnum)/h(CPnum);     px(CPnum-1) = ( x(CPnum-1)-k(CPnum-1)*px(CPnum) )/h(CPnum-1);          for i=CPnum-2:-1:1         px(i) = ( x(i)-c*px(i+1)-k(i)*px(CPnum) )/h(i);     end     px(CPnum+1) = px(1);     px(CPnum+2) = px(2);     end

---

% 插值计算三次周期性b样条曲线 function p = Cubic_Bsp(u,x)     b(1) = ((1-u)^3)/6;     b(2) = (3*u^3-6*u^2+4)/6;     b(3) = (-3*u^3+3*u^2+3*u+1)/6;     b(4) = (u^3)/6;

    p = x*b'; end

---

%test.m 测试

clear;clc x = [-1;-1;1;1;0.2]; y = [1;-0.5;-1;1;0.8]; z = zeros(6,1); z(3) = 1;

CPnum = size(x,1);

%x,y,z坐标分别计算B样条曲线控制点

px = LU_B1(CPnum, x); py = LU_B1(CPnum, y); pz = LU_B1(CPnum, z);

%首尾相连，曲线闭合

x(CPnum+1) = x(1); y(CPnum+1) = y(1); z(CPnum+1) = z(1);

figure; plot3(x,y,z,'ro-'); hold on; for i=1:CPnum%用这个循环     c=num2str(i);     c=[' ',c];     text(x(i),y(i),z(i),c)     %text(x(i),y(i),c) end

plot3(px,py,pz,'g-'); for i=1:CPnum+1%用这个循环     c=num2str(i);     c=[' ',c];     text(px(i),py(i),pz(i),c)     plot3(px(i),py(i),pz(i),'*'); end

axis equal

px(CPnum+3) = px(3); py(CPnum+3) = py(3); pz(CPnum+3) = pz(3);

px(CPnum+4) = px(4); py(CPnum+4) = py(4); pz(CPnum+4) = pz(4);

%两点之间对10个点进行插值计算

nP = 10; delta = 1.0/nP; for j=1:CPnum     u = 0.0;     for i=1:nP                  p = px';         Xi = p(j:j+3);         xx(i+(j-1)*nP) = Cubic_Bsp(u,Xi);           p = py';         Yi = p(j:j+3);         yy(i+(j-1)*nP) = Cubic_Bsp(u,Yi);                          p = pz';         Zi = p(j:j+3);         zz(i+(j-1)*nP) = Cubic_Bsp(u,Zi);                          u = u + delta;     end end xx(end+1) = x(1); yy(end+1) = y(1); zz(end+1) = z(1); plot3(xx,yy,zz,'b--')

 

结果如图所示，红色为三维空间5个不共面点

绿色为B样条曲线控制点

蓝色为B样条拟合曲线