Polyline Fitting of Planar Points under Min-Sum Criteria
Boris Aronov∗

Tetsuo Asano†

Naoki Katoh‡

Kurt Mehlhorn§

Takeshi Tokuyama¶

Abstract
Fitting a curve of a certain type to a given set of points in the plane is a basic problem instatistics and has numerous applications. We consider fitting a polyline with k joints underthe min-sum criteria with respect to L1 - and L2 -metrics, which are more appropriate measuresthan uniform and Hausdorff metrics in statistical context. We present efficient algorithms forthe 1-joint versions of the problem and fully polynomial-time approximation schemes for thegeneral k-joint versions.

1

Introduction

Curve fitting aims to approximate a given set of points in the plane by a curve of a certain type.
This is a fundamental problem in statistics, and has numerous applications. In particular, it is abasic operation in regression analysis. Linear regression approximates a point set by a line, while
non-linear regression approximates it by a non-linear function from a given family.
In this paper, we consider the case where the points are fitted by a polygonal curve (polyline)with k joints, see Figure 1. This is often referred to as polygonal approximation or polygonal
fitting problem. It is used widely. For example, it is commonly employed in scientific andbusiness analysis to represent a data set by a polyline with a small number of joints. The best
representation is the polyline minimizing the error of approximation. Error is either defined as
the maximum (vertical) distance of any input point from the polyline (min-max-optimization) or
the sum of vertical distances (min-sum-approximation).
Min-max-approximation by a polyline is well studied. In one popular formulation one minimizes the maximum of the vertical distance (called the uniform metric or Chebyschev error
function) from the points to the curve. Hakimi and Schmeichel gave a O(n2 log n) time algorithm for this problem [11]; the time complexity was later improved to O(n2 ) [27] and then toO(n log n) [10]. Another popular approach is to minimize the Hausdorff measure that is the
maximum of the Euclidean distances between the points and the output curve. This problem
can also be solved in polynomial time [23]. These problems are closely related to curve simplification, in which the input is a polyline with n edges rather than a set of n points; this question
arises in geographic information systems (see the survey [28]) and has received much attention incomputational geometry [5, 12, 15, 22].
The minimize-the-maximum (min-max ) formulation is useful in pattern recognition applications. However, in applications to statistics, its serious deficiency is its extreme sensitivity to∗

Polytechnic University, Brooklyn, NY 11201, USA; http://cis.poly.edu/˜aronov. Supported in part by NSF
ITR Grant CCR-00-81964. Part of the work was carried out while B.A. was visiting JAIST, Universitat Politècnica
de Catalunya, and MPII.
†
Japan Advanced Institute of Science and Technology, Tatsunokuchi, 923-1292, Japan; t-asano@jaist.ac.jp. Partially supported by the Ministry of Education, Science, Sports and Culture, Grant-in-Aid for Scientific Research (B).
‡Kyoto University, Kyoto, 606-8501, Japan; naoki@archi.kyoto-u.ac.jp.
§
Max-Planck-Institut für Informatik, D-66123, Saarbrücken, Germany; mehlhorn@mpi-sb.mpg.de.
¶Tohoku University, Sendai, 980-8579, Japan; tokuyama@dais.is.tohoku.ac.jp.

1

v3
v4
v1

e1

e4
e3

e2
v2

e5

Figure 1: A 4-joint polyline fitting a set of points.
the presence and location of outliers. Even a single outlier can drastically change the output,
while outliers, real or imagined, are common in statistical data. For this reason, minimize-thesum-of-errors (min-sum) methods are considerably more popular in statistics: The most basic
one is the least-squares method that minimizes the sum of the squares of vertical distances between the input points and the output curve. In this paper, we call it the L2 -fitting problem(the term least-squares fitting is commonly used as well), and regression line in statistics usuallyrefers to the L2 -fitting line. If the output curve is either a straight line or a low-degree algebraic
curve, it is quite easy to compute the optimal L2 fitting. Another criterion is L1 -minimization,in which we minimize the sum of vertical distances from the candidate curve to the points being
fitted. L1 fitting is more resilient to outliers than L2 fitting; however, it is usually more expensive (or complicated) to compute the optimal solution. For linear regression, the L1 -optimalline can be computed in linear time [13], but it requires sophisticated computational techniques.
Several other formulations have been proposed for further reducing the effect of outliers on thelinear regression. Repeated median regression is a well-known example, and efficient solutions are
known for several other criteria [16]. Takeshi says: Inserted a reference suggested by a referee. ←−Boris says: ok, slightly reworded Linear regression considering Euclidean or Manhattan distance ←−between points and the line has also been studied [21].
In this paper, we focus on the L1 - and L2 -fitting problems when the desired curve is a kjoint polyline; in other words, it is a continuous piecewise-linear x-monotone curve with k + 1
linear components. We assume that a coordinate system is fixed, and the input points are
sorted with respect to their x-coordinate values. To the authors’ knowledge, the computational
complexity of the optimal k-joint problem under either of these minimization criteria has notbeen previously investigated. More specifically, it seems that an efficient solution of the L1 -fitting
problem extending the result of Imai et al. [13] is theoretically challenging even for the 1-joint
problem.
In this paper, we begin by considering the 1-joint problem. We give algorithms of complexity
O(n) and Õ(n4/3 ) time for the L2 and L1 criteria, respectively.1 The L2 -fitting algorithm is simple
and practical, whereas the L1 -fitting algorithm depends on using a semi-dynamic range searchdata structure and parametric search. For general k, we present two approximation schemes. Let
zopt be the minimum fitting error for a k-joint polyline and let ε be a positive constant. We givea polynomial-time approximation scheme (PTAS) to compute a b(1 + ε)kc-joint fitting whose
error is at most zopt and we describe a fully polynomial-time approximation scheme (FPTAS) to
compute a k-joint polyline with (1 + ε)zopt fitting error, and consequently show that the problems
cannot be strongly NP-hard, although their NP-hardness remains open. Boris says: I asked
before: “Can we do 1+eps approximation of BOTH?” We do not have it in the paper now, but
1

We write f (n) = Õ(g(n)) if there exists an absolute constant c ≥ 0 such that f (n) = O(g(n) logc n).

2perhaps we should?
←−
Intuitively, why are the problems we consider in this paper more difficult than some related
questions? We have mentioned that the uniform metric fitting problem can be solved efficiently.
The key point is that its corresponding decision problem to determine whether there exists a kjoint polyline with a uniform error less than a given value w is geometrically a stabbing problem.
The k-joint path must go through n vertical line segments of length 2w centered at the input
points, and one can continuously move a feasible polyline so that each link becomes extremal ingeometric sense, that is, goes through a pair of endpoints of vertical segments. Hence, we candesign a polynomial-time algorithm to find the optimal path based on dynamic programming.
The L1 - and L2 -fitting problems seem to be more subtle: We do not have reformulations of thecorresponding decision problems in terms of stabbing.

2

Preliminaries

