source file of the GNU LilyPond music typesetter
- (c) 2006 Joe Neeman <joeneeman@gmail.com>
+ (c) 2006--2009 Joe Neeman <joeneeman@gmail.com>
*/
#include "skyline.hh"
+#include <deque>
+#include <cstdio>
#include "ly-smobs.icc"
but the distance routine does.
*/
-#define EPS 1e-10
+/* 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 inline bool
-approx_equal (Real x, Real y)
-{
- return abs (x - y) < EPS || (isinf (x) && isinf (y) && ((x > 0) == (y > 0)));
-}
-
-static inline bool
-approx_greater_than (Real x, Real y)
-{
- return x > y + EPS;
-}
-
-static inline bool
-approx_less_than (Real x, Real y)
-{
- return x < y - EPS;
-}
-
-static inline bool
-approx_less_equal (Real x, Real y)
-{
- return x <= y + EPS;
-}
-
-static inline bool
-approx_greater_equal (Real x, Real y)
+static void
+print_buildings (list<Building> const &b)
{
- return x >= y - EPS;
+ for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
+ i->print ();
}
void
Skyline::print () const
{
- for (list<Building>::const_iterator i = buildings_.begin ();
- i != buildings_.end (); i++)
- {
- (*i).print ();
- }
+ print_buildings (buildings_);
}
-bool
-Skyline::is_legal_skyline () const
+void
+Skyline::print_points () const
{
- list<Building>::const_iterator i;
- Real last_x = -infinity_f;
- for (i = buildings_.begin (); i != buildings_.end (); i++)
- {
- if (i->iv_[LEFT] != last_x)
- return false;
- last_x = i->iv_[RIGHT];
- if (isinf (i->iv_.length ()) && i->height_[LEFT] != i->height_[RIGHT])
- return false;
- }
- return last_x == infinity_f;
+ vector<Offset> ps (to_points (X_AXIS));
+
+ for (vsize i = 0; i < ps.size (); i++)
+ printf ("(%f,%f)%s" , ps[i][X_AXIS], ps[i][Y_AXIS],
+ (i%2)==1 ? "\n" : " ");
}
Building::Building (Real start, Real start_height, Real end_height, Real end)
- : iv_ (start, end)
{
- height_[LEFT] = start_height;
- height_[RIGHT] = end_height;
-
if (isinf (start) || isinf (end))
assert (start_height == end_height);
- precompute ();
+ end_ = end;
+ precompute (start, start_height, end_height, end);
}
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];
- iv_ = b[horizon_axis];
- iv_.widen (horizon_padding + EPS);
- height_[LEFT] = height;
- height_[RIGHT] = height;
-
- if (sane ())
- precompute ();
+ end_ = end;
+ precompute (start, height, height, end);
}
void
-Building::precompute ()
+Building::precompute (Real start, Real start_height, Real end_height, Real end)
{
- slope_ = (height_[RIGHT] - height_[LEFT]) / (iv_.length());
- if (height_[LEFT] == height_[RIGHT]) /* in case they're both infinity */
+ 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;
assert (!isinf (slope_) && !isnan (slope_));
- if (isinf (iv_[START]))
+ if (isinf (start))
{
- assert (slope_ == 0);
- y_intercept_ = height_[LEFT];
+ assert (start_height == end_height);
+ y_intercept_ = start_height;
}
else
- y_intercept_ = height_[LEFT] - slope_ * iv_[START];
+ y_intercept_ = start_height - slope_ * start;
}
Real
Building::height (Real x) const
{
- if (isinf (x))
- return (x > 0) ? height_[RIGHT] : height_[LEFT];
- return slope_*x + y_intercept_;
+ return isinf (x) ? y_intercept_ : slope_*x + y_intercept_;
}
void
Building::print () const
{
- printf ("X[%f,%f] -> Y[%f,%f]\n",
- iv_[LEFT], iv_[RIGHT],
- height_[LEFT], height_[RIGHT]);
+ printf ("%f x + %f ends at %f\n", slope_, y_intercept_, end_);
}
Real
-Building::intersection (Building const &other) const
+Building::intersection_x (Building const &other) const
{
- return (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
+ Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
+ return isnan (ret) ? -infinity_f : ret;
}
void
Building::leading_part (Real chop)
{
- assert (chop > iv_[LEFT] && chop <= iv_[RIGHT] && !approx_equal (chop, iv_[LEFT]));
- iv_[RIGHT] = chop;
- height_[RIGHT] = height (chop);
+ assert (chop <= end_);
+ end_ = chop;
}
Building
-Building::sloped_neighbour (Real horizon_padding, Direction d) const
+Building::sloped_neighbour (Real start, Real horizon_padding, Direction d) const
{
- Real left = iv_[d];
- Real right = iv_[d] + d * horizon_padding;
- Real left_height = height_[d];
- Real right_height = height_[d] - horizon_padding;
+ 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);
return Building (left, left_height, right_height, right);
}
-bool
-Building::sane () const
+static Real
+first_intersection (Building const &b, list<Building> *const s, Real start_x)
{
- return approx_less_than (iv_[LEFT], iv_[RIGHT])
- && !isinf (height_[RIGHT])
- && !isinf (height_[LEFT]);
-}
+ while (!s->empty () && start_x < b.end_)
+ {
+ Building c = s->front ();
+ if (c.conceals (b, start_x))
+ return start_x;
-static void
-skyline_trailing_part (list<Building> *sky, Real x)
-{
- if (approx_equal (x, sky->front ().iv_[RIGHT]))
- sky->pop_front ();
- else
- assert (x < sky->front ().iv_[RIGHT]);
+ Real i = b.intersection_x (c);
+ if (i > start_x && i <= b.end_ && i <= c.end_)
+ return i;
- if (!sky->empty ())
- {
- sky->front ().iv_[LEFT] = x;
- sky->front ().height_[LEFT] = sky->front ().height (x);
+ start_x = c.end_;
+ if (b.end_ > c.end_)
+ s->pop_front ();
}
+ return b.end_;
}
bool
-Building::conceals_beginning (Building const &other) const
+Building::conceals (Building const &other, Real x) const
{
- if (approx_equal (intersection (other), iv_[LEFT]) || approx_equal (height_[LEFT], other.height_[LEFT]))
- return slope_ > other.slope_;
- return height_[LEFT] > other.height_[LEFT];
-}
+ if (slope_ == other.slope_)
+ return y_intercept_ > other.y_intercept_;
-bool
-Building::conceals (Building const &other) const
-{
- assert (iv_[LEFT] <= other.iv_[LEFT]);
- return (iv_[RIGHT] >= other.iv_[RIGHT])
- && approx_greater_equal (height (other.iv_[LEFT]), other.height_[LEFT])
- && approx_greater_equal (height (other.iv_[RIGHT]), other.height_[RIGHT]);
+ /* 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
Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
list<Building> *const result)
{
+ if (s1->empty () || s2->empty ())
+ {
+ programming_error ("tried to merge an empty skyline");
+ return;
+ }
+
+ Real x = -infinity_f;
while (!s1->empty ())
{
- if (s2->front ().conceals_beginning (s1->front ()))
+ if (s2->front ().conceals (s1->front (), x))
swap (s1, s2);
Building b = s1->front ();
- while (!s2->empty () && b.conceals (s2->front ()))
- s2->pop_front ();
+ Real end = first_intersection (b, s2, x);
+
if (s2->empty ())
{
result->push_front (b);
break;
}
- /* s2 either intersects with b or it ends after b */
- Real end = infinity_f;
- Real s2_start_height = s2->front ().height_[LEFT];
- Real s2_end_height = s2->front ().height_[RIGHT];
- Real s1_start_height = b.height (s2->front ().iv_[LEFT]);
- Real s1_end_height = b.height (s2->front ().iv_[RIGHT]);
- if (approx_greater_than (s2_start_height, s1_start_height))
- end = s2->front ().iv_[LEFT];
- else if (approx_greater_than (s2_end_height, s1_end_height))
- end = b.intersection (s2->front ());
- end = min (end, b.iv_[RIGHT]);
-
- b.leading_part (end);
- result->push_front (b);
-
- skyline_trailing_part (s1, end);
- skyline_trailing_part (s2, end);
+ /* 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 ();
+
+ x = end;
}
result->reverse ();
}
}
static void
-single_skyline (Building b, Real horizon_padding, list<Building> *const ret)
+single_skyline (Building b, Real start, Real horizon_padding, list<Building> *const ret)
{
- b.iv_.widen (horizon_padding);
-
- if (!isinf (b.iv_[RIGHT]))
- ret->push_front (Building (b.iv_[RIGHT], -infinity_f,
+ 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 (horizon_padding > 0 && !isinf (b.iv_.length ()))
- ret->push_front (b.sloped_neighbour (horizon_padding, RIGHT));
+ if (sloped_neighbours)
+ ret->push_front (b.sloped_neighbour (start, horizon_padding, RIGHT));
- if (b.iv_[RIGHT] > b.iv_[LEFT])
+ if (b.end_ > start + EPS)
ret->push_front (b);
- if (horizon_padding > 0 && !isinf (b.iv_.length ()))
- ret->push_front (b.sloped_neighbour (horizon_padding, LEFT));
- if (!isinf (b.iv_[LEFT]))
+ 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, b.iv_[LEFT]));
+ -infinity_f, start - horizon_padding));
}
-void
-Skyline::internal_build_skyline (list<Building> *buildings, list<Building> *const result)
+/* 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];
+
+ if (iv[LEFT] - horizon_padding < last_end)
+ {
+ i++;
+ continue;
+ }
+
+ if (iv[LEFT] - horizon_padding > last_end + EPS)
+ result.push_front (Building (last_end, -infinity_f, -infinity_f, iv[LEFT] - 2*horizon_padding));
+
+ 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));
+
+ 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 = buildings->size ();
+ vsize size = boxes->size ();
if (size == 0)
{
- empty_skyline (result);
- return;
+ list<Building> result;
+ empty_skyline (&result);
+ return result;
}
else if (size == 1)
{
- single_skyline (buildings->front (), 0, result);
- return;
+ list<Building> result;
+ single_skyline (Building (boxes->front (), horizon_padding, horizon_axis, sky),
+ boxes->front ()[horizon_axis][LEFT], horizon_padding, &result);
+ return result;
}
- list<Building> right_half;
- list<Building>::iterator i = buildings->begin ();
-
- for (vsize s = 0; s < size/2; s++)
- i++;
- right_half.splice (right_half.end (), *buildings, i, buildings->end ());
+ deque<list<Building> > partials;
+ boxes->sort (LessThanBox (horizon_axis));
+ while (!boxes->empty ())
+ partials.push_back (non_overlapping_skyline (boxes, horizon_padding, horizon_axis, sky));
- list<Building> right;
- list<Building> left;
- internal_build_skyline (&right_half, &right);
- internal_build_skyline (buildings, &left);
- internal_merge_skyline (&right, &left, result);
+ /* 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_);
+ 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++)
{
*/
Skyline::Skyline (vector<Box> const &boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
{
- list<Building> bldgs;
+ list<Box> filtered_boxes;
sky_ = sky;
+ Axis vert_axis = other_axis (horizon_axis);
for (vsize i = 0; i < boxes.size (); i++)
{
- Building front (boxes[i], horizon_padding, horizon_axis, sky);
- if (front.sane ())
- {
- bldgs.push_front (front);
- if (horizon_padding > 0 && !isinf (front.iv_.length ()))
- {
- bldgs.push_front (front.sloped_neighbour (horizon_padding, LEFT));
- bldgs.push_front (front.sloped_neighbour (horizon_padding, RIGHT));
- }
- }
+ 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]);
}
- internal_build_skyline (&bldgs, &buildings_);
- assert (is_legal_skyline ());
+ 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, 0, horizon_axis, sky);
- single_skyline (front, horizon_padding, &buildings_);
+ Building front (b, horizon_padding, horizon_axis, sky);
+ single_skyline (front, b[horizon_axis][LEFT], horizon_padding, &buildings_);
}
void
list<Building> my_bld;
my_bld.splice (my_bld.begin (), buildings_);
internal_merge_skyline (&other_bld, &my_bld, &buildings_);
- assert (is_legal_skyline ());
}
void
list<Building> other_bld;
list<Building> my_bld;
+ if (isnan (b[other_axis (a)][LEFT])
+ || isnan (b[other_axis (a)][RIGHT]))
+ {
+ programming_error ("insane box for skyline");
+ return;
+ }
+
+ /* 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, 0, a, sky_), horizon_padding, &other_bld);
+ single_skyline (Building (b, horizon_padding, a, sky_), b[a][LEFT], horizon_padding, &other_bld);
internal_merge_skyline (&other_bld, &my_bld, &buildings_);
- assert (is_legal_skyline ());
}
void
{
list<Building>::iterator end = buildings_.end ();
for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
- {
- i->height_[LEFT] += sky_ * r;
- i->height_[RIGHT] += sky_ * r;
- i->y_intercept_ += sky_ * r;
- }
- assert (is_legal_skyline ());
+ i->y_intercept_ += sky_ * r;
}
void
-Skyline::shift (Real r)
+Skyline::shift (Real s)
{
list<Building>::iterator end = buildings_.end ();
for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
{
- i->iv_[LEFT] += r;
- i->iv_[RIGHT] += r;
+ i->end_ += s;
+ i->y_intercept_ -= s * i->slope_;
}
}
list<Building>::const_iterator j = other.buildings_.begin ();
Real dist = -infinity_f;
+ Real start = -infinity_f;
while (i != buildings_.end () && j != other.buildings_.end ())
{
- Interval iv = intersection (i->iv_, j->iv_);
- dist = max (dist, max (i->height (iv[LEFT]) + j->height (iv[LEFT]),
- i->height (iv[RIGHT]) + j->height (iv[RIGHT])));
- if (i->iv_[RIGHT] <= j->iv_[RIGHT])
+ 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;
}
list<Building>::const_iterator i;
for (i = buildings_.begin (); i != buildings_.end (); i++)
{
- if (i->iv_[RIGHT] >= airplane)
+ if (i->end_ >= airplane)
return sky_ * i->height (airplane);
}
Skyline::set_minimum_height (Real h)
{
Skyline s (sky_);
- s.buildings_.front ().height_[LEFT] = h * sky_;
- s.buildings_.front ().height_[RIGHT] = h * sky_;
s.buildings_.front ().y_intercept_ = h * sky_;
merge (s);
}
vector<Offset>
-Skyline::to_points () const
+Skyline::to_points (Axis horizon_axis) const
{
vector<Offset> out;
+ Real start = -infinity_f;
for (list<Building>::const_iterator i (buildings_.begin ());
i != buildings_.end (); i++)
{
- if (!isinf (i->iv_[LEFT]) && !isinf (i->height_[LEFT]))
- out.push_back (Offset (i->iv_[LEFT], sky_ * i->height_[LEFT]));
- if (!isinf (i->iv_[RIGHT]) && !isinf (i->height_[RIGHT]))
- out.push_back (Offset (i->iv_[RIGHT], sky_ * i->height_[RIGHT]));
+ out.push_back (Offset (start, sky_ * i->height (start)));
+ out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
+ start = i->end_;
}
- return out;
-}
-
-Skyline_pair::Skyline_pair ()
- : skylines_ (Skyline (DOWN), Skyline (UP))
-{
-}
-
-Skyline_pair::Skyline_pair (vector<Box> const &boxes, Real padding, Axis a)
- : skylines_ (Skyline (boxes, padding, a, DOWN), Skyline (boxes, padding, a, UP))
-{
-}
-
-Skyline_pair::Skyline_pair (Box const &b, Real padding, Axis a)
- : skylines_ (Skyline (b, padding, a, DOWN), Skyline (b, padding, a, UP))
-{
-}
-void
-Skyline_pair::raise (Real r)
-{
- skylines_[UP].raise (r);
- skylines_[DOWN].raise (r);
-}
-
-void
-Skyline_pair::shift (Real r)
-{
- skylines_[UP].shift (r);
- skylines_[DOWN].shift (r);
-}
+ if (horizon_axis == Y_AXIS)
+ for (vsize i = 0; i < out.size (); i++)
+ out[i] = out[i].swapped ();
-void
-Skyline_pair::insert (Box const &b, Real padding, Axis a)
-{
- skylines_[UP].insert (b, padding, a);
- skylines_[DOWN].insert (b, padding, a);
-}
-
-void
-Skyline_pair::merge (Skyline_pair const &other)
-{
- skylines_[UP].merge (other[UP]);
- skylines_[DOWN].merge (other[DOWN]);
+ return out;
}
-Skyline&
-Skyline_pair::operator [] (Direction d)
+bool
+Skyline::is_empty () const
{
- return skylines_[d];
+ Building b = buildings_.front ();
+ return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
}
-Skyline const&
-Skyline_pair::operator [] (Direction d) const
-{
- return skylines_[d];
-}
/****************************************************************/
IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
-IMPLEMENT_SIMPLE_SMOBS (Skyline_pair);
-IMPLEMENT_TYPE_P (Skyline_pair, "ly:skyline-pair?");
-IMPLEMENT_DEFAULT_EQUAL_P (Skyline_pair);
-
SCM
Skyline::mark_smob (SCM)
{
+ ASSERT_LIVE_IS_ALLOWED ();
return SCM_EOL;
}
return 1;
}
-
-SCM
-Skyline_pair::mark_smob (SCM)
-{
- return SCM_EOL;
-}
-
-int
-Skyline_pair::print_smob (SCM s, SCM port, scm_print_state *)
-{
- Skyline_pair *r = (Skyline_pair *) SCM_CELL_WORD_1 (s);
- (void) r;
-
- scm_puts ("#<Skyline-pair>", port);
- return 1;
-}