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