1.介绍
 
本项目使用ubuntu16.04，ubuntu的其他版本应该也是可以的；
 
点云库PCL某某这本书其实是以官方的英文手册翻译成的，英文链接见以下地址：
 
https://pcl.readthedocs.io/projects/tutorials/en/latest/region_growing_segmentation.html#region-growing-segmentation
 
但是想做相关的算法开发，只看书中的内容还是不够明白，只能找来官方文档读了，以下基本是翻译再加一些我自己的理解，如果有不清楚的还是建议自己去啃英文说明。
 
使用pcl:: regiong类实现区域增长算法。该算法的目的是在光滑性约束条件下合并相近的点，实现区域的分割，适用于不同平/曲面的分割。该算法的输出是不同曲面单元的集合，其中每个单元可以认为是一个光滑表面的点的总成。该算法基于点之间的法线角度比较完成划分。
 
2.算法实现
 
因为生长点是根据曲率进行寻找的，选择曲率最小的点，所以首先根据曲率将点进行分类。由于选择的曲率是最小的点，所以生长点附近一般均为平面或者相对较平坦的面。
 
将生长点添加进种子的集合。然后对种子里的每个生长点进行生长运算。
 
种子点的运算过程如下：
 
测试每个临近点与种子点法向量的夹角，如果角度小于临界值，则添加进当前的种子区域；
测试每个临近点的曲率，如果曲率小于临界值，则将该点添加进当前区域的种子列表；
第三步的意思其实是在找到曲率小于临界值的点后添加进种子列表，原来的种子点就要删除了，但是这一个地方有一个问题，是否需要原有的种子周围临近点全部搜索完毕，这个应该是要搜索完毕的；但是如果发现多个曲率小于临界值的点，是否需要都添加？我的理解是都需要添加；
注：曲面的一阶微分的分量包括法向量，曲面的二阶微分的分量包括曲率；从某种意义上说法向量的夹角是类似与曲率的。但是为什么区域增长算法通过夹角生成区域，而通过曲率生成种子列表？这个后面会有描述。
 
待种子列表中为空时，说明该区域生长完毕；
 
这里贴出算法伪代码：
 
 
 
 
 
 
 
关于伪代码这里就不解释了，只要花时间认真去看，相信是可以理解的。其中的\是去除的意思。
 
代码如下：
 
#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/filters/passthrough.h>
#include <pcl/segmentation/region_growing.h>

int
main (int argc, char** argv)
{
//加载pcd文件
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
  if ( pcl::io::loadPCDFile <pcl::PointXYZ> ("region_growing_tutorial.pcd", *cloud) == -1)
  {
    std::cout << "Cloud reading failed." << std::endl;
    return (-1);
  }
//通过NormalEstimation类计算法线
  pcl::search::Search<pcl::PointXYZ>::Ptr tree (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 (50);
  normal_estimator.compute (*normals);
//在z轴方向滤波，剔除0~1.0范围之外的点
// PassThrough是滤波类，得到的indices是各点的索引值
  pcl::IndicesPtr indices (new std::vector <int>);
  pcl::PassThrough<pcl::PointXYZ> pass;
  pass.setInputCloud (cloud);
  pass.setFilterFieldName ("z");
  pass.setFilterLimits (0.0, 1.0);
  pass.filter (*indices);
//设置一个簇中包含点的最大值门限与最小值门限，通过区域生长算法生成的所有的簇，小于50个点或者大//于1000000个点都被丢弃，不包含在簇内
  pcl::RegionGrowing<pcl::PointXYZ, pcl::Normal> reg;
  reg.setMinClusterSize (50);
  reg.setMaxClusterSize (1000000);
//算法需要k近邻搜索，搜索近邻点的数量限制为30
  reg.setSearchMethod (tree);
  reg.setNumberOfNeighbours (30);
  reg.setInputCloud (cloud);
  //reg.setIndices (indices);
  reg.setInputNormals (normals);
//此处两个初始化的值是最关键的，直接影响分割结果的好坏；
// setSmoothnessThreshold是法线的夹角门限
// setCurvatureThreshold是曲率门限，对于满足发现夹角的点，如果也满足曲率门限
//则将会被添加进种子列表内；
  reg.setSmoothnessThreshold (3.0 / 180.0 * M_PI);
  reg.setCurvatureThreshold (1.0);
//下面执行区域生长算法，返回一个簇的数组
  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." << std::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++;
    if (counter % 10 == 0)
      std::cout << std::endl;
  }
  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);
}
 
