lily/skyline.cc

   1 /* skyline.cc -- implement the Skyline class
   2
   3    source file of the GNU LilyPond music typesetter
   4
   5    (c) 2006--2007 Joe Neeman <joeneeman@gmail.com>
   6 */
   7
   8 #include "skyline.hh"
   9 #include <deque>
  10
  11 #include "ly-smobs.icc"
  12
  13 /* A skyline is a sequence of non-overlapping buildings: something like
  14    this:
  15                    _______
  16                   |       \                                 ________
  17                   |        \                       ________/        \
  18         /\        |          \                    /                  \
  19        /  --------             \                 /                    \
  20       /                          \              /                      \
  21      /                             ------------/                        ----
  22    --
  23    Each building has a starting position, and ending position, a starting
  24    height and an ending height.
  25
  26    The following invariants are observed:
  27     - the start of the first building is at -infinity
  28     - the end of the last building is at infinity
  29     - if a building has infinite length (ie. the first and last buildings),
  30       then its starting height and ending height are equal
  31     - the end of one building is the same as the beginning of the next
  32       building
  33
  34    We also allow skylines to point down (the structure is exactly the same,
  35    but we think of the part above the line as being filled with mass and the
  36    part below as being empty). ::distance finds the minimum distance between
  37    an UP skyline and a DOWN skyline.
  38
  39    Note that we store DOWN skylines upside-down. That is, in order to compare
  40    a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
  41    This means that the merging routine doesn't need to be aware of direction,
  42    but the distance routine does.
  43 */
  44
  45 /*
  46   FIXME:
  47
  48   * Consider to use
  49
  50   typedef list<Skyline_point> Skyline;
  51   struct Skyline_point
  52   {
  53     Real x;
  54     Drul_array<Real> ys;
  55   };
  56
  57   this is a cleaner representation, as it doesn't duplicate the X, and
  58   doesn't need bogus buildings at infinity  --hwn.
  59
  60
  61   * All the messing around with EPS is very fishy.  There are no
  62   complicated numerical algorithms involved, so EPS should not be
  63   necessary.
  64
  65   --hwn
  66
  67
  68  */
  69
  70 #define EPS 1e-10
  71
  72 static inline bool
  73 approx_equal (Real x, Real y)
  74 {
  75   return abs (x - y) < EPS || (isinf (x) && isinf (y) && ((x > 0) == (y > 0)));
  76 }
  77
  78 static inline bool
  79 approx_greater_than (Real x, Real y)
  80 {
  81   return x > y + EPS;
  82 }
  83
  84 static inline bool
  85 approx_less_than (Real x, Real y)
  86 {
  87   return x < y - EPS;
  88 }
  89
  90 static inline bool
  91 approx_less_equal (Real x, Real y)
  92 {
  93   return x <= y + EPS;
  94 }
  95
  96 static inline bool
  97 approx_greater_equal (Real x, Real y)
  98 {
  99   return x >= y - EPS;
 100 }
 101
 102 void
 103 Skyline::print () const
 104 {
 105   for (list<Building>::const_iterator i = buildings_.begin ();
 106        i != buildings_.end (); i++)
 107     {
 108       (*i).print ();
 109     }
 110 }
 111
 112 bool
 113 is_legal_skyline (list<Building> const &buildings)
 114 {
 115   list<Building>::const_iterator i;
 116   Real last_x = -infinity_f;
 117   for (i = buildings.begin (); i != buildings.end (); i++)
 118     {
 119       if (i->iv_[LEFT] != last_x)
 120         return false;
 121       last_x = i->iv_[RIGHT];
 122       if (isinf (i->iv_.length ()) && i->height_[LEFT] != i->height_[RIGHT])
 123         return false;
 124     }
 125   return last_x == infinity_f;
 126 }
 127
 128 Building::Building (Real start, Real start_height, Real end_height, Real end)
 129   : iv_ (start, end)
 130 {
 131   height_[LEFT] = start_height;
 132   height_[RIGHT] = end_height;
 133
 134   if (isinf (start) || isinf (end))
 135     assert (start_height == end_height);
 136
 137   precompute ();
 138 }
 139
 140 Building::Building (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
 141 {
 142   Real height = sky * b[other_axis (horizon_axis)][sky];
 143
 144   iv_ = b[horizon_axis];
 145   iv_.widen (horizon_padding + EPS);
 146   height_[LEFT] = height;
 147   height_[RIGHT] = height;
 148
 149   if (sane ())
 150     precompute ();
 151 }
 152
 153 void
 154 Building::precompute ()
 155 {
 156   slope_ = (height_[RIGHT] - height_[LEFT]) / (iv_.length());
 157   if (height_[LEFT] == height_[RIGHT]) /* in case they're both infinity */
 158     slope_ = 0;
 159
 160   assert (!isinf (slope_) && !isnan (slope_));
 161
 162   if (isinf (iv_[START]))
 163     {
 164       assert (slope_ == 0);
 165       y_intercept_ = height_[LEFT];
 166     }
 167   else
 168     y_intercept_ = height_[LEFT] - slope_ * iv_[START];
 169 }
 170
 171 Real
 172 Building::height (Real x) const
 173 {
 174   if (isinf (x))
 175     return (x > 0) ? height_[RIGHT] : height_[LEFT];
 176   return slope_*x + y_intercept_;
 177 }
 178
 179 void
 180 Building::print () const
 181 {
 182   printf ("X[%f,%f] -> Y[%f,%f]\n",
 183           iv_[LEFT], iv_[RIGHT],
 184           height_[LEFT], height_[RIGHT]);
 185 }
 186
 187 Real
 188 Building::intersection_x (Building const &other) const
 189 {
 190   return (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
 191 }
 192
 193 void
 194 Building::leading_part (Real chop)
 195 {
 196   assert (chop > iv_[LEFT] && chop <= iv_[RIGHT] && !approx_equal (chop, iv_[LEFT]));
 197   iv_[RIGHT] = chop;
 198   height_[RIGHT] = height (chop);
 199 }
 200
 201 Building
 202 Building::sloped_neighbour (Real horizon_padding, Direction d) const
 203 {
 204   Real left = iv_[d];
 205   Real right = iv_[d] + d * horizon_padding;
 206   Real left_height = height_[d];
 207   Real right_height = height_[d] - horizon_padding;
 208   if (d == LEFT)
 209     {
 210       swap (left, right);
 211       swap (left_height, right_height);
 212     }
 213   return Building (left, left_height, right_height, right);
 214 }
 215
 216 bool
 217 Building::sane () const
 218 {
 219   return approx_less_than (iv_[LEFT], iv_[RIGHT])
 220     && !isinf (height_[RIGHT])
 221     && !isinf (height_[LEFT]);
 222 }
 223
 224 static void
 225 skyline_trailing_part (list<Building> *sky, Real x)
 226 {
 227   if (approx_equal (x, sky->front ().iv_[RIGHT]))
 228     sky->pop_front ();
 229   else
 230     assert (x < sky->front ().iv_[RIGHT]);
 231
 232   if (!sky->empty ())
 233     {
 234       sky->front ().iv_[LEFT] = x;
 235       sky->front ().height_[LEFT] = sky->front ().height (x);
 236     }
 237 }
 238
 239 bool
 240 Building::conceals_beginning (Building const &other) const
 241 {
 242   bool w = false;
 243   Real h = other.height (iv_[LEFT]);
 244   if (approx_equal (height_[LEFT], h))
 245     w = slope_ > other.slope_;
 246   else if (height_[LEFT] > h)
 247     w = true;
 248   else
 249     w = false;
 250
 251   return w;
 252 }
 253
 254 bool
 255 Building::conceals (Building const &other) const
 256 {
 257   assert (iv_[LEFT] <= other.iv_[LEFT]);
 258   return (iv_[RIGHT] >= other.iv_[RIGHT])
 259     && approx_greater_equal (height (other.iv_[LEFT]), other.height_[LEFT])
 260     && approx_greater_equal (height (other.iv_[RIGHT]), other.height_[RIGHT]);
 261 }
 262
 263 void
 264 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
 265                                  list<Building> *const result)
 266 {
 267   while (!s1->empty ())
 268     {
 269       if (s2->front ().conceals_beginning (s1->front ()))
 270         swap (s1, s2);
 271
 272       Building b = s1->front ();
 273       while (!s2->empty () && b.conceals (s2->front ()))
 274         s2->pop_front ();
 275       if (s2->empty ())
 276         {
 277           result->push_front (b);
 278           break;
 279         }
 280
 281       /* s2 either intersects with b or it ends after b */
 282       Real end = infinity_f;
 283       Real s2_start_height = s2->front ().height_[LEFT];
 284       Real s2_end_height = s2->front ().height_[RIGHT];
 285       Real s1_start_height = b.height (s2->front ().iv_[LEFT]);
 286       Real s1_end_height = b.height (s2->front ().iv_[RIGHT]);
 287       if (approx_greater_than (s2_start_height, s1_start_height))
 288         end = s2->front ().iv_[LEFT];
 289       else if (approx_greater_than (s2_end_height, s1_end_height))
 290         end = b.intersection_x (s2->front ());
 291       end = min (end, b.iv_[RIGHT]);
 292
 293       b.leading_part (end);
 294       result->push_front (b);
 295
 296       skyline_trailing_part (s1, end);
 297       skyline_trailing_part (s2, end);
 298     }
 299   result->reverse ();
 300 }
 301
 302 static void
 303 empty_skyline (list<Building> *const ret)
 304 {
 305   ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
 306 }
 307
 308 static void
 309 single_skyline (Building b, Real horizon_padding, list<Building> *const ret)
 310 {
 311   b.iv_.widen (horizon_padding);
 312
 313   if (!isinf (b.iv_[RIGHT]))
 314     ret->push_front (Building (b.iv_[RIGHT], -infinity_f,
 315                                -infinity_f, infinity_f));
 316   if (horizon_padding > 0 && !isinf (b.iv_.length ()))
 317     ret->push_front (b.sloped_neighbour (horizon_padding, RIGHT));
 318
 319   if (b.iv_[RIGHT] > b.iv_[LEFT])
 320     ret->push_front (b);
 321
 322   if (horizon_padding > 0 && !isinf (b.iv_.length ()))
 323     ret->push_front (b.sloped_neighbour (horizon_padding, LEFT));
 324   if (!isinf (b.iv_[LEFT]))
 325     ret->push_front (Building (-infinity_f, -infinity_f,
 326                                -infinity_f, b.iv_[LEFT]));
 327 }
 328
 329 /* remove a non-overlapping set of buildings from BUILDINGS and build a skyline
 330    out of them */
 331 static list<Building>
 332 non_overlapping_skyline (list<Building> *const buildings)
 333 {
 334   list<Building> result;
 335   Real last_end = -infinity_f;
 336   list<Building>::iterator i = buildings->begin ();
 337   while (i != buildings->end ())
 338     {
 339       if (approx_less_than (i->iv_[LEFT], last_end))
 340         {
 341           i++;
 342           continue;
 343         }
 344
 345       if (approx_greater_than (i->iv_[LEFT], last_end))
 346         result.push_back (Building (last_end, -infinity_f, -infinity_f, i->iv_[LEFT]));
 347       else
 348         i->iv_[LEFT] = last_end;
 349
 350       last_end = i->iv_[RIGHT];
 351       list<Building>::iterator j = i;
 352       i++;
 353       result.splice (result.end (), *buildings, j);
 354     }
 355   if (last_end < infinity_f)
 356     result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
 357   assert (is_legal_skyline (result));
 358   return result;
 359 }
 360
 361 list<Building>
 362 Skyline::internal_build_skyline (list<Building> *buildings)
 363 {
 364   vsize size = buildings->size ();
 365
 366   if (size == 0)
 367     {
 368       list<Building> result;
 369       empty_skyline (&result);
 370       return result;
 371     }
 372   else if (size == 1)
 373     {
 374       list<Building> result;
 375       single_skyline (buildings->front (), 0, &result);
 376       return result;
 377     }
 378
 379   deque<list<Building> > partials;
 380   buildings->sort ();
 381   while (!buildings->empty ())
 382     partials.push_back (non_overlapping_skyline (buildings));
 383
 384   /* we'd like to say while (partials->size () > 1) but that's O(n).
 385      Instead, we exit in the middle of the loop */
 386   while (!partials.empty ())
 387     {
 388       list<Building> merged;
 389       list<Building> one = partials.front ();
 390       partials.pop_front ();
 391       if (partials.empty ())
 392         return one;
 393
 394       list<Building> two = partials.front ();
 395       partials.pop_front ();
 396       internal_merge_skyline (&one, &two, &merged);
 397       partials.push_back (merged);
 398     }
 399   assert (0);
 400   return list<Building> ();
 401 }
 402
 403 Skyline::Skyline ()
 404 {
 405   sky_ = UP;
 406   empty_skyline (&buildings_);
 407 }
 408
 409 Skyline::Skyline (Skyline const &src)
 410 {
 411   sky_ = src.sky_;
 412
 413   for (list<Building>::const_iterator i = src.buildings_.begin ();
 414        i != src.buildings_.end (); i++)
 415     {
 416       buildings_.push_back (Building ((*i)));
 417     }
 418 }
 419
 420 Skyline::Skyline (Direction sky)
 421 {
 422   sky_ = sky;
 423   empty_skyline (&buildings_);
 424 }
 425
 426 /*
 427   build skyline from a set of boxes. If horizon_padding > 0, expand all the boxes
 428   by that amount and add 45-degree sloped boxes to the edges of each box (of
 429   width horizon_padding). That is, the total amount of horizontal expansion is
 430   horizon_padding*4, half of which is sloped and half of which is flat.
 431
 432   Boxes should have fatness in the horizon_axis (after they are expanded by
 433   horizon_padding), otherwise they are ignored.
 434  */
 435 Skyline::Skyline (vector<Box> const &boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 436 {
 437   list<Building> bldgs;
 438   sky_ = sky;
 439
 440   for (vsize i = 0; i < boxes.size (); i++)
 441     {
 442       Building front (boxes[i], horizon_padding, horizon_axis, sky);
 443       if (front.sane ())
 444         {
 445           bldgs.push_front (front);
 446           if (horizon_padding > 0 && !isinf (front.iv_.length ()))
 447             {
 448               bldgs.push_front (front.sloped_neighbour (horizon_padding, LEFT));
 449               bldgs.push_front (front.sloped_neighbour (horizon_padding, RIGHT));
 450             }
 451         }
 452     }
 453
 454   buildings_ = internal_build_skyline (&bldgs);
 455   assert (is_legal_skyline (buildings_));
 456 }
 457
 458 Skyline::Skyline (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
 459 {
 460   sky_ = sky;
 461   Building front (b, 0, horizon_axis, sky);
 462   single_skyline (front, horizon_padding, &buildings_);
 463 }
 464
 465 void
 466 Skyline::merge (Skyline const &other)
 467 {
 468   assert (sky_ == other.sky_);
 469
 470   list<Building> other_bld (other.buildings_);
 471   list<Building> my_bld;
 472   my_bld.splice (my_bld.begin (), buildings_);
 473   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 474   assert (is_legal_skyline (buildings_));
 475 }
 476
 477 void
 478 Skyline::insert (Box const &b, Real horizon_padding, Axis a)
 479 {
 480   list<Building> other_bld;
 481   list<Building> my_bld;
 482
 483   my_bld.splice (my_bld.begin (), buildings_);
 484   single_skyline (Building (b, 0, a, sky_), horizon_padding, &other_bld);
 485   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 486   assert (is_legal_skyline (buildings_));
 487 }
 488
 489 void
 490 Skyline::raise (Real r)
 491 {
 492   list<Building>::iterator end = buildings_.end ();
 493   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 494     {
 495       i->height_[LEFT] += sky_ * r;
 496       i->height_[RIGHT] += sky_ * r;
 497       i->y_intercept_ += sky_ * r;
 498     }
 499   assert (is_legal_skyline (buildings_));
 500 }
 501
 502 void
 503 Skyline::shift (Real r)
 504 {
 505   list<Building>::iterator end = buildings_.end ();
 506   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 507     {
 508       i->iv_[LEFT] += r;
 509       i->iv_[RIGHT] += r;
 510     }
 511 }
 512
 513 Real
 514 Skyline::distance (Skyline const &other) const
 515 {
 516   assert (sky_ == -other.sky_);
 517   list<Building>::const_iterator i = buildings_.begin ();
 518   list<Building>::const_iterator j = other.buildings_.begin ();
 519
 520   Real dist = -infinity_f;
 521   while (i != buildings_.end () && j != other.buildings_.end ())
 522     {
 523       Interval iv = intersection (i->iv_, j->iv_);
 524       dist = max (dist, max (i->height (iv[LEFT]) + j->height (iv[LEFT]),
 525                              i->height (iv[RIGHT]) + j->height (iv[RIGHT])));
 526       if (i->iv_[RIGHT] <= j->iv_[RIGHT])
 527         i++;
 528       else
 529         j++;
 530     }
 531   return dist;
 532 }
 533
 534 Real
 535 Skyline::height (Real airplane) const
 536 {
 537   assert (!isinf (airplane));
 538
 539   list<Building>::const_iterator i;
 540   for (i = buildings_.begin (); i != buildings_.end (); i++)
 541     {
 542       if (i->iv_[RIGHT] >= airplane)
 543         return sky_ * i->height (airplane);
 544     }
 545
 546   assert (0);
 547   return 0;
 548 }
 549
 550 Real
 551 Skyline::max_height () const
 552 {
 553   Skyline s (-sky_);
 554   s.set_minimum_height (0);
 555   return sky_ * distance (s);
 556 }
 557
 558 void
 559 Skyline::set_minimum_height (Real h)
 560 {
 561   Skyline s (sky_);
 562   s.buildings_.front ().height_[LEFT] = h * sky_;
 563   s.buildings_.front ().height_[RIGHT] = h * sky_;
 564   s.buildings_.front ().y_intercept_ = h * sky_;
 565   merge (s);
 566 }
 567
 568
 569 vector<Offset>
 570 Skyline::to_points () const
 571 {
 572   vector<Offset> out;
 573
 574   for (list<Building>::const_iterator i (buildings_.begin ());
 575        i != buildings_.end (); i++)
 576     {
 577       if (!isinf (i->iv_[LEFT]) && !isinf (i->height_[LEFT]))
 578         out.push_back (Offset (i->iv_[LEFT], sky_ * i->height_[LEFT]));
 579       if (!isinf (i->iv_[RIGHT]) && !isinf (i->height_[RIGHT]))
 580         out.push_back (Offset (i->iv_[RIGHT], sky_ * i->height_[RIGHT]));
 581     }
 582   return out;
 583 }
 584
 585 Skyline_pair::Skyline_pair ()
 586   : skylines_ (Skyline (DOWN), Skyline (UP))
 587 {
 588 }
 589
 590 Skyline_pair::Skyline_pair (vector<Box> const &boxes, Real padding, Axis a)
 591   : skylines_ (Skyline (boxes, padding, a, DOWN), Skyline (boxes, padding, a, UP))
 592 {
 593 }
 594
 595 Skyline_pair::Skyline_pair (Box const &b, Real padding, Axis a)
 596   : skylines_ (Skyline (b, padding, a, DOWN), Skyline (b, padding, a, UP))
 597 {
 598 }
 599
 600 void
 601 Skyline_pair::raise (Real r)
 602 {
 603   skylines_[UP].raise (r);
 604   skylines_[DOWN].raise (r);
 605 }
 606
 607 void
 608 Skyline_pair::shift (Real r)
 609 {
 610   skylines_[UP].shift (r);
 611   skylines_[DOWN].shift (r);
 612 }
 613
 614 void
 615 Skyline_pair::insert (Box const &b, Real padding, Axis a)
 616 {
 617   skylines_[UP].insert (b, padding, a);
 618   skylines_[DOWN].insert (b, padding, a);
 619 }
 620
 621 void
 622 Skyline_pair::merge (Skyline_pair const &other)
 623 {
 624   skylines_[UP].merge (other[UP]);
 625   skylines_[DOWN].merge (other[DOWN]);
 626 }
 627
 628 Skyline&
 629 Skyline_pair::operator [] (Direction d)
 630 {
 631   return skylines_[d];
 632 }
 633
 634 Skyline const&
 635 Skyline_pair::operator [] (Direction d) const
 636 {
 637   return skylines_[d];
 638 }
 639
 640 /****************************************************************/
 641
 642
 643 IMPLEMENT_SIMPLE_SMOBS (Skyline);
 644 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
 645 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
 646
 647 IMPLEMENT_SIMPLE_SMOBS (Skyline_pair);
 648 IMPLEMENT_TYPE_P (Skyline_pair, "ly:skyline-pair?");
 649 IMPLEMENT_DEFAULT_EQUAL_P (Skyline_pair);
 650
 651 SCM
 652 Skyline::mark_smob (SCM)
 653 {
 654   ASSERT_LIVE_IS_ALLOWED();
 655   return SCM_EOL;
 656 }
 657
 658 int
 659 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
 660 {
 661   Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
 662   (void) r;
 663
 664   scm_puts ("#<Skyline>", port);
 665
 666   return 1;
 667 }
 668
 669 SCM
 670 Skyline_pair::mark_smob (SCM)
 671 {
 672   return SCM_EOL;
 673 }
 674
 675 int
 676 Skyline_pair::print_smob (SCM s, SCM port, scm_print_state *)
 677 {
 678   Skyline_pair *r = (Skyline_pair *) SCM_CELL_WORD_1 (s);
 679   (void) r;
 680
 681   scm_puts ("#<Skyline-pair>", port);
 682   return 1;
 683 }