Binocular stereo matching algorithm based on MST cost aggregation

1.   Introduction

Image sensors are widely used in automatic and sensing devices to realize object detection and recognition, thanks to the fast development of imaging and computer technologies. Stereo matching is of great importance for computer vision appositions, such as autonomous driving ^[1,2], target detection and recognition ^[3,4]. It also plays important role in the fields of robot navigation ^[5], and space detection ^[6,7]. Various stereo matching algorithms have been developed so far. While most algorithms are cable of meeting the matching quality, it tends to have the problem of high cost of computer time. Real time applications demand enhanced matching quality and reduced processing time. Therefore, applicable methods are needed to meet and balance these two seemly contradictory requirements for real time processing.

Stereo matching algorithms can be classified into global and local algorithms, normally. The local algorithms tend to have lower computational complexity and lower accuracy comparing with global algorithms. Recently developed algorithms target to find balance between quality and computation time. Among those, the non-global stereo matching method based on Minimum Spanning Tree (MST) is proved to be efficient and advantageous ^[8]. The method however still suffers the problem of high computation time when deals with real time applications. Parallel processors and graphics processor (GPU) are great help to increase the computational speed. Especially parallel processor, having the advantage of low power requirements, are of many research concerns. There are some promising results have been obtained in the development of parallel-sequential matching algorithms.

In this paper, an advanced method is proposed to build the MST based stereo matching algorithm by Yang on parallel processors, which leads to reduced computation time at the same maintain the original algorithm's advantage of high accuracy.

2.   MST algorithm and cost aggregation

The weighted connected graph $G = (V, E)$ is supposed to have $n$ vertices, $\left| {V = n} \right|$. Then, the spanning tree of G is a minimal connected subetaaph ${G}^{\mathrm{\rm{'}}} = ({V}^{\mathrm{\rm{'}}}, {E}^{\mathrm{\rm{'}}})$ of G. where, ${V}^{\mathrm{\rm{'}}}$ satisfies ${V}^{\mathrm{\rm{'}}} = V$ and $E\mathrm{\rm{'}}$ satisfies ${E}^{\mathrm{\rm{'}}}\subseteq E\bigwedge \left|{E}^{\mathrm{\rm{'}}}\right| = n-1$.

If $G\mathrm{\text{'}}$ is the MST of $G$, the ${E}^{\mathrm{\text{'}}}$ satisfies ${E}_{0}\subseteq E$ for any ${E}_{0}$. When $\left|{E}_{0}\right| = n-1$, the equation (1) will always hold.

where, $W(e)$ represents the weight of an edge $e$.

The MST of a connected graph can be calculated by Prim algorithm or Kruskal algorithm ^[9]. The difference between these two algorithms is that the Prim algorithm starts from an origin, and continuously expands its range to select the edge with the smallest weight; The Kruskal algorithm looks at the whole, and constantly selects the smallest weight edge from the whole.

In most local stereo matching algorithms, the matching cost of a pixel $p$ is calculated within a certain domain window of $p$. In other words, the matching cost is determined by the pixels which are located inside the domain widow. The pixels which are outside of the domain widow have no influxes on the matching cost.

In order to make the calculation of matching cost global, that is, all pixels in the image can affect the matching cost of the pixel, image $I$ can be regarded as an undirected weighted connected graph $G = (V, E)$. The vertex set $V$ is all pixels of $I$ and edge set $E$ is the relationship between adjacent pixels of $I$. Therefore, $G$ is a 4-connected grid graph. The weight of the edge $(u, v)$ of $G$ is defined as the gradient of the pixel gray value as shown in equation (2).

According to the structure of MST, for the two points $p$ and $q$ in the image, if $q$ is far away from $p$ and the color differences between $q$ and $p$ is large, the number of path hops between $q$ and $p$ in MST is large and the influence on the calculation of matching cost of $p$ is weak.

Therefore, the distance $D\left(p, q\right)$ between the two points $p$ and $q$ in the MST can represent not only the spatial difference of pixels, but also the intensity difference of pixel values between them. The definition of similarity $S\left(p, q\right)$ between $p$ and $q$ is shown in equation (3).

where, $\sigma$ is the adjustment parameter of similarity. Therefore, under MST, bilateral filtering function used in cost agammaegation stage can be extended to equation (4).

