2016-12-14 11:31:09 jia20003 阅读数 10714

图像处理之角点检测与亚像素角点定位

 

角点是图像中亮度变化最强地方反映了图像的本质特征,提取图像中的角点可以有效提高图像处理速度与精准度。所以对于整张图像来说特别重要,角点检测与提取的越准确图像处理与分析结果就越接近真实。同时角点检测对真实环境下的对象识别、对象匹配都起到决定性作用。Harris角点检测是图像处理中角点提取的经典算法之一,应用范围广发,在经典的SIFT特征提取算法中Harris角点检测起到关键作用。通常对角点检测算法都有如下要求:


1. 基于灰度图像、能够自动调整运行稳定,检测出角点的数目。

2. 对噪声不敏感、有一定的噪声抑制,有较强的角点角点检测能力。

3. 准确性够高,能够正确发现角点位置

4. 算法尽可能的快与运行时间短


Harris角点检测基本上满足了上述四点要求,所以被广发应用,除了Harris角点检测,另外一种常见的角点检测算法-Shi-Tomasi角点检测也得到了广发应用,OpenCV中对这两种算法均有实现API可以调用。关于Harris角点检测原理可以看我之前写的博文:

http://blog.csdn.net/jia20003/article/details/16908661

关于Shi-Tomasi角点检测,与Harris角点检测唯一不同就是在计算角点响应值R上面。


然后根据输入的阈值T大于该阈值的R对应像素点即为图像中角点位置坐标。此刻坐标往往都是整数出现,而在真实的世界中坐标多数时候都不是整数,假设我们计算出来的角点位置P(34, 189)而实际上准确角点位置是P(34.278, 189.706)这样带小数的位置,而这样的准确位置寻找过程就叫做子像素定位或者亚像素定位。这一步在SURF与SIFT算法中都有应用而且非常重要。常见的亚像素级别精准定位方法有三类:

1. 基于插值方法

2. 基于几何矩寻找方法

3. 拟合方法 - 比较常用

拟合方法中根据使用的公式不同可以分为高斯曲面拟合与多项式拟合等等。以高斯拟合为例


这样就求出了亚像素的位置。使用亚像素位置进行计算得到结果将更加准确,对图像特征提取、匹配结果效果显著。OpenCV中已经对角点检测实现了亚像素级别的API可以调用。

代码演示

OpenCV亚像素角点检测例子:

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace cv;
using namespace std;

Mat src, gray_src;
int max_corners = 10;
int max_trackbar = 30;
const char* output_title = "subpxiel-result";
void GoodFeature2Track_Demo(int, void*);
int main(int argc, char** argv) {
	src = imread("D:/vcprojects/images/home.jpg");
	if (src.empty()) {
		printf("could not load image...\n");
		return -1;
	}
	cvtColor(src, gray_src, COLOR_BGR2GRAY);
	namedWindow("input", CV_WINDOW_AUTOSIZE);
	namedWindow(output_title, CV_WINDOW_AUTOSIZE);
	imshow("input", src);

	createTrackbar("Corners:", output_title, &max_corners, max_trackbar, GoodFeature2Track_Demo);
	GoodFeature2Track_Demo(0, 0);

	waitKey(0);
	return 0;
}

void GoodFeature2Track_Demo(int, void*) {
	if (max_corners < 1) {
		max_corners = 1;
	}
	vector<Point2f> corners;
	double qualityLevel = 0.01;
	double minDistance = 10;
	int blockSize = 3;
	double k = 0.04;
	goodFeaturesToTrack(gray_src, corners, max_corners, qualityLevel, minDistance, Mat(), blockSize, false, k);
	cout << "number of corners : " << corners.size() << endl;
	Mat copy = src.clone();
	for (size_t t = 0; t < corners.size(); t++) {
		circle(copy, corners[t], 4, Scalar(255, 0, 0), 2, 8, 0);
	}
	imshow(output_title, copy);

	// locate corner point on sub pixel level
	Size winSize = Size(5, 5);
	Size zerozone = Size(-1, -1);
	TermCriteria criteria = TermCriteria(TermCriteria::EPS + TermCriteria::MAX_ITER, 40, 0.001);
	cornerSubPix(gray_src, corners, winSize, zerozone, criteria);
	for (size_t t = 0; t < corners.size(); t++) {
		cout << (t+1) << ".point[x, y]=" << corners[t].x << "," << corners[t].y << endl;
	}

	return;
}

原图如下:


运行结果:


转载请注明来自【jia20003】的博客!


2018-11-20 00:18:48 Godsolve 阅读数 299

实现Harris角点检测算法,并与OpenCV的cornerHarris函数的结果进行比较。

特征点在图像中一般有具体的坐标,并具有某些数学特征,如局部最大或最小灰度、以及某些梯度特征等。角点可以简单的认为是两条边的交点。如下图所示:

在这里插入图片描述
在各个方向上移动小窗口,如果在所有方向上移动,窗口内灰度都发生变化,则认为是角点;如果任何方向都不变化,则是均匀区域;如果灰度只在一个方向上变化,则可能是图像边缘。
而Harris角点检测算法就是基于图像的这种特点来进行计算的。该算法的主要步骤为:

  1. 计算图像I(x,y)I(x,y)在XX方向和YY方向的梯度
    在这里插入图片描述
  2. 计算图像两个方向梯度的乘积在这里插入图片描述
  3. 使用窗口高斯函数分别对I_x2、I_y2、IxIy进行高斯加权,生成矩阵M。
    在这里插入图片描述
    在这里插入图片描述
  4. 计算每个像素的Harris响应值R,并设定一阈值T,对小于阈值T的R置零。
    在这里插入图片描述
    在这里插入图片描述
  5. 在一个固定窗口大小的邻域内(5×55×5)进行非极大值抑制,局部极大值点即为图像中的角点。
    在这里插入图片描述

运行结果为:
在这里插入图片描述
改了一下响应函数的参数值:
在这里插入图片描述
但是发现这种方法做出来的效果并不是很好,虽然找出来了一些角点,但是还有很多角点没有找到,而且找到了很多错误的或者在当前阈值下不应该出现的角点。

之后找到了opencv中自带的角点检测函数,并用其对相同的图片进行了处理,结果如下:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

后来又查了些资料,使用cornerEigenValsAndVecs()函数和minMaxLoc()函数来重新写了一下角点检测算法,测试结果还不错:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

角点检测(Corner Detection)也称为特征点检测,是图像处理和计算机视觉中用来获取图像局部特征点的一类方法,广泛应用于运动检测、图像匹配、视频跟踪、三维建模以及目标识别等领域中。
在实验过程中,也发现了Harris角点的一些性质:
参数α对角点检测的影响:增大α的值,将减小角点响应值R,减少被检测角点的数量;减小α的值,将增大角点响应值R,增加被检测角点的数量。
此外,该方法还可以扩展到多尺度Harris角点检测,但是在此次实验中没有实现,可以考虑在之后加以实现。


同时也欢迎各位关注我的微信的公众号 南木的下午茶

在这里插入图片描述


代码:

// CVE8.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//

#include "pch.h"
#include <iostream>
#include <string>
#include <list>
#include <vector>
#include <map>
#include <cmath>
#include <stack>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <stdio.h>
using namespace std;
using namespace cv;

/*
RGB转换成灰度图像的一个常用公式是:
Gray = R*0.299 + G*0.587 + B*0.114
*/
//******************灰度转换函数*************************  
//第一个参数image输入的彩色RGB图像的引用;  
//第二个参数imageGray是转换后输出的灰度图像的引用;  
//*******************************************************
void ConvertRGB2GRAY(const Mat &image, Mat &imageGray);

//******************Sobel卷积因子计算X、Y方向梯度和梯度方向角********************  
//第一个参数imageSourc原始灰度图像;  
//第二个参数imageSobelX是X方向梯度图像;  
//第三个参数imageSobelY是Y方向梯度图像;  
//第四个参数pointDrection是梯度方向角数组指针  
//*************************************************************  
void SobelGradDirction(Mat &imageSource, Mat &imageSobelX, Mat &imageSobelY);

//******************计算Sobel的X方向梯度幅值的平方*************************  
//第一个参数imageGradX是X方向梯度图像;    
//第二个参数SobelAmpXX是输出的X方向梯度图像的平方  
//*************************************************************  
void SobelXX(const Mat imageGradX, Mat_<float> &SobelAmpXX);

//******************计算Sobel的Y方向梯度幅值的平方*************************    
//第一个参数imageGradY是Y方向梯度图像;  
//第二个参数SobelAmpXX是输出的Y方向梯度图像的平方  
//*************************************************************  
void SobelYY(const Mat imageGradY, Mat_<float> &SobelAmpYY);

