该边界提取采用PCL库里边的方法，基于法线估计来实现的边界检测与提取：

首先从原始点云上计算出法线，再由法线结合数据估计出边界。（这样解释还是特别抽像吧）

------------法线求解：（平面的法线是垂直于它的单位向量。在点云的表面的法线被定义为垂直于与点云表面相切的平面的向量。法线提供了关于曲面的曲率信息）

对点云数据集的每个点的法线估计，可以看作是对表面法线的近似推断。（因此该表面的判断就是你寻找的周围几个点或者半径内几个点组成的平面，就是下述代码中reforn这个参数，该参数的设置一般设置为分辨率的10倍时，效果较好，主要是对于法线估计。邻域半径选择太小了，噪声较大，估计的法线就容易出错，而搜索邻域半径设置的太大估计速度就比较慢。）

----------求解原理

确定表面一点法线的问题近似于估计表面的一个相切面法线的问题，因此转换过来以后就变成一个最小二乘法平面拟合估计问题。

平面方程用法线式表示为:

为平面上点处法向量的方向余弦，|p|为原点到平面的距离。

     (此处a,b,c能够构成该平面的一个法向量n，后边的约束条件三方和为1只是控制了法向量大小，并没有改变方向，是为了下边利用点到平面的距离公式时，控制分母为1简化计算)

求平面方程即转化为求四个参数。









将上式整理为：                        (由该式子联系求特征值特征向量可以想到，为特征值，即为对应的特征向量)

即转化到了求解矩阵A的特征值与特征向量的问题，矩阵A即为n个点的协方差矩阵。即为该矩阵的一个特征向量。



//PCL中的NormalEstimation
#include <pcl/point_types.h>
#include <pcl/features/normal_3d.h>

{
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);

  ... read, pass in or create a point cloud ...

  // Create the normal estimation class, and pass the input dataset to it
  pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne;
  ne.setInputCloud (cloud);

  // Create an empty kdtree representation, and pass it to the normal estimation object.
  // Its content will be filled inside the object, based on the given input dataset (as no other search surface is given).
  pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ> ());
  ne.setSearchMethod (tree);

  // Output datasets
  pcl::PointCloud<pcl::Normal>::Ptr cloud_normals (new pcl::PointCloud<pcl::Normal>);

  // Use all neighbors in a sphere of radius 3cm
  ne.setRadiusSearch (0.03);

  // Compute the features
  ne.compute (*cloud_normals);

  // cloud_normals->points.size () should have the same size as the input cloud->points.size ()*
}

OpenMP加速法线估计

PCL提供了表面法线估计的加速实现，基于OpenMP使用多核/多线程来加速计算。 该类的名称是pcl :: NormalEstimationOMP，其API与单线程pcl :: NormalEstimation 100％兼容。 在具有8个内核的系统上，一般计算时间可以加快6-8倍。

include <pcl/point_types.h>
#include <pcl/features/normal_3d_omp.h>

{
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);

  ... read, pass in or create a point cloud ...

  // Create the normal estimation class, and pass the input dataset to it
  pcl::NormalEstimationOMP<pcl::PointXYZ, pcl::Normal> ne;
  ne.setNumberOfThreads(12);  // 手动设置线程数，否则提示错误
  ne.setInputCloud (cloud);

  // Create an empty kdtree representation, and pass it to the normal estimation object.
  // Its content will be filled inside the object, based on the given input dataset (as no other search surface is given).
  pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ> ());
  ne.setSearchMethod (tree);

  // Output datasets
  pcl::PointCloud<pcl::Normal>::Ptr cloud_normals (new pcl::PointCloud<pcl::Normal>);

  // Use all neighbors in a sphere of radius 3cm
  ne.setRadiusSearch (0.03);

  // Compute the features
  ne.compute (*cloud_normals);

  // cloud_normals->points.size () should have the same size as the input cloud->points.size ()*
}


上述完整的介绍了法线的估计过程，自认为只要有稍微的高数基础应该都能看懂的哈哈

