图像的局部特征研究--Hu不变矩

研究问题：最近需要对九宫格里面的九幅图像做出识别，判断哪一个格子里面的图片才是最特殊的，而其他图片是相似的，即其他八个格子里的图可以看作是一图作旋转，增添，缩放，镜像等操作，就像下面这样：

                    

图1.右下角的小人是要找的目标图像，其他的小车可看作有其中一幅旋转平移等操作得到

我的思路是先对这个图片进行分割，分割成9幅小图，再对每一幅图片求特征表示及描述，根据这些图的特点，选择了一个比较简单的衡量图像特征的区域描述子--Hu矩。

下面详细叙述Hu矩的不变性：

在黎曼积分意义下，将二维（p+q）阶矩定义为：

                                                                 （1）

可以看作是一个密度分布函数，在这里我们研究的对象是一幅图，其分布函数区间为 [ 0, 255]。

相应的（p+q）阶中心矩则定义如下：

                                                              （2）

其中

                                                                                                            （3）

将积分符号换成求和符号会比较助于理解，但是为了保持证明的完整性，仍然保留积分符号考虑整个实数域。

①中心距对于的平移具有不变性：

假如新的坐标为：                      ，其中α，β是常数                                                                                  （4）

易知和的中心距是相同的。

②矩对于缩放具有不变性：

         ，α是个常数                                                                                                                                       （5）

可以看作是(x, y)分别乘于系数α得到，对于每一个α系数由公式（2），有：

                                                                                                                                                                         （6）

因为α是个常数，那么变换前后的中心距有这样的关系：

                                                                                                                                                                         （7）

最后可以得到：

                                                                                            （8）

这也可以称为相似矩不变性。

③矩对于旋转具有不变性：

                                                                                                                                              （9）

旋转矩阵的模为1。

                                                                                                                                       （10）

将（9）和（10）跟公式（2）结合，也可以得到和的关系，这可以称为正交不变性。

上述这些只是数学上的简单说明，论证的是图像矩的不变性，对于实际的应用，更常用的是用二阶矩和三阶矩构造出的7个不变矩组：

首先，对中心矩进行归一化：

， 其中，p+q = 2,3......。                                                                                                                            （11）

那么，存在7个不变矩组为：

                                                                                   （12）

这些矩组对于平移、尺度变化、镜像（内部为负号）和旋转是不变的。

为了解决一开始提出的问题，就简单地写了一个程序：

