Matlab 绘制空间任意方向圆柱体
一 系统函数
Matlab本身自带了绘制圆柱体函数,如下:
[X,Y,Z] = cylinder
[X,Y,Z] = cylinder(r)
[X,Y,Z] = cylinder(r,n)
cylinder(axes_handle,...)
cylinder(...)
说明:
- cylinder 生成单位圆柱的 x、y 和 z 坐标。可使用 surf 或 mesh 绘制圆柱形对象,或不提供输出参数直接绘制它。
- [X,Y,Z] = cylinder 返回半径等于 1 的圆柱的 x、y 和 z 坐标。该圆柱绕其周长有 20 个等距点。
- [X,Y,Z] = cylinder® 使用 r 定义剖面曲线以返回圆柱的 x、y 和 z 坐标。cylinder 将 r 中的每个元素视为沿圆柱单位高度的等距高度的半径。该圆柱绕其周长有 20 个等距点。
- [X,Y,Z] = cylinder(r,n) 基于向量 r 定义的剖面曲线返回圆柱的 x、y 和 z 坐标。该圆柱绕其周长有 n 个等距点。
- cylinder(axes_handle,…) 将图形绘制到带有句柄 axes_handle 的坐标区中,而不是当前坐标区 (gca) 中。
- 不带任何输出参数的 cylinder(…) 使用 surf 绘制圆柱。
(备注:以上来源于Matlab帮助文档)
可见系统自带函数只能生成垂直于XY平面的圆柱体,要想绘制任意方向的圆柱体则需要在此基础上写程序实现。
二 由空间任意两点绘制空间圆柱体
方法一
- 思路:首先由系统自带函数
cylinder
生成标准圆柱体,该圆柱体轴线垂直于XY平面,底面原点位于坐标系原点,半径为r,高度为1。然后在此基础上通过旋转、拉伸变换得到想要的任意方向的圆柱体。 - 实现:
首先将cylinder
生成的圆柱体坐标的z坐标扩大L=norm(u1-u2)倍,然后将x,y,z坐标绕坐标y轴旋转angle(z),再将x,y,z坐标绕坐标z轴旋转angle(x),得到最终的圆柱体坐标,最后用surf
进行绘制。 - 代码
以下代码转至:https://blog.csdn.net/wushang_1314/article/details/91147550
function plotcylinder(u1,u2,color_a,r,alpha)
% 根据空间两点画圆柱
% u1,u2 ——空间两点
% color ——颜色
% r ——半决
% alpha ——透明度
L=norm(u1-u2);
ROD=u2-u1;
[X,Y,Z]=cylinder(r,100);
x1=X*0;y1=Y*0;z1=Z*0;
Z=L*Z-L/2;
ROD_midpoint=(u1+u2)/2;
x=ROD_midpoint(1);
y=ROD_midpoint(2);
z=ROD_midpoint(3);
a=[1 0 0];b=[0 1 0];c=[0 0 1];
if(ROD(2)==0||ROD(1)==0)
if(ROD(2)==0) % 在XZ平面
angel=acos(dot(ROD,c)/norm(ROD)/norm(c));
if(ROD(1)<0)
angel=-angel; %%
end
A2=[cos(angel) 0 sin(angel);0 1 0; -sin(angel) 0 cos(angel)]; % 绕Y轴旋转
for i=1:length(X(1,:))
u=[X(1,i) Y(1,i) Z(1,i)]';
u1=A2*u;
x1(1,i)=u1(1);y1(1,i)=u1(2);z1(1,i)=u1(3);
u=[X(1,i) Y(1,i) Z(2,i)]';
u1=A2*u;
x1(2,i)=u1(1);y1(2,i)=u1(2);z1(2,i)=u1(3);
end
end
if(ROD(1)==0) % 在YZ平面
angel=acos(dot(ROD,c)/norm(ROD)/norm(c));
angel=2*pi-angel; % 此处maybe要增加类似上面的判定是否变为负数
A1=[1 0 0;0 cos(angel) -sin(angel);0 sin(angel) cos(angel)]; % 绕X轴旋转
for i=1:length(X(1,:))
u=[X(1,i) Y(1,i) Z(1,i)]';
u1=A1*u;
x1(1,i)=u1(1);y1(1,i)=u1(2);z1(1,i)=u1(3);
u=[X(1,i) Y(1,i) Z(2,i)]';
u1=A1*u;
x1(2,i)=u1(1);y1(2,i)=u1(2);z1(2,i)=u1(3);
end
end
else
% 先绕Z轴旋转,在绕Y轴旋转
angel=acos(dot(ROD,c)/norm(ROD)/norm(c));
A2=[cos(angel) 0 sin(angel);0 1 0; -sin(angel) 0 cos(angel)];
angel=acos(dot(ROD,a)/norm(ROD)/norm(a));
if(ROD(2)<0)
angel=2*pi-angel;
end
A3=[cos(angel) -sin(angel) 0;sin(angel) cos(angel) 0;0 0 1];
A=A3*A2;
for i=1:length(X(1,:))
u=[X(1,i) Y(1,i) Z(1,i)]';
u1=A*u;
x1(1,i)=u1(1);y1(1,i)=u1(2);z1(1,i)=u1(3);
u=[X(1,i) Y(1,i) Z(2,i)]';
u1=A*u;
x1(2,i)=u1(1);y1(2,i)=u1(2);z1(2,i)=u1(3);
end
end
fill3(x1(1,:)+x,y1(1,:)+y,z1(1,:)+z,color_a,'EdgeColor','none','FaceAlpha',alpha) %底面
hold on
fill3(x1(2,:)+x,y1(2,:)+y,z1(2,:)+z,color_a,'EdgeColor','none','FaceAlpha',alpha) %顶面
hold on
surf(x1+x,y1+y,z1+z,'facecolor',color_a,'edgecolor','none','FaceAlpha',alpha) %圆柱表面
% axis equal
% xlabel('x')
% ylabel('y')
% zlabel('z')
set(gcf,'color','w') %gcf是当前图窗的句柄,用图窗句柄可查询和修改图窗属性
- 缺陷
但是笔者在使用过程中发现某些情况会出错,如下:
方法二
通过观察[X,Y,Z] = cylinder(r,n)
生成的圆柱体数据类型时发现,X,Y,Z均为一个2行n列的矩阵,用于存储圆柱底面和顶面两个圆的空间数据点,然后通过surf(X,Y,Z,'facecolor',color,'edgecolor','none','FaceAlpha',alpha)
进行绘制。
即只需求出空间圆柱底面与顶面圆的方程即可绘制圆柱体,具体方法如下:
- 思路:
已知空间两点,即可分别求出两个圆的方程,通过空间一点与法向即可求出空间圆的方程,方法可由以下博客得出: http://blog.sina.com.cn/s/blog_6496e38e0102vi7e.html
- 实现:
function myplotcylinder(u1,u2,color,r,alpha)
% 根据空间两点画圆柱
% u1,u2 ——空间两点
% color ——颜色
% r ——半径
% alpha ——透明度
n=u2-u1;
theta=(0:2*pi/100:2*pi)'; %theta角从0到2*pi
a=cross(n,[1 0 0]); %n与i叉乘,求取a向量
if ~any(a) %如果a为零向量,将n与j叉乘 %any检测矩阵中是否有非零元素,如果有,则返回1,否则,返回0
a=cross(n,[0 1 0]);
end
b=cross(n,a); %求取b向量
a=a/norm(a); %单位化a向量
b=b/norm(b); %单位化b向量
%a,b是满足既垂直于n,又互相垂直的任意单位向量
%底面圆方程
c1=u1(1)*ones(size(theta,1),1);
c2=u1(2)*ones(size(theta,1),1);
c3=u1(3)*ones(size(theta,1),1);
x1=c1+r*a(1)*cos(theta)+r*b(1)*sin(theta); %圆上各点的x坐标
y1=c2+r*a(2)*cos(theta)+r*b(2)*sin(theta); %圆上各点的y坐标
z1=c3+r*a(3)*cos(theta)+r*b(3)*sin(theta); %圆上各点的z坐标
%顶面圆方程
c1=u2(1)*ones(size(theta,1),1);
c2=u2(2)*ones(size(theta,1),1);
c3=u2(3)*ones(size(theta,1),1);
x2=c1+r*a(1)*cos(theta)+r*b(1)*sin(theta); %圆上各点的x坐标
y2=c2+r*a(2)*cos(theta)+r*b(2)*sin(theta); %圆上各点的y坐标
z2=c3+r*a(3)*cos(theta)+r*b(3)*sin(theta); %圆上各点的z坐标
X(1,:)=x1(:);
X(2,:)=x2(:);
Y(1,:)=y1(:);
Y(2,:)=y2(:);
Z(1,:)=z1(:);
Z(2,:)=z2(:);
fill3(X(1,:),Y(1,:),Z(1,:),color,'EdgeColor','none','FaceAlpha',alpha) %底面
hold on
fill3(X(2,:),Y(2,:),Z(2,:),color,'EdgeColor','none','FaceAlpha',alpha) %顶面
hold on
surf(X,Y,Z,'facecolor',color,'edgecolor','none','FaceAlpha',alpha) %圆柱表面
<完>
后记:笔者才疏学浅,如有错误,望指出。