下边接着开始，怎么根据求的法线来找边界呢？？?

基于法线完成的边界估计主要是利用各个法线方向之间的夹角来做的判断（所以有个设置角度的阈值参数）。（此处还是没有太明白怎么根据法线夹角确定哪些点是边界点，如有特别明白的可以下边评论留言。）

对于边界的估计就是这个函数boundEst.setRadiusSearch(re)，参数re也设置为分辨率（此处的分辨率指的是点云的密度）的10倍，太小则内部的很多点就都当成边界点了。最后一个参数是边界判断时的角度阈值，默认值为PI/2，此处设置为PI/4，用户也可以根据需要进行更改。

#include <iostream>
#include <pcl/console/parse.h>
#include <pcl/filters/extract_indices.h>
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <boost/thread/thread.hpp>
#include <pcl/features/boundary.h>
#include <math.h>
#include <boost/make_shared.hpp>
#include <pcl/point_cloud.h>
#include <pcl/visualization/cloud_viewer.h>
#include <pcl/visualization/range_image_visualizer.h>
#include <pcl/features/normal_3d.h>
#include <pcl/filters/covariance_sampling.h>
#include <pcl/filters/normal_space.h>
#include <pcl/kdtree/kdtree_flann.h>
#include <pcl/io/ply_io.h>
#include <pcl/filters/statistical_outlier_removal.h>
int estimateBorders(pcl::PointCloud<pcl::PointXYZI>::Ptr &cloud, float re, float reforn)
{
	pcl::PointCloud<pcl::Boundary> boundaries; //保存边界估计结果
	pcl::BoundaryEstimation<pcl::PointXYZI, pcl::Normal, pcl::Boundary> boundEst; //定义一个进行边界特征估计的对象
	pcl::NormalEstimation<pcl::PointXYZI, pcl::Normal> normEst; //定义一个法线估计的对象
	pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>); //保存法线估计的结果
	pcl::PointCloud<pcl::PointXYZI>::Ptr cloud_boundary(new pcl::PointCloud<pcl::PointXYZI>);
	normEst.setInputCloud(pcl::PointCloud<pcl::PointXYZI>::Ptr(cloud));
	normEst.setRadiusSearch(reforn); //设置法线估计的半径//normEst.setKSearch(10);//表示计算点云法向量时，搜索的点云个数
	normEst.compute(*normals); //将法线估计结果保存至normals
    //输出法线的个数
    std:cout << "reforn: " << reforn << std::endl;
    std::cerr << "normals: " << normals->size() << std::endl;

	boundEst.setInputCloud(cloud); //设置输入的点云
	boundEst.setInputNormals(normals); //设置边界估计的法线，因为边界估计依赖于法线
	boundEst.setRadiusSearch(re); //设置边界估计所需要的半径,//这里的Threadshold为一个浮点值，可取点云模型密度的10倍
	boundEst.setAngleThreshold(M_PI / 4); //边界估计时的角度阈值M_PI / 4  并计算k邻域点的法线夹角,若大于阈值则为边界特征点
	boundEst.setSearchMethod(pcl::search::KdTree<pcl::PointXYZI>::Ptr(new pcl::search::KdTree<pcl::PointXYZI>)); //设置搜索方式KdTree
	boundEst.compute(boundaries); //将边界估计结果保存在boundaries

	std::cerr << "AngleThreshold: " << M_PI / 4 << std::endl;
	//输出边界点的个数
	std::cerr << "boundaries: " << boundaries.points.size() << std::endl;
	//存储估计为边界的点云数据，将边界结果保存为pcl::PointXYZ类型
	for (int i = 0; i < cloud->points.size(); i++)
	{
		if (boundaries[i].boundary_point > 0)
		{
			cloud_boundary->push_back(cloud->points[i]);
		}
		
	}
	pcl::PCDWriter writer;
	std::stringstream ss;
	ss << "boundary" << ".pcd";
	writer.write<pcl::PointXYZI>(ss.str(), *cloud_boundary, false);
	//可视化显示原始点云与边界提取结果
	boost::shared_ptr<pcl::visualization::PCLVisualizer> MView(new pcl::visualization::PCLVisualizer("边界提取"));

	int v1(0);
	MView->createViewPort(0.0, 0.0, 0.5, 1.0, v1);
	MView->setBackgroundColor(0.3, 0.3, 0.3, v1);
	MView->addText("Raw point clouds", 10, 10, "v1_text", v1);
	int v2(0);
	MView->createViewPort(0.5, 0.0, 1, 1.0, v2);
	MView->setBackgroundColor(0.5, 0.5, 0.5, v2);
	MView->addText("Boudary point clouds", 80, 80, "v2_text", v2);
	MView->addPointCloud<pcl::PointXYZI>(cloud, "sample cloud", v1);
	MView->addPointCloud<pcl::PointXYZI>(cloud_boundary, "cloud_boundary", v2);
	MView->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1, 0, 0, "sample cloud", v1);
	MView->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 0, 1, 0, "cloud_boundary", v2);
	MView->addCoordinateSystem(1.0);
	MView->initCameraParameters();
	MView->spin();
	return 0;
}
int
main(int argc, char** argv)
{
	srand(time(NULL));
	float re, reforn;
	re = std::atof(argv[2]);
	reforn = std::atof(argv[3]);
	pcl::PointCloud<pcl::PointXYZI>::Ptr cloud_src(new pcl::PointCloud<pcl::PointXYZI>);
	pcl::io::loadPCDFile(argv[1], *cloud_src);
	 //创建滤波器对象
	pcl::PointCloud<pcl::PointXYZI>::Ptr cloud_filtered(new pcl::PointCloud<pcl::PointXYZI>);
	pcl::StatisticalOutlierRemoval<pcl::PointXYZI> sor;
	sor.setInputCloud(cloud_src);
	sor.setMeanK(100);//寻找每个点的50个最近邻点
	sor.setStddevMulThresh(3.0);//一个点的最近邻距离超过全局平均距离的一个标准差以上，就会舍弃
	sor.filter(*cloud_filtered);
	std::cout << "cloud_src: " << cloud_src->points.size() << std::endl;
	std::cout << "cloud_filtered: " << cloud_filtered->points.size() << std::endl;
	estimateBorders(cloud_src, re, reforn);
	return 0;
}