//******************计算Sobel的XY方向梯度幅值的乘积*************************    
//第一个参数imageGradX是X方向梯度图像;
//第二个参数imageGradY是Y方向梯度图像;
//第二个参数SobelAmpXY是输出的XY方向梯度图像 
//*************************************************************  
void SobelXY(const Mat imageGradX, const Mat imageGradY, Mat_<float> &SobelAmpXY);

//****************计算一维高斯的权值数组*****************
//第一个参数size是代表的卷积核的边长的大小
//第二个参数sigma表示的是sigma的大小
//*******************************************************
double *getOneGuassionArray(int size, double sigma);

//****************高斯滤波函数的实现*****************
//第一个参数srcImage是代表的输入的原图
//第二个参数dst表示的是输出的图
//第三个参数size表示的是卷积核的边长的大小
//*******************************************************
void MyGaussianBlur(Mat_<float> &srcImage, Mat_<float> &dst, int size);

//****计算局部特涨结果矩阵M的特征值和响应函数H = (A*B - C) - k*(A+B)^2******
//M
//A  C
//C  B
//Tr(M)=a+b=A+B
//Det(M)=a*b=A*B-C^2
//计算输出响应函数的值得矩阵
//****************************************************************************
void harrisResponse(Mat_<float> &GaussXX, Mat_<float> &GaussYY, Mat_<float> &GaussXY, Mat_<float> &resultData, float k);


//***********非极大值抑制和满足阈值及某邻域内的局部极大值为角点**************
//第一个参数是响应函数的矩阵
//第二个参数是输入的灰度图像
//第三个参数表示的是输出的角点检测到的结果图
void LocalMaxValue(Mat_<float> &resultData, Mat &srcGray, Mat &ResultImage, int kSize);

int main()
{
	const Mat srcImage = imread("E:/C++/CVE8/img.jpg");
	if (!srcImage.data)
	{
		printf("could not load image...\n");
		return -1;
	}
	imshow("srcImage", srcImage);
	Mat srcGray;
	ConvertRGB2GRAY(srcImage, srcGray);
	Mat imageSobelX;
	Mat imageSobelY;
	Mat resultImage;
	Mat_<float> imageSobelXX;
	Mat_<float> imageSobelYY;
	Mat_<float> imageSobelXY;
	Mat_<float> GaussianXX;
	Mat_<float> GaussianYY;
	Mat_<float> GaussianXY;
	Mat_<float> HarrisRespond;
	//计算Soble的XY梯度
	SobelGradDirction(srcGray, imageSobelX, imageSobelY);
	//计算X方向的梯度的平方
	SobelXX(imageSobelX, imageSobelXX);
	SobelYY(imageSobelY, imageSobelYY);
	SobelXY(imageSobelX, imageSobelY, imageSobelXY);
	//计算高斯模糊XX YY XY
	MyGaussianBlur(imageSobelXX, GaussianXX, 3);
	MyGaussianBlur(imageSobelYY, GaussianYY, 3);
	MyGaussianBlur(imageSobelXY, GaussianXY, 3);
	harrisResponse(GaussianXX, GaussianYY, GaussianXY, HarrisRespond, 0.02);
	LocalMaxValue(HarrisRespond, srcGray, resultImage, 3);
	imshow("imageSobelX", imageSobelX);
	imshow("imageSobelY", imageSobelY);
	imshow("resultImage", resultImage);
	waitKey(0);
	return 0;
}
void ConvertRGB2GRAY(const Mat &image, Mat &imageGray)
{
	if (!image.data || image.channels() != 3)
	{
		return;
	}
	//创建一张单通道的灰度图像
	imageGray = Mat::zeros(image.size(), CV_8UC1);
	//取出存储图像像素的数组的指针
	uchar *pointImage = image.data;
	uchar *pointImageGray = imageGray.data;
	//取出图像每行所占的字节数
	size_t stepImage = image.step;
	size_t stepImageGray = imageGray.step;
	for (int i = 0; i < imageGray.rows; i++)
	{
		for (int j = 0; j < imageGray.cols; j++)
		{
			pointImageGray[i*stepImageGray + j] = (uchar)(0.114*pointImage[i*stepImage + 3 * j] + 0.587*pointImage[i*stepImage + 3 * j + 1] + 0.299*pointImage[i*stepImage + 3 * j + 2]);
		}
	}
}


//存储梯度膜长
void SobelGradDirction(Mat &imageSource, Mat &imageSobelX, Mat &imageSobelY)
{
	imageSobelX = Mat::zeros(imageSource.size(), CV_32SC1);
	imageSobelY = Mat::zeros(imageSource.size(), CV_32SC1);
	//取出原图和X和Y梯度图的数组的首地址
	uchar *P = imageSource.data;
	uchar *PX = imageSobelX.data;
	uchar *PY = imageSobelY.data;

	//取出每行所占据的字节数
	int step = imageSource.step;
	int stepXY = imageSobelX.step;

	int index = 0;//梯度方向角的索引
	for (int i = 1; i < imageSource.rows - 1; ++i)
	{
		for (int j = 1; j < imageSource.cols - 1; ++j)
		{
			//通过指针遍历图像上每一个像素   
			double gradY = P[(i + 1)*step + j - 1] + P[(i + 1)*step + j] * 2 + P[(i + 1)*step + j + 1] - P[(i - 1)*step + j - 1] - P[(i - 1)*step + j] * 2 - P[(i - 1)*step + j + 1];
			PY[i*stepXY + j * (stepXY / step)] = abs(gradY);

			double gradX = P[(i - 1)*step + j + 1] + P[i*step + j + 1] * 2 + P[(i + 1)*step + j + 1] - P[(i - 1)*step + j - 1] - P[i*step + j - 1] * 2 - P[(i + 1)*step + j - 1];
			PX[i*stepXY + j * (stepXY / step)] = abs(gradX);
		}
	}
	//将梯度数组转换成8位无符号整型
	convertScaleAbs(imageSobelX, imageSobelX);
	convertScaleAbs(imageSobelY, imageSobelY);
}


void SobelXX(const Mat imageGradX, Mat_<float> &SobelAmpXX)
{
	SobelAmpXX = Mat_<float>(imageGradX.size(), CV_32FC1);
	for (int i = 0; i < SobelAmpXX.rows; i++)
	{
		for (int j = 0; j < SobelAmpXX.cols; j++)
		{
			SobelAmpXX.at<float>(i, j) = imageGradX.at<uchar>(i, j)*imageGradX.at<uchar>(i, j);
		}
	}
	//convertScaleAbs(SobelAmpXX, SobelAmpXX);
}

void SobelYY(const Mat imageGradY, Mat_<float> &SobelAmpYY)
{
	SobelAmpYY = Mat_<float>(imageGradY.size(), CV_32FC1);
	for (int i = 0; i < SobelAmpYY.rows; i++)
	{
		for (int j = 0; j < SobelAmpYY.cols; j++)
		{
			SobelAmpYY.at<float>(i, j) = imageGradY.at<uchar>(i, j)*imageGradY.at<uchar>(i, j);
		}
	}
	//convertScaleAbs(SobelAmpYY, SobelAmpYY);
}

void SobelXY(const Mat imageGradX, const Mat imageGradY, Mat_<float> &SobelAmpXY)
{
	SobelAmpXY = Mat_<float>(imageGradX.size(), CV_32FC1);
	for (int i = 0; i < SobelAmpXY.rows; i++)
	{
		for (int j = 0; j < SobelAmpXY.cols; j++)
		{
			SobelAmpXY.at<float>(i, j) = imageGradX.at<uchar>(i, j)*imageGradY.at<uchar>(i, j);
		}
	}
	//convertScaleAbs(SobelAmpXY, SobelAmpXY);
}



//计算一维高斯的权值数组
double *getOneGuassionArray(int size, double sigma)
{
	double sum = 0.0;
	//定义高斯核半径
	int kerR = size / 2;

	//建立一个size大小的动态一维数组
	double *arr = new double[size];
	for (int i = 0; i < size; i++)
	{

		// 高斯函数前的常数可以不用计算,会在归一化的过程中给消去
		arr[i] = exp(-((i - kerR)*(i - kerR)) / (2 * sigma*sigma));
		sum += arr[i];//将所有的值进行相加

	}
	//进行归一化	
	for (int i = 0; i < size; i++)
	{
		arr[i] /= sum;
		cout << arr[i] << endl;
	}
	return arr;
}

