IEEE TVCG 9(1) Jan 2003

Computing and Rendering Point Set Surfaces
Marc Alexa

Johannes Behr

Daniel Cohen-Or

Shachar Fleishman

David Levin

Claudio T. Silva

TU Darmstadt

ZGDV Darmstadt

Tel Aviv University

Tel Aviv University

Tel Aviv University

AT&T Labs

Abstract
We advocate the use of point sets to represent shapes. We provide a definition of a smooth manifold surface from a set of points
close to the original surface. The definition is based on local maps
from differential geometry, which are approximated by the methodof moving least squares (MLS). The computation of points on thesurface is local, which results in an out-of-core technique that canhandle any point set.
We show that the approximation error is bounded and present
tools to increase or decrease the density of the points, thus, allowing
an adjustment of the spacing among the points to control the error.
To display the point set surface, we introduce a novel point rendering technique. The idea is to evaluate the local maps according
to the image resolution. This results in high quality shading effects
and smooth silhouettes at interactive frame rates.
CR Categories:
I.3.3 [Computer Graphics]: Picture/Image
Generation—Display algorithms; I.3.5 [Computer Graphics]:
Computational Geometry and Object Modeling—Curve, surface,solid, and object representations;
Keywords: surface representation and reconstruction, moving
least squares, point sample rendering, 3D acquisition

1

Introduction

Point sets are receiving a growing amount of attention as a representation of models in computer graphics. One reason for this is the
emergence of affordable and accurate scanning devices generatinga dense point set, which is an initial representation of the physicalmodel [34]. Another reason is that highly detailed surfaces requirea large number of small primitives, which contribute to less thana pixel when displayed, so that points become an effective display
primitive [41, 46].
A point-based representation should be as small as possible
while conveying the shape, in the sense that the point set is neither
noisy nor redundant. In [1] we have presented tools to adjust thedensity of points so that a smooth surface can be well-reconstructed.
Figure 1 shows a point set with varying density. Our approach is
motivated by differential geometry and aims at minimizing the geometric error of the approximation. This is done by locally approximating the surface with polynomials using moving least squares(MLS). Here, we include proofs and explanations of the underly-

Figure 1: A point set representing a statue of an angel. The densityof points and, thus, the accuracy of the shape representation arechanging (intentionally) along the vertical direction.

ing mathematics, detail a more robust way to compute points on thesurface, and derive bounds on the approximation error.
We understand the generation of points on the surface of a shape
as a sampling process. The number of points is adjusted by either
up-sampling or down-sampling the representation. Given a data setof points P = {pi } (possibly acquired by a 3D scanning device),
we define a smooth surface SP (MLS surface) based on the input
points (the definition of the surface is given in Section 3). We suggest replacing the points P defining SP with a reduced set R = {ri }defining an MLS surface SR which approximates SP . This generalparadigm is illustrated in 2D in Figure 2: Points P, depicted in purple, define a curve SP (also in purple). SP is re-sampled with points
ri ∈ SP (red points). This typically lighter point set called the representation points now defines the red curve SR which approximates

ri

pi
SP
a)

and its interaction with the real-world object being acquired. Fora non-trivial object, it is necessary to perform multiple scans, eachin its own coordinate system, and to register the scans [6]. In general, areas of the objects are likely to be covered by several samples
from scans performed from different positions. One can think of
the output of the registration as a thick point set.
A common approach is to generate a triangulated surface modelover the thick point set. There are several efficient triangulation
techniques, such as [2, 3, 5, 7, 17]. One of the shortcomings of thisapproach is that the triangulated model is likely to be rough, containing bumps and other kinds of undesirable features, such as holes
and tunnels, and be non-manifold. Further processing of the triangulated models, such as smoothing [51, 14] or manifold conversion [20], becomes necessary. The prominent difficulty is that the
point set might not actually interpolate a smooth surface. We callconsolidation the process of “massaging” the point set into a surface. Some techniques, such as Hoppe et al [23], Curless and Levoy[12], and Wheeler et al [55] consolidate their sampled data by using an implicit representation based on a distance function definedon a volumetric grid. In [23], the distances are taken as the signed
distance to a locally defined tangent plan. This technique needs
further processing [24, 22] to generate a smooth surface. Curlessand Levoy [12] use the structure of the range scans and essentiallyscan convert each range surface into the volume, properly weighting the multiply scanned areas. Their technique is robust to noise
and is able to take relative confidence of the samples into account.
The work of Wheeler et al [55] computes the signed distance to aconsensus surface defined by weighted averaging of the different
scans. One of the nice properties of the volumetric approach is that
it is possible to prove under certain conditions that the output is aleast-square fit of the input points (see [12]).
The volumetric sign-distance techniques described above are related to a new field in computer graphics called Volume Graphics
pioneered by Kaufman and colleagues [26, 54, 50] which aims toaccurately define how to deal with volumetric data directly, and answer questions related to the proper way to convert between surface
and volume representations.
It is also possible to consolidate the point set by performing
weighted averaging directly on the data points. In [53], model triangulation is performed first, then averaging is performed in areas
which overlap. In [49], the data points are first averaged, weighted
by a confidence in each measurement, and then triangulated.
Another approach to define surfaces from the data points is toperform some type of surface fitting [18], such as fitting a polynomial [31] or an algebraic surface [42] to the data. In general, it
is necessary to know the intrinsic topology of the data and (sometimes) have a parametrization before surface fitting can be applied.
Since this is a non trivial task, Krishnamurthy and Levoy [28] haveproposed a semi-automatic technique for fitting smooth surfaces todense polygon meshes created by Curless and Levoy [12]. Anotherform of surface fitting algorithms couples some form of high-level
model recognition with a fitting process [44].
The process of sampling (or resampling) surfaces has been studied in different settings. For instance, surface simplification algorithms [9] sample surfaces in different ways to optimize rendering
performance. Related to our work are algorithms which use particle
systems for sampling surfaces. Turk [52] proposes a technique forcomputing level of details of triangular surfaces by first randomly
spreading points on a triangular surface, then optimizing their positions by letting each point repel their neighbors. He uses an approximation of surface curvature to weight the number of points
which should be placed in a given area of the surface. A related
approach is to use physically-based particle systems to sample an
implicit surface [56, 13]. Crossno and Angel [11] describe a system
for sampling isosurfaces, where they use the curvature to automatically modulate the repulsive forces.

SP
b)

ri
SR

SR
c)

d)

SP

