OpenCV實驗（8）：基于區域的視覺立體匹配

一. 實驗要求

二. 實驗思路

1. 想到啥說啥

實驗要求參考文獻《Obtaining Depth Maps From Color Images By Region Based Stereo Matching Algorithms》，而且實驗給的圖也是這個論文中的，所以該實驗可以看做是對論文的復現，對于這最后一個實驗，如果有更好的想法或者發現的問題可以評論或私信我，
繼續參考學長的博客https://blog.csdn.net/qq_41748260/article/details/103992462
論文中立體匹配演算法的思想直接參考學長的博客就可以了，這里就不重復了，

2. 具體程序

根據"a) Global Error Energy Minimization by Smoothing Functions"，中的steps1、2、3就實作了視差圖的計算，再根據"Filtering Unreliable Disparity Estimation By Average Error Thresholding Mechanism"計算具有可靠差異的視差圖，接著根據"Depth Map Generation From Disparity Map"，從視差圖生成深度圖即可，

發現的一些問題和更具體的已經在代碼注釋中說明了，

三. 實驗效果

四. 實驗代碼

下面這個是完整的代碼，代碼注釋中給出了完整詳細的解釋

import numpy as np
import cv2 as cv
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
from mpl_toolkits.mplot3d import Axes3D


def disparity_GEEMBSF(imgLeft, imgRight, windowSize=(3, 3), dMax=30, alpha=1):
    """
    論文中的
    "a) Global Error Energy Minimization by Smoothing Functions"
    函式名的命名方式：disparity 表示視差，后面的一串為該方法標題縮寫

    :param imgLeft: 左圖
    :param imgRight: 右圖
    :param windowSize: (n, m)，視窗的尺寸為 n x m
    :param dMax: 中值濾波迭代次數
    :param alpha: 閾值系數
    :return: 視差圖，平均誤差能量矩陣，函式運行時間
    """
    timeBegin = cv.getTickCount()  # 記錄開始時間

    n, m = windowSize  # 視窗大小
    rows, cols, channels = imgLeft.shape  # 該實驗中 imgLeft 和 imgRight 的 shape 是一樣的，rows=185,cols=231,channels=3
    # 觀察到論文中和實驗要求中所給的左右原始圖是185(行)x231(列)像素的，而結果圖中大約是185(行)x190(列)的，
    # 我們的結果中imgDisparity會看到右邊緣有大量黑色的區域，如果將其去掉，那么會更美觀也會與論文中的結果更加符合
    # 因此在前面我們將cols=190
    cols = 190
    errorEnergyMatrixD = np.zeros((rows, cols, dMax), dtype=np.float)  # 誤差能量矩陣（共dMax層），方便后續計算
    errorEnergyMatrixAvgD = np.zeros((rows, cols, dMax), dtype=np.float)  # 平均誤差能量矩陣（共dMax層），方便后續計算
    imgDisparity = np.zeros((rows, cols), dtype=np.uint8)  # 具有可靠差異的視差圖，將作為結果回傳

    # 計算誤差能量矩陣 errorEnergyMatrix，平均誤差能量矩陣 errorEnergyMatrixAvg
    # 發現的問題：公式(1)在求和計算時，x從i到i+n，y從j到j+m，求和符號是包括上界和下界的，
    # 那么視窗大小就變為了(n+1, m+1)，這與論文中提到的視窗大小為(n, m)是不符的，這一點使我疑惑，
    # 我在處理的時候，計算的是x從i到i+n-1，y從j到j+m-1，

    # 先padding，這樣方便使用numpy加速計算
    imgLeftPlus = cv.copyMakeBorder(imgLeft, 0, 0, n-1, m-1+dMax, borderType=cv.BORDER_REPLICATE)
    imgRightPlus = cv.copyMakeBorder(imgRight, 0, 0, n-1, m-1, borderType=cv.BORDER_REPLICATE)
    # 迭代 dMax 次
    for d in range(dMax):
        # 對整個影像進行遍歷
        for i in range(rows):
            for j in range(cols):
                # 對于每個 (i, j, d) 根據公式(1)計算誤差能量矩陣
                errorEnergy = (imgLeftPlus[i:i+n, j+d:j+m+d, ...] - imgRightPlus[i:i+n, j:j+m, ...]) ** 2
                errorEnergyMatrixD[i, j, d] = np.sum(errorEnergy) / (3 * n * m)
        # 對 errorEnergyMatrix 進行遍歷
        for i in range(rows):
            for j in range(cols):
                # 對于每個 (i, j, d) 根據公式(2)計算平均誤差能量矩陣
                errorEnergyMatrixAvgD[i, j, d] = np.sum(errorEnergyMatrixD[i:i+n, j:j+m, d]) / (n * m)
        # 論文中說到了（公式(1)下方）
        # For a predetermined disparity
        # search range (w), every e(i, j, d) matrix respect to disparity is smoothed by applying
        # averaging filter many times. (See Figure 1.b)
        # 也就是對于每個e(i, j, d)，進行多次平均濾波，在這里我選擇執行3次，
        for k in range(3):
            for i in range(rows):
                for j in range(cols):
                    # 對于每個 (i, j, d) 根據公式(2)計算平均誤差能量矩陣
                    # 下面i+n和j+m越了界也是沒有問題的，切片會正常計算
                    errorEnergyMatrixAvgD[i, j, d] = np.sum(errorEnergyMatrixAvgD[i:i + n, j:j + m, d]) / (n * m)

    errorEnergyMatrixAvg = np.min(errorEnergyMatrixAvgD, axis=2)  # 平均誤差能量矩陣（最終的，只有一層）
    imgDisparity[:, :] = np.argmin(errorEnergyMatrixAvgD, axis=2)  # 視差圖
    imgOrignal = imgDisparity.copy()  # 保留一份，并作為結果回傳
    # cv.imwrite('../images/out/disparity.png', imgDisparity)

    # cv.imshow('imgDisparity1', imgDisparity)
    # imgDisparityHist = cv.equalizeHist(imgDisparity)  # imgDisparity 太黑了，直方圖均衡化方便觀察
    # cv.imshow('imgDisparity and imgDisparityHist111', np.hstack([imgDisparity, imgDisparityHist]))  # 水平排列兩幅影像進行顯示
    # cv.imshow('imgDisparityHistNon', imgDisparityHist)
    # cv.imwrite('../images/out/imgDisparityHist.png', imgDisparityHist)

    # 下面的部分我們實作論文中的：（包含公式 5、6、7、8、9）
    # 可靠差異的視差圖
    # "Filtering Unreliable Disparity Estimation By Average Error Thresholding Mechanism"
    Ve = alpha * np.mean(imgDisparity)  # 計算Ve
    temp = errorEnergyMatrixAvg.copy()
    temp[temp > Ve] = 0
    temp[temp != 0] = 1
    temp = temp.astype(np.int)
    imgDisparity = np.multiply(imgDisparity, temp).astype(np.uint8)  # 大于Ve的設定為0
    # Sd = np.sum(temp)  # 計算Sd
    # Rd = 1 / (np.sum(np.multiply(imgDisparity, temp).astype(np.float)) * Sd)  # 計算Rd

    timeEnd = cv.getTickCount()  # 記錄結束時間
    time = (timeEnd - timeBegin) / cv.getTickFrequency()  # 計算總時間

    return imgOrignal, imgDisparity, errorEnergyMatrixAvg, time


