lily/skyline.cc

   1 /*
   2   This file is part of LilyPond, the GNU music typesetter.
   3
   4   Copyright (C) 2006--2012 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
  59 /* If we start including very thin buildings, numerical accuracy errors can
  60    arise. Therefore, we ignore all buildings that are less than epsilon wide. */
  61 #define EPS 1e-5
  62
  63 static void
  64 print_buildings (list<Building> const &b)
  65 {
  66   for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
  67     i->print ();
  68 }
  69
  70 void
  71 Skyline::print () const
  72 {
  73   print_buildings (buildings_);
  74 }
  75
  76 void
  77 Skyline::print_points () const
  78 {
  79   vector<Offset> ps (to_points (X_AXIS));
  80
  81   for (vsize i = 0; i < ps.size (); i++)
  82     printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
  83             (i % 2) == 1 ? "\n" : " ");
  84 }
  85
  86 Building::Building (Real start, Real start_height, Real end_height, Real end)
  87 {
  88   if (isinf (start) || isinf (end))
  89     assert (start_height == end_height);
  90
  91   start_ = start;
  92   end_ = end;
  93   precompute (start, start_height, end_height, end);
  94 }
  95
  96 Building::Building (Box const &b, Axis horizon_axis, Direction sky)
  97 {
  98   Real start = b[horizon_axis][LEFT];
  99   Real end = b[horizon_axis][RIGHT];
 100   Real height = sky * b[other_axis (horizon_axis)][sky];
 101
 102   start_ = start;
 103   end_ = end;
 104   precompute (start, height, height, end);
 105 }
 106
 107 void
 108 Building::precompute (Real start, Real start_height, Real end_height, Real end)
 109 {
 110   slope_ = 0.0; /* if they were both infinite, we would get nan, not 0, from the prev line */
 111   if (start_height != end_height)
 112     slope_ = (end_height - start_height) / (end - start);
 113
 114   assert (!isinf (slope_) && !isnan (slope_));
 115
 116   if (isinf (start))
 117     {
 118       assert (start_height == end_height);
 119       y_intercept_ = start_height;
 120     }
 121   else
 122     y_intercept_ = start_height - slope_ * start;
 123 }
 124
 125 inline Real
 126 Building::height (Real x) const
 127 {
 128   return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
 129 }
 130
 131 void
 132 Building::print () const
 133 {
 134   printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
 135 }
 136
 137 inline Real
 138 Building::intersection_x (Building const &other) const
 139 {
 140   Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
 141   return isnan (ret) ? -infinity_f : ret;
 142 }
 143
 144 void
 145 Building::leading_part (Real chop)
 146 {
 147   assert (chop <= end_);
 148   end_ = chop;
 149 }
 150
 151 // Returns a shift s such that (x + s, y) intersects the roof of
 152 // this building.  If no such shift exists, returns infinity_f.
 153 Real
 154 Building::shift_to_intersect (Real x, Real y) const
 155 {
 156   // Solve for s: y = (x + s)*m + b
 157   Real ret = (y - y_intercept_ - slope_ * x) / slope_;
 158
 159   if (ret >= start_ && ret <= end_ && !isinf (ret))
 160     return ret;
 161   return infinity_f;
 162 }
 163
 164 // Returns the interval of horizontal shifts for which this
 165 // building (pointing up) overlaps the other building (pointing down).
 166 Interval
 167 Building::overlapping_shift_interval (Building const &other) const
 168 {
 169   Interval iv;
 170
 171   // If one building is empty, there will never be an overlap.
 172   if (y_intercept_ == -infinity_f || other.y_intercept_ == -infinity_f)
 173     return iv;
 174
 175   // There are two kinds of interesting positions:
 176   // - when the horizontal extents of the buildings just touch
 177   // - when an endpoint of one building intersects the roof of the other.
 178   // The interval we are looking for is the smallest one that
 179   // contains all of the interesting points.
 180
 181
 182   Real my_y1 = height (start_);
 183   Real my_y2 = height (end_);
 184   Real his_y1 = -other.height (other.start_); // "-" because OTHER points down
 185   Real his_y2 = -other.height (other.end_);
 186
 187   // If both buildings are infinite in the same direction,
 188   // the return value is either empty or full.
 189   if ((isinf (start_) && isinf (other.start_))
 190       || (isinf (end_) && isinf (other.end_)))
 191     return (y_intercept_ > other.y_intercept_)
 192            ? Interval (-infinity_f, infinity_f) : Interval ();
 193
 194   // ...when the horizontal extents of the buildings just touch...
 195   if (my_y1 >= his_y2)
 196     iv.add_point (other.end_ - start_);
 197   if (my_y2 >= his_y1)
 198     iv.add_point (other.start_ - end_);
 199
 200   // ...when an endpoint of one building intersects the roof of the other.
 201   Real p1 = shift_to_intersect (other.start_, his_y1);
 202   Real p2 = shift_to_intersect (other.end_, his_y2);
 203   // "-my_y1" because OTHER points down:
 204   Real p3 = other.shift_to_intersect (start_, -my_y1);
 205   Real p4 = other.shift_to_intersect (end_, -my_y2);
 206   if (!isinf (p1))
 207     iv.add_point (p1);
 208   if (!isinf (p2))
 209     iv.add_point (p2);
 210   if (!isinf (p3))
 211     iv.add_point (p3);
 212   if (!isinf (p4))
 213     iv.add_point (p4);
 214
 215   return iv;
 216 }
 217
 218 static Real
 219 first_intersection (Building const &b, list<Building> *const s, Real start_x)
 220 {
 221   while (!s->empty () && start_x < b.end_)
 222     {
 223       Building c = s->front ();
 224
 225       // conceals and intersection_x involve multiplication and
 226       // division. Avoid that, if we can.
 227       if (c.y_intercept_ == -infinity_f)
 228         {
 229           if (c.end_ > b.end_)
 230             return b.end_;
 231           start_x = c.end_;
 232           s->pop_front ();
 233           continue;
 234         }
 235
 236       if (c.conceals (b, start_x))
 237         return start_x;
 238
 239       Real i = b.intersection_x (c);
 240       if (i > start_x && i <= b.end_ && i <= c.end_)
 241         return i;
 242
 243       start_x = c.end_;
 244       if (b.end_ > c.end_)
 245         s->pop_front ();
 246     }
 247   return b.end_;
 248 }
 249
 250 bool
 251 Building::conceals (Building const &other, Real x) const
 252 {
 253   if (slope_ == other.slope_)
 254     return y_intercept_ > other.y_intercept_;
 255
 256   /* their slopes were not equal, so there is an intersection point */
 257   Real i = intersection_x (other);
 258   return (i <= x && slope_ > other.slope_)
 259          || (i > x && slope_ < other.slope_);
 260 }
 261
 262 // Remove redundant empty buildings from the skyline.
 263 // If there are two adjacent empty buildings, they can be
 264 // turned into one.
 265 void
 266 Skyline::normalize ()
 267 {
 268   bool last_empty = false;
 269   list<Building>::iterator i;
 270   for (i = buildings_.begin (); i != buildings_.end (); i++)
 271     {
 272       if (last_empty && i->y_intercept_ == -infinity_f)
 273         {
 274           list<Building>::iterator last = i;
 275           last--;
 276           last->end_ = i->end_;
 277           buildings_.erase (i);
 278           i = last;
 279         }
 280       last_empty = (i->y_intercept_ == -infinity_f);
 281     }
 282
 283   assert (buildings_.front ().start_ == -infinity_f);
 284   assert (buildings_.back ().end_ == infinity_f);
 285 }
 286
 287 void
 288 Skyline::deholify ()
 289 {
 290   // Since a skyline should always be normalized, we can
 291   // assume that there are never two adjacent empty buildings.
 292   // That is, if center is empty then left and right are not.
 293   list<Building>::iterator left = buildings_.begin ();
 294   list<Building>::iterator center = buildings_.begin ();
 295   list<Building>::iterator right;
 296
 297   for (right = buildings_.begin (); right != buildings_.end (); right++)
 298     {
 299       if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
 300         {
 301           Real p1 = left->height (left->end_);
 302           Real p2 = right->height (right->start_);
 303           *center = Building (center->start_, p1, p2, center->end_);
 304
 305           left = center;
 306           center = right;
 307         }
 308     }
 309 }
 310
 311 void
 312 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
 313                                  list<Building> *const result) const
 314 {
 315   if (s1->empty () || s2->empty ())
 316     {
 317       programming_error ("tried to merge an empty skyline");
 318       return;
 319     }
 320
 321   Real x = -infinity_f;
 322   Real last_end = -infinity_f;
 323   while (!s1->empty ())
 324     {
 325       if (s2->front ().conceals (s1->front (), x))
 326         swap (s1, s2);
 327
 328       Building b = s1->front ();
 329       Building c = s2->front ();
 330
 331       // Optimization: if the other skyline is empty at this point,
 332       // we can avoid testing some intersections. Just grab as many
 333       // buildings from s1 as we can, and shove them onto the output.
 334       if (c.y_intercept_ == -infinity_f
 335           && c.end_ >= b.end_)
 336         {
 337           list<Building>::iterator i = s1->begin ();
 338           i++;
 339           while (i != s1->end () && i->end_ <= c.end_)
 340             i++;
 341
 342           s1->front ().start_ = x;
 343           result->splice (result->end (), *s1, s1->begin (), i);
 344           x = result->back ().end_;
 345           last_end = x;
 346           continue;
 347         }
 348
 349       Real end = first_intersection (b, s2, x);
 350       if (s2->empty ())
 351         {
 352           b.start_ = last_end;
 353           result->push_back (b);
 354           break;
 355         }
 356
 357       /* only include buildings wider than epsilon */
 358       if (end > x + EPS)
 359         {
 360           b.leading_part (end);
 361           b.start_ = last_end;
 362           last_end = b.end_;
 363           result->push_back (b);
 364         }
 365
 366       if (end >= s1->front ().end_)
 367         s1->pop_front ();
 368
 369       x = end;
 370     }
 371 }
 372
 373 static void
 374 empty_skyline (list<Building> *const ret)
 375 {
 376   ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
 377 }
 378
 379 /*
 380   Given Building 'b', build a skyline containing only that building.
 381 */
 382 static void
 383 single_skyline (Building b, list<Building> *const ret)
 384 {
 385   if (b.end_ > b.start_ + EPS)
 386     {
 387       ret->push_back (Building (-infinity_f, -infinity_f,
 388                                 -infinity_f, b.start_));
 389       ret->push_back (b);
 390       ret->push_back (Building (b.end_, -infinity_f,
 391                                 -infinity_f, infinity_f));
 392     }
 393   else
 394     {
 395       empty_skyline (ret);
 396     }
 397 }
 398
 399 /* remove a non-overlapping set of boxes from BOXES and build a skyline
 400    out of them */
 401 static list<Building>
 402 non_overlapping_skyline (list<Building> *const buildings)
 403 {
 404   list<Building> result;
 405   Real last_end = -infinity_f;
 406   Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
 407   list<Building>::iterator i = buildings->begin ();
 408   while (i != buildings->end ())
 409     {
 410       Real x1 = i->start_;
 411       Real y1 = i->height (i->start_);
 412       Real x2 = i->end_;
 413       Real y2 = i->height (i->end_);
 414
 415       // Drop buildings that will obviously have no effect.
 416       if (last_building.height (x1) >= y1
 417           && last_building.end_ >= x2
 418           && last_building.height (x2) >= y2)
 419         {
 420           list<Building>::iterator j = i++;
 421           buildings->erase (j);
 422           continue;
 423         }
 424
 425       if (x1 < last_end)
 426         {
 427           i++;
 428           continue;
 429         }
 430
 431       if (x1 > last_end + EPS)
 432         result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
 433
 434       result.push_back (*i);
 435       last_building = *i;
 436       last_end = i->end_;
 437
 438       list<Building>::iterator j = i++;
 439       buildings->erase (j);
 440     }
 441
 442   if (last_end < infinity_f)
 443     result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
 444   return result;
 445 }
 446
 447 class LessThanBuilding
 448 {
 449 public:
 450   bool operator () (Building const &b1, Building const &b2)
 451   {
 452     return b1.start_ < b2.start_
 453            || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
 454   }
 455 };
 456
 457 /**
 458    BUILDINGS is a list of buildings, but they could be overlapping
 459    and in any order.  The returned list of buildings is ordered and non-overlapping.
 460 */
 461 list<Building>
 462 Skyline::internal_build_skyline (list<Building> *buildings) const
 463 {
 464   vsize size = buildings->size ();
 465
 466   if (size == 0)
 467     {
 468       list<Building> result;
 469       empty_skyline (&result);
 470       return result;
 471     }
 472   else if (size == 1)
 473     {
 474       list<Building> result;
 475       single_skyline (buildings->front (), &result);
 476       return result;
 477     }
 478
 479   deque<list<Building> > partials;
 480   buildings->sort (LessThanBuilding ());
 481   while (!buildings->empty ())
 482     partials.push_back (non_overlapping_skyline (buildings));
 483
 484   /* we'd like to say while (partials->size () > 1) but that's O (n).
 485      Instead, we exit in the middle of the loop */
 486   while (!partials.empty ())
 487     {
 488       list<Building> merged;
 489       list<Building> one = partials.front ();
 490       partials.pop_front ();
 491       if (partials.empty ())
 492         return one;
 493
 494       list<Building> two = partials.front ();
 495       partials.pop_front ();
 496       internal_merge_skyline (&one, &two, &merged);
 497       partials.push_back (merged);
 498     }
 499   assert (0);
 500   return list<Building> ();
 501 }
 502
 503 Skyline::Skyline ()
 504 {
 505   sky_ = UP;
 506   empty_skyline (&buildings_);
 507 }
 508
 509 Skyline::Skyline (Skyline const &src)
 510 {
 511   sky_ = src.sky_;
 512
 513   /* doesn't a list's copy constructor do this? -- jneem */
 514   for (list<Building>::const_iterator i = src.buildings_.begin ();
 515        i != src.buildings_.end (); i++)
 516     {
 517       buildings_.push_back (Building ((*i)));
 518     }
 519 }
 520
 521 Skyline::Skyline (Direction sky)
 522 {
 523   sky_ = sky;
 524   empty_skyline (&buildings_);
 525 }
 526
 527 /*
 528   Build skyline from a set of boxes.
 529
 530   Boxes should have fatness in the horizon_axis, otherwise they are ignored.
 531  */
 532 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
 533 {
 534   list<Building> buildings;
 535   sky_ = sky;
 536
 537   Axis vert_axis = other_axis (horizon_axis);
 538   for (vsize i = 0; i < boxes.size (); i++)
 539     {
 540       Interval iv = boxes[i][horizon_axis];
 541       if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
 542         buildings.push_front (Building (boxes[i], horizon_axis, sky));
 543     }
 544
 545   buildings_ = internal_build_skyline (&buildings);
 546   normalize ();
 547 }
 548
 549 /*
 550   build skyline from a set of line segments.
 551
 552   Buildings should have fatness in the horizon_axis, otherwise they are ignored.
 553  */
 554 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
 555 {
 556   list<Building> buildings;
 557   sky_ = sky;
 558
 559   for (vsize i = 0; i < segments.size (); i++)
 560     {
 561       Drul_array<Offset> const &seg = segments[i];
 562       Offset left = seg[LEFT];
 563       Offset right = seg[RIGHT];
 564       if (left[horizon_axis] > right[horizon_axis])
 565         swap (left, right);
 566
 567       Real x1 = left[horizon_axis];
 568       Real x2 = right[horizon_axis];
 569       Real y1 = left[other_axis (horizon_axis)] * sky;
 570       Real y2 = right[other_axis (horizon_axis)] * sky;
 571
 572       if (x1 + EPS < x2)
 573         buildings.push_back (Building (x1, y1, y2, x2));
 574     }
 575
 576   buildings_ = internal_build_skyline (&buildings);
 577   normalize ();
 578 }
 579
 580 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
 581 {
 582   sky_ = sky;
 583
 584   deque<Skyline> partials;
 585   for (vsize i = 0; i < skypairs.size (); i++)
 586     partials.push_back (Skyline ((skypairs[i])[sky]));
 587
 588   while (partials.size () > 1)
 589     {
 590       Skyline one = partials.front ();
 591       partials.pop_front ();
 592       Skyline two = partials.front ();
 593       partials.pop_front ();
 594
 595       one.merge (two);
 596       partials.push_back (one);
 597     }
 598
 599   if (partials.size ())
 600     buildings_.swap (partials.front ().buildings_);
 601   else
 602     buildings_.clear ();
 603 }
 604
 605 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
 606 {
 607   sky_ = sky;
 608   Building front (b, horizon_axis, sky);
 609   single_skyline (front, &buildings_);
 610 }
 611
 612 void
 613 Skyline::merge (Skyline const &other)
 614 {
 615   assert (sky_ == other.sky_);
 616
 617   if (other.is_empty ())
 618     return;
 619
 620   if (is_empty ())
 621     {
 622       buildings_ = other.buildings_;
 623       return;
 624     }
 625
 626   list<Building> other_bld (other.buildings_);
 627   list<Building> my_bld;
 628   my_bld.splice (my_bld.begin (), buildings_);
 629   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 630   normalize ();
 631 }
 632
 633 void
 634 Skyline::insert (Box const &b, Axis a)
 635 {
 636   list<Building> other_bld;
 637   list<Building> my_bld;
 638
 639   if (isnan (b[other_axis (a)][LEFT])
 640       || isnan (b[other_axis (a)][RIGHT]))
 641     {
 642       programming_error ("insane box for skyline");
 643       return;
 644     }
 645
 646   /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
 647   Interval iv = b[a];
 648   if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
 649     return;
 650
 651   my_bld.splice (my_bld.begin (), buildings_);
 652   single_skyline (Building (b, a, sky_), &other_bld);
 653   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 654   normalize ();
 655 }
 656
 657 void
 658 Skyline::raise (Real r)
 659 {
 660   list<Building>::iterator end = buildings_.end ();
 661   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 662     i->y_intercept_ += sky_ * r;
 663 }
 664
 665 void
 666 Skyline::shift (Real s)
 667 {
 668   list<Building>::iterator end = buildings_.end ();
 669   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 670     {
 671       i->start_ += s;
 672       i->end_ += s;
 673       i->y_intercept_ -= s * i->slope_;
 674     }
 675 }
 676
 677 Real
 678 Skyline::distance (Skyline const &other, Real horizon_padding) const
 679 {
 680   Real dummy;
 681   return internal_distance (other, horizon_padding, &dummy);
 682 }
 683
 684 Real
 685 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
 686 {
 687   Real touch;
 688   internal_distance (other, horizon_padding, &touch);
 689   return touch;
 690 }
 691
 692 Real
 693 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
 694 {
 695   if (horizon_padding == 0.0)
 696     return internal_distance (other, touch_point);
 697
 698   // Note that it is not necessary to build a padded version of other,
 699   // because the same effect can be achieved just by doubling horizon_padding.
 700   Skyline padded_this = padded (horizon_padding);
 701   return padded_this.internal_distance (other, touch_point);
 702 }
 703
 704 Skyline
 705 Skyline::padded (Real horizon_padding) const
 706 {
 707   list<Building> pad_buildings;
 708   for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
 709     {
 710       if (i->start_ > -infinity_f)
 711         {
 712           Real height = i->height (i->start_);
 713           if (height > -infinity_f)
 714             {
 715               // Add the sloped building that pads the left side of the current building.
 716               Real start = i->start_ - 2 * horizon_padding;
 717               Real end = i->start_ - horizon_padding;
 718               pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
 719
 720               // Add the flat building that pads the left side of the current building.
 721               start = i->start_ - horizon_padding;
 722               end = i->start_;
 723               pad_buildings.push_back (Building (start, height, height, end));
 724             }
 725         }
 726
 727       if (i->end_ < infinity_f)
 728         {
 729           Real height = i->height (i->end_);
 730           if (height > -infinity_f)
 731             {
 732               // Add the flat building that pads the right side of the current building.
 733               Real start = i->end_;
 734               Real end = start + horizon_padding;
 735               pad_buildings.push_back (Building (start, height, height, end));
 736
 737               // Add the sloped building that pads the right side of the current building.
 738               start = end;
 739               end += horizon_padding;
 740               pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
 741             }
 742         }
 743     }
 744
 745   // The buildings may be overlapping, so resolve that.
 746   list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
 747
 748   // Merge the padding with the original, to make a new skyline.
 749   Skyline padded (sky_);
 750   list<Building> my_buildings = buildings_;
 751   padded.buildings_.clear ();
 752   internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
 753   padded.normalize ();
 754
 755   return padded;
 756 }
 757
 758 Real
 759 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
 760 {
 761   assert (sky_ == -other.sky_);
 762
 763   list<Building>::const_iterator i = buildings_.begin ();
 764   list<Building>::const_iterator j = other.buildings_.begin ();
 765
 766   Real dist = -infinity_f;
 767   Real start = -infinity_f;
 768   Real touch = -infinity_f;
 769   while (i != buildings_.end () && j != other.buildings_.end ())
 770     {
 771       Real end = min (i->end_, j->end_);
 772       Real start_dist = i->height (start) + j->height (start);
 773       Real end_dist = i->height (end) + j->height (end);
 774       dist = max (dist, max (start_dist, end_dist));
 775
 776       if (end_dist == dist)
 777         touch = end;
 778       else if (start_dist == dist)
 779         touch = start;
 780
 781       if (i->end_ <= j->end_)
 782         i++;
 783       else
 784         j++;
 785       start = end;
 786     }
 787
 788   *touch_point = touch;
 789   return dist;
 790 }
 791
 792 // changes the direction that the skyline is pointing
 793 void
 794 Skyline::invert ()
 795 {
 796   list<Building>::iterator i;
 797   for (i = buildings_.begin (); i != buildings_.end (); i++)
 798     if (!isinf (i->y_intercept_))
 799       {
 800         i->y_intercept_ *= -1;
 801         i->slope_ *= -1;
 802       }
 803
 804   sky_ = -sky_;
 805 }
 806
 807 Real
 808 Skyline::height (Real airplane) const
 809 {
 810   assert (!isinf (airplane));
 811
 812   list<Building>::const_iterator i;
 813   for (i = buildings_.begin (); i != buildings_.end (); i++)
 814     {
 815       if (i->end_ >= airplane)
 816         return sky_ * i->height (airplane);
 817     }
 818
 819   assert (0);
 820   return 0;
 821 }
 822
 823 Real
 824 Skyline::max_height () const
 825 {
 826   Real ret = -infinity_f;
 827
 828   list<Building>::const_iterator i;
 829   for (i = buildings_.begin (); i != buildings_.end (); ++i)
 830     {
 831       ret = max (ret, i->height (i->start_));
 832       ret = max (ret, i->height (i->end_));
 833     }
 834
 835   return sky_ * ret;
 836 }
 837
 838 Real
 839 Skyline::max_height_position () const
 840 {
 841   Skyline s (-sky_);
 842   s.set_minimum_height (0);
 843   return touching_point (s);
 844 }
 845
 846 void
 847 Skyline::set_minimum_height (Real h)
 848 {
 849   Skyline s (sky_);
 850   s.buildings_.front ().y_intercept_ = h * sky_;
 851   merge (s);
 852 }
 853
 854 vector<Offset>
 855 Skyline::to_points (Axis horizon_axis) const
 856 {
 857   vector<Offset> out;
 858
 859   Real start = -infinity_f;
 860   for (list<Building>::const_iterator i (buildings_.begin ());
 861        i != buildings_.end (); i++)
 862     {
 863       out.push_back (Offset (start, sky_ * i->height (start)));
 864       out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
 865       start = i->end_;
 866     }
 867
 868   if (horizon_axis == Y_AXIS)
 869     for (vsize i = 0; i < out.size (); i++)
 870       out[i] = out[i].swapped ();
 871
 872   return out;
 873 }
 874
 875 // Returns the smallest (non-negative) shift in the given
 876 // direction which will result in THIS and OTHER not overlapping.
 877 // Warning: this function is O(n^2 log n). Use sparingly.
 878 Real
 879 Skyline::smallest_shift (Skyline const &other,
 880                          Direction d,
 881                          Real horizon_padding,
 882                          Real vertical_padding) const
 883 {
 884   // If one or both of the paddings is zero, this can
 885   // be optimized...
 886   Skyline padded_me = padded (horizon_padding);
 887   padded_me.raise (vertical_padding);
 888
 889   list<Building>::const_iterator i = padded_me.buildings_.begin ();
 890   list<Building>::const_iterator j = other.buildings_.begin ();
 891   list<Building>::const_iterator i_end = padded_me.buildings_.end ();
 892   list<Building>::const_iterator j_end = other.buildings_.end ();
 893
 894   // Find all shifts that are not allowed.
 895   vector<Interval> forbidden_shifts;
 896   for (; i != i_end; ++i)
 897     if (i->y_intercept_ != -infinity_f)
 898       for (j = other.buildings_.begin (); j != j_end; ++j)
 899         {
 900           Interval iv = i->overlapping_shift_interval (*j);
 901           if (!iv.is_empty ())
 902             forbidden_shifts.push_back (iv);
 903         }
 904
 905   // Now comes the trick: we want to find the smallest point
 906   // that is not in the union of forbidden_shifts. We can represent
 907   // the union of forbidden_shifts as a skyline, where a point is
 908   // allowed if it has height -infinity_f and forbidden otherwise.
 909   vector<Box> boxes;
 910   for (vector<Interval>::iterator k = forbidden_shifts.begin ();
 911        k != forbidden_shifts.end (); ++k)
 912     boxes.push_back (Box (*k, Interval (-1, 0)));
 913   Skyline s (boxes, X_AXIS, UP);
 914
 915   // Find the smallest (ie. closest to zero, in the appropriate direction)
 916   // coordinate where the height of s is -infinity_f.
 917   Real last_good_point = -infinity_f;
 918   for (i = s.buildings_.begin (); i != s.buildings_.end (); ++i)
 919     {
 920       if (d == LEFT && i->start_ > 0)
 921         return last_good_point;
 922
 923       if (i->y_intercept_ == -infinity_f)
 924         {
 925           if (i->start_ <= 0 && i->end_ >= 0)
 926             return 0;
 927           if (d == RIGHT && i->start_ >= 0)
 928             return i->start_;
 929
 930           last_good_point = i->end_;
 931         }
 932     }
 933
 934   return infinity_f * d;
 935 }
 936
 937 Real
 938 Skyline::left () const
 939 {
 940   for (list<Building>::const_iterator i (buildings_.begin ());
 941        i != buildings_.end (); i++)
 942     if (i->y_intercept_ > -infinity_f)
 943       return i->start_;
 944
 945   return infinity_f;
 946 }
 947
 948 Real
 949 Skyline::right () const
 950 {
 951   for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
 952        i != buildings_.rend (); ++i)
 953     if (i->y_intercept_ > -infinity_f)
 954       return i->end_;
 955
 956   return -infinity_f;
 957 }
 958
 959 bool
 960 Skyline::is_empty () const
 961 {
 962   if (!buildings_.size ())
 963     return true;
 964   Building b = buildings_.front ();
 965   return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
 966 }
 967
 968 bool
 969 Skyline::is_singleton () const
 970 {
 971   return buildings_.size () == 3;
 972 }
 973
 974 void
 975 Skyline::clear ()
 976 {
 977   buildings_.clear ();
 978   empty_skyline (&buildings_);
 979 }
 980
 981 /****************************************************************/
 982
 983 IMPLEMENT_SIMPLE_SMOBS (Skyline);
 984 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
 985 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
 986
 987 SCM
 988 Skyline::mark_smob (SCM s)
 989 {
 990   ASSERT_LIVE_IS_ALLOWED (s);
 991   return SCM_EOL;
 992 }
 993
 994 int
 995 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
 996 {
 997   Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
 998   (void) r;
 999
1000   scm_puts ("#<Skyline>", port);
1001
1002   return 1;
1003 }
1004
1005 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
1006 SCM
1007 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
1008 {
1009   LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
1010
1011   Real horizon_padding = 0;
1012   if (horizon_padding_scm != SCM_UNDEFINED)
1013     {
1014       LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
1015       horizon_padding = scm_to_double (horizon_padding_scm);
1016     }
1017
1018   Skyline *skyline = Skyline::unsmob (skyline_scm);
1019   Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
1020   return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
1021 }
1022
1023 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
1024 SCM
1025 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
1026 {
1027   LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
1028
1029   Real horizon_padding = 0;
1030   if (horizon_padding_scm != SCM_UNDEFINED)
1031     {
1032       LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
1033       horizon_padding = scm_to_double (horizon_padding_scm);
1034     }
1035
1036   Skyline *skyline = Skyline::unsmob (skyline_scm);
1037   Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
1038   return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
1039 }
1040
1041 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
1042 SCM
1043 Skyline::get_max_height (SCM skyline_scm)
1044 {
1045   return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
1046 }
1047
1048 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
1049 SCM
1050 Skyline::get_max_height_position (SCM skyline_scm)
1051 {
1052   return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
1053 }
1054
1055 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
1056 SCM
1057 Skyline::get_height (SCM skyline_scm, SCM x_scm)
1058 {
1059   Real x = robust_scm2double (x_scm, 0.0);
1060   return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
1061 }