void MyGaussianBlur(Mat_<float> &srcImage, Mat_<float> &dst, int size)
{
	CV_Assert(srcImage.channels() == 1 || srcImage.channels() == 3); // 只处理单通道或者三通道图像
	int kerR = size / 2;
	dst = srcImage.clone();
	int channels = dst.channels();
	double* arr;
	arr = getOneGuassionArray(size, 1);//先求出高斯数组

									   //遍历图像 水平方向的卷积
	for (int i = kerR; i < dst.rows - kerR; i++)
	{
		for (int j = kerR; j < dst.cols - kerR; j++)
		{
			float GuassionSum[3] = { 0 };
			//滑窗搜索完成高斯核平滑
			for (int k = -kerR; k <= kerR; k++)
			{

				if (channels == 1)//如果只是单通道
				{
					GuassionSum[0] += arr[kerR + k] * dst.at<float>(i, j + k);//行不变,列变换,先做水平方向的卷积
				}
				else if (channels == 3)//如果是三通道的情况
				{
					Vec3f bgr = dst.at<Vec3f>(i, j + k);
					auto a = arr[kerR + k];
					GuassionSum[0] += a * bgr[0];
					GuassionSum[1] += a * bgr[1];
					GuassionSum[2] += a * bgr[2];
				}
			}
			for (int k = 0; k < channels; k++)
			{
				if (GuassionSum[k] < 0)
					GuassionSum[k] = 0;
				else if (GuassionSum[k] > 255)
					GuassionSum[k] = 255;
			}
			if (channels == 1)
				dst.at<float>(i, j) = static_cast<float>(GuassionSum[0]);
			else if (channels == 3)
			{
				Vec3f bgr = { static_cast<float>(GuassionSum[0]), static_cast<float>(GuassionSum[1]), static_cast<float>(GuassionSum[2]) };
				dst.at<Vec3f>(i, j) = bgr;
			}

		}
	}

	//竖直方向
	for (int i = kerR; i < dst.rows - kerR; i++)
	{
		for (int j = kerR; j < dst.cols - kerR; j++)
		{
			float GuassionSum[3] = { 0 };
			//滑窗搜索完成高斯核平滑
			for (int k = -kerR; k <= kerR; k++)
			{

				if (channels == 1)//如果只是单通道
				{
					GuassionSum[0] += arr[kerR + k] * dst.at<float>(i + k, j);//行变,列不换,再做竖直方向的卷积
				}
				else if (channels == 3)//如果是三通道的情况
				{
					Vec3f bgr = dst.at<Vec3f>(i + k, j);
					auto a = arr[kerR + k];
					GuassionSum[0] += a * bgr[0];
					GuassionSum[1] += a * bgr[1];
					GuassionSum[2] += a * bgr[2];
				}
			}
			for (int k = 0; k < channels; k++)
			{
				if (GuassionSum[k] < 0)
					GuassionSum[k] = 0;
				else if (GuassionSum[k] > 255)
					GuassionSum[k] = 255;
			}
			if (channels == 1)
				dst.at<float>(i, j) = static_cast<float>(GuassionSum[0]);
			else if (channels == 3)
			{
				Vec3f bgr = { static_cast<float>(GuassionSum[0]), static_cast<float>(GuassionSum[1]), static_cast<float>(GuassionSum[2]) };
				dst.at<Vec3f>(i, j) = bgr;
			}

		}
	}
	delete[] arr;
}

void harrisResponse(Mat_<float> &GaussXX, Mat_<float> &GaussYY, Mat_<float> &GaussXY, Mat_<float> &resultData, float k)
{
	//创建一张响应函数输出的矩阵
	resultData = Mat_<float>(GaussXX.size(), CV_32FC1);
	for (int i = 0; i < resultData.rows; i++)
	{
		for (int j = 0; j < resultData.cols; j++)
		{
			float a = GaussXX.at<float>(i, j);
			float b = GaussYY.at<float>(i, j);
			float c = GaussXY.at<float>(i, j);
			resultData.at<float>(i, j) = a * b - c * c - k * (a + b)*(a + b);
		}
	}
}


//非极大值抑制
void LocalMaxValue(Mat_<float> &resultData, Mat &srcGray, Mat &ResultImage, int kSize)
{
	int r = kSize / 2;
	ResultImage = srcGray.clone();
	for (int i = r; i < ResultImage.rows - r; i++)
	{
		for (int j = r; j < ResultImage.cols - r; j++)
		{
			if (resultData.at<float>(i, j) > resultData.at<float>(i - 1, j - 1) &&
				resultData.at<float>(i, j) > resultData.at<float>(i - 1, j) &&
				resultData.at<float>(i, j) > resultData.at<float>(i - 1, j - 1) &&
				resultData.at<float>(i, j) > resultData.at<float>(i - 1, j + 1) &&
				resultData.at<float>(i, j) > resultData.at<float>(i, j - 1) &&
				resultData.at<float>(i, j) > resultData.at<float>(i, j + 1) &&
				resultData.at<float>(i, j) > resultData.at<float>(i + 1, j - 1) &&
				resultData.at<float>(i, j) > resultData.at<float>(i + 1, j) &&
				resultData.at<float>(i, j) > resultData.at<float>(i + 1, j + 1))
			{
				if ((int)resultData.at<float>(i, j) > 18000)
				{
					circle(ResultImage, Point(i, j), 5, Scalar(0, 0, 255), 2, 4, 0);
				}
			}

		}
	}
}

#include <opencv2/opencv.hpp>
#include <iostream>
#include <math.h>

using namespace cv;
using namespace std;

// 定义全局变量
const string harris_winName = "自定义角点检测";
Mat src_img, gray_img;                       // src_img表示原图, gray_img表示灰度图
Mat harris_dst_img, harris_response_img;     // harris_dst_img存储自相关矩阵M的特征值和特征向量,harris_response_img存储响应函数的结果

double min_respense_value;			  // 响应函数的结果矩阵中的最小值
double max_respense_value;			  // 响应函数的结果矩阵中的最大值

int qualityValue = 30;
int max_qualityValue = 100;              // 通过qualityValue/max_qualityValue的结果作为qualitylevel来计算阈值
RNG  random_number_generator;             // 定义一个随机数发生器
void self_defining_Harris_Demo(int, void*);      //TrackBar回调函数声明

// 主函数
int main()
{
	src_img = imread("E:/C++/CVE8/img.jpg");
	if (src_img.empty())
	{
		printf("could not load the image...\n");
		return -1;
	}
	namedWindow("原图", CV_WINDOW_AUTOSIZE);
	imshow("原图", src_img);
	cvtColor(src_img, gray_img, COLOR_BGR2GRAY);      //将彩色图转化为灰度图

	// 计算特征值
	int blockSize = 3;
	int ksize = 3;
	double k = 0.04;
	harris_dst_img = Mat::zeros(src_img.size(), CV_32FC(6));
	// 目标图像harris_dst_img存储自相关矩阵M的特征值和特征向量,
	// 并将它们以(λ1, λ2, x1, y1, x2, y2)的形式存储。其中λ1, λ2是M未经过排序的特征值;
	// x1, y1是对应于λ1的特征向量;x2, y2是对应于λ2的特征向量。
	// 因此目标矩阵为6通道,即 CV_32FC(6)的矩阵。

	harris_response_img = Mat::zeros(src_img.size(), CV_32FC1);
	// harris_response_img用来存储通过每个像素值所对应的自相关矩阵所计算得到的响应值

	cornerEigenValsAndVecs(gray_img, harris_dst_img, blockSize, ksize, 4);
	// 该函数用来计算每个像素值对应的自相关矩阵的特征值和特征向量

	// 计算响应函数值
	for (int row = 0; row < harris_dst_img.rows; ++row)
	{
		for (int col = 0; col < harris_dst_img.cols; ++col)
		{
			double eigenvalue1 = harris_dst_img.at<Vec6f>(row, col)[0];     // 获取特征值1
			double eigenvalue2 = harris_dst_img.at<Vec6f>(row, col)[1];		// 获取特征值2
			harris_response_img.at<float>(row, col) = eigenvalue1 * eigenvalue2 - k * pow((eigenvalue1 + eigenvalue2), 2);
			// 通过响应公式R=λ1*λ2 - k*(λ1+λ2)*(λ1+λ2)来计算每个像素对应的响应值
		}
	}
	minMaxLoc(harris_response_img, &min_respense_value, &max_respense_value, 0, 0, Mat());   // 寻找响应矩阵中的最小值和最大值
	namedWindow(harris_winName, CV_WINDOW_AUTOSIZE);
	createTrackbar("Quality Value:", harris_winName, &qualityValue, max_qualityValue, self_defining_Harris_Demo);    //创建TrackBar
	self_defining_Harris_Demo(0, 0);

	waitKey(0);
	return 0;
}


