]> git.donarmstrong.com Git - lilypond.git/blobdiff - lily/skyline.cc
Merge remote-tracking branch 'origin/translation'
[lilypond.git] / lily / skyline.cc
index c6ec8c4229fdb4eb4692f96d096d2309d83a63af..33d2823d47e0d752934558f75ea74e90e85ead41 100644 (file)
@@ -1,23 +1,38 @@
-/* skyline.cc -- implement the Skyline class
+/*
+  This file is part of LilyPond, the GNU music typesetter.
+
+  Copyright (C) 2006--2014 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 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 skyline is a sequence of non-overlapping buildings: something like
    this:
                    _______
-                  /       \                                ________
-                 /         \                      ________/        \
-        /\                \                    /                  \
-       /  -----             \                 /                    \
-      /                         \              /                      \
-     /                            ------------/                        ----
+                  |       \                                 ________
+                  |        \                       ________/        \
+        /\        |          \                    /                  \
+       /  --------             \                 /                    \
+      /                          \              /                      \
+     /                             ------------/                        ----
    --
    Each building has a starting position, and ending position, a starting
    height and an ending height.
    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.
-*/
 
-#define EPS 1e-10
+   From 2007 through 2012, buildings of width less than EPS were discarded,
+   citing numerical accuracy concerns.  We remember that floating point
+   comparisons of nearly-equal values can be affected by rounding error.
+   Also, some target machines use the x87 floating point unit, which provides
+   extended precision for intermediate results held in registers. On this type
+   of hardware comparisons such as
+     double c = 1.0/3.0; boolean compare = (c == 1.0/3.0)
+   could go either way because the 1.0/3.0 is allowed to be kept
+   higher precision than the variable 'c'.
+   Alert to these considerations, we now accept buildings of zero-width.
+*/
 
-static inline bool
-equal (Real x, Real y)
+static void
+print_buildings (list<Building> const &b)
 {
-  return abs (x - y) < EPS || (isinf (x) && isinf (y) && ((x > 0) == (y > 0)));
+  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
-Skyline::is_legal_skyline () const
+void
+Skyline::print_points () const
 {
-  list<Building>::const_iterator i;
-  Real last_x = -infinity_f;
-  Real last_h = -infinity_f;
-  for (i = buildings_.begin (); i != buildings_.end (); i++)
-    {
-      if (i->iv_[LEFT] != last_x)
-       return false;
-      if (i != buildings_.begin () && !equal (i->start_height_, last_h))
-       return false;
-      last_x = i->iv_[RIGHT];
-      last_h = i->end_height_;
-    }
-  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, Real max_slope)
-  : iv_ (start, end)
+Building::Building (Real start, Real start_height, Real end_height, Real end)
 {
-  start_height_ = start_height;
-  end_height_ = end_height;
+  if (isinf (start) || isinf (end))
+    assert (start_height == end_height);
 
-  if (isinf (start))
-    assert (isinf (start_height) || start_height == end_height);
-  if (isinf (end))
-    assert (isinf (end_height) || start_height == end_height);
+  start_ = start;
+  end_ = end;
+  precompute (start, start_height, end_height, end);
+}
 
-  precompute (max_slope);
+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];
+
+  start_ = start;
+  end_ = end;
+  precompute (start, height, height, end);
 }
 
 void
-Building::precompute (Real max_slope)
-{
-  slope_ = (end_height_ - start_height_) / (iv_.length());
-  if (start_height_ == end_height_)
-    slope_ = 0;
-  if (isinf (slope_) || isnan (slope_))
-    slope_ = max_slope * (start_height_ < end_height_ ? 1 : -1);
-
-#if 0
-  /*
-    this check is sensitive to roundoff errors when converting to/from
-    sequences of points.
-   */
-  assert (abs (slope_) <= max_slope + EPS);
-#endif
-  
-  if (isinf (iv_[START]))
+Building::precompute (Real start, Real start_height, Real end_height, Real end)
+{
+  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 (start))
     {
-      if (isinf (iv_[STOP]))
-       zero_height_ = start_height_;
-      else
-       zero_height_ = end_height_ - slope_ * iv_[STOP];
+      assert (start_height == end_height);
+      y_intercept_ = start_height;
+    }
+  else if (fabs(slope_) > 1e6)
+    // too steep to be stored in slope-intercept form, given round-off error
+    {
+      slope_ = 0.0;
+      y_intercept_ = max(start_height, end_height);
     }
   else
-    zero_height_ = start_height_ - slope_ * iv_[START];
+    y_intercept_ = start_height - slope_ * start;
 }
 
