Stereo Matching by Hierarchical, Microcanonical Annealing

STEREO MATCHING BY HIERARCHICAL, MICROCANONICAL ANNEALING

Stephen T. Barnard

Artificial Intelligence Center SRI International, Menlo Park, C~lifornia 94025

April 16, 1987

h and lR are the left and right image intensities, or some sim­ ple function of them (such as a Laplacian). Subscripts i and j range over all sites in the left image lattice. D is the disparity map. The first term is the absolute difference in intensity between corresponding points (where the correspondences are defined by the current state of the. system), and the second is proportional to the local spatial variation in the current state (the magni­ tude of the gradient of the disparity map). The constant >. is used to balance the terms. 1 This equation expresses two com­ peting constraints: first, the image intensities of corresponding points should be more-or-less equal, and second, the disparity map should be more-or-less continuous.

The Hamiltonian,

ABSTRACT

An improved stochastic stereo-matching algorithm is pre­ sented. It incorporates two substantial modifications to an ear­ lier version: a new variation of simulated annealing that is faster, simpler, and more controllable than the conventional "heat­ bath" version, and a hierarchical, coarse-to-fine-resolution con­ trol structure. The Hamiltonian used in the original model is minimized, but far more efficiently. The basis of microcanoni­ cal annealing is the Creutz algorithm . Unlike its counterpart, the familiar Metropolis algorithm, the Creutz algorithm simu­ lates a thermally isolated system at equilibrium. The hierarchi­ cal control structure, together with a Brownian state-transition function, tracks ground states across scale, beginning with small, coarsely coded levels. Results are shown for a 512 x 512 pair with 50 pixels of disparity.

I INTRODUCTION

The paper describes a stochastic optimization approach to stereo matching. An early version is described in [1). Two major extensions have been made that permit substantially im­ proved performance: First, a new variety of simulated annealing is employed which is more efficient, more easily implemented, and offers more control over the annealing process than the standard form. Secondly, the system is extended to operate over several levels of resolution, exploiting relatively quickly computed min­ ima at lower levels of resolution to initialize its state at higher levels. This method leads both to more efficiency and to the ability to deal with much larger ranges of disparities.

The basic representation remains unchanged. The state of the system encodes a dense map of discrete horizontal disparities , defined over the left image, which specify corresponding points in the right image .. In the improved design, this state is relative to a particular level of resolution.

The energy of a lattice site is composed of two terms, each of which expresses a constraint important in stereo matching:

Support for this work was provided by the Defense Advanced Research Projects Agency under contracts DCA 76-85-C-0004 and MDA903-86-C0084.

is unchanged; therefore, the ground states that we seek to deter­ mine or to approximate are the same. The new design is strictly concerned with improved performance. The techniques used to obtain it are not restricted to stereo matching, and should be con­ sidered for any simulated-annealing approach to low-level vision problems.

II MICROCANONICAL ANNEALIN.Q.

The conventional simulated annealing algorithm uses an adaptation of the Metropolis algorithm to bring a system to equi­ librium at decreasing temperatures. The Metropolis algorithm [2) defines a Markov process that generates a sequence of states, such that the probability of occurrence of any particular state is proportional to its Boltzman weight,

