SD-MVS: Segmentation-Driven Deformation Multi-View Stereo with Spherical Refinement and EM optimization

Some recent methods (Wang et al. 2021; Yuesong Wang et al. 2023) attempt to leverage patch deformation to improve matching cost or propagation scheme. As shown in Fig. 1, due to their insufficiency in exploiting edge information, they often cross boundary and reference areas with discontinuous depths, thereby yielding unsatisfactory results, especially when confronting with scenarios characterized by extensive similar colors and occlusions like forests and farmlands. Simultaneously, superpixel-based segmentation approaches (Romanoni and Matteucci 2019) also struggle in precisely recognizing certain critical edges within these scenarios. They also lack instance semantic information to broaden receptive field, thereby meet pixelwise characteristic.

SAM segmentation can mitigate this issue as it separates different instances to extract subtle edges information while neglecting robust illumination disturbance. Consequently, we can leverage instance segmentation to better exploit and further introduce edge information into patch deformation. Specifically, we perform instance segmentation using SAM for input image $I_{i}$ to generate masks for diverse instances, denoted as $\mathcal{F}$. Hence we have $M=\mathcal{F}(I_{i})$, where $M$ is an image mask whose size is consistent with $I_{i}$.

For each pixel $p$, we compute the bilateral weighted adaption of normalized cross correlation score (NCC) (Schönberger et al. 2016) between reference images $I_{i}$ and source image $I_{j}$, which can be calculated as follows:

$$\rho\left(p,\mathrm{W}_{p}^{i}\right)=\frac{cov\left(\mathrm{W}_{p}^{i},\mathrm{W}_{p}^{j}\right)}{\sqrt{cov\left(\mathrm{W}_{p}^{i},\mathrm{W}_{p}^{i}\right)cov\left(\mathrm{W}_{p}^{j},\mathrm{W}_{p}^{j}\right)}}$$

(1)

where $cov$ is weighted covariance, $\mathrm{W}_{p}^{i}$ and $\mathrm{W}_{p}^{j}$ are respectively the corresponding images patches on image $I_{i}$ and $I_{j}$.

The goal of minimizing the matching cost is to obtain the optimal matching depths via the computation of color differences. However, when objects with varying depths exhibit similar colors, they are susceptible to generating matching inaccuracies, as shown in Fig. 4(a). Therefore, we introduce patch deformation to compute matching cost upon the sample patch $\mathrm{W}$ intersecting with different instances.

Specifically, we first measure the distances from the corresponding central pixel $p$ to the left, right, lower and upper boundaries of $M$, denoted respectively as $d_{l}$, $d_{r}$, $d_{d}$, and $d_{u}$. Then we can deform the shape of $\mathrm{W}$ to match these boundaries. The new shape of deformed patch can be defined as:

$$\left[\frac{d_{l}+d_{r}}{d_{l}+d_{r}+d_{d}+d_{u}}L,\frac{d_{d}+d_{u}}{d_{l}+d_{r}+d_{d}+d_{u}}L\right]$$

(2)

where $L$ denotes the side length of the square patch before patch deformation. Additionally, we reposition the patch’s center by adding an offset:

$$\Delta o(p)=\left(\frac{d_{l}-d_{r}}{d_{l}+d_{r}}L_{h},\frac{d_{d}-d_{u}}{d_{u}+d_{d}}L_{v}\right)$$

(3)

where $L_{h}$ and $L_{v}$ are respectively the horizontal and vertical length of deformed patch. The new center of the sample patch now becomes $p+\Delta o(p)$.

Both patch deformation and center offset allow pixels positioned at boundary regions to orient their patches more intensively towards the center of its own instance. Enhancing the receptive field for homogenous pixels in such approach can yield more robust results, consequently reducing potential errors in estimation. Note that considering the runtime, we restrict the number of calculations for each window such that the number of calculations after deformation never surpasses the initial total number $(L/2)^{2}$.

After SAM-based instance segmentation, pixels within the same instance typically exhibit similar depths, whereas noticeable depth discontinuities frequently arise at the boundaries between instances. Considering that propagation involves updating potential depths and normals within the surrounding area for each pixel, depth discontinuities will inevitably impact propagation. Consequently, we leverage patch deformation to adaptively alter the propagation scheme.

