图像形状特征（七）--Zernike矩



原文：http://blog.csdn.NET/wrj19860202/article/details/6334275

Zernike在1934年引入了一组定义在单位圆 上的复值函数集{ }，{ }具有完备性和正交性，使得它可以表示定义在单位圆盘内的任何平方可积函数。其定义为：

表示原点到点 的矢量长度； 表示矢量 与 轴逆时针方向的夹角。

是实值径向多项式：

称为Zernike多项式。

Zernike多项式满足正交性：

其中

为克罗内克符号，

是 的共轭多项式。

由于Zernike多项式的正交完备性，所以在单位圆内的任何图像 都可以唯一的用下面式子来展开：

式子中 就是Zernike矩，其定义为：

注意式子中 和 采用的是不同的坐标系（ 采用直角坐标，而 采用的极坐标系，在计算的时候要进行坐标转换）

对于离散的数字图像，可将积分形式改为累加形式：

我们在计算一副图像的Zernike矩时，必须将图像的中心移到坐标的原点，将图像的像素点映射到单位圆内，由于Zernike矩具有旋转不变性，我们可以将 作为图像的不变特征，其中图像的低频特征有p值小的 提取，高频特征由p值高的 提取。从上面可以看出，Zernike矩可以构造任意高阶矩。

由于Zernike矩只具有旋转不变性，不具有平移和尺度不变性，所以要提前对图像进行归一化，我们采用标准矩的方法来归一化一副图像，标准矩定义为：

，

由标准矩我们可以得到图像的"重心"，

我们将图像的"重心"移动到单位圆的圆心（即坐标的原点），便解决了平移问题。

我们知道 表征了图像的"面积"，归一图像的尺度无非就是把他们的大小变为一致的，（这里的大小指的是图像目标物的大小，不是整幅图像的大小，"面积"也是目标物的"面积"）。

所以，对图像进行变换 就可以达到图像尺寸一致的目的。

综合上面结果，对图像进行 变换，最终图像 的Zernike矩就是平移，尺寸和旋转不变的。

Zernike 不变矩相比 Hu 不变矩识别效果会好一些，因为他描述了图像更多的细节内容，特别是高阶矩，但是由于 Zernike 不变矩计算时间比较长，所以出现了很多快速的算法，大家可以 google 一下。

用 Zernike 不变矩来识别手势轮廓，识别率大约在 40%~50% 之间，跟 Hu 不变矩一样， Zernike 不变矩一般用来描述目标物形状占优势的图像，不适合用来描述纹理丰富的图像，对于纹理图像，识别率一般在 20%~30% 左右，很不占优势。

[cpp] view plain copy print ?

