ransac算法思想

Ransac： Random Sample Consensus, 隨機抽樣一致性。RANSAC算法在1981年由Fischler和Bolles首次提出。

Ransac是一種通過使用觀測到的數據點來估計數學模型參數的迭代方法。其中數據點包括內點(inlier)，外點(outlier)。outlier對模型的估計沒有價值，因此該方法也可以叫做outlier檢測方法。其中inliners指樣本中正確的的樣本點，常稱為內點，內群點，正確樣本，outliners指樣本中的噪聲點(如測量時誤差很大的點，數據值不准確的點等)，常稱為外點，離群點，噪聲點。

所以通俗點將，就是我們有一組觀測數據，但是數據點中有很多噪聲點，如果我們直接用模型來擬合數據，模型肯定會受噪聲點影響，而且噪聲點比例越大，擬合出的模型也會越不准確。Ransac就是用來解決樣本中的外點問題(噪聲)，最多可處理50%的外點情況。

1. Ransac的假設和思想

再深入點討論，有一組觀測數據，要建立一個模型來擬合數據，我們第一步要做的，肯定是用一個標准對數據進行篩選，去除噪聲點，讓數據盡量干凈和准確。但是如果我們沒有一個合適的篩選標准，該怎么辦呢？我們可以假設：觀測數據中除了外點(噪聲)，肯定存在內點(准確點)，而且內點的個數夠我們用來擬合一個模型(比如擬合一條直線，至少要兩個數據點)。有了這個假設，我們可以從觀測數據中隨機挑選出n個數據點，用這n個點來擬合模型，然后對這個模型進行評價。如果通過迭代，反復重復這個過程，只要迭代次數過大，我們隨機挑選出來的n個樣本點是有可能全部是內點的，再配合上模型評價，就能找到最優的擬合模型。這就是ransac的主要假設和思想，官方敘述如下：

• RANSAC的基本假設是 “內群”數據可以通過幾組模型參數來敘述其數據分布，而“離群”數據則 是不適合模型化的數據。 數據會受噪聲影響, 噪聲指的是離群：例如從極端的噪聲, 錯誤解釋有關數據的測量, 或不正確的假設。RANSAC假定，給定一組（通常很小的）內群，存在一個 程序，這個程序可以估算最佳解釋或最適用於這一數據模型的參數。

宏觀點來看，Ransac是一種思想，一種利用迭代來增加概率的方法，即：

• RANSAC是一種思想，一個求解已知模型的參數的框架。它不限定某一特定的問題，可以是計算機視覺的問題，同樣也可以是統計數學，甚至可以是經濟學領域的模型參數估計問題。
• RANSAC是一種迭代的方法，用來在一組包含離群的被觀測數據中估算出數學模型的參數。 RANSAC 是一個非確定性算法，在某種意義上說，它會產生一個在一定概率下合理的結果，其允許使用更多次的迭代來使其概率增加。

2. Ransac基本流程和參數設置

2.1 Ransac步驟

RANSAC算法的輸入:

• 1.一組觀測數據（往往含有較大的噪聲或無效點
• 2.一個用於解釋觀測數據的參數化模型
• 3.一些可信的參數

迭代過程：

1. 在數據中隨機選擇n個點設定為內群
2. 計算適合內群的模型，如線性直線模型y=ax+b
3. 把其它剛才沒選到的點帶入剛才建立的模型中，計算是否為內群點
4. 記下內群數量
5. 重復以上步驟, 迭代k次
6. 比較哪次計算中內群數量最多,內群最多的那次所建的模型就是我們所要求的解

注意：不同問題對應的數學模型不同，因此在計算模型參數時方法必定不同，RANSAC的作用不在於 計算模型參數。（這導致ransac的缺點在於要求數學模型已知）

2.2 Ransac參數設置

在上面的步驟中，有兩個參數比較重要，需要提前確定，即：

• 一開始的時候我們要隨機選擇多少點(n)
• 以及要重復迭代多少次(k)

假設每個點是真正內群的概率為 w：

\[w = 內群的數目/（內群數目+外群數目） \]

通常我們不知道 w 是多少, \(w^n\)是所選擇的n個點都是內群的機率, \(1-w^n\)是所選擇的n個點至少有一個不是內群的機率, \((1 − w^n)^k\)是表示重復 k 次都沒有全部的n個點都是內群的機率, 假設算法跑 k 次以后成功的機率是p，那么

\[1 − p = (1 − w^n)^k\\ p = 1 − (1 − w^n)^k \]