The adaptive checkerboard propagation scheme (Xu and Tao 2019) is conducted by introducing the optimal hypotheses from four near and four far search domains, as illustrated in Fig. 5 (a). However, his search domain between two adjacent diagonal directions is too dense, which leads to an imbalanced search space density and a risk of selecting redundant values. Hence we modify its oblique direction into a straight line extending to the corner of each patch.

Subsequently, we propose patch deformation on propagation via SAM, which adjusts the propagation patch shape and direction for each pixel. As illustrated in Fig. 5 (b), we adapt the propagation directions according to the shape of the surrounding mask. Specifically, denoting $l_{l}$, $l_{r}$, $l_{d}$, and $l_{u}$ as the length from the central pixel $p$ to the left, right, lower and upper edges of the patch, respectively, we obtain:

$$l_{u}=\frac{d_{u}}{d_{u}+d_{d}}L_{v},l_{l}=\frac{d_{l}}{d_{l}+d_{r}}L_{h}$$

(4)

Both $l_{r}$ and $l_{d}$ can be obtained similarly. Therefore, the directions and lengths of slanted branch $l_{ul}$ is given by:

$$l_{ul}=\sqrt{{l_{u}}^{2}+{l_{l}}^{2}},\alpha_{ur}=\mathrm{arc}\tan\left(\frac{l_{u}}{l_{l}}\right)$$

(5)

where $l_{ul}$ refers to the length of the up-right branch, and $\alpha_{ur}$ represents the angle between the upward branch and the up-right branch. Corresponding lengths and directions of other three slanted branches can be obtained similarly.

Having adjusted all directions and lengths, we encounter another challenge: the searching domain for each branch is unbalanced. Since the process of selecting a pixel with the minimal cost is essentially a spatial neighborhood search, an imbalance will emerge due to the different length of branches. The search along a shorter branch is suffered from unreliable results due to its minor search domain.

To address this, we accordingly modify the searching strategy in the propagation scheme, as shown in Fig. 5 (c). Specifically, we employ eight different colors to depict separate search domains on the eight directions centered on $p$. Instead of taking the central pixel $p$ as the dividing point, we use the midpoint of the sum of the lengths of two opposite branches to divide the search domain. In experiments, pixels with the same color are grouped into the same domain, with CUDA operators balance the load of searching for minima within each color-specific region. Therefore, our proposed strategy ensures load-balance across all directions and allows for faster convergence.

Many conventional methods adopt cascading architectures by sequentially loading different scales of images into GPU, as shown in Fig. 6 (a). This may result in a time-consuming performance due to the limited transfer speed between CPU and GPU. Therefore, we draw inspiration from mipmap (Williams 1983) in computer graphics, a technique to load different scales of images in parallel at once, to replace the previous cascading architecture into our proposed parallel architecture.

Specifically, we first perform image downsampling in the CPU. Subsequently, multi-scale images are assembled and loaded together into the GPU, as depicted in Fig. 6 (b). Then multi-scale images are processed together through matching cost, propagation and refinement in the GPU. Finally, all predicted depth images are transferred back into the CPU. Denoting the maximum memory consumption of ACMMP cascading architectures as $\sigma$, and the number of memory read operations as $k$, this technique enables us to load all scales of images in the GPU memory at a reasonable cost of $\frac{4}{3}\sigma$ instead of sequentially loading images, thereby eliminating the need for $k-1$ additional memory read operations.

Based on this architerture, we further introduce multi-scale consistency on matching cost and propagation. Regarding matching cost, we first apply SAM segmentation on the $k$-th level downsampled image. Based on segmentation results, we construct deformed patch and further compute $k$-th level matching cost, denoted as $c_{k}$. Therefore, the multi-scale matching cost is given by:

$$C_{ms}=\frac{\sum_{k}{c_{k}}}{k}\,\$$

(6)

Concerning with propagation, the multi-scale consistency aggregates the search domain for all scales in each direction, yielding a total of eight distinct search domains. Conclusively, eight values with the lowest cost within each domain are chose as new hypothesis for further computation.