Figure 2: An illustration of the paradigm: The possibly noisy or redundant point set (purple points) defines a manifold (purple curve).
This manifold is sampled with (red) representation points. The representation points define a different manifold (red curve). The spacing of representation points depends on the desired accuracy of theapproximation.

SP .
Compared to the earlier version [1] we give more and new details
on the (computational) properties of the surface definition:
Smooth manifold: The surface defined by the point set is a 2manifold and expected to be C∞ smooth, given that the points
are sufficiently close to the smooth surface being represented.
Bounded sampling error: Let SR be defined by the set of representation points {ri } ⊂ SP . The representation has boundederror ε, if d(SP , SR ) < ε, where d(·, ·) is the Hausdorff distance.
Local computation: For computing a point on the surface only alocal neighborhood of that point is required. This results ina small memory footprint, which depends only on the anticipated feature size and not the number of points (in contrast to
several other implicit surface definitions, e.g. those based onradial basis functions).
In addition to giving more implementation details for the rendering method from [1], we give a justification that builds on theerror bounds introduced here. This connection to the bounded error
substantiates the claimed properties of our rendering scheme:High quality: Since SR is a smooth surface, proper re-samplingleads to smooth silhouettes and normals resulting in superior
rendering quality at interactive frame rates.
Single step procedure: Re-sampling respects screen space resolution and guarantees sufficient sampling, i.e. no holes have to
be filled in a post-processing step.

2

Related work

Consolidation
Recent technological and algorithmic advances have improved theprocess of automatic acquisition of 3D models. Acquiring the geometry of an object starts with data acquisition, usually performed
with a range scanner. This raw data contains errors (e.g., line-ofsight error [21, 47]) mainly due to noise intrinsic to the sensor used2

Lee [30] uses a moving-least squares approach to the reconstruction of curves from unorganized and noisy points. He proposesa solution for reconstructing two and three-dimensional curves bythinning the point sets. Although his approach resembles the oneused here (and based on theory developed in [33]), his projection
procedure is different, and requires several iterations to converge toa clean point set (i.e., it is not actually a projection, but more of aconverging smoothing step).
An alternative to our point-based modeling mechanism is proposed by Linsen [36]. His work is based on extending the k neighborhood of a point to a “fan” by using an angle criterion (i.e., aneighborhood should cover the full 360 degrees around). Using thissimple scheme, Linsen proposes a variety of operations for point
clouds, including rendering, smoothing, and some modeling operations.

g
n

fi

r

pi

q

H
Figure 3: The MLS projection procedure: First, a local reference
domain H for the purple point r is generated. The projection of r
onto H defines its origin q (the red point). Then, a local polynomial
approximation g to the heights fi of points pi over H is computed.
In both cases, the weight for each of the pi is a function of the
distance to q (the red point). The projection of r onto g (the blue
point) is the result of the MLS projection procedure.

Point sample rendering
Following the pioneering work of Levoy and Whitted [35], several
researchers have recently proposed using “points” as the basic rendering primitive, instead of traditional rendering primitives, suchas triangulated models. One of the main reasons for this trend isthat in complex models the triangle size is decreasing to pixel resolution. This is particularly true for real-world objects acquired as“textured” point clouds [39].
Grossman and Dally [19] presented techniques for converting
geometric models into point-sampled data sets, and algorithms for
efficiently rendering the point sets. Their technique addresses several fundamental issues, including the sampling rate of conversion
from triangles to points, and several rendering issues, including
handling “gaps” in the rendered images and efficient visibility data
structures. The Surfels technique of Pfister et al [41] builds andimproves on this earlier work. They present alternative techniques
for the sampling of the triangle mesh, including visibility testing,texture filtering, and shading.
Rusinkiewicz and Levoy [46] introduce a technique which uses
a hierarchy of spheres of different radii to model a high-resolution
model. Their technique uses small spheres to model the vertices
at the highest resolution, and a set of bounding spheres to modelintermediate levels. Together with each spherical sample, they also
save other associated data, such as normals. Their system is capableof time-critical rendering, as it adapts the depth of tree traversal tothe available time for rendering a given frame. In [45], they show
how their system can be used for streaming models over a network.
Kalaiah and Varshney introduced an effective method for rendering point primitives that requires the computation of the principal
curvatures and a local coordinate frame for each point [25]. Theirapproach renders a surface as a collection of local neighborhoods,
and it is similar to our rendering technique proposed later, although
they do not use dynamic level of detail in their system.
All the above techniques account for local illumination. Schaufler and Jensen [48] propose to compute global illumination effects
directly on point-sampled geometry by a ray tracing technique. The
actual intersection point is computed, based on a local approximation of the surface, assuming a uniform sampling of the surface.
Point-based rendering suffers from the limited resolution of the
fixed number of sample points representing the model. At somedistance, the screen space resolution is high relative to the point
samples, which causes under-sampling. Tackling this problem byinterpolating the surface points in image space is limited. A better
approach is to re-sample the surface during rendering at the desired resolution in object-space, guaranteeing that sampling density
is sufficient with respect to the image space resolution.
Hybrid polygon-point approaches have been proposed. Cohenet al [10] introduces a simplification technique which transitions
triangles into (possibly multiple) points for faster rendering. Theirsystem uses an extension of DeFloriani et al’s Multi-Triangulation
data structure [15, 16]. A similar system has been developed byChen and Nguyen [8] as an extension of QSplat.

3

Defining the surface - projecting

Our approach relies on the idea that the given point set implicitly
defines a surface. We build upon recent work by Levin [33]. The
main idea is the definition of a projection procedure, which projects
any point near the point set onto the surface. Then, the MLS surface
is defined as the points projecting onto themselves. In the following, we explain the projection procedure, prove the projection andmanifold properties, motivate the smoothness conjecture, and, givedetails on how to efficiently compute the projection.

3.1

The projection procedure

Let points pi ∈ R3 , i ∈ {1, . . . , N}, be sampled from a surface S (possibly with a measurement noise). The goal is to project a pointr ∈ R3 near S onto a two-dimensional surface SP that approximates
the pi . The following procedure is motivated by differential geometry, namely that the surface can be locally approximated by afunction.
1. Reference domain: Find a local reference domain (plane) forr (see Figure 3). The local plane H = {x|hn, xi − D = 0, x ∈R3 }, n ∈ R3 , knk = 1 is computed so as to minimize a local
weighted sum of square distances of the points pi to the plane.
The weights attached to pi are defined as the function of the
distance of pi to the projection of r on the plane H, rather thanthe distance to r. Assume q is the projection of r onto H, thenH is found by locally minimizingN∑ hn, pi i − D2θ kpi − qk



