立體匹配-----NCC視差匹配

目錄

一、立體匹配算法

　　　1.立體匹配算法分類

二、NCC 視差匹配方法

           1.原理

 　　　2.NCC計算公式

　　　3.算法流程

　　　4.代碼實現

　　    5.不同場景運行

三、結論

四、遇到的問題及解決方法

一、立體匹配算法

1.立體匹配算法分類

 在立體匹配中，匹配問題可以看成是尋找兩組數據相關程度的過程。根據采用圖像表示的基元不同，立體匹配算法有多種分類。

  ①根據算法運行時約束的作用范圍：分為局部匹配算法和全局匹配算法。

  ②基於生成的視差圖：可分為稠密匹配和稀疏匹配。稠密匹配：是基於生成的視差圖，對於所有像素都能生成確定視差值，稱為稠密匹配。稀疏匹配：只選擇關鍵像素點[通常為角點或者邊緣點]計算視差值的方法稱為稀疏匹配，該算法計算速度較快，但后續還需要通過插值算法計算缺失像素點的視差值，因此應用場景上有很大限制。

1.1全局匹配算法

  全局（半全局）立體匹配算法主要是采用了全局的優化理論方法估計視差，建立一個全局能量函數，其包含一個數據項和平滑項，通過最小化全局能量函數得到最優的視差值。其中，圖割、置信傳播、動態規划,粒子群算法、遺傳算法等優化算法都是常用的求解能量最小化的方法。（以后作補充，也可參考我之前的粒子群算法、遺傳算法文章）

1.2局部匹配算法

（1）算法介紹：局部立體匹配算法又稱基於窗口的方法或基於支持區域的方法。算法對參考圖像中的每個像素計算一個合適大小、形狀和權重的窗口，然后對這個窗口內的視差值進行加權平均。理想的支持窗口應該完全覆蓋弱紋理區域，並在窗口內深度連續。與全局立體匹配算法相似，通過優化一個代價函數的方法計算最佳視差。但是，在局部立體匹配算法的能量函數中，只有基於局部區域的約束數據項，沒有平滑項。局部匹配算法僅利用某一點鄰域的灰度、顏色、梯度等信息進行計算匹配代價，計算復雜度較低，大多實時的立體匹配算法都屬於局部立體匹配的范疇，但局部立體匹配算法對低紋理區域、重復紋理區域、視差不連續和遮擋區域匹配效果不理想。

（2）基本原理：在參考圖像中選擇一個點，選擇該點鄰域內一個支持窗口，然后依據一定的相似性判斷准則，在待匹配圖像中尋找與支持窗口最相似的子窗口，該子窗口所對應的像素點即為對應的匹配點。

（3）不同場景運用：固定窗口代價聚合使用固定大小和形狀的窗口作為代價聚合的基元，通常是一個矩形，並假設支持窗口內的其它像素點與待匹配點具有相同的視差。固定窗口法精度不高，但易實現、耗時短，在一些對實時性要求極高的場合得到了應用。

  基於雙邊濾波的代價聚合算法仍然使用固定大小和形狀的窗口，但窗口內的元素權重不同，權重由目標圖像在該窗口內像素與窗口中心的灰度差和距離計算。基於雙邊濾波的代價聚合算法精度高，但計算復雜，實時性差，算法性能隨窗口尺寸指數增加。

  基於分割的代價聚合算法的主要思想是：預先將作為參考圖像的左圖進行分割，對於支持窗口內與窗口中心處於同一分割的像素，對應的權值取1，否則為一個遠小於1的正數。但是圖像分割是一個非常耗時的操作，同樣無法在實時性要求較高的場合使用。

  基於十字的代價聚合算法（Cross-based Cost Aggregation,CBCA）的支持窗口形狀並不確定，會根據匹配點鄰域的灰度值而改變，具體實現后面的更新中會介紹。該方法可以使用GPU並行計算，具有較好的實時性，現廣泛應用於各種算法的代價聚合步驟。

（4）優點：算法成熟、計算簡單、速度快，能進行圖像實時處理，匹配精度較高，計算量比全局匹配算法更小。

（5）缺點：對噪聲敏感，弱紋理區易造成誤匹配

二、NCC 視差匹配方法

1.原理：對於原始的圖像內任意一個像素 點(px, py )構建一個n x n的鄰域作為匹配窗口。然后對於目標相素位置(px + d, py )同樣構建一個n x n大小的匹配窗口，對兩個窗口進行相似度度量,注意這里的有一個取值范圍。對於兩幅圖像來說，在進行NCC計算之前要對圖像處理，也就是將兩幀圖像校正到水平位置，即光心處於同一 水平線上，此時極線是水平的,否則匹配過程只能在傾斜的極線方向上完成,這將消耗更多的計算資源。

