一、Fast算法
1、基本原理
Fast特征點檢測feature2D原理是在圓周上按順時針方向從1到16的順序對圓周像素點進行編號。如果在圓周上有N個連續的像素的亮度都比圓心像素的亮度Ip加上閾值t還要亮，或者比圓心像素的亮度減去閾值還要暗，則圓心像素被稱為角點。
算法核心：利用周圍像素比較的信息可以得到特征點，簡單、高效。
FAST特征檢測算法來源於corner的定義，基於特征點周圍的像素灰度值。檢測候選特征點周圍一圈的像素值，如果候選區域內像素點足夠多且與候選點灰度值差值足夠大，則認為一個特征點。所以思路是：構建差值窗口，閾值選擇（點足夠多）
其中I(x)為圓周上任意一點的灰度，I(p)為圓心的灰度，Ed為灰度值差的閾值，如果N大於給定閾值，一般為周圍圓圈點的四分之三，則認為p是一個特征點。
在原理基礎上，為了提高運算速度，在計算時采用額外加速方法。
在點周圍每隔90度的四個點，如果有3個和候選點的灰度值值足夠大才認為此候選點為特征點候選點。如果不滿足此條件直接丟棄。程序中采用半徑為3，共有16（N）個周圍像素需要比較。FAST_9,FAST_10就是表示周圍像素個數。
2、算法流程
根據算法原理在設計程序時大體流程為：
1、設置閾值：用於比較是否周圍像素點和候選點的差值是否足夠大，閾值選擇很重要，也是一個缺陷
2、構建移動窗口：程序中設計為半徑為3，大約16個像素組成的區域，與中心點像素比較
3、候選像素與構建的周圍區域比較：算法采用先與圖中位置法檢查在位置1，9，5和13四個位置的像素，首先檢測位置1和位置9，如果它們都比閾值暗或比閾值亮，再檢測位置5和位置13。如果$P$是一個角點，那么上述四個像素點中至少有3個應該必須都大於$I_p+t$或者小於$I_p-t$，因為若是一個角點，超過四分之三圓的部分應該滿足判斷條件。如果滿足，則檢測圓內所有點。如果不滿足直接舍棄
4、對角點進行非極大值抑制，得到角點輸出。
以上方法還是有不夠魯棒的地方，但可以通過機器學習和非極大值抑制的方法來增強魯棒性。
1、計算得分函數，它的值V是特征點與其圓周上16個像素點的絕對差值中所有連續10個像素中的最小值的最大值，而且該值還要大於閾值t；
2、在3×3的特征點鄰域內（而不是圖像鄰域），比較V；
3、剔除掉非極大值的特征點
3、算法性質
- 通過周圍區域判斷四個角上點不能拒絕許多的候選點；
- 檢測出來的角點不是最優的，這是因為它的效率取決於問題的排序與角點的分布；
- 對於角點分析的結果被丟棄了；
- 多個特征點容易擠在一起
- 閾值選擇對結果有很大影響
二、算法源碼
注：此源碼來自於opencv，在此基礎進行分析和理解。
template<int patternSize> void FAST_t(InputArray _img, std::vector<KeyPoint>& keypoints, int threshold, bool nonmax_suppression) { Mat img = _img.getMat(); //提取出輸入圖像矩陣 //K為圓周連續像素的個數 //N用於循環圓周的像素點，因為要首尾連接，所以N要比實際圓周像素數量多K+1個 const int K = patternSize/2, N = patternSize + K + 1; #if CV_SSE2 const int quarterPatternSize = patternSize/4; (void)quarterPatternSize; #endif int i, j, k, pixel[25]; //找到圓周像素點相對於圓心的偏移量 makeOffsets(pixel, (int)img.step, patternSize); //特征點向量清零 keypoints.clear(); //保證閾值不大於255，不小於0 threshold = std::min(std::max(threshold, 0), 255); #if CV_SSE2 __m128i delta = _mm_set1_epi8(-128), t = _mm_set1_epi8((char)threshold), K16 = _mm_set1_epi8((char)K); (void)K16; (void)delta; (void)t; #endif // threshold_tab為閾值列表，在進行閾值比較的時候，只需查該表即可 uchar threshold_tab[512]; /*為閾值列表賦值，該表分為三段：第一段從threshold_tab[0]至threshold_tab[255 - threshold]，值為1，落在該區域的值表示滿足角點判斷條件2；第二段從threshold_tab[255 – threshold]至threshold_tab[255 + threshold]，值為0，落在該區域的值表示不是角點；第三段從threshold_tab[255 + threshold]至threshold_tab[511]，值為2，落在該區域的值表示滿足角點判斷條件1*/ for( i = -255; i <= 255; i++ ) threshold_tab[i+255] = (uchar)(i < -threshold ? 1 : i > threshold ? 2 : 0); //開辟一段內存空間 AutoBuffer<uchar> _buf((img.cols+16)*3*(sizeof(int) + sizeof(uchar)) + 128); uchar* buf[3]; /*buf[0、buf[1]和buf[2]分別表示圖像的前一行、當前行和后一行。因為在非極大值抑制的步驟2中，是要在3×3的角點鄰域內進行比較，因此需要三行的圖像數據。因為只有得到了當前行的數據，所以對於上一行來說，才湊夠了連續三行的數據，因此輸出的非極大值抑制的結果是上一行數據的處理結果*/ buf[0] = _buf; buf[1] = buf[0] + img.cols; buf[2] = buf[1] + img.cols; //cpbuf存儲角點的坐標位置，也是需要連續三行的數據 int* cpbuf[3]; cpbuf[0] = (int*)alignPtr(buf[2] + img.cols, sizeof(int)) + 1; cpbuf[1] = cpbuf[0] + img.cols + 1; cpbuf[2] = cpbuf[1] + img.cols + 1; memset(buf[0], 0, img.cols*3); //buf數組內存清零 //遍歷整幅圖像像素，尋找角點 //由於圓的半徑為3個像素，因此圖像的四周邊界都留出3個像素的寬度 for(i = 3; i < img.rows-2; i++) { //得到圖像行的首地址指針 const uchar* ptr = img.ptr<uchar>(i) + 3; //得到buf的某個數組，用於存儲當前行的得分函數的值V uchar* curr = buf[(i - 3)%3]; //得到cpbuf的某個數組，用於存儲當前行的角點坐標位置 int* cornerpos = cpbuf[(i - 3)%3]; memset(curr, 0, img.cols); //清零 int ncorners = 0; //檢測到的角點數量 if( i < img.rows - 3 ) { //每一行都留出3個像素的寬度 j = 3; #if CV_SSE2 if( patternSize == 16 ) { for(; j < img.cols - 16 - 3; j += 16, ptr += 16) { __m128i m0, m1; __m128i v0 = _mm_loadu_si128((const __m128i*)ptr); __m128i v1 = _mm_xor_si128(_mm_subs_epu8(v0, t), delta); v0 = _mm_xor_si128(_mm_adds_epu8(v0, t), delta); __m128i x0 = _mm_sub_epi8(_mm_loadu_si128((const __m128i*)(ptr + pixel[0])), delta); __m128i x1 = _mm_sub_epi8(_mm_loadu_si128((const __m128i*)(ptr + pixel[quarterPatternSize])), delta); __m128i x2 = _mm_sub_epi8(_mm_loadu_si128((const __m128i*)(ptr + pixel[2*quarterPatternSize])), delta); __m128i x3 = _mm_sub_epi8(_mm_loadu_si128((const __m128i*)(ptr + pixel[3*quarterPatternSize])), delta); m0 = _mm_and_si128(_mm_cmpgt_epi8(x0, v0), _mm_cmpgt_epi8(x1, v0)); m1 = _mm_and_si128(_mm_cmpgt_epi8(v1, x0), _mm_cmpgt_epi8(v1, x1)); m0 = _mm_or_si128(m0, _mm_and_si128(_mm_cmpgt_epi8(x1, v0), _mm_cmpgt_epi8(x2, v0))); m1 = _mm_or_si128(m1, _mm_and_si128(_mm_cmpgt_epi8(v1, x1), _mm_cmpgt_epi8(v1, x2))); m0 = _mm_or_si128(m0, _mm_and_si128(_mm_cmpgt_epi8(x2, v0), _mm_cmpgt_epi8(x3, v0))); m1 = _mm_or_si128(m1, _mm_and_si128(_mm_cmpgt_epi8(v1, x2), _mm_cmpgt_epi8(v1, x3))); m0 = _mm_or_si128(m0, _mm_and_si128(_mm_cmpgt_epi8(x3, v0), _mm_cmpgt_epi8(x0, v0))); m1 = _mm_or_si128(m1, _mm_and_si128(_mm_cmpgt_epi8(v1, x3), _mm_cmpgt_epi8(v1, x0))); m0 = _mm_or_si128(m0, m1); int mask = _mm_movemask_epi8(m0); if( mask == 0 ) continue; if( (mask & 255) == 0 ) { j -= 8; ptr -= 8; continue; } __m128i c0 = _mm_setzero_si128(), c1 = c0, max0 = c0, max1 = c0; for( k = 0; k < N; k++ ) { __m128i x = _mm_xor_si128(_mm_loadu_si128((const __m128i*)(ptr + pixel[k])), delta); m0 = _mm_cmpgt_epi8(x, v0); m1 = _mm_cmpgt_epi8(v1, x); c0 = _mm_and_si128(_mm_sub_epi8(c0, m0), m0); c1 = _mm_and_si128(_mm_sub_epi8(c1, m1), m1); max0 = _mm_max_epu8(max0, c0); max1 = _mm_max_epu8(max1, c1); } max0 = _mm_max_epu8(max0, max1); int m = _mm_movemask_epi8(_mm_cmpgt_epi8(max0, K16)); for( k = 0; m > 0 && k < 16; k++, m >>= 1 ) if(m & 1) { cornerpos[ncorners++] = j+k; if(nonmax_suppression) curr[j+k] = (uchar)cornerScore<patternSize>(ptr+k, pixel, threshold); } } } #endif for( ; j < img.cols - 3; j++, ptr++ ) { //當前像素的灰度值 int v = ptr[0]; //由當前像素的灰度值，確定其在閾值列表中的位置 const uchar* tab = &threshold_tab[0] - v + 255; //pixel[0]表示圓周上編號為0的像素相對於圓心坐標的偏移量 //ptr[pixel[0]表示圓周上編號為0的像素值 //tab[ptr[pixel[0]]]表示相對於當前像素（即圓心）圓周上編號為0的像素值在閾值列表threshold_tab中所查詢得到的值，如果為1，說明I0 < Ip - t，如果為2，說明I0 > Ip + t，如果為0，說明 Ip – t < I0 < Ip + t。因此通過tab，就可以得到當前像素是否滿足角點條件。 //編號為0和8（即直徑在圓周上的兩個像素點）在列表中的值相或后得到d。d=0說明編號為0和8的值都是0；d=1說明編號為0和8的值至少有一個為1，而另一個不能為2；d=2說明編號為0和8的值至少有一個為2，而另一個不能為1；d=3說明編號為0和8的值有一個為1，另一個為2。只可能有這四種情況。 int d = tab[ptr[pixel[0]]] | tab[ptr[pixel[8]]]; //d=0說明圓周上不可能有連續12個像素滿足角點條件，因此當前值一定不是角點，所以退出此次循環，進入下一次循環 if( d == 0 ) continue; //繼續進行其他直徑上兩個像素點的判斷 d &= tab[ptr[pixel[2]]] | tab[ptr[pixel[10]]]; d &= tab[ptr[pixel[4]]] | tab[ptr[pixel[12]]]; d &= tab[ptr[pixel[6]]] | tab[ptr[pixel[14]]]; //d=0說明上述d中至少有一個d為0，所以肯定不是角點；另一種情況是一個d為2，而另一個d為1，相與后也為0，這說明一個是滿足角點條件1，而另一個滿足角點條件2，所以肯定也不會有連續12個像素滿足同一個角點條件的，因此也一定不是角點。 if( d == 0 ) continue; //繼續判斷圓周上剩余的像素點 d &= tab[ptr[pixel[1]]] | tab[ptr[pixel[9]]]; d &= tab[ptr[pixel[3]]] | tab[ptr[pixel[11]]]; d &= tab[ptr[pixel[5]]] | tab[ptr[pixel[13]]]; d &= tab[ptr[pixel[7]]] | tab[ptr[pixel[15]]]; //如果滿足if條件，則說明有可能滿足角點條件2 if( d & 1 ) { //vt為真正的角點條件，即Ip – t，count為連續像素的計數值 int vt = v - threshold, count = 0; //遍歷整個圓周 for( k = 0; k < N; k++ ) { int x = ptr[pixel[k]]; //提取出圓周上的像素值 if(x < vt) //如果滿足條件2 { //連續計數，並判斷是否大於K（K為圓周像素的一半） if( ++count > K ) { //進入該if語句，說明已經得到一個角點 //保存該點的位置，並把當前行的角點數加1 cornerpos[ncorners++] = j; //進行非極大值抑制的第一步，計算得分函數 if(nonmax_suppression) curr[j] = (uchar)cornerScore<patternSize>(ptr, pixel, threshold); break; //退出循環 } } else count = 0; //連續像素的計數值清零 } } //如果滿足if條件，則說明有可能滿足角點條件1 if( d & 2 ) { //vt為真正的角點條件，即Ip + t，count為連續像素的計數值 int vt = v + threshold, count = 0; //遍歷整個圓周 for( k = 0; k < N; k++ ) { int x = ptr[pixel[k]]; //提取出圓周上的像素值 if(x > vt) //如果滿足條件1 { //連續計數，並判斷是否大於K（K為圓周像素的一半） if( ++count > K ) { //進入該if語句，說明已經得到一個角點 //保存該點的位置，並把當前行的角點數加1 cornerpos[ncorners++] = j; //進行非極大值抑制的第一步，計算得分函數 if(nonmax_suppression) curr[j] = (uchar)cornerScore<patternSize>(ptr, pixel, threshold); break; //退出循環 } } else count = 0; //連續像素的計數值清零 } } } } //保存當前行所檢測到的角點數 cornerpos[-1] = ncorners; //i=3說明只僅僅計算了一行的數據，還不能進行非極大值抑制的第二步，所以不進行下面代碼的操作，直接進入下一次循環 if( i == 3 ) continue; //以下代碼是進行非極大值抑制的第二步，即在3×3的角點鄰域內對得分函數的值進行非極大值抑制。因為經過上面代碼的計算，已經得到了當前行的數據，所以可以進行上一行的非極大值抑制。因此下面的代碼進行的是上一行的非極大值抑制。 //提取出上一行和上兩行的圖像像素 const uchar* prev = buf[(i - 4 + 3)%3]; const uchar* pprev = buf[(i - 5 + 3)%3]; //提取出上一行所檢測到的角點位置 cornerpos = cpbuf[(i - 4 + 3)%3]; //提取出上一行的角點數 ncorners = cornerpos[-1]; //在上一行內遍歷整個檢測到的角點 for( k = 0; k < ncorners; k++ ) { j = cornerpos[k]; //得到角點的位置 int score = prev[j]; //得到該角點的得分函數值 //在3×3的角點鄰域內，計算當前角點是否為最大值，如果是則壓入特性值向量中 if( !nonmax_suppression || (score > prev[j+1] && score > prev[j-1] && score > pprev[j-1] && score > pprev[j] && score > pprev[j+1] && score > curr[j-1] && score > curr[j] && score > curr[j+1]) ) { keypoints.push_back(KeyPoint((float)j, (float)(i-1), 7.f, -1, (float)score)); } } } }
在該函數內，對閾值列表理解起來可能有一定的難度，下面我們舉一個具體的例子來進行講解。設我們選取的閾值threshold為30，則根據
for( i = -255; i <= 255; i++ )
threshold_tab[i+255] = (uchar)(i < -threshold ? 1 : i > threshold? 2 : 0);
我們可以從-255到255一共分為3段：-255～-30，-30～30，30～255。由於數組的序號不能小於0，因此在給threshold_tab數組賦值上，序號要加上255，這樣區間就變為：0～225，225～285，285～510，而這三個區間對應的值分別為1，0和2。設我們當前像素值為40，則根據
const uchar* tab = &threshold_tab[0] -v + 255;
tab的指針指向threshold_tab[215]處，因為255-40=215。這樣在圓周像素與當前像素進行比較時，使用的是threshold_tab[215]以后的值。例如圓周上編號為0的像素值為5，則該值在閾值列表中的位置是threshold_tab[215 + 5]，是threshold_tab[220]。它在閾值列表中的第一段，即threshold_tab[220] = 1，說明編號為0的像素滿足角點條件2。我們來驗證一下：5 < 40 – 30，確實滿足條件2；如果圓周上編號為1的像素值為80，則該值在閾值列表中的位置是threshold_tab[295]（即215 + 80 = 295），而它在閾值列表中的第三段，即threshold_tab[295] = 2，因此它滿足角點條件1，即80 > 40 + 30；而如果圓周上編號為2的像素值為45，則threshold_tab[260] = 0，它不滿足角點條件，即40 – 30 < 45 < 40 + 30。
三、opencv函數解析
1、測試函數
void main() { Mat src; src = imread("D:/Demo.jpg"); // vector of keyPoints std::vector<KeyPoint> keyPoints; // construction of the fast feature detector object FastFeatureDetector fast(40); // 檢測的閾值為40 // feature point detection fast.detect(src,keyPoints); drawKeypoints(src, keyPoints, src, Scalar::all(-1), DrawMatchesFlags::DRAW_OVER_OUTIMG); imshow("FAST feature", src); cvWaitKey(0); }
2、函數解釋
在OpenCV中，當patternSize為16時，用以下數組表示這16個點相對於圓心的坐標：
static const int offsets16[][2] =
{
{0, 3}, { 1, 3}, { 2, 2}, { 3, 1}, { 3, 0}, { 3, -1}, { 2, -2}, { 1, -3},
{0, -3}, {-1, -3}, {-2, -2}, {-3, -1}, {-3, 0}, {-3, 1}, {-2, 2}, {-1, 3}
};
OpenCV用函數來計算圓周上的點相對於圓心坐標在原圖像中的位置：
void makeOffsets(int pixel[25], int rowStride, int patternSize) { //分別定義三個數組，用於表示patternSize為16，12和8時，圓周像素對於圓心的相對坐標位置 static const int offsets16[][2] = { {0, 3}, { 1, 3}, { 2, 2}, { 3, 1}, { 3, 0}, { 3, -1}, { 2, -2}, { 1, -3}, {0, -3}, {-1, -3}, {-2, -2}, {-3, -1}, {-3, 0}, {-3, 1}, {-2, 2}, {-1, 3} }; static const int offsets12[][2] = { {0, 2}, { 1, 2}, { 2, 1}, { 2, 0}, { 2, -1}, { 1, -2}, {0, -2}, {-1, -2}, {-2, -1}, {-2, 0}, {-2, 1}, {-1, 2} }; static const int offsets8[][2] = { {0, 1}, { 1, 1}, { 1, 0}, { 1, -1}, {0, -1}, {-1, -1}, {-1, 0}, {-1, 1} }; //根據patternSize值，得到具體應用上面定義的哪個數組 const int (*offsets)[2] = patternSize == 16 ? offsets16 : patternSize == 12 ? offsets12 : patternSize == 8 ? offsets8 : 0; CV_Assert(pixel && offsets); int k = 0; //代入輸入圖像每行的像素個數，得到圓周像素的絕對坐標位置 for( ; k < patternSize; k++ ) pixel[k] = offsets[k][0] + offsets[k][1] * rowStride; //由於要計算連續的像素，因此要循環的多列出一些值 for( ; k < 25; k++ ) pixel[k] = pixel[k - patternSize]; }
template<> int cornerScore<16>(const uchar* ptr, const int pixel[], int threshold) { const int K = 8, N = K*3 + 1; //v為當前像素值 int k, v = ptr[0]; short d[N]; //計算當前像素值與其圓周像素值之間的差值 for( k = 0; k < N; k++ ) d[k] = (short)(v - ptr[pixel[k]]); #if CV_SSE2 __m128i q0 = _mm_set1_epi16(-1000), q1 = _mm_set1_epi16(1000); for( k = 0; k < 16; k += 8 ) { __m128i v0 = _mm_loadu_si128((__m128i*)(d+k+1)); __m128i v1 = _mm_loadu_si128((__m128i*)(d+k+2)); __m128i a = _mm_min_epi16(v0, v1); __m128i b = _mm_max_epi16(v0, v1); v0 = _mm_loadu_si128((__m128i*)(d+k+3)); a = _mm_min_epi16(a, v0); b = _mm_max_epi16(b, v0); v0 = _mm_loadu_si128((__m128i*)(d+k+4)); a = _mm_min_epi16(a, v0); b = _mm_max_epi16(b, v0); v0 = _mm_loadu_si128((__m128i*)(d+k+5)); a = _mm_min_epi16(a, v0); b = _mm_max_epi16(b, v0); v0 = _mm_loadu_si128((__m128i*)(d+k+6)); a = _mm_min_epi16(a, v0); b = _mm_max_epi16(b, v0); v0 = _mm_loadu_si128((__m128i*)(d+k+7)); a = _mm_min_epi16(a, v0); b = _mm_max_epi16(b, v0); v0 = _mm_loadu_si128((__m128i*)(d+k+8)); a = _mm_min_epi16(a, v0); b = _mm_max_epi16(b, v0); v0 = _mm_loadu_si128((__m128i*)(d+k)); q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); v0 = _mm_loadu_si128((__m128i*)(d+k+9)); q0 = _mm_max_epi16(q0, _mm_min_epi16(a, v0)); q1 = _mm_min_epi16(q1, _mm_max_epi16(b, v0)); } q0 = _mm_max_epi16(q0, _mm_sub_epi16(_mm_setzero_si128(), q1)); q0 = _mm_max_epi16(q0, _mm_unpackhi_epi64(q0, q0)); q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 4)); q0 = _mm_max_epi16(q0, _mm_srli_si128(q0, 2)); threshold = (short)_mm_cvtsi128_si32(q0) - 1; #else //a0為閾值 int a0 = threshold; //滿足角點條件2時，更新閾值 for( k = 0; k < 16; k += 2 ) { //a為d[k+1]，d[k+2]和d[k+3]中的最小值 int a = std::min((int)d[k+1], (int)d[k+2]); a = std::min(a, (int)d[k+3]); //如果a小於閾值，則進行下一次循環 if( a <= a0 ) continue; //更新閾值 //a為從d[k+1]到d[k+8]中的最小值 a = std::min(a, (int)d[k+4]); a = std::min(a, (int)d[k+5]); a = std::min(a, (int)d[k+6]); a = std::min(a, (int)d[k+7]); a = std::min(a, (int)d[k+8]); //從d[k]到d[k+9]中的最小值與a0比較，哪個大，哪個作為新的閾值 a0 = std::max(a0, std::min(a, (int)d[k])); a0 = std::max(a0, std::min(a, (int)d[k+9])); } //滿足角點條件1時，更新閾值 int b0 = -a0; for( k = 0; k < 16; k += 2 ) { int b = std::max((int)d[k+1], (int)d[k+2]); b = std::max(b, (int)d[k+3]); b = std::max(b, (int)d[k+4]); b = std::max(b, (int)d[k+5]); if( b >= b0 ) continue; b = std::max(b, (int)d[k+6]); b = std::max(b, (int)d[k+7]); b = std::max(b, (int)d[k+8]); b0 = std::min(b0, std::max(b, (int)d[k])); b0 = std::min(b0, std::max(b, (int)d[k+9])); } threshold = -b0-1; #endif #if VERIFY_CORNERS testCorner(ptr, pixel, K, N, threshold); #endif //更新后的閾值作為輸出 return threshold; }