非结构化点云的区域增长：验证失败 pcl
Region growing for unstructured pointclouds: failed verification pcl
我正在 python 中针对非结构化点云开发区域增长算法。我尝试根据 http://pointclouds.org/documentation/tutorials/region_growing_segmentation.php#region-growing-segmentation 和 SEGMENTATION OF POINT CLOUDS USING SMOOTHNESS CONSTRAINT.
实现算法
使用相同的点云，我无法使用 pcl 实现重现结果，尽管所有步骤和参数似乎都相同。我的实现产生了 217 个区域，而 pcl 的实现产生了 24 个（过度分割）。可以找到点云：
https://github.com/scrox/region-growing/blob/master/test.txt
是什么导致了这种显着差异？我也欢迎我的 python 代码的性能改进。
import math
import numpy as np
from sklearn.neighbors import KDTree
unique_rows=np.loadtxt("test.txt")
tree = KDTree(unique_rows, leaf_size=2)
dist,nn_glob = tree.query(unique_rows[:len(unique_rows)], k=30)
def normalsestimation(pointcloud,nn_glob,VP=[0,0,0]):
ViewPoint = np.array(VP)
normals = np.empty((np.shape(pointcloud)))
curv = np.empty((len(pointcloud),1))
for index in range(len(pointcloud)):
nn_loc = pointcloud[nn_glob[index]]
COV = np.cov(nn_loc,rowvar=False)
eigval, eigvec = np.linalg.eig(COV)
idx = np.argsort(eigval)
nor = eigvec[:,idx][:,0]
if nor.dot((ViewPoint-pointcloud[index,:])) > 0:
normals[index] = nor
else:
normals[index] = - nor
curv[index] = eigval[idx][0] / np.sum(eigval)
return normals,curv
def regiongrowing(pointcloud,nn_glob,theta_th = 'auto', cur_th = 'auto'):
normals,curvature = normalsestimation(pointcloud,nn_glob=nn_glob)
order = curvature[:,0].argsort().tolist()
region = []
if theta_th == 'auto':
theta_th = 15.0 / 180.0 * math.pi # in radians
if cur_th == 'auto':
cur_th = np.percentile(curvature,98)
while len(order) > 0:
region_cur = []
seed_cur = []
poi_min = order[0] #poi_order[0]
region_cur.append(poi_min)
seed_cur.append(poi_min)
order.remove(poi_min)
for i in range(len(seed_cur)):
nn_loc = nn_glob[seed_cur[i]]
for j in range(len(nn_loc)):
nn_cur = nn_loc[j]
if all([nn_cur in order , np.arccos(np.abs(np.dot(normals[seed_cur[i]],normals[nn_cur])))<theta_th]):
region_cur.append(nn_cur)
order.remove(nn_cur)
if curvature[nn_cur] < cur_th:
seed_cur.append(nn_cur)
region.append(region_cur)
return region
region = regiongrowing(unique_rows,nn_glob)
这里是 pcl 代码：
#include <iostream>
#include <vector>
#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/search/search.h>
#include <pcl/search/kdtree.h>
#include <pcl/features/normal_3d.h>
#include <pcl/visualization/cloud_viewer.h>
#include <pcl/segmentation/region_growing.h>
int
main (int argc, char** argv)
{
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
if ( pcl::io::loadPCDFile <pcl::PointXYZ> ("test.pcd", *cloud) == -1)
{
std::cout << "Cloud reading failed." << std::endl;
return (-1);
}
pcl::search::Search<pcl::PointXYZ>::Ptr tree = boost::shared_ptr<pcl::search::Search<pcl::PointXYZ> > (new pcl::search::KdTree<pcl::PointXYZ>);
pcl::PointCloud <pcl::Normal>::Ptr normals (new pcl::PointCloud <pcl::Normal>);
pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> normal_estimator;
normal_estimator.setSearchMethod (tree);
normal_estimator.setInputCloud (cloud);
normal_estimator.setKSearch (30);
normal_estimator.compute (*normals);
pcl::RegionGrowing<pcl::PointXYZ, pcl::Normal> reg;
reg.setSearchMethod (tree);
reg.setNumberOfNeighbours (30);
reg.setInputCloud (cloud);
reg.setInputNormals (normals);
reg.setSmoothnessThreshold (15.0 / 180.0 * M_PI);
reg.setCurvatureThreshold (0.088614);
std::vector <pcl::PointIndices> clusters;
reg.extract (clusters);
std::cout << "Number of clusters is equal to " << clusters.size () << std::endl;
std::cout << "First cluster has " << clusters[0].indices.size () << " points." << endl;
std::cout << "These are the indices of the points of the initial" <<
std::endl << "cloud that belong to the first cluster:" << std::endl;
int counter = 0;
while (counter < clusters[0].indices.size ())
{
std::cout << clusters[0].indices[counter] << ", ";
counter++;
}
std::cout << std::endl;
pcl::PointCloud <pcl::PointXYZRGB>::Ptr colored_cloud = reg.getColoredCloud ();
pcl::visualization::CloudViewer viewer ("Cluster viewer");
viewer.showCloud(colored_cloud);
while (!viewer.wasStopped ())
{
}
return (0);
}
我认为这与您 def regiongrowing(pointcloud,nn_glob,theta_th = 'auto', cur_th = 'auto'): 中的 for i in range(len(seed_cur)): 有关。我相信当 python 使用 for 循环时，它实际上并没有遍历循环开始后添加到列表中的任何元素。所以
if curvature[nn_cur] < cur_th:
seed_cur.append(nn_cur)
它不会继续循环 nn_cur。我使用 while 循环稍微修改了您的代码。我的结果是 21 个区域。与 PCL 中的情况不完全相同，但更接近。我的函数是regiongrowing1，供参考。
import math
import numpy as np
from sklearn.neighbors import KDTree
unique_rows=np.loadtxt("test.txt")
tree = KDTree(unique_rows, leaf_size=2)
dist,nn_glob = tree.query(unique_rows[:len(unique_rows)], k=30)
def normalsestimation(pointcloud,nn_glob,VP=[0,0,0]):
ViewPoint = np.array(VP)
normals = np.empty((np.shape(pointcloud)))
curv = np.empty((len(pointcloud),1))
for index in range(len(pointcloud)):
nn_loc = pointcloud[nn_glob[index]]
COV = np.cov(nn_loc,rowvar=False)
eigval, eigvec = np.linalg.eig(COV)
idx = np.argsort(eigval)
nor = eigvec[:,idx][:,0]
if nor.dot((ViewPoint-pointcloud[index,:])) > 0:
normals[index] = nor
else:
normals[index] = - nor
curv[index] = eigval[idx][0] / np.sum(eigval)
return normals,curv
#seed_count = 0
#while seed_count < len(current_seeds)
def regiongrowing(pointcloud,nn_glob,theta_th = 'auto', cur_th = 'auto'):
normals,curvature = normalsestimation(pointcloud,nn_glob=nn_glob)
order = curvature[:,0].argsort().tolist()
region = []
if theta_th == 'auto':
theta_th = 15.0 / 180.0 * math.pi # in radians
if cur_th == 'auto':
cur_th = np.percentile(curvature,98)
while len(order) > 0:
region_cur = []
seed_cur = []
poi_min = order[0] #poi_order[0]
region_cur.append(poi_min)
seed_cur.append(poi_min)
order.remove(poi_min)
for i in range(len(seed_cur)):
nn_loc = nn_glob[seed_cur[i]]
for j in range(len(nn_loc)):
nn_cur = nn_loc[j]
if all([nn_cur in order , np.arccos(np.abs(np.dot(normals[seed_cur[i]],normals[nn_cur])))<theta_th]):
region_cur.append(nn_cur)
order.remove(nn_cur)
if curvature[nn_cur] < cur_th:
seed_cur.append(nn_cur)
region.append(region_cur)
return region
def regiongrowing1(pointcloud,nn_glob,theta_th = 'auto', cur_th = 'auto'):
normals,curvature = normalsestimation(pointcloud,nn_glob=nn_glob)
order = curvature[:,0].argsort().tolist()
region = []
if theta_th == 'auto':
theta_th = 15.0 / 180.0 * math.pi # in radians
if cur_th == 'auto':
cur_th = np.percentile(curvature,98)
while len(order) > 0:
region_cur = []
seed_cur = []
poi_min = order[0] #poi_order[0]
region_cur.append(poi_min)
seedval = 0
seed_cur.append(poi_min)
order.remove(poi_min)
# for i in range(len(seed_cur)):#change to while loop
while seedval < len(seed_cur):
nn_loc = nn_glob[seed_cur[seedval]]
for j in range(len(nn_loc)):
nn_cur = nn_loc[j]
if all([nn_cur in order , np.arccos(np.abs(np.dot(normals[seed_cur[seedval]],normals[nn_cur])))<theta_th]):
region_cur.append(nn_cur)
order.remove(nn_cur)
if curvature[nn_cur] < cur_th:
seed_cur.append(nn_cur)
seedval+=1
region.append(region_cur)
return region
region= regiongrowing(unique_rows,nn_glob)
region1 = regiongrowing1(unique_rows,nn_glob)
尝试一下，看看你得到了什么。
我正在尝试 运行 像您的代码一样使用 python 为点云增长的区域，但它有一个错误。
C:\Users\youss\AppData\Local\Temp\ipykernel_1300869360685.py:23: RuntimeWarning: arccos 中遇到无效值
如果全部（[nn_cur 按顺序排列，np.arccos（np.abs（np.dot（法线[seed_cur[seedval]]，法线[nn_cur] )))
我正在 python 中针对非结构化点云开发区域增长算法。我尝试根据 http://pointclouds.org/documentation/tutorials/region_growing_segmentation.php#region-growing-segmentation 和 SEGMENTATION OF POINT CLOUDS USING SMOOTHNESS CONSTRAINT.
实现算法
使用相同的点云，我无法使用 pcl 实现重现结果，尽管所有步骤和参数似乎都相同。我的实现产生了 217 个区域，而 pcl 的实现产生了 24 个（过度分割）。可以找到点云：
https://github.com/scrox/region-growing/blob/master/test.txt
是什么导致了这种显着差异？我也欢迎我的 python 代码的性能改进。
import math
import numpy as np
from sklearn.neighbors import KDTree
unique_rows=np.loadtxt("test.txt")
tree = KDTree(unique_rows, leaf_size=2)
dist,nn_glob = tree.query(unique_rows[:len(unique_rows)], k=30)
def normalsestimation(pointcloud,nn_glob,VP=[0,0,0]):
ViewPoint = np.array(VP)
normals = np.empty((np.shape(pointcloud)))
curv = np.empty((len(pointcloud),1))
for index in range(len(pointcloud)):
nn_loc = pointcloud[nn_glob[index]]
COV = np.cov(nn_loc,rowvar=False)
eigval, eigvec = np.linalg.eig(COV)
idx = np.argsort(eigval)
nor = eigvec[:,idx][:,0]
if nor.dot((ViewPoint-pointcloud[index,:])) > 0:
normals[index] = nor
else:
normals[index] = - nor
curv[index] = eigval[idx][0] / np.sum(eigval)
return normals,curv
def regiongrowing(pointcloud,nn_glob,theta_th = 'auto', cur_th = 'auto'):
normals,curvature = normalsestimation(pointcloud,nn_glob=nn_glob)
order = curvature[:,0].argsort().tolist()
region = []
if theta_th == 'auto':
theta_th = 15.0 / 180.0 * math.pi # in radians
if cur_th == 'auto':
cur_th = np.percentile(curvature,98)
while len(order) > 0:
region_cur = []
seed_cur = []
poi_min = order[0] #poi_order[0]
region_cur.append(poi_min)
seed_cur.append(poi_min)
order.remove(poi_min)
for i in range(len(seed_cur)):
nn_loc = nn_glob[seed_cur[i]]
for j in range(len(nn_loc)):
nn_cur = nn_loc[j]
if all([nn_cur in order , np.arccos(np.abs(np.dot(normals[seed_cur[i]],normals[nn_cur])))<theta_th]):
region_cur.append(nn_cur)
order.remove(nn_cur)
if curvature[nn_cur] < cur_th:
seed_cur.append(nn_cur)
region.append(region_cur)
return region
region = regiongrowing(unique_rows,nn_glob)
这里是 pcl 代码：
#include <iostream>
#include <vector>
#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/search/search.h>
#include <pcl/search/kdtree.h>
#include <pcl/features/normal_3d.h>
#include <pcl/visualization/cloud_viewer.h>
#include <pcl/segmentation/region_growing.h>
int
main (int argc, char** argv)
{
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
if ( pcl::io::loadPCDFile <pcl::PointXYZ> ("test.pcd", *cloud) == -1)
{
std::cout << "Cloud reading failed." << std::endl;
return (-1);
}
pcl::search::Search<pcl::PointXYZ>::Ptr tree = boost::shared_ptr<pcl::search::Search<pcl::PointXYZ> > (new pcl::search::KdTree<pcl::PointXYZ>);
pcl::PointCloud <pcl::Normal>::Ptr normals (new pcl::PointCloud <pcl::Normal>);
pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> normal_estimator;
normal_estimator.setSearchMethod (tree);
normal_estimator.setInputCloud (cloud);
normal_estimator.setKSearch (30);
normal_estimator.compute (*normals);
pcl::RegionGrowing<pcl::PointXYZ, pcl::Normal> reg;
reg.setSearchMethod (tree);
reg.setNumberOfNeighbours (30);
reg.setInputCloud (cloud);
reg.setInputNormals (normals);
reg.setSmoothnessThreshold (15.0 / 180.0 * M_PI);
reg.setCurvatureThreshold (0.088614);
std::vector <pcl::PointIndices> clusters;
reg.extract (clusters);
std::cout << "Number of clusters is equal to " << clusters.size () << std::endl;
std::cout << "First cluster has " << clusters[0].indices.size () << " points." << endl;
std::cout << "These are the indices of the points of the initial" <<
std::endl << "cloud that belong to the first cluster:" << std::endl;
int counter = 0;
while (counter < clusters[0].indices.size ())
{
std::cout << clusters[0].indices[counter] << ", ";
counter++;
}
std::cout << std::endl;
pcl::PointCloud <pcl::PointXYZRGB>::Ptr colored_cloud = reg.getColoredCloud ();
pcl::visualization::CloudViewer viewer ("Cluster viewer");
viewer.showCloud(colored_cloud);
while (!viewer.wasStopped ())
{
}
return (0);
}
我认为这与您 def regiongrowing(pointcloud,nn_glob,theta_th = 'auto', cur_th = 'auto'): 中的 for i in range(len(seed_cur)): 有关。我相信当 python 使用 for 循环时，它实际上并没有遍历循环开始后添加到列表中的任何元素。所以
if curvature[nn_cur] < cur_th:
seed_cur.append(nn_cur)
它不会继续循环 nn_cur。我使用 while 循环稍微修改了您的代码。我的结果是 21 个区域。与 PCL 中的情况不完全相同，但更接近。我的函数是regiongrowing1，供参考。
import math
import numpy as np
from sklearn.neighbors import KDTree
unique_rows=np.loadtxt("test.txt")
tree = KDTree(unique_rows, leaf_size=2)
dist,nn_glob = tree.query(unique_rows[:len(unique_rows)], k=30)
def normalsestimation(pointcloud,nn_glob,VP=[0,0,0]):
ViewPoint = np.array(VP)
normals = np.empty((np.shape(pointcloud)))
curv = np.empty((len(pointcloud),1))
for index in range(len(pointcloud)):
nn_loc = pointcloud[nn_glob[index]]
COV = np.cov(nn_loc,rowvar=False)
eigval, eigvec = np.linalg.eig(COV)
idx = np.argsort(eigval)
nor = eigvec[:,idx][:,0]
if nor.dot((ViewPoint-pointcloud[index,:])) > 0:
normals[index] = nor
else:
normals[index] = - nor
curv[index] = eigval[idx][0] / np.sum(eigval)
return normals,curv
#seed_count = 0
#while seed_count < len(current_seeds)
def regiongrowing(pointcloud,nn_glob,theta_th = 'auto', cur_th = 'auto'):
normals,curvature = normalsestimation(pointcloud,nn_glob=nn_glob)
order = curvature[:,0].argsort().tolist()
region = []
if theta_th == 'auto':
theta_th = 15.0 / 180.0 * math.pi # in radians
if cur_th == 'auto':
cur_th = np.percentile(curvature,98)
while len(order) > 0:
region_cur = []
seed_cur = []
poi_min = order[0] #poi_order[0]
region_cur.append(poi_min)
seed_cur.append(poi_min)
order.remove(poi_min)
for i in range(len(seed_cur)):
nn_loc = nn_glob[seed_cur[i]]
for j in range(len(nn_loc)):
nn_cur = nn_loc[j]
if all([nn_cur in order , np.arccos(np.abs(np.dot(normals[seed_cur[i]],normals[nn_cur])))<theta_th]):
region_cur.append(nn_cur)
order.remove(nn_cur)
if curvature[nn_cur] < cur_th:
seed_cur.append(nn_cur)
region.append(region_cur)
return region
def regiongrowing1(pointcloud,nn_glob,theta_th = 'auto', cur_th = 'auto'):
normals,curvature = normalsestimation(pointcloud,nn_glob=nn_glob)
order = curvature[:,0].argsort().tolist()
region = []
if theta_th == 'auto':
theta_th = 15.0 / 180.0 * math.pi # in radians
if cur_th == 'auto':
cur_th = np.percentile(curvature,98)
while len(order) > 0:
region_cur = []
seed_cur = []
poi_min = order[0] #poi_order[0]
region_cur.append(poi_min)
seedval = 0
seed_cur.append(poi_min)
order.remove(poi_min)
# for i in range(len(seed_cur)):#change to while loop
while seedval < len(seed_cur):
nn_loc = nn_glob[seed_cur[seedval]]
for j in range(len(nn_loc)):
nn_cur = nn_loc[j]
if all([nn_cur in order , np.arccos(np.abs(np.dot(normals[seed_cur[seedval]],normals[nn_cur])))<theta_th]):
region_cur.append(nn_cur)
order.remove(nn_cur)
if curvature[nn_cur] < cur_th:
seed_cur.append(nn_cur)
seedval+=1
region.append(region_cur)
return region
region= regiongrowing(unique_rows,nn_glob)
region1 = regiongrowing1(unique_rows,nn_glob)
尝试一下，看看你得到了什么。
我正在尝试 运行 像您的代码一样使用 python 为点云增长的区域，但它有一个错误。
C:\Users\youss\AppData\Local\Temp\ipykernel_1300869360685.py:23: RuntimeWarning: arccos 中遇到无效值
如果全部（[nn_cur 按顺序排列，np.arccos（np.abs（np.dot（法线[seed_cur[seedval]]，法线[nn_cur] )))