]> git.donarmstrong.com Git - lilypond.git/blobdiff - lily/skyline.cc
Adds an explanatory comment to skyline.cc.
[lilypond.git] / lily / skyline.cc
index 5eedbe25373f39f53da02f671c7b00f6826c5074..3f52d82e9040fd65b104bd98593c034265fbd19d 100644 (file)
@@ -1,12 +1,26 @@
-/* skyline.cc -- implement the Skyline class
+/*
+  This file is part of LilyPond, the GNU music typesetter.
+
+  Copyright (C) 2006--2012 Joe Neeman <joeneeman@gmail.com>
+
+  LilyPond is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  LilyPond is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
 
-   source file of the GNU LilyPond music typesetter
-   (c) 2006--2007 Joe Neeman <joeneeman@gmail.com>
+  You should have received a copy of the GNU General Public License
+  along with LilyPond.  If not, see <http://www.gnu.org/licenses/>.
 */
 
 #include "skyline.hh"
+#include "skyline-pair.hh"
 #include <deque>
+#include <cstdio>
 
 #include "ly-smobs.icc"
 
    a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
    This means that the merging routine doesn't need to be aware of direction,
    but the distance routine does.
-*/
-
-/*
-  FIXME:
-
-  * Consider to use
-
-  typedef list<Skyline_point> Skyline;
-  struct Skyline_point
-  {
-    Real x;
-    Drul_array<Real> ys; 
-  };
-
-  this is a cleaner representation, as it doesn't duplicate the X, and
-  doesn't need bogus buildings at infinity  --hwn.
-
-
-  * All the messing around with EPS is very fishy.  There are no
-  complicated numerical algorithms involved, so EPS should not be
-  necessary.
-
-  --hwn
-  
-  
- */
 
-#define EPS 1e-10
+   Be careful about numerical accuracy. When dealing with extremely small values,
+   computation errors may arise due to the use of floating point arithmetic.
+   For example, if left and right have equal values to start with, in C++
+   they may not receive the same value after
 
-static inline bool
-approx_equal (Real x, Real y)
-{
-  return abs (x - y) < EPS || (isinf (x) && isinf (y) && ((x > 0) == (y > 0)));
-}
+   left = left*factor + offset;
+   right = right*factor + offset;
 
-static inline bool
-approx_greater_than (Real x, Real y)
-{
-  return x > y + EPS;
-}
-
-static inline bool
-approx_less_than (Real x, Real y)
-{
-  return x < y - EPS;
-}
+   Which is very unfortunate. Maybe using GCC compiler options to disallow
+   extended precision for intermediate results and/or the choice to store
+   intermediates with less than full precision would retain some kind of
+   deterministic behavior that way.
 
-static inline bool
-approx_less_equal (Real x, Real y)
-{
-  return x <= y + EPS;
-}
+   Anyway, it seems that accepting extremely narrow building in skylines
+   doesn't cause accuracy problems to us, so we allow arbitrarily small buildings.
+   However, as Keith pointed out, problems may appear if more than one operation
+   is done before storing the result, and/or there are different code paths
+   for doing the operations to the different ends of an interval.
+*/
 
-static inline bool
-approx_greater_equal (Real x, Real y)
+static void
+print_buildings (list<Building> const &b)
 {
-  return x >= y - EPS;
+  for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
+    i->print ();
 }
 
 void
 Skyline::print () const
 {
-  for (list<Building>::const_iterator i = buildings_.begin ();
-       i != buildings_.end (); i++)
-    {
-      (*i).print ();
-    }
+  print_buildings (buildings_);
 }
 
-bool
-is_legal_skyline (list<Building> const &buildings)
+void
+Skyline::print_points () const
 {
-  list<Building>::const_iterator i;
-  Real last_x = -infinity_f;
-  for (i = buildings.begin (); i != buildings.end (); i++)
-    {
-      if (i->iv_[LEFT] != last_x)
-       return false;
-      last_x = i->iv_[RIGHT];
-      if (isinf (i->iv_.length ()) && i->height_[LEFT] != i->height_[RIGHT])
-       return false;
-    }
-  return last_x == infinity_f;
+  vector<Offset> ps (to_points (X_AXIS));
+
+  for (vsize i = 0; i < ps.size (); i++)
+    printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
+            (i % 2) == 1 ? "\n" : " ");
 }
 
 Building::Building (Real start, Real start_height, Real end_height, Real end)
-  : iv_ (start, end)
 {
-  height_[LEFT] = start_height;
-  height_[RIGHT] = end_height;
-
   if (isinf (start) || isinf (end))
     assert (start_height == end_height);
 
-  precompute ();
+  start_ = start;
+  end_ = end;
+  precompute (start, start_height, end_height, end);
 }
 
