一、理论出处
本次实现三次Bezier可展曲面的理论出处为：
作者： Chih-Hsing Chu, Carlo H. Sequin
文章名： Developable Bezier patches: properties and design
出处： Computer-Aided Design 34 (2002) 511—527
二、可展的相关理论及算法
具体图示：
具体理论：
四、 matlab程序实现
function Developable_Cubic_surface
% date : 2020.10.16 ---- 202.10.17
% 第三篇论文的编程实现
% Developable Bezier patches: properties and design
% 双三次的可展曲面设计
% author ： mw_1422102031
clear;
clc;
%A曲线的控制顶点
A0 = [0,0,0.1]; A1 = [1,2,2]; A2 = [2,1,1]; A3 = [3,0,-3];
ax = [A0(1),A1(1),A2(1),A3(1)];
ay = [A0(2),A1(2),A2(2),A3(2)];
az = [A0(3),A1(3),A2(3),A3(3)];
[U,V] = meshgrid(0:0.05:1,0:0.05:1);
%给定c0，c3的方向向量
c0 = [0,1,1];
c3 = [0,3,-1];
%由控制顶点计算a1,a2,a3
a1 = []; a2 = []; a3 = [];
a1(1) = ax(2) - ax(1); a1(2) = ay(2) - ay(1); a1(3) = az(2) - az(1);
a2(1) = ax(3) - ax(2); a2(2) = ay(3) - ay(2); a2(3) = az(3) - az(2);
a3(1) = ax(4) - ax(3); a3(2) = ay(4) - ay(3); a3(3) = az(4) - az(3);
%计算求解过程中的中间量
d11 = -det([a2;c0;a1])*det([a3;c0;c3]);
d = det([a2;c0;a1])*det([a2;a3;c3]) - det([a1;c0;a3])*det([a3;a1;c3]);
d1 = d11 / d;
d21 = det([a1;c0;c3])*det([a3;a1;c3]);
d2 = d21 / d;
d31 = det([a1;c0;a3])*det([a3;c0;c3]);
d3 = d31 / d;
d41 = -det([a1;c0;c3])*det([a2;a3;c3]);
d4 = d41 / d;
% 计算非线性方程组的系数
f1 = det([a1;c0;a3]);
f2 = 3*det([a1;c0;c3]) + 3*d2*det([a1;c0;a3]);
f3 = 3*d1* det([a1;c0;a3]);
f4 = 4* det([a2;c0;c3]) + 4*d2* det([a2;c0;a3]) + d4* det([a3;c0;a1]);
f5 = 4*d1* det([a2;c0;a3]) + d3* det([a3;c0;a1]);
g1 = d1* det([a1;a3;c3]) + 4*d3* det([a2;a1;c3]) + 4* det([a2;c0;c3]);
g2 = d2* det([a1;a3;c3]) + 4*d4* det([a2;a1;c3]);
g3 = det([a3;c0;c3]);
g4 = 3* det([a3;c0;c3]) + 3*d3* det([a3;a1;c3]);
g5 = 3*d4* det([a3;a1;c3]);
h1 = det([a1;c0;c3]);
h2 = det([a2;c0;c3]);
h3 = 3*det([a2;c0;c3]) + 3*d3* det([a2;a1;c3]) + 3*d2*det([a2;c0;a3]) + ...
3*d1*d4*det([a2;a1;a3]) + 3*d2*d3*det([a2;a1;a3]);
h4 = 3*d4*det([a2;a1;c3]) + 3*d2*d4*det([a2;a1;a3]);
h5 = 3*d1*det([a2;a1;a3]) + 3*d1*d3*det([a2;a1;a3]);
h6 = det([a3;c0;c3]);
k = 1;
syms x y z
f(x,y,z) = f1*x*y + f2*k*z + f3*z^2 + f4*k*x + f5*x*z;
g(x,y,z) = g1*y*z + g2*k*y + g3*x*y + g4*k*z + g5*k^2;
h(x,y,z) = h1*y*z + h2*x*y + h3*k*z + h4*k^2 + h5*z^2 + h6*k*x;
eqns = [f(x,y,z)==0,g(x,y,z)==0,h(x,y,z)==0];
[x,y,z] = vpasolve(eqns,x,y,z); %会出现多个解
% 得到的的非线性方程组的解不唯一，成为了列向量的形式
A = [x,y,z];
%求解非线性方程组得到的向量A(:,i) = (α,β,μ),i = 1,2,3,...； α,β,μ对应论文中的未知量
%确定c0，c3 的具体向量值，求解c1，c2
%β = beta , α = alpha ， μ = mu ， η = eta = k, ρ = rho ，δ = delta
x = A(:,2);
%因为求出来的解是虚数，我将虚数取模了
x = [norm(x(1),Inf),norm(x(2),Inf),norm(x(3),Inf)];
c0_bar = x(1)*c0;
c3_bar = x(2)*c3;
rho = d1*x(3) + d2*k;
delta = d3*x(3) + d4*k;
c1 = x(3)*c0 + delta*a1;
c2 = k*c3 + rho*a3;
%计算B曲线的的控制顶点
B0 = A0 + c0_bar; B1 = A1 + c1; B2 = A2 + c2; B3 = A3 + c3_bar;
bx = [B0(1),B1(1),B2(1),B3(1)];
by = [B0(2),B1(2),B2(2),B3(2)];
bz = [B0(3),B1(3),B2(3),B3(3)];
%A_bezier曲线
x = V;
n=length(ax)-1;
Ax = 0;
Ay = 0;
Az = 0;
for i = 1:n+1
Ax = Ax + ax(i)*B(x,n,i-1) ; %将控制顶点与Bernstein基函数相乘得到bezier曲线
Ay = Ay + ay(i)*B(x,n,i-1) ;
Az = Az + az(i)*B(x,n,i-1) ;
end
figure(1)
plot3(ax,ay,az,'k',ax,ay,az,'m*') %控制多边形
hold on
plot3(Ax,Ay,Az,'r') %bezier曲线A(v)
%B_bezier曲线
x = V;
n=length(bx)-1;
Bx = 0; By = 0; Bz = 0;
for i = 1:n+1
Bx = Bx + bx(i)*B(x,n,i-1) ; %将控制顶点与Bernstein基函数相乘得到bezier曲线
By = By + by(i)*B(x,n,i-1) ;
Bz = Bz + bz(i)*B(x,n,i-1) ;
end
figure(2)
plot3(bx,by,bz,'k',bx,by,bz,'m*') %控制多边形
hold on
plot3(Bx,By,Bz,'r')
%将两条边界bezier曲线带入直纹面方程
X = (1-U).*Ax + U.*Bx ;
Y = (1-U).*Ay + U.*By ;
Z = (1-U).*Az + U.*Bz ;
%两种方式画图
figure(3)
mesh(X,Y,Z)
end
结果显示：
具体的理论在我之前的博客里面有：https://blog.csdn.net/mw_1422102031/article/details/108932741
超链接
五、 matlab程序实现（中间编程出错误的版本）
function X_Developable_Cubic_surface
%用一般牛顿法求解非线性方程出现错误，可能涉及到精度问题，具体原因不找了
% date : 2020.10.16 ---- 202.10.17
% 第三篇论文的编程实现
% Developable Bezier patches: properties and design
% 双三次的可展曲面设计
% author ： mw_1422102031
clear;
clc;
%A曲线的控制顶点
A0 = [0,0,0]; A1 = [1,2,2]; A2 = [2,1,1]; A3 = [3,0,-1];
ax = [A0(1),A1(1),A2(1),A3(1)];
ay = [A0(2),A1(2),A2(2),A3(2)];
az = [A0(3),A1(3),A2(3),A3(3)];
[U,V] = meshgrid(0:0.02:1,0:0.02:1);
%给定c0，c3的方向向量
c0 = 0.01*[0,1,1];
c3 = 0.01*[0,3,-1];
%由控制顶点计算a1,a2,a3
a1 = []; a2 = []; a3 = [];
a1(1) = ax(2) - ax(1); a1(2) = ay(2) - ay(1); a1(3) = az(2) - az(1);
a2(1) = ax(3) - ax(2); a2(2) = ay(3) - ay(2); a2(3) = az(3) - az(2);
a3(1) = ax(4) - ax(3); a3(2) = ay(4) - ay(3); a3(3) = az(4) - az(3);
%计算求解过程中的中间量
d11 = -det([a2;c0;a1])*det([a3;c0;c3]);
d = det([a2;c0;a1])*det([a2;a3;c3]) - det([a1;c0;a3])*det([a3;a1;c3]);
d1 = d11 / d;
d21 = det([a1;c0;c3])*det([a3;a1;c3]);
d2 = d21 / d;
d31 = det([a1;c0;a3])*det([a3;c0;c3]);
d3 = d31 / d;
d41 = -det([a1;c0;c3])*det([a2;a3;c3]);
d4 = d41 / d;
% 计算非线性方程组的系数
f1 = det([a1;c0;a3]);
f2 = 3*det([a1;c0;c3]) + 3*d2*det([a1;c0;a3]);
f3 = 3*d1* det([a1;c0;a3]);
f4 = 4* det([a2;c0;c3]) + 4*d2* det([a2;c0;a3]) + d4* det([a3;c0;a1]);
f5 = 4*d1* det([a2;c0;a3]) + d3* det([a3;c0;a1]);
g1 = d1* det([a1;a3;c3]) + 4*d3* det([a2;a1;c3]) + 4* det([a2;c0;c3]);
g2 = d2* det([a1;a3;c3]) + 4*d4* det([a2;a1;c3]);
g3 = det([a3;c0;c3]);
g4 = 3* det([a3;c0;c3]) + 3*d3* det([a3;a1;c3]);
g5 = 3*d4* det([a3;a1;c3]);
h1 = det([a1;c0;c3]);
h2 = det([a2;c0;c3]);
h3 = 3*det([a2;c0;c3]) + 3*d3* det([a2;a1;c3]) + 3*d2*det([a2;c0;a3]) + ...
3*d1*d4*det([a2;a1;a3]) + 3*d2*d3*det([a2;a1;a3]);
h4 = 3*d4*det([a2;a1;c3]) + 3*d2*d4*det([a2;a1;a3]);
h5 = 3*d1*det([a2;a1;a3]) + 3*d1*d3*det([a2;a1;a3]);
h6 = det([a3;c0;c3]);
%给出非线性方程组
k = 1.5;
syms x y z
f(x,y,z) = f1*x*y + f2*k*z + f3*z^2 + f4*k*x + f5*x*z;
g(x,y,z) = g1*y*z + g2*k*y + g3*x*y + g4*k*z + g5*k^2;
h(x,y,z) = h1*y*z + h2*x*y + h3*k*z + h4*k^2 + h5*z^2 + h6*k*x;
%求这三个方程的偏导数
f_1(x,y,z) = diff(f(x,y,z),x); f_2(x,y,z) = diff(f(x,y,z),y); f_3(x,y,z) = diff(f(x,y,z),z);
g_1(x,y,z) = diff(g(x,y,z),x); g_2(x,y,z) = diff(g(x,y,z),y); g_3(x,y,z) = diff(g(x,y,z),z);
h_1(x,y,z) = diff(h(x,y,z),x); h_2(x,y,z) = diff(h(x,y,z),2); h_3(x,y,z) = diff(h(x,y,z),z);
%%一般牛顿法求解非线性方程组
x0 = [0.1, 0.1, -0.1]';
eps1 = 1e-8 ;
eps2 = 1e-8 ;
KMAX = 20;
x = x0;
k = 0;
for k = 0:KMAX
%求解F(x)
m = zeros(size(x));
m(1) = f(x(1),x(2),x(3)) ;
m(2) = g(x(1),x(2),x(3)) ;
m(3) = h(x(1),x(2),x(3)) ;
b = m;
norm_b = norm(b,Inf);
if norm_b<eps1
break
end
%求解Jac(x)
A = [];
A = [f_1(x(1),x(2),x(3)), f_2(x(1),x(2),x(3)), f_3(x(1),x(2),x(3));...
g_1(x(1),x(2),x(3)), g_2(x(1),x(2),x(3)), g_3(x(1),x(2),x(3));...
h_1(x(1),x(2),x(3)), h_2(x(1),x(2),x(3)), h_3(x(1),x(2),x(3))];
dx = -A\b;
x = x+dx;
norm_dx = norm(dx,Inf);
if norm_dx < eps2
break
end
end
x
%求解非线性方程组得到的向量x = (α,β,μ)； α,β,μ对应论文中的参数
%确定c0，c3 的具体向量值，求解c1，c2
%β = beta , α = alpha ， μ = mu ， η = eta = k, ρ = rho ，δ = delta
c0_bar = x(1)*c0;
c3_bar = x(2)*c3;
rho = d1*x(3) + d2*k;
delta = d3*x(3) + d4*k;
c1 = x(3)*c0 + delta*a1;
c2 = k*c3 + rho*a3;
%计算B曲线的的控制顶点
B0 = A0 + c0_bar; B1 = A1 + c1; B2 = A2 + c2; B3 = A3 + c3_bar;
bx = [B0(1),B1(1),B2(1),B3(1)];
by = [B0(2),B1(2),B2(2),B3(2)];
bz = [B0(3),B1(3),B2(3),B3(3)];
%A_bezier曲线
x = V;
n=length(ax)-1;
Ax = 0;
Ay = 0;
Az = 0;
for i = 1:n+1
Ax = Ax + ax(i)*B(x,n,i-1) ; %将控制顶点与Bernstein基函数相乘得到bezier曲线
Ay = Ay + ay(i)*B(x,n,i-1) ;
Az = Az + az(i)*B(x,n,i-1) ;
end
figure(1)
plot3(ax,ay,az,'k',ax,ay,az,'m*') %控制多边形
hold on
plot3(Ax,Ay,Az,'r') %bezier曲线A(v)
%B_bezier曲线
x = V;
n=length(bx)-1;
Bx = 0; By = 0; Bz = 0;
for i = 1:n+1
Bx = Bx + bx(i)*B(x,n,i-1) ; %将控制顶点与Bernstein基函数相乘得到bezier曲线
By = By + by(i)*B(x,n,i-1) ;
Bz = Bz + bz(i)*B(x,n,i-1) ;
end
figure(2)
plot3(bx,by,bz,'k',bx,by,bz,'m*') %控制多边形
hold on
plot3(Bx,By,Bz,'r')
%将两条边界bezier曲线带入直纹面方程
X = (1-U).*Ax + U.*Bx ;
Y = (1-U).*Ay + U.*By ;
Z = (1-U).*Az + U.*Bz ;
%两种方式画图
figure(3)
mesh(X,Y,Z)
end
结果就不展示了，本身这个编程就是失败的。
这个我用一般牛顿法去求解非线性方程组，中间迭代的时候一直出错，并且得到得结果是错的。然后用matlab中自带的求解非线性方程组的方法把这个问题解决的。
2855
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
