lily/skyline.cc

   1 /* skyline.cc -- implement the Skyline class
   2
   3    source file of the GNU LilyPond music typesetter
   4
   5    (c) 2006 Joe Neeman <joeneeman@gmail.com>
   6 */
   7
   8 #include "skyline.hh"
   9
  10 #include "ly-smobs.icc"
  11
  12 /* A skyline is a sequence of non-overlapping buildings: something like
  13    this:
  14                    _______
  15                   /       \                                ________
  16                  /         \                      ________/        \
  17         /\      /           \                    /                  \
  18        /  -----/              \                 /                    \
  19       /                         \              /                      \
  20      /                            ------------/                        ----
  21    --
  22    Each building has a starting position, and ending position, a starting
  23    height and an ending height.
  24
  25    The following invariants are observed:
  26     - the start of the first building is at -infinity
  27     - the end of the last building is at infinity
  28     - if a building has infinite length (ie. the first and last buildings),
  29       then its starting height and ending height are equal
  30     - the end of one building is the same as the beginning of the next
  31       building
  32
  33    We also allow skylines to point down (the structure is exactly the same,
  34    but we think of the part above the line as being filled with mass and the
  35    part below as being empty). ::distance finds the minimum distance between
  36    an UP skyline and a DOWN skyline.
  37
  38    Note that we store DOWN skylines upside-down. That is, in order to compare
  39    a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
  40    This means that the merging routine doesn't need to be aware of direction,
  41    but the distance routine does.
  42 */
  43
  44 #define EPS 1e-10
  45
  46 static inline bool
  47 equal (Real x, Real y)
  48 {
  49   return abs (x - y) < EPS || (isinf (x) && isinf (y) && ((x > 0) == (y > 0)));
  50 }
  51
  52 void
  53 Skyline::print () const
  54 {
  55   for (list<Building>::const_iterator i = buildings_.begin ();
  56        i != buildings_.end (); i++)
  57     {
  58       (*i).print ();
  59     }
  60 }
  61
  62 bool
  63 Skyline::is_legal_skyline () const
  64 {
  65   list<Building>::const_iterator i;
  66   Real last_x = -infinity_f;
  67   Real last_h = -infinity_f;
  68   for (i = buildings_.begin (); i != buildings_.end (); i++)
  69     {
  70       if (i->iv_[LEFT] != last_x)
  71         return false;
  72       if (i != buildings_.begin () && !equal (i->height_[LEFT], last_h))
  73         return false;
  74       last_x = i->iv_[RIGHT];
  75       last_h = i->height_[RIGHT];
  76     }
  77   return last_x == infinity_f;
  78 }
  79
  80 Building::Building (Real start, Real start_height, Real end_height,
  81                     Real end, Real max_slope)
  82   : iv_ (start, end)
  83 {
  84   height_[LEFT] = start_height;
  85   height_[RIGHT] = end_height;
  86
  87   if (isinf (start))
  88     assert (isinf (start_height) || start_height == end_height);
  89   if (isinf (end))
  90     assert (isinf (end_height) || start_height == end_height);
  91
  92   precompute (max_slope);
  93 }
  94
  95 void
  96 Building::precompute (Real max_slope)
  97 {
  98   slope_ = (height_[RIGHT] - height_[LEFT]) / (iv_.length());
  99   if (height_[LEFT] == height_[RIGHT])
 100     slope_ = 0;
 101   if (isinf (slope_) || isnan (slope_))
 102     slope_ = max_slope * (height_[LEFT] < height_[RIGHT] ? 1 : -1);
 103
 104 #if 0
 105   /*
 106     this check is sensitive to roundoff errors when converting to/from
 107     sequences of points.
 108    */
 109   assert (abs (slope_) <= max_slope + EPS);
 110 #endif
 111
 112   if (isinf (iv_[START]))
 113     {
 114       if (isinf (iv_[STOP]))
 115         zero_height_ = height_[LEFT];
 116       else
 117         zero_height_ = height_[RIGHT] - slope_ * iv_[STOP];
 118     }
 119   else
 120     zero_height_ = height_[LEFT] - slope_ * iv_[START];
 121 }
 122
 123 Real
 124 Building::height (Real x) const
 125 {
 126   if (isinf (x))
 127     return (x > 0) ? height_[RIGHT] : height_[LEFT];
 128   return slope_*x + zero_height_;
 129 }
 130
 131 void
 132 Building::print () const
 133 {
 134   printf ("X[%f,%f] -> Y[%f,%f]\n",
 135           iv_[LEFT], iv_[RIGHT],
 136           height_[LEFT], height_[RIGHT]);
 137 }
 138
 139 Real
 140 Building::intersection (Building const &other) const
 141 {
 142   return (zero_height_ - other.zero_height_) / (other.slope_ - slope_);
 143 }
 144
 145 void
 146 Building::leading_part (Real chop, Real h)
 147 {
 148   assert (chop > iv_[LEFT] && chop <= iv_[RIGHT] && !equal (chop, iv_[LEFT]));
 149   assert (equal (h, height (chop)));
 150   iv_[RIGHT] = chop;
 151   height_[RIGHT] = h;
 152 }
 153
 154 static void
 155 skyline_trailing_part (list<Building> *sky, Real x)
 156 {
 157   if (equal (x, sky->front ().iv_[RIGHT]))
 158     sky->pop_front ();
 159   else
 160     assert (x < sky->front ().iv_[RIGHT]);
 161
 162   if (!sky->empty ())
 163     {
 164       sky->front ().iv_[LEFT] = x;
 165       sky->front ().height_[LEFT] = sky->front ().height (x);
 166     }
 167 }
 168
 169 bool
 170 Building::obstructs (Building const &other) const
 171 {
 172   if (equal (intersection (other), iv_[LEFT]) || equal (height_[LEFT], other.height_[LEFT]))
 173     return slope_ > other.slope_ || (slope_ == other.slope_ && zero_height_ > other.zero_height_);
 174   return height_[LEFT] > other.height_[LEFT];
 175 }
 176
 177 void
 178 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
 179                                  list<Building> *const result)
 180 {
 181   while (!s1->empty ())
 182     {
 183       if (s2->front ().obstructs (s1->front ()))
 184         swap (s1, s2);
 185
 186       Building b = s1->front ();
 187       while (s2->front ().iv_[RIGHT] < b.iv_[RIGHT]
 188              && s2->front ().height_[RIGHT] <= b.height (s2->front ().iv_[RIGHT]) + EPS)
 189         s2->pop_front ();
 190
 191       /* the front of s2 either intersects with b or it ends after b */
 192       Real end = infinity_f;
 193       Real s2_end_height = s2->front ().height_[RIGHT];
 194       Real s1_end_height = b.height (s2->front ().iv_[RIGHT]);
 195       if (s2_end_height > s1_end_height + EPS)
 196         end = b.intersection (s2->front ());
 197       end = min (end, b.iv_[RIGHT]);
 198       Real height = b.height (end);
 199
 200       b.leading_part (end, height);
 201       result->push_front (b);
 202
 203       skyline_trailing_part (s1, end);
 204       if (!s1->empty ())
 205         s1->front ().height_[LEFT] = height;
 206       skyline_trailing_part (s2, end);
 207     }
 208   result->reverse ();
 209 }
 210
 211 static void
 212 empty_skyline (list<Building> *const ret)
 213 {
 214   ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f, 0));
 215 }
 216
 217 static void
 218 single_skyline (Building const &b, list<Building> *const ret, Real max_slope)
 219 {
 220   if (!isinf (b.iv_[RIGHT]))
 221     ret->push_front (Building (b.iv_[RIGHT], b.height_[RIGHT],
 222                                -infinity_f, infinity_f, max_slope));
 223   ret->push_front (b);
 224   if (!isinf (b.iv_[LEFT]))
 225     ret->push_front (Building (-infinity_f, -infinity_f,
 226                                b.height_[LEFT], b.iv_[LEFT], max_slope));
 227 }
 228
 229 void
 230 Skyline::internal_build_skyline (list<Building> *buildings, list<Building> *const result)
 231 {
 232   vsize size = buildings->size ();
 233
 234   if (size == 0)
 235     {
 236       empty_skyline (result);
 237       return;
 238     }
 239   else if (size == 1)
 240     {
 241       single_skyline (buildings->front (), result, max_slope_);
 242       return;
 243     }
 244
 245   list<Building> right_half;
 246   list<Building>::iterator i = buildings->begin ();
 247
 248   for (vsize s = 0; s < size/2; s++)
 249     i++;
 250   right_half.splice (right_half.end (), *buildings, i, buildings->end ());
 251
 252   list<Building> right;
 253   list<Building> left;
 254   internal_build_skyline (&right_half, &right);
 255   internal_build_skyline (buildings, &left);
 256   internal_merge_skyline (&right, &left, result);
 257 }
 258
 259 Skyline::Skyline ()
 260 {
 261   max_slope_ = 2;
 262   sky_ = UP;
 263   empty_skyline (&buildings_);
 264
 265
 266 }
 267
 268 Skyline::Skyline (Skyline const &src)
 269 {
 270   max_slope_ = src.max_slope_;
 271   sky_ = src.sky_;
 272
 273   for (list<Building>::const_iterator i = src.buildings_.begin ();
 274        i != src.buildings_.end (); i++)
 275     {
 276       buildings_.push_back (Building ((*i)));
 277     }
 278 }
 279
 280 Skyline::Skyline (Direction sky)
 281 {
 282   max_slope_ = 2;
 283   sky_ = sky;
 284   empty_skyline (&buildings_);
 285 }
 286
 287 /*
 288   build skyline from a set of boxes.
 289
 290   Boxes should have fatness in the horizon_axis, otherwise they are ignored.
 291  */
 292 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
 293 {
 294   list<Building> bldgs;
 295   sky_ = sky;
 296   max_slope_ = 2;
 297
 298   for (vsize i = 0; i < boxes.size (); i++)
 299     {
 300       Interval iv = boxes[i][horizon_axis];
 301       Real height = sky * boxes[i][other_axis (horizon_axis)][sky];
 302       if (!iv.is_empty () && !isinf (height) && !equal (iv[LEFT], iv[RIGHT]))
 303         bldgs.push_front (Building (iv[LEFT], height, height, iv[RIGHT],
 304                                     max_slope_));
 305     }
 306
 307   internal_build_skyline (&bldgs, &buildings_);
 308   assert (is_legal_skyline ());
 309 }
 310
 311 Skyline::Skyline (vector<Offset> const &points, Real max_slope, Direction sky)
 312 {
 313   sky_ = sky;
 314   max_slope_ = max_slope;
 315
 316   for (vsize i = 1; i < points.size (); i++)
 317     {
 318       buildings_.push_back (Building (points[i-1][X_AXIS], sky * points[i-1][Y_AXIS],
 319                                       sky * points[i][Y_AXIS],
 320
 321                                       points[i][X_AXIS],
 322                                       max_slope));
 323
 324     }
 325
 326   assert (is_legal_skyline ());
 327 }
 328
 329 void
 330 Skyline::merge (Skyline const &other)
 331 {
 332   assert (sky_ == other.sky_);
 333
 334   list<Building> other_bld (other.buildings_);
 335   list<Building> my_bld;
 336   my_bld.splice (my_bld.begin (), buildings_);
 337   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 338   assert (is_legal_skyline ());
 339 }
 340
 341 void
 342 Skyline::insert (Box const &b, Axis a)
 343 {
 344   list<Building> other_bld;
 345   list<Building> my_bld;
 346   Interval iv = b[a];
 347   Real height = sky_ * b[other_axis (a)][sky_];
 348
 349   assert (!iv.is_empty ());
 350
 351   my_bld.splice (my_bld.begin (), buildings_);
 352   single_skyline (Building (iv[LEFT], height, height, iv[RIGHT], max_slope_), &other_bld, max_slope_);
 353   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
 354   assert (is_legal_skyline ());
 355 }
 356
 357 void
 358 Skyline::raise (Real r)
 359 {
 360   list<Building>::iterator end = buildings_.end ();
 361   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
 362     {
 363       i->height_[LEFT] += sky_ * r;
 364       i->height_[RIGHT] += sky_ * r;
 365       i->zero_height_ += sky_ * r;
 366     }
 367   assert (is_legal_skyline ());
 368 }
 369
 370 Real
 371 Skyline::distance (Skyline const &other) const
 372 {
 373   assert (sky_ == -other.sky_);
 374   list<Building>::const_iterator i = buildings_.begin ();
 375   list<Building>::const_iterator j = other.buildings_.begin ();
 376
 377   Real dist = -infinity_f;
 378   while (i != buildings_.end () && j != other.buildings_.end ())
 379     {
 380       Interval iv = intersection (i->iv_, j->iv_);
 381       dist = max (dist, max (i->height (iv[LEFT]) + j->height (iv[LEFT]),
 382                              i->height (iv[RIGHT]) + j->height (iv[RIGHT])));
 383       if (i->iv_[RIGHT] <= j->iv_[RIGHT])
 384         i++;
 385       else
 386         j++;
 387     }
 388   return dist;
 389 }
 390
 391 Real
 392 Skyline::height (Real airplane) const
 393 {
 394   assert (!isinf (airplane));
 395
 396   list<Building>::const_iterator i;
 397   for (i = buildings_.begin (); i != buildings_.end (); i++)
 398     {
 399       if (i->iv_[RIGHT] >= airplane)
 400         return sky_ * i->height (airplane);
 401     }
 402
 403   assert (0);
 404   return 0;
 405 }
 406
 407 Real
 408 Skyline::max_height () const
 409 {
 410   Skyline s (-sky_);
 411   s.set_minimum_height (0);
 412   return sky_ * distance (s);
 413 }
 414
 415 void
 416 Skyline::set_minimum_height (Real h)
 417 {
 418   Skyline s (sky_);
 419   s.buildings_.front ().height_[LEFT] = h * sky_;
 420   s.buildings_.front ().height_[RIGHT] = h * sky_;
 421   s.buildings_.front ().zero_height_ = h * sky_;
 422   merge (s);
 423 }
 424
 425
 426 vector<Offset>
 427 Skyline::to_points () const
 428 {
 429   vector<Offset> out;
 430
 431   bool first = true;
 432   for (list<Building>::const_iterator i (buildings_.begin ());
 433        i != buildings_.end (); i++)
 434     {
 435       if (first)
 436         out.push_back (Offset ((*i).iv_[LEFT], sky_ * (*i).height_[LEFT]));
 437
 438       first = false;
 439       out.push_back (Offset ((*i).iv_[RIGHT], sky_ * (*i).height_[RIGHT]));
 440     }
 441
 442   return out;
 443 }
 444
 445 /****************************************************************/
 446
 447
 448 IMPLEMENT_SIMPLE_SMOBS (Skyline);
 449 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
 450 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
 451
 452 SCM
 453 Skyline::mark_smob (SCM)
 454 {
 455   return SCM_EOL;
 456 }
 457
 458 int
 459 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
 460 {
 461   Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
 462   (void) r;
 463
 464   scm_puts ("#<Skyline>", port);
 465
 466   return 1;
 467 }