lily/skyline.cc

   1 /*
   2   This file is part of LilyPond, the GNU music typesetter.
   3
   4   Copyright (C) 2006--2011 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 <deque>
  22 #include <cstdio>
  23
  24 #include "ly-smobs.icc"
  25
  26 /* A skyline is a sequence of non-overlapping buildings: something like
  27    this:
  28                    _______
  29                   |       \                                 ________
  30                   |        \                       ________/        \
  31         /\        |          \                    /                  \
  32        /  --------             \                 /                    \
  33       /                          \              /                      \
  34      /                             ------------/                        ----
  35    --
  36    Each building has a starting position, and ending position, a starting
  37    height and an ending height.
  38
  39    The following invariants are observed:
  40     - the start of the first building is at -infinity
  41     - the end of the last building is at infinity
  42     - if a building has infinite length (ie. the first and last buildings),
  43       then its starting height and ending height are equal
  44     - the end of one building is the same as the beginning of the next
  45       building
  46
  47    We also allow skylines to point down (the structure is exactly the same,
  48    but we think of the part above the line as being filled with mass and the
  49    part below as being empty). ::distance finds the minimum distance between
  50    an UP skyline and a DOWN skyline.
  51
  52    Note that we store DOWN skylines upside-down. That is, in order to compare
  53    a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
  54    This means that the merging routine doesn't need to be aware of direction,
  55    but the distance routine does.
  56 */
  57
  58 /* If we start including very thin buildings, numerical accuracy errors can
  59    arise. Therefore, we ignore all buildings that are less than epsilon wide. */
  60 #define EPS 1e-5
  61
  62 static void
  63 print_buildings (list<Building> const &b)
  64 {
  65   for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
  66     i->print ();
  67 }
  68
  69 void
  70 Skyline::print () const
  71 {
  72   print_buildings (buildings_);
  73 }
  74
  75 void
  76 Skyline::print_points () const
  77 {
  78   vector<Offset> ps (to_points (X_AXIS));
  79
  80   for (vsize i = 0; i < ps.size (); i++)
  81     printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
  82             (i % 2) == 1 ? "\n" : " ");
  83 }
  84
  85 Building::Building (Real start, Real start_height, Real end_height, Real end)
  86 {
  87   if (isinf (start) || isinf (end))
  88     assert (start_height == end_height);
  89
  90   end_ = end;
  91   precompute (start, start_height, end_height, end);
  92 }
  93
  94 Building::Building (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
  95 {
  96   Real start = b[horizon_axis][LEFT] - horizon_padding;
  97   Real end = b[horizon_axis][RIGHT] + horizon_padding;
  98   Real height = sky * b[other_axis (horizon_axis)][sky];
  99
 100   end_ = end;
 101   precompute (start, height, height, end);
 102 }
 103
 104 void
 105 Building::precompute (Real start, Real start_height, Real end_height, Real end)
 106 {
 107   slope_ = (end_height - start_height) / (end - start);
 108   if (start_height == end_height) /* if they were both infinite, we would get nan, not 0, from the prev line */
 109     slope_ = 0;
 110
 111   assert (!isinf (slope_) && !isnan (slope_));
 112
 113   if (isinf (start))
 114     {
 115       assert (start_height == end_height);
 116       y_intercept_ = start_height;
 117     }
 118   else
 119     y_intercept_ = start_height - slope_ * start;
 120 }
 121
 122 Real
 123 Building::height (Real x) const
 124 {
 125   return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
 126 }
 127
 128 void
 129 Building::print () const
 130 {
 131   printf ("%f x + %f ends at %f\n", slope_, y_intercept_, end_);
 132 }
 133
 134 Real
 135 Building::intersection_x (Building const &other) const
 136 {
 137   Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
 138   return isnan (ret) ? -infinity_f : ret;
 139 }
 140
 141 void
 142 Building::leading_part (Real chop)
 143 {
 144   assert (chop <= end_);
 145   end_ = chop;
 146 }
 147
 148 Building
 149 Building::sloped_neighbour (Real start, Real horizon_padding, Direction d) const
 150 {
 151   Real x = (d == LEFT) ? start : end_;
 152   Real left = x;
 153   Real right = x + d * horizon_padding;
 154   Real left_height = height (x);
 155   Real right_height = left_height - horizon_padding;
 156   if (d == LEFT)
 157     {
 158       swap (left, right);
 159       swap (left_height, right_height);
 160     }
 161   return Building (left, left_height, right_height, right);
 162 }
 163
 164 static Real
 165 first_intersection (Building const &b, list<Building> *const s, Real start_x)
 166 {
 167   while (!s->empty () && start_x < b.end_)
 168     {
 169       Building c = s->front ();
 170       if (c.conceals (b, start_x))
 171         return start_x;
 172
 173       Real i = b.intersection_x (c);
 174       if (i > start_x && i <= b.end_ && i <= c.end_)
 175         return i;
 176
 177       start_x = c.end_;
 178       if (b.end_ > c.end_)
 179         s->pop_front ();
 180     }
 181   return b.end_;
 182 }
 183
 184 bool
 185 Building::conceals (Building const &other, Real x) const
 186 {
 187   if (slope_ == other.slope_)
 188     return y_intercept_ > other.y_intercept_;
 189
 190   /* their slopes were not equal, so there is an intersection point */
 191   Real i = intersection_x (other);
 192   return (i <= x && slope_ > other.slope_)
 193          || (i > x && slope_ < other.slope_);
 194 }
 195
 196 void
 197 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
 198                                  list<Building> *const result)
 199 {
 200   if (s1->empty () || s2->empty ())
 201     {
 202       programming_error ("tried to merge an empty skyline");
 203       return;
 204     }
 205
 206   Real x = -infinity_f;
 207   while (!s1->empty ())
 208     {
 209       if (s2->front ().conceals (s1->front (), x))
 210         swap (s1, s2);
 211
 212       Building b = s1->front ();
 213       Real end = first_intersection (b, s2, x);
 214
 215       if (s2->empty ())
 216         {
 217           result->push_front (b);
 218           break;
 219         }
 220
 221       /* only include buildings wider than epsilon */
 222       if (end > x + EPS)
 223         {
 224           b.leading_part (end);
 225           result->push_front (b);
 226         }
 227
 228       if (end >= s1->front ().end_)
 229         s1->pop_front ();
 230
 231       x = end;
 232     }
 233   result->reverse ();
 234 }
 235
 236 static void
 237 empty_skyline (list<Building> *const ret)
 238 {
 239   ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
 240 }
 241
 242 /*
 243   Given Building 'b' with starting wall location 'start', extend each side
 244   with a sloped roofline of width 'horizon_padding'; put the skyline in 'ret'
 245 */
 246 static void
 247 single_skyline (Building b, Real start, Real horizon_padding, list<Building> *const ret)
 248 {
 249   bool sloped_neighbours = horizon_padding > 0 && !isinf (start) && !isinf (b.end_);
 250   if (!isinf (b.end_))
 251     ret->push_front (Building (b.end_ + horizon_padding, -infinity_f,
 252                                -infinity_f, infinity_f));
 253   if (sloped_neighbours)
 254     ret->push_front (b.sloped_neighbour (start, horizon_padding, RIGHT));
 255
 256   if (b.end_ > start + EPS)
 257     ret->push_front (b);
 258
 259   if (sloped_neighbours)
 260     ret->push_front (b.sloped_neighbour (start, horizon_padding, LEFT));
 261
 262   if (!isinf (start))
 263     ret->push_front (Building (-infinity_f, -infinity_f,
 264                                -infinity_f, start - horizon_padding));
 265 }
 266
 267 /* remove a non-overlapping set of boxes from BOXES and build a skyline
 268    out of them */
 269 static list<Building>
 270 non_overlapping_skyline (list<Box> *const boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 271 {
 272   list<Building> result;
 273   Real last_end = -infinity_f;
 274   list<Box>::iterator i = boxes->begin ();
 275   while (i != boxes->end ())
 276     {
 277       Interval iv = (*i)[horizon_axis];
 278
 279       if (iv[LEFT] - horizon_padding < last_end)
 280         {
 281           i++;
 282           continue;
 283         }
 284
 285       if (iv[LEFT] - horizon_padding > last_end + EPS)
 286         result.push_front (Building (last_end, -infinity_f, -infinity_f, iv[LEFT] - 2 * horizon_padding));
 287
 288       Building b (*i, horizon_padding, horizon_axis, sky);
 289       bool sloped_neighbours = horizon_padding > 0 && !isinf (iv.length ());
 290       if (sloped_neighbours)
 291         result.push_front (b.sloped_neighbour (iv[LEFT] - horizon_padding, horizon_padding, LEFT));
 292       result.push_front (b);
 293       if (sloped_neighbours)
 294         result.push_front (b.sloped_neighbour (iv[LEFT] - horizon_padding, horizon_padding, RIGHT));
 295
 296       list<Box>::iterator j = i++;
 297       boxes->erase (j);
 298       last_end = result.front ().end_;
 299     }
 300   if (last_end < infinity_f)
 301     result.push_front (Building (last_end, -infinity_f, -infinity_f, infinity_f));
 302   result.reverse ();
 303   return result;
 304 }
 305
 306 class LessThanBox
 307 {
 308   Axis a_;
 309
 310 public:
 311   LessThanBox (Axis a)
 312   {
 313     a_ = a;
 314   }
 315
 316   bool operator () (Box const &b1, Box const &b2)
 317   {
 318     return b1[a_][LEFT] < b2[a_][LEFT];
 319   }
 320 };
 321
 322 list<Building>
 323 Skyline::internal_build_skyline (list<Box> *boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 324 {
 325   vsize size = boxes->size ();
 326
 327   if (size == 0)
 328     {
 329       list<Building> result;
 330       empty_skyline (&result);
 331       return result;
 332     }
 333   else if (size == 1)
 334     {
 335       list<Building> result;
 336       single_skyline (Building (boxes->front (), horizon_padding, horizon_axis, sky),
 337                       boxes->front ()[horizon_axis][LEFT] - horizon_padding,
 338                       horizon_padding, &result);
 339       return result;
 340     }
 341
 342   deque<list<Building> > partials;
 343   boxes->sort (LessThanBox (horizon_axis));
 344   while (!boxes->empty ())
 345     partials.push_back (non_overlapping_skyline (boxes, horizon_padding, horizon_axis, sky));
 346
 347   /* we'd like to say while (partials->size () > 1) but that's O (n).
 348      Instead, we exit in the middle of the loop */
 349   while (!partials.empty ())
 350     {
 351       list<Building> merged;
 352       list<Building> one = partials.front ();
 353       partials.pop_front ();
 354       if (partials.empty ())
 355         return one;
 356
 357       list<Building> two = partials.front ();
 358       partials.pop_front ();
 359       internal_merge_skyline (&one, &two, &merged);
 360       partials.push_back (merged);
 361     }
 362   assert (0);
 363   return list<Building> ();
 364 }
 365
 366 Skyline::Skyline ()
 367 {
 368   sky_ = UP;
 369   empty_skyline (&buildings_);
 370 }
 371
 372 Skyline::Skyline (Skyline const &src)
 373 {
 374   sky_ = src.sky_;
 375
 376   /* doesn't a list's copy constructor do this? -- jneem */
 377   for (list<Building>::const_iterator i = src.buildings_.begin ();
 378        i != src.buildings_.end (); i++)
 379     {
 380       buildings_.push_back (Building ((*i)));
 381     }
 382 }
 383
 384 Skyline::Skyline (Direction sky)
 385 {
 386   sky_ = sky;
 387   empty_skyline (&buildings_);
 388 }
 389
 390 /*
 391   build padded skyline from an existing skyline with padding
 392   added to it.
 393 */
 394
 395 Skyline::Skyline (Skyline const &src, Real horizon_padding, Axis a)
 396 {
 397   /*
 398      We extract boxes from the skyline, then build a new skyline from
 399      the boxes.
 400      A box is created for every horizontal portion of the skyline
 401      Because skylines are defined positive, and then inverted if they
 402      are to be down-facing, we create the new skyline in the UP
 403      direction, then give it the down direction if needed.
 404   */
 405   Real start = -infinity_f;
 406   list<Box> boxes;
 407
 408   // establish a baseline box
 409   boxes.push_back (Box (Interval (-infinity_f, infinity_f),
 410                         Interval (0, 0)));
 411   list<Building>::const_iterator end = src.buildings_.end ();
 412   for (list<Building>::const_iterator i = src.buildings_.begin (); i != end; start = i->end_, i++)
 413     if ((i->slope_ == 0) && !isinf (i->y_intercept_))
 414       boxes.push_back (Box (Interval (start, i->end_),
 415                             Interval (-infinity_f, i->y_intercept_)));
 416   buildings_ = internal_build_skyline (&boxes, horizon_padding, X_AXIS, UP);
 417   sky_ = src.sky_;
 418 }
 419
 420 /*
 421   build skyline from a set of boxes. If horizon_padding > 0, expand all the boxes
 422   by that amount and add 45-degree sloped boxes to the edges of each box (of
 423   width horizon_padding). That is, the total amount of horizontal expansion is
 424   horizon_padding*4, half of which is sloped and half of which is flat.
 425
 426   Boxes should have fatness in the horizon_axis (after they are expanded by
 427   horizon_padding), otherwise they are ignored.
 428  */
 429 Skyline::Skyline (vector<Box> const &boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 430 {
 431   list<Box> filtered_boxes;
 432   sky_ = sky;
 433
 434   Axis vert_axis = other_axis (horizon_axis);
 435   for (vsize i = 0; i < boxes.size (); i++)
 436     {
 437       Interval iv = boxes[i][horizon_axis];
 438       iv.widen (horizon_padding);
 439       if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
 440         filtered_boxes.push_front (boxes[i]);
 441     }
 442
 443   buildings_ = internal_build_skyline (&filtered_boxes, horizon_padding, horizon_axis, sky);
 444 }
 445
 446 Skyline::Skyline (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
 447 {
 448   sky_ = sky;
 449   Building front (b, horizon_padding, horizon_axis, sky);
 450   single_skyline (front, b[horizon_axis][LEFT] - horizon_padding,
 451                   horizon_padding, &buildings_);
 452 }
 453
 454 void
 455 Skyline::merge (Skyline const &other)
 456 {
 457   assert (sky_ == other.sky_);
 458
 459   list<Building> other_bld (other.buildings_);
 460   list<Building> my_bld;
 461   my_bld.splice (my_bld.begin (), buildings_);
 462   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 463 }
 464
 465 void
 466 Skyline::insert (Box const &b, Real horizon_padding, Axis a)
 467 {
 468   list<Building> other_bld;
 469   list<Building> my_bld;
 470
 471   if (isnan (b[other_axis (a)][LEFT])
 472       || isnan (b[other_axis (a)][RIGHT]))
 473     {
 474       programming_error ("insane box for skyline");
 475       return;
 476     }
 477
 478   /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
 479   Interval iv = b[a];
 480   iv.widen (horizon_padding);
 481   if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
 482     return;
 483
 484   my_bld.splice (my_bld.begin (), buildings_);
 485   single_skyline (Building (b, horizon_padding, a, sky_), b[a][LEFT] - horizon_padding,
 486                   horizon_padding, &other_bld);
 487   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 488 }
 489
 490 void
 491 Skyline::raise (Real r)
 492 {
 493   list<Building>::iterator end = buildings_.end ();
 494   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 495     i->y_intercept_ += sky_ * r;
 496 }
 497
 498 void
 499 Skyline::shift (Real s)
 500 {
 501   list<Building>::iterator end = buildings_.end ();
 502   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 503     {
 504       i->end_ += s;
 505       i->y_intercept_ -= s * i->slope_;
 506     }
 507 }
 508
 509 Real
 510 Skyline::distance (Skyline const &other, Real horizon_padding) const
 511 {
 512   assert (sky_ == -other.sky_);
 513
 514   Skyline const *padded_this = this;
 515   Skyline const *padded_other = &other;
 516
 517   /*
 518     For systems, padding is not added at creation time.  Padding is
 519     added to AxisGroup objects when outside-staff objects are added.
 520     Thus, when we want to place systems with horizontal padding,
 521     we do it at distance calculation time.
 522   */
 523   if (horizon_padding != 0.0)
 524     {
 525       padded_this = new Skyline (*padded_this, horizon_padding, X_AXIS);
 526       padded_other = new Skyline (*padded_other, horizon_padding, X_AXIS);
 527     }
 528
 529   list<Building>::const_iterator i = padded_this->buildings_.begin ();
 530   list<Building>::const_iterator j = padded_other->buildings_.begin ();
 531
 532   Real dist = -infinity_f;
 533   Real start = -infinity_f;
 534   while (i != padded_this->buildings_.end () && j != padded_other->buildings_.end ())
 535     {
 536       Real end = min (i->end_, j->end_);
 537       Real start_dist = i->height (start) + j->height (start);
 538       Real end_dist = i->height (end) + j->height (end);
 539       dist = max (dist, max (start_dist, end_dist));
 540       if (i->end_ <= j->end_)
 541         i++;
 542       else
 543         j++;
 544       start = end;
 545     }
 546   return dist;
 547 }
 548
 549 Real
 550 Skyline::height (Real airplane) const
 551 {
 552   assert (!isinf (airplane));
 553
 554   list<Building>::const_iterator i;
 555   for (i = buildings_.begin (); i != buildings_.end (); i++)
 556     {
 557       if (i->end_ >= airplane)
 558         return sky_ * i->height (airplane);
 559     }
 560
 561   assert (0);
 562   return 0;
 563 }
 564
 565 Real
 566 Skyline::max_height () const
 567 {
 568   Skyline s (-sky_);
 569   s.set_minimum_height (0);
 570   return sky_ * distance (s);
 571 }
 572
 573 void
 574 Skyline::set_minimum_height (Real h)
 575 {
 576   Skyline s (sky_);
 577   s.buildings_.front ().y_intercept_ = h * sky_;
 578   merge (s);
 579 }
 580
 581 vector<Offset>
 582 Skyline::to_points (Axis horizon_axis) const
 583 {
 584   vector<Offset> out;
 585
 586   Real start = -infinity_f;
 587   for (list<Building>::const_iterator i (buildings_.begin ());
 588        i != buildings_.end (); i++)
 589     {
 590       out.push_back (Offset (start, sky_ * i->height (start)));
 591       out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
 592       start = i->end_;
 593     }
 594
 595   if (horizon_axis == Y_AXIS)
 596     for (vsize i = 0; i < out.size (); i++)
 597       out[i] = out[i].swapped ();
 598
 599   return out;
 600 }
 601
 602 bool
 603 Skyline::is_empty () const
 604 {
 605   Building b = buildings_.front ();
 606   return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
 607 }
 608
 609 void
 610 Skyline::clear ()
 611 {
 612   buildings_.clear ();
 613   empty_skyline (&buildings_);
 614 }
 615
 616 /****************************************************************/
 617
 618 IMPLEMENT_SIMPLE_SMOBS (Skyline);
 619 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
 620 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
 621
 622 SCM
 623 Skyline::mark_smob (SCM)
 624 {
 625   ASSERT_LIVE_IS_ALLOWED ();
 626   return SCM_EOL;
 627 }
 628
 629 int
 630 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
 631 {
 632   Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
 633   (void) r;
 634
 635   scm_puts ("#<Skyline>", port);
 636
 637   return 1;
 638 }