Some Experiments in Applying Inductive Inference Principles to Surface Reconstruction

1. Introduction

In its simplest form, the Minimum Description-L<?ngth (MDL) principle states that induction can be accomph~hed by finding the shortest description for a ~t c_>f o~servattons in a suitable language, where the descnptton mcludes a theory that accounts for the o~servations. It has been ~ro~ed mathematically that choosmg the shortest descnptiOn enables one to eventually converge on an "appropriate" theory, given a sufficient number of observations, wh~re appropriateness is judged in terms of the ability to predict observations [Pednault 1988; Barron 1985; Barron and Cover 1983· Rissanen 1978]. This convergence property gives the MDL principle tremendous potential for application to vision and other induction problems, since it guarantees robust inferences for sufficient amounts of data, even in the presence of noise.

The MDL principle, though, is computationally intractable in its most general form. Its convergence properties have been studied with respect to languages that are Turing-equivalent and thus allow any co'!lp~tab~e function to be represented. To actually apply the pnnctpl~ m this general context, one would have to solve the haltmg problem of Turing machines, which is clearly impossible.

However, as I have argued elsewhere [Ped~ult.1984]~ the halting problem does not doom the MDL pnnctple, smce one can limit the range of theories to a tractable subset, and/or employ approximation techniques that ~ttempt to fi~d theories that are as close to the optimum as IS computationally feasible. This pragmatic approach enables one to construct efficient algorithms by tailoring the optimization to specific properties of. interest. This l~ad~ to the following methodology for applymg the MDL pnnctple to computer vision and other induction problems:

considered in the examples. Accordin~ to the principle, the theory that best accounts for a collection of observations is the one that yields the shortest description [Pednault 1988; Barron 1985; Barron and Cover 1983; Rissanen 1978]. A description consists of a machine-readable representation of the theory plus an encoding of the observations with respect to the theory. When applied to surface reconstruction the "theory" is the function that defines the surface and the "observations" are the data points. These points are encoded in terms of the difference between the points and the reconstructed surface. The optimal reconstruction is then the function S that minimizes the sum

where l(S) is the length in bits of a machine-readable de~ription of S, z1, .... ~n are the data points, and /(zt.·· .. zn I S) ts the number of bits needed to encode the difference between these points and S. The combination of a function and an encoding of the difference constitutes an exact representation of the data points. The optimal surface is the one that minimizes the length of this representation.

In general, S will be a piecewise-continuous function that .can be deco!Dposed. into a collection of regions and c~ntmuou.s functions, with one function per region. l(S) wtll thus mc~rporate the number of bits needed to specify both the reg10n boundaries and the functions within each regi?~· If these fun~tions form a collection of parametric famihes, the function associated with a region can be described; br f~st specifying the family it belongs to and then specifymg its parameters. The number of bits needed to supply thi~ information is included in l(S). Note that l(S) wdl thus mcrease as the number of regions and function parameters increase. We can therefore view l(S) as measuring the complexity of the reconstructed surface.

To encode the data points, the difference between their values and S is analyzed statistically by viewing the difference as a random process. This random process, together with the function S, induces a probability dis~bution p on the data points. If we employ the Shannon coding technique developed in information theory [Gallager 1968]. we can use p to encode the points using no more than

bits. l(zt ..... zn IS) can thus be thought of as a degree-of-fit term. If the S reproduces the data points exactly, p(z1, ··· ,zn) will equal one and /(z1, .. · ,zn I S) will be zero. As the fit ?egrades, p(zt.····zn) will decrease and l(z1, ... ,zn IS) will mcrease.

Sine~ th~ de &Tee of fit can vary from region to region, we are not JUStified m assuming that the probability distribution p remams constant. To allow it to vary, the distribution m!JS~ be specif!ed alo.ng with the interpolating function wtthm each reg10n. This can be accomplished by treating p as a member of a parametric family of distributions and then encoding its parameters. The number of bits needed to specify this information is included in l(S).

The optimum surface is the one that achieves the minimal total coding length over all partitionings of the points into regions and all possible assignments of surfaces and noise models to each region.

3. Polynomial Surface Reconstruction

For t~e experiments reported in this paper, piecewise­ polynomtal functions were used in conjunction with a Gaussian noise model to describe surfaces and the random displacement of the data points. The algorithms developed actually allow any linear combination of basis functions. Polynomials were used in the experiments because of their familiarity in order to facilitate the presentation of the results. A_ Gaussian noise model was selected, since optimal polynom~al coefficie.nts can then be determined using a computationally-efficient least-squares algorithm. Nonlinear families of interpolating functions can also be used within the framework, as can other noise models but different techniques must then be employed to fmd optimal function parameters (the least-squares algorithm presumes a Gaussian noise model and linear combination of basis functions) .