我們可以通過P反算得到抽取次數K：\(K=log(1-P)/log(1-w^n)。\)

所以如果希望成功機率高：

• 當n不變時，k越大，則p越大；
• 當w不變時，n越大，所需的k就越大。
• 通常w未知，所以n選小一點比較好。

實際設置經驗

可以先設置好P和n，然后K設置為一個很大的數，然后在迭代的過程中，每次迭代完都可以計算一次w，然后對K進行動態修改。 opencv的findHomography函數中ransac就是采用的這一策略：

H = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5)  # 采用Ransac方法計算單應性矩陣

參考：https://blog.csdn.net/laobai1015/article/details/51683076/

3. Ransac應用

在計算機視覺方面，Ransac最主要的應用主要包括直線擬合，平面擬合，單應性矩陣擬合等。這里簡單介紹下直線擬合和透視矩陣擬合。

3.1 Ransac直線擬合

假設我們知道兩個變量X與Y之間呈線性關系，Y=aX+b，我們想確定參數a與b的具體值。通過實驗， 可以得到一組X與Y的測試值。雖然理論上兩個未知數的方程只需要兩組值即可確認，但由於系統誤差的原因，任意取兩點算出的a與b的值都不盡相同。通常情況下，我們可以采用最小二乘法擬合出直線方程（最小二乘法：通過計算最小均方差關於參數a、b的偏導數為零時的值）。但是最小二乘法只適合於誤差較小的情況，如果測量數據中外點較多，誤差很大，就需要采用Ransac算法。

在下圖中分別是Ransac和最小二乘法擬合的直線，可以看出兩者的差別。直接采用最小二乘法擬合直線，直線會受離群點影響，偏離中間的內群數據點，而ransac只是挑選部分數據點進行直線擬合，所以不會受離群點影響。

下面是一段示例代碼，比較了ransac和最小二乘法擬合直線的區別：

#include<iostream>
#include<opencv2/opencv.hpp>
#include<vector>
#include<set>
using namespace std;

bool ransacFitLine(vector<cv::Point2d>& points, cv::Vec4f& bestParam, vector<cv::Point2d>& inlinerPoints, int initInlinerNums, int iter, double thresh, int EndInlinerNums){
/*
    @param points:  所有的樣本點.
    @param bestParam:  擬合最好的直線方程
    @param inlinerPoints:  擬合最好時的內點集合
    @param initInlinerNums： 用來擬合直線的初始化內點個數
    @param iter： 最大迭代次數
    @param thresh： 閾值，用來判斷某個點是否在擬合的直線上
    @param EndInlinerNums： 閾值，內點個數超過閾值時停止迭代
 */

    set<int> randomIndex;
    vector<cv::Point2d> randomPoints;
    vector<cv::Point2d> tempInliners;
    int n = points.size();
    double minErr=DBL_MAX;
    bool flag = false;

    while(iter>=0){

        randomIndex.clear();
        randomPoints.clear();
        tempInliners.clear();

        // 隨機選取initInlinerNums個點，擬合直線，其他點作為測試點
        for(int i=0; i<initInlinerNums; i++){
            int index = rand()%n;
            randomIndex.insert(index);
            randomPoints.push_back(points[index]);
        }

        // 采用隨機點擬合直線
        cv::Vec4f lineParam;
        cv::fitLine(randomPoints, lineParam, cv::DIST_L2, 0, 0.01, 0.01);
        double linek = lineParam[1]/lineParam[0];
        
        // 計算測試點是否也在擬合直線上, 在直線上的點當作inliner
        for(int i=0; i<n; i++){
            if(randomIndex.find(i)==randomIndex.end()){
                double err = linek*(points[i].x - lineParam[2]) + lineParam[3] - points[i].y;
                err = err*err; // 差值的平方作為誤差
                if (err < thresh){
                    tempInliners.push_back(points[i]);
                }
            }
        }
        
        // 當前inliner內點個數超過閾值時
        if(tempInliners.size()>EndInlinerNums){
            //采用所有的內點，擬合直線並計算平均誤差
            tempInliners.insert(tempInliners.end(), randomPoints.begin(), randomPoints.end());
            cv::fitLine(randomPoints, lineParam, cv::DIST_L2, 0, 0.01, 0.01);
            double linek = lineParam[1]/lineParam[0];
            double sumErr = 0;
            for(int i=0; i<tempInliners.size(); i++){
                double err = linek*(points[i].x - lineParam[2]) + lineParam[3] - points[i].y;
                sumErr += err*err; // 差值的平方作為誤差
            }
            sumErr = sumErr/tempInliners.size();  // 誤差平方和的平均值
            if(minErr>sumErr){  // 記錄最小的平均誤差，
                minErr = sumErr;
                bestParam = lineParam;
                inlinerPoints.assign(tempInliners.begin(), tempInliners.end()); 
                flag = true;
            }
        }
        iter--;
    }
    if(!flag) cout << "did't meet fit acceptance criteria" << endl;

    return flag;
}

