在前面我們已經陸續介紹了許多特徵檢測運算元，我們可以根據影象區域性的自相關函式求得Harris角點，後面又提到了兩種十分優秀的特徵點以及他們的描述方法SIFT特徵和SURF特徵。SURF特徵是為了提高運算效率對SIFT特徵的一種近似，雖然在有些實驗環境中已經達到了實時，但是我們實踐工程應用中，特徵點的提取與匹配只是整個應用演算法中的一部分，所以我們對於特徵點的提取必須有更高的要求，從這一點來看前面介紹的的那些特徵點方法都不可取。

一 FAST演算法原理

為了解決這個問題，Edward Rosten和Tom Drummond在2006年發表的“Machine learning for high-speed corner detection”文章中提出了一種FAST特徵，並在2010年對這篇論文作了小幅度的修改後重新發表。FAST的全稱為Features From Accelerated Segment Test。

Rosten等人將FAST角點定義為：若某畫素點與其周圍領域內足夠多的畫素點處於不同的區域，則該畫素點可能為角點。也就是某些屬性與眾不同，考慮灰度影象，即若該點的灰度值比其周圍領域內足夠多的畫素點的灰度值大或者小，則該點可能為角點。

二 FAST演算法步驟

1. 從圖片中選取一個畫素$P$，下面我們將判斷它是不是一個特徵點。我們首先把它的亮度值設為$I_p$;
2. 設定一個合適的閾值$t$;
3. 考慮以該畫素點為中心的一個半徑等於3畫素的離散化的Bresenham圓，這個圓的邊界上有16個畫素；

4.現在，如果在這個大小為16個畫素的圓上有$n$個連續的畫素點，他們的畫素值要麼都比$I_p+t$大，要麼都比$I_p-t$小，那麼他就是一個角點。$n$的值可以設定為12或者9，實驗證明選擇9可能會有更好的效果。

上面的演算法中，對於影象中的每一個點，我們都要去遍歷其鄰域圓上的16個點的畫素，效率較低。我們下面提出了一種高效的測試（high-speed test）來快速排除一大部分非角點的畫素。該方法僅僅檢查在位置1，9，5和13四個位置的畫素，首先檢測位置1和位置9，如果它們都比閾值暗或比閾值亮，再檢測位置5和位置13。如果$p$是一個角點，那麼上述四個畫素點中至少有3個應該必須都大於$I_p+t$或者小於$I_p-t$，因為若是一個角點，超過四分之三圓的部分應該滿足判斷條件。如果不滿足，那麼$p$不可能是一個角點。對於所有點做上面這一部分初步的檢測後，符合條件的將成為候選的角點，我們再對候選的角點，做完整的測試，即檢測圓上的所有點。

上面的演算法效率實際上是很高的，但是有點一些缺點：

• 當我們設定$n<12$時就不能使用快速演算法來過濾非角點的點；
• 檢測出來的角點不是最優的，這是因為它的效率取決於問題的排序與角點的分佈；
• 對於角點分析的結果被丟棄了；
• 多個特徵點容易擠在一起。

三 使用機器學習做一個角點分類器

1. 首先選取你進行角點提取的應用場景下很多張的測試圖片;
2. 使用 FAST 演算法找出每幅影象的特徵點;
3. 對每一個特徵點，將其周圍的 16 個畫素儲存構成一個向量。對於步驟二中得到的角點，把他儲存在$P$中;
4. 對於影象上的每一個畫素點，它周圍鄰域圓上位置為$x$，$x\in\{1,2,...,16\}$的點表示為$p\to{x}$，可以用下面的判斷公式將該點$p\to{x}$分為三類:$$S_{p\to{x}}=\begin{cases}d,I_{p\to{x}}≤I_p-t \quad(darker)\\ s, I_p-t < I_{p\to{x}}<I_p+t \quad(similar)\\ b, I_p+t≤I_{p\to{x}} \quad(brighter)\end{cases}$$
5. 我們任意16個位置中的一個位置$x$，可以把集合$P$分為三個部分$P_d,P_s,P_b$，其中$P_d$的定義如下，$P_s$和$P_b$的定義與其類似：$$P_b=\{p\in{P}:S_{p\to{x}}=b\}$$ 換句話說，對於給定的位置$x$，它都是可以把所有影象中的點分為三類，第一類$P_d$包含了所有位置$x$處的畫素在閾值$t$下暗於中心畫素，$P_b$包含了所有位置$x$處的畫素在閾值$t$亮於中心元素；
6. 定義一個新的布林變數$K_p$，如果$p$是一個角點，那些$K_p$為真，否則為假；
7. 使用ID3演算法（決策樹分類器）來查詢每一個子集；
8. 遞迴計算所有的子集直至$K_p$的熵為0；
9. 被建立的決策樹就用於於其他圖片的FAST檢測。

 四 非極大值抑制

