Piecewise-linear interpolation between polygonal slices.
In this paper we present a new technique for piecewise-linear surface reconstruction from a series of parallel polygonal cross-sections. This is an important problem in medical imaging, surface reconstruction from topographic data, and other applications. We reduce the problem, as in most previous works, to a series of problems of piecewise-linear interpolation between each pair of successive slices. Our algorithm uses a partial curve matching technique for matching parts of the contours, an optimal triangulation of 3-D polygons for resolving the unmatched parts, and a minimum spanning tree heuristic for interpolating between non simply connected regions. Unlike previous attempts at solving this problem, our algorithm seems to handle successfully any kind of data. It allows multiple contours in each slice, with any hierarchy of contour nesting, and avoids the introduction of counter-intuitive bridges between contours, proposed in some earlier papers to handle interpolation between multiply connected regions. Experimental results on various complex examples, involving actual medical imaging data, are presented, and show the good and robust performance of our algorithm.
Introduction
The problem of reconstructing the boundary of a solid
object from a series of parallel planar cross-sections has
attracted much attention in the literature during the
past two decades. The main motivation for this problem
comes from medical imaging applications, where
cross-sections of human organs, such as bones, tumors
and tissues, are obtained by CT (Computerized Tomog-
raphy) or MRI (Magnetic Resonance Imaging) appa-
Work on this paper by both authors has been supported by
a grant from the G.I.F., the German-Israeli Foundation for Scientific
Research and Development. Work on this paper by the
second author has also been supported by National Science Foundation
Grant CCR-91-22103, and by grants from the U.S.-Israeli
Binational Science Foundation, and the Israel Science Fund administered
by the Israeli Academy of Sciences.
y School of Mathematical Sciences, Tel-Aviv University, Tel-Aviv
69978, Israel, and Algotec Systems Ltd., Raanana, Israel.
z School of Mathematical Sciences, Tel-Aviv University, Tel-Aviv
69978, Israel, and Courant Institute of Mathematical Sci-
ences, New York University, New York, NY 10012, USA
rata. These cross-sections, hereafter called slices, are
the basis for interpolating the boundary surface of the
organ. The interpolated object can then be displayed
in graphics applications, or (more recently) even manufactured
by an NC (Numerically Controlled) or an RP
(Rapid Prototyping) machine. Another motivation for
this problem is the non-destructive digitization of ob-
jects: after an object is scanned by an echo-graphic or
an X-ray apparatus, the obtained slices are used for the
reconstruction of the original object. Yet another motivation
is the reconstruction of a 3-dimensional model
of a terrain from topographic elevation contours.
Many solutions were suggested for the pure raster in-
terpolation. These usually handle two raster images,
where each pixel is either white or black, or assigned
a grey-level taken from a fixed range. The interpolation
produces one or more intermediate raster images,
which smoothly and locally turn the first image into the
second one. Then, the bounding surface is detected using
other methods, such as edge detection techniques,
for identifying places of transition from the inside to
the outside of the object. In the grey level case, these
methods include some thresholding mechanism which
decides which levels are 'inside' the object and which
are not. Cline et al. [6] attempted to convert directly
the voxel data into a polyhedral surface, suggesting the
marching cubes technique, which produced very small
triangles whose size was roughly the same as that of the
input voxels.
Many other solutions, including the approach taken
in this paper, assume that the interpolation is preceded
by an edge-detection process, which is invoked for each
of the slices. Thus, each slice is then assumed to be represented
by a hierarchy of non-crossing contours, each
being a closed simple Jordan curve, which represent the
boundaries between "material" and "non-material" ar-
eas; in general, the depth and breadth of this hierarchy
is not restricted, and a contour may enclose any number
of other contours, which themselves may enclose
other contours, and so on. In practice, each contour is
given as a discrete circular sequence of points along it,
and we can thus regard it as a simple closed polygonal
curve, whose vertices are the given points. Finally, we
may also assume that the exterior, unbounded region in
each planar slice represents "non-material" (the model
is assumed to be bounded).
Thus the problem that we face is: Given a series of
parallel planar slices (which we will assume to be parallel
to the xy-plane), each consisting of a collection
of non-crossing, but possibly nested, closed and simple
polygonal curves, with the above properties, we want
to reconstruct a polyhedral solid model whose cross
sections along the given planes coincide with the input
slices. A natural simplification of the problem, also
taken in most earlier works, is to consider only a single
pair of successive parallel slices, and to construct a solid
model within the layer delimited by the planes of the
slices, which interpolates between the given slices. The
union (or, rather, concatenation) of these models will
give us a solution model for the full problem.
Before continuing, we should remark that the solution
is not uniquely defined, and the measure of 'goodness'

Figure

