近似平面点云一般特征——粗糙度

粗糙度

由前文《点云的凹凸性检验》

点云凹凸性检验（2/2）_三尺流流的博客-CSDN博客简单描绘点云里的“沟壑”和“山峰”。有的杂点请自行处理。https://blog.csdn.net/qq_55433334/article/details/125255479?spm=1001.2014.3001.5502的实验结果，我们可以看到：即便是在光滑的凸起表面，也存在凹陷的部分杂点。这些凹陷与凸起交错的部分表示：这部分的点云虽然是描述的是凸起体，但是仍然是表面较粗糙的凸起体。

问：

 现有一近似平面的粗糙点云，需要以此对应平面为基准，使用数值定量描述此平面点云的表面粗糙度。

粗糙度表征量

描述平面粗糙度时，需要先确定基准平面。这里推荐大佬博主M＆Q的最小二乘法：

【MQ笔记】超简单的最小二乘法拟合平面（Python）_M＆Q的博客-CSDN博客_最小二乘法拟合平面

在我们用二乘法进行平面拟合的时候，即存在处于平面上方的点也有处于平面下方的点，这种起伏的大小决定了平面是否是粗糙的。Sa为最老牌的三维表面粗糙度参数之一,此粗糙度描述的是凹凸不平的模型体积除以底面积得到的平均高度。

由公式：

可得到复杂模型相对于基准平面的算术平均粗糙度，我们目前暂时以这个值作为平面是否光滑的表征量。但是点云是一种“数字信号”，而非“模拟信号”直接调用平均粗糙度公式时我们是缺少可以积分的函数Z（x，y）的。

此时，依据点云对此公式做一点点修改：

逻辑流

1.导入近似平面的点云。

2.最小二乘法拟合粗糙度基准平面。

2.计算原点云到基准平面点云投影的距离。

3.调用公式计算点云距离的均值。

完整代码

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


def roughness_projection(point_list, plane_list):

    sa_list = []
    for i in range(len(plane_list)):
        A, B, C, D = plane_list[i]
        f_a, f_b, f_c, f_d = -A, -B, -C, -D

        x = point_list[i][:, 0]
        y = point_list[i][:, 1]
        z = point_list[i][:, 2]


        distance = []

        xp = []
        yp = []
        zp = []
        for j in range(len(x)):
            xp.append(
                ((f_b ** 2 + f_c ** 2) * x[j] - f_a * (f_b * y[j] + f_c * z[j] + f_d)) / (f_a ** 2 + f_b ** 2 + f_c ** 2))
            yp.append(
                ((f_a ** 2 + f_c ** 2) * y[j] - f_b * (f_a * x[j] + f_c * z[j] + f_d)) / (f_a ** 2 + f_b ** 2 + f_c ** 2))
            zp.append(
                ((f_a ** 2 + f_b ** 2) * z[j] - f_c * (f_a * x[j] + f_b * y[j] + f_d)) / (f_a ** 2 + f_b ** 2 + f_c ** 2))
            distance.append(((x[j] - xp[j]) ** 2 + (y[j] - yp[j]) ** 2 + (z[j] - zp[j]) ** 2) ** (0.5))
        sa = mean(distance)
        sa_list.append(sa)
    return sa_list

def roughness_view(point_list, plane_list):

    fig1 = plt.figure()
    ax1 = fig1.add_subplot(111, projection='3d')
    ax1.set_xlabel("x")
    ax1.set_ylabel("y")
    ax1.set_zlabel("z")

    for i in range(len(plane_list)):
        x = point_list[i][:, 0]
        y = point_list[i][:, 1]
        z = point_list[i][:, 2]
        A, B, C, D = plane_list[i]
        a, b, c, d = A, B, -C, D

        ax1.scatter(x, y, z, color='gray', marker='.')

        x_p = np.linspace(min(x),max(x), 100)
        y_p = np.linspace(min(y),max(y), 100)
        x_p, y_p = np.meshgrid(x_p, y_p)
        z_p =a * x_p + b * y_p + d
        ax1.plot_wireframe(x_p, y_p, z_p, rstride=10, cstride=10)

    plt.show()
    return None



if __name__ == '__main__':
    # 拟合平面：z = 0.559 * x + 0.164 * y + -11.822
    plane_list = [[0.559, 0.164, -1, -11.822]]
    point_list = []
    pcd = o3d.io.read_point_cloud("D:\1.pcd")
    pcd = np.asarray(pcd.points)
    point_list.append(pcd)

    sa_list = roughness_projection(point_list, plane_list)
    print("sa算术平均粗糙度：",sa_list)
    roughness_view(point_list, plane_list)

 

 即为此近似平面点云的算术平均粗糙度。