图像中关键点提取和pcl中点云的关键点提取

2.图像中关键点提取和pcl中点云的关键点提取
2.1二维图像的Harris角点提取
经典的Harris角点提取方法实质上就是对Moravec算子的改良和优化。从图像分析的角度来定义角点可以有以下两种定义：
a.角点可以是两个边缘的的角点；
b.角点是邻域内具有两个主方向的特征点；
Harris角点检测原理主要参考：Harris角点提取和Harris角点检测原理
Harris角点检测使用了图像窗口模板在待检测图像上在各个方向上平移，从而探测角点，图像窗口为：[u,v],产生灰度变化的自相关函数如下：

其中：w(u,v)是窗口函数，I(∆x+u,∆y+v)为平移后的灰度图像，I(u,v)为灰度图像。
窗口函数（权重矩阵）可以是平坦的，也可以是高斯的（通常为高斯滤波器）如下图：


二次项函数本质上是一个椭圆函数，椭圆的长半径和短半径由λ_1 、λ_2(即特征值)决定的，方向由M(x,y)的特征向量决定，椭圆方程为：[∆x,∆y]M(x,y)[(∆x，∆y)]T=1,（T表示转置）椭圆示意图如下：

椭圆函数特征值与图像中的角点、直线（边缘）和平面之间的关系如下图所示。共可分为三种情况：

a.一个特征值大，另一个特征值小，自相关函数值在某一方向上大，在其他方向上小，则表示检测到边；
b.两个特征值都小，且近似相等；自相关函数数值在各个方向上都小，则表示检测到图像中的平面（平坦部分）；
c.两个特征值都大，且近似相等，自相关函数在所有方向都增大。表示检测到了角点。
二次项函数可以求得特征值，但是Harris方法不必计算具体的特征值，而是计算角点响应值R来判断角点。R的计算公式为：

Harris角点检测实现步骤大致如下：
1.计算图像I(x,y)在X和Y两个方向的梯度Ix 、Iy;
2.计算图像两个方向的梯度乘积；
3.使用高斯函数对
进行高斯加权，生成矩阵M的元素A、B和C;
4.计算每个像素的Harris响应值R,并对小于某一阈值t的R置为0；
5.在3X3或5X5的邻域内进行非极大值抑制，局部最大值即为图像的角点
Harris角点检测特点：
a.增大α的值，将减小角点响应值R，降低角点检测的灵敏性，减少被检测角点数量，减小α值，将增大角点响应值R,增加角点检测的灵敏性，增加被检测角点数量；
b.Harris角点检测对亮度和对比度的变化不敏感；
c.Harris角点检测具有旋转不变性；
d.Harris角点检测具有尺度不变性
2.2Harris角点实例分析：

#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <iostream>
#include "opencv2/imgproc/types_c.h"
using namespace cv;
using namespace std;


Mat image;
Mat imageGray;
int thresh = 200;
int MaxThresh = 255;

void Trackbar(int, void*);  //阈值控制

int main()
{
	image = imread("rectified.JPG");
	cvtColor(image, imageGray, CV_RGB2GRAY);
	GaussianBlur(imageGray, imageGray, Size(5, 5), 1); // 滤波
	namedWindow("Corner Detected", WINDOW_NORMAL);

	createTrackbar("threshold：", "Corner Detected", &thresh, MaxThresh, Trackbar);
	imshow("Corner Detected", image);
	Trackbar(0, 0);
	waitKey();
	return 0;
}

void Trackbar(int, void*)
{
	Mat dst, dst8u, dstshow, imageSource;
	dst = Mat::zeros(image.size(), CV_32FC1);
	imageSource = image.clone();
	cornerHarris(imageGray, dst, 3, 3, 0.04, BORDER_DEFAULT);
	normalize(dst, dst8u, 0, 255, CV_MINMAX);  //归一化
	convertScaleAbs(dst8u, dstshow);
	imshow("dst", dstshow);  //dst显示
	for (int i = 0; i < image.rows; i++)
	{
		for (int j = 0; j < image.cols; j++)
		{
			if (dstshow.at<uchar>(i, j) > thresh)  //阈值判断，大于设定的阈值则显示出角点
			{
				circle(imageSource, Point(j, i), 2, Scalar(0, 0, 255), 2); //标注角点
			}
		}
	}
	imshow("Corner Detected", imageSource);
}

