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 (X_AXIS)
490 && !boxes[i].is_empty (Y_AXIS))
491 buildings.push_front (Building (boxes[i], horizon_axis, sky));
493 buildings_ = internal_build_skyline (&buildings);
498 build skyline from a set of line segments.
500 Segments can be articulated from left to right or right to left.
501 In the case of the latter, they will be stored internally as left to right.
503 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
505 list<Building> buildings;
508 for (vsize i = 0; i < segments.size (); i++)
510 Drul_array<Offset> const &seg = segments[i];
511 Offset left = seg[LEFT];
512 Offset right = seg[RIGHT];
513 if (left[horizon_axis] > right[horizon_axis])
516 Real x1 = left[horizon_axis];
517 Real x2 = right[horizon_axis];
518 Real y1 = left[other_axis (horizon_axis)] * sky;
519 Real y2 = right[other_axis (horizon_axis)] * sky;
522 buildings.push_back (Building (x1, y1, y2, x2));
525 buildings_ = internal_build_skyline (&buildings);
529 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
533 deque<Skyline> partials;
534 for (vsize i = 0; i < skypairs.size (); i++)
535 partials.push_back (Skyline ((skypairs[i])[sky]));
537 while (partials.size () > 1)
539 Skyline one = partials.front ();
540 partials.pop_front ();
541 Skyline two = partials.front ();
542 partials.pop_front ();
545 partials.push_back (one);
548 if (partials.size ())
549 buildings_.swap (partials.front ().buildings_);
554 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
557 if (!b.is_empty (X_AXIS) && !b.is_empty (Y_AXIS))
559 Building front (b, horizon_axis, sky);
560 single_skyline (front, &buildings_);
566 Skyline::merge (Skyline const &other)
568 assert (sky_ == other.sky_);
570 if (other.is_empty ())
575 buildings_ = other.buildings_;
579 list<Building> other_bld (other.buildings_);
580 list<Building> my_bld;
581 my_bld.splice (my_bld.begin (), buildings_);
582 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
587 Skyline::insert (Box const &b, Axis a)
589 list<Building> other_bld;
590 list<Building> my_bld;
592 if (isnan (b[other_axis (a)][LEFT])
593 || isnan (b[other_axis (a)][RIGHT]))
595 programming_error ("insane box for skyline");
599 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
600 if (b.is_empty (X_AXIS) || b.is_empty (Y_AXIS))
603 my_bld.splice (my_bld.begin (), buildings_);
604 single_skyline (Building (b, a, sky_), &other_bld);
605 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
610 Skyline::raise (Real r)
612 list<Building>::iterator end = buildings_.end ();
613 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
614 i->y_intercept_ += sky_ * r;
618 Skyline::shift (Real s)
620 list<Building>::iterator end = buildings_.end ();
621 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
625 i->y_intercept_ -= s * i->slope_;
630 Skyline::distance (Skyline const &other, Real horizon_padding) const
633 return internal_distance (other, horizon_padding, &dummy);
637 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
640 internal_distance (other, horizon_padding, &touch);
645 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
647 if (horizon_padding == 0.0)
648 return internal_distance (other, touch_point);
650 // Note that it is not necessary to build a padded version of other,
651 // because the same effect can be achieved just by doubling horizon_padding.
652 Skyline padded_this = padded (horizon_padding);
653 return padded_this.internal_distance (other, touch_point);
657 Skyline::padded (Real horizon_padding) const
659 if (horizon_padding < 0.0)
660 warning ("Cannot have negative horizon padding. Junking.");
662 if (horizon_padding <= 0.0)
665 list<Building> pad_buildings;
666 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
668 if (i->start_ > -infinity_f)
670 Real height = i->height (i->start_);
671 if (height > -infinity_f)
673 // Add the sloped building that pads the left side of the current building.
674 Real start = i->start_ - 2 * horizon_padding;
675 Real end = i->start_ - horizon_padding;
676 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
678 // Add the flat building that pads the left side of the current building.
679 start = i->start_ - horizon_padding;
681 pad_buildings.push_back (Building (start, height, height, end));
685 if (i->end_ < infinity_f)
687 Real height = i->height (i->end_);
688 if (height > -infinity_f)
690 // Add the flat building that pads the right side of the current building.
691 Real start = i->end_;
692 Real end = start + horizon_padding;
693 pad_buildings.push_back (Building (start, height, height, end));
695 // Add the sloped building that pads the right side of the current building.
697 end += horizon_padding;
698 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
703 // The buildings may be overlapping, so resolve that.
704 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
706 // Merge the padding with the original, to make a new skyline.
707 Skyline padded (sky_);
708 list<Building> my_buildings = buildings_;
709 padded.buildings_.clear ();
710 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
717 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
719 assert (sky_ == -other.sky_);
721 list<Building>::const_iterator i = buildings_.begin ();
722 list<Building>::const_iterator j = other.buildings_.begin ();
724 Real dist = -infinity_f;
725 Real start = -infinity_f;
726 Real touch = -infinity_f;
727 while (i != buildings_.end () && j != other.buildings_.end ())
729 Real end = min (i->end_, j->end_);
730 Real start_dist = i->height (start) + j->height (start);
731 Real end_dist = i->height (end) + j->height (end);
732 dist = max (dist, max (start_dist, end_dist));
734 if (end_dist == dist)
736 else if (start_dist == dist)
739 if (i->end_ <= j->end_)
746 *touch_point = touch;
751 Skyline::height (Real airplane) const
753 assert (!isinf (airplane));
755 list<Building>::const_iterator i;
756 for (i = buildings_.begin (); i != buildings_.end (); i++)
758 if (i->end_ >= airplane)
759 return sky_ * i->height (airplane);
767 Skyline::max_height () const
769 Real ret = -infinity_f;
771 list<Building>::const_iterator i;
772 for (i = buildings_.begin (); i != buildings_.end (); ++i)
774 ret = max (ret, i->height (i->start_));
775 ret = max (ret, i->height (i->end_));
782 Skyline::direction () const
788 Skyline::left () const
790 for (list<Building>::const_iterator i (buildings_.begin ());
791 i != buildings_.end (); i++)
792 if (i->y_intercept_ > -infinity_f)
799 Skyline::right () const
801 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
802 i != buildings_.rend (); ++i)
803 if (i->y_intercept_ > -infinity_f)
810 Skyline::max_height_position () const
813 s.set_minimum_height (0);
814 return touching_point (s);
818 Skyline::set_minimum_height (Real h)
821 s.buildings_.front ().y_intercept_ = h * sky_;
826 Skyline::to_points (Axis horizon_axis) const
830 Real start = -infinity_f;
831 for (list<Building>::const_iterator i (buildings_.begin ());
832 i != buildings_.end (); i++)
834 out.push_back (Offset (start, sky_ * i->height (start)));
835 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
839 if (horizon_axis == Y_AXIS)
840 for (vsize i = 0; i < out.size (); i++)
841 out[i] = out[i].swapped ();
847 Skyline::is_empty () const
849 if (!buildings_.size ())
851 Building b = buildings_.front ();
852 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
859 empty_skyline (&buildings_);
862 /****************************************************************/
864 IMPLEMENT_SIMPLE_SMOBS (Skyline);
865 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
866 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
869 Skyline::mark_smob (SCM s)
871 ASSERT_LIVE_IS_ALLOWED (s);
876 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
878 Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
881 scm_puts ("#<Skyline>", port);
886 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
888 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
890 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
892 Real horizon_padding = 0;
893 if (horizon_padding_scm != SCM_UNDEFINED)
895 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
896 horizon_padding = scm_to_double (horizon_padding_scm);
899 Skyline *skyline = Skyline::unsmob (skyline_scm);
900 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
901 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
904 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
906 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
908 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
910 Real horizon_padding = 0;
911 if (horizon_padding_scm != SCM_UNDEFINED)
913 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
914 horizon_padding = scm_to_double (horizon_padding_scm);
917 Skyline *skyline = Skyline::unsmob (skyline_scm);
918 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
919 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
922 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
924 Skyline::get_max_height (SCM skyline_scm)
926 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
929 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
931 Skyline::get_max_height_position (SCM skyline_scm)
933 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
936 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
938 Skyline::get_height (SCM skyline_scm, SCM x_scm)
940 Real x = robust_scm2double (x_scm, 0.0);
941 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
944 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
946 "Return whether @var{sky} is empty.")
948 Skyline *s = Skyline::unsmob (sky);
949 LY_ASSERT_SMOB (Skyline, sky, 1);
950 return scm_from_bool (s->is_empty ());