(1)

i=1

where θ is a smooth, monotone decreasing function, which is
positive on the whole space. By setting q = r + tn for somet ∈ R, (1) can be rewritten as:N∑ n, pi − r − tn

i=1

3

2



θ kpi − r − tnk(2)We define the operator Q(r) = q = r + tn as the local minimum of Eq. 2 with smallest t and the local tangent plane Hnear r accordingly. The local reference domain is then given
by an orthonormal coordinate system on H so that q is theorigin of this system.

Gaussiand2

θ (d) = e− h2where h is a fixed parameter reflecting the anticipated spacing between neighboring points. By changing h the surface can be tuned
to smooth out features of size < h in S. More specifically, a smallvalue for h causes the Gaussian to decay faster and the approximation is more local. Conversely, large values for h result in a moreglobal approximation, smoothing out sharp or oscillatory featuresof the surface. Figure 4 illustrates the effect of different h values.

2. Local map: The reference domain for r is used to computea local bivariate polynomial approximation to the surface in aneighborhood of r (see Figure 3). Let qi be the projection ofpi onto H, and fi the height of pi over H, i.e fi = n · (pi − q).
Compute the coefficients of a polynomial approximation g sothat the weighted least squares error

3.3N∑ g(xi , yi ) − fi2

θ kpi − qk(3)Computing the projection

We explain how to efficiently compute the projection and what values should be chosen for the polynomial degree and h. Furthermore,
we discuss trade-offs between accuracy and speed.
Step 1 of the projection procedure is a non-linear optimization
problem. Usually, Eq. 2 will have more than one local minimum.
By definition, the local minimum with smallest t has to be chosen,
which means the plane should be close to r. For minimizing (2) wehave to use some iterative scheme, which descends toward the next
local minimum. Without any additional information we start witht = 0 and first approximate
 the normal. Note that in this case theweights θi = θ kpi − rk are fixed. Let B = {b jk }, B ∈ R3×3 thematrix of weighted covariances be defined by

i=1

is minimized. Here (xi , yi ) is the representation of qi in a local coordinate system in H. Note that, again, the distancesused for the weight function are from the projection of r
onto H. The projection P of r onto SP is defined by thepolynomial value at the origin, i.e. P(r) = q + g(0, 0)n =r + (t + g(0, 0))n.

3.2

(4)

Properties of the projection procedure

In our application scenario the projection property (i.e.
P(P(r)) = P(r)) is extremely important: We want to use
the above procedure to compute points exactly on the surface. The
implication of using a projection operator is that the set of points
we project and the order of the points we project do not change thesurface.
From Eq. 2 it is clear that if (t, n) is a minimizer for r then (s, n)
is a minimizer for r + (s −t)n. Assuming the minimization is global
in a neighborhood U ∈ R3 of q then Q is a projection operation inU, because the pair (0, n) is a minimizer for r − tn = q. Further,also q + g(0, 0)n = r + (t + g(0, 0))n is a minimizer and, thus, also
P a projection in U.
We like to stress that the projection property results from using
distances to q rather than r. If θ depended on r the procedure would
not be a projection because the values of θ would change along thenormal direction.
The surface SP is formally defined as the subset of all points
in R3 that project onto themselves. A simple counting argumentshows that this subset is a two-parameter family and, thus, SP a
2-manifold. Eq. 2 essentially has 6 degrees of freedom: 3 for r,
2 for n, and 1 for t. On the other hand there are 4 independentnecessary conditions: For a local minimum the partial derivativesof the normal and t have to be zero and for r to be on the surfacet = 0 is necessary. This leaves 6 − 4 = 2 free parameters for thesurface. It is clear from simple examples that the parameter family
includes manifolds which cannot be expressed as functions over R2 .
The particular charm of this surface definition is that it avoids
piecewise parameterizations. No subset of the surface is parameterized over a planar piece but every single point has its own support
plane. This avoids the common problems of piecewise parameterizations for shapes, e.g. parameterization dependence, distortions in
the parameterization, and continuity issues along the boundaries ofpieces.
However, in this approach we have the problem of proving continuity for any point on the surface because its neighbors have different support planes in general. The intuition for SP ∈ C∞ is, ofcourse, that (1) interpreted as a function R6 → R is C∞ so that alsoa particular kernel of its gradient is C∞ .
The approximation of single points is mainly dictated by the radial weight function θ . The weight function suggested in [33] is ab jk = ∑ θi pi − r j pi − rk .
ijkThen, the minimization problem (2) can be rewritten in bilinearform
min nT Bn
(5)
knk=1

and the solution of the minimization problem is given as the Eigenvector of B that corresponds to the smallest Eigenvalue.
If a normal is computed (or known in advance) it is fixed and thefunction is minimized with respect to t. This is a non-linear minimization problem in one dimension. In general, it is not solvablewith deterministic numerical methods because the number of minima is only bounded by the number N of points pi . However, in allpractical cases we have found that for t ∈ [−h/2, h/2] there is only
one local minimum. This is intuitively clear as h is connected to thefeature size and features smaller than h are smoothed out, i.e. features with distance smaller than h do not provide several minima.
For this reason we make the assumption that any point r to be projected is at most h/2 away from its projection. In all our practicalapplications this is naturally satisfied.
Using these prerequisites the local minimum is bracketed with±h/2. The partial derivativeN

2 ∑ n, pi − r − tni=1n, pi − r − tn
1+
h2

2

!
2

ekpi −r−tnk /h2(6)

can be easily evaluated together with the function itself. The derivative is exploited in a simple iterative minimization scheme as explained in [43, Chapter 10.3].
Once t 6= 0 fixing t and minimization with respect to the normal
direction is also a non-linear problem because q = r + tn changes
and, thus, also the weights change. The search space can be visualized as the tangent planes of a sphere with center point r and radiust. However, in practice we have found the normal (or q) to change
only slightly so that we approximate the sphere locally around thecurrent value of r + tn as the current plane defined t and n. On thisplane a conjugate gradient scheme [43] is used to minimize among

4

