Using particles to sample and control implicit surfaces.
We present a new particle-based approach to sampling and controlling implicit surfaces. A simple constraint locks a set of particles onto a surface while the particles and the surface move. We use the constraint to make surfaces follow particles, and to make particles follow surfaces. We implement control points for direct manipulation by specifying particle motions, then solving for surface motion that maintains the constraint. For sampling and rendering, we run the constraint in the order direction, creating floater particles that roam freely over the surface. Local repulsion is used to make floaters spread evenly across the surface. By varying the radius of repulsion adaptively, and fissioning or killing particles based on the local density, we can achieve good sampling distributions very rapidly, and maintain them even in the face of rapid and extreme deformations and changes in surface topology.
Introduction
Implicit surfaces have proven to be useful for modeling, ani-
mation, and visualization. One appeal of implicit models is that
new surfaces can be created by adding or otherwise combining the
functions that define them, producing a variety of subtle and interesting
shape effects. Another is their role in the visualization of
volume data. In addition, the implicit representation lends itself to
such calculations as ray/surface intersection and inside/outside test.
However, implicit surfaces suffer from two serious drawbacks: first,
although well suited to ray tracing, they are not easily rendered at interactive
speeds, reflecting the underlying problem that it is difficult
to sample them systematically. This is particularly a problem if we
wish to render time-varying surfaces in real time, which is vital for
interactive sculpting. Second, the shapes of implicit surfaces have
proven to be more difficult to specify and control than those of their
parametric counterparts.
Mail to the authors should be addressed to the Department of Computer
Science, Carnegie Mellon University, 5000Forbes Ave, Pittsburgh PA 15213,
USA. Email should be addressed to Andrew Witkin as aw@cs.cmu.edu, and
to Paul Heckbert as ph@cs.cmu.edu.
c
ACM. Reprinted from Computer Graphics, Proc. SIGGRPAH '94.
In this paper, we present a new particle-based approach to sampling
and shape control of implicit surfaces that addresses these
problems. At the heart of our approach is a simple constraint that
locks a collection of particles onto an implicit surface while both
the particles and the surface move. We can use the constraint to
make the surface follow the particles, or to make the particles follow
the surface. Our formulation is differential: we specify and solve
for velocities rather than positions, and the behavior of the system
is governed by differential equations that integrate these velocities
over time.
We control surface shape by moving particles interactively, solving
for surface motion that keeps the particles on the surface. This
technique lets us pin down the surface at some points while interactively
dragging others. These particles act as control points for
direct manipulation of the surface.
For sampling and rendering, we run the constraint in the other
direction, creating particles that may roam freely over the surface,
but are compelled to follow it as it moves. We call these particles
floaters. Our starting point is the idea that uniform sampling density
can be achieved by making the particles repel each other. This
approach was used by Turk [29] to resample polygon meshes, and
by Figueiredo et al. [12] to sample implicit surfaces.
Simple repulsion can work quite well for stationary surfaces, but
only if a reasonably good initial sampling is available. If large-scale
non-uniformities exist, convergence can be extremely slow
for even moderate sampling densities. We eliminate the need for
a good starting point, and dramatically accelerate convergence, by
employing an iterative "fissioning" approach, in which we start with
a small number of particles and a large radius of repulsion, allow
them to equilibrate, then split each particle in two, reducing the
radius of repulsion. We repeat the process until the desired sampling
density is reached. Each level thus inherits a distribution that is
already uniform at large scale, requiring just a few iterations to iron
out the local irregularities.
Global fissioning still fails to handle surfaces that move and de-
form, since irregularities can arise after the density becomes high.
We introduce a local adaptive repulsion algorithm, in which the repulsion
radius and particle birth and death are regulated based on
local density. This method is fast enough to maintain good sampling
even in the face of rapid and extreme surface motion and deformation

The remainder of the paper is organized as follows: we begin
by discussing previous related work. Then we introduce the basic
constraint mechanism that attaches particles to surfaces. Next we
describe the use of particles for surface shape control. We then
explain our adaptive repulsion sampling algorithm. After describing
the implementation and results, we conclude with a discussion of
future work.
Previous Work
Related work can be divided into two categories: sampling methods
and control methods.
2.1 Sampling Methods
Related research on surface sampling includes both particle-based
sampling techniques and polygonization techniques for implicit surfaces

