一、前言
本文主要实现了对stl文件的基本操作。首先我们识别了stl文件，读取其中的关键参数；然后我们对圆台进行了基本的操作，包括移动和旋转；最后我们对圆台进行切片，并读取切片的凸包和周长。
本文需要了解基础的三维知识，以及与numpy-stl相关的知识，可参考以下文献：
-
https://blog.csdn.net/ucsheep/article/details/121463494
-
https://blog.csdn.net/fuhanghang/article/details/84145833
-
-
https://pypi.org/project/numpy-stl/
-
https://numpy-stl.readthedocs.io/en/latest/usage.html
-
https://numpy-stl.readthedocs.io/en/latest/stl.html
-
https://www.codenong.com/7291b83ab659ca629bca/
-
https://w.wol.ph/2015/01/28/readingwriting-3d-stl-files-numpy-stl/
-
https://zhuanlan.zhihu.com/p/631144380
二、构建一个空心圆台
使用三维软件构建空心圆台（不用在乎尺寸，可随意绘制，要求底面必须与XY平面平行），圆台的上面和下面如图所示：
将模型保存为STL格式文件。文件可在本文绑定资源中下载。
三、安装numpy-stl，并读取文件
在conda环境中使用以下命令安装：
pip install numpy-stl使用自带的方法导入做好的空心圆台：
import numpy as np
from stl import mesh
mesh_ = mesh.Mesh.from_file('Cone.stl')
# print('所有点的x坐标：', mesh_.x)
# print('所有点的y坐标：', mesh_.y)
# print('所有点的z坐标：', mesh_.z)
# print('所有点的v0：', mesh_.v0)
# print('所有点的v1：', mesh_.v1)
# print('所有点的v2：', mesh_.v2)
# print('所有点的normals：', mesh_.normals)
print('所有点的points：', mesh_.points)
"""
所有点的points： [[ 4.5384525e+01 -1.2160616e+01 1.1811921e-04 ... 3.3223778e+01
-3.3223621e+01 1.1811921e-04]
[ 3.3223778e+01 -3.3223621e+01 1.1811921e-04 ... 1.2160756e+01
-4.5384342e+01 1.1811921e-04]
[ 4.5384525e+01 -1.2160616e+01 1.1811921e-04 ... 1.2160756e+01
-4.5384342e+01 1.1811921e-04]
...
[-1.1030294e+01 -2.9554515e+00 2.9312531e+01 ... -9.6316395e+00
-5.5607967e+00 2.9709482e+01]
[-1.0742725e+01 -2.8783977e+00 2.9709482e+01 ... -1.1419408e+01
1.7084386e-04 2.9312531e+01]
[-1.1419408e+01 1.7084386e-04 2.9312531e+01 ... -1.0742725e+01
-2.8783977e+00 2.9709482e+01]]
"""导入文件会被读取为mesh格式的数据，相关的参数在代码中已有体现，如果需要细致了解，参考源代码Base.py。其中，读取points会得到数组包含数组，基本的格式为：
[[x1, y1, z1, x2, y2, z2, x3, y3, z3,]...[]]其中最里面的数组包含9个数，分别是三角形网格的三个顶点。
四、渲染显示
为了方便在后台观察模型，而不是去打开别的软件，使用以下代码用matplotlib绘制模型：
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.tri as mtri
def show(mesh_):
# 绘制三维图形
fig = plt.figure(dpi=200)
ax = fig.add_subplot(111, projection='3d')
for mesh in mesh_:
vertices = np.array([mesh[0:3], mesh[3:6], mesh[6:9]]) # 这里定义了一个三角锥体
try:
ax.plot_trisurf(vertices[:, 0], vertices[:, 1], vertices[:, 2], linewidth=0.2, antialiased=True)
except:
print('vertices:', vertices)
# 设置坐标轴标签
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.title('3D Grid Triangles Visualization')
plt.show()
show(mesh_)我们主要是绘制了三角形在三维空间中的位置，结果如下所示：
五、移动旋转
使用库自带的方法函数就可以实现平移和旋转：
#平移
mesh_.translate(np.array([0, 30, 0]))#在y的方向平移
#旋转
mesh_.rotate([0.0, 1.0, 0.0], np.radians(180))#绕y轴旋转90度
show(mesh_)我们将模型翻转了180°，结果如下所示：
六、 切片
我们想要实现使用移动的XY平面对模型进行切割，首先考虑平面的标准法向量为（0, 0, 1)，才能在z方向自由移动。我们从模型的最低点开始逐渐往上进行切片，直到最顶端。切片实际上是集合操作，没有现成的方法。所以，需要记住以下几个推论：当一个平面与三角网格相交，至少有一个交点，最多有两个交点。在同一个位置，只且只有一个网格与平面相交。
由以上两个推论：我们就可以根据两个点和平面的z值，计算两点之间的直线与平面之间的交点（高中几何知识，不在此做推理，详细请看下文代码）。现展示切片后结果：
七、求凸包和周长
当我们切圆台的时候，会得到一大堆点，但是我们并不需要全部的点，我们只需要整个图像的外包线，也就是周长。可以使用scipy中的ConvexHull计算凸包，函数定义如下，供之后调用：
from scipy.spatial import ConvexHull
import numpy as np
#定义计算凸包和周长的函数
def get_hull_and_content(points):
# 计算凸包
hull = ConvexHull(points)
# 提取凸包的顶点坐标
convex_hull_points = points[hull.vertices]
# 将最后一个顶点复制到列表的开头以闭合凸包
convex_hull_points = np.vstack([convex_hull_points, convex_hull_points[0]])
# 计算凸包的周长
perimeter = np.sum(np.linalg.norm(convex_hull_points[1:] - convex_hull_points[:-1], axis=1))
# print("凸包的周长为:", perimeter)
return perimeter八、切片求周长代码、可视化周长变化
整体的代码如下所示，代码的逻辑在其中有标注。
delta = 0.1
min_z = np.min(mesh_.z)
max_z = np.max(mesh_.z)
print('min_z:', min_z, 'max_z:', max_z)#整个模型的最低点和最高点
Content = []#当前的index情况下，所有的凸包的周长
index = min_z
while True:
# print('start index:', index)
Cross = []#当前的index情况下，所有的交点
for mesh in mesh_:
#给出当前的最低和最高的z
# print('mesh:', mesh)
z_min = np.min([mesh[2], mesh[5], mesh[8]])
z_max = np.max([mesh[2], mesh[5], mesh[8]])
if z_min <= index and z_max >= index:
# print('z_min:', z_min, 'z_max:', z_max, 'index:', index)
#如果当前的三角形与当前的平面有交点, 直接计算交点
if index == mesh[2]:
#如果当前的平面与三角形的一个点重合
if not any(np.array_equal([mesh[0], mesh[1], mesh[2]], p) for p in Cross):
Cross.append([mesh[0], mesh[1], mesh[2]])
elif index == mesh[5]:
if not any(np.array_equal([mesh[3], mesh[4], mesh[5]], p) for p in Cross):
Cross.append([mesh[3], mesh[4], mesh[5]])
elif index == mesh[8]:
if not any(np.array_equal([mesh[6], mesh[7], mesh[8]], p) for p in Cross):
Cross.append([mesh[6], mesh[7], mesh[8]])
else:
#计算交点
#计算三角形的三个点,因为是网格，所以会有一半的重复点
#print('三角形的三个点')
p1 = np.array([mesh[0], mesh[1], mesh[2]])
p2 = np.array([mesh[3], mesh[4], mesh[5]])
p3 = np.array([mesh[6], mesh[7], mesh[8]])
Lower = []
Upper = []
for p in [p1, p2, p3]:
if p[2] < index:
Lower.append(p)
else:
Upper.append(p)
if len(Lower) == 2:
#如果有两个点在下面，一个点在上面,则交换位置，包证有两个点在上面，一个点在下面
R = Lower
Lower = Upper
Upper = R
for m in Upper:
x1 = m[0]
y1 = m[1]
z1 = m[2]
x3 = Lower[0][0]
y3 = Lower[0][1]
z3 = Lower[0][2]
S = (index - z3)/(z1 - z3)
x2 = x3 + S*(x1 - x3)
y2 = y3 + S*(y1 - y3)
if not any(np.array_equal([x2, y2, index], p) for p in Cross):
Cross.append([x2, y2, index])
#绘制散点图
X=[]
Y=[]
Points = []
for p in Cross:
X.append(p[0])
Y.append(p[1])
Points.append((p[0], p[1]))
# print('Points:', len(Points))
#plt.scatter(X, Y)
#plt.show()
# print(Points)
content = get_hull_and_content(np.array(Points))
Content.append(content)
# print('Cross Number:', len(Cross))
# print('end index:', index)
if index == max_z:
print('退出循环：',index, max_z)
break
else:
index = index+delta
if index > max_z:
index = max_z
print('Content:', Content)为了检验我们的切片是否正确，我们直接可视化所有的周长：
plt.figure(dpi=200)
# 生成对应的x坐标
x = list(range(1, len(Content) + 1))
# 绘制折线图
plt.plot(x, Content)
plt.xlabel('X Label') # 设置x轴标签
plt.ylabel('Y Label') # 设置y轴标签
plt.title('Line Plot') # 设置图标题
plt.show()结果如下所示，经过我们的计算，周长线性增长是必然，同时也说明了我们的方法是合理的。
九、后记
本文只是实现了一个简单的切片工作，如果需要进行更复杂的操作，需要读者根据自己的需求自定义。
627
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
