Harris算子以及未来的规划...

Harris算子以及未来规划...

找工作的时候，被问到了harris，面试官问梯度矩阵特征值在边缘区域的分布情况。当时答的感觉也没啥问题，但面试官就是在我说出答案后特别诧异.....回来后挺难过.....后来又问了相机标定的流程，我只敢写个相机模型，没怎么答上来相机标定怎么计算。。感觉coding还是不够深入，找到工作后，发现自己很多不足的地方。打算赶紧抢救下。。

1.Harris 实现，这个应该不难。

2.多尺度的harris ，Hessian，DoG，HarrisLaplace，HessianLaplace，看看论文和VLFeat 实现下。

3.相机标定，在看机器人视觉测量与控制这本书，近几天不忙应该能coding一个basic相机标定的算法。

4.再把三大特征看完 HoG，Harr，LBP

4.1看下HDR的论文，实现一下算法过程，此外把scaling 图像缩放看一下

5.项目上面的话，线结构光的可以改进下激光平面标定的算法，用PNP算一下

6.把PNP的视频看完，还有几篇论文。

7.做一下知识管理，找工作后发现好的知识管理和高效率很重要。

8.未来：关注物体检测（FastRCNN和YOLO）人脸识别，以及人脸特征点。以及GIMP的源码。写点linux上跑的图像处理软件。

9.学习层面的话，先把操作系统，算法设计看完，计算机网络也恶补一下，开始入嵌入式操作系统的坑，最好是把安卓移植过程了解清楚，以及安卓的一些坑。

10.再长远的话学java，入安卓的坑了，写一点在安卓上跑的app。

以下是Harris C++源码：

参考过https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harris.cpp

在处理非极大抑制上做的蛮好。

#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <math.h>
using namespace std;
using namespace cv;
class CmpHarris
{
public:
	CmpHarris(const Mat& img, int threshold = 0.01, float _lamada = 0.04) : thresh(threshold), lamada(_lamada) {
		if (img.channels() == 3)
		{
			cvtColor(img, image, CV_RGB2GRAY);
			image.convertTo(image, CV_32FC1);
		}
		else
		{
			image = img.clone();
			image.convertTo(image, CV_32FC1);
		}
		gradx2 = cv::Mat::zeros(img.size(), CV_32FC1);
		grady2 = cv::Mat::zeros(img.size(), CV_32FC1);
		gradxy = cv::Mat::zeros(img.size(), CV_32FC1);

	};
	
	void CmpHarris::getHarrisRespose(vector<Point2d> &response,Mat& res);
	~CmpHarris()
	{
		
	}
private:
	Mat image;
	Mat gradx2;
	Mat grady2;
	Mat gradxy;
	float thresh;
	float lamada;
	int cmpSobel(Mat& image, Mat& gardx2, Mat& gardy2, Mat& gardxy);
	int Harris(vector<Point2d> &response,Mat& res, Mat& gradx2, Mat& grady2, Mat & gradxy, float thresh, float lamada);
};
int CmpHarris::cmpSobel(Mat& image, Mat& gardx2, Mat& gardy2, Mat& gradxy)
{
	int width = image.cols;
	int height = image.rows;
	Mat gradX, gradY;
	gradX = Mat::zeros(gradx2.size(), CV_32FC1);
	gradY = Mat::zeros(grady2.size(), CV_32FC1);
	for (int i = 1; i < height - 1; i++)
	{
		float* data = image.ptr<float>(i);
		float* gradxd = gradX.ptr<float>(i);
		float* gradyd = gradY.ptr<float>(i);
		float* updata = image.ptr<float>(i - 1);
		float* downdata = image.ptr<float>(i + 1);
		for (int j = 1; j < width - 1; j++)
		{
			float gradx, grady;
			float a00, a01, a02, a10, a12, a20, a21, a22;
			a00 = updata[j - 1]; a01 = updata[j]; a02 = updata[j + 1];
			a10 = data[j - 1]; a12 = data[j + 1];
			a20 = downdata[j - 1]; a21 = downdata[j]; a22 = downdata[j + 1];
			/*gradx = -a00 - a10 * 2 - a20 + a02 + 2 * a12 + a22;
			grady = -a00 - a01 * 2 - a02 + a20 + 2 * a21 + a22;*/
			gradx = -a00 - a10 * 2 - a20 + a02 + 2 * a12 + a22;
			grady = -a00 - a01 * 2 - a02 + a20 + 2 * a21 + a22;
			gradxd[j] = gradx;
			gradyd[j] = grady;
		}
	}
	gradx2 = gradX.mul(gradX);
	grady2 = gradY.mul(gradY);
	gradxy = gradX.mul(gradY);
	return 1;
}
int CmpHarris::Harris(vector<Point2d> &response,Mat& res, Mat& gradx2, Mat& grady2, Mat & gradxy,float thresh, float lamada)
{
	//首先进行行卷积

	Mat gaussKernel = getGaussianKernel(7, 2);
	Mat HarrisResponse;
	HarrisResponse = Mat::zeros(gradx2.size(),CV_32FC1);
	filter2D(gradx2, gradx2, CV_32FC1, gaussKernel);
	filter2D(grady2, grady2, CV_32FC1, gaussKernel);
	filter2D(gradxy, gradxy, CV_32FC1, gaussKernel);
	for (int i = 0; i < gradx2.rows; i++)
	{
		float* dataRes = HarrisResponse.ptr<float>(i);
		for (int j = 0; j <  gradx2.cols; j++)
		{
			float sumGradx2, sumGrady2, sumGradxy;
			sumGradx2 = sumGrady2 = sumGradxy = 0;
			sumGradx2 = gradx2.at<float>(i, j);
			sumGrady2 = grady2.at<float>(i, j);
			sumGradxy = gradxy.at<float>(i, j);
			float det_m = sumGradx2 *sumGrady2 - sumGradxy  * sumGradxy;
			float trace_m = sumGradx2 + sumGrady2;
			dataRes[j] = det_m - lamada * trace_m * trace_m;
		}
	}
	//NMS +　aaaaaaaaaaaaaa
	double maxResponse = 0;
	minMaxLoc(HarrisResponse, NULL, &maxResponse,NULL,NULL);
	Mat dilated;
	Mat localMax;

	//膨胀就是求局部最大值的操作，核B与图形卷积，
	//即计算核B覆盖的区域的像素点的最大值，
	//并把这个最大值赋值给参考点指定的像素
	dilate(HarrisResponse, dilated, Mat());

	compare(HarrisResponse, dilated, localMax, CMP_EQ);

	Mat cornerMap;
	double qualityLevel = 0.01;
	double threshs = qualityLevel * maxResponse;
	cornerMap = HarrisResponse > threshs;
	bitwise_and(cornerMap, localMax, res);
	for (int i = 0; i < res.rows; i++)
	{
		uchar* p = res.ptr<uchar>(i);
		for (int j = 0; j < res.cols; j++)
		{
			if (p[j])
			{
				response.push_back(Point2d(i, j));
			}
		}
	 }
	/*imshow("cv1", dilated);
	imshow("c2", HarrisResponse);
	imshow("cv3", localMax);
	imshow("cv4", res);
	waitKey(0);
	*/

	return 1;
}

