Real-time stereo matching architecture based on 2D MRF model: a memory-efficient systolic array

There is a growing need in computer vision applications for stereopsis, requiring not only accurate distance but also fast and compact physical implementation. Global energy minimization techniques provide remarkably precise results. But they suffer from huge computational complexity. One of the main challenges is to parallelize the iterative computation, solving the memory access problem between the big external memory and the massive processors. Remarkable memory saving can be obtained with our memory reduction scheme, and our new architecture is a systolic array. If we expand it into N ’ s multiple chips in a cascaded manner, we can cope with various ranges of image resolutions. We have realized it using the FPGA technology. Our architecture records 19 times smaller memory than the global minimization technique, which is a principal step toward real-time chip implementation of the various iterative image processing algorithms with tiny and distributed memory resources like optical flow, image restoration, etc.


Introduction
The stereo matching problem is to find the corresponding points in a pair of images portraying the same scene.The underlying principle is that two cameras separated by a baseline capture slightly dissimilar views of the same scene.Finding the corresponding pairs is known to be the most challenging step in the binocular stereo problem.
As shown in Table 1, the conventional methods can be categorized into the local and global methods [1].The unit, million disparity estimations per second (MDE/s), is the product of the number of pixels, disparity levels, and frame-rate and therefore, stands for the overall computational speed.Note that the global methods have the low throughput due to their small number of processors.
The local method, typically window correlation and dynamic programming (DP) methods, examines subimages only to obtain local minima as solutions.Inherently, this method needs relatively small operations and memory, making it the popular approach in real-time DSP systems [2,3] and parallel VLSI chips [4][5][6][7].The local method can be easily realized in the massive parallel structure as shown in Table 1.Nevertheless, there are many situations where this method may fail: the occlusion, uniform texture, ambiguity of the low texture, etc.Even further, the window method tends to yield blurred results around the object boundary.
In contrast, the global method, typically graph cut [8,9] and BP [10][11][12], deals with whole images, resulting in the global minima, analogously to the approximated global minimum principle.This approach has the advantage of low error rate but tends to need huge computational loads and memory resources.Recently, some researchers realized BP using PC aided by specialized parallel processors on GPU graphic card [13].As described in Table 1, the so-called real-time BP can yield reasonable results only for the small throughput (MDE/s).Unfortunately, the specialized GPU relies upon high speed clocks and a small number of processors, which cannot be regarded as fully parallel architecture.Thus, it has the throughput limitation.Nevertheless, this system is successfully used in the realtime computer vision area [14].There is no full parallel system that has fast computational power (MDE/s) for the high resolution images or the fast frame rates.Further, there is no genuine compact hardware dedicated to the global stereo matching in real time.Most of the existing systems are impractical in terms of size, power requirement, and expense and are not suitable for compact applications like robot vision.
If a massive parallel architecture is realized as shown in Figure 1 then the computational time may be reduced drastically.However, this global matching architecture is not workable simply because of the enormous data bus bandwidth between the processors and the big external memory resource.In an effort to avoid this bottleneck, the memories must be evenly distributed throughout the processors so that each processor may access its own memory unhindered by the others.This distributed approach also raises problems when the number of processors is excessively large and the memory size is too big, making the VLSI implementation a formidable task.Therefore, we need to use distributed internal memories of small size, which can be easily accessed by many processors simultaneously.
Consider the one chip solution with a systolic array and efficient memory configuration.To avoid the huge memory, we tried to implement the BP on the FPGA by reducing the memory size [15], which is similar to the hierarchical iteration sequence [16].In this paper, we use IF scheme [16] for our architecture and make it 2 times smaller than IF considering the message propagation direction, as we will call "Fast belief propagation (FBP)".Based on this method, we built a full parallel architecture that is efficient in memory usage as well as equivalent to the original belief propagation (BP) method in terms of accuracy.
For a real-time application with small and compact hardware, GPU-and CPU-based system is not good due to their bulky size.We used this architecture to build a stereo vision chip and observed the expected performance-realtime and small memory for high precision depth images.
The remainder of this paper is organized as follows.Section 2 explains the background of the belief propagation.Section 3 defines a layer structure and explains an FBP sequence.A new iteration filter algorithm considering iteration directions is described in Section 4. For a VLSI realization, Section 5 suggests a parallel architecture and its memory complexity.Experiments are presented in Section 6. Section 7 draws conclusions on our newly developed architecture.