-Real 
+inline Real
 Building::height (Real x) const
 {
-  if (isinf (x))
-    return (x > 0) ? end_height_ : start_height_;
-  return slope_*x + zero_height_;
+  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],
-         start_height_, end_height_);
+  printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
 }
 
-Real
-Building::intersection (Building const &other) const
+inline Real
+Building::intersection_x (Building const &other) const
 {
-  return (zero_height_ - other.zero_height_) / (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, Real h)
+Building::leading_part (Real chop)
 {
-  assert (chop > iv_[LEFT] && chop <= iv_[RIGHT] && !equal (chop, iv_[LEFT]));
-  assert (equal (h, height (chop)));
-  iv_[RIGHT] = chop;
-  end_height_ = h;
+  assert (chop <= end_);
+  end_ = chop;
 }
 
-static void
-skyline_trailing_part (list<Building> *sky, Real x)
+// 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
 {
-  if (equal (x, sky->front ().iv_[RIGHT]))
-    sky->pop_front ();
-  else
-    assert (x < sky->front ().iv_[RIGHT]);
+  // Solve for s: y = (x + s)*m + b
+  Real ret = (y - y_intercept_ - slope_ * x) / slope_;
 
-  if (!sky->empty ())
+  if (ret >= start_ && ret <= end_ && !isinf (ret))
+    return ret;
+  return infinity_f;
+}
+
+static Real
+first_intersection (Building const &b, list<Building> *s, Real start_x)
+/* Return the first x >= start_x where skyline s above Building b.
+ * Removes buildings from s that are concealed by b. */
+{
+  while (!s->empty () && start_x < b.end_)
     {
-      sky->front ().iv_[LEFT] = x;
-      sky->front ().start_height_ = sky->front ().height (x);
+      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 b.end_;
 }
 
 bool
-Building::obstructs (Building const &other) const
+Building::conceals (Building const &other, Real x) const
+{
+  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_);
+}
+
+// Remove redundant empty buildings from the skyline.
+// If there are two adjacent empty buildings, they can be
+// turned into one.
+void
+Skyline::normalize ()
 {
-  if (equal (intersection (other), iv_[LEFT]) || equal (start_height_, other.start_height_))
-    return slope_ > other.slope_ || (slope_ == other.slope_ && zero_height_ > other.zero_height_);
-  return start_height_ > other.start_height_;
+  bool last_empty = false;
+  list<Building>::iterator i;
+
+  for (i = buildings_.begin (); i != buildings_.end (); i++)
+    {
+      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);
+    }
+
+  assert (buildings_.front ().start_ == -infinity_f);
+  assert (buildings_.back ().end_ == infinity_f);
+}
+
+void
+Skyline::deholify ()
+{
+  // 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 ().obstructs (s1->front ()))
-       swap (s1, s2);
+      if (s2->front ().conceals (s1->front (), x))
+        swap (s1, s2);
 
       Building b = s1->front ();
-      while (s2->front ().iv_[RIGHT] < b.iv_[RIGHT]
-            && s2->front ().end_height_ <= b.height (s2->front ().iv_[RIGHT]) + EPS)
-       s2->pop_front ();
-
-      /* the front of s2 either intersects with b or it ends after b */
-      Real end = infinity_f;
-      Real s2_end_height = s2->front ().end_height_;
-      Real s1_end_height = b.height (s2->front ().iv_[RIGHT]);
-      if (s2_end_height > s1_end_height + EPS)
-       end = b.intersection (s2->front ());
-      end = min (end, b.iv_[RIGHT]);
-      Real height = b.height (end);
-
-      b.leading_part (end, height);
-      result->push_front (b);
-
-      skyline_trailing_part (s1, end);
-      if (!s1->empty ())
-       s1->front ().start_height_ = height;
-      skyline_trailing_part (s2, end);
+      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;
+        }
+      // first_intersection() removes buildings from s2 if b hides them
+      Real end = first_intersection (b, s2, x);
+      if (s2->empty ())
+        {
+          b.start_ = last_end;
+          result->push_back (b);
+          break;
+        }
+
+      // Should be (end > x), during ver2.19.  end == x happens fairly often,
+      // and we do not need to keep vertical segments within a skyline.
+      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 ();
+      // Should add during ver2.19 (to avoid an endless loop
+      // when  merging identical skylines with a vertical segment)
+      // if (end >= s2->front().end_) s2->pop_front();
+
+      x = end;
     }
-  result->reverse ();
 }
 
 static void
 empty_skyline (list<Building> *const ret)
 {
-  ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f, 0));
+  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 const &b, list<Building> *const ret, Real max_slope)
+single_skyline (Building b, list<Building> *const ret)
 {
-  if (!isinf (b.iv_[RIGHT]))
-    ret->push_front (Building (b.iv_[RIGHT], b.end_height_,
-                              -infinity_f, infinity_f, max_slope));
-  ret->push_front (b);
-  if (!isinf (b.iv_[LEFT]))
-    ret->push_front (Building (-infinity_f, -infinity_f,
-                              b.start_height_, b.iv_[LEFT], max_slope));
+  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));
 }
 
