基于区域增长(RegionGrowing)的分割算法——参照pcl源码

Posted 2020-11-24 xiaoniubenrecord-6161

• 不想做只会调API的程序员，进而重写了pcl::RegionGrowing类（然而还是基于了pcl的数据结构，哎，学习有点迷茫）。
• 分割，顾名思义，按照一定的规律将点云分成一类类的。方便于接下来对点云的操作处理。不同的应用方向会用到不同的分割方法。本篇介绍的基于区域增长的算法，最终达到的理想效果是 点云按照估计的曲率聚类，但本人做了一些小的demo示例之后其实对于实际的应用还是一头雾水（也可能是比较菜吧 /emoji sad/）然而该算法的思想却是可以给做别的提供一些思路。
• 算法思想：首先根据点的曲率进行由小到大的排序，然后将曲率最小的点设为初始种子点，这样算法就会从场景中最平滑的区域开始，可以减少分割数。然后对于满足阈值（曲率和相邻法线夹角）的点继续加入seed队列来判断，直到一个seed队列全部判断清空，这意味着一块区域分割完毕。按照曲率顺序选择未判断的点开始新seed队列，重复以上步骤，直到最后一个点判度完毕。

 

#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <cmath>

#include <pcl-1.8/pcl/common/common.h>
#include <pcl-1.8/pcl/console/print.h>
#include <pcl-1.8/pcl/search/kdtree.h>
#include <pcl-1.8/pcl/PCLPointCloud2.h>
#include <pcl-1.8/pcl/point_types.h>
#include <pcl-1.8/pcl/PointIndices.h>
#include <pcl-1.8/pcl/features/normal_3d.h>
#include <pcl-1.8/pcl/pcl_exports.h>
#include <pcl-1.8/pcl/io/pcd_io.h>
#include <pcl-1.8/pcl/visualization/pcl_visualizer.h>

#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/Geometry>

using namespace std;
typedef pcl::PointXYZ PointT;

