不变矩
图像算法:不变矩
一 工具:VC6.0+OpenCV1.0
语言:CPP
二 原理:http://blog.csdn.net/byxdaz/archive/2009/12/28/5089448.aspx
三 说明:
本源码参考了http://blog.csdn.net/dadaadao/archive/2011/01/04/6114705.aspx,对里面的一些错误(比较多)进行了修正,同时增加了匹配功能,即可以用来实现不变矩匹配,但在匹配中总发现有点点蹊跷,反复查看程序未发现问题,暂时放在博客,如果大家发现了错误,请不吝赐教!
四 源码:
1 /*===============================================//
2 功能:不变矩匹配
3 时间:3/28/2011 SkySeraph HQU
4 参考:
5 //===============================================*/
6 #include "iostream"
7 usingnamespace std;
8
9 #include "cv.h"
10 #include "highgui.h"
11
12 #include "math.h"
13
14 #pragma comment(lib,"highgui.lib")
15 #pragma comment(lib,"cv.lib")
16 #pragma comment(lib,"cvaux.lib")
17 #pragma comment(lib,"cxcore.lib")
18
19 constchar* filename ="D:\\My Documents\\My Pictures\\Images\\1.bmp";
20 constchar* filename2 ="D:\\My Documents\\My Pictures\\Images\\2.bmp";
21
22 /*=============================================*/
23 double M[7] = {0}; //HU不变矩
24
25 bool HuMoment(IplImage* img)
26 {
27
28 int bmpWidth = img->width;
29 int bmpHeight = img->height;
30 int bmpStep = img->widthStep;
31 int bmpChannels = img->nChannels;
32 uchar*pBmpBuf = (uchar*)img->imageData;
33
34 double m00=0,m11=0,m20=0,m02=0,m30=0,m03=0,m12=0,m21=0; //中心矩
35 double x0=0,y0=0; //计算中心距时所使用的临时变量(x-x')
36 double u20=0,u02=0,u11=0,u30=0,u03=0,u12=0,u21=0;//规范化后的中心矩
37 //double M[7]; //HU不变矩
38 double t1=0,t2=0,t3=0,t4=0,t5=0;//临时变量,
39 //double Center_x=0,Center_y=0;//重心
40 int Center_x=0,Center_y=0;//重心
41 int i,j; //循环变量
42
43 // 获得图像的区域重心
44 double s10=0,s01=0,s00=0; //0阶矩和1阶矩 //注:二值图像的0阶矩表示面积
45 for(j=0;j<bmpHeight;j++)//y
46 {
47 for(i=0;i<bmpWidth;i++)//x
48 {
49 s10+=i*pBmpBuf[j*bmpStep+i];
50 s01+=j*pBmpBuf[j*bmpStep+i];
51 s00+=pBmpBuf[j*bmpStep+i];
52 }
53 }
54 Center_x=(int)(s10/s00+0.5);
55 Center_y=(int)(s01/s00+0.5);
56
57 // 计算二阶、三阶矩
58 m00=s00;
59 for(j=0;j<bmpHeight;j++)
60 {
61 for(i=0;i<bmpWidth;i++)//x
62 {
63 x0=(i-Center_x);
64 y0=(j-Center_y);
65 m11+=x0*y0*pBmpBuf[j*bmpStep+i];
66 m20+=x0*x0*pBmpBuf[j*bmpStep+i];
67 m02+=y0*y0*pBmpBuf[j*bmpStep+i];
68 m03+=y0*y0*y0*pBmpBuf[j*bmpStep+i];
69 m30+=x0*x0*x0*pBmpBuf[j*bmpStep+i];
70 m12+=x0*y0*y0*pBmpBuf[j*bmpStep+i];
71 m21+=x0*x0*y0*pBmpBuf[j*bmpStep+i];
72 }
73 }
74
75 // 计算规范化后的中心矩
76 u20=m20/pow(m00,2);
77 u02=m02/pow(m00,2);
78 u11=m11/pow(m00,2);
79 u30=m30/pow(m00,2.5);
80 u03=m03/pow(m00,2.5);
81 u12=m12/pow(m00,2.5);
82 u21=m21/pow(m00,2.5);
83
84 // 计算中间变量。
85 t1=(u20-u02);
86 t2=(u30-3*u12);
87 t3=(3*u21-u03);
88 t4=(u30+u12);
89 t5=(u21+u03);
90
91 // 计算不变矩
92 M[0]=u20+u02;
93 M[1]=t1*t1+4*u11*u11;
94 M[2]=t2*t2+t3*t3;
95 M[3]=t4*t4+t5*t5;
96 M[4]=t2*t4*(t4*t4-3*t5*t5)+t3*t5*(3*t4*t4-t5*t5);
97 M[5]=t1*(t4*t4-t5*t5)+4*u11*t4*t5;
98 M[6]=t3*t4*(t4*t4-3*t5*t5)-t2*t5*(3*t4*t4-t5*t5);
99
100
101 /*cout<<M[0]<<endl;//<<二"<<M[0]<<"三"<<M[0]<<"四"<<M[0]<<"五"<<M[0]<<"六"<<M[0]<<"七"<<M[0]<<endl;
102 cout<<M[1]<<endl;
103 cout<<M[2]<<endl;
104 cout<<M[3]<<endl;
105 cout<<M[4]<<endl;
106 cout<<M[5]<<endl;
107 cout<<M[6]<<endl;
108 cout<<endl;*/
109 returntrue;
110 }
111
112
113 int main(char argc,char** argv)
114 {
115 int i;
116 double Sa[7] = {0},Ta[7] ={0};
117
118 ///*源图像
119 IplImage*img = cvLoadImage(filename,0);//灰度
120 HuMoment(img);
121 for(i=0;i<7;i++)
122 {
123 Sa[i] = M[i];
124 M[i] =0;
125 }
126 cout<<Sa[0]<<endl;
127 cout<<Sa[1]<<endl;
128 cout<<Sa[2]<<endl;
129 cout<<Sa[3]<<endl;
130 cout<<Sa[4]<<endl;
131 cout<<Sa[5]<<endl;
132 cout<<Sa[6]<<endl;
133 cout<<endl;
134 //*/
135
136
137 ///*模板图
138 IplImage*tpl = cvLoadImage(filename2,0);//灰度
139 HuMoment(tpl);
140 for(i=0;i<7;i++)
141 {
142 Ta[i] = M[i];
143 M[i] =0;
144 }
145 cout<<Ta[0]<<endl;
146 cout<<Ta[1]<<endl;
147 cout<<Ta[2]<<endl;
148 cout<<Ta[3]<<endl;
149 cout<<Ta[4]<<endl;
150 cout<<Ta[5]<<endl;
151 cout<<Ta[6]<<endl;
152 cout<<endl;
153
154
155 // 计算相似度
156 double dbR =0; //相似度
157 double dSigmaST =0;
158 double dSigmaS =0;
159 double dSigmaT =0;
160 double temp =0;
161
162 for(i=0;i<7;i++)
163 {
164 temp = Sa[i]*Ta[i];
165 dSigmaST+=temp;
166 dSigmaS+=pow(Sa[i],2);
167 dSigmaT+=pow(Ta[i],2);
168 }
169 dbR = dSigmaST/(sqrt(dSigmaS)*sqrt(dSigmaT));
170 printf("%lf\n",dbR);
171 //cout<<dbR<<endl;
172
173 cvReleaseImage(&img);
174 cvReleaseImage(&tpl);
175
176 return0;
177 }
2 功能:不变矩匹配
3 时间:3/28/2011 SkySeraph HQU
4 参考:
5 //===============================================*/
6 #include "iostream"
7 usingnamespace std;
8
9 #include "cv.h"
10 #include "highgui.h"
11
12 #include "math.h"
13
14 #pragma comment(lib,"highgui.lib")
15 #pragma comment(lib,"cv.lib")
16 #pragma comment(lib,"cvaux.lib")
17 #pragma comment(lib,"cxcore.lib")
18
19 constchar* filename ="D:\\My Documents\\My Pictures\\Images\\1.bmp";
20 constchar* filename2 ="D:\\My Documents\\My Pictures\\Images\\2.bmp";
21
22 /*=============================================*/
23 double M[7] = {0}; //HU不变矩
24
25 bool HuMoment(IplImage* img)
26 {
27
28 int bmpWidth = img->width;
29 int bmpHeight = img->height;
30 int bmpStep = img->widthStep;
31 int bmpChannels = img->nChannels;
32 uchar*pBmpBuf = (uchar*)img->imageData;
33
34 double m00=0,m11=0,m20=0,m02=0,m30=0,m03=0,m12=0,m21=0; //中心矩
35 double x0=0,y0=0; //计算中心距时所使用的临时变量(x-x')
36 double u20=0,u02=0,u11=0,u30=0,u03=0,u12=0,u21=0;//规范化后的中心矩
37 //double M[7]; //HU不变矩
38 double t1=0,t2=0,t3=0,t4=0,t5=0;//临时变量,
39 //double Center_x=0,Center_y=0;//重心
40 int Center_x=0,Center_y=0;//重心
41 int i,j; //循环变量
42
43 // 获得图像的区域重心
44 double s10=0,s01=0,s00=0; //0阶矩和1阶矩 //注:二值图像的0阶矩表示面积
45 for(j=0;j<bmpHeight;j++)//y
46 {
47 for(i=0;i<bmpWidth;i++)//x
48 {
49 s10+=i*pBmpBuf[j*bmpStep+i];
50 s01+=j*pBmpBuf[j*bmpStep+i];
51 s00+=pBmpBuf[j*bmpStep+i];
52 }
53 }
54 Center_x=(int)(s10/s00+0.5);
55 Center_y=(int)(s01/s00+0.5);
56
57 // 计算二阶、三阶矩
58 m00=s00;
59 for(j=0;j<bmpHeight;j++)
60 {
61 for(i=0;i<bmpWidth;i++)//x
62 {
63 x0=(i-Center_x);
64 y0=(j-Center_y);
65 m11+=x0*y0*pBmpBuf[j*bmpStep+i];
66 m20+=x0*x0*pBmpBuf[j*bmpStep+i];
67 m02+=y0*y0*pBmpBuf[j*bmpStep+i];
68 m03+=y0*y0*y0*pBmpBuf[j*bmpStep+i];
69 m30+=x0*x0*x0*pBmpBuf[j*bmpStep+i];
70 m12+=x0*y0*y0*pBmpBuf[j*bmpStep+i];
71 m21+=x0*x0*y0*pBmpBuf[j*bmpStep+i];
72 }
73 }
74
75 // 计算规范化后的中心矩
76 u20=m20/pow(m00,2);
77 u02=m02/pow(m00,2);
78 u11=m11/pow(m00,2);
79 u30=m30/pow(m00,2.5);
80 u03=m03/pow(m00,2.5);
81 u12=m12/pow(m00,2.5);
82 u21=m21/pow(m00,2.5);
83
84 // 计算中间变量。
85 t1=(u20-u02);
86 t2=(u30-3*u12);
87 t3=(3*u21-u03);
88 t4=(u30+u12);
89 t5=(u21+u03);
90
91 // 计算不变矩
92 M[0]=u20+u02;
93 M[1]=t1*t1+4*u11*u11;
94 M[2]=t2*t2+t3*t3;
95 M[3]=t4*t4+t5*t5;
96 M[4]=t2*t4*(t4*t4-3*t5*t5)+t3*t5*(3*t4*t4-t5*t5);
97 M[5]=t1*(t4*t4-t5*t5)+4*u11*t4*t5;
98 M[6]=t3*t4*(t4*t4-3*t5*t5)-t2*t5*(3*t4*t4-t5*t5);
99
100
101 /*cout<<M[0]<<endl;//<<二"<<M[0]<<"三"<<M[0]<<"四"<<M[0]<<"五"<<M[0]<<"六"<<M[0]<<"七"<<M[0]<<endl;
102 cout<<M[1]<<endl;
103 cout<<M[2]<<endl;
104 cout<<M[3]<<endl;
105 cout<<M[4]<<endl;
106 cout<<M[5]<<endl;
107 cout<<M[6]<<endl;
108 cout<<endl;*/
109 returntrue;
110 }
111
112
113 int main(char argc,char** argv)
114 {
115 int i;
116 double Sa[7] = {0},Ta[7] ={0};
117
118 ///*源图像
119 IplImage*img = cvLoadImage(filename,0);//灰度
120 HuMoment(img);
121 for(i=0;i<7;i++)
122 {
123 Sa[i] = M[i];
124 M[i] =0;
125 }
126 cout<<Sa[0]<<endl;
127 cout<<Sa[1]<<endl;
128 cout<<Sa[2]<<endl;
129 cout<<Sa[3]<<endl;
130 cout<<Sa[4]<<endl;
131 cout<<Sa[5]<<endl;
132 cout<<Sa[6]<<endl;
133 cout<<endl;
134 //*/
135
136
137 ///*模板图
138 IplImage*tpl = cvLoadImage(filename2,0);//灰度
139 HuMoment(tpl);
140 for(i=0;i<7;i++)
141 {
142 Ta[i] = M[i];
143 M[i] =0;
144 }
145 cout<<Ta[0]<<endl;
146 cout<<Ta[1]<<endl;
147 cout<<Ta[2]<<endl;
148 cout<<Ta[3]<<endl;
149 cout<<Ta[4]<<endl;
150 cout<<Ta[5]<<endl;
151 cout<<Ta[6]<<endl;
152 cout<<endl;
153
154
155 // 计算相似度
156 double dbR =0; //相似度
157 double dSigmaST =0;
158 double dSigmaS =0;
159 double dSigmaT =0;
160 double temp =0;
161
162 for(i=0;i<7;i++)
163 {
164 temp = Sa[i]*Ta[i];
165 dSigmaST+=temp;
166 dSigmaS+=pow(Sa[i],2);
167 dSigmaT+=pow(Ta[i],2);
168 }
169 dbR = dSigmaST/(sqrt(dSigmaS)*sqrt(dSigmaT));
170 printf("%lf\n",dbR);
171 //cout<<dbR<<endl;
172
173 cvReleaseImage(&img);
174 cvReleaseImage(&tpl);
175
176 return0;
177 }