2 This file is part of LilyPond, the GNU music typesetter.
4 Copyright (C) 2006--2015 Joe Neeman <joeneeman@gmail.com>
6 LilyPond is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
11 LilyPond is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with LilyPond. If not, see <http://www.gnu.org/licenses/>.
21 #include "skyline-pair.hh"
29 /* A skyline is a sequence of non-overlapping buildings: something like
39 Each building has a starting position, and ending position, a starting
40 height and an ending height.
42 The following invariants are observed:
43 - the start of the first building is at -infinity
44 - the end of the last building is at infinity
45 - if a building has infinite length (ie. the first and last buildings),
46 then its starting height and ending height are equal
47 - the end of one building is the same as the beginning of the next
50 We also allow skylines to point down (the structure is exactly the same,
51 but we think of the part above the line as being filled with mass and the
52 part below as being empty). ::distance finds the minimum distance between
53 an UP skyline and a DOWN skyline.
55 Note that we store DOWN skylines upside-down. That is, in order to compare
56 a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
57 This means that the merging routine doesn't need to be aware of direction,
58 but the distance routine does.
60 From 2007 through 2012, buildings of width less than EPS were discarded,
61 citing numerical accuracy concerns. We remember that floating point
62 comparisons of nearly-equal values can be affected by rounding error.
63 Also, some target machines use the x87 floating point unit, which provides
64 extended precision for intermediate results held in registers. On this type
65 of hardware comparisons such as
66 double c = 1.0/3.0; boolean compare = (c == 1.0/3.0)
67 could go either way because the 1.0/3.0 is allowed to be kept
68 higher precision than the variable 'c'.
69 Alert to these considerations, we now accept buildings of zero-width.
73 print_buildings (list<Building> const &b)
75 for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
80 Skyline::print () const
82 print_buildings (buildings_);
86 Skyline::print_points () const
88 vector<Offset> ps (to_points (X_AXIS));
90 for (vsize i = 0; i < ps.size (); i++)
91 printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
92 (i % 2) == 1 ? "\n" : " ");
95 Building::Building (Real start, Real start_height, Real end_height, Real end)
97 if (isinf (start) || isinf (end))
98 assert (start_height == end_height);
102 precompute (start, start_height, end_height, end);
105 Building::Building (Box const &b, Axis horizon_axis, Direction sky)
107 Real start = b[horizon_axis][LEFT];
108 Real end = b[horizon_axis][RIGHT];
109 Real height = sky * b[other_axis (horizon_axis)][sky];
113 precompute (start, height, height, end);
117 Building::precompute (Real start, Real start_height, Real end_height, Real end)
119 slope_ = 0.0; /* if they were both infinite, we would get nan, not 0, from the prev line */
120 if (start_height != end_height)
121 slope_ = (end_height - start_height) / (end - start);
123 assert (!isinf (slope_) && !isnan (slope_));
127 assert (start_height == end_height);
128 y_intercept_ = start_height;
130 else if (fabs (slope_) > 1e6)
131 // too steep to be stored in slope-intercept form, given round-off error
134 y_intercept_ = std::max (start_height, end_height);
137 y_intercept_ = start_height - slope_ * start;
141 Building::height (Real x) const
143 return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
147 Building::print () const
149 printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
153 Building::intersection_x (Building const &other) const
155 Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
156 return isnan (ret) ? -infinity_f : ret;
159 // Returns a shift s such that (x + s, y) intersects the roof of
160 // this building. If no such shift exists, returns infinity_f.
162 Building::shift_to_intersect (Real x, Real y) const
164 // Solve for s: y = (x + s)*m + b
165 Real ret = (y - y_intercept_ - slope_ * x) / slope_;
167 if (ret >= start_ && ret <= end_ && !isinf (ret))
173 Building::above (Building const &other, Real x) const
175 return (isinf (y_intercept_) || isinf (other.y_intercept_) || isinf (x))
176 ? y_intercept_ > other.y_intercept_
177 : (slope_ - other.slope_) * x + y_intercept_ > other.y_intercept_;
180 // Remove redundant empty buildings from the skyline.
181 // If there are two adjacent empty buildings, they can be
184 Skyline::normalize ()
186 bool last_empty = false;
187 list<Building>::iterator i;
189 for (i = buildings_.begin (); i != buildings_.end (); i++)
191 if (last_empty && i->y_intercept_ == -infinity_f)
193 list<Building>::iterator last = i;
195 last->end_ = i->end_;
196 buildings_.erase (i);
199 last_empty = (i->y_intercept_ == -infinity_f);
202 assert (buildings_.front ().start_ == -infinity_f);
203 assert (buildings_.back ().end_ == infinity_f);
207 Skyline::internal_merge_skyline (list<Building> *sb, list<Building> *sc,
208 list<Building> *const result) const
210 if (sb->empty () || sc->empty ())
212 programming_error ("tried to merge an empty skyline");
216 Building b = sb->front ();
217 for (; !sc->empty (); sc->pop_front ())
219 /* Building b is continuing from the previous pass through the loop.
220 Building c is newly-considered, and starts no earlier than b started.
221 The comments draw b as if its roof had zero slope ----.
222 with dashes where b lies above c.
223 The roof of c could rise / or fall \ through the roof of b,
224 or the vertical sides | of c could intersect the roof of b. */
225 Building c = sc->front ();
226 if (b.end_ < c.end_) /* finish with b */
228 if (b.end_ <= b.start_) /* we are already finished with b */
230 else if (c.above (b, c.start_)) /* -| . | */
234 if (m.end_ > m.start_)
235 result->push_back (m);
236 if (b.above (c, b.end_)) /* -|\--. */
239 n.end_ = b.start_ = b.intersection_x (c);
240 result->push_back (n);
241 result->push_back (b);
246 if (c.above (b, b.end_)) /* ---/ . | */
247 b.end_ = b.intersection_x (c);
250 result->push_back (b);
252 /* 'c' continues further, so move it into 'b' for the next pass. */
256 else /* b.end_ > c.end_ so finish with c */
258 if (c.above (b, c.start_)) /* -| |---. */
262 if (m.end_ > m.start_)
263 result->push_back (m);
264 if (b.above (c, c.end_)) /* -| \---. */
265 c.end_ = b.intersection_x (c);
267 else if (c.above (b, c.end_)) /* ---/|--. */
270 c.start_ = m.end_ = b.intersection_x (c);
271 result->push_back (m);
273 else /* c is completely hidden by b */
275 result->push_back (c);
279 if (b.end_ > b.start_)
280 result->push_back (b);
284 empty_skyline (list<Building> *const ret)
286 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
290 Given Building 'b', build a skyline containing only that building.
293 single_skyline (Building b, list<Building> *const ret)
295 assert (b.end_ >= b.start_);
297 if (b.start_ != -infinity_f)
298 ret->push_back (Building (-infinity_f, -infinity_f,
299 -infinity_f, b.start_));
301 if (b.end_ != infinity_f)
302 ret->push_back (Building (b.end_, -infinity_f,
303 -infinity_f, infinity_f));
306 /* remove a non-overlapping set of boxes from BOXES and build a skyline
308 static list<Building>
309 non_overlapping_skyline (list<Building> *const buildings)
311 list<Building> result;
312 Real last_end = -infinity_f;
313 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
314 list<Building>::iterator i = buildings->begin ();
315 while (i != buildings->end ())
318 Real y1 = i->height (i->start_);
320 Real y2 = i->height (i->end_);
322 // Drop buildings that will obviously have no effect.
323 if (last_building.height (x1) >= y1
324 && last_building.end_ >= x2
325 && last_building.height (x2) >= y2)
327 list<Building>::iterator j = i++;
328 buildings->erase (j);
338 // Insert empty Buildings into any gaps. (TODO: is this needed? -KOH)
340 result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
342 result.push_back (*i);
346 list<Building>::iterator j = i++;
347 buildings->erase (j);
350 if (last_end < infinity_f)
351 result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
355 class LessThanBuilding
358 bool operator () (Building const &b1, Building const &b2)
360 return b1.start_ < b2.start_
361 || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
366 BUILDINGS is a list of buildings, but they could be overlapping
367 and in any order. The returned list of buildings is ordered and non-overlapping.
370 Skyline::internal_build_skyline (list<Building> *buildings) const
372 vsize size = buildings->size ();
376 list<Building> result;
377 empty_skyline (&result);
382 list<Building> result;
383 single_skyline (buildings->front (), &result);
387 deque<list<Building> > partials;
388 buildings->sort (LessThanBuilding ());
389 while (!buildings->empty ())
390 partials.push_back (non_overlapping_skyline (buildings));
392 /* we'd like to say while (partials->size () > 1) but that's O (n).
393 Instead, we exit in the middle of the loop */
394 while (!partials.empty ())
396 list<Building> merged;
397 list<Building> one = partials.front ();
398 partials.pop_front ();
399 if (partials.empty ())
402 list<Building> two = partials.front ();
403 partials.pop_front ();
404 internal_merge_skyline (&one, &two, &merged);
405 partials.push_back (merged);
408 return list<Building> ();
414 empty_skyline (&buildings_);
417 Skyline::Skyline (Direction sky)
420 empty_skyline (&buildings_);
424 Build skyline from a set of boxes.
426 Boxes should be non-empty on both axes. Otherwise, they will be ignored
428 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
430 list<Building> buildings;
433 for (vsize i = 0; i < boxes.size (); i++)
434 if (!boxes[i].is_empty (X_AXIS)
435 && !boxes[i].is_empty (Y_AXIS))
436 buildings.push_front (Building (boxes[i], horizon_axis, sky));
438 buildings_ = internal_build_skyline (&buildings);
443 build skyline from a set of line segments.
445 Segments can be articulated from left to right or right to left.
446 In the case of the latter, they will be stored internally as left to right.
448 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
450 list<Building> buildings;
453 for (vsize i = 0; i < segments.size (); i++)
455 Drul_array<Offset> const &seg = segments[i];
456 Offset left = seg[LEFT];
457 Offset right = seg[RIGHT];
458 if (left[horizon_axis] > right[horizon_axis])
459 std::swap (left, right);
461 Real x1 = left[horizon_axis];
462 Real x2 = right[horizon_axis];
463 Real y1 = left[other_axis (horizon_axis)] * sky;
464 Real y2 = right[other_axis (horizon_axis)] * sky;
467 buildings.push_back (Building (x1, y1, y2, x2));
470 buildings_ = internal_build_skyline (&buildings);
474 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
478 deque<Skyline> partials;
479 for (vsize i = 0; i < skypairs.size (); i++)
480 partials.push_back (Skyline ((skypairs[i])[sky]));
482 while (partials.size () > 1)
484 Skyline one = partials.front ();
485 partials.pop_front ();
486 Skyline two = partials.front ();
487 partials.pop_front ();
490 partials.push_back (one);
493 if (partials.size ())
494 buildings_.swap (partials.front ().buildings_);
499 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
502 if (!b.is_empty (X_AXIS) && !b.is_empty (Y_AXIS))
504 Building front (b, horizon_axis, sky);
505 single_skyline (front, &buildings_);
511 Skyline::merge (Skyline const &other)
513 assert (sky_ == other.sky_);
515 if (other.is_empty ())
520 buildings_ = other.buildings_;
524 list<Building> other_bld (other.buildings_);
525 list<Building> my_bld;
526 my_bld.splice (my_bld.begin (), buildings_);
527 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
532 Skyline::insert (Box const &b, Axis a)
534 list<Building> other_bld;
535 list<Building> my_bld;
537 if (isnan (b[other_axis (a)][LEFT])
538 || isnan (b[other_axis (a)][RIGHT]))
540 programming_error ("insane box for skyline");
544 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
545 if (b.is_empty (X_AXIS) || b.is_empty (Y_AXIS))
548 my_bld.splice (my_bld.begin (), buildings_);
549 single_skyline (Building (b, a, sky_), &other_bld);
550 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
555 Skyline::raise (Real r)
557 list<Building>::iterator end = buildings_.end ();
558 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
559 i->y_intercept_ += sky_ * r;
563 Skyline::shift (Real s)
565 list<Building>::iterator end = buildings_.end ();
566 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
570 i->y_intercept_ -= s * i->slope_;
575 Skyline::distance (Skyline const &other, Real horizon_padding) const
578 return internal_distance (other, horizon_padding, &dummy);
582 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
585 internal_distance (other, horizon_padding, &touch);
590 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
592 if (horizon_padding == 0.0)
593 return internal_distance (other, touch_point);
595 // Note that it is not necessary to build a padded version of other,
596 // because the same effect can be achieved just by doubling horizon_padding.
597 Skyline padded_this = padded (horizon_padding);
598 return padded_this.internal_distance (other, touch_point);
602 Skyline::padded (Real horizon_padding) const
604 if (horizon_padding < 0.0)
605 warning ("Cannot have negative horizon padding. Junking.");
607 if (horizon_padding <= 0.0)
610 list<Building> pad_buildings;
611 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
613 if (i->start_ > -infinity_f)
615 Real height = i->height (i->start_);
616 if (height > -infinity_f)
618 // Add the sloped building that pads the left side of the current building.
619 Real start = i->start_ - 2 * horizon_padding;
620 Real end = i->start_ - horizon_padding;
621 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
623 // Add the flat building that pads the left side of the current building.
624 start = i->start_ - horizon_padding;
626 pad_buildings.push_back (Building (start, height, height, end));
630 if (i->end_ < infinity_f)
632 Real height = i->height (i->end_);
633 if (height > -infinity_f)
635 // Add the flat building that pads the right side of the current building.
636 Real start = i->end_;
637 Real end = start + horizon_padding;
638 pad_buildings.push_back (Building (start, height, height, end));
640 // Add the sloped building that pads the right side of the current building.
642 end += horizon_padding;
643 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
648 // The buildings may be overlapping, so resolve that.
649 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
651 // Merge the padding with the original, to make a new skyline.
652 Skyline padded (sky_);
653 list<Building> my_buildings = buildings_;
654 padded.buildings_.clear ();
655 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
662 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
664 assert (sky_ == -other.sky_);
666 list<Building>::const_iterator i = buildings_.begin ();
667 list<Building>::const_iterator j = other.buildings_.begin ();
669 Real dist = -infinity_f;
670 Real start = -infinity_f;
671 Real touch = -infinity_f;
672 while (i != buildings_.end () && j != other.buildings_.end ())
674 Real end = std::min (i->end_, j->end_);
675 Real start_dist = i->height (start) + j->height (start);
676 Real end_dist = i->height (end) + j->height (end);
677 dist = std::max (dist, std::max (start_dist, end_dist));
679 if (end_dist == dist)
681 else if (start_dist == dist)
684 if (i->end_ <= j->end_)
691 *touch_point = touch;
696 Skyline::height (Real airplane) const
698 assert (!isinf (airplane));
700 list<Building>::const_iterator i;
701 for (i = buildings_.begin (); i != buildings_.end (); i++)
703 if (i->end_ >= airplane)
704 return sky_ * i->height (airplane);
712 Skyline::max_height () const
714 Real ret = -infinity_f;
716 list<Building>::const_iterator i;
717 for (i = buildings_.begin (); i != buildings_.end (); ++i)
719 ret = std::max (ret, i->height (i->start_));
720 ret = std::max (ret, i->height (i->end_));
727 Skyline::direction () const
733 Skyline::left () const
735 for (list<Building>::const_iterator i (buildings_.begin ());
736 i != buildings_.end (); i++)
737 if (i->y_intercept_ > -infinity_f)
744 Skyline::right () const
746 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
747 i != buildings_.rend (); ++i)
748 if (i->y_intercept_ > -infinity_f)
755 Skyline::max_height_position () const
758 s.set_minimum_height (0);
759 return touching_point (s);
763 Skyline::set_minimum_height (Real h)
766 s.buildings_.front ().y_intercept_ = h * sky_;
771 Skyline::to_points (Axis horizon_axis) const
775 Real start = -infinity_f;
776 for (list<Building>::const_iterator i (buildings_.begin ());
777 i != buildings_.end (); i++)
779 out.push_back (Offset (start, sky_ * i->height (start)));
780 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
784 if (horizon_axis == Y_AXIS)
785 for (vsize i = 0; i < out.size (); i++)
786 out[i] = out[i].swapped ();
792 Skyline::is_empty () const
794 if (!buildings_.size ())
796 Building b = buildings_.front ();
797 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
804 empty_skyline (&buildings_);
807 /****************************************************************/
809 const char Skyline::type_p_name_[] = "ly:skyline?";
811 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
813 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
815 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
817 Real horizon_padding = 0;
818 if (!SCM_UNBNDP (horizon_padding_scm))
820 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
821 horizon_padding = scm_to_double (horizon_padding_scm);
824 Skyline *skyline = unsmob<Skyline> (skyline_scm);
825 Skyline *other_skyline = unsmob<Skyline> (other_skyline_scm);
826 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
829 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
831 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
833 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
835 Real horizon_padding = 0;
836 if (!SCM_UNBNDP (horizon_padding_scm))
838 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
839 horizon_padding = scm_to_double (horizon_padding_scm);
842 Skyline *skyline = unsmob<Skyline> (skyline_scm);
843 Skyline *other_skyline = unsmob<Skyline> (other_skyline_scm);
844 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
847 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
849 Skyline::get_max_height (SCM skyline_scm)
851 return scm_from_double (unsmob<Skyline> (skyline_scm)->max_height ());
854 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
856 Skyline::get_max_height_position (SCM skyline_scm)
858 return scm_from_double (unsmob<Skyline> (skyline_scm)->max_height_position ());
861 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
863 Skyline::get_height (SCM skyline_scm, SCM x_scm)
865 Real x = robust_scm2double (x_scm, 0.0);
866 return scm_from_double (unsmob<Skyline> (skyline_scm)->height (x));
869 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
871 "Return whether @var{sky} is empty.")
873 Skyline *s = unsmob<Skyline> (sky);
874 LY_ASSERT_SMOB (Skyline, sky, 1);
875 return scm_from_bool (s->is_empty ());