Group Title: Department of Computer and Information Science and Engineering Technical Reports
Title: A Fast solution to the surface reconstruction problem
Full Citation
Permanent Link:
 Material Information
Title: A Fast solution to the surface reconstruction problem
Series Title: Department of Computer and Information Science and Engineering Technical Reports
Physical Description: Book
Language: English
Creator: Vemuri, Baba C.
Lai, S. H.
Publisher: Department of Computer and Information Sciences, University of Florida
Place of Publication: Gainesville, Fla.
Copyright Date: 1993
 Record Information
Bibliographic ID: UF00095177
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.


This item has the following downloads:

199396 ( PDF )

Full Text

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,


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
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


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
K, = Ko+ U,V, (8)


D -I
-I 0
K0 ..0
Ko =

0 -I D
Us = U1 U2 Uk]
V, = [ V1 V2 .* Vk]

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)

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)

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)


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)
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)

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).


E15 -"


dashed error tolerance 106
sold error tolerance 10

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

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

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)

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

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,
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.


[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.

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs