一、前言
本文是在我前一篇博文的基础上进一步扩展，除了对实体进行切片，求周长，还新增了旋转实体的操作。
（源代码）基于numpy-stl操作stl文件——读取圆台z轴截面的周长
1. 模型准备
创建一个不规则的实体模型，并保存为stl：
模型与xy平面平行，由两个回转的药丸形状构成。
2. 模型渲染
为了方便可视化模型，我们使用以下代码在notebook里面渲染：
from stl import mesh
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.tri as mtri
import time
mesh_ = mesh.Mesh.from_file('test.stl')
def show(mesh_):
# 绘制三维图形
time1 = time.time()
fig = plt.figure(dpi=300)
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()
time2 = time.time()
print('caculate time:', time2 - time1)
show(mesh_)结果如下所示：
二、分析
我们想要对这个平躺着的模型进行切片处理，最开始方便办法就是将模型立起来。但是模型的轴向向量并不可知，如何计算得到模型的轴向向量是本文难题。得到模型轴向向量之后，我们可以使用轴向向量进行坐标变换，将模型的轴与z轴对齐。由此便可以直接套用我之前一篇的博客的工作，得到切片的周长。
1. 如何获取轴向向量
首先我们观察模型可知，模型与xy平面平行。那么，当我们使用一个垂直于xy平面且过原点的平面y=kx于模型相交时，会得到一个椭圆的面。 当且仅当，这个椭圆的周长最小时，这个切平面的法向量就是模型的轴向向量。
2. 代码实现获取轴向向量
在上一节中我们已经说明了几何原理，这一节我们使用代码实现这个想法。
2.1 导入相关库
from scipy.spatial import ConvexHull
import numpy as np
from stl import mesh
import matplotlib.pyplot as plt
import time2.2 计算空间三角形与和y=kx平面的交点
def check_intersection_with_plane(A, B, C, k):
def line_plane_intersection(P1, P2, k):
x1, y1, z1 = P1
x2, y2, z2 = P2
# 计算 t 的分母和分子
denominator = k * (x2 - x1) - (y2 - y1)
if denominator == 0:
# 如果分母为0，说明直线平行于平面，没有交点或是无穷多个交点
return None
t = (y1 - k * x1) / denominator
# 判断 t 是否在 [0, 1] 范围内，确定交点是否在线段上
if 0 <= t <= 1:
# 计算交点的坐标
x = x1 + t * (x2 - x1)
y = y1 + t * (y2 - y1)
z = z1 + t * (z2 - z1)
return (x, y, z)
else:
return None
# 检查三角形的每一条边是否与平面 y = kx 相交
edges = [(A, B), (B, C), (C, A)]
intersection_points = []
for P1, P2 in edges:
intersection = line_plane_intersection(P1, P2, k)
if intersection:
intersection_points.append(intersection)
if len(intersection_points) > 2:
# print(len(intersection_points))
intersection_points=[]
return intersection_points if intersection_points else None以上代码计算了，如果空间中的三个点构造的三角形与平面y=kx相交，则返回交点。由我前一篇博文的推论可知，交点只会有一点或者两点。具体实现原理是高中几何知识，不在此赘述。
2.3 计算空间中的点的凸包及周长
#计算凸包
def compute_convex_hull(points):
hull = ConvexHull(points)
return hull
def calculate_hull_perimeter(hull):
perimeter = 0
for simplex in hull.simplices:
# 每个 simplex 包含凸包中的两个点的索引
point1 = hull.points[simplex[0]]
point2 = hull.points[simplex[1]]
# 计算两个点之间的欧几里得距离
edge_length = np.linalg.norm(point1 - point2)
perimeter += edge_length
return perimeter以上两个函数计算了一个交点集合的凸包，求出的周长是某一k值下切到的。
2.4 计算单次切片的周长
#切片计算周长
def slides(mesh_,k):
Cross = []
# 遍历每个三角形
for triangle in mesh_.points:
A = np.array([triangle[0], triangle[1], triangle[2]])
B = np.array([triangle[3], triangle[4], triangle[5]])
C = np.array([triangle[6], triangle[7], triangle[8]])
has_intersection = check_intersection_with_plane(A, B, C, k)
if has_intersection:
for h in has_intersection:
if not any(np.array_equal(h, p) for p in Cross):
Cross.append(h)
# print(f"Number of intersection points: {len(Cross)}")
# 创建一个三维图形对象
plt.figure()
M=[]
Z=[]
# 绘制散点图
if Cross:
X, Y, Z = zip(*Cross)
x = np.sqrt(np.array(X)**2 + np.array(Y)**2)
# 根据 X 的正负确定投影后的正负
M = np.sign(X) * x
plt.scatter(M, Z, c='r', s=1)
plt.show()
result = []
for i in range(len(M)):
result.append((M[i],Z[i]))
# return result
# 计算三维凸包
hull = compute_convex_hull(result)
# 计算凸包的周长
perimeter = calculate_hull_perimeter(hull)
return perimeter这个函数输入的mesh：也就是我们模型的网格坐标，和k斜率：平面位置。这个函数能够返回在当前k值下，与模型相交的点的凸包的周长，如下所示：
2.5 计算轴向向量
def uniform_split(a, b, k):
return np.linspace(a, b, k)
def get_normal(mesh_):
# 示例用法
a = -1
b = 1
Num = 10
Min_k = 0
while True:
X=uniform_split(a, b, Num)
Content = []
for k in X:
# if k != 0:
# print(f"y = {k}x")
Content.append(slides(mesh_,k))
plt.figure(dpi=200)
# 绘制折线图
plt.plot(X, Content)
plt.xlabel('k') # 设置x轴标签
plt.ylabel('content') # 设置y轴标签
plt.title('Line Plot') # 设置图标题
plt.show()
print('Stop:',abs(Min_k-X[np.argmin(Content)]) < 0.0001)
if abs(Min_k-X[np.argmin(Content)]) < 0.0001 :
break
else:
print('Min K:',Min_k,'————》 Min X；',X[np.argmin(Content)],'delta:',abs(X[np.argmin(Content)]-Min_k),'Num:',Num)
delta= X[np.argmin(Content)]-Min_k
Min_k = X[np.argmin(Content)]
a = Min_k - 3*delta
b = Min_k + 3*delta
Num +=10
print('Now a:',a,'b:',b,'k:',Num)
return [-Min_k,1,0] 这段代码实现了我们如何去搜索最优的平面，达到切面的周长最小。我们首先设定k在-1和1之间（依据经验设定，不同的模型应该不同），分为10段，初始k为0。在每一次计算得到切面周长的集合之后，我们更新k的值，并将上下限限制在k的正负三倍delta之间，同时分段数加10。更新上下限和增加分段数，能够加速算法的收索，快速达到最优解。当delta小于0.0001时，结束搜索，并返回平面的法向量，这个法向量便是模型的轴向量。从下图可以看出，最终的搜索结果为：k = 0.72650124左右。（旋转完成是因为，函数都被集成了，所以直接运行下去，可以不看）
3. 对模型进行旋转
import numpy as np
from stl import mesh
def rotation_(my_mesh):
# 定义计算旋转矩阵的函数
def rotation_matrix_from_vectors(vec1, vec2):
""" Find the rotation matrix that aligns vec1 to vec2
:param vec1: A 3d "source" vector
:param vec2: A 3d "destination" vector
:return mat: A transform matrix (3x3)
"""
a, b = (vec1 / np.linalg.norm(vec1)).reshape(3), (vec2 / np.linalg.norm(vec2)).reshape(3)
v = np.cross(a, b)
c = np.dot(a, b)
s = np.linalg.norm(v)
kmat = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
return np.eye(3) + kmat + kmat @ kmat * ((1 - c) / (s ** 2))
normal = get_normal(my_mesh)
# 定义原始方向向量和目标方向向量（单位向量）
initial_vector = np.array(normal) #
target_vector = np.array([0,0,1]) #z轴方向
target_unit_vector = target_vector / np.linalg.norm(target_vector)
# 计算旋转矩阵
R = rotation_matrix_from_vectors(initial_vector, target_unit_vector)
# 旋转所有顶点
my_mesh.vectors = np.dot(my_mesh.vectors.reshape(-1, 3), R).reshape(-1, 3, 3)
# 保存旋转后的 STL 文件
my_mesh.save('my_rotated_model_surface.stl')
print('旋转完成')
return my_mesh以上这个函数会根据给出的mesh网格，自动计算轴向量，然后使用库方法惊醒旋转，保存为my_rotated_model_surface.stl文件。结果如下所示：
三、切片求周长
我们已经得到了想要的标准模型，接下来使用我上一篇博文的方法，进行切分：
from scipy.spatial import ConvexHull
import numpy as np
import matplotlib.pyplot as plt
#定义计算凸包和周长的函数
def get_hull_and_content(points):
# 计算凸包
try:
hull = ConvexHull(points)
except:
print('points:', points)
print('凸包计算失败')
return 0
# 提取凸包的顶点坐标
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
# 提取切片周长
def slide_mesh(mesh_,delta):
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
time1 = time.time()
x_for_plot = []
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.axis('equal')
# plt.show()
# print(Points)
if len(Points) > 2:
#至少三个点，形成一个平面
content = get_hull_and_content(np.array(Points))
Content.append(content)
x_for_plot.append(index)
# print('Cross Number:', len(Cross))
# print('end index:', index)
if index == max_z:
print('退出循环：',index, max_z)
break
else:
time2 = time.time()
print('当前index:', index, '已经运行时间:', time2 - time1)
index = index+delta
if index > max_z:
index = max_z
print('Content:', Content)
#delta = 10 , 60s
#delta = 5, 123s
return x_for_plot, Content四、结果
接下来的代码将上文的函数组合使用：
# 读取 STL 文件
my_mesh = mesh.Mesh.from_file('test.stl')
my_mesh = rotation_(my_mesh)#旋转模型
delta = 5 #切片的间隔,高度差
x_for_plot, Content = slide_mesh(my_mesh,delta=delta)#计算切片周长
#旋转后的文件名字为：my_rotated_model_surface.stl
# 绘制周长的变化图
plt.figure(dpi=200)
x_max = np.max(x_for_plot)
x_min = np.min(x_for_plot)
X = [x - x_min for x in x_for_plot]
# 绘制折线图
plt.plot(X, Content)
plt.xlabel('location') # 设置x轴标签
plt.ylabel('content') # 设置y轴标签
plt.title('Line Plot') # 设置图标题
plt.show()
#总用时： 192s按照我们给定的切片间隔，将输出周长随高度变化的折线图：
从结果中可以看出，我们的算法基本上解决了这个问题。
扫一扫
447
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
