2 This file is part of LilyPond, the GNU music typesetter.
4 Copyright (C) 2006--2014 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"
25 /* A skyline is a sequence of non-overlapping buildings: something like
35 Each building has a starting position, and ending position, a starting
36 height and an ending height.
38 The following invariants are observed:
39 - the start of the first building is at -infinity
40 - the end of the last building is at infinity
41 - if a building has infinite length (ie. the first and last buildings),
42 then its starting height and ending height are equal
43 - the end of one building is the same as the beginning of the next
46 We also allow skylines to point down (the structure is exactly the same,
47 but we think of the part above the line as being filled with mass and the
48 part below as being empty). ::distance finds the minimum distance between
49 an UP skyline and a DOWN skyline.
51 Note that we store DOWN skylines upside-down. That is, in order to compare
52 a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
53 This means that the merging routine doesn't need to be aware of direction,
54 but the distance routine does.
56 From 2007 through 2012, buildings of width less than EPS were discarded,
57 citing numerical accuracy concerns. We remember that floating point
58 comparisons of nearly-equal values can be affected by rounding error.
59 Also, some target machines use the x87 floating point unit, which provides
60 extended precision for intermediate results held in registers. On this type
61 of hardware comparisons such as
62 double c = 1.0/3.0; boolean compare = (c == 1.0/3.0)
63 could go either way because the 1.0/3.0 is allowed to be kept
64 higher precision than the variable 'c'.
65 Alert to these considerations, we now accept buildings of zero-width.
69 print_buildings (list<Building> const &b)
71 for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
76 Skyline::print () const
78 print_buildings (buildings_);
82 Skyline::print_points () const
84 vector<Offset> ps (to_points (X_AXIS));
86 for (vsize i = 0; i < ps.size (); i++)
87 printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
88 (i % 2) == 1 ? "\n" : " ");
91 Building::Building (Real start, Real start_height, Real end_height, Real end)
93 if (isinf (start) || isinf (end))
94 assert (start_height == end_height);
98 precompute (start, start_height, end_height, end);
101 Building::Building (Box const &b, Axis horizon_axis, Direction sky)
103 Real start = b[horizon_axis][LEFT];
104 Real end = b[horizon_axis][RIGHT];
105 Real height = sky * b[other_axis (horizon_axis)][sky];
109 precompute (start, height, height, end);
113 Building::precompute (Real start, Real start_height, Real end_height, Real end)
115 slope_ = 0.0; /* if they were both infinite, we would get nan, not 0, from the prev line */
116 if (start_height != end_height)
117 slope_ = (end_height - start_height) / (end - start);
119 assert (!isinf (slope_) && !isnan (slope_));
123 assert (start_height == end_height);
124 y_intercept_ = start_height;
126 else if (fabs (slope_) > 1e6)
127 // too steep to be stored in slope-intercept form, given round-off error
130 y_intercept_ = max (start_height, end_height);
133 y_intercept_ = start_height - slope_ * start;
137 Building::height (Real x) const
139 return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
143 Building::print () const
145 printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
149 Building::intersection_x (Building const &other) const
151 Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
152 return isnan (ret) ? -infinity_f : ret;
156 Building::leading_part (Real chop)
158 assert (chop <= end_);
162 // Returns a shift s such that (x + s, y) intersects the roof of
163 // this building. If no such shift exists, returns infinity_f.
165 Building::shift_to_intersect (Real x, Real y) const
167 // Solve for s: y = (x + s)*m + b
168 Real ret = (y - y_intercept_ - slope_ * x) / slope_;
170 if (ret >= start_ && ret <= end_ && !isinf (ret))
176 first_intersection (Building const &b, list<Building> *s, Real start_x)
177 /* Return the first x (start_x <= x <= b.end)
178 * where skyline s is clear of Building b.
179 * Removes buildings from s that are fully concealed by b. */
181 for ( ; !s->empty (); s->pop_front ())
183 Building c = s->front ();
185 start_x = max (start_x, c.start_);
186 if (start_x >= b.end_)
189 // height and intersection_x involve multiplication and
190 // division. Avoid that, if we can.
191 if (c.y_intercept_ != -infinity_f)
193 if (c.height (start_x) > b.height (start_x))
196 if (c.height (b.end_) > b.height (b.end_))
198 Real i = b.intersection_x (c);
199 if (i <= b.end_ && i <= c.end_)
200 return max (start_x, i);
203 if (c.end_ >= b.end_)
204 return b.end_; // leave this c on the list s
209 // Remove redundant empty buildings from the skyline.
210 // If there are two adjacent empty buildings, they can be
213 Skyline::normalize ()
215 bool last_empty = false;
216 list<Building>::iterator i;
218 for (i = buildings_.begin (); i != buildings_.end (); i++)
220 if (last_empty && i->y_intercept_ == -infinity_f)
222 list<Building>::iterator last = i;
224 last->end_ = i->end_;
225 buildings_.erase (i);
228 last_empty = (i->y_intercept_ == -infinity_f);
231 assert (buildings_.front ().start_ == -infinity_f);
232 assert (buildings_.back ().end_ == infinity_f);
236 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
237 list<Building> *const result) const
239 if (s1->empty () && s2->empty ())
240 programming_error ("tried to merge two empty skylines");
242 bool s2_pending = false;
243 Real last_end = -infinity_f;
244 while (!s1->empty () && !s2->empty ())
246 Building b = s1->front ();
247 Building c = s2->front ();
248 if (s2_pending // the head of s2 starts before the head of s1
249 || c.height (last_end) > b.height (last_end))
254 // Now b, the head of s1, starts earlier or lies above s2.
256 // Optimization: if the other skyline is empty at this point,
257 // we can avoid testing some intersections. Just grab as many
258 // buildings from s1 as we can, and shove them onto the output.
259 if (c.y_intercept_ == -infinity_f
262 list<Building>::iterator i = s1->begin ();
264 while (i != s1->end () && i->end_ <= c.end_)
267 s1->front ().start_ = last_end;
268 result->splice (result->end (), *s1, s1->begin (), i);
269 last_end = result->back ().end_;
272 // first_intersection() removes buildings from s2 if b hides them
273 Real x = first_intersection (b, s2, last_end);
279 b.leading_part (x); // chops b, leaving the part up to x
282 result->push_back (b);
287 if (b.end_ <= c.end_)
289 // Any trailing part of b is concealed by c
291 // Consider c next for inclusion in the merged skyline
296 // the trailing part of c is fully exposed, goes directly to merged
300 result->push_back (c);
303 if (last_end < infinity_f)
307 s1->front ().start_ = last_end;
308 result->splice (result->end (), *s1);
310 else if (!s2->empty ())
312 s2->front ().start_ = last_end;
313 result->splice (result->end (), *s2);
319 empty_skyline (list<Building> *const ret)
321 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
325 Given Building 'b', build a skyline containing only that building.
328 single_skyline (Building b, list<Building> *const ret)
330 assert (b.end_ >= b.start_);
332 if (b.start_ != -infinity_f)
333 ret->push_back (Building (-infinity_f, -infinity_f,
334 -infinity_f, b.start_));
336 if (b.end_ != infinity_f)
337 ret->push_back (Building (b.end_, -infinity_f,
338 -infinity_f, infinity_f));
341 /* remove a non-overlapping set of boxes from BOXES and build a skyline
343 static list<Building>
344 non_overlapping_skyline (list<Building> *const buildings)
346 list<Building> result;
347 Real last_end = -infinity_f;
348 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
349 list<Building>::iterator i = buildings->begin ();
350 while (i != buildings->end ())
353 Real y1 = i->height (i->start_);
355 Real y2 = i->height (i->end_);
357 // Drop buildings that will obviously have no effect.
358 if (last_building.height (x1) >= y1
359 && last_building.end_ >= x2
360 && last_building.height (x2) >= y2)
362 list<Building>::iterator j = i++;
363 buildings->erase (j);
373 // Insert empty Buildings into any gaps. (TODO: is this needed? -KOH)
375 result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
377 result.push_back (*i);
381 list<Building>::iterator j = i++;
382 buildings->erase (j);
385 if (last_end < infinity_f)
386 result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
390 class LessThanBuilding
393 bool operator () (Building const &b1, Building const &b2)
395 return b1.start_ < b2.start_
396 || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
401 BUILDINGS is a list of buildings, but they could be overlapping
402 and in any order. The returned list of buildings is ordered and non-overlapping.
405 Skyline::internal_build_skyline (list<Building> *buildings) const
407 vsize size = buildings->size ();
411 list<Building> result;
412 empty_skyline (&result);
417 list<Building> result;
418 single_skyline (buildings->front (), &result);
422 deque<list<Building> > partials;
423 buildings->sort (LessThanBuilding ());
424 while (!buildings->empty ())
425 partials.push_back (non_overlapping_skyline (buildings));
427 /* we'd like to say while (partials->size () > 1) but that's O (n).
428 Instead, we exit in the middle of the loop */
429 while (!partials.empty ())
431 list<Building> merged;
432 list<Building> one = partials.front ();
433 partials.pop_front ();
434 if (partials.empty ())
437 list<Building> two = partials.front ();
438 partials.pop_front ();
439 internal_merge_skyline (&one, &two, &merged);
440 partials.push_back (merged);
443 return list<Building> ();
449 empty_skyline (&buildings_);
452 Skyline::Skyline (Direction sky)
455 empty_skyline (&buildings_);
459 Build skyline from a set of boxes.
461 Boxes should be non-empty on both axes. Otherwise, they will be ignored
463 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
465 list<Building> buildings;
468 for (vsize i = 0; i < boxes.size (); i++)
469 if (!boxes[i].is_empty (X_AXIS)
470 && !boxes[i].is_empty (Y_AXIS))
471 buildings.push_front (Building (boxes[i], horizon_axis, sky));
473 buildings_ = internal_build_skyline (&buildings);
478 build skyline from a set of line segments.
480 Segments can be articulated from left to right or right to left.
481 In the case of the latter, they will be stored internally as left to right.
483 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
485 list<Building> buildings;
488 for (vsize i = 0; i < segments.size (); i++)
490 Drul_array<Offset> const &seg = segments[i];
491 Offset left = seg[LEFT];
492 Offset right = seg[RIGHT];
493 if (left[horizon_axis] > right[horizon_axis])
496 Real x1 = left[horizon_axis];
497 Real x2 = right[horizon_axis];
498 Real y1 = left[other_axis (horizon_axis)] * sky;
499 Real y2 = right[other_axis (horizon_axis)] * sky;
502 buildings.push_back (Building (x1, y1, y2, x2));
505 buildings_ = internal_build_skyline (&buildings);
509 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
513 deque<Skyline> partials;
514 for (vsize i = 0; i < skypairs.size (); i++)
515 partials.push_back (Skyline ((skypairs[i])[sky]));
517 while (partials.size () > 1)
519 Skyline one = partials.front ();
520 partials.pop_front ();
521 Skyline two = partials.front ();
522 partials.pop_front ();
525 partials.push_back (one);
528 if (partials.size ())
529 buildings_.swap (partials.front ().buildings_);
534 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
537 if (!b.is_empty (X_AXIS) && !b.is_empty (Y_AXIS))
539 Building front (b, horizon_axis, sky);
540 single_skyline (front, &buildings_);
546 Skyline::merge (Skyline const &other)
548 assert (sky_ == other.sky_);
550 if (other.is_empty ())
555 buildings_ = other.buildings_;
559 list<Building> other_bld (other.buildings_);
560 list<Building> my_bld;
561 my_bld.splice (my_bld.begin (), buildings_);
562 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
567 Skyline::insert (Box const &b, Axis a)
569 list<Building> other_bld;
570 list<Building> my_bld;
572 if (isnan (b[other_axis (a)][LEFT])
573 || isnan (b[other_axis (a)][RIGHT]))
575 programming_error ("insane box for skyline");
579 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
580 if (b.is_empty (X_AXIS) || b.is_empty (Y_AXIS))
583 my_bld.splice (my_bld.begin (), buildings_);
584 single_skyline (Building (b, a, sky_), &other_bld);
585 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
590 Skyline::raise (Real r)
592 list<Building>::iterator end = buildings_.end ();
593 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
594 i->y_intercept_ += sky_ * r;
598 Skyline::shift (Real s)
600 list<Building>::iterator end = buildings_.end ();
601 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
605 i->y_intercept_ -= s * i->slope_;
610 Skyline::distance (Skyline const &other, Real horizon_padding) const
613 return internal_distance (other, horizon_padding, &dummy);
617 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
620 internal_distance (other, horizon_padding, &touch);
625 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
627 if (horizon_padding == 0.0)
628 return internal_distance (other, touch_point);
630 // Note that it is not necessary to build a padded version of other,
631 // because the same effect can be achieved just by doubling horizon_padding.
632 Skyline padded_this = padded (horizon_padding);
633 return padded_this.internal_distance (other, touch_point);
637 Skyline::padded (Real horizon_padding) const
639 if (horizon_padding < 0.0)
640 warning ("Cannot have negative horizon padding. Junking.");
642 if (horizon_padding <= 0.0)
645 list<Building> pad_buildings;
646 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
648 if (i->start_ > -infinity_f)
650 Real height = i->height (i->start_);
651 if (height > -infinity_f)
653 // Add the sloped building that pads the left side of the current building.
654 Real start = i->start_ - 2 * horizon_padding;
655 Real end = i->start_ - horizon_padding;
656 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
658 // Add the flat building that pads the left side of the current building.
659 start = i->start_ - horizon_padding;
661 pad_buildings.push_back (Building (start, height, height, end));
665 if (i->end_ < infinity_f)
667 Real height = i->height (i->end_);
668 if (height > -infinity_f)
670 // Add the flat building that pads the right side of the current building.
671 Real start = i->end_;
672 Real end = start + horizon_padding;
673 pad_buildings.push_back (Building (start, height, height, end));
675 // Add the sloped building that pads the right side of the current building.
677 end += horizon_padding;
678 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
683 // The buildings may be overlapping, so resolve that.
684 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
686 // Merge the padding with the original, to make a new skyline.
687 Skyline padded (sky_);
688 list<Building> my_buildings = buildings_;
689 padded.buildings_.clear ();
690 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
697 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
699 assert (sky_ == -other.sky_);
701 list<Building>::const_iterator i = buildings_.begin ();
702 list<Building>::const_iterator j = other.buildings_.begin ();
704 Real dist = -infinity_f;
705 Real start = -infinity_f;
706 Real touch = -infinity_f;
707 while (i != buildings_.end () && j != other.buildings_.end ())
709 Real end = min (i->end_, j->end_);
710 Real start_dist = i->height (start) + j->height (start);
711 Real end_dist = i->height (end) + j->height (end);
712 dist = max (dist, max (start_dist, end_dist));
714 if (end_dist == dist)
716 else if (start_dist == dist)
719 if (i->end_ <= j->end_)
726 *touch_point = touch;
731 Skyline::height (Real airplane) const
733 assert (!isinf (airplane));
735 list<Building>::const_iterator i;
736 for (i = buildings_.begin (); i != buildings_.end (); i++)
738 if (i->end_ >= airplane)
739 return sky_ * i->height (airplane);
747 Skyline::max_height () const
749 Real ret = -infinity_f;
751 list<Building>::const_iterator i;
752 for (i = buildings_.begin (); i != buildings_.end (); ++i)
754 ret = max (ret, i->height (i->start_));
755 ret = max (ret, i->height (i->end_));
762 Skyline::direction () const
768 Skyline::left () const
770 for (list<Building>::const_iterator i (buildings_.begin ());
771 i != buildings_.end (); i++)
772 if (i->y_intercept_ > -infinity_f)
779 Skyline::right () const
781 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
782 i != buildings_.rend (); ++i)
783 if (i->y_intercept_ > -infinity_f)
790 Skyline::max_height_position () const
793 s.set_minimum_height (0);
794 return touching_point (s);
798 Skyline::set_minimum_height (Real h)
801 s.buildings_.front ().y_intercept_ = h * sky_;
806 Skyline::to_points (Axis horizon_axis) const
810 Real start = -infinity_f;
811 for (list<Building>::const_iterator i (buildings_.begin ());
812 i != buildings_.end (); i++)
814 out.push_back (Offset (start, sky_ * i->height (start)));
815 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
819 if (horizon_axis == Y_AXIS)
820 for (vsize i = 0; i < out.size (); i++)
821 out[i] = out[i].swapped ();
827 Skyline::is_empty () const
829 if (!buildings_.size ())
831 Building b = buildings_.front ();
832 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
839 empty_skyline (&buildings_);
842 /****************************************************************/
844 const char Skyline::type_p_name_[] = "ly:skyline?";
846 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
848 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
850 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
852 Real horizon_padding = 0;
853 if (horizon_padding_scm != SCM_UNDEFINED)
855 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
856 horizon_padding = scm_to_double (horizon_padding_scm);
859 Skyline *skyline = Skyline::unsmob (skyline_scm);
860 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
861 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
864 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
866 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
868 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
870 Real horizon_padding = 0;
871 if (horizon_padding_scm != SCM_UNDEFINED)
873 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
874 horizon_padding = scm_to_double (horizon_padding_scm);
877 Skyline *skyline = Skyline::unsmob (skyline_scm);
878 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
879 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
882 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
884 Skyline::get_max_height (SCM skyline_scm)
886 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
889 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
891 Skyline::get_max_height_position (SCM skyline_scm)
893 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
896 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
898 Skyline::get_height (SCM skyline_scm, SCM x_scm)
900 Real x = robust_scm2double (x_scm, 0.0);
901 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
904 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
906 "Return whether @var{sky} is empty.")
908 Skyline *s = Skyline::unsmob (sky);
909 LY_ASSERT_SMOB (Skyline, sky, 1);
910 return scm_from_bool (s->is_empty ());