lily/skyline.cc

   1 /*
   2   This file is part of LilyPond, the GNU music typesetter.
   3
   4   Copyright (C) 2006--2010 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 skyline from a set of boxes. If horizon_padding > 0, expand all the boxes
 387   by that amount and add 45-degree sloped boxes to the edges of each box (of
 388   width horizon_padding). That is, the total amount of horizontal expansion is
 389   horizon_padding*4, half of which is sloped and half of which is flat.
 390
 391   Boxes should have fatness in the horizon_axis (after they are expanded by
 392   horizon_padding), otherwise they are ignored.
 393  */
 394 Skyline::Skyline (vector<Box> const &boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 395 {
 396   list<Box> filtered_boxes;
 397   sky_ = sky;
 398
 399   Axis vert_axis = other_axis (horizon_axis);
 400   for (vsize i = 0; i < boxes.size (); i++)
 401     {
 402       Interval iv = boxes[i][horizon_axis];
 403       iv.widen (horizon_padding);
 404       if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
 405         filtered_boxes.push_front (boxes[i]);
 406     }
 407
 408   buildings_ = internal_build_skyline (&filtered_boxes, horizon_padding, horizon_axis, sky);
 409 }
 410
 411 Skyline::Skyline (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
 412 {
 413   sky_ = sky;
 414   Building front (b, horizon_padding, horizon_axis, sky);
 415   single_skyline (front, b[horizon_axis][LEFT], horizon_padding, &buildings_);
 416 }
 417
 418 void
 419 Skyline::merge (Skyline const &other)
 420 {
 421   assert (sky_ == other.sky_);
 422
 423   list<Building> other_bld (other.buildings_);
 424   list<Building> my_bld;
 425   my_bld.splice (my_bld.begin (), buildings_);
 426   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 427 }
 428
 429 void
 430 Skyline::insert (Box const &b, Real horizon_padding, Axis a)
 431 {
 432   list<Building> other_bld;
 433   list<Building> my_bld;
 434
 435   if (isnan (b[other_axis (a)][LEFT])
 436       || isnan (b[other_axis (a)][RIGHT]))
 437     {
 438       programming_error ("insane box for skyline");
 439       return;
 440     }
 441
 442   /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
 443   Interval iv = b[a];
 444   iv.widen (horizon_padding);
 445   if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
 446     return;
 447
 448   my_bld.splice (my_bld.begin (), buildings_);
 449   single_skyline (Building (b, horizon_padding, a, sky_), b[a][LEFT], horizon_padding, &other_bld);
 450   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 451 }
 452
 453 void
 454 Skyline::raise (Real r)
 455 {
 456   list<Building>::iterator end = buildings_.end ();
 457   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 458     i->y_intercept_ += sky_ * r;
 459 }
 460
 461 void
 462 Skyline::shift (Real s)
 463 {
 464   list<Building>::iterator end = buildings_.end ();
 465   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 466     {
 467       i->end_ += s;
 468       i->y_intercept_ -= s * i->slope_;
 469     }
 470 }
 471
 472 Real
 473 Skyline::distance (Skyline const &other) const
 474 {
 475   assert (sky_ == -other.sky_);
 476   list<Building>::const_iterator i = buildings_.begin ();
 477   list<Building>::const_iterator j = other.buildings_.begin ();
 478
 479   Real dist = -infinity_f;
 480   Real start = -infinity_f;
 481   while (i != buildings_.end () && j != other.buildings_.end ())
 482     {
 483       Real end = min (i->end_, j->end_);
 484       Real start_dist = i->height (start) + j->height (start);
 485       Real end_dist = i->height (end) + j->height (end);
 486       dist = max (dist, max (start_dist, end_dist));
 487       if (i->end_ <= j->end_)
 488         i++;
 489       else
 490         j++;
 491       start = end;
 492     }
 493   return dist;
 494 }
 495
 496 Real
 497 Skyline::height (Real airplane) const
 498 {
 499   assert (!isinf (airplane));
 500
 501   list<Building>::const_iterator i;
 502   for (i = buildings_.begin (); i != buildings_.end (); i++)
 503     {
 504       if (i->end_ >= airplane)
 505         return sky_ * i->height (airplane);
 506     }
 507
 508   assert (0);
 509   return 0;
 510 }
 511
 512 Real
 513 Skyline::max_height () const
 514 {
 515   Skyline s (-sky_);
 516   s.set_minimum_height (0);
 517   return sky_ * distance (s);
 518 }
 519
 520 void
 521 Skyline::set_minimum_height (Real h)
 522 {
 523   Skyline s (sky_);
 524   s.buildings_.front ().y_intercept_ = h * sky_;
 525   merge (s);
 526 }
 527
 528
 529 vector<Offset>
 530 Skyline::to_points (Axis horizon_axis) const
 531 {
 532   vector<Offset> out;
 533
 534   Real start = -infinity_f;
 535   for (list<Building>::const_iterator i (buildings_.begin ());
 536        i != buildings_.end (); i++)
 537     {
 538       out.push_back (Offset (start, sky_ * i->height (start)));
 539       out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
 540       start = i->end_;
 541     }
 542
 543   if (horizon_axis == Y_AXIS)
 544     for (vsize i = 0; i < out.size (); i++)
 545       out[i] = out[i].swapped ();
 546
 547   return out;
 548 }
 549
 550 bool
 551 Skyline::is_empty () const
 552 {
 553   Building b = buildings_.front ();
 554   return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
 555 }
 556
 557 void
 558 Skyline::clear ()
 559 {
 560   buildings_.clear ();
 561   empty_skyline (&buildings_);
 562 }
 563
 564 /****************************************************************/
 565
 566
 567 IMPLEMENT_SIMPLE_SMOBS (Skyline);
 568 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
 569 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
 570
 571 SCM
 572 Skyline::mark_smob (SCM)
 573 {
 574   ASSERT_LIVE_IS_ALLOWED ();
 575   return SCM_EOL;
 576 }
 577
 578 int
 579 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
 580 {
 581   Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
 582   (void) r;
 583
 584   scm_puts ("#<Skyline>", port);
 585
 586   return 1;
 587 }