import numpy as np
import open3d as o3d
import matplotlib.pyplot as plt
import os
from matplotlib.gridspec import GridSpec
from tornado.gen import Return
# 设置中文字体为宋体
plt.rcParams["font.family"] = ["SimSun"]
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
class DirectionalCompressionPointCloudProcessor:
def __init__(self, random_seed=42, output_dir="output_models",
compression_factor=0.5, compression_direction=np.array([0, 0, 1])):
"""初始化点云处理器,支持沿任意方向向量压缩"""
self.random_seed = random_seed
self.output_dir = output_dir # 输出目录
self.compression_factor = compression_factor # 压缩因子,1.0表示不压缩,值越小压缩程度越高
self.compression_direction = self._normalize_vector(compression_direction) # 压缩方向向量(已归一化)
np.random.seed(random_seed)
# 创建输出目录(如果不存在)
os.makedirs(self.output_dir, exist_ok=True)
# 初始化数据结构
self.data = None
self.points = None
self.labels = None
self.unique_labels = {}
# 可视化几何体存储
self.geometries = []
self.geometries_pcd = []
# 处理参数
self.cluster_index = 0
self.min_points_threshold = 20
self.small_cluster_threshold = 10
self.outlier_std_ratio = 1.2
# 平滑处理参数
self.smoothing_iterations = 5 # 平滑迭代次数
self.smoothing_lambda = 0.5 # 平滑强度参数
# 法线计算参数 - 使用固定KNN参数
self.knn_param = 20 # 固定KNN参数
self.consistent_knn = 10 # 法线一致性检查的邻点数
# 可视化参数
self.visualization_pause_time = 0 # 可视化暂停时间,0表示无限等待
self.coord_offset_ratio = 0.3 # 坐标系偏移比例
self.coord_position_choice = 0 # 0:左下, 1:右下, 2:左上, 3:右上
# 三视图参数
self.viewpoint_size = 100 # 三视图点大小
self.save_views = True # 是否保存三视图
def _normalize_vector(self, vector):
"""归一化向量"""
norm = np.linalg.norm(vector)
if norm == 0:
return np.array([0, 0, 1]) # 默认Z轴方向
return vector / norm
def set_compression_direction(self, direction):
"""设置压缩方向向量"""
self.compression_direction = self._normalize_vector(direction)
print(f"已设置压缩方向向量: {self.compression_direction}")
def load_data(self, file_path):
"""加载点云数据"""
try:
self.data = np.genfromtxt(
file_path,
delimiter=',',
skip_header=1
)
self.points = self.data[:, :3]
self.labels = self.data[:, 3].astype(int)
for label in np.unique(self.labels):
self.unique_labels[label] = self.points[self.labels == label]
# 计算点云整体边界框,用于放置坐标系
self._calculate_overall_bounding_box()
print(f"成功加载数据,共 {len(self.points)} 个点,{len(self.unique_labels)} 个类簇")
return True
except Exception as e:
print(f"数据加载失败: {str(e)}")
return False
def _calculate_overall_bounding_box(self):
"""计算整个点云的边界框,用于确定坐标系位置"""
if self.points is not None and len(self.points) > 0:
self.overall_min = np.min(self.points, axis=0)
self.overall_max = np.max(self.points, axis=0)
self.overall_center = np.mean(self.points, axis=0)
self.overall_size = self.overall_max - self.overall_min
else:
self.overall_min = np.array([0, 0, 0])
self.overall_max = np.array([100, 100, 100])
self.overall_center = np.array([50, 50, 50])
self.overall_size = np.array([100, 100, 100])
def _filter_outliers(self, pcd):
"""过滤离群点"""
cl, ind = pcd.remove_statistical_outlier(
nb_neighbors=15,
std_ratio=self.outlier_std_ratio
)
return cl, ind
def _orient_normals_externally(self, pcd):
"""将法线方向统一向外"""
points = np.asarray(pcd.points)
normals = np.asarray(pcd.normals)
# 使用包围盒中心方法
bbox = pcd.get_axis_aligned_bounding_box()
center = bbox.get_center()
to_center = points - center
dot_products = np.sum(normals * to_center, axis=1)
flip_mask = dot_products > 0 # 指向中心的需要翻转
normals[flip_mask] *= -1
pcd.normals = o3d.utility.Vector3dVector(normals)
return pcd
def _create_colored_coordinate_frame(self, size=1.0):
"""创建自定义颜色的坐标系:x轴蓝色,y轴绿色,z轴红色"""
frame = o3d.geometry.TriangleMesh.create_coordinate_frame(size=size)
# 设置坐标轴颜色
frame.vertex_colors[0] = [0, 0, 1] # x轴顶点1 - 蓝色
frame.vertex_colors[1] = [0, 0, 1] # x轴顶点2 - 蓝色
frame.vertex_colors[2] = [0, 1, 0] # y轴顶点1 - 绿色
frame.vertex_colors[3] = [0, 1, 0] # y轴顶点2 - 绿色
frame.vertex_colors[4] = [1, 0, 0] # z轴顶点1 - 红色
frame.vertex_colors[5] = [1, 0, 0] # z轴顶点2 - 红色
# 根据选择计算不同的坐标系位置
offset = self.overall_size * self.coord_offset_ratio
if self.coord_position_choice == 0: # 左下
coord_position = self.overall_min - offset
elif self.coord_position_choice == 1: # 右下
coord_position = np.array([self.overall_max[0] + offset[0],
self.overall_min[1] - offset[1],
self.overall_min[2] - offset[2]])
elif self.coord_position_choice == 2: # 左上
coord_position = np.array([self.overall_min[0] - offset[0],
self.overall_max[1] + offset[1],
self.overall_min[2] - offset[2]])
else: # 右上
coord_position = self.overall_max + offset
# 平移坐标系到计算位置
frame.translate(coord_position)
return frame
def _laplacian_smoothing(self, mesh, iterations=5, lambda_factor=0.5):
"""对网格进行拉普拉斯平滑处理,减少尖锐边缘和凹陷"""
if not mesh:
return mesh
# 兼容旧版本Open3D,不使用clone()方法
# 创建网格的深拷贝以避免修改原始数据
smoothed_mesh = o3d.geometry.TriangleMesh()
smoothed_mesh.vertices = o3d.utility.Vector3dVector(np.asarray(mesh.vertices))
smoothed_mesh.triangles = o3d.utility.Vector3iVector(np.asarray(mesh.triangles))
smoothed_mesh.vertex_normals = o3d.utility.Vector3dVector(np.asarray(mesh.vertex_normals))
smoothed_mesh.triangle_normals = o3d.utility.Vector3dVector(np.asarray(mesh.triangle_normals))
# 获取顶点和邻接信息
vertices = np.asarray(smoothed_mesh.vertices)
adjacency_list = smoothed_mesh.adjacency_list
# 执行多轮平滑
for _ in range(iterations):
new_vertices = vertices.copy()
for i in range(len(vertices)):
# 确保 i 在 adjacency_list 的有效范围内
if i < len(adjacency_list) and adjacency_list[i]:
neighbors = adjacency_list[i]
# 计算相邻顶点的平均值
neighbor_avg = np.mean(vertices[neighbors], axis=0)
# 应用拉普拉斯平滑公式
new_vertices[i] = vertices[i] + lambda_factor * (neighbor_avg - vertices[i])
# 更新顶点位置
smoothed_mesh.vertices = o3d.utility.Vector3dVector(new_vertices)
# 重新计算法线
smoothed_mesh.compute_vertex_normals()
return smoothed_mesh
def _directional_compression(self, mesh, compression_factor=None, use_max_edge=True):
"""
基于单一边缘平面的方向压缩算法:
将方向向量作为法向量,以最边缘的一个顶点确定投影平面,
所有顶点都将相对于这个单一平面进行比例移动
参数:
compression_factor: 压缩因子,0-1之间的值表示压缩比例
use_max_edge: 如果为True,使用正方向最边缘的平面;如果为False,使用负方向最边缘的平面
"""
if not mesh:
return mesh
# 如果未指定压缩因子,使用类中定义的默认值
if compression_factor is None:
compression_factor = self.compression_factor
# 确保压缩因子在有效范围内
if compression_factor <= 0 or compression_factor > 1:
raise ValueError("压缩因子必须是0到1之间的值(不包括0,包括1)")
# 获取网格顶点
vertices = np.asarray(mesh.vertices)
# 确保压缩方向向量已归一化
plane_normal = self.compression_direction / np.linalg.norm(self.compression_direction)
# 计算所有顶点在压缩方向上的投影
projections = np.dot(vertices, plane_normal)
# 找到投影的最大值和最小值,确定沿压缩方向的边缘位置
min_proj = np.min(projections)
max_proj = np.max(projections)
original_thickness = max_proj - min_proj
if original_thickness <= 0: # 避免厚度为零的情况
return mesh
# 选择单一边缘平面作为投影基准
if use_max_edge:
# 使用正方向最边缘的平面
edge_proj = max_proj
# 平面方程:plane_normal · x + d = 0,其中d = -edge_proj
d = -edge_proj
# 计算每个顶点到边缘平面的带符号距离(负值表示在平面内侧)
distances = np.dot(vertices, plane_normal) + d
else:
# 使用负方向最边缘的平面
edge_proj = min_proj
d = -edge_proj
distances = np.dot(vertices, plane_normal) + d
# 处理顶点位置
new_vertices = []
for i, vertex in enumerate(vertices):
dist = distances[i]
# 应用压缩因子:新距离 = 原始距离 × 压缩因子
new_dist = dist * compression_factor
# 计算需要移动的距离
move_distance = new_dist - dist
# 计算移动向量
move_vector = move_distance * plane_normal
# 计算新顶点位置
new_vertex = vertex + move_vector
new_vertices.append(new_vertex)
# 更新顶点
mesh.vertices = o3d.utility.Vector3dVector(new_vertices)
# 平滑处理
mesh = self._laplacian_smoothing(
mesh,
iterations=self.smoothing_iterations,
lambda_factor=self.smoothing_lambda
)
# 重新计算法线
mesh.compute_vertex_normals()
return mesh
def _save_3d_model(self, mesh, cluster_index, label):
"""保存三维模型为多种格式"""
if not mesh:
return
# 构建文件名(包含类簇信息)
base_filename = f"Ore_cluster_{label}_sub_{cluster_index}"
# 保存为PLY格式(包含颜色信息)
ply_path = os.path.join(self.output_dir, f"{base_filename}_DirCompress.ply")
o3d.io.write_triangle_mesh(ply_path, mesh)
print(f"已保存PLY模型: {ply_path}")
# 保存为STL格式(广泛用于3D打印)
stl_path = os.path.join(self.output_dir, f"{base_filename}_DirCompress.stl")
o3d.io.write_triangle_mesh(stl_path, mesh)
print(f"已保存STL模型: {stl_path}")
# 保存为OBJ格式(在许多3D软件中兼容)
obj_path = os.path.join(self.output_dir, f"{base_filename}_DirCompress.obj")
o3d.io.write_triangle_mesh(obj_path, mesh)
print(f"已保存OBJ模型: {obj_path}")
def _process_sub_cluster(self, sub_cluster, cluster_index, label):
"""处理单个子簇并确保生成的网格为闭合体模型"""
# 添加微小噪声打破共面性
noise = np.random.normal(0, 1e-6, sub_cluster.shape)
sub_cluster += noise
# 创建点云对象
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(sub_cluster)
# 过滤点数太少的子簇
if len(pcd.points) < self.small_cluster_threshold:
print(f"子簇点数不足 {self.small_cluster_threshold},跳过处理")
return None, None
# 不过滤离群点
filtered_pcd = pcd
print(f"使用KNN参数计算法向量: {self.knn_param}个最近邻")
# 计算法向量 - 使用固定KNN参数
filtered_pcd.estimate_normals(
o3d.geometry.KDTreeSearchParamKNN(knn=self.knn_param),
fast_normal_computation=False # 禁用快速计算,提高精度
)
# 法线方向一致化
filtered_pcd.orient_normals_consistent_tangent_plane(k=self.consistent_knn)
# 设置颜色
color = plt.get_cmap("Accent")(cluster_index % 8)[:3]
filtered_pcd.paint_uniform_color(color)
# 泊松曲面重建
# try:
with o3d.utility.VerbosityContextManager(
o3d.utility.VerbosityLevel.Debug
) as cm:
mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(
filtered_pcd,
depth=5,
width= 0,
scale=1.1,
linear_fit=True
)
# 后处理:确保网格为闭合流形(兼容性更强)
if mesh is not None:
# 清理网格
mesh.remove_degenerate_triangles()
mesh.remove_duplicated_triangles()
mesh.remove_duplicated_vertices()
mesh.remove_unreferenced_vertices()
# 沿指定方向压缩模型
mesh = self._directional_compression(mesh)
# mesh = self.compress_mesh(mesh)
# 计算法线(确保可视化正常)
mesh.compute_vertex_normals()
mesh.paint_uniform_color(color)
# 保存三维模型
self._save_3d_model(mesh, cluster_index, label)
return filtered_pcd, mesh
# except Exception as e:
# print(f"泊松重建失败: {str(e)}")
# return filtered_pcd, None
def process_clusters(self):
"""处理所有类簇"""
if not self.unique_labels:
print("没有可处理的类簇数据")
return
# 遍历每个类簇
for label, cluster_points in self.unique_labels.items():
if label == -1: # 跳过噪声点
continue
if len(cluster_points) < self.min_points_threshold:
print(f"类簇 {label} 点数不足 {self.min_points_threshold},跳过处理")
continue
print(f"处理类簇 {label},包含 {len(cluster_points)} 个点")
current_clusters = self._get_sub_clusters(label, cluster_points)
for sub_cluster in current_clusters:
# 传递label参数用于保存文件
pcd, mesh = self._process_sub_cluster(sub_cluster, self.cluster_index, label)
if pcd:
self.geometries_pcd.append(pcd)
if mesh:
self.geometries.append(mesh)
self.cluster_index += 1
print(f"共处理 {self.cluster_index} 个有效类簇")
def _get_sub_clusters(self, label, cluster_points):
"""获取子簇,直接返回原始类簇"""
return [cluster_points]
def visualize_orthographic_views(self):
"""生成并显示三视图(正视图、侧视图、俯视图)"""
if self.points is None or len(self.points) == 0:
print("没有点云数据可生成三视图")
return
# 创建图形和子图
fig = plt.figure(figsize=(15, 10))
gs = GridSpec(2, 2, figure=fig)
# 正视图 (X-Y平面)
ax1 = fig.add_subplot(gs[0, 0])
# 侧视图 (Y-Z平面)
ax2 = fig.add_subplot(gs[0, 1])
# 俯视图 (X-Z平面)
ax3 = fig.add_subplot(gs[1, 0])
# 设置标题
ax1.set_title('正视图 (X-Y平面)')
ax2.set_title('侧视图 (Y-Z平面)')
ax3.set_title('俯视图 (X-Z平面)')
# 设置坐标轴标签
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax2.set_xlabel('Y')
ax2.set_ylabel('Z')
ax3.set_xlabel('X')
ax3.set_ylabel('Z')
# 确保各视图比例一致
ax1.set_aspect('equal')
ax2.set_aspect('equal')
ax3.set_aspect('equal')
# 获取颜色映射
cmap = plt.get_cmap("Accent")
# 按类别绘制点
for i, (label, points) in enumerate(self.unique_labels.items()):
if label == -1: # 跳过噪声点
continue
color = cmap(i % 8)[:3]
# 正视图 (X-Y)
ax1.scatter(points[:, 0], points[:, 1], c=[color], s=self.viewpoint_size, alpha=0.6,
label=f'类别 {label}')
# 侧视图 (Y-Z)
ax2.scatter(points[:, 1], points[:, 2], c=[color], s=self.viewpoint_size, alpha=0.6)
# 俯视图 (X-Z)
ax3.scatter(points[:, 0], points[:, 2], c=[color], s=self.viewpoint_size, alpha=0.6)
# 添加图例
ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
# 调整布局
plt.tight_layout()
# 保存三视图(如果启用)
if self.save_views:
views_path = os.path.join(self.output_dir, "point_cloud_orthographic_views.png")
plt.savefig(views_path, dpi=300, bbox_inches='tight')
print(f"已保存三视图: {views_path}")
# 显示三视图
plt.show()
def visualize_results(self):
"""可视化处理结果,包括3D视图和三视图"""
# 先显示三视图
self.visualize_orthographic_views()
if not self.geometries_pcd and not self.geometries:
print("没有可可视化的3D结果")
return
# 创建自定义颜色的坐标系
coord_frame_size = np.max(self.overall_size) * 0.1
colored_coord_frame = self._create_colored_coordinate_frame(size=coord_frame_size)
# 可视化点云(显示法线)
if self.geometries_pcd:
print("显示点云与法线方向 (按Q键关闭窗口)")
vis = o3d.visualization.Visualizer()
vis.create_window(window_name="点云与法线方向", width=1024, height=768)
# 添加自定义坐标系
vis.add_geometry(colored_coord_frame)
# 添加点云
for pcd in self.geometries_pcd:
vis.add_geometry(pcd)
# 设置渲染参数
opt = vis.get_render_option()
opt.point_size = 3
opt.line_width = 1
opt.show_coordinate_frame = False
# 运行可视化
vis.run()
vis.destroy_window()
# 可视化点云和重建网格
if self.geometries:
print("显示点云和泊松重建网格 (按Q键关闭窗口)")
vis = o3d.visualization.Visualizer()
vis.create_window(window_name="点云和泊松重建网格", width=1024, height=768)
# 添加自定义坐标系
vis.add_geometry(colored_coord_frame)
# 添加几何体
for geom in self.geometries + self.geometries_pcd:
vis.add_geometry(geom)
opt = vis.get_render_option()
opt.mesh_show_back_face = True
opt.point_size = 2
opt.show_coordinate_frame = False
vis.run()
vis.destroy_window()
def compress_mesh(self, mesh):
# 检查网格是否有效
if not mesh.has_vertices():
raise ValueError("STL 文件中没有顶点数据。")
# 2. 定义压缩方向
compression_direction = np.array([-0.9848, -0.1736, 0.0])
compression_direction = compression_direction / np.linalg.norm(compression_direction)
# 3. 定义压缩比例
scale_factor = 0.04 # 压缩比例
# 4. 获取顶点数据
vertices = np.asarray(mesh.vertices)
# 新增:计算顶点中心点,用于后续平移校正
center = np.mean(vertices, axis=0)
# 新增:将顶点平移到以原点为中心
vertices_centered = vertices - center
# 5. 对顶点进行投影与压缩
dot_products = np.dot(vertices_centered, compression_direction)
projections = np.outer(dot_products, compression_direction)
orthogonal_components = vertices_centered - projections
compressed_projections = projections * scale_factor
new_vertices_centered = orthogonal_components + compressed_projections
#将顶点平移回原来的位置
new_vertices = new_vertices_centered + center
# 6. 更新网格顶点
mesh.vertices = o3d.utility.Vector3dVector(new_vertices)
# 恢复法线计算(很重要)
mesh.compute_vertex_normals()
return mesh
if __name__ == "__main__":
# 定义压缩方向向量
compression_direction = np.array([-0.9848, -0.1736, 0]) # 短轴方向
# 创建处理器实例,指定输出目录和压缩参数
processor = DirectionalCompressionPointCloudProcessor(
output_dir="3d_models_New_output",
compression_factor= 0.04, # 压缩因子,0.3表示压缩到原来的30%
compression_direction=compression_direction
)
# 加载数据(请替换为实际文件路径)
file_path = 'midpoints_with_cluster_index_Ore2Part_filtered_10_602030.csv'
if processor.load_data(file_path):
# 可以在处理前动态更改压缩方向
# processor.set_compression_direction(np.array([0, 1, 0])) # 例如改为Y轴方向
# 处理类簇(处理过程中会自动保存模型)
processor.process_clusters()
# 可视化结果(现在包括三视图)
processor.visualize_results()
else:
print("无法继续处理,数据加载失败")
如何保证压缩后的网格模型尽量减少三角网交叉现象