Figure 4: The effect of different values for parameter h. A point set representing an Aphrodite statue defines an MLS surface. The left sideshows an MLS surface resulting from a small value and reveals a surface micro-structure. The right side shows a larger value for h, smoothingout small features or noise.
all q on the plane. The main idea is to fix a subspace for minimization in which n cannot vanish so that the constraint knk = 1 is always satisfied. The use of a simple linear subspace makes the computation of partial derivatives efficient and, thus, conjugate gradientmethods applicable. Clearly, this search space effectively changes
t resulting in a theoretically worse convergence behavior. In practice, the difference between the sphere and the plane is small for theregion of interest and the effect is not noticeable.
Using these pieces the overall implementation of computing the
support plane looks like this:coefficients of the polynomial leads to a system of linear equationsof size equal to the number of coefficients, e.g. 10 for a third degree
polynomial.
Through experimentation, we found that high degree polynomials are likely to oscillate. Polynomials of degree 3 to 4 have proved
to be very useful as they produce good fits of the neighborhood, do
not oscillate, and are quickly computed.

3.4

Data structures and trade-offs

The most time-consuming step in computing the projection of apoint r is collecting the coefficients from each of the pi (i.e. computing the sum). Implemented naively, this process takes O(N)time, where N is the number of points. We exploit the effect of
the quickly decreasing weight function in two ways:

Initial normal estimate in r: The normal in a point might begiven as part of the input (e.g. estimate from range images)
or could be computed when refining a point set surface (e.g.
from close points in the already processed point set). If nonormal is available it is computed using the Eigenvector of
the matrix of weighted co-variances B.

1. In a certain distance from r the weight function is effectively
zero. We call this the neglect distance dn , which depends
solely on h. A regular grid with cell size 2dn is used to partition the point set. For each projection operation a maximumof 8 cells is needed. This results in a very small memory footprint and yields a simple and effective out of core implementation, which makes the storage requirements of this approach
independent of the total size of the point set.

Iterative non-linear minimization: The following two steps arerepeated as long as any of the parameters changes more thana predefined ε:
1. Minimize along t, where the minimum is initially bracketed by t = ±h/2.
2. Minimize q on the current plane H : (t, n) using conjugate gradients. The new value for q = r + tn leads to
new values for the pair (t, n).

2. Using only the selected cells the terms are collected using ahierarchical method inspired by solutions to the N-body problem [4]. The basic observation is, that a cluster of points far
from r can be combined into one point. To exploit this idea,
each cell is organized as an Octree. Leaf nodes contain the pi ;inner nodes contain information about the number of points
in the subtree and their centroid. Then, terms are collected

The last pair (t, n) defines the resulting support plane H.
The second step of the projection procedure is a standard linear
least squares problem: Once the plane H is computed, the weights
θi = θ (kpi − qk) are known. The gradient of (3) over the unknown

5

from the nodes of the Octree. If the node’s dimension is muchsmaller than its distance to r, the centroid is used for computing the coefficients; otherwise the subtree is traversed. In
addition, whole nodes can be neglected if their distance to r islarger than dn .
The idea of neglecting points could be also made independent of
numerical issues by using a compactly supported weight function.
However, the weight function has to be smooth. An example forsuch an alternative is θ (x) = 2x3 − 3x2 + 1.
A simple way to trade accuracy for speed is to assume that theplane H passes through the point to be projected. This assumption
is reasonable for input points, which are expected to be close to
the surface they define (e.g. input that has been smoothed). This
simplification saves the cost of the iterative minimization scheme.

3.5

(a)

Results

Actual timings for the projection procedure depend heavily on thefeature size h. On a standard Pentium PC the points of the bunny
were projected at a rate of 1500-3500 points per second. This meansthe 36K points of the bunny are projected onto the surface they
define themselves in 10-30 seconds. Smaller values for h lead tofaster projection since the neighborhoods and, thus, the number ofpoints taken into account are smaller.
As has been stressed before, the memory requirements mainly
depend on the local feature size h. As long as h is small with respect
to the total diameter of the model the use of main memory of the
projection procedure is negligible.

4

(c)
(b)
Figure 5: Points acquired by range scanning devices or vertices
from a processed mesh typically have uneven sampling density on
the surface (a). The sampling techniques discussed here allow toevenly re-sample the object (b) to ensure a sufficient density forfurther processing steps, for example rendering (c).

Approximation error

Consider again the setting depicted in Figure 2. Input points {pi }define a surface SP , which are then represented by a set of points{ri } ∈ SP . However, the set {ri } defines a surface SR which approximates SP . Naturally, we would like to have an upper bound on the
distance between SR and SP .
From differential geometry, we know that a smooth surface canbe locally represented as a function over a local coordinate system. We also know that in approximating a bivariate function fby a polynomial g of total degree m, the approximation error iskg − f k ≤ M · hm+1 [32]. The constant M involves the (m + 1)-thderivatives of f , i.e. M ∈ O(k f (m+1) k). In the case of surface approximation, since SP is infinitely smooth, there exists a constant
Mm+1 , involving the (m + 1)-th derivatives of SP , such thatkSP − SR k ≤ Mm+1 hm+1

normal ni , a corresponding orthonormal system {ui , vi , ni }, and an
origin at qi .
Gi = {qi + x · ui + y · vi + gi (x, y) · ni | (x, y) ∈ [−h, h]2 }.

These error bounds explain nicely why curvature (see e.g.
[17]) is an important criterion for piecewise linear approximations(meshes): The error of piecewise linear functions depends linearlyon the second derivatives and the spacing of points. This means,more points are needed where the curvature is higher.
However, when the surface is approximated with higher orderpolynomials, curvature is irrelevant to the approximation error. Using cubic polynomials, the approximation error depends on the
fourth order derivatives of the surface. Note that our visual system cannot sense smoothness beyond second order [38]. From thatpoint of view, the sampling density locally depends on an “invisible” criterion. We have found it to be sufficient to fix a spacing hof points ri on SP . A smaller h will cause the error to be smaller.
However, using higher order polynomials to approximate SP causes
the error to decrease faster when h is reduced.

(7)

where SR is computed using polynomials of degree m.
Note that this error bound holds for the MLS surface SR , which
is obtained by the projection procedure applied to the data set R,
using the same parameter h as for SP . Moreover, the same type ofapproximation order holds for piecewise approximation: Assume
each point ri defines a support plane so that SP is a function overa patch [−h, h]2 around ri on the plane and, further, that the points{ri } leave no hole of radius more than h on the surface. Then,the above arguments hold and the approximation error of the unionof local, non-conforming, polynomial patches Gi around points ri
approximates SP as
SP −[Gi ≤ Mm+1 h

m+1

.

(9)

5

Generating the representation point set

A given point set might have erroneous point locations (i.e. isnoisy), may contain too many points (i.e. is redundant) or notenough points (i.e. is under-sampled).
The problem of noise is handled by projecting the points onto
the MLS surface they define. The result of the projection procedure is a thin point set. Redundancy is avoided by decimating the(8)Here, Gi is the polynomial patch defined by an MLS polynomial gi for the point ri , with corresponding reference plane Hi with

6

a)