A k-joint polyline is an alternating sequence P = (e1 , v1 , e2 , v2 , . . . , ek , vk , ek+1 ) of line segments(links) and joint vertices (joints), where es and es+1 share the endpoint vs , for s = 1, 2, . . . , k,and e1 and ek+1 are infinite rays. We denote the link es on line y = as x − bs by (as , bs ) if theinterval of the values of x corresponding to the link is understood. A joint vs is represented by
the pair (us , vs ) of its coordinate values. Thus, the connectivity and monotonicity of the polyline
can be guaranteed by requiring that vs = as us − bs = as+1 us − bs+1 , for s = 1, 2, . . . , k + 1, andu1 < . . . < uk . Boris says: Technically, this is wrong, as the last equality only works up to s = k
not s = k + 1.
←−
We now formulate the problem of fitting a k-joint polyline to an n-point set. Given a setof points S = {p1 = (x1 , y1 ), p2 = (x2 , y2 ), . . . , pn = (xn , yn )} with x1 < x2 < . . . < xn
and an integer k, and setting u0 = −∞ and uk+1 = ∞ for convenience, find a polyline P =((a1 , b1 ), (u1 , v1 ), (a2 , b2 ), (u2 , v2 ), . . . , (uk , vk ), (ak+1 , bk+1 )) minimizing one of the following threequantities for L1 -, L2 -, and uniform metric fitting, respectively:
L1 :

k+1
X

X

|as xi − bs − yi |,(1)(as xi − bs − yi )2 ,(2)

s=1 us−1 <xi ≤us

L2 :

k+1
X

X

s=1 us−1 <xi ≤us

Uniform metric :

max



max

s=1,...,k+1 us−1 ≤xi ≤us|as xi − bs − yi | .

(3)

For k = 0, the problems are linear regression problems. The L2 -linear
Pn regressionPisn wellknownP
as the Gaussian
least-squares
method.
Once
wecomputeA=ni=1 xi , Bn =i=1 yi ,
Pn
Pn
n
2
2
Cn = i=1 xi , Dn = i=1 xi , and En = i=1 xi yi in linear time, we can construct an optimal
fitting line y = ax − b by considering the partial derivatives of the objective function and solving
a 2 × 2 system of linear equations. The linear regression problem with respect to the uniform
error is equivalent to finding a pair of parallel lines at the minimum vertical distance that containall the given points between them. This can be done by applying the rotating caliper method thatcomputes antipodal pairs of points on the convex hull of the point set. For an x-sorted point set
this can be done in O(n) time [24]. The L1 -linear regression problem is more involved; however,a linear-time algorithm has been devised by Imai et al. [13] based on Megiddo’s prune-and-searchparadigm.

3

3

Fitting a 1-joint polyline

We consider the problem of fitting a 1-joint polyline to a set of points. We proceed in two steps.
We first assume that the joint vertex lies in a fixed interval [xq , xq+1 ] and later eliminate thisassumption. Let S1 (q) = {p1 , p2 , . . . , pq } and S2 (q) = {pq+1 , . . . , pn }. Our objective polyline
consists of two links lying on lines `1 : y = a1 x − b1 and `2 : y = a2 x − b2 , respectively. We call atuple (a1 , b1 , a2 , b2 ) feasible if the two lines y = a1 x − b1 and y = a2 x − b2 meet at a point whose2x-coordinate u = ab11 −b−a2 lies in the interval [xq , xq+1 ]. Our goal here is to find a feasible tuple(a1 , b1 , a2 , b2 ) representing a 1-joint polyline minimizing
q
X
i=1
q
X
i=1

n
X

|a1 xi − b1 − yi | +|a2 xi − b2 − yi | and(4)(a2 xi − b2 − yi )2 ,

(5)

i=q+1
n
X

(a1 xi − b1 − yi )2 +i=q+1

for L1 - and L2 -fitting, respectively. Minimizing (4) is equivalent to, provided a1 6= a2 , minimizingPn
i=1 wi subject to−wi ≤ a1 xi − b1 − yi ≤ wi , for i ≤ q,
−wi ≤ a2 xi − b2 − yi ≤ wi , for i ≥ q + 1,
b1 − b2≤ xq+1 ,xq ≤a1 − a2and(6)where the last line represents the feasibility condition.
Lemma 3.1. For either L1 - or L2 -fitting criterion, the 1-joint problem for a fixed q reduces tosolving two convex programming problems.
Proof. Disregarding the feasibility constraint, the problem is clearly a quadratic programming
problem for the L2 case and a linear programming problem for the L1 case. The feasibility
constraint requiring that the lines `1 and `2 meet in the strip between xq and xq+1 can be
expressed by different linear constraints depending on whether a1 ≤ a2 . Thus, we can decompose
the (L1 or L2 ) problem into two subproblems. If a1 ≤ a2 , the lines meet in the strip if and only
if `1 is not below `2 at xq and is not above it at xq+1 . Thus, the additional constraint becomesxq (a2 − a1 ) ≤ b2 − b1 ≤ xq+1 (a2 − a1 ) .

(7)

In the opposite case, the additional constraint is
xq+1 (a2 − a1 ) ≤ b2 − b1 ≤ xq (a2 − a1 ) .

(8)

Clearly, each subproblem is a convex programming problem, as claimed.
From the above lemma, it is clear that the optimal 1-joint polyline can be computed by using
linear/quadratic programming. However, we aim to design combinatorial algorithms for these
problems. Indeed, we can classify the solution into two types: (a) An inequality in (7), (8) holdswith equality. (b) All of the inequalities in (7), (8) are strict. We call the solution fixed in
the former case and free otherwise. From the form of the expressions in (7), (8) we deduce thefollowing simple observation.
Lemma 3.2. If the solution is fixed, the joint is located on either of the two vertical lines x = xq ,x = xq+1 .
4

If the joint is on the line x = xq+1 , we can regard it as a solution for the partition into
S1 (q + 1) = S1 (q) ∪ {pq+1 } and S2 (q + 1) = S2 (q) \ {pq+1 }. Thus, for each partition, we essentiallyneed to solve two subproblems: (1) the free problem and (2) the fixed problem where the joint ison the vertical line x = xq . This leads to the following generic algorithm: For each partition of S
into two intervals S1 and S2 , we first consider the free problem ignoring the feasibility constraint,
and check whether the resulting solution is feasible, i.e., we verify that the intersection pointlies in the strip between pq and pq+1 . If it is feasible, it is the best solution for the partition.
Otherwise, we consider the fixed solution adding the constraint that the joint lie on x = xq , andreport the solution for the partition. After processing all n − 1 possible partitions, we report thesolution with the smallest error.
If it takes O(f (n)) time to process a subproblem for each partition, the total time complexity
is O(nf (n)). For efficiency, we design a dynamic algorithm to process each partition so that f (n)is reduced in the amortized sense.

3.1

The L2 1-joint problem