2.NCC計算公式：

 

 

 

求出的NCC(p,d)值在[−1,1]之間

3.算法流程：

 ①采集圖像：通過標定好的雙目相機采集圖像，當然也可以用兩個單目相機來組合成雙目相機。 

 ②極線校正：校正的目的是使兩幀圖像極線處於水平方向，或者說是使兩幀圖像的光心處於同一水平線上。通過校正極線可以方便后續的NCC操作。 

 ③特征匹配：這里便是我們利用NCC做匹配的步驟，匹配方法如上所述，右視圖中與左視圖待測像素同一水平線上相關性最高的即為最優匹配。完成匹配后，需要記錄其視差d，即待測像素水平方向xl與匹配像素水平方向xr之間的差值d=xr−xl，最終可以得到一個與原始圖像尺寸相同的視差圖D。

 ④深度恢復：通過上述匹配結果得到的視差圖D，我們可以很簡單的利用相似三角形反推出以左視圖為參考系的深度圖。計算原理如下圖所示：

 

 

 

如圖，Tx為雙目相機基線，f為相機焦距，這些可以通過相機標定步驟得到。而xr - xl就是視差D。

通過公式 z = f * Tx / d可以很簡單地得到以左視圖為參考系的深度圖了。

4.代碼實現：

（1）源代碼

# -*- coding: utf-8 -*-
from PIL import Image
from pylab import *
import cv2
from numpy import *
from numpy.ma import array
from scipy.ndimage import filters
def plane_sweep_ncc(im_l,im_r,start,steps,wid):
    """ 使用歸一化的互相關計算視差圖像 """
    m,n = im_l.shape
    # 保存不同求和值的數組
    mean_l = zeros((m,n))
    mean_r = zeros((m,n))
    s = zeros((m,n))
    s_l = zeros((m,n))
    s_r = zeros((m,n))
    # 保存深度平面的數組
    dmaps = zeros((m,n,steps))
    # 計算圖像塊的平均值
    filters.uniform_filter(im_l,wid,mean_l)
    filters.uniform_filter(im_r,wid,mean_r)
    # 歸一化圖像
    norm_l = im_l - mean_l
    norm_r = im_r - mean_r
    # 嘗試不同的視差
    for displ in range(steps):
        # 將左邊圖像移動到右邊，計算加和
        filters.uniform_filter(np.roll(norm_l, -displ - start) * norm_r, wid, s) # 和歸一化
        filters.uniform_filter(np.roll(norm_l, -displ - start) * np.roll(norm_l, -displ - start), wid, s_l)
        filters.uniform_filter(norm_r*norm_r,wid,s_r) # 和反歸一化
        # 保存 ncc 的分數
        dmaps[:,:,displ] = s / sqrt(s_l * s_r)
        # 為每個像素選取最佳深度
    return np.argmax(dmaps, axis=2)

def plane_sweep_gauss(im_l,im_r,start,steps,wid):
 """ 使用帶有高斯加權周邊的歸一化互相關計算視差圖像 """
 m,n = im_l.shape
 # 保存不同加和的數組
 mean_l = zeros((m,n))
 mean_r = zeros((m,n))
 s = zeros((m,n))
 s_l = zeros((m,n))
 s_r = zeros((m,n))
 # 保存深度平面的數組
 dmaps = zeros((m,n,steps))
 # 計算平均值
 filters.gaussian_filter(im_l,wid,0,mean_l)
 filters.gaussian_filter(im_r,wid,0,mean_r)
 # 歸一化圖像
 norm_l = im_l - mean_l
 norm_r = im_r - mean_r
 # 不同的視差
 for displ in range(steps):
     # 將左邊圖像移動到右邊，計算加和
     filters.gaussian_filter(np.roll(norm_l, -displ - start) * norm_r, wid, 0, s) # 和歸一化
     filters.gaussian_filter(np.roll(norm_l, -displ - start) * np.roll(norm_l, -displ - start), wid, 0, s_l)
     filters.gaussian_filter(norm_r*norm_r,wid,0,s_r) # 和反歸一化
     # 保存 ncc 的分數
     dmaps[:,:,displ] = s / np.sqrt(s_l * s_r)
 # 為每個像素選取最佳深度
 return np.argmax(dmaps, axis=2)

im_l = array(Image.open(r'C:/Users/LE/PycharmProjects/untitled/NCC/3.jpg').convert('L'), 'f')
im_r = array(Image.open(r'C:/Users/LE/PycharmProjects/untitled/NCC/4.jpg').convert('L'),'f')
# 開始偏移，並設置步長
steps = 12
start = 4
# ncc 的寬度
wid = 9
res = plane_sweep_ncc(im_l,im_r,start,steps,wid)
import scipy.misc
scipy.misc.imsave('depth.png',res)
show()