Turk used repelling particles on surfaces to uniformly resample
a static surface [28] and to simplify a polygonization by reducing
the number of polygons [29]. Hoppe et al. also explored mesh
simplification, framing it as an optimization problem with penalties
for geometric error, number of samples, and edge length [16]. Their
method did not restrict the points to a surface, however, as Turk's
and ours do.
Szeliski and Tonnesen used oriented particles to model surfaces
[27]. Their technique allowed the user to move the particles interac-
tively, employing short-range repulsion and long-range attraction to
keep the particles from clumping or flying apart. The system generated
a surface by connecting neighboring particles appropriately,
but it did not manipulate a high level representation for a surface,
such as a parametric patch or an implicit function, as ours does. The
output of their system was a triangulation. Their system bears a superficial
resemblance to ours because we both use disks to visualize
the surface, but in other respects our techniques are quite different.
An implicit surface, also called an iso-surface, is the set of points
x that satisfy are typically defined by
starting with simple building block functions and by creating new
implicit functions using the sum, min, or max of simpler functions.
When the building blocks are polynomials in x , y, and z, the resulting
surfaces are called algebraic surfaces, and when the building blocks
are spherical Gaussian functions, the surfaces are called blobbies [8],
"soft objects", or "metaballs". The use of sums of implicit functions
allows blend surfaces to be created [24], and the use of min and max
yields the union and intersection of solid objects.
Rendering an implicit surface is often difficult. If a ray tracer
is used, intersecting a ray with an implicit surface reduces to one-dimensional
root-finding, but this can be very slow for a complex
implicit function [8]. To exploit the speed of graphics hardware,
we would prefer to render using a z-buffer algorithm. This requires
converting the implicit surface into polygons or other parametric
surfaces.
Most existing methods for polygonizing implicit surfaces sub-divide
space into a uniform grid of cubical or tetrahedral voxels,
sample the function at the vertices of the voxels, and then, for each
voxel whose vertices are not all in or all out, generate polygon(s) approximating
the surface [33,18,21,10]. This approach is often called
the marching cubes algorithm. Improvements on this algorithm use
adaptive subdivision based on curvature [9]. Unfortunately, all of
these algorithms will miss small surface features if the initial sampling
grid is too coarse, except Snyder's, which uses interval arithmetic
to guarantee that the topology of the polygonization matches
the topology of the real surface [26].
These polygonization algorithms were designed for static sur-
to polygonize a changing surface with them would require
beginning from scratch each time. The algorithm of Jevans et al. is
an exception. It re-polygonizes only those voxels that change [17].
Physically-based approaches to the polygonization of implicit
surfaces were pioneered by Figueiredo et al. [12]. One of the two
methods they describe starts with particles randomly scattered in
3-D space, subjects them to forces that pull them to the surface
(an idea proposed in [11]), and uses repulsion between particles to
distribute them uniformly over the surface. Their technique uses
penalty methods, however, which lead to stiff differential equations
whose solution is generally either slow to repel into a nice pattern, or
inaccurate at staying on the surface. Once the particles have reached
equilibrium, a polygonization is found using Delaunay triangulation.
Their work resembles ours most closely, but our simulation method
differs from theirs, and our technique supports interactive control
of surfaces and incremental sampling of changing surfaces, while
theirs does not.
2.2 Control Methods
One of the principal disadvantages of implicit modeling relative
to parametric modeling is the difficulty of controlling the shape of
an implicit surface [11]. The effect of the parameters of an implicit
surface is often non-intuitive.
With algebraic surfaces, for instance, it is hard to predict the surface
shape given its coefficients. Modeling is further complicated by
the global nature of an algebraic surface's polynomial basis func-
tions, which prevent local shape control. For these reasons and
others, piecewise algebraic surfaces have recently become popular
[25]. Piecewise algebraic surfaces are typically defined by a
weighted sum of Bernstein polynomials over a lattice of tetrahe-
dra. Least squares methods for fitting surfaces to a set of points are
available both for standard algebraic surfaces [22] and for piecewise
algebraic surfaces [1]. Pratt's algorithm can fit a surface with m parameters
to n points (n > m) in time O((n+m)m 2 ). These methods
are limited to algebraic surfaces, however.
Blobby models employ local basis functions, so they are often
more intuitive to work with than algebraic surfaces [8]. In an interactive
blobby modeling system, a user might use dials or sliders to
adjust the position and radius of each blobby center [7], but arriving
at a desired surface is a matter of guesswork, and the real time
display is typically just a wireframe, with a higher quality rendering
requiring off-line ray tracing or polygonization. Some recent work
has fit blobby models to a set of surface points, but the method is
quite slow, one example requiring days of computer time to fit 2900
control points using 1200 parameters [20]. Direct manipulation of a
blobby surface at interactive speeds has remained an open problem.
The differential methods we use to constrain the motion of particles
and surfaces are rooted in classical mechanics (see, e.g. [15]
for a discussion of mechanical constraints and constraint forces) and
are closely related to constraint methods used in physically based
modeling for computer graphics [5,2,3,32,31,4]. Allied methods
have also been used for interactive geometric modeling [30,14].
3 The Particle/Surface Constraint
In this section we derive the basic machinery that allows us to
attach moving particles to moving surfaces. First we derive a basic
constraint on particle and surface velocities that establishes, then
maintains contact as the system evolves over time. We then pose two
related problems: solve for particle velocities given time derivatives
of the surface parameters, and solve for surface derivatives given
particle velocities. Since the problem will generally be undercon-
strained, we express it as a constrained optimization.
Notation: We use boldface to denote vectors, and italics for
scalars. Subscripts denote partial differentiation. Superscript i or j
denote the i th or j th member of a collection of objects. E.g. p i is
the i th in a collection of vectors, and F x is the derivative of scalar
F with respect to vector x, hence a vector. Superscripts other than
have their usual meaning as exponents, e.g. |x - c| 2 or e -x 2 .
A dot, as in -
q, denotes a derivative with respect to time.
3.1 The Basic Constraint
We represent the moving implicit surface by F(x,
where x is position in space, and q(t) is a vector of m time-varying
shape parameters. For example, an implicit sphere could be defined
by
center c and radius r . The parameter
vector q would then be the 4-vector [c x , c y , c z , r ].
The condition that a collection of n moving particles lie on the
surface is
where p (t) is the trajectory of the i th particle. In order for this
condition to be met from some initial time t 0 onward, it suffices
that equation 1 is satisfied at t 0 , and that the time derivative -
thereafter. Since we want to manipulate velocities rather than
positions, we obtain an expression for -
F using the chain rule:
q, (2)
F , F x , and Fq evaluated at p i . By
setting -
F i to zero in equation 2, we obtain n linear constraints on the
s and on -
q. In principle, if we began with a valid state and ensured
that these conditions were met at every instant thereafter, we would
be guaranteed that the particles remained on the surface. In practice,
we might not have valid initial conditions, and numerical integration
errors would cause drift over time. We cure these problems using
a feedback term [6], setting -
is a feedback
constant. This yields the set of n linear constraint equations
3.2 Constrained Optimization
We employ these constraints in two ways: first, in order to use
particles to move the surface, we solve for -
q given the -
's. Second,
to use mutually repelling particles to sample the surface, we solve
for the -
q. In either case, we generally wish to solve underconstrained
systems. To do so we minimize a quadratic function
of -
q, subject to the constraints. The objective function we
use here is
G =2
| -
where P and Q are known desired values for -
These desired values can be used in a variety of ways. Setting P to
zero minimizes particle velocities. Setting Q to zero minimizes the
In unconstrained optimization we require that the gradient of
the objective function vanish. At a constrained minimum, we require
instead that the gradient of the objective function be a linear
combination of the gradients of the constraint functions [13]. This
condition ensures that no further local improvement can be made
without violating the constraints. In the case of a point constrained
to a surface, this condition is easily visualized: the gradient of the
objective function must lie normal to the surface, so that its orthogonal
projection onto the tangent plane vanishes. Though harder to
visualize, the idea is the same in higher dimensions.
1 Although we do not give the derivation here, a straightforward and useful
generalization is to allow error to be measured using an arbitrary symmetric
positive-definite metric tensor, e.g. ( -
In particular, it
is possible to automatically compute a sensitivity matrix, analogous to the
mass matrix in mechanics, that compensates for scale differences among the
components of Fq (see [31].)
The classical method of Lagrange multipliers [13] solves constrained
optimization problems by adding to the gradient of the objective
a linear combination of constraint gradients, with unknown
coefficients. One then solves simultaneously for the original un-
knowns, and for the coefficients. In the case of linear constraints
and a quadratic objective, this is a linear problem.
The two problems we wish to solve-obtaining -
q, and -
given -
to minimize the same objective subject to the same
constraints, differing only in regard to the knowns and unknowns.
Even so, the solutions will turn out to be quite different because of
the structure of C j 's dependencies on -
q. We next consider
each problem in turn.
3.3 Floaters
In solving for the -
's, the requirement that the gradient of the
objective be a linear combination of the constraint gradients is expressed
by
for some value of the unknown coefficients # i . The summation over
drops out because C j cannot dependon -
In addition
we require that the constraints be met, i.e. that C i
Equation 4 allows us to express the -
's in terms of the unknown
's. Substituting for -
in equation 3 gives
We may solve for each # i independently. Doing so yields
x
Substituting into equation 4 yields
x
which is the particle velocity that solves the constrained optimization
problem. Notice that in the case that the surface is not moving and
the constraints are met, so that
this reduces to
which is just the orthogonal projection of P i onto the surface's tangent
plane at p i .
3.4 Control Points
We follow the same procedure in solving for -
q, except that derivatives
of C j and G are taken with respect to -
q. The condition that the
gradient of the objective be a linear combination of the constraint
gradients is
This time, the sum does not vanish, because every C j generally
depends on -
q.
We next use equation 6 to substitute for -
q in equation 3:
Rearranging gives us the n - n matrix equation to be solved for
Note that element (i, j ) of the matrix is just the dot product F i
q .
Having solved for the # j 's, we then solve for -
using equation
q . (8)

Summary

In this section we have given the solutions to two very closely
related problems:
. Given the instantaneous surface motion -
q, solve for particle
velocities -
deviation from desired velocities
subject to the constraint that the particles stay on the sur-
Each particle's constrained velocity may be computed
independently.
. Given the particle velocities -
solve for the implicit function
time derivative -
q that minimizes deviation from a desired time
derivative Q, again, subject to the constraint that the particles
must remain on the surface. Calculating -
q entails the solution
of an n - n linear system, where n is the number of particles.
We combine these methods by maintaining two populations of
particles: control points and floaters. Control points are moved
explicitly by the user, and -
q is calculated to make the surface follow
them. In contrast, floaters' velocities are calculated to make them
follow the surface, once -
q has been computed.
Adaptive Sampling
In this section we address the problem of sampling implicit sur-
faces, building on the floater mechanism that we presented in the
previous section. Good sampling is a requirement both for quick
rendering and for the evaluation of integrals such as surface area or
volume.
Our primary goal is to obtain sampling distributions that are either
(a) uniform, with user-specified density, (b) or non-uniform, with
density based on local criteria such as surface curvature. We wish
to reach the specified distribution quickly from a few seed points
(ideally, only one per connected component) and to maintain a good
distribution as the surface moves and deforms. To support interactive
sculpting, we must be able to update at least a few hundred sample
points at 10Hz or better. Additional goals are that the particles should
move as little as possible in response to surface motion, and that only
basic and generic information about the function F be required. It
should not be necessary to supply a surface parameterization.
The starting point for our approach is the idea, introduced by
Turk [28] and by Figueiredo et al. [12], that particles can be made
to spread out to uniform density by local repulsion, relying on the
finiteness of the surface to limit growth. Simple repulsion can do a
good job at ironing out local irregularities given a reasonably good
initial sampling (as in Turk's application to resampling of a polygon
mesh) but is extremely slow to converge if the initial sampling is
irregular at large scale, and fails completely to track surface motions
and deformations.
After describing our basic repulsion scheme, we introduce the
idea of global fissioning: we start the sampling process with a very
small number of particles but a very large radius of interaction, coming
close to equilibrium in just a few iterations. We then fission each
particle, imposing random displacements that are smaller than the
interaction radius. At the same time, we scale the interaction radius
to a smaller value. We now have a new starting point, locally
irregular but with nearly uniform large-scale structure. A few iterations
suffice to smooth out the small irregularities and reach a new
equilibrium. The scaling and fissioning process is repeated until the
target sampling density is reached.
Global fissioning still fails to handle surface motion: should new
nonuniformities be introduced after the fissioning process termi-
nates, the system suffers all of the shortcomings of simple fixed-
scale repulsion. So, for example, the sudden introduction of a bulge
in the surface can create a gaping hole in the sampling pattern that
will be repaired extremely slowly, if at all. Intuitively, we would
like particles at the edge of such voids to "feel" the reduction of
density, expand their radii of interaction to quickly fill the hole, then
begin fissioning to restore full density. On the other hand, if density
becomes too high, we would like particles to die off until the desired
density is restored. We will conclude the section by describing
a fast and robust adaptive repulsion scheme that provides just this
behavior, meeting all of our goals.
4.1 Simple Repulsion
As a windowed density measure, we employ a simple Gaussian
energy function based on distances between particles in 3-D. We
define the energy of particle i due to particle j to be:
is the vector between particles, # is a global
repulsion amplitude parameter, and # , called the global repulsion
radius, is the standard deviation of the Gaussian. The repulsion
radius controls the range of the repulsion "force." Note that
The energy of particle i in its current position is defined as:
Ultimately, we would like to reach the global minimum of each
by varying the particle positions on the surface. Finding the
global minimum is impractical, but we can find a local minimum by
gradient descent: each particle moves in the direction that reduces its
energy fastest. We therefore choose each particle's desired velocity
to be negatively proportional to the gradient of energy with respect
to its position:
The formulas for energy and desired velocity have been carefully
chosen here so that "energy" is unitless, while desired velocity is
proportional to distance. This guarantees that the sampling pattern
computed by this simple repulsion method scales with a surface.
If desired particle velocities are set in this way, and constrained
particle velocities are computed with equation 5, particles repel, but
their behavior is highly dependent on the parameter # . The slope of a
c
Gaussian peaks at distances of -# and it is near zero at much smaller
or much greater distances. When the distance between particles is
not between .03# and 3# , for instance, the repulsion is below 7%
of its peak. If # is chosen too small then particles will (nearly) stop
spreading when their separation is about 3# , and if # is chosen too
big then distant particles will repel more than nearby ones, and the
resulting sampling pattern will be poor. The best value for # is about
.3 p (surface area)/(number of particles).
4.2 Global Fissioning
If a surface is seeded with several floater particles, and an initial
value of # can be found that causes these particles to disperse, then
the sampling can be repeatedly refined by allowing the particles
to reach equilibrium, then simultaneously fissioning each particle
into two, giving the new particles a small random displacement, and
simultaneously dividing # by # 2. The particles are considered to be
at equilibrium when their net forces, and hence their speeds, get low.
With this global fissioning scheme, early generations will spread out
sparsely, and succeeding generations will fill in more densely.
Simple repulsion with global fissioning is acceptable for maintaining
a good distribution on a very slowly changing surface, but the
population is always a power of two, and particles do not redistribute
quickly in response to rapid surface changes. Global fissioning fails
to adapt to changes in a surface adequately, as mentioned earlier.
4.3 Adaptive Repulsion
To develop a more adaptive repulsion scheme, we employ an analogy
to a population of organisms distributing itself uniformly across
an area. Specifically, imagine a population of pioneers spreading
West and colonizing America. In order to settle the entire country
as quickly as possible, a good rule is for each male-female pair to
spread out as much as possible away from their neighbors, until the
encroachment on them is roughly equal in all directions, and only
then to homestead and have children. If the encroachment from
neighbors is low, then each pair can claim more land (be greedier),
but when neighbors are pressing in, each pair must relinquish land.
Early pioneers travel great distances and claim huge tracts of land,
while later generations move less and divide up successively smaller
shares until the desired density is achieved.
These ideas can be applied to particle behavior. To achieve uniform
densities quickly, and maintain them as the surface moves or
deforms, we will allow each particle to have its own repulsion radius
and to decide independently when it should fission or die.
A particle's radius should grow when all of the forces on it are small
and it should shrink when the forces on it are big. For a particle
near equilibrium, birth and death occur when the density is too low
or too high, respectively. We now quantify these principles.
Similar to the simple repulsion scheme, we define the energy of
particle i due to particle j as:
Note that the global parameter # has been replaced by the local
parameter
The energy at particle i is defined as:
The repulsion force and desired velocity is again proportional to the
gradient of energy with respect to position:
The time-varying repulsion radii will be controlled differentially.
We want the radius to grow when the energy is too low and to
shrink when the energy is too high. This can be done indirectly by
controlling the energies.
As stated earlier, our energy measure is scale-invariant. That is,
if all surfaces and samples are scaled (p i and # ), the
constant. Therefore, to ensure that neighboring particles repel each
other, we can simply drive all of their energies to a global desired
energy level, -
To arrive at a value for -
, we consider an ideal
hexagonal close-packing, which is the best uniform sampling pattern
for a planar surface. In this configuration, all # i should be equal,
and the distance between nearest neighbors should be roughly 2#
to guarantee strong repulsion forces. Since each particle has six
nearest neighbors in this configuration, the desired energy should be
roughly -
The portion of a particle's repulsion energy that is directly affected
by a change in its own repulsion radius is:
To keep D i near the desired value, we use the linear feedback equation

E)
where # is the feedback constant.
The change to the repulsion radius of a particle that will yield
this change in energy can be derived with the chain rule: -
neglecting the latter terms, thus:
The rule aboveworks fine for particles that are exerting some force on
their neighbors, but it causes infinite radius change when a particle
is alone in a sparsely sampled region of a surface (or is the first
particle), where D i
In such cases we want the radius
to grow, but not catastrophically, so we modify equation 11:
for some #. The change in energy with respect to a change in radius
is:
Using equations 9, 12, 10, and 13 to control particle positions and
repulsion radii will do a good job of moving particles into sparse
regions quickly, but their radii might become very large, and hence
the density might remain too low.
4.4 Adaptive Fission/Death
To achieve uniform density it is necessary that large-radius particles
fission. Likewise, particles that are overcrowded should be
considered for death.
We use the following criteria to control birth and death of parti-
cles: A particle is fissioned iff:

Figure

1: This sequence illustrates the adaptive repulsion and fis-
sioning mechanism. The topmost image shows a deliberately poor
sampling of a blobby cylinder produced using simple repulsion: the
cylinder was rapidly stretched, leaving the sample points behind.
The remaining images, from top to bottom, show the recovery of
good sampling when adaptive repulsion is enabled. The particles
at the frontier increase their radii of repulsion, rapidly filling the
voids. As the particles slow down, they fission, restoring the desired
sampling density. This process takes about four seconds on an
SGI Crimson.
. the particle is near equilibrium, | -
| < # , and
. either the particle's repulsion radius is huge (# max ), or
it is adequately energized and its radius is above the desired
radius (D i > # -
Fission splits a single particle in two. The two particles are given
initial radii of # i /
# 2 and a desired velocity that is a random direction
scaled by a fraction of # i . A particle dies iff:
. the particle is near equilibrium, | -
| < # , and
. the particle's repulsion radius is too small, # i < # -
# , and
. the following biased randomized test succeeds: R > # i
where R is a uniform random number between 0 and 1.
The death criteria are made stochastic to prevent mass suicide in
overcrowded regions.
This combination of adaptive repulsion, fissioning, and death
is much more responsive to changes in the surface shape than the
simple repulsion scheme.
5 Implementation and Results
The techniques described above have been implemented in about
3700 lines of C++ code. Particular implicit function classes are
derived from a generic implicit function base class. Adding a new
implicit function to the system is easy, requiring only the implementation
of functions F , F x , Fq , and bounding box. Each of these
except Fq is standard in any system employing implicit functions.
For example, we define the blobby sphere implicit function to be
the sum of Gaussians of the distance to each of k center points [8].
The parameter vector q consists of 4k parameters: a bias b plus
four parameters for each sphere (a center 3-vector c i and standard
deviation s ). Thus,
If we define
then the functions needed by the system are
where
If we assume that g i beyond a radius of 3s , then a conservative
bounding box for blobby spheres is the bounding box of
non-blobby spheres with centers c i and radii 3s i .
We have also implemented spheres and blobby cylinders. A
blobby cylinder function is defined to be the sum of Gaussians of
the distance to each of several line segments. A system of k blobby
cylinders has 7k parameters: a bias plus seven parameters for
each cylinder (two endpoints and a standard deviation).
It is often useful to freeze some of these parameters to a fixed
value so that they will not be modified during interaction. This is
c
done simply by leaving them out of the q and Fq vectors. To get
blobs of equal radii, for instance, one would omit all s .
The system starts up with a single floater positioned arbitrarily
in the bounding box of the surface and then begins the physical
simulation by repeating the following differential step:
. The user interface sets desired control point velocities P i . Stationary
control points of course have zero desired velocity,
while control points being dragged by the user have desired
velocities that are calculated as a function of cursor position.
. Set Q, the desired values for the time derivatives of the surface
parameters. These are typically set to zero to minimize parametric
change in the surface, but they could also be calculated
to attract the surface toward a default shape.
. Compute the actual surface parameter changes, -
q, as constrained
by the control point velocities, using equations 7 and
8.
. Compute repulsion forces between floaters to set their desired
velocities using equation 9.
. Compute actual floater velocities, as constrained by the
already-computed surface time derivatives, using equation 5.
(When the gradient F x is near zero, however, the surface is locally
ill-defined, and it is best to leave such floaters motionless,
. Compute the change to floater repulsion radii, -
using equations
12, 10, and 13.
. Update the positions of the control points and floaters using
Euler's method, that is: p
similar formulas to update the surface parameters q from -
q,
and the floater repulsion radii # i from -
# .
. Test each floater for possible fission/death.
. Redisplay the floaters and control points as disks tangent to
the surface, with normal given by F x and (for floaters) radius
proportional to # i .
Using the mouse, the user can pick a control point and move it in
a plane perpendicular to the view direction. Pulling control point i
sets the desired control point velocity P i . Since the velocities feed
into the constrained optimization solution, which in turn feeds into
a numerical differential equation solver, some care must be taken to
ensure that control point motions are reasonably smooth and well
behaved, which they might not be if positions were set directly by
polling the pointing device. A simple solution which works well is to
make the velocity of the dragged particle proportional to the vector
from the point to the 3-D cursor position. This in effect provides
spring coupling between the cursor and the control point. Although
the control point can lag behind the cursor as a result, performance
is brisk enough that the lag is barely noticeable. Similar dragging
schemes are described in [14,31]. The user can also create and delete
control points and adjust the desired repulsion radius -
# through a
slider.
The matrix in equation 7 is symmetric and in general positive
definite. It thus lends itself to solution by Cholesky decomposition
[23], which is easy to implement, robust and efficient. However, the
matrix can become singular if inconsistent or redundant constraints
are applied, that is if the number of constraints exceeds m, or if
some of the F i
q 's are linear combinations of others. While the
former condition is easy to detect by counting, the latter is not.
The problem of singularities can be circumvented by using a least-squares
technique, or singular value decomposition [23].
The system is fast enough to run at interactive rates. Let m be
the number of degrees of freedom of the implicit surface, let n be
the number of control points, and let r be the number of floaters.
The most expensive parts of the algorithm are the computation of
the n -n matrix of equation 7, which has cost O(mn 2 ), the solution
of the linear system, which has cost O(n 3 ), the computation of repulsion
forces between all pairs of floaters, which currently has cost
O(r 2 ), and the display of the floaters, which has cost O(r) (with a
large constant). Our current system does not handle overconstrained
thus the total asymptotic cost of the algorithm
is O(mn 2
We have run simulations as complex as
Above floaters, the O(r 2 ) repulsion cost has dominated,
but this could easily be optimized using spatial data structures. For
smaller numbers of floaters (r < 150), our system runs at interactive
rates (10 Hz or faster on a Silicon Graphics workstation with 100
MHz processor).
The following parameter settings are recommended (where d is
surface diameter):
coefficients to keep
particles from drifting off sur-
face, and keep particles ener-
gized, respectively
desired energy
prevent divide-by-zero
desired repulsion radius (user-
controllable)
that this changes over time)
speed (multiple of
.2 fraction of -
fraction of -
# , for death
Most of these parameters can be set once and forgotten. The only
parameter that a user would typically need to control is the desired
repulsion radius, -
# .
Overall, the method meets our goals, it is fast, and it has proven
to be very robust. It has recovered from even violent user interaction
causing very rapid shape change. The adaptive sampling, fission,
and death techniques seem to be well tuned and to work well to-
gether, as we have not seen the system oscillate, diverge, or die with
the current parameter settings. During interaction, -
# is the only
parameter that needs to be varied.
Another result of this work is that we have discovered that implicit
surfaces are slippery: when you attempt to move them using control
points they often slip out of your grasp.
6 Conclusions
In this paper we have presented a new particle-based method for
sampling and control of implicit surfaces. It is capable of supporting
real-time rendering and direct manipulation of surfaces. Our
control method is not limited to algebraic surfaces as many previous
techniques are; it allows fast control of general procedural implicit
functions through control points on the surface. We have presented
a dynamic sampling and rendering method for implicit surfaces that
samples a changing surface more quickly than existing methods.
The use of constraint methods allows particles to follow the surface
as it changes, and to do this more rapidly and accurately than with
penalty methods. Our algorithms for adaptive repulsion, fission, and
death of particles are capable of generating good sampling patterns
much more quickly than earlier repulsion schemes, and they sample
the surface well even during rapid shape changes.

Figure

2: This sequence illustrates the construction of a shape composed
of blobby cylinders. The shape was created by direct manipulation
of control points using the mouse. In the topmost image, all
three cylinder primitives are superimposed. Each subsequent image
represents the result of a single mouse motion.
There are a number of directions for future research.
We intend to investigate other uses for the samplings we obtain.
One of these is the calculation of surface integrals for area, volume,
or surface fairness measures such as those described in [19,30].
Another is the creation of polygon meshes.
To polygonize a surface within the framework presented here it
is necessary to infer topology from the sample points. This is more
difficult than finding a polygonization from a set of samples on a
grid in 3-D, as in marching cubes algorithms, where an approximate
topology is suggested by the signs of the samples and by the topology
of the grid itself. Delaunay triangulation in 2-D or 3-D is one
possible way to extract topology [12,27]. A more robust alternative
would employ Lipschitz conditions and interval arithmetic [26]. To
preserve the basic advantages of our method, we would require a
polygonization algorithm that allows efficient dynamic updates as
the surface changes.
Although we developed it to sample implicit surfaces, our adaptive
repulsion scheme can be applied to meshing or sampling of
parametric surfaces as well: each floater would be defined by its
position in the surface's 2-D parameter space, rather than position
in 3-D space.
Several performance and numerical issues remain to be addressed.
As we tackle more complex models, we could exploit sparsity in
F's dependence on q. Notably, with local bases such as blobby
models, the dependence of F on faraway elements is negligible.
An additional numerical issue is the handling of singular constraint
matrices, due to overdetermined or dependent constraints. Excellent
results can be obtained using least-squares techniques.
An additional area of investigation is the use of local criteria,
notably surface curvature, to control sampling density. Surface curvature
can be measured directly, at the cost of taking additional
derivatives of F . Since this places a considerable extra burden on
the implementor of implicit primitives, an alternative is to estimate
curvature at each floater based on positions and normals of nearby
points. Having established a desired density at each point, based
on curvature or any other criterion, relatively simple modifications
to the adaptive repulsion scheme will yield the desired nonuniform
density. Another possible density criterion is the user's focus of
interest, e.g. the neighborhood of a control point being dragged.
Finally, there is room for considerable further work in interactive
sculpting of implicit surfaces. Dragging one control point at a time
can be somewhat limiting given the slippery behavior of the surface.
However, the basic control-point machinery developed here could be
used to build more complex sculpting tools that influence multiple
surface points in coordinated ways.

Acknowledgements

The authors wish to thank Scott Draves and Sebastian Grassia for
their contributions to this work. This research was supported in part
by a Science andTechnology Center Grant from the National Science
Foundation, #BIR-8920118, by an NSF High Performance Computing
and Communications Grant, #BIR-9217091, by the Engineering
Design Research Center, an NSF Engineering Research Center at
Carnegie Mellon University, by Apple Computer, Inc, and by an
equipment grant from Silicon Graphics, Inc. The second author
was supported by NSF Young Investigator Award #CCR-9357763.



--R


Analytical methods for dynamic simulation of non-penetrating rigid bodies
Curved surfaces and coherence for non-penetrating rigid body simulation
Dynamic simulation of non-penetrating flexible bodies
A modeling system based on dynamic constaints.
Stabilization of constraints and integrals of motion in dynamical systems.
Practical uses for implicit surfaces in an- imation
A generalization of algebraic surface drawing.
Polygonization of implicit surfaces.
An implicit surface polygonizer.
Interactive techniques for implicit modeling.
Luiz Henrique de Figueiredo
Practical Optimization.

Classical Mechanics.
Mesh optimization.
Speeding up 3-D animation for simulation
Marching cubes: A high resolution 3D surface reconstruction algorithm.
Functional minimization for fair surface design.
"blobby model"
An evaluation of implicit surface tilers.
Direct least-squares fitting of algebraic sur- faces

A constructive geometry for computer graphics.
Piecewise algebraic surface patches.
Generative Modeling for Computer Graphics and CAD.
Surface modeling with oriented particle systems.
Generating textures on arbitrary surfaces using reaction-diffusion

Variational surface mod- eling
Interactive dynamics.
AndrewWitkin and
Data structure for soft objects.
--TR
Direct least-squares fitting of algebraic surfaces
Marching cubes: A high resolution 3D surface construction algorithm
Polygonization of implicit surfaces
Analytical methods for dynamic simulation of non-penetrating rigid bodies
Interactive dynamics
Interactive techniques for implicit modeling
Curved surfaces and coherence for non-penetrating rigid body simulation
Volumetric shape description of range data using MYAMPERSANDldquo;Blobby ModelMYAMPERSANDrdquo;
Generating textures on arbitrary surfaces using reaction-diffusion
Generative modeling for computer graphics and CAD
Re-tiling polygonal surfaces
Variational surface modeling
Functional optimization for fair surface design
Surface modeling with oriented particle systems
Dynamic simulation of non-penetrating flexible bodies
Through-the-lens camera control
Physically-based methods for polygonization of implicit surfaces
Higher-order interpolation and least-squares approximation using implicit algebraic surfaces
Mesh optimization
An implicit surface polygonizer
A Generalization of Algebraic Surface Drawing
Fast animation and control of nonrigid structures
An Evaluation of Implicit Surface Tilers

--CTR
Tasso Karkanis , A. James Stewart, Curvature-Dependent Triangulation of Implicit Surfaces, IEEE Computer Graphics and Applications, v.21 n.2, p.60-69, March 2001
Hans Khling Pedersen, Decorating implicit surfaces, Proceedings of the 22nd annual conference on Computer graphics and interactive techniques, p.291-300, September 1995
Gary Yngve , Greg Turk, Robust Creation of Implicit Surfaces from Polygonal Meshes, IEEE Transactions on Visualization and Computer Graphics, v.8 n.4, p.346-359, October 2002
Yoshitomo Jo , Masafumi Oka , Akinori Kimura , Kyoko Hasegawa , Ayumu Saitoh , Susumu Nakata , Akihiro Shibata , Satoshi Tanaka, Technical Section: Stochastic visualization of intersection curves of implicit surfaces, Computers and Graphics, v.31 n.2, p.230-242, April, 2007
Patricia Crossno , Edward Angel, Isosurface extraction using particle systems, Proceedings of the 8th conference on Visualization '97, p.495-ff., October 18-24, 1997, Phoenix, Arizona, United States
Marek Teichmann , Michael Capps, Surface reconstruction with anisotropic density-scaled alpha shapes, Proceedings of the conference on Visualization '98, p.67-72, October 18-23, 1998, Research Triangle Park, North Carolina, United States
Jerzy Karczmarczuk, Teaching of image synthesis in functional style, Proceedings of the 2005 workshop on Functional and declarative programming in education, September 25-25, 2005, Tallinn, Estonia
Nina Amenta , Marshall Bern , Manolis Kamvysselis, A new Voronoi-based surface reconstruction algorithm, Proceedings of the 25th annual conference on Computer graphics and interactive techniques, p.415-421, July 1998
Hans Khling Pedersen, A framework for interactive texturing on curved surfaces, Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, p.295-302, August 1996
John C. Hart, Using the CW-complex to represent the topological structure of implicit surfaces and solids, ACM SIGGRAPH 2005 Courses, July 31-August
Max Froumentin , Eric Varlet, Dynamic implicit surface tesselation, Proceedings of the ACM symposium on Virtual reality software and technology, p.79-86, September 1997, Lausanne, Switzerland
John C. Hart , Ed Bachta , Wojciech Jarosz , Terry Fleury, Using particles to sample and control more complex implicit surfaces, ACM SIGGRAPH 2005 Courses, July 31-August
Patricia Crossno , Edward Angel, Visual debugging of visualization software: a case study for particle systems, Proceedings of the conference on Visualization '99: celebrating ten years, p.417-420, October 1999, San Francisco, California, United States
Bryan S. Morse , Terry S. Yoo , Penny Rheingans , David T. Chen , K. R. Subramanian, Interpolating implicit surfaces from scattered surface data using compactly supported radial basis functions, ACM SIGGRAPH 2005 Courses, July 31-August
Barbara J. Meier, Painterly rendering for animation, Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, p.477-484, August 1996
Jing Hua , Hong Qin, Haptics-Based Dynamic Implicit Solid Modeling, IEEE Transactions on Visualization and Computer Graphics, v.10 n.5, p.574-586, September 2004
Matthias Mller , Simon Schirm , Matthias Teschner, Interactive blood simulation for virtual surgery based on smoothed particle hydrodynamics, Technology and Health Care, v.12 n.1, p.25-31, February 2004
Frederic Triquet , Laurent Grisoni , Philippe Meseure , Christophe Chaillou, Realtime visualization of implicit objects with contact control, Proceedings of the 1st international conference on Computer graphics and interactive techniques in Australasia and South East Asia, February 11-14, 2003, Melbourne, Australia
Julie Dorsey , Hans Khling Pedersen , Pat Hanrahan, Flow and changes in appearance, ACM SIGGRAPH 2005 Courses, July 31-August
Fausto Bernardini , Joshua Mittleman , Holly Rushmeier , Cludio Silva , Gabriel Taubin, The Ball-Pivoting Algorithm for Surface Reconstruction, IEEE Transactions on Visualization and Computer Graphics, v.5 n.4, p.349-359, October 1999
Julie Dorsey , Hans Khling Pedersen , Pat Hanrahan, Flow and changes in appearance, Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, p.411-420, August 1996
Tasso Karkanis , A. James Stewart, Curvature-Dependent Triangulation of Implicit Surfaces, IEEE Computer Graphics and Applications, v.21 n.2, p.60-69, March 2001
Julie Dorsey , Hans Kohling Pederseny , Pat Hanrahan, Flow and changes in appearance, ACM SIGGRAPH 2006 Courses, July 30-August 03, 2006, Boston, Massachusetts
Greg Turk , David Banks, Image-guided streamline placement, Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, p.453-460, August 1996
Matthias Mller , David Charypar , Markus Gross, Particle-based fluid simulation for interactive applications, Proceedings of the ACM SIGGRAPH/Eurographics symposium on Computer animation, July 26-27, 2003, San Diego, California
Hans Pedersen , Karan Singh, Organic labyrinths and mazes, Proceedings of the 4th international symposium on Non-photorealistic animation and rendering, June 05-07, 2006, Annecy, France
Siu-Wing Cheng , Tamal K. Dey , Edgar A. Ramos , Tathagata Ray, Sampling and meshing a surface with guaranteed topology and geometry, Proceedings of the twentieth annual symposium on Computational geometry, June 08-11, 2004, Brooklyn, New York, USA
Lee Markosian , Jonathan M. Cohen , Thomas Crulli , John Hughes, Skin: a constructive approach to modeling free-form shapes, Proceedings of the 26th annual conference on Computer graphics and interactive techniques, p.393-400, July 1999
Aravind Kalaiah , Amitabh Varshney, Modeling and Rendering of Points with Local Geometry, IEEE Transactions on Visualization and Computer Graphics, v.9 n.1, p.30-42, January
J. D. Boissonnat , S. Oudot, Provably good surface sampling and approximation, Proceedings of the Eurographics/ACM SIGGRAPH symposium on Geometry processing, June 23-25, 2003, Aachen, Germany
D. Stahl , N. Ezquerra , G. Turk, Bag-of-particles as a deformable model, Proceedings of the symposium on Data Visualisation 2002, May 27-29, 2002, Barcelona, Spain
Robert Blanding , Cole Brooking , Mark Ganter , Duane Storti, A skeletal-based solid editor, Proceedings of the fifth ACM symposium on Solid modeling and applications, p.141-150, June 08-11, 1999, Ann Arbor, Michigan, United States
Greg Turk, Texture synthesis on surfaces, Proceedings of the 28th annual conference on Computer graphics and interactive techniques, p.347-354, August 2001
B. Crespin, Dynamic triangulation of variational implicit surfaces using incremental Delaunay tetrahedralization, Proceedings of the 2002 IEEE symposium on Volume visualization and graphics, October 28-29, 2002, Boston, Massachusetts
Gershon Elber, Line Art Illustrations of Parametric and Implicit Forms, IEEE Transactions on Visualization and Computer Graphics, v.4 n.1, p.71-81, January 1998
Antoine Bouthors , Matthieu Nesme, Twinned meshes for dynamic triangulation of implicit surfaces, Proceedings of Graphics Interface 2007, May 28-30, 2007, Montreal, Canada
Mark Pauly , Richard Keiser , Leif P. Kobbelt , Markus Gross, Shape modeling with point-sampled geometry, ACM Transactions on Graphics (TOG), v.22 n.3, July
Mark Pauly , Markus Gross , Leif P. Kobbelt, Efficient simplification of point-sampled surfaces, Proceedings of the conference on Visualization '02, October 27-November 01, 2002, Boston, Massachusetts
Ron J. Balsys , Kevin G. Suffern, Point based rendering of non-manifold surfaces with contours, Proceedings of the 2nd international conference on Computer graphics and interactive techniques in Australasia and South East Asia, June 15-18, 2004, Singapore
Greg Turk , James F. O'Brien, Modelling with implicit surfaces that interpolate, ACM SIGGRAPH 2005 Courses, July 31-August
Wen Y. Su , John C. Hart, A programmable particle system framework for shape modeling, ACM SIGGRAPH 2005 Courses, July 31-August
Greg Turk , James F. O'brien, Modelling with implicit surfaces that interpolate, ACM Transactions on Graphics (TOG), v.21 n.4, p.855-873, October 2002
Yu-Kun Lai , Qian-Yi Zhou , Shi-Min Hu , Ralph R. Martin, Feature sensitive mesh segmentation, Proceedings of the 2006 ACM symposium on Solid and physical modeling, June 06-08, 2006, Cardiff, Wales, United Kingdom
Victor J. Milenkovic, Position-based physics: simulating the motion of many highly interacting spheres and polyhedra, Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, p.129-136, August 1996
Bryan S. Morse , Weiming Liu , Terry S. Yoo , Kalpathi Subramanian, Active contours using a constraint-based implicit representation, ACM SIGGRAPH 2005 Courses, July 31-August
Jatin Chhugani , Subodh Kumar, Budget sampling of parametric surface patches, Proceedings of the symposium on Interactive 3D graphics, April 27-30, 2003, Monterey, California
Barton T. Stander , John C. Hart, Guaranteeing the topology of an implicit surface polygonization for interactive modeling, ACM SIGGRAPH 2005 Courses, July 31-August
Barton T. Stander , John C. Hart, Guaranteeing the topology of an implicit surface polygonization for interactive modeling, Proceedings of the 24th annual conference on Computer graphics and interactive techniques, p.279-286, August 1997
M. Mller , R. Keiser , A. Nealen , M. Pauly , M. Gross , M. Alexa, Point based animation of elastic, plastic and melting objects, Proceedings of the 2004 ACM SIGGRAPH/Eurographics symposium on Computer animation, August 27-29, 2004, Grenoble, France
Kurt W. Fleischer , David H. Laidlaw , Bena L. Currin , Alan H. Barr, Cellular texture generation, Proceedings of the 22nd annual conference on Computer graphics and interactive techniques, p.239-248, September 1995
Ireneusz Tobor , Patrick Reuter , Christophe Schlick, Reconstructing multi-scale variational partition of unity implicit surfaces with attributes, Graphical Models, v.68 n.1, p.25-41, January 2006
Haixia Du, Interactive shape design using volumetric implicit PDEs, Proceedings of the eighth ACM symposium on Solid modeling and applications, June 16-20, 2003, Seattle, Washington, USA
Patrick Reuter , Ireneusz Tobor , Christophe Schlick , Sbastien Dedieu, Point-based modelling and rendering using radial basis functions, Proceedings of the 1st international conference on Computer graphics and interactive techniques in Australasia and South East Asia, February 11-14, 2003, Melbourne, Australia
Welch , Andrew Witkin, Free-form shape design using triangulated surfaces, Proceedings of the 21st annual conference on Computer graphics and interactive techniques, p.247-256, July 1994
Nathan Bell , Yizhou Yu , Peter J. Mucha, Particle-based simulation of granular materials, Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, July 29-31, 2005, Los Angeles, California
Haixia Du , Hong Qin, A shape design system using volumetric implicit PDEs, ACM SIGGRAPH 2005 Courses, July 31-August
Marie-Paule Cani-Gascuel , Mathieu Desbrun, Animation of Deformable Models Using Implicit Surfaces, IEEE Transactions on Visualization and Computer Graphics, v.3 n.1, p.39-50, January 1997
Marc Alexa , Johannes Behr , Daniel Cohen-Or , Shachar Fleishman , David Levin , Claudio T. Silva, Point set surfaces, Proceedings of the conference on Visualization '01, October 21-26, 2001, San Diego, California
Marc Alexa , Johannes Behr , Daniel Cohen-Or , Shachar Fleishman , David Levin , Claudio T. Silva, Computing and Rendering Point Set Surfaces, IEEE Transactions on Visualization and Computer Graphics, v.9 n.1, p.3-15, January
Ronald P. Fedkiw , Guillermo Sapiro , Chi-Wang Shu, Shock capturing, level sets, and PDE based methods in computer vision and image processing: a review of Osher's contributions, Journal of Computational Physics, v.185 n.2, p.309-341, March
constrained optimization;physically based modeling;interaction;adaptive sampling