//  回调函数实现
void self_defining_Harris_Demo(int, void*)
{
	if (qualityValue < 10)
	{
		qualityValue = 10;       // 控制qualitylevel的下限值
	}
	Mat result_img = src_img.clone();    // 输出图像
	float threshold_value = min_respense_value + (((double)qualityValue) / max_qualityValue)*(max_respense_value - min_respense_value);
	for (int row = 0; row < result_img.rows; row++)
	{
		for (int col = 0; col < result_img.cols; col++)
		{
			float resp_value = harris_response_img.at<float>(row, col);
			if (resp_value > threshold_value)
			{
				circle(result_img, Point(col, row), 2, Scalar(random_number_generator.uniform(0, 255),
					random_number_generator.uniform(0, 255), random_number_generator.uniform(0, 255)), 2, 8, 0);
			}
		}
	}
	imshow(harris_winName, result_img);
}
// CVE8.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//

#include "pch.h"
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <iostream>

using namespace cv;
using namespace std;

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

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

int main()
{
	image = imread("E:/C++/CVE8/img.jpg");
	cvtColor(image, imageGray, CV_RGB2GRAY);
	GaussianBlur(imageGray, imageGray, Size(5, 5), 1); // 滤波
	namedWindow("Corner Detected");
	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);
}
#include "pch.h"
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <math.h>
#include <iostream>

using namespace cv;
using namespace std;

int thresh = 130;
int max_count = 255;
Mat img, img_gray;
const char* output_title = "Harris Corner Dectction Result";
void Harris_Demo(int, void *);

int main(int argv, char** argc) {
	img = imread("E:/C++/CVE8/图片1.png");
	if (img.empty()) {
		printf("colud not load image...");
		return -1;
	}
	namedWindow("input image", CV_WINDOW_AUTOSIZE);
	imshow("input image", img);
	//以上是图像处理的标准开头
	namedWindow(output_title, CV_WINDOW_AUTOSIZE);
	cvtColor(img, img_gray, CV_BGR2GRAY);

	createTrackbar("Threshold", output_title, &thresh, max_count, Harris_Demo);
	Harris_Demo(0, 0);

	waitKey(0);
	return 0;
}

void Harris_Demo(int, void *) {

	Mat dst, norm_dst, normScaleDst;
	dst = Mat::zeros(img_gray.size(), CV_32FC1);
	//harris角点核心函数
	int blockSize = 3;
	int ksize = 3;
	int k = 0.04;

	cornerHarris(img_gray, dst, blockSize, ksize, k, BORDER_DEFAULT);
	//上述输出的取值范围并不是0-255 需要按照最大最小值进行归一化
	normalize(dst, norm_dst, 0, 255, NORM_MINMAX, CV_32FC1, Mat());
	convertScaleAbs(norm_dst, normScaleDst);

	Mat resultImg = img.clone();
	//用彩色来显示
	for (int row = 0; row < resultImg.rows; row++) {
		//定义每一行的指针
		uchar* currentRow = normScaleDst.ptr(row);
		for (int col = 0; col < resultImg.cols; col++) {
			int value = (int)*currentRow;
			if (value > thresh) {
				circle(resultImg, Point(col, row), 2, Scalar(0, 0, 255), 2, 8, 0);
			}
			currentRow++;
		}
	}

	imshow(output_title, resultImg);
}
2019-08-29 22:01:35 a18612039484 阅读数 99

OpenCV基础

opencv是一个计算机视觉库,用于图像处理。

常规图像处理

API:

import numpy as np
import cv2 as cv

# 读取图片
original = cv.imread('')

# 显示图片
cv.imshow('title',original)
cv.waitkey() # 不wait直接闪退,waitkey靠点击退出

# 剪裁图片
h,w = original.shape[:2]
l,t = int(w/4),int(h/4)
r,b = int(w*3/4),int(h*3/4)
cropped = original[t:b,l:r]

# 图像缩放
scaled = cv.resize(original,(int(w/4),int(h/4)),
    interpolation = cv.INTER_LINEAR)

# 图像保存
cv.imwrite('path/name',img)

边缘检测

思想:通过识别亮度梯度变化最大的像素点从而检测出物体的边缘
API:

import cv2 as cv
original = cv.imread('',cv.IMREAD_GRAYSCALE)
# 索贝尔边缘识别
s_img = cv.Sobel(origianl,cv.CV_64F,1,0,ksize = 5)
cv.imshow('sobel',s_img)
# 拉普拉斯边缘识别
l_img = cv.Laplacian(original,cv.CV_64F)
cv.imshow('Laplacian',l_img)
# 卡尼边缘识别
c_img = cv.Canny(original,50,240)
cv.imshow('canny',c_img)

cv.waitkey()

亮度提升

直方图均衡化实现亮度提升,更有利于边缘识别

import cv2 as cv
original = cv.imread('')
# 彩色图转为灰度图
gray = cv.cvtColor(original,cv.COLOR_BGR2GRAY)
# 直方图均衡化
equalized_gray = cv.equalizeHist(gray)
cv.imshow('Equalized Gray',equalized_gray)

# 提升亮度
yuv = cv.cvtColot(original,cv.COLOR_BGR2YUV)
yuv[...,0] = cv.equalizeHist(yuv[...,0])
equalize_color = cv.cvtColor(yuv,cv.COLOR_YUV2BGR)
cv.imshow('Equalized Color',equalized_color)
cv.waitKey()

角点检测

角点:平直棱线的交汇点
API:

# 灰度化原图
gray = cv.cvtColor(original,cv.COLOR_BGR2GRAY)
# Harris 角点检测器
# 边缘水平方向、垂直方向颜色值改变超过阈值7、5时即为边缘
# 边缘线方向改变超过阈值0.04弧度即为一个角点。
corners = cv.cornerHarris(gray,7,5,0.04)

# 绘制角点
mixture = original.copy()
mixture[corners>corners.max()*0.01]=[0,0,255]
cv.imshow('Corner',mixture)
cv.waitKey()

图像识别

特征点检测

常用的特征点检测方法:STAR检测和SIFT检测
特征点检测结合了边缘检测和角点检测

STAR特征点检测

API:

import cv2 as cv
original = cv.imread('')
gray = cv.cvtColor(original,cv.COLOR_BGR_2GRAY)

# 创建star特征点检测器
star = cv.xfeatures2d.StarDetector_create()
# 检测出gray图像所有的特征点
keypoints = star.detect(gray)
mixture = origianl.copy()
cv.drawKeypoints(
    original,keypoints,mixture,
    flags = cv.DRAN_MATCHES_FLAGS_DRAN_RICH_KEYPOINTS
)
cv.imshow('Mixture',mixture)
cv.waitKey()

SIFT特征点监测

API:

import cv2 as cv
original = cv.imread('')
gray = cv.cvtColor(original,cv.COLOR_BGR_2GRAY)

# 创建sift特征点检测器
star = cv.xfeatures2d.SIFT_create()
# 检测出gray图像所有的特征点
keypoints = star.detect(gray)
mixture = origianl.copy()
cv.drawKeypoints(
    original,keypoints,mixture,
    flags = cv.DRAN_MATCHES_FLAGS_DRAN_RICH_KEYPOINTS
)
cv.imshow('Mixture',mixture)
cv.waitKey()

特征值矩阵

特征值矩阵记录图像的每个特征点的梯度信息,提取特征值矩阵,只要有足够多的的样本,就可以做隐马尔可夫模型。类似于mfcc(梅尔频率倒谱系数–语音识别)
API:

import cv2 as cv
original = cv.imread('')
gray = cv.cvtColor(original,cv.COLOR_BGR2GRAY)
sift = cv.xfeatures2d.SIFT_create()
keypoints = sift.detect(gray)
# 创建特征值矩阵desc
_,desc = sift.compute(gray,keypoints)

物体识别案例

API:

import numpy as np
import cv2 as cv
import hmmlearn.hmm as hl

def search_files(directory):
    #详见CSDN ghcjasongo 工具三
    pass

# 获取训练集
train_objects = search_files('')
train_x,train_y = [],[]
for label,filenames in train_objects.items():
    descs = np.array[]
    for filename in filenames:
        image = cv.imread(filename)
        gray = cv.cvtColor(image,cv.COLOR_BGR2GRAY)
        # 范围缩放,使特征描述矩阵样本数量一致
        h,w = gray.shape[:2]
        f = 200/min(h,w)
        gray = cv.resize(gray,None,fx=f,fy=y)

        # 提取sift特征值点
        sift = cv.xfeatures2d.SIFT_create()
        keypoints = sift.detect(gray)
        _,desc = sift.compute(gray,keypoints)
        if len(descs) == 0:
            descs == desc
        else:
            descs = np.append(descs,desc,axis = 0)
    train_x.append(descs)
    train_y.appedn(label)

# 训练隐马尔可夫模型
models = {}
for descs,label in zip(train_x,train_y):
    model = hl.GaussianHMM(
        n_components = 4,covariance_type = 'diag',n_iter = 100
    )
    model[label] = model.fit(descs)