def depthGeneration(imgDisparity, f=30, T=20):
    """
    實作論文中的"Depth Map Generation From Disparity Map"
    根據視差圖，實作深度圖

    :param imgDisparity: 具有可靠差異的視差圖
    :param f: 焦距
    :param T: 間距
    :return: 深度圖
    """
    # 實作公式(4)
    rows, cols = imgDisparity.shape
    imgDepth = np.zeros((rows, cols), dtype=np.uint8)
    for i in range(rows):
        for j in range(cols):
            if imgDisparity[i, j] == 0:
                imgDepth[i, j] = 0
            else:
                imgDepth[i, j] = f * T // imgDisparity[i, j]
    # imgDepthHist = cv.equalizeHist(imgDepth)
    # cv.imshow('imgDepthHist', imgDepthHist)

    return imgDepth


imgLeft = cv.imread('../images/images5/view1m.png')
imgRight = cv.imread('../images/images5/view5m.png')
imgOrignal, imgDisparity, errorEnergyMatrixAvg, time = disparity_GEEMBSF(imgLeft, imgRight)
print(time)  # 列印函式運行時間
imgDepth = depthGeneration(imgDisparity)

# imgDisparityHist = cv.equalizeHist(imgDisparity)  # imgDisparity 太黑了，直方圖均衡化方便觀察
# cv.imshow('imgDisparityHistReal', imgDisparityHist)

# 下面繪制三幅影像
# 視差圖影像
plt.figure()
plt.imshow(imgOrignal, cmap='gray', vmin=0, vmax=40)
plt.title('disparity figure')  # 設定標題
ax = plt.gca()  # 回傳坐標軸實體
x_major_locator = MultipleLocator(20)  # 刻度間隔為20
y_major_locator = MultipleLocator(20)
ax.xaxis.set_major_locator(x_major_locator)  # 設定坐標間隔
ax.yaxis.set_major_locator(y_major_locator)
plt.colorbar()  # 設定colorBar
plt.show()

# 深度圖影像
plt.figure()
plt.imshow(imgDepth, cmap='gray', vmin=0, vmax=120)
plt.title('depth figure (cm)')  # 設定標題
ax = plt.gca()  # 回傳坐標軸實體
x_major_locator = MultipleLocator(20)  # 刻度間隔為20
y_major_locator = MultipleLocator(20)
ax.xaxis.set_major_locator(x_major_locator)  # 設定坐標間隔
ax.yaxis.set_major_locator(y_major_locator)
plt.colorbar()  # 設定colorBar
plt.show()

# 3D影像
# 這里又使我疑惑，實驗要求中說是“將深度圖用3D視圖展示”，但實驗結果中3D視圖的z軸是0~40也就是視差范圍
# 在實驗結果中的深度圖顯示的是cm，所以深度圖的3D視圖的縱坐標也應該是厘米，而不是視差
# 但為了能和實驗結果的圖匹配，我將視差圖進行了3D展示
fig = plt.figure()
ax = Axes3D(fig)
rows, cols = imgDisparity.shape
x = np.arange(cols)
y = np.arange(rows)
x, y = np.meshgrid(x, y)
z = imgDisparity
surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap=plt.get_cmap('gray'),
                       vmin=0, vmax=40)
ax.set_title('3D figure')
ax.set_xlabel('x')
ax.set_ylabel('y')

ax.set_zticks(np.arange(0, 41, 10))
ax.set_xticks(np.arange(0, 185, 30))
ax.set_yticks(np.arange(0, 190, 30))
plt.savefig('3D figure.png', bbox_inches='tight')
plt.show()

cv.waitKey(0)
cv.destroyAllWindows()