int main(int argc,char** argv)
{
    pcl::PointCloud<PointT>::Ptr cloud(new pcl::PointCloud<PointT>());
    pcl::io::loadPCDFile(argv[1],*cloud);
    //pcl::io::loadPLYFile(argv[1],*cloud); //load ply format;

    //normal estimation;
    pcl::search::KdTree<PointT>::Ptr tree(new pcl::search::KdTree<PointT>());
    pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>());
    pcl::NormalEstimation<PointT,pcl::Normal> ne;
    tree->setInputCloud(cloud);
    ne.setInputCloud(cloud);
    ne.setSearchMethod(tree);
    ne.setKSearch(20);
    ne.compute(*normals);

    int k = 20;
    int points_num = cloud->points.size();
    vector<int> k_nebor_index;
    vector<float> k_nebor_index_dis;
    vector<vector<int>> point_k_index;
    point_k_index.resize(points_num,k_nebor_index);
    for (size_t i = 0; i < points_num; i++)
    {
       if (tree->nearestKSearch(cloud->points[i],k,k_nebor_index,k_nebor_index_dis))
       {
          point_k_index[i].swap(k_nebor_index);
       }
       else
       {
           PCL_ERROR("WARNING:FALSE NEARERTKSEARCH");
       }
    }

    vector<pair<float,int>> vec_curvature;
    for (size_t i = 0; i < points_num; i++)
    {
        vec_curvature.push_back(make_pair(normals->points[i].curvature,i));
    }
    sort(vec_curvature.begin(),vec_curvature.end());
    float curvature_threshold = 0.1;
    float normal_threshold = cosf(10.0/180.0*M_PI);
    int seed_orginal = vec_curvature[0].second;
    int counter_0 = 0;
    int segment_laber(0);
    vector<int> segmen_num;
    vector<int> point_laber;
    point_laber.resize(points_num,-1);
    while(counter_0 < points_num)
    {   
        queue<int> seed;
        seed.push(seed_orginal);
        point_laber[seed_orginal] = segment_laber;
        int counter_1(1);
        while(!seed.empty())
        {
            int curr_seed = seed.front();
            seed.pop();
            int curr_nebor_index(0);
            while (curr_nebor_index<k)
            {
                bool is_a_seed = false;
                int cur_point_idx = point_k_index[curr_seed][curr_nebor_index];
                if (point_laber[cur_point_idx] != -1)
                {
                    curr_nebor_index++;
                    continue;
                }
                Eigen::Map<Eigen::Vector3f> vec_curr_point(static_cast<float*>(normals->points[curr_seed].normal));
                Eigen::Map<Eigen::Vector3f> vec_nebor_point (static_cast<float*>(normals->points[cur_point_idx].normal));
                float dot_normals = fabsf(vec_curr_point.dot(vec_nebor_point));
                if(dot_normals<normal_threshold)
                {
                    is_a_seed = false;
                }
                else if(normals->points[cur_point_idx].curvature>curvature_threshold)
                {
                    is_a_seed = false;
                }
                else
                {
                    is_a_seed = true;
                }                
                if (!is_a_seed)
                {
                    curr_nebor_index++;
                    continue;
                }
                point_laber[cur_point_idx] = segment_laber;
                counter_1++;
                if (is_a_seed)
                {
                    seed.push(cur_point_idx);
                }
                curr_nebor_index++;  
            }      
        }
        segment_laber++; 
        counter_0 +=counter_1;
        segmen_num.push_back(counter_1);
        for (size_t i = 0; i < points_num; i++)
        {
            int index_curvature = vec_curvature[i].second;
            if(point_laber[index_curvature] == -1)
            {
                seed_orginal = index_curvature;
                break;
            }           
        }
    }
    cout<<"seg_num:"<<segmen_num.size()<<endl;
    //summary of segmentation results
    vector<pcl::PointIndices> cluster;
    pcl::PointIndices segment_points;
    int seg_num = segmen_num.size();
    cluster.resize(seg_num,segment_points);
    for(size_t i_seg = 0;i_seg<seg_num;i_seg++)
    {
        cluster[i_seg].indices.resize(segmen_num[i_seg],0);
    }
    vector<int> counter_2;
    counter_2.resize(seg_num,0);
    for(size_t i_point =0 ;i_point<points_num;i_point++)
    {
        int seg_idx = point_laber[i_point];
        int nebor_idx = counter_2[seg_idx];
        cluster[seg_idx].indices[nebor_idx] = i_point;
        counter_2[seg_idx] +=1;
    }

    //Remove outline points
    vector<pcl::PointIndices> clusters;
    int minNum = 100; 
    int maxNum = 100000;
    for(size_t i_cluster = 0;i_cluster<seg_num;i_cluster++)
    {
        if(cluster[i_cluster].indices.size() < maxNum && cluster[i_cluster].indices.size() > minNum)
        {
            clusters.push_back(cluster[i_cluster]);
        }
    }

    // visualization
    boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer(new pcl::visualization::PCLVisualizer("RegionGrowing Viewer"));
    pcl::PointCloud<pcl::PointXYZRGB>::Ptr cloud_color(new pcl::PointCloud<pcl::PointXYZRGB>());
    srand(time(nullptr));
    vector<unsigned char> color;
    for (size_t i = 0; i < clusters.size(); i++)
    {
        color.push_back(static_cast<unsigned char>(rand()%256));
        color.push_back(static_cast<unsigned char>(rand()%256));
        color.push_back(static_cast<unsigned char>(rand()%256));
    }
    int color_index(0);
    for(size_t i = 0;i<clusters.size();i++)
    {
        for(size_t j = 0; j<clusters[i].indices.size();j++)
        {
            pcl::PointXYZRGB n_points;
            n_points.x = cloud->points[clusters[i].indices[j]].x;
            n_points.y = cloud->points[clusters[i].indices[j]].y;
            n_points.z = cloud->points[clusters[i].indices[j]].z;
            n_points.r = color[3*color_index];
            n_points.g = color[3*color_index+1];
            n_points.b = color[3*color_index+2];
            cloud_color->push_back(n_points);
        }
        color_index++;
    }
    viewer->addPointCloud(cloud_color);
    viewer->spin();
    return 0;
}

CMakeLists.txt:

cmake_minimum_required (VERSION 2.8)
project(my_regionGrowing)

set(CMAKE_BUILD_TYPE "Release")
find_package(PCL 1.8 REQUIRED)
include_directories(my_regionGrowing ${PCL_INCLUDE_DIRS})
link_libraries(${PCL_LIBRARY_DIRS})
add_definitions(${PCL_DEFINITIONS})

add_executable(RegionGrowing main.cpp)
target_link_libraries(RegionGrowing ${PCL_LIBRARIES})

示例demo:

hand.pcd

颜色部分为随即赋值，代码部分还进行了点的去除，我们将聚类中小于100 和 大于 100000的点去除了，还有阈值部分，也调了几次参数，但是对于阈值的设置还是很模糊的。（/emoji sad/) 然而这样实现一遍还是比调API要开心很多，因为本身代码能力也比较差，也不知道关于此算法有没有复杂度更低，更优美简洁的实现方法了。对于我的学习方法是否正确有效，本人也比较迷茫。希望看到这篇文章的大佬能够给予些指导，也希望能够和大家一起进步。对于文中披漏部分希望大家能够指正。#参照PCL-1.8/pcl/segmentation/impl/region_growing.hpp/