We show how to construct an optimal L2 -fitting 1-joint polyline in linear time. We process
the partitions (S1 (q), S2 (q)) starting from= n − 1, in order.
We maintain
Pq thePq q = 1 to q PPqq2
2, D =
sums, variances,andcovariances
A
=
x
,
B
=
y
,
C
=
x
q
q
q
q
i=1 yi ,
i=1 i
i=1 i
i=1 i
Pq
and Eq = i=1 xi yi incrementally, at constant amortized cost. They also
Pprovide us with thecorresponding values for S2 (q) if we precompute those values for S, i.e., ni=q+1 xi = An − Aqetc.
For the free case, the objective function is separable,in the sense that the optimal solution
Pq
2
can
be
identified
byfinding(a,b)
minimizing(ax
1 1
i=1 1 i − b1 − yi ) and (a2 , b2 ) minimizingPn2
j=q+1 (a2 xj − b2 − yj ) independently. Each can be computed in O(1) time from the values of
Aq , . . . , Eq as explained in section 2. The feasibility check of the solution is done in O(1) timeby computing the intersection point of the corresponding pair of lines. It remains to solve thesubproblems with the additional constraint that the joint is at x = xq . Put
f (a1 , b1 , a2 , b2 ) =q
X(a1 xi − b1 − yi )2 +

i=1

n
X

(a2 xj − b2 − yj )2 ,

(9)

j=q+1

g(a1 , b1 , a2 , b2 ) = a1 xq − b1 − a2 xq + b2 , and(10)

L(a1 , b1 , a2 , b2 ) = f (a1 , b1 , a2 , b2 ) − λg(a1 , b1 , a2 , b2 ),(11)

so that f (·) is the function to be minimized and the joint constraint can be expressed as g(·) = 0.
Then, by the Kuhn-Tucker condition the optimal solution Zopt = (a01 , b01 , a02 , b02 ) describing a best
L2 -fitting 1-joint polyline for a fixed value of q has to satisfy
∂L
∂L
∂L
∂L
=
=
=
= 0,
∂a1 Zopt∂b1 Zopt∂a2 Zopt∂b2 Zopt(12)g(Zopt ) = 0.

(13)

and
This gives us a set of five linear equations that must be satisfied by the optimal parameter valuesof a1 , b1 , a2 , b2 and the Lagrange multiplier λ. The coefficients can be expressed in terms of
xq , Aq , . . . , Eq , and this system can be solved in constant time for each q. Thus, we have thefollowing:
Theorem 3.3. L2 -optimal 1-joint fitting can be computed in linear time.

5

3.2
3.2.1

The L1 1-joint problem
Semi-dynamic L1 linear regression

We start with the problem of computing the optimal linear L1 -fittingP(i.e., linear regression) of
the input point set, i.e., we seek the line `opt : y = ax − b minimizing ni=1 |axi − b − yi |.
Boris says: In the next two sentences the word “problem” appears 4 times. A bit too much. ←−The difficulty with the L1 -fitting problem is that, written in linear programming terms (as in (6)),
it has n + 2 variables, in contrast to the least-squares case where the problem is directly solved
as a bivariate problem. Nonetheless, the problem can be solved by a brute-force combinatorial
algorithm in O(n3 ) time, since there are O(n2 ) possible linear dissections of the point set which
can be enumerated in Θ(n2 ) worst-case time by constructing the dual arrangement, and one cancompute the optimal line in linear time once the dissection by the line is given (this algorithm
can be easily sped up to constant or near-constant amortized time per dissection). Moreover,by Lemma 3.4, the optimal line bisects S into two equal-size subsets; in other words, the line
is a halving line. Using this fact, Imai et al. [13] devised an optimal linear-time algorithm forcomputing `opt based on the multidimensional prune-and-search paradigm. In order to design anefficient algorithm for the 1-joint fitting problem, we consider a semi-dynamic version of the L1
linear regression for a point set P with low amortized time complexity, where we dynamicallymaintain P with insertions and deletions under an assumption that P is always a subset of a
fixed universe S of size n that is given from the outset. (In fact, for our application, it is sufficient
to be able to start with P = ∅ and handle only insertions, and to start with P = S and handle
only deletions. Moreover, the order of insertions and deletions is known in advance. The data
structure we describe below is more general.)Consider the dual space, with pi = (xi , yi ) transformed to the dual line Y = fi (X) wherefi (X) = xi X − yi . The line y = ax − b is transformed to the point (a, b) in the dual space. The
kth level of the arrangement A = A(S ∗ ) of the set S ∗ of dual lines is the trajectory of the kth
largest value among fi (X).2 We call the dn/2eth level the median level.
Lemma 3.4 (Imai et al. [13]). If the optimal L1 -fitting line is given by y = aopt x − bopt , its dual
point (aopt , bopt ) is on the median level if n is odd, and between the n2 th level and the ( n2 + 1)thlevel if Boris says: Shouldn’t “( n2 + 1)th” be “( n2 + 1)st”? n is even.
←−In case the optimal slope is not unique, let aopt denote the smallest such slope. Now, given
X-value t, consider the point (t, fi (t)) for each i = 1, 2, . . . , n, and let F (t) be the sum of thebn/2c largest values in {fi (t) : i = 1, 2, . . . , n} and G(t) be the sum of the bn/2c smallest values inthe same set. Put H(t) = F (t) − G(t). H(t) gives the L1 fitting error of the dual line of any point(t, y) on the median level (or between the two median levels if n is even). Thus, by Lemma 3.4,
H(t) is minimized at t = aopt .
Lemma 3.5. F (t) is a convex function, while G(t) is concave. As a consequence, H(t) is also
convex. At t = aopt the slope of H(t) changes from negative to non-negative.
Proof. The convexity follows directly from the fact that, in any line arrangement, the portionsof the lines lying on or below (resp. on or above) any fixed level k can be decomposed into k
non-overlapping concave (resp. convex) chains; see, for example, [3].
Suppose a fixed universe S ∗ of lines is given. We need a data structure that maintains a
subset P ∗ ⊆ S ∗ and supports the following operations on P ∗ :Median-location query For a query value t, return the point on the bn/2cth highest line at
X = t.
2

We use an asterisk to denote geometric dual of a point, line, or a set of lines/points.

6

