图像的局部特征研究--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矩作为特征描述子的用途比较大，但是其他地方不考虑直接用这一算法。