Review of belief propagation
The basic concept of belief propagation (BP) is to find iteratively the maximum a posteriori (MAP) solution on a 2-D Markov random field (MRF).All the parameters and variables are defined on the 2-D graph Figure 2 (we use the notation from [10]).P: a set of nodes on 2-D MRF, which in fact corresponds to pixels on an image.D: a set of hidden states stored in the nodes.p P: a node that is located on the coordinate p = (p 0 , p 1 ).d p D: a hidden state at p. g l , g r : left and right images of N 0 by N 1 size.Also, N E denotes the edge set and therefore, (p, q) N E for an edge between two nodes p and q.
With the help of these notations, the pairwise MRF energy model can be defined as determining the estimate d, given an energy function E(•): D(d p ) is the data cost for the node p having the state d p .Similarly, V (d p , d q ) is the edge cost for a pair of neighbor nodes p and q having states d p and d q , respectively.
We assume a condition of parallel optics without the loss of generality.Then, stereo matching simply involves finding a point (p 0 , p 1 + d p ) in the right image which corresponds to a point (p 0 , p 1 ) in the left image.Thus, the hidden state d p represents the offset between the corresponding pixels, as is called disparity.
At each state d p , the data cost constrained by the left and right images is defined as where C d and K d are a weighting factor and upper bound of the cost, respectively.This upper bound is useful in making the data cost robust to occlusions and artifacts that may violate the common assumptions that the ambient brightness must be uniform.
Also, the disparity should vary smoothly almost everywhere except at some places like object boundaries.In order to allow this discontinuity, we keep the edge cost V (d p , d q ) constant whenever the difference becomes larger than the predefined parameter K d : where C v and K v are similarly defined as the constant.
Finding the state d with minimum energy in Equation 1 amounts to the estimation problem with MAP.As is well known, the approximated MAP solution d can be estimated using the following BP update [10]: N(p)\q is the neighbors of node p excluding q, a is the normalization value, and S is the state size.This equation expresses the following mechanism.The message m l pq (d q ) at node p is updated at time l and then sent to the neighbor node q.After L iterations, the expected dp at each node can be decided with Equation 7.
Let us explain the hierarchical BP in brief.It is based on the iteration scheme in multiple different scale levels.Between the levels, 2 × 2 scale change is considered to aid the coarse-to-fine iteration.According to this scheme, we need to over-sample the message and data costs in the coarse level to obtain the cost for the finer level.In this paper, and  D k p k denote the iteration number, the iteration time index, the node, message, and data cost in the M/2 k by N/2 k hierarchical graph of the scale level k [0, K -1], respectively.Here, K -1 means the coarsest level.As shown in Figure 3, the data cost at k is calculated from the data cost at k -1 by the summation over a 2 × 2 block.At the scale level 0, the data cost If the memory complexity at each node is B bits, the overall memory size is

