lily/skyline.cc

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