Slope-sum query For a query point p = (t, y), return the sum of the slopes of lines below p at
X = t.
Height-sum query For a query value p = (t, y), return the sum of the Y -coordinates of the lines
below p at X = t. The height-sum query is reduced to a slope-sum query plus a constantterm-sum query that reports the sum of the constant terms of the equations representing
lines.
Update A line in S ∗ is added to or removed from P ∗ .
Suppose a data structure supporting such queries on a set P ∗ ⊆ S ∗ of lines in O(τ (n)) time
is available, where n = |S ∗ |. Then we can query the slopes of F and G at t, and hence computethe slope of H at t in O(τ (n)) time.
Takeshi says: I add some explanation for the case where we encounter a break point. We indeed
need not consider slope if we apply perturbation, but I do not change in that way. Boris says: ←−Slightly rephrased←−
It may happen that the piecewise linear function H has a break point at t. In order to handle
this case, we apply symbolic perturbation and examine the slopes of H at t +  and t −  for an
infinitesimally small  > 0 to determine whether H increases to the right of t, to the left of t, orin both directions. This perturbation is only activated when a linear decision required in queries
needs a tie-break. Because of convexity of H, we have the following:
Lemma 3.6. Given t, we can decide whether t < aopt , t > aopt , or t = aopt in O(τ (n)) time.
Thus, we can perform binary search to find aopt . We show below how to make this searchstrongly polynomial. Once we know aopt , we determine bopt by the median-location query att = aopt .
3.2.2

Semi-dynamic data structure for the queries