角点显示如下，可以鼠标拖动滑动块显示出不同的阈值下的角点数量
thresh=200时：

thresh=149时：

2.3点云的关键点提取
3DHarris点云关键点提取借鉴了二维图像的关键点提取思想，不过3DHarris关键点检测使用的是点云表面法向量的信息，这区别于2DHarris角点检测使用的图像梯度信息。
3DHarris角点提取步骤：

1. 求解离散点云的法线，构造法线协方差矩阵，法线协方差矩阵如下：

式中n为当前点邻域点个数。
2. 根据角点响应函数计算点云中每个点的角点响应值，计算公式为：

3. 根据角点响应阈值判断当前点是否为角点，若R3d大于给定的阈值，同时为局部极大值点则判定该点为角点，否则不是角点。
注意：Harris输出的点云必须有强度信息（I）的pcl::PointXYZI,因为评估值保存在（I）里面，因此后续关键点的保存和可视化中不能直接使用Harris输出的pcl::PointXYZI格式点云，需要通过关键点的索引从原始点云中重新获取。
点云Harris角点信息提取代码实现如下：
主要参考：郭浩，苏伟，朱德海，王可.点云库pcl从入门到精通「M」.北京:机械工业出版社，2019:182-190

#include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/visualization/cloud_viewer.h>
#include <pcl/visualization/pcl_visualizer.h>
#include <pcl/keypoints/harris_3d.h>
#include <boost/thread/thread.hpp>
#include <stdlib.h>
#include <iostream>

using namespace std;
using namespace std;
int main(int argc, char** argv)
{

	// 创建点云对象指针
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);

	pcl::io::loadPCDFile<pcl::PointXYZ>("rabbit.pcd", *cloud);
	boost::shared_ptr<pcl::visualization::PCLVisualizer> viewer(new pcl::visualization::PCLVisualizer("3D Viewer"));
	viewer->addPointCloud<pcl::PointXYZ>(cloud, "sample cloud1");

	//harris detect

	pcl::HarrisKeypoint3D<pcl::PointXYZ, pcl::PointXYZI, pcl::Normal> harris;
	harris.setInputCloud(cloud);//设置输入点云 指针
	harris.setNonMaxSupression(true);
	harris.setRadius(0.6f);//　块体半径
	harris.setThreshold(0.01f);//数量阈值
	//新建的点云必须初始化，清零，否则指针会越界
	//注意Harris的输出点云必须是有强度(I)信息的 pcl::PointXYZI，因为评估值保存在I分量里
	pcl::PointCloud<pcl::PointXYZI>::Ptr cloud_out_ptr(new pcl::PointCloud<pcl::PointXYZI>);

	// 计算特征点
	harris.compute(*cloud_out_ptr);
	// 关键点
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_harris_ptr(new pcl::PointCloud<pcl::PointXYZ>);
	int size = cloud_out_ptr->size();
	pcl::PointXYZ point;
	//可视化结果不支持XYZI格式点云，所有又要导回XYZ格式。。。。
	for (int i = 0; i < size; ++i)
	{
		point.x = cloud_out_ptr->at(i).x;
		point.y = cloud_out_ptr->at(i).y;
		point.z = cloud_out_ptr->at(i).z;
		cloud_harris_ptr->push_back(point);
	}
	pcl::visualization::PointCloudColorHandlerCustom<pcl::PointXYZ> harris_color_handler(cloud_harris_ptr, 0, 255, 0);//第一个参数类型为　指针
	viewer->addPointCloud<pcl::PointXYZ>(cloud_harris_ptr, harris_color_handler, "harris");//第一个参数类型为　指针
	viewer->setPointCloudRenderingProperties(pcl::visualization::PCL_VISUALIZER_POINT_SIZE, 5, "harris");

	while (!viewer->wasStopped())
	{
		viewer->spinOnce(100);

	}
	system("pause");
	return 0;
}

总结:它的参数设置主要在设定搜索半径和阈值上面.
显示效果如图：

注：3D点云的关键提取算子除了有3DHarris还有3DSIFT、NARF(Normal Aligned Radial Feature)中文名：法向对准径向特征，这些算子仍然使用的是点云的法向量信息。