b)

Figure 6: Noisy input points (green points) are projected onto theirsmooth MLS curve (orange points). The figures in (a) show thepoint sets and a close-up view. The decimation process is shown in(b). Points are color-coded as in Figure 7.
Figure 8: The up-sampling process: Points are added at verticesof the Voronoi diagram. In each step, the vertex with the largest
empty circle is chosen. The process is repeated until the radiusof the largest circle is smaller than a specified bound. The wavy
torus originally consisting of 800 points has been up-sampled to
20K points.

point set, taking care that it persists to be a good approximation of
the MLS surface defined by the original point set. In the case ofunder-sampling, the input point set needs to be up-sampled. In thefollowing sections, we show techniques to remove and add points.
Figure 5 illustrates the idea of re-sampling at the example of the
Buddha statue. Using the techniques presented below it is possible
to re-sample the geometry to be evenly sampled on the surface.

5.1

in the sparsest set, all of the points have a high error, that is, theirremoval will cause a large error. As the decimation proceeds, fewer
points remain and their importance grows and the error associatedwith them is larger. Figure 7 shows the decimation process in 3Dwith corresponding renderings of the point sets.

Down-sampling

Given a point set, the decimation process repeatedly removes thepoint that contributes the smallest amount of information to theshape. The contribution of a point to the shape or the error of thesampling is dictated by the definition of the shape. If the point set isreconstructed by means of a triangulation, criteria from the specifictriangulation algorithm should control the re-sampling. Criteria include the distance of points on the surface [23], curvature [17], or
distance from the medial axis of the shape [2]. For a direct displayof the point set on a screen homogeneous distribution of the points
over the surface is required [19, 41]. Here, we derive a criterion
motivated by our definition of the surface.
The contribution of a projected point pi to the surface SP canbe estimated by comparing SP with SP−{p } . Computing the Hausi
dorff distance between both surfaces is expensive and not suitable
for an iterative removal of points of a large point set. Instead, we
approximate the contribution of pi by its distance from its projection onto the surface SP−{p } . Thus, we estimate the difference ofi
SP and SP−{p } by projecting pi onto SP−{p } (projecting pi underiithe assumption it was not part of P).
The contribution values of all points are inserted into a priority
queue. At each step of the decimation process, the point with thesmallest error is removed from the point set and from the priority
queue. After the removal of a point, the error values of nearby
points have to be recalculated since they might have been affectedby the removal. This process is repeated until the desired numberof points is reached or the contributions of all points exceed someprespecified bound.
Figure 6 illustrates our decimation process applied on the set of
red points depicted in (a). First, the red points are projected on the
MLS surface to yield the blue points. A close-up view over a part ofthe points shows the relation between the input (red) points and the
projected points. In (b), we show three snapshots of the decimation
process, where points are colored according to their error value;
blue represents low error and red represents high error. Note that

5.2

Up-sampling

In some cases, the density of the point set might not be sufficient for
the intended usage (e.g. direct point rendering or piecewise reconstructions). Motivated by the error bounds presented in Section 4
we try to decrease the spacing among points. Additional points
should be placed (and then projected to the MLS surface) where
the spacing among points is larger then a specified bound.
The basic idea of our approach is to compute Voronoi diagramson the MLS surface and add points at vertices of this diagram. Notethat the vertices of the Voronoi diagram are exactly those points on
the surface with maximum distance to several of the existing points.
This idea is related to Lloyd’s method [37], i.e techniques using
Voronoi diagrams to achieve a certain distribution of points [40].
However, computing the Voronoi diagram on the MLS surface isexcessive and local approximations are used instead. More specifically, our technique works as follows: In each step, one of the existing points is selected randomly. A local linear approximation isbuilt and nearby points are projected onto this plane. The Voronoi
diagram of these points is computed. Each Voronoi vertex is thecenter of a circle that touches three or more of the points without
including any point. The circle with largest radius is chosen andits center is projected to the MLS surface. The process is repeatediteratively until the radius of the largest circle is less than a userspecified threshold. (see Figure 8). At the end of the process, thedensity of points is locally near-uniform on the surface. Figure 8
shows the original sparse point set containing 800 points, and thesame object after adding 20K points over its MLS surface.

7

a)

b)

c)

d)

e)

Figure 7: The point set representing an Aphrodite statue is projected onto a smooth MLS-surface. After removing redundant points, a setof 37K points represents the statue (a). The corresponding rendering is shown in (b). The point set is decimated using point removal. An
intermediate stage of the reduction process is shown in (c). Note that the points are color-coded with respect to their importance. Blue points
do not contribute much to the shape and might be removed; red points are important for the definition of the shape. The final point set in (e)contains 20K points. The corresponding rendering is depicted in (d) and is visually close to the one in (b).

6Rendering

The challenge of our interactive point rendering approach is to use
the representation points and (when necessary) create new points
by sampling the implicitly defined surface at a resolution sufficient
to conform to the screen space resolution (see Figure 9 for an illustration of that approach).
Usually, the representation points are not sufficient to render theobject in screen space. In some regions, it is not necessary to render
all points as they are occluded, backfacing, or have higher density
than needed. However, typically, points are not dense enough tobe projected directly as a single pixel and more points need to begenerated by interpolation in object space.

6.1

hFigure 10: The patch size of a polynomial: Points inside a ballof radius h around the red point r are projected onto the support
plane of the red point. The patch size is defined as the bounding
disk (in local coordinates) of the projections. Note that using a diskof radius h would lead to unpleasant effects in some cases, as thepolynomial might leave the ball of radius h.

Culling and view dependency