We show how to realize semi-dynamic median-location- and sum-queries. As a preliminary step,
we describe a semi-dynamic data structure for vertical ray queries, i.e., queries of the form: Givena vertical upward ray starting at (t, z) determine the number of lines in P ∗ intersected by theray, the sum of their slopes, and the sum of their constant terms. A dual line Y = xi X − yiis above (t, z) iff the primal point (xi , yi ) is above the line y = tx − z. Thus our queries are
reduced to half-space queries in the primal plane. We use the partition-tree data structure of√
Matoušek [4, 17, 19]. It supports half-space queries on sets with n points in time O( n), linear
space, and preprocessing time O(n log n).
We build a partition tree T (S) on the set S of points dual to the lines in S ∗ (in fact, these
are the points to which a line is being fitted). A standard construction proceeds as follows: With
each node v of the partition tree we associate a point set S(v) ⊆ S and a triangle ∆(v) ⊃ S(v),where S(v) ⊂ S(parent(v)) at any node v other than the root and S(v) = S P
at the root. In
addition wePalso store at v the size |S(v)| of S(v) and the sums ξ(S(v)) =pi ∈S(v) xi andχ(S(v)) = pi ∈S(v) yi of the slopes and constant terms of the corresponding dual lines. Since
the point sets S(v), over all children v of a node w in the tree, by definition of a partition tree,
partition the set S(w), and |S(v)| is at most a fraction of |S(w)|, this tree has linear size andlogarithmic depth. For our purposes, we modify the partition tree to obtain a new tree T (S, P )where the same ∆(v) as in T (S) is associated with every node v, but v stores P (v) = S(v) ∩ P ,ξ(P (v)) and χ(P (v)) instead of the corresponding values for S(v). This data structure enables
us to execute the half-plane range query in P , and thus the vertical ray query in P ∗ .
Our data structure is semi-dynamic. When P changes, with a point p being added or removed,
what we need to update is just values |P (v)|, ξ(P (v)), and χ(P (v)) for each node v where p is7relevant. Since the sets S(v) for all nodes v at a fixed level of the partition tree form a partition ofS, only one node must be updated at each level; to facilitate the update one might associate witheach point p ∈ S a list of length O(log n) containing the nodes v of the tree with p ∈ S(v). Boris
says: Doesn’t this decision increase the space to Θ(n log n). Thus, the update can be performed in ←−O(log n) time. This ends the description of the semi-dynamic vertical ray query data structure.
Our sum-queries can be done by using the vertical ray query.
We next turn to the median-location query data structure. For a given t, let m(t) = (t, y(t))
be the intersection of the vertical line X = t and the median level of the dual arrangement A(S ∗ ).
We can use the vertical query data structure to compare any given η with y(t). We perform avertical ray query to find the number of lines above (t, η). If it is less than bn/2c, y(t) < η;otherwise y(t) ≥ η. This suggests computing y(t) by some kind of binary search. If we had thesorted list of intersections between the vertical line X = t and the lines in S ∗ available, we could
perform a binary search on L by using O(log n) ray queries. However, it takes O(n log n) time
to compute the list, which is too expensive since we aim for a sublinear query time. Instead, weconstruct a data structure which can simulate the binary search without explicitly computing thesorted list.
Lemma 3.7. We can construct a randomized data structure in time O(n log n) such that, given√t, we can compute y(t) in O( n log n) time. The query time bound holds for every vertical line
X = t with high probability.
Proof. Takeshi says: I add some remarks on dynamization. Boris says: Rephrased. Recall that
we start with an underlying fixed set S of points and aim to build a semi-dynamic data structure
for answering queries on a subset P ∗ of the set S ∗ of dual lines.
We fix a small constant ε > 0, and randomly select cn1−ε lines from Ψ0 = S ∗ , to have aset Ψ1 of lines, where c is a suitable constant. From the results of Clarkson and Shor [9], if theconstant c is sufficiently large, with high probability every vertical segment intersecting no lineof Ψ1 intersects at most nε log n lines of S ∗ . In other words, Ψ1 is the dual of an (nε−1 log n)-net
εn, 1}-net of Ψ∗i . Thus,of S. Similarly, we construct Ψi+1 from Ψi such that Ψ∗i+1 is a min{ n |Ψlogi|εwe have a filtration Ψ0 ⊃ Ψ1 ⊃ . . . ⊃ Ψk , and |Ψk | ≤ n . The number k of layers is a function of
ε and c only, so the construction takes O(n) time.
Additionally, we construct a dual range-searching data structure for Ψi such that for a query√
vertical interval I we can report all lines in Ψi meeting I in O( n + K) time, where K is thenumber of reported lines. In primal space a vertical interval corresponds to a strip bounded by
two parallel lines and hence we may use partition trees as described above to implement reporting
queries. The preprocessing time is O(n log n).
Now, our algorithm for finding y(t) is as follows: Given t, we first compute all the intersections
between X = t and the lines of Ψk , sort them to have a list y (1) , y (2) , . . . y (m) of y-values wherem ≤ n . By using the semi-dynamic vertical ray-query data structure for P ∗ , we can decide
whether y(t) > y (i) or Boris says: I changed y(t) ≤ y (i+1) to y(t) ≤ y i , is this ok? y(t) ≤ y i√in O( n) time. Accordingly, we can perform binary search to find the unique interval Ik =√(y (i) , y (i+1) ] containing y(t) in O( n log n) time.
The interior of the vertical interval Ik containing y(t) is crossed by no line of Ψk . By using√
the dual range-searching data structure, we extract, in time O( n + K), the set of K lines inΨk−1 intersecting Ik ; K = O(nε log n) with high probability. Boris says: Added... Starting with
Ik and the extracted set of K lines, we compute Ik−1 ⊆ Ik not intersected by any line of Ψk−1 .
Proceeding recursively, we eventually obtain y(t), since at the last level of the filtration we arriveat an interval I0 containing y(t), with no line of Ψ0 = S ∗ crossing its interior. The total time is√√O(nε log2 n + n log n) = O( n log n).

8

←−
←−

←−

←−

√
We have described an O( n log n) time realization of the semi-dynamic query data structure,√i.e., τ (n) = O( n log n) in this data structure.
3.2.3

Determining aopt

We finally come to the strongly polynomial method for determining aopt via parametric search [25].
We remark that we do not employ parametric search to compute y(t) for a fixed t, since it is notalways possible to use parametric search in a nested fashion, and there are technical difficulties
in applying multi-dimensional parametric search paradigm [20, 26] to our problem.
Parametric search identifies a real number aopt . The search has two ingredients:• A decision procedure D(t) of one real parameter. The procedure tests whether t < aopt intime O(TD ).
• A master program M (t) which takes as input a real parameter t ∈ (−∞, ∞). M (t) is aparallel program which takes time TM with p processors for any fixed value of t. We assume
that M (t) proceeds in parallel rounds, as follows: Each processor performs some calculationand then makes a decision that depends on whether or not the input value t is larger thana just computed threshold value which does not depend on t (more generally, there may
be a constant number of such threshold values per processor per round and the decision
depends on where t lies among them). Once the decisions are made by all processors, thenext parallel round begins. Thus, in each round of the parallel program, we have at most p
threshold values that subdivide (−∞, ∞) into at most p + 1 subintervals, and all decisions
in the round are determined by the subinterval containing t. We assume subintervals to beopen on the left and closed on the right and we require that aopt is generated as a thresholdvalue when M (t) is run with t = aopt . In this case, the final interval computed by theprocedure described below has aopt as its right boundary.
Parametric search simulates the execution of M (t) for t = aopt . We maintain an interval I(open on the left and closed on the right) containing aopt and a set Q of real numbers. Initially,
I = (−∞, ∞) and Q = ∅. We simulate M (t), one round at a time. In each round, we generate thethreshold value of each processor and, if it is contained in I, insert it into Q. Then, we computethe median m of Q and compare m with aopt by invoking D(m). If m ≥ aopt , we replace I by
I ∩ (−∞, m], otherwise by I ∩ (m, ∞). The elements of Q not contained in the updated I are
removed from Q. We iterate this process until I has no element of Q in its interior. Then, we candetermine all decisions in the current round and proceed to the next round. Thus, we complete
the simulation in O(pTM + TM TD log p) time.
The decision procedure: In our case aopt is the slope of the optimal L1 -line. The decision
procedure proceeds as follows. Let t be the input. Boris says: Slight notational confusion: Earlier
we used Ik , Ik−1 to denote the vertical interval we keep narrowing down (and its new shorter version).
Here we call them s, s0 and use I to denote the interval of values for aopt . What do we do? If it
causes not too much trouble, I suggest we leave I alone (no subscripts!) and replace s by Ik and s0
by Ik−1 or would this be too hard to follow?
←−
1. Determine the y-coordinate y(t) of the median level of the dual arrangement at t using
the filtration (Ψi ) of the dual arrangement. At level i of the filtration, we have a verticalsegment s delimited by two lines of Ψi and not properly intersected by any line of Ψi .
Boris says: Replaced Ψi+1 by Ψi−1 throughout, as they are numbered from the largest to
the smallest.
←−9(a) First determine the lines of Ψi−1 intersecting s. This is a range query, where the range
is the strip defined by the two primal parallel lines dual to the endpoints of s. The√
running time is O( n + K) where K = Õ(nε ) is the number of lines intersected. Whatkinds of comparisons are we making here? We have a line tx + b, for given b, and wecheck whether this line intersects a triangle. This can be rewritten as a comparison oft with a real number.
(b) Once we have determined the intersected lines, we sort the intersections and performbinary search on them. For each point, we need to count the number of lines above it.
This is again a range query. So with log n “number of lines above a point” queries, wehave reduced s to s0 delimited by two adjacent lines of Ψi−1 .
The comparison of two lines at the vertical line t, or of at + b with a0 t + b0 , is equivalent
to comparing t with a real number.
2. Once we have found y(t), we use a single slope query and a single constant term query to
determine the slope of H(t) at t. The sign of this slope determines whether t < aopt .
The Master Program: The master program is simply the parallel version of the decision
procedure. Takeshi says: I gave minor change. Boris says: slightly reworded It uses O(τ (n)) ←−processors, where τ (n) is the sequential query time. We first present the version with τ (n) = ←−√O( n log n), and refine our description below. The range queries are easily parallelized because
they amount to walking down at most log n paths in the partition tree. Thus the parallel timeof a range query is O(log n). Step 1(b) amounts to a parallel sort (time O(log n)) followed bya binary search (time O(log n)). Each step of the binary search requires a range query. Step 2
is also a single range query. Thus the parallel time is O(log2 n) with O(τ (n)) processors. When
executed with t = aopt , the master program generates aopt as a threshold value since two duallines intersect at aopt .
Some Intuition: It is instructive to see what the parametric search does geometrically. In
each Ψi it determines a trapezoid T with two vertical walls (sides) and two non-vertical edges(top and bottom). The non-vertical edges lie on adjacent lines in Ψi . The median level of thefull arrangement intersects the sides and avoids the top and bottom. In addition, aopt lies in thex-span I of T . The interval I is the interval maintained in the parametric search.
In step 1(a), we narrow T by moving in its vertical walls, which ensures that no lines of Ψi−1
enters the reduced trapezoid through its top and bottom. In other words, all vertical segments
connecting the top to the bottom of the trapezoid intersect the same set of lines of Ψi−1 .
Next we come to sorting the intersections (step 1(b)) with the generic vertical segment spanned
by T . Step 1(b) will further move in the vertical walls until the trapezoid contains no intersections
between lines in Ψi−1 .
At this point the lines in Ψi−1 are sorted and we can perform binary search on them. Eachsearch step will move one of the non-vertical walls. Also in each search step we do a range searchand this may further move in the vertical walls (since the median level must stay within theshrinking trapezoid).
This completes our description of the parametric-search-based procedure for determining aopt .
3.2.4

Time-space tradeoff

To speed up the query time τ (n) and thus the overall algorithm, we generalize the data structure
to allow it to use super-linear storage based on Matoušek’s construction [18]. If we can use O(m)space for n < m < n2 , we first select r = O(m/n) points from S and construct a dual cutting,i.e., a decomposition of the dual plane into cells, such that each cell C is intersected by at most10

n/r lines dual to points of S; the number of cells required is O(r2 ) and the computation time isO(nr).
Let S(C) be the set of lines intersecting C. We construct a point-location data structure on
the cutting. For each cell C, we store the cumulative statistics (the sum of slopes etc.) for theset of lines
ppassing below C, and construct the partition tree for S(C). The query time of eachtree is Õ( n/r). When P changes, we need to update the data stored in each of the O(r2 ) cellsof the cutting, and also the O(r) partition trees corresponding to sets containing the updated
point. Thus, update time is O(r2 + r log n).
Update time can be sped up by not storing the statistics for each cell explicitly, but rather
retrieving them when needed at a cost of O(log r). This reduces the time needed for an updateto O(r log r) as shown below. Thus, if we set r = n1/3 , the update time and query time isτ (n) = Õ(n1/3 ). The space and preprocessing time is Õ(n4/3 ). The parallel time complexity is
not affected by the space-time trade-off.
Now we explain how to reduce the update time to O(r log r). The issue is that each of O(r2 )cells in the cutting has an attached data structure (partition tree) and some additional constantsize cumulative statistics (total number of lines passing below/above the cell and sums of theirslopes and constant terms). Let ` be the line inserted into or deleted from P in the current
update. The partition trees are easy to update when a line is inserted/removed and there are
only O(r) of them to modify; indeed, only the partition trees stored at cells intersected by ` areaffected, and there are O(r) such cells because of the zone theorem [7] for arrangements.
On the other hand, the line ` contributes to the global statistics of O(r2 ) cells; thus we needΘ(r2 ) time to update them directly. Instead, we store the information implicitly, so that it canbe computed as needed in time O(log r).
Consider a possibly non-simple spanning path σ in the cell adjacency graph of the cutting
obtained, for example, by tracing along a spanning tree of the graph. The length of the path isO(r2 ). Number the cells in the order of their occurrence along σ; note that a cell might appear on
the path more than once—only the first occurrence contains real data and subsequent occurrence
are treated as dummy.
Build a balanced binary tree B on top of σ. Each internal node v corresponds to a collectionof consecutive nodes on σ. Let the above count of v be the number of lines in P which both (1)pass above all the cells corresponding to the collection, and (2) do not pass above all the cells
corresponding to the collection for the parent node of v. We call this number the above count ofv. Symmetrically, we define the below count. The tree B is static, but each node v stores with it
the above and below counts that are maintained dynamically.
Initialize all the counts to 0 when P = ∅. Now consider adding a line ` to P (deletion ishandled symmetrically). Line ` intersects O(r) cells of the cutting. We can afford to computewhich cells these are, explicitly (and we have to do it, to update the partition tree information).
This gives m = O(r) cell numbers. Sort them to obtain c0 = −∞, c1 , c2 , ..., cm , cm+1 = +∞.
Each interval (ci , ci+1 ), if non-empty (i.e. if ci+1 > ci + 1), corresponds to a connected subpathof σ lying completely to one side of `. We can test in constant time on which side it lies bychecking one of its cells. Now decompose this interval into a logarithmic number of canonicalones (corresponding to nodes of B) and increment above/below counts at the O(log r) nodes.
Repeat for each of the O(r) subpaths. The time for an update of global counts is thus reducedto O(r log r), as claimed.
Takeshi says: Added some details on paralle algorithm. Boris says: I am confused: Why do ←−we need any of this at all? Don’t we need SOME master algorithm that makes a decision at the value
we are looking for, and it need not be related to the decision algorithm? Why can’t we just use the oldmaster algorithm? Or is the number of processors too large? Perhaps we should explain something? ←−Since we apply parametric searching later, we need to have a parallel algorithm for the range
query with the time-space tradeoff. Note that we do not need parallel algorithms for construction
11and update of the data structure. For simplicity, we consider the counting halfplane range query.
We first find in the arrangement of r lines the cell C containing the dual point of the line defining
the halfplane. This is point location, and can be done in O(log n) time by a single processor. The
retrieval of the cumulative statistics of C can be also done in O(log n) time by apsingle processor.
Now, we do range searching in S(C) in time O(log n) by using O(τ (n)) = Õ( n/r) processors
as discussed in Section 3.2.3. Thus, the parallel range query is done in O(log n) time by using
O(τ (n)) processors. Naturally, median-location, slope-sum, and height-sum queries can be donein polylogarithmic time (indeed, O(log2 n) time) by using O(τ (n)) processors.
3.2.5

Algorithm for L1 1-joint fitting

Takeshi says: This subsubsecion has been considerably changed.
←−
Finally, we describe the algorithm to find the L1 -optimal 1-joint polyline fitting a set S of npoints in the plane. Recall that there are two different types of solutions:Type 1 There is an index q such that the 1-joint polyline consists of the optimal L1 -fitting lineof S1 (q) = {p1 , p2 , . . . , pq } and that of S2 (q) = {pq+1 , pq+2 , . . . , pn }.
Type 2 There is an index q such that the joint lies on the vertical line x = xq .
If the optimal solution is of type 1, we compute an optimal L1 -fitting line for S1 (q) and S2 (q)separately, for every q = 1, 2, . . . , n, by using the semi-dynamic algorithm with S as the universe.
Takeshi says: Minor change is given. Boris says: ok Thus, the time complexity is Õ(nτ (n)) ←−for examining all q = 1, 2, . . . , n. In particular, if we use quasi-linear space Õ(n), the time ←−complexity is Õ(n1.5 ), and if we use O(n4/3 ) space, the time complexity is Õ(n4/3 ).
Takeshi says: I gave major revision below. Otherwise, the optimal solution is of type 2. If ←−the index q is given, our objective function is0

f (a, a , z) =qX|a(xi − xq ) − yi + z| +

i=1

n
X

|a0 (xi − xq ) − yi + z)|.

