Dynamic Subdivision Surfaces For Natural Terrain Modeling
C ,1 !,,, ,,l 1, ,; Il, i I
Dept. of CISE
University of Florida
Gainesville, FL 32611
Hong Qin
Dept. of CS
State University of "' York
Stony Brook, NY 11794
Baba C. Vemuri
Dept. of CISE
University of Florida
Gainesville, FL 32611
Abstract
We propose a novel technique for r,.'!!,.. '. ".A re
alistic terrain models from elevation data by using a
new I, ...... ', finite element methodbased subdivision
,., f,,,, model for .,. f',,, reconstruction and a variant
of the successive random addition method for fractal
., .. ...i.. 1!., subdivision *,f...' m odel is
based on the '1.,ii., [1,, subdivision scheme, a popular
technique for generating smooth ,..,.' .... of arbitrary
",1."l~;.,,. and is embedded in a i,'. v based modeling
paradigm. I !, initialized model I, f,*! ... under the in
fluence of the .'!!,. ,', ,I data forces, at the same time
Gaussian noise is added at 1I'T[7 T, . levels of subdivi
sion so that a naturallooking ,. f .. is reconstructed.
We introduce a noise addition technique in the butter
fly subdivision scheme to generate fractals and describe
a novel technique for locally parameterizing the sur
face generated by the subdivision algorithm. We also
introduce ,.1.',.. 1 quantities required to develop the
1,,..... ', subdivision ,.',f.i.. model and derive the gov
erning ,., '....... ,'[T, .'l equation using Lagrangian
mechanics and the finite element method. Our exper
iments with real terrain data demonstrate the fT .....
of the proposed technique for modeling natural terrains
from elevation data.
1 Introduction
Natural terrain modeling from elevation data is an
important application area in computer vision and
graphics. Roughness is an essential characteristic of
the natural terrains and hence traditional surface re
construction methods using smoothness constraints do
not yield the desired results. Fournier et al. [1] first
proposed a random displacement technique for ii
thesizing fractal surfaces which was later modified by
Saupe [2] in his successive random addition scheme of
generating fractals. Yokoya et al. [3] improved these
schemes by adding data constraints. Szeliski and Ter
zopoulos [4] proposed constrained fractal surfaces us
ing a Gibbs sampler algorithm which was later im
proved by Vemuri et al. [5]. Arakawa and Krotkov
[6] refined the original Gibbs sampler technique by re
defining the temperature (control) parameter to ob
tain better control of roughness in the fitted (con
strained) fractal surface. However, the (CPU) execu
tion times reported in their work are very high, thus
making their scheme unattractive for i! i1ii applica
tions.
All the techniques for natural terrain modeling
mentioned above usually needs a grid of very large
size to model realistic terrains, especially with irregu
larly spaced data. Riimelin [7] developed an interest
ing fractal interpolation algorithm which can generate
interpolating fractal surfaces for regularly or irregu
larly spaced data. However, this scheme is not com
putationally practical for large problems as pointed
out in [5].
In this paper, we propose a novel technique for nat
ural terrain modeling using a new dynamic subdivi
sion surface model in . 1 ,i. L, 1I ii with a variant of
the successive random addition method for fractal sur
face ill,,  Subdivision surfaces are widely used in
computer graphics for modeling surfaces of (known)
arbitrary geometry/I', 1,I _. To model a surface of
a known l' I ," 1' first an arbitrary polygonal mesh
of the same genus is chosen as the initial mesh a.k.a.
control mesh. Next, a fixed set of rules are applied re
cursively to refine the initial mesh. After each step of
subdivision, a finer I" .1 _. .1 mesh with more vertices,
edges and faces is obtained. In the limit, a smooth
surface of the same genus as of the initial mesh is
obtained. It 11 i be noted that the limit ,..f.'.. is
described by the same degrees of freedom as in the ini
tial mesh, thereby offering a very compact represen
tation of a potentially complicated shape. In [8], we
embedded a specific subdivision scheme in a I'1! 
based modeling paradigm and the resulting dynamic
subdivision surface model was used to recover com
plex smooth shapes of arbitrary 1 .1 **1 ,_ from large
range and volume data sets using very few degrees of
freedom.
The dynamic subdivision surface model proposed
in [8] relies on the nature of the specific subdivision
rules used. The approach taken in this paper is much
more general in the sense that it can be used with
...., r I1 of subdivision schemes. However, we choose
the butterfly scheme [9] to demonstrate the concept.
This scheme has the added advantage of an interpola
tory subdivision scheme where the limit surface inter
polates the initial mesh whereas in [8], the limit sur
face only approximates the initial mesh. Once we em
bed the chosen subdivision surface model into i'1! i 
based paradigm, the initialized model deforms under
the influence of synthesized forces to fit the elevation
data via the principle of !i. i minimization. The in
terpolatory feature of the butterfly subdivision scheme
reduces the computational overhead involved in prop
agating the external forces from the limit surface to
the initial control mesh. Recalling the fact that the
limit surface in ;io subdivision scheme is a function
of the degrees of freedom of the initial mesh, once an
approximate surface is recovered from the data, the
model adopts a new initial mesh which is obtained via
a subdivision of the original initial mesh. Note that
this new initial mesh and the original initial mesh have
the same limit surface, but the new initial mesh has
more degrees of freedom thereby facilitating the recov
ery of the local features in the surface. This process is
continued till a prescribed error criteria for fitting the
data points is achieved.
The limit surface of .;i subdivision scheme is
smooth, whereas the natural terrains are not. The
dynamic subdivision surface model, without i~! mod
ifications, recovers a smooth surface from elevation
data. To recover a natural terrain, we perturb the
new vertex positions of the dynamic model at vari
ous levels of subdivision in a similar fashion as in the
successive random addition method [2]. The variance
of the added Gaussian noise at different subdivision
levels controls the roughness of the final reconstructed
natural terrain. The proposed model can  ill. i ..
natural terrains from elevation data using very few
degrees of freedom. The effectiveness of the proposed
technique is illustrated by successful experiments with
real elevation data sets in a later section.
2 Formulation
I ,,.i 1: The smoothing effect of the subdivision pro
cess on the triangles of the initial mesh.
The butterfly subdivision scheme [9], like 11ii i
other subdivision schemes used in computer graph
ics applications, starts with an initial triangular mesh
which is also known as control mesh. The vertices
of the control mesh are known as control points. In
each step of subdivision, the initial (control) mesh is
refined through the transformation of each triangular
face into a patch with four triangular faces as illus
trated in the Fig.1. After one step of refinement, the
new mesh in the finer level retains the vertices of each
triangular face in the previous level and hence, in
terpolates the coarser mesh in the previous level. In
addition, every edge in each triangular face is spilt by
adding a new vertex whose position is obtained by an
;.! ll combination of the neighboring vertex positions
in the coarser level as shown in I i, 2' i The name,
butterfly subdivision, originated from the "I ii 1 111,l '
like configuration of the contributing vertices. The
weighting factors for different contributing vertex po
sitions are shown in Fig.2(b). The vertex ej+1 in the
Vi 8. 
4 12 2
(a) (b)
Si I, 2: (a) The contributing vertices in the jth level
for the vertex in the j+lth level corresponding to the
edge between vI and v2; (b) the weighing factors for
different vertices.
j + 1th level of subdivision, corresponding to the edge
connecting vertices vi and v. at level j, is obtained by
j+1 = jvv v
12 = 0.5(vi+v2)+2. 1 +vw +v +v V ),
where 0 < w < 1, and vJ denotes the position of the
ith vertex at the jth level.
2.1 "Nat rural" Surface Generation
The original butterfly subdivision scheme produces
a smooth surface in the limit. However, we are I i i 
to recover natural terrains, which are rough. There
fore, we need to modify the butterfly subdivision
scheme such that the limit surface looks like "i Ii 1I I!"
terrain. We adopt a similar technique as in the succes
sive random addition method [2] for generating frac
tal surfaces. In this successive random addition tech
nique, an equally spaced rectangular grid is refined by
interpolating the midpoints of each rectangular cell
(thereby dividing each rectangular cell into four rect
angular cells) and then all grid positions are perturbed
by addition of a Gaussian noise. This process is car
ried out recursively to obtain a fractal surface whose
roughness is controlled by the variance of the added
Gaussian noise at [1!.1 i. !! refinement levels. We have
. ..., I'fi, ,I the ',.!,i, fl/., subdivision scheme by perturbing
the vertex positions at various levels of subdivision by
addition of a Gaussian noise whose variance controls
the roughness of the resulting limit r ..,. f. This pro
cess of adding Gaussian noise is similar to that of the
successive random addition method mentioned above,
the difference being the vertex positions obtained us
ing butterfly subdivision rules are perturbed instead of
grid points obtained through midpoint interpolation.
2.2 Local Parameterization
The limit surface defined by the recursive subdivi
sion process is of arbitrary I1. 1. I. _ where a global
parameterization i!! not be possible. However, we
can locally parameterize the limit surface over the
domain defined by the initial mesh. The idea is to
track i,!! arbitrary point on the initial mesh across
the meshes obtained via the subdivision process with
Gaussian noise addition, so that a correspondence can
be established between the point being tracked in the
initial mesh and its mapping on the limit surface. We
note that the limit surface can be represented by same
number of triangular patches as that of the triangular
faces in the initial mesh (see Fig.l). Therefore, the
limit surface s can be expressed as
s = Sk, (1)
k=l
where n is the number of triangular faces in the initial
mesh and sk is the triangular patch in the limit surface
corresponding to the kth triangular face in the initial
mesh.
I !,.1. 3: Tracking a point x through various levels of
subdivision.
We are now ready to describe the parameterization
of the limit surface over the initial mesh. The process
is best explained through an example. We choose a
simple planar mesh shown in I i, ; as the initial mesh.
An arbitrary point x inside the triangular face abc is
tracked over the meshes obtained through subdivision.
It ii be noted that .,./. point inside the triangular
patch in the limit '.. f .. corresponding to the face abc
in the initial mesh depends (/Il, on the vertices in the
initial mesh which are within the 2neighborhood of
the vertices a, b and c due to the local nature of the
subdivision process.
In the rest of the formulation, superscripts are used
to indicate the subdivision level. For example, vwi
denotes the collection of vertices at level j which con
trols the patch in the limit surface corresponding to
the triangular face uvw at the jth level of subdivi
sion. Let vb be the collection of vertices in the
initial mesh which are within the 2neighborhood of
the vertices a,b and c (marked black in Fig.3). Let
the number of such vertices be r. Then, the vector
vbO, which is the concatenation of the (x, y, z) posi
tions for all the r vertices, is of dimension 3r. These r
vertices control the triangular patch in the limit sur
face corresponding to the triangular face abc in the
initial mesh. Now, there exists (3r x 3r) subdivi
sion matrices (Aabc)t, (Aabc)1, (Aabc)r and (Aabc)m
such that v1bf (1v0 le 0
such that v = (Aabc)Vbc V = (Aabc)vbc
vfe = (Aabc)rVbc, and vdef (Abc)Vabc, where
the subscripts t, 1, r and m stand for top, left, right
and middle respectively (indicating the relative posi
tion of the new triangle with respect to the original tri
angle), and V v fe and vcef are the concate
nation of the (x, y, z) positions for the vertices in the
2neighborhood of the corresponding triangle in the
newly obtained subdivided mesh. The new vertices
in this level of subdivision are marked green in Fig.3.
The 2neighborhood configuration of the vertices in
the newly obtained triangles is exactly the same as
that of the original triangle, hence local subdivision
matrices are square and the vector dimensions are the
same. C il out one more level of subdivision,
along with the old vertices, we get a new set of ver
tices which are marked in magenta in I i, ; Adopting
a similar approach, we obtain v g = (Abed)tV ed
(Aabc)tVbed, Vhg = (Abed) ed = (Aabc)led, V ih
(Abed),Vd = (Aabc)V, ed, Vhi = (Abed)Vd ed
(Aabr),lVed.
Let x be a point with '. i !Iri, coordinates
(aobc, fObc>, n ) inside the triangular face abc. When
the initial mesh is subdivided, x becomes a point in
side the triangular face bed with '. ! i.1 i i coordi
nates (a( ed' l ed .) Another level of subdivision
takes x inside the triangular face I 1,' with '. i. I i i
coordinates (a dg, 0g2i, gi). Let sb denote the jth
level approximation of the triangular patch Sabc in the
limit surface corresponding to the triangular face abc
in the initial mesh. Now vO can be written as
r r r
abc bx, cx,...,ay,by, cy,...,a ,b,,c ...]
where the subscripts x, y and z indicate the x, y and
z coordinates respectively of the corresponding vertex
position. The expressions for v ed and v2 can also
be written in a similar manner. Next, we construct
the matrix Bacb as follows:
r r
abe, bc f bc0 0, ...., 0, O,...,0, O,..., 0
0, ... a0, fbc Obc, Ybc, 0, . ., 0, ..., O
r r r
O,.. ., ,... bc, 0 ...
oabe
The matrices Bcd and Bgi can be constructed in a
similar fashion. We can now write s bc(), Sbc(x),
and s bc(x) as
sobc(x)
s c(x)
Sabc(x)
b, b(X)Vbc,
B ,Vik(Aabc)1VA bc,
B ,i(Aabc)t(Aabc)1V .bc.
Proceeding in a similar fashion, we can derive an
expression for s b,(x), the jth level approximation
of Sabc(x). It can be shown that the expression for
s bc(x) can be simplified and written as
sabc(x) = B ,xi(Abc)vbc
= Bbc(X)Vbc, (2)
where x is inside the triangular face uvw at level j and
B'bc,(x) = B xi(A'bc). I ii! I11,, we can complete
the local parameterization process by writing
Sbc(x) = (lim Bbc (x))b = BbXVbc (3)
j= Bb a cbc. a
In the above equation, Babc is the collection of ba
sis functions at the vertices of v bc. It II also be
noted that the butterfly subdivision scheme is a sta
tionary subdivision process, and hence new vertex po
sitions are obtained by ;.ihl. combinations of nearby
vertices. This guarantees that each row of the matri
ces (Aabc)t, (Aabc)l, (Aabc)r and (Aabc)m sums to one.
The largest eigenvalue of such matrices is 1 and there
fore the limit in Eqn.3 exists. Now, if we assume that
the triangular face abc is the kth face in the initial
mesh, then Eqn.3 can be rewritten as
sk (x) = Bk (x, %" = Bk(x)Akp, (4)
where p is the concatenation of the (x,y,z) positions
of all the vertices in the initial mesh and the matrix
Ak, when postmultiplied by p, selects the vertices
v controlling the kth triangular patch in the limit
surface. If there are t vertices in the initial mesh and
r of them control the kth patch, then p is a vector of
dimension 3t, Ak is a (3r x 3t) matrix and Bk(x) is a
(3 x 3r) matrix.
Combining Eqn.1 and Eqn.4, we get
s(x) = ( Bk(x)Ak)p = J(x)p, (5)
k=l
where the (3 x 3t) matrix J is the collection of basis
functions for the corresponding vertices in the initial
mesh. The vector p is also known as the degrees of
freedom vector of the limit surface s.
We now treat the vertex positions in the initial
mesh defining the limit surface s as a function of time
in order to develop the new dynamic butterfly subdi
vision model. The velocity of the surface model can
be expressed as s(x, p) = J(x)p, where an overstruck
dot denotes a time derivative and x e So, So being
the domain defined by the initial mesh.
3 Finite Element Implementation
We have already pointed out in Section 2 that the
limit surface obtained by the recursive application
of the butterfly subdivision rules along with Gaus
sian noise addition can be represented by triangular
patches. We consider each patch in the limit surface
as an element. The number of such patches is equal
to the number of triangular faces in the initial mesh
as mentioned earlier. The governing motion equation
of the dynamic model is given by
Mpi + Dp + Kp = fp,
where fp is the generalized force vector and M, D, and
K are the mass, damping and stiffness matrices of the
model. We use the same example as in Section 2 (re
fer Fig.3) to develop the mass, damping and stiffness
matrices for the element corresponding to the trian
gular face abc in Fig.3. The derivations hold for ;io
element in general.
The mass matrix for the element given by Sab,
(Eqn.3) can be written as
N = JSb t(x)B b((x)Bab,(x)dx.
JxESabc
However, in this paper, we develop a simplified
method for the elemental mass matrix derivation. The
triangular patch in the limit surface corresponding to
the face abc in the initial mesh is approximated by
a triangular mesh with 4J faces obtained after j lev
els of subdivision of the original triangular face abc
(each subdivision step splits one triangular face into
4 triangular faces). In addition, we choose a discrete
mass .1. ii function which has nonzero values only
at the vertex positions of the jth subdivision level
mesh. Then the mass matrix can be expressed as
k
p(vi){BI bc(v)} {BI b ( )},
i=1
where k is the number of vertices in the triangular
mesh with 4J faces. This approximation has been
found to be very effective and !. I for implemen
tation purposes. The elemental damping matrix can
be obtained in an exactly similar fashion.
We assign the internal. i. each element in the
limit surface, thereby defining the internal !. i of
the subdivision surface model. We take a similar ap
proach as in the derivation of the elemental mass and
damping matrix and assign the internal !. i to a j
th level approximation of the element. In this paper,
we assign spring !!. to the discretized model as its
internal !!. i For the example used throughout the
paper, this !!. i at the jth level of approximation
can be written as
, abc
0 v vQ
1= {jb}T (Kjbc) vjb,
where ki is the spring constant, vi and vJ, the 1th
and mth vertex in the jth level mesh, are in the 1
neighborhood of each other, Q is the domain defined
by all such vertex pairs, is the natural length of
the spring connected between vi and vi, and v3
n' abc
is the concatenation of the (x,y,z) positions of all the
vertices in the jth subdivision level of the triangular
face abc in the initial mesh. It iii be noted that
the entries in (KJb,) depend on the distance between
the connected vertices and hence, (KEb,), unlike other
elemental matrices, is a function of time which needs
to be recomputed in each time step.
The expression for the internal ii. can be
rewritten as
1 v( IT}(Kj.{vab}
aEbc 2 v abc (Kjabc) Vabc I
1 f (A'bc) Vab } T (K bc) (A' bc) Vab,
SVabT{(Ab)T(Kj)(Ab) }IVab, ,
where (Ajb,) and v
where (Abc) and vabc are same as in Eqn.2. There
fore, the expression for the elemental stiffness matrix
is given by
Kab, = (Ab )T (Kjb)(Ab).
The generalized force vector fp in Eqn.6 represents
the net effect of all externally applied forces. We at
tach springs from each point in the elevation data set
to the nearest point on the model and the exerted
point forces deform the initialized model to fit the
given data set.
3.1 Model Subdivision.
The initialized model grows dynamically accord
ing to the equation of motion (Eqn.6) and when an
approximate fit to the given elevation data set is
achieved, the number of control vertices can be in
creased by replacing the original initial mesh by a
new initial mesh obtained by 1i''l i!i a single step of
the '..,i. ., fl., subdivision algorithm. This increases the
number of degrees of freedom to represent the same
limit surface and a new equilibrium position for the
model with a better fit to the given data set can be
achieved. The error of fit criteria for the discrete data
is based on distance between the data points and the
points on the limit surface where the corresponding
springs are attached.
4 Results
We present two natural terrain synthesis results in
this section. The initialized model is deformed by ap
I'1 ii, spring forces on its limit surface from the dis
crete data points. At each time step, every control
vertex position is perturbed by adding a random noise
drawn from a Gaussian distribution. The variance of
the Gaussian distribution determines the roughness of
the !ill i .1 surface. In both the experiments, the
initialized butterfly subdivision model has an initial
(control) mesh with 98 triangular faces and 68 con
trol vertices. The 'natural" looking limit surface of
the initialized model is deformed by the forces exerted
from the discrete elevation data points. When an ap
proximate fitting is obtained, the model is subdivided
to obtain a closer fit using more degrees of freedom
(control vertices) of the new initial mesh. The fit
ted surface has 1568 triangular faces and 841 control
vertices in both the experiments. It !1i ~ be noted
that !1. i of same ,i1,1i i natural terrains using
the existing techniques requires a large number of grid
points (of the order of 10') [4, 5, 6, 7] and hence the
proposed technique provides a more compact represen
tation of the 0ill. i .1 terrain. The elevation data
values are scaled to fit an unit cube and the variance of
added noise is 104 for the  i 1.1 i , 1 fractal surface
depicted in Fig.4. The corresponding value of noise
variance for fractal surfaces depicted in I i, ., ,, (b)
and (c) are 106, 105 and 104 respectively. In the
first experiment, 4096 elevation data points are used
whereas the second data set comprised of 3099 eleva
tion values. The error in fitting is approximately one
percent in both the experiments.
5 Conclusions
In this paper, we have presented a novel technique
for synthesizing natural terrains from elevation data
by using a finite element methodbased dynamic but
terfly subdivision surface model in . .i !!i, I i i with
a modified successive random addition technique. We
have modified the butterfly subdivision scheme and
presented a local parameterization of the subdivision
scheme. We have also incorporated the advantages of
freeform deformable models in the modified butter
fly subdivision scheme and have developed a new I, 
of finite element based on the butterfly subdivision
scheme. The In  of our model is shown via exper
iments indicating a promising future of the proposed
model in natural terrain modeling.
References
[1] A. Fournier, D. Fussel, and L. Carpenter, "Com
puter rendering of stochastic models," Communi
cations of the ACM, vol. 25, no. 6, pp. 371384,
1982.
[2] D. Saupe, "Algorithms for random fractals," in
I !. Science of Fractal Images, H.O. Peitgen and
D. Saupe, Eds., pp. 71 136. Springer Verlag,
1988.
[3] N. Yokoya, K. Yamamoto, and N. Funakubo,
"Ii I Ibased ;,!! 1 i and interpolation of 3d
natural surfaces and their application to terrain
modeling," Computer Vision, Graphics and Im
age j., ....... vol. 46, pp. 284302, 1989.
[4] R. Szeliski and D. Terzopoulos, "1 ..!!! splines to
fractals," in Proceedings of AC f. SIGGRAPH'89,
August 1989, pp. 5160.
[5] B.C. Vemuri, C. Mandal, and S.H. Lai, "A fast
Gibbs sampler for synthesizing constrained frac
tals," IEEE I,......; ,. on Visualization and
Computer Graphics, vol. 3, no. 4, pp. 337 351,
October December 1997.
(a) (b) (c)
I o', 4: (a) discrete elevation data set (4096 points),
(b) fitted dynamic butterfly subdivision surface model
with 841 vertices (without noise addition), and (c) fit
ted dynamic subdivision surface model with 841 ver
tices when Gaussian noise is added.
(a) (b) (c)
1 i,oi' 5: Synthesized natural terrain of different
roughness using the dynamic butterfly subdivision
surface model with 841 vertices from a data set of
3099 elevation values.
[6] K. Arakawa and E. Krotkov, "i ,. I i! surface re
construction for modeling natural terrain," in Pro
ceedings of the IEEE C,.f' ..... on Computer Vi
sion and Pattern Recognition, 1993.
[7] W. Riimelin, "I i, I interpolation of random
fields of fractional brownian motion," in Frac
tal Geometry and Computer Graphics, J.L. En
carnacao, H.O. Peitgen, G. Sakas and G. Englert,
Eds., pp. 122 132. SpringerVerlag, 1992.
[8] C. Mandal, B.C. Vemuri, and H. Qin, "S .p. re
covery using dynamic subdivision surfaces," in
Proceedings of the International C.. f,. .... on
Computer Vision, DB..iil, India, January 1998,
pp. 805 810.
[9] N. D D. Levin, and J.A. Gregory, "A butterfly
subdivision scheme for surface interpolation with
tension control," .I' .1 l. .I..... .. on Graphics,
vol. 9, no. 2, pp. 160 169, April 1990.