where, $C_d^A(p)$ represents the agammaegation cost of pixels $p$ at parallax $d$, ${C_d}(q)$ represents the matching cost of pixels $q$ at parallax $d$, $c(p)$ represents the node connected with $p$.

The similarity of two pixels in an image depends on the distance between the two nodes in the MST. Therefore, the MST can be organized into a tree structure. The leaf to root (L2R) cost agammaegation process is used to calculate the cost of each pixel affected by the nodes in its subtree. The root to leaf (R2L) cost agammaegation process is used to calculate the cost of each pixel affected by nodes other than those in its subtree in MST. After L2R and R2L processes, the cost of each pixel affected by all other pixels can be obtained. The final matching costs were agammaegated over each influencing pixel and then over each pixel included.

3.   An improved BFS algorithm

Breadth first search (BFS) is a search algorithm in graph structure ^[10]. The algorithm starts from a source node, traverses the child nodes of the node in turn, and then traverses the child nodes of these child nodes, and so on, until traversing the complete graph. Therefore, BFS algorithm needs to use an open closed table to record the traversed nodes. The open records the nodes to be traversed, and the closed records the traversed nodes. Similar to the BFS algorithm, Level Synchronous Parallel BFS (LSP-BFS) algorithm maintains three node sets: visited node set $V$, current level node set $C$, and next level node set $N$. The algorithm takes out the elements in $C$ and places the adjacent nodes in $N$ in parallel. This process is iterated until $N$ is empty set. Before the next iteration process, let $C = N$, $N = \mathrm{\varnothing }$.

Since the process of cost agammaegation algorithm based on MST is carried out on the MST of the reference image, it is necessary to reduce algorithm to tree structure BFS algorithm. The BFS algorithm is mainly used for the general graphs which may have loops. The tree structure is an acyclic graph. The nodes in the tree structure have clear partial sequence relations. When a node is accessed, the parent node of the node must have been accessed, and the child node of the node must not have been accessed. Therefore, it doesn't need to check whether the adjacent points have already been accessed in the tree structure BFS algorithm. It can save the closed container.

3.1.   LSP-BFS algorithm for tree structure

BFS algorithm on tree structure can omit closed container. Similarly, LSP-BFS algorithm on tree structure can omit the visited node set $V$. In the extended traversal range, it only needs all the child nodes to join the next level node set $N$. It is not necessary to determine whether the child node has been traversed. The pseudo code of the improved LSP-BFS algorithm is shown in algorithm 1.

It can be seen that LSP-BFS algorithm has two characteristics:

(1) Each layer requires additional synchronization. Synchronization is a time-consuming operation. If the MST height is too high, more synchronization operations are required.

(2) The number of parallel operations in a BFS layer depends on the number of nodes in that layer. The BFS layer with more nodes is more parallel, which can take better advantage to the performance of parallel devices. There are four images in Figure 1. The four images in are converted into MST. The vertex corresponding to the upper left pixel as the root node. A tree is constructed with the root node. The results of the number of nodes in the tree hierarchy are shown in Table 1. The node hierarchy distribution is shown in Figure 2. The statistical method is to sum the number of nodes per 50 layers. In the distribution map, the abscissa is the serial number of the hierarchy (numbering every 50 levels), and the ordinate is the number of nodes in the hierarchy (summing the nodes of every 50 layers). It can be seen from Figure 2 that the number of nodes in higher level and lower level is less, and the number of nodes in middle level is more. Therefore, in general, the node distribution of tree structure is sparse in the higher level and lower level, and the middle level is more denser.

According to these two characteristics, in order to give full play to the performance of parallel devices, the parallel BFS algorithm needs to reduce the height of the tree, and increases the number of nodes in each layer of the tree. As shown in Figure 3, tree (a) and tree (b) originate from the same acyclic graph, but the root node of tree (a) is ${v}_{1}$, the root node of tree (b) is ${v}_{2}$. The height of tree (a) and (b) are 5 and 4, respectively. The node density in the second layer of tree (a) is significantly lower than that in the tree (b). Therefore, for the same spanning tree, when the root nodes are different, the lower the height of the tree, the higher the nodes density in the middle level.