# 测试集数据整理同训练集
# 测试数据
pred_y = []
for descs in test_x:  # 遍历3次
    # 验证每个模型对当前desc的匹配度得分
    best_score, best_label = None, None
    for label, model in models.items():
        score = model.score(descs)
        if (best_score is None) or (best_score < score):
            best_score = score
            best_label = label
    pred_y.append(best_label)

print(test_y)
print(pred_y)

人脸识别

视频捕捉

API:


import cv2 as cv
# 获取视频捕捉设备
video_capture = cv.VideoCapture(0)
while True# 读取一帧
    frame = video_capture.read()[1]
    cv.imshow('VideoCapture',frame)
    if cv.waitKey(33) == 27:
        break
video_capture.release()
cv.destroyAllWindows()

人脸定位

哈尔级联人脸定位
API:

import cv2 as cv
# 哈尔级联人脸定位
face_d = cv.CascadeClassifier('haar/face.xml')
eye_d = cv.CascadeClassifier('haar/eye.xml')
nose_d = cv.CascadeClassifier('haar/nose.xml')
# 创建人脸捕捉器
vc = cv.VideoCapture(0)
while True:
    frame = vc.read[1]
    # 1.3为最小人脸尺寸,最多抓5张脸
    faces = face_d.detectMultiScale(frame,1.3,5)
    for l,t,w,h in faces:
        a,b = int(w/2),int(h/2)
        cv.ellipse(frame,(l+a,t+b),(a,b),0,0,360,(255,0,255),2)
        # 绘制椭圆
        #cv.ellipse(
        #face, # 图像
        #(l + a, t + b), # 椭圆心
        #(a, b), # 半径
        #0, # 椭圆旋转角度
        #0, 360, # 起始角, 终止角
        #(255, 0, 255), # 颜色
        #2 # 线宽
        #)
        face = frame[t:t+h,l:l+w]
        eyes = eye_d.detectMultiScale(face,1.3,5)
        for l,t,w,h in eyes:
            a,b = a,b = int(w/2),int(h/2)
            cv.ellipse(face,(l+a,t+b),(a,b),0,0,360,(255,0,255),2)
        noses = nose_d.detectMultiScale(face,1.3,5)
        for l,t,w,h in noses:
            a,b = a,b = int(w/2),int(h/2)
            cv.ellipse(face,(l+a,t+b),(a,b),0,0,360,(255,0,255),2)
    cv.imshow('VideoCapture',frame)
    if cv.waitKey(33) == 27:
        break
vc.release()
cv.destroyAllWindows()

人脸识别

opencv下的lbph(局部二值模式直方图)

import os
import numpy as np
import cv2 as cv
import sklearn.preprocessing as sp

# 创建哈尔级联人脸定位器
fd = cv.CascadeClassifier('haar/face.xml')
# 详见工具三
def search_faces(directory):
    pass

train_faces = search_faces(
    'faces/training')
# 创建标签编码器
codec = sp.LabelEncoder()
codec.fit(list(train_faces.keys()))
train_x, train_y = [], []
for label, filenames in train_faces.items():
    for filename in filenames:
        image = cv.imread(filename)
        # 对人脸进行灰度处理
        gray = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
        # 创建人脸识别
        faces = fd.detectMultiScale(gray, 1.1, 2,
                                    minSize=(100, 100))
        for l, t, w, h in faces:
            train_x.append(
                gray[t:t + h, l:l + w])
            train_y.append(
                codec.transform([label])[0])
train_y = np.array(train_y)

# !!!!!!!!!!!!!!!!!!
# 局部二值模式直方图人脸识别分类器
model = cv.face.LBPHFaceRecognizer_create()
# 训练模型
model.train(train_x, train_y)

# 整理测试集数据
test_faces = search_faces(
    'faces/testing')
test_x, test_y, test_z = [], [], []
for label, filenames in test_faces.items():
    for filename in filenames:
        image = cv.imread(filename)
        gray = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
        faces = fd.detectMultiScale(gray, 1.1, 2,
                                    minSize=(100, 100))
        for l, t, w, h in faces:
            test_x.append(
                gray[t:t + h, l:l + w])
            test_y.append(
                codec.transform([label])[0])
            a, b = int(w / 2), int(h / 2)
            # 绘制椭圆
            cv.ellipse(image, (l + a, t + b),
                       (a, b), 0, 0, 360,
                       (255, 0, 255), 2)
            test_z.append(image)
test_y = np.array(test_y)

# 预测
pred_test_y = []
for face in test_x:
    pred_code = model.predict(face)[0]
    pred_test_y.append(pred_code)

escape = False
while not escape:
    for code, pred_code, image in zip(
            test_y, pred_test_y, test_z):
        label, pred_label = \
            codec.inverse_transform([code, pred_code])
        text = '{} {} {}'.format(
            label,
            '==' if code == pred_code else '!=',
            pred_label)
        # 将文本写在图片上
        cv.putText(image, text, (10, 60),
                   cv.FONT_HERSHEY_SIMPLEX, 2,
                   (255, 255, 255), 6)
        cv.imshow('Recognizing...', image)
        if cv.waitKey(1000) == 27:
            escape = True
            break
2016-08-22 12:23:57 qq_20823641 阅读数 5178

一、角点

    图像处理和与计算机视觉领域,兴趣点(interest points),或称作关键点(keypoints)、特征点(feature points) 被大量用于解决物体识别,图像识别、图像匹配、视觉跟踪、三维重建等一系列的问题。我们不再观察整幅图,而是选择某些特殊的点,然后对他们进行局部有的放矢的分析。如果能检测到足够多的这种点,同时他们的区分度很高,并且可以精确定位稳定的特征,那么这个方法就有使用价值。

图像特征类型可以被分为如下三种:

        <1>边缘                   <2>角点 (感兴趣关键点)                 <3>斑点(Blobs)(感兴趣区域)

      其中,角点是个很特殊的存在。他们在图像中可以轻易地定位,同时,他们在人造物体场景,比如门、窗、桌等出随处可见。因为角点位于两条边缘的交点处,代表了两个边缘变化的方向上的点,,所以他们是可以精确定位的二维特征,甚至可以达到亚像素的精度。且其图像梯度有很高的变化,这种变化是可以用来帮助检测角点的。需要注意的是,角点与位于相同强度区域上的点不同,与物体轮廓上的点也不同,因为轮廓点难以在相同的其他物体上精确定位。

       角点检测(Corner Detection)是计算机视觉系统中用来获得图像特征的一种方法,广泛应用于运动检测、图像匹配、视频跟踪、三维建模和目标识别等领域中。也称为特征点检测。

       角点通常被定义为两条边的交点,更严格的说,角点的局部邻域应该具有两个不同区域的不同方向的边界。而实际应用中,大多数所谓的角点检测方法检测的是拥有特定特征的图像点,而不仅仅是“角点”。这些特征点在图像中有具体的坐标,并具有某些数学特征,如局部最大或最小灰度、某些梯度特征等。

现有的角点检测算法并不是都十分的健壮。很多方法都要求有大量的训练集和冗余数据来防止或减少错误特征的出现。另外,角点检测方法的一个很重要的评价标准是其对多幅图像中相同或相似特征的检测能力,并且能够应对光照变化、图像旋转等图像变化。

 如果某一点在任意方向的一个微小变动都会引起灰度很大的变化,那么我们就把它称之为角点

首先我们来看三幅图片理解什么是角点:

我们在图片以某像素点为中心,取一窗口,当窗口向各个方向移动时,其内部灰度值变化不是很明显,则该点即处在平坦区域(如左边图);当其内部灰度值只在几个固定的方向上变化较为明显,那么该点则处在边缘区域(如图中间部分);当向各个方向移动,其变化都是很明显,则该点为角点(如图右)


