]> git.donarmstrong.com Git - lilypond.git/blob - lily/skyline.cc
0250fc07f4743ba44e7bc1ac8af23ed6057aa435
[lilypond.git] / lily / skyline.cc
1 /*
2   This file is part of LilyPond, the GNU music typesetter.
3
4   Copyright (C) 2006--2012 Joe Neeman <joeneeman@gmail.com>
5
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.
10
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.
15
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/>.
18 */
19
20 #include "skyline.hh"
21 #include <deque>
22 #include <cstdio>
23
24 #include "ly-smobs.icc"
25
26 /* A skyline is a sequence of non-overlapping buildings: something like
27    this:
28                    _______
29                   |       \                                 ________
30                   |        \                       ________/        \
31         /\        |          \                    /                  \
32        /  --------             \                 /                    \
33       /                          \              /                      \
34      /                             ------------/                        ----
35    --
36    Each building has a starting position, and ending position, a starting
37    height and an ending height.
38
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
45       building
46
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.
51
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.
56 */
57
58 /* If we start including very thin buildings, numerical accuracy errors can
59    arise. Therefore, we ignore all buildings that are less than epsilon wide. */
60 #define EPS 1e-5
61
62 static void
63 print_buildings (list<Building> const &b)
64 {
65   for (list<Building>::const_iterator i = b.begin (); i != b.end (); i++)
66     i->print ();
67 }
68
69 void
70 Skyline::print () const
71 {
72   print_buildings (buildings_);
73 }
74
75 void
76 Skyline::print_points () const
77 {
78   vector<Offset> ps (to_points (X_AXIS));
79
80   for (vsize i = 0; i < ps.size (); i++)
81     printf ("(%f,%f)%s", ps[i][X_AXIS], ps[i][Y_AXIS],
82             (i % 2) == 1 ? "\n" : " ");
83 }
84
85 Building::Building (Real start, Real start_height, Real end_height, Real end)
86 {
87   if (isinf (start) || isinf (end))
88     assert (start_height == end_height);
89
90   end_ = end;
91   precompute (start, start_height, end_height, end);
92 }
93
94 Building::Building (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
95 {
96   Real start = b[horizon_axis][LEFT] - horizon_padding;
97   Real end = b[horizon_axis][RIGHT] + horizon_padding;
98   Real height = sky * b[other_axis (horizon_axis)][sky];
99
100   end_ = end;
101   precompute (start, height, height, end);
102 }
103
104 void
105 Building::precompute (Real start, Real start_height, Real end_height, Real end)
106 {
107   slope_ = (end_height - start_height) / (end - start);
108   if (start_height == end_height) /* if they were both infinite, we would get nan, not 0, from the prev line */
109     slope_ = 0;
110
111   assert (!isinf (slope_) && !isnan (slope_));
112
113   if (isinf (start))
114     {
115       assert (start_height == end_height);
116       y_intercept_ = start_height;
117     }
118   else
119     y_intercept_ = start_height - slope_ * start;
120 }
121
122 Real
123 Building::height (Real x) const
124 {
125   return isinf (x) ? y_intercept_ : slope_ * x + y_intercept_;
126 }
127
128 void
129 Building::print () const
130 {
131   printf ("%f x + %f ends at %f\n", slope_, y_intercept_, end_);
132 }
133
134 Real
135 Building::intersection_x (Building const &other) const
136 {
137   Real ret = (y_intercept_ - other.y_intercept_) / (other.slope_ - slope_);
138   return isnan (ret) ? -infinity_f : ret;
139 }
140
141 void
142 Building::leading_part (Real chop)
143 {
144   assert (chop <= end_);
145   end_ = chop;
146 }
147
148 Building
149 Building::sloped_neighbour (Real start, Real horizon_padding, Direction d) const
150 {
151   Real x = (d == LEFT) ? start : end_;
152   Real left = x;
153   Real right = x + d * horizon_padding;
154   Real left_height = height (x);
155   Real right_height = left_height - horizon_padding;
156   if (d == LEFT)
157     {
158       swap (left, right);
159       swap (left_height, right_height);
160     }
161   return Building (left, left_height, right_height, right);
162 }
163
164 static Real
165 first_intersection (Building const &b, list<Building> *const s, Real start_x)
166 {
167   while (!s->empty () && start_x < b.end_)
168     {
169       Building c = s->front ();
170       if (c.conceals (b, start_x))
171         return start_x;
172
173       Real i = b.intersection_x (c);
174       if (i > start_x && i <= b.end_ && i <= c.end_)
175         return i;
176
177       start_x = c.end_;
178       if (b.end_ > c.end_)
179         s->pop_front ();
180     }
181   return b.end_;
182 }
183
184 bool
185 Building::conceals (Building const &other, Real x) const
186 {
187   if (slope_ == other.slope_)
188     return y_intercept_ > other.y_intercept_;
189
190   /* their slopes were not equal, so there is an intersection point */
191   Real i = intersection_x (other);
192   return (i <= x && slope_ > other.slope_)
193          || (i > x && slope_ < other.slope_);
194 }
195
196 void
197 Skyline::internal_merge_skyline (list<Building> *s1, list<Building> *s2,
198                                  list<Building> *const result)
199 {
200   if (s1->empty () || s2->empty ())
201     {
202       programming_error ("tried to merge an empty skyline");
203       return;
204     }
205
206   Real x = -infinity_f;
207   while (!s1->empty ())
208     {
209       if (s2->front ().conceals (s1->front (), x))
210         swap (s1, s2);
211
212       Building b = s1->front ();
213       Real end = first_intersection (b, s2, x);
214
215       if (s2->empty ())
216         {
217           result->push_front (b);
218           break;
219         }
220
221       /* only include buildings wider than epsilon */
222       if (end > x + EPS)
223         {
224           b.leading_part (end);
225           result->push_front (b);
226         }
227
228       if (end >= s1->front ().end_)
229         s1->pop_front ();
230
231       x = end;
232     }
233   result->reverse ();
234 }
235
236 static void
237 empty_skyline (list<Building> *const ret)
238 {
239   ret->push_front (Building (-infinity_f, -infinity_f, -infinity_f, infinity_f));
240 }
241
242 /*
243   Given Building 'b' with starting wall location 'start', extend each side
244   with a sloped roofline of width 'horizon_padding'; put the skyline in 'ret'
245 */
246 static void
247 single_skyline (Building b, Real start, Real horizon_padding, list<Building> *const ret)
248 {
249   bool sloped_neighbours = horizon_padding > 0 && !isinf (start) && !isinf (b.end_);
250   if (!isinf (b.end_))
251     ret->push_front (Building (b.end_ + horizon_padding, -infinity_f,
252                                -infinity_f, infinity_f));
253   if (sloped_neighbours)
254     ret->push_front (b.sloped_neighbour (start, horizon_padding, RIGHT));
255
256   if (b.end_ > start + EPS)
257     ret->push_front (b);
258
259   if (sloped_neighbours)
260     ret->push_front (b.sloped_neighbour (start, horizon_padding, LEFT));
261
262   if (!isinf (start))
263     ret->push_front (Building (-infinity_f, -infinity_f,
264                                -infinity_f, start - horizon_padding));
265 }
266
267 /* remove a non-overlapping set of boxes from BOXES and build a skyline
268    out of them */
269 static list<Building>
270 non_overlapping_skyline (list<Box> *const boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
271 {
272   list<Building> result;
273   Real last_end = -infinity_f;
274   list<Box>::iterator i = boxes->begin ();
275   while (i != boxes->end ())
276     {
277       Interval iv = (*i)[horizon_axis];
278
279       if (iv[LEFT] - horizon_padding < last_end)
280         {
281           i++;
282           continue;
283         }
284
285       if (iv[LEFT] - horizon_padding > last_end + EPS)
286         result.push_front (Building (last_end, -infinity_f, -infinity_f, iv[LEFT] - 2 * horizon_padding));
287
288       Building b (*i, horizon_padding, horizon_axis, sky);
289       bool sloped_neighbours = horizon_padding > 0 && !isinf (iv.length ());
290       if (sloped_neighbours)
291         result.push_front (b.sloped_neighbour (iv[LEFT] - horizon_padding, horizon_padding, LEFT));
292       result.push_front (b);
293       if (sloped_neighbours)
294         result.push_front (b.sloped_neighbour (iv[LEFT] - horizon_padding, horizon_padding, RIGHT));
295
296       list<Box>::iterator j = i++;
297       boxes->erase (j);
298       last_end = result.front ().end_;
299     }
300   if (last_end < infinity_f)
301     result.push_front (Building (last_end, -infinity_f, -infinity_f, infinity_f));
302   result.reverse ();
303   return result;
304 }
305
306 class LessThanBox
307 {
308   Axis a_;
309
310 public:
311   LessThanBox (Axis a)
312   {
313     a_ = a;
314   }
315
316   bool operator () (Box const &b1, Box const &b2)
317   {
318     return b1[a_][LEFT] < b2[a_][LEFT];
319   }
320 };
321
322 list<Building>
323 Skyline::internal_build_skyline (list<Box> *boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
324 {
325   vsize size = boxes->size ();
326
327   if (size == 0)
328     {
329       list<Building> result;
330       empty_skyline (&result);
331       return result;
332     }
333   else if (size == 1)
334     {
335       list<Building> result;
336       single_skyline (Building (boxes->front (), horizon_padding, horizon_axis, sky),
337                       boxes->front ()[horizon_axis][LEFT] - horizon_padding,
338                       horizon_padding, &result);
339       return result;
340     }
341
342   deque<list<Building> > partials;
343   boxes->sort (LessThanBox (horizon_axis));
344   while (!boxes->empty ())
345     partials.push_back (non_overlapping_skyline (boxes, horizon_padding, horizon_axis, sky));
346
347   /* we'd like to say while (partials->size () > 1) but that's O (n).
348      Instead, we exit in the middle of the loop */
349   while (!partials.empty ())
350     {
351       list<Building> merged;
352       list<Building> one = partials.front ();
353       partials.pop_front ();
354       if (partials.empty ())
355         return one;
356
357       list<Building> two = partials.front ();
358       partials.pop_front ();
359       internal_merge_skyline (&one, &two, &merged);
360       partials.push_back (merged);
361     }
362   assert (0);
363   return list<Building> ();
364 }
365
366 Skyline::Skyline ()
367 {
368   sky_ = UP;
369   empty_skyline (&buildings_);
370 }
371
372 Skyline::Skyline (Skyline const &src)
373 {
374   sky_ = src.sky_;
375
376   /* doesn't a list's copy constructor do this? -- jneem */
377   for (list<Building>::const_iterator i = src.buildings_.begin ();
378        i != src.buildings_.end (); i++)
379     {
380       buildings_.push_back (Building ((*i)));
381     }
382 }
383
384 Skyline::Skyline (Direction sky)
385 {
386   sky_ = sky;
387   empty_skyline (&buildings_);
388 }
389
390 /*
391   build padded skyline from an existing skyline with padding
392   added to it.
393 */
394
395 Skyline::Skyline (Skyline const &src, Real horizon_padding, Axis /*a*/)
396 {
397   /*
398      We extract boxes from the skyline, then build a new skyline from
399      the boxes.
400      A box is created for every horizontal portion of the skyline
401      Because skylines are defined positive, and then inverted if they
402      are to be down-facing, we create the new skyline in the UP
403      direction, then give it the down direction if needed.
404   */
405   Real start = -infinity_f;
406   list<Box> boxes;
407
408   // establish a baseline box
409   // FIXME: This has hardcoded logic, assuming a == X_AXIS!
410   boxes.push_back (Box (Interval (-infinity_f, infinity_f),
411                         Interval (0, 0)));
412   list<Building>::const_iterator end = src.buildings_.end ();
413   for (list<Building>::const_iterator i = src.buildings_.begin (); i != end; start = i->end_, i++)
414     if ((i->slope_ == 0) && !isinf (i->y_intercept_))
415       boxes.push_back (Box (Interval (start, i->end_),
416                             Interval (-infinity_f, i->y_intercept_)));
417   buildings_ = internal_build_skyline (&boxes, horizon_padding, X_AXIS, UP);
418   sky_ = src.sky_;
419 }
420
421 /*
422   build skyline from a set of boxes. If horizon_padding > 0, expand all the boxes
423   by that amount and add 45-degree sloped boxes to the edges of each box (of
424   width horizon_padding). That is, the total amount of horizontal expansion is
425   horizon_padding*4, half of which is sloped and half of which is flat.
426
427   Boxes should have fatness in the horizon_axis (after they are expanded by
428   horizon_padding), otherwise they are ignored.
429  */
430 Skyline::Skyline (vector<Box> const &boxes, Real horizon_padding, Axis horizon_axis, Direction sky)
431 {
432   list<Box> filtered_boxes;
433   sky_ = sky;
434
435   Axis vert_axis = other_axis (horizon_axis);
436   for (vsize i = 0; i < boxes.size (); i++)
437     {
438       Interval iv = boxes[i][horizon_axis];
439       iv.widen (horizon_padding);
440       if (iv.length () > EPS && !boxes[i][vert_axis].is_empty ())
441         filtered_boxes.push_front (boxes[i]);
442     }
443
444   buildings_ = internal_build_skyline (&filtered_boxes, horizon_padding, horizon_axis, sky);
445 }
446
447 Skyline::Skyline (Box const &b, Real horizon_padding, Axis horizon_axis, Direction sky)
448 {
449   sky_ = sky;
450   Building front (b, horizon_padding, horizon_axis, sky);
451   single_skyline (front, b[horizon_axis][LEFT] - horizon_padding,
452                   horizon_padding, &buildings_);
453 }
454
455 void
456 Skyline::merge (Skyline const &other)
457 {
458   assert (sky_ == other.sky_);
459
460   list<Building> other_bld (other.buildings_);
461   list<Building> my_bld;
462   my_bld.splice (my_bld.begin (), buildings_);
463   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
464 }
465
466 void
467 Skyline::insert (Box const &b, Real horizon_padding, Axis a)
468 {
469   list<Building> other_bld;
470   list<Building> my_bld;
471
472   if (isnan (b[other_axis (a)][LEFT])
473       || isnan (b[other_axis (a)][RIGHT]))
474     {
475       programming_error ("insane box for skyline");
476       return;
477     }
478
479   /* do the same filtering as in Skyline (vector<Box> const&, etc.) */
480   Interval iv = b[a];
481   iv.widen (horizon_padding);
482   if (iv.length () <= EPS || b[other_axis (a)].is_empty ())
483     return;
484
485   my_bld.splice (my_bld.begin (), buildings_);
486   single_skyline (Building (b, horizon_padding, a, sky_), b[a][LEFT] - horizon_padding,
487                   horizon_padding, &other_bld);
488   internal_merge_skyline (&other_bld, &my_bld, &buildings_);
489 }
490
491 void
492 Skyline::raise (Real r)
493 {
494   list<Building>::iterator end = buildings_.end ();
495   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
496     i->y_intercept_ += sky_ * r;
497 }
498
499 void
500 Skyline::shift (Real s)
501 {
502   list<Building>::iterator end = buildings_.end ();
503   for (list<Building>::iterator i = buildings_.begin (); i != end; i++)
504     {
505       i->end_ += s;
506       i->y_intercept_ -= s * i->slope_;
507     }
508 }
509
510 Real
511 Skyline::distance (Skyline const &other, Real horizon_padding) const
512 {
513   Real dummy;
514   return internal_distance (other, horizon_padding, &dummy);
515 }
516
517 Real
518 Skyline::touching_point (Skyline const &other, Real horizon_padding) const
519 {
520   Real touch;
521   internal_distance (other, horizon_padding, &touch);
522   return touch;
523 }
524
525 Real
526 Skyline::internal_distance (Skyline const &other, Real horizon_padding, Real *touch_point) const
527 {
528   assert (sky_ == -other.sky_);
529
530   Skyline const *padded_this = this;
531   Skyline const *padded_other = &other;
532   bool created_tmp_skylines = false;
533
534   /*
535     For systems, padding is not added at creation time.  Padding is
536     added to AxisGroup objects when outside-staff objects are added.
537     Thus, when we want to place systems with horizontal padding,
538     we do it at distance calculation time.
539   */
540   if (horizon_padding != 0.0)
541     {
542       padded_this = new Skyline (*padded_this, horizon_padding, X_AXIS);
543       padded_other = new Skyline (*padded_other, horizon_padding, X_AXIS);
544       created_tmp_skylines = true;
545     }
546
547   list<Building>::const_iterator i = padded_this->buildings_.begin ();
548   list<Building>::const_iterator j = padded_other->buildings_.begin ();
549
550   Real dist = -infinity_f;
551   Real start = -infinity_f;
552   Real touch = -infinity_f;
553   while (i != padded_this->buildings_.end () && j != padded_other->buildings_.end ())
554     {
555       Real end = min (i->end_, j->end_);
556       Real start_dist = i->height (start) + j->height (start);
557       Real end_dist = i->height (end) + j->height (end);
558       dist = max (dist, max (start_dist, end_dist));
559
560       if (end_dist == dist)
561         touch = end;
562       else if (start_dist == dist)
563         touch = start;
564
565       if (i->end_ <= j->end_)
566         i++;
567       else
568         j++;
569       start = end;
570     }
571
572   if (created_tmp_skylines)
573     {
574       delete padded_this;
575       delete padded_other;
576     }
577
578   *touch_point = touch;
579   return dist;
580 }
581
582 Real
583 Skyline::height (Real airplane) const
584 {
585   assert (!isinf (airplane));
586
587   list<Building>::const_iterator i;
588   for (i = buildings_.begin (); i != buildings_.end (); i++)
589     {
590       if (i->end_ >= airplane)
591         return sky_ * i->height (airplane);
592     }
593
594   assert (0);
595   return 0;
596 }
597
598 Real
599 Skyline::max_height () const
600 {
601   Skyline s (-sky_);
602   s.set_minimum_height (0);
603   return sky_ * distance (s);
604 }
605
606 Real
607 Skyline::max_height_position () const
608 {
609   Skyline s (-sky_);
610   s.set_minimum_height (0);
611   return touching_point (s);
612 }
613
614 void
615 Skyline::set_minimum_height (Real h)
616 {
617   Skyline s (sky_);
618   s.buildings_.front ().y_intercept_ = h * sky_;
619   merge (s);
620 }
621
622 vector<Offset>
623 Skyline::to_points (Axis horizon_axis) const
624 {
625   vector<Offset> out;
626
627   Real start = -infinity_f;
628   for (list<Building>::const_iterator i (buildings_.begin ());
629        i != buildings_.end (); i++)
630     {
631       out.push_back (Offset (start, sky_ * i->height (start)));
632       out.push_back (Offset (i->end_, sky_ * i->height (i->end_)));
633       start = i->end_;
634     }
635
636   if (horizon_axis == Y_AXIS)
637     for (vsize i = 0; i < out.size (); i++)
638       out[i] = out[i].swapped ();
639
640   return out;
641 }
642
643 bool
644 Skyline::is_empty () const
645 {
646   Building b = buildings_.front ();
647   return b.end_ == infinity_f && b.y_intercept_ == -infinity_f;
648 }
649
650 void
651 Skyline::clear ()
652 {
653   buildings_.clear ();
654   empty_skyline (&buildings_);
655 }
656
657 /****************************************************************/
658
659 IMPLEMENT_SIMPLE_SMOBS (Skyline);
660 IMPLEMENT_TYPE_P (Skyline, "ly:skyline?");
661 IMPLEMENT_DEFAULT_EQUAL_P (Skyline);
662
663 SCM
664 Skyline::mark_smob (SCM s)
665 {
666   ASSERT_LIVE_IS_ALLOWED (s);
667   return SCM_EOL;
668 }
669
670 int
671 Skyline::print_smob (SCM s, SCM port, scm_print_state *)
672 {
673   Skyline *r = (Skyline *) SCM_CELL_WORD_1 (s);
674   (void) r;
675
676   scm_puts ("#<Skyline>", port);
677
678   return 1;
679 }
680
681 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_touching_point, 3, 1, "")
682 SCM
683 Skyline::get_touching_point (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
684 {
685   LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
686
687   Real horizon_padding = 0;
688   if (horizon_padding_scm != SCM_UNDEFINED)
689     {
690       LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
691       horizon_padding = scm_to_double (horizon_padding_scm);
692     }
693
694   Skyline *skyline = Skyline::unsmob (skyline_scm);
695   Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
696   return scm_from_double (skyline->touching_point (*other_skyline, horizon_padding));
697 }
698
699 MAKE_SCHEME_CALLBACK_WITH_OPTARGS (Skyline, get_distance, 3, 1, "")
700 SCM
701 Skyline::get_distance (SCM skyline_scm, SCM other_skyline_scm, SCM horizon_padding_scm)
702 {
703   LY_ASSERT_SMOB (Skyline, other_skyline_scm, 1);
704
705   Real horizon_padding = 0;
706   if (horizon_padding_scm != SCM_UNDEFINED)
707     {
708       LY_ASSERT_TYPE (scm_is_number, horizon_padding_scm, 3);
709       horizon_padding = scm_to_double (horizon_padding_scm);
710     }
711
712   Skyline *skyline = Skyline::unsmob (skyline_scm);
713   Skyline *other_skyline = Skyline::unsmob (other_skyline_scm);
714   return scm_from_double (skyline->distance (*other_skyline, horizon_padding));
715 }
716
717 MAKE_SCHEME_CALLBACK (Skyline, get_max_height, 1)
718 SCM
719 Skyline::get_max_height (SCM skyline_scm)
720 {
721   return scm_from_double (Skyline::unsmob (skyline_scm)->max_height ());
722 }
723
724 MAKE_SCHEME_CALLBACK (Skyline, get_max_height_position, 1)
725 SCM
726 Skyline::get_max_height_position (SCM skyline_scm)
727 {
728   return scm_from_double (Skyline::unsmob (skyline_scm)->max_height_position ());
729 }
730
731 MAKE_SCHEME_CALLBACK (Skyline, get_height, 2)
732 SCM
733 Skyline::get_height (SCM skyline_scm, SCM x_scm)
734 {
735   Real x = robust_scm2double (x_scm, 0.0);
736   return scm_from_double (Skyline::unsmob (skyline_scm)->height (x));
737 }