The structure of our rendering system is similar to QSplat [46]. The
input points are arranged into a bounding sphere hierarchy. For eachnode, we store a position, a radius, a normal coverage, and optionally a color. The leaf nodes additionally store the orientation of
the support plane and the coefficients of the associated polynomial.
The hierarchy is used to cull the nodes with the view frustum and to
apply a hierarchical backface culling [29]. Note that culling is important for our approach since the cost of rendering the leaf nodes(evaluating the polynomials) is relatively high compared to simpler
primitives. Moreover, if the traversal reaches a node with an extent
that projects to a size of less than a pixel, this node is simply projected to the frame-buffer without traversing its subtree. When thetraversal reaches a leaf node and the extent of its bounding sphere
projects to more than one pixel in screen space, additional points
have to be generated.
The radius of a bounding sphere in leaf nodes is simply set to theanticipated feature size h. The bounding spheres for inner nodes
naturally bound all of their subtree’s bounding spheres. To computethe size of a bounding sphere in pixel space the radius is scaledby the model-view transform and then divided by the distance ofcenter and eye point in z-direction to accommodate the projective
transform.

6.2

It is clear that the union of these polynomial patches is not acontinous surface. However, the Hausdorff-error of this approximation is not worse than the error of the surface computed by projecting every point using the operator P (see Section 4). Since theHausdorff-error in object space limits the screen space error the error bound can be used to make the approximation error invisble in
the resulting image.
However, to conform with the requirements formulated in Section 4 the point set and the associated polyonmials are required tobe near-uniform on the surface. It might be necessary to first process a given point set with the up-sampling methods presented in
Section 5. This way, we ensure that the local, non-conforming (i.e.
overlapping or intersecting) polynomials are a good approximation
to the surface inside a patch [−h, h]2 around a point and, thus, the resulting image shows a smooth surface. However, most dense pointsets can be readily displayed with the approach presented here. For
example, Figure 11 shows several renderings of the original Stanford Bunny data.
It is critical to properly define the extent of a polynomial patch on
the supporting plane, such that neighboring patches are guaranteed
to overlap (to avoid holes) but do not overlap more than necessary.
Since no inter-point connectivity information is available, it is unclear which points are immediate neighbors of a given point on thesurface.
To compute the extent of a polynomial patch associated to thepoint pi on the support plane H all points inside a ball of radius

Sampling additional points

The basic idea is to generate a grid of points sufficient to cover theextent of a leaf node. However, projecting the grid points using
the method described in Section 3 is too slow in interactive applications. The main idea of our interactive rendering approach is tosample the polynomials associated with the representation points
rather than really projecting the points.

8

Figure 9: Up-sampling could be used generate enough points on the surface to conform with the resolution of the image to be rendered. The
right image shows a close-up rendered with splats.
h around the projection q = Q(pi ) are collected. These points areprojected to the support plane H, which leads to local (u, v) coordinates for each projected point. The extent is defined by a circle
around q that encloses all projected points. More specifically, assume q has the local coordinate (0,0), the radius of the circle isgiven by the largest 2-norm of all local (u, v) coordinates of projected points. Since the spacing of points is expected to be less thanh, patches of neighboring points are guaranteed to overlap.
Note that using a constant extent (e.g. a disk of radius h on the
support plane) can lead to errors, as the polynomial g over H mightleave the ball of radius h, in which a good approximation of thepoint set is expected. Figure 10 illustrates the computation of thepatch sizes.
The grid spacing d should be so computed that neighboring
points have a screen space distance of less than a pixel. Thus, thegrid spacing depends on the orientation of the polyonmial patchwith respect to the screen. Since the normals change on each polyonmial patch we rather use a simple heuristic to conservatively estimate d: The grid spacing d is computed so that a grid perpendicular
to the viewing direction is sufficiently sampled in image space.
If the grid is, indeed, perpendicular to the viewing direction, thesampling is also correct on the polynomial. If the grid is not perpendicular to the viewing direction, the projected area might be oversampled or undersampled depending on the orientation of the support plane and the derivatives of the polyonmial. Note, however,that sufficient sampling is guaranteed if the derivatives are less than
1. In practice we have rarely encountered higher derivatives so we
decided not to evaluate the maximum derivatives of all polynomial
patches. However, this could be done in a preprocess and the density could be adjusted accordingly.
Upon the view-dependent grid spacing d, the polynomials areevaluated by a forward difference approach, where the polynomialis scanned across its extent in its local u, v parametric space. The
affine map transforming from support plane coordinates to world
coordinates is factored into polynomial evaluation, thus, generating points in world coordinates. These points are then fed into thegraphics pipeline to be projected to the screen.
Surprisingly, we have found that quad meshes are processed

faster by current graphics hardware than a set of points. For thisreason we use quad meshes to represent the polynomial patches.
This has the additional advantage that during lazy evaluation (seebelow) no holes occur.

6.3

Grid pyramids

The time-critical factor is the view-dependent evaluation of thepoints on the polynomial. Optimally, these are recomputed whenever the projected screen space size changes. To accelerate the rendering process, we store a grid pyramid with various resolutions per
representation point. Initially, the pyramid levels are created, but nogrid is actually evaluated. When a specific grid resolution is needed,
the system creates and stores the level that slightly oversamples thepolynomial for a specific resolution, such that small changes in theviewing position do not result in new evaluations.
To enhance the interactivity of our approach we use lazy evaluation during rotating or zooming. Once the viewer stops moving, aproper grid is chosen from the pyramid.

6.4

Results

We have tested our approach on a variety of point sets. Figure 11
shows the renderings of the Stanford Bunny. In (a), the original
point set is shown. Splatting (b) is not leading to good results, because the model is not sampled densely enough. The traditional
Gouraud-shaded mesh in (c) and (g) is compared to our approachin (d) and (h). Note the accuracy of the highlights. The nonconforming local polynomial patches are shown color-coded in (e)and (f). An example of a environment mapping to demonstrate thenormal continuity is given in Figure 12. Note that silhouettes andnormals are smooth, which leads to less distortions on the boundary
and in the reflections.
The frame rates we achieve are mainly dictated by the numberof visible representation points (i.e. graph traversal time) and thenumber of pixels to be filled. All models depicted in the paper
are displayed at more than 5 frames per second in a 5122 screen
window (see the accompanying video for more information). The

9

a)

b)

c)

d)

e)

f)

g)

h)

Figure 11: The Stanford Bunny: The points defining the bunny are depicted in (a) (some points are culled). Points are splatted in (b) tosatisfy screen space resolution. Note the difference of a piecewise linear mesh over the points (c) and close-up in (g) to the rendering of
non-conforming polynomial patches (d) and (h). The patches are color-coded in (e) and (f).

number of representation points ranges from 1000 (for the torus)
to 900K (for the angel statue). Tests are performed on a PC with
GeForce2 graphics board.

7

Conclusion

