【論文筆記】PatchMatch Stereo - Stereo Matching with Slanted Support Windows

—— patch match應用到立體匹配的一篇論文，晦澀難懂，硬是硬着頭皮看了幾遍

NOTE

1. 隨機搜索的最優平面和視差的關系？

2. “estimate an individual 3D plane at each pixel”是對每一個窗口的中心點都找一個最優平面？

3. 隨機值初始化的理論基礎？

4. 在全局法中的作用

5. 實際實現理解

answer：

1. 當空間中的點x, y均確定，法向量也確定之后，最優平面和最優視差有着一一對應的關系，注意這個平面是定義在(x, y, 視差)空間下的平面
2. 對，每一個點都會有一個平面，在搜索和傳播的時候通過鄰域視差相似性以及左右視差一致性來進行傳播，而是否傳播的條件就是對比傳播前后該點的代價值是否減小。這時候就涉及到了代價聚合的思想，在搜索傳播的過程中已經完成代價聚合。
3. TODO
4. 見With Global Methods部分
5. 見第6部分：代碼理解

前言

局部法窗口內視差為常數的假設

傳統局部法求代價聚合的方法都是假設一個窗口內的視差為常數，在包含邊緣的窗口以及斜面的窗口，此情況便不成立了；

對於包含邊緣的窗口這類場景，已經有BFCA, GFCA等保邊平滑濾波應用到代價聚合來處理這個問題，避免代價值跨邊緣聚合；

亞像素精度的實現

通常有兩種方法：

1. disparity refine中用各種fitting來實現亞像素增強
2. 在匹配過程中就將此問題考慮進去（這時候，看上去就不得用WTA了吧）

Abstract

主要思想：

對局部窗口中的每一個點，均估計一個最優平面；

最優平面對應着最優視差，最優視差對應着最低左右匹配代價值；

最優平面的估計方法：

• spatial propagation：planes propagate among the center point and its neighbor points
• view propagation: planes propagate among the left view and the right view
• temporal propagation: planes propagate among the preceding and consecutive frames

貢獻點：

• 將Patch Match應用到立體匹配中，用於解決包含邊緣的局部窗口以及包含斜面的窗口內視差不符合常數這個假設的問題；
• 在匹配過程中，同步搜索最優平面以實現亞像素精度視差匹配
• 可擴展至全局法立體匹配中，優化其效果；

貢獻點1的體現可以下面這張圖來體現：

常規的局部法匹配，通常都是用figure 1(a)來假設窗口內的視差值，和實際視差並不符，figure 1(b)是本文的思想，對比即可知。

性能與效果：

性能就不說了，真的是有點慢，n多個全圖代價計算與視差更新循環...

效果：目前還未進行大規模圖集測試，作者說是當時Middlebury上排第2的，局部法當中效果最好的。

Future:

關於future的第一要務，肯定是加速，作者也說要嘗試在GPU上實現real-time了；

另一方面，還想嘗試應用到光流領域。

研究背景

算法

模型

對每個像素點\(p\), 尋找一個最優平面\(f_p\), 一旦\(f_p\)找到，則這個點的視差便可由下式計算得到：

\[d_p = a_{f_p}p_x + b_{f_p}p_y + c_{f_p} \tag{1} \]

最優平面的定義:

\[f_p = \arg\min_{f\in\mathcal{F}}m(p, f) \tag{2} \]

值得注意：\(\mathcal{F}\)的取值是無窮的，\(m(p, f)\)就是常規的代價聚合公式：

\[m(p, f) = \sum_{q \in W_p}{w(p, q) \cdot \rho(q, q - (a_{f_p}p_x + b_{f_p}p_y + c_{f_p}))} \tag{3} \]

其中的\(w(p, q) = e^{- \frac{|I_p - I_q|}{\gamma}}\), 就是常規的顏色相似度度量，而\(\rho(q, q - (a_{f_p}p_x + b_{f_p}p_y + c_{f_p}))\)為匹配代價SAD + gradient，其的計算方式如下：