-Building::Building (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
+Building::Building (Box const &b, Axis horizon_axis, Direction sky)
 {
+  Real start = b[horizon_axis][LEFT];
+  Real end = b[horizon_axis][RIGHT];
   Real height = sky * b[other_axis (horizon_axis)][sky];
 
-  iv_ = b[horizon_axis];
-  iv_.widen (horizon_padding + EPS);
-  height_[LEFT] = height;
-  height_[RIGHT] = height;
-
-  if (sane ())
-    precompute ();
+  start_ = start;
+  end_ = end;
+  precompute (start, height, height, end);
 }
 
 void
-Building::precompute ()
+Building::precompute (Real start, Real start_height, Real end_height, Real end)
 {
-  slope_ = (height_[RIGHT] - height_[LEFT]) / (iv_.length());
-  if (height_[LEFT] == height_[RIGHT]) /* in case they're both infinity */
-    slope_ = 0;
+  slope_ = 0.0; /* if they were both infinite, we would get nan, not 0, from the prev line */
+  if (start_height != end_height)
+    slope_ = (end_height - start_height) / (end - start);
 
   assert (!isinf (slope_) && !isnan (slope_));
 
-  if (isinf (iv_[START]))
+  if (isinf (start))
     {
-      assert (slope_ == 0);
-      y_intercept_ = height_[LEFT];
+      assert (start_height == end_height);
+      y_intercept_ = start_height;
     }
   else
-    y_intercept_ = height_[LEFT] - slope_ * iv_[START];
+    y_intercept_ = start_height - slope_ * start;
 }
 
-Real 
+inline Real
 Building::height (Real x) const
 {
-  if (isinf (x))
-    return (x > 0) ? height_[RIGHT] : height_[LEFT];
-  return slope_*x + y_intercept_;
+  return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
 }
 
 void
 Building::print () const
 {
-  printf ("X[%f,%f] -> Y[%f,%f]\n",
-         iv_[LEFT], iv_[RIGHT],
-         height_[LEFT], height_[RIGHT]);
+  printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
 }
 
-Real
+inline Real
 Building::intersection_x (Building const &other) const
 {
-  return (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
+  Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
+  return isnan (ret) ? -infinity_f : ret;
 }
 
 void
 Building::leading_part (Real chop)
 {
-  assert (chop > iv_[LEFT] && chop <= iv_[RIGHT] && !approx_equal (chop, iv_[LEFT]));
-  iv_[RIGHT] = chop;
-  height_[RIGHT] = height (chop);
+  assert (chop <= end_);
+  end_ = chop;
 }
 
-Building
-Building::sloped_neighbour (Real horizon_padding, Direction d) const
+// Returns a shift s such that (x + s, y) intersects the roof of
+// this building.  If no such shift exists, returns infinity_f.
+Real
+Building::shift_to_intersect (Real x, Real y) const
 {
-  Real left = iv_[d];
-  Real right = iv_[d] + d * horizon_padding;
-  Real left_height = height_[d];
-  Real right_height = height_[d] - horizon_padding;
-  if (d == LEFT)
+  // Solve for s: y = (x + s)*m + b
+  Real ret = (y - y_intercept_ - slope_ * x) / slope_;
+
+  if (ret >= start_ && ret <= end_ && !isinf (ret))
+    return ret;
+  return infinity_f;
+}
+
+static Real
+first_intersection (Building const &b, list<Building> *const s, Real start_x)
+{
+  while (!s->empty () && start_x < b.end_)
     {
-      swap (left, right);
-      swap (left_height, right_height);
+      Building c = s->front ();
+
+      // conceals and intersection_x involve multiplication and
+      // division. Avoid that, if we can.
+      if (c.y_intercept_ == -infinity_f)
+        {
+          if (c.end_ > b.end_)
+            return b.end_;
+          start_x = c.end_;
+          s->pop_front ();
+          continue;
+        }
+
+      if (c.conceals (b, start_x))
+        return start_x;
+
+      Real i = b.intersection_x (c);
+      if (i > start_x && i <= b.end_ && i <= c.end_)
+        return i;
+
+      start_x = c.end_;
+      if (b.end_ > c.end_)
+        s->pop_front ();
     }
-  return Building (left, left_height, right_height, right);
+  return b.end_;
 }
 
 bool
-Building::sane () const
+Building::conceals (Building const &other, Real x) const
 {
-  return approx_less_than (iv_[LEFT], iv_[RIGHT])
-    && !isinf (height_[RIGHT])
-    && !isinf (height_[LEFT]);
+  if (slope_ == other.slope_)
+    return y_intercept_ > other.y_intercept_;
+
+  /* their slopes were not equal, so there is an intersection point */
+  Real i = intersection_x (other);
+  return (i <= x && slope_ > other.slope_)
+         || (i > x && slope_ < other.slope_);
 }
 
-static void
-skyline_trailing_part (list<Building> *sky, Real x)
+// Remove redundant empty buildings from the skyline.
+// If there are two adjacent empty buildings, they can be
+// turned into one.
+void
+Skyline::normalize ()
 {
-  if (approx_equal (x, sky->front ().iv_[RIGHT]))
-    sky->pop_front ();
-  else
-    assert (x < sky->front ().iv_[RIGHT]);
+  bool last_empty = false;
+  list<Building>::iterator i;
 
-  if (!sky->empty ())
+  for (i = buildings_.begin (); i != buildings_.end (); i++)
     {
-      sky->front ().iv_[LEFT] = x;
-      sky->front ().height_[LEFT] = sky->front ().height (x);
+      if (last_empty && i->y_intercept_ == -infinity_f)
+        {
+          list<Building>::iterator last = i;
+          last--;
+          last->end_ = i->end_;
+          buildings_.erase (i);
+          i = last;
+        }
+      last_empty = (i->y_intercept_ == -infinity_f);
     }
-}
 
-bool
-Building::conceals_beginning (Building const &other) const
-{
-  bool w = false;
-  Real h = other.height (iv_[LEFT]);
-  if (approx_equal (height_[LEFT], h))
-    w = slope_ > other.slope_;    
-  else if (height_[LEFT] > h) 
-    w = true;
-  else 
-    w = false;
-
-  return w;
+  assert (buildings_.front ().start_ == -infinity_f);
+  assert (buildings_.back ().end_ == infinity_f);
 }
 
-bool
-Building::conceals (Building const &other) const
+void
+Skyline::deholify ()
 {
-  assert (iv_[LEFT] <= other.iv_[LEFT]);
-  return (iv_[RIGHT] >= other.iv_[RIGHT])
-    && approx_greater_equal (height (other.iv_[LEFT]), other.height_[LEFT])
-    && approx_greater_equal (height (other.iv_[RIGHT]), other.height_[RIGHT]);
+  // Since a skyline should always be normalized, we can
+  // assume that there are never two adjacent empty buildings.
+  // That is, if center is empty then left and right are not.
+  list<Building>::iterator left = buildings_.begin ();
+  list<Building>::iterator center = buildings_.begin ();
+  list<Building>::iterator right;
+
+  for (right = buildings_.begin (); right != buildings_.end (); right++)
+    {
+      if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
+        {
+          Real p1 = left->height (left->end_);
+          Real p2 = right->height (right->start_);
+          *center = Building (center->start_, p1, p2, center->end_);
+
+          left = center;
+          center = right;
+        }
+    }
 }
 
 void
 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
-                                list<Building> *const result)
+                                 list<Building> *const result) const
 {
+  if (s1->empty () || s2->empty ())
+    {
+      programming_error ("tried to merge an empty skyline");
+      return;
+    }
+
+  Real x = -infinity_f;
+  Real last_end = -infinity_f;
   while (!s1->empty ())
     {
-      if (s2->front ().conceals_beginning (s1->front ()))
-       swap (s1, s2);
+      if (s2->front ().conceals (s1->front (), x))
+        swap (s1, s2);
 
       Building b = s1->front ();
-      while (!s2->empty () && b.conceals (s2->front ()))
-       s2->pop_front ();
+      Building c = s2->front ();
+
+      // Optimization: if the other skyline is empty at this point,
+      // we can avoid testing some intersections. Just grab as many
+      // buildings from s1 as we can, and shove them onto the output.
+      if (c.y_intercept_ == -infinity_f
+          && c.end_ >= b.end_)
+        {
+          list<Building>::iterator i = s1->begin ();
+          i++;
+          while (i != s1->end () && i->end_ <= c.end_)
+            i++;
+
+          s1->front ().start_ = x;
+          result->splice (result->end (), *s1, s1->begin (), i);
+          x = result->back ().end_;
+          last_end = x;
+          continue;
+        }
+
+      Real end = first_intersection (b, s2, x);
       if (s2->empty ())
-       {
-         result->push_front (b);
-         break;
-       }
-
-      /* s2 either intersects with b or it ends after b */
-      Real end = infinity_f;
-      Real s2_start_height = s2->front ().height_[LEFT];
-      Real s2_end_height = s2->front ().height_[RIGHT];
-      Real s1_start_height = b.height (s2->front ().iv_[LEFT]);
-      Real s1_end_height = b.height (s2->front ().iv_[RIGHT]);
-      if (approx_greater_than (s2_start_height, s1_start_height))
-       end = s2->front ().iv_[LEFT];
-      else if (approx_greater_than (s2_end_height, s1_end_height))
-       end = b.intersection_x (s2->front ());
-      end = min (end, b.iv_[RIGHT]);
-
-      b.leading_part (end);
-      result->push_front (b);
-
-      skyline_trailing_part (s1, end);
-      skyline_trailing_part (s2, end);
+        {
+          b.start_ = last_end;
+          result->push_back (b);
+          break;
+        }
+
+      if (end >= x)
+        {
+          b.leading_part (end);
+          b.start_ = last_end;
+          last_end = b.end_;
+          result->push_back (b);
+        }
+
+      if (end >= s1->front ().end_)
+        s1->pop_front ();
+
+      x = end;
     }
-  result->reverse ();
 }
 
 static void