如何解決從鄰近的位置選取了多個特徵點的問題，我們可以使用Non-Maximal Suppression來解決。

1. 為每一個檢測到的特徵點計算它的響應大小(score function)$V$，這裡$V$定義為點$p$和它周圍16個畫素點的絕對偏差之和;
2. 考慮兩個相鄰的特徵點，並比較它們的$V$值;
3. $V$值較低的點會被刪除。

五 OpenCV庫FAST特徵檢測

# -*- coding: utf-8 -*-
"""
Created on Mon Aug 27 16:09:48 2018

@author: lenovo
"""

'''
FAST角點檢測
'''
import cv2

'''1、載入圖片'''
img1 = cv2.imread('./image/match1.jpg')
img1 = cv2.resize(img1,dsize=(600,400))
image1 = img1.copy()


'''2、提取特徵點'''
#建立一個FAST物件，傳入閾值t  可以處理RGB色彩空間影象  
fast = cv2.FastFeatureDetector_create(threshold=50)
keypoints1 = fast.detect(image1,None)
#在影象上繪製關鍵點
image1 = cv2.drawKeypoints(image=image1,keypoints = keypoints1,outImage=image1,color=(255,0,255),flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)

#輸出預設引數
print("Threshold: ", fast.getThreshold())
print("nonmaxSuppression: ", fast.getNonmaxSuppression())
print("neighborhood: ", fast.getType())
print("Total Keypoints with nonmaxSuppression: ", len(keypoints1))

#顯示影象
cv2.imshow('fast_keypoints1',image1)
cv2.waitKey(20)

