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

1 篇文章 0 订阅

% 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--')

• 8
点赞
• 7
评论
• 50
收藏
• 一键三连
• 扫一扫，分享海报

11-02

09-09 4257
03-16 1万+
04-23
09-05
08-19
12-20
11-14 1672
03-28