@@ -305,61 +337,87 @@ empty_skyline (list<Building> *const ret)
   ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
 }
 
+/*
+  Given Building 'b', build a skyline containing only that building.
+*/
 static void
-single_skyline (Building b, Real horizon_padding, list<Building> *const ret)
-{
-  b.iv_.widen (horizon_padding);
-  
-  if (!isinf (b.iv_[RIGHT]))
-    ret->push_front (Building (b.iv_[RIGHT], -infinity_f,
-                              -infinity_f, infinity_f));
-  if (horizon_padding > 0 && !isinf (b.iv_.length ()))
-    ret->push_front (b.sloped_neighbour (horizon_padding, RIGHT));
-  
-  if (b.iv_[RIGHT] > b.iv_[LEFT])
-    ret->push_front (b);
-
-  if (horizon_padding > 0 && !isinf (b.iv_.length ()))
-    ret->push_front (b.sloped_neighbour (horizon_padding, LEFT));
-  if (!isinf (b.iv_[LEFT]))
-    ret->push_front (Building (-infinity_f, -infinity_f,
-                              -infinity_f, b.iv_[LEFT]));
-}
-
-/* remove a non-overlapping set of buildings from BUILDINGS and build a skyline
+single_skyline (Building b, list<Building> *const ret)
+{
+  assert (b.end_ >= b.start_);
+
+  if (b.start_ != -infinity_f)
+    ret->push_back (Building (-infinity_f, -infinity_f,
+                              -infinity_f, b.start_));
+  ret->push_back (b);
+  if (b.end_ != infinity_f)
+    ret->push_back (Building (b.end_, -infinity_f,
+                              -infinity_f, infinity_f));
+}
+
+/* remove a non-overlapping set of boxes from BOXES and build a skyline
    out of them */
 static list<Building>
 non_overlapping_skyline (list<Building> *const buildings)
 {
   list<Building> result;
   Real last_end = -infinity_f;
+  Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
   list<Building>::iterator i = buildings->begin ();
   while (i != buildings->end ())
     {
-      if (approx_less_than (i->iv_[LEFT], last_end))
-       {
-         i++;
-         continue;
-       }
-
-      if (approx_greater_than (i->iv_[LEFT], last_end))
-       result.push_back (Building (last_end, -infinity_f, -infinity_f, i->iv_[LEFT]));
-      else
-       i->iv_[LEFT] = last_end;
-
-      last_end = i->iv_[RIGHT];
-      list<Building>::iterator j = i;
-      i++;
-      result.splice (result.end (), *buildings, j);
+      Real x1 = i->start_;
+      Real y1 = i->height (i->start_);
+      Real x2 = i->end_;
+      Real y2 = i->height (i->end_);
+
+      // Drop buildings that will obviously have no effect.
+      if (last_building.height (x1) >= y1
+          && last_building.end_ >= x2
+          && last_building.height (x2) >= y2)
+        {
+          list<Building>::iterator j = i++;
+          buildings->erase (j);
+          continue;
+        }
+
+      if (x1 < last_end)
+        {
+          i++;
+          continue;
+        }
+
+      if (x1 >= last_end)
+        result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
+
+      result.push_back (*i);
+      last_building = *i;
+      last_end = i->end_;
+
+      list<Building>::iterator j = i++;
+      buildings->erase (j);
     }
+
   if (last_end < infinity_f)
     result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
-  assert (is_legal_skyline (result));
   return result;
 }
 
