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);
238 // Since a skyline should always be normalized, we can
239 // assume that there are never two adjacent empty buildings.
240 // That is, if center is empty then left and right are not.
241 list<Building>::iterator left = buildings_.begin ();
242 list<Building>::iterator center = buildings_.begin ();
243 list<Building>::iterator right;
245 for (right = buildings_.begin (); right != buildings_.end (); right++)
247 if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
249 printf ("We are here");
250 exit (17); // not-reached
251 Real p1 = left->height (left->end_);
252 Real p2 = right->height (right->start_);
253 *center = Building (center->start_, p1, p2, center->end_);
262 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
263 list<Building> *const result) const
265 if (s1->empty () && s2->empty ())
266 programming_error ("tried to merge two empty skylines");
268 bool s2_pending = false;
269 Real last_end = -infinity_f;
270 while (!s1->empty () && !s2->empty ())
272 Building b = s1->front ();
273 Building c = s2->front ();
274 if (s2_pending // the head of s2 starts before the head of s1
275 || c.height (last_end) > b.height (last_end))
280 // Now b, the head of s1, starts earlier or lies above s2.
282 // Optimization: if the other skyline is empty at this point,
283 // we can avoid testing some intersections. Just grab as many
284 // buildings from s1 as we can, and shove them onto the output.
285 if (c.y_intercept_ == -infinity_f
288 list<Building>::iterator i = s1->begin ();
290 while (i != s1->end () && i->end_ <= c.end_)
293 s1->front ().start_ = last_end;
294 result->splice (result->end (), *s1, s1->begin (), i);
295 last_end = result->back ().end_;
298 // first_intersection() removes buildings from s2 if b hides them
299 Real x = first_intersection (b, s2, last_end);
305 b.leading_part (x); // chops b, leaving the part up to x
308 result->push_back (b);
313 if (b.end_ <= c.end_)
315 // Any trailing part of b is concealed by c
317 // Consider c next for inclusion in the merged skyline
322 // the trailing part of c is fully exposed, goes directly to merged
326 result->push_back (c);
329 if (last_end < infinity_f)
333 s1->front ().start_ = last_end;
334 result->splice (result->end (), *s1);
336 else if (!s2->empty ())
338 s2->front ().start_ = last_end;
339 result->splice (result->end (), *s2);
345 empty_skyline (list<Building> *const ret)
347 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
351 Given Building 'b', build a skyline containing only that building.
354 single_skyline (Building b, list<Building> *const ret)
356 assert (b.end_ >= b.start_);
358 if (b.start_ != -infinity_f)
359 ret->push_back (Building (-infinity_f, -infinity_f,
360 -infinity_f, b.start_));
362 if (b.end_ != infinity_f)
363 ret->push_back (Building (b.end_, -infinity_f,
364 -infinity_f, infinity_f));
367 /* remove a non-overlapping set of boxes from BOXES and build a skyline
369 static list<Building>
370 non_overlapping_skyline (list<Building> *const buildings)
372 list<Building> result;
373 Real last_end = -infinity_f;
374 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
375 list<Building>::iterator i = buildings->begin ();
376 while (i != buildings->end ())
379 Real y1 = i->height (i->start_);
381 Real y2 = i->height (i->end_);
383 // Drop buildings that will obviously have no effect.
384 if (last_building.height (x1) >= y1
385 && last_building.end_ >= x2
386 && last_building.height (x2) >= y2)
388 list<Building>::iterator j = i++;
389 buildings->erase (j);
399 // Insert empty Buildings into any gaps. (TODO: is this needed? -KOH)
401 result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
403 result.push_back (*i);
407 list<Building>::iterator j = i++;
408 buildings->erase (j);
411 if (last_end < infinity_f)
412 result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
416 class LessThanBuilding
419 bool operator () (Building const &b1, Building const &b2)
421 return b1.start_ < b2.start_
422 || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
427 BUILDINGS is a list of buildings, but they could be overlapping
428 and in any order. The returned list of buildings is ordered and non-overlapping.
431 Skyline::internal_build_skyline (list<Building> *buildings) const
433 vsize size = buildings->size ();
437 list<Building> result;
438 empty_skyline (&result);
443 list<Building> result;
444 single_skyline (buildings->front (), &result);
448 deque<list<Building> > partials;
449 buildings->sort (LessThanBuilding ());
450 while (!buildings->empty ())
451 partials.push_back (non_overlapping_skyline (buildings));
453 /* we'd like to say while (partials->size () > 1) but that's O (n).
454 Instead, we exit in the middle of the loop */
455 while (!partials.empty ())
457 list<Building> merged;
458 list<Building> one = partials.front ();
459 partials.pop_front ();
460 if (partials.empty ())
463 list<Building> two = partials.front ();
464 partials.pop_front ();
465 internal_merge_skyline (&one, &two, &merged);
466 partials.push_back (merged);
469 return list<Building> ();
475 empty_skyline (&buildings_);
478 Skyline::Skyline (Direction sky)
481 empty_skyline (&buildings_);
485 Build skyline from a set of boxes.
487 Boxes should be non-empty on both axes. Otherwise, they will be ignored
489 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
491 list<Building> buildings;
494 for (vsize i = 0; i < boxes.size (); i++)
495 if (!boxes[i].is_empty (X_AXIS)
496 && !boxes[i].is_empty (Y_AXIS))
497 buildings.push_front (Building (boxes[i], horizon_axis, sky));
499 buildings_ = internal_build_skyline (&buildings);
504 build skyline from a set of line segments.
506 Segments can be articulated from left to right or right to left.
507 In the case of the latter, they will be stored internally as left to right.
509 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
511 list<Building> buildings;
514 for (vsize i = 0; i < segments.size (); i++)
516 Drul_array<Offset> const &seg = segments[i];
517 Offset left = seg[LEFT];
518 Offset right = seg[RIGHT];
519 if (left[horizon_axis] > right[horizon_axis])
522 Real x1 = left[horizon_axis];
523 Real x2 = right[horizon_axis];
524 Real y1 = left[other_axis (horizon_axis)] * sky;
525 Real y2 = right[other_axis (horizon_axis)] * sky;
528 buildings.push_back (Building (x1, y1, y2, x2));
531 buildings_ = internal_build_skyline (&buildings);
535 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
539 deque<Skyline> partials;
540 for (vsize i = 0; i < skypairs.size (); i++)
541 partials.push_back (Skyline ((skypairs[i])[sky]));
543 while (partials.size () > 1)
545 Skyline one = partials.front ();
546 partials.pop_front ();
547 Skyline two = partials.front ();
548 partials.pop_front ();
551 partials.push_back (one);
554 if (partials.size ())
555 buildings_.swap (partials.front ().buildings_);
560 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
563 if (!b.is_empty (X_AXIS) && !b.is_empty (Y_AXIS))
565 Building front (b, horizon_axis, sky);
566 single_skyline (front, &buildings_);
572 Skyline::merge (Skyline const &other)
574 assert (sky_ == other.sky_);
576 if (other.is_empty ())
581 buildings_ = other.buildings_;
585 list<Building> other_bld (other.buildings_);
586 list<Building> my_bld;
587 my_bld.splice (my_bld.begin (), buildings_);
588 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
593 Skyline::insert (Box const &b, Axis a)
595 list<Building> other_bld;
596 list<Building> my_bld;
598 if (isnan (b[other_axis (a)][LEFT])
599 || isnan (b[other_axis (a)][RIGHT]))
601 programming_error ("insane box for skyline");
605 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
606 if (b.is_empty (X_AXIS) || b.is_empty (Y_AXIS))
609 my_bld.splice (my_bld.begin (), buildings_);
610 single_skyline (Building (b, a, sky_), &other_bld);
611 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
616 Skyline::raise (Real r)
618 list<Building>::iterator end = buildings_.end ();
619 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
620 i->y_intercept_ += sky_ * r;
624 Skyline::shift (Real s)
626 list<Building>::iterator end = buildings_.end ();
627 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
631 i->y_intercept_ -= s * i->slope_;
636 Skyline::distance (Skyline const &other, Real horizon_padding) const
639 return internal_distance (other, horizon_padding, &dummy);
643 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
646 internal_distance (other, horizon_padding, &touch);
651 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
653 if (horizon_padding == 0.0)
654 return internal_distance (other, touch_point);
656 // Note that it is not necessary to build a padded version of other,
657 // because the same effect can be achieved just by doubling horizon_padding.
658 Skyline padded_this = padded (horizon_padding);
659 return padded_this.internal_distance (other, touch_point);
663 Skyline::padded (Real horizon_padding) const
665 if (horizon_padding < 0.0)
666 warning ("Cannot have negative horizon padding. Junking.");
668 if (horizon_padding <= 0.0)
671 list<Building> pad_buildings;
672 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
674 if (i->start_ > -infinity_f)
676 Real height = i->height (i->start_);
677 if (height > -infinity_f)
679 // Add the sloped building that pads the left side of the current building.
680 Real start = i->start_ - 2 * horizon_padding;
681 Real end = i->start_ - horizon_padding;
682 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
684 // Add the flat building that pads the left side of the current building.
685 start = i->start_ - horizon_padding;
687 pad_buildings.push_back (Building (start, height, height, end));
691 if (i->end_ < infinity_f)
693 Real height = i->height (i->end_);
694 if (height > -infinity_f)
696 // Add the flat building that pads the right side of the current building.
697 Real start = i->end_;
698 Real end = start + horizon_padding;
699 pad_buildings.push_back (Building (start, height, height, end));
701 // Add the sloped building that pads the right side of the current building.
703 end += horizon_padding;
704 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
709 // The buildings may be overlapping, so resolve that.
710 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
712 // Merge the padding with the original, to make a new skyline.
713 Skyline padded (sky_);
714 list<Building> my_buildings = buildings_;
715 padded.buildings_.clear ();
716 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
723 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
725 assert (sky_ == -other.sky_);
727 list<Building>::const_iterator i = buildings_.begin ();
728 list<Building>::const_iterator j = other.buildings_.begin ();
730 Real dist = -infinity_f;
731 Real start = -infinity_f;
732 Real touch = -infinity_f;
733 while (i != buildings_.end () && j != other.buildings_.end ())
735 Real end = min (i->end_, j->end_);
736 Real start_dist = i->height (start) + j->height (start);
737 Real end_dist = i->height (end) + j->height (end);
738 dist = max (dist, max (start_dist, end_dist));
740 if (end_dist == dist)
742 else if (start_dist == dist)
745 if (i->end_ <= j->end_)
752 *touch_point = touch;
757 Skyline::height (Real airplane) const
759 assert (!isinf (airplane));
761 list<Building>::const_iterator i;
762 for (i = buildings_.begin (); i != buildings_.end (); i++)
764 if (i->end_ >= airplane)
765 return sky_ * i->height (airplane);
773 Skyline::max_height () const
775 Real ret = -infinity_f;
777 list<Building>::const_iterator i;
778 for (i = buildings_.begin (); i != buildings_.end (); ++i)
780 ret = max (ret, i->height (i->start_));
781 ret = max (ret, i->height (i->end_));
788 Skyline::direction () const
794 Skyline::left () const
796 for (list<Building>::const_iterator i (buildings_.begin ());
797 i != buildings_.end (); i++)
798 if (i->y_intercept_ > -infinity_f)
805 Skyline::right () const
807 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
808 i != buildings_.rend (); ++i)
809 if (i->y_intercept_ > -infinity_f)
816 Skyline::max_height_position () const
819 s.set_minimum_height (0);
820 return touching_point (s);
824 Skyline::set_minimum_height (Real h)
827 s.buildings_.front ().y_intercept_ = h * sky_;
832 Skyline::to_points (Axis horizon_axis) const
836 Real start = -infinity_f;
837 for (list<Building>::const_iterator i (buildings_.begin ());
838 i != buildings_.end (); i++)
840 out.push_back (Offset (start, sky_ * i->height (start)));
841 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
845 if (horizon_axis == Y_AXIS)
846 for (vsize i = 0; i < out.size (); i++)
847 out[i] = out[i].swapped ();
853 Skyline::is_empty () const
855 if (!buildings_.size ())
857 Building b = buildings_.front ();
858 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
865 empty_skyline (&buildings_);
868 /****************************************************************/
870 const char Skyline::type_p_name_[] = "ly:skyline?";
872 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
874 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
876 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
878 Real horizon_padding = 0;
879 if (horizon_padding_scm != SCM_UNDEFINED)
881 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
882 horizon_padding = scm_to_double (horizon_padding_scm);
885 Skyline *skyline = Skyline::unsmob (skyline_scm);
886 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
887 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
890 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
892 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
894 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
896 Real horizon_padding = 0;
897 if (horizon_padding_scm != SCM_UNDEFINED)
899 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
900 horizon_padding = scm_to_double (horizon_padding_scm);
903 Skyline *skyline = Skyline::unsmob (skyline_scm);
904 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
905 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
908 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
910 Skyline::get_max_height (SCM skyline_scm)
912 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
915 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
917 Skyline::get_max_height_position (SCM skyline_scm)
919 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
922 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
924 Skyline::get_height (SCM skyline_scm, SCM x_scm)
926 Real x = robust_scm2double (x_scm, 0.0);
927 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
930 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
932 "Return whether @var{sky} is empty.")
934 Skyline *s = Skyline::unsmob (sky);
935 LY_ASSERT_SMOB (Skyline, sky, 1);
936 return scm_from_bool (s->is_empty ());