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 Be careful about numerical accuracy. When dealing with extremely small values,
59 computation errors may arise due to the use of floating point arithmetic.
60 For example, if left and right have equal values to start with, in C++
61 they may not receive the same value after
63 left = left*factor + offset;
64 right = right*factor + offset;
66 Which is very unfortunate. Maybe using GCC compiler options to disallow
67 extended precision for intermediate results and/or the choice to store
68 intermediates with less than full precision would retain some kind of
69 deterministic behavior that way.
71 Anyway, it seems that accepting extremely narrow building in skylines
72 doesn't cause accuracy problems to us, so we allow arbitrarily small buildings.
73 However, as Keith pointed out, problems may appear if more than one operation
74 is done before storing the result, and/or there are different code paths
75 for doing the operations to the different ends of an interval.
79 print_buildings (list<Building> const &b)
81 for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
86 Skyline::print () const
88 print_buildings (buildings_);
92 Skyline::print_points () const
94 vector<Offset> ps (to_points (X_AXIS));
96 for (vsize i = 0; i < ps.size (); i++)
97 printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
98 (i % 2) == 1 ? "\n" : " ");
101 Building::Building (Real start, Real start_height, Real end_height, Real end)
103 if (isinf (start) || isinf (end))
104 assert (start_height == end_height);
108 precompute (start, start_height, end_height, end);
111 Building::Building (Box const &b, Axis horizon_axis, Direction sky)
113 Real start = b[horizon_axis][LEFT];
114 Real end = b[horizon_axis][RIGHT];
115 Real height = sky * b[other_axis (horizon_axis)][sky];
119 precompute (start, height, height, end);
123 Building::precompute (Real start, Real start_height, Real end_height, Real end)
125 slope_ = 0.0; /* if they were both infinite, we would get nan, not 0, from the prev line */
126 if (start_height != end_height)
127 slope_ = (end_height - start_height) / (end - start);
129 assert (!isinf (slope_) && !isnan (slope_));
133 assert (start_height == end_height);
134 y_intercept_ = start_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;
160 Building::leading_part (Real chop)
162 assert (chop <= end_);
166 // Returns a shift s such that (x + s, y) intersects the roof of
167 // this building. If no such shift exists, returns infinity_f.
169 Building::shift_to_intersect (Real x, Real y) const
171 // Solve for s: y = (x + s)*m + b
172 Real ret = (y - y_intercept_ - slope_ * x) / slope_;
174 if (ret >= start_ && ret <= end_ && !isinf (ret))
180 first_intersection (Building const &b, list<Building> *const s, Real start_x)
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_;
311 Real end = first_intersection (b, s2, x);
315 result->push_back (b);
321 b.leading_part (end);
324 result->push_back (b);
327 if (end >= s1->front ().end_)
335 empty_skyline (list<Building> *const ret)
337 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
341 Given Building 'b', build a skyline containing only that building.
344 single_skyline (Building b, list<Building> *const ret)
346 assert (b.end_ >= b.start_);
348 if (b.start_ != -infinity_f)
349 ret->push_back (Building (-infinity_f, -infinity_f,
350 -infinity_f, b.start_));
352 if (b.end_ != infinity_f)
353 ret->push_back (Building (b.end_, -infinity_f,
354 -infinity_f, infinity_f));
357 /* remove a non-overlapping set of boxes from BOXES and build a skyline
359 static list<Building>
360 non_overlapping_skyline (list<Building> *const buildings)
362 list<Building> result;
363 Real last_end = -infinity_f;
364 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
365 list<Building>::iterator i = buildings->begin ();
366 while (i != buildings->end ())
369 Real y1 = i->height (i->start_);
371 Real y2 = i->height (i->end_);
373 // Drop buildings that will obviously have no effect.
374 if (last_building.height (x1) >= y1
375 && last_building.end_ >= x2
376 && last_building.height (x2) >= y2)
378 list<Building>::iterator j = i++;
379 buildings->erase (j);
390 result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
392 result.push_back (*i);
396 list<Building>::iterator j = i++;
397 buildings->erase (j);
400 if (last_end < infinity_f)
401 result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
405 class LessThanBuilding
408 bool operator () (Building const &b1, Building const &b2)
410 return b1.start_ < b2.start_
411 || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
416 BUILDINGS is a list of buildings, but they could be overlapping
417 and in any order. The returned list of buildings is ordered and non-overlapping.
420 Skyline::internal_build_skyline (list<Building> *buildings) const
422 vsize size = buildings->size ();
426 list<Building> result;
427 empty_skyline (&result);
432 list<Building> result;
433 single_skyline (buildings->front (), &result);
437 deque<list<Building> > partials;
438 buildings->sort (LessThanBuilding ());
439 while (!buildings->empty ())
440 partials.push_back (non_overlapping_skyline (buildings));
442 /* we'd like to say while (partials->size () > 1) but that's O (n).
443 Instead, we exit in the middle of the loop */
444 while (!partials.empty ())
446 list<Building> merged;
447 list<Building> one = partials.front ();
448 partials.pop_front ();
449 if (partials.empty ())
452 list<Building> two = partials.front ();
453 partials.pop_front ();
454 internal_merge_skyline (&one, &two, &merged);
455 partials.push_back (merged);
458 return list<Building> ();
464 empty_skyline (&buildings_);
467 Skyline::Skyline (Skyline const &src)
471 /* doesn't a list's copy constructor do this? -- jneem */
472 for (list<Building>::const_iterator i = src.buildings_.begin ();
473 i != src.buildings_.end (); i++)
475 buildings_.push_back (Building ((*i)));
479 Skyline::Skyline (Direction sky)
482 empty_skyline (&buildings_);
486 Build skyline from a set of boxes.
488 Boxes should be non-empty on both axes. Otherwise, they will be ignored
490 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
492 list<Building> buildings;
495 for (vsize i = 0; i < boxes.size (); i++)
496 if (!boxes[i].is_empty ())
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 Building front (b, horizon_axis, sky);
564 single_skyline (front, &buildings_);
569 Skyline::merge (Skyline const &other)
571 assert (sky_ == other.sky_);
573 if (other.is_empty ())
578 buildings_ = other.buildings_;
582 list<Building> other_bld (other.buildings_);
583 list<Building> my_bld;
584 my_bld.splice (my_bld.begin (), buildings_);
585 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
590 Skyline::insert (Box const &b, Axis a)
592 list<Building> other_bld;
593 list<Building> my_bld;
595 if (isnan (b[other_axis (a)][LEFT])
596 || isnan (b[other_axis (a)][RIGHT]))
598 programming_error ("insane box for skyline");
602 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
606 my_bld.splice (my_bld.begin (), buildings_);
607 single_skyline (Building (b, a, sky_), &other_bld);
608 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
613 Skyline::raise (Real r)
615 list<Building>::iterator end = buildings_.end ();
616 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
617 i->y_intercept_ += sky_ * r;
621 Skyline::shift (Real s)
623 list<Building>::iterator end = buildings_.end ();
624 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
628 i->y_intercept_ -= s * i->slope_;
633 Skyline::distance (Skyline const &other, Real horizon_padding) const
636 return internal_distance (other, horizon_padding, &dummy);
640 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
643 internal_distance (other, horizon_padding, &touch);
648 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
650 if (horizon_padding == 0.0)
651 return internal_distance (other, touch_point);
653 // Note that it is not necessary to build a padded version of other,
654 // because the same effect can be achieved just by doubling horizon_padding.
655 Skyline padded_this = padded (horizon_padding);
656 return padded_this.internal_distance (other, touch_point);
660 Skyline::padded (Real horizon_padding) const
662 if (horizon_padding < 0.0)
663 warning ("Cannot have negative horizon padding. Junking.");
665 if (horizon_padding <= 0.0)
668 list<Building> pad_buildings;
669 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
671 if (i->start_ > -infinity_f)
673 Real height = i->height (i->start_);
674 if (height > -infinity_f)
676 // Add the sloped building that pads the left side of the current building.
677 Real start = i->start_ - 2 * horizon_padding;
678 Real end = i->start_ - horizon_padding;
679 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
681 // Add the flat building that pads the left side of the current building.
682 start = i->start_ - horizon_padding;
684 pad_buildings.push_back (Building (start, height, height, end));
688 if (i->end_ < infinity_f)
690 Real height = i->height (i->end_);
691 if (height > -infinity_f)
693 // Add the flat building that pads the right side of the current building.
694 Real start = i->end_;
695 Real end = start + horizon_padding;
696 pad_buildings.push_back (Building (start, height, height, end));
698 // Add the sloped building that pads the right side of the current building.
700 end += horizon_padding;
701 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
706 // The buildings may be overlapping, so resolve that.
707 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
709 // Merge the padding with the original, to make a new skyline.
710 Skyline padded (sky_);
711 list<Building> my_buildings = buildings_;
712 padded.buildings_.clear ();
713 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
720 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
722 assert (sky_ == -other.sky_);
724 list<Building>::const_iterator i = buildings_.begin ();
725 list<Building>::const_iterator j = other.buildings_.begin ();
727 Real dist = -infinity_f;
728 Real start = -infinity_f;
729 Real touch = -infinity_f;
730 while (i != buildings_.end () && j != other.buildings_.end ())
732 Real end = min (i->end_, j->end_);
733 Real start_dist = i->height (start) + j->height (start);
734 Real end_dist = i->height (end) + j->height (end);
735 dist = max (dist, max (start_dist, end_dist));
737 if (end_dist == dist)
739 else if (start_dist == dist)
742 if (i->end_ <= j->end_)
749 *touch_point = touch;
754 Skyline::height (Real airplane) const
756 assert (!isinf (airplane));
758 list<Building>::const_iterator i;
759 for (i = buildings_.begin (); i != buildings_.end (); i++)
761 if (i->end_ >= airplane)
762 return sky_ * i->height (airplane);
770 Skyline::max_height () const
772 Real ret = -infinity_f;
774 list<Building>::const_iterator i;
775 for (i = buildings_.begin (); i != buildings_.end (); ++i)
777 ret = max (ret, i->height (i->start_));
778 ret = max (ret, i->height (i->end_));
785 Skyline::direction () const
791 Skyline::left () const
793 for (list<Building>::const_iterator i (buildings_.begin ());
794 i != buildings_.end (); i++)
795 if (i->y_intercept_ > -infinity_f)
802 Skyline::right () const
804 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
805 i != buildings_.rend (); ++i)
806 if (i->y_intercept_ > -infinity_f)
813 Skyline::max_height_position () const
816 s.set_minimum_height (0);
817 return touching_point (s);
821 Skyline::set_minimum_height (Real h)
824 s.buildings_.front ().y_intercept_ = h * sky_;
829 Skyline::to_points (Axis horizon_axis) const
833 Real start = -infinity_f;
834 for (list<Building>::const_iterator i (buildings_.begin ());
835 i != buildings_.end (); i++)
837 out.push_back (Offset (start, sky_ * i->height (start)));
838 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
842 if (horizon_axis == Y_AXIS)
843 for (vsize i = 0; i < out.size (); i++)
844 out[i] = out[i].swapped ();
850 Skyline::is_empty () const
852 if (!buildings_.size ())
854 Building b = buildings_.front ();
855 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
862 empty_skyline (&buildings_);
865 /****************************************************************/
867 IMPLEMENT_SIMPLE_SMOBS (Skyline);
868 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
869 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
872 Skyline::mark_smob (SCM s)
874 ASSERT_LIVE_IS_ALLOWED (s);
879 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
881 Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
884 scm_puts ("#<Skyline>", port);
889 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
891 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
893 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
895 Real horizon_padding = 0;
896 if (horizon_padding_scm != SCM_UNDEFINED)
898 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
899 horizon_padding = scm_to_double (horizon_padding_scm);
902 Skyline *skyline = Skyline::unsmob (skyline_scm);
903 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
904 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
907 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
909 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
911 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
913 Real horizon_padding = 0;
914 if (horizon_padding_scm != SCM_UNDEFINED)
916 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
917 horizon_padding = scm_to_double (horizon_padding_scm);
920 Skyline *skyline = Skyline::unsmob (skyline_scm);
921 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
922 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
925 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
927 Skyline::get_max_height (SCM skyline_scm)
929 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
932 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
934 Skyline::get_max_height_position (SCM skyline_scm)
936 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
939 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
941 Skyline::get_height (SCM skyline_scm, SCM x_scm)
943 Real x = robust_scm2double (x_scm, 0.0);
944 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
947 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
949 "Return whether @var{sky} is empty.")
951 Skyline *s = Skyline::unsmob (sky);
952 LY_ASSERT_SMOB (Skyline, sky, 1);
953 return scm_from_bool (s->is_empty ());