另外,关于角点的具体描述可以有几种:

  1. ·   一阶导数(即灰度的梯度)的局部最大所对应的像素点;
  2. ·  两条及两条以上边缘的交点;
  3. ·  图像中梯度值和梯度方向的变化速率都很高的点;
  4. ·  角点处的一阶导数最大,二阶导数为零,指示物体边缘变化不连续的方向。

 二、moravvec角点

           Moravec 在1981年提出Moravec角点检测算子[1],并将它应用于立体匹配。

          首先, 计算每个像素点的兴趣值, 即以该像素点为中心, 取一个w*w(如:5×5)的方形窗口, 计算0度、45度、90度、135度四个方向灰度差的平方和, 取其中的最小值作为该像素点的兴趣值.E就是像素的变化值。Moravec算子对四个方向进行加权求和来确定变化的大小,然和设定阈值,来确定到底是边还是角点。

                                                   

                                         图   以3×3为例 黑色窗口为I(x,y) 红色窗口为I(x+u,y+v)

             其中四种移位 (u,v) = (1,0), (1,1), (0,1), (-1, 1).w(x,y)为方形二值窗口,若像素点在窗口内,则取值为1, 否则为0。

      moravec角点检测步骤:

      (1)对于每一个像素点,计算在E(u,v),在我们的算法中,(u,v)的取值是((1,0), (1,1),(0,1), (-1, 1).当然,你自己可以改成(1,0),(1,1),(0,1),(-1,1),(-1,0),(-1,-1),(0,-1),(1,-1)8种情况

    (2)计算最小值对每个位置minValue = min{E(u,v)},其中(u,v) = (1,0),(1,1), (0,1), (-1, 1).

      (3)对每个位置minValue 进行判断,是不是大于设定阈值,如果是大于设定阈值,接着判断是不是局部极大,在判断角点的时候,必须判断每个方向的patch的变化。

moravec角点检测主要有两个缺点

  • 不具有旋转不变性
  • 对边缘点的反应比较强烈

<span style="font-size:18px;">#include <opencv2/opencv.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
using namespace std;
using namespace cv;
// MoravecCorners角点检测
cv::Mat MoravecCorners(cv::Mat srcImage, int kSize, int threshold)
{
	cv::Mat resMorMat = srcImage.clone();
	int r = kSize / 2;
	const int nRows = srcImage.rows;
	const int nCols = srcImage.cols;
	int nConut = 0;
	CvPoint *pPoint = new CvPoint[nRows*nCols];
	for (int i = r; i < srcImage.rows-r; i++)
	{
		for (int j = r; j<srcImage.cols-r; j++)
		{
			int wV1, wV2, wV3, wV4;
			wV1 = wV2 = wV3 = wV4 = 0;
			for (int k = -r; k < r; k++)
				wV1 += (srcImage.at<uchar>(i,j+k)-
				srcImage.at<uchar>(i,j+k+1))*(srcImage.at
				<uchar>(i,j+k)-srcImage.at<uchar>(i,j+k+1));
			for (int k = -r; k < r; k++)
				wV2 += (srcImage.at<uchar>(i+k,j)-
				srcImage.at<uchar>(i+k+1,j))*(srcImage.at
				<uchar>(i+k,j)-srcImage.at<uchar>(i+k+1,j));
			for (int k = -r; k < r; k++)
				wV3 += (srcImage.at<uchar>(i+k,j+k)-
				srcImage.at<uchar>(i+k+1,j+k+1))*(srcImage.at
				<uchar>(i+k,j+k)-srcImage.at<uchar>(i+k+1,j+k+1));
			for (int k = -r; k < r; k++)
				wV4 += (srcImage.at<uchar>(i+k,j-k)-
				srcImage.at<uchar>(i+k+1,j-k-1))*(srcImage.at
				<uchar>(i+k,j-k)-srcImage.at<uchar>(i+k+1,j-k-1));
			int value = min(min(wV1,wV2), min(wV3,wV4));
			if (value > threshold)
			{
				pPoint[nConut] = cvPoint(j,i);
				nConut++;
			}
		}
	}
	for (int i = 0; i < nConut; i++)
		cv::circle(resMorMat, pPoint[i], 5, cv::Scalar(255,0,0));
    return resMorMat;
}

int main()
{
	cv::Mat srcImage = imread("lena.jpg",0);
	if (!srcImage.data)
		return -1;
	cv::Mat resMorMat =  MoravecCorners(srcImage, 5,10000);
	cv::imshow("srcImage", srcImage);
	cv::imshow("resMorMat",resMorMat);
	cv::waitKey(0);
	return 0;
}</span>

 


 三、harris角点检测

       在harris的角点检测中,使用的是高斯窗口,所以w(x,y)表示的是高斯窗口中的权重。此时 当u和v取两组相互垂直的值时,E(u,v)都有较大值的点。

 

<1>计算图像I(x,y)在x和y两个方向的梯度Ix,Iy  


     

结果解释:

  • 角点:最直观的印象就是在水平、竖直两个方向上变化均较大的点,即Ix、Iy都较大
  • 边缘:仅在水平、或者仅在竖直方向有较大的变化量,即Ix和Iy只有其一较大
  • 平坦地区:在水平、竖直方向的变化量均较小,即Ix、Iy都较小
或者说,r1 r2是特征值
  • r1,r2都很小,对应于图像中的平滑区域
  • r1,r2都很大,对应于图像中的角点
  • r1,r2一个很大,一个很小,对应于图像中的边缘

      Harris角点检测最直观的解释是:在任意两个相互垂直的方向上,都有较大变化的点。

      在moravec角点检测中,w(x,y)的取值是二元的,在窗口内部就取值为1,在窗口外部就取值为0,在harris的角点检测中,使用的是高斯窗口,所以w(x,y)表示的是高斯窗口中的权重。此时 当u和v取两组相互垂直的值时,E(u,v)都有较大值的。

Harris角点检测算法有诸多优点:

  •  旋转不变性,椭圆转过一定角度但是其形状保持不变(特征值保持不变)
  •  对于图像灰度的仿射变化具有部分的不变性,由于仅仅使用了图像的一介导数,对于图像灰度平移变化不变;对于图像灰度尺度变化不变

当然Harris也有许多不完善的地方:它对尺度很敏感,不具备几何尺度不变性。

 

      cornerHarris 函数用于在OpenCV中运行Harris角点检测算子处理图像。和cornerMinEigenVal( )以及cornerEigenValsAndVecs( )函数类似,cornerHarris 函数对于每一个像素(x,y)在邻域内,计算2x2梯度的协方差矩阵,接着它计算如下式子:

 

即可以找出输出图中的局部最大值,即找出了角点。

<span style="font-size:18px;">C++: void cornerHarris(InputArraysrc,OutputArray dst, int blockSize, int ksize, double k,intborderType=BORDER_DEFAULT )
第一个参数,InputArray类型的src,输入图像,即源图像,填Mat类的对象即可,且需为单通道8位或者浮点型图像。
第二个参数,OutputArray类型的dst,函数调用后的运算结果存在这里,即这个参数用于存放Harris角点检测的输出结果,和源图片有一样的尺寸和类型。
第三个参数,int类型的blockSize,表示邻域的大小,更多的详细信息在cornerEigenValsAndVecs()中有讲到。
第四个参数,int类型的ksize,表示Sobel()算子的孔径大小。
第五个参数,double类型的k,Harris参数。
第六个参数,int类型的borderType,图像像素的边界模式,注意它有默认值BORDER_DEFAULT。更详细的解释,参考borderInterpolate( )函数。</span>

<span style="font-size:18px;">#include <opencv2/opencv.hpp>  
#include "opencv2/highgui/highgui.hpp"  
#include "opencv2/imgproc/imgproc.hpp"  
using namespace cv;  
using namespace std;  
#define WINDOW_NAME1 "【程序窗口1】"        
#define WINDOW_NAME2 "【程序窗口2】"        
Mat g_srcImage, g_srcImage1,g_grayImage;  
int thresh = 30; 
int max_thresh = 175; 
void on_CornerHarris( int, void* );//回调函数  
int main( int argc, char** argv )  
{  
    g_srcImage = imread( "1.jpg", 1 );  
     if(!g_srcImage.data ) { printf("读取图片错误,请确定目录下是否有imread函数指定的图片存在~! \n"); return false; }    
     imshow("原始图",g_srcImage);  
    g_srcImage1=g_srcImage.clone( );  
    cvtColor( g_srcImage1, g_grayImage, CV_BGR2GRAY );  
    namedWindow( WINDOW_NAME1, CV_WINDOW_AUTOSIZE );  
    createTrackbar( "阈值: ", WINDOW_NAME1, &thresh, max_thresh, on_CornerHarris );  
    on_CornerHarris( 0, 0 );  
    waitKey(0);  
    return(0);  
}  
  
void on_CornerHarris( int, void* )  
{  

    Mat dstImage;
    Mat normImage;
    Mat scaledImage;
    dstImage = Mat::zeros( g_srcImage.size(), CV_32FC1 );  
    g_srcImage1=g_srcImage.clone( );  
    cornerHarris( g_grayImage, dstImage, 2, 3, 0.04, BORDER_DEFAULT );  
    normalize( dstImage, normImage, 0, 255, NORM_MINMAX, CV_32FC1, Mat() );  
    convertScaleAbs( normImage, scaledImage );
    for( int j = 0; j < normImage.rows ; j++ )  
    { for( int i = 0; i < normImage.cols; i++ )  
    {  
        if( (int) normImage.at<float>(j,i) > thresh+80 )  
        {  
            circle( g_srcImage1, Point( i, j ), 5,  Scalar(10,10,255), 2, 8, 0 );  
            circle( scaledImage, Point( i, j ), 5,  Scalar(0,10,255), 2, 8, 0 );  
        }  
    }  
    }  
    imshow( WINDOW_NAME1, g_srcImage1 );  
    imshow( WINDOW_NAME2, scaledImage );  
  
}  
</span>


 

四、知识补充

     关于矩阵知识的一点补充:好长时间没看过线性代数的话,这一段比较难理解。可以看到M是实对称矩阵,这里简单温习一下实对称矩阵和二次型的一些知识点吧。

1. 关于特征值和特征向量: 

特征值的特征向量的概念忘了就自己查吧,这里只说关键的。对于实对称矩阵M(设阶数为n),则一定有n个实特征值,每个特征值对应一组特征向量(这组向量中所有向量共线),不同特征值对应的特征向量间相互正交;(注意这里说的是实对称矩阵,不是所有的矩阵都满足这些条件)

2. 关于对角化:

对角化是指存在一个正交矩阵Q,使得  Q’MQ 能成为一个对角阵(只有对角元素非0),其中Q’是Q的转置(同时也是Q的逆,因为正交矩阵的转置就是其逆)。一个矩阵对角化后得到新矩阵的行列式和矩阵的迹(对角元素之和)均与原矩阵相同。如果M是n阶实对称矩阵,则Q中的第 j 列就是第 j 个特征值对应的一个特征向量(不同列的特征向量两两正交)。

3. 关于二次型:

对于一个n元二次多项式,f(x1,x2....xn)= ∑ ( aij*xi*xj ) ,其中 i 和 j 的求和区间均为 [1,n] ,

可将其各次的系数 aij 写成一个n*n矩阵M,由于 aij 和 aji 的对称等价关系,一般将 aij 和 aji 设为一样的值,均为xi*xj 的系数的二分之一。这样,矩阵M就是实对称矩阵了。即二次型的矩阵默认都是实对称矩阵

4. 关于二次型的标准化(正交变换法):

二次型的标准化是指通过构造一个n阶可逆矩阵 C,使得向量 ( x1,x2...xn ) = C * (y1,y2...yn),把n维向量 x 变换成n维向量 y ,并代入f(x1,x2....xn) 后得到 g(y1,y2...yn),而后者的表达式中的二次项中不包含任何交叉二次项 yi*yj(全部都是平方项 yi^2),也即表达式g的二次型矩阵N是对角阵。用公式表示一下 f 和 g ,(下面的表达式中 x 和 y都代表向量,x' 和 y' 代表转置)

f = x' * M * x ;

g = f = x' * M * x = (Cy)' * M * (Cy) = y'* (C'MC) * y = y' * N * y  ;