i=q+1

The variable z corresponds to the y-coordinate value of the joint, and a and a0 are slopes of two
halflines to the left and right of the joint, respectively. Boris says: Rephrased. As a function in ←−three variables a, a0 , z, the function f is convex as a sum of convex terms. We want to find thetriple (a, a0 , z) minimizing f .
Boris says: (1) I removed Lemma 3.8 proving convexity of fopt . (2) why do we need the next
sentence? The restriction of f obtained by fixing the variable z is also a convex function in a and ←−a0 . Putf (a, a0 , z).
fopt (z) = max
0a,aIt is convex, as the maximum of a set of convex functions.
Now, suppose (xq , ηopt ) is the location of the optimal joint. For any value η of z, if we cancompute f opt(η), f opt(η +), f opt(η −) for an  > 0, we can decide whether η > ηopt or η < ηopt ,or η −  ≤ ηopt ≤ η + . Indeed, the first case and the second case are when f opt(η − ) < f opt(η)
and f opt(η + ) < f opt(η), respectively, and otherwise the third case occurs. This enables us toapply binary searching and parametric searching for computing ηopt .
For each q, we guess the y-coordinate value η of the joint vertex (xq , η) and compute f opt(η).
For the purpose, we compute the best line, in the sense of L1 fitting, approximating each of
S1 (q) and S2 (q) going through the (for now, fixed) joint by using almost the same strategy as insection 3.2.1.
We focus on the best line approximating S1 (q). It suffices to determine the slope of this line.
In the dual space, we just need to compute a point p = (a(p), b(p)) on the line Y = xq X − η such
12

