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

  • 不想做只会调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/

原文地址：https://www.cnblogs.com/xiaoniubenrecord-6161/p/12186576.html

时间： 2024-10-11 20:39:45

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

PCL源码剖析之MarchingCubes算法

原文:http://blog.csdn.net/lming_08/article/details/19432877 MarchingCubes算法简介 MarchingCubes(移动立方体)算法是目前三围数据场等值面生成中最常用的方法.它实际上是一个分而治之的方法,把等值面的抽取分布于每个体素中进行.对于每个被处理的体素,以三角面片逼近其内部的等值面片.每个体素是一个小立方体,构造三角面片的处理过程对每个体素都“扫描”一遍,就好像一个处理器在这些体素上移动一样,由此得名移动立方体算法. MC算

OpenCV学习笔记（27）KAZE 算法原理与源码分析（一）非线性扩散滤波

http://blog.csdn.net/chenyusiyuan/article/details/8710462 OpenCV学习笔记(27)KAZE 算法原理与源码分析(一)非线性扩散滤波 2013-03-23 17:44 16963人阅读 评论(28) 收藏 举报 分类: 机器视觉(34) 版权声明:本文为博主原创文章,未经博主允许不得转载. 目录(?)[+] KAZE系列笔记: OpenCV学习笔记(27)KAZE 算法原理与源码分析(一)非线性扩散滤波 OpenCV学习笔记(28)KA

iOS下使用SHA1WithRSA算法加签源码

首先了解一下几个相关概念,以方便后面遇到的问题的解决: RSA算法:1977年由Ron Rivest.Adi Shamirh和LenAdleman发明的,RSA就是取自他们三个人的名字.算法基于一个数论:将两个大素数相乘非常容易,但要对这个乘积的结果进行因式分解却非常困难,因此可以把乘积公开作为公钥.该算法能够抵抗目前已知的所有密码攻击.RSA算法是一种非对称算法,算法需要一对密钥,使用其中一个加密,需要使用另外一个才能解密.我们在进行RSA加密通讯时,就把公钥放在客户端,私钥留在服务器. PE

基于Linux平台下网络病毒Caem.c源码及解析

Came.c型病毒在这里主要修改了用户的密码,同时对用户的终端设备进行了监视.希望与大家共同交流 转载请注明出处:http://blog.csdn.net/u010484477     O(∩_∩)O谢谢 #define HOME "/" #define TIOCSCTTY 0x540E #define TIOCGWINSZ 0x5413 #define TIOCSWINSZ 0x5414 #define ECHAR 0x1d #define PORT 39617 #define BU

基于Windows Socket的安全通信(C++实现,附源码)

先了解一下Socket的相关函数原型 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 //加载套接字库 int PASCAL FAR WSAStartup(WORD wVersionRequired, LPWSADATA lpWSAData); //释放套接字库资源 int PASCAL FAR WSACleanup(void); //创建套接字 SOCKET PASCAL FAR socket (int af,int type,int pr

Spark技术内幕：Master基于ZooKeeper的High Availability（HA）源码实现

如果Spark的部署方式选择Standalone,一个采用Master/Slaves的典型架构,那么Master是有SPOF(单点故障,Single Point of Failure).Spark可以选用ZooKeeper来实现HA. ZooKeeper提供了一个Leader Election机制,利用这个机制可以保证虽然集群存在多个Master但是只有一个是Active的,其他的都是Standby,当Active的Master出现故障时,另外的一个Standby Master会被选举出来.由于

【开源下载】基于winform的xml菜单编辑器（c#源码)

xml编辑器源码 最近帮朋友做了一个档案管理系统,客户端能够把文件上传到服务器,也能够从服务器下载,支持多用户.通讯框架使用的networkcomms v3框架. 这个档案管理系统中用到了树形目录,使用人员需要随时调整左侧的目录,考虑到使用数据库的比较繁琐,就想到了一个方法,即可以在客户端编辑左侧的那个目录,保存成一个xml文件.修改完成后需要的话可以把这个xml文件上传到服务器,其他人员可以从服务器加载这个xml文件.虽然简单,但也比较好的满足了朋友的需求.今天刚好有时间,把左侧目录的编辑页面

基于数据波动性的分割算法

我们常见的分割算法有很多种,比如能量法,包络线法之类的,但这些算法难以实现实时分割,今天我给大家分享一个原创的分割算法,是在以前项目中用过的,这两天加以优化,最中整理了一个MATLAB版本的,给大家分享一下. 算法的原理简单介绍一下: 这里给出了一段肌音信号(已经分割好了),是用加速度传感器在手上采集的,每次完成一次动作,就会产生一个数据波动,如果我们需要分析这样一段信号的特征,需要先将这些信号分割出来.   我们在分析时,主要任务时提取出信号帧起始点对应的数据索引,这里我写了一个分割函数: f

机器学习IB1算法的weka源码详细解析（1NN）

机器学习的1NN最近邻算法,在weka里叫IB1,是因为Instance Base  1 ,也就是只基于一个最近邻的实例的惰性学习算法. 下面总结一下,weka中对IB1源码的学习总结. 首先需要把 weka-src.jar 引入编译路径,否则无法跟踪源码. 1)读取data数据,完成 IB1 分类器的调用,结果预测评估.为了后面的跟踪. try { File file = new File("F:\\tools/lib/data/contact-lenses.arff"); Arff