右键项目->属性->在命令行输入参数：如下所示：（以上需要三个参数）



结果如下：

                       

s

由Edelsbrunner H提出的alpha shapes算法是一种简单、有效的快速提取边界点算法。其克服了点云边界点形状影响的缺点，可快速准确提取边界点，其原理如下：

如下图所示，对于任意形状的平面点云，若一个半径为a的圆，绕其进行滚动。若滚动圆半径a足够小时，则点云中每一点均为边界点；若适当增大到一定程度，其只在边界点上进行滚动，其滚动的轨迹为点云边界。







效果图如下：

      

参考文献：

[1]平面点云边界提取算法研究[D].长沙理工大学,2017.

文章目录
1.用霍夫变换提取直线
2.边界跟踪


图像的基本特征之一是直线。一般物体平面图像的轮廓可近似为直线及弧线的组合，对物体轮廓的检测与识别可以转化为对这些基元的检测与提取。另外，在运动图像分析和估计领域也可以采用直线对应法实现刚体旋转量和位移量的测量，所以直线检测对图像算法的研究具有重要意义。

边缘是一个局部的概念，一个区域的边界是一个整体的概念，边界追踪是一种串行的图像分割技术。图像由于噪声以及光照不均匀等原因，边缘点可能是不连续的，边界追踪可以将其变为有意义的信息。下面分别介绍直线的提取与边界追踪。

1.用霍夫变换提取直线
Hough变换是一种利用图像的全局特征将特定形状的边缘连接起来，形成连续平滑边缘的一种方法。它通过将原图像上的点映射到用于累加的参数空间，实现对已知解析式曲线的识别。由于利用了图像的全局特性，所以受噪声和边界间断的影响较小，鲁棒性能好。Hough变换常用来对图像中的直线进行识别。

