opencv實現一種改進的Fast特徵檢測算法

時間 2019-11-19

標籤 opencv 實現一種改進 fast 特徵檢測算法简体版

原文原文鏈接

引言html

　　以前瞭解了Fast算法以後使用opencv本身實現了下，具體見http://www.cnblogs.com/Wiley-hiking/p/6898049.html。不過算法也有缺點，主要就是對邊緣點和噪點的抗干擾能力不強。在《基於FAST改進的快速角點探測算法》文章中做者提出一種改進的Fast角點算法，提升算法的穩定性和抗干擾能力。本身讀完後使用opencv實現了該算法，這裏將學習過程進行一個記錄。ios

1.原始Fast檢測算法算法

有關原始Fast檢測算法，本身寫的小結在http://www.cnblogs.com/Wiley-hiking/p/6898049.html。數組

2.改進的Fast檢測算法ide

文章中指出，原始Fast檢測算法優勢是計算量小，缺點是（1）算法閾值須要人工設定，適應性很差（2）原始算法會提取獲得不少邊緣點或者局部非極大值點。針對上述問題，做者提出改進的Fast檢測算法，主要步驟包括：函數

（1）根據圖像性質自適應肯定Fast檢測算法中的閾值；做者提出使用圖像中像素灰度值前i個最大值和前i個最小值差值和，乘以一個係數結果做爲閾值。學習

（2）使用（1）獲得的閾值進行Fast檢測角點獲得候選點ui

（3）對於（2）獲得的候選點，利用Hessian矩陣去除邊緣點spa

（4）使用拉普拉斯極值排除局部非極大值點.net

　　如今詳細說明各步驟的具體實現。

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 }