The proposed fast belief propagation sequence
In this section, we propose our FBP algorithm and architecture that enable us to run the BP on the FPGA with tiny distributed RAMs and show the remarkable memory reduction.It is 2 times smaller than the Iteration Filter's memory reduction scheme [16].Before entering this section, I recommend for readers to understand the Iteration Filter scheme [16] that is wholly different from the normal iteration sequence and shows the amazing memory reduction effect.We redesign the Iteration Filter algorithm and implement it on the FPGA.
If we consider a separate layer for each iteration, then we can build a stack of layers.In this structure, the iteration can be represented as the upward propagation.Thus, Figure 4 can be redrawn as Figure 5. From this interpretation, we are considering the 2D graph with the iteration as the 3D layer graph (p 0 , p 1 , l) with the propagation.Let us define message and data cost sets at each node and layer l as: From these definitions, we can simplify the message update function in Equation 5 as: where (N(p), l -1) and M((N(p), l -1)) = {M(u, l -1)| u N(p)} represent the neighbor nodes and their message costs in the buffer, respectively.
As an initialization stage, each node p observes the input to obtain the data cost D(p, 0).Afterward, in every iteration l, each node calculates the new message M(p, l) according to the update function f(•) and after then stores it as M(p, l -1) in the buffer.
Let Q(l) and M(Q(l)) denote the set of nodes in lth layer and its message cost set, respectively.Then, M(Q (l)) can be updated from M(Q(l -1)) and D(Q(l -1)) in the buffer: Consider a new FBP computing order based on the IF scheme.Note that Q(p 0l, l) forms a linear array of M nodes on the p 1 axis in the lth layer.If we collect all the layers of Q(p 0l, l) in terms of p 0 then Q(p 0 ) forms a planar array of LM nodes: with the notation Q(p 0l, l) and Q(p 0 ), we can build an efficient computation order.We will call this memory-efficient BP sequence, FBP.The cost of Q(p 0 ) is updated from the buffer of the message M (Q(p 0 -1)), M(Q(p 0 -2)), and data cost D(Q(p 0 -1)) as described in Algorithm 1.As shown in Figure 6, our memory resource consists of local and layer buffers.The layer buffer stores all the layers' costs of Q(p 0 -1) and Q(p 0 -2).The local buffer holds only one layer's costs on Q(p 0 , l -1).
Algorithm 1: FBP algorithm For ℓp 0 in the lth iteration layer profile, each node at (p 0l, p 1 ) and the lth layer can be updated from the node at N(p 0l, p 1 ) and the (l -1)th layer.Thus, as shown in Figure 7 and Equation 17, the nodes at Q(p 0 , l) can be computed from Q(p 0 , l -1), Q(p 0 -1, l -1), and Q(p 0 -2, l -1).
Q(p 0 , l) and Q(p 0 , l -1) belong to Q(p0).Hence, given the layer buffer Q(p0 -2) and Q(p0 -1) and the local buffer Q(p 0 , l -1), the costs in Q(p 0 , l) are updated at each layer l recursively, which sequence is described in Figure 6a, b, and 6c.That is, given M(Q(p 0 -1)), M(Q (p 0 -2)), and D(Q(p 0 -1)), we can calculate M(Q(p 0 )).The new costs in local buffer should be stored in the layer buffer to process the next set Q(p 0 + 1) in the next time.This sequence shifts the layer buffer to the p 0 axis direction.Then, for p 0 from 0 to N + L -1, we can obtain the final iterated message M(Q(p 0 , L)).For the example, as shown in Figure 6b, and 6c, the location of the buffer is changed from Q(p 0 = 5) to Q(p 0 = 6) by our sequence.
In the hierarchical case, as shown in Figure 6d, we can construct the hierarchical layer structure by considering the hierarchical iterations.At each level, we can follow the FBP sequence at each level only if considering two by two scale changes between levels.Please refer to [16] for the detailed hierarchical memory reduction scheme of IF.
If we use the notation B as BP memory complexity at each node and consider the nodes of L k by M/2 k size in Q k (•), we need two layer buffers of the BL k M/2 k size and one local buffer of BM/2 k size at each level k.Thus, compared with the hierarchical BP, the overall memory size can be reduced from ) bits by adopting the iteration filter scheme to our VLSI sequence.This can be shown as follows.
If we approximately consider the total memory as the 0th level, the reduction rate amounts to N/(2L 0 + 1) times when 2L 0 ≪ N. In summary, the update sequence must be effective whenever N, one of the image size components is big, and L 0 , the iteration number, is small.

New iteration sequence considering the iteration direction
Let us consider the message propagation direction for the further memory reduction.As shown at the definition of M(p, l) in Equation 9, we assumed that the messages of all the directions are stored in the buffer.However, due to the message propagation direction information, we can reduce the memory resource 2 times smaller.Among the neighbor messages M(N(p), l -1), only m l−1 rp (d p ) for r N (p) is necessary for updating M (p, l).In Figure 8, let us denote the message propagation direction as Δ = p -N(p).The needed messages for the update are the ones that are propagated from neighboring node N(p) to p. Except for the message of the direction Δ = [+1 0] that is propagated from local buffer, all the other messages are being loaded from the layer buffer.This is summarized at the access column part of Table 2. But, in the data cost case, as shown in Figure 9, we do not need to consider the propagation direction and simply read D(Q(p 0 -1, l -1)) in the layer buffer Q(p 0 -1) for D(Q(p 0 , l)) because D(Q(p 0 , l)) is equal to D(Q(p 0 -1, l -1)) like Equation 12.
As explained in the FBP algorithm, at each update time, the location of the buffer is shifted to p 0 axis being updated by the new cost.The newly updated messages and data cost in the local buffer should be stored in the layer buffer for the processing of the next Q(p 0 + 1).Thus, if the messages from all possible directions be saved in the local buffer, then some messages can be transferred to Q(p 0 -1, l -1).At the same time, some old costs in Q(p 0 -1, l -1) are moved to Q(p 0 -2, l -1) in a similar way.With this scheme, the number of propagation directions to be stored at the buffer is described at the store(Δ) part in Table 2.
From the definition in Equations 15 and 16, the number of nodes is LM for both Q(p 0 -2) and Q(p 0 -1) and M for Q(p 0 -(l -1), l -1).Table 2 shows the required number of messages and data costs at each node.The number of states is S, and the number of bits for the message cost and data cost is B m and B D , respectively.Then, by multiplying all the parts, we can calculate the memory size of the buffer as shown in Table 3.
If B = 4B m S + B D S, then we can obtain as follows: If you compare Equations 20 and 22, the value a is changed from two to one.Therefore, due to the propagation direction of BP, we can obtain 2 times smaller memory than the iteration filter [16].

