Surface Approximation and Geometric Partitions.
Motivated by applications in computer graphics, visualization, and scientific computation, we study the computational complexity of the following problem: given a set S of n points sampled from a bivariate function f(x,y) and an input parameter $\eps > 0$, compute a piecewise-linear function $\Sigma(x,y)$ of minimum complexity (that is, an xy-monotone polyhedral surface, with a minimum number of vertices, edges, or faces) such that $| \Sigma(x_p, y_p) \; - \; z_p | \:\:\leq\:\: \eps$ for all $(x_p, y_p, z_p) \in S$. We give hardness evidence for this problem, by showing that a closely related problem is NP-hard. The main result of our paper is a polynomial-time approximation algorithm that computes a piecewise-linear surface of size O(Ko log Ko), where Ko is the complexity of an optimal surface satisfying the constraints of the problem.The technique developed in our paper is more general and applies to several other problems that deal with partitioning of points (or other objects) subject to certain geometric constraints.  For instance, we get the same approximation bound for the following problem arising in machine learning: given n "red" and m "blue" points in the plane, find a minimum number of pairwise disjoint triangles such that each blue point is covered by some triangle and no red point lies in any of the triangles.
Introduction
In scientific computation, visualization, and computer graphics, the modeling and construction
of surfaces is an important area. A small sample of some recent papers [2, 3, 5, 7, 10,
13, 20, 21] on this topic gives an indication of the scope and importance of this problem.
The first author has been supported by National Science Foundation Grant CCR-93-01259 and an NYI
award.
Rather than delve into any specific problem studied in these papers, we focus on a general,
abstract problem that seems to underlie them all.
In many scientific and computer graphics applications, computation takes place over
a surface in three dimensions. The surface is generally modeled by piecewise linear (or,
sometimes piecewise cubic) patches, whose vertices lie either on or in the close vicinity of
the actual surface. In order to ensure that all local features of the surface are captured,
algorithms for an automatic generation of these models generally sample at a dense set of
regularly spaced points. Demands for real-time speed and reasonable performance, however,
require the models to have as small a combinatorial complexity as possible. A common
technique employed to reduce the complexity of the model is to somehow "thin" the surface
by deleting vertices with relatively "flat" neighborhoods. Only ad hoc and heuristic methods
are known for this key step. Most of the thinning methods follow a set of local rules (such
as deleting edges or vertices whose incident faces are almost coplanar), which are applied
in an arbitrary order until they are no longer applicable. Not surprisingly, these methods
come with no performance guarantee, and generally no quantitative statement can be made
about the surface approximation computed by them.
In this paper, we address the complexity issues of the surface approximation problem
for surfaces that are xy-monotone. These surfaces represent the graphs of bivariate functions
f(x; y), and they arise quite naturally in many scientific applications. One possible
approach for handling arbitrary surfaces is to break them into monotone pieces, and apply
our algorithm individually on each piece. Let us formally define the main problem studied
in our paper.
Let f be a function of two variables x and y, and let S be a set of n points sampled
from f . A piecewise linear function \Sigma is called an "-approximation of f , for " ? 0, if
for every point Given S and ", the surface approximation problem is
to compute a piecewise linear function \Sigma that "-approximates f with a minimum number
of breakpoints. The breakpoints of \Sigma can be arbitrary points of R 3 , and they are not
necessarily points of S. In many applications, f is generally a function hypothesized to
fit the observed data-the function \Sigma is a computationally efficient substitute for f . The
parameter " is used to achieve a complexity-quality tradeoff - smaller the ", higher the
fidelity of the approximation. (The graph of a piecewise linear function of two variables is
also called a polyhedral terrain in computational geometry literature; the breakpoints of the
function are the vertices of the terrain.)
The state of theoretical knowledge on the surface approximation problem appears to be
rather slim. The provable performance bounds are known only for convex surfaces. For this
special case, an O(n 3 presented by Mitchell and Suri [16] for computing
an "-approximation of a convex polytope of n vertices in R 3 . Their algorithm produces
an approximation of size O(K o log n), where K o is the size of an optimal "-approximation.
Extending their work, Clarkson [6] has proposed an O(K randomized
algorithm for computing an approximation of size O(K o log K can be an
arbitrarily small positive number.
In this paper, we study the approximation problem for surfaces that correspond to
graphs of bivariate functions. We show that it is NP-Hard to decide if a surface can be
"-approximated using k vertices (or facets). The main result of our paper, however, is a
polynomial-time algorithm for computing an "-approximating surface of a guaranteed qual-
ity. If an optimal "-approximating surface of f has K our algorithm produces
a surface with O(K vertices. Observe that we are dealing with two approximation
measures here: ", which measures the absolute z difference between f and the "simplified"
surface \Sigma, and the ratio between the sizes of the optimal surface and the output of our
algorithm. For the lack of better terminology, we use the term "approximation" for both
measures. Notice though that " is an input (user-specified) parameter, while log K o is the
approximation guarantee provided by the analysis of our algorithm.
The key to our approximation method is an algorithm for partitioning a set of points
in the plane by pairwise disjoint triangles. This is an instance of the geometric set cover
problem, with an additional disjointness constraint on the covering objects (triangles).
Observe that the disjointness condition on covering objects precludes the well-known greedy
method for set covering [11, 14]; in fact, we can show that a greedy solution has size
in the worst-case. Let us now reformulate our surface approximation problem as a
constrained geometric partitioning problem.
p denote the orthogonal projection of a point p 2 R 3 onto the xy-plane z = 0. In
general, for any set A ae R 3 , we use -
A to denote the orthogonal projection of A onto the
xy-plane. Then, in order to get an "-approximation of f , it suffices to find a set of triangles
in 3-space such that (i) the projections of these triangles on plane z = 0 are pairwise disjoint
and they cover the projected set of points -
S, and (ii) the vertical distance between a triangle
and any point of S directly above/below it is at most ". Our polynomial-time algorithm
produces a family of O(K We stitch together these triangles to
produce the desired surface \Sigma. The "stitching" process introduces at most a constant factor
more triangles.
The geometric partition framework also includes several extensions and generalizations
of the basic surface approximation problem. For instance, we can formulate a stronger
version of the problem by replacing each sample point by a horizontal triangle (or, any
polygon). Specifically, we are given a family of horizontal triangles (or polygons) in 3-
space, whose projections on the xy-plane are pairwise disjoint. We want a piecewise linear,
"-approximating surface whose maximum vertical distance from any point on the triangles
is ". Our approximation algorithm works equally well for this variant of the problem-this
variant addresses the case when some local features of the surface are known in detail;
unfortunately, our method works only for horizontal features.
Finally, let us mention the planar bichromatic partition problem that is of independent
interest in the machine learning literature: Given a set R of 'red' points and another set
B of 'blue' points in the plane, find a minimum number of pairwise disjoint triangles such
that each blue point lies in a triangle and no red point lies in any of the triangles. Our
algorithm gives a solution with O(K
The running time of our algorithms, though polynomial, is quite high, and at the moment
has only theoretical value. These being some of the first results in this area, however, we
expect that the theoretical time complexity of these problems would improve with further
work. Perhaps, some of the ideas in our paper may also shed light on the theoretical
performance of some of the "practical" algorithms that are used in the trade.
A Proof of NP-Hardness
We show that both the planar bichromatic partition problem and the surface approximation
problem are NP-Hard , by a reduction from the planar 3-SAT. We do not know whether
they are in NP since the coordinates of the triangles in the solution may be quite large. We
recall that the 3-SAT problem consists of n variables x clauses
each with three literals C
x k . The problem is to decide
whether the boolean formula
has a truth assignment. An instance of 3-SAT is called planar if its variable-clause graph
is planar. In other words, F is an instance of the planar 3-SAT if the graph G(F E)
is planar (see [12]), where V and E are defined as follows:
appears in C
Theorem 2.1 The planar bichromatic partition problem is NP-Hard.
Proof: Our construction is similar to the one used by Fowler et al. [9], who prove the
intractability of certain planar geometric covering problems (without the disjointness con-
dition); see also [4, 8] for similar constructions. We first describe our construction for the
bichromatic partition problem. To simplify the proof, our construction allows three or more
points to lie on a line-the construction can be modified easily to remove these degeneracies.
Let F be a boolean formula, and let E) be a straight-line planar embedding of
the graph G(F ). We construct an instance of the red-blue point partition problem whose
solution determines whether F is satisfiable.
a 11
(iv)

Figure

1: (i) An instance of planar 3-SAT:
instance of bichromatic partition, (iii) Details of P 1 and C 1 , only some of the red
points lying near P 1 and C 1 are shown, (iv) Two possible coverings of blue points on P 2 .
We start by placing a 'blue' point at each clause node C m. Let x i be a
variable node, and let e i1 il be the edges incident to it. In the plane embedding of
G, the edges e ij form a "star" (see Figure 1 (i)). We replace this star by its Minkowski sum
with a disk of radius ffi , for a sufficiently small ffi ? 0. Before performing the Minkowski sum,
however, we shrink all the edges of the star at x i by 2ffi, so that the "star-shaped polygons"
meeting at a clause node do not overlap (see Figure 1 (ii)). Let P i denote the star-shaped
polygon corresponding to x i . In the polygon P i , corresponding to each edge e ij , there is a
tube, consisting of two copies of e ij , each translated by distance ffi, plus a circular arc s ij
near the clause node C j .
We place an even number of (say 2k i ) 'blue' points on the boundary of P i , as follows.
We put two points a ij and b ij on the circular arc s ij near its tip. If C j contains the literal
six points on the straight-line portion of P i 's boundary, three each on translated
copies of the edge e ij . On each copy, we move the middle point slightly inwards so as to
replace the original edge of P i by a path of length two. On the other hand, if C j contains
the literal -
four points on straight-line portion of P i 's boundary, two each on
translated copies of the edge e ij . Thus, the number of blue points added for edge e ij is
either six of eight. (2k i is the total number of points put along P i .) Let B denote the set
of all blue points placed in this way, and let
Finally, we scatter a large (but polynomially bounded) number of 'red' points so that (i)
any segment connecting two blue points that are not adjacent along the boundary of some
contains a red point, and (ii) any triangle with three blue points as its vertices contains
at least one red point unless the triangle is defined by a

Figure

(iii).) Let R be a set of red points satisfying the above two properties.
We claim that the set of blue points B can be covered by k disjoint triangles, none of
which contains any red point, if and only if the formula F has a truth assignment. Our
proof is similar to the one in Fowler et al. [9]; we only sketch the main ideas. The red points
act as enforcers, ensuring that only those blue points that are adjacent on the boundary of
a P i can be covered by a single triangle. Thus, the minimum number of triangles needed
to cover all the points on P i is k i . Further, there are precisely two ways to cover these
points using k i triangles- in one covering, a ij and b ij are covered by a single triangle for
those clauses only in which x and, in the other covering, a ij and b ij are covered by
a single triangle for those clauses only in which -
Figure 4 (iv). We regard the
first covering as setting x and the second covering as setting x
Suppose 1. For any clause C j that contains x i , the points a ij and b ij are covered
by a single triangle, and we can cover the clause point corresponding to C j by the same
triangle. The same holds if x and the clause C j contains -
In other words, the clause
point of C j can be covered for free if C j is satisfied. Thus, the set of blue points B can be
covered by k triangles if and only if the clause point for each clause C j is covered for free,
that is, the formula F has a truth assignment. This completes our proof of NP-Hardness
of the planar bichromatic partition problem. 2
Remark: The preceding construction is degenerate in that most of the red points lie on
segments connecting two blue points. There are several ways to remove these collinearities;
we briefly describe one of them. For each polygon P i , replace every blue point b on P i by
two blue points b placed very close to b. (We do not make copies of 'clause points' C j ,
m.) For every pair of blue points b l that we did not want to cover by a single
triangle in the original construction, we place a red point in the convex hull of b 0
l .
If there are 4k i blue points on the boundary of P i , they can be covered by k i triangles, and
there are exactly two ways to cover these blue points by k i triangles, as earlier. Following
a similar, but more involved, argument, we can prove that the set of all blue points can be
covered by
triangles if and only if F is satisfiable.
Theorem 2.2 The 3-dimensional surface approximation problem is NP-Hard.
Proof: Our construction is similar in spirit to the one for the bichromatic partition problem,
albeit slightly more complex in detail. We use points of three colors: red, white and black.
The 'white' points lie on the plane z = 0, the `black' points lie on the plane z = 2A, and
the 'red' points lie between A is a sufficiently large constant.
To maintain a connection with the previous construction, the black and white points play
the role of blue points, while the red points play the role of enforcers as before, restricting
the choice of "legal" triangles that can cover the black or white points. We will describe
the construction in the xy-plane, which represents the orthogonal projection of the actual
construction. The actual construction is obtained simply by vertically translating each
point to its correct plane, depending on its color.
We start out again by putting a 'black' point at each clause node C j . Then, for each
variable x i , we construct the "star-shaped" polygon P i ; this part is identical to the previous
construction. We replace each of the two straight-line edges of P i by "concave chains," bent
inward, and also make a small "dent" at the tip of the circular arc s ij , as shown in Figure 2.
We place 12 points on each arm of P i , alternating in color black and white, as follows. At
the tip of the circular arc s ij , we put a white point c ij at the outer endpoint of the dent
and a black point d ij at the inner endpoint of the dent (Figure 2 (ii)). The rest of the
construction is shown in Figure 2 (i) - we put two more points a ij on the circular
arc and 4 points ff l
on each of the two concave chains. The two points
surrounding namely, a ij and b ij , are such that any segment connecting them to any
point on the two concave chains lies inside P i . Next, corresponding to each edge e ij of the
graph G(F ), we put a 'white' point c 0
ij on the segment joining c ij and the clause point C j ,
very close to c ij such that
(See

Figure

2 (iii).) This condition says that, in the final construction when the black
and white points have been translated to their correct z-plane, the vertical distance between
ij and the segment C no more than "-recall that " is the input measure of
approximation.
This completes the placement of white and black points. The only remaining part of
the construction is the placement of 'red' points, which we now describe.
a 22
a
d 21
c 0c 0d 22
c 22
e
c 0(iii)
c 0b21
c 21
ff 1a 21

Figure

2: Placing points on the polygon P 2 corresponding to Figure 1: (i) Modified P 2
and points on P 2 , (ii) points and triangles in the neighborhood of c points and
triangles near C 1 .
We add a set of triangles, each containing a large (but polynomially bounded) number of
red points-the role of these triangles is to restrict the choice of legal triangles that can cover
black/white points. The set of triangles associated with P i is labeled T i . The construction
of T i is detailed in Figure 2 (i). Specifically, for an edge e ij , if x then we put a small
triangle that intersects the segment b ij c 0
ij but not b ij c ij . On the other hand, if -
then we put a small triangle that intersects a ij c 0
ij but not a ij c ij . Next, we put a small
number of triangles inside P i , near its concave chains, so that at most three consecutive
points along P i may be covered by one triangle without intersecting any triangle of T i . We
ensure that one of these triangles intersects the triangle 4ff 1
so that fff 1
cannot be covered by a single triangle. We also place three triangles near each clause C j ,
each containing a large number of red points; see Figure 2 (iii). Finally, we translate black
and white points in the z-direction, as described earlier. Let f- be the set of all
'red' triangles. We move all points in - i vertically to the plane z
There are
two ways to cover the points of P i with 2k i legal (non-intersecting) triangles-one in which
are covered by a single triangle, and the one in which b ij are covered by a
triangle. These coverings are associated with the true and false settings of the variable x i .
Let P denote the set of all points constructed, let t denote the total number of 'red'
triangles, and let
We claim that there exist a polyhedral terrain with
3(k+m+t) vertices that "-approximates P if and only if F has a truth assignment, provided
that " is sufficiently small-recall that m is the number of clauses in F . The claim follows
from the observations that it is always better to cover all red points lying in a horizontal
triangle - i by - i itself, and that a clause C j requires one triangle to cover its points if and
only if one of the literals in C j is set true; otherwise it requires two triangles. (For instance,
if C j contains the literal x i and x i is set true, then the triangle a ij can be enlarged
slightly to cover c 0
ij . The remaining three points for the clause C j can be covered by one
additional triangle.) The rest of the argument is the same as for the bichromatic partition
problem. Finally, we can perturb the points slightly so that no four of them are coplanar.In the remainder of the paper, we develop our approximation algorithms.
3 A Canonical Trapezoidal Partition
We introduce an abstract geometric partitioning problem in the plane, which captures
the essence of both the surface approximation problem as well as the bichromatic points
partition problem. The abstract problem deals with trapezoidal partitions under a boolean
constraint function C satisfying the "subset restriction" property. More precisely, let C be
a boolean function from compact, connected subsets of the plane to f0; 1g satisfying the
following property:
For technical reasons, we choose to work with "trapezoids" instead of triangles, where
the top and bottom edges of each trapezoid are parallel to the X-axis. The trapezoids and
triangles are equivalent for the purpose of approximation-each triangle can be decomposed
into two trapezoids, and each trapezoid can be decomposed into two triangles.
Given a set of n points P in the plane, a family of trapezoids is
called a valid trapezoidal partition (a trapezoidal partition for brevity) of P with respect to
a boolean constraint function C if the following conditions hold:
(ii) \Delta covers all the points: P ae
(iii) The trapezoids in \Delta have pairwise disjoint interiors.
We can cast our bichromatic partition problem in this abstract framework by setting
set of 'blue' points) and, for a trapezoid - ae R 2 , defining only
if - is empty of red points, that is, - In the surface approximation problem, we
set
(the orthogonal projection of S on the plane z = 0) and a trapezoid - ae R 2 has
only if - can be vertically lifted to a planar trapezoid - in R 3 so that the
vertical distance between -
- and any point of S directly above/below it is at most ".
The space of optimal solutions for our abstract problem is potentially infinite-the
vertices of the triangles in our problem can be anywhere in the plane. For our approximation
results, however, we show that a restricted choice of trapezoids suffices.
Given a set of n points P in the plane, let L(P ) denote the set consisting of the following
lines: the horizontal lines passing through a point of P , and the lines passing through two
points of P . Thus, jL(P We will call the lines of L(P ) the canonical lines
determined by P . We say that a trapezoid \Delta ae R 2 is canonical if all of its edges belong to
lines in L(P ). A trapezoidal partition \Delta is canonical if all of its trapezoids are canonical.
The following lemma shows that by limiting ourselves to canonical trapezoidal partitions
only, we sacrifice at most a constant (multiplicative) factor in our approximation.

Figure

3: A canonical trapezoidal partition
Lemma 3.1 Any trapezoidal partition of P with k trapezoids can be transformed into a
canonical trapezoidal partition of P with at most 4k trapezoids.
Proof: We give a construction for transforming each trapezoid \Delta 2 \Delta into four trapezoids
4, with pairwise disjoint interiors, so that \Delta i together cover all the
points in P " \Delta. By (3.1), the new set of trapezoids is a valid trapezoidal partition of P .
Our construction works as follows.
Consider the convex hull of the points P "\Delta. If the convex hull itself is a trapezoid,
we return that trapezoid. Otherwise, let '; denote the left, right, top and bottom edges
of \Delta, as shown in Figure 4 (i). We perform the following four steps, which constitute our
transformation algorithm.
r
(iv)
r

Figure

4: Illustration of the canonicalization.
(i) We shrink the trapezoid \Delta by translating each of its four bounding edges towards the
interior, until it meets a point of P \Delta . Let \Delta 0 ' \Delta denote the smaller trapezoid thus
obtained respectively, denote a point of P \Delta lying on
the left, right, top, and bottom edge of \Delta 0 ; we break ties arbitrarily if more than one
point lies on an edge.
(ii) We partition \Delta 0 into two trapezoids, \Delta L and \Delta R , by drawing the line segment p u
as shown in (Figure 4 (ii).
(iii) We next partition \Delta L into two trapezoids \Delta LU and \Delta LB , by drawing the maximal
horizontal segment through p ' . Let p 0
' denote right endpoint of this segment. Similarly,
we partition \Delta R into \Delta RU and \Delta RB , lying respectively above and below the horizontal
line segment p r p 0
r .
(iv) We rotate the line supporting the left boundary of \Delta LU around the point p ' in clock-wise
direction until it meets a point of the set denote the intersection
of this line and the top edge of \Delta LU . We set

Figure

4 (iv)). (If
a triangle, which we regard as a degenerate trapezoid; e.g. \Delta 4 in

Figure

4 (iii).) The top and bottom edges of \Delta 1 contain p u and p ' , respectively, the
left edge contains p ' and q ' , and the right edge is determined by the pair of points
p u and p b . Thus, the trapezoid \Delta 1 is in canonical position. The three remaining
trapezoids are constructed similarly.
In the above construction, if any of the four trapezoids \Delta i does not cover any point
of P \Delta , then we can discard it. Thus, an arbitrary trapezoid of the partition \Delta can be
transformed into at most four canonical trapezoids. This completes the proof of the lemma.4 Greedy Algorithms
At this point, we can obtain a weak approximation result using the canonical trapezoidal
partition. Roughly speaking, we can use the greedy set covering heuristic [6, 14], ignoring the
disjointness constraint, and then refine the heuristic output to produce disjoint trapezoids.
Unfortunately, the last step can increase the complexity of the solution quite significantly.
Theorem 4.1 Given a set P of n points in the plane and a boolean constraint function C
satisfying (3.1) that can be evaluated in polynomial time, we can compute an O(K
log K
size trapezoidal partition of P respecting C in polynomial time, where K o is the size of an
optimal trapezoidal partition.
Proof: Consider the set \Xi of all valid, canonical trapezoids in the plane-the set \Xi has
O(n 6 ) trapezoids. We form the geometric set-system
X can be computed by testing each \Delta 2 \Xi whether it is valid. We compute a set cover of
X using the greedy algorithm [11, 14] in polynomial time. The greedy algorithm returns a
set R consisting of O(K necessarily disjoint. In order to produce a
disjoint cover, we first compute the arrangement A(R) of the plane induced by R. Then,
we decompose each face of A(R) into trapezoids by drawing a horizontal segment through
each vertex until the segment hits an edge of the arrangement. The resulting partition is
a trapezoidal partition of P . The number of trapezoids in the partition is O((K
- since the arrangement A(R) has this size. The total running time of the algorithm is
polynomial. 2
Remark: (i) The canonical form of trapezoids is used only to construct a finite family of
trapezoids to search for an approximate solution. A direct application of the definition in
the previous subsection gives a family of O(n 6 ) canonical trapezoids. By using a slightly
different canonical form, we can reduce the size of canonical triangles to O(n 4 ). In another
paper [1], we present an near-quadratic time algorithm for finding an approximation of size
(ii) One can show that the number of trapezoids produced by the above algorithm is \Omega\Gamma K
in the worst case.
A Recursively Separable Partition 13
5 A Recursively Separable Partition
We now develop our main approximation algorithm. The algorithm is based on dynamic
programming, and it depends on two key ideas-a recursively separable partition and a
compliant partition. These partitions are further specializations of the canonical trapezoidal
partition introduced in the previous section, and they are central to our algorithm's
performance.
A trapezoidal partition \Delta is called recursively separable if the following holds:
ffl \Delta consists of a single trapezoid, or
ffl there exists a line ' not intersecting the interior of any trapezoid in \Delta such that (i)
are both nonempty, where ' are the two
half-planes defined by ', and (ii) each of \Delta recursively separable.
The following key lemma gives an upper bound on the penalty incurred by our approximation
algorithm if only recursively separable trapezoidal partitions are used.
Lemma 5.1 Let P be a finite set of points in the plane and let \Delta be a trapezoidal partition
of P with k trapezoids. There exists a recursively separable partition \Delta ? of P with O(k log
trapezoids. In addition, each separating line is either a horizontal line passing through a
vertex of \Delta or a line supporting an edge of a trapezoid in \Delta.
Proof: We present a recursive algorithm for computing \Delta ? . Our algorithm is similar to
the binary space partition algorithm proposed by Paterson and Yao [17]. We assume that
the boundaries of the trapezoids in \Delta are also pairwise disjoint-this assumption is only to
simplify our proof.
At each recursive step of the algorithm, the subproblem under consideration lies in a
trapezoid T . (This containing trapezoid may degenerate to a triangle, or it may even be
unbounded.) The top and bottom edges of T (if they exist) pass through the vertices of \Delta,
while the left and right edges (if they exist) are portions of edges of \Delta. Initially T is a set to
an appropriately large trapezoid containing the family \Delta. Let \Delta T denote the trapezoidal
partition of P " T obtained by intersecting \Delta with T , and let V T be the set of vertices of
lying in the interior of T . An edge of \Delta T cannot intersect the left or right edge of T ,
because they are portions of the edge of T . Therefore, each edge of \Delta T either lies in the
interior of T , or intersects only the top and bottom edges of T .
If stop. Otherwise, we proceed as follows. If there is
a trapezoid 4 2 \Delta T that completely crosses T (that is, its vertices lie on the top and
bottom edges of T ), then we do the following. If 4 is the leftmost trapezoid of \Delta T , then
we partition T into two trapezoids by drawing a line through the right edge of \Delta,
so that T 1 contains 4 and T 2 contains the remaining trapezoids of \Delta T . If 4 is not the
leftmost trapezoid of \Delta T , then we partition T into by drawing a line through the left
edge of \Delta.
If every trapezoid \Delta in \Delta T has at least one vertex in the interior of T , and so the
previous condition is not met, then we choose a point with a median y-coordinate.
We partition T into trapezoids T by drawing a horizontal line ' v passing through v.
Each trapezoid partitioned into two trapezoids by adding the
segment . At the end of this dividing step, let \Delta 1 and \Delta 2 be the set of trapezoids of
that lie in T 1 and T 2 , respectively. We recursively refine \Delta 1 and \Delta 2 into separable partitions
respectively, and return \Delta ?
2 . This completes the description of the
algorithm.
We now prove that \Delta   satisfies the properties claimed in the lemma. It is clear that \Delta
is recursively separable and that each separating line of \Delta   either supports an edge of \Delta or
it is horizontal. To bound the size of \Delta   , we charge each trapezoid of \Delta   to its bottom-left
vertex. Each such vertex is either a bottom-left vertex of a trapezoid of \Delta, or it is an
intersection point of a left edge of a trapezoid of \Delta with the extension of a horizontal edge
of another trapezoid of \Delta. There are only k vertices of the first type, so it suffices to bound
the number of vertices of the second type. Since the algorithm extends a horizontal edge
of a trapezoid of \Delta T only if every trapezoid of \Delta T has at least one vertex in the interior
of T , and we always extend a horizontal edge with a median y-coordinate, it is easily seen
that the number of vertices of the second type is O(k log k). This completes the proof. 2
Remark: Given a family \Delta of k disjoint orthogonal rectangles partitioning P , we can find
a set of O(k) recursively separable rectangles that forms a rectangular partition of P -this
uses the orthogonal binary space partition algorithm of Paterson and Yao [18].
6 An Approximation Algorithm
Lemma 5.1 applies to any trapezoidal partition of P . In particular, if we start with a
canonical trapezoidal partition \Delta, then the output partition \Delta ? is both canonical and
recursively separable, and each separating line in \Delta ? belongs to the family of canonical
lines L(P ). For the lack of a better term, we call a trapezoidal partition of P that satisfies
these conditions a compliant partition. Lemmas 3.1 and 5.1 together imply the following
useful theorem.
Theorem 6.1 Let P be a set of n points in the plane and let C be a boolean constraint
function satisfying the condition (3.1). If there is a trapezoidal partition of P respecting C
with k trapezoids, then there is a compliant partition of P also respecting C with O(k log
trapezoids.
In the remainder of this section, we give a polynomial-time algorithm, using dynamic
programming, for constructing an optimal compliant partition. By Theorem 6.1, this partition
has O(K trapezoids. Recall that the set consists of all canonical
lines determined by P .
Consider a subset of points R ' P , and a canonical trapezoid \Delta containing R. Let
oe(R; \Delta) denote the size of an optimal compliant partition of R in \Delta; the size of a partition
is the number of trapezoids in the partition. Theorem 6.1 gives the following recursive
definition of oe:
min
where the minimum is over all those lines ' 2 L that are either horizontal and intersect
\Delta, or intersect both the top and bottom edges of \Delta; ' denote the positive and
negative half-planes induced by . The goal of our dynamic
programming algorithm is to compute oe(P; T ), for some canonical trapezoid T enclosing
all the points P . We now describe how the dynamic programming systematically computes
the required partial answers.
Every canonical trapezoid \Delta in the plane can be described (uniquely) by a 6-tuple
consisting of integers between 1 and n. The first two numbers fix two
points p i and p j through which the lines containing the top and bottom edges of \Delta pass;
the second pair fixes the points p k 1
through which the line containing the left edge of
passes; and the third pair fixes the points p l 1
through which the line containing the
right edge of \Delta passes. (In case of ties, we may choose the points closest to the corners of
\Delta.) We use the notation \Delta(i; for the trapezoid associated with the 6-tuple
If the 6-tuple does not produce a trapezoid, then \Delta(i;
undefined.
We use the abbreviated notation
to denote the size of an optimal compliant partition for the points contained
in \Delta(i;
The quantity oe(i; undefined if the trapezoid \Delta(i;
If the points in P are sorted in increasing order of their y-coordinates, then \Delta(i;
is defined only for i - j. Our dynamic programming algorithm computes the oe values as
follows.
If C (\Delta(i;
Otherwise,
min
where the last minimum varies over all pairs of points such that the line passing
through them intersects both the top and the bottom edge of \Delta(i;
If \Xi denotes the set of all canonical trapezoids, then )-each 6-tuple is
associated with at most one unique trapezoid. If Q(n) denotes the time to decide whether
for an arbitrary trapezoid \Delta, then we can initially compute all trapezoids for
which C these trapezoids, we initially set
For all the remaining trapezoids in \Xi, we use the recursive formula presented above to
compute their oe. Computing oe for a trapezoid requires computing the minimum of O(n 2 )
quantities. Thus the total running time of the algorithm is O(n 8 ). The following theorem
states the main result of our paper.
Theorem 6.2 Given a set P of n points in the plane and a boolean constraint C satisfying
condition (3.1), we can compute a geometric partition of P with respect to C using
the number of trapezoids in an optimal partition. Our
algorithm runs in worst-case time O(n 8 +n 6 Q(n)), where Q(n) is the time to decide whether
any subset R ' P .
Remark: By computing oe's in a more clever order and exploiting certain geometric properties
of a geometric partition, the time complexity of the above algorithm can be improved
by one order of magnitude. This minor improvement, however, doesn't seem worth the
effort needed to explain it.
Theorem 6.2 immediately implies polynomial-time approximation algorithms for the
surface approximation and the planar bichromatic partition problem. In the case of the
surface approximation problem, deciding C (\Delta) for a trapezoid \Delta requires checking whether
there is a plane - in R 3 such that the vertical distance between - and the points covered by
\Delta is at most ". This problem can be solved in linear time using the fixed-dimensional linear
programming algorithm of Megiddo [15]. more practical algorithm, running in time
O(n log n), is the following. Let A ' P denote the set of points covered by the trapezoid
\Delta. For a point respectively, denote the point p translated vertically
up and down by ". Let A
only if sets A + and A \Gamma can be separated by a plane. The two sets are separable if their
convex hulls are disjoint. This can be tested in O(n log n) time-for instance, see the book
by Preparata and Shamos [19].)
Theorem 6.3 Given a set S of n points in R 3 and a parameter " ? 0, in polynomial time
we can compute an "-approximate polyhedral terrain \Sigma with O(K
K o is the number of vertices in an optimal terrain. Our algorithm runs in O(n 8 ) worse-case
time.
In the planar bichromatic partition problem, deciding whether C checking
whether \Delta contains any point from the red set R. This can clearly be done in O(n)
time. Actually, with a preprocessing requiring O(n 2 log O(1) n) time, this test can be made in
O(log n) time for any trapezoid \Delta. Our main purpose, however, is to show the polynomial
time for the approximation algorithm.
Theorem 6.4 Given a set R of 'red' points and another set B of `blue' points in the plane,
we can find in polynomial time a set of O(K disjoint triangles that cover B but do
not contain any red point; K o is the number of triangles in an optimal solution.
Remark: In view of the remark following Theorem 6.1, given a set R of 'red' points and
another set B of 'blue' points in the plane, we can find in polynomial time a set of O(K
disjoint orthogonal rectangles that cover B but do not contain any red point. In this
case, the time complexity improves by a few orders of magnitude, because there are only
rectangle and each rectangle is subdivided into two rectangles by drawing
a horizontal or vertical line passing through one of the points. Omitting all the details, the
running time in this case is O(n 5 ).
Extensions
We can extend our algorithm to a slightly stronger form of surface approximation. In the
basic problem, the implicit function (surface) is represented by a set of sample points S.
What if the sample consists of two-dimensional compact, connected pieces? In this section,
we show an extension of our algorithm that deals with the case when the sample consists of
a set T of n horizontal triangles with pairwise disjoint xy-projection. (Since any polygon
can be decomposed into triangles, this case also handles polygons.) Our goal is to compute
a polyhedral terrain \Sigma, so that the vertical distance between any point in T
is at most ". We produce a terrain \Sigma with O(K is the
number of vertices in an optimal such surface.
be the input set of n horizontal triangles in R 3 with the property
that their vertical projections on the plane are pairwise disjoint. We will consistently
use the following notational convention: for an object s 2 R 3 , -
s denotes its orthogonal
projection on the plane z = 0, and for a subset g ' -
denotes the portion of s in R 3 such
that g. Abusing the notation slightly, we say that a set \Xi of trapezoids (or triangles)
in R 3 "-approximates T within a region Q ' R 2 if the vertical distance between T and \Xi
in Q is at most " and the vertical projections of trapezoids of \Xi are disjoint on z = 0.
Let S denote the set of vertices of the triangles in T , and let -
S be their orthogonal
projection on z = 0. We set
S, as the set of points in our abstract problem. The
constraint function is defined as follows. Given a trapezoid \Delta 2 R 3 , we have C
and only if the vertical distance between \Delta and any point in
directly above/below
\Delta is at most ". (Thus, while the point set P includes only the vertices of T , the constraint
set takes into consideration the whole triangles.) The constraint C satisfies (3.1), and it can
be computed in polynomial time.
It is also clear that the size of an optimal trapezoidal partition of P with respect to
C is a lower bound on the size of a similar partition for T , the set of triangles. We first
apply Theorem 6.3 to obtain a family \Delta of O(k log trapezoids that "-approximates P
with respect to C; clearly k - K . The next step of our algorithm is to extend \Delta to a
polyhedral terrain that "-approximates the triangles of T . Care must to be exercised in this
step if one wants to add only O(k log new trapezoids. In the second step, we work with
the projection of T and \Delta in the plane
(i) (ii) (iii)

Figure

5: (i) -
T and -
\Delta, (ii) R 1 and G; (iii) R 2 and Q i 's.
\Deltag and -
g. Let R be the set of connected components
of the closure of
\Delta. That is, R is the portion of
T lying in the
common exterior of -
\Delta, as shown in Figure 5 (i). R is a collection of simple polygons, each
of which is contained in a triangle of -
T . Since the corners of the triangles of -
are covered
by -
\Delta, the vertices of all polygons in R lie on the boundary of -
\Delta, and each edge of R is
contained in an edge of -
\Delta or of -
T . Let R 1 ' R be the subset of polygons that touch at
least three edges of trapezoids in -
\Delta, and let R
For each polygon P i 2 R 1 , we compute a set of triangles that "-approximate T within
. For a vertex lying on the boundary of a trapezoid -
\Delta, let -
v denote the point
on \Delta whose xy-projection is v. Let -
T i be the triangle containing P i . We triangulate P i ,
and, for each triangle 4abc in the triangulation, we pick 4-a - b-c. Since T i is parallel to the
xy-plane, it can be proved that the maximum vertical distance between 4-a - b-c and T i is ".
We repeat this step for all polygons in R 1 . The number of triangles generated in this step
is proportional to the number of vertices in the polygons of R 1 , which we bound in the
following lemma.
Lemma 7.1 The total number of vertices in the polygons of R 1 is O(k log k).
Proof: Each vertex of a polygon in R 1 is either (i) a vertex of a trapezoid in -
\Delta, or (ii) an
intersection point of an edge of -
\Delta with an edge of a triangle in T . There are only O(k log
vertices of type (i), so it remains to bound the number of vertices of type (ii).
We construct an undirected graph be the
set of edges in -
\Delta. 1 To avoid confusion, we will call the edges of \Gamma segments and those of E
arcs. For each segment fl i , we place a point i close to fl i , inside the trapezoid bounded by
. The set of resulting points forms the node set V . If there is an edge pq of a polygon in
R 1 such that we add the arc (i; j) to E; see Figure 5 (ii). It is easily seen
that G is a planar graph, and that log k). Fix a pair of segments
that (1; be the set of edges in R 1 , sorted either left to
right or top to bottom, as the case may be, that are incident to fl 1 and fl 2 . Let jE 12
Assume that for every 1 - lies on fl 1 and q i lies on fl 2 . The number of vertices
of type (ii) is obviously 2
We call two edges equivalent if the
interior of the convex hull of p i q i and p j q j does not intersect any trapezoid of -
\Delta. This
equivalence relation partitions E 12 into equivalent classes, each consisting of a contiguous
subsequence of E 12 . Let - ij denote the number of equivalence classes in E ij .
1: There are at most two edges in each equivalence class of E 12 .
Proof: Assume for the sake of a contradiction that three edges
belong to an equivalence class. Further, assume that the triangle -
T bounded by p i q i
lies below p i q i (see Figure 6 (i)). Since the quadrilateral Q defined by p i q i and q
does not contain any trapezoid of -
\Delta, p i+1 q i+1 is also an edge of -
T . But then Q
is a connected component of R and it touches only two edges of -
\Delta, thereby implying that
is not an edge of a polygon of R 1 , a contradiction. 2
Thus,
1 The segments of \Gamma may overlap, because the trapezoids of -
can touch each other. If a segment fl i of
\Gamma is an edge of two trapezoids, then no edge of R1 can be incident to fl i .

Figure

Edges in an equivalent class of E 12 , and (ii) edges in different equivalent classes
Next, we bound the quantity
12 be two consecutive
equivalent classes of E 12 , let p i q i be the bottom edge of E j
12 , and let p i+1 q i+1 be the top
edge of E j+1
12 . The quadrilateral Q
contains at least one trapezoid -
of -
\Delta.
We call the triangle edges of Q, and p i the trapezoidal edges of Q.
The triangle edges of Q
are adjacent in E 12 . Let
be
the set of resulting quadrilaterals. Since
suffices to bound the
number of quadrilaterals in Q.
Consider the planar subdivision induced by Q and call it A(Q). For each bounded face
Q(f) be the smallest quadrilateral of Q that contains f . Since the boundaries
of quadrilaterals do not cross, Q(f) is well defined.
2: Every face f of A(Q) can be uniquely associated with a trapezoid
\Delta such
that
Proof: The claim is obviously true for the unbounded face, so assume that f is a bounded
does not contain any other quadrilateral of Q, then so by
definition f contains a trapezoid of -
\Delta. If Q does contain another quadrilateral
of Q, let Q j be the largest trapezoid that lies inside Q i -that is, @Q j is a part of @f . If
none of the trapezoidal edges of Q j lies in the interior of Q i , then the trapezoidal edges of
both lie on the same segments of \Gamma, say, . Consequently, the triangle edges of
both belong to E 12 , which is impossible because then the triangle edges of Q i are
not adjacent in E 12 . Hence, one of the trapezoidal edges of Q j lies in the interior of Q i . Let
\Delta be the trapezoid bounded by this edge. Since the triangle edges of Q j lie outside -
the interior of -
\Delta does not intersect any edge of R 1 , -
lies in f . This completes the proof
of Claim 2. 2
By Claim 2, the number of faces in A(Q) is at most j -
log k). This completes
the proof of the lemma. 2
Next, we partition the polygons of R 2 into equivalence classes in the same way as we
Closing Remarks 21
partitioned the edges of E 12 in the proof of Lemma 7.1. That is, we call two polygons
endpoints lie on the same pair of edges in -
\Delta, and (ii) the
interior of the convex hull of does not intersect any trapezoid of -
\Delta. Using the same
argument as in the proof of the above lemma, the following lemma can be established.
Lemma 7.2 The edges of R 2 can be partitioned into O(k log equivalence classes.
For each equivalence class be the convex hull of E i -observe that Q i
is a convex quadrilateral, as illustrated in Figure 5 (iii). Each quadrilateral Q i can be "-
approximated using at most three triangles in R 3 in the same way as we approximated each
polygon P i of R 1 . By Lemma 7.2, the total number of triangles created in this step is also
O(k log k).
Putting together these pieces, we obtain the following lemma.
Lemma 7.3 The family of trapezoids \Delta can be supplemented with O(k log additional
trapezoids in R 3 so that all the triangles of T are "-approximated. The orthogonal projection
of all the trapezoids on the plane disjoint.
The area not covered by the projection of trapezoids found in the preceding lemma, of
course, can be approximated without any regards to the triangles of T . The final surface
has O(K trapezoids and it "-approximates the family of triangles T . We finish with
a statement of our main theorem in this section.
Theorem 7.4 Given a set of n horizontal triangles in R 3 , with pairwise disjoint projection
on the plane z = 0, and a parameter " ? 0, we can compute in polynomial time a "-
approximate polyhedral terrain of size O(K o log K is the size of an optimal
"-approximate terrain.
8 Closing Remarks
We presented an approximation technique for certain geometric covering problems with
a disjointness constraint. Our algorithm achieves a logarithmic performance guarantee
on the size of the cover, thus matching the bound achieved by the "greedy set cover"
heuristic for arbitrary sets and no disjointness constraint. Applications of our result include
polynomial time algorithms to approximate a monotone, polyhedral surface in 3-space, and
to approximate the disjoint cover by triangles of red-blue points. We also proved that these
problems are NP-Hard .
The surface approximation problem is an important problem in visualization and computer
graphics. The state of theoretical knowledge on this problem appears to be rather
slim. Except for the convex surfaces, no approximation algorithms with good performance
guarantees are known [6, 16]. For the approximation of convex polytopes, it turns out that
one does not need disjoint covering, and therefore the greedy set cover heuristic works.
We conclude by mentioning some open problems. An obvious open problem is to reduce
the running time of our algorithm for it to be of any practical value. Finding efficient
heuristics with good performance guarantees seems hard for most of the geometric partitioning
problems, and requires further work. A second problem of great practical interest
is to "-approximate general polyhedra-this problem arises in many real applications of
computer modeling. To the best of our knowledge, the latter problem remains open even
for the special case where one wants to find a minimum-vertex polyhedral surface that lies
between two monotone surfaces. The extension of our algorithm presented in Section 7 does
not work because we do not know how to handle the last fill-in stage.



--R

Fast greedy algorithms for geometric covering and other problems
An algorithm for piecewise linear approximation of an implicitly defined two-dimensional surfaces
An algorithm for piecewise linear approximation of an implicitly defined manifold
Decision trees for geometric models
Polygonization of implicit surfaces
Algorithms for polytope covering and approximation
Simplification of objects rendered by polygonal approximations
Several hardness results on problems of point separation and line stabbing
Optimal packing and covering in the plane are NP-complete
Piecewise linear approximations of digitized space curves with applications Scientific Visualization of Physical Phenomena pp.
Approximation algorithms for combinatorial problems
Computing 11
A high resolution 3D surface construction algo- rithm
On the ratio of optimal integral and fractional cover

Separation and approximation of polyhedral surfaces
Efficient binary space partitions for hidden-surface removal and solid modeling
Optimal binary space partitions for orthogonal objects
Computational Geometry: An Introduction
Decimation of triangle meshes
An adaptive subdivision method for surface fitting from sampled data
On some link distance problems in a simple polygon

--TR

--CTR
Gabriel Peyr , Stphane Mallat, Surface compression with geometric bandelets, ACM Transactions on Graphics (TOG), v.24 n.3, July 2005
David Cohen-Steiner , Pierre Alliez , Mathieu Desbrun, Variational shape approximation, ACM Transactions on Graphics (TOG), v.23 n.3, August 2004
Pankaj K. Agarwal , Boris Aronov , Vladlen Koltun, Efficient algorithms for bichromatic separability, ACM Transactions on Algorithms (TALG), v.2 n.2, p.209-227, April 2006
Pankaj K. Agarwal , Boris Aronov , Vladlen Koltun, Efficient algorithms for bichromatic separability, Proceedings of the fifteenth annual ACM-SIAM symposium on Discrete algorithms, January 11-14, 2004, New Orleans, Louisiana
Michiel Smid , Rahul Ray , Ulrich Wendt , Katharina Lange, Computing large planar regions in terrains, with an application to fracture surfaces, Discrete Applied Mathematics, v.139 n.1-3, p.253-264,
Mark de Berg , Micha Streppel, Approximate range searching using binary space partitions, Computational Geometry: Theory and Applications, v.33 n.3, p.139-151, February 2006
Bernard Chazelle , C. Seshadhri, Online geometric reconstruction, Proceedings of the twenty-second annual symposium on Computational geometry, June 05-07, 2006, Sedona, Arizona, USA
Joseph S. B. Mitchell , Joseph O'Rourke, Computational geometry, ACM SIGACT News, v.32 n.3, p.63-72, 09/01/2001
simplification;visualization;approximation algorithms;levels of detail;dynamic programming;machine learning;terrains