OpenCV实现基于Zernike矩的亚像素边缘检测

   在做物体检测时，由于成本和应用场合的限制，不能够一味地增加相机的分辨率，或者已经用了分辨率很高的相机，但是视野范围很大，仍然无法实现很高的精度，这时就要考虑亚像素技术，亚像素技术就是在两个像素点之间进行进一步的细分，从而得到亚像素级别的边缘点的坐标（也就是float类型的坐标），一般来说，现有的技术可以做到2细分、4细分，甚至很牛的能做到更高，通过亚像素边缘检测技术的使用，可以节约成本，提高识别精度。

    常见的亚像素技术包括灰度矩、Hu矩、空间矩、Zernike矩等，通过查阅多篇论文，综合比较精度、时间和实现难度，选择Zernike矩作为项目中的亚像素边缘检测算法。基于Zernike矩的边缘检测原理，下一篇文章详细再总结一下，这里大体把步骤说一下，并使用OpenCV实现。

    大体步骤就是，首先确定使用的模板大小，然后根据Zernike矩的公式求出各个模板系数，例如我选择的是7*7模板（模板越大精度越高，但是运算时间越长），要计算M00，M11R，M11I，M20，M31R，M31I，M40七个Zernike模板。第二步是对待检测图像进行预处理，包括滤波二值化等，也可以在进行一次Canny边缘检测。第三步是将预处理的图像跟7个Zernike矩模板分别进行卷积，得到7个Zernike矩。第四步是把上一步的矩乘上角度校正系数（这是因为利用Zernike的旋转不变性，Zernike模型把边缘都简化成了x=n的直线，这里要调整回来）。第五步是计算距离参数l和灰度差参数k，根据k和l判断该点是否为边缘点。以下是基于opencv的实现。

[cpp]  view plain  copy

