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))
164 // Returns the interval of horizontal shifts for which this
165 // building (pointing up) overlaps the other building (pointing down).
167 Building::overlapping_shift_interval (Building const &other) const
171 // If one building is empty, there will never be an overlap.
172 if (y_intercept_ == -infinity_f || other.y_intercept_ == -infinity_f)
175 // There are two kinds of interesting positions:
176 // - when the horizontal extents of the buildings just touch
177 // - when an endpoint of one building intersects the roof of the other.
178 // The interval we are looking for is the smallest one that
179 // contains all of the interesting points.
182 Real my_y1 = height (start_);
183 Real my_y2 = height (end_);
184 Real his_y1 = -other.height (other.start_); // "-" because OTHER points down
185 Real his_y2 = -other.height (other.end_);
187 // If both buildings are infinite in the same direction,
188 // the return value is either empty or full.
189 if ((isinf (start_) && isinf (other.start_))
190 || (isinf (end_) && isinf (other.end_)))
191 return (y_intercept_ > other.y_intercept_)
192 ? Interval (-infinity_f, infinity_f) : Interval ();
194 // ...when the horizontal extents of the buildings just touch...
196 iv.add_point (other.end_ - start_);
198 iv.add_point (other.start_ - end_);
200 // ...when an endpoint of one building intersects the roof of the other.
201 Real p1 = shift_to_intersect (other.start_, his_y1);
202 Real p2 = shift_to_intersect (other.end_, his_y2);
203 // "-my_y1" because OTHER points down:
204 Real p3 = other.shift_to_intersect (start_, -my_y1);
205 Real p4 = other.shift_to_intersect (end_, -my_y2);
219 first_intersection (Building const &b, list<Building> *const s, Real start_x)
221 while (!s->empty () && start_x < b.end_)
223 Building c = s->front ();
225 // conceals and intersection_x involve multiplication and
226 // division. Avoid that, if we can.
227 if (c.y_intercept_ == -infinity_f)
236 if (c.conceals (b, start_x))
239 Real i = b.intersection_x (c);
240 if (i > start_x && i <= b.end_ && i <= c.end_)
251 Building::conceals (Building const &other, Real x) const
253 if (slope_ == other.slope_)
254 return y_intercept_ > other.y_intercept_;
256 /* their slopes were not equal, so there is an intersection point */
257 Real i = intersection_x (other);
258 return (i <= x && slope_ > other.slope_)
259 || (i > x && slope_ < other.slope_);
262 // Remove redundant empty buildings from the skyline.
263 // If there are two adjacent empty buildings, they can be
266 Skyline::normalize ()
268 bool last_empty = false;
269 list<Building>::iterator i;
270 for (i = buildings_.begin (); i != buildings_.end (); i++)
272 if (last_empty && i->y_intercept_ == -infinity_f)
274 list<Building>::iterator last = i;
276 last->end_ = i->end_;
277 buildings_.erase (i);
280 last_empty = (i->y_intercept_ == -infinity_f);
283 assert (buildings_.front ().start_ == -infinity_f);
284 assert (buildings_.back ().end_ == infinity_f);
290 // Since a skyline should always be normalized, we can
291 // assume that there are never two adjacent empty buildings.
292 // That is, if center is empty then left and right are not.
293 list<Building>::iterator left = buildings_.begin ();
294 list<Building>::iterator center = buildings_.begin ();
295 list<Building>::iterator right;
297 for (right = buildings_.begin (); right != buildings_.end (); right++)
299 if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
301 Real p1 = left->height (left->end_);
302 Real p2 = right->height (right->start_);
303 *center = Building (center->start_, p1, p2, center->end_);
312 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
313 list<Building> *const result) const
315 if (s1->empty () || s2->empty ())
317 programming_error ("tried to merge an empty skyline");
321 Real x = -infinity_f;
322 Real last_end = -infinity_f;
323 while (!s1->empty ())
325 if (s2->front ().conceals (s1->front (), x))
328 Building b = s1->front ();
329 Building c = s2->front ();
331 // Optimization: if the other skyline is empty at this point,
332 // we can avoid testing some intersections. Just grab as many
333 // buildings from s1 as we can, and shove them onto the output.
334 if (c.y_intercept_ == -infinity_f
337 list<Building>::iterator i = s1->begin ();
339 while (i != s1->end () && i->end_ <= c.end_)
342 s1->front ().start_ = x;
343 result->splice (result->end (), *s1, s1->begin (), i);
344 x = result->back ().end_;
349 Real end = first_intersection (b, s2, x);
353 result->push_back (b);
357 /* only include buildings wider than epsilon */
360 b.leading_part (end);
363 result->push_back (b);
366 if (end >= s1->front ().end_)
374 empty_skyline (list<Building> *const ret)
376 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
380 Given Building 'b', build a skyline containing only that building.
383 single_skyline (Building b, list<Building> *const ret)
385 if (b.end_ > b.start_ + EPS)
387 ret->push_back (Building (-infinity_f, -infinity_f,
388 -infinity_f, b.start_));
390 ret->push_back (Building (b.end_, -infinity_f,
391 -infinity_f, infinity_f));
399 /* remove a non-overlapping set of boxes from BOXES and build a skyline
401 static list<Building>
402 non_overlapping_skyline (list<Building> *const buildings)
404 list<Building> result;
405 Real last_end = -infinity_f;
406 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
407 list<Building>::iterator i = buildings->begin ();
408 while (i != buildings->end ())
411 Real y1 = i->height (i->start_);
413 Real y2 = i->height (i->end_);
415 // Drop buildings that will obviously have no effect.
416 if (last_building.height (x1) >= y1
417 && last_building.end_ >= x2
418 && last_building.height (x2) >= y2)
420 list<Building>::iterator j = i++;
421 buildings->erase (j);
431 if (x1 > last_end + EPS)
432 result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
434 result.push_back (*i);
438 list<Building>::iterator j = i++;
439 buildings->erase (j);
442 if (last_end < infinity_f)
443 result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
447 class LessThanBuilding
450 bool operator () (Building const &b1, Building const &b2)
452 return b1.start_ < b2.start_
453 || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
458 BUILDINGS is a list of buildings, but they could be overlapping
459 and in any order. The returned list of buildings is ordered and non-overlapping.
462 Skyline::internal_build_skyline (list<Building> *buildings) const
464 vsize size = buildings->size ();
468 list<Building> result;
469 empty_skyline (&result);
474 list<Building> result;
475 single_skyline (buildings->front (), &result);
479 deque<list<Building> > partials;
480 buildings->sort (LessThanBuilding ());
481 while (!buildings->empty ())
482 partials.push_back (non_overlapping_skyline (buildings));
484 /* we'd like to say while (partials->size () > 1) but that's O (n).
485 Instead, we exit in the middle of the loop */
486 while (!partials.empty ())
488 list<Building> merged;
489 list<Building> one = partials.front ();
490 partials.pop_front ();
491 if (partials.empty ())
494 list<Building> two = partials.front ();
495 partials.pop_front ();
496 internal_merge_skyline (&one, &two, &merged);
497 partials.push_back (merged);
500 return list<Building> ();
506 empty_skyline (&buildings_);
509 Skyline::Skyline (Skyline const &src)
513 /* doesn't a list's copy constructor do this? -- jneem */
514 for (list<Building>::const_iterator i = src.buildings_.begin ();
515 i != src.buildings_.end (); i++)
517 buildings_.push_back (Building ((*i)));
521 Skyline::Skyline (Direction sky)
524 empty_skyline (&buildings_);
528 Build skyline from a set of boxes.
530 Boxes should have fatness in the horizon_axis, otherwise they are ignored.
532 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
534 list<Building> buildings;
537 Axis vert_axis = other_axis (horizon_axis);
538 for (vsize i = 0; i < boxes.size (); i++)
540 Interval iv = boxes[i][horizon_axis];
541 if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
542 buildings.push_front (Building (boxes[i], horizon_axis, sky));
545 buildings_ = internal_build_skyline (&buildings);
550 build skyline from a set of line segments.
552 Buildings should have fatness in the horizon_axis, otherwise they are ignored.
554 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
556 list<Building> buildings;
559 for (vsize i = 0; i < segments.size (); i++)
561 Drul_array<Offset> const &seg = segments[i];
562 Offset left = seg[LEFT];
563 Offset right = seg[RIGHT];
564 if (left[horizon_axis] > right[horizon_axis])
567 Real x1 = left[horizon_axis];
568 Real x2 = right[horizon_axis];
569 Real y1 = left[other_axis (horizon_axis)] * sky;
570 Real y2 = right[other_axis (horizon_axis)] * sky;
573 buildings.push_back (Building (x1, y1, y2, x2));
576 buildings_ = internal_build_skyline (&buildings);
580 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
584 deque<Skyline> partials;
585 for (vsize i = 0; i < skypairs.size (); i++)
586 partials.push_back (Skyline ((skypairs[i])[sky]));
588 while (partials.size () > 1)
590 Skyline one = partials.front ();
591 partials.pop_front ();
592 Skyline two = partials.front ();
593 partials.pop_front ();
596 partials.push_back (one);
599 if (partials.size ())
600 buildings_.swap (partials.front ().buildings_);
605 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
608 Building front (b, horizon_axis, sky);
609 single_skyline (front, &buildings_);
613 Skyline::merge (Skyline const &other)
615 assert (sky_ == other.sky_);
617 if (other.is_empty ())
622 buildings_ = other.buildings_;
626 list<Building> other_bld (other.buildings_);
627 list<Building> my_bld;
628 my_bld.splice (my_bld.begin (), buildings_);
629 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
634 Skyline::insert (Box const &b, Axis a)
636 list<Building> other_bld;
637 list<Building> my_bld;
639 if (isnan (b[other_axis (a)][LEFT])
640 || isnan (b[other_axis (a)][RIGHT]))
642 programming_error ("insane box for skyline");
646 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
648 if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
651 my_bld.splice (my_bld.begin (), buildings_);
652 single_skyline (Building (b, a, sky_), &other_bld);
653 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
658 Skyline::raise (Real r)
660 list<Building>::iterator end = buildings_.end ();
661 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
662 i->y_intercept_ += sky_ * r;
666 Skyline::shift (Real s)
668 list<Building>::iterator end = buildings_.end ();
669 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
673 i->y_intercept_ -= s * i->slope_;
678 Skyline::distance (Skyline const &other, Real horizon_padding) const
681 return internal_distance (other, horizon_padding, &dummy);
685 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
688 internal_distance (other, horizon_padding, &touch);
693 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
695 if (horizon_padding == 0.0)
696 return internal_distance (other, touch_point);
698 // Note that it is not necessary to build a padded version of other,
699 // because the same effect can be achieved just by doubling horizon_padding.
700 Skyline padded_this = padded (horizon_padding);
701 return padded_this.internal_distance (other, touch_point);
705 Skyline::padded (Real horizon_padding) const
707 list<Building> pad_buildings;
708 for (list<Building>::const_iterator i = buildings_.begin (); i != buildings_.end (); ++i)
710 if (i->start_ > -infinity_f)
712 Real height = i->height (i->start_);
713 if (height > -infinity_f)
715 // Add the sloped building that pads the left side of the current building.
716 Real start = i->start_ - 2 * horizon_padding;
717 Real end = i->start_ - horizon_padding;
718 pad_buildings.push_back (Building (start, height - horizon_padding, height, end));
720 // Add the flat building that pads the left side of the current building.
721 start = i->start_ - horizon_padding;
723 pad_buildings.push_back (Building (start, height, height, end));
727 if (i->end_ < infinity_f)
729 Real height = i->height (i->end_);
730 if (height > -infinity_f)
732 // Add the flat building that pads the right side of the current building.
733 Real start = i->end_;
734 Real end = start + horizon_padding;
735 pad_buildings.push_back (Building (start, height, height, end));
737 // Add the sloped building that pads the right side of the current building.
739 end += horizon_padding;
740 pad_buildings.push_back (Building (start, height, height - horizon_padding, end));
745 // The buildings may be overlapping, so resolve that.
746 list<Building> pad_skyline = internal_build_skyline (&pad_buildings);
748 // Merge the padding with the original, to make a new skyline.
749 Skyline padded (sky_);
750 list<Building> my_buildings = buildings_;
751 padded.buildings_.clear ();
752 internal_merge_skyline (&pad_skyline, &my_buildings, &padded.buildings_);
759 Skyline::internal_distance (Skyline const &other, Real *touch_point) const
761 assert (sky_ == -other.sky_);
763 list<Building>::const_iterator i = buildings_.begin ();
764 list<Building>::const_iterator j = other.buildings_.begin ();
766 Real dist = -infinity_f;
767 Real start = -infinity_f;
768 Real touch = -infinity_f;
769 while (i != buildings_.end () && j != other.buildings_.end ())
771 Real end = min (i->end_, j->end_);
772 Real start_dist = i->height (start) + j->height (start);
773 Real end_dist = i->height (end) + j->height (end);
774 dist = max (dist, max (start_dist, end_dist));
776 if (end_dist == dist)
778 else if (start_dist == dist)
781 if (i->end_ <= j->end_)
788 *touch_point = touch;
793 Skyline::height (Real airplane) const
795 assert (!isinf (airplane));
797 list<Building>::const_iterator i;
798 for (i = buildings_.begin (); i != buildings_.end (); i++)
800 if (i->end_ >= airplane)
801 return sky_ * i->height (airplane);
809 Skyline::max_height () const
811 Real ret = -infinity_f;
813 list<Building>::const_iterator i;
814 for (i = buildings_.begin (); i != buildings_.end (); ++i)
816 ret = max (ret, i->height (i->start_));
817 ret = max (ret, i->height (i->end_));
824 Skyline::max_height_position () const
827 s.set_minimum_height (0);
828 return touching_point (s);
832 Skyline::set_minimum_height (Real h)
835 s.buildings_.front ().y_intercept_ = h * sky_;
840 Skyline::to_points (Axis horizon_axis) const
844 Real start = -infinity_f;
845 for (list<Building>::const_iterator i (buildings_.begin ());
846 i != buildings_.end (); i++)
848 out.push_back (Offset (start, sky_ * i->height (start)));
849 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
853 if (horizon_axis == Y_AXIS)
854 for (vsize i = 0; i < out.size (); i++)
855 out[i] = out[i].swapped ();
861 Skyline::is_empty () const
863 if (!buildings_.size ())
865 Building b = buildings_.front ();
866 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
873 empty_skyline (&buildings_);
876 /****************************************************************/
878 IMPLEMENT_SIMPLE_SMOBS (Skyline);
879 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
880 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
883 Skyline::mark_smob (SCM s)
885 ASSERT_LIVE_IS_ALLOWED (s);
890 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
892 Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
895 scm_puts ("#<Skyline>", port);
900 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
902 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
904 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
906 Real horizon_padding = 0;
907 if (horizon_padding_scm != SCM_UNDEFINED)
909 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
910 horizon_padding = scm_to_double (horizon_padding_scm);
913 Skyline *skyline = Skyline::unsmob (skyline_scm);
914 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
915 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
918 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
920 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
922 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
924 Real horizon_padding = 0;
925 if (horizon_padding_scm != SCM_UNDEFINED)
927 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
928 horizon_padding = scm_to_double (horizon_padding_scm);
931 Skyline *skyline = Skyline::unsmob (skyline_scm);
932 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
933 return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
936 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
938 Skyline::get_max_height (SCM skyline_scm)
940 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
943 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
945 Skyline::get_max_height_position (SCM skyline_scm)
947 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
950 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
952 Skyline::get_height (SCM skyline_scm, SCM x_scm)
954 Real x = robust_scm2double (x_scm, 0.0);
955 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));