图像上任意直线区域都可以一一对应参数空间中的点，而图像上的任意像素都同时存在于很多直线区域之上。将图像上的直线区域想象为容器，把特定像素想象成放在容器中的棋子，只不过在这里，每个棋子都可以同时存在于多个容器中。那么Hough变换可以理解为依次检查图像上的每个棋子（特定像素），对于每个棋子，找到所有包含它的容器（平面上的直线区域），并为每个容器的计数器加1，这样就可以统计出每个容器（平面上的直线区域）所包含的棋子（特定像素数量）。当图像上某个直线区域包含的特定像素足够多时，就可以认为这个直线区域表示的直线存在。

在MATLAB中，通过Hough变换在图像中检测直线需要以下3个步骤。

（1）利用hough()函数执行霍夫变换，得到霍夫矩阵；

（2）利用houghpeaks（）函数在霍夫矩阵中寻找峰值点；

（3）利用houghlines（）函数在前两步结果的基础上得到原二值图像中的直线信息。这些函数的调用方法为
[H,theta,rho]=hough(BW,param1,vall,param2,val2)

其中，BW是边缘检测后的图像；param1、val1及param2、val2为可选参数；H是变换得到的Hough矩阵；theta和rho为分别对应于Hough矩阵每一列和每一行的θ和ρ值组成的向量。
peaks=houghpeaks(H,numpeaks,param1,val1,param2,val2)

其中，H是由hough（）函数得到的Hough矩阵；numpeaks是要寻找的峰值数目，默认为1；peaks是一个Q×2的矩阵，每行的两个元素分别为某一峰值点在Hough矩阵中的行、列索引，Q为找到的峰值点的数目。
lines=houghlines(BW,theta,rho,peaks,param1,val1,param2,val2)

其中，BW是边缘检测后的图像；theta和rho是Hough矩阵每一列和每一行的θ和ρ值组成的向量，由hough函数返回；peaks是一个包含峰值点信息的Q×2的矩阵，由houghpeaks函数返回；lines是一个结构体数组，数组长度是找到的直线条数。

对图像进行Hough变换提取直线
**clear all
close all
RGB=imread('peppers.png');
I=rgb2gray(RGB);
BW=edge(I,'Canny');
[H,T,R]=hough(BW,'RhoResolution',0.5,'Theta',-90:0.5:89.5);%计算霍夫变换
subplot(121)
imshow(RGB)
title('原始图像')
subplot(122)
imshow(imadjust(mat2gray(H)),'XData',T,'YData',R,...
    'InitialMagnification','fit');
title('Hough变换')
xlabel('\theta轴'),ylabel('\rho轴');
axis on
axis normal
hold on
colormap(hot)



将一幅图像进行Hough变换，标出峰值位置
clear all
close all
I=imread('eight.tif');
subplot(121)
imshow(I)
title('原始图像')
BW=edge(imrotate(I,50,'crop'),'prewitt');
[H,T,R]=hough(BW);
P=houghpeaks(H,2);
subplot(122)
imshow(H,[],'XData',T,'YData',R,'InitialMagnification','fit');%显示Hough变换
title('Hough变换')
xlabel('\theta轴'),ylabel('\rho轴')
axis on,axis normal
hold on
plot(T(P(:,2)),R(P(:,1)),'s','color','white');%标记颜色为白色



对图像进行Hough变换，并标出最长的线段
clear all
close all
I=imread('blobs.png');
subplot(131)
imshow(I)
title('原始图像')
rotI=imrotate(I,45,'crop');
BW=edge(rotI,'sobel');
[H,T,R]=hough(BW);
subplot(132)
imshow(H,[],'XData',T,'YData',R,...
    'InitialMagnification','fit')