<pre name="code" class="cpp">void calcHu(Mat image)
{
	int bmpWidth = image.cols;
	int bmpHeight = image.rows;
	int bmpStep = image.step; 
	int bmpChannels = image.channels();
	uchar* pBmpBuf = image.data;

	double m00=0,m11=0,m20=0,m02=0,m30=0,m03=0,m12=0,m21=0;  //中心矩 
	double x0=0,y0=0;                                        //计算中心距时所使用的临时变量（x-x'） 
	double u20=0,u02=0,u11=0,u30=0,u03=0,u12=0,u21=0;        //规范化后的中心矩
	double t1=0,t2=0,t3=0,t4=0,t5=0;//临时变量
	int Center_x=0,Center_y=0;//重心 
	int i,j;            //循环变量

	//  获得图像的区域重心(普通矩)
	double s10=0,s01=0,s00=0;  //0阶矩和1阶矩  
	for(j=0;j<bmpHeight;j++)//y
	{
		for(i=0;i<bmpWidth;i++)//x
		{
			s10+=i*pBmpBuf[j*bmpStep+i];
			s01+=j*pBmpBuf[j*bmpStep+i];
			s00+=pBmpBuf[j*bmpStep+i];
		}
	}
	Center_x=(int)(s10/s00+0.5);
	Center_y=(int)(s01/s00+0.5);

	//  计算二阶、三阶矩(中心矩)
	m00=s00; 
	for(j=0;j<bmpHeight;j++) 
	{
		for(i=0;i<bmpWidth;i++)//x 
		{ 
			x0=(i-Center_x); 
			y0=(j-Center_y); 
			m11+=x0*y0*pBmpBuf[j*bmpStep+i]; 
			m20+=x0*x0*pBmpBuf[j*bmpStep+i]; 
			m02+=y0*y0*pBmpBuf[j*bmpStep+i]; 
			m03+=y0*y0*y0*pBmpBuf[j*bmpStep+i];
			m30+=x0*x0*x0*pBmpBuf[j*bmpStep+i]; 
			m12+=x0*y0*y0*pBmpBuf[j*bmpStep+i]; 
			m21+=x0*x0*y0*pBmpBuf[j*bmpStep+i]; 
		} 
	} 

	//  计算规范化后的中心矩: mij/pow(m00,((i+j+2)/2)
	u20=m20/pow(m00,2); 
	u02=m02/pow(m00,2); 
	u11=m11/pow(m00,2);
	u30=m30/pow(m00,2.5); 
	u03=m03/pow(m00,2.5);
	u12=m12/pow(m00,2.5); 
	u21=m21/pow(m00,2.5);

	//  计算中间变量
	t1=(u20-u02); 
	t2=(u30-3*u12); 
	t3=(3*u21-u03); 
	t4=(u30+u12);
	t5=(u21+u03);

	//  计算不变矩 
	M[0]=u20+u02; 
	M[1]=t1*t1+4*u11*u11; 
	M[2]=t2*t2+t3*t3; 
	M[3]=t4*t4+t5*t5;
	M[4]=t2*t4*(t4*t4-3*t5*t5)+t3*t5*(3*t4*t4-t5*t5); 
	M[5]=t1*(t4*t4-t5*t5)+4*u11*t4*t5;
	M[6]=t3*t4*(t4*t4-3*t5*t5)-t2*t5*(3*t4*t4-t5*t5);
}

完整的程序如下：

#include <opencv2/opencv.hpp>
#include <ctime>
using namespace std;
using namespace cv;

double M[7] = {0}; 

vector<Mat> cutImage(Mat image)
{
	vector<Mat> roi;
	int rwidth = image.cols/3;
	int rheight = image.rows/3;

	for (int y = 0;y < image.rows-10;y +=rheight)
	{
		for (int x = 0;x <  image.cols-10;x += rwidth)
		{
			Rect rect(x, y, rwidth, rheight);
			rect &= Rect(0, 0, image.cols, image.rows);;
			roi.push_back(image(rect));
		}
	}
	return roi;
}