\[\rho(q, q') = (1 - \alpha) \cdot \min(|I_q - I_q'|, \tau_{col}) + \alpha \cdot \min{(|\nabla{I_q} - \nabla{I_q'}|, \tau_{grad})}\tag{4} \]

With Patch Match

前一部分已經很清楚的敘述了整個匹配的思想，遺留問題是對於每個像素點\(p\)，尋找一個最優平面\(f_p\)，這個平面要怎么尋找呢？

作者用patch match的思想來解決，他自己也稱：

We show that an ideal algorithm to solve this problem is PatchMatch [1] that we extend to find an approximate nearest neighbor according to a plane

隨機值初始化的理論基礎——大數定理？

We find the plane for a region by initializing each pixel with a random plane. The hope is that after this random initialization at least one pixel of the region carries a plane that is close to the correct one.

基於patch match展開的一系列迭代優化：

Random Initialization

對於點\((x_0, y_0)\)，隨機賦值\(minDisp < z0 < maxDisp\)，即有\(P = (x_0, y_0, z_0)\)在平面上, 設平面法向量為 \(\vec{n} = (n_x, n_y, n_z)\), 則公式1中的系數有：

\[a_f := -\frac{n_x}{n_z}, b_f := - \frac{n_y}{n_z}; c_f := \frac{n_xx_0 + n_yy_0 + n_zz_0}{n_z} \]

上式由點法式方程結合式1可推導得到，點法式方程如下：

\[n_x \cdot (x - x_0) + n_y \cdot (y - y_0) + n_z \cdot (z - z_0) = 0 \tag{5} \]

而上式中的 \(z\) 即為傳說中的視差 \(d_p\)

NOTE:

1. 可以通過強制\(\vec{n} = (0, 0, 1)\)來將窗口設置為fronto-parallel
2. 可以通過強制視差為整數來關閉亞像素優化

Iteration

作者一共進行了3次迭代，奇數次迭代從左上角像素遍歷到右下角像素，偶數次迭代從右下角遍歷到左上角；

每次迭代需要依次對左視角和右視角執行以下4步。。。

相當於以下4步要全圖做6遍，這速度怎么可能快得起來。。。

1. Spatial Propagation

   主要思想：相鄰像素點間很可能具有相似的平面

   偽代碼： check \(m(p, f_q) < m(p, f_p)\)，如果是，則將當前點的平面更新為鄰域點\(q\)的平面\(f_q\)

   Note: 在奇數次迭代，check當前點左邊和上邊兩個點的平面，在偶數次迭代，check當前點右邊和下邊兩個點的平面

2. View Propagation

   主要思想：一個像素點和其在另一視角匹配點間很可能具有相似的平面

   We check all pixels of the second view that have our current pixel p as a matching point according to their current plane.

   偽代碼： check \(m(p, f_p') < m(p, f_p)\)，如果是，則將當前點的平面更新為另一視角匹配點\(p'\)的平面\(f_p'\)

3. Temporal Propagation

   ——主要用於視頻處理

   主要思想：在視頻連續幀間，同一位置的像素點很可能具有相似的平面

   偽代碼：check \(m(p, f_p') < m(p, f_p)\)，如果是，則將當前點的平面更新為另一幀中對應位置點\(p'\)的平面\(f_p'\)

4. Plane Refinement

   ——看上去邏輯和patch match那個搜索有關系

   設\(z_0\)在空間中變化范圍為：\([-\Delta_{z_0}^{max}, \Delta_{z_0}^{max}]\),法向量\(\vec{n}\)的變化范圍為 \([-\Delta_{n}^{max}, \Delta_{n}^{max}]\), 設\(z_0' := z_0 + \Delta_{z_0}\), \(\vec{n'} := unit(\vec{n} + \vec{\Delta_n})\), 通過由點\(P' = (x_0, y_0, z_0')\)和法向量\(\vec{n'}\)來確定的新平面\(f_p'\)；

   check \(m(p, f_p') < m(p, f_p)\)，如果是，則將當前點的平面更新為平面\(f_p'\)。

   這部分的refine是迭代的，看上去就是一個二分搜索，同時對\(z\)和法向量\(\vec{n}\)的搜索：

   剛開始\(\Delta_{z_0}^{max} = maxdisp / 2\), \(\Delta_{n}^{max}=1\), 接着以\(\Delta_{z_0}^{max} := \Delta_{z_0}^{max} / 2, \Delta_{n}^{max} := \Delta_{n}^{max} / 2\)幅度減小搜索半徑；

   作者在注釋部分提出，闊以用梯度下降替代這一步。

Post-processing

整個后處理流程比較常規：LR check出可靠區域和不可靠區域(通常是遮擋區域和誤匹配區域)，對不可靠區域用可靠區域的值進行填充，最后使用權重中值濾波來refine，主要是為了去除一些水平條紋

• LR Check
• Invalidated Pixel Filling
  • 用與此不可靠點最相近的可靠點左右視差最小值來進行填充，而非直接填充常數
• Weighted Median Filter
  • 對填充區域使用WMF

With Gloabal Methods

全局法中使用的基於局部窗口計算的自適應加權匹配代價具有獲取視差不連續區域的特性，但也繼承了"fronto-parallel"的問題；

作者指出，只要窗口半徑為1，就不存在這個問題，也就是對每個像素點計算了一個匹配代價，在補充材料里提出這種匹配代價闊以embed為能量函數的數據項，用於全局法優化；

數據項[8] + 平滑項[13] + \(\alpha\)-expansion[14]優化，思想如下：

代碼理解

// spatial propagation
/*********************************************************************
1. 設當前點對應的平面為舊平面
2. 左鄰域點對應的平面為新平面（設當前為奇數次迭代）
3. 注意p_neigb = planes[cpv](ny, nx)傳入plane_match_cost函數中計算cost值：
   實質就是將上一次(包含隨機初始化)鄰域點的視差作為當前點視差傳入plane_match_cost
   計算匹配代價；對應着Iteration部分的spatial propagation的思想“相鄰像素點間很
   可能具有相似的'平面' <=> '視差'”
*********************************************************************/
for(auto it = offsets.begin(); it < offsets.end(); ++it)
{
    std::pair<int, int> ofs = *it;

    int ny = y + ofs.first;
    int nx = x + ofs.second;

    if(!inside(nx, ny, 0, 0, cols, rows))
        continue;

    Plane p_neigb = planes[cpv](ny, nx);
    float new_cost = plane_match_cost(p_neigb, x, y, WINDOW_SIZE, cpv);

    if(new_cost < old_cost)
    {
        old_plane = p_neigb;
        old_cost = new_cost;
    }
}

// plane refinement
// 這部分比較直觀，就是空間上的二分隨機搜索，由隨機平面法向量和隨機視差確定，具體見理論部分
int sign = (cpv == 0) ? -1 : 1;	// -1 processing left, +1 processing right

float max_dz  = max_delta_z;
float max_dn = max_delta_n;

Plane& old_plane = planes[cpv](y, x);
float& old_cost = costs[cpv](y, x);

while(max_dz >= end_dz)
{
    // Searching a random plane starting from the actual one
    std::random_device rd;
    std::mt19937 gen(rd());

    std::uniform_real_distribution<> rand_z(-max_dz, +max_dz);
    std::uniform_real_distribution<> rand_n(-max_dn, +max_dn);

    float z = old_plane[0] * x + old_plane[1] * y + old_plane[2];
    float delta_z = rand_z(gen);
    cv::Vec3f new_point(x, y, z + delta_z);

    cv::Vec3f n = old_plane.getNormal();
    cv::Vec3f delta_n(rand_n(gen), rand_n(gen), rand_n(gen));
    cv::Vec3f new_normal = n + delta_n;
    cv::normalize(new_normal, new_normal);

    // test the new plane
    Plane new_plane(new_point, new_normal);
    float new_cost = plane_match_cost(new_plane, x, y, WINDOW_SIZE, cpv);

    if(new_cost < old_cost)
    {
        old_plane = new_plane;
        old_cost = new_cost;
    }

    max_dz /= 2.0f;
    max_dn /= 2.0f;
}

// view propagation，這部分邏輯有點繞，實際上就是左右視差一致性？
/*********************************************************************
1. 13行提取當前視角下當前點對應的平面
2. 19行將當前點視角變換后到對應視角的點計算出來(簡單理解為根據當前點平面查到當前點視差，
   然后找到對應視角找對應的點，將(x + sign * 視差, y, 視差)構造為新平面new_plane
3. 而對應視角下的點(x + sign * 視差, y)在上一次計算后會有自己的一個cost:old_cost
   對應着Iteration部分的view propagation的思想“左右視角下，匹配的兩個點具有相似的
   '平面' <=> '視差'”
*********************************************************************/
int sign = (cpv == 0) ? -1 : 1;	// -1 processing left, +1 processing right

// current plane
Plane view_plane = planes[cpv](y, x);

// computing matching point in other view		
// reparameterized corresopndent plane in other view
int mx, my;
Plane new_plane = view_plane.viewTransform(x, y, sign, mx, my);

if(!inside(mx, my, 0, 0, views[0].cols, views[0].rows))
    return;

// check if this reparameterized plane is better in the other view
float& old_cost = costs[1-cpv](my, mx);
float  new_cost = plane_match_cost(new_plane, mx, my, WINDOW_SIZE, 1-cpv);

if(new_cost < old_cost)
{
    planes[1-cpv](my, mx) = new_plane;
    old_cost = new_cost;
}