title('Hough变换')
xlabel('\theta轴'),ylabel('\rho轴')
axis on,axis normal
hold on
P=houghpeaks(H,5,'threshold',ceil(0.3*max(H(:))));%寻找极值点
x=T(P(:,2));
y=R(P(:,1));
plot(x,y,'s','color','white');
lines=houghlines(BW,T,R,P,'FillGap',5,'MinLength',7);%找出对应的直线边缘
subplot(133)
imshow(rotI)
title('检测线段')
hold on
max_len=0;
for k=1:length(lines)
    xy=[lines(k).point1;lines(k).point2];
    plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');
    %标记直线边缘对应的起点
    plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','blue')
    plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red')
    len=norm(lines(k).point1-lines(k).point2);%计算直线边缘长度
    if(len>max_len)
        max_len=len;
        xy_long=xy;
    end
end
plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color','b')


2.边界跟踪
边界追踪是指从图像的一个边缘点出发，然后根据某种判别准则搜索下一个边缘点以此跟踪目标边界。边界跟踪的算法步骤如下：

（1）确定边界的起始搜索点，起始点的选择很关键，对某些图像，选择不同的起始点会导致不同的结果

（2）确定合适边界判别准则和搜索准则，判别准则用于判断一个点是不是边界点，搜索准则用于指导如何搜索下一个边缘点；

（3）确定搜索终止的条件

在MATLAB中，bwtraceboundary函数采用基于曲线跟踪的策略，给定搜索起始点和搜索方向及其返回该起始点的一条边界；bwboundaries函数用于获取二值图像中对象的轮廓。这两种函数的用法如下所示：
B=bwtraceboundary(BW,P,fstep)
B=bwtraceboundary(BW,P,fstep,conn)
B=bwtraceboundary(...,N,dir)

其中，非零像素表示对象，零像素构成背景；bwtraceboundary返回值为B，一个Q×2的矩阵，Q是区域边界像素的数量；B保存有边界像素的行、列坐标；参数P是一个指定行、列坐标的二元矢量，表示对象边界上开始跟踪的那个点；fstep表示初始查找方向，用于寻找对象中与P相连的下一个像素，fstep是寻找下一个相连像素点的初始方向；dir是寻找边界的方向，即顺时针还是逆时针。
B=bwboundaries(BW,conn)

其中，B是一个P×1数组，P为对象个数，每个对象是Q×2的矩阵，对应于对象轮廓像素的坐标。

对图像进行边界跟踪提取。
clear all
close all
RGB=imread('saturn.png');
figure
subplot(131)
imshow(RGB)
title('原始图像')
I=rgb2gray(RGB);
threshold=graythresh(I);
BW=im2bw(I,threshold);
subplot(132)
imshow(BW)
title('二值图像')
dim=size(BW);
col=round(dim(2)/2)-90;%计算起始点的列坐标
row=find(BW(:,col),1);%计算起始点行坐标
connectivity=8;
num_points=180;
contour=bwtraceboundary(BW,[row col],'N',connectivity,num_points);%提取边界
subplot(133)
imshow(RGB)
hold on
plot(contour(:,2),contour(:,1),'g','LineWidth',2)
title('结果图像')



对图像进行边界跟踪提取
clear all
close all
I=imread('blobs.png');
figure(1)
subplot(121)
imshow(I)
title('原始图像')
B=bwboundaries(I);
D=B{1,1};
subplot(122)
plot(D(:,2),D(:,1))
set(gca,'YDir','reverse')
title('边界标记后的图像')



对图像中不同的区域表示不同的颜色
clear all
close all
I=imread('tape.png');
figure(1)
subplot(131)
imshow(I)
title('原始图像')
BW=im2bw(I,graythresh(I));%生成二值图像
subplot(132)
imshow(BW)
title('二值图像')
[B,L]=bwboundaries(BW,'noholes');%提取边界，并返回边界元胞数组B和区域标志数组L
subplot(133)
imshow(label2rgb(L,@jet,[.5,.5,.5]))%以不同的颜色标志不同的区域
title('彩色标机图像')
hold on
for k=1:length(B)
    boundary=B{k};
    plot(boundary(:,2),boundary(:,1),'w','LineWidth',2)%在图像上叠画边界
end