需要将region_growing_tutorial.pcd替换成你想分析的pcd文件，或者自己去网站下载相应的pcd文件；
 
3.编译与执行
 
CMakeList.txt内容如下：
 
cmake_minimum_required(VERSION 2.8 FATAL_ERROR)

project(region_growing_segmentation)

find_package(PCL 1.5 REQUIRED)

include_directories(${PCL_INCLUDE_DIRS})
link_directories(${PCL_LIBRARY_DIRS})
add_definitions(${PCL_DEFINITIONS})

add_executable (region_growing_segmentation region_growing_segmentation.cpp)
target_link_libraries (region_growing_segmentation ${PCL_LIBRARIES})
 
将CMakeList.txt与cpp文件放在同一个文件夹，然后依次执行
 
cmake CMakeList.txt
 
make
 
./region_growing_segmentation
 
执行之后就显示出了分割后的图，不同区域通过不同颜色分割，其中的红色点是为分区的点，
 
  reg.setMinClusterSize (50);
 
  reg.setMaxClusterSize (1000000);
 
红色点区域上述规划的小于50或者大于1000000的情况
 
4.问题总结
 
关于依据曲率选择种子的问题
 
根据官方的文档可以看出，是在通过法线夹角筛选后的点中选择满足曲率要求的种子点的。同时，种子点也是加进种子列表的，其原来的点并没有删除，只有当原来的点的规定半径区域的点运算完毕后才会将其删除，这一点官方文档描述的不是很清楚，但是伪代码写的是对的。
 
关于初始种子的选取问题
 
另一篇文章有描述，是将所有点按照曲率值进行排序，找到最小曲率点，将其添加到种子点集，链接如下：
 
https://cloud.tencent.com/developer/article/1475904
 
 

 
 
原文链接：Region growing segmentation
 
 
本教程使用到的点云数据：source files
 
在本篇教程中，我们将学习如何使用由pcl::RegionGrowing类实现的区域生长算法。该算法的目的是合并在平滑约束条件下足够接近的点。因此，该算法的输出数据结构是由聚类组成的数组，其中每个聚类都是被认为是同一光滑表面的一部分的点的集合。该算法的工作原理（光滑度的计算）是基于两点法线之间的角度比较。
 
目录
 
基本原理
 
算法伪代码
 
输入
 
初始化
 
算法实现
 
程序代码
 
代码分析
 
实验结果
 
 其他数据结果
 

基本原理
 
首先，它根据点的曲率值对点进行排序。需要这样做是因为区域从曲率最小的点开始增长。这样做的原因是曲率最小的点位于平坦区域(从最平坦区域生长可以减少段的总数)。
 
我们有了分类过的云。直到云中没有未标记点时，算法选取曲率值最小的点，开始区域的增长。这个过程如下所示：
 
选中的点被添加到名为种子的集合中。
对于每一个种子点，找到它的邻近点： 
  
算出每个相邻点的法线和当前种子点的法线之间的角度，如果角度小于阈值，则将当前点添加到当前区域。
然后计算每个邻居点的曲率值，如果曲率小于阈值，那么这个点被添加到种子中。
将当前的种子从种子列表中移除。
如果种子列表变成空的，这意味着该区域生长已完成，继续重复上述过程。
 
 
 
区域生长算法： 将具有相似性的点云集合起来构成区域。
 首先对每个需要分割的区域找出一个种子点作为生长的起点，然后将种子点周围邻域中与种子有相同或 相似性质的点合并到种子像素所在的区域中。而新的点继续作为种子向四周生长，直到再没有满足条件 的像素可以包括进来，一个区域就生长而成了。
 算法流程：
 
 
计算法线normal 和曲率curvatures，依据曲率升序 排序；
选择曲率最低的为初始种子点，种子周围的临近点和 种子点云相比较；
法线的方向是否足够相近（法线夹角足够 r p y）， 法线夹角阈值；
曲率是否足够小（表面处在同一个弯曲程度），区域 差值阈值；
如果满足2，3则该点可用做种子点; 6. 如果只满足2，则归类而不做种子;
 
算法伪代码
 
输入
 
Point cloud = 
Point normals = 
Points curvatures = 
Neighbour finding function 
Curvature threshold 
Angle threshold 
初始化
 
区域列表置为空： Region list 
 