+class LessThanBuilding
+{
+public:
+  bool operator () (Building const &b1, Building const &b2)
+  {
+    return b1.start_ < b2.start_
+           || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
+  }
+};
+
+/**
+   BUILDINGS is a list of buildings, but they could be overlapping
+   and in any order.  The returned list of buildings is ordered and non-overlapping.
+*/
 list<Building>
-Skyline::internal_build_skyline (list<Building> *buildings)
+Skyline::internal_build_skyline (list<Building> *buildings) const
 {
   vsize size = buildings->size ();
 
@@ -372,16 +430,16 @@ Skyline::internal_build_skyline (list<Building> *buildings)
   else if (size == 1)
     {
       list<Building> result;
-      single_skyline (buildings->front (), 0, &result);
+      single_skyline (buildings->front (), &result);
       return result;
     }
 
   deque<list<Building> > partials;
-  buildings->sort ();
+  buildings->sort (LessThanBuilding ());
   while (!buildings->empty ())
     partials.push_back (non_overlapping_skyline (buildings));
 
-  /* we'd like to say while (partials->size () > 1) but that's O(n).
+  /* we'd like to say while (partials->size () > 1) but that's O (n).
      Instead, we exit in the middle of the loop */
   while (!partials.empty ())
     {
@@ -389,7 +447,7 @@ Skyline::internal_build_skyline (list<Building> *buildings)
       list<Building> one = partials.front ();
       partials.pop_front ();
       if (partials.empty ())
-       return one;
+        return one;
 
       list<Building> two = partials.front ();
       partials.pop_front ();
@@ -403,13 +461,14 @@ Skyline::internal_build_skyline (list<Building> *buildings)
 Skyline::Skyline ()
 {
   sky_ = UP;
-  empty_skyline (&buildings_);  
+  empty_skyline (&buildings_);
 }
 
 Skyline::Skyline (Skyline const &src)
 {
   sky_ = src.sky_;
+
+  /* doesn't a list's copy constructor do this? -- jneem */
   for (list<Building>::const_iterator i = src.buildings_.begin ();
        i != src.buildings_.end (); i++)
     {
@@ -424,42 +483,86 @@ Skyline::Skyline (Direction sky)
 }
 
 /*
-  build skyline from a set of boxes. If horizon_padding > 0, expand all the boxes
-  by that amount and add 45-degree sloped boxes to the edges of each box (of
-  width horizon_padding). That is, the total amount of horizontal expansion is
-  horizon_padding*4, half of which is sloped and half of which is flat.
+  Build skyline from a set of boxes.
 
-  Boxes should have fatness in the horizon_axis (after they are expanded by
-  horizon_padding), otherwise they are ignored.
+  Boxes should be non-empty on both axes.  Otherwise, they will be ignored
  */
-Skyline::Skyline (vector<Box> const &boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
+Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
 {
-  list<Building> bldgs;
+  list<Building> buildings;
   sky_ = sky;
 
   for (vsize i = 0; i < boxes.size (); i++)
+    if (!boxes[i].is_empty ())
+      buildings.push_front (Building (boxes[i], horizon_axis, sky));
+
+  buildings_ = internal_build_skyline (&buildings);
+  normalize ();
+}
+
+/*
+  build skyline from a set of line segments.
+
+  Segments can be articulated from left to right or right to left.
+  In the case of the latter, they will be stored internally as left to right.
+ */
+Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
+{
+  list<Building> buildings;
+  sky_ = sky;
+
+  for (vsize i = 0; i < segments.size (); i++)
     {
-      Building front (boxes[i], horizon_padding, horizon_axis, sky);
-      if (front.sane ())
-       {
-         bldgs.push_front (front);
-         if (horizon_padding > 0 && !isinf (front.iv_.length ()))
-           {
-             bldgs.push_front (front.sloped_neighbour (horizon_padding, LEFT));
-             bldgs.push_front (front.sloped_neighbour (horizon_padding, RIGHT));
-           }
-       }
+      Drul_array<Offset> const &seg = segments[i];
+      Offset left = seg[LEFT];
+      Offset right = seg[RIGHT];
+      if (left[horizon_axis] > right[horizon_axis])
+        swap (left, right);
+
+      Real x1 = left[horizon_axis];
+      Real x2 = right[horizon_axis];
+      Real y1 = left[other_axis (horizon_axis)] * sky;
+      Real y2 = right[other_axis (horizon_axis)] * sky;
+
+      if (x1 <= x2)
+        buildings.push_back (Building (x1, y1, y2, x2));
     }
-  
-  buildings_ = internal_build_skyline (&bldgs);
-  assert (is_legal_skyline (buildings_));
+
+  buildings_ = internal_build_skyline (&buildings);
+  normalize ();
 }
 
-Skyline::Skyline (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
+Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
 {
   sky_ = sky;
-  Building front (b, 0, horizon_axis, sky);
-  single_skyline (front, horizon_padding, &buildings_);
+
+  deque<Skyline> partials;
+  for (vsize i = 0; i < skypairs.size (); i++)
+    partials.push_back (Skyline ((skypairs[i])[sky]));
+
+  while (partials.size () > 1)
+    {
+      Skyline one = partials.front ();
+      partials.pop_front ();
+      Skyline two = partials.front ();
+      partials.pop_front ();
+
+      one.merge (two);
+      partials.push_back (one);
+    }
+
+  if (partials.size ())
+    buildings_.swap (partials.front ().buildings_);
+  else
+    buildings_.clear ();
+}
+
+Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
+{
+  sky_ = sky;
+  Building front (b, horizon_axis, sky);
+  single_skyline (front, &buildings_);
+  normalize ();
 }
 
 void
@@ -467,23 +570,43 @@ Skyline::merge (Skyline const &other)
 {
   assert (sky_ == other.sky_);
 
+  if (other.is_empty ())
+    return;
+
+  if (is_empty ())
+    {
+      buildings_ = other.buildings_;
+      return;
+    }
+
   list<Building> other_bld (other.buildings_);
   list<Building> my_bld;
   my_bld.splice (my_bld.begin (), buildings_);
   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
-  assert (is_legal_skyline (buildings_));
+  normalize ();
 }
 
 void
-Skyline::insert (Box const &b, Real horizon_padding, Axis a)
+Skyline::insert (Box const &b, Axis a)
 {
   list<Building> other_bld;
   list<Building> my_bld;
 
+  if (isnan (b[other_axis (a)][LEFT])
+      || isnan (b[other_axis (a)][RIGHT]))
+    {
+      programming_error ("insane box for skyline");
+      return;
+    }
+
+  /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
+  if (b.is_empty ())
+    return;
+
   my_bld.splice (my_bld.begin (), buildings_);
-  single_skyline (Building (b, 0, a, sky_), horizon_padding, &other_bld);
+  single_skyline (Building (b, a, sky_), &other_bld);
   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
-  assert (is_legal_skyline (buildings_));
+  normalize ();
 }
 
 void
@@ -491,43 +614,139 @@ Skyline::raise (Real r)
 {
   list<Building>::iterator end = buildings_.end ();
   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
-    {
-      i->height_[LEFT] += sky_ * r;
-      i->height_[RIGHT] += sky_ * r;
-      i->y_intercept_ += sky_ * r;
-    }
-  assert (is_legal_skyline (buildings_));
+    i->y_intercept_ += sky_ * r;
 }
 
 void
-Skyline::shift (Real r)
+Skyline::shift (Real s)
 {
   list<Building>::iterator end = buildings_.end ();
   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
     {
-      i->iv_[LEFT] += r;
-      i->iv_[RIGHT] += r;
+      i->start_ += s;
+      i->end_ += s;
+      i->y_intercept_ -= s * i->slope_;
     }
 }
 
 Real
-Skyline::distance (Skyline const &other) const
+Skyline::distance (Skyline const &other, Real horizon_padding) const
+{
+  Real dummy;
+  return internal_distance (other, horizon_padding, &dummy);
+}
+
+Real
+Skyline::touching_point (Skyline const &other, Real horizon_padding) const
+{
+  Real touch;
+  internal_distance (other, horizon_padding, &touch);
+  return touch;
+}
+
+Real
+Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
+{
+  if (horizon_padding == 0.0)
+    return internal_distance (other, touch_point);
+
+  // Note that it is not necessary to build a padded version of other,
+  // because the same effect can be achieved just by doubling horizon_padding.
+  Skyline padded_this = padded (horizon_padding);
+  return padded_this.internal_distance (other, touch_point);
+}
+
+Skyline
+Skyline::padded (Real horizon_padding) const
+{
+  if (horizon_padding < 0.0)
+    warning ("Cannot have negative horizon padding.  Junking.");
+
+  if (horizon_padding <= 0.0)
+    return *this;
+
+  list<Building> pad_buildings;
+  for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
+    {
+      if (i->start_ > -infinity_f)
+        {
+          Real height = i->height (i->start_);
+          if (height > -infinity_f)
+            {
+              // Add the sloped building that pads the left side of the current building.
+              Real start = i->start_ - 2 * horizon_padding;
+              Real end = i->start_ - horizon_padding;
+              pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
+
+              // Add the flat building that pads the left side of the current building.
+              start = i->start_ - horizon_padding;
+              end = i->start_;
+              pad_buildings.push_back (Building (start, height, height, end));
+            }
+        }
+
+      if (i->end_ < infinity_f)
+        {
+          Real height = i->height (i->end_);
+          if (height > -infinity_f)
+            {
+              // Add the flat building that pads the right side of the current building.
+              Real start = i->end_;
+              Real end = start + horizon_padding;
+              pad_buildings.push_back (Building (start, height, height, end));
+
+              // Add the sloped building that pads the right side of the current building.
+              start = end;
+              end += horizon_padding;
+              pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
+            }
+        }
+    }
+
+  // The buildings may be overlapping, so resolve that.
+  list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
+
+  // Merge the padding with the original, to make a new skyline.
+  Skyline padded (sky_);
+  list<Building> my_buildings = buildings_;
+  padded.buildings_.clear ();
+  internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
+  padded.normalize ();
+
+  return padded;
+}
+
+Real
+Skyline::internal_distance (Skyline const &other, Real *touch_point) const
 {
   assert (sky_ == -other.sky_);
+
   list<Building>::const_iterator i = buildings_.begin ();
   list<Building>::const_iterator j = other.buildings_.begin ();
 
   Real dist = -infinity_f;
+  Real start = -infinity_f;
+  Real touch = -infinity_f;
   while (i != buildings_.end () && j != other.buildings_.end ())
     {
-      Interval iv = intersection (i->iv_, j->iv_);
-      dist = max (dist, max (i->height (iv[LEFT]) + j->height (iv[LEFT]),
-                            i->height (iv[RIGHT]) + j->height (iv[RIGHT])));
-      if (i->iv_[RIGHT] <= j->iv_[RIGHT])
-       i++;
+      Real end = min (i->end_, j->end_);
+      Real start_dist = i->height (start) + j->height (start);
+      Real end_dist = i->height (end) + j->height (end);
+      dist = max (dist, max (start_dist, end_dist));
+
+      if (end_dist == dist)
+        touch = end;
+      else if (start_dist == dist)
+        touch = start;
+
+      if (i->end_ <= j->end_)
+        i++;
       else
-       j++;
+        j++;
+      start = end;
     }
+
+  *touch_point = touch;
   return dist;
 }
 
@@ -539,8 +758,8 @@ Skyline::height (Real airplane) const
   list<Building>::const_iterator i;
   for (i = buildings_.begin (); i != buildings_.end (); i++)
     {
-      if (i->iv_[RIGHT] >= airplane)
-       return sky_ * i->height (airplane);
+      if (i->end_ >= airplane)
+        return sky_ * i->height (airplane);
     }
 
   assert (0);
@@ -550,121 +769,109 @@ Skyline::height (Real airplane) const
 Real
 Skyline::max_height () const
 {
-  Skyline s (-sky_);
-  s.set_minimum_height (0);
-  return sky_ * distance (s);
+  Real ret = -infinity_f;
+
+  list<Building>::const_iterator i;
+  for (i = buildings_.begin (); i != buildings_.end (); ++i)
+    {
+      ret = max (ret, i->height (i->start_));
+      ret = max (ret, i->height (i->end_));
+    }
+
+  return sky_ * ret;
 }
 
-void
-Skyline::set_minimum_height (Real h)
+Direction
+Skyline::direction () const
 {
-  Skyline s (sky_);
-  s.buildings_.front ().height_[LEFT] = h * sky_;
-  s.buildings_.front ().height_[RIGHT] = h * sky_;
-  s.buildings_.front ().y_intercept_ = h * sky_;
-  merge (s);
+  return sky_;
 }
 
-
-vector<Offset>
-Skyline::to_points () const
+Real
+Skyline::left () const
 {
-  vector<Offset> out;
-
   for (list<Building>::const_iterator i (buildings_.begin ());
        i != buildings_.end (); i++)
-    {
-      if (!isinf (i->iv_[LEFT]) && !isinf (i->height_[LEFT]))
-       out.push_back (Offset (i->iv_[LEFT], sky_ * i->height_[LEFT]));
-      if (!isinf (i->iv_[RIGHT]) && !isinf (i->height_[RIGHT]))
-       out.push_back (Offset (i->iv_[RIGHT], sky_ * i->height_[RIGHT]));
-    }
-  return out;
-}
+    if (i->y_intercept_ > -infinity_f)
+      return i->start_;
 
-bool
-Skyline::is_empty () const
-{
-  return buildings_.empty ();
+  return infinity_f;
 }
 
-Skyline_pair::Skyline_pair ()
-  : skylines_ (Skyline (DOWN), Skyline (UP))
+Real
+Skyline::right () const
 {
-}
+  for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
+       i != buildings_.rend (); ++i)
+    if (i->y_intercept_ > -infinity_f)
+      return i->end_;
 
-Skyline_pair::Skyline_pair (vector<Box> const &boxes, Real padding, Axis a)
-  : skylines_ (Skyline (boxes, padding, a, DOWN), Skyline (boxes, padding, a, UP))
-{
+  return -infinity_f;
 }
 
-Skyline_pair::Skyline_pair (Box const &b, Real padding, Axis a)
-  : skylines_ (Skyline (b, padding, a, DOWN), Skyline (b, padding, a, UP))
+Real
+Skyline::max_height_position () const
 {
+  Skyline s (-sky_);
+  s.set_minimum_height (0);
+  return touching_point (s);
 }
 
 void
-Skyline_pair::raise (Real r)
+Skyline::set_minimum_height (Real h)
 {
-  skylines_[UP].raise (r);
-  skylines_[DOWN].raise (r);
+  Skyline s (sky_);
+  s.buildings_.front ().y_intercept_ = h * sky_;
+  merge (s);
 }
 
-void
-Skyline_pair::shift (Real r)
+vector<Offset>
+Skyline::to_points (Axis horizon_axis) const
 {
-  skylines_[UP].shift (r);
-  skylines_[DOWN].shift (r);
-}
+  vector<Offset> out;
 
-void
-Skyline_pair::insert (Box const &b, Real padding, Axis a)
-{
-  skylines_[UP].insert (b, padding, a);
-  skylines_[DOWN].insert (b, padding, a);
-}
+  Real start = -infinity_f;
+  for (list<Building>::const_iterator i (buildings_.begin ());
+       i != buildings_.end (); i++)
+    {
+      out.push_back (Offset (start, sky_ * i->height (start)));
+      out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
+      start = i->end_;
+    }
 
-void
-Skyline_pair::merge (Skyline_pair const &other)
-{
-  skylines_[UP].merge (other[UP]);
-  skylines_[DOWN].merge (other[DOWN]);
-}
+  if (horizon_axis == Y_AXIS)
+    for (vsize i = 0; i < out.size (); i++)
+      out[i] = out[i].swapped ();
 
-bool
-Skyline_pair::is_empty () const
-{
-  return skylines_[UP].is_empty ()
-    && skylines_[DOWN].is_empty ();
+  return out;
 }
 
-Skyline&
-Skyline_pair::operator [] (Direction d)
+bool
+Skyline::is_empty () const
 {
-  return skylines_[d];
+  if (!buildings_.size ())
+    return true;
+  Building b = buildings_.front ();
+  return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
 }
 
-Skyline const&
-Skyline_pair::operator [] (Direction d) const
+void
+Skyline::clear ()
 {
-  return skylines_[d];
+  buildings_.clear ();
+  empty_skyline (&buildings_);
 }
 
 /****************************************************************/
 
-
 IMPLEMENT_SIMPLE_SMOBS (Skyline);
 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
 
-IMPLEMENT_SIMPLE_SMOBS (Skyline_pair);
-IMPLEMENT_TYPE_P (Skyline_pair, "ly:skyline-pair?");
-IMPLEMENT_DEFAULT_EQUAL_P (Skyline_pair);
-
 SCM
-Skyline::mark_smob (SCM)
+Skyline::mark_smob (SCM s)
 {
-  ASSERT_LIVE_IS_ALLOWED();
+  ASSERT_LIVE_IS_ALLOWED (s);
   return SCM_EOL;
 }
 
@@ -679,18 +886,69 @@ Skyline::print_smob (SCM s, SCM port, scm_print_state *)
   return 1;
 }
 
+MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
 SCM
-Skyline_pair::mark_smob (SCM)
+Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
 {
-  return SCM_EOL;
+  LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
+
+  Real horizon_padding = 0;
+  if (horizon_padding_scm != SCM_UNDEFINED)
+    {
+      LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
+      horizon_padding = scm_to_double (horizon_padding_scm);
+    }
+
+  Skyline *skyline = Skyline::unsmob (skyline_scm);
+  Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
+  return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
 }
 
-int
-Skyline_pair::print_smob (SCM s, SCM port, scm_print_state *)
+MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
+SCM
+Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
 {
-  Skyline_pair *r = (Skyline_pair *) SCM_CELL_WORD_1 (s);
-  (void) r;
+  LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
 
-  scm_puts ("#<Skyline-pair>", port);
-  return 1;
+  Real horizon_padding = 0;
+  if (horizon_padding_scm != SCM_UNDEFINED)
+    {
+      LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
+      horizon_padding = scm_to_double (horizon_padding_scm);
+    }
+
+  Skyline *skyline = Skyline::unsmob (skyline_scm);
+  Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
+  return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
+}
+
+MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
+SCM
+Skyline::get_max_height (SCM skyline_scm)
+{
+  return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
+}
+
+MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
+SCM
+Skyline::get_max_height_position (SCM skyline_scm)
+{
+  return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
+}
+
+MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
+SCM
+Skyline::get_height (SCM skyline_scm, SCM x_scm)
+{
+  Real x = robust_scm2double (x_scm, 0.0);
+  return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
+}
+
+LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
+           1, 0, 0, (SCM sky),
+           "Return whether @var{sky} is empty.")
+{
+  Skyline *s = Skyline::unsmob (sky);
+  LY_ASSERT_SMOB (Skyline, sky, 1);
+  return scm_from_bool (s->is_empty ());
 }