因此 C‘MC = N。正交变换法,就是直接将M对角化得到N,而N中对角线的元素就是M的特征值。正交变换法中得到的 C 正好是一个正交矩阵,其每一列都是两两正交的单位向量,因此 C 的作用仅仅是将坐标轴旋转(不会有放缩)。

http://blog.csdn.net/newthinker_wei/article/details/45603583

http://www.360doc.com/content/15/1212/23/20007814_519967668.shtml

http://blog.csdn.net/crzy_sparrow/article/details/7391511

http://blog.csdn.net/lu597203933/article/details/15088485

http://blog.csdn.net/poem_qianmo/article/details/29356187


图像识别算法交流 QQ群:145076161,欢迎图像识别与图像算法,共同学习与交流

2016-03-15 15:28:41 qq_14821023 阅读数 5626

在本周的计算机视觉与模式识别作业中,给定输入图像是两张普通A4打印纸,上面可能有手写笔记或者打印内容但是拍照时角度不正。要求输出:
1. 图像的边缘;
2. 计算 A4纸边缘的各直线方程;
3. 提取A4纸的4个角点。
这里写图片描述
这里写图片描述


作业要求的是使用C++的CImg库,与OpenCV和MATLAB相比,CImg还是有点不太方便。

  • 图像边缘的提取
    这个任务比较简单,找一个像Prewitt算子、拉普拉斯算子或是其他算子,对输入图像进行卷积计算就可以了,原理是计算图像像素的变化,而变化剧烈的地方就是图像边缘出现的地方。
    这里我使用了Sobel算子来进行边缘提取并二值化。
    int sobelX[3][3] = { { -1,0,1 },{ -2,0,2 },{ -1,0,1 } };
    int sobelY[3][3] = { { 1,2,1 },{ 0,0,0 },{ -1,-2,-1 } };
    Sobel算子就是分别求图像X、Y方向的梯度,再求梯度平方和,小于某个阈值置为0,反之为255。
    最后结果如下:
    图1
    图2
  • Harris角点的检测
    主要时间都花费在这部分的实现上面了。这里推荐一篇文章,里面详解了Harris角点检测算法的原理、步骤,还给出了OpenCV实现的代码。
    http://www.cnblogs.com/ronny/p/4009425.html
    还有一些有帮助的文章:
    http://blog.csdn.net/l_inyi/article/details/8915116
    http://blog.csdn.net/berguiliu/article/details/24985825
    关键步骤如下:
    这里写图片描述
    而我自己则是在理解该文章所写算法的基础上,使用CImg库按步骤地搜索网上的资料,完成每一部分之后,检测出了角点。
    在实现过程中遇到的问题主要有:
    • 在求x、y方向梯度的时候,一开始使用的是一位的模板,导致检测到的角点过多,容易混淆。后来经过实验,改成了二维的3X3模板,效果较好;
    • 在第三步使用高斯加权的时候,我以为是对各梯度矩阵使用高斯分布的公式:这里写图片描述
      角点检测的结果:这里写图片描述
      明显不太对。
      后来查阅资料才认识到应该使用高斯低通滤波窗口对Ixx、Iyy、Ixy进行卷积。更正之后的结果:
      这里写图片描述这里写图片描述
  • 求边缘直线方程
    在Harris角点检测所得的图片当中,可以看到四个角的角点是独立出来的,能够比较容易地提取它们的坐标。因此可以利用Harris角点检测过程中计算的result数组来找到四个顶点的坐标,有了坐标之后四条直线的方程也就知道了。
    最后求得的方程:
    角点和方程

差点忘了贴代码:

// CV_HW2.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include <iostream>
#include <algorithm>
#include "CImg.h"
using namespace cimg_library;
using namespace std;

/*

CImg<> edgeDetect(CImg<> src) {
    int mask[3][3] = { {0,-1,0},{-1,4,-1},{0,-1,0} };
    int width = src.width();
    int height = src.height();
    CImg<double> newImg(width, height, 1, 3);
    newImg.fill(0);
    for (int i = 1; i < width - 1; i++) {
        for (int j = 1; j < height - 1; j++) {

            for (int k = 0; k < 3; k++) {
                double sum = 0.0;
                for (int row = 0; row < 3; row++) {
                    for (int col = 0; col < 3; col++) {
                        sum += mask[row][col] * src(i - 1 + row, j - 1 + col, k);
                    }
                }
                if (sum < 0)
                    sum = 0;
                else if (sum > 255)
                    sum = 255;
                newImg(i, j, k) = sum;
            }
        }
    }
    newImg.display();
    return newImg;
}
*/

CImg<> edgeDetect(CImg<> src) {
    int sobelX[3][3] = { { -1,0,1 },{ -2,0,2 },{ -1,0,1 } };
    int sobelY[3][3] = { { 1,2,1 },{ 0,0,0 },{ -1,-2,-1 } };
    int width = src.width();
    int height = src.height();
    CImg<double> G(width, height, 1, 3);
    CImg<double> newImg(width, height, 1, 3);
    newImg.fill(0);
    for (int i = 1; i < width - 1; i++) {
        for (int j = 1; j < height - 1; j++) {
            double sumX = 0.0;
            double sumY = 0.0;
            for (int row = 0; row < 3; row++) {
                for (int col = 0; col < 3; col++) {
                    sumX += sobelX[row][col] * src(i - 1 + row, j - 1 + col, 0);
                    sumY += sobelY[row][col] * src(i - 1 + row, j - 1 + col, 0);
                }
            }

            newImg(i, j, 0) = sqrt(sumX*sumX + sumY*sumY) > 85 ? 255 : 0;
            newImg(i, j, 1) = newImg(i, j, 0);
            newImg(i, j, 2) = newImg(i, j, 0);
        }
    }
    newImg.display();
    //newImg.dilate(3.5);
    return newImg;
}

void print(double k, double b) {
    cout << "y = " << k << " * " << "x";
    if (b > 0) {
        cout << "+" << b;
    }
    else if (b < 0) {
        cout << b;
    }
    cout << endl;
}