可用的点云列表：Available points list 
 
算法实现
 
 
程序代码
 
#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/filters/passthrough.h>
#include <pcl/segmentation/region_growing.h>
#include <pcl/console/print.h>
#include <pcl/console/parse.h>
#include <pcl/console/time.h>
#include <windows.h>
#include <stdio.h>
#include <psapi.h>
void PrintMemoryInfo( )
{
	HANDLE hProcess;
	PROCESS_MEMORY_COUNTERS pmc;

	hProcess=GetCurrentProcess();
	printf( "\nProcess ID: %u\n", hProcess );

	// Print information about the memory usage of the process.
	//输出进程使用的内存信息
   
	if (NULL == hProcess)
		return;

	if ( GetProcessMemoryInfo( hProcess, &pmc, sizeof(pmc)) )
	{
		printf( "\tPageFaultCount: 0x%08X\n", pmc.PageFaultCount );
		printf( "\tPeakWorkingSetSize: 0x%08X\n", 
				  pmc.PeakWorkingSetSize );
		printf( "\tWorkingSetSize: 0x%08X\n", pmc.WorkingSetSize );
		printf( "\tQuotaPeakPagedPoolUsage: 0x%08X\n", 
				  pmc.QuotaPeakPagedPoolUsage );
		printf( "\tQuotaPagedPoolUsage: 0x%08X\n", 
				  pmc.QuotaPagedPoolUsage );
		printf( "\tQuotaPeakNonPagedPoolUsage: 0x%08X\n", 
				  pmc.QuotaPeakNonPagedPoolUsage );
		printf( "\tQuotaNonPagedPoolUsage: 0x%08X\n", 
				  pmc.QuotaNonPagedPoolUsage );
		printf( "\tPagefileUsage: 0x%08X\n", pmc.PagefileUsage ); 
		printf( "\tPeakPagefileUsage: 0x%08X\n", 
				  pmc.PeakPagefileUsage );
	}

	CloseHandle( hProcess );
}

using namespace pcl::console;
int
main (int argc, char** argv)
{

	if(argc<2)
	{
		std::cout<<".exe xx.pcd -kn 50 -bc 0 -fc 10.0 -nc 0 -st 30 -ct 0.05"<<endl;

		return 0;
	}//如果输入参数小于1个，输出提示
	time_t start,end,diff[5],option;
	start = time(0); 
	int KN_normal=50; //设置默认输入参数
	bool Bool_Cuting=false;//设置默认输入参数
	float far_cuting=10,near_cuting=0,SmoothnessThreshold=30.0,CurvatureThreshold=0.05;//设置默认输入参数
	parse_argument (argc, argv, "-kn", KN_normal);
	parse_argument (argc, argv, "-bc", Bool_Cuting);
	parse_argument (argc, argv, "-fc", far_cuting);
	parse_argument (argc, argv, "-nc", near_cuting);
	parse_argument (argc, argv, "-st", SmoothnessThreshold);
	parse_argument (argc, argv, "-ct", CurvatureThreshold);//设置输入参数方式
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
	if ( pcl::io::loadPCDFile <pcl::PointXYZ> (argv[1], *cloud) == -1)
	{
		std::cout << "Cloud reading failed." << std::endl;
		return (-1);
	}// 加载输入点云数据
	end = time(0); 
	diff[0] = difftime (end, start); 
	PCL_INFO ("\Loading pcd file takes(seconds): %d\n", diff[0]); 
	//Noraml estimation step(1 parameter)
	pcl::search::Search<pcl::PointXYZ>::Ptr tree = boost::shared_ptr<pcl::search::Search<pcl::PointXYZ> > (new pcl::search::KdTree<pcl::PointXYZ>);//创建一个指向kd树搜索对象的共享指针
	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 (KN_normal);// 设置用于法向量估计的k近邻数目
	normal_estimator.compute (*normals);//计算并输出法向量
	end = time(0); 
	diff[1] = difftime (end, start)-diff[0]; 
	PCL_INFO ("\Estimating normal takes(seconds): %d\n", diff[1]); 
	//optional step: cutting the part are far from the orignal point in Z direction.2 parameters
	pcl::IndicesPtr indices (new std::vector <int>);//创建一组索引
	if(Bool_Cuting)//判断是否需要直通滤波
	{

		pcl::PassThrough<pcl::PointXYZ> pass;//设置直通滤波器对象
		pass.setInputCloud (cloud);//设置输入点云
		pass.setFilterFieldName ("z");//设置指定过滤的维度
		pass.setFilterLimits (near_cuting, far_cuting);//设置指定纬度过滤的范围
		pass.filter (*indices);//执行滤波，保存滤波结果在上述索引中
	}


	// 区域生长算法的5个参数
	pcl::RegionGrowing<pcl::PointXYZ, pcl::Normal> reg;//创建区域生长分割对象
	reg.setMinClusterSize (50);//设置一个聚类需要的最小点数
	reg.setMaxClusterSize (1000000);//设置一个聚类需要的最大点数
	reg.setSearchMethod (tree);//设置搜索方法
	reg.setNumberOfNeighbours (30);//设置搜索的临近点数目
	reg.setInputCloud (cloud);//设置输入点云
	if(Bool_Cuting)reg.setIndices (indices);//通过输入参数设置，确定是否输入点云索引
	reg.setInputNormals (normals);//设置输入点云的法向量
	reg.setSmoothnessThreshold (SmoothnessThreshold / 180.0 * M_PI);//设置平滑阈值
	reg.setCurvatureThreshold (CurvatureThreshold);//设置曲率阈值

	std::vector <pcl::PointIndices> clusters;
	reg.extract (clusters);//获取聚类的结果，分割结果保存在点云索引的向量中。
	end = time(0); 
	diff[2] = difftime (end, start)-diff[0]-diff[1]; 
	PCL_INFO ("\Region growing takes(seconds): %d\n", diff[2]); 

	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++;
	if (counter % 10 == 0)
	std::cout << std::endl;
	}
	std::cout << std::endl;
	*/
	PrintMemoryInfo();
	pcl::PointCloud <pcl::PointXYZRGB>::Ptr colored_cloud = reg.getColoredCloud ();
	pcl::visualization::CloudViewer viewer ("区域增长分割方法");
	viewer.showCloud(colored_cloud);
	while (!viewer.wasStopped ())
	{
	}//进行可视化

	return (0);
}
 
