http://www.cnblogs.com/polly333/p/5417353.html
Fast特征檢測
一、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;
}