#關閉非極大值抑制
fast.setNonmaxSuppression(0)
keypoints1 = fast.detect(image1,None)
print("Total Keypoints without nonmaxSuppression: ", len(keypoints1))
image1 = cv2.drawKeypoints(image=image1,keypoints = keypoints1,outImage=image1,color=(255,0,255),flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv2.imshow('fast_keypoints1 nms',image1)

cv2.waitKey(0)
cv2.destroyAllWindows()

程式執行效果如下：

可以看到經過非極大值抑制之後，特徵點從2092個降低到了1106個。如果你修改閾值$t$，你會發現$t$越大，檢測到的特徵點越小。

如果你還記得我們之前介紹SIFT特徵和SURF特徵，我們忽略演算法引數的影響，從總體上來看，你會發現FAST的特徵點數量遠遠多於前著。這是受多方面元素影響，一方面是受演算法本身影響，這兩個演算法是完全不同的；另一方面FAST特徵對噪聲比較敏感，從圖片上我們也可以觀察到，比如廣場上許多的的噪聲點。

除了上面我所說的這些，FAST演算法還有以下需要改進的地方：

• 由於FAST演算法依賴於一個閾值$t$，因此演算法還需要人為干涉；
• FAST演算法不產生多尺度特徵而且FAST特徵點沒有方向資訊，這樣就會失去旋轉不變性；

後面我會介紹ORB演算法，ORB將基於FAST關鍵點檢測的技術和基於BRIFE描述符的技術相結合，ORB演算法解決了上面我所講述到的缺點，可以用來替代SIFT和SURF演算法，與兩者相比，ORB擁有更快的速度。

六 自己實現FAST特徵檢測

由於FAST演算法比較簡單，因此我們可以按照我們之前所講述的步驟，自己去實現它，程式碼如下：

# -*- coding: utf-8 -*-
"""
Created on Mon Aug 27 20:22:51 2018

@author: lenovo
"""
import numpy as np
import cv2
from matplotlib import pyplot as plt

'''
自己實現FAST角點檢測演算法:不依賴OpenCV庫
參考程式碼：https://github.com/tbliu/FAST
'''
def rgb2gray(image):
    '''
    轉換圖片空間RGB->gray
    
    args:
        image：輸入RGB圖片資料        
    return：
        返回灰度圖片
    '''
    rows,cols = image.shape[:2]
    grayscale = np.zeros((rows,cols),dtype=np.uint8)
    for row in range(0,rows):
        for col in range(0,cols):
            red,green,blue = image[row][col]
            gray = int(0.3*red+0.59*green+0.11*blue)
            grayscale[row][col] = gray
    return grayscale

def bgr2gray(image):
    '''
    轉換圖片空間BGR->gray
    
    args:
        image：輸入BGR圖片資料        
    return：
        返回灰度圖片
    '''
    rows,cols = image.shape[:2]
    grayscale = image.copy()for row in range(0,rows):
        for col in range(0,cols):
            blue,green,red = image[row][col]
            gray = int(0.3*red+0.59*green+0.11*blue)
            grayscale[row][col] = gray
    return grayscale

def medianBlur(image,ksize=3,):
    '''
    中值濾波，去除椒鹽噪聲
    
    args:
        image：輸入圖片資料,要求為灰度圖片
        ksize：濾波視窗大小        
    return：
        中值濾波之後的圖片
    '''
    rows,cols = image.shape[:2]
    #輸入校驗
    half = ksize//2
    startSearchRow = half
    endSearchRow = rows-half-1
    startSearchCol = half
    endSearchCol = cols-half-1
    dst = np.zeros((rows,cols),dtype=np.uint8)
    #中值濾波
    for y in range(startSearchRow,endSearchRow):
        for x in range(startSearchCol,endSearchCol):
            window = []
            for i in range(y-half,y+half+1):
                for j in range(x-half,x+half+1):
                    window.append(image[i][j])
            #取中間值         
            window = np.sort(window,axis=None)                   
            if len(window)%2 == 1:            
                medianValue = window[len(window)//2]
            else:
                medianValue = int((window[len(window)//2]+window[len(window)//2+1])/2)
            dst[y][x] = medianValue
    return dst        
    
def circle(row,col):
    '''
    對於圖片上一畫素點位置(row,col)，獲取其鄰域圓上16個畫素點座標，圓由16個畫素點組成
    
    args:
        row：行座標 注意row要大於等於3
        col：列座標 注意col要大於等於3       
    '''
    if row < 3 or col < 3:
        return
    point1 = (row-3, col)
    point2 = (row-3, col+1)
    point3 = (row-2, col+2)
    point4 = (row-1, col+3)
    point5 = (row, col+3)
    point6 = (row+1, col+3)
    point7 = (row+2, col+2)
    point8 = (row+3, col+1)
    point9 = (row+3, col)
    point10 = (row+3, col-1)
    point11 = (row+2, col-2)
    point12 = (row+1, col-3)
    point13 = (row, col-3)
    point14 = (row-1, col-3)
    point15 = (row-2, col-2)
    point16 = (row-3, col-1)
    
    return [point1, point2,point3,point4,point5,point6,point7,point8,point9,point10,point11,point12, point13,point14,point15,point16]

def is_corner(image,row,col,threshold):
    '''
    檢測影象位置(row,col)處畫素點是不是角點
    如果圓上有12個連續的點滿足閾值條件，那麼它就是一個角點
    
    方法：
        如果位置1和9它的畫素值比閾值暗或比閾值亮，則檢測位置5和位置15
        如果這些畫素符合標準，請檢查畫素5和13是否相符
        如果滿足有3個位置滿足閾值條件，則它是一個角點
        重複迴圈函式返回的每個點如果沒有滿足閾值，則不是一個角落
        
        注意：這裡我們簡化了論文章中的角點檢測過程，會造成一些誤差
    
    args:
        image：輸入圖片資料,要求為灰度圖片
        row：行座標 注意row要大於等於3
        col：列座標 注意col要大於等於3 
        threshold：閾值        
    return : 
        返回True或者False
    '''
    #校驗
    rows,cols = image.shape[:2]
    if row < 3 or col < 3 :
        return False    
    if row >= rows-3 or col >= cols-3:
        return False    
    intensity = int(image[row][col])
    ROI = circle(row,col)
    #獲取位置1,9,5,13的畫素值
    row1, col1 = ROI[0]
    row9, col9 = ROI[8]
    row5, col5 = ROI[4]
    row13, col13 = ROI[12]
    intensity1 = int(image[row1][col1])
    intensity9 = int(image[row9][col9])
    intensity5 = int(image[row5][col5])
    intensity13 = int(image[row13][col13])
    #統計上面4個位置中滿足  畫素值  >  intensity + threshold點的個數
    countMore = 0
    #統計上面4個位置中滿足 畫素值  < intensity - threshold點的個數
    countLess = 0
    if intensity1 - intensity > threshold:
        countMore += 1 
    elif intensity1 + threshold < intensity:
        countLess += 1
    if intensity9 - intensity > threshold:
        countMore += 1
    elif intensity9 + threshold < intensity:
        countLess += 1
    if intensity5 - intensity > threshold:
        countMore += 1
    elif intensity5 + threshold < intensity:
        countLess += 1
    if intensity13 - intensity > threshold:
        countMore += 1
    elif intensity13 + threshold < intensity:
        countLess += 1
        
    return countMore >= 3 or countLess>=3

def areAdjacent(point1, point2):
    """
    通過尤拉距離來確定兩個點是否相鄰,如果它們在彼此的四個畫素內，則兩個點相鄰
    
    args:
        point1:畫素點1的位置
        point2:畫素點2的位置
     return : 
        返回True或者False
    """
    row1, col1 = point1
    row2, col2 = point2
    xDist = row1 - row2
    yDist = col1 - col2
    return (xDist ** 2 + yDist ** 2) ** 0.5 <= 4

def calculateScore(image,point):
    """ 
    計算非極大值抑制的分數
    為每一個檢測到的特徵點計算它的響應大小，得分V定義為p和它周圍16個畫素點的絕對偏差之和
    考慮兩個相鄰的特徵點，並比較它們的V,V值較小的點移除
    
    args:
        image：輸入圖片資料,要求為灰度圖片
        point: 角點座標
    """
    col, row  = point
    intensity = int(image[row][col])
    ROI = circle(row,col)
    values = []
    for p in ROI:
        values.append(int(image[p]))
    score = 0
    for value in values:
        score += abs(intensity - value)
    return score
    
def suppress(image, corners):
    '''
    Performs non-maximal suppression on the list of corners.
    For adjacent corners, discard the one with the smallest score.
    Otherwise do nothing

    Since we iterate through all the pixels in the image in order, any adjacent 
    corner points should be next to each other in the list of all corners

    Non-maximal suppression throws away adjacent corners which are the same point in real life
    
    args:
        image： is a numpy array of intensity values. NOTE: Image must be grayscale
        corners ： a list of (x,y) tuples   where x is the column index,and y is the row index
    '''
    i = 1
    #由於相鄰的角點在corners列表中彼此相鄰，所以我們寫成下面形式
    while i < len(corners):
        currPoint = corners[i]
        prevPoint = corners[i - 1]
        #判斷兩個角點是否相鄰
        if areAdjacent(prevPoint, currPoint):
            #計算非極大值抑制的分數
            currScore = calculateScore(image, currPoint)
            prevScore = calculateScore(image, prevPoint)
            #移除較小分數的點
            if (currScore > prevScore):
                del(corners[i - 1])
            else:
                del(corners[i])
        else:
            i += 1
            continue
    return


def detect(image, threshold=50,nonMaximalSuppress=True):
    '''
    corners = detect(image, threshold) performs the detection
    on the image and returns the corners as a list of (x,y) tuples
    where x is the column index, and y is the row index

    Nonmaximal suppression is implemented by default. 


    args： 
        image： is a numpy array of intensity values. NOTE: Image must be grayscale
        threshold：threshold is an int used to filter out non-corners. 
    return：
        returns the corners as a list of (x,y) tuples   where x is the column index,
        and y is the row index
    '''    
    corners = []
    rows,cols = image.shape[:2]
    #中值濾波
    image = medianBlur(image,3)
    cv2.imshow('medianBlur',image)   
    cv2.waitKey(20)
    #開始搜尋角點
    for row in range(rows):
        for col in range(cols):
            if is_corner(image, row, col, threshold):
                corners.append((col, row))  
    #非極大值抑制
    if nonMaximalSuppress:
        suppress(image, corners) 
    return corners;

def test():
    image = cv2.imread('./image/match1.jpg')
    image = cv2.resize(image,dsize=(600,400))
    imgray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    corners = detect(imgray)
    print('檢測到的角點個數為：',len(corners))     
    for point in corners:
        cv2.circle(image,point,1,(0,255,0),1)
    cv2.imshow('FAST',image)   
    cv2.waitKey(0)
    cv2.destroyAllWindows()
    
if __name__=='__main__':
    test()

執行之後效果如下：

我們在FAST演算法檢測之前使用了中值濾波，為了去除噪聲的影響，你也可以嘗試使用高斯濾波或者均指濾波等濾波手段，這裡就不在介紹，由於我們在實現濾波的時候並沒有對影象邊緣進行填充，因此在圖片四周並不會模糊。

右圖為檢測到的角點，大概有888個畫素點，這是使用了非極大值抑制之後的效果,與我們使用OpenCV庫檢測到的1106個近似相等.

參考文章：