1. 理论
角点是一个局部概念。窗口在局部范围移动,如果采集到的灰度值没有显著变化,那么该区域不存在角点;如果在某一方向上有显著变化,那么该区域可能存在一条直线;如果在各个方向都有变化,那么该区域可能存在角点。
Harris以当前待检点为原点,调整x和y方向的偏移值获取相邻像素与该像素的相似度:,其中w(u,v)反映了每个像素的权重,Harris使用高斯分布对其赋值。
根据泰勒公式,有,其中I_x和I_y分别为x和y方向的偏导。这样,相似度公式可以整理为,
其中。显然,这是个二次项函数,而二次项函数从几何学上来说是一个椭圆函数,由M(x,y)决定的椭圆函数的两个特征值反映了其扁率和尺寸。Harris算法就是根据这两个特征值来确定当前点是否为角点的(图1-1)。
图1-1. 椭圆特征值与图像特征的关系
但是,计算矩阵特征值比较麻烦,在这里用矩阵M的行列式和迹来反映特征值的情况:,其中,,显然,当两个特征值较大时R取较大值。
至此,可以归纳出Harris检测角点的算法步骤:
(1)计算图像在X和Y方向上的梯度
(2)计算图像两个方向的点乘结果
(3)高斯函数对其进行加权,得到矩阵M
(4)计算每个像素的Harris响应值R,使用阈值t的置为0
(5)在3*3或5*5邻域内进行非最大值抑制(取窗口内最大值对应的点作为角点)
2. Harris角点的性质
(1)由于微分算子的特性(邻域灰度值分布特点),图像整体上的亮度、对比度变化并不会改变Harris角点的位置,但由于阈值 t 的关系,可能会影响Harris角点的检测数量
(2)由于二阶矩阵的特性(椭圆方向变化并不影响特征值本身),Harris角点具有旋转不变性
(3)因为窗口尺寸是一定的,因此Harris角点不具备尺度不变性(角点位置放大可能被检测为边缘或曲线)
3. 代码
#define PI 3.14159
#define MAXPOINTS 50000
typedef struct myPoint
{
int x;
int y;
}_point;
typedef struct mySize
{
int height;
int width;
}_size;
template<typename T1, typename T2>
void convolution(T1* img, T2** out, _size imgSize, T2* kernel, _size kerSize);
void dotMultiply(int* src1, int* src2, int** dst, _size srcSize);
void gaussian(int* src, float** out, _size imgSize, _size gaussSize, float sita);
void strengthR(float* A, float* B, float* C, float** R, float* maxR, _size imgSize, float alpha);
void depressR(float* srcR, _point* dstR, int* numDst, _size imgSize, _size depressWin, float threshold);
void _cornerHarris(unsigned char* src, _size srcSize, _point* corners, int* numCorner, float alpha);
void main()
{
cv::Mat img = cv::imread("../file/butterfly.jpg", 0);
unsigned char* imgBuffer = new unsigned char[img.rows*img.cols];
for(int i=0; i<img.rows; i++)
{
uchar* ptr = img.ptr<uchar>(i);
for(int j=0; j<img.cols; j++)
imgBuffer[i*img.cols+j] = ptr[j];
}
_size srcSize;
srcSize.height = img.rows;
srcSize.width = img.cols;
_point corners[MAXPOINTS];
int numCorners = 0;
_cornerHarris(imgBuffer, srcSize, corners, &numCorners, 0.04);
cv::Mat colorImg;
cv::cvtColor(img, colorImg, CV_GRAY2BGR);
for(int i=0; i<numCorners; i++)
{
cv::Point center(corners[i].x, corners[i].y);
cv::circle(colorImg, center, 2, cv::Scalar(0,0,255));
}
cv::namedWindow("show");
cv::imshow("show", colorImg);
cv::waitKey(0);
}
template<typename T1, typename T2>
void convolution(T1* img, T2** out, _size imgSize, T2* kernel, _size kerSize)
{
for(int i=0; i<imgSize.height; i++)
{
for(int j=0; j<imgSize.width; j++)
{
int count = 0;
T2 sumValue = 0;
for(int m=i-kerSize.height; m<=i+kerSize.height; m++)
{
for(int n=j-kerSize.width; n<=j+kerSize.width; n++)
{
if(m>=0 && m<imgSize.height && n>=0 && n<imgSize.width)
sumValue += T2(img[m*imgSize.width+n])*kernel[count];
count++;
}
}
(*out)[i*imgSize.width+j] = sumValue;
}
}
}
void dotMultiply(int* src1, int* src2, int** dst, _size srcSize)
{
for(int i=0; i<srcSize.height; i++)
{
for(int j=0; j<srcSize.width; j++)
(*dst)[i*srcSize.width+j] = src1[i*srcSize.width+j]*src2[i*srcSize.width+j];
}
}
void gaussian(int* src, float** out, _size imgSize, _size gaussSize, float sita)
{
int sizeKernel = (2*gaussSize.height+1) * (2*gaussSize.width+1);
float* gaussKernel = new float[sizeKernel];
for(int i=-gaussSize.height; i<=gaussSize.height; i++)
{
for(int j=-gaussSize.width; j<=gaussSize.width; j++)
{
float tmp = -1*(i*i+j*j)/(2*sita*sita);
gaussKernel[(i+gaussSize.height)*(2*gaussSize.width+1)+(j+gaussSize.width)] =
exp(tmp)/(2*PI*sita*sita);
}
}
convolution(src, out, imgSize, gaussKernel, gaussSize);
}
void strengthR(float* A, float* B, float* C, float** R, float* maxR, _size imgSize, float alpha)
{
*maxR = 0;
for(int i=0; i<imgSize.height; i++)
{
for(int j=0; j<imgSize.width; j++)
{
float detM = A[i*imgSize.width+j]*B[i*imgSize.width+j] - C[i*imgSize.width+j]*C[i*imgSize.width+j];
float traceM = A[i*imgSize.width+j] + B[i*imgSize.width+j];
float result = detM - alpha*(traceM*traceM);
if(result > *maxR)
*maxR = result;
(*R)[i*imgSize.width+j] = result;
}
}
}
void depressR(float* srcR, _point* dstR, int* numDst, _size imgSize, _size depressWin, float threshold)
{
*numDst = 0;
for(int i=0; i<imgSize.height; i++)
{
for(int j=0; j<imgSize.width; j++)
{
float flagValue = srcR[i*imgSize.width+j]<threshold?0:srcR[i*imgSize.width+j];
int numPoint = 0, numPass = 0;
for(int m=i-depressWin.height; m<=i+depressWin.height; m++)
{
for(int n=j-depressWin.width; n<=j+depressWin.width; n++)
{
if(m>=0 && m<imgSize.height && n>=0 && n<imgSize.width)
{
float compareValue = srcR[m*imgSize.width+n]<threshold?0:srcR[m*imgSize.width+n];
if(flagValue > compareValue)
numPass ++;
numPoint ++;
}
}
}
if(numPoint == numPass+1)
{
_point corner;
corner.x = j;
corner.y = i;
dstR[(*numDst)++] = corner;
}
}
}
}
void _cornerHarris(unsigned char* src, _size srcSize, _point* corners, int* numCorner, float alpha)
{
int kernel[3] = {-1, 0, 1};
_size sizeX, sizeY;
sizeX.height = 0; sizeX.width = 1;
sizeY.height = 1; sizeY.width = 0;
int* gradientX = new int[srcSize.height*srcSize.width];
int* gradientY = new int[srcSize.height*srcSize.width];
convolution(src, &gradientX, srcSize, kernel, sizeX);
convolution(src, &gradientY, srcSize, kernel, sizeY);
int* gradientXX = new int[srcSize.height*srcSize.width];
int* gradientYY = new int[srcSize.height*srcSize.width];
int* gradientXY = new int[srcSize.height*srcSize.width];
dotMultiply(gradientX, gradientX, &gradientXX, srcSize);
dotMultiply(gradientY, gradientY, &gradientYY, srcSize);
dotMultiply(gradientX, gradientY, &gradientXY, srcSize);
delete[] gradientX;
delete[] gradientY;
_size gaussSize;
gaussSize.height = 1; gaussSize.width = 1;
float sita = 1;
float* M_A = new float[srcSize.height*srcSize.width];
float* M_B = new float[srcSize.height*srcSize.width];
float* M_C = new float[srcSize.height*srcSize.width];
gaussian(gradientXX, &M_A, srcSize, gaussSize, sita);
gaussian(gradientYY, &M_B, srcSize, gaussSize, sita);
gaussian(gradientXY, &M_C, srcSize, gaussSize, sita);
delete[] gradientXX;
delete[] gradientYY;
delete[] gradientXY;
float maxR = 0;
float* R = new float[srcSize.height*srcSize.width];
strengthR(M_A, M_B, M_C, &R, &maxR, srcSize, alpha);
delete[] M_A;
delete[] M_B;
delete[] M_C;
float threshold = 0.1;
threshold *= maxR;
_size depressSize;
depressSize.height = 1; depressSize.width = 1;
depressR(R, corners, numCorner, srcSize, depressSize, threshold);
delete[] R;
}
效果图(左:原图,中:旋转30度,右:放大2倍)