Systolic VLSI architecture
Our architecture has four hierarchical levels.This level affects the iteration times.The higher hierarchical levels make iteration times smaller because the message can be converged faster in the coarse level.In our FBP architecture, it makes the memory size much smaller because our memory resource is dependent on iteration times.The HFBP algorithm can be easily realized with a systolic array architecture.As depicted in Figure 10, it consists of identical PE groups with nearest neighbor communication.In our implementation, it has a total of 20 PE groups.The PE group is divided into eight identical PEs as shown in Figure 11.Therefore, it amount to 160 PEs for processing a pair of 160 × 240 images.Figure 12 represents the local and layer buffer assignment for each PE k = 1,...,7 in the PE group.Thus, the 8/2 k number of PEs in the group is activated at level k due to the scale-down of the hierarchical structure.
As shown in Figure 11, the PE group consists of two parts.The first part is the data cost module that computes the initial costs using the left and right scan lines of the images.The other group is for updating the message and data cost.The pixel data from the left and right cameras enter into the PE group and each PE computes the data cost and the new message using the old messages from neighboring PEs and its own buffers.Figure 13 shows the data cost module that calculates the hierarchical data costs along the levels 0 to 3. In Figure 13b, the left and right scan lines are first stored in the registers, and then the right scan line registers are shifted by state d to compute D p (d p ) according to Equation 3.For each state, the data cost D p (d) at level 0 is obtained by taking the absolute difference of the left and right pixel values.On the other hand, B in Figure 13c is used for computing the higher level data cost D k p k (d).For the level k's cost, the previous level k-1 data costs are summed up and then accumulated over 2 k scan lines.This is equivalent to applying the summation of the 2 k × 2 k window for the hierarchical data cost; each data cost is used by the PE at each level.Data costs at each level, computed in the data cost module, are processed and saved in the corresponding PEs and buffers.See Figure 13.As described in Figure 12, the multiplexer (MUX) selects the messages and data costs at each level from which new messages and data costs can be updated and saved at the local buffer.Meanwhile, the old costs in this buffer are shifted into the layer buffer.In the four scale levels, 4-to-1 message multiplexer (MUX) is used.
For S number of states, the time complexity O(S) is needed to update one message at each node by forward, backward, and normalization operations [10].Normally, it needs 3S steps.As explained in Equation 9, four messages that are propagated to neighbor nodes need to be computed at each node.To compute these messages, our system needs only 6S clocks due to the pipeline structure.See Figure 14.
Since (M/2 k ) nodes are handled by (M/2 k ) processors in parallel on p k 1 axis, the total required clocks are reduced from As a whole, each PE calculates the messages in parallel by accessing the local buffer or the layer buffer which is located in the neighboring PEs or PE groups.

Experimental results
Our new architecture has been tested by both a simulation and FPGA realization.