P
that qi=1 |a(p)xi − b(p) − yi | is minimized. We observe that the above expression is convex as afunction of a, and hence θ(p) = θ+ (p) − θ− (p) is monotone and changes the sign at p, where θ+ (p)(θ− (p)) is the sum of slopes of lines above p (resp. below p). Thus, we can apply binary searchby using slope-sum query, and this binary search can be performed in O(log n) steps by using thefiltration as described in Lemma 3.7. Thus, we can compute f opt(η) in O(τ (n) log2 n) time byusing our semi-dynamic query structures, and in polylogarithmic time using O(τ (n)) processors
in parallel.
Now we have f opt(η) in hand for any given η, and we can apply binary search for computing
the optimal value ηopt . In order to construct a strongly polynomial algorithm, we apply parametricsearch. Given η, our algorithm runs in polylogarithmic time using O(τ (n)) processors. Thus, theparametric search paradigm [25] is applicable here. Therefore, for a fixed q, the second case of
the problem can be handled in Õ(τ (n)) time.
Since there are n candidates of q, the overall time complexity for examining for all q is Õ(τ (n)).
By considering the space-time tradeoff on τ (n), we have the following:
Theorem 3.8. The optimal L1 -fitting 1-joint polyline is computed in Õ(n1.5 ) randomized timeusing quasi-linear space, and Õ(n4/3 ) randomized time using O(n4/3 ) space.

4

Fitting a k-joint polyline

The k-joint fitting problem is polynomial-time solvable for fixed k. We describe the algorithm ina non-deterministic fashion. We guess the partition of x1 , . . . , xn into k intervals each of whichcorresponds to a line segment of the polyline. Also, we guess whether each joint is free or fixed.
We decompose the problem at the free joints and obtain a set of subproblems. In each subproblem,
we add the linear constraints corresponding to the fixed condition (i.e., each joint is located ona guessed vertical line). Thus, each subproblem is a convex programming problem: a linear
program for L1 , and a quadratic program for L2 . We solve each subproblem separately to obtainthe solution of the whole problem. Note that this strategy works because of the convexity of each
subproblem. There are O((3n)k ) different choices of the guesses, thus we can be replace guessing
by a brute-force search to have a polynomial-time deterministic algorithm if k is a constant.
For a general k, we do not know whether the problem is in the class P or not. Thus, we would
like to consider approximation algorithms. One possible approach is to relax the requirementthat number of joints be exactly k. We can design a PTAS for it.
Theorem 4.1. Let zopt be the optimal L1 (or L2 ) error of a k-joint fitting. Then, for any constantε > 0, we can compute a b(1 + ε)kc-joint fitting whose error is at most zopt in polynomial time.
Proof. We ignore continuity and approximate the points by using a piecewise-linear (not necessarily continuous) function with k linear pieces. This can be done by preparing the optimal linear
regression for each subinterval of consecutive points of S, and then applying dynamic programming. We can restore continuity by inserting at most k steep (nearly vertical) line segments. The
resulting polyline has at most 2k joints and error at most zopt . We can improve 2k to b 3k
2 c byapplying the 1-joint algorithm instead of linear regression algorithm, and further improve it tob(1 + ε)kc by using the r-joint algorithm mentioned above for r = dε−1 e.
Another approach is to keep the number of joints at k and approximate the fitting error. We
give a FPTAS for it. We only discuss the L1 case, since the L2 case is analogous. Let zopt be theoptimal L1 -error, and we aim to find a k-joint polyline whose error is at most (1 + ε)zopt . We
remark that if zopt = 0, our solution is exactly the same as the solution for the uniform metricfitting problem, and thus we may assume zopt > 0. Recall that the uniform metric fitting problem
can be solved in O(n log n) time [10]. The following is a trivial but crucial observation:
13

Lemma 4.2. Let z∞ be the optimal error for the uniform metric k-joint fitting problem. Then,z∞ ≤ zopt ≤ nz∞ .
Proof. The sum of the errors in the uniform-metric–optimal polyline is at most nz∞ . Hencenz∞ ≥ zopt . On the other hand, every k-joint polyline has a data point in S such that thevertical distance to the polyline is at least z∞ , so zopt ≥ z∞ .
Our strategy is as follows: We call the n vertical lines through our input points the columnlines. We give a set of portal points on each column line, and call a k-joint polyline a tame polyline
if each of its links satisfies the condition that the line containing the link goes through a pair of
portal points.
On each column line, the distance between its data point and the intersection point with the
optimal polyline is at most zopt , thus at most nz∞ . Thus, on the ith column line, we place theportals in the vertical range [yi − nz∞ , yi + nz∞ ]. The portal points are placed symmetricallyabove and below yi . The jth portal above yi is located at the y-value yi + (1 + 2ε )j−1 δ, where∞ε
and j = 1, 2, . . . , M . We choose the largest M satisfying (1 + 2ε )M δ ≤ nz∞ , and henceδ = z2nM = O(ε−1 log(n + ε−1 )). We also put portals at heights yi and yi ± nz∞ . In this way the numberof portals in any column is at most 2M + 3.
We call a closed interval between adjacent portals in a column a prime interval.
Lemma 4.3. There exists a tame polyline whose L1 error is at most (1 + ε)zopt .
Proof. We start from the optimal polyline `opt , and deform it to obtain a tame polyline. We
proceed sequentially, left to right. Consider the line containing the leftmost link of `opt . We
continuously move the line to a tame line without crossing any portal point during the movement;
if the line started off passing through a portal point, we rotate it around it; if the line started off
passing through two, it is already tame. The right joint of the current link is accordingly moved
to the intersection of the new line and the line containing the right neighbor link. It may happen
that during this transformation a joint crosses a column line. However, the intersection pointsof the original and the deformed polylines with a column line are located in a common prime
interval. We repeat this operation, proceeding from left to right, to obtain a tame polyline.
Now, consider the change of vertical distances between a point pi and the two polylines. The
polylines go through the same prime interval of the column line through pi . An index i is calleda near-index if the polylines goes through a prime interval containing pi in its closure; otherwise
it is called a far-index. For the near-indices, the summation of the errors caused by the newεpolyline is bounded by nδ = z∞
2 . For each far-index, the errors caused by the new polyline at
the column is bounded by 1 + 2ε times the one caused by the old (i.e. optimal) polyline. Thus,
the combined error of the new polyline for all the far indices is at most (1 + 2ε )zopt . In total, theεerror of the new polyline is bounded by (1 + 2ε )zopt + z∞
2 ≤ (1 + ε)zopt .
Thus, it suffices to compute the optimal tame polyline. There are M n portals, and thus
N = O(M 2 n2 ) lines going through a pair of portals. Let L be the set of these lines. We design
a dynamic programming algorithm. For the ith column, for each line ` ∈ L and each m ≤ k,
we record the approximation error of the best m-joint tame polyline up to the current column
whose (rightmost) link covering pi is on `. When we proceed to the (i + 1)th column, eachapproximation error is updated. If there is an intersection between lines ` and `0 in the interval(xi , xi+1 ], we consider the polylines that have the intersection as a possible joint. This can bedone by copying the data for ` to `0 and vice versa incrementing the joint number by one, andthen keeping the smaller of the current and the new (copied) error for each of the pairs (`, m)and (`0 , m) for m = 1, 2, . . . , k. Then, we add the distance from pi+1 to each polyline. Finally,
we select the minimum error at the nth column, and retrieve the polyline by backtracking.

