opencv實現一種改進的Fast特征檢測算法
阿新 • • 發佈：2017-05-26
sheng 特征檢測 local 溫習 現在 ble map 閾值 lag 
引言
之前了解了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 }View Code
4.算法效果
圖1是原始Fast檢測出來的特征點，點數較多；圖2是圖1經過非極大值抑制後的結果，觀察可發現局部非極大值點被剔除；圖3是圖2利用Hessian矩陣剔除邊緣點後的結果，在圖像下方邊緣處效果比較明顯。
5.參考文獻
[1]《基於FAST改進的快速角點探測算法》 燕鵬等
opencv實現一種改進的Fast特征檢測算法