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