-void
-Skyline::internal_build_skyline (list<Building> *buildings, list<Building> *const result)
+/* 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 ())
+    {
+      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;
+        }
+
+      // Insert empty Buildings into any gaps. (TODO: is this needed? -KOH)
+      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));
+  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) const
 {
   vsize size = buildings->size ();
 
   if (size == 0)
     {
-      empty_skyline (result);
-      return;
+      list<Building> result;
+      empty_skyline (&result);
+      return result;
     }
   else if (size == 1)
     {
-      single_skyline (buildings->front (), result, max_slope_);
-      return;
+      list<Building> result;
+      single_skyline (buildings->front (), &result);
+      return result;
     }
 
-  list<Building> right_half;
-  list<Building>::iterator i = buildings->begin ();
-
-  for (vsize s = 0; s < size/2; s++)
-    i++;
-  right_half.splice (right_half.end (), *buildings, i, buildings->end ());
+  deque<list<Building> > partials;
+  buildings->sort (LessThanBuilding ());
+  while (!buildings->empty ())
+    partials.push_back (non_overlapping_skyline (buildings));
 
-  list<Building> right;
-  list<Building> left;
-  internal_build_skyline (&right_half, &right);
-  internal_build_skyline (buildings, &left);
-  internal_merge_skyline (&right, &left, result);
+  /* 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 ())
+    {
+      list<Building> merged;
+      list<Building> one = partials.front ();
+      partials.pop_front ();
+      if (partials.empty ())
+        return one;
+
+      list<Building> two = partials.front ();
+      partials.pop_front ();
+      internal_merge_skyline (&one, &two, &merged);
+      partials.push_back (merged);
+    }
+  assert (0);
+  return list<Building> ();
 }
 
 Skyline::Skyline ()
 {
-  max_slope_ = 2;
   sky_ = UP;
   empty_skyline (&buildings_);
-
-  
-}
-
-Skyline::Skyline (Skyline const &src)
-{
-  max_slope_ = src.max_slope_;
-  sky_ = src.sky_;
-  for (list<Building>::const_iterator i = src.buildings_.begin ();
-       i != src.buildings_.end (); i++)
-    {
-      buildings_.push_back (Building ((*i)));
-    }
 }
 
 Skyline::Skyline (Direction sky)
 {
-  max_slope_ = 2;
   sky_ = sky;
   empty_skyline (&buildings_);
 }
 
 /*
-  build skyline from a set of boxes.
+  Build skyline from a set of boxes.
 
-  Boxes should have fatness in the horizon_axis, otherwise they are ignored.
+  Boxes should be non-empty on both axes.  Otherwise, they will be ignored
  */
 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
 {
-  list<Building> bldgs;
+  list<Building> buildings;
   sky_ = sky;
-  max_slope_ = 2;
 
   for (vsize i = 0; i < boxes.size (); i++)
+    if (!boxes[i].is_empty (X_AXIS)
+        && !boxes[i].is_empty (Y_AXIS))
+      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++)
     {
-      Interval iv = boxes[i][horizon_axis];
-      Real height = sky * boxes[i][other_axis (horizon_axis)][sky];
-      if (!iv.is_empty () && !isinf (height) && !equal (iv[LEFT], iv[RIGHT]))
-       bldgs.push_front (Building (iv[LEFT], height, height, iv[RIGHT],
-                                   max_slope_));
+      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));
     }
-  
-  internal_build_skyline (&bldgs, &buildings_);
-  assert (is_legal_skyline ());
+
+  buildings_ = internal_build_skyline (&buildings);
+  normalize ();
 }
 
-Skyline::Skyline (vector<Offset> const &points, Real max_slope, Direction sky)
+Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
 {
   sky_ = sky;
-  max_slope_ = max_slope;
 
-  for (vsize i = 1; i < points.size (); i++)
-    {
-      buildings_.push_back (Building (points[i-1][X_AXIS], sky * points[i-1][Y_AXIS],
-                                     sky * points[i][Y_AXIS],
+  deque<Skyline> partials;
+  for (vsize i = 0; i < skypairs.size (); i++)
+    partials.push_back (Skyline ((skypairs[i])[sky]));
 
-                                     points[i][X_AXIS], 
-                                     max_slope));
+  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);
     }
 
