FAST特徵點檢測的matlab實現
FAST特徵點檢測的matlab原始碼實現
1. 簡介
Features From Accelerate Segment Test (FAST) 是一種常用的特徵點檢測方法。它具有速度快，特徵點數量可控的優點。由FAST衍生出來的Oriented FAST and Rotated Brief (ORB)描述子也常常被應用於實時性較強的應用中。OPENCV中有關於FAST實現的原始碼，在matlab中也可以直接使用detectFASTFeatures()函式提取FAST特徵點。但是matlab中的我們無法看到FAST的實現方式，也無法對其進行優化改進。本文提供了一個關於FAST的matlab版本，幫助那些喜歡matlab的朋友們學習參考。為了讓整個博文具有完整性，我會簡單地介紹一下FAST的實現原理，便於與後面matlab的原始碼對應，方便讀者的理解與使用，如果熟悉FAST原理的小夥伴可以直接跳過原理部分講解，閱讀原始碼，也可以點選下面的連結下載我寫的FAST的matlab原始碼。
2. FAST的原理介紹
和大部分特徵點提取方法一樣，FAST特徵點提取主要分為兩個步驟：
- 特徵點檢測
- 極值點抑制
2.1 特徵點檢測
從原理的角度講，FAST並不像如Harris，SIFT，SURF那樣具有非常嚴格的數學背景，它更像是一種基於主觀經驗判斷的方法。其核心思想就是找到一些比周圍點的淺或者深的點作為特徵點，如圖1所示
圖1. FAST特徵點提取示意圖
FAST特徵點檢測的方法，就是比較檢測點的灰度值P0與檢測點附近的點的灰度值P1,P2,…,P16
剛才提到的FAST特徵點提取是要找到相比周圍點淺或者深的點，為了減小影象噪聲的干擾，FAST作者提出的方法是將P1,P2,…,P16看成一個迴圈列表，如果他們中有連續n個點大於P0 + threshold或者連續n個點小於P0 - threshold，則認為該檢測點是“可能”的特徵點（後面要進行極值點抑制，被抑制的點會從特徵點的候選名單裡剔除，為了便於表達，下文中我就不嚴格區別特徵點的candidate和真正的特徵點）。這裡的threshold主要還是為了防止由於影象噪聲引起的檢測誤差，可以是一個相對於P0
圖2. 連續9個及以上的點大於待檢測點。該檢測點為特徵點
假設圖2中，檢測點的灰度值P0 = 100,用橙色表示，threshold = 10，n = 9，紅色的點代表灰度值大於P0 + threshold (110)的點，黃色代表灰度值小於P0 - threshold(90)的點，灰色代表灰度值介於90到110之間的點。從圖中可以看到由於P13，P14，P15，P16，P1，P2，P3，P4，P5，P6這10個連續的點都大於P0 + threshold (110)，並且連續的點數大於等於n = 9，因此該檢測點為特徵點。同樣的方法可以得到圖3中有連續的9個點（P1到P9）小於P0 - threshold(90)，因此該檢測點也為特徵點。
圖3. 連續9個及以上的點小於待檢測點。該檢測點為特徵點
下面給出一個反例。如圖4所示，雖然P0周圍的點都滿足|Pi-P0|>thershold,但是沒有連續的n個點大於P0+threshold或者沒有連續的n個點小於P0-threshold，因此該檢測點不是特徵點。需要指出的是作者的原文的表述中有一定的奇異性。“n” contiguous pixels out of the 16 need to be either above or below P0 by the threshold, if the pixel needs to be detected as an interest point .[1] 這裡容易誤解成有n個連續的點滿足 |Pi-P0|>thershold，可以認為是特徵點。但是我在仔細研讀趙春江解讀的opencv的FAST原始碼後[2]，發現作者的本意是要求這n個連續的點必須都大於P0+threshold，或者都小於P0-threshold，才滿足特徵點判定條件。（當然特徵點檢測演算法是人定義的，並沒有嚴格的對與錯之分，這裡我只是儘可能地還原作者的本意。）
圖4. 非特徵點的案例
這裡有一個問題：為什麼要與圓周上的點的灰度值做比較？而不是直接與檢測點相鄰的每一個點進行比較呢？我認為可能是FAST的作者希望能予以特徵點一些相對冗餘的包容度（這一點恰好與基於一階導數差分或者二階導數差分判斷的特徵點方法相反）。影象特徵點其實並不是一個非常精確的概念，隨著光影，攝像頭焦距的變化，特徵點很可能會在原本位置附近“浮動”，因此跳開檢測點相鄰區域，只比較圓周上的點，可以忽略特徵點這種浮動帶來的影響。如果相鄰區域出現多個特徵點可以通過極值點抑制的方法來剔除。
上述過程是用於判定待定點是否為特徵點的方法，根據FAST特徵點的定義判定一個檢測點是否為特徵點，需要依次檢視周圍16個點的連續9個點是否大於P0 + threshold以及小於P0 - threshold的情況，每一個檢測點最糟情況是進行2169次判斷。如果這樣的操作需要遍歷全圖每一個點，顯然計算量也不小，FAST的作者提出了一個“偷懶”的方法。即優先判斷上下左右四個點P1，P5，P9，P13。並且作者還引入了一個閾值列表來加速該過程的判斷，具體細節可以參考趙春江解讀的opencv的FAST原始碼[2]。不過需要指出的是作者的這種優化思想主要是針對CPU逐點計算的優化思想（大部分點不是特徵點），對於FPGA，GPU等並行加速設計時，if語句的判斷操作並不適用。下面的matlab實現中，為了讓程式碼看起來更容易理解，我只注重原理的實現，並不追求時效性。
2.2 極值點抑制
極值點抑制是指把上述過程中提取出來可能的特徵點，通過一定的邏輯進行篩選。留下那些更具有代表性的點。Harris，SIFT，和SURF的極值點抑制是利用Hessian矩陣的行列式計算特徵點邊界特性，邊界特性強的點其特徵相似度較大，會被剔除。
FAST方法的特徵提取相對更加顯獲一些。FAST的作者把“更具代表性的點”理解為在一定區域內相比與周圍環境更淺或者更深的點。
如何讓演算法找到這些灰度更淺或者更深的點？首先就必須給這些特徵點進行打分，其次在一定的區域內進行篩選（分數最高的留下，分數低的淘汰）。
第一，FAST的打分方法是改變threshold，找到最大的一個threshold使該點恰好滿足特徵點判定的條件，即滿足n個以上連續的點大於P0 + threshold，或者P0-threshold*。這個threshold即為該點的得分值V = threshold。如果影象的畫素是UINT型別，V一定為圓周上某個與特徵點差值的絕對值減一，即V=|Pk-P0|-1，k為1到16中的某一個點。演算法實現過程中只需要遍歷這16個點找到滿足條件的最大的thershold*即可。
圖5. 特徵點打分
假設特徵點的灰度值P0 = 100，threshold = 10，當threhsold* = 11時，灰度值為111的點不再滿足大於P0 + threshold* 的條件。但是還是有9個連續的點滿足大於P0 + threshold*。當threshold* = 14時，該檢測點恰好不滿足特徵點的判定條件，因此threshold*的最大值為V = 14 -1 = 13。我們可以遍歷整張圖片如果該點不是特徵點V(i,j) = 0，如果該點為特徵點V(i,j) >= threshold。
第二，極值點抑制方法。剛才的方法可以得到一個特徵點打分V(i,j)的map，遍歷每一個可能的特徵點，如果該特徵點的V不是在以它為中心3x3的區域內最大值，則刪除該點，如果是則保留。
3. matlab原始碼實現
matlab的原始碼可以通過以下的連結下載：
https://download.csdn.net/download/qq_35721810/10867848
檔案中有一張png的測試圖片和三個.m檔案。如下圖所示
圖5 檔案說明
matlabFast.m是使用matlab中自帶的detectFASTFeatures()函式提取的FAST特徵點，可以用於對比實驗。myFAST.m是FAST特徵點檢測的實現。testMyFAST.m是FAST檢測的一個用例。讀者可以直接執行testMyFAST.m看到FAST原始碼的實現結果。
下面我主要針對myFAST.m進行解讀。
myFAST是一個實現FAST特徵點檢測的類，這個類中有四個properties
properties
numInChain = 12;
threshold = 0.3;
img;
featureMap;
end
numInChain代表連續點的數量，即上文中提到的n。threshold代表判斷特徵點超出的閾值，這裡我使用的是比例。img是影象，下面初始化的時候，我會把原始影象轉成灰度圖並儲存在img中。featureMap就是上文提到的V，它是一個與原圖一樣大小的map，如果該點不是特徵點featureMap(i, j) = 0，如果該點是特徵點featureMap(i, j) 即為該點的打分值。
methods中的第一個函式為myFAST的建構函式。
function obj = myFAST(img_)
if ndims(img_) == 3
obj.img = rgb2gray(img_);
else
obj.img = img_;
end
end
img_為原始影象，建構函式把原始影象轉為灰度圖，並存在類的img中。
函式operate()是FAST實現的主要函式。
% FAST的實現
function [corners] = operate(obj)
corners = []; % 特徵點
tempCorners = []; % 特徵點的candidate
[m, n] = size(obj.img);
obj.featureMap = zeros(m, n); %特徵點的打分值，如果不是特徵點打分為0
for i = 4:m-4
for j = 4:n-4
% p0為檢測點
p0 = obj.img(i,j);
% 依次提取當前點半徑為3的Bresenham圓圓周上的點，檢測點正上方為p（1）
p = [];
p(end + 1) = obj.img(i - 3 ,j ); %p1
p(end + 1) = obj.img(i - 3 ,j + 1); %p2
p(end + 1) = obj.img(i - 2 ,j + 2); %p3
p(end + 1) = obj.img(i - 1 ,j + 3); %p4
p(end + 1) = obj.img(i ,j + 3); %p5
p(end + 1) = obj.img(i + 1 ,j + 3); %p6
p(end + 1) = obj.img(i + 2 ,j + 2); %p7
p(end + 1) = obj.img(i + 3 ,j + 1); %p8
p(end + 1) = obj.img(i + 3 ,j ); %p9
p(end + 1) = obj.img(i + 3 ,j - 1); %p10
p(end + 1) = obj.img(i + 2 ,j - 2); %p11
p(end + 1) = obj.img(i + 1 ,j - 3); %p12
p(end + 1) = obj.img(i ,j - 3); %p13
p(end + 1) = obj.img(i - 1 ,j - 3); %p14
p(end + 1) = obj.img(i - 2 ,j - 2); %p15
p(end + 1) = obj.img(i - 3 ,j - 1); %p16
thr = obj.threshold*p0;
%判斷該點是否為特徵點的candidate
%featureFlag = 1 代表周圍有n個以上的點大於p0 + thr
%featureFlag = 2 代表周圍有n個以上的點小於p0 - thr
featureFlag = FastFeaturePointDetect(p0,p,thr,obj.numInChain);
if featureFlag == 1 || featureFlag == 2
%如果是特徵點，計算特徵點的打分值
val = ExtremeSuppression(p0,p,thr, obj.numInChain,featureFlag);
tempCorners(end + 1).x = j;
tempCorners(end).y = i;
tempCorners(end).val = val;
obj.featureMap(i,j) = val;
end
end
end
%極值點抑制
for k = 1:length(tempCorners)
j = tempCorners(k).x;
i = tempCorners(k).y;
tempList(1) = obj.featureMap(i - 1,j - 1);
tempList(2) = obj.featureMap(i - 1,j );
tempList(3) = obj.featureMap(i - 1,j + 1);
tempList(4) = obj.featureMap(i ,j - 1);
tempList(5) = obj.featureMap(i ,j + 1);
tempList(6) = obj.featureMap(i - 1,j - 1);
tempList(7) = obj.featureMap(i - 1,j );
tempList(8) = obj.featureMap(i - 1,j + 1);
%如果特徵點的打分值是附近3*3區域裡最大的則將tempCorners（k）付給corners
if tempCorners(k).val > max(tempList)
corners(end+1).x = tempCorners(k).x;
corners(end).y = tempCorners(k).y;
corners(end).val = tempCorners(k).val;
end
end
end
end
end
函式operate()中有兩個大迴圈，第一個迴圈會遍歷全圖（由於Bresenham圓半徑為3，因此中心點從（4，4）開始），依次判斷圖中每一個畫素是否為特徵點的candidate，如果是，則將該點的座標儲存下來，並且為該特徵點打分，並且記錄在featureMap中。第二個迴圈遍歷每一個特徵點的candidate，判斷是否是3x3區域中打分值最大的。如果是，則保留存入corners裡，如果不是，則刪除。
第一個迴圈中有兩個主要的函式: FastFeaturePointDetect()和ExtremeSuppression()。FastFeaturePointDetect()是用來判斷該點是否為特徵點；ExtremeSuppression()是用來給特徵點打分的（這裡名命為ExtremeSuppression極值點抑制，可能寫程式碼時有欠考慮，但其實這個函式的作用只是為特徵點打分，抑制是operate()的第二個大迴圈實現的）
%判斷該點是否為特徵點的candidate
function flag = FastFeaturePointDetect(p0,p,thr, numInChain)
flag = 0;
% 遍歷16個點，看是否有numInChain個點大於p0 + thr
% 迴圈列表的index = mod(k-1,length(p))+1
for i = 1:16
if p(i) - p0 > thr
counter = 1;
for j = 1:numInChain - 1
k = i + j;
if p(mod(k-1,length(p))+1) - p0 > thr
counter = counter + 1;
else
break;
end
end
if counter == numInChain
flag = 1;
return
end
else
continue;
end
end
% 遍歷16個點，看是否有numInChain個點小於p0 - thr
for i = 1:16
if p0 - p(i) > thr
counter = 1;
for j = 1:numInChain - 1
k = i + j;
if p0 - p(mod(k-1,length(p))+1) > thr
counter = counter + 1;
else
break;
end
end
if counter == numInChain
flag = 2;
return
end
else
continue;
end
end
end
下面一次給大家講解上述的兩個函式。首先FastFeaturePointDetect(p0,p,thr, numInChain)的輸入p0代表待檢測點，即上文中的P0，p代表p0周圍的16個點，thr代表threshold，numInChain代表n個連續的點數。這個函式中一共有兩個二個迴圈，第一個迴圈用來判斷是否有numInChain個點大於p0 + thr；第二個迴圈用來判斷是否有numInChain個點小於p0 - thr。這裡為了讓實現更加簡單易懂，我並沒有嚴格按照FAST作者推薦的“偷懶”方法去實現（我猜用matlab的朋友們並不關心演算法的執行速度）。這裡有一個使用迴圈列表的小技巧，在matlab中迴圈列表的index可以根據index = mod(k-1,length( p ))+1來實現，用一個例子說明，假設length( p ) = 16，k為17，k此時已經超出了p的長度，可以計算出mod(17 - 1,16)+1 = 1，因此index正好指向p的第一元素。這裡的featureFlag為該函式的輸出結果，如果結果為1代表有numInChain個點大於p0 + thr，如果結果為2代表有numInChain個點小於p0 - thr，如果結果為0代表該點不是特徵點。
%給特徵點打分
function [val] = ExtremeSuppression(p0, p,thr, numInChain,featureFlag)
N = length(p);
if featureFlag == 1
d = p - double(p0);
else
d = double(p0) - p;
end
a0 = thr;
%打分的結果一定是d（1）到d(16)中的一個，遍歷16次找出最大的滿足條件的d(k)
for i = 1:16
a = min([d(mod(i,N)+1),d(mod(i+1,N)+1),d(mod(i+2,N)+1)]);
if a < a0
continue;
end
for j = 4:numInChain-1
a = min(a,d(mod(i+j-1,N)+1));
end
a0 = max(a0,min(a, d(mod(i-1,N)+1) ));
a0 = max(a0, min(a, d(mod(i-1,N)+1) ));
end
val = a0 -1;
end
特徵點的打分用ExtremeSuppression(p0, p,thr, numInChain,featureFlag)函式實現。它的基本思想就是打分值一定是這16個點中其中一個，因此遍歷這16個點，如果這個點之後的連續9個點中最小的點，大於thr，則將其記錄下來。滿足上面條件最大的點的灰度值減一，即為所求結果。
4. 結果展示
最後，是我的演算法和matlab自帶函式的效果對比圖，總體效果上看還是差不多的，但是個人感覺matlab裡detectFASTFeatures()的特徵點提取還是比我自己實現的效果要好一點吧，可能裡面的極值點抑制與我的方法有一些差異。
圖6. myFAST的實現結果
圖7. matlab函式的實現結果
[1]: Features from Accelerated Segment Test, Deepak Geetha Viswanathan
[2]: https://blog.csdn.net/zhaocj/article/details/40301561