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 /*
 422   build skyline from a set of boxes. If horizon_padding > 0, expand all the boxes
 423   by that amount and add 45-degree sloped boxes to the edges of each box (of
 424   width horizon_padding). That is, the total amount of horizontal expansion is
 425   horizon_padding*4, half of which is sloped and half of which is flat.
 426
 427   Boxes should have fatness in the horizon_axis (after they are expanded by
 428   horizon_padding), otherwise they are ignored.
 429  */
 430 Skyline::Skyline (vector<Box> const &boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 431 {
 432   list<Box> filtered_boxes;
 433   sky_ = sky;
 434
 435   Axis vert_axis = other_axis (horizon_axis);
 436   for (vsize i = 0; i < boxes.size (); i++)
 437     {
 438       Interval iv = boxes[i][horizon_axis];
 439       iv.widen (horizon_padding);
 440       if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
 441         filtered_boxes.push_front (boxes[i]);
 442     }
 443
 444   buildings_ = internal_build_skyline (&filtered_boxes, horizon_padding, horizon_axis, sky);
 445 }
 446
 447 Skyline::Skyline (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
 448 {
 449   sky_ = sky;
 450   Building front (b, horizon_padding, horizon_axis, sky);
 451   single_skyline (front, b[horizon_axis][LEFT] - horizon_padding,
 452                   horizon_padding, &buildings_);
 453 }
 454
 455 void
 456 Skyline::merge (Skyline const &other)
 457 {
 458   assert (sky_ == other.sky_);
 459
 460   list<Building> other_bld (other.buildings_);
 461   list<Building> my_bld;
 462   my_bld.splice (my_bld.begin (), buildings_);
 463   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 464 }
 465
 466 void
 467 Skyline::insert (Box const &b, Real horizon_padding, Axis a)
 468 {
 469   list<Building> other_bld;
 470   list<Building> my_bld;
 471
 472   if (isnan (b[other_axis (a)][LEFT])
 473       || isnan (b[other_axis (a)][RIGHT]))
 474     {
 475       programming_error ("insane box for skyline");
 476       return;
 477     }
 478
 479   /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
 480   Interval iv = b[a];
 481   iv.widen (horizon_padding);
 482   if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
 483     return;
 484
 485   my_bld.splice (my_bld.begin (), buildings_);
 486   single_skyline (Building (b, horizon_padding, a, sky_), b[a][LEFT] - horizon_padding,
 487                   horizon_padding, &other_bld);
 488   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 489 }
 490
 491 void
 492 Skyline::raise (Real r)
 493 {
 494   list<Building>::iterator end = buildings_.end ();
 495   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 496     i->y_intercept_ += sky_ * r;
 497 }
 498
 499 void
 500 Skyline::shift (Real s)
 501 {
 502   list<Building>::iterator end = buildings_.end ();
 503   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 504     {
 505       i->end_ += s;
 506       i->y_intercept_ -= s * i->slope_;
 507     }
 508 }
 509
 510 Real
 511 Skyline::distance (Skyline const &other, Real horizon_padding) const
 512 {
 513   assert (sky_ == -other.sky_);
 514
 515   Skyline const *padded_this = this;
 516   Skyline const *padded_other = &other;
 517
 518   /*
 519     For systems, padding is not added at creation time.  Padding is
 520     added to AxisGroup objects when outside-staff objects are added.
 521     Thus, when we want to place systems with horizontal padding,
 522     we do it at distance calculation time.
 523   */
 524   if (horizon_padding != 0.0)
 525     {
 526       padded_this = new Skyline (*padded_this, horizon_padding, X_AXIS);
 527       padded_other = new Skyline (*padded_other, horizon_padding, X_AXIS);
 528     }
 529
 530   list<Building>::const_iterator i = padded_this->buildings_.begin ();
 531   list<Building>::const_iterator j = padded_other->buildings_.begin ();
 532
 533   Real dist = -infinity_f;
 534   Real start = -infinity_f;
 535   while (i != padded_this->buildings_.end () && j != padded_other->buildings_.end ())
 536     {
 537       Real end = min (i->end_, j->end_);
 538       Real start_dist = i->height (start) + j->height (start);
 539       Real end_dist = i->height (end) + j->height (end);
 540       dist = max (dist, max (start_dist, end_dist));
 541       if (i->end_ <= j->end_)
 542         i++;
 543       else
 544         j++;
 545       start = end;
 546     }
 547   return dist;
 548 }
 549
 550 Real
 551 Skyline::height (Real airplane) const
 552 {
 553   assert (!isinf (airplane));
 554
 555   list<Building>::const_iterator i;
 556   for (i = buildings_.begin (); i != buildings_.end (); i++)
 557     {
 558       if (i->end_ >= airplane)
 559         return sky_ * i->height (airplane);
 560     }
 561
 562   assert (0);
 563   return 0;
 564 }
 565
 566 Real
 567 Skyline::max_height () const
 568 {
 569   Skyline s (-sky_);
 570   s.set_minimum_height (0);
 571   return sky_ * distance (s);
 572 }
 573
 574 void
 575 Skyline::set_minimum_height (Real h)
 576 {
 577   Skyline s (sky_);
 578   s.buildings_.front ().y_intercept_ = h * sky_;
 579   merge (s);
 580 }
 581
 582
 583 vector<Offset>
 584 Skyline::to_points (Axis horizon_axis) const
 585 {
 586   vector<Offset> out;
 587
 588   Real start = -infinity_f;
 589   for (list<Building>::const_iterator i (buildings_.begin ());
 590        i != buildings_.end (); i++)
 591     {
 592       out.push_back (Offset (start, sky_ * i->height (start)));
 593       out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
 594       start = i->end_;
 595     }
 596
 597   if (horizon_axis == Y_AXIS)
 598     for (vsize i = 0; i < out.size (); i++)
 599       out[i] = out[i].swapped ();
 600
 601   return out;
 602 }
 603
 604 bool
 605 Skyline::is_empty () const
 606 {
 607   Building b = buildings_.front ();
 608   return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
 609 }
 610
 611 void
 612 Skyline::clear ()
 613 {
 614   buildings_.clear ();
 615   empty_skyline (&buildings_);
 616 }
 617
 618 /****************************************************************/
 619
 620
 621 IMPLEMENT_SIMPLE_SMOBS (Skyline);
 622 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
 623 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
 624
 625 SCM
 626 Skyline::mark_smob (SCM)
 627 {
 628   ASSERT_LIVE_IS_ALLOWED ();
 629   return SCM_EOL;
 630 }
 631
 632 int
 633 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
 634 {
 635   Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
 636   (void) r;
 637
 638   scm_puts ("#<Skyline>", port);
 639
 640   return 1;
 641 }