A Fast Solution to the Surface Reconstruction Problem

B. C. Vemuril and S. H. Lai2

1Department of Computer & Information Sciences

2Department of Electrical Engineering

University of Florida, Gainesville, FL 32611

email:, Iii;'-. II1,a.cis.ufl.edu

Abstract

Surface reconstruction from range data acquired using a vareity of sources has been a very active area of

research in computational vision over the past decade. Generalized splines have emerged as the single most

popular approximation tool to this end. In this paper we present a new and fast algorithm for solving the surface

reconstruction problem using the membrane spline which yields a Co surface. Our algorithm for dense data

constraints takes O(N) computational time, where N is the number of nodes in the discretization. For sparse

data constraints, the algorithm requires O(log N/Tr2) iterations with each iteration taking O(N) time. The number

of iterations in this case depends on a prespecified error tolerance. We demonstrate the algorithm performance

on synthesized sparse range data.

1 Introduction

The problem of reconstructing and representing three-dimensional shapes has received an enormous amount of

attention in vision research for the past decade. 3D shapes may be described either by the bounding surfaces or by

the enclosed volumes. In this paper we will be concerned only with the surface based representations constructed

from depth/range data. Depth or range data may be obtained in numerous ways eg., from stereo, from shading,

from motion, or directly from a laser ranging device. In all these cases with the exception of the last technique, the

data obtained are sparse due to the nature of the solution methods used. For eg., in depth from stereo techniques,

the range data are specified only at a sparse set of locations eg., along object boundaries. The "filling in" may be

achieved by using a surface interpolation technique. There are numerous surface interpolation methods in numerical

analysis [6] and computer vision literature [1]. Among these techniques, the most popular ones in vision have been

those using a variational principle formulation [9, 10, 1]. This formulation leads to what are known as energy-based

surface models which may be used for interpolating sparse data, to smooth noisy depth data, and also to integrate

data from multiple sensors or view points.

We will now briefly present a variational principle formulation of the surface reconstruction problem [9]. The

solution to the general surface reconstruction problem is defined by a configuration of minimal energy E of a physical

system comprised of the thin plate patches, the membrane strips, and the springs. To solve the variational principle

an explicit expression must be obtained for the potential energy functional E : V F- R. V is the admissible space

of the variational principle. The admissible space V consists of functions which represent possible plate-membrane

configurations, a (Sobolev) space of bounded-energy functions. The precise expression for the variational principle

is given by

(v) = p(., Y){r(X, )(v + 2v + V) + [1 r(X, y)](v + v )} dxdy, (1)

where v(x, y) is the admissible deflection function and vx, vy its partial derivatives assumed to be small. p(x, y) and

r(x, y) are called continuity control functions, they constitute an explicit representation of depth and orientation

discontinuities, respectively. The functions range over the interval [0, 1] and need not vary continuously. Setting the

continuity control functions p and r is tantamount to prior knowledge of the location of the discontinuities. Such

knowledge may be obtained from registered multi-sensor data.

The above energy expression serves as a stabilizer in the overall variational principle for the surface reconstruction

problem. To the stabilizer, we add data constraints via what are known as penalty terms. For the surface recon-

struction from depth data, the following penalty term which measures the discrepancy between the surface and data

weighted by the uncertainty in the data may be used,

d(v) = i (v(x, yi) d)2 (2)

Where, di are the depth data points specified in the domain 2 and ci are the uncertainty associated with the data.

The total energy is

C = Ad + CE, (3)

Where A is the regularization parameter that controls the amount of smoothing performed. The goal is to find a u

that minimizes the total potential energy C(v) in equation 3.

To compute a numerical solution to the above minimization problem, we first discretize the functionals E, and Sd

using finite element techniques [9]. For constant p(x, y) and r(x, y) the energy function is a quadratic form. The

energy due to the data compatibility term in discrete form becomes,

Ed(x, d) = d(x d)K( -d). (4)

Where x is the discretized surface, d are the data points, and Kd is a diagonal matrix (for uncorrelated noise in the

data) containing the uncertainty oi associated with the data points. The smoothness energy in discrete form is

E,(x) = xT Asx (5)

