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 static void
 243 single_skyline (Building b, Real start, Real horizon_padding, list<Building> *const ret)
 244 {
 245   bool sloped_neighbours = horizon_padding > 0 && !isinf (start) && !isinf (b.end_);
 246   if (!isinf (b.end_))
 247     ret->push_front (Building (b.end_ + horizon_padding, -infinity_f,
 248                                -infinity_f, infinity_f));
 249   if (sloped_neighbours)
 250     ret->push_front (b.sloped_neighbour (start, horizon_padding, RIGHT));
 251
 252   if (b.end_ > start + EPS)
 253     ret->push_front (b);
 254
 255   if (sloped_neighbours)
 256     ret->push_front (b.sloped_neighbour (start, horizon_padding, LEFT));
 257
 258   if (!isinf (start))
 259     ret->push_front (Building (-infinity_f, -infinity_f,
 260                                -infinity_f, start - horizon_padding));
 261 }
 262
 263 /* remove a non-overlapping set of boxes from BOXES and build a skyline
 264    out of them */
 265 static list<Building>
 266 non_overlapping_skyline (list<Box> *const boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 267 {
 268   list<Building> result;
 269   Real last_end = -infinity_f;
 270   list<Box>::iterator i = boxes->begin ();
 271   while (i != boxes->end ())
 272     {
 273       Interval iv = (*i)[horizon_axis];
 274
 275       if (iv[LEFT] - horizon_padding < last_end)
 276         {
 277           i++;
 278           continue;
 279         }
 280
 281       if (iv[LEFT] - horizon_padding > last_end + EPS)
 282         result.push_front (Building (last_end, -infinity_f, -infinity_f, iv[LEFT] - 2*horizon_padding));
 283
 284       Building b (*i, horizon_padding, horizon_axis, sky);
 285       bool sloped_neighbours = horizon_padding > 0 && !isinf (iv.length ());
 286       if (sloped_neighbours)
 287         result.push_front (b.sloped_neighbour (iv[LEFT] - horizon_padding, horizon_padding, LEFT));
 288       result.push_front (b);
 289       if (sloped_neighbours)
 290         result.push_front (b.sloped_neighbour (iv[LEFT] - horizon_padding, horizon_padding, RIGHT));
 291
 292       list<Box>::iterator j = i++;
 293       boxes->erase (j);
 294       last_end = result.front ().end_;
 295     }
 296   if (last_end < infinity_f)
 297     result.push_front (Building (last_end, -infinity_f, -infinity_f, infinity_f));
 298   result.reverse ();
 299   return result;
 300 }
 301
 302 class LessThanBox
 303 {
 304   Axis a_;
 305
 306 public:
 307   LessThanBox (Axis a)
 308   {
 309     a_ = a;
 310   }
 311
 312   bool operator() (Box const &b1, Box const &b2)
 313   {
 314     return b1[a_][LEFT] < b2[a_][LEFT];
 315   }
 316 };
 317
 318 list<Building>
 319 Skyline::internal_build_skyline (list<Box> *boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 320 {
 321   vsize size = boxes->size ();
 322
 323   if (size == 0)
 324     {
 325       list<Building> result;
 326       empty_skyline (&result);
 327       return result;
 328     }
 329   else if (size == 1)
 330     {
 331       list<Building> result;
 332       single_skyline (Building (boxes->front (), horizon_padding, horizon_axis, sky),
 333                       boxes->front ()[horizon_axis][LEFT], horizon_padding, &result);
 334       return result;
 335     }
 336
 337   deque<list<Building> > partials;
 338   boxes->sort (LessThanBox (horizon_axis));
 339   while (!boxes->empty ())
 340     partials.push_back (non_overlapping_skyline (boxes, horizon_padding, horizon_axis, sky));
 341
 342   /* we'd like to say while (partials->size () > 1) but that's O (n).
 343      Instead, we exit in the middle of the loop */
 344   while (!partials.empty ())
 345     {
 346       list<Building> merged;
 347       list<Building> one = partials.front ();
 348       partials.pop_front ();
 349       if (partials.empty ())
 350         return one;
 351
 352       list<Building> two = partials.front ();
 353       partials.pop_front ();
 354       internal_merge_skyline (&one, &two, &merged);
 355       partials.push_back (merged);
 356     }
 357   assert (0);
 358   return list<Building> ();
 359 }
 360
 361 Skyline::Skyline ()
 362 {
 363   sky_ = UP;
 364   empty_skyline (&buildings_);
 365 }
 366
 367 Skyline::Skyline (Skyline const &src)
 368 {
 369   sky_ = src.sky_;
 370
 371   /* doesn't a list's copy constructor do this? -- jneem */
 372   for (list<Building>::const_iterator i = src.buildings_.begin ();
 373        i != src.buildings_.end (); i++)
 374     {
 375       buildings_.push_back (Building ((*i)));
 376     }
 377 }
 378
 379 Skyline::Skyline (Direction sky)
 380 {
 381   sky_ = sky;
 382   empty_skyline (&buildings_);
 383 }
 384
 385 /*
 386   build padded skyline from an existing skyline with padding
 387   added to it.
 388 */
 389
 390 Skyline::Skyline (Skyline const &src, Real horizon_padding, Axis a)
 391 {
 392   /*
 393      We extract boxes from the skyline, then build a new skyline from
 394      the boxes.
 395      A box is created for every horizontal portion of the skyline
 396      Because skylines are defined positive, and then inverted if they
 397      are to be down-facing, we create the new skyline in the UP
 398      direction, then give it the down direction if needed.
 399   */
 400   Real start = -infinity_f;
 401   list<Box> boxes;
 402
 403   // establish a baseline box
 404   boxes.push_back (Box (Interval (-infinity_f, infinity_f),
 405                         Interval (0, 0)));
 406   list<Building>::const_iterator end = src.buildings_.end ();
 407   for (list<Building>::const_iterator i = src.buildings_.begin (); i != end; start=i->end_, i++ )
 408     if ((i->slope_ == 0) && !isinf (i->y_intercept_))
 409       boxes.push_back (Box (Interval (start, i->end_),
 410                             Interval (-infinity_f , i->y_intercept_)));
 411   buildings_ = internal_build_skyline (&boxes, horizon_padding, X_AXIS, UP);
 412   sky_ = src.sky_;
 413 }
 414
 415
 416 /*
 417   build skyline from a set of boxes. If horizon_padding > 0, expand all the boxes
 418   by that amount and add 45-degree sloped boxes to the edges of each box (of
 419   width horizon_padding). That is, the total amount of horizontal expansion is
 420   horizon_padding*4, half of which is sloped and half of which is flat.
 421
 422   Boxes should have fatness in the horizon_axis (after they are expanded by
 423   horizon_padding), otherwise they are ignored.
 424  */
 425 Skyline::Skyline (vector<Box> const &boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 426 {
 427   list<Box> filtered_boxes;
 428   sky_ = sky;
 429
 430   Axis vert_axis = other_axis (horizon_axis);
 431   for (vsize i = 0; i < boxes.size (); i++)
 432     {
 433       Interval iv = boxes[i][horizon_axis];
 434       iv.widen (horizon_padding);
 435       if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
 436         filtered_boxes.push_front (boxes[i]);
 437     }
 438
 439   buildings_ = internal_build_skyline (&filtered_boxes, horizon_padding, horizon_axis, sky);
 440 }
 441
 442 Skyline::Skyline (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
 443 {
 444   sky_ = sky;
 445   Building front (b, horizon_padding, horizon_axis, sky);
 446   single_skyline (front, b[horizon_axis][LEFT], horizon_padding, &buildings_);
 447 }
 448
 449 void
 450 Skyline::merge (Skyline const &other)
 451 {
 452   assert (sky_ == other.sky_);
 453
 454   list<Building> other_bld (other.buildings_);
 455   list<Building> my_bld;
 456   my_bld.splice (my_bld.begin (), buildings_);
 457   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 458 }
 459
 460 void
 461 Skyline::insert (Box const &b, Real horizon_padding, Axis a)
 462 {
 463   list<Building> other_bld;
 464   list<Building> my_bld;
 465
 466   if (isnan (b[other_axis (a)][LEFT])
 467       || isnan (b[other_axis (a)][RIGHT]))
 468     {
 469       programming_error ("insane box for skyline");
 470       return;
 471     }
 472
 473   /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
 474   Interval iv = b[a];
 475   iv.widen (horizon_padding);
 476   if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
 477     return;
 478
 479   my_bld.splice (my_bld.begin (), buildings_);
 480   single_skyline (Building (b, horizon_padding, a, sky_), b[a][LEFT], horizon_padding, &other_bld);
 481   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 482 }
 483
 484 void
 485 Skyline::raise (Real r)
 486 {
 487   list<Building>::iterator end = buildings_.end ();
 488   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 489     i->y_intercept_ += sky_ * r;
 490 }
 491
 492 void
 493 Skyline::shift (Real s)
 494 {
 495   list<Building>::iterator end = buildings_.end ();
 496   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 497     {
 498       i->end_ += s;
 499       i->y_intercept_ -= s * i->slope_;
 500     }
 501 }
 502
 503 Real
 504 Skyline::distance (Skyline const &other, Real horizon_padding) const
 505 {
 506   assert (sky_ == -other.sky_);
 507
 508   Skyline const *padded_this = this;
 509   Skyline const *padded_other = &other;
 510
 511   /*
 512     For systems, padding is not added at creation time.  Padding is
 513     added to AxisGroup objects when outside-staff objects are added.
 514     Thus, when we want to place systems with horizontal padding,
 515     we do it at distance calculation time.
 516   */
 517   if (horizon_padding != 0.0)
 518     {
 519       padded_this = new Skyline (*padded_this, horizon_padding, X_AXIS);
 520       padded_other = new Skyline (*padded_other, horizon_padding, X_AXIS);
 521     }
 522
 523   list<Building>::const_iterator i = padded_this->buildings_.begin ();
 524   list<Building>::const_iterator j = padded_other->buildings_.begin ();
 525
 526   Real dist = -infinity_f;
 527   Real start = -infinity_f;
 528   while (i != padded_this->buildings_.end () && j != padded_other->buildings_.end ())
 529     {
 530       Real end = min (i->end_, j->end_);
 531       Real start_dist = i->height (start) + j->height (start);
 532       Real end_dist = i->height (end) + j->height (end);
 533       dist = max (dist, max (start_dist, end_dist));
 534       if (i->end_ <= j->end_)
 535         i++;
 536       else
 537         j++;
 538       start = end;
 539     }
 540   return dist;
 541 }
 542
 543 Real
 544 Skyline::height (Real airplane) const
 545 {
 546   assert (!isinf (airplane));
 547
 548   list<Building>::const_iterator i;
 549   for (i = buildings_.begin (); i != buildings_.end (); i++)
 550     {
 551       if (i->end_ >= airplane)
 552         return sky_ * i->height (airplane);
 553     }
 554
 555   assert (0);
 556   return 0;
 557 }
 558
 559 Real
 560 Skyline::max_height () const
 561 {
 562   Skyline s (-sky_);
 563   s.set_minimum_height (0);
 564   return sky_ * distance (s);
 565 }
 566
 567 void
 568 Skyline::set_minimum_height (Real h)
 569 {
 570   Skyline s (sky_);
 571   s.buildings_.front ().y_intercept_ = h * sky_;
 572   merge (s);
 573 }
 574
 575
 576 vector<Offset>
 577 Skyline::to_points (Axis horizon_axis) const
 578 {
 579   vector<Offset> out;
 580
 581   Real start = -infinity_f;
 582   for (list<Building>::const_iterator i (buildings_.begin ());
 583        i != buildings_.end (); i++)
 584     {
 585       out.push_back (Offset (start, sky_ * i->height (start)));
 586       out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
 587       start = i->end_;
 588     }
 589
 590   if (horizon_axis == Y_AXIS)
 591     for (vsize i = 0; i < out.size (); i++)
 592       out[i] = out[i].swapped ();
 593
 594   return out;
 595 }
 596
 597 bool
 598 Skyline::is_empty () const
 599 {
 600   Building b = buildings_.front ();
 601   return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
 602 }
 603
 604 void
 605 Skyline::clear ()
 606 {
 607   buildings_.clear ();
 608   empty_skyline (&buildings_);
 609 }
 610
 611 /****************************************************************/
 612
 613
 614 IMPLEMENT_SIMPLE_SMOBS (Skyline);
 615 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
 616 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
 617
 618 SCM
 619 Skyline::mark_smob (SCM)
 620 {
 621   ASSERT_LIVE_IS_ALLOWED ();
 622   return SCM_EOL;
 623 }
 624
 625 int
 626 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
 627 {
 628   Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
 629   (void) r;
 630
 631   scm_puts ("#<Skyline>", port);
 632
 633   return 1;
 634 }