（超级详细）基于numpy-stl，对不规则柱状实体旋转并求沿轴切片的周长_numpy stl mesh教程

一、前言

本文是在我前一篇博文的基础上进一步扩展，除了对实体进行切片，求周长，还新增了旋转实体的操作。

（源代码）基于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 time

2.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

按照我们给定的切片间隔，将输出周长随高度变化的折线图：

从结果中可以看出，我们的算法基本上解决了这个问题。