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;
129 y_intercept_ = start_height - slope_ * start;
133 Building::height (Real x) const
135 return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
139 Building::print () const
141 printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
145 Building::intersection_x (Building const &other) const
147 Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
148 return isnan (ret) ? -infinity_f : ret;
152 Building::leading_part (Real chop)
154 assert (chop <= end_);
158 // Returns a shift s such that (x + s, y) intersects the roof of
159 // this building. If no such shift exists, returns infinity_f.
161 Building::shift_to_intersect (Real x, Real y) const
163 // Solve for s: y = (x + s)*m + b
164 Real ret = (y - y_intercept_ - slope_ * x) / slope_;
166 if (ret >= start_ && ret <= end_ && !isinf (ret))
172 first_intersection (Building const &b, list<Building> *const s, Real start_x)
174 while (!s->empty () && start_x < b.end_)
176 Building c = s->front ();
178 // conceals and intersection_x involve multiplication and
179 // division. Avoid that, if we can.
180 if (c.y_intercept_ == -infinity_f)
189 if (c.conceals (b, start_x))
192 Real i = b.intersection_x (c);
193 if (i > start_x && i <= b.end_ && i <= c.end_)
204 Building::conceals (Building const &other, Real x) const
206 if (slope_ == other.slope_)
207 return y_intercept_ > other.y_intercept_;
209 /* their slopes were not equal, so there is an intersection point */
210 Real i = intersection_x (other);
211 return (i <= x && slope_ > other.slope_)
212 || (i > x && slope_ < other.slope_);
215 // Remove redundant empty buildings from the skyline.
216 // If there are two adjacent empty buildings, they can be
219 Skyline::normalize ()
221 bool last_empty = false;
222 list<Building>::iterator i;
224 for (i = buildings_.begin (); i != buildings_.end (); i++)
226 if (last_empty && i->y_intercept_ == -infinity_f)
228 list<Building>::iterator last = i;
230 last->end_ = i->end_;
231 buildings_.erase (i);
234 last_empty = (i->y_intercept_ == -infinity_f);
237 assert (buildings_.front ().start_ == -infinity_f);
238 assert (buildings_.back ().end_ == infinity_f);
244 // Since a skyline should always be normalized, we can
245 // assume that there are never two adjacent empty buildings.
246 // That is, if center is empty then left and right are not.
247 list<Building>::iterator left = buildings_.begin ();
248 list<Building>::iterator center = buildings_.begin ();
249 list<Building>::iterator right;
251 for (right = buildings_.begin (); right != buildings_.end (); right++)
253 if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
255 Real p1 = left->height (left->end_);
256 Real p2 = right->height (right->start_);
257 *center = Building (center->start_, p1, p2, center->end_);
266 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
267 list<Building> *const result) const
269 if (s1->empty () || s2->empty ())
271 programming_error ("tried to merge an empty skyline");
275 Real x = -infinity_f;
276 Real last_end = -infinity_f;
277 while (!s1->empty ())
279 if (s2->front ().conceals (s1->front (), x))
282 Building b = s1->front ();
283 Building c = s2->front ();
285 // Optimization: if the other skyline is empty at this point,
286 // we can avoid testing some intersections. Just grab as many
287 // buildings from s1 as we can, and shove them onto the output.
288 if (c.y_intercept_ == -infinity_f
291 list<Building>::iterator i = s1->begin ();
293 while (i != s1->end () && i->end_ <= c.end_)
296 s1->front ().start_ = x;
297 result->splice (result->end (), *s1, s1->begin (), i);
298 x = result->back ().end_;
303 Real end = first_intersection (b, s2, x);
307 result->push_back (b);
313 b.leading_part (end);
316 result->push_back (b);
319 if (end >= s1->front ().end_)
327 empty_skyline (list<Building> *const ret)
329 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
333 Given Building 'b', build a skyline containing only that building.
336 single_skyline (Building b, list<Building> *const ret)
338 assert (b.end_ >= b.start_);
340 if (b.start_ != -infinity_f)
341 ret->push_back (Building (-infinity_f, -infinity_f,
342 -infinity_f, b.start_));
344 if (b.end_ != infinity_f)
345 ret->push_back (Building (b.end_, -infinity_f,
346 -infinity_f, infinity_f));
349 /* remove a non-overlapping set of boxes from BOXES and build a skyline
351 static list<Building>
352 non_overlapping_skyline (list<Building> *const buildings)
354 list<Building> result;
355 Real last_end = -infinity_f;
356 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
357 list<Building>::iterator i = buildings->begin ();
358 while (i != buildings->end ())
361 Real y1 = i->height (i->start_);
363 Real y2 = i->height (i->end_);
365 // Drop buildings that will obviously have no effect.
366 if (last_building.height (x1) >= y1
367 && last_building.end_ >= x2
368 && last_building.height (x2) >= y2)
370 list<Building>::iterator j = i++;
371 buildings->erase (j);
381 // Insert empty Buildings into any gaps. (TODO: is this needed? -KOH)
383 result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
385 result.push_back (*i);
389 list<Building>::iterator j = i++;
390 buildings->erase (j);
393 if (last_end < infinity_f)
394 result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
398 class LessThanBuilding
401 bool operator () (Building const &b1, Building const &b2)
403 return b1.start_ < b2.start_
404 || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
409 BUILDINGS is a list of buildings, but they could be overlapping
410 and in any order. The returned list of buildings is ordered and non-overlapping.
413 Skyline::internal_build_skyline (list<Building> *buildings) const
415 vsize size = buildings->size ();
419 list<Building> result;
420 empty_skyline (&result);
425 list<Building> result;
426 single_skyline (buildings->front (), &result);
430 deque<list<Building> > partials;
431 buildings->sort (LessThanBuilding ());
432 while (!buildings->empty ())
433 partials.push_back (non_overlapping_skyline (buildings));
435 /* we'd like to say while (partials->size () > 1) but that's O (n).
436 Instead, we exit in the middle of the loop */
437 while (!partials.empty ())
439 list<Building> merged;
440 list<Building> one = partials.front ();
441 partials.pop_front ();
442 if (partials.empty ())
445 list<Building> two = partials.front ();
446 partials.pop_front ();
447 internal_merge_skyline (&one, &two, &merged);
448 partials.push_back (merged);
451 return list<Building> ();
457 empty_skyline (&buildings_);
460 Skyline::Skyline (Skyline const &src)
464 /* doesn't a list's copy constructor do this? -- jneem */
465 for (list<Building>::const_iterator i = src.buildings_.begin ();
466 i != src.buildings_.end (); i++)
468 buildings_.push_back (Building ((*i)));
472 Skyline::Skyline (Direction sky)
475 empty_skyline (&buildings_);
479 Build skyline from a set of boxes.
481 Boxes should be non-empty on both axes. Otherwise, they will be ignored
483 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
485 list<Building> buildings;
488 for (vsize i = 0; i < boxes.size (); i++)
489 if (!boxes[i].is_empty ())
490 buildings.push_front (Building (boxes[i], horizon_axis, sky));
492 buildings_ = internal_build_skyline (&buildings);
497 build skyline from a set of line segments.
499 Segments can be articulated from left to right or right to left.
500 In the case of the latter, they will be stored internally as left to right.
502 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
504 list<Building> buildings;
507 for (vsize i = 0; i < segments.size (); i++)
509 Drul_array<Offset> const &seg = segments[i];
510 Offset left = seg[LEFT];
511 Offset right = seg[RIGHT];
512 if (left[horizon_axis] > right[horizon_axis])
515 Real x1 = left[horizon_axis];
516 Real x2 = right[horizon_axis];
517 Real y1 = left[other_axis (horizon_axis)] * sky;
518 Real y2 = right[other_axis (horizon_axis)] * sky;
521 buildings.push_back (Building (x1, y1, y2, x2));
524 buildings_ = internal_build_skyline (&buildings);
528 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
532 deque<Skyline> partials;
533 for (vsize i = 0; i < skypairs.size (); i++)
534 partials.push_back (Skyline ((skypairs[i])[sky]));
536 while (partials.size () > 1)
538 Skyline one = partials.front ();
539 partials.pop_front ();
540 Skyline two = partials.front ();
541 partials.pop_front ();
544 partials.push_back (one);
547 if (partials.size ())
548 buildings_.swap (partials.front ().buildings_);
553 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
556 Building front (b, horizon_axis, sky);
557 single_skyline (front, &buildings_);
562 Skyline::merge (Skyline const &other)
564 assert (sky_ == other.sky_);
566 if (other.is_empty ())
571 buildings_ = other.buildings_;
575 list<Building> other_bld (other.buildings_);
576 list<Building> my_bld;
577 my_bld.splice (my_bld.begin (), buildings_);
578 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
583 Skyline::insert (Box const &b, Axis a)
585 list<Building> other_bld;
586 list<Building> my_bld;
588 if (isnan (b[other_axis (a)][LEFT])
589 || isnan (b[other_axis (a)][RIGHT]))
591 programming_error ("insane box for skyline");
595 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
599 my_bld.splice (my_bld.begin (), buildings_);
600 single_skyline (Building (b, a, sky_), &other_bld);
601 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
606 Skyline::raise (Real r)
608 list<Building>::iterator end = buildings_.end ();
609 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
610 i->y_intercept_ += sky_ * r;
614 Skyline::shift (Real s)
616 list<Building>::iterator end = buildings_.end ();
617 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
621 i->y_intercept_ -= s * i->slope_;
626 Skyline::distance (Skyline const &other, Real horizon_padding) const
629 return internal_distance (other, horizon_padding, &dummy);
633 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
636 internal_distance (other, horizon_padding, &touch);
641 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
643 if (horizon_padding == 0.0)
644 return internal_distance (other, touch_point);
646 // Note that it is not necessary to build a padded version of other,
647 // because the same effect can be achieved just by doubling horizon_padding.
648 Skyline padded_this = padded (horizon_padding);
649 return padded_this.internal_distance (other, touch_point);
653 Skyline::padded (Real horizon_padding) const
655 if (horizon_padding < 0.0)
656 warning ("Cannot have negative horizon padding. Junking.");
658 if (horizon_padding <= 0.0)
661 list<Building> pad_buildings;
662 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
664 if (i->start_ > -infinity_f)
666 Real height = i->height (i->start_);
667 if (height > -infinity_f)
669 // Add the sloped building that pads the left side of the current building.
670 Real start = i->start_ - 2 * horizon_padding;
671 Real end = i->start_ - horizon_padding;
672 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
674 // Add the flat building that pads the left side of the current building.
675 start = i->start_ - horizon_padding;
677 pad_buildings.push_back (Building (start, height, height, end));
681 if (i->end_ < infinity_f)
683 Real height = i->height (i->end_);
684 if (height > -infinity_f)
686 // Add the flat building that pads the right side of the current building.
687 Real start = i->end_;
688 Real end = start + horizon_padding;
689 pad_buildings.push_back (Building (start, height, height, end));
691 // Add the sloped building that pads the right side of the current building.
693 end += horizon_padding;
694 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
699 // The buildings may be overlapping, so resolve that.
700 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
702 // Merge the padding with the original, to make a new skyline.
703 Skyline padded (sky_);
704 list<Building> my_buildings = buildings_;
705 padded.buildings_.clear ();
706 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
713 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
715 assert (sky_ == -other.sky_);
717 list<Building>::const_iterator i = buildings_.begin ();
718 list<Building>::const_iterator j = other.buildings_.begin ();
720 Real dist = -infinity_f;
721 Real start = -infinity_f;
722 Real touch = -infinity_f;
723 while (i != buildings_.end () && j != other.buildings_.end ())
725 Real end = min (i->end_, j->end_);
726 Real start_dist = i->height (start) + j->height (start);
727 Real end_dist = i->height (end) + j->height (end);
728 dist = max (dist, max (start_dist, end_dist));
730 if (end_dist == dist)
732 else if (start_dist == dist)
735 if (i->end_ <= j->end_)
742 *touch_point = touch;
747 Skyline::height (Real airplane) const
749 assert (!isinf (airplane));
751 list<Building>::const_iterator i;
752 for (i = buildings_.begin (); i != buildings_.end (); i++)
754 if (i->end_ >= airplane)
755 return sky_ * i->height (airplane);
763 Skyline::max_height () const
765 Real ret = -infinity_f;
767 list<Building>::const_iterator i;
768 for (i = buildings_.begin (); i != buildings_.end (); ++i)
770 ret = max (ret, i->height (i->start_));
771 ret = max (ret, i->height (i->end_));
778 Skyline::direction () const
784 Skyline::left () const
786 for (list<Building>::const_iterator i (buildings_.begin ());
787 i != buildings_.end (); i++)
788 if (i->y_intercept_ > -infinity_f)
795 Skyline::right () const
797 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
798 i != buildings_.rend (); ++i)
799 if (i->y_intercept_ > -infinity_f)
806 Skyline::max_height_position () const
809 s.set_minimum_height (0);
810 return touching_point (s);
814 Skyline::set_minimum_height (Real h)
817 s.buildings_.front ().y_intercept_ = h * sky_;
822 Skyline::to_points (Axis horizon_axis) const
826 Real start = -infinity_f;
827 for (list<Building>::const_iterator i (buildings_.begin ());
828 i != buildings_.end (); i++)
830 out.push_back (Offset (start, sky_ * i->height (start)));
831 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
835 if (horizon_axis == Y_AXIS)
836 for (vsize i = 0; i < out.size (); i++)
837 out[i] = out[i].swapped ();
843 Skyline::is_empty () const
845 if (!buildings_.size ())
847 Building b = buildings_.front ();
848 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
855 empty_skyline (&buildings_);
858 /****************************************************************/
860 IMPLEMENT_SIMPLE_SMOBS (Skyline);
861 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
862 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
865 Skyline::mark_smob (SCM s)
867 ASSERT_LIVE_IS_ALLOWED (s);
872 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
874 Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
877 scm_puts ("#<Skyline>", port);
882 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
884 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
886 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
888 Real horizon_padding = 0;
889 if (horizon_padding_scm != SCM_UNDEFINED)
891 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
892 horizon_padding = scm_to_double (horizon_padding_scm);
895 Skyline *skyline = Skyline::unsmob (skyline_scm);
896 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
897 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
900 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
902 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
904 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
906 Real horizon_padding = 0;
907 if (horizon_padding_scm != SCM_UNDEFINED)
909 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
910 horizon_padding = scm_to_double (horizon_padding_scm);
913 Skyline *skyline = Skyline::unsmob (skyline_scm);
914 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
915 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
918 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
920 Skyline::get_max_height (SCM skyline_scm)
922 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
925 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
927 Skyline::get_max_height_position (SCM skyline_scm)
929 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
932 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
934 Skyline::get_height (SCM skyline_scm, SCM x_scm)
936 Real x = robust_scm2double (x_scm, 0.0);
937 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
940 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
942 "Return whether @var{sky} is empty.")
944 Skyline *s = Skyline::unsmob (sky);
945 LY_ASSERT_SMOB (Skyline, sky, 1);
946 return scm_from_bool (s->is_empty ());