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.
59 /* If we start including very thin buildings, numerical accuracy errors can
60 arise. Therefore, we ignore all buildings that are less than epsilon wide. */
64 print_buildings (list<Building> const &b)
66 for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
71 Skyline::print () const
73 print_buildings (buildings_);
77 Skyline::print_points () const
79 vector<Offset> ps (to_points (X_AXIS));
81 for (vsize i = 0; i < ps.size (); i++)
82 printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
83 (i % 2) == 1 ? "\n" : " ");
86 Building::Building (Real start, Real start_height, Real end_height, Real end)
88 if (isinf (start) || isinf (end))
89 assert (start_height == end_height);
93 precompute (start, start_height, end_height, end);
96 Building::Building (Box const &b, Axis horizon_axis, Direction sky)
98 Real start = b[horizon_axis][LEFT];
99 Real end = b[horizon_axis][RIGHT];
100 Real height = sky * b[other_axis (horizon_axis)][sky];
104 precompute (start, height, height, end);
108 Building::precompute (Real start, Real start_height, Real end_height, Real end)
110 slope_ = 0.0; /* if they were both infinite, we would get nan, not 0, from the prev line */
111 if (start_height != end_height)
112 slope_ = (end_height - start_height) / (end - start);
114 assert (!isinf (slope_) && !isnan (slope_));
118 assert (start_height == end_height);
119 y_intercept_ = start_height;
122 y_intercept_ = start_height - slope_ * start;
126 Building::height (Real x) const
128 return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
132 Building::print () const
134 printf ("%f x + %f from %f to %f\n", slope_, y_intercept_, start_, end_);
138 Building::intersection_x (Building const &other) const
140 Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
141 return isnan (ret) ? -infinity_f : ret;
145 Building::leading_part (Real chop)
147 assert (chop <= end_);
151 // Returns a shift s such that (x + s, y) intersects the roof of
152 // this building. If no such shift exists, returns infinity_f.
154 Building::shift_to_intersect (Real x, Real y) const
156 // Solve for s: y = (x + s)*m + b
157 Real ret = (y - y_intercept_ - slope_ * x) / slope_;
159 if (ret >= start_ && ret <= end_ && !isinf (ret))
165 first_intersection (Building const &b, list<Building> *const s, Real start_x)
167 while (!s->empty () && start_x < b.end_)
169 Building c = s->front ();
171 // conceals and intersection_x involve multiplication and
172 // division. Avoid that, if we can.
173 if (c.y_intercept_ == -infinity_f)
182 if (c.conceals (b, start_x))
185 Real i = b.intersection_x (c);
186 if (i > start_x && i <= b.end_ && i <= c.end_)
197 Building::conceals (Building const &other, Real x) const
199 if (slope_ == other.slope_)
200 return y_intercept_ > other.y_intercept_;
202 /* their slopes were not equal, so there is an intersection point */
203 Real i = intersection_x (other);
204 return (i <= x && slope_ > other.slope_)
205 || (i > x && slope_ < other.slope_);
208 // Remove redundant empty buildings from the skyline.
209 // If there are two adjacent empty buildings, they can be
212 Skyline::normalize ()
214 bool last_empty = false;
215 list<Building>::iterator i;
217 for (i = buildings_.begin (); i != buildings_.end (); i++)
219 if (last_empty && i->y_intercept_ == -infinity_f)
221 list<Building>::iterator last = i;
223 last->end_ = i->end_;
224 buildings_.erase (i);
227 last_empty = (i->y_intercept_ == -infinity_f);
230 assert (buildings_.front ().start_ == -infinity_f);
231 assert (buildings_.back ().end_ == infinity_f);
237 // Since a skyline should always be normalized, we can
238 // assume that there are never two adjacent empty buildings.
239 // That is, if center is empty then left and right are not.
240 list<Building>::iterator left = buildings_.begin ();
241 list<Building>::iterator center = buildings_.begin ();
242 list<Building>::iterator right;
244 for (right = buildings_.begin (); right != buildings_.end (); right++)
246 if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
248 Real p1 = left->height (left->end_);
249 Real p2 = right->height (right->start_);
250 *center = Building (center->start_, p1, p2, center->end_);
259 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
260 list<Building> *const result) const
262 if (s1->empty () || s2->empty ())
264 programming_error ("tried to merge an empty skyline");
268 Real x = -infinity_f;
269 Real last_end = -infinity_f;
270 while (!s1->empty ())
272 if (s2->front ().conceals (s1->front (), x))
275 Building b = s1->front ();
276 Building c = s2->front ();
278 // Optimization: if the other skyline is empty at this point,
279 // we can avoid testing some intersections. Just grab as many
280 // buildings from s1 as we can, and shove them onto the output.
281 if (c.y_intercept_ == -infinity_f
284 list<Building>::iterator i = s1->begin ();
286 while (i != s1->end () && i->end_ <= c.end_)
289 s1->front ().start_ = x;
290 result->splice (result->end (), *s1, s1->begin (), i);
291 x = result->back ().end_;
296 Real end = first_intersection (b, s2, x);
300 result->push_back (b);
304 /* only include buildings wider than epsilon */
307 b.leading_part (end);
310 result->push_back (b);
313 if (end >= s1->front ().end_)
321 empty_skyline (list<Building> *const ret)
323 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
327 Given Building 'b', build a skyline containing only that building.
330 single_skyline (Building b, list<Building> *const ret)
332 if (b.end_ > b.start_ + EPS)
334 if (b.start_ != -infinity_f)
335 ret->push_back (Building (-infinity_f, -infinity_f,
336 -infinity_f, b.start_));
338 if (b.end_ != infinity_f)
339 ret->push_back (Building (b.end_, -infinity_f,
340 -infinity_f, infinity_f));
348 /* remove a non-overlapping set of boxes from BOXES and build a skyline
350 static list<Building>
351 non_overlapping_skyline (list<Building> *const buildings)
353 list<Building> result;
354 Real last_end = -infinity_f;
355 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
356 list<Building>::iterator i = buildings->begin ();
357 while (i != buildings->end ())
360 Real y1 = i->height (i->start_);
362 Real y2 = i->height (i->end_);
364 // Drop buildings that will obviously have no effect.
365 if (last_building.height (x1) >= y1
366 && last_building.end_ >= x2
367 && last_building.height (x2) >= y2)
369 list<Building>::iterator j = i++;
370 buildings->erase (j);
380 if (x1 > last_end + EPS)
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.height (b1.start_) > b2.height (b1.start_));
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 have fatness in the horizon_axis, otherwise they are ignored.
481 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
483 list<Building> buildings;
486 Axis vert_axis = other_axis (horizon_axis);
487 for (vsize i = 0; i < boxes.size (); i++)
489 Interval iv = boxes[i][horizon_axis];
490 if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
491 buildings.push_front (Building (boxes[i], horizon_axis, sky));
494 buildings_ = internal_build_skyline (&buildings);
499 build skyline from a set of line segments.
501 Buildings should have fatness in the horizon_axis, otherwise they are ignored.
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 Building front (b, horizon_axis, sky);
558 single_skyline (front, &buildings_);
563 Skyline::merge (Skyline const &other)
565 assert (sky_ == other.sky_);
567 if (other.is_empty ())
572 buildings_ = other.buildings_;
576 list<Building> other_bld (other.buildings_);
577 list<Building> my_bld;
578 my_bld.splice (my_bld.begin (), buildings_);
579 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
584 Skyline::insert (Box const &b, Axis a)
586 list<Building> other_bld;
587 list<Building> my_bld;
589 if (isnan (b[other_axis (a)][LEFT])
590 || isnan (b[other_axis (a)][RIGHT]))
592 programming_error ("insane box for skyline");
596 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
598 if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
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->y_intercept_ += 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++)
623 i->y_intercept_ -= s * i->slope_;
628 Skyline::distance (Skyline const &other, Real horizon_padding) const
631 return internal_distance (other, horizon_padding, &dummy);
635 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
638 internal_distance (other, horizon_padding, &touch);
643 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
645 if (horizon_padding == 0.0)
646 return internal_distance (other, touch_point);
648 // Note that it is not necessary to build a padded version of other,
649 // because the same effect can be achieved just by doubling horizon_padding.
650 Skyline padded_this = padded (horizon_padding);
651 return padded_this.internal_distance (other, touch_point);
655 Skyline::padded (Real horizon_padding) const
657 list<Building> pad_buildings;
658 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
660 if (i->start_ > -infinity_f)
662 Real height = i->height (i->start_);
663 if (height > -infinity_f)
665 // Add the sloped building that pads the left side of the current building.
666 Real start = i->start_ - 2 * horizon_padding;
667 Real end = i->start_ - horizon_padding;
668 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
670 // Add the flat building that pads the left side of the current building.
671 start = i->start_ - horizon_padding;
673 pad_buildings.push_back (Building (start, height, height, end));
677 if (i->end_ < infinity_f)
679 Real height = i->height (i->end_);
680 if (height > -infinity_f)
682 // Add the flat building that pads the right side of the current building.
683 Real start = i->end_;
684 Real end = start + horizon_padding;
685 pad_buildings.push_back (Building (start, height, height, end));
687 // Add the sloped building that pads the right side of the current building.
689 end += horizon_padding;
690 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
695 // The buildings may be overlapping, so resolve that.
696 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
698 // Merge the padding with the original, to make a new skyline.
699 Skyline padded (sky_);
700 list<Building> my_buildings = buildings_;
701 padded.buildings_.clear ();
702 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
709 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
711 assert (sky_ == -other.sky_);
713 list<Building>::const_iterator i = buildings_.begin ();
714 list<Building>::const_iterator j = other.buildings_.begin ();
716 Real dist = -infinity_f;
717 Real start = -infinity_f;
718 Real touch = -infinity_f;
719 while (i != buildings_.end () && j != other.buildings_.end ())
721 Real end = min (i->end_, j->end_);
722 Real start_dist = i->height (start) + j->height (start);
723 Real end_dist = i->height (end) + j->height (end);
724 dist = max (dist, max (start_dist, end_dist));
726 if (end_dist == dist)
728 else if (start_dist == dist)
731 if (i->end_ <= j->end_)
738 *touch_point = touch;
743 Skyline::height (Real airplane) const
745 assert (!isinf (airplane));
747 list<Building>::const_iterator i;
748 for (i = buildings_.begin (); i != buildings_.end (); i++)
750 if (i->end_ >= airplane)
751 return sky_ * i->height (airplane);
759 Skyline::max_height () const
761 Real ret = -infinity_f;
763 list<Building>::const_iterator i;
764 for (i = buildings_.begin (); i != buildings_.end (); ++i)
766 ret = max (ret, i->height (i->start_));
767 ret = max (ret, i->height (i->end_));
774 Skyline::direction () const
780 Skyline::left () const
782 for (list<Building>::const_iterator i (buildings_.begin ());
783 i != buildings_.end (); i++)
784 if (i->y_intercept_ > -infinity_f)
791 Skyline::right () const
793 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
794 i != buildings_.rend (); ++i)
795 if (i->y_intercept_ > -infinity_f)
802 Skyline::max_height_position () const
805 s.set_minimum_height (0);
806 return touching_point (s);
810 Skyline::set_minimum_height (Real h)
813 s.buildings_.front ().y_intercept_ = h * sky_;
818 Skyline::to_points (Axis horizon_axis) const
822 Real start = -infinity_f;
823 for (list<Building>::const_iterator i (buildings_.begin ());
824 i != buildings_.end (); i++)
826 out.push_back (Offset (start, sky_ * i->height (start)));
827 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
831 if (horizon_axis == Y_AXIS)
832 for (vsize i = 0; i < out.size (); i++)
833 out[i] = out[i].swapped ();
839 Skyline::is_empty () const
841 if (!buildings_.size ())
843 Building b = buildings_.front ();
844 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
851 empty_skyline (&buildings_);
854 /****************************************************************/
856 IMPLEMENT_SIMPLE_SMOBS (Skyline);
857 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
858 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
861 Skyline::mark_smob (SCM s)
863 ASSERT_LIVE_IS_ALLOWED (s);
868 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
870 Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
873 scm_puts ("#<Skyline>", port);
878 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
880 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
882 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
884 Real horizon_padding = 0;
885 if (horizon_padding_scm != SCM_UNDEFINED)
887 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
888 horizon_padding = scm_to_double (horizon_padding_scm);
891 Skyline *skyline = Skyline::unsmob (skyline_scm);
892 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
893 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
896 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
898 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
900 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
902 Real horizon_padding = 0;
903 if (horizon_padding_scm != SCM_UNDEFINED)
905 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
906 horizon_padding = scm_to_double (horizon_padding_scm);
909 Skyline *skyline = Skyline::unsmob (skyline_scm);
910 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
911 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
914 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
916 Skyline::get_max_height (SCM skyline_scm)
918 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
921 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
923 Skyline::get_max_height_position (SCM skyline_scm)
925 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
928 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
930 Skyline::get_height (SCM skyline_scm, SCM x_scm)
932 Real x = robust_scm2double (x_scm, 0.0);
933 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
936 LY_DEFINE (ly_skyline_empty_p, "ly:skyline-empty?",
938 "Return whether @var{sky} is empty.")
940 Skyline *s = Skyline::unsmob (sky);
941 LY_ASSERT_SMOB (Skyline, sky, 1);
942 return scm_from_bool (s->is_empty ());