where K, is a very large (n2xn2), sparse and banded matrix called the stiffness matrix. The resulting energy function

is a quadratic in x

E(x)= x TKx b+c (6)

with K = AK, + Kd and b = Kdd. The minimum of this energy function u is found by solving the large sparse

linear system

Kx = b (7)

Another way of approaching the solution to the minimization of the energy C is using calculus of variations [9].

The necessary condition for the minimum of the energy functional is a partial differential equation called the Euler-

Lagrange equation. Discretization of this PDE leads to a system of algebraic equations as given in 7. The large

(n2, n2) sparse linear system in equation 7 can be solved using iterative techniques and the best known algorithms

either multi-grid method [8] or hierarchical conjugate gradient [7] take O(NlogN) time where, N = n2. In this

paper, we develop a fast algorithm to solve the above linear system with the plate terms disabled. Our algorithm

takes O(N) time for dense data constraints and O((N log N)/7r2) time for sparse data constraints. Our approach is

based on converting this linear system into a Lyapunov matrix equation and then solving it using an iterative scheme

called the alternating direction implicit (ADI) method. We will describe the details of our technique in a subsequent

section.

The rest of the paper is organized as follows, first we formulate the problem and derive the necessary results

in section 2. Numerical solution technique is then described in section 3. We present some experimental results

depicting the performance of our algorithm in section 4.

2 Problem Formulation

In this section, we will pose and formulate the problem of solving the algebraic system of equations resulting from

the discretization of the Poisson equation or the minimization of the discrete version of the constrained variational

principle discussed in the previous section.

The surface reconstruction problem can be posed as follows: Given a dense or sparse set of depth data, find a

surface which passes close to the data and is smooth.

The stabilizer discussed in the previous section imposes smoothness constraints on the surface reconstruction

problem leading to a Co surface i.e., the membrane spline or a C1 surface namely, the thin plate spline or a

combination thereof called the thin-plate-membrane spline. As mentioned earlier, the numerical solution to the

surface reconstruction problem involves solving a linear system Kx = b. The matrix K is very large and sparse

but is not lightly banded. This is primarily due to the fact that a vector has been used to represent a 2D array

which leads to representing neighboring interactions in 2D by interactions between very distant elements in 1D. No

existing algorithms for such a block banded matrix lead to an O(N) solution. The K matrix is a sum of two matrices

K, and Kd with K, containing the membrane molecules defined on the interior and boundary of f as well as at

discontinuities (see figure 1). The Kd matrix is a diagonal matrix containing non-zero weights. We will now describe

X

Figure 1: Computational molecules, (a) in the interior of f and (b) at the boundary of f

the derivation to transform the above linear system to a Lyapunov matrix equation.

2.1 The Lyapunov Matrix Equation

To transform the linear system Kx = b into a Lyapunov matrix equation, we will first split the matrix K, into two

components.

K, = Ko+ U,V, (8)

where

D -I

-I 0

K0 ..0

Ko =

-I

0 -I D

Us = U1 U2 Uk]

V, = [ V1 V2 .* Vk]

and

4 -1

-1 4 -1

D= ".. E rn.

-1 4 -1

-1 4

Each of the column vectors ui and vi C Rn2" 1 contain only one or two nonzero entry mapped to the boundary of the

discretized f or the discontinuity locations. The matrix Ko is a well-structured matrix, the contents of which were

