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;
216 for (i = buildings_.begin (); i != buildings_.end (); i++)
218 if (last_empty && i->y_intercept_ == -infinity_f)
220 list<Building>::iterator last = i;
222 last->end_ = i->end_;
223 buildings_.erase (i);
226 last_empty = (i->y_intercept_ == -infinity_f);
229 assert (buildings_.front ().start_ == -infinity_f);
230 assert (buildings_.back ().end_ == infinity_f);
236 // Since a skyline should always be normalized, we can
237 // assume that there are never two adjacent empty buildings.
238 // That is, if center is empty then left and right are not.
239 list<Building>::iterator left = buildings_.begin ();
240 list<Building>::iterator center = buildings_.begin ();
241 list<Building>::iterator right;
243 for (right = buildings_.begin (); right != buildings_.end (); right++)
245 if (center != buildings_.begin () && center->y_intercept_ == -infinity_f)
247 Real p1 = left->height (left->end_);
248 Real p2 = right->height (right->start_);
249 *center = Building (center->start_, p1, p2, center->end_);
258 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
259 list<Building> *const result) const
261 if (s1->empty () || s2->empty ())
263 programming_error ("tried to merge an empty skyline");
267 Real x = -infinity_f;
268 Real last_end = -infinity_f;
269 while (!s1->empty ())
271 if (s2->front ().conceals (s1->front (), x))
274 Building b = s1->front ();
275 Building c = s2->front ();
277 // Optimization: if the other skyline is empty at this point,
278 // we can avoid testing some intersections. Just grab as many
279 // buildings from s1 as we can, and shove them onto the output.
280 if (c.y_intercept_ == -infinity_f
283 list<Building>::iterator i = s1->begin ();
285 while (i != s1->end () && i->end_ <= c.end_)
288 s1->front ().start_ = x;
289 result->splice (result->end (), *s1, s1->begin (), i);
290 x = result->back ().end_;
295 Real end = first_intersection (b, s2, x);
299 result->push_back (b);
303 /* only include buildings wider than epsilon */
306 b.leading_part (end);
309 result->push_back (b);
312 if (end >= s1->front ().end_)
320 empty_skyline (list<Building> *const ret)
322 ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
326 Given Building 'b', build a skyline containing only that building.
329 single_skyline (Building b, list<Building> *const ret)
331 if (b.end_ > b.start_ + EPS)
333 ret->push_back (Building (-infinity_f, -infinity_f,
334 -infinity_f, b.start_));
336 ret->push_back (Building (b.end_, -infinity_f,
337 -infinity_f, infinity_f));
345 /* remove a non-overlapping set of boxes from BOXES and build a skyline
347 static list<Building>
348 non_overlapping_skyline (list<Building> *const buildings)
350 list<Building> result;
351 Real last_end = -infinity_f;
352 Building last_building (-infinity_f, -infinity_f, -infinity_f, infinity_f);
353 list<Building>::iterator i = buildings->begin ();
354 while (i != buildings->end ())
357 Real y1 = i->height (i->start_);
359 Real y2 = i->height (i->end_);
361 // Drop buildings that will obviously have no effect.
362 if (last_building.height (x1) >= y1
363 && last_building.end_ >= x2
364 && last_building.height (x2) >= y2)
366 list<Building>::iterator j = i++;
367 buildings->erase (j);
377 if (x1 > last_end + EPS)
378 result.push_back (Building (last_end, -infinity_f, -infinity_f, x1));
380 result.push_back (*i);
384 list<Building>::iterator j = i++;
385 buildings->erase (j);
388 if (last_end < infinity_f)
389 result.push_back (Building (last_end, -infinity_f, -infinity_f, infinity_f));
393 class LessThanBuilding
396 bool operator () (Building const &b1, Building const &b2)
398 return b1.start_ < b2.start_
399 || (b1.start_ == b2.start_ && b1.height (b1.start_) > b2.height (b1.start_));
404 BUILDINGS is a list of buildings, but they could be overlapping
405 and in any order. The returned list of buildings is ordered and non-overlapping.
408 Skyline::internal_build_skyline (list<Building> *buildings) const
410 vsize size = buildings->size ();
414 list<Building> result;
415 empty_skyline (&result);
420 list<Building> result;
421 single_skyline (buildings->front (), &result);
425 deque<list<Building> > partials;
426 buildings->sort (LessThanBuilding ());
427 while (!buildings->empty ())
428 partials.push_back (non_overlapping_skyline (buildings));
430 /* we'd like to say while (partials->size () > 1) but that's O (n).
431 Instead, we exit in the middle of the loop */
432 while (!partials.empty ())
434 list<Building> merged;
435 list<Building> one = partials.front ();
436 partials.pop_front ();
437 if (partials.empty ())
440 list<Building> two = partials.front ();
441 partials.pop_front ();
442 internal_merge_skyline (&one, &two, &merged);
443 partials.push_back (merged);
446 return list<Building> ();
452 empty_skyline (&buildings_);
455 Skyline::Skyline (Skyline const &src)
459 /* doesn't a list's copy constructor do this? -- jneem */
460 for (list<Building>::const_iterator i = src.buildings_.begin ();
461 i != src.buildings_.end (); i++)
463 buildings_.push_back (Building ((*i)));
467 Skyline::Skyline (Direction sky)
470 empty_skyline (&buildings_);
474 Build skyline from a set of boxes.
476 Boxes should have fatness in the horizon_axis, otherwise they are ignored.
478 Skyline::Skyline (vector<Box> const &boxes, Axis horizon_axis, Direction sky)
480 list<Building> buildings;
483 Axis vert_axis = other_axis (horizon_axis);
484 for (vsize i = 0; i < boxes.size (); i++)
486 Interval iv = boxes[i][horizon_axis];
487 if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
488 buildings.push_front (Building (boxes[i], horizon_axis, sky));
491 buildings_ = internal_build_skyline (&buildings);
496 build skyline from a set of line segments.
498 Buildings should have fatness in the horizon_axis, otherwise they are ignored.
500 Skyline::Skyline (vector<Drul_array<Offset> > const &segments, Axis horizon_axis, Direction sky)
502 list<Building> buildings;
505 for (vsize i = 0; i < segments.size (); i++)
507 Drul_array<Offset> const &seg = segments[i];
508 Offset left = seg[LEFT];
509 Offset right = seg[RIGHT];
510 if (left[horizon_axis] > right[horizon_axis])
513 Real x1 = left[horizon_axis];
514 Real x2 = right[horizon_axis];
515 Real y1 = left[other_axis (horizon_axis)] * sky;
516 Real y2 = right[other_axis (horizon_axis)] * sky;
519 buildings.push_back (Building (x1, y1, y2, x2));
522 buildings_ = internal_build_skyline (&buildings);
526 Skyline::Skyline (vector<Skyline_pair> const &skypairs, Direction sky)
530 deque<Skyline> partials;
531 for (vsize i = 0; i < skypairs.size (); i++)
532 partials.push_back (Skyline ((skypairs[i])[sky]));
534 while (partials.size () > 1)
536 Skyline one = partials.front ();
537 partials.pop_front ();
538 Skyline two = partials.front ();
539 partials.pop_front ();
542 partials.push_back (one);
545 if (partials.size ())
546 buildings_.swap (partials.front ().buildings_);
551 Skyline::Skyline (Box const &b, Axis horizon_axis, Direction sky)
554 Building front (b, horizon_axis, sky);
555 single_skyline (front, &buildings_);
559 Skyline::merge (Skyline const &other)
561 assert (sky_ == other.sky_);
563 if (other.is_empty ())
568 buildings_ = other.buildings_;
572 list<Building> other_bld (other.buildings_);
573 list<Building> my_bld;
574 my_bld.splice (my_bld.begin (), buildings_);
575 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
580 Skyline::insert (Box const &b, Axis a)
582 list<Building> other_bld;
583 list<Building> my_bld;
585 if (isnan (b[other_axis (a)][LEFT])
586 || isnan (b[other_axis (a)][RIGHT]))
588 programming_error ("insane box for skyline");
592 /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
594 if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
597 my_bld.splice (my_bld.begin (), buildings_);
598 single_skyline (Building (b, a, sky_), &other_bld);
599 internal_merge_skyline (&other_bld, &my_bld, &buildings_);
604 Skyline::raise (Real r)
606 list<Building>::iterator end = buildings_.end ();
607 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
608 i->y_intercept_ += sky_ * r;
612 Skyline::shift (Real s)
614 list<Building>::iterator end = buildings_.end ();
615 for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
619 i->y_intercept_ -= s * i->slope_;
624 Skyline::distance (Skyline const &other, Real horizon_padding) const
627 return internal_distance (other, horizon_padding, &dummy);
631 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
634 internal_distance (other, horizon_padding, &touch);
639 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
641 if (horizon_padding == 0.0)
642 return internal_distance (other, touch_point);
644 // Note that it is not necessary to build a padded version of other,
645 // because the same effect can be achieved just by doubling horizon_padding.
646 Skyline padded_this = padded (horizon_padding);
647 return padded_this.internal_distance (other, touch_point);
651 Skyline::padded (Real horizon_padding) const
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::left () const
772 for (list<Building>::const_iterator i (buildings_.begin ());
773 i != buildings_.end (); i++)
774 if (i->y_intercept_ > -infinity_f)
781 Skyline::right () const
783 for (list<Building>::const_reverse_iterator i (buildings_.rbegin ());
784 i != buildings_.rend (); ++i)
785 if (i->y_intercept_ > -infinity_f)
792 Skyline::max_height_position () const
795 s.set_minimum_height (0);
796 return touching_point (s);
800 Skyline::set_minimum_height (Real h)
803 s.buildings_.front ().y_intercept_ = h * sky_;
808 Skyline::to_points (Axis horizon_axis) const
812 Real start = -infinity_f;
813 for (list<Building>::const_iterator i (buildings_.begin ());
814 i != buildings_.end (); i++)
816 out.push_back (Offset (start, sky_ * i->height (start)));
817 out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
821 if (horizon_axis == Y_AXIS)
822 for (vsize i = 0; i < out.size (); i++)
823 out[i] = out[i].swapped ();
829 Skyline::is_empty () const
831 if (!buildings_.size ())
833 Building b = buildings_.front ();
834 return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
841 empty_skyline (&buildings_);
844 /****************************************************************/
846 IMPLEMENT_SIMPLE_SMOBS (Skyline);
847 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
848 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
851 Skyline::mark_smob (SCM s)
853 ASSERT_LIVE_IS_ALLOWED (s);
858 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
860 Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
863 scm_puts ("#<Skyline>", port);
868 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
870 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
872 LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
874 Real horizon_padding = 0;
875 if (horizon_padding_scm != SCM_UNDEFINED)
877 LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
878 horizon_padding = scm_to_double (horizon_padding_scm);
881 Skyline *skyline = Skyline::unsmob (skyline_scm);
882 Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
883 return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
886 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
888 Skyline::get_distance (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->distance (*other_skyline, horizon_padding));
904 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
906 Skyline::get_max_height (SCM skyline_scm)
908 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
911 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
913 Skyline::get_max_height_position (SCM skyline_scm)
915 return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
918 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
920 Skyline::get_height (SCM skyline_scm, SCM x_scm)
922 Real x = robust_scm2double (x_scm, 0.0);
923 return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));