1: Contour association and tiling
of a proposed solution is rather subjective and intu-
itive. Of course, if each of the two slices consists of a
single contour, and these contours roughly 'sit above
each other', then we expect the solution to be a single
'drum-like' polytope whose boundary consists of a circular
sequence of triangles 'wrapping around' the two
contours; see Figure 1(a). However, even in the simple
case of one contour in each slice, if the xy-projections of
the two contours are far away from each other, it is not
clear which is a better solution: to construct a highly
slanted 'pipe' that connects between these contours (as
in

Figure

1(b)), or to regard the lower contour as the
top cover of some pillar-like solid, and the upper contour
as the bottom cover of another pillar-like solid (as
in

Figure

1(c)). The choice of a solution can become
much more arbitrary in more involved cases.
We first briefly review the fairly extensive literature
on this problem. Most of the earlier works only studied
the variant where each slice contains only one con-
tour. In the sequel we denote it as the one-to-one case,
as opposed to the one-to-many and the many-to-many
cases. These studies either sought a global optimization
of some objective function, or settled with a local tiling-
advancing rule, after the tiling starting points at the two
contours were somehow determined, e.g. the closest pair
of vertices between the contours.
In the full version we review these works in more de-
tail, but since they are all inferior to our solution, in
that they cannot handle arbitrary data, we omit this
review here. These works are by Keppel [19], Fuchs et
al. [9], Sloan and Painter [28], Cook et al. [7], Christiansen
and Sederberg [5], Shantz [25], Batnitzky et al.
[2], Sloan and Hrechanyk [26], Ganapathy and Dennehy
[10], and by Wang and Aggarwal [29]. A good survey on
all the works cited so far is given by Sloan and Painter
[27]. They decompose each method into its building
blocks, and compared the various methods accordingly.
They also describe a testbed for evaluating and comparing
these techniques.
Zyda, Jones and Hogan [32] made an attempt to handle
the many-to-many case, but their method could not
handle any branching cases, and produced unsatisfactory
results for partially overlapping contours, as in Figure
4. Other limitations involved specific geometries.
However, they suggested solutions for these limitations,
which might help in certain cases, but which required
some interaction with the user.
Boissonnat [3] presented a totally different approach.
He constructed the Delaunay triangulation for each
slice, projected one triangulation onto the other, and
obtained a collection of tetrahedra, aiming to maximize
the sum of their volumes. This was a considerable step
towards handling the case where each slice has multiple
contours. Boissonnat mentioned three typical examples
where his standard method failed to produce good re-
sults, thus requiring special treatment. The first example
contained two overlapping contours but with considerable
differences in their geometry. The second example
consisted of two similar contours, where one of
them also contained a hole polygon. And the last example
showed a branching problem but without contour
overlaps. Our algorithm handles successfully (and eas-
ily) all these 'bad' examples. Boissonnat suggested a
correction scheme, which either changed the geometry
of one of the slices, or constructed one or two intermediate
slices between the original ones. These idea were
further developed by Boissonnat and Geiger ([4], [11]).
There have been a few other recent works that also
attempted to handle the more general cases; these are
by Kehtarnavaz and De Figueiredo [17], Kehtarnavaz,
Simar and De Figueiredo [18], Ekoule, Peyrin and Odet
[8], Meyers, Skinner and Sloan [22], and by Welzl and
Wolfers [30]. However, all these treatments suffered
from a variety of problems, treated only certain restricted
cases, and worked only in favorable situations
where the two slices closely resemble each other.
Finally, Gitlin, O'Rourke and Subramanian [12] prove
that it is not always possible to find an interpolating
polyhedron between two polygons, that lie in two parallel
planes. That is, any attempted interpolation produces
a self-intersecting surface. This result holds only
when the interpolating triangles are all assumed to connect
between the two polygons vertices, i.e. to have two
vertices taken from one polygon and the third from the
other polygon. We do not assume this in our approach,
and we indeed get a non-intersecting interpolation when
we apply our algorithm on their example.
Two comprehensive reviews of many of the works on
reconstructing 3D objects from cross-sections are given
by Schumaker [23] and by Hagen, M-uller and Nielson
[14]. The first review is not restricted to piecewise-linear
constructions based on polygonal slices, but also refers
to parametric and tensor representations of contours,
and describes volumetric and surface approaches.
We propose a new approach to the interpolation prob-
lem. Our algorithm handles well slices with multiple
contours, and does not rely on any resemblance between
the slices. We accept slices which contain any number
of contours, arbitrarily different in their geometries.
The xy-projections of contours of different slices may
arbitrarily overlap; we do not make any distinction in
the treatment of contours which fully overlap, partially
overlap or do not overlap at all. Many of the previous
works, such as [19], [9], [26] and [10], either prohibit
the creation of triangles in the same slice, or specifically
define steps where this action is allowed. We do
not make a distinction between triangles which connect
the two slices and those which totally lie within a single
slice. Generally speaking, our algorithm is based on
an analysis of the resemblance between the two slices.
Thus, we separately treat similar contour portions that
are matched between the two slices, and then treat the
remaining portions that do not match. We refrain from
creating artificial bridges between contours that might
conflict with the geometry of the other slice, as already
noted in [25] and others. In the only case where we construct
such bridges, they are guaranteed not to conflict
with the geometry of the other slice, i.e. not to intersect
the projection of any other contour. We do not have
to introduce intermediate slices. To recap, our algorithm
appears to overcome all the technical difficulties
that hampered the performance of previous solutions,
it treats data in full generality, and the extensive experimentation
that we have conducted indicates that it
performs very well on complicated large-size data. We