-  assert (is_legal_skyline ());
+  if (partials.size ())
+    buildings_.swap (partials.front ().buildings_);
+  else
+    buildings_.clear ();
+}
+
+Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
+{
+  sky_ = sky;
+  if (!b.is_empty (X_AXIS) && !b.is_empty (Y_AXIS))
+    {
+      Building front (b, horizon_axis, sky);
+      single_skyline (front, &buildings_);
+      normalize ();
+    }
 }
 
 void
@@ -331,11 +568,20 @@ 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 ());
+  normalize ();
 }
 
 void
@@ -343,13 +589,22 @@ Skyline::insert (Box const &b, Axis a)
 {
   list<Building> other_bld;
   list<Building> my_bld;
-  Interval iv = b[a];
-  Real height = sky_ * b[other_axis (a)][sky_];
+
+  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 (X_AXIS) || b.is_empty (Y_AXIS))
+    return;
 
   my_bld.splice (my_bld.begin (), buildings_);
-  single_skyline (Building (iv[LEFT], height, height, iv[RIGHT], max_slope_), &other_bld, max_slope_);
+  single_skyline (Building (b, a, sky_), &other_bld);
   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
-  assert (is_legal_skyline ());
+  normalize ();
 }
 
 void
@@ -357,32 +612,139 @@ Skyline::raise (Real r)
 {
   list<Building>::iterator end = buildings_.end ();
   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
+    i->y_intercept_ += sky_ * r;
+}
+
+void
+Skyline::shift (Real s)
+{
+  list<Building>::iterator end = buildings_.end ();
+  for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
+    {
+      i->start_ += s;
+      i->end_ += s;
+      i->y_intercept_ -= s * i->slope_;
+    }
+}
+
+Real
+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)
     {
-      i->start_height_ += sky_ * r;
-      i->end_height_ += sky_ * r;
-      i->zero_height_ += sky_ * r;
+      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));
+            }
+        }
     }
-  assert (is_legal_skyline ());
+
+  // 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::distance (Skyline const &other) const
+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;
 }
 
@@ -394,8 +756,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);
@@ -404,52 +766,110 @@ Skyline::height (Real airplane) const
 
 Real
 Skyline::max_height () const
+{
+  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;
+}
+
+Direction
+Skyline::direction () const
+{
+  return sky_;
+}
+
+Real
+Skyline::left () const
+{
+  for (list<Building>::const_iterator i (buildings_.begin ());
+       i != buildings_.end (); i++)
+    if (i->y_intercept_ > -infinity_f)
+      return i->start_;
+
+  return infinity_f;
+}
+
+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_;
+
+  return -infinity_f;
+}
+
+Real
+Skyline::max_height_position () const
 {
   Skyline s (-sky_);
   s.set_minimum_height (0);
-  return sky_ * distance (s);
+  return touching_point (s);
 }
 
 void
 Skyline::set_minimum_height (Real h)
 {
   Skyline s (sky_);
-  s.buildings_.front ().start_height_ = h * sky_;
-  s.buildings_.front ().end_height_ = h * sky_;
-  s.buildings_.front ().zero_height_ = h * sky_;
+  s.buildings_.front ().y_intercept_ = h * sky_;
   merge (s);
 }
 
-
 vector<Offset>
-Skyline::to_points () const
+Skyline::to_points (Axis horizon_axis) const
 {
   vector<Offset> out;
 
-  bool first = true;
+  Real start = -infinity_f;
   for (list<Building>::const_iterator i (buildings_.begin ());
        i != buildings_.end (); i++)
     {
-      if (first)
-       out.push_back (Offset ((*i).iv_[LEFT], sky_ * (*i).start_height_));
-
-      first = false;
-      out.push_back (Offset ((*i).iv_[RIGHT], sky_ * (*i).end_height_));
+      out.push_back (Offset (start, sky_ * i->height (start)));
+      out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
+      start = i->end_;
     }
 
+  if (horizon_axis == Y_AXIS)
+    for (vsize i = 0; i < out.size (); i++)
+      out[i] = out[i].swapped ();
+
   return out;
 }
 
-/****************************************************************/
+bool
+Skyline::is_empty () const
+{
+  if (!buildings_.size ())
+    return true;
+  Building b = buildings_.front ();
+  return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
+}
 
+void
+Skyline::clear ()
+{
+  buildings_.clear ();
+  empty_skyline (&buildings_);
+}
+
+/****************************************************************/
 
 IMPLEMENT_SIMPLE_SMOBS (Skyline);
 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
 
 SCM
-Skyline::mark_smob (SCM)
+Skyline::mark_smob (SCM s)
 {
+  ASSERT_LIVE_IS_ALLOWED (s);
   return SCM_EOL;
 }
 
@@ -463,3 +883,70 @@ 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::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
+{
+  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));
+}
+
+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)
+{
+  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->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 ());
+}