Software simulation
First, we verify our VLSI algorithm using the Middlebury data set with a software simulation.In the previous sections, we presented a new architecture which is  equivalent to HBP in terms of input-output relationship and which is a systolic array with a small memory space.Hence, it is suitable for VLSI implementation.
The requirement for both memory resource and computation time is only dependent on the layer number L k .Therefore, it is reasonable to analyze the performance in terms of iterations as well as various images.We specify the accuracy using the following equation.where d is the estimated disparity, d True is the true disparity, Pm is the area except for the occlusion part, and N is the pixel number in its area.This error means the rate where the disparity error is larger than 1.
For fair comparison, the same parameters are used throughout the experiments: and K d = 60.Figures 15 and 16 are the results of the Middlebury test images.In Figure 15, four levels are used both for HBP and HFBP.The layer number at each level is assigned as (8,8,8,8) from coarse-to-fine scale levels.With the same iterations, HFBP and HBP show the same lower error results.
Figure 16 shows the relationship between the iteration layers and FBP's average memory reduction rates when compared with HBP, where the same iteration times, (L, L, L, L), are applied for each layer.Due to the hierarchical scheme, the iteration converged around 28 iterations and yielded 0.8% maximum error.The remarkable result, though, is the memory reduction, which is around 32 times.In fact, even less memory is possible for a higher error rate.Thus, this architecture makes the performance scalable between the space and accuracy.
Table 4 compares our FBP FPGA with other real-time systems in terms of error.It is evident that our method shows almost the same error as Real-time BP.Here, real-time BP is also based on the HBP algorithm [ and known for the lowest error among real-time systems.
If we use Equation 22, the total buffer size becomes 3.3 Mb, which is 19 times smaller than HBP's 62 Mb.Also, for processing one frame image, the 160 PEs need 0.6 MHz clocks.This speed amounts to 18.8 MHz clocks processing 15 frames in 1 s.In order to achieve maximum 36.8MDE/s throughput for a 160 × 480 image, only a 18.8 MHz system clock is necessary ideally.Tables 5 and 6 show the computational performance between our new system and other systems.The local matching is effectively implemented as the pipeline and parallel structure since it does not need to access the huge memory size iteratively.GPU is the SIMD processor with a high speed core clock and external memory clock.Even if it is not a full parallel structure, it operates in real time due to the high clock speed and small number of parallel processors.But, our system is the fully parallel and can operate at the much slower 25 MHz clock speed.Furthermore, our system has one chip solution that consumes less memory resources inside the FPGA and can easily be parallelized to multiple chips due to the systolic array architecture.This simple and regular architecture is suitable for VLSI implementation.In addition, the semi-global matching [17] needs two frames' latency times, but our FBP has the latency time below one frame due to the processing sequence like the filter.For a higher resolution solution, we need to increase the computational power.It is possible by simply cascading several chips together in proportion to the image size or increasing the clock speed.
It has been observed that the FPGA, incorporating 160 PEs, operates at a 25 MHz clock rate.For convenience, more specifications are summarized in Table 7. Ideally, to store the local and layer buffers, our necessary memory size is around 3.3 Mb.But, in the real implementation, we used 395 internal block RAMs in FPGA, which amount to 7.1 Mb.Incidentally, assigning each buffer to Block RAMs may result in unused leak memory, that is waste, that can be avoided in full ASICs.The new architecture is implemented in FPGA as shown in Figure 17.Here, Figure 17a is a block diagram and Figure 17b is a photo of the actual board.As can be seen, two cameras supply a pair of video streams and two FPGAs perform preprocessing and our FBP algorithm.The disparity map forms a stream from FPGA to a grabber through Camlink cables.From the video RAM on the grabber board, the PC reads the disparity data and converts it to a gray scale image for the observation.Figure 18 shows the typical video output of the FPGA.

Conclusions
In this paper, a new architecture for the global stereo matching algorithm has been presented.The key idea is to rearrange the computation order in BP to obtain a parallel and memory-efficient structure.As the results show, our system spends 19 times less memory than the ordinary BP.The memory space can be negotiated with the iteration number.The architecture is also scalable in terms of image size; the regular structure can be easily expanded by cascading identical modules.
When applied to binocular stereo vision, this architecture shows the ability to process stereo matching in real time.Experimental results confirm that this array architecture easily provides high throughput with low clock speed where small iterations are guaranteed by the hierarchical iteration scheme.
In the future, we plan to realize this architecture with a small and compact ASIC chip.Beyond the programmable chips, we can simply expect a real-time chip with higher resolution and the lowest error rate with huge PE numbers.Unlike the bulky GPU and CPU systems, making the complex stereo matching system with a compact chip may lead to many real-time vision applications.
Furthermore, if we change the message and data cost model, our memory-efficient architecture can be considered to other BP-based motion estimation and image restoration [10].The combined effort of parallel processing and efficient memory usage makes a chance to implement a compact VLSI chip.Furthermore, more general iterative algorithms can be considered, which communicate only neighbor pixels in the image, such as GBP typical cut [18].As explained in [16], if we apply the IF scheme to these algorithms, we can reduce their memory resources to a tiny   size.Thus, if they have simple update logics for the iteration, then full parallel VLSI architectures may be realizable.

Figure 1
Figure 1 Alternative architectures for parallel algorithms.(a) Parallel processors with a global memory.(b) Massive systolic array processors.

Figure 2 1 Figure 3
Figure 2 A 2-D regular graph which corresponds to a 2-D image.

Figure 10
Figure 10 Systolic array architecture of FBP.

Figure 13 Figure 14
Figure 13 Architecture for hierarchical data cost module.(a) Hierarchical summation.(b) Data cost module A. (c) Summation module B.
(a) Overall system (b) Hardware board

Figure 17
Figure 17 The overall hardware system.(a) Overall system.(b) Hardware board.

Table 1
Comparison of several real-time stereo systems

Table 3
FBP buffer size 10]Figure12Activated PE and hierarchical buffer assignment at each level in the PE group.

Table 4 Disparity
error comparison of several real-time methods (%)

Table 5
Comparisons of computation time between the real-time systems

Table 6
Comparisons of hardware spec.between the realtime systems

Table 7
Additional hardware specifications used in our system