void calcHu(Mat image)
{
	int bmpWidth = image.cols;
	int bmpHeight = image.rows;
	int bmpStep = image.step; 
	int bmpChannels = image.channels();
	uchar* pBmpBuf = image.data;

	double m00=0,m11=0,m20=0,m02=0,m30=0,m03=0,m12=0,m21=0;  //中心矩 
	double x0=0,y0=0;                                        //计算中心距时所使用的临时变量（x-x'） 
	double u20=0,u02=0,u11=0,u30=0,u03=0,u12=0,u21=0;        //规范化后的中心矩
	double t1=0,t2=0,t3=0,t4=0,t5=0;//临时变量
	int Center_x=0,Center_y=0;//重心 
	int i,j;            //循环变量

	//  获得图像的区域重心(普通矩)
	double s10=0,s01=0,s00=0;  //0阶矩和1阶矩  
	for(j=0;j<bmpHeight;j++)//y
	{
		for(i=0;i<bmpWidth;i++)//x
		{
			s10+=i*pBmpBuf[j*bmpStep+i];
			s01+=j*pBmpBuf[j*bmpStep+i];
			s00+=pBmpBuf[j*bmpStep+i];
		}
	}
	Center_x=(int)(s10/s00+0.5);
	Center_y=(int)(s01/s00+0.5);

	//  计算二阶、三阶矩(中心矩)
	m00=s00; 
	for(j=0;j<bmpHeight;j++) 
	{
		for(i=0;i<bmpWidth;i++)//x 
		{ 
			x0=(i-Center_x); 
			y0=(j-Center_y); 
			m11+=x0*y0*pBmpBuf[j*bmpStep+i]; 
			m20+=x0*x0*pBmpBuf[j*bmpStep+i]; 
			m02+=y0*y0*pBmpBuf[j*bmpStep+i]; 
			m03+=y0*y0*y0*pBmpBuf[j*bmpStep+i];
			m30+=x0*x0*x0*pBmpBuf[j*bmpStep+i]; 
			m12+=x0*y0*y0*pBmpBuf[j*bmpStep+i]; 
			m21+=x0*x0*y0*pBmpBuf[j*bmpStep+i]; 
		} 
	} 

	//  计算规范化后的中心矩: mij/pow(m00,((i+j+2)/2)
	u20=m20/pow(m00,2); 
	u02=m02/pow(m00,2); 
	u11=m11/pow(m00,2);
	u30=m30/pow(m00,2.5); 
	u03=m03/pow(m00,2.5);
	u12=m12/pow(m00,2.5); 
	u21=m21/pow(m00,2.5);

	//  计算中间变量
	t1=(u20-u02); 
	t2=(u30-3*u12); 
	t3=(3*u21-u03); 
	t4=(u30+u12);
	t5=(u21+u03);

	//  计算不变矩 
	M[0]=u20+u02; 
	M[1]=t1*t1+4*u11*u11; 
	M[2]=t2*t2+t3*t3; 
	M[3]=t4*t4+t5*t5;
	M[4]=t2*t4*(t4*t4-3*t5*t5)+t3*t5*(3*t4*t4-t5*t5); 
	M[5]=t1*(t4*t4-t5*t5)+4*u11*t4*t5;
	M[6]=t3*t4*(t4*t4-3*t5*t5)-t2*t5*(3*t4*t4-t5*t5);
}

int compareHu(double mo[9][7])
{
	int no = 0;
	double sum = 0, min = 100, max = 0;
	
	for (int i = 0;i<9;i++)
	{
		sum += mo[i][0];
		if(mo[i][0]>max)
			max = mo[i][0];
		if(mo[i][0]<min)
			min = mo[i][0];
	}

	sum /= 9;
	if(sum - min > max - sum) max = min;

	for (int i = 0;i<9;i++)
		if (mo[i][0]==max){ no = i;break;}
	
	//cout<<max<<endl<<no<<endl;
	return no;
}

void drawCross(int n, Mat image)
{
	int centerx = 0, centery = 0, widstep = 0, heistep = 0;

	widstep = image.cols/6;
	heistep = image.rows/6;
	
	centerx = n%3 * widstep * 2 + widstep;
	centery = n/3 * heistep * 2 + heistep;
	
	Scalar color(0, 0, 255);
	line(image, Point(centerx-20, centery), Point(centerx+20, centery), color, 2);
	line(image, Point(centerx, centery-10), Point(centerx, centery+10), color, 2);
}

int main()
{   
	Mat src, binary;
	vector <Mat> srcRoi;
	int count = 12;
	double moment[9][7] = {0};

	while (--count)
	{
		char imageName[10];
		sprintf_s(imageName, "%d.jpg", count);
		src = imread(imageName, 1);
		Canny(src, binary, 50, 100);

		srcRoi = cutImage(binary);

		for (int i = 0;i < 9;i++)
		{
			calcHu(srcRoi[i]);

			for (int j = 0;j<7;j++)
				moment[i][j] = M[j];	
		}

		int no = compareHu(moment);
		drawCross(no, src);
		imshow(imageName, src);
		waitKey(0);
	}

	return 0;
}

总结：Hu不变矩是一个比较简单而且有效的特征描述子，需要注意的是，要计算的图片最好是处理过的边缘图或者轮廓图，这能提高鲁棒性，但是，相比其他匹配算法如SURF来说，Hu矩的局限也很大，它受图片的质量影响比较大，可以说，Hu矩作为特征描述子的用途比较大，但是其他地方不考虑直接用这一算法。