代码分析
 
然后设置最小和最大集群大小。这意味着在分割完成后，所有点小于最小值(或大于最大值)的聚类将被丢弃。最小值和最大值的默认值分别为1和“尽可能多”。
 
算法在内部结构中需要K最近邻搜索，所以这里是提供搜索方法并设置邻居数量的地方。之后，它接收到必须分割的点云、点下标和法线。
 
  pcl::RegionGrowing<pcl::PointXYZ, pcl::Normal> reg;
  reg.setMinClusterSize (50);
  reg.setMaxClusterSize (1000000);
  reg.setSearchMethod (tree);
  reg.setNumberOfNeighbours (30);
  reg.setInputCloud (cloud);
  //reg.setIndices (indices);
  reg.setInputNormals (normals);
 
这两行是算法初始化中最重要的部分，因为它们负责上述的平滑约束。第一种方法以弧度为单位设置角度，作为法向偏差的允许范围。如果点之间的法线偏差小于平滑阈值，那么建议它们在同一簇中(新的点-被测试的点-将被添加到簇中)。第二个是曲率阈值。如果两点有一个小的法向偏差，那么它们之间的曲率差被测试。
 
 reg.setSmoothnessThreshold (3.0 / 180.0 * M_PI);
  reg.setCurvatureThreshold (1.0);
 
RegionGrowing类提供了一个返回彩色云的方法，其中每个集群都有自己的颜色。因此，在这部分代码中，实例化pcl::visualization::CloudViewer以查看分割的结果——相同颜色的云 
 
pcl::PointCloud <pcl::PointXYZRGB>::Ptr colored_cloud = reg.getColoredCloud ();
  pcl::visualization::CloudViewer viewer ("Cluster viewer");
  viewer.showCloud(colored_cloud);
  while (!viewer.wasStopped ())
  {
  }

  return (0);
}
 
实验结果
 
原始点云：
 
 
 使用区域生长算法分割后的点云（每种颜色代表一个聚类）：
 
 
 在最后一张图片中，你可以看到彩色云有许多红点。这意味着这些点属于被拒绝的聚类，因为它们有太多/太少的点。
 
使用命令行进行输入： 
 
 
 
