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