【点云】生成有凹凸的平面

文章目录

• 前言
• 高斯函数
• • 原理
  • 代码
  • 保存
• 测试
• • 测试1 ：领域曲率代码
  • 测试2：高斯曲率代码
• 加上噪点
• • 测试1
  • 测试2
  • 改进

前言

尝试用一些数据生成有凹凸面的点云。
我们姑且把z轴当成有凹凸的缺陷，x轴和y轴共同组成一个平面。

高斯函数

原理

高斯函数wiki中，我们得知

其中，σ为标准差，用来控制“钟形”的宽度。

根据wiki中下面的举例sigma_X = 1;sigma_Y = 2;可以看出，σx=σy时，高斯的水平集是个圆，σx不等于σy时，高斯的水平集是个椭圆。（可以这样想：在平面上的公式，圆和椭圆的区别）

继续往下看：

代码

再根据wiki中下面的 Octave 代码，仿写出python代码：

# 导入 numpy 和 open3d 库import numpy as npimport open3d as o3d# 定义高斯函数的参数A = 1x0 = 0y0 = 0sigma_X = 1sigma_Y = 2# 生成 X 和 Y 的坐标网格X, Y = np.meshgrid(np.arange(-5, 5.1, 0.1), np.arange(-5, 5.1, 0.1))#-5到5，步长为0.1# 创建 open3d 点云对象pcd = o3d.geometry.PointCloud()# 循环旋转角度for theta in np.arange(0, np.pi, np.pi / 100): # 0 到 π，步长为 π / 100#可以改变这个值# 计算高斯函数的系数a = np.cos(theta) ** 2 / (2 * sigma_X ** 2) + np.sin(theta) ** 2 / (2 * sigma_Y ** 2)b = np.sin(2 * theta) / (4 * sigma_X ** 2) - np.sin(2 * theta) / (4 * sigma_Y ** 2)c = np.sin(theta) ** 2 / (2 * sigma_X ** 2) + np.cos(theta) ** 2 / (2 * sigma_Y ** 2)# 计算 Z 的坐标Z = A * np.exp(-(a * (X - x0) ** 2 + 2 * b * (X - x0) * (Y - y0) + c * (Y - y0) ** 2))# 将 X, Y, Z 合并为点云矩阵，形状为 (n, 3)points = np.stack((X, Y, Z), axis=-1)points = points.reshape(-1, 3)# 更新点云的坐标pcd.points = o3d.utility.Vector3dVector(points)# 添加坐标coord = o3d.geometry.TriangleMesh.create_coordinate_frame(size=1, origin=[0, 0, 0])#x红色，y绿色，z蓝色# 可视化点云o3d.visualization.draw_geometries([pcd, coord])

得到

若我们改为sigma_X = 1，sigma_Y = 1，则

发现中间确实为圆，与上述猜想一致。

若我们需要凹陷的缺陷，则改为A=-1即可。

保存

# 保存点云o3d.io.write_point_cloud("flaw.pcd",pcd ) 

点云大小如下：

测试

【最详解】如何进行点云的凹凸缺陷检测（opene3D）
拿出之前写的凹凸检测代码开始测试，首先测试上述这种无噪音的。
记得一定要根据点云的大小改radius = 0.5 #邻域半径,否则一点效果也没有

测试1 ：领域曲率代码

参数如下：

结果如下：

意外的还算不错。

测试2：高斯曲率代码

也是改了radius =0.5。

果然结果还是这个更好。

加上噪点

在之前代码的基础上更改如下，并改成椭圆形缺陷。

结果：

测试1

测试2

发现在针对椭圆形的凹凸缺陷都不够灵敏了。

改进

一开始以为在使用邻近搜索中，用的方法不太好，用的是在球内的点搜索，或许换个方法就可以了。–2024.2.17
后发现修改为如下，依旧不太好。

k, idx, _ = kdtree.search_knn_vector_3d(cloud.points[i], num_knn) 

后发现在打印出数据中，数据太小，而定义为平面的数据宽泛太大，于是在测试2中改为如下

limit_max = 1e-3for i in range(len(curvatures)):if -limit_max<curvatures[i][0] < limit_max and -limit_max<curvatures[i][1] <limit_max: #平坦np.asarray(pcd.colors)[i] = [0, 0, 0]#黑elif -limit_max<curvatures[i][0] < limit_max and curvatures[i][1] >limit_max:#凸np.asarray(pcd.colors)[i] = [1, 0, 0]#红elif -limit_max<curvatures[i][0] < limit_max and -limit_max<curvatures[i][1] <limit_max: #凹np.asarray(pcd.colors)[i] = [0, 1, 0]#绿elif curvatures[i][0] < -limit_max and curvatures[i][1] >limit_max: #鞍形脊 大部分凸，少部分凹np.asarray(pcd.colors)[i] = [0, 0, 1]#蓝elif curvatures[i][0] < -limit_max and curvatures[i][1] <-limit_max: #鞍形谷 大部分凹，少部分凸np.asarray(pcd.colors)[i] = [0, 1, 1]#青elif curvatures[i][0] > limit_max and curvatures[i][1] >limit_max: #峰 np.asarray(pcd.colors)[i] = [1, 0, 1]#紫elif curvatures[i][0] > limit_max and curvatures[i][1] <-limit_max: #阱np.asarray(pcd.colors)[i] = [1, 1, 0]#黄

结果便好很多了。