int main(int argc, char* argv[]){
    cv::RNG rng(200);//seed為200

    vector<cv::Point2d> points;

    // 生成500個隨機樣本點
    int nSamples = 500;   // 生成樣本數據個數
    double k = 1 + rng.gaussian(1);  // 0-2之間的隨機斜率
    cout << k << endl;
    for(int i=0; i<nSamples; i++){
        double x = rng.uniform(40, 600);
        double y = k*x + rng.uniform(-30, 30);  // 加入隨機噪聲
        x += rng.uniform(-30, 30);             // 加入隨機噪聲
        points.push_back(cv::Point2d(x, y));
    }

    // 加入outlier
    for(int j=0; j< 100; j++){
        double x = rng.uniform(200, 600);
        double y = rng.uniform(10, 300);
        points.push_back(cv::Point2d(x, y));
    }


    cv::Mat img = cv::Mat(720, 720, CV_8UC3, cv::Scalar(255,255,255));

    //繪制point
    for(int i=0; i<points.size(); i++){
        cv::circle(img, points[i], 3, cv::Scalar(0, 0, 0), 2);
    }  
    
     // 最小二乘法
    if(1){
        cv::Vec4f lineParam;
        fitLine(points, lineParam, cv::DIST_L2, 0.0, 0.01, 0.01);
        //獲取點斜式的點和斜率, y=k(x-x0) + y0
        cv::Point point0;
        point0.x = lineParam[2];
        point0.y = lineParam[3];
        double linek = lineParam[1] / lineParam[0];
        

        //繪制最小二乘法擬合的直線
        cv::Point2d point1, point2;
        point1.x=10;
        point1.y = linek*(point1.x-point0.x) + point0.y;
        point2.x = 600;
        point2.y = linek*(point2.x-point0.x) + point0.y;
        cv::line(img, point1, point2, cv::Scalar(0, 255, 0), 4);
    }    

    // ransac擬合直線
    if(1){   // 為了避免命名沖突，所以放在if局部范圍內
        cv::Vec4f lineParam;
        vector<cv::Point2d> inlinerPoints;
        bool ret = ransacFitLine(points, lineParam, inlinerPoints, 50, 1000, 7e3, 300);
        if(ret){
            cout << ret << endl;
            cout << lineParam << endl;
            cout << inlinerPoints.size() << endl;
            
            //獲取點斜式的點和斜率, y=k(x-x0) + y0
            cv::Point point0;
            point0.x = lineParam[2];
            point0.y = lineParam[3];
            double linek = lineParam[1] / lineParam[0];
            

            //繪制擬合的直線
            cv::Point2d point1, point2;
            point1.x=10;
            point1.y = linek*(point1.x-point0.x) + point0.y;
            point2.x = 600;
            point2.y = linek*(point2.x-point0.x) + point0.y;
            cv::line(img, point1, point2, cv::Scalar(0, 0, 255), 4);
        }else{
            cout << "Fitting failed, try to change parameters for ransac!!!" << endl;
        }
        
    }

    cv::imshow("img", img);
    cv::waitKey(0);
    cv::destroyAllWindows();
    return 0;
}

代碼結果如下：

3.2 Ransac單應性矩陣擬合

關於單應矩陣理解：https://zhuanlan.zhihu.com/p/49435367

在進行圖像匹配，全景拼接等時，常會用到單應性矩陣。例如在做全景拼接時，對於同一個場景，相機在不同角度拍攝了兩張照片，一般先尋找兩幅圖片的匹配特征點，然后通過匹配特征點的對應關系計算出一個矩陣，這個矩陣就是單應性矩陣，利用這個矩陣就能將兩張圖片組合在一起。所以單應矩陣描述的是兩組坐標點的關系，即：

\[\begin{bmatrix} x_0 \\ y_0 \\ 1 \end{bmatrix} =H\begin{bmatrix} X'_0 \\ Y'_0 \\ 1 \end{bmatrix} \\ \]

上面H即為單應性矩陣，\((x_0, y_0)\)為第一幅圖片上的坐標點，\((X'_0, Y'_0)\)為第二幅圖片上的坐標點，兩個點是一組匹配點。H一般為3x3的矩陣，形式如下：

