2 This file is part of LilyPond, the GNU music typesetter.
4 Copyright (C) 2006--2012 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 #include "ly-smobs.icc"
27 /* A skyline is a sequence of non-overlapping buildings: something like
37 Each building has a starting position, and ending position, a starting
38 height and an ending height.
40 The following invariants are observed:
41 - the start of the first building is at -infinity
42 - the end of the last building is at infinity
43 - if a building has infinite length (ie. the first and last buildings),
44 then its starting height and ending height are equal
45 - the end of one building is the same as the beginning of the next
48 We also allow skylines to point down (the structure is exactly the same,
49 but we think of the part above the line as being filled with mass and the
50 part below as being empty). ::distance finds the minimum distance between
51 an UP skyline and a DOWN skyline.
53 Note that we store DOWN skylines upside-down. That is, in order to compare
54 a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
55 This means that the merging routine doesn't need to be aware of direction,
56 but the distance routine does.
58 From 2007 through 2012, buildings of width less than EPS were discarded,
59 citing numerical accuracy concerns. We remember that floating point
60 comparisons of nearly-equal values can be affected by rounding error.
61 Also, some target machines use the x87 floating point unit, which provides
62 extended precision for intermediate results held in registers. On this type
63 of hardware comparisons such as
64 double c = 1.0/3.0; boolean compare = (c == 1.0/3.0)
65 could go either way because the 1.0/3.0 is allowed to be kept
66 higher precision than the variable 'c'.
67 Alert to these considerations, we now accept buildings of zero-width.
71 print_buildings (list<Building> const &b)
73 for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
78 Skyline::print () const
80 print_buildings (buildings_);
84 Skyline::print_points () const
86 vector<Offset> ps (to_points (X_AXIS));
88 for (vsize i = 0; i < ps.size (); i++)
89 printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
90 (i % 2) == 1 ? "\n" : " ");
93 Building::Building (Real start, Real start_height, Real end_height, Real end)
95 if (isinf (start) || isinf (end))
96 assert (start_height == end_height);
100 precompute (start, start_height, end_height, end);
103 Building::Building (Box const &b, Axis horizon_axis, Direction sky)
105 Real start = b[horizon_axis][LEFT];
106 Real end = b[horizon_axis][RIGHT];
107 Real height = sky * b[other_axis (horizon_axis)][sky];
111 precompute (start, height, height, end);
115 Building::precompute (Real start, Real start_height, Real end_height, Real end)
117 slope_ = 0.0; /* if they were both infinite, we would get nan, not 0, from the prev line */
118 if (start_height != end_height)
119 slope_ = (end_height - start_height) / (end - start);
121 assert (!isinf (slope_) && !isnan (slope_));
125 assert (start_height == end_height);
126 y_intercept_ = start_height;
128 else if (fabs(slope_) > 1e6)
129 // too steep to be stored in slope-intercept form, given round-off error
132 y_intercept_ = max(start_height, end_height);
135 y_intercept_ = start_height - slope_ * start;
139 Building::height (Real x) const
141 return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
145 Building::print () const
147 printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
151 Building::intersection_x (Building const &other) const
153 Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
154 return isnan (ret) ? -infinity_f : ret;
158 Building::leading_part (Real chop)
160 assert (chop <= end_);
164 // Returns a shift s such that (x + s, y) intersects the roof of
165 // this building. If no such shift exists, returns infinity_f.
167 Building::shift_to_intersect (Real x, Real y) const
169 // Solve for s: y = (x + s)*m + b
170 Real ret = (y - y_intercept_ - slope_ * x) / slope_;
172 if (ret >= start_ && ret <= end_ && !isinf (ret))
178 first_intersection (Building const &b, list<Building> *s, Real start_x)
179 /* Return the first x >= start_x where skyline s above Building b.
180 * Removes buildings from s that are concealed by b. */
182 while (!s->empty () && start_x < b.end_)
184 Building c = s->front ();
186 // conceals and intersection_x involve multiplication and
187 // division. Avoid that, if we can.
188 if (c.y_intercept_ == -infinity_f)
197 if (c.conceals (b, start_x))
200 Real i = b.intersection_x (c);
201 if (i > start_x && i <= b.end_ && i <= c.end_)
212 Building::conceals (Building const &other, Real x) const
214 if (slope_ == other.slope_)
215 return y_intercept_ > other.y_intercept_;
217 /* their slopes were not equal, so there is an intersection point */
218 Real i = intersection_x (other);
219 return (i <= x && slope_ > other.slope_)
220 || (i > x && slope_ < other.slope_);
223 // Remove redundant empty buildings from the skyline.
224 // If there are two adjacent empty buildings, they can be
227 Skyline::normalize ()
229 bool last_empty = false;
230 list<Building>::iterator i;
232 for (i = buildings_.begin (); i != buildings_.end (); i++)
234 if (last_empty && i->y_intercept_ == -infinity_f)
236 list<Building>::iterator last = i;
238 last->end_ = i->end_;
239 buildings_.erase (i);
242 last_empty = (i->y_intercept_ == -infinity_f);
245 assert (buildings_.front ().start_ == -infinity_f);
246 assert (buildings_.back ().end_ == infinity_f);
252 // Since a skyline should always be normalized, we can
253 // assume that there are never two adjacent empty buildings.
254 // That is, if center is empty then left and right are not.
255 list<Building>::iterator left = buildings_.begin ();
256 list<Building>::iterator center = buildings_.begin ();
257 list<Building>::iterator right;
259 for (right = buildings_.begin (); right != buildings_.end (); right++)
261 if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
263 Real p1 = left->height (left->end_);
264 Real p2 = right->height (right->start_);
265 *center = Building (center->start_, p1, p2, center->end_);
274 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
275 list<Building> *const result) const
277 if (s1->empty () || s2->empty ())
279 programming_error ("tried to merge an empty skyline");
283 Real x = -infinity_f;
284 Real last_end = -infinity_f;
285 while (!s1->empty ())
287 if (s2->front ().conceals (s1->front (), x))
290 Building b = s1->front ();
291 Building c = s2->front ();
293 // Optimization: if the other skyline is empty at this point,
294 // we can avoid testing some intersections. Just grab as many
295 // buildings from s1 as we can, and shove them onto the output.
296 if (c.y_intercept_ == -infinity_f
299 list<Building>::iterator i = s1->begin ();
301 while (i != s1->end () && i->end_ <= c.end_)
304 s1->front ().start_ = x;
305 result->splice (result->end (), *s1, s1->begin (), i);
306 x = result->back ().end_;
310 // first_intersection() removes buildings from s2 if b hides them
311 Real end = first_intersection (b, s2, x);
315 result->push_back (b);
319 // Should be (end > x), during ver2.19. end == x happens fairly often,
320 // and we do not need to keep vertical segments within a skyline.
323 b.leading_part (end);
326 result->push_back (b);
329 if (end >= s1->front ().end_)
331 // Should add during ver2.19 (to avoid an endless loop
332 // when merging identical skylines with a vertical segment)
333 // if (end >= s2->front().end_) s2->pop_front();
340 empty_skyline (list<Building> *const ret)
342 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
346 Given Building 'b', build a skyline containing only that building.
349 single_skyline (Building b, list<Building> *const ret)
351 assert (b.end_ >= b.start_);
353 if (b.start_ != -infinity_f)
354 ret->push_back (Building (-infinity_f, -infinity_f,
355 -infinity_f, b.start_));
357 if (b.end_ != infinity_f)
358 ret->push_back (Building (b.end_, -infinity_f,
359 -infinity_f, infinity_f));
362 /* remove a non-overlapping set of boxes from BOXES and build a skyline
364 static list<Building>
365 non_overlapping_skyline (list<Building> *const buildings)
367 list<Building> result;
368 Real last_end = -infinity_f;
369 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
370 list<Building>::iterator i = buildings->begin ();
371 while (i != buildings->end ())
374 Real y1 = i->height (i->start_);
376 Real y2 = i->height (i->end_);
378 // Drop buildings that will obviously have no effect.
379 if (last_building.height (x1) >= y1
380 && last_building.end_ >= x2
381 && last_building.height (x2) >= y2)
383 list<Building>::iterator j = i++;
384 buildings->erase (j);
394 // Insert empty Buildings into any gaps. (TODO: is this needed? -KOH)
396 result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
398 result.push_back (*i);
402 list<Building>::iterator j = i++;
403 buildings->erase (j);
406 if (last_end < infinity_f)
407 result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
411 class LessThanBuilding
414 bool operator () (Building const &b1, Building const &b2)
416 return b1.start_ < b2.start_
417 || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
422 BUILDINGS is a list of buildings, but they could be overlapping
423 and in any order. The returned list of buildings is ordered and non-overlapping.
426 Skyline::internal_build_skyline (list<Building> *buildings) const
428 vsize size = buildings->size ();
432 list<Building> result;
433 empty_skyline (&result);
438 list<Building> result;
439 single_skyline (buildings->front (), &result);
443 deque<list<Building> > partials;
444 buildings->sort (LessThanBuilding ());
445 while (!buildings->empty ())
446 partials.push_back (non_overlapping_skyline (buildings));
448 /* we'd like to say while (partials->size () > 1) but that's O (n).
449 Instead, we exit in the middle of the loop */
450 while (!partials.empty ())
452 list<Building> merged;
453 list<Building> one = partials.front ();
454 partials.pop_front ();
455 if (partials.empty ())
458 list<Building> two = partials.front ();
459 partials.pop_front ();
460 internal_merge_skyline (&one, &two, &merged);
461 partials.push_back (merged);
464 return list<Building> ();
470 empty_skyline (&buildings_);
473 Skyline::Skyline (Direction sky)
476 empty_skyline (&buildings_);
480 Build skyline from a set of boxes.
482 Boxes should be non-empty on both axes. Otherwise, they will be ignored
484 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
486 list<Building> buildings;
489 for (vsize i = 0; i < boxes.size (); i++)
490 if (!boxes[i].is_empty (X_AXIS)
491 && !boxes[i].is_empty (Y_AXIS))
492 buildings.push_front (Building (boxes[i], horizon_axis, sky));
494 buildings_ = internal_build_skyline (&buildings);
499 build skyline from a set of line segments.
501 Segments can be articulated from left to right or right to left.
502 In the case of the latter, they will be stored internally as left to right.
504 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
506 list<Building> buildings;
509 for (vsize i = 0; i < segments.size (); i++)
511 Drul_array<Offset> const &seg = segments[i];
512 Offset left = seg[LEFT];
513 Offset right = seg[RIGHT];
514 if (left[horizon_axis] > right[horizon_axis])
517 Real x1 = left[horizon_axis];
518 Real x2 = right[horizon_axis];
519 Real y1 = left[other_axis (horizon_axis)] * sky;
520 Real y2 = right[other_axis (horizon_axis)] * sky;
523 buildings.push_back (Building (x1, y1, y2, x2));
526 buildings_ = internal_build_skyline (&buildings);
530 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
534 deque<Skyline> partials;
535 for (vsize i = 0; i < skypairs.size (); i++)
536 partials.push_back (Skyline ((skypairs[i])[sky]));
538 while (partials.size () > 1)
540 Skyline one = partials.front ();
541 partials.pop_front ();
542 Skyline two = partials.front ();
543 partials.pop_front ();
546 partials.push_back (one);
549 if (partials.size ())
550 buildings_.swap (partials.front ().buildings_);
555 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
558 if (!b.is_empty (X_AXIS) && !b.is_empty (Y_AXIS))
560 Building front (b, horizon_axis, sky);
561 single_skyline (front, &buildings_);
567 Skyline::merge (Skyline const &other)
569 assert (sky_ == other.sky_);
571 if (other.is_empty ())
576 buildings_ = other.buildings_;
580 list<Building> other_bld (other.buildings_);
581 list<Building> my_bld;
582 my_bld.splice (my_bld.begin (), buildings_);
583 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
588 Skyline::insert (Box const &b, Axis a)
590 list<Building> other_bld;
591 list<Building> my_bld;
593 if (isnan (b[other_axis (a)][LEFT])
594 || isnan (b[other_axis (a)][RIGHT]))
596 programming_error ("insane box for skyline");
600 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
601 if (b.is_empty (X_AXIS) || b.is_empty (Y_AXIS))
604 my_bld.splice (my_bld.begin (), buildings_);
605 single_skyline (Building (b, a, sky_), &other_bld);
606 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
611 Skyline::raise (Real r)
613 list<Building>::iterator end = buildings_.end ();
614 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
615 i->y_intercept_ += sky_ * r;
619 Skyline::shift (Real s)
621 list<Building>::iterator end = buildings_.end ();
622 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
626 i->y_intercept_ -= s * i->slope_;
631 Skyline::distance (Skyline const &other, Real horizon_padding) const
634 return internal_distance (other, horizon_padding, &dummy);
638 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
641 internal_distance (other, horizon_padding, &touch);
646 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
648 if (horizon_padding == 0.0)
649 return internal_distance (other, touch_point);
651 // Note that it is not necessary to build a padded version of other,
652 // because the same effect can be achieved just by doubling horizon_padding.
653 Skyline padded_this = padded (horizon_padding);
654 return padded_this.internal_distance (other, touch_point);
658 Skyline::padded (Real horizon_padding) const
660 if (horizon_padding < 0.0)
661 warning ("Cannot have negative horizon padding. Junking.");
663 if (horizon_padding <= 0.0)
666 list<Building> pad_buildings;
667 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
669 if (i->start_ > -infinity_f)
671 Real height = i->height (i->start_);
672 if (height > -infinity_f)
674 // Add the sloped building that pads the left side of the current building.
675 Real start = i->start_ - 2 * horizon_padding;
676 Real end = i->start_ - horizon_padding;
677 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
679 // Add the flat building that pads the left side of the current building.
680 start = i->start_ - horizon_padding;
682 pad_buildings.push_back (Building (start, height, height, end));
686 if (i->end_ < infinity_f)
688 Real height = i->height (i->end_);
689 if (height > -infinity_f)
691 // Add the flat building that pads the right side of the current building.
692 Real start = i->end_;
693 Real end = start + horizon_padding;
694 pad_buildings.push_back (Building (start, height, height, end));
696 // Add the sloped building that pads the right side of the current building.
698 end += horizon_padding;
699 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
704 // The buildings may be overlapping, so resolve that.
705 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
707 // Merge the padding with the original, to make a new skyline.
708 Skyline padded (sky_);
709 list<Building> my_buildings = buildings_;
710 padded.buildings_.clear ();
711 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
718 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
720 assert (sky_ == -other.sky_);
722 list<Building>::const_iterator i = buildings_.begin ();
723 list<Building>::const_iterator j = other.buildings_.begin ();
725 Real dist = -infinity_f;
726 Real start = -infinity_f;
727 Real touch = -infinity_f;
728 while (i != buildings_.end () && j != other.buildings_.end ())
730 Real end = min (i->end_, j->end_);
731 Real start_dist = i->height (start) + j->height (start);
732 Real end_dist = i->height (end) + j->height (end);
733 dist = max (dist, max (start_dist, end_dist));
735 if (end_dist == dist)
737 else if (start_dist == dist)
740 if (i->end_ <= j->end_)
747 *touch_point = touch;
752 Skyline::height (Real airplane) const
754 assert (!isinf (airplane));
756 list<Building>::const_iterator i;
757 for (i = buildings_.begin (); i != buildings_.end (); i++)
759 if (i->end_ >= airplane)
760 return sky_ * i->height (airplane);
768 Skyline::max_height () const
770 Real ret = -infinity_f;
772 list<Building>::const_iterator i;
773 for (i = buildings_.begin (); i != buildings_.end (); ++i)
775 ret = max (ret, i->height (i->start_));
776 ret = max (ret, i->height (i->end_));
783 Skyline::direction () const
789 Skyline::left () const
791 for (list<Building>::const_iterator i (buildings_.begin ());
792 i != buildings_.end (); i++)
793 if (i->y_intercept_ > -infinity_f)
800 Skyline::right () const
802 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
803 i != buildings_.rend (); ++i)
804 if (i->y_intercept_ > -infinity_f)
811 Skyline::max_height_position () const
814 s.set_minimum_height (0);
815 return touching_point (s);
819 Skyline::set_minimum_height (Real h)
822 s.buildings_.front ().y_intercept_ = h * sky_;
827 Skyline::to_points (Axis horizon_axis) const
831 Real start = -infinity_f;
832 for (list<Building>::const_iterator i (buildings_.begin ());
833 i != buildings_.end (); i++)
835 out.push_back (Offset (start, sky_ * i->height (start)));
836 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
840 if (horizon_axis == Y_AXIS)
841 for (vsize i = 0; i < out.size (); i++)
842 out[i] = out[i].swapped ();
848 Skyline::is_empty () const
850 if (!buildings_.size ())
852 Building b = buildings_.front ();
853 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
860 empty_skyline (&buildings_);
863 /****************************************************************/
865 IMPLEMENT_SIMPLE_SMOBS (Skyline);
866 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
867 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
870 Skyline::mark_smob (SCM s)
872 ASSERT_LIVE_IS_ALLOWED (s);
877 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
879 Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
882 scm_puts ("#<Skyline>", port);
887 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
889 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
891 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
893 Real horizon_padding = 0;
894 if (horizon_padding_scm != SCM_UNDEFINED)
896 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
897 horizon_padding = scm_to_double (horizon_padding_scm);
900 Skyline *skyline = Skyline::unsmob (skyline_scm);
901 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
902 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
905 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
907 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
909 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
911 Real horizon_padding = 0;
912 if (horizon_padding_scm != SCM_UNDEFINED)
914 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
915 horizon_padding = scm_to_double (horizon_padding_scm);
918 Skyline *skyline = Skyline::unsmob (skyline_scm);
919 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
920 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
923 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
925 Skyline::get_max_height (SCM skyline_scm)
927 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
930 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
932 Skyline::get_max_height_position (SCM skyline_scm)
934 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
937 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
939 Skyline::get_height (SCM skyline_scm, SCM x_scm)
941 Real x = robust_scm2double (x_scm, 0.0);
942 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
945 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
947 "Return whether @var{sky} is empty.")
949 Skyline *s = Skyline::unsmob (sky);
950 LY_ASSERT_SMOB (Skyline, sky, 1);
951 return scm_from_bool (s->is_empty ());