discussed in the last section. The matrix UV (= iv contains entries indicating the difference between

the matrices K, and K0.

The special matrix Ko can be decomposed as the sum of two tensor products of matrices A and I, i.e.,

Ko =A I + I A, (9)

with

2 -1

-1 2 -1

A = .. .. C lRex,

-1 2 -1

-1 2

and D is the tensor (Kronecker) product.

By using the special structure of the matrix Ko, we can rewrite any linear system of the form Koz = f as the

following Lyapunov matrix equation

AZ + ZA = F, (10)

where Z and F are the original (n x n) matrices corresponding to their concatenated n2 x 1 vectors z and f respectively.

Therefore, instead of solving the original large linear system Kx = b in equation 7, our approach is to solve the

corresponding Lyapunov matrix equation with a modified right-hand side matrix B. Notice that the matrix A in

equation 10 is very tightly banded (tridiagonal). In addition to being tridiagonal, the matrix A is also SPD and

Toeplitz. We can use the ADI (Alternating Direction Implicit) method to solve this Lyapunov matrix equation in

O(N) operations for dense data and in O(NloglN/72) for the sparse data constraints. We will describe he ADI

method in the next section.

2.2 Modification of the Right-hand Side

Our method is based on the ADI technique tailored for the Lyapunov matrix equation with a tridiagonal, SPD and

Toeplitz matrix A. In the following, we are going to show how the original Kx = b problem can be converted to an

equivalent linear system form Kox = b' by modifying its right-hand side. The solution of Kox = b' may be obtained

by solving the corresponding Lyapunov matrix equation containing the modified right-hand side.

The derivation is mainly based on the Sherman-Morrison-Woodbury formula [3]

(A + UVT)-1 = A-1 A- U(I + VA-U)-VTA-1, (11)

where A E nxn, U E Rnxm V E Rmxn, and I E l xm is the identity matrix.

By decomposing every nonzero entry of the diagonal matrix Kd into the outer product of u$ and v' and putting

them into U and V as follows:

U = [ ... Uk U ...Uk

V = Avi ... Avk Vj ... vk

and using equations (3) and (4), we have

K = AK + AUV + Kd, (12)

= AKo + UVT. (13)

By substituting the Sherman-Morrison-Woodbury formula into x = K-lb, we get the following result

x = Ko,(b Uy)/A, (14)

with

y = (I+ A- VTKolU)- (A -VTKKob). (15)

Therefore, the right-hand side of Kox = b' is given by b' = (b Uy)/A. The vector y E R(k+k')x1 participates in

modifying the values of the b vector at locations corresponding to the data constraints, discontinuities and boundary

of f. Computation of vector y is nontrivial and will be the focus of discussion in a subsequent section.

3 Numerical Solution

Our numerical algorithm can be described as follows:

1. Compute c := A-1VTKo b.

2. Form the matrix E := I + A -VTKoU.

3. Solve Ey = c and modify the right-hand side to b' = (b Uy)/A.

4. Solve AX + XA = B' where B' is the matrix form of b'.

In steps (1) and (4), the ADI method is used to solve the Lyapunov matrix equations with different right-hand side B

(the matrix form of vector b) and B'. After computing xo = Kolb using the ADI method, the vector c is obtained

by multiplying xo with VT whose rows serve as sampling vectors as each of them contain only single nonzero entries.

The linear system Ey = c is nonsymmetric and needs to be solved efficiently. In this paper, we use the biconjugate

gradient (BCG) technique [] in a hierarchical basis [7], to solve this problem. In the following, the details of the ADI

method for the Lyapunov equation, the formation of matrix E, and the solution to the nonsymmetric linear system

Ey = c will be given.

3.1 ADI Method for the Lyapunov Equation (Stages (1) and (4) of the algorithm)

The ADI method for solving a Lyapunov matrix equation is described in [5]. For any Lyapunov matrix equation of

the form AX + XA = B, the ADI method involves the following steps in each iteration j = 1, 2,..., J,

(A+p3I)XN_ = B-Xji(A-pjI) (16)

Xj(A+pjI) = B-(A-pjI)Xj_ (17)

The nice properties of matrix A help to advance the ADI iterations very efficiently in each step of the technique.

Since A is SPD and tridiagonal, we can use Cholesky decomposition for the tridiagonal SPD matrix (A + p I) to

compute the updated solution. This solution update requires only O(N) time per iteration.

We now examine the convergence of the ADI method for the problems in stages (1) and (4) of our algorithm. Let

the initial X = 0 (zero matrix), and AXt = Xt-X*, where X* is the true solution of the Lyapunov matrix equation,

A,, A2, . ., and qi, q2, qn be the eigenvalues and the eigenvectors of A respectively. For our particular matrix

A, the left and right eigenvalues and eigenvectors are the same, and the eigenvectors are orthogonal, i.e.,

A = QDQ, (18)

where

Q = [qi q2 .. qn],

D = diag[A1, A2,.*, An]

sin k7rl

_2 sin 2-

qk + 1

k fnk Jet

and Ak = 2 2 cos k- for k = 1, . ., n. Notice that Q is the matrix corresponding to the discrete sine transform

and Q = QT = Q-1

By taking the tensor product of the eigen vectors, qk ql, k, I = 1, . ., n, we can generate an orthogonal basis for

Rnxn and express the true solution X* in this basis as follows:

X* = ak(qk qj). (19)

k=l 1=1

Subtracting each of the equations 16 and 17 from the Lyapunov equation AX + XA = B, and expressing X in the

above tensor product basis, we have the following result for the error after t iterations,

AX, = [n tA-3 )A pj akl(qk q(). (20)

k=1 l=1 j=

The error AXt is bounded by

|AXt < p X*l, (21)

where

ti (Ak pj )(A PJ- pj

Pt = max I (-)

k,l j=1 Ak + Pj A, + Pj

Obviously, the error bound depends on the choice of the ADI parameters pj. Choosing the optimal set of ADI

parameters pj's leads to the ADI minimax parameter problem. Given the bounds for the eigenvalues of the SPD

matrix A as, Vk, 0 < a < Ak < b, we have,

Pt < min max fl( )2. (22)

l,...,t O

The right-hand side is the minimax ADI parameter problem. The solution to this minimax problem gives the

optimum set of ADI parameters. Computing the optimal choice of the ADI parameters is described in [5].

The minimax analysis also provides the solution to the number of iterations needed in ADI to achieve a specified

error tolerance c in X such that I|AXj| < c lX, || This number of iterations is given by

ln4In T4

J = [ k (23)

72

where

1

k1 = I

m + V/m_'2

m = -(K2(A) + ),

22A K2(A)

[z] denotes the smallest integer larger than z, and K2 is the ratio of the largest and smallest eigenvalues, i.e.,the

2-condition number of A.

In the dense and sparse data cases, we have different A matrices in the Lyapunov matrix equations. In the

following, we give the analysis of numbers of iterations needed for each of these different cases.

3.1.1 Dense Data Case

The data compatibility matrix Kd for dense data case is the identity matrix. We can add this identity matrix to the

matrix Ko to form the Lyapunov matrix equation with different matrix A i.e., it is the original A plus the scaled

identity matrix -I. This doesn't change the eigenvectors of A, but its eigenvalues are shifted by Therefore, all

the eigenvalues of A are bounded between and 4 + independent of the size of the matrix A. From equation

23, we can see the number of iterations needed will be bounded by a constant number for a given error ratio. Hence,

the number of iterations needed in ADI is independent of the size of the problem for the dense data case. Thus, the

computational time for the dense data case is O(N).

20-

E15 -"

10-

dashed error tolerance 106

sold error tolerance 10

... . ........ .... .... .... .... .... .. . 00

n

Figure 2: The number of iterations needed in the ADI method with respect to the size n of the matrix A E Rnxn

3.1.2 Sparse Data Case

If we put the largest and smallest eigenvalues, 2 2 cos n- and 2 2 cos-, into equation 23, we can get the

following approximation for the number of iterations,

4 16 log n

J In (log 2 + ). (24)

C 7 2o7 2/2

In figure 2, we plot the number of iterations as a function of the input matrix size (nx n). We can see that the number

of iterations required is very small for practical ranges of matrix sizes typically encountered in vision problems.

The as is application of minmax analysis to the ADI method for our problem leads to a very loose upper bound

on the number of iterations. This is due to the fact that no information about the spectrum of the membrane

interpolant is incorporated into the analysis. However, we do have the a priori knowledge about the spectrum of

the exact solution, i.e. the coefficients akl's in equation 19. Taking this information into account leads to a different

minimax problem which will yield a different choice of ADI parameters. Thus, we can get a tighter bound for the

required number of iteration than that given in equation 23. We leave this improvement as the focus of our future

research.

3.2 Formation of the Matrix E (Stage (2))

Each entry (i, j) in matrix E can be written as

E(i,) = A- TKo- u + (i j). (25)

Recall that every vector vi and uj contain only one or two nonzero elements. This means that vfK- lu is the

sampling at the locations vi of the impulse response of the linear system Ko 1 or the Lyapunov system with the impulse

at the locations uj. Unfortunately, this system is not shift-invariant and this makes the problem of efficiently forming

the matrix E more difficult.

Assume the nonzero elements of vi and uj are located at (p,, qi) and (pj, qj), respectively. We can find an analytic

form of the solution at (p,, qi) to the Lyapunov matrix equation AX + XA = B, with B(k, 1) = 6(k pj, qj) by

substituting equation 18 for A and solving for X. This leads to,

X(p, q) = fn(pi pj, qi j) + fn (pi + pj, i + qj) fn (pi pj, qi + qj) fn(pi + pj, qi qj), (26)

where

fn n(p, Cos kpr Cos lq

q 8 2 cos k cos 1l2

Sk=l 1=1 con+s 1 cosn+1

The function fn(p, q), for 0 < p, q < n + 1, should be precomputed and stored as an ((n + 2) x (n + 2)) lookup table.

Then each entry in the matrix E can be efficiently calculated from the above equation and the look-up table.

3.3 Solution to Ey = c (Stage 3)

The linear system Ey = c is nonsymmetric and indefinite. For the dense data case, the size of the matrix E is k x k,

where k is the number of discontinuities. For the sparse data case, The size of the matrix E is (k+k') x (k+k'), where

k is the number of boundary points and discontinuity locations and k' is the number of data points. A reasonable

assumption is (k + k') O(n). Consequently, it is important to have a fast algorithm to compute the solution y to

this nonsymmetric linear system.

The conjugate-gradient(CG) method is a very powerful iterative scheme to solve the symmetric positive definite

linear system. A natural generalization of the CG-type method to solve the general nonsymmetric linear system is

called the biconjugate gradient(BCG) method. The BCG algorithm [4] is as follows.

0. Choose yo, and set qo = ro = c Eyo; Choose ro : 0, and set qo = io;

1. Compute a k = r_ rk- ,f = qkEqk-1, and ak = O k/af;

2. Update Yk = Yk-1 + akqk-i;

3. Set rk = rk-1 kEqk-1, and rk = ik-1 akET'qk-1;

4. Compute fN = rk, and Ak = 3Na;

5. Set qk =rk + kqk-1, and qk = rk + kk-l;

6. If rk w 0 or rk w 0, stop; else go to 1.

Like CG, The operations and storage requirement per step in BCG is constant. The convergence behavior is similar to

the CG method except that the BCG algorithm is susceptible to breakdowns and numerical instabilities. To be more

specific, division by 0 (or a number very close to 0) may occur in computing ak and Ak. These two different breakdowns

may be avoided by using the look-ahead Lanczos algorithm or the quasi-minimal residual(QMR) algorithm [2].

As in the CG-algorithm, preconditioning techniques can speed up the convergence of the BCG algorithm. In this

paper, we use the technique similar to the hierarchical conjugate gradient(HCG) method, which involves smoothing

the residual vector in each step. However, in BCG, we have to smooth two residual vectors rk and rk in each step.

The smoothing operation is the same as that in HCG. We call this preconditioned BCG method the hierarchical

biconjugate gradient(HBCG) method. In the BCG algorithm, each iteration requires that Eqk-1 and ETqk 1 be

computed. Since the matrix E is a dense matrix, the computational cost for the direct matrix-vector products in

each iteration is may be too much. After examining the matrix E, we find that most of its entries away from the

diagonal are very close to 0 and the overall structure of E is very smooth. Therefore, we can approximate the matrix

and vector with very low-order polynomials and use the polynomial integration to approximate the inner product in

the off-diagonal insignificant region. This leads to a significantly fast scheme in forming the required matrix vector

products.

4 Experimental Result

present one example of the algorithm performance on a 64 64node surface reconstruction problem. The input

40 60 404 60

20 20

0 0 0 0

(c) (b)

data set is very sparse and contains 15 data points randomly distributed in the plane as shown in figure 3.3. The

discontinuity locations are specified along the line between coordinates (1,32) and (30,32).

6 ,0,

40

70,

20

o50

60 60

40 40 60

40 40

20120 20

() (d)

Figure fir Surface reconst find ruction example; (a) original sparse data i.e., RHS of Kx = b, (b) modifipresented in section 3. The

i.e., RHS of Kox = b', (c) solution after 1 iteration of ADI method and (d) after 16 iterations.

Our algorithm was tested on several data sets (both dense and sparse) of varying input sizes. In this section, we

present one example of the algorithm performance on a 64 x 64-node surface reconstruction problem. The input

data set is very sparse and contains 15 data points randomly distributed in the plane as shown in figure 3.3. The

discontinuity locations are specified along the line between coordinates (1,32) and (30,32).

The first step is to find the modified right-hand side by using the numerical method presented in section 3. The

modified right-hand side B' is shown in figure 3.3. This depicts the changes that are made to the original data before

applying the ADI method. Note the changes along the boundary, and to the data values within the boundary. The

ADI iterations are then employed to obtain the interpolated surface. From equation 23, it is evident that just 16

iterations can guarantee the error tolerance c to be within 10-6. We depict the surfaces after 1 and 16 ADI iterations

in figure 3.3.

In figure 3.3, we can see that just one ADI iteration can produce the global shape of the exact surface. This is

because we advance the ADI iterations with the parameters pj's starting with a small value and proceeding toward

large values. The smaller parameters are used to recover the low frequency detail of the surface and the large values

refine the shape with high frequency components. Consequently, it is not surprising that the global shape was

recovered in one iteration.

This phenomena of global to local shape recovery is quite different from the other iterative numerical methods.

Most of the other numerical methods do not exhibit this type of spectral characteristic during the iterations. In

addition to the fast convergence rate, another advantage of the ADI method is that we can make use of the spectral

properties of the prior surface model (in our case the membrane spline) to speed up the convergence rate. This can

be achieved via adaptation of the ADI parameters to the spectrum of this model constrained by the given data.

5 Conclusion

In this paper, we presented a new and fast algorithm for solving the surface reconstruction problem using a membrane

spline interpolant. The problem is formulated as the solution to a Lyapunov matrix equation with an appropriate

right-hand side. We derive this right-hand side and present fast numerical methods to compute it. To solve the

Lyapunov matrix equation, we use the ADI algorithm and obtain convergence (for a prespecified tolerance) in very

few iterations for different sizes of the A matrix. In addition, we can speed up the convergence rate of the ADI

method by making use of the observation that the ADI parameters are related to different amounts of error reduction

in different spectral ranges. We will address this issue in our future research.

References

[1] R.M. Bolle and B.C. Vemuri. "On three-dimensional surface reconstruction methods,". IEEE Trans. Pattern

Anal. Machine Intel., 13(1):1-13, 1991.

[2] R.W. Freund and N.M. Nachtigal. "QMR: a quasi-minimal residual method for non-Hermitian matrices,".

Technical Report 90.51, RIACS, NASA Ames Research Center, Moffett Field, 1990.

[3] G.H. Golub and C.F. Van Loan. Matrix Computations. The Johns Hopkins University Press, 2nd edition, 1989.

[4] C. Lanczos. "An iteration method for the solution of the eigenvalue problem of linear differential and integral

operators,". J. Res. Nail. Bur. Stand., 4'. ..-282, 1950.

[5] A. Lu and E. L. Wachspress. "Solution of Lyapunov equations by alternating direction implicit iteration,"

Computers Math. Applic., 21(9):43-58, 1991.

[6] L.L. Schumaker. "Fitting surfaces to scattered data,". In G.G. Lorentz, C.K. Chui, and L.L. Schumaker, editors,

Approximation Theory II, pages 203-'',.;. New York: Academic, 1976.

[7] R. Szeliski. "Fast surface interpolation using hierarchical basis functions,". IEEE Trans. Pattern Anal. Machine

Intell., 12(6):513-528, 1990.

[8] D. Terzopoulos. II ....i analysis using multigrid relaxation methods,". IEEE Trans. Pattern Anal. Machine

Intell., 8:129-139, 1986.

[9] D. Terzopoulos. "The computation of visible-surface representations,". IEEE Trans. Pattern Anal. Machine

Intell., 10(4):417-438, 1988.

[10] B.C. Vemuri, A. Mitiche, and J.K. Aggarwal. "Curvature-based representation of objects from range data,".

Image and Vision Comput., 4(2):107-114, 1986.