\[H = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{bmatrix} \]

這個矩陣有9個未知數(自由度為8)，為了求解矩陣我們一般令\(a_{33}=1\)， 則一組匹配點可以列如下兩個方程：

若有四組匹配點，就能列出8個方程，便能計算出單應矩陣H， 如下圖所示：

opencv的getPerspectiveTransform函數就是通過四組匹配點來計算單應性矩陣(投影矩陣/透視矩陣)，如下：

　H = cv2.getPerspectiveTransform(src, dst)   # 返回3*3的矩陣
 參數：    
      src：原圖像中的四組坐標，如 np.float32([[56,65],[368,52],[28,387],[389,390]])
      dst: 轉換后的對應四組坐標，如np.float32([[0,0],[300,0],[0,300],[300,300]])

如果我們只有四組匹配點，而且能夠保證匹配點都是正確的，的確可以采用上述方法求解單應矩陣，但是如果我們有100組匹配點，而且有些匹配點可能是錯誤的，該怎么求解一個擬合最好的矩陣呢?

仔細觀察上面八個方程的形式，發現是一個矩陣求解的問題：

\[AH=B \]

其中，A是8x8的矩陣，H是8x1的矩陣，B也是8x1的矩陣。若有100組匹配點，則A是200x8的矩陣，B是200x1的矩陣。

聯想到上面最小二乘法擬合直線的問題：\(y=wx\)，那么在這里x就是矩陣A的每一行(8x1的向量)，\(w\)就是矩陣H，\(y\)

就是矩陣B的每一行(1x1), 如下圖所示：

所以100組匹配點就是對應200條數據，采用最小二乘法進行擬合(不算是直線了，算平面擬合了)。

有了求解方法，另一個問題就是：100組匹配點中可能存在錯誤的匹配，很明顯錯誤的匹配就是所謂的outlier了。很自然的能想到：采用ransac代替上面的最小二乘法即可。

上述就是采用Ransac擬合單應性矩陣的思想，在opencv中findHomography函數實現了上述過程，如下所示：

#findHomography(srcPoints,dstPoints,method=0,ransacReprojThreshold=3,mask=noArray(),maxIters=2000,confidence=0.995)	

H = cv2.findHomography(src_pts, dst_pts, method=cv2.RANSAC, ransacReprojThreshold=5)  # 生成變換矩陣
參數：
	src_pts：原圖像中的坐標(大於等於四組)
	dst_pts: 轉換后的對應坐標(大於等於四組)
	method：求解單應矩陣的方法，取值如下：
			0：最小二乘法
			RANSAC：RANSAC擬合方法
			LMEDS： Least-Median robust method
			RHO：PROSAC-based robust method
   ransacReprojThreshold：ransac中判斷是否是內點的閾值

下面看一個全景拼接的示例，全景拼接的常規流程 如下：

• 1、針對某個場景拍攝多張/序列圖像

• 2、通過匹配特征（sift匹配）計算下一張圖像與上一張圖像之間的匹配特征點

• 3、對匹配特征點進行篩選，剔除掉部分匹配點。采用ransanc擬合匹配點之間的單應性矩陣

• 3、利用單應性矩陣進行圖像映射，將下一張圖像疊加到上一張圖像的坐標系中

• 4、變換后的融合/合成
  

下面是一個代碼示例：

# coding: utf-8
import numpy as np
import cv2

left_img = cv2.imdecode(np.fromfile(r"D:\筆記\cpp\image_processing\圖像拼接\left.jpg", dtype=np.uint8), 1)
left_img = cv2.resize(left_img, (600, 400))
right_img = cv2.imdecode(np.fromfile(r"D:\筆記\cpp\image_processing\圖像拼接\right.jpg", dtype=np.uint8), 1)
right_img = cv2.resize(right_img, (600, 400))
left_gray = cv2.cvtColor(left_img, cv2.COLOR_BGR2GRAY)
right_gray = cv2.cvtColor(right_img, cv2.COLOR_BGR2GRAY)

hessian = 300
surf = cv2.xfeatures2d.SIFT_create(hessian) # 將Hessian Threshold設置為400,閾值越大能檢測的特征就越少
kp1, des1 = surf.detectAndCompute(left_gray, None)  # 查找關鍵點和描述符
kp2, des2 = surf.detectAndCompute(right_gray, None)

# kp1s = np.float32([kp.pt for kp in kp1])
# kp2s = np.float32([kp.pt for kp in kp2])