In differential geometry, a smooth surface is characterized by theexistence of smooth local maps at any point. In this work we use
this as a framework to approximate a smooth surface defined bya set of points and we introduced new techniques to resample thesurface to generate an adequate representation of the surface.
To render such surfaces, the surface is covered by a finite number, as small as possible, of non-conforming, overlapping, polynomial patches. We showed that the error of these approximations isbounded and dependent on the spacing among points. Thus, it ispossible to provide a point set representation that conforms with a
specified tolerance.
Our paradigm for representing surfaces advocates the use of apoint set (without connectivity) as a representation of shapes. This
representation is universal in the sense that it is used from the beginning (i.e. acquisition) to the end (i.e. rendering) of a graphicalrepresentation’s life cycle. Moreover, we believe that this work is
a step towards rendering with higher order polynomials. Note that
we have used simple primitives (points) to achieve this goal. Thisadmits to the current trend of integrating high quality rendering features into graphics hardware.
It would be interesting to integrate our approach with combinatorial methods such as the one of Amenta et al [2]. This would combine topological guarantees with the additional precision of higher
order approximations and the possibility of smoothing out noise or
small features.
Using different values for h, it is easy to generate more smoothor more detailed versions of a surface from one point set (see, for

Figure 12: Comparison of mesh rendering with our technique withenvironment mapping. The left column shows renderings of a meshconsisting of 1000 vertices. The right column shows our techniqueusing the vertices as input points. The environment maps emphasize the improved normal and boundary continuity.

10

example, Figure 4). A set of different versions could be used as asmooth-to-detailed hierarchy and would allow for multiresolution
modeling [27]. Of course, h is not necessarily a global parameterand could be adapted to the local feature size. Varying h has several implications and utility in handling point sets (see [30] for anice introduction to the issues in two dimensions), such as properly accounting for differences in sampling rate and levels of noise
during the acquisition process. Also, non radial functions might benecessary to properly account for sharp features in the models.

[20] A. P. Gueziec, G. Taubin, F. Lazarus, and W. Horn. Converting sets of polygons
to manifold surfaces by cutting and stitching. IEEE Visualization ’98, pages
383–390, October 1998.
[21] P. Hebert, D. Laurendeau, and D. Poussart. Scene reconstruction and description: Geometric primitive extraction from multiple view scattered data. In IEEE
Computer Vision and Pattern Recognition 1993, pages 286–292, 1993.
[22] H. Hoppe, T. DeRose, T. Duchamp, M. Halstead, H. Jin, J. McDonald,
J. Schweitzer, and W. Stuetzle. Piecewise smooth surface reconstruction. Proceedings of SIGGRAPH 94, pages 295–302, July 1994.
[23] H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, and W. Stuetzle. Surface
reconstruction from unorganized points. Computer Graphics (Proceedings of
SIGGRAPH 92), 26(2):71–78, July 1992.
[24] H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, and W. Stuetzle. Mesh optimization. Proceedings of SIGGRAPH 93, pages 19–26, August 1993.

Acknowledgements[25] A. Kalaiah and A. Varshney. Differential point rendering. In Rendering Techniques 2001: 12th Eurographics Workshop on Rendering, pages 139–150. Eurographics, June 2001. ISBN 3-211-83709-4.

We gratefully acknowledge the helpful comments of Markus Gross,Jörg Peters, Ulrich Reif, and several anonymous reviewers.
This work was supported by a grant from the Israeli Ministry of
Science, a grant from GIF (German Israeli Foundation) and a grantfrom the Israeli Academy of Sciences (center of excellence). The
bunny model is courtesy of the Stanford Computer Graphics Laboratory. The angel statue in Figure 1 was scanned by Peter Neugebauer at Fraunhofer IGD in Darmstadt, Germany using a structured
light scanner and the QTSculptor system.

[26] A. Kaufman, D. Cohen, and R. Yagel. Volume graphics. IEEE Computer,26(7):51–64, July 1993.
[27] L. P. Kobbelt, T. Bareuther, and H.-P. Seidel. Multiresolution shape deformations for meshes with dynamic vertex connectivity. Computer Graphics Forum,
19(3):249–259, Aug. 2000.
[28] V. Krishnamurthy and M. Levoy. Fitting smooth surfaces to dense polygon
meshes. Proceedings of SIGGRAPH 96, pages 313–324, August 1996.
[29] S. Kumar, D. Manocha, W. Garrett, and M. Lin. Hierarchical back-face computation. In X. Pueyo and P. Schröder, editors, Eurographics Rendering Workshop
1996, pages 235–244, New York City, NY, June 1996. Eurographics, Springer
Wien.

References[30] I.-K. Lee. Curve reconstruction from unorganized points. Computer Aided Geometric Design, 17(2):161–177, February 2000.

[1] M. Alexa, J. Behr, D. Cohen-Or, S. Fleishman, D. Levin, and C. T. Silva. Point
set surfaces. In IEEE Visualization 2001, pages 21–28, October 2001. ISBN
0-7803-7200-x.

[31] Z. Lei, M. M. Blane, and D. B. Cooper. 3L fitting of higher degree implicit polynomials. In Proceedings of Third IEEE Workshop on Applications of Computer
Vision, Sarasota, Florida, USA, Dec. 1996.

[2] N. Amenta, M. Bern, and M. Kamvysselis. A new voronoi-based surface reconstruction algorithm. Proceedings of SIGGRAPH 98, pages 415–422, July 1998.

[32] D. Levin. The approximation power of moving least-squares. Mathematics ofComputation, 67(224), 1998.

[3] C. L. Bajaj, F. Bernardini, and G. Xu. Automatic reconstruction of surfaces
and scalar fields from 3d scans. Proceedings of SIGGRAPH 95, pages 109–118,
August 1995.

[33] D. Levin. Mesh-independent surface interpolation. Advances in Computational
Mathematics, in press, 2001.

[4] J. Barnes and P. Hut. A hierarchical O(N log N) force calculation algorithm.
Nature, 324:446–449, Dec. 1986. Technical report.

[34] M. Levoy, K. Pulli, B. Curless, S. Rusinkiewicz, D. Koller, L. Pereira, M. Ginzton, S. Anderson, J. Davis, J. Ginsberg, J. Shade, and D. Fulk. The digitalmichelangelo project: 3d scanning of large statues. Proceedings of SIGGRAPH
2000, pages 131–144, July 2000.

[5] F. Bernardini, J. Mittleman, H. Rushmeier, C. Silva, and G. Taubin. The ballpivoting algorithm for surface reconstruction. IEEE Transactions on Visualization and Computer Graphics, 5(4):349–359, October - December 1999.

