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_));
124 assert (start_height == end_height);
126 start_height_ = start_height;
129 // Because of the way slope is calculated as a ratio of usually small
130 // differences, its precision is not sufficient for extrapolating
131 // significantly from the original interval. For that reason, we
132 // don't store the y0 value that would allow more precalculation by using
133 // y = x * slope + y0
135 // y = (x - xstart) * slope + ystart
138 Building::height (Real x) const
140 return (isinf (x) || isinf (start_)) ? start_height_
141 : slope_ * (x - start_) + start_height_;
145 Building::print () const
147 printf ("%f (x - x0) + %f from %f to %f\n", slope_, start_height_, start_, end_);
151 Building::intersection_x (Building const &other) const
154 - (other.start_height_ - start_height_
155 - other.slope_ * (other.start_ - start_))
156 /(other.slope_ - slope_);
157 /* A numerically symmetric version would be
158 x = (x1+x0)/2 - (y1-y0 - (s1+s0)(x1-x0)/2)/(s1-s0)
160 return isnan (ret) ? -infinity_f : ret;
164 Building::leading_part (Real chop)
166 assert (chop <= end_);
171 first_intersection (Building const &b, list<Building> *const s, Real start_x)
173 while (!s->empty () && start_x < b.end_)
175 Building c = s->front ();
177 // conceals and intersection_x involve multiplication and
178 // division. Avoid that, if we can.
179 if (c.start_height_ == -infinity_f)
188 if (c.conceals (b, start_x))
191 Real i = b.intersection_x (c);
192 if (i > start_x && i <= b.end_ && i <= c.end_)
202 // conceals returns true if *this is strictly above other at x
205 Building::conceals (Building const &other, Real x) const
207 return height (x) > other.height (x);
210 // Remove redundant empty buildings from the skyline.
211 // If there are two adjacent empty buildings, they can be
214 Skyline::normalize ()
216 bool last_empty = false;
217 list<Building>::iterator i;
219 for (i = buildings_.begin (); i != buildings_.end (); i++)
221 if (last_empty && i->start_height_ == -infinity_f)
223 list<Building>::iterator last = i;
225 last->end_ = i->end_;
226 buildings_.erase (i);
229 last_empty = (i->start_height_ == -infinity_f);
232 assert (buildings_.front ().start_ == -infinity_f);
233 assert (buildings_.back ().end_ == infinity_f);
239 // Since a skyline should always be normalized, we can
240 // assume that there are never two adjacent empty buildings.
241 // That is, if center is empty then left and right are not.
242 list<Building>::iterator left = buildings_.begin ();
243 list<Building>::iterator center = buildings_.begin ();
244 list<Building>::iterator right;
246 for (right = buildings_.begin (); right != buildings_.end (); right++)
248 if (center != buildings_.begin () && center->start_height_ == -infinity_f)
250 Real p1 = left->height (left->end_);
251 Real p2 = right->start_height_;
252 *center = Building (center->start_, p1, p2, center->end_);
261 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
262 list<Building> *const result) const
264 if (s1->empty () || s2->empty ())
266 programming_error ("tried to merge an empty skyline");
270 Real x = -infinity_f;
271 Real last_end = -infinity_f;
272 while (!s1->empty ())
274 if (s2->front ().conceals (s1->front (), x))
277 Building b = s1->front ();
278 Building c = s2->front ();
280 // Optimization: if the other skyline is empty at this point,
281 // we can avoid testing some intersections. Just grab as many
282 // buildings from s1 as we can, and shove them onto the output.
283 if (c.start_height_ == -infinity_f
286 list<Building>::iterator i = s1->begin ();
288 while (i != s1->end () && i->end_ <= c.end_)
291 s1->front ().start_height_ = s1->front ().height (x);
292 s1->front ().start_ = x;
293 result->splice (result->end (), *s1, s1->begin (), i);
294 x = result->back ().end_;
299 Real end = first_intersection (b, s2, x);
302 b.start_height_ = b.height (last_end);
304 result->push_back (b);
310 b.leading_part (end);
311 b.start_height_ = b.height (last_end);
314 result->push_back (b);
317 if (end >= s1->front ().end_)
325 empty_skyline (list<Building> *const ret)
327 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
331 Given Building 'b', build a skyline containing only that building.
334 single_skyline (Building b, list<Building> *const ret)
336 assert (b.end_ >= b.start_);
338 if (b.start_ != -infinity_f)
339 ret->push_back (Building (-infinity_f, -infinity_f,
340 -infinity_f, b.start_));
342 if (b.end_ != infinity_f)
343 ret->push_back (Building (b.end_, -infinity_f,
344 -infinity_f, infinity_f));
347 /* remove a non-overlapping set of boxes from BOXES and build a skyline
349 static list<Building>
350 non_overlapping_skyline (list<Building> *const buildings)
352 list<Building> result;
353 Real last_end = -infinity_f;
354 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
355 list<Building>::iterator i = buildings->begin ();
356 while (i != buildings->end ())
359 Real y1 = i->start_height_;
361 Real y2 = i->height (i->end_);
363 // Drop buildings that will obviously have no effect.
364 if (last_building.height (x1) >= y1
365 && last_building.end_ >= x2
366 && last_building.height (x2) >= y2)
368 list<Building>::iterator j = i++;
369 buildings->erase (j);
379 // Insert empty Buildings into any gaps. (TODO: is this needed? -KOH)
381 result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
383 result.push_back (*i);
387 list<Building>::iterator j = i++;
388 buildings->erase (j);
391 if (last_end < infinity_f)
392 result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
396 class LessThanBuilding
399 bool operator () (Building const &b1, Building const &b2)
401 return b1.start_ < b2.start_
402 || (b1.start_ == b2.start_ && b1.start_height_ > b2.start_height_);
407 BUILDINGS is a list of buildings, but they could be overlapping
408 and in any order. The returned list of buildings is ordered and non-overlapping.
411 Skyline::internal_build_skyline (list<Building> *buildings) const
413 vsize size = buildings->size ();
417 list<Building> result;
418 empty_skyline (&result);
423 list<Building> result;
424 single_skyline (buildings->front (), &result);
428 deque<list<Building> > partials;
429 buildings->sort (LessThanBuilding ());
430 while (!buildings->empty ())
431 partials.push_back (non_overlapping_skyline (buildings));
433 /* we'd like to say while (partials->size () > 1) but that's O (n).
434 Instead, we exit in the middle of the loop */
435 while (!partials.empty ())
437 list<Building> merged;
438 list<Building> one = partials.front ();
439 partials.pop_front ();
440 if (partials.empty ())
443 list<Building> two = partials.front ();
444 partials.pop_front ();
445 internal_merge_skyline (&one, &two, &merged);
446 partials.push_back (merged);
449 return list<Building> ();
455 empty_skyline (&buildings_);
458 Skyline::Skyline (Skyline const &src)
462 /* doesn't a list's copy constructor do this? -- jneem */
463 for (list<Building>::const_iterator i = src.buildings_.begin ();
464 i != src.buildings_.end (); i++)
466 buildings_.push_back (Building ((*i)));
470 Skyline::Skyline (Direction sky)
473 empty_skyline (&buildings_);
477 Build skyline from a set of boxes.
479 Boxes should be non-empty on both axes. Otherwise, they will be ignored
481 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
483 list<Building> buildings;
486 for (vsize i = 0; i < boxes.size (); i++)
487 if (!boxes[i].is_empty (X_AXIS)
488 && !boxes[i].is_empty (Y_AXIS))
489 buildings.push_front (Building (boxes[i], horizon_axis, sky));
491 buildings_ = internal_build_skyline (&buildings);
496 build skyline from a set of line segments.
498 Segments can be articulated from left to right or right to left.
499 In the case of the latter, they will be stored internally as left to right.
501 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
503 list<Building> buildings;
506 for (vsize i = 0; i < segments.size (); i++)
508 Drul_array<Offset> const &seg = segments[i];
509 Offset left = seg[LEFT];
510 Offset right = seg[RIGHT];
511 if (left[horizon_axis] > right[horizon_axis])
514 Real x1 = left[horizon_axis];
515 Real x2 = right[horizon_axis];
516 Real y1 = left[other_axis (horizon_axis)] * sky;
517 Real y2 = right[other_axis (horizon_axis)] * sky;
520 buildings.push_back (Building (x1, y1, y2, x2));
523 buildings_ = internal_build_skyline (&buildings);
527 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
531 deque<Skyline> partials;
532 for (vsize i = 0; i < skypairs.size (); i++)
533 partials.push_back (Skyline ((skypairs[i])[sky]));
535 while (partials.size () > 1)
537 Skyline one = partials.front ();
538 partials.pop_front ();
539 Skyline two = partials.front ();
540 partials.pop_front ();
543 partials.push_back (one);
546 if (partials.size ())
547 buildings_.swap (partials.front ().buildings_);
552 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
555 if (!b.is_empty (X_AXIS) && !b.is_empty (Y_AXIS))
557 Building front (b, horizon_axis, sky);
558 single_skyline (front, &buildings_);
564 Skyline::merge (Skyline const &other)
566 assert (sky_ == other.sky_);
568 if (other.is_empty ())
573 buildings_ = other.buildings_;
577 list<Building> other_bld (other.buildings_);
578 list<Building> my_bld;
579 my_bld.splice (my_bld.begin (), buildings_);
580 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
585 Skyline::insert (Box const &b, Axis a)
587 list<Building> other_bld;
588 list<Building> my_bld;
590 if (isnan (b[other_axis (a)][LEFT])
591 || isnan (b[other_axis (a)][RIGHT]))
593 programming_error ("insane box for skyline");
597 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
598 if (b.is_empty (X_AXIS) || b.is_empty (Y_AXIS))
601 my_bld.splice (my_bld.begin (), buildings_);
602 single_skyline (Building (b, a, sky_), &other_bld);
603 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
608 Skyline::raise (Real r)
610 list<Building>::iterator end = buildings_.end ();
611 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
612 i->start_height_ += sky_ * r;
616 Skyline::shift (Real s)
618 list<Building>::iterator end = buildings_.end ();
619 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
627 Skyline::distance (Skyline const &other, Real horizon_padding) const
630 return internal_distance (other, horizon_padding, &dummy);
634 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
637 internal_distance (other, horizon_padding, &touch);
642 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
644 if (horizon_padding == 0.0)
645 return internal_distance (other, touch_point);
647 // Note that it is not necessary to build a padded version of other,
648 // because the same effect can be achieved just by doubling horizon_padding.
649 Skyline padded_this = padded (horizon_padding);
650 return padded_this.internal_distance (other, touch_point);
654 Skyline::padded (Real horizon_padding) const
656 if (horizon_padding < 0.0)
657 warning ("Cannot have negative horizon padding. Junking.");
659 if (horizon_padding <= 0.0)
662 list<Building> pad_buildings;
663 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
665 if (i->start_ > -infinity_f)
667 Real height = i->start_height_;
668 if (height > -infinity_f)
670 // Add the sloped building that pads the left side of the current building.
671 Real start = i->start_ - 2 * horizon_padding;
672 Real end = i->start_ - horizon_padding;
673 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
675 // Add the flat building that pads the left side of the current building.
676 start = i->start_ - horizon_padding;
678 pad_buildings.push_back (Building (start, height, height, end));
682 if (i->end_ < infinity_f)
684 Real height = i->height (i->end_);
685 if (height > -infinity_f)
687 // Add the flat building that pads the right side of the current building.
688 Real start = i->end_;
689 Real end = start + horizon_padding;
690 pad_buildings.push_back (Building (start, height, height, end));
692 // Add the sloped building that pads the right side of the current building.
694 end += horizon_padding;
695 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
700 // The buildings may be overlapping, so resolve that.
701 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
703 // Merge the padding with the original, to make a new skyline.
704 Skyline padded (sky_);
705 list<Building> my_buildings = buildings_;
706 padded.buildings_.clear ();
707 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
714 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
716 assert (sky_ == -other.sky_);
718 list<Building>::const_iterator i = buildings_.begin ();
719 list<Building>::const_iterator j = other.buildings_.begin ();
721 Real dist = -infinity_f;
722 Real start = -infinity_f;
723 Real touch = -infinity_f;
724 while (i != buildings_.end () && j != other.buildings_.end ())
726 Real end = min (i->end_, j->end_);
727 Real start_dist = i->height (start) + j->height (start);
728 Real end_dist = i->height (end) + j->height (end);
729 dist = max (dist, max (start_dist, end_dist));
731 if (end_dist == dist)
733 else if (start_dist == dist)
736 if (i->end_ <= j->end_)
743 *touch_point = touch;
748 Skyline::height (Real airplane) const
750 assert (!isinf (airplane));
752 list<Building>::const_iterator i;
753 for (i = buildings_.begin (); i != buildings_.end (); i++)
755 if (i->end_ >= airplane)
756 return sky_ * i->height (airplane);
764 Skyline::max_height () const
766 Real ret = -infinity_f;
768 list<Building>::const_iterator i;
769 for (i = buildings_.begin (); i != buildings_.end (); ++i)
771 ret = max (ret, i->start_height_);
772 ret = max (ret, i->height (i->end_));
779 Skyline::direction () const
785 Skyline::left () const
787 for (list<Building>::const_iterator i (buildings_.begin ());
788 i != buildings_.end (); i++)
789 if (i->start_height_ > -infinity_f)
796 Skyline::right () const
798 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
799 i != buildings_.rend (); ++i)
800 if (i->start_height_ > -infinity_f)
807 Skyline::max_height_position () const
810 s.set_minimum_height (0);
811 return touching_point (s);
815 Skyline::set_minimum_height (Real h)
818 s.buildings_.front ().start_height_ = h * sky_;
823 Skyline::to_points (Axis horizon_axis) const
827 Real start = -infinity_f;
828 for (list<Building>::const_iterator i (buildings_.begin ());
829 i != buildings_.end (); i++)
831 out.push_back (Offset (start, sky_ * i->height (start)));
832 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
836 if (horizon_axis == Y_AXIS)
837 for (vsize i = 0; i < out.size (); i++)
838 out[i] = out[i].swapped ();
844 Skyline::is_empty () const
846 if (!buildings_.size ())
848 Building b = buildings_.front ();
849 return b.end_ == infinity_f && b.start_height_ == -infinity_f;
856 empty_skyline (&buildings_);
859 /****************************************************************/
861 IMPLEMENT_SIMPLE_SMOBS (Skyline);
862 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
863 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
866 Skyline::mark_smob (SCM s)
868 ASSERT_LIVE_IS_ALLOWED (s);
873 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
875 Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
878 scm_puts ("#<Skyline>", port);
883 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
885 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
887 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
889 Real horizon_padding = 0;
890 if (horizon_padding_scm != SCM_UNDEFINED)
892 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
893 horizon_padding = scm_to_double (horizon_padding_scm);
896 Skyline *skyline = Skyline::unsmob (skyline_scm);
897 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
898 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
901 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
903 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
905 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
907 Real horizon_padding = 0;
908 if (horizon_padding_scm != SCM_UNDEFINED)
910 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
911 horizon_padding = scm_to_double (horizon_padding_scm);
914 Skyline *skyline = Skyline::unsmob (skyline_scm);
915 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
916 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
919 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
921 Skyline::get_max_height (SCM skyline_scm)
923 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
926 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
928 Skyline::get_max_height_position (SCM skyline_scm)
930 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
933 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
935 Skyline::get_height (SCM skyline_scm, SCM x_scm)
937 Real x = robust_scm2double (x_scm, 0.0);
938 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
941 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
943 "Return whether @var{sky} is empty.")
945 Skyline *s = Skyline::unsmob (sky);
946 LY_ASSERT_SMOB (Skyline, sky, 1);
947 return scm_from_bool (s->is_empty ());