void CmpHarris::getHarrisRespose(vector<Point2d> &response,Mat& res)
{
	cmpSobel(this->image, this->gradx2, this->grady2, this->gradxy);
	Harris(response,res, this->gradx2, this->grady2, this->gradxy, this->thresh, this->lamada);
}

int main(int argc, char **argv) {
	cv::Mat im;
	Mat image, res, res_norm, res_norm_scale;
	vector<Point2d> HarrisRespose;
	im = cv::imread("D://2.jpg");
	CmpHarris harris(im, 0.01, 0.04);
	harris.getHarrisRespose(HarrisRespose,res);
	if (!HarrisRespose.empty())
	{
		for (int i = 0; i < HarrisRespose.size(); i++)
		{
			Point FeatPoint(HarrisRespose[i].y, HarrisRespose[i].x);
			circle(im, FeatPoint, 1, CV_RGB(0, 255, 0), 2, 0, 0);
		}
	}
	namedWindow("harris");
	imshow("harris", im);
	waitKey(0);
	return 0;
}

上图是C++版本生成harris角点

下面是matlab版本的：

% in_image-待检测的rgb图像数组
% a--角点参数响应，取值范围：0.04~0.06
% [posr，posc]-角点坐标
in_image = imread('D:\\2.jpg');
a = 0.04;
in_image=rgb2gray(in_image);
I=double(in_image);
%%%%计算xy方向梯度%%%%%

fx=[-1,0,1];%x方向梯度模板
Ix=filter2(fx,I);%x方向滤波
fy=[-1;0;1];%y方向梯度模板(注意是分号)
Iy=filter2(fy,I);
%%%%计算两个方向梯度的乘积%%%%%
Ix2=Ix.^2;
Iy2=Iy.^2;
Ixy=Ix.*Iy;
%%%%使用高斯加权函数对梯度乘积进行加权%%%%
%产生一个7*7的高斯窗函数，sigma值为2
h=fspecial('gaussian',[7,7],2);
IX2=filter2(h,Ix2);
IY2=filter2(h,Iy2);
IXY=filter2(h,Ixy);
%%%%%计算每个像元的Harris响应值%%%%%
[height,width]=size(I);
R=zeros(height,width);
%像素(i,j)处的Harris响应值
for i=1:height
    for j=1:width
        M=[IX2(i,j) IXY(i,j);IXY(i,j) IY2(i,j)];
        R(i,j)=det(M)-a*(trace(M))^2;
    end
end

%%%%%去掉小阈值的Harris值%%%%%
Rmax=max(max(R));
%阈值
t=0.01*Rmax;
for i=1:height
    for j=1:width
        if R(i,j)<t
            R(i,j)=0;
        end
    end
end
figure
imshow(R);
hold on
%%%%%进行3*3领域非极大值抑制%%%%%%%%%
corner_peaks=imregionalmax(R);
%imregionalmax对二维图片，采用8领域（默认，也可指定）查找极值，三维图片采用26领域
%极值置为1，其余置为0
num=sum(sum(corner_peaks));
%%%%%%显示所提取的Harris角点%%%%
[posr,posc]=find(corner_peaks==1);
figure
imshow(in_image);
hold on
for i=1:length(posr)
    plot(posc(i),posr(i),'r+');
end

上图是matlab生成的harris角点