[35] M. Levoy and T. Whitted. The use of points as a display primitive. Tr 85-022,
University of North Carolina at Chapel Hill, 1985.

[6] P. Besl and N. McKay. A method for registration of 3-d shapes. IEEE
Transaltions on Pattern Analysis and Machine Intelligence, 14(2):239–256,
February 1992.

[36] L. Linsen. Point cloud representation. Technical report, Fakultät fuer Informatik,
Universität Karlsruhe, 2001.

[7] J.-D. Boissonnat. Geometric structues for three-dimensional shape representation. ACM Transactions on Graphics, 3(4):266–286, October 1984.

[37] S. P. Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28:128–137, 1982.

[8] B. Chen and M. X. Nguyen. Pop: a hybrid point and polygon rendering system
for large data. In IEEE Visualization 2001, pages 45–52, October 2001. ISBN
0-7803-7200-x.

[38] D. Marr. Vision: A Computational Investigation into the Human Representation
and Processing of Visual Information. W. H. Freeman, Sept. 1983.
[39] D. K. McAllister, L. F. Nyland, V. Popescu, A. Lastra, and C. McCue. Real-timerendering of real-world environments. Eurographics Rendering Workshop 1999,
June 1999.

[9] P. Cignoni, C. Montani, and R. Scopigno. A comparison of mesh simplification
algorithms. Computers & Graphics, 22(1):37–54, February 1998.
[10] J. D. Cohen, D. G. Aliaga, and W. Zhang. Hybrid simplification: combining
multi-resolution polygon and point rendering. In IEEE Visualization 2001, pages
37–44, October 2001. ISBN 0-7803-7200-x.

[40] A. Okabe, B. Boots, and K. Sugihara. Spatial Tesselations — Concepts andApplications of Voronoi Diagrams. Wiley, Chichester, 1992.

[11] P. Crossno and E. Angel. Isosurface extraction using particle systems. IEEE
Visualization ’97, pages 495–498, November 1997.

[41] H. Pfister, M. Zwicker, J. van Baar, and M. Gross. Surfels: Surface elements
as rendering primitives. Proceedings of SIGGRAPH 2000, pages 335–342, July
2000.

[12] B. Curless and M. Levoy. A volumetric method for building complex models
from range images. Proceedings of SIGGRAPH 96, pages 303–312, August
1996.

[42] V. Pratt. Direct least-squares fitting of algebraic surfaces. Computer Graphics(Proceedings of SIGGRAPH 87), 21(4):145–152, July 1987.

[13] L. H. de Figueiredo, J. de Miranda Gomes, D. Terzopoulos, and L. Velho.
Physically-based methods for polygonization of implicit surfaces. Graphics Interface ’92, pages 250–257, May 1992.

[43] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical
Recipes in C: The Art of Scientific Computing (2nd ed.). Cambridge University
Press, 1992.

[14] M. Desbrun, M. Meyer, P. Schröder, and A. H. Barr. Implicit fairing of irregularmeshes using diffusion and curvature flow. Proceedings of SIGGRAPH 99, pages
317–324, August 1999.

[44] R. Ramamoorthi and J. Arvo. Creating generative models from range images.
Proceedings of SIGGRAPH 99, pages 195–204, August 1999.
[45] S. Runsinkiewicz and M. Levoy. Streaming qsplat: A viewer for networked
visualization of large, dense models. In to appear in Proceedings of Symposiumon Interactive 3D Graphics. ACM, 2001.

[15] L. D. Floriani, P. Magillo, and E. Puppo. Building and traversing a surface atvariable resolution. IEEE Visualization ’97, pages 103–110, November 1997.
ISBN 0-58113-011-2.

[46] S. Rusinkiewicz and M. Levoy. Qsplat: A multiresolution point rendering system
for large meshes. Proceedings of SIGGRAPH 2000, pages 343–352, July 2000.

[16] L. D. Floriani, P. Magillo, and E. Puppo. Efficient implementation of multitriangulations. IEEE Visualization ’98, pages 43–50, October 1998. ISBN 08186-9176-X.

[47] M. Rutishauser, M. Stricker, and M. Trobina. Merging range images of arbitrarily
shaped objects. In IEEE Computer Vision and Pattern Recognition 1994, pages
573–580, 1994.

[17] M. Gopi, S. Krishnan, and C. T. Silva. Surface reconstruction based on lower
dimensional localized delaunay triangulation. Computer Graphics Forum, 19(3),August 2000.

[48] G. Schaufler and H. W. Jensen. Ray tracing point sampled geometry. Rendering
Techniques 2000: 11th Eurographics Workshop on Rendering, pages 319–328,
June 2000.

[18] A. Goshtasby and W. D. O’Neill. Surface fitting to scattered data by a sum ofgaussians. Computer Aided Geometric Design, 10(2):143–156, April 1993.

[49] M. Soucy and D. Laurendeau. Surface modeling from dynamic integration ofmultiple range views. In 1992 International Conference on Pattern Recognition,
pages I:449–452, 1992.

[19] J. P. Grossman and W. J. Dally. Point sample rendering. Eurographics Rendering
Workshop 1998, pages 181–192, June 1998.

11

[50] M. Sramek and A. E. Kaufman. Alias-free voxelization of geometric objects.
IEEE Transactions on Visualization and Computer Graphics, 5(3):251–267, July
- September 1999. ISSN 1077-2626.
[51] G. Taubin. A signal processing approach to fair surface design. Proceedings of
SIGGRAPH 95, pages 351–358, August 1995.
[52] G. Turk. Re-tiling polygonal surfaces. Computer Graphics (Proceedings of
SIGGRAPH 92), 26(2):55–64, July 1992.
[53] G. Turk and M. Levoy. Zippered polygon meshes from range images. Proceedings of SIGGRAPH 94, pages 311–318, July 1994.
[54] S. W. Wang and A. E. Kaufman. Volume sculpting. 1995 Symposium on Interactive 3D Graphics, pages 151–156, April 1995. ISBN 0-89791-736-7.
[55] M. Wheeler, Y. Sato, and K. Ikeuchi. Consensus surfaces for modeling 3d objects
from multiple range images. In IEEE International Conference on Computer
Vision 1998, pages 917–924, 1998.
[56] A. P. Witkin and P. S. Heckbert. Using particles to sample and control implicit
surfaces. Proceedings of SIGGRAPH 94, pages 269–278, July 1994.

12