-/*
- skyline.cc -- implement Skyline_entry and funcs.
+/* skyline.cc -- implement the Skyline class
- source file of the GNU LilyPond music typesetter
-
- (c) 2002--2005 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;
}