python - 加速 3D 点云中的 RANSAC 平面检测
问题描述
现在我正在使用 RANSAC 对 3D 点云数据进行平面分割。我的代码的基本流程图是:
- 选择 3 个随机点然后创建一个候选平面
- 检查到平面一定距离阈值内的所有其他点。如果有很多点低于阈值,那么我会保存飞机。
- 如果覆盖点低于阈值,则选择另一个 RANSAC 点(回到 no 1)
令我沮丧的是,如何加快进程?想象一下,如果我有 1 个 mio 点,那么我必须检查所有指向我的候选平面的点,如果不符合阈值,那么我选择另一个随机点并一遍又一遍地检查 1 个 mio 点。对于我的代码的任何建议甚至更正,我将非常感谢
这是我的代码,数据集可以在这里下载。对于这个示例,我使用了一小部分真实数据。
import numpy as np
import pandas as pd
import math
from scipy.spatial import cKDTree
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import random
data = pd.read_csv('data_sample.txt', usecols=[0,1,2], header=None, delimiter=' ')
df_points = pd.DataFrame(data)
%matplotlib qt
#Change parameter
distance_threshold = 0.06
minimum_point_threshold = 400
sampling_resolution = 0.06
#search nearest neighbor
kdtree = cKDTree(df_points)
temp_pair_point = []
list_pair_point = []
list_candidates = []
#Select random points
for h in df_points.index:
pair_point = kdtree.query_ball_point(np.array(df_points.iloc[h]),r=sampling_resolution, return_sorted=False)
temp_pair_point.append(pair_point)
for k in range(0,len(temp_pair_point)):
pairsize = len(temp_pair_point[k])
if pairsize > 10:
list_pair_point.append(temp_pair_point[k])
for m in range(0,len(list_pair_point)):
candidates = random.choices(list_pair_point[m], k=3)
list_candidates.append(candidates)
#select random point from nearest neighbor
list_sample_point = []
for i in range(0,len(list_candidates)):
list_test_distance = []
p1 = np.array(data.iloc[list_candidates[i][0]])
p2 = np.array(data.iloc[list_candidates[i][1]])
p3 = np.array(data.iloc[list_candidates[i][2]])
sample = [p1,p2,p3]
#create plane
x1,y1,z1 =np.ravel(p1)
x2,y2,z2 =np.ravel(p2)
x3,y3,z3 =np.ravel(p3)
m_u = [x2-x1,y2-y1,z2-z1]
m_v = [x3-x2,y3-y2,z3-z2]
m_normal = np.cross(m_u,m_v)
m_pos = p1
m_dist = np.dot(-m_normal, m_pos)
for j in data.index:
test_point = np.array(data.iloc[j])
test_distance = np.dot(-m_normal,test_point)
#print('Distance ={}'.format(test_distance))
if abs(test_distance) < distance_threshold:
list_test_distance.append(test_distance)
if len(list_test_distance) > minimum_point_threshold:
list_sample_point.append(sample)
#Check candidate plane
list_u = []
list_v = []
list_normal = []
list_distance = []
list_min_max = []
for i in range(0,len(list_sample_point)):
psample = list_sample_point[i]
psample1 = np.array(psample[0])
psample2 = np.array(psample[1])
psample3 = np.array(psample[2])
min_max = [float(min(np.transpose(psample)[0])), float(max(np.transpose(psample)[0])),float(min(np.transpose(psample)[1])), float(max(np.transpose(psample)[1])),float(min(np.transpose(psample)[2])), float(max(np.transpose(psample)[2]))]
list_min_max.append(min_max)
u = [psample2[0]-psample1[0], psample2[1]-psample1[1], psample2[2]-psample1[2]]
v = [psample3[0]-psample2[0], psample3[1]-psample2[1], psample3[2]-psample2[2]]
normal = np.cross(u,v)
d = np.dot(-normal,psample1)
list_u.append(u)
list_v.append(v)
list_normal.append(normal)
list_distance.append(d)
#Plot candidate plane
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
#ax.scatter3D(xs=plot_x, ys=plot_y, zs=plot_z, s=30, c='red')
ax.set_xlim3d(-25,-20)
ax.set_ylim3d(5,10)
ax.set_xlim3d(5,10)
ax.scatter3D(xs=df_points[0], ys=df_points[1], zs=df_points[2], s=1, c='green')
for j in range(0,len(list_normal)):
if list_normal[j][2] != 0:
print('Normal plane to Z axis')
xx, yy = np.meshgrid(np.arange(list_min_max[j][0],list_min_max[j][1]+(list_min_max[j][1]-list_min_max[j][0]),(list_min_max[j][1]-list_min_max[j][0])), np.arange(list_min_max[j][2],list_min_max[j][3]+(list_min_max[j][3]-list_min_max[j][2]),(list_min_max[j][3]-list_min_max[j][2])), sparse=True)
z = (-list_normal[j][0] * xx - list_normal[j][1] * yy -list_distance[j]) / list_normal[j][2]
ax.plot_surface(xx, yy, z, color=np.random.rand(3,), alpha=0.5)
elif list_normal[j][0] != 0:
print('Normal plane to X axis')
yy, zz, = np.meshgrid(np.arange(list_min_max[j][2],list_min_max[j][3]+(list_min_max[j][3]-list_min_max[j][2]),(list_min_max[j][3]-list_min_max[j][2])), np.arange(list_min_max[j][4],list_min_max[j][5]+(list_min_max[j][5]-list_min_max[j][4]),(list_min_max[j][5]-list_min_max[j][4])), sparse=True)
x = (-list_normal[j][1] * yy - list_normal[j][2] * zz -list_distance[j]) / list_normal[j][0]
ax.plot_surface(x, yy, zz, color=np.random.rand(3,), alpha=0.5)
elif list_normal[j][1] !=0:
print('Normal plane to Y axis')
xx, zz, = np.meshgrid(np.arange(list_min_max[j][0],list_min_max[j][1]+(list_min_max[j][1]-list_min_max[j][0]),(list_min_max[j][1]-list_min_max[j][0])), np.arange(list_min_max[j][4],list_min_max[j][5]+(list_min_max[j][5]-list_min_max[j][4]),(list_min_max[j][5]-list_min_max[j][4])), sparse=True)
y = (-list_normal[j][0] * xx - list_normal[j][2] * zz -list_distance[j]) / list_normal[j][1]
ax.plot_surface(xx, y, zz, color=np.random.rand(3,), alpha=0.5)
解决方案
首先,RANSAC 并不是在每个点上都运行。使用您的代码,它会从您拥有的 1945 个点中计算出 1729 个平面。通常,在启动 RANSAC 之前,您会计算要运行的迭代次数。这取决于预期平面的内点和异常点的数量:
k = log(1-p)/log(1-(1-e)^s)
p 是您想要的概率,即 RANSAC 会导致一个离群值自由平面,例如 99.99% 的概率。e 是异常值的百分比,例如 80% 异常值e = 0.8
。这导致 1147 次迭代。但是通常你应该能够假设一个较低的异常概率,假设是 20%,因为选择了邻居(就像你一样)现在突然只有 13 次迭代来检测一个概率为 99.99% 的平面,所以只有 13(! ) 而不是 1729(!)。如果您想找到 100 架飞机,它仍然可以将您的迭代次数减少约 25%。
其次,如果您仍然想在每个点上运行它,您可以固定代码的某些部分:首先,您使用两个 for 循环来处理可以写成一个的东西,并且您继续计算test_distance
即使len(list_test_distance) > minimum_point_threshold
可能已经是真的。通过将代码更改为:
#select random point from nearest neighbor
list_sample_point = []
list_u = []
list_v = []
list_normal = []
list_distance = []
list_min_max = []
for i in range(0,len(list_candidates)):
list_test_distance = []
p1 = np.array(data.iloc[list_candidates[i][0]])
p2 = np.array(data.iloc[list_candidates[i][1]])
p3 = np.array(data.iloc[list_candidates[i][2]])
sample = [p1,p2,p3]
#create plane
x1,y1,z1 =np.ravel(p1)
x2,y2,z2 =np.ravel(p2)
x3,y3,z3 =np.ravel(p3)
m_u = [x2-x1,y2-y1,z2-z1]
m_v = [x3-x2,y3-y2,z3-z2]
m_normal = np.cross(m_u,m_v)
m_pos = p1
m_dist = np.dot(-m_normal, m_pos)
for j in data.index:
test_point = np.array(data.iloc[j])
test_distance = np.dot(-m_normal,test_point)
#print('Distance ={}'.format(test_distance))
if abs(test_distance) < distance_threshold:
list_test_distance.append(test_distance)
if len(list_test_distance) > minimum_point_threshold:
#list_sample_point.append(sample) I think you don't need this one any more as well
psample = sample
psample1 = np.array(psample[0])
psample2 = np.array(psample[1])
psample3 = np.array(psample[2])
min_max = [float(min(np.transpose(psample)[0])), float(max(np.transpose(psample)[0])),float(min(np.transpose(psample)[1])), float(max(np.transpose(psample)[1])),float(min(np.transpose(psample)[2])), float(max(np.transpose(psample)[2]))]
list_min_max.append(min_max)
u = [psample2[0]-psample1[0], psample2[1]-psample1[1], psample2[2]-psample1[2]]
v = [psample3[0]-psample2[0], psample3[1]-psample2[1], psample3[2]-psample2[2]]
normal = np.cross(u,v)
d = np.dot(-normal,psample1)
list_u.append(u)
list_v.append(v)
list_normal.append(normal)
list_distance.append(d)
break
我将笔记本电脑上的运行时间从 7 分 11 秒更改为 2 分 35 秒,使其速度提高了 2.5 倍以上。结果应该仍然是相同的(纠正我,如果我错了)
推荐阅读
- laravel - 登录后laravel 7 csrf令牌更改
- python - SendBird 错误 400403 无效值:JSON 正文
- python - 定义和求解无冗余方程(Python)
- scala - 使用功能构造的组合批处理
- flutter - Flutter - Ontap InkWell 在 Stack 中不起作用
- c# - C# 应用程序是否可以读取其他应用程序发送的通知并相应地执行操作?
- swift - videoGravity resizeAspectFill 不起作用 - Swift - AVFoundation
- javascript - MapBox 通过 props 设置多个标记 react
- dns - 如何在没有网页的情况下拥有域名?
- sql - oracle中如何通过不同角色在其他列上获取不同列的值