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] - 2*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   if (isnan (b[other_axis (a)][LEFT])
 423       || isnan (b[other_axis (a)][RIGHT]))
 424     {
 425       programming_error ("insane box for skyline");
 426       return;
 427     }
 428
 429   /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
 430   Interval iv = b[a];
 431   iv.widen (horizon_padding);
 432   if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
 433     return;
 434
 435   my_bld.splice (my_bld.begin (), buildings_);
 436   single_skyline (Building (b, horizon_padding, a, sky_), b[a][LEFT], horizon_padding, &other_bld);
 437   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 438 }
 439
 440 void
 441 Skyline::raise (Real r)
 442 {
 443   list<Building>::iterator end = buildings_.end ();
 444   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 445     i->y_intercept_ += sky_ * r;
 446 }
 447
 448 void
 449 Skyline::shift (Real s)
 450 {
 451   list<Building>::iterator end = buildings_.end ();
 452   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 453     {
 454       i->end_ += s;
 455       i->y_intercept_ -= s * i->slope_;
 456     }
 457 }
 458
 459 Real
 460 Skyline::distance (Skyline const &other) const
 461 {
 462   assert (sky_ == -other.sky_);
 463   list<Building>::const_iterator i = buildings_.begin ();
 464   list<Building>::const_iterator j = other.buildings_.begin ();
 465
 466   Real dist = -infinity_f;
 467   Real start = -infinity_f;
 468   while (i != buildings_.end () && j != other.buildings_.end ())
 469     {
 470       Real end = min (i->end_, j->end_);
 471       Real start_dist = i->height (start) + j->height (start);
 472       Real end_dist = i->height (end) + j->height (end);
 473       dist = max (dist, max (start_dist, end_dist));
 474       if (i->end_ <= j->end_)
 475         i++;
 476       else
 477         j++;
 478       start = end;
 479     }
 480   return dist;
 481 }
 482
 483 Real
 484 Skyline::height (Real airplane) const
 485 {
 486   assert (!isinf (airplane));
 487
 488   list<Building>::const_iterator i;
 489   for (i = buildings_.begin (); i != buildings_.end (); i++)
 490     {
 491       if (i->end_ >= airplane)
 492         return sky_ * i->height (airplane);
 493     }
 494
 495   assert (0);
 496   return 0;
 497 }
 498
 499 Real
 500 Skyline::max_height () const
 501 {
 502   Skyline s (-sky_);
 503   s.set_minimum_height (0);
 504   return sky_ * distance (s);
 505 }
 506
 507 void
 508 Skyline::set_minimum_height (Real h)
 509 {
 510   Skyline s (sky_);
 511   s.buildings_.front ().y_intercept_ = h * sky_;
 512   merge (s);
 513 }
 514
 515
 516 vector<Offset>
 517 Skyline::to_points (Axis horizon_axis) const
 518 {
 519   vector<Offset> out;
 520
 521   Real start = -infinity_f;
 522   for (list<Building>::const_iterator i (buildings_.begin ());
 523        i != buildings_.end (); i++)
 524     {
 525       out.push_back (Offset (start, sky_ * i->height (start)));
 526       out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
 527       start = i->end_;
 528     }
 529
 530   if (horizon_axis == Y_AXIS)
 531     for (vsize i = 0; i < out.size (); i++)
 532       out[i] = out[i].swapped ();
 533
 534   return out;
 535 }
 536
 537 bool
 538 Skyline::is_empty () const
 539 {
 540   Building b = buildings_.front ();
 541   return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
 542 }
 543
 544
 545 /****************************************************************/
 546
 547
 548 IMPLEMENT_SIMPLE_SMOBS (Skyline);
 549 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
 550 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
 551
 552 SCM
 553 Skyline::mark_smob (SCM)
 554 {
 555   ASSERT_LIVE_IS_ALLOWED ();
 556   return SCM_EOL;
 557 }
 558
 559 int
 560 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
 561 {
 562   Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
 563   (void) r;
 564
 565   scm_puts ("#<Skyline>", port);
 566
 567   return 1;
 568 }