]> git.donarmstrong.com Git - lilypond.git/blobdiff - lily/skyline.cc
Run `make grand-replace'.
[lilypond.git] / lily / skyline.cc
index 1d7782008dae7df5f0ccbd28184cafa569ff323b..6e64a2c96b6983ba4fac58ebb05026536c058c1e 100644 (file)
-/*   
-  skyline.cc -- implement Skyline_entry and funcs.
+/* skyline.cc -- implement the Skyline class
 
-  source file of the GNU LilyPond music typesetter
-
-  (c) 2002--2004 Han-Wen Nienhuys <hanwen@cs.uu.nl>
+   source file of the GNU LilyPond music typesetter
+   (c) 2006--2008 Joe Neeman <joeneeman@gmail.com>
 */
 
-#include "skyline.hh" 
+#include "skyline.hh"
+#include <deque>
+
+#include "ly-smobs.icc"
+
+/* A skyline is a sequence of non-overlapping buildings: something like
+   this:
+                   _______
+                  |       \                                 ________
+                  |        \                       ________/        \
+        /\        |          \                    /                  \
+       /  --------             \                 /                    \
+      /                          \              /                      \
+     /                             ------------/                        ----
+   --
+   Each building has a starting position, and ending position, a starting
+   height and an ending height.
+
+   The following invariants are observed:
+    - the start of the first building is at -infinity
+    - the end of the last building is at infinity
+    - if a building has infinite length (ie. the first and last buildings),
+      then its starting height and ending height are equal
+    - the end of one building is the same as the beginning of the next
+      building
+
+   We also allow skylines to point down (the structure is exactly the same,
+   but we think of the part above the line as being filled with mass and the
+   part below as being empty). ::distance finds the minimum distance between
+   an UP skyline and a DOWN skyline.
+
+   Note that we store DOWN skylines upside-down. That is, in order to compare
+   a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
+   This means that the merging routine doesn't need to be aware of direction,
+   but the distance routine does.
+*/
 
-/*
-  A skyline is a shape of the form:
+/* If we start including very thin buildings, numerical accuracy errors can
+   arise. Therefore, we ignore all buildings that are less than epsilon wide. */
+#define EPS 1e-5
 
+static void
+print_buildings (list<Building> const &b)
+{
+  for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
+    i->print ();
+}
 
-                  ----
-                  |  |
-         ---------|  |
-         |           |
-          |          |
-          |          |______
-  --------|                 |___
-  
+void
+Skyline::print () const
+{
+  print_buildings (buildings_);
+}
 
+void
+Skyline::print_points () const
+{
+  vector<Offset> ps (to_points (X_AXIS));
 
-  This file deals with building such skyline structure, and computing
-  the minimum distance between two opposing skylines.
-  
-   
-  Invariants for a skyline:
+  for (vsize i = 0; i < ps.size (); i++)
+    printf ("(%f,%f)%s" , ps[i][X_AXIS], ps[i][Y_AXIS],
+           (i%2)==1 ? "\n" : " ");
+}
 
-  skyline[...].width_ forms a partition of the real interval, where
-  the segments are adjacent, and ascending. Hence we have
-  
-  skyline.top ().width_[RIGHT] = inf
-  skyline[0].width_[LEFT] = -inf
-  
- */
+Building::Building (Real start, Real start_height, Real end_height, Real end)
+{
+  if (isinf (start) || isinf (end))
+    assert (start_height == end_height);
 
+  end_ = end;
+  precompute (start, start_height, end_height, end);
+}
 
-const Real EPS = 1e-12;  
+Building::Building (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
+{
+  Real start = b[horizon_axis][LEFT] - horizon_padding;
+  Real end = b[horizon_axis][RIGHT] + horizon_padding;
+  Real height = sky * b[other_axis (horizon_axis)][sky];
 
-/*
-  TODO: avoid unnecessary fragmentation.
+  end_ = end;
+  precompute (start, height, height, end);
+}
 
-  This is O (n^2), searching and insertion.  Could be O (n log n) with
-  binsearch.
-*/
 void
-insert_extent_into_skyline (Array<Skyline_entry> *line, Box b, Axis line_axis,
-                           Direction d)
+Building::precompute (Real start, Real start_height, Real end_height, Real end)
 {
-  Interval extent = b[line_axis];
-  if (extent.is_empty ())
-    return;
-  
-  Real stick_out = b[other_axis (line_axis)][d];
+  slope_ = (end_height - start_height) / (end - start);
+  if (start_height == end_height) /* if they were both infinite, we would get nan, not 0, from the prev line */
+    slope_ = 0;
 
-  /*
-    Intersect each segment of LINE with EXTENT, and if non-empty, insert relevant segments. 
-   */
-  for (int i = line->size (); i--;)
+  assert (!isinf (slope_) && !isnan (slope_));
+
+  if (isinf (start))
     {
-      Interval w = line->elem (i).width_;
-      w.intersect (extent);
+      assert (start_height == end_height);
+      y_intercept_ = start_height;
+    }
+  else
+    y_intercept_ = start_height - slope_ * start;
+}
 
-      if (extent[LEFT] >= w[RIGHT])
-       break;
-      
-      Real my_height = line->elem (i).height_;
+Real 
+Building::height (Real x) const
+{
+  return isinf (x) ? y_intercept_ : slope_*x + y_intercept_;
+}
 
-      if (!w.is_empty () &&
-         w.length () > EPS
-         && d* (my_height - stick_out) < 0)
-       {
-         Interval e1 (line->elem (i).width_[LEFT], extent[LEFT]);
-         Interval e3 (extent[RIGHT], line->elem (i).width_[RIGHT]);
+void
+Building::print () const
+{
+  printf ("%f x + %f ends at %f\n", slope_, y_intercept_, end_);
+}
 
-         if (!e3.is_empty () && e3.length () > EPS)
-           line->insert (Skyline_entry (e3, my_height), i+1);
+Real
+Building::intersection_x (Building const &other) const
+{
+  Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
+  return isnan (ret) ? -infinity_f : ret;
+}
 
-         line->elem_ref (i).height_ = stick_out;
-         line->elem_ref (i).width_ = w;
-         if (!e1.is_empty () && e1.length () > EPS)
-           line->insert (Skyline_entry (e1, my_height), i );
-       }
+void
+Building::leading_part (Real chop)
+{
+  assert (chop <= end_);
+  end_ = chop;
+}
 
+Building
+Building::sloped_neighbour (Real start, Real horizon_padding, Direction d) const
+{
+  Real x = (d == LEFT) ? start : end_;
+  Real left = x;
+  Real right = x + d * horizon_padding;
+  Real left_height = height (x);
+  Real right_height = left_height - horizon_padding;
+  if (d == LEFT)
+    {
+      swap (left, right);
+      swap (left_height, right_height);
+    }
+  return Building (left, left_height, right_height, right);
+}
+
+static Real
+first_intersection (Building const &b, list<Building> *const s, Real start_x)
+{
+  while (!s->empty () && start_x < b.end_)
+    {
+      Building c = s->front ();
+      if (c.conceals (b, start_x))
+       return start_x;
+
+      Real i = b.intersection_x (c);
+      if (i > start_x && i <= b.end_ && i <= c.end_)
+       return i;
 
+      start_x = c.end_;
+      if (b.end_ > c.end_)
+       s->pop_front ();
     }
+  return b.end_;
+}
+
+bool
+Building::conceals (Building const &other, Real x) const
+{
+  if (slope_ == other.slope_)
+    return y_intercept_ > other.y_intercept_;
+
+  /* their slopes were not equal, so there is an intersection point */
+  Real i = intersection_x (other);
+  return (i <= x && slope_ > other.slope_)
+    || (i > x && slope_ < other.slope_);
 }
 
 void
-merge_skyline (Array<Skyline_entry> * a1,
-              Array<Skyline_entry> const  & a2,
-              Direction dir)
+Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
+                                list<Building> *const result)
 {
-  for (int i = 0; i < a2.size (); i++)
+  if (s1->empty () || s2->empty ())
     {
-      Box b;
-      b[X_AXIS] = a2[i].width_;
-      b[Y_AXIS][dir] = a2[i].height_;
-      b[Y_AXIS][-dir] = dir * infinity_f ;
+      programming_error ("tried to merge an empty skyline");
+      return;
+    }
+
+  Real x = -infinity_f;
+  while (!s1->empty ())
+    {
+      if (s2->front ().conceals (s1->front (), x))
+       swap (s1, s2);
+
+      Building b = s1->front ();
+      Real end = first_intersection (b, s2, x);
+
+      if (s2->empty ())
+       {
+         result->push_front (b);
+         break;
+       }
+
+      /* only include buildings wider than epsilon */
+      if (end > x + EPS)
+       {
+         b.leading_part (end);
+         result->push_front (b);
+       }
+
+      if (end >= s1->front ().end_)
+       s1->pop_front ();
 
-      insert_extent_into_skyline (a1, b, X_AXIS, dir);
+      x = end;
     }
+  result->reverse ();
 }
 
+static void
+empty_skyline (list<Building> *const ret)
+{
+  ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
+}
 
-Array<Skyline_entry>
-empty_skyline (Direction d)
+static void
+single_skyline (Building b, Real start, Real horizon_padding, list<Building> *const ret)
 {
-  Array<Skyline_entry> skyline;
+  bool sloped_neighbours = horizon_padding > 0 && !isinf (start) && !isinf (b.end_);
+  if (!isinf (b.end_))
+    ret->push_front (Building (b.end_ + horizon_padding, -infinity_f,
+                              -infinity_f, infinity_f));
+  if (sloped_neighbours)
+    ret->push_front (b.sloped_neighbour (start, horizon_padding, RIGHT));
+  
+  if (b.end_ > start + EPS)
+    ret->push_front (b);
 
-  Interval i;
-  i.set_empty ();
-  i.swap ();
-  Skyline_entry e;
-  e.width_ = i;
-  e.height_ = -d * infinity_f; 
-  skyline.push (e);
-  return skyline;
+  if (sloped_neighbours)
+    ret->push_front (b.sloped_neighbour (start, horizon_padding, LEFT));
+
+  if (!isinf (start))
+    ret->push_front (Building (-infinity_f, -infinity_f,
+                              -infinity_f, start - horizon_padding));
 }
 
-Array<Skyline_entry>
-extents_to_skyline (Array<Box> const &extents, Axis a, Direction d)
+/* remove a non-overlapping set of boxes from BOXES and build a skyline
+   out of them */
+static list<Building>
+non_overlapping_skyline (list<Box> *const boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 {
+  list<Building> result;
+  Real last_end = -infinity_f;
+  list<Box>::iterator i = boxes->begin ();
+  while (i != boxes->end ())
+    {
+      Interval iv = (*i)[horizon_axis];
 
-  Array<Skyline_entry> skyline = empty_skyline (d);
+      if (iv[LEFT] - horizon_padding < last_end)
+       {
+         i++;
+         continue;
+       }
 
-  /*
-    This makes a cubic algorithm (array  insertion is O (n),
-    searching the array dumbly is O (n), and for n items, we get O (n^3).)
+      if (iv[LEFT] - horizon_padding > last_end + EPS)
+       result.push_front (Building (last_end, -infinity_f, -infinity_f, iv[LEFT] - 2*horizon_padding));
 
-    We could do a lot better (n log (n), using a balanced tree) but
-    that seems overkill for now.
-   */
-  for (int j = extents.size (); j--; )
-    insert_extent_into_skyline (&skyline, extents[j], a, d);
+      Building b (*i, horizon_padding, horizon_axis, sky);
+      bool sloped_neighbours = horizon_padding > 0 && !isinf (iv.length ());
+      if (sloped_neighbours)
+       result.push_front (b.sloped_neighbour (iv[LEFT] - horizon_padding, horizon_padding, LEFT));
+      result.push_front (b);
+      if (sloped_neighbours)
+       result.push_front (b.sloped_neighbour (iv[LEFT] - horizon_padding, horizon_padding, RIGHT));
 
-  return skyline;
+      list<Box>::iterator j = i++;
+      boxes->erase (j);
+      last_end = result.front ().end_;
+    }
+  if (last_end < infinity_f)
+    result.push_front (Building (last_end, -infinity_f, -infinity_f, infinity_f));
+  result.reverse ();
+  return result;
 }
 
+class LessThanBox
+{
+  Axis a_;
+
+public:
+  LessThanBox (Axis a)
+  {
+    a_ = a;
+  }
+
+  bool operator() (Box const &b1, Box const &b2)
+  {
+    return b1[a_][LEFT] < b2[a_][LEFT];
+  }
+};
+
+list<Building>
+Skyline::internal_build_skyline (list<Box> *boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
+{
+  vsize size = boxes->size ();
+
+  if (size == 0)
+    {
+      list<Building> result;
+      empty_skyline (&result);
+      return result;
+    }
+  else if (size == 1)
+    {
+      list<Building> result;
+      single_skyline (Building (boxes->front (), horizon_padding, horizon_axis, sky),
+                     boxes->front ()[horizon_axis][LEFT], horizon_axis, &result);
+      return result;
+    }
 
+  deque<list<Building> > partials;
+  boxes->sort (LessThanBox (horizon_axis));
+  while (!boxes->empty ())
+    partials.push_back (non_overlapping_skyline (boxes, horizon_padding, horizon_axis, sky));
+
+  /* we'd like to say while (partials->size () > 1) but that's O (n).
+     Instead, we exit in the middle of the loop */
+  while (!partials.empty ())
+    {
+      list<Building> merged;
+      list<Building> one = partials.front ();
+      partials.pop_front ();
+      if (partials.empty ())
+       return one;
+
+      list<Building> two = partials.front ();
+      partials.pop_front ();
+      internal_merge_skyline (&one, &two, &merged);
+      partials.push_back (merged);
+    }
+  assert (0);
+  return list<Building> ();
+}
+
+Skyline::Skyline ()
+{
+  sky_ = UP;
+  empty_skyline (&buildings_);
+}
+
+Skyline::Skyline (Skyline const &src)
+{
+  sky_ = src.sky_;
+  /* doesn't a list's copy constructor do this? -- jneem */
+  for (list<Building>::const_iterator i = src.buildings_.begin ();
+       i != src.buildings_.end (); i++)
+    {
+      buildings_.push_back (Building ((*i)));
+    }
+}
+
+Skyline::Skyline (Direction sky)
+{
+  sky_ = sky;
+  empty_skyline (&buildings_);
+}
 
 /*
-  minimum distance that can be achieved between baselines. "Clouds" is
-  a skyline pointing down.
+  build skyline from a set of boxes. If horizon_padding > 0, expand all the boxes
+  by that amount and add 45-degree sloped boxes to the edges of each box (of
+  width horizon_padding). That is, the total amount of horizontal expansion is
+  horizon_padding*4, half of which is sloped and half of which is flat.
 
-  This is an O (n) algorithm.
+  Boxes should have fatness in the horizon_axis (after they are expanded by
+  horizon_padding), otherwise they are ignored.
  */
-Real
-skyline_meshing_distance (Array<Skyline_entry> const &buildings,
-                         Array<Skyline_entry> const &clouds)
+Skyline::Skyline (vector<Box> const &boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
 {
-  int i = buildings.size () -1;
-  int j  = clouds.size () -1;
+  list<Box> filtered_boxes;
+  sky_ = sky;
 
-  Real distance = - infinity_f;
+  Axis vert_axis = other_axis (horizon_axis);
+  for (vsize i = 0; i < boxes.size (); i++)
+    {
+      Interval iv = boxes[i][horizon_axis];
+      iv.widen (horizon_padding);
+      if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
+       filtered_boxes.push_front (boxes[i]);
+    }
   
-  while (i > 0 || j > 0)
+  buildings_ = internal_build_skyline (&filtered_boxes, horizon_padding, horizon_axis, sky);
+}
+
+Skyline::Skyline (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
+{
+  sky_ = sky;
+  Building front (b, horizon_padding, horizon_axis, sky);
+  single_skyline (front, b[horizon_axis][LEFT], horizon_padding, &buildings_);
+}
+
+void
+Skyline::merge (Skyline const &other)
+{
+  assert (sky_ == other.sky_);
+
+  list<Building> other_bld (other.buildings_);
+  list<Building> my_bld;
+  my_bld.splice (my_bld.begin (), buildings_);
+  internal_merge_skyline (&other_bld, &my_bld, &buildings_);
+}
+
+void
+Skyline::insert (Box const &b, Real horizon_padding, Axis a)
+{
+  list<Building> other_bld;
+  list<Building> my_bld;
+
+  if (isnan (b[other_axis (a)][LEFT])
+      || isnan (b[other_axis (a)][RIGHT]))
     {
-      Interval w = buildings[i].width_;
-      w.intersect (clouds[j].width_);
-      
-      if (!w.is_empty ())
-       distance = distance >? (buildings[i].height_ - clouds[j].height_);
+      programming_error ("insane box for skyline");
+      return;
+    }
 
-      if (i>0 && buildings[i].width_[LEFT] >=  clouds[j].width_[LEFT])
-       {
-         i--;
-       }
-      else if (j > 0 && buildings[i].width_[LEFT] <=  clouds[j].width_[LEFT])
-       {
-         j--;
-       }       
+  /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
+  Interval iv = b[a];
+  iv.widen (horizon_padding);
+  if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
+    return;
+
+  my_bld.splice (my_bld.begin (), buildings_);
+  single_skyline (Building (b, horizon_padding, a, sky_), b[a][LEFT], horizon_padding, &other_bld);
+  internal_merge_skyline (&other_bld, &my_bld, &buildings_);
+}
+
+void
+Skyline::raise (Real r)
+{
+  list<Building>::iterator end = buildings_.end ();
+  for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
+    i->y_intercept_ += sky_ * r;
+}
+
+void
+Skyline::shift (Real s)
+{
+  list<Building>::iterator end = buildings_.end ();
+  for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
+    {
+      i->end_ += s;
+      i->y_intercept_ -= s * i->slope_;
     }
+}
+
+Real
+Skyline::distance (Skyline const &other) const
+{
+  assert (sky_ == -other.sky_);
+  list<Building>::const_iterator i = buildings_.begin ();
+  list<Building>::const_iterator j = other.buildings_.begin ();
 
-  return distance;
+  Real dist = -infinity_f;
+  Real start = -infinity_f;
+  while (i != buildings_.end () && j != other.buildings_.end ())
+    {
+      Real end = min (i->end_, j->end_);
+      Real start_dist = i->height (start) + j->height (start);
+      Real end_dist = i->height (end) + j->height (end);
+      dist = max (dist, max (start_dist, end_dist));
+      if (i->end_ <= j->end_)
+       i++;
+      else
+       j++;
+      start = end;
+    }
+  return dist;
 }
 
-Skyline_entry::Skyline_entry ()
+Real
+Skyline::height (Real airplane) const
 {
-  height_ = 0.0;
+  assert (!isinf (airplane));
+
+  list<Building>::const_iterator i;
+  for (i = buildings_.begin (); i != buildings_.end (); i++)
+    {
+      if (i->end_ >= airplane)
+       return sky_ * i->height (airplane);
+    }
+
+  assert (0);
+  return 0;
 }
 
-Skyline_entry::Skyline_entry (Interval i, Real r)
+Real
+Skyline::max_height () const
 {
-  width_ = i;
-  height_ = r;
-  
+  Skyline s (-sky_);
+  s.set_minimum_height (0);
+  return sky_ * distance (s);
 }
 
 void
-heighten_skyline (Array<Skyline_entry> *buildings, Real ground)
+Skyline::set_minimum_height (Real h)
+{
+  Skyline s (sky_);
+  s.buildings_.front ().y_intercept_ = h * sky_;
+  merge (s);
+}
+
+
+vector<Offset>
+Skyline::to_points (Axis horizon_axis) const
 {
-  for (int i = 0; i < buildings->size (); i++)
-    buildings->elem_ref (i).height_ += ground; 
+  vector<Offset> out;
+
+  Real start = -infinity_f;
+  for (list<Building>::const_iterator i (buildings_.begin ());
+       i != buildings_.end (); i++)
+    {
+      out.push_back (Offset (start, sky_ * i->height (start)));
+      out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
+      start = i->end_;
+    }
+
+  if (horizon_axis == Y_AXIS)
+    for (vsize i = 0; i < out.size (); i++)
+      out[i] = out[i].swapped ();
+
+  return out;
+}
+
+bool
+Skyline::is_empty () const
+{
+  Building b = buildings_.front ();
+  return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
+}
+
+
+/****************************************************************/
+
+
+IMPLEMENT_SIMPLE_SMOBS (Skyline);
+IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
+IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
+
+SCM
+Skyline::mark_smob (SCM)
+{
+  ASSERT_LIVE_IS_ALLOWED ();
+  return SCM_EOL;
+}
+
+int
+Skyline::print_smob (SCM s, SCM port, scm_print_state *)
+{
+  Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
+  (void) r;
+
+  scm_puts ("#<Skyline>", port);
+
+  return 1;
 }