During the patch-matching phase, we consider not only the multi-scale matching cost $C_{ms}$, but also the reprojection error $C_{rp}$ and the projection color gradient error $C_{pc}$. $C_{rp}$ proposed in ACMMP validates depth estimation from geometric consistency. $C_{pc}$ measures color consistency between current pixel $p_{i}$ in reference image $I_{i}$ and its corresponding pixel $p_{j}$ in source images $I_{j}$:

$$C_{pc}=\max\left\{\left\|\nabla I_{j}\left(p_{j}\right)-\nabla I_{i}\left(p_{i}\right)\right\|,\tau\right\}$$

(7)

where $\nabla$ represents the Laplacian Operator, $p_{j}$ denotes pixel in image $I_{j}$ the projected by pixel $p_{i}$ in $I_{i}$, and $\tau$ is the truncation threshold to robustify the cost against outliers. With these terms, our the aggregated costs $C_{ag}$ can be given by:

$$C_{ag}=w_{ms}C_{ms}+w_{rp}C_{rp}+w_{pc}C_{pc}$$

(8)

where $w_{ms}$, $w_{rp}$, and $w_{pc}$ respectively represent the aggregation weights of each component.

As shown in Fig. 7, given the normalized normal, we first randomly choose two orthogonal vectors, $e_{1}$ and $e_{2}$, perpendicular to the normal $n$ as the perturbation direction. We then use the angles $\theta_{1}$ and $\theta_{2}$ as the degree of rotation for iterative refinement. The normal first undergoes a counterclockwise rotation by $\theta_{1}$ degrees around $e_{1}$ as the rotation axis. Subsequently, the normal is further rotated counterclockwise by $\theta_{2}$ degrees around $e_{2}$ as the rotation axis. According to Rodrigues’ rotation formula, the ultimate updated normal $n^{\prime\prime}$ is given by:

$$\begin{cases}n^{\prime}=cos\theta_{1}\cdot n+sin\theta_{1}(e_{1}\times n)\\ n^{\prime\prime}=cos\theta_{2}\cdot n^{\prime}+sin\theta_{2}(e_{2}\times n^{\prime})\end{cases}$$

(10)

This is analogous to sliding a vertex directed by the normal on the surface of a sphere, which ensures the preservation of normalization for the normal vector both before and after rotation. By finding two orthogonal bases perpendicular to the normal for refinement, it can be ensured that perturbations in each direction are equivalent. This approach aligns more closely with the geometric essence of the normal, which is defined on a sphere rather than individual axes in the $xyz$ coordinate system. As a result, our approach boosts the robustness and stability during the refinement process.

We also utilize gradient descent in our method. The primary merit of gradient descent lies in its ability to logically restrict the search space to the vicinity of probable solutions. Denoting the number of total iterations as $N_{max}$, the rotation angle $\theta$ for the $i$-th round is randomly selected from range $\left[0,5*2^{N_{max}-i}\right]$. After one round of refinement for depth $d$ and normal $n$, we determine the new direction for local perturbations $e_{1}^{\prime}$ and $e_{2}^{\prime}$ based on the result of the previous search. As such, we get:

$$\begin{cases}e_{1}^{\prime}\leftarrow n^{\prime\prime}-n\\ e_{2}^{\prime}\leftarrow e_{1}^{\prime}\times n^{\prime\prime}\end{cases}$$

(11)

Here, $e_{1}^{\prime}$ is aligned with the vector sum of the previous round’s perturbation, while $e_{2}^{\prime}$ is a vector perpendicular to both $n^{\prime}$ and $e_{1}^{\prime}$, as shown in Fig. 7(c). The primary merit of gradient descent lies in its ability to restrict the search domain of neighbourhood solutions. Each round of search takes place on the orthogonal plane defined by the previous search direction and the current normal direction, thereby enabling faster convergence to the optimal solution.

ACMMP employs a fixed interval for local perturbations on depth, while static perturbation range cannot adapt well to locally varying scene depth. Addressing this, we introduce pixelwise depth search interval chosen within the deformed patch.

