如何解决平面噪声点云边缘锐化多次运行结果不一致问题？

我有一个平面的带噪声点云，想要估算其尺寸，因此需要锐化边缘。我采用基于直方图均值的统计方法过滤边缘点，具体实现代码如下：

import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt

def filter_cloud_borders(pcl:o3d.geometry.PointCloud, visu:bool=True) -> o3d.geometry.PointCloud:
    #
    min_bbox = pcl.get_minimal_oriented_bounding_box()
    # Box Pose
    box_pose = np.eye(4)
    box_pose[:3, :3] = min_bbox.R
    box_pose[:3, 3] = min_bbox.center
    # align it to center
    pcl.transform(np.linalg.inv(box_pose))
    pcl_pts = np.asarray(pcl.points)

    # compute the histogram of X and Y values separately
    hist_x, bins_x = np.histogram(pcl_pts[:, 0], bins=500)
    hist_y, bins_y = np.histogram(pcl_pts[:, 1], bins=500)

    # compute the mean of the histogram values for X and Y separately
    mean_hist_x = np.mean(hist_x)
    mean_hist_y = np.mean(hist_y)

    # find indices where histogram values are greater than the mean for X and Y
    ## TODO: adjust these limits properly to avoid over cropping
    crop_factor = 0.9
    indices_greater_than_mean_x = np.where(hist_x > crop_factor*mean_hist_x)[0]
    indices_greater_than_mean_y = np.where(hist_y > crop_factor*mean_hist_y)[0]

    # Filter point cloud points based on X and Y values
    filtered_points = pcl_pts[
        (pcl_pts[:, 0] > bins_x[indices_greater_than_mean_x[0]]) & 
        (pcl_pts[:, 0] < bins_x[indices_greater_than_mean_x[-1]]) &
        (pcl_pts[:, 1] > bins_y[indices_greater_than_mean_y[0]]) &
        (pcl_pts[:, 1] < bins_y[indices_greater_than_mean_y[-1]])
    ]

    # apply transformation to the filtered points
    transformed_points = np.hstack((filtered_points, np.ones((len(filtered_points), 1))))  # Convert points to homogeneous coordinates
    transformed_points = np.dot(box_pose, transformed_points.T).T[:, :3]  # Apply transformation

    new_cloud = o3d.geometry.PointCloud()
    new_cloud.points = o3d.utility.Vector3dVector(transformed_points)
    new_cloud.estimate_normals()
    #
    # pcl.transform(box_pose)
    #
    if(visu):
        plt.plot(hist_x)
        plt.axhline(y=mean_hist_x, color="g", linestyle="-")
        plt.axvline(x=indices_greater_than_mean_x[0], color="m", linestyle="-")
        plt.axvline(x=indices_greater_than_mean_x[-1], color="r", linestyle="-")
        plt.show()

        plt.plot(hist_y)
        plt.axhline(y=mean_hist_y, color="g", linestyle="-")
        plt.axvline(x=indices_greater_than_mean_y[0], color="m", linestyle="-")
        plt.axvline(x=indices_greater_than_mean_y[-1], color="r", linestyle="-")
        plt.show()
        # Origin
        origin = o3d.geometry.TriangleMesh.create_coordinate_frame(
            size=0.1, origin=[0, 0, 0])
        # colors
        min_bbox.color = [0.3, 0.4, 0.5]
        new_cloud.paint_uniform_color([0.5, 0.3, 0.7])
        ## Visualization
        o3d.visualization.draw_geometries([origin, pcl, new_cloud, min_bbox])
    
    return new_cloud

但我发现同一数据集多次运行时，得到的直方图分布并不一致！我怀疑这是因为每次运行时最小方向包围盒（尤其是方向）不同导致的，但不确定。请问该如何解决这个问题？

你的怀疑是正确的：Open3D的get_minimal_oriented_bounding_box()处理对称或近似对称的平面点云时，确实可能返回不同朝向的包围盒，导致点云对齐后的X/Y轴方向随机变化，最终让直方图分布每次都不一样。

解决核心是固定点云对齐后的坐标系，确保每次运行的基准一致，以下是两种可靠的实现方法：

方法1：基于平面法线固定坐标系（推荐） #

既然是平面点云，先估算平面法线并将其固定为Z轴，再将点云对齐到这个固定坐标系，后续用轴对齐包围盒替代方向包围盒，彻底避免朝向随机问题：

import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt

