lily/skyline.cc

   1 /*
   2   This file is part of LilyPond, the GNU music typesetter.
   3
   4   Copyright (C) 2006--2014 Joe Neeman <joeneeman@gmail.com>
   5
   6   LilyPond is free software: you can redistribute it and/or modify
   7   it under the terms of the GNU General Public License as published by
   8   the Free Software Foundation, either version 3 of the License, or
   9   (at your option) any later version.
  10
  11   LilyPond is distributed in the hope that it will be useful,
  12   but WITHOUT ANY WARRANTY; without even the implied warranty of
  13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14   GNU General Public License for more details.
  15
  16   You should have received a copy of the GNU General Public License
  17   along with LilyPond.  If not, see <http://www.gnu.org/licenses/>.
  18 */
  19
  20 #include "skyline.hh"
  21 #include "skyline-pair.hh"
  22 #include <deque>
  23 #include <cstdio>
  24
  25
  26 /* A skyline is a sequence of non-overlapping buildings: something like
  27    this:
  28                    _______
  29                   |       \                                 ________
  30                   |        \                       ________/        \
  31         /\        |          \                    /                  \
  32        /  --------             \                 /                    \
  33       /                          \              /                      \
  34      /                             ------------/                        ----
  35    --
  36    Each building has a starting position, and ending position, a starting
  37    height and an ending height.
  38
  39    The following invariants are observed:
  40     - the start of the first building is at -infinity
  41     - the end of the last building is at infinity
  42     - if a building has infinite length (ie. the first and last buildings),
  43       then its starting height and ending height are equal
  44     - the end of one building is the same as the beginning of the next
  45       building
  46
  47    We also allow skylines to point down (the structure is exactly the same,
  48    but we think of the part above the line as being filled with mass and the
  49    part below as being empty). ::distance finds the minimum distance between
  50    an UP skyline and a DOWN skyline.
  51
  52    Note that we store DOWN skylines upside-down. That is, in order to compare
  53    a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
  54    This means that the merging routine doesn't need to be aware of direction,
  55    but the distance routine does.
  56
  57    From 2007 through 2012, buildings of width less than EPS were discarded,
  58    citing numerical accuracy concerns.  We remember that floating point
  59    comparisons of nearly-equal values can be affected by rounding error.
  60    Also, some target machines use the x87 floating point unit, which provides
  61    extended precision for intermediate results held in registers. On this type
  62    of hardware comparisons such as
  63      double c = 1.0/3.0; boolean compare = (c == 1.0/3.0)
  64    could go either way because the 1.0/3.0 is allowed to be kept
  65    higher precision than the variable 'c'.
  66    Alert to these considerations, we now accept buildings of zero-width.
  67 */
  68
  69 static void
  70 print_buildings (list<Building> const &b)
  71 {
  72   for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
  73     i->print ();
  74 }
  75
  76 void
  77 Skyline::print () const
  78 {
  79   print_buildings (buildings_);
  80 }
  81
  82 void
  83 Skyline::print_points () const
  84 {
  85   vector<Offset> ps (to_points (X_AXIS));
  86
  87   for (vsize i = 0; i < ps.size (); i++)
  88     printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
  89             (i % 2) == 1 ? "\n" : " ");
  90 }
  91
  92 Building::Building (Real start, Real start_height, Real end_height, Real end)
  93 {
  94   if (isinf (start) || isinf (end))
  95     assert (start_height == end_height);
  96
  97   start_ = start;
  98   end_ = end;
  99   precompute (start, start_height, end_height, end);
 100 }
 101
 102 Building::Building (Box const &b, Axis horizon_axis, Direction sky)
 103 {
 104   Real start = b[horizon_axis][LEFT];
 105   Real end = b[horizon_axis][RIGHT];
 106   Real height = sky * b[other_axis (horizon_axis)][sky];
 107
 108   start_ = start;
 109   end_ = end;
 110   precompute (start, height, height, end);
 111 }
 112
 113 void
 114 Building::precompute (Real start, Real start_height, Real end_height, Real end)
 115 {
 116   slope_ = 0.0; /* if they were both infinite, we would get nan, not 0, from the prev line */
 117   if (start_height != end_height)
 118     slope_ = (end_height - start_height) / (end - start);
 119
 120   assert (!isinf (slope_) && !isnan (slope_));
 121
 122   if (isinf (start))
 123     {
 124       assert (start_height == end_height);
 125       y_intercept_ = start_height;
 126     }
 127   else if (fabs(slope_) > 1e6)
 128     // too steep to be stored in slope-intercept form, given round-off error
 129     {
 130       slope_ = 0.0;
 131       y_intercept_ = max(start_height, end_height);
 132     }
 133   else
 134     y_intercept_ = start_height - slope_ * start;
 135 }
 136
 137 inline Real
 138 Building::height (Real x) const
 139 {
 140   return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
 141 }
 142
 143 void
 144 Building::print () const
 145 {
 146   printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
 147 }
 148
 149 inline Real
 150 Building::intersection_x (Building const &other) const
 151 {
 152   Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
 153   return isnan (ret) ? -infinity_f : ret;
 154 }
 155
 156 void
 157 Building::leading_part (Real chop)
 158 {
 159   assert (chop <= end_);
 160   end_ = chop;
 161 }
 162
 163 // Returns a shift s such that (x + s, y) intersects the roof of
 164 // this building.  If no such shift exists, returns infinity_f.
 165 Real
 166 Building::shift_to_intersect (Real x, Real y) const
 167 {
 168   // Solve for s: y = (x + s)*m + b
 169   Real ret = (y - y_intercept_ - slope_ * x) / slope_;
 170
 171   if (ret >= start_ && ret <= end_ && !isinf (ret))
 172     return ret;
 173   return infinity_f;
 174 }
 175
 176 static Real
 177 first_intersection (Building const &b, list<Building> *s, Real start_x)
 178 /* Return the first x >= start_x where skyline s above Building b.
 179  * Removes buildings from s that are concealed by b. */
 180 {
 181   while (!s->empty () && start_x < b.end_)
 182     {
 183       Building c = s->front ();
 184
 185       // conceals and intersection_x involve multiplication and
 186       // division. Avoid that, if we can.
 187       if (c.y_intercept_ == -infinity_f)
 188         {
 189           if (c.end_ > b.end_)
 190             return b.end_;
 191           start_x = c.end_;
 192           s->pop_front ();
 193           continue;
 194         }
 195
 196       if (c.conceals (b, start_x))
 197         return start_x;
 198
 199       Real i = b.intersection_x (c);
 200       if (i > start_x && i <= b.end_ && i <= c.end_)
 201         return i;
 202
 203       start_x = c.end_;
 204       if (b.end_ > c.end_)
 205         s->pop_front ();
 206     }
 207   return b.end_;
 208 }
 209
 210 bool
 211 Building::conceals (Building const &other, Real x) const
 212 {
 213   if (slope_ == other.slope_)
 214     return y_intercept_ > other.y_intercept_;
 215
 216   /* their slopes were not equal, so there is an intersection point */
 217   Real i = intersection_x (other);
 218   return (i <= x && slope_ > other.slope_)
 219          || (i > x && slope_ < other.slope_);
 220 }
 221
 222 // Remove redundant empty buildings from the skyline.
 223 // If there are two adjacent empty buildings, they can be
 224 // turned into one.
 225 void
 226 Skyline::normalize ()
 227 {
 228   bool last_empty = false;
 229   list<Building>::iterator i;
 230
 231   for (i = buildings_.begin (); i != buildings_.end (); i++)
 232     {
 233       if (last_empty && i->y_intercept_ == -infinity_f)
 234         {
 235           list<Building>::iterator last = i;
 236           last--;
 237           last->end_ = i->end_;
 238           buildings_.erase (i);
 239           i = last;
 240         }
 241       last_empty = (i->y_intercept_ == -infinity_f);
 242     }
 243
 244   assert (buildings_.front ().start_ == -infinity_f);
 245   assert (buildings_.back ().end_ == infinity_f);
 246 }
 247
 248 void
 249 Skyline::deholify ()
 250 {
 251   // Since a skyline should always be normalized, we can
 252   // assume that there are never two adjacent empty buildings.
 253   // That is, if center is empty then left and right are not.
 254   list<Building>::iterator left = buildings_.begin ();
 255   list<Building>::iterator center = buildings_.begin ();
 256   list<Building>::iterator right;
 257
 258   for (right = buildings_.begin (); right != buildings_.end (); right++)
 259     {
 260       if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
 261         {
 262           Real p1 = left->height (left->end_);
 263           Real p2 = right->height (right->start_);
 264           *center = Building (center->start_, p1, p2, center->end_);
 265
 266           left = center;
 267           center = right;
 268         }
 269     }
 270 }
 271
 272 void
 273 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
 274                                  list<Building> *const result) const
 275 {
 276   if (s1->empty () || s2->empty ())
 277     {
 278       programming_error ("tried to merge an empty skyline");
 279       return;
 280     }
 281
 282   Real x = -infinity_f;
 283   Real last_end = -infinity_f;
 284   while (!s1->empty ())
 285     {
 286       if (s2->front ().conceals (s1->front (), x))
 287         swap (s1, s2);
 288
 289       Building b = s1->front ();
 290       Building c = s2->front ();
 291
 292       // Optimization: if the other skyline is empty at this point,
 293       // we can avoid testing some intersections. Just grab as many
 294       // buildings from s1 as we can, and shove them onto the output.
 295       if (c.y_intercept_ == -infinity_f
 296           && c.end_ >= b.end_)
 297         {
 298           list<Building>::iterator i = s1->begin ();
 299           i++;
 300           while (i != s1->end () && i->end_ <= c.end_)
 301             i++;
 302
 303           s1->front ().start_ = x;
 304           result->splice (result->end (), *s1, s1->begin (), i);
 305           x = result->back ().end_;
 306           last_end = x;
 307           continue;
 308         }
 309       // first_intersection() removes buildings from s2 if b hides them
 310       Real end = first_intersection (b, s2, x);
 311       if (s2->empty ())
 312         {
 313           b.start_ = last_end;
 314           result->push_back (b);
 315           break;
 316         }
 317
 318       // Should be (end > x), during ver2.19.  end == x happens fairly often,
 319       // and we do not need to keep vertical segments within a skyline.
 320       if (end >= x)
 321         {
 322           b.leading_part (end);
 323           b.start_ = last_end;
 324           last_end = b.end_;
 325           result->push_back (b);
 326         }
 327
 328       if (end >= s1->front ().end_)
 329         s1->pop_front ();
 330       // Should add during ver2.19 (to avoid an endless loop
 331       // when  merging identical skylines with a vertical segment)
 332       // if (end >= s2->front().end_) s2->pop_front();
 333
 334       x = end;
 335     }
 336 }
 337
 338 static void
 339 empty_skyline (list<Building> *const ret)
 340 {
 341   ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
 342 }
 343
 344 /*
 345   Given Building 'b', build a skyline containing only that building.
 346 */
 347 static void
 348 single_skyline (Building b, list<Building> *const ret)
 349 {
 350   assert (b.end_ >= b.start_);
 351
 352   if (b.start_ != -infinity_f)
 353     ret->push_back (Building (-infinity_f, -infinity_f,
 354                               -infinity_f, b.start_));
 355   ret->push_back (b);
 356   if (b.end_ != infinity_f)
 357     ret->push_back (Building (b.end_, -infinity_f,
 358                               -infinity_f, infinity_f));
 359 }
 360
 361 /* remove a non-overlapping set of boxes from BOXES and build a skyline
 362    out of them */
 363 static list<Building>
 364 non_overlapping_skyline (list<Building> *const buildings)
 365 {
 366   list<Building> result;
 367   Real last_end = -infinity_f;
 368   Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
 369   list<Building>::iterator i = buildings->begin ();
 370   while (i != buildings->end ())
 371     {
 372       Real x1 = i->start_;
 373       Real y1 = i->height (i->start_);
 374       Real x2 = i->end_;
 375       Real y2 = i->height (i->end_);
 376
 377       // Drop buildings that will obviously have no effect.
 378       if (last_building.height (x1) >= y1
 379           && last_building.end_ >= x2
 380           && last_building.height (x2) >= y2)
 381         {
 382           list<Building>::iterator j = i++;
 383           buildings->erase (j);
 384           continue;
 385         }
 386
 387       if (x1 < last_end)
 388         {
 389           i++;
 390           continue;
 391         }
 392
 393       // Insert empty Buildings into any gaps. (TODO: is this needed? -KOH)
 394       if (x1 > last_end)
 395         result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
 396
 397       result.push_back (*i);
 398       last_building = *i;
 399       last_end = i->end_;
 400
 401       list<Building>::iterator j = i++;
 402       buildings->erase (j);
 403     }
 404
 405   if (last_end < infinity_f)
 406     result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
 407   return result;
 408 }
 409
 410 class LessThanBuilding
 411 {
 412 public:
 413   bool operator () (Building const &b1, Building const &b2)
 414   {
 415     return b1.start_ < b2.start_
 416            || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
 417   }
 418 };
 419
 420 /**
 421    BUILDINGS is a list of buildings, but they could be overlapping
 422    and in any order.  The returned list of buildings is ordered and non-overlapping.
 423 */
 424 list<Building>
 425 Skyline::internal_build_skyline (list<Building> *buildings) const
 426 {
 427   vsize size = buildings->size ();
 428
 429   if (size == 0)
 430     {
 431       list<Building> result;
 432       empty_skyline (&result);
 433       return result;
 434     }
 435   else if (size == 1)
 436     {
 437       list<Building> result;
 438       single_skyline (buildings->front (), &result);
 439       return result;
 440     }
 441
 442   deque<list<Building> > partials;
 443   buildings->sort (LessThanBuilding ());
 444   while (!buildings->empty ())
 445     partials.push_back (non_overlapping_skyline (buildings));
 446
 447   /* we'd like to say while (partials->size () > 1) but that's O (n).
 448      Instead, we exit in the middle of the loop */
 449   while (!partials.empty ())
 450     {
 451       list<Building> merged;
 452       list<Building> one = partials.front ();
 453       partials.pop_front ();
 454       if (partials.empty ())
 455         return one;
 456
 457       list<Building> two = partials.front ();
 458       partials.pop_front ();
 459       internal_merge_skyline (&one, &two, &merged);
 460       partials.push_back (merged);
 461     }
 462   assert (0);
 463   return list<Building> ();
 464 }
 465
 466 Skyline::Skyline ()
 467 {
 468   sky_ = UP;
 469   empty_skyline (&buildings_);
 470 }
 471
 472 Skyline::Skyline (Direction sky)
 473 {
 474   sky_ = sky;
 475   empty_skyline (&buildings_);
 476 }
 477
 478 /*
 479   Build skyline from a set of boxes.
 480
 481   Boxes should be non-empty on both axes.  Otherwise, they will be ignored
 482  */
 483 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
 484 {
 485   list<Building> buildings;
 486   sky_ = sky;
 487
 488   for (vsize i = 0; i < boxes.size (); i++)
 489     if (!boxes[i].is_empty (X_AXIS)
 490         && !boxes[i].is_empty (Y_AXIS))
 491       buildings.push_front (Building (boxes[i], horizon_axis, sky));
 492
 493   buildings_ = internal_build_skyline (&buildings);
 494   normalize ();
 495 }
 496
 497 /*
 498   build skyline from a set of line segments.
 499
 500   Segments can be articulated from left to right or right to left.
 501   In the case of the latter, they will be stored internally as left to right.
 502  */
 503 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
 504 {
 505   list<Building> buildings;
 506   sky_ = sky;
 507
 508   for (vsize i = 0; i < segments.size (); i++)
 509     {
 510       Drul_array<Offset> const &seg = segments[i];
 511       Offset left = seg[LEFT];
 512       Offset right = seg[RIGHT];
 513       if (left[horizon_axis] > right[horizon_axis])
 514         swap (left, right);
 515
 516       Real x1 = left[horizon_axis];
 517       Real x2 = right[horizon_axis];
 518       Real y1 = left[other_axis (horizon_axis)] * sky;
 519       Real y2 = right[other_axis (horizon_axis)] * sky;
 520
 521       if (x1 <= x2)
 522         buildings.push_back (Building (x1, y1, y2, x2));
 523     }
 524
 525   buildings_ = internal_build_skyline (&buildings);
 526   normalize ();
 527 }
 528
 529 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
 530 {
 531   sky_ = sky;
 532
 533   deque<Skyline> partials;
 534   for (vsize i = 0; i < skypairs.size (); i++)
 535     partials.push_back (Skyline ((skypairs[i])[sky]));
 536
 537   while (partials.size () > 1)
 538     {
 539       Skyline one = partials.front ();
 540       partials.pop_front ();
 541       Skyline two = partials.front ();
 542       partials.pop_front ();
 543
 544       one.merge (two);
 545       partials.push_back (one);
 546     }
 547
 548   if (partials.size ())
 549     buildings_.swap (partials.front ().buildings_);
 550   else
 551     buildings_.clear ();
 552 }
 553
 554 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
 555 {
 556   sky_ = sky;
 557   if (!b.is_empty (X_AXIS) && !b.is_empty (Y_AXIS))
 558     {
 559       Building front (b, horizon_axis, sky);
 560       single_skyline (front, &buildings_);
 561       normalize ();
 562     }
 563 }
 564
 565 void
 566 Skyline::merge (Skyline const &other)
 567 {
 568   assert (sky_ == other.sky_);
 569
 570   if (other.is_empty ())
 571     return;
 572
 573   if (is_empty ())
 574     {
 575       buildings_ = other.buildings_;
 576       return;
 577     }
 578
 579   list<Building> other_bld (other.buildings_);
 580   list<Building> my_bld;
 581   my_bld.splice (my_bld.begin (), buildings_);
 582   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 583   normalize ();
 584 }
 585
 586 void
 587 Skyline::insert (Box const &b, Axis a)
 588 {
 589   list<Building> other_bld;
 590   list<Building> my_bld;
 591
 592   if (isnan (b[other_axis (a)][LEFT])
 593       || isnan (b[other_axis (a)][RIGHT]))
 594     {
 595       programming_error ("insane box for skyline");
 596       return;
 597     }
 598
 599   /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
 600   if (b.is_empty (X_AXIS) || b.is_empty (Y_AXIS))
 601     return;
 602
 603   my_bld.splice (my_bld.begin (), buildings_);
 604   single_skyline (Building (b, a, sky_), &other_bld);
 605   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 606   normalize ();
 607 }
 608
 609 void
 610 Skyline::raise (Real r)
 611 {
 612   list<Building>::iterator end = buildings_.end ();
 613   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 614     i->y_intercept_ += sky_ * r;
 615 }
 616
 617 void
 618 Skyline::shift (Real s)
 619 {
 620   list<Building>::iterator end = buildings_.end ();
 621   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 622     {
 623       i->start_ += s;
 624       i->end_ += s;
 625       i->y_intercept_ -= s * i->slope_;
 626     }
 627 }
 628
 629 Real
 630 Skyline::distance (Skyline const &other, Real horizon_padding) const
 631 {
 632   Real dummy;
 633   return internal_distance (other, horizon_padding, &dummy);
 634 }
 635
 636 Real
 637 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
 638 {
 639   Real touch;
 640   internal_distance (other, horizon_padding, &touch);
 641   return touch;
 642 }
 643
 644 Real
 645 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
 646 {
 647   if (horizon_padding == 0.0)
 648     return internal_distance (other, touch_point);
 649
 650   // Note that it is not necessary to build a padded version of other,
 651   // because the same effect can be achieved just by doubling horizon_padding.
 652   Skyline padded_this = padded (horizon_padding);
 653   return padded_this.internal_distance (other, touch_point);
 654 }
 655
 656 Skyline
 657 Skyline::padded (Real horizon_padding) const
 658 {
 659   if (horizon_padding < 0.0)
 660     warning ("Cannot have negative horizon padding.  Junking.");
 661
 662   if (horizon_padding <= 0.0)
 663     return *this;
 664
 665   list<Building> pad_buildings;
 666   for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
 667     {
 668       if (i->start_ > -infinity_f)
 669         {
 670           Real height = i->height (i->start_);
 671           if (height > -infinity_f)
 672             {
 673               // Add the sloped building that pads the left side of the current building.
 674               Real start = i->start_ - 2 * horizon_padding;
 675               Real end = i->start_ - horizon_padding;
 676               pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
 677
 678               // Add the flat building that pads the left side of the current building.
 679               start = i->start_ - horizon_padding;
 680               end = i->start_;
 681               pad_buildings.push_back (Building (start, height, height, end));
 682             }
 683         }
 684
 685       if (i->end_ < infinity_f)
 686         {
 687           Real height = i->height (i->end_);
 688           if (height > -infinity_f)
 689             {
 690               // Add the flat building that pads the right side of the current building.
 691               Real start = i->end_;
 692               Real end = start + horizon_padding;
 693               pad_buildings.push_back (Building (start, height, height, end));
 694
 695               // Add the sloped building that pads the right side of the current building.
 696               start = end;
 697               end += horizon_padding;
 698               pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
 699             }
 700         }
 701     }
 702
 703   // The buildings may be overlapping, so resolve that.
 704   list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
 705
 706   // Merge the padding with the original, to make a new skyline.
 707   Skyline padded (sky_);
 708   list<Building> my_buildings = buildings_;
 709   padded.buildings_.clear ();
 710   internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
 711   padded.normalize ();
 712
 713   return padded;
 714 }
 715
 716 Real
 717 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
 718 {
 719   assert (sky_ == -other.sky_);
 720
 721   list<Building>::const_iterator i = buildings_.begin ();
 722   list<Building>::const_iterator j = other.buildings_.begin ();
 723
 724   Real dist = -infinity_f;
 725   Real start = -infinity_f;
 726   Real touch = -infinity_f;
 727   while (i != buildings_.end () && j != other.buildings_.end ())
 728     {
 729       Real end = min (i->end_, j->end_);
 730       Real start_dist = i->height (start) + j->height (start);
 731       Real end_dist = i->height (end) + j->height (end);
 732       dist = max (dist, max (start_dist, end_dist));
 733
 734       if (end_dist == dist)
 735         touch = end;
 736       else if (start_dist == dist)
 737         touch = start;
 738
 739       if (i->end_ <= j->end_)
 740         i++;
 741       else
 742         j++;
 743       start = end;
 744     }
 745
 746   *touch_point = touch;
 747   return dist;
 748 }
 749
 750 Real
 751 Skyline::height (Real airplane) const
 752 {
 753   assert (!isinf (airplane));
 754
 755   list<Building>::const_iterator i;
 756   for (i = buildings_.begin (); i != buildings_.end (); i++)
 757     {
 758       if (i->end_ >= airplane)
 759         return sky_ * i->height (airplane);
 760     }
 761
 762   assert (0);
 763   return 0;
 764 }
 765
 766 Real
 767 Skyline::max_height () const
 768 {
 769   Real ret = -infinity_f;
 770
 771   list<Building>::const_iterator i;
 772   for (i = buildings_.begin (); i != buildings_.end (); ++i)
 773     {
 774       ret = max (ret, i->height (i->start_));
 775       ret = max (ret, i->height (i->end_));
 776     }
 777
 778   return sky_ * ret;
 779 }
 780
 781 Direction
 782 Skyline::direction () const
 783 {
 784   return sky_;
 785 }
 786
 787 Real
 788 Skyline::left () const
 789 {
 790   for (list<Building>::const_iterator i (buildings_.begin ());
 791        i != buildings_.end (); i++)
 792     if (i->y_intercept_ > -infinity_f)
 793       return i->start_;
 794
 795   return infinity_f;
 796 }
 797
 798 Real
 799 Skyline::right () const
 800 {
 801   for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
 802        i != buildings_.rend (); ++i)
 803     if (i->y_intercept_ > -infinity_f)
 804       return i->end_;
 805
 806   return -infinity_f;
 807 }
 808
 809 Real
 810 Skyline::max_height_position () const
 811 {
 812   Skyline s (-sky_);
 813   s.set_minimum_height (0);
 814   return touching_point (s);
 815 }
 816
 817 void
 818 Skyline::set_minimum_height (Real h)
 819 {
 820   Skyline s (sky_);
 821   s.buildings_.front ().y_intercept_ = h * sky_;
 822   merge (s);
 823 }
 824
 825 vector<Offset>
 826 Skyline::to_points (Axis horizon_axis) const
 827 {
 828   vector<Offset> out;
 829
 830   Real start = -infinity_f;
 831   for (list<Building>::const_iterator i (buildings_.begin ());
 832        i != buildings_.end (); i++)
 833     {
 834       out.push_back (Offset (start, sky_ * i->height (start)));
 835       out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
 836       start = i->end_;
 837     }
 838
 839   if (horizon_axis == Y_AXIS)
 840     for (vsize i = 0; i < out.size (); i++)
 841       out[i] = out[i].swapped ();
 842
 843   return out;
 844 }
 845
 846 bool
 847 Skyline::is_empty () const
 848 {
 849   if (!buildings_.size ())
 850     return true;
 851   Building b = buildings_.front ();
 852   return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
 853 }
 854
 855 void
 856 Skyline::clear ()
 857 {
 858   buildings_.clear ();
 859   empty_skyline (&buildings_);
 860 }
 861
 862 /****************************************************************/
 863
 864 const char Skyline::type_p_name_[] = "ly:skyline?";
 865
 866 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
 867 SCM
 868 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
 869 {
 870   LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
 871
 872   Real horizon_padding = 0;
 873   if (horizon_padding_scm != SCM_UNDEFINED)
 874     {
 875       LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
 876       horizon_padding = scm_to_double (horizon_padding_scm);
 877     }
 878
 879   Skyline *skyline = Skyline::unsmob (skyline_scm);
 880   Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
 881   return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
 882 }
 883
 884 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
 885 SCM
 886 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
 887 {
 888   LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
 889
 890   Real horizon_padding = 0;
 891   if (horizon_padding_scm != SCM_UNDEFINED)
 892     {
 893       LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
 894       horizon_padding = scm_to_double (horizon_padding_scm);
 895     }
 896
 897   Skyline *skyline = Skyline::unsmob (skyline_scm);
 898   Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
 899   return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
 900 }
 901
 902 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
 903 SCM
 904 Skyline::get_max_height (SCM skyline_scm)
 905 {
 906   return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
 907 }
 908
 909 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
 910 SCM
 911 Skyline::get_max_height_position (SCM skyline_scm)
 912 {
 913   return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
 914 }
 915
 916 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
 917 SCM
 918 Skyline::get_height (SCM skyline_scm, SCM x_scm)
 919 {
 920   Real x = robust_scm2double (x_scm, 0.0);
 921   return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
 922 }
 923
 924 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
 925            1, 0, 0, (SCM sky),
 926            "Return whether @var{sky} is empty.")
 927 {
 928   Skyline *s = Skyline::unsmob (sky);
 929   LY_ASSERT_SMOB (Skyline, sky, 1);
 930   return scm_from_bool (s->is_empty ());
 931 }