Specifically, for each pixel, we extract the depth values of all pixels encompassed by its deformed patch, and choose the maximal and minimal values from this set as depth boundary for perturbations. Additionally, considering our iterative refinement strategy, during the $i$-th iteration, the pixelwise search interval is chosen within the deformed patch gained from $i$-th downsampled image, thereby narrowing the perturbation interval to yield more accurate hypothesis.

By fixing $w_{ms}$, $w_{rp}$, and $w_{pc}$, we can optimize the aggregated cost $C_{ag}$, formulated as:

$$\underset{C_{ms},C_{rp},C_{pc}}{\min}C_{ag}=w_{ms}C_{ms}+w_{rp}C_{rp}+w_{pc}C_{pc}$$

(12)

After optimization, we can get the optimal depth estimation under current hyperparameters.

By fixing $C_{ms}$,$C_{rp}$ and $C_{pc}$, we can optimize $w_{ms}$, $w_{rp}$ and $w_{pc}$, defined by:

	$\displaystyle\underset{w_{ms},w_{ms},w_{pc}}{\min}$	$\displaystyle C_{ag}=w_{ms}C_{ms}+w_{rp}C_{rp}+w_{pc}C_{pc},$		(13)
	$\displaystyle s.t.$	$\displaystyle w_{ms}+w_{rp}+w_{pc}=1,$
		$\displaystyle w_{ms},w_{rp},w_{pc}>\eta$

All hyperparameters are required to exceed a minimal value $\eta$, and we implement a normalization constraint ensuring that their sum equals $1$ to mitigate significant variances. Following the E-step optimization, we can alternatively optimize the hyperparameters and feed them back into the E-step for the next round of aggregated cost optimization.

Since it may be challenging to obtain the analytical solution to the optimization problem in M-step, we will use numerical optimization methods such as Newton’s method (Qi and Sun 1993) to obtain the optimal solutions for $w_{ms}$, $w_{rp}$, and $w_{pc}$. A comprehensive formula derivation of the optimization can be found in supplementary material.

In practical situations, there might be partial pixels with depth estimation errors when all pixels are selected. Hence, we only select pixels where SIFT features can be matched between different images, and then calculate the aggregate cost between the pixels corresponding to these features.

In terms of matching cost, we respectively remove patch deformation (w/o. Adp. Cost), multi-scale consistency (w/o. Mul. Cost) and both of them (w/. ACM. Cost). Since w/. ACM. Cost has neither deformable nor multi-scale, it produces the worst results. w/o. Mul. Cost slightly outperformed w/o. Adp. Cost, yet both are inferior to SD-MVS, implying that patch deformation contribute more than multi-scale consistency.

In terms of propagation, we respectively remove patch deformation (w/o. Adp. Pro.), multi-scale consistency (w/o. Mul. Pro.) and apply propagation scheme from ACMMP (w/. ACM. Pro.). Given that patches in ACMMP do not deform in accordance with the patch, its performance fell short of expectations. Both w/o. Adp. Pro. and w/o. Mul. Pro. delivered similar results, yet fell short in comparison to SD-MVS, indicating that both patch deformation and multi-scale consistency on propagation are equally crucial.

In terms of refinement, we respectively remove refinement (w/o. Ref.), exchange the refinement module into Gipuma (Galliani, Lasinger, and Schindler 2015) (w/. Gip. Ref.) and switch the refinement module into ACMMP (w/. ACM. Ref.). As observed, the absence of refinement significantly diminishes the results. However, introducing Gipuma refinement brings about noticeable progress, with further advancements achieved after adopting ACMMP refinement. Nonetheless, both refinement methods are worse than SD-MVS, proving the necessity of spherical gradient refinement.

We conduct two experiments (w/o. EM A and w/o. EM B) by removing EM-based Optimization and respectively setting $(w_{ms},w_{rp},w_{pc})$ to $(1,0.5,0.5)$ and $(1,0.2,0.2)$. The results highlight the impact of hyperparameter settings on the final results. Furthermore, their inferior performances compared to SD-MVS evidences the importance of automatic parameter tuning by the proposed EM-based Optimization.