void line(double x1, double y1, double x2, double y2) {
    double k = (y1 - y2) / (x1 - x2);
    double b =y1 - k*x1;
    print(k, b);
}

CImg<> pointDetect(CImg<> src) {
    int width = src.width();
    int height = src.height();

    double upLeftX = 0.0;
    double upLeftY = 0.0;
    double downLeftX = 0.0;
    double downLeftY = 0.0;
    double upRightX = 0.0;
    double upRightY = 0.0;
    double downRightX = 0.0;
    double downRightY = 0.0;
    //求左上顶点

    for (int j = 100; j < height - 50; j++) {
        for (int i = 100; i < width; i++) {

            if (upLeftX == 0 && upLeftY == 0 && src(i, j, 0) == 1) {
                upLeftX = i;
                upLeftY = j;
                break;
            }
        }
        if (upLeftX || upLeftY) {
            break;
        }
    }
    //求右上顶点

    for (int i = width - 50; i >= 0; i--) {
        for (int j = 100; j < height; j++) {
            if (upRightX == 0 && upRightY == 0 && src(i, j,0) == 1) {
                upRightX = i;
                upRightY = j;
                break;
            }
        }
        if (upRightX || upRightY) {
            break;
        }
    }
    //求左下顶点

    for (int j = height - 110; j >= 0; j--) {
        for (int i = width-100; i >= 0; i--) {
            if (downLeftX == 0 && downLeftY == 0 && src(i, j, 0) == 1) {
                downLeftX = i;
                downLeftY = j;
                break;
            }
        }
        if (downLeftX || downLeftY) {
            break;
        }
    }
    //求右下顶点
    for (int j = downLeftY-1; j >= 0; j--) {
        for (int i = width - 50; i >= 0; i--) {
            if (downRightX == 0 && downRightY == 0 && src(i, j, 0) == 1
                && abs((i - downLeftX) - (upRightX - upLeftX)) < 60) {
                downRightX = i;
                downRightY = j;
                break;
            }
        }
        if (downRightX && downRightY) {
            break;
        }
    }

    cout << "图像的四个角点坐标分别是:" << endl;
    cout << "左上角upLeft:(" << upLeftX << "," << upLeftY  << ")"<< endl;
    cout << "右上角upRight:(" << upRightX << " ," << upRightY << ")" << endl;
    cout << "右下角downRight:(" << downRightX << " ," << downRightY << ")" << endl;
    cout << "左下角downLeft:(" << downLeftX << " ," << downLeftY << ")" << endl;

    cout << "上方直线方程为:";
    line(upLeftX, upLeftY, upRightX, upRightY);

    cout << "左方直线方程为:";
    line(upLeftX, upLeftY, downLeftX, downLeftY);

    cout << "右方直线方程为:";
    line(upRightX, upRightY, downRightX, downRightY);

    cout << "下方直线方程为:";
    line(downLeftX, downLeftY, downRightX, downRightY);
    //src.dilate(3);
    //src.display();
    return src;
}


CImg<> harrisCorners(CImg<> src) {
    int width = src.width();
    int height = src.height();

    int horizontal[3][3] = { { 5,0, -5 },{ 8,0,-8 }, { 5,0,-5 } };
    int vertical[3][3] = { { 5,8, 5 },{ 0,0,0 },{ -5,-8,-5 } };

    CImg<double> Ix(width, height, 1, 3);
    CImg<double> Iy(width, height, 1, 3);

    CImg<double> Ixx(width, height, 1, 3);
    CImg<double> Iyy(width, height, 1, 3);
    CImg<double> Ixy(width, height, 1, 3);

    //x、y方向梯度
    for (int i = 1; i < width; i++) {
        for (int j = 1; j < height - 1; j++) {
            for (int k = 0; k < 3; k++) {
                double sumX = 0.0;
                double sumY = 0.0;
                for (int row = 0; row < 3; row++) {
                    for (int col = 0; col < 3; col++) {
                        sumX += src(row + i - 1, col + j - 1, k)*horizontal[row][col];
                        sumY += src(row + i - 1, col + j - 1, k)*vertical[row][col];
                    }
                }
                Ix(i, j, k) = sumX;
                Iy(i, j, k) = sumY;
            }
        }
    }

    //两个方向梯度的乘积
    for (int i = 0; i < width; i++) {
        for (int j = 0; j < height; j++) {
            for (int k = 0; k < 3; k++) {
                Ixx(i, j, k) = Ix(i, j, k)*Ix(i, j, k);
                Iyy(i, j, k) = Iy(i, j, k)*Iy(i, j, k);
                Ixy(i, j, k) = Ix(i, j, k)*Iy(i, j, k);
            }
        }
    }
    //高斯加权
    CImg<double> R(width, height, 1, 3);
    double pi = 3.1415926;
    double alpha = 0.04;
    double maxR[3] = { INT_MIN };
    double threshHold[3];
    int gauss[3][3] = { {1,2,1},{2,4,2},{1,2,1} };
    for (int i = 1; i < width-1; i++) {
        for (int j = 1; j < height-1; j++) {
            for (int k = 0; k < 3; k++) {
                double A = 0.0;
                double B = 0.0;
                double C = 0.0;
                for (int row = 0; row < 3; row++) {
                    for (int col = 0; col < 3; col++) {
                        A += Ixx(row + i - 1, col + j - 1, k)*gauss[row][col];
                        B += Ixy(row + i - 1, col + j - 1, k)*gauss[row][col];
                        C += Iyy(row + i - 1, col + j - 1, k)*gauss[row][col];
                    }
                }
                A /= 16.0;
                B /= 16.0;
                C /= 16.0;
                //double A = (Ixx(i - 1, j - 1, k) + Ixx(i - 1, j + 1, k) + Ixx(i + 1, j - 1, k) + Ixx(i + 1, j + 1, k) + (Ixx(i - 1, j, k) + Ixx(i, j - 1, k) + Ixx(i + 1, j, k) + Ixx(i, j + 1, k)) * 2.0 + Ixx(i, j, k) * 4) / 16;
                //double A = exp(-Ixx(i, j, k)*Ixx(i, j, k) / 8.0) / (2.0 * sqrt(2.0 * pi));
                //double B = exp(-Ixy(i, j, k)*Ixy(i, j, k) / 8.0) / (2.0 * sqrt(2.0 * pi));
                //double C = exp(-Iyy(i, j, k)*Iyy(i, j, k) / 8.0) / (2.0 * sqrt(2.0 * pi));

                //Harris响应值R
                R(i, j, k) = (A*C - B*B) - alpha*(A + C)*(A + C);
                if (maxR[k] < R(i, j, k)) {
                    maxR[k] = R(i, j, k);
                    threshHold[k] = 0.01*maxR[k];
                }
            }
        }
    }
    //小于阈值的R置为0
    for (int i = 0; i < width; i++) {
        for (int j = 0; j < height; j++) {
            for (int k = 0; k < 3; k++) {
                if (R(i, j, k) < threshHold[k]) {
                    R(i, j, k) = 0;
                }
            }
        }
    }
    CImg<double> result(width, height, 1, 3);
    result.fill(0);
    for (int i = 1; i < width - 1; i++) {
        for (int j = 1; j < height - 1; j++) {
            for (int k = 0; k < 3; k++) {
                //局部非最大值抑制
                if (R(i, j, k) > threshHold[k] && R(i, j, k) > R(i - 1, j - 1, k) && R(i, j, k) > R(i - 1, j, k) && R(i, j, k) > R(i - 1, j + 1, k)
                    && R(i, j, k) > R(i, j - 1, k) && R(i, j, k) > R(i, j + 1, k) && R(i, j, k) > R(i + 1, j - 1, k) && R(i, j, k) > R(i + 1, j, k)
                    && R(i, j, k) > R(i + 1, j + 1, k)) {
                    result(i, j, k) = 1;
                }
            }
        }
    }
    pointDetect(result);
    unsigned char color[] = { 255,255,0 };
    for (int i = 0; i < width; i++) {
        for (int j = 0; j < height; j++) {
            if (result(i, j, 0) == 1) {
                src.draw_circle(i, j, 5, color);
            }
        }
    }
    src.display();
    return src;
}

int main()
{
    CImg<> src1("1.bmp");
    CImg<> tmp = edgeDetect(src1);
    tmp.save_bmp("edge1.bmp");
    CImg<> tmp2 = harrisCorners(src1);
    tmp2.save_bmp("harris1.bmp");

    CImg<> src("2.bmp");
    CImg<> tmp3 = edgeDetect(src);
    tmp3.save_bmp("edge2.bmp");
    CImg<> tmp5 = harrisCorners(src);
    tmp5.save_bmp("harris2.bmp");
    return 0;
}

opencv图像角点检测

阅读数 5052

没有更多推荐了,返回首页