.(\.polynomial-time algorithm has been developed to fmd optimal reconstructions for one dimensional surfaces (i.e., curves). Finding an optimal reconstruction is computationally infeasible for a two dimensional surface since an exponential number of regions must be examined i~ order to find the optimum. For a multidimensional surface, one has no choi_ce but to use approximation techniques to find reco!lstructions !hat are as close to the optimum as is ~mputation~ly feasib~e. One approach I am investigating lS to use a ~enes of optimal 1D reconstructions to guide the reconstruction of 2D surfaces. Early experiments with this appro~imation technique have proved quite promising.

With lD surface reconstruction, regions become intervals on a line. Polynomials are used to describe the surface (curv~) withi!l each regi~n. To obtain an efficient algorithm, each mtervalis encoded mdependently and the resulting codes are then co.ncatenated to form the encoding of the overall surface. This allows the total coding length to be minimized using dy~am~c programming techniques, which produces a polynomial-time algorithm.

The data poin~ are assumed to take on integer values from 0 to 255. ThlS corresponds to range data quantized to 8 bits accuracy. Because of the discrete nature of the data the noise model employed is a bounded discrete �aus~ian distribution having the probability mass function

where z represents the value of a data point, and.

The parameter ~ represents the true height of the surface at a point, while y determines the distribution of errors when measuring o. Note that o and y2 only loosely correspond to the m~an and variance of the distribution. They are not numencally equal to these statistics because the distribution is bounded and discrete.

The underlying surface is assumed to be sampled at uniformly spaced increments so that z; represents the value of the i'th sample. The polynomials are then expressed as functions of i. The measurement errors for the samples within an interpolation interval are assumed to be statistically independent. The joint distribution is then the product of the marginal distributions for the points in the interval. The parameter y2 is a assumed to be constant over an interyal, while o is given by the interpolating poly~omial. . Thus, . if a0+ ··· +a,.im is the polynomial ass~~ated. Wi~ tJ:Ie mterval containing points i 1 through it> the JOmt distnbuuon for these points is given by

The number of bits needed to encode points i1 through h relative to the polynomial is given by the negative log of this expression:

S; 1 h represents an encoding of the polynomial function, the n~mber of points in the interval,. and the noise p~eter y . An upper bound on l(z;l' ... ,z)l I S;y~) can be obtamed by replacing the~(· .. ,y) term by Lmins p~5,y)]. Minimizing this upper bound with respect to the polynomial coefficients reduces to a least-squares fit problem and produces near­ optimal coefficients. Minimization with respect to y2 can be accomplished via a Newton-Raphson method.

The total coding length is the sum of the coding lengths for each of the intervals. Thus, if i" ... ,i1 are the starting points of the intervals and h .... Jt are the ending points, where i 1 = 1,it = n, and ik+t = 1 + h for 1 ~ k < n, then the total coding length is given by

The number of intervals and their endpoints are determined by minimizing this sum using dynamic programming techniques.

The quantity l(Si ·J is equal to the number of bits needed to encode the wid:K of the interval (i.e., h- ik + 1), the order of the polynomial, the polynomial coefficients, and the noise parameter y2. The interval width and the order of the polynomial are integers and are encoded using the Rissanen­ Cover-Elias universal prior for integers [Rissanen 1983; Leung-Yan-Cheong and Cover 1978; Elias 1975]. If m is the order of the polynomial, the number of bits needed to encode m and the interval width is given by log; (m + 1) + c and log;(j~- ik + 1) + c , respectively, where c = 2.865 and where log2 (x) = log2(x) + log2log2(x) + ... up to but not including the first negative term.

Whereas the interval width and the order of the polynomial are integers, y2 and the polynomial coefficients are real-valued and must therefore be quantized. Rissanen [1983] has proposed a quantization scheme that would minimize the worst-case coding length assuming uniform quantization steps. However, this scheme was found to produce somewhat inefficient codes in the experiments reported here, mainly because the data points are already quantized to 8 bits. A nonuniform quantization scheme was therefore developed instead to achieve more compact encodings.

4. Experimental Results

Figures 1 and 2 illustrate the performance of the reconstruction algorithm using synthetic data. In these figures, Gaussian white noise of various standard deviations is added to two different piecewise polynomial curves. The noisy curves and their resulting reconstructions are plotted side by side for comparison. The underlying curve in each figure consists of three polynomial regions of orders 0, 1, and 2, respectively. The curve in Figure 1 is continuous, while in Figure 2 tt is discontinuous. In Figures la and 2a, the curves are sampled at 48 points and the standard deviation of the noise is increased from 0 to 32 in a logarithmic fashion. In Figures 1 b and 2b, the standard deviation of the noise is held constant at 32 and the number of data points is increased from 48 to 768 by factors of two in Figure 1 b, and from 48 to 384 in Figure 2b.

Figure 3 presents the results of applying the algorithm to the rows of a 128 by 128 range image of a mechanical part. Figure 3a shows the original image, while Figure 3b shows the reconstruction. In Figure 3c, a row and a column from the original image are plotted along side their reconstructions. The positions of the row and column are indicated by horizontal and vertical marks along the borders of Figures 3a and 3b.

In all of these examples, the number of intervals, their locations, and the order of the polynomials were determined entirely by minimizing the total encoding length of the data points. The initial implementation allowed up to 15th order polynomials. Problems with numerical instability, however, required that the polynomials be limited to at most 5th order in the examples presented. Further work is being done to resolve these stability problems.

S. Discussion

The experimental results demonstrate that the convergence properties of the MDL principle can be preserved in practical applications. Figures lb and 2b clearly show that an exact surface reconstruction is obtained asymptotically when the true surface is a member of the subset considered (i.e., polynomials) and one is minimizing the description length over this subset. In both cases, the structures of the underlying curves are recovered at the highest data samplings. The use of a high noise level illustrates that convergence occurs independent of the amount of noise. The noise level does affect the rate of convergence (i.e., the higher the noise, the slower the convergence), but not the eventual convergence. Note that convergence is much faster for the discontinuous curve.

Figures 1 and 2 together demonstrate that the degree to which the reconstruction approximates the true surface adjusts automatically to the number of data points and to the amount of error in the measurements when the number of points is insufficient for exact reconstruction. Coarser approximations are made when the measurement errors are large and/or the number of points is small, and finer approximations are made when the errors are small and/or the number of points large. The result is a graceful decrease in the accuracy of the reconstruction as measurement error increases or the number of data points decreases. This compensation occurs automatically without the manual adjustment of parameters. Note that the decrease in accuracy is more gradual for the discontinuous curve.

The mechanical part shown in Figure 3 does not have polynomial surfaces, yet the reconstruction produced by the algorithm is quite accurate. Thus, even when the true surface lies outside of the subset considered, a good approximation to the surface can still be obtained. The principal requirement is for the subset to have sufficient latitude for an adequate approximation. In addition, the noise model must adequately reflect the statistical properties of the measurement errors, which does happen to be the case for the mechanical part.

The surfaces in Figure 3 are actually two dimensional, not one dimensional as is assumed by the algorithm. However, finding a minimal encoding of a 2D surface is computationally infeasible, since an exponential number of regions would have be examined to find the optimum. I am investigating an approach to 2D reconstruction that utilizes a series of 1D reconstructions determined through exact minimal encodings. The intuition is that high-quality 1D reconstructions based on polynomials of sufficient degree will provide tight constraints on the 2D surface. The quality of the reconstructions obtained on the mechanical part strongly supports the viability of this approach. It also provides evidence that one need not necessarily find the absolute shortest description in order to obtain a good approximation to the true surface-a description that is short enough can do just as well.

6. Conclusions

The strength of the MDL principle lies in its convergence properties [Pednault 1988; Barron 1985; Barron and Cover 1983; Rissanen 1978]. The benefits of these properties provide strong motivation for discovering ways of overcoming the computational barriers posed by the principle in its most general form. The results presented here demonstrate that these barriers can be overcome in specific applications by limiting the range of admissible theories to a tractable subset and/or by employing approximation techniques that attempt to find theories that are as close to the optimum as is computationally feasible.

The approach that I am exploring in applying the MDL principle to surface reconstruction illustrates one way of dealing with these computational issues. Different means are also being explored by other researchers in computer vision. Segen [1980, 1985, 1988] has applied the principle to several problems in vision and perception using various approximation techniques. Smith and Wolf [1984] have considered an approach closely related to the one presented here in applying the MDL principle to curve matching for image correspondence. Leclerc [1988] is investigating a distributed 2D approximation to the MDL principle for the purpose of image partitioning. Pentland [1988] has developed an algorithm motivated by the MDL principle for the reconstruction of su~rquadrics from binary images. Though different approaches to constructing approximation algorithms are taken in each instance, they all fall within the basic paradigm outlined above and in Section 1.

This previous work focused primarily on the practical application of the MDL principle without considering its convergence properties. The results presented in this paper demonstrate that convergence is one of the key aspects of the MDL principle. Since convergence is not guaranteed when approximation techniques are employed, it should be included as one of the tests for inappropriate behavior when following the methodology described in Section 1. The variety of applications previously considered, together with the convergence properties examined here, demonstrate the flexibility and viability of the MDL principle. Undoubtedly, the future will see many more applications of this principle.

Acknowledgements

I wish to thank Steve Zucker for introducing me to computer vision years ago, Tom Cover for having introduced me to Kolmogorov complexity and the minimal description-length principle, and Steve, Tom, Peter Cheeseman, Mike Georgeff, Yvan Leclerc, John Mohammed, Jakob Segen, and many others for the stimulating and challenging (if not heated) discussions on these topics over the years. John Gabbe and Larry O'Gorman provided helpful comments on earlier drafts. Again, I thank my faithful editors Corinne and Marla Babcock. Gerard Medioni at the University of Southern California kindly provided the range image that appears in Figure 3.