1. Mat M00 = (Mat_<double>(7, 7) <<  
2.     0, 0.0287, 0.0686, 0.0807, 0.0686, 0.0287, 0,  
3.     0.0287, 0.0815, 0.0816, 0.0816, 0.0816, 0.0815, 0.0287,  
4.     0.0686, 0.0816, 0.0816, 0.0816, 0.0816, 0.0816, 0.0686,  
5.     0.0807, 0.0816, 0.0816, 0.0816, 0.0816, 0.0816, 0.0807,  
6.     0.0686, 0.0816, 0.0816, 0.0816, 0.0816, 0.0816, 0.0686,  
7.     0.0287, 0.0815, 0.0816, 0.0816, 0.0816, 0.0815, 0.0287,  
8.     0, 0.0287, 0.0686, 0.0807, 0.0686, 0.0287, 0);  
9.   
10. Mat M11R = (Mat_<double>(7, 7) <<  
11.     0, -0.015, -0.019, 0, 0.019, 0.015, 0,  
12.     -0.0224, -0.0466, -0.0233, 0, 0.0233, 0.0466, 0.0224,  
13.     -0.0573, -0.0466, -0.0233, 0, 0.0233, 0.0466, 0.0573,  
14.     -0.069, -0.0466, -0.0233, 0, 0.0233, 0.0466, 0.069,  
15.     -0.0573, -0.0466, -0.0233, 0, 0.0233, 0.0466, 0.0573,  
16.     -0.0224, -0.0466, -0.0233, 0, 0.0233, 0.0466, 0.0224,  
17.     0, -0.015, -0.019, 0, 0.019, 0.015, 0);  
18.   
19. Mat M11I = (Mat_<double>(7, 7) <<  
20.     0, -0.0224, -0.0573, -0.069, -0.0573, -0.0224, 0,  
21.     -0.015, -0.0466, -0.0466, -0.0466, -0.0466, -0.0466, -0.015,  
22.     -0.019, -0.0233, -0.0233, -0.0233, -0.0233, -0.0233, -0.019,  
23.     0, 0, 0, 0, 0, 0, 0,  
24.     0.019, 0.0233, 0.0233, 0.0233, 0.0233, 0.0233, 0.019,  
25.     0.015, 0.0466, 0.0466, 0.0466, 0.0466, 0.0466, 0.015,  
26.     0, 0.0224, 0.0573, 0.069, 0.0573, 0.0224, 0);  
27.   
28. Mat M20 = (Mat_<double>(7, 7) <<  
29.     0, 0.0225, 0.0394, 0.0396, 0.0394, 0.0225, 0,  
30.     0.0225, 0.0271, -0.0128, -0.0261, -0.0128, 0.0271, 0.0225,  
31.     0.0394, -0.0128, -0.0528, -0.0661, -0.0528, -0.0128, 0.0394,  
32.     0.0396, -0.0261, -0.0661, -0.0794, -0.0661, -0.0261, 0.0396,  
33.     0.0394, -0.0128, -0.0528, -0.0661, -0.0528, -0.0128, 0.0394,  
34.     0.0225, 0.0271, -0.0128, -0.0261, -0.0128, 0.0271, 0.0225,  
35.     0, 0.0225, 0.0394, 0.0396, 0.0394, 0.0225, 0);  
36.   
37. Mat M31R = (Mat_<double>(7, 7) <<  
38.     0, -0.0103, -0.0073, 0, 0.0073, 0.0103, 0,  
39.     -0.0153, -0.0018, 0.0162, 0, -0.0162, 0.0018, 0.0153,  
40.     -0.0223, 0.0324, 0.0333, 0, -0.0333, -0.0324, 0.0223,  
41.     -0.0190, 0.0438, 0.0390, 0, -0.0390, -0.0438, 0.0190,  
42.     -0.0223, 0.0324, 0.0333, 0, -0.0333, -0.0324, 0.0223,  
43.     -0.0153, -0.0018, 0.0162, 0, -0.0162, 0.0018, 0.0153,  
44.     0, -0.0103, -0.0073, 0, 0.0073, 0.0103, 0);  
45.   
46. Mat M31I = (Mat_<double>(7, 7) <<  
47.     0, -0.0153, -0.0223, -0.019, -0.0223, -0.0153, 0,  
48.     -0.0103, -0.0018, 0.0324, 0.0438, 0.0324, -0.0018, -0.0103,  
49.     -0.0073, 0.0162, 0.0333, 0.039, 0.0333, 0.0162, -0.0073,  
50.     0, 0, 0, 0, 0, 0, 0,  
51.     0.0073, -0.0162, -0.0333, -0.039, -0.0333, -0.0162, 0.0073,  
52.     0.0103, 0.0018, -0.0324, -0.0438, -0.0324, 0.0018, 0.0103,  
53.     0, 0.0153, 0.0223, 0.0190, 0.0223, 0.0153, 0);  
54.   
55. Mat M40 = (Mat_<double>(7, 7) <<  
56.     0, 0.013, 0.0056, -0.0018, 0.0056, 0.013, 0,  
57.     0.0130, -0.0186, -0.0323, -0.0239, -0.0323, -0.0186, 0.0130,  
58.     0.0056, -0.0323, 0.0125, 0.0406, 0.0125, -0.0323, 0.0056,  
59.     -0.0018, -0.0239, 0.0406, 0.0751, 0.0406, -0.0239, -0.0018,  
60.     0.0056, -0.0323, 0.0125, 0.0406, 0.0125, -0.0323, 0.0056,  
61.     0.0130, -0.0186, -0.0323, -0.0239, -0.0323, -0.0186, 0.0130,  
62.     0, 0.013, 0.0056, -0.0018, 0.0056, 0.013, 0);  
63.   
64. vector<Point2d> SubEdgePoints;  
65.   
66. Mat ZerImageM00;  
67. filter2D(smoothImage, ZerImageM00, CV_64F, M00, Point(-1, -1), 0, BORDER_DEFAULT);  
68.   
69. Mat ZerImageM11R;  
70. filter2D(smoothImage, ZerImageM11R, CV_64F, M11R, Point(-1, -1), 0, BORDER_DEFAULT);  
71.   
72. Mat ZerImageM11I;  
73. filter2D(smoothImage, ZerImageM11I, CV_64F, M11I, Point(-1, -1), 0, BORDER_DEFAULT);  
74.   
75.   
76. Mat ZerImageM20;  
77. filter2D(smoothImage, ZerImageM20, CV_64F, M20, Point(-1, -1), 0, BORDER_DEFAULT);  
78.   
79.   
80. Mat ZerImageM31R;  
81. filter2D(smoothImage, ZerImageM31R, CV_64F, M31R, Point(-1, -1), 0, BORDER_DEFAULT);  
82.   
83. Mat ZerImageM31I;  
84. filter2D(smoothImage, ZerImageM31I, CV_64F, M31I, Point(-1, -1), 0, BORDER_DEFAULT);  
85.   
86. Mat ZerImageM40;  
87. filter2D(smoothImage, ZerImageM40, CV_64F, M40, Point(-1, -1), 0, BORDER_DEFAULT);  
88.   
89. int row_number = smoothImageAgain.rows;  
90. int col_number = smoothImageAgain.cols;  
91.   
92. for (int i = 0; i < row_number; i++)  
93. {  
94.     for (int j = 0; j <col_number; j++)  
95.     {  
96.         if (ZerImageM00.at<double>(i, j) == 0)  
97.         {  
98.             continue;  
99.         }  
100.           
101.         double theta_temporary = atan2(ZerImageM31I.at<double>(i, j), ZerImageM31R.at<double>(i, j));  
102.   
103.         //compute Z11'/Z31'  
104.         double rotated_z11 = 0.0;  
105.         rotated_z11 = sin(theta_temporary)*(ZerImageM11I.at<double>(i, j)) + cos(theta_temporary)*(ZerImageM11R.at<double>(i, j));  
106.         double rotated_z31 = 0.0;  
107.         rotated_z31 = sin(theta_temporary)*(ZerImageM31I.at<double>(i, j)) + cos(theta_temporary)*(ZerImageM31R.at<double>(i, j));  
108.   
109.         //compute l  
110.         double l_method1 = sqrt((5 * ZerImageM40.at<double>(i, j) + 3 * ZerImageM20.at<double>(i, j)) / (8 * ZerImageM20.at<double>(i, j)));  
111.         double l_method2 = sqrt((5 * rotated_z31 + rotated_z11) / (6 * rotated_z11));  
112.         double l = (l_method1 + l_method2) / 2;  
113.           
114.         //compute k/h  
115.         double k, h;          
116.         k = 3 * rotated_z11 / 2 / pow((1 - l_method2*l_method2), 1.5);  
117.         h = (ZerImageM00.at<double>(i, j) - k*PI / 2 + k*asin(l_method2) + k*l_method2*sqrt(1 - l_method2*l_method2)) / PI;  
118.   
119.         //judge the edge  
120.         double k_value = 20.0;  
121.         double l_value = sqrt(2) / g_N;  
122.           
123.         double absl = abs(l_method2 - l_method1);  
124.         if (k >= k_value && absl <= l_value)  
125.         {  
126.             Point2d point_temporary;  
127.             point_temporary.x = j + g_N*l*cos(theta_temporary) / 2;  
128.             point_temporary.y = i + g_N*l*sin(theta_temporary) / 2;  
129.             SubEdgePoints.push_back(point_temporary);  
130.         }  
131.         else  
132.         {  
133.             continue;  
134.         }  
135.     }  
136. }  
137. //show edge  
138. for (size_t i = 0; i < SubEdgePoints.size(); i++)  
139. {  
140. Point center_forshow(cvRound(SubEdgePoints[i].x), cvRound(SubEdgePoints[i].y));  
141. circle(SubImage, center_forshow, 1, Scalar(0, 97, 255), 1, 8, 0);  
142. }  
143. imshow("亚像素边缘", SubImage);