D:\PCLProject\pcl-project\pcl_segmentation\4_region_growing_segmentation\cmake_bin\Release>region_growing_segmentation.exe pig1.pcd -kn 50 -bc 0 -fc 10.0 -nc 0 -st 30 -ct 0.05
 Loading pcd file takes(seconds): 0
 Estimating normal takes(seconds): 1
 Region growing takes(seconds): 0
 Number of clusters is equal to 141
 First cluster has 136600 points.
 These are the indices of the points of the initial
 cloud that belong to the first cluster:
 
 
Process ID: 4294967295
         PageFaultCount: 0x00004A9E
         PeakWorkingSetSize: 0x03A69000
         WorkingSetSize: 0x03A68000
         QuotaPeakPagedPoolUsage: 0x0002DD40
         QuotaPagedPoolUsage: 0x0002D5B8
         QuotaPeakNonPagedPoolUsage: 0x00003508
         QuotaNonPagedPoolUsage: 0x00003480
         PagefileUsage: 0x03427000
         PeakPagefileUsage: 0x03427000
 
 
 其他数据结果
 
原始点云：
 
 
 
 处理之后的点云：
 
 
 

这一篇是补充前一篇的。相比前一篇的第二种实现方法，这里使用了均值、标准差随迭代更新的策略，且根据区域的性质，对三通道取不同的权值。修正后的程序为：
 
clear all;close all;clc;
img=imread('parrot.png');
subplot(3,1,1),
imshow(img);
[m,n,~]=size(img);
seeds=[400,200,1;230,600,2;440,670,3];  %第3列为label
labelNum=3;
newSeeds=seeds;
hold on;
plot(seeds(:,2)',seeds(:,1)','gs','linewidth',1);
title('原始图像及种子位置');
flag=1;
ycc=double(rgb2ycbcr(img));
Y=ycc(:,:,1);
Cb=ycc(:,:,2);
Cr=ycc(:,:,3);
img=double(img);
R=img(:,:,1);
G=img(:,:,2);
B=img(:,:,3);
mask=zeros(m,n);
for i=1:size(seeds,1)
    mask(seeds(i,1),seeds(i,2))=1;
end

avgs=zeros(labelNum,3);%初始均值
stds=repmat([30,10,10],[size(seeds,1),1]); %初始标准差
w=[5,5,5;3,2.2,2.8;3,3,3.5];  %各区域三通道权重

while(flag)
    temp=[];
    
    for i=1:labelNum
        index=find(seeds(:,3)==i);
        selected=seeds(index,1:2);
        indexS=m*(selected(:,2)-1)+selected(:,1);
        yS=Y(indexS);
        CbS=Cb(indexS);
        CrS=Cr(indexS);
        rgbSamp=[yS,CbS,CrS];
        avgs(i,:)=sum(rgbSamp,1)/length(index);
        if size(rgbSamp,1)>1
            stds(i,:)=std(rgbSamp,1);
        end
    end
    stds=max(5,stds);
    for i=1:size(newSeeds,1)
        rowS=newSeeds(i,1);
        colS=newSeeds(i,2);
        label=newSeeds(i,3);
        rlow=max(1,rowS-1);
        rhigh=min(m,rowS+1);
        clow=max(1,colS-1);
        chigh=min(n,colS+1);
        for row=rlow:rhigh
            for col=clow:chigh
                if mask(row,col)==1  %已经添加过的排除掉
                    continue;
                end
                diffycc=reshape(ycc(row,col,:),1,size(ycc,3))-avgs(label,:);
                judge=abs(diffycc)<w(label,:).*stds(label,:);
                if sum(judge)==3 % R、G、B、Y均在3倍均方根范围内
                    temp=[temp;row,col,label];
                    mask(row,col)=1;
                end
            end
        end
    end
    newSeeds=temp;
    seeds=[seeds;newSeeds];
    if isempty(temp)
        flag=0;
    end
end
imlabel=zeros(m,n);
for i=1:size(seeds,1)
    imlabel(seeds(i,1),seeds(i,2))=seeds(i,3);
end
mask=cat(3,mask,mask,mask);
sgImg=mask.*img;
imlabel=0.5*mat2gray(imlabel)+0.5;
subplot(3,1,2);
imshow(uint8(sgImg));
title('分割结果');
subplot(3,1,3);
imshow(imlabel);
title('种子点生长获得的掩模');
 
分割结果：
 
 
 

  
  
 
 
可以看到，这种方法比第一种方法获得了更好的分割效果。 
 
 

转载：      有梦想的田园犬 https://blog.csdn.net/AmbitiousRuralDog/article/details/80267519