When the number of nodes in the same layer is large enough, LSP-BFS can play the largest role. Therefore, in order to make LSP-BFS algorithm can be applied to solve the agammaegation cost of MST, it is necessary to reduce the height of tree. If the middle node ${v}_{r}$ of the longest path of acyclic graph is taken as the root node of MST, the MST with the lowest height can be obtained. Therefore, it is necessary to find the longest path of acyclic graph. A simple method is to use BFS twice to find the longest path. The basic process of the algorithm is as follows: first, let any point ${v}_{1}$ of the graph be starting point, use BFS to find the point ${v}_{2}$ which is farthest from ${v}_{1}$. Then from the point ${v}_{2}$ start, use BFS to find the point ${v}_{3}$ which is farthest from ${v}_{2}$. The path of ${v}_{2}$ to ${v}_{3}$ is the longest path of the graph.

This algorithm can find the longest path, but it is need to traverse all of the longest paths to find ${v}_{r}$. This process takes extra time. In order to solve this problem, this paper proposes a pruning method. The node ${v}_{r}$ can be found out in twice level traversal time by using the pruning method.

In an acyclic graph, nodes with only one adjacent node are called branch nodes. The node ${v}_{r}$ can be found by deleting these branch nodes. For a given acyclic graph $G = \left(V, E\right)$ and a source node $s$, the process of the algorithm is as follows:

(1) Starting from the source node $s$, all the branch nodes are searched by BFS, then all the branch nodes are put into the branch node set $N$ of the next layer;

(2) Let the branch node set $V = N$, $N = \mathrm{\varnothing }$. For each node $v$ in $V$, if the adjacent node of node $v$ is not in the next level branch node set $N$, then $v$ is added to $N$ and deleted from $G$.

(3) Step (2) is performed continuously until the number of nodes in $G$ is 1 or 2. If the number of remaining nodes is 1, the node is ${\upsilon _r}$. If the number of nodes is 2, any one of them can be regarded as ${\upsilon _r}$.

The pseudo code of pruning method is shown in algorithm 2. In this method, the time complexity of finding ${v}_{r}$ is $O\left(2l\right) = O\left(3l\right) = O\left(l\right)$. where $l$ is the longest path length of $G$. Because $G$ is traversed twice, the time complexity of pruning method is $O\left(2l\right) = O\left(l\right)$. The time complexity of the algorithm for finding the longest path through BFS twice is also $O\left(2l\right) = O\left(l\right)$. Therefore, the performance of the pruning algorithm is equivalent to that of BFS twice method. It does not need to find ${v}_{r}$ on the longest path, so the pruning method can reduce one traversal time, and the algorithm design is more simpler.

For the test samples in Figure 1, after reducing the height of the tree by using pruning method, the statistics of nodes for each layer are shown in Figure 4 and table 2. In the Figure 4, the black dotted line represents the node distribution of each layer when the tree height is not reduced; the red solid line represents the node distribution of each layer when the tree height is reduced. It can be seen from Figure 4 that as the tree height decreases, the density of nodes for the middle layer also increases.

It can be seen from table 2, the tree height of test examples (a) and (b) is reduced by about 45%, and the tree height of test examples (c) and (d) is reduced by about 25%.

3.2.   Improved L2R cost aggregation algorithm

The acyclic graph $G$ and a source node $s$ are obtained for a certain image. The improved L2R cost agammaegation algorithm can maintain a queue $Q$ and a linked list $L$, where $L$ is the storing linked list. The process of the algorithm is as follows:

(1) Put $s$ in the team;

(2) Use the variable $l$ to save the length of $Q$ and insert a new node $N$ at the end of $L$;

(3) Set an element from $Q$ as $e$, insert the $e$ into the tail of $N$, and queue all the child nodes of $e$;

(4) Perform step (3) $l$ times;

(5) Continue to perform steps (2) to (4) until $Q$ is empty;

(6) Traverse $L$ from the tail to the head. When each node of $L$ is traversed, the agammaegate cost of all nodes in the node is calculated. This step can be calculated in parallel.

Obviously, the nodes in the linked list L store the nodes of each layer of $G$. After traversing $G$, the length of $L$ is the height of the tree with $s$ as the root node. Variable $l$ stores the number of nodes in each layer. The pseudo code of the algorithm is shown in algorithm 3.