def filter_cloud_borders(pcl: o3d.geometry.PointCloud, visu: bool = True) -> o3d.geometry.PointCloud:
    # 1. 估算平面法线并固定坐标系
    pcl.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30))
    # 取平均法线作为平面基准方向（归一化确保一致性）
    avg_normal = np.mean(np.asarray(pcl.normals), axis=0)
    avg_normal = avg_normal / np.linalg.norm(avg_normal)
    
    # 构建固定坐标系：Z轴为法线，X/Y轴正交化生成
    z_axis = avg_normal
    # 选一个与Z轴不共线的向量作为初始X轴，确保正交
    x_axis = np.array([1, 0, 0]) if abs(np.dot(z_axis, np.array([1,0,0]))) < 0.9 else np.array([0,1,0])
    x_axis = x_axis - np.dot(x_axis, z_axis) * z_axis  # 格拉姆-施密特正交化
    x_axis = x_axis / np.linalg.norm(x_axis)
    y_axis = np.cross(z_axis, x_axis)  # Y轴由叉乘生成
    
    # 构建变换矩阵：将点云中心移到原点并对齐到固定坐标系
    center = np.mean(np.asarray(pcl.points), axis=0)
    transform = np.eye(4)
    transform[:3, 0] = x_axis
    transform[:3, 1] = y_axis
    transform[:3, 2] = z_axis
    transform[:3, 3] = center
    inv_transform = np.linalg.inv(transform)
    
    # 对齐点云到固定坐标系
    pcl_aligned = pcl.transform(inv_transform)
    pcl_pts = np.asarray(pcl_aligned.points)

    # 2. 原直方图过滤逻辑保持不变
    hist_x, bins_x = np.histogram(pcl_pts[:, 0], bins=500)
    hist_y, bins_y = np.histogram(pcl_pts[:, 1], bins=500)

    mean_hist_x = np.mean(hist_x)
    mean_hist_y = np.mean(hist_y)

    crop_factor = 0.9
    indices_greater_than_mean_x = np.where(hist_x > crop_factor * mean_hist_x)[0]
    indices_greater_than_mean_y = np.where(hist_y > crop_factor * mean_hist_y)[0]

    filtered_points = pcl_pts[
        (pcl_pts[:, 0] > bins_x[indices_greater_than_mean_x[0]]) & 
        (pcl_pts[:, 0] < bins_x[indices_greater_than_mean_x[-1]]) &
        (pcl_pts[:, 1] > bins_y[indices_greater_than_mean_y[0]]) &
        (pcl_pts[:, 1] < bins_y[indices_greater_than_mean_y[-1]])
    ]

    # 3. 将过滤后的点云转换回原坐标系
    transformed_points = np.hstack((filtered_points, np.ones((len(filtered_points), 1))))
    transformed_points = np.dot(transform, transformed_points.T).T[:, :3]

    new_cloud = o3d.geometry.PointCloud()
    new_cloud.points = o3d.utility.Vector3dVector(transformed_points)
    new_cloud.estimate_normals()

    if visu:
        plt.plot(hist_x)
        plt.axhline(y=mean_hist_x, color="g", linestyle="-")
        plt.axvline(x=indices_greater_than_mean_x[0], color="m", linestyle="-")
        plt.axvline(x=indices_greater_than_mean_x[-1], color="r", linestyle="-")
        plt.show()

        plt.plot(hist_y)
        plt.axhline(y=mean_hist_y, color="g", linestyle="-")
        plt.axvline(x=indices_greater_than_mean_y[0], color="m", linestyle="-")
        plt.axvline(x=indices_greater_than_mean_y[-1], color="r", linestyle="-")
        plt.show()

        origin = o3d.geometry.TriangleMesh.create_coordinate_frame(size=0.1, origin=[0, 0, 0])
        aligned_bbox = pcl_aligned.get_axis_aligned_bounding_box()
        aligned_bbox.color = [0.3, 0.4, 0.5]
        new_cloud.paint_uniform_color([0.5, 0.3, 0.7])
        o3d.visualization.draw_geometries([origin, pcl, new_cloud, aligned_bbox.transform(transform)])
    
    return new_cloud

方法2：强制固定最小包围盒朝向 #

如果必须使用最小方向包围盒，可以通过保存初始旋转矩阵作为参考，后续将新计算的包围盒与参考对齐：

1. 第一次运行时，保存min_bbox.R作为基准旋转矩阵
2. 后续运行时，通过四元数匹配将新的旋转矩阵对齐到基准方向

不过这种方法可靠性低于方法1，因为对称点云的最小包围盒主轴本身存在歧义。

额外优化建议 #

• 直方图的bins参数建议根据点云实际尺寸动态计算（比如按每毫米一个bin），避免坐标系缩放导致分辨率不稳定
• 用中位数替代均值作为过滤阈值，减少噪声对直方图统计的干扰

内容的提问来源于stack exchange，提问作者bhomaidan90