The Visual Computer manuscript No.
(will be inserted by the editor)
Huaiping Yang · Bert J¨uttler
Evolution of T-spline Level Sets for Meshing Non-uniformly
Sampled and Incomplete Data
Abstract Given a large s e t of unorg anized point sample
data, we propose a new framework for co mputing a tri-
angular mesh representing an approximating piecewise
smooth surface. The data may be non–uniformly dis-
tributed, noisy, and they may contain holes. This frame-
work is based on the combination of two types of surface
representations: triangular meshes, and T-spline level sets,
which are implicit surfaces defined by refinable spline
functions allowing T-junctions. Our metho d contains thre e
main steps. Firstly, we construct an implicit representa-
tion of a smooth (C
2
in our case) surface, by using an
evolution process of T-spline level s e ts, such that the
implicit surface captures the topology and outline of the
object to be reconstructed. The initial mesh with high
quality is obta ined through the marching triangulation
of the implicit surface. Secondly, we project each data
point to the initial mesh, and get a scalar displacement
field. Detailed features will be captured by the displaced
mesh. Finally, we present an additional evolution pro-
cess, which combines data-driven velocities and feature-
preserving bilateral filters, in order to reproduce sharp
features. We also show that various shape constraints,
such as distance field constraints, range constraints and
volume constraints can be naturally added to our frame-
work, which is helpful to obtain a desired reconstruction
result, especially when the given data contains noise and
inaccuracies.
Keywords Mesh reconstruction · Point cloud ·
Displacement maps · T-spline · Level sets
1 Introduction
We consider the problem of surface reconstruction and
approximation from unorg anized data points. Due to the
H. Yang, B. J¨uttler
Institute of A pplied Geometry, Johannes K epler University
4040 Linz, Austria
E-mail: yang.huaiping@jku.at, bert.juettler@jku.at
widespread use of 3D scanning devices for shap e acquisi-
tion, this problem has an increasing number of applica-
tions in computer graphics, computer aided design, com-
puter vision and image processing. Depending on the
area of the application, the reconstructed s urface may
have a explicit representation (e.g. meshes) o r an im-
plicit representation (e.g. level sets). Among the various
approaches, these two represe ntations may complement
each other. On the one hand, the implicit representa-
tions [51] offer advantages such as the non-existence of
the parametrization problem, repairing capabilities of in-
complete data and simple operations of shape editing,
but they are hard to model sharp features [39]. On the
other hand, the explicit representations can easily ha n-
dle sharp features, but have difficulties when processing
topology changes.
1.1 Our Work
We propose a hybrid model for sur face reconstruction by
combining two types of r e presentations: an implicit T-
spline level set and a mesh. Given a set of unorganized
and noisy data points without normals as input, we want
to reconstruct a mesh surface which approximates the
data. We develop a three-phas e algorithm (cf. Fig. 1) to
perform this reconstruction:
1. Initial mesh generation (Fig. 1 (c)). In the first phase,
we use an evolution pro c e ss (Fig. 1 (a) and (b)) to
create an implicit representation, which is defined as
the zero level set of a C
2
T-spline scalar function. The
obtained T-spline level set (with correct topology) is
to serve as a smooth base sur fa c e S
0
for the displace-
ment mapping. A high-quality initial mes h (with ac-
curate normals) is generated from the implicit func-
tion by using the marching triangulation [26] method.
2. Displacement mapping (Fig. 1 (d)). In the sec ond
phase, we produce a smooth scalar displac e ment field,
which is computed by projecting data points to the
initial mesh. Small geometric features are then con-
structed by the displace ment mapping of the mesh
2 Huaiping Yang, Bert J¨uttler
(a) (b) (c) (d) (e)
Fig. 1 Mesh reconstruction of the Rocker-Arm model. The figure shows the data points and the T-mesh in (a), an interme-
diate T-spline level set during the evolution in (b), the final T-spline level set (the initial mesh) in (c), the displaced mesh
in (d) , and the final mesh with sharp features in (e).
along the normal direction, which is guided by the
smooth gradient vector field of the implicit function.
3. Recovering sharp features (Fig. 1 (e)). In the third
phase, we use an additional evolution process, which
combines data-driven velo c ities and feature-pre serving
bilateral filters, in order to better represent sharp fea-
tures.
This paper is an extensio n to our work [55] presented
at Shape Modeling International 2007. In this extended
version, we concentrate more on the pro cess generating
the base surfa c e , i.e., we show how to formulate the evo-
lution of T-spline level sets. More specifically, we dis-
cuss how to combine different shape constraints, such
as distance field constraints, ra nge constr aints and vol-
ume constraints, into the framework such that a more
robust and effective evolution process ca n be established
that produces the desired result. The distance field con-
straints help us to avoid additional branches and singu-
larities o f the level sets, without having to use costly
re-initialization steps. The range constraints allow us
to specify regions lying inside or outside of the recon-
structed s urface. A suitable volume constraint can be
used to stop the level sets from entering the holes in the
data. We note that these cons traints have been used in
our previous work [18] in 2D for dual evolution of planar
B-spline curves and T-spline level sets.
Our method combines two types of representations:
the implicit T-spline level set and the mesh. This com-
bination strategy makes our method benefit from the
advantages of both representations. On the one hand,
the evolution of T-spline level sets is able to capture the
complex topology of noisy data. On the other hand, the
mesh re presentation helps to produce detailed a nd sharp
features. Compared with other existing approaches for
similar purposes, our method a lso shows the following
advantages:
1) The use of two complementary representations can
improve and speed up some geometr ic computations in
the algorithm. For example, with the help of implicit rep-
resentation, a high-q uality initial mesh can be obtained,
and the projection of the data points to the base surface
can be efficiently computed by Newton iteration.
2) We us e an evolution process to recover the sharp
features. The evolution is governed by a combination of
two terms: a data-driven velocity and a bilateral filtering.
This evolution process can produce sharp features, which
are faithful to the given data.
3) By incorporating the r ange constraints and volume
constraints, we are able to conveniently exploit the a pri-
ori knowledge about the geometric properties of the ob-
ject to be r e constructed. This is especially helpful when
the given data contains noise and inaccuracies.
1.2 Related Work
Numerous approaches have been proposed to compute
a surface approximation of a given set of unorganized
points. Depending on the area of the application, differ-
ent representations have been used, such as tria ngular
meshes [3,7,15,38,2,41,47], subdivision surfaces [28,4 4,
12], parametric spline surfaces [16,23], discretized level
sets [42,57,8], sca lar s pline functions [45,33], radial basis
functions [10,40] and point set surfaces [1,43,19]. It is be-
yond the scope of this paper to give a detailed overview
of all the existing work. [29] gives an excellent survey
of the previous work on surface ex traction from point
clouds.
In this paper, we suggest a new reconstruction algo-
rithm combining two types of repres e ntations: an implicit
T-spline level set and a mesh. T he proposed method r e -
lies on two main tools: displacement mapping [14] and
bilateral filtering [48,50]. In the re mainder of this sec-
tion, we will describe some related work about these two
tools.
Evolution of T-spline Level Sets for Meshing Non-uniformly Sampled and Incomplete Data 3
Given a smooth base surface S
0
, a displaced surface S
can be generated by a scalar field (a displacement map),
which specifies the displacement values along the norma l
directions of S
0
. The use of displacement maps is quite
popular for geometric modeling purpos e s. They are used
in high end rendering systems [14,5], to capture the fine
detail of a 3D photography model [36], for geometric sim-
plification with appearance-preserving [13], for building
semi-regular multiresolution meshes from an arbitrary
connectivity input mesh[25 ] and fo r multiresolution mesh
deformations [35 ]. In order to avoid c racks between ad-
jacent triangles of a mesh, the interpolated norma l is
used [2 4] to displace the surface, B- spline surfaces are
fitted [36] to the mesh before the displacement mapping,
displaced subdivis ion surfaces [37] are suggested which is
based on the butterfly subdivision scheme. There a re also
existing approaches for reconstructing a displaced sub-
division sur fa ce directly from a given set of points [31].
Most recently, the author in [56] presents a displace d
surface representation based on a manifold structure.
Extracting sharp featur e s from 3D data is impor-
tant [49], but difficult due to the feature-insensitive sam-
pling and the noise of the given data. Many approaches
have been proposed to addres s this problem [27,30,53,
6]. As a non-itera tive scheme for edge-preserving smooth-
ing, the bilateral filter is used in [20] and [32] to denoise
a given surface. The author in [52] presents a robust gen-
eral approach conducting bilateral filters to govern the
sharping of triangular meshes. Also , the authors in [4]
conduct the bilateral filter in the reconstruction o f sur-
faces from scattered data.
The remainder of the paper is organized as follows.
The next sec tion describes how to formulate the evo-
lution process of T-spline level sets, in o rder to gener-
ate the smooth base surface for the object reconstruc-
tion. Section 3 discuss e s the co mbination of different
constraints into the evolution process. In particular, it
is shown how to incorporate distance field constraints,
range constraints and volume cons traints such that the
T-spline level sets evolution will be more robust and ef-
fective when dealing with non-uniformly sampled and
incomplete data. Section 4 pre sents how to reconstruct
geometric details and sharp features of the object by de-
riving a mesh repre sentation from the smo oth implicit
surface. After presenting some experimental results in
Section 5, we conclude this paper and discuss future
work.
2 Evolution of T-spline Level Set for Base
Surface Generation
In this section, we describe how to generate the base
surface through the evolution of T-spline level sets. The
input of our algorithm is a set of unorganized (maybe
noisy and defected) data points (p
k
)
k=1,2,...,n
, which are
scattered over an unknown piecewise smooth surface S
p
.
The base surface represented S
0
by T-spline level sets
will capture the topology and the outline of the surface
S
p
.
2.1 Definition of T-spline Level Sets
The T-spline level set Γ (f) is defined as the zero set of
a trivariate T-spline function f over some domain D,
Γ (f) = { (x, y, z) ∈ D ⊂ R
3
| f (x, y, z) = 0 }, (1)
where
f(x, y, z) =
P
n
i=1
c
i
B
i
(x, y, z)
P
n
i=1
B
i
(x, y, z)
(2)
are T-spline [46] functions in 3D, with the real c oeffi-
cients (control points) c
i
, i = 1, 2, ..., n. For cubic T-
splines, the basis functions are
B
i
(x, y, z) = N
3
i0
(x)N
3
i0
(y)N
3
i0
(z), (3)
where N
3
i0
(x), N
3
i0
(y) and N
3
i0
(z) are certain cubic B-
splines, whose associated knot vectors are determined
by the T-spline control grids (T-mesh).
In order to simplify the notation, we use x to repre-
sent the point x = (x, y, z) in 3D, and gather the control
coefficients (in a suitable ordering) in a column vector c.
The T-spline basis functions form another column vector
b = [b
1
, b
2
, ..., b
n
]
⊤
, where
b
i
=
B
i
(x)
P
n
i=1
B
i
(x)
, i = 1, 2, ..., n.
Now the T-spline function f(x) = b(x)
⊤
c.
T-splines [46] are generalizations of tensor product B-
splines. T-spline control grids permit T-junctions, which
provides a valuable property of local refinement. As shown
in Figure 2, much less (decreased by 71%) control points
are needed to represent the same level set (a horse) by
trivaria te T-splines than that by trivariate tensor prod-
uct B-splines. Since a T-spline function is piecewise ra-
tional, the T-spline level sets are piecewise algebraic sur -
faces (C
2
in our case).
The T-spline control grid (T-mesh) (cf. Fig. 1 (a)) is
constructed through octree subdivision of the function
domain D , as follows.
1. Set the initial T-mesh to be an axis-aligned bounding
box of the cubic domain which contains all the data
points.
2. For e ach cell co ntaining more than n
0
data points
(n
0
is a user-specified constant value), subdivide it
by applying the octree subdivision.
3. Repeat step 2 until a user-s pecified thr eshold (e.g., a
maximum level of subdivision) is r eached or no cell
contains more than n
0
data points any more.
In this way, the distribution of T-spline control points is
made adaptive to the geometry of the object to be recon-
structed, which usually leads to a sparse repre sentation.
4 Huaiping Yang, Bert J¨uttler
(a) B-splines with 5148 control p oints (b) T-splines with 1484 control points
Fig. 2 A horse implicitly defined by trivariate B-spline (T-spline) scalar functions.
2.2 Evolution of T-spline level s ets
Consider a T-spline level set Γ (f) defined as the zero set
of a time-dependent function f(x, τ), where
f(x, τ) = b(x)
⊤
c(τ), (4)
with some time parameter τ. It will be subject to the
evolution process
∂x
∂τ
· n(x, τ) = v(x, τ) · n(x, τ), x ∈ Γ (f), (5)
which is driven by a (possibly time-dependent) vector
field v. Note that we only consider the normal compo-
nent (v · n) o f the velocity along the level set, since the
tangent veloc ity does not change the shape at all. The
unit normal vector n is given by
n(x, τ ) =
∇f(x, τ)
|∇f(x, τ)|
. (6)
During the evolution, the definition of T-spline level
sets
f(x, τ) ≡ 0, x ∈ Γ (f), (7)
implies
∂f(x, τ)
∂τ
+ ∇f(x, τ) ·
∂x
∂τ
= 0, x ∈ Γ (f ). (8)
Combining (5), (6) and (8), we get the evolution equation
of T-spline level sets under the vector field v,
∂f(x, τ)
∂τ
= −v(x, τ) · ∇f(x, τ), x ∈ Γ (f). (9)
In our method, we always start the evolution of T-
spline level sets from an initial level set, which contains
all data points inside. The initial level set function f can
be found by approximating it to the signed distance field
of a bounding sphere or a rough offset of the surface to
be reconstructed [57].
2.3 Evolution Speed Function
The speed function v plays a a key role in the evolution
of T-spline level s ets. In order to capture the base surface
of the given data points, we use a similar speed function
as that proposed by Caselles et al. [11],
v = e(d)(γ + κ)n − (1 − e(d))(∇d · n)n (10)
where γ is a constant velocity (which is also known as a
balloon force), κ = div(∇f/|∇f |) is the mean curvature
of the level set surface, d is the unsigned distance func-
tion of the data points, and e(d) = 1 − e
−ηd
2
is the edge
detector function. η is pre-defined, and its value is af-
fected by the size of the data range. In our exper imental
setting, all data points are contained in a unit bounding
box (−1 ≤ x, y, z ≤ 1), and we usually set η = 1.
The spee d function (10) is a linear combination of
two par ts. In the first part, the mean curvatur e term
κ makes the level set smooth, and the ba lloon force γ is
used to increase the speed of the evolution (especially for
capturing narrow concave boundaries). The second part
attracts the level set to the detected boundary edges.
The edge detector function e(d) is used as weighting co-
efficients to balance the influence of the two parts.
2.4 Discretization of the Evolution Equation
The evolution equation (9) under the speed function (10)
is discretized by uniformly sampling a set of points x
j
,
j = 1, . . . , N
0
on the T-spline level set. We use the
marching triangulation method to generate the uniform
samples, which will be described later in Section 4.1.
Since usually the number of sample p oints N
0
is much
larger than the number of T-s pline control coefficients
n, the equation ca n not be exa c tly sa tisfied. We use a
least-squar e s approach to choose the time deriva tive of
Evolution of T-spline Level Sets for Meshing Non-uniformly Sampled and Incomplete Data 5
the function f by solving
E =
N
0
X
j=1
(
˙
f(x
j
, τ) + v(x
j
, τ) · ∇f (x
j
, τ) )
2
→ Min (11)
where
˙
f(x
j
, τ) = ∂f (x
j
, τ)/∂τ is the time derivative of f.
In order to prevent the problem from being ill-pos e d,
we a dd a regulariz ation term to the above function E to
be minimized,
E + λ
c
|
˙
c|
2
→ Min, (12)
where
˙
c is the time derivative of the T-spline control
coefficients, and λ
c
is a pre-scribed small constant [17].
Since in our ca se
˙
f(x, τ) = b(x)
⊤
˙
c(τ) (cf. Equa-
tion (4)), the solution to (12) can be found by solving
a sparse linear system of equations.
Then we generate the updated control coe fficients
c(τ + ∆τ) = c(τ) +
˙
c∆τ. (13)
by using an explicit Euler step ∆τ. ∆τ is chosen as
∆τ = min(1, {
h
|v(x
j
, τ) · n(x
j
, τ)|
}
0≤j≤N
0
) (14)
where h is a user-defined constant to indicate the ma x-
imum allowed evolution step size for each sample point
on the T-spline level set.
3 Constraints for T-spline level sets evolution
In this section, we present thre e different constra ints in-
cluding distance field constra ints, range constraints, and
volume constraints, which can be applied to improve the
robustness and effectiveness of T-spline level sets evolu-
tion.
3.1 Distance Field Constraint
The distance field constraint is used to restore the signed
distance property of the T-spline level set function during
the evolution, without having to use the c ostly level set
reinitialization proc e dure [54].
Since an ideal signed distance function φ satisfies
|∇φ| = 1 everywhere in the domain of D, we pro pose
the following constraint ter m
S
0
=
Z
D
(
∂|∇f(x, τ )|
∂τ
+|∇f(x, τ)|−1 )
2
dV → Min, (15)
where dV is the volume element. This constra int acts as
a penalty function to p e nalize the deviation of f from a
signed distance function.
In practice, we realize tha t it is usua lly sufficient to
maintain the signed distance property of f only in a small
neighborhood o f the zero level set (like the narrow-band
(a) (b) (c)
Fig. 3 The distance field constraint (DFC) for T-spline level
sets in 2D. The figure shows the data points with T-mesh in
(a), the resulting level set function without DFC in (b), and
the resulting level set function with DFC in (c).
level set method). Thus we propose to use a much more
efficient version of the distance field constraint,
S =
N
0
X
j=1
(
∂|∇f(x
j
, τ)|
∂τ
+ |∇f(x
j
, τ)| − 1)
2
, (16)
by only considering the sample points x
j
, j = 1, . . . , N
0
on the zero level set. Combining the distance field con-
straint into the evolution equation, we ge t the new ob-
jective function to b e minimized
E + λ
c
|
˙
c|
2
+ λ
s
S → Min. (17)
We usually cho ose the weight λ
s
= 0.1 in our experi-
ments. More details about the influence and choice of λ
s
are describ e d in [54].
Fig. 3 illustrates an example of using the distance
field constraint. As shown in (b), without the using of
distance field constraints, the level set function bec omes
very steep (or flat) in s ome regions, which may cause
difficulties in maintaining the numerical stability during
the evolution. This situation can be cured by applying
the distance field constraint, as shown in (c).
3.2 Range Constraint
One typical a dvantage of the implicit representation of a
surface is fast inside/outside distinction. In our ca se, we
assume that the T- spline level set function f(x) > 0 re p-
resents the region outside the surface, and f (x) < 0 in-
side the surface. The r ange constraints are used to spec-
ify regions which should lie inside/outside the surface to
be reconstructed. More specifically, we co nsider a set o f
points {x
k
}
k=1,..., N
1
which should no t lie outside the zero
level set, and another set of points {y
k
}
k=1,...,N
2
which
should not lie inside the zero level set,
f(x
k
) ≤ 0, k = 1, . . . , N
1
.
f(y
k
) ≥ 0, k = 1, . . . , N
2
.
(18)
The above range constraints can be achieved by adding
corresponding penalty ter ms to the objective function (17)
such that the time derivative satisfies
˙
f(x
k
) < 0 if f(x
k
) >