14

There are O(N 2 ) intersections of lines, and it takes O(k) time for each intersection forcopying and updating. This requires O(N 2 k) work and dominates the running time. Since
N = O(n2 M 2 ) = O(n2 ε−2 log2 (n + ε−1 )), we have the following:
Theorem 4.4. A (1+ε)-approximation, i.e., a k-joint polyline with error 1+ε times the optimal,
for each of the L1 and L2 k-joint problems can be computed in O(kn4 ε−4 log4 (n + ε−1 )) time.

5

Concluding remarks

A major open problem is to determine the complexity class of the k-joint problem for L1 - andL2 -fitting. The corresponding L1 or L2 polyline approximation problem where the input is a
curve is also interesting.
We remark that the curve simplification problem under the L1 -measure is to minimize thearea between input and output polylines. In the restricted case where the vertex set or the setof lines supporting edges of the output polyline is required to be a subset of that of the input
polyline, the problem is reduced to the k-link shortest path problem in a graph. In particular, if
the input polyline is convex, this problem is related to matrix searching (see [2]). However, for thegeneral case the authors are not aware of an efficient algorithm, and it is an interesting research
problem. A related question for which polynomial-time algorithms have been constructed for L1
and L2 measures is approximating an x-monotone curve by a k-peak curve, i.e., a curve with at
most k local maxima [6]. The k-joint problem seems to be more difficult than the k-peak problem
because controlling continuity is a challenge in the former case.
Acknowledgment:on convexity.

The authors would like to thank Jirı́ Matoušek for a stimulating discussion

References[1] P. K. Agarwal, S. Hal-Peled, N.H. Mustafa, and Y. Wang, “Near-linear time approximation algorithm
for curve simplification,” Proc. 2002 European Symp. Algorithms, LNCS 2461, (2002) pp.29–41.
[2] A. Aggarwal, B. Schieber, and T. Tokuyama, “Finding a minimum-weight k-link path in graphs withthe concave Monge property and applications,” Discrete Comput. Geom., 12 (1994) 263–280.
[3] P. Agarwal, B. Aronov, T. Chan, M. Sharir, “On levels in arrangements of lines, segments, planes,
and triangles,” Discrete Comput. Geom., 19 (1998) 315–331.
[4] P. Agarwal and J. Matoušek, “Ray shooting and parametric search,” SIAM J. Comput., 22 (1993)
794–806.
[5] P. K. Agarwal and K. R. Varadarajan, “Efficient algorithms for approximating polygonal chains,”
Discrete Comput. Geom., 23 (2000) 273–291.
[6] J. Chun, K. Sadakane, and T. Tokuyama, “Linear time algorithm for approximating a curve by asingle-peaked curve,” Proc. 14th Internat. Symp. Algorithms Comput. (ISAAC 2003), LNCS 2906,
2003, pp. 6–16.
[7] H. Edelsbrunner, Algorithms in Combinatorial Geometry, ETACS Monographs on Theoretical Computer Science 10, Springer-Verlag, 1987.
[8] D. Eu and G.T. Toussaint, “On approximating polygonal curves in two and three dimensions,”
CVGIP: Graphical Models And Image processing 56 (1994) 231–246.
[9] K. L. Clarkson and P. W. Shor, “Application of random sampling in computational geometry,” Discrete Comput. Geom., 4 (1989) 423–432.
[10] M. Goodrich, “Efficient piecewise-linear function approximation using the uniform metric,” Discrete
Comput. Geom., 14 (1995) 445–462.

15[11] S. Hakimi and E. Schmeichel, “Fitting polygonal functions to a set of points in the plane,” GraphicalModels and Image Processing, 53 (1991) 132–136.
[12] H. Imai and M. Iri: “Polygonal approximations of a curve - Formulations and algorithms,” Computational Morphology, Elsevier Science Publishers B.V. (North Holland), 1988, 71–86.
[13] H. Imai, K. Kato, and P. Yamamoto: “A linear-time algorithm for linear L1 approximation of points,”
Algorithmica, 4 (1989) 77–96.
[14] N. Katoh and T. Tokuyama, “Notes on computing peaks in k-levels and parametric spanning trees,”
Proc. 17th ACM Symp. on Computational Geometry, 2001, pp. 241–248.
[15] Y. Kurozumi and W.A. Davis: “Polygonal approximation by the minimax method,” Computer Graphics and Image Processing, 19 (1982) 248–264.
[16] S. Langerman and W. Steiger, “Optimization in arrangements.” Proc. Symp. Theor. Aspects Computer
Science (STACS2003), LNCS 2607, 2003, pp. 50–61.
[17] J. Matoušek, “Efficient partition trees,” Discrete Comput. Geom., 8 (1992) 315–334.
[18] J. Matoušek, “Range searching with efficient hierarchical cutting,” Proc. 8th ACM Symp. on Comput.
Geom., (1992) 276–287.
[19] J. Matoušek, “Geometric range searching,” ACM Computing Surveys, 26 (1994) 421–461.
[20] J. Matoušek and O. Schwarzkopf, “Linear optimization queries,” Proc. 8th Annual ACM Symp. Comput. Geom., 1992, pp. 16–25.
[21] N. Megiddo and A. Tamir, “Finding least-distances lines,” SIAM J. Algebraic Discrete Methods, 4(1983) pp. 207–211.
[22] A. Melkman and J. O’Rourke, “On polygonal chain approximation,” Computational Morphology,
Elsevier Science Publishers B.V. (North Holland), 1988, 87–95.
[23] J. O’Rourke and G. Toussaint, “Pattern recognition,” Chapter 43 of Handbook of Discrete and Computational Geometry (eds. J. Goodman and J. O’Rourke), CRC Press, 1997.
[24] F. P. Preparata and M. I. Shamos, Computational Geometry, an Introduction, Springer-Verlag, New
York, 1985.
[25] J. Salowe, “Parametric search,” Chapter 37 of Handbook of Discrete and Computational Geometry(eds. J. Goodman and J. O’Rourke), CRC Press, 1997.
[26] T. Tokuyama, “Minimax parametric optimization problems and multi-dimensional parametric searching,” Proc. 33rd Symp. Theory Comput. (2001), pp. 75–83.
[27] D. P. Wang, N. F. Huang, H. S. Chao, and R. C. T. Lee, “Plane sweep algorithms for polygonal
approximation problems with applications,” Proc. 4th Internat. Symp. Algorithms Comput. (ISAAC
2003), LNCS 762, 1993, pp. 515–522.
[28] Robert Weibel, “Generalization of spatial data: principles and selected algorithms,” Algorithmic
Foundations of Geographic Information Systems, LNCS 1340, 1997, pp. 99–152.

16