2 This file is part of LilyPond, the GNU music typesetter.
4 Copyright (C) 2006--2014 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"
26 /* A skyline is a sequence of non-overlapping buildings: something like
36 Each building has a starting position, and ending position, a starting
37 height and an ending height.
39 The following invariants are observed:
40 - the start of the first building is at -infinity
41 - the end of the last building is at infinity
42 - if a building has infinite length (ie. the first and last buildings),
43 then its starting height and ending height are equal
44 - the end of one building is the same as the beginning of the next
47 We also allow skylines to point down (the structure is exactly the same,
48 but we think of the part above the line as being filled with mass and the
49 part below as being empty). ::distance finds the minimum distance between
50 an UP skyline and a DOWN skyline.
52 Note that we store DOWN skylines upside-down. That is, in order to compare
53 a DOWN skyline with an UP skyline, we need to flip the DOWN skyline first.
54 This means that the merging routine doesn't need to be aware of direction,
55 but the distance routine does.
57 From 2007 through 2012, buildings of width less than EPS were discarded,
58 citing numerical accuracy concerns. We remember that floating point
59 comparisons of nearly-equal values can be affected by rounding error.
60 Also, some target machines use the x87 floating point unit, which provides
61 extended precision for intermediate results held in registers. On this type
62 of hardware comparisons such as
63 double c = 1.0/3.0; boolean compare = (c == 1.0/3.0)
64 could go either way because the 1.0/3.0 is allowed to be kept
65 higher precision than the variable 'c'.
66 Alert to these considerations, we now accept buildings of zero-width.
70 print_buildings (list<Building> const &b)
72 for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
77 Skyline::print () const
79 print_buildings (buildings_);
83 Skyline::print_points () const
85 vector<Offset> ps (to_points (X_AXIS));
87 for (vsize i = 0; i < ps.size (); i++)
88 printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
89 (i % 2) == 1 ? "\n" : " ");
92 Building::Building (Real start, Real start_height, Real end_height, Real end)
94 if (isinf (start) || isinf (end))
95 assert (start_height == end_height);
99 precompute (start, start_height, end_height, end);
102 Building::Building (Box const &b, Axis horizon_axis, Direction sky)
104 Real start = b[horizon_axis][LEFT];
105 Real end = b[horizon_axis][RIGHT];
106 Real height = sky * b[other_axis (horizon_axis)][sky];
110 precompute (start, height, height, end);
114 Building::precompute (Real start, Real start_height, Real end_height, Real end)
116 slope_ = 0.0; /* if they were both infinite, we would get nan, not 0, from the prev line */
117 if (start_height != end_height)
118 slope_ = (end_height - start_height) / (end - start);
120 assert (!isinf (slope_) && !isnan (slope_));
124 assert (start_height == end_height);
125 y_intercept_ = start_height;
127 else if (fabs(slope_) > 1e6)
128 // too steep to be stored in slope-intercept form, given round-off error
131 y_intercept_ = max(start_height, end_height);
134 y_intercept_ = start_height - slope_ * start;
138 Building::height (Real x) const
140 return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
144 Building::print () const
146 printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
150 Building::intersection_x (Building const &other) const
152 Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
153 return isnan (ret) ? -infinity_f : ret;
157 Building::leading_part (Real chop)
159 assert (chop <= end_);
163 // Returns a shift s such that (x + s, y) intersects the roof of
164 // this building. If no such shift exists, returns infinity_f.
166 Building::shift_to_intersect (Real x, Real y) const
168 // Solve for s: y = (x + s)*m + b
169 Real ret = (y - y_intercept_ - slope_ * x) / slope_;
171 if (ret >= start_ && ret <= end_ && !isinf (ret))
177 first_intersection (Building const &b, list<Building> *s, Real start_x)
178 /* Return the first x >= start_x where skyline s above Building b.
179 * Removes buildings from s that are concealed by b. */
181 while (!s->empty () && start_x < b.end_)
183 Building c = s->front ();
185 // conceals and intersection_x involve multiplication and
186 // division. Avoid that, if we can.
187 if (c.y_intercept_ == -infinity_f)
196 if (c.conceals (b, start_x))
199 Real i = b.intersection_x (c);
200 if (i > start_x && i <= b.end_ && i <= c.end_)
211 Building::conceals (Building const &other, Real x) const
213 if (slope_ == other.slope_)
214 return y_intercept_ > other.y_intercept_;
216 /* their slopes were not equal, so there is an intersection point */
217 Real i = intersection_x (other);
218 return (i <= x && slope_ > other.slope_)
219 || (i > x && slope_ < other.slope_);
222 // Remove redundant empty buildings from the skyline.
223 // If there are two adjacent empty buildings, they can be
226 Skyline::normalize ()
228 bool last_empty = false;
229 list<Building>::iterator i;
231 for (i = buildings_.begin (); i != buildings_.end (); i++)
233 if (last_empty && i->y_intercept_ == -infinity_f)
235 list<Building>::iterator last = i;
237 last->end_ = i->end_;
238 buildings_.erase (i);
241 last_empty = (i->y_intercept_ == -infinity_f);
244 assert (buildings_.front ().start_ == -infinity_f);
245 assert (buildings_.back ().end_ == infinity_f);
251 // Since a skyline should always be normalized, we can
252 // assume that there are never two adjacent empty buildings.
253 // That is, if center is empty then left and right are not.
254 list<Building>::iterator left = buildings_.begin ();
255 list<Building>::iterator center = buildings_.begin ();
256 list<Building>::iterator right;
258 for (right = buildings_.begin (); right != buildings_.end (); right++)
260 if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
262 Real p1 = left->height (left->end_);
263 Real p2 = right->height (right->start_);
264 *center = Building (center->start_, p1, p2, center->end_);
273 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
274 list<Building> *const result) const
276 if (s1->empty () || s2->empty ())
278 programming_error ("tried to merge an empty skyline");
282 Real x = -infinity_f;
283 Real last_end = -infinity_f;
284 while (!s1->empty ())
286 if (s2->front ().conceals (s1->front (), x))
289 Building b = s1->front ();
290 Building c = s2->front ();
292 // Optimization: if the other skyline is empty at this point,
293 // we can avoid testing some intersections. Just grab as many
294 // buildings from s1 as we can, and shove them onto the output.
295 if (c.y_intercept_ == -infinity_f
298 list<Building>::iterator i = s1->begin ();
300 while (i != s1->end () && i->end_ <= c.end_)
303 s1->front ().start_ = x;
304 result->splice (result->end (), *s1, s1->begin (), i);
305 x = result->back ().end_;
309 // first_intersection() removes buildings from s2 if b hides them
310 Real end = first_intersection (b, s2, x);
314 result->push_back (b);
318 // Should be (end > x), during ver2.19. end == x happens fairly often,
319 // and we do not need to keep vertical segments within a skyline.
322 b.leading_part (end);
325 result->push_back (b);
328 if (end >= s1->front ().end_)
330 // Should add during ver2.19 (to avoid an endless loop
331 // when merging identical skylines with a vertical segment)
332 // if (end >= s2->front().end_) s2->pop_front();
339 empty_skyline (list<Building> *const ret)
341 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
345 Given Building 'b', build a skyline containing only that building.
348 single_skyline (Building b, list<Building> *const ret)
350 assert (b.end_ >= b.start_);
352 if (b.start_ != -infinity_f)
353 ret->push_back (Building (-infinity_f, -infinity_f,
354 -infinity_f, b.start_));
356 if (b.end_ != infinity_f)
357 ret->push_back (Building (b.end_, -infinity_f,
358 -infinity_f, infinity_f));
361 /* remove a non-overlapping set of boxes from BOXES and build a skyline
363 static list<Building>
364 non_overlapping_skyline (list<Building> *const buildings)
366 list<Building> result;
367 Real last_end = -infinity_f;
368 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
369 list<Building>::iterator i = buildings->begin ();
370 while (i != buildings->end ())
373 Real y1 = i->height (i->start_);
375 Real y2 = i->height (i->end_);
377 // Drop buildings that will obviously have no effect.
378 if (last_building.height (x1) >= y1
379 && last_building.end_ >= x2
380 && last_building.height (x2) >= y2)
382 list<Building>::iterator j = i++;
383 buildings->erase (j);
393 // Insert empty Buildings into any gaps. (TODO: is this needed? -KOH)
395 result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
397 result.push_back (*i);
401 list<Building>::iterator j = i++;
402 buildings->erase (j);
405 if (last_end < infinity_f)
406 result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
410 class LessThanBuilding
413 bool operator () (Building const &b1, Building const &b2)
415 return b1.start_ < b2.start_
416 || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
421 BUILDINGS is a list of buildings, but they could be overlapping
422 and in any order. The returned list of buildings is ordered and non-overlapping.
425 Skyline::internal_build_skyline (list<Building> *buildings) const
427 vsize size = buildings->size ();
431 list<Building> result;
432 empty_skyline (&result);
437 list<Building> result;
438 single_skyline (buildings->front (), &result);
442 deque<list<Building> > partials;
443 buildings->sort (LessThanBuilding ());
444 while (!buildings->empty ())
445 partials.push_back (non_overlapping_skyline (buildings));
447 /* we'd like to say while (partials->size () > 1) but that's O (n).
448 Instead, we exit in the middle of the loop */
449 while (!partials.empty ())
451 list<Building> merged;
452 list<Building> one = partials.front ();
453 partials.pop_front ();
454 if (partials.empty ())
457 list<Building> two = partials.front ();
458 partials.pop_front ();
459 internal_merge_skyline (&one, &two, &merged);
460 partials.push_back (merged);
463 return list<Building> ();
469 empty_skyline (&buildings_);
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 const char Skyline::type_p_name_[] = "ly:skyline?";
866 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
868 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
870 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
872 Real horizon_padding = 0;
873 if (horizon_padding_scm != SCM_UNDEFINED)
875 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
876 horizon_padding = scm_to_double (horizon_padding_scm);
879 Skyline *skyline = Skyline::unsmob (skyline_scm);
880 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
881 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
884 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
886 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
888 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
890 Real horizon_padding = 0;
891 if (horizon_padding_scm != SCM_UNDEFINED)
893 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
894 horizon_padding = scm_to_double (horizon_padding_scm);
897 Skyline *skyline = Skyline::unsmob (skyline_scm);
898 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
899 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
902 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
904 Skyline::get_max_height (SCM skyline_scm)
906 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
909 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
911 Skyline::get_max_height_position (SCM skyline_scm)
913 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
916 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
918 Skyline::get_height (SCM skyline_scm, SCM x_scm)
920 Real x = robust_scm2double (x_scm, 0.0);
921 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
924 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
926 "Return whether @var{sky} is empty.")
928 Skyline *s = Skyline::unsmob (sky);
929 LY_ASSERT_SMOB (Skyline, sky, 1);
930 return scm_from_bool (s->is_empty ());