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