引言
之前了解了Fast算法之后使用opencv自己实现了下，具体见http://www.cnblogs.com/Wiley-hiking/p/6898049.html。不过算法也有缺点，主要就是对边缘点和噪点的抗干扰能力不强。在《基于FAST改进的快速角点探测算法》文章中作者提出一种改进的Fast角点算法，提高算法的稳定性和抗干扰能力。自己读完后使用opencv实现了该算法，这里将学习过程进行一个记录。
1.原始Fast检测算法
有关原始Fast检测算法，自己写的小结在http://www.cnblogs.com/Wiley-hiking/p/6898049.html。
2.改进的Fast检测算法
文章中指出，原始Fast检测算法优点是计算量小，缺点是（1）算法阈值需要人工设定，适应性不好（2）原始算法会提取得到很多边缘点或者局部非极大值点。针对上述问题，作者提出改进的Fast检测算法，主要步骤包括：
（1）根据图像性质自适应确定Fast检测算法中的阈值；作者提出使用图像中像素灰度值前i个最大值和前i个最小值差值和，乘以一个系数结果作为阈值。
（2）使用（1）得到的阈值进行Fast检测角点得到候选点
（3）对于（2）得到的候选点，利用Hessian矩阵去除边缘点
（4）使用拉普拉斯极值排除局部非极大值点
现在详细说明各步骤的具体实现。
2.1自适应确定阈值
这里需要得到图像像素中前10个最大值和前10个最小值，直观感受就是——这是一个排序问题呀。从网上搜索一些有关排序的文章温习下，确定使用下述方法实现
（1）得到前10个最大值；首先读取图像中前10个像素灰度值保存到数组并进行降序排列（升序排列也可以，使用改进的冒泡法实现，参考http://blog.csdn.net/cbs612537/article/details/8294960/）；之后依次读取剩下的像素，每读入一个像素，做如下操作：将该像素灰度值与前面10个像素组合并重新降序排列并将最小值剔除，得到新的降序排列数组（仍然包含10个像素）；重复上述步骤直至读取完所有像素，最终的数组包含的10个像素灰度值就是前10个最大值
（2）得到前10个最小值；步骤与上述类似，只不过将新的像素灰度值与前面10个像素组合重新排序后剔除最大值。
这里我在将新像素合入数组再剔除最大值（最小值）的过程等效为新元素插入的过程。由于插入前数组就是有序的，我就使用了二分法进行插入并更新数组，提高算法效率。得到前10个最大值和前10个最小值之后按照上述方法得到Fast特征检测中需要用到的threshold。
2.2Fast特征检测
带入2.1中得到的threshold进行计算，具体参考http://www.cnblogs.com/Wiley-hiking/p/6898049.html
2.3使用hessian矩阵去除边缘点
有关hessian矩阵的介绍和应用，参考http://blog.csdn.net/lwzkiller/article/details/55050275。我们这里只需要知道，角点的Hessian矩阵两个特征值比较相近；而利用矩阵的特性我们可以在不求得特征值的情况下计算两个特征值的比值，这里直接给出应用结论如下。
而我们是在离散的灰度值点上计算hessian矩阵，二阶导数的计算相对简化。有关离散函数二阶导数计算方面参考http://blog.csdn.net/xiaofengsheng/article/details/6023368
2.4使用拉普拉斯极值排除非局部极大值点
由于Fast算法中可以利用膨胀与比较的组合实现非极大值的抑制，这里我就没有使用该步骤
3.opencv代码实现
开发环境是vs2012+opencv2.4.13，源代码贴出来
1 #include <iostream>
2 #include <core/core.hpp>
3 #include <highgui/highgui.hpp>
4 #include <imgproc/imgproc.hpp>
5 #include <features2d/features2d.hpp>
6 #include <stdlib.h>
7
8 using namespace std;
9 using namespace cv;
10
11
12 int getSum(uchar *p, int length)
13 {
14 int sum = 0;
15 for(int i=0;i<length; i++)
16 {
17 sum += *(p+i);
18 }
19 return sum;
20 }
21
22 void bubbleSortUp(int *p)
23 {
24 uchar flag = 1;
25 for(int i=0; i < 10 && flag; i++)
26 {
27 flag = 0;
28 for(int j=9; j > i; j--)
29 {
30 if(p[j] < p[j-1])
31 {
32 uchar tmp = p[j];
33 p[j] = p[j-1];
34 p[j-1] = tmp;
35 flag = 1;
36 }
37 }
38 }
39 }
40
41 void bubbleSortDown(int *p)
42 {
43 uchar flag = 1;
44 for(int i=0; i < 10 && flag; i++)
45 {
46 flag = 0;
47 for(int j=9; j > i; j--)
48 {
49 if(p[j] > p[j-1])
50 {
51 uchar tmp = p[j];
52 p[j] = p[j-1];
53 p[j-1] = tmp;
54 flag = 1;
55 }
56 }
57 }
58 }
59
60 void printArray(int *p, int len)
61 {
62 for(int i=0; i<len; i++)
63 {
64 cout<<p[i]<<" ";
65 }
66 cout<<endl;
67 }
68
69 void binaryInsertUp(int *p, uchar value)
70 {
71 int left = 0;
72 int right = 9;
73 int mid;
74 if(value < p[0])
75 return;
76 if(value > p[9])
77 {
78 for(int i=0; i<9; i++)
79 {
80 p[i] = p[i+1];
81
82 }
83 p[9] = value;
84 return;
85 }
86
87 while(left < right)
88 {
89 mid = (left + right)/2;
90 if(mid == left || mid == right)
91 break;
92 if(value < p[mid])
93 right = mid;
94 else
95 left = mid;
96 }
97 for(int i = 0; i < left; i++)
98 {
99 p[i] = p[i+1];
100 }
101 p[left] = value;
102 }
103
104 void binaryInsertDown(int *p, uchar value)
105 {
106 int left = 0;
107 int right = 9;
108 int mid;
109 if(value > p[0])
110 return;
111 if(value < p[9])
112 {
113 for(int i=0; i<9; i++)
114 {
115 p[i] = p[i+1];
116
117 }
118 p[9] = value;
119 return;
120 }
121
122 while(left < right)
123 {
124 mid = (left + right)/2;
125 if(mid == left || mid == right)
126 break;
127 if(value < p[mid])
128 left = mid;
129 else
130 right = mid;
131 }
132 for(int i = 0; i < left; i++)
133 {
134 p[i] = p[i+1];
135 }
136 p[left] = value;
137 }
138
139 int main(int argc, char* argv[])
140 {
141 /* 1.读入图像 */
142 Mat image = imread("../church01.jpg", 0);
143 if(!image.data)
144 return 0;
145
146 namedWindow("Original Image");
147 imshow("Original Image", image);
148
149 Mat fastImg(image.size(), CV_8U, Scalar(0));//用于保存Fast检测的候选点
150 Mat fastScore(image.size(), CV_32F, Scalar(0));//用于计算候选点score
151 vector<Point> points;
152 int rows, cols, threshold;
153 rows = image.rows;
154 cols = image.cols;
155 threshold = 50;
156 int maxValues[10], minValues[10];
157
158 /* 2.计算Fast算法中的阈值参数 */
159 Mat_<uchar>::const_iterator it = image.begin<uchar>();
160 int count = 0;
161 while(count < 10)
162 {
163 maxValues[count] = *it;
164 minValues[count] = *it;
165 count++;
166 it++;
167 }
168 bubbleSortUp(maxValues);
169 bubbleSortDown(minValues);
170
171 while(it != image.end<uchar>())
172 {
173 int value = *it;
174 binaryInsertUp(maxValues, value);
175 binaryInsertDown(minValues, value);
176 it++;
177 }
178 printArray(maxValues, 10);
179 printArray(minValues, 10);
180
181 int diff = 0;
182 for(int i = 0; i < 10; i++)
183 {
184 diff += maxValues[i] - minValues[i];
185 }
186 threshold = (int)(0.15f*diff/10.0);
187
188 /* 3.使用Fast算法检测得到候选点 */
189 for(int x = 3; x < rows - 3; x++)
190 {
191 for(int y = 3; y < cols - 3; y++)
192 {
193 uchar delta[16] = {0};
194 uchar diff[16] = {0};
195 delta[0] = (diff[0] = abs(image.at<uchar>(x,y) - image.at<uchar>(x, y-3))) > threshold;
196 delta[8] = (diff[8] = abs(image.at<uchar>(x,y) - image.at<uchar>(x, y+3))) > threshold;
197 if(getSum(delta, 16) == 0)
198 continue;
199 else
200 {
201 delta[12] = (diff[12] = abs(image.at<uchar>(x,y) - image.at<uchar>(x-3, y))) > threshold;
202 delta[4] = (diff[4] = abs(image.at<uchar>(x,y) - image.at<uchar>(x+3, y))) > threshold;
203 if(getSum(delta, 16) < 3)
204 continue;
205
206 else
207 {
208 delta[1] = (diff[1] = abs(image.at<uchar>(x,y) - image.at<uchar>(x+1, y-3))) > threshold;
209 delta[2] = (diff[2] = abs(image.at<uchar>(x,y) - image.at<uchar>(x+2, y-2))) > threshold;
210 delta[3] = (diff[3] = abs(image.at<uchar>(x,y) - image.at<uchar>(x+3, y-1))) > threshold;
211
212 delta[5] = (diff[5] = abs(image.at<uchar>(x,y) - image.at<uchar>(x+3, y+1))) > threshold;
213 delta[6] = (diff[6] = abs(image.at<uchar>(x,y) - image.at<uchar>(x+2, y+2))) > threshold;
214 delta[7] = (diff[7] = abs(image.at<uchar>(x,y) - image.at<uchar>(x+1, y+3))) > threshold;
215
216 delta[9] = (diff[9] = abs(image.at<uchar>(x,y) - image.at<uchar>(x-1, y+3))) > threshold;
217 delta[10] = (diff[10] = abs(image.at<uchar>(x,y) - image.at<uchar>(x-2, y+2))) > threshold;
218 delta[11] = (diff[11] = abs(image.at<uchar>(x,y) - image.at<uchar>(x-3, y+1))) > threshold;
219
220 delta[13] = (diff[13] = abs(image.at<uchar>(x,y) - image.at<uchar>(x-3, y-1))) > threshold;
221 delta[14] = (diff[14] = abs(image.at<uchar>(x,y) - image.at<uchar>(x-2, y-2))) > threshold;
222 delta[15] = (diff[15] = abs(image.at<uchar>(x,y) - image.at<uchar>(x-1, y-3))) > threshold;
223
224 if(getSum(delta, 16) >= 12)
225 {
226 points.push_back(Point(y,x));
227 fastScore.at<float>(x,y) = getSum(diff, 16);
228 fastImg.at<uchar>(x,y) = 255;
229 }
230 }
231 }
232 }
233
234 }
235
236 vector<Point>::const_iterator itp = points.begin();
237 while(itp != points.end())
238 {
239 circle(image, *itp, 3, Scalar(255), 1);
240 ++itp;
241 }
242 namedWindow("Fast Image");
243 imshow("Fast Image", image);//Fast检测候选点在原图中标记
244
245 /* 4.局部非极大值抑制 */
246 image = imread("../church01.jpg", 0);
247 Mat dilated(fastScore.size(), CV_32F, Scalar(0));
248 Mat localMax;
249 //Mat element(7, 7, CV_8U, Scalar(1));
250 dilate(fastScore, dilated, Mat());
251 compare(fastScore, dilated, localMax, CMP_EQ);
252 bitwise_and(fastImg, localMax, fastImg);
253
254 for(int x = 0;x < fastImg.cols; x++)
255 {
256 for(int y = 0; y < fastImg.rows; y++)
257 {
258 if(fastImg.at<uchar>(y,x))
259 {
260 circle(image, Point(x,y), 3, Scalar(255), 1);
261 }
262 }
263 }
264 namedWindow("Fast Image2");
265 imshow("Fast Image2", image);//局部非极大值抑制后候选点在原图中标记
266
267 /* 5.计算Hessian矩阵去除边缘点和不稳定点 */
268 image = imread("../church01.jpg", 0);
269 Mat cornerStrength(fastScore.size(), CV_64F, Scalar(0));
270 for(int x = 0;x < fastImg.rows; x++)
271 {
272 for(int y = 0; y < fastImg.cols; y++)
273 {
274 if(fastImg.at<uchar>(x,y))
275 {
276 if((x>0) && (x<fastImg.rows-1)
277 && (y>0) && (y<fastImg.cols-1))
278 {
279 double Ixx = image.at<uchar>(x+1,y) + image.at<uchar>(x-1,y)-2*image.at<uchar>(x,y);
280 double Iyy = image.at<uchar>(x,y+1) + image.at<uchar>(x,y-1)-2*image.at<uchar>(x,y);
281 double Ixy = image.at<uchar>(x+1,y+1) + image.at<uchar>(x,y)-image.at<uchar>(x,y+1)-image.at<uchar>(x+1,y);
282 cornerStrength.at<double>(x,y) = (Ixx+Iyy)*(Ixx+Iyy)/(Ixx*Iyy-Ixy*Ixy);
283
284 }
285 }
286 }
287 }
288 int t = 10;
289 Mat cornerMap;
290 cornerMap = cornerStrength > t;
291 image = imread("../church01.jpg", 0);
292
293 for(int x = 0;x < cornerMap.cols; x++)
294 {
295 for(int y = 0; y < cornerMap.rows; y++)
296 {
297 if(cornerMap.at<uchar>(y,x))
298 {
299 circle(image, Point(x,y), 3, Scalar(255), 1);
300
301 }
302 }
303 }
304 namedWindow("improvedFast");
305 imshow("improvedFast", image);//最终检测得到的角点在原图中标记
306
307 waitKey();
308 return 0;
309 }
4.算法效果
图1是原始Fast检测出来的特征点，点数较多；图2是图1经过非极大值抑制后的结果，观察可发现局部非极大值点被剔除；图3是图2利用Hessian矩阵剔除边缘点后的结果，在图像下方边缘处效果比较明显。
5.参考文献
[1]《基于FAST改进的快速角点探测算法》 燕鹏等
时间： 2024-08-06 11:52:52