1. C++代码如下：  
2.   
3. /*计算一行的像素个数 
4. imwidth:图像宽度 
5. deep:图像深度（8位灰度图为1，24位彩色图为3） 
6. */  
7. #define  bpl(imwidth, deep) ((imwidth*deep*8+31)/32*4)  
8. /*获取像素值 
9. psrcBmp:图像数据指针 
10. nsrcBmpWidth:图像宽度，以像素为单位 
11. x,y:像素点 
12. deep:图像的位数深度，（1表示8位的灰度图，3表示24位的RGB位图） 
13. */  
14. COLORREF J_getpixel( const BYTE *psrcBmp, const int nsrcBmpWidth, const int x, const int y, int deep = 3)  
15. {  
16.     if (deep == 3)  
17.     {  
18.         return RGB(*(psrcBmp + x*3 + y*bpl(nsrcBmpWidth, deep) + 2 ) ,   
19.             *(psrcBmp + x*3 + y*bpl(nsrcBmpWidth, deep) + 1 ) ,   
20.             *(psrcBmp + x*3 + y*bpl(nsrcBmpWidth, deep) +0 ));  
21.     }  
22.     else if (deep == 1)  
23.     {  
24.         return *(psrcBmp + x + y*bpl(nsrcBmpWidth, deep));  
25.     }  
26. }  
27.   
28. //获取标准矩(只支持8位灰度图)  
29. void GetStdMoment(BYTE *psrcBmp ,   
30.                  int nsrcBmpWidth,  
31.                  int nsrcBmpHeight,  
32.                  double *m)  
33. {  
34.     for ( int p = 0 ; p < 2 ; p++ )  
35.         for ( int q = 0 ; q < 2 ; q++ )  
36.         {  
37.             if( p == 1 && q == 1)  
38.                 break;  
39.             for ( int y = 0 ; y < nsrcBmpHeight ; y++ )  
40.                 for ( int x = 0 ; x < nsrcBmpWidth ; x++ )  
41.                     m[p*2+q] += (pow( (double)x , p ) * pow( (double)y , q ) * J_getpixel(psrcBmp , nsrcBmpWidth , x ,y, 1));  
42.         }  
43. }  
44.   
45. //阶乘  
46. double Factorial( int n )  
47. {  
48.     if( n < 0 )  
49.         return -1;  
50.   
51.     double m = 1;  
52.     for(int i = 2 ; i <= n ; i++)  
53.     {  
54.         m *= i;  
55.     }  
56.     return m;  
57. }  
58.   
59. //阶乘数，计算好方便用，提高速度  
60. double factorials[11] = {1 , 1 , 2 , 6 , 24 , 120 , 720 , 5040 , 40320 , 362880 , 39916800};  
61.   
62. //把图像映射到单位圆，获取像素极坐标半径  
63. double GetRadii(int nsrcBmpWidth,  
64.        int nsrcBmpHeight,  
65.        int x0,  
66.        int y0,  
67.        int x,  
68.        int y)  
69. {  
70.     double lefttop = sqrt(((double)0 - x0)*(0 - x0) + (0 - y0)*(0 - y0));  
71.     double righttop = sqrt(((double)nsrcBmpWidth - 1 - x0)*(nsrcBmpWidth - 1 - x0) + (0 - y0)*(0 - y0));  
72.     double leftbottom = sqrt(((double)0 - x0)*(0 - x0) + (nsrcBmpHeight - 1 - y0)*(nsrcBmpHeight - 1 - y0));  
73.     double rightbottom = sqrt(((double)nsrcBmpWidth - 1 - x0)*(nsrcBmpWidth - 1 - x0) + (nsrcBmpHeight - 1 - y0)*(nsrcBmpHeight - 1 - y0));  
74.   
75.     double maxRadii = lefttop;  
76.     maxRadii < righttop ? righttop : maxRadii;  
77.     maxRadii < leftbottom ? leftbottom : maxRadii;  
78.     maxRadii < rightbottom ? rightbottom : maxRadii;  
79.   
80.     double Radii = sqrt(((double)x - x0)*(x - x0) + (y - y0)*(y - y0))/maxRadii;  
81.     if(Radii > 1)  
82.     {  
83.         Radii = 1;  
84.     }  
85.     return Radii;  
86. }  
87.   
88. //把图像映射到单位圆，获取像素极坐标角度  
89. double GetAngle(int nsrcBmpWidth,  
90.                 int nsrcBmpHeight,  
91.                 int x,  
92.                 int y)  
93. {  
94.     double o;  
95.   
96.     double dia = sqrt((double)nsrcBmpWidth*nsrcBmpWidth + nsrcBmpHeight*nsrcBmpHeight);  
97.     int x0 = nsrcBmpWidth / 2;  
98.     int y0 = nsrcBmpHeight / 2;  
99.     double x_unity = (x - x0)/(dia/2);   
100.     double y_unity = (y - y0)/(dia/2);  
101.   
102.     if( x_unity == 0 && y_unity >= 0 )  
103.         o=pi/2;  
104.     else if( x_unity ==0 && y_unity <0)  
105.         o=1.5*pi;  
106.     else  
107.         o=atan( y_unity / x_unity );  
108.     if(o*y<0)    //第三象限  
109.         o=o+pi;  
110.   
111.     return o;  
112. }  
113.   
114. //Zernike不变矩  
115. J_GetZernikeMoment(BYTE *psrcBmp ,   
116.                 int nsrcBmpWidth,  
117.                 int nsrcBmpHeight,  
118.                 double *Ze )  
119. {  
120.     double R[count][count] = {0.0};  
121.     double V[count][count] = {0.0};  
122.   
123.     double M[4] = {0.0};  
124.     GetStdMoment(psrcBmp , nsrcBmpWidth , nsrcBmpHeight , M);  
125.     int x0 = (int)(M[2]/M[0]+0.5);  
126.     int y0 = (int)(M[1]/M[0]+0.5);  
127.   
128.     for(int n = 0 ; n < count ; n++)  
129.     {  
130.         for (int m = 0 ; m < count ; m++)  
131.         {  
132.             //优化算法，只计算以下介数  
133.   
134.             if( (n == 1 && m == 0) ||  
135.                 (n == 1 && m == 1) ||  
136.                 (n == 2 && m == 0) ||  
137.                 (n == 2 && m == 1) ||  
138.                 (n == 2 && m == 2) ||  
139.                 (n == 3 && m == 0) ||  
140.                 (n == 3 && m == 1) ||  
141.                 (n == 3 && m == 2) ||  
142.                 (n == 3 && m == 3) ||  
143.                 (n == 4 && m == 0) ||  
144.                 (n == 4 && m == 1) ||  
145.                 (n == 4 && m == 2) ||  
146.                 (n == 4 && m == 3) ||  
147.                 (n == 4 && m == 4))  
148.   
149.             {  
150.                 for(int y = 0 ; y < nsrcBmpHeight ; y++)  
151.                 {  
152.                     for (int x = 0 ; x < nsrcBmpWidth ; x++)  
153.                     {  
154.                         for(int s = 0 ; (s <= (n - m)/2 ) && n >= m ; s++)  
155.                         {  
156.                             R[n][m] += pow( -1.0, s )  
157.                                     * ( n - s > 10 ? Factorial( n - s ) : factorials[ n - s ] )  
158.                                     * pow( GetRadii( nsrcBmpWidth, nsrcBmpHeight, x0, y0, x, y ), n - 2 * s )  
159.                                     / ( ( s > 10 ? Factorial( s ) : factorials[ s ] )  
160.                                     * ( ( n + m ) / 2 - s > 10 ? Factorial( ( n + m ) / 2 - s ) : factorials[ ( n + m ) / 2 - s ] )  
161.                                     * ( ( n - m ) / 2 - s > 10 ? Factorial( ( n - m ) / 2 - s ) : factorials[ ( n - m ) / 2 - s ] ) );  
162.                         }  
163.                         Ze[ n * count + m ] += R[ n ][ m ]  
164.                                             * J_getpixel( psrcBmp, nsrcBmpWidth, x ,y, 1)  
165.                                             * cos( m * GetAngle( nsrcBmpWidth, nsrcBmpHeight, x, y) );//实部  
166.   
167.                         V[n][m] += R[ n ][ m ]   
168.                                 * J_getpixel( psrcBmp, nsrcBmpWidth, x, y, 1)  
169.                                 * sin( m * GetAngle( nsrcBmpWidth, nsrcBmpHeight, x, y ) );//虚部  
170.   
171.                         R[n][m] = 0.0;  
172.                     }  
173.                 }  
174.                 *(Ze+n*count + m) = sqrt( (*(Ze+n*count + m))*(*(Ze+n*count + m)) + V[n][m]*V[n][m] )*(n+1)/pi/M[0];  
175.             }  
176.         }  
177.     }  
178. }