where E( S) is the. energy of state S and (J is the inverse tem­ perature of the system. This ·process generates ·samples from the canonical ensemble; that is, the system is considered to be immersed in a heat bath with a controllable temperature. An­ nealing is accomplished by imposing a schedule for reducing the temperature, 1/ (J, so as to keep the system close to equilibrium.

Creutz has described an alternative technique that simulates the microcanonical ensemble [3). In this method, the total energy remains constant (for some fixed point in the schedule). Instead of simulating a system immersed in a heat bath, the Creutz al­ gorithm simulates a thermally isolated system in which energy is conserved. It performs a random walk through state space, constraining states to a surface of constant energy. The simplest way to accomplish this is to augment the representation with an additional degree of freedom, called a demon, that carries a vari­ able amount of energy, Ev. The total energy of the system is now:

'The performance of the system is not very sensitive to the value chosen for>. . In the example of Section IV, as well as all three examples in [1), >. = 5 was used.

Normally, the demon is constrained to have nonnegative energy, although the possibility of giving it negative energy is useful in annealing, as will be discussed below.

In the Metropolis algorithm a potential new state S' is chosen randomly, and is accepted or rejected based on the change in energy:

If t::,E is negative, the new state is accepted; otherwise, it is accepted with probability exp( -{3t::,E).

The Creutz algorithm is quite similar. If t::,E is negative, the new state is accepted, and the demon energy is increased (Ev <­ Ev - flE). If t::,E is nonnegative, however, acceptance of the new state is contingent upon Ev: if t::,E < Ev the change is accepted, and the demon energy is decreased (Ev <-- Ev- t::,E); otherwise, the new state is rejected. Clearly, the total energy of the system remains constant .

The Creutz algorithm has several advantages. Unlike the Metropolis algorithm, it does not require the evaluation of tran­ scendental functions. Of course, in practice these functions can be stored as tables, but we would like our algorithm to be adapt­ able to fine-grained parallel processors such as the Connection Machine. The small amount of local memory in such machines makes lookup tables unattractive. The Creutz algorithm can eas­ ily be implemented with only integer arithmetic- again, a signif­ icant advantage for fine-grained parallel processors and for VLSI implementation. Experiments indicate that the Creutz method can be programmed to run an order of magnitude faster than the conventional Metropolis method for discrete systems [4]. A further important advantage is that microcanonical simulation does not require high-quality random numbers.

In conventional simulated annealing, we control the process by· specifying the temperature. In the microcanonical version, however, temperature is not a control parameter; it is a statistical feature of the system. In fact, standard arguments can be used to show that at equilibrium the demon energies have a Boltzman distribution:

If the demon energies are non-negative, the inverse temperature can be determined from the mean value:

Control of microcanonical annealing is accomplished by peri­ odically removing.energy from the system. The method used to generate the results in Section IV is as follows: ·

1. Assume that the process begins in a random state So of high energy with respect to the ground state. Call this Eo = E(So) . The initial demon energy is zero.
2. Remove a fixed proportion of this energy, 8E, by reducing the demon energy (Ev <-- Ev - 8E). The results of Sec­ tion IV were obtained with 8E = Eo/300.
3. Run the Creutz algorithm until the system reaches equilib­ rium. A reasonable test for equilibrium is to consider the rate of accepted moves to states of higher energy. If the system is large enough, this rate will increase steadily until it approximates the rate of moves to states of lower energy. We terminate the process for a particular energy level when the rate of accepted moves to higher energy states decreases (measured over N site visits) .
4. Repeat steps (2) and (3) until no further improvement is observed.

This procedure has only one free parameter: the ratio 8E / Ev. While Ev < 0, the Creutz method operates as a "greedy" algo, rithm, accepting only moves to lower energy states and increasing Ev. When Ev becomes positive (which happens quickly if 8E is small), the demon begins to exchange energy between lattice po­ sitions. As the ground state is approached, Ev remains negative because it cannot absorb more energy from the lattice.

The algorithm described above is sequential, and therefore is not suitable for parallel processing. Each local state transition can change Ev , which in turn can affect the next transition. For­ tunately, the technique can be modified to a parallel one by using a separate demon for each lattice site. (As we add demons , the technique moves toward a canonical-ensemble simulation. In fact, if the number of demons is very large compared to the number of sites, the technique specializes to the Metropolis algorithm [3].) Preliminary experiments with one demon per site indicate good performance, although the results of Section IV were generated with the single-demon algorithm.

III HIERARCHICAL ANNEALING

In the original model [1], annealing was performed only at the level of resolution of the stereo images. In some cases, the images were first bandpass-filtered to remove low-frequency com­ ponents- in effect, a simple photometric correction for inconsis­ tent sensor gains or film development. Lower and upper bounds on disparity were specified in advance, and all state transitions were considered with equal probability. As the range of dispar­ ity became large, this scheme required much larger amounts of computation. If we have m permissible disparities and N pix­ els, the size of the state space is mN. If we double the range of disparity, the size of the state space increases by a factor of 2N. Since N is a rather large nurn'ber (2 18 in the example shown in Section IV) , we see that the state space grows explosively with increasing disparity range.

A natural extension of the method is to adopt a hierarchical, coarse-to-fine control structure. At a coarse level of resolution, the number of lattice sites (i.e., the number of pixels) and the range of disparity are small; therefore the size ofthe state space is relatively small. 2 We should be able to compute an approximate ground state quickly, and then use it to initialize the annealing process at the next, finer level of resolution.

· Although the state space may be large in absolute terms.

The Laplacian pyramid, originally developed as a compact image-coding technique [5], offers an efficient representation for hierarchical annealing. In a Laplacian pyramid, an im­ age is transformed into a sequence of bandpass-filtered copies, J0 ,II, I 2 , .. . ,In, each of which is smaller than its predecessor by a factor of 1/2 in linear dimension (a factor of 1/4 in area), with the center frequency of the passband reduced by one octave. This transform can be computed efficiently by recursively applying a small generating kernel to create a Gaussian (low-passed) pyra­ mid, and then differencing successive low-passed images to con­ struct the Laplacian pyramid.

After constructing Laplacian pyramids from the original stereo images, disparity is reduced by a factor of 1/2 in successive lev­ els. Therefore, at some level, disparity is small everywhere. For typical stereo images, we can take this to be level n - 4. [For example, if the original images were a power of 2 in linear di­ mension, the Laplacian images at level n - 4 would be 16 x 16 pixels. Disparities in the range of 0 to 63 pixels in a pair of 512 x 512 images would be reduced to the range of 0 to 1 pixel, with truncation, at the (n- 4)th level.] We shall start annealing at this level, find an approximate ground state, and then expand the solution to the next level. To make this coarse-to-fine strategy work, however, we need two further modifications. We must use a different process for generating state transitions, and we must specify how a low-resolution result is used to start the annealing process at the next-higher level.

In the original model, the probability of choosing a new dis­ parity for consideration as a new state was uniformly distributed over the prior range of disparities:

This is not compatible with our intention of guiding the process with lower-resolution results, however. Because we are now as­ suming the disparity of a lattice site to be close to its correct value, a more effective generating process is to restrict the dis­ parities to increase or decrease by one pixel:

with the further restriction that the disparities are not allowed to specify corresponding points outside the boundary of the right image. In this scheme, the system undergoes Brownian motion through state space. An additional feature of this state-transition function is that it is no longer necessary to specify bounds on disparity in advance.

Expanding a low-resolution result to the next level is slightly more complicated. Obviously, one should begin by simply dou­ bling the size of the low-resolution lattice and doubling the dis­ parity values. Having done this, however, we find that the new state has an artificially low energy because every odd disparity value is "unoccupied," and the new map is therefore more uni­ form than it should be. A spurious symmetry is imposed on the new state that is solely due to the quantization of the previous result, which is likely to place the system near a metastable state (a local minimum) from which it cannot recover. Fortunately, there is an easy solution to this problem: Destroy this symmetry by adding heat. One effective way is as follows:

1. Compute the energy ES of the initial state at level k. (This
2. state has been determined by doubling the result at level k + 1.)
3. Add a fixed proportion of this energy, 6.E", by increasing the demon temperature (Ev +- Ev + 6.E") . The results of Section IV were obtained with 6.E" = ES/10.
4. Run the Creutz algorithm until the system reaches equilib­ rium.
5. Repeat steps (2) and (3) until the number of accepted higher­ energy states exceeds the number of rejected higher-energy states.

IV RESULTS

The images in Figure 1 are an aerial stereo pair (512 x 512 pixels each) covering part of Martin Marietta's test site for the Autonomous Land Vehicle Project near Denver. The terrain is dominated by a long, steep "hogback" (running diagonally across the middle of the images), th.at separates two broad valleys. The highest terrain is in the upper left. Several roads and a few buildings may be seen. Disparity ranges from zero to 50 pixels.

The results of applying the microcanonical, hierarchical an­ nealing algorithm to these data are illustrated in Figure 2, which shows the approximate ground state at the highest level of res­ olution. Disparity values are displayed in 15 pseudocolors with wraparound.

The largest problem solved by the previous version of the sys­ tem was on the order of a 256 x 256 pair with 15 levels of dispar­ ity. Solving this problem required about 12 hours on a Symbolics 3600. The result shown here is a 512 x 512 pair with 50 levels of disparity, and was computed in about the same amount of time. Consider the size of the state space of the two problems. The state space of the current problem 3 exceeds the size of the old problem by a factor of more than 10 36000 ·

V CONCLUSIONS

Two major improvements to a stochastic stereo-matching system have been described: first , a simulation of the micro­ canonical ensemble, as opposed to the canonical ensemble of con­ ventional simulated annealing; and second, the extension to a hi­ erarchical control structure based on Laplacian pyramids. Both

techniques are rather general in nature and can probably be used in other simulated-annealing applications. Together, they permit solutions of stereo-matching problems that are far beyond the competence of the original system.

We are currently evaluating the quantitative increase in per­ formance attributable to each new technique. Qualitatively, the use of microcanonical annealing yields perhaps an order of mag­ nitude increase in efficiency. A potentially more important ben­ efit is that its computation is simpler than the computation of standard annealing, and is therefore more readily implemented in fine-grained parallel systems. The hierarchical control struc­ ture, combined with the Brownian state-transition function, con­ tributes most of the increased performance.

3 Considering the range of disparity to be fixed at 50 pixels.