# 顯示特征點
img_with_drawKeyPoint_left = cv2.drawKeypoints(left_gray, kp1, np.array([]), flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv2.imshow("img_with_drawKeyPoint_left", img_with_drawKeyPoint_left)

img_with_drawKeyPoint_right = cv2.drawKeypoints(right_gray, kp2, np.array([]), flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
cv2.imshow("img_with_drawKeyPoint_right", img_with_drawKeyPoint_right)



'''
BFMatcher簡稱暴力匹配，意思就是嘗試所有可能匹配，實現最佳匹配。

FlannBasedMatcher簡稱最近鄰近似匹配。
是一種近似匹配方法，並不追求完美！，因此速度更快。
可以調整FlannBasedMatcher參數改變匹配精度或改變算法速度。
參考：https://blog.csdn.net/claroja/article/details/83411108
'''
FLANN_INDEX_KDTREE = 0  # 建立FLANN匹配器的參數

indexParams = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)  # 配置索引，密度樹的數量為5
searchParams = dict(checks=50)  # 指定遞歸次數
# FlannBasedMatcher：是目前最快的特征匹配算法（最近鄰搜索）
flann = cv2.FlannBasedMatcher(indexParams, searchParams)  # 建立匹配器

# 參考https://docs.opencv.org/master/db/d39/classcv_1_1DescriptorMatcher.html#a378f35c9b1a5dfa4022839a45cdf0e89
'''
int queryIdx –>是測試圖像的特征點描述符（descriptor）的下標，同時也是描述符對應特征點（keypoint)的下標。

int trainIdx –> 是樣本圖像的特征點描述符的下標，同樣也是相應的特征點的下標。

int imgIdx –>當樣本是多張圖像的話有用。

float distance –>代表這一對匹配的特征點描述符（本質是向量）的歐氏距離，數值越小也就說明兩個特征點越相像。
'''

matches = flann.knnMatch(des1, des2, k=2)  # 得出匹配的關鍵點


good = []
# 提取優秀的特征點
for m, n in matches:
    if m.distance < 0.7 * n.distance:  # 如果第一個鄰近距離比第二個鄰近距離的0.7倍小，則保留
        good.append(m)

src_pts = np.array([kp1[m.queryIdx].pt for m in good])  # 查詢圖像的特征描述子索引
dst_pts = np.array([kp2[m.trainIdx].pt for m in good])  # 訓練(模板)圖像的特征描述子索引


# findHomography參考https://docs.opencv.org/master/d9/d0c/group__calib3d.html#ga4abc2ece9fab9398f2e560d53c8c9780
# 單應矩陣:https://www.cnblogs.com/wangguchangqing/p/8287585.html
H = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5)  # 生成變換矩陣

h1, w1 = left_gray.shape[:2]
h2, w2 = right_gray.shape[:2]
shift = np.array([[1.0, 0, w1], [0, 1.0, 0], [0, 0, 1.0]])
# 點積
M = np.dot(shift, H[0])  # 獲取左邊圖像到右邊圖像的投影映射關系, 向右移動w1?

dst = cv2.warpPerspective(left_img, M, (w1+w2, max(h1, h2)))  # 透視變換，新圖像可容納完整的兩幅圖
cv2.imshow('left_img', dst)  # 顯示，第一幅圖已在標准位置
dst[0:h2, w1:w1+w2] = right_img  # 將第二幅圖放在右側
# cv2.imwrite('tiled.jpg',dst_corners)
cv2.imshow('total_img', dst)
cv2.imshow('leftgray', left_img)
cv2.imshow('rightgray', right_img)
cv2.waitKey(0)
cv2.destroyAllWindows()

效果如下：

參考：

https://blog.csdn.net/hzwwpgmwy/article/details/83578694

https://zhuanlan.zhihu.com/p/82793501

4 Ransac優缺點

優點：

它能魯棒的估計模型參數。例如，它能從包含大量局外點的數據集中估計出高精度的參數。

缺點：

1. 它計算參數的迭代次數沒有上限；如果設置迭代次數的上限，得到的結果可能不是最優的結果，甚至可能得到錯誤的結果。
2. RANSAC只有一定的概率得到可信的模型，概率與迭代次數成正比。
3. 它要求設置跟問題相關的閥值。
4. RANSAC只能從特定的數據集中估計出一個模型，如果存在兩個（或多個）模型，RANSAC不能找到別的模型。
5. 要求數學模型已知

參考：https://zhuanlan.zhihu.com/p/36301702/
https://blog.csdn.net/leonardohaig/article/details/104570965