4.   Experimental results and analysis

All the algorithms are implemented and tested on Windows 10 operating system computer by using Visual C + + and NVIDIA CUDA. The hardware conditions of the computer are Intel (R) core (TM) i5-6300HQ CPU @ 2.30 GHz, NVIDIA GeForce GTX 960m graphics card and 8.00 GB memory. Middlebury stereo vision test set is widely used test set in academia. it includes test sets for binocular vision, multi vision and other technologies. Cones, Teddy, Tsukuba and Venus binocular vision test sets are used in this paper, as shown in Figure 5. The first row images are the reference images of the test set (the image of the left imaging system); the second row images are the target images of the test set (the image of the right imaging system); the third row images are the real disparity images. The parameters of the test set are shown in Table 3.

The process of parallel cost agammaegation based on MST is as follows: (1) median filtering; (2) calculating the weight of edges; (3) reordering by the weight; (4) generating MST by the Kruskal algorithm; (5) pruning and constructing tree; (6) L2R process; (7) R2L process.

Image median filtering is a nonlinear filtering method. The calculation formula for each pixel of image median filtering is shown in equation (5).

where, $k$ is the size of filtering window, the value of $k$ is odd. When the window is too large, the image will become blurred. The parameter selection in this paper starts from the minimum window (k = 1) and gradually grows.

The test results using image median filter are shown in Figure 6 and table 4. In the Figure 6, the images in the first row are real disparity images, and the images in the rows 2, 3, 4 and 5 are the result of setting $k = 1$, $k = 3$, $k = 5$, and $k = 7$, respectively. It can be seen from Figure 6 and table 4 that the image median filtering has a great influence on the results. When $k = 3$, the effect of image median filtering is better. Equation (6) and equation (7) are selected to test the edge weight. The test results are shown in Figure 7 and table 5. In the Figure 7, the images in the rows 1, 2 and 3 are the real disparity images, the result images of a function ${W}_{1}({p}_{1}, {p}_{2})$ is selected and the result images of a function ${W}_{2}({p}_{1}, {p}_{2})$ is selected, respectively. It can be seen from Figure 6 and Figure 7 that the effect is the better when the function ${W}_{2}({p}_{1}, {p}_{2})$ is selected.

When calculating MST by using Kruskal algorithm, it is need to sort all the edges in descending order of weight. The time complexity of common sorting algorithms is $O\left(nlog\left(n\right)\right)$. However, the weight of the edge will not be greater than 255 in this paper, so histogram sorting can be used. The time complexity of this method is $O\left(1\right)$.

The time test results of the algorithm are shown in table 6. The numbers in brackets in the table 6 represent the running speed ratio. The running speed ratio of MST algorithm using GPU is compared with the corresponding algorithm using CPU. The running speed ratio of the algorithm proposed in this paper is compared with the corresponding MST algorithm using GPU. It can be seen from table 6 that the MST algorithm can be accelerated by about 17 times by using GPU, and the optimization algorithm proposed in this paper can increase the speed by about 20% on this basis. Therefore, the optimization algorithm proposed in this paper can speed up about 20 times compared with the algorithm using CPU, and it can meet the requirements of real-time processing.

5.   Conclusion

An advanced method is developed in this paper to improve the computation efficiency for stereo matching. The proposed method is to reduce the height of the tree used for MST algorithm in parallel processing, hence dramatically speed up the computation. An agammaegation of matching cost is also developed simultaneously in order to achieve the hierarchical synchronous parallel computing. The cost is based on Leaf-to-Root agammaegation.

The simulation results show that proposed parallel processing stereo matching algorithm using reduced high of spanning tree and L2R agammaegation matching cost enhance the computation speed. It is a matter of 20 times speed-up as proved. This method therefore provides a promising approach to real-time stereo processing, which can of great importance of Binocular stereo matching.

Acknowledgments

This work was supported by the Natural Science Foundation of Jiangsu Province (BK20191012); Scientific Research Foundation of Nanjing Institute of Technology (JCYJ201822, CKJB201803 and CXY201932).

Conflict of interest

The authors declare there is no conflict of interest.