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