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 (Direction sky)
463 empty_skyline (&buildings_);
467 Build skyline from a set of boxes.
469 Boxes should be non-empty on both axes. Otherwise, they will be ignored
471 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
473 list<Building> buildings;
476 for (vsize i = 0; i < boxes.size (); i++)
477 if (!boxes[i].is_empty (X_AXIS)
478 && !boxes[i].is_empty (Y_AXIS))
479 buildings.push_front (Building (boxes[i], horizon_axis, sky));
481 buildings_ = internal_build_skyline (&buildings);
486 build skyline from a set of line segments.
488 Segments can be articulated from left to right or right to left.
489 In the case of the latter, they will be stored internally as left to right.
491 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
493 list<Building> buildings;
496 for (vsize i = 0; i < segments.size (); i++)
498 Drul_array<Offset> const &seg = segments[i];
499 Offset left = seg[LEFT];
500 Offset right = seg[RIGHT];
501 if (left[horizon_axis] > right[horizon_axis])
504 Real x1 = left[horizon_axis];
505 Real x2 = right[horizon_axis];
506 Real y1 = left[other_axis (horizon_axis)] * sky;
507 Real y2 = right[other_axis (horizon_axis)] * sky;
510 buildings.push_back (Building (x1, y1, y2, x2));
513 buildings_ = internal_build_skyline (&buildings);
517 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
521 deque<Skyline> partials;
522 for (vsize i = 0; i < skypairs.size (); i++)
523 partials.push_back (Skyline ((skypairs[i])[sky]));
525 while (partials.size () > 1)
527 Skyline one = partials.front ();
528 partials.pop_front ();
529 Skyline two = partials.front ();
530 partials.pop_front ();
533 partials.push_back (one);
536 if (partials.size ())
537 buildings_.swap (partials.front ().buildings_);
542 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
545 if (!b.is_empty (X_AXIS) && !b.is_empty (Y_AXIS))
547 Building front (b, horizon_axis, sky);
548 single_skyline (front, &buildings_);
554 Skyline::merge (Skyline const &other)
556 assert (sky_ == other.sky_);
558 if (other.is_empty ())
563 buildings_ = other.buildings_;
567 list<Building> other_bld (other.buildings_);
568 list<Building> my_bld;
569 my_bld.splice (my_bld.begin (), buildings_);
570 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
575 Skyline::insert (Box const &b, Axis a)
577 list<Building> other_bld;
578 list<Building> my_bld;
580 if (isnan (b[other_axis (a)][LEFT])
581 || isnan (b[other_axis (a)][RIGHT]))
583 programming_error ("insane box for skyline");
587 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
588 if (b.is_empty (X_AXIS) || b.is_empty (Y_AXIS))
591 my_bld.splice (my_bld.begin (), buildings_);
592 single_skyline (Building (b, a, sky_), &other_bld);
593 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
598 Skyline::raise (Real r)
600 list<Building>::iterator end = buildings_.end ();
601 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
602 i->y_intercept_ += sky_ * r;
606 Skyline::shift (Real s)
608 list<Building>::iterator end = buildings_.end ();
609 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
613 i->y_intercept_ -= s * i->slope_;
618 Skyline::distance (Skyline const &other, Real horizon_padding) const
621 return internal_distance (other, horizon_padding, &dummy);
625 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
628 internal_distance (other, horizon_padding, &touch);
633 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
635 if (horizon_padding == 0.0)
636 return internal_distance (other, touch_point);
638 // Note that it is not necessary to build a padded version of other,
639 // because the same effect can be achieved just by doubling horizon_padding.
640 Skyline padded_this = padded (horizon_padding);
641 return padded_this.internal_distance (other, touch_point);
645 Skyline::padded (Real horizon_padding) const
647 if (horizon_padding < 0.0)
648 warning ("Cannot have negative horizon padding. Junking.");
650 if (horizon_padding <= 0.0)
653 list<Building> pad_buildings;
654 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
656 if (i->start_ > -infinity_f)
658 Real height = i->height (i->start_);
659 if (height > -infinity_f)
661 // Add the sloped building that pads the left side of the current building.
662 Real start = i->start_ - 2 * horizon_padding;
663 Real end = i->start_ - horizon_padding;
664 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
666 // Add the flat building that pads the left side of the current building.
667 start = i->start_ - horizon_padding;
669 pad_buildings.push_back (Building (start, height, height, end));
673 if (i->end_ < infinity_f)
675 Real height = i->height (i->end_);
676 if (height > -infinity_f)
678 // Add the flat building that pads the right side of the current building.
679 Real start = i->end_;
680 Real end = start + horizon_padding;
681 pad_buildings.push_back (Building (start, height, height, end));
683 // Add the sloped building that pads the right side of the current building.
685 end += horizon_padding;
686 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
691 // The buildings may be overlapping, so resolve that.
692 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
694 // Merge the padding with the original, to make a new skyline.
695 Skyline padded (sky_);
696 list<Building> my_buildings = buildings_;
697 padded.buildings_.clear ();
698 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
705 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
707 assert (sky_ == -other.sky_);
709 list<Building>::const_iterator i = buildings_.begin ();
710 list<Building>::const_iterator j = other.buildings_.begin ();
712 Real dist = -infinity_f;
713 Real start = -infinity_f;
714 Real touch = -infinity_f;
715 while (i != buildings_.end () && j != other.buildings_.end ())
717 Real end = min (i->end_, j->end_);
718 Real start_dist = i->height (start) + j->height (start);
719 Real end_dist = i->height (end) + j->height (end);
720 dist = max (dist, max (start_dist, end_dist));
722 if (end_dist == dist)
724 else if (start_dist == dist)
727 if (i->end_ <= j->end_)
734 *touch_point = touch;
739 Skyline::height (Real airplane) const
741 assert (!isinf (airplane));
743 list<Building>::const_iterator i;
744 for (i = buildings_.begin (); i != buildings_.end (); i++)
746 if (i->end_ >= airplane)
747 return sky_ * i->height (airplane);
755 Skyline::max_height () const
757 Real ret = -infinity_f;
759 list<Building>::const_iterator i;
760 for (i = buildings_.begin (); i != buildings_.end (); ++i)
762 ret = max (ret, i->height (i->start_));
763 ret = max (ret, i->height (i->end_));
770 Skyline::direction () const
776 Skyline::left () const
778 for (list<Building>::const_iterator i (buildings_.begin ());
779 i != buildings_.end (); i++)
780 if (i->y_intercept_ > -infinity_f)
787 Skyline::right () const
789 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
790 i != buildings_.rend (); ++i)
791 if (i->y_intercept_ > -infinity_f)
798 Skyline::max_height_position () const
801 s.set_minimum_height (0);
802 return touching_point (s);
806 Skyline::set_minimum_height (Real h)
809 s.buildings_.front ().y_intercept_ = h * sky_;
814 Skyline::to_points (Axis horizon_axis) const
818 Real start = -infinity_f;
819 for (list<Building>::const_iterator i (buildings_.begin ());
820 i != buildings_.end (); i++)
822 out.push_back (Offset (start, sky_ * i->height (start)));
823 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
827 if (horizon_axis == Y_AXIS)
828 for (vsize i = 0; i < out.size (); i++)
829 out[i] = out[i].swapped ();
835 Skyline::is_empty () const
837 if (!buildings_.size ())
839 Building b = buildings_.front ();
840 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
847 empty_skyline (&buildings_);
850 /****************************************************************/
852 IMPLEMENT_SIMPLE_SMOBS (Skyline);
853 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
854 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
857 Skyline::mark_smob (SCM s)
859 ASSERT_LIVE_IS_ALLOWED (s);
864 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
866 Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
869 scm_puts ("#<Skyline>", port);
874 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
876 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
878 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
880 Real horizon_padding = 0;
881 if (horizon_padding_scm != SCM_UNDEFINED)
883 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
884 horizon_padding = scm_to_double (horizon_padding_scm);
887 Skyline *skyline = Skyline::unsmob (skyline_scm);
888 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
889 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
892 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
894 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
896 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
898 Real horizon_padding = 0;
899 if (horizon_padding_scm != SCM_UNDEFINED)
901 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
902 horizon_padding = scm_to_double (horizon_padding_scm);
905 Skyline *skyline = Skyline::unsmob (skyline_scm);
906 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
907 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
910 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
912 Skyline::get_max_height (SCM skyline_scm)
914 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
917 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
919 Skyline::get_max_height_position (SCM skyline_scm)
921 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
924 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
926 Skyline::get_height (SCM skyline_scm, SCM x_scm)
928 Real x = robust_scm2double (x_scm, 0.0);
929 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
932 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
934 "Return whether @var{sky} is empty.")
936 Skyline *s = Skyline::unsmob (sky);
937 LY_ASSERT_SMOB (Skyline, sky, 1);
938 return scm_from_bool (s->is_empty ());