Figure

2: Matching contour portions
regard our algorithm as a significant step in the solution
of this problem. For an illustration of the performance
of our algorithm on real-life examples, see Figures 17
and 18.
Here is a brief overview of our algorithm. We first
match similar contour portions between the two slices
(e.g. the upper arcs pq and p 0 q 0 in Figure 2). Then we
'stitch' (or `tile') each pair of matched contour portions
by a sequence of adjacent triangles forming a 'band' between
the portions. With some care, if we take the union
of the original contours and new bands, and cancel out
duplicate edges, we obtain a collection of closed spatial
polygonal curves, each of which may be composed of
pieces of contours on both slices and of some edges of
the connecting triangles. Moreover, our matching procedure
essentially guarantees that the xy-projections of
these curves are pairwise disjoint, although they may
be nested within each other. If no nesting occurs, we
simply triangulate each of these spatial polygons, using
a simple dynamic programming approach, which
roughly aims to minimize the total area of the triangula-
tion. If nesting occurs, we take one polygon P with all
polygons whose xy-projections are directly
nested within that of P , and apply a minimum spanning
tree procedure that introduces edges connecting
between these polygons and yielding an xy-projection
which is simply connected, so we can then proceed to
triangulate the resulting polygonal curve, as above. See

Figure

3 for an illustration. More details of all these
steps will be given later in the paper.
For the purpose of identifying matching portions of
the contours we use a partial curve matching technique,
which was first suggested by Kalvin et al. [16] and
by Schwartz and Sharir [24]. This technique, which
uses the so-called Geometric Hashing method, originally
solved the following curve matching problem: Given two
curves in the plane, such that one is a (slight deformation
of a) proper subcurve of the other, find the translation
and rotation of the subcurve that yields the best
least-squares fit to the appropriate portion of the longer
curve.
This technique was extended and used in computer
vision for automatic identification of partially obscured
objects in two or three dimensions, an important problem
in robotics applications of computer vision, which
has attracted much attention; see Hong and Wolfson
[15], Wolfson [31], and Kishon, Hastie and Wolfson [20].
A simplified variant of this technique has recently
been used by Barequet and Sharir [1] for the totally
different problem of detecting and repairing gaps in the
boundary of a polyhedron, a problem which often arises
in the creation of polyhedral approximations of CAD
models. In the variant used in [1], as well as the one
used in this paper, no motion of one curve relative to
the other is allowed, so the technique becomes considerably
simpler.
The paper is organized as follows. In Section 2 we give
a more precise definition of the problem and present
an overview of the algorithm (more detailed than the
one given above). The later sections describe in detail
certain phases of the algorithm. Section 3 describes the
matching of contour portions, and Section 4 describes
the actual surface reconstruction. In Section 5 we briefly
analyze the complexity of the algorithm, and present
experimental results. We end in Section 6 with some
concluding remarks.
2 Overview of the Algorithm
The input to the algorithm, as described in the Intro-
duction, is a pair of parallel planar slices parallel to the
xy-plane, each slice consisting of a list of closed and
simple polygonal contours, which do not intersect each
other. The containment hierarchy of the contours may
be omitted; in this case we compute it ourselves, using
a simple line-sweep technique. Contours of even nesting
level (starting at 0, which labels external contours)
are such that their interior, in a sufficiently small neighborhood
of the contour, is the "material", and contours
of odd level are those whose interior, sufficiently near
them, is the "non-material". We orient each contour so
that, when viewing the contour from above, the material
lies to the right of the contour (thus all even-level
contours are oriented in the clockwise direction, when
viewed from above, and all odd-level contours are oriented
in the counterclockwise direction). If necessary,
we re-orient the contours in these consistent directions.
We remark that we need to compute the contour hierarchy
only to obtain the consistent orientation of contours;
the hierarchy itself is not used in the algorithm.
Our proposed algorithm consists of the following
steps:
1. Data acquisition:
(a) Orient all the contours in each slice in consistent
directions, as explained above. If the input does not
include this information, compute the contour nesting
hierarchy in each slice, and use it to obtain the desired
orientations.
2. Matching contour portions:
(a) Discretize each contour polygon into a cyclic sequence
of vertices, so that the arc length between each
pair of consecutive vertices is equal to some (small)
given parameter.
(b) Vote for contour matches. Each pair of distinct vertices
of the discretized contours, one on each slice, whose
mutual horizontal distance is below some threshold pa-
rameter, contributes one vote. The vote is for the match
between these two contours with the appropriate shift,
which maps one of these vertices to the other.
(c) Transform the resulting votes into a collection of
candidates of partial contour matches.
3. Reconstructing the surface:
(a) Stitch together each pair of contour portions that
have been matched in the above step, by adding triangles
which connect between these portions. The new
triangles are oriented consistently with the contours.
(b) Combine the remaining contour edges into spatial
cycles (denoted as clefts), obtained by taking the union
of the contour edges in both slices and of the edges of
the stitching triangles, and by canceling out duplicate,
oppositely-oriented edges. When projected onto the xy-
plane, these cycles do not intersect each other. Find
their nesting hierarchy, and, for each cleft C, construct,
if necessary, a system of straight 'bridges' that connect
between C and its holes (immediate children in the hi-
erarchy), so as to turn them into a single cycle, which
now replaces the cleft C.
(c) Triangulate the resulting 3-D clefts, using a 3-D min-
Figure

3: The different steps of our algorithm
imum area triangulation technique.
The various steps of the algorithm are illustrated in

Figure

3. Figure 3(a) shows a pair of slices in a branching
situation. The lower slice contains one contour (de-
noted by S 1 ), and the upper slice contains three contours
(denoted by S 2 ). Figure 3(b) shows the tiling of
the three matches found between these two slices. Figure
3(c) shows the remaining clefts, which form a shallow
hierarchy of nested polygons. Figure 3(d) shows
the clefts after the hole elimination step, and their
minimum-area triangulations are shown in Figure 3(e).
The final surface reconstruction is shown in Figure 3(f).
Matching Contour Portions
For lack of space, we omit here a detailed description
of the data acquisition step, which was already briefly
described above, and is anyway fairly straightforward.
3.1 Contour Discretization
Each contour polygon is refined and discretized into a
cyclic sequence of points. This is done by choosing some
sufficiently small arc length parameter s, and by generating
equally-spaced points, at distance s apart from
each other (along the polygon boundary). We need to
choose s so that it is much smaller than (a) the length of
any original edge of any contour, and (b) the minimum
straight distance between any pair of contour points
which lie on different contours (on the same slice) or
lie on the same contour and their distance along the
contour is sufficiently large.
3.2 Voting for Contour Matches
We first try to match pairs of portions of contours,
where each portion in a pair belongs to a different
slice. These matches aim to detect regions of similarity
between the boundaries of the interpolating object
along the two slices. Naturally, contour portions
which are similar in the two slices must have similar sequences
of footprints. Thus, our next goal is to search
for pairs of sufficiently long subsequences that closely
match each other. In our approach, two subsequences
are said to closely
match each other, if, for some chosen parameter " ? 0,
the number of indices k for which kp
is sufficiently close to '. (Here the norm k \Delta k is the
Euclidean distance between the xy-projections of the
points.) We perform the following voting process, where
votes are given to good point-to-point matches.
The contours are given as cyclic ordered sequences
of vertices, and we break each of them arbitrarily to
make it into a linear sequence. All the vertices of
the lower slice are preprocessed for range-searching, so
that, for each vertex v of the upper slice, we can efficiently
locate all the vertices of the lower slice that
lie in some "-neighborhood of (the xy-projections of) v.
We have used a simple heuristic projection method for
the range searching, that reduces the search to a pair
of 1-dimensional searches; this may be inefficient in the
worst case, but works very well in practice.
The positions along a contour sequence c, whose
length is ' c , are numbered from 0 to ' c \Gamma 1. Assume
that the querying vertex v is in position i of contour
sequence c 1 . Then, each vertex retrieved by the query,
which is in position j in contour sequence c 2 , contributes
a vote for the match between contours c 1 and c 2 with a
shift equal to (j \Gamma i) (mod ' c2 ).
Obviously, matches between long portions of contours
are reflected by a large number of votes for the appropriate
shift between the matching contours. We thus
collect the shifts with a large number of votes, identify
(roughly) the starting and ending points of the corresponding
match, and possibly also remove (small) portions
of the matches, near their endpoints, which are
also shared by other matches.
Note that this procedure will only match contour portions
that are consistently oriented. This is important,
since it guarantees that the material region lies on the
same side of the two contours, thus implying that if
we stitch the gap between these two contours, we get
a proper portion of the boundary of a possible interpolating
solid. Also note that the " parameter for the
range-searching queries is not a function of the input.
It is rather our a priori estimation of the physical size
of the difference between similar contour portions of the
two slices.
In some cases small portions of the contours are included
in more than one candidate match. This usually
happens at the connection between two different
matches, involving the same contour on one slice and
different contours on the other. We simply eliminate
those portions common to more than one significant
candidate match.
3.3 Accepting Match Candidates
Each match is given a score. The score may also reflect,
in addition to the number of votes for the appropriate
shift, other quality measures (such as the closeness of
the vertices on the two matching contour portions). Our
setting of the scoring function is described in Section 5.
A crucial property of the matches, which we have to
achieve when tuning the parameters which control the
detection of matches, is that every intersection between
the xy-projections of any pair of contours, one from each
slice, will lead to a match between the corresponding

Figure

4: Intersecting contours with no long matching
portions
contours, although it might be very short. This will ensure
that there will be at least one match between each
pair of overlapping contours (namely, contours with intersecting
xy-projections), even if the relative amount
of overlap is very small. The reason for this is given
shortly below.
We achieve this by accepting very short match candidates
(even of length 2 or 3), as long as their quality
fits our bounds. We have to make sure that the discretization
parameter is sufficiently small with respect
to the voting threshold and with respect to the smallest
contour feature size, as described above, in order not
to miss contour intersections. An extreme example is
shown in Figure 4, where four short matches are indeed
found between the two contours.
4 Reconstructing the Surface
4.1 Stitching the Matches
Each match consists of two directed polygonal chains,
whose xy-projections are very close to each other. We
arbitrarily choose one end of the match, and 'merge'
the two chains as if they were sorted lists of numbers.
In each step of the merge we have pointers to the current
vertices, u and v, in the two chains, and make a
decision as to which chain should be advanced, say v advances
to a new vertex w. Then we add the new triangle
4uvw to the boundary of the constructed polyhedron,
and advance the current vertex (from v to w) along the
appropriate chain. When we reach the last vertex of
one chain, we may further advance only the other one.
This process terminates when we reach the last vertices
of both chains.
The triangles that we create are oriented consistently
with the contours. To achieve this, we invert the orientation
of all contour polygons of one slice, say the lower
one, and orient each triangle so that the edge it shares
with a contour is oriented oppositely to the contour;
see

Figure

5 for an illustration of the directions of the
tiling triangles. This orientation guarantees that the
reconstructed solid boundary is oriented in a consistent
manner.
Several advancing rules were examined, and the following
simplest one proved itself the best. Assume
that the current vertices of the two chains are v 1
j , which are followed by v 1
Then, if jv 1
j+1 j, we advance
the first chain; otherwise we advance the second
chain 1 . This bears close resemblance to the merging
That is, we advance so that the newly added triangle has
smaller perimeter; actually, for program efficiency, we eventually
used the squares of the distances, with equally good results.

Figure

5: Tiling a match

Figure

match and the clefts near an intersection
of contours
of two sorted lists, and turns out to produce reasonably
looking triangulated boundary patches between the
matched contour portions. Figure 5 shows such a tri-
angulation. Alternative advancing rules were described
by Christiansen and Sederberg [5], and by Ganapathy
and Dennehy [10].
4.2 Filling the Clefts
After tiling the matching contour portions, we remain
with the unmatched portions. These, combined with
the extreme edges of the tiling sequences, form a collection
of closed 3-D polygons, which we refer to as clefts.
Finding these clefts is straightforward. First, recall that
we have already inverted the orientation of all contour
polygons of the lower slice. Consider the contour polygons
of the two slices (the inverted lower slice and the
upper slice), as well as the tiling triangles, as the formal
sum of their directed edges, and add up all these
polygons, with the convention that ~e
these cancellations, the resulting sum consists of all the
desired clefts. Since each polygon that participates in
the formal sum is a directed cycle, the resulting sum
is easily seen to represent a collection of pairwise edge-disjoint
directed cycles.
This collection of 3-D polygonal cleft cycles has the
property that, when projected onto the xy-plane, no
two cycles intersect. This is because we have already
detected all the projected contour intersections in the
matching phase, and because the extreme tiling edges
in the triangulation of a match are likely to degenerate
or almost degenerate in the projection. This situation
is illustrated in Figure 6. The only way this property
can be violated is when an extreme tiling edge crosses
another edge of a contour. For this to happen, either
the extreme edge must be highly slanted, or the contour
must have very sharp turns. This can be avoided
by an appropriate fine tuning of the discretization and
matching threshold parameters, and in any case we did
not face such a situation in all of our comprehensive
experiments (which involved fairly complex and rather
'adversary' data). To recap, while in rare, worst-case

Figure

7: orientations
scenarios, projected cleft cycles might intersect, we will
assume in what follows that this does not occur.
The xy-projections of these cleft polygons might
again form a hierarchy of polygon nesting. We check
this possibility by invoking again the same line-sweeping
procedure used in the data acquisition step. This time,
all the polygons are guaranteed to be correctly oriented.

Figure

7 illustrates this situation. Figure 7(a) shows a
branching case, where the lower slice contains one con-
tour, and the upper slice contains two contours. After
tiling the detected match, inverting the lower slice,
and canceling out opposite edge occurrences, two nested
cleft cycles remain, already oriented consistently, as
shown as Figure 7(b).
Let C be a cleft whose immediate children in the hierarchy
are We define an undirected weighted
complete graph G, so that each vertex of G is one of
these cycles, and the weight of an edge connecting
two cycles is the minimum distance between the
xy-projections of these cycles. (We assume that our
contour discretization is dense enough, so that the minimum
vertex-to-vertex distance, which is what we actually
computed, serves as a sufficiently good approximation
of the actual minimum distance.) We now compute
a minimum spanning tree T of G, and form a bridge between
each pair of cycles connected by an edge of T ; the
bridge is a straight segment connecting the two nearest
vertices on these cycles. It is easily verified that the
xy-projections of the bridges do not cross each other,
and also do not cross the xy-projection of any cycle.
We create two oppositely-oriented copies of each bridge
and add them to the given cycles. This eliminates the
and replaces the whole configuration
by a single composite, self-touching but otherwise
simple polygonal cycle.
We emphasize that bridge assembly was the most significant
obstacle in previous works, which used this tool
for reducing a branching situation to the simple one-
to-one case. That was because a bridge in one slice
could conflict with the geometry of the other slice, by
having its projection intersect a contour, as shown in

Figure

8(b). As noted, we do not face this problem.
We note that this procedure is required only in complicated
cases (few of them are presented in Section 5);
in most practical instances clefts do not tend to be
nested. In any case, after this hole elimination step,
we are left with a collection of closed polygonal cleft cy-
cles, with the property that their xy-projections enclose
pairwise-disjoint regions.
Our next goal is to triangulate each cleft cycle. Since
many such triangulations are possible, we seek a triangulation
which minimizes the total area of the triangles.
More precisely, we want to solve the following problem:
Given a 3-dimensional closed polygonal curve P , and
an objective function F defined on all triangles (called
weight in the sequel), find the triangulation of P (i.e., a

Figure

8: Bridges in simple branching cases
collection of triangles spanned by the vertices of P , so
that each edge of P is incident to exactly one triangle,
and all other triangle edges are incident to two triangles
which minimizes the total sum of F over its
triangles.
For this purpose, we closely follow the dynamic programming
technique of Klincsek [21] for finding a polygon
triangulation in the plane, which minimizes the total
sum of edge lengths. Let
be the given polygon. Let W i;j (0 -
denote the weight of the best triangulation of the polygonal
curve (v We then apply a
straightforward dynamic programming approach, which
initializes each W i;i+2 to F(v
computes each W (i; k), for
We omit here
the easy further details. The actual value of F(u; v; w)
that we have used is a weighted average of several terms,
the most significant one being the area of the triangle;
see Section 5 for more details.
Finally, if the interpolation involves the uppermost or
lowermost slice in the given sequence, we also have to
add to the reconstructed object boundary the regions
of material on that slice, so as to 'close' the volume of
that object.
5 Experimental Results
We have implemented the whole algorithm on a Digital
DECstation 5000/240 and on a Sun SparcStation
II in C. We have experimented with the algorithm
on several data files obtained by CT or MRI scan-
ners, and obtained very good results in practically all
cases. The input usually consisted of about fifty to
one hundred cross-sections, from 0.5 mm to 2.0 mm
apart. The tuning of the parameters (discretization
length and size of neighborhood in the geometric hash-
ing) was very robust, and large variation of these parameters
produced nearly identical results. We usually
used 1.0 mm as the discretization parameter, and
3.0 mm for the voting threshold (for human organs
whose global size was between 5 to 20 cm in all di-
mensions). We allowed up to two successive point mis-matches
along a match. A point-to-point match contributed
the amount of 1=(d + 0:1) to the match score,
where d was the xy distance between the two points.
We considered only match candidates which received
4 votes or more and whose scores were above 15:0.
(Our discretization was sufficiently dense so that this
choice still captured all contour overlaps.) The objective
function F for the cleft-triangulation was taken to

Figure

9: A synthetic example

Figure

10: A synthetic branching example
be 0:85A A is the area of the
triangle, P is its perimeter, and R is the ratio between
the largest and the smallest of its three edges. All these
parameters were user defined, but modifying them did
not achieve any better results.
We can measure the complexity of the algorithm in
terms of four variables: k, the total number of input
points along the contour edges of two consecutive slices,
n, the total number of points after the arc length discretization
step, -, the number of clefts, and h, the
size of the biggest cleft (after the hole-elimination step).
Usually, k is considerably smaller than n. The number
- of clefts could be in the worst case as large as \Theta(k 2 ),
but in practice - is smaller than, or at least comparable
with the number of contours c (which is typically much
smaller than k).
Due to lack of space, we do not give here a formal
analysis of the complexity of the algorithm. Such
an analysis produces terms that may be very large;
for example, the cleft triangulation step (based on dynamic
programming) takes O(-h 3 ) time, which is, the-
oretically, and in complex situations also pragmatically,
the most expensive portion of our algorithm. The entire
algorithm runs on practical instances in average
Here are some specific examples of the performance
of the algorithm:

Figure

9 shows a simple case, similar to the first problematic
example of Boissonnat [3]. Each slice in this
example contains exactly one contour. The tiled match
appears in white, whereas the triangulated cleft appears
in black.

Figure

shows a synthetic branching example,
where the lower slice contains two contours, whereas
the upper slice contains only one contour. Two long
matches were found and tiled by the white trian-
gles. The remaining cleft appears in between the two
matches, and its triangulation appears in black.

Figure

11 shows a more complicated synthetic exam-
ple, where the lower slice contains one contour, and the
upper slice contains three contours. Figure 11(a) shows
a top view of this situation, whereas Figure 11(b) shows
an isometric view of it. Figure 11(c) shows the surface
reconstruction. The tiling of the single match appears
in white. The remaining cleft consists of one cycle which
encloses two other cycles. Two bridge constructions
compose the three polygons into a single one, and its
(a)
(b) Before reconstruction
(c) After reconstruction

Figure

11: A synthetic complicated example
(a) Before reconstruction
(b) After reconstruction

Figure

12: A simple case
triangulation appears in black.
The next several figures show the typical performance
of our algorithm, as observed from its execution on a series
of cross-sections of a human jaw bone. Figure 12
shows the reconstructed surface between two slices. The
tiles of the matches appear in white, and the triangulations
of the two remaining clefts appear in black. Figure
13 shows two slices with a branching region. The
match tiles appear in white, whereas the cleft triangulations
appear in black. Figure 14 presents a consider-
(a) Before reconstruction
(b) After reconstruction

Figure

13: A simple branching case
(a) Before reconstruction
(b) After reconstruction

Figure

14: A complicated branching case
(a) Before reconstruction
(b) After reconstruction

Figure

15: A multiple branching case
(a) Before reconstruction
(b) After reconstruction

Figure

case
ably more complicated situation, where the reconstruction
turns out to be 'intuitively correct'. Figure 14(a)
shows that the lower slice contains a contour with two
hole contours, and the upper slice contains only one
contour, which mostly lies above the "material" region
of the lower slice. The surface reconstruction is shown
in

Figure

14(b). The three match tiles appear in white,
and the cleft triangulations appear in black. The reader
may verify that the two "non-material pillars" represented
by the two hole contours in the lower slice were
connected to the unbounded "non-material" region in
the upper slice. Similar complex examples are given in

Figures

15 and 16.
The reconstruction of the whole jaw bone, whose in-
Figure

17: A fully reconstructed human jaw bone
(a) (b)
(c) (d)

Figure

put consisted of 96 slices, is shown in Figure 17. The
result was a valid polyhedral description, which contained
about 60,000 triangles. The reconstructed jaw
contained 209 3-D cavities (fully enclosed in the outer
3-D reconstructed boundary).
Our experimentations were mostly performed on
medical imaging data. However, as we mentioned in
the introduction, the reconstruction problem has other
applications as well. One such application is the re-construction
of a terrain from elevation contour data in
topographic maps. This application is much simpler,
because of the xy-monotonicity of terrains. For exam-
ple, contours cannot overlap in this case (although they
can be nested). Figure 18 shows the reconstruction of
a terrain in the Zikhron-Ya'akov area in Israel. Figures
show the elevation contours (a top view and
an isometric view). The data contained seventeen levels
which are 10 meters apart, starting from 20 and ending
at meters above see level. These levels consisted of
Model Synthetic Jaw Bone
Slices 2 96
Layers 1 95
Contours:
Total 3.0 417.0
Per Slice 1.5 4.3
Vertices:
Total 143.0 17,646.0
Per Contour 47.7 42.3
Matches:
Total 2.0 330.0
Per Layer 2.0 3.5
Clefts:
Total 1.0 275.0
Per Layer 1.0 2.9
Time (Seconds):
All:
Stitching 0.07 12.50
Triangulation 0.13 37.12
Total
Per Layer 0.20 0.52

Table

1: Performance of the algorithm
contours made of 5,329 edges. Figures 18(c,d) show
the full reconstruction of the terrain (again, a top view
and an isometric view).
Finally, we measured the performance of our implementation
on two of the examples described above. All
the time measurements were taken on a Digital DECstation
5000/240. Table 1 summarizes the performance
of the algorithm on the simple branching case shown in

Figure

10, and on the whole jaw bone shown in Figure
17. Note that the data of the jaw bone was accumulated
from 95 layer reconstructions. Therefore, some
of the data is averaged per layer. We note that in medical
data, at least the data with which we experimented,
successive slices tend to differ a lot in their geometries.
Therefore, many clefts are created and their (relatively
time-consuming) triangulations affect the total running
time of the algorithm.
6 Conclusion
We have proposed in this paper an algorithm for solving
the practical problem of polyhedral interpolation between
parallel polygonal slices. This problem has many
medical and geographic applications, and also appears
to be a fairly basic and interesting problem in computer
graphics and solid modeling.
Our method produces a relatively smooth boundary,
due to the contour discretization. In situations where
the addition of new vertices is not desired, e.g. due to
data explosion, our system uses the discretization only
for the matching step, but the tiling itself and the following
minimum-area triangulation are performed on
the contours containing only the original points.
We feel that our technique reconstructed the boundary
of various organs in an intuitively appealing manner;
one might say that our algorithm demonstrated some
'understanding' of the underlying problem. The results
were more than adequate even in extreme cases of tiling
between two seemingly totally different slices.
We plan to continue the experimentation with our al-
gorithm, to test its performance limits and see if there
are data instances on which the algorithm might not
perform well, thus requiring further calibrations and enhancements
(in view of our experimentation so far, we
doubt that anything really problematic will arise).

Acknowledgment

We wish to thank Haim Wolfson for helpful discussions
concerning the geometric hashing technique, Emo Welzl
and Barbara Wolfers for helpful discussions concerning
the solid reconstruction problem, and for providing us
with data files and with a comprehensive bibliography
on that problem, and Jean-Daniel Boissonnat and Bernhard
Geiger for supplying us with data files.



--R

Filling gaps in the boundary of a polyhedron

Shape reconstruction from planar cross sec- tions
Three dimensional reconstruction of complex shapes based on the Delaunay triangula- tion
Conversion of complex contour line definitions into polygonal element mo- saics
Two algorithms for the three-dimensional reconstruction of tomograms

A triangulation algorithm from arbitrary shaped multiple planar contours
Optimal surface reconstruction from planar contours
A new general triangulation method for planar contours
Construction et utilisation des mod'eles d'organes en vue de l'assistance au diagnostic et aux interventions chirur
On reconstructing polyhedra from parallel slices
Primitives for the manipulation of general subdivisions and the computation of Voronoi diagrams

An improved model-based matching method using footprints
boundary matching using footprints
A framework for surface reconstruction from 3D contours
A syntactic/semantic technique for surface reconstruction from cross-sectional contours
Approximating complex surfaces by triangulation of contour lines

Minimal triangulations of polygonal do- mains
Surfaces from con- tours: the correspondence and branching problems
Reconstructing 3D objects from cross- sections
of partially obscured objects in two and three dimensions by matching noisy characteristic curves
Surface definition for branching contour-defined objects
Surface reconstruction from sparse data
From contours to surfaces: testbed and initial results
Pessimal guesses may be opti- mal: A counterintuitive search result
Surface reconstructionand representation of 3D scenes
Surface reconstruction between simple polygons via angle criteria
IEEE Transactions on Pattern Analysis and Machine Intelligence
Surface construction from planar contours
--TR
Surface reconstruction and representation of 3-D scenes
From contours to surfaces: testbed and initial results
Two-dimensional, model-based, boundary matching using footprints
of partially obscured objects in two and three dimensions by matching noisy characteristic
A syntactic/semantic technique for surface reconstruction from cross-sectional contours
A framework for surface reconstruction from 3D contours
Shape reconstruction from planar cross sections
Pessimal Guesses may be Optimal
On curve matching
A triangulation algorithm from arbitrary shaped multiple planar contours
Primitives for the manipulation of general subdivisions and the computation of Voronoi
Optimal surface reconstruction from planar contours
Surface Reconstruction Between Simple Polygons via Angle Criteria

--CTR
Takeo Igarashi , Satoshi Matsuoka , Hidehiko Tanaka, Teddy: a sketching interface for 3D freeform design, Proceedings of the 26th annual conference on Computer graphics and interactive techniques, p.409-416, July 1999
Matthew T. Dickerson , Scott A. McElfresh , Mark Montague, New algorithms and empirical findings on minimum weight triangulation heuristics (extended abstract), Proceedings of the eleventh annual symposium on Computational geometry, p.238-247, June 05-07, 1995, Vancouver, British Columbia, Canada
Takeo Igarashi , Satoshi Matsuoka , Hidehiko Tanaka, Teddy: a sketching interface for 3D freeform design, ACM SIGGRAPH 2006 Courses, July 30-August 03, 2006, Boston, Massachusetts
Gill Barequet , Daniel Shapiro , Ayellet Tal, History consideration in reconstructing polyhedral surfaces from parallel slices, Proceedings of the 7th conference on Visualization '96, p.149-ff., October 28-29, 1996, San Francisco, California, United States
Daniel Cohen-Or , Amira Solomovic , David Levin, Three-dimensional distance field metamorphosis, ACM Transactions on Graphics (TOG), v.17 n.2, p.116-141, April 1998
James D. Fix , Richard E. Ladner, Multiresolution banded refinement to accelerate surface reconstruction from polygons, Proceedings of the fourteenth annual symposium on Computational geometry, p.240-248, June 07-10, 1998, Minneapolis, Minnesota, United States
F. Hurtado , M. Noy , J. Urrutia, Flipping edges in triangulations, Proceedings of the twelfth annual symposium on Computational geometry, p.214-223, May 24-26, 1996, Philadelphia, Pennsylvania, United States
Nina Amenta , Marshall Bern, Surface reconstruction by Voronoi filtering, Proceedings of the fourteenth annual symposium on Computational geometry, p.39-48, June 07-10, 1998, Minneapolis, Minnesota, United States
slice interpolation;branching surfaces;surface fitting;polyhedra;geometric hashing;curve matching;tiling;dynamic programming;surface reconstruction;triangulation