lily/beam-quanting.cc

   1 /*
   2   beam-quanting.cc -- implement Beam quanting functions
   3
   4   source file of the GNU LilyPond music typesetter
   5
   6   (c) 1997--2006 Han-Wen Nienhuys <hanwen@xs4all.nl>
   7   Jan Nieuwenhuizen <janneke@gnu.org>
   8 */
   9
  10 #include "beam.hh"
  11
  12 #include <algorithm>
  13 using namespace std;
  14
  15 #include "grob.hh"
  16 #include "align-interface.hh"
  17 #include "international.hh"
  18 #include "output-def.hh"
  19 #include "pointer-group-interface.hh"
  20 #include "staff-symbol-referencer.hh"
  21 #include "stem.hh"
  22 #include "warn.hh"
  23 #include "main.hh"
  24
  25 Real
  26 get_detail (SCM alist, SCM sym, Real def)
  27 {
  28   SCM entry = scm_assq (sym, alist);
  29
  30   if (scm_is_pair (entry))
  31     return robust_scm2double (scm_cdr (entry), def);
  32   return def;
  33 }
  34
  35 void
  36 Beam_quant_parameters::fill (Grob *him)
  37 {
  38   SCM details = him->get_property ("details");
  39
  40   INTER_QUANT_PENALTY = get_detail (details, ly_symbol2scm ("inter-quant-penalty"), 1000.0);
  41   SECONDARY_BEAM_DEMERIT = get_detail (details, ly_symbol2scm ("secondary-beam-demerit"), 10.0);
  42   STEM_LENGTH_DEMERIT_FACTOR = get_detail (details, ly_symbol2scm ("stem-length-demerit-factor"), 5);
  43   REGION_SIZE = get_detail (details, ly_symbol2scm ("region-size"), 2);
  44   BEAM_EPS = get_detail (details, ly_symbol2scm ("beam-eps"), 1e-3);
  45
  46   // possibly ridiculous, but too short stems just won't do
  47   STEM_LENGTH_LIMIT_PENALTY = get_detail (details, ly_symbol2scm ("stem-length-limit-penalty"), 5000);
  48   DAMPING_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("damping-direction-penalty"), 800);
  49   MUSICAL_DIRECTION_FACTOR = get_detail (details, ly_symbol2scm ("musical-direction-factor"), 400);
  50   IDEAL_SLOPE_FACTOR = get_detail (details, ly_symbol2scm ("ideal-slope-factor"), 10);
  51   ROUND_TO_ZERO_SLOPE = get_detail (details, ly_symbol2scm ("round-to-zero-slope"), 0.02);
  52 }
  53
  54 static Real
  55 shrink_extra_weight (Real x, Real fac)
  56 {
  57   return fabs (x) * ((x < 0) ? fac : 1.0);
  58 }
  59
  60 struct Quant_score
  61 {
  62   Real yl;
  63   Real yr;
  64   Real demerits;
  65
  66 #if DEBUG_BEAM_SCORING
  67   string score_card_;
  68 #endif
  69 };
  70
  71 /*
  72   TODO:
  73
  74   - Make all demerits customisable
  75
  76   - One sensible check per demerit (what's this --hwn)
  77
  78   - Add demerits for quants per se, as to forbid a specific quant
  79   entirely
  80 */
  81
  82 int
  83 best_quant_score_idx (vector<Quant_score> const &qscores)
  84 {
  85   Real best = 1e6;
  86   int best_idx = -1;
  87   for (vsize i = qscores.size (); i--;)
  88     {
  89       if (qscores[i].demerits < best)
  90         {
  91           best = qscores [i].demerits;
  92           best_idx = i;
  93         }
  94     }
  95
  96   return best_idx;
  97 }
  98
  99 MAKE_SCHEME_CALLBACK (Beam, quanting, 2);
 100 SCM
 101 Beam::quanting (SCM smob, SCM posns)
 102 {
 103   Grob *me = unsmob_grob (smob);
 104
 105   Beam_quant_parameters parameters;
 106   parameters.fill (me);
 107
 108   Real yl = scm_to_double (scm_car (posns));
 109   Real yr = scm_to_double (scm_cdr (posns));
 110
 111   /*
 112     Calculations are relative to a unit-scaled staff, i.e. the quants are
 113     divided by the current staff_space.
 114
 115   */
 116   Real ss = Staff_symbol_referencer::staff_space (me);
 117   Real thickness = Beam::get_thickness (me) / ss;
 118   Real slt = Staff_symbol_referencer::line_thickness (me) / ss;
 119
 120   Real dy_mus = robust_scm2double (me->get_property ("least-squares-dy"), 0);
 121   Real straddle = 0.0;
 122   Real sit = (thickness - slt) / 2;
 123   Real inter = 0.5;
 124   Real hang = 1.0 - (thickness - slt) / 2;
 125   Real quants [] = {straddle, sit, inter, hang };
 126
 127   int num_quants = int (sizeof (quants) / sizeof (Real));
 128   vector<Real> quantsl;
 129   vector<Real> quantsr;
 130
 131   /*
 132     going to REGION_SIZE == 2, yields another 0.6 second with
 133     wtk1-fugue2.
 134
 135
 136     (result indexes between 70 and 575)  ? --hwn.
 137
 138   */
 139
 140   /*
 141     Do stem computations.  These depend on YL and YR linearly, so we can
 142     precompute for every stem 2 factors.
 143   */
 144   vector<Grob*> stems
 145     = extract_grob_array (me, "stems");
 146   vector<Stem_info> stem_infos;
 147   vector<Real> base_lengths;
 148   vector<Real> stem_xposns;
 149
 150   Drul_array<bool> dirs_found (0, 0);
 151   Grob *common[2];
 152   for (int a = 2; a--;)
 153     common[a] = common_refpoint_of_array (stems, me, Axis (a));
 154
 155   Grob *fvs = first_visible_stem (me);
 156   Grob *lvs = last_visible_stem (me);
 157   Real xl = fvs ? fvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
 158   Real xr = fvs ? lvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
 159
 160   /*
 161     We store some info to quickly interpolate.
 162
 163     Sometimes my head is screwed on backwards.  The stemlength are
 164     AFFINE linear in YL and YR. If YL == YR == 0, then we might have
 165     stem_y != 0.0, when we're cross staff.
 166
 167   */
 168   for (vsize i = 0; i < stems.size (); i++)
 169     {
 170       Grob *s = stems[i];
 171
 172       Stem_info si (Stem::get_stem_info (s));
 173       si.scale (1 / ss);
 174       stem_infos.push_back (si);
 175       dirs_found[stem_infos.back ().dir_] = true;
 176
 177       bool f = to_boolean (s->get_property ("french-beaming"))
 178         && s != lvs && s != fvs;
 179
 180       base_lengths.push_back (calc_stem_y (me, s, common, xl, xr,
 181                                       Interval (0, 0), f) / ss);
 182       stem_xposns.push_back (s->relative_coordinate (common[X_AXIS], X_AXIS));
 183     }
 184
 185   bool xstaff = false;
 186   if (lvs && fvs)
 187     {
 188       Grob *commony = fvs->common_refpoint (lvs, Y_AXIS);
 189       xstaff = Align_interface::has_interface (commony);
 190     }
 191
 192   Direction ldir = Direction (stem_infos[0].dir_);
 193   Direction rdir = Direction (stem_infos.back ().dir_);
 194   bool is_knee = dirs_found[LEFT] && dirs_found[RIGHT];
 195
 196   int region_size = (int) parameters.REGION_SIZE;
 197
 198   /*
 199     Knees are harder, lets try some more possibilities for knees.
 200   */
 201   if (is_knee)
 202     region_size += 2;
 203
 204   /*
 205     Asymetry ? should run to <= region_size ?
 206   */
 207   for (int i = -region_size; i < region_size; i++)
 208     for (int j = 0; j < num_quants; j++)
 209       {
 210         quantsl.push_back (i + quants[j] + int (yl));
 211         quantsr.push_back (i + quants[j] + int (yr));
 212       }
 213
 214   vector<Quant_score> qscores;
 215
 216   for (vsize l = 0; l < quantsl.size (); l++)
 217     for (vsize r = 0; r < quantsr.size (); r++)
 218       {
 219         Quant_score qs;
 220         qs.yl = quantsl[l];
 221         qs.yr = quantsr[r];
 222         qs.demerits = 0.0;
 223
 224         qscores.push_back (qs);
 225       }
 226
 227   /* This is a longish function, but we don't separate this out into
 228      neat modular separate subfunctions, as the subfunctions would be
 229      called for many values of YL, YR. By precomputing various
 230      parameters outside of the loop, we can save a lot of time. */
 231   for (vsize i = qscores.size (); i--;)
 232     {
 233       Real d = score_slopes_dy (qscores[i].yl, qscores[i].yr,
 234                                 dy_mus, yr- yl,
 235                                 xr - xl,
 236                                 xstaff, &parameters);
 237       qscores[i].demerits += d;
 238
 239 #if DEBUG_BEAM_SCORING
 240       qscores[i].score_card_ += to_string ("S%.2f", d);
 241 #endif
 242     }
 243
 244   Real rad = Staff_symbol_referencer::staff_radius (me);
 245   Drul_array<int> edge_beam_counts
 246     (Stem::beam_multiplicity (stems[0]).length () + 1,
 247      Stem::beam_multiplicity (stems.back ()).length () + 1);
 248
 249   Real beam_translation = get_beam_translation (me) / ss;
 250
 251   Real reasonable_score = (is_knee) ? 200000 : 100;
 252   for (vsize i = qscores.size (); i--;)
 253     if (qscores[i].demerits < reasonable_score)
 254       {
 255         Real d = score_forbidden_quants (qscores[i].yl, qscores[i].yr,
 256                                          rad, slt, thickness, beam_translation,
 257                                          edge_beam_counts, ldir, rdir, &parameters);
 258         qscores[i].demerits += d;
 259
 260 #if DEBUG_BEAM_SCORING
 261         qscores[i].score_card_ += to_string (" F %.2f", d);
 262 #endif
 263       }
 264
 265   for (vsize i = qscores.size (); i--;)
 266     if (qscores[i].demerits < reasonable_score)
 267       {
 268         Real d = score_stem_lengths (stems, stem_infos,
 269                                      base_lengths, stem_xposns,
 270                                      xl, xr,
 271                                      is_knee,
 272                                      qscores[i].yl, qscores[i].yr, &parameters);
 273         qscores[i].demerits += d;
 274
 275 #if DEBUG_BEAM_SCORING
 276         qscores[i].score_card_ += to_string (" L %.2f", d);
 277 #endif
 278       }
 279
 280   int best_idx = best_quant_score_idx (qscores);
 281
 282 #if DEBUG_BEAM_SCORING
 283   SCM inspect_quants = me->get_property ("inspect-quants");
 284   if (to_boolean (me->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-scoring")))
 285       && scm_is_pair (inspect_quants))
 286     {
 287       Drul_array<Real> ins = ly_scm2interval (inspect_quants);
 288
 289       Real mindist = 1e6;
 290       for (vsize i = 0; i < qscores.size (); i++)
 291         {
 292           Real d = fabs (qscores[i].yl- ins[LEFT]) + fabs (qscores[i].yr - ins[RIGHT]);
 293           if (d < mindist)
 294             {
 295               best_idx = i;
 296               mindist = d;
 297             }
 298         }
 299       if (mindist > 1e5)
 300         programming_error ("can't find quant");
 301     }
 302 #endif
 303
 304   Interval final_positions;
 305   if (best_idx < 0)
 306     {
 307       warning (_ ("no feasible beam position"));
 308       final_positions = Interval (0, 0);
 309     }
 310   else
 311     {
 312       final_positions = Drul_array<Real> (qscores[best_idx].yl,
 313                                           qscores[best_idx].yr);
 314     }
 315
 316 #if DEBUG_BEAM_SCORING
 317   if (best_idx >= 0
 318       && to_boolean (me->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-scoring"))))
 319     {
 320       qscores[best_idx].score_card_ += to_string ("i%d", best_idx);
 321
 322       // debug quanting
 323       me->set_property ("quant-score",
 324                         scm_makfrom0str (qscores[best_idx].score_card_.c_str ()));
 325     }
 326 #endif
 327
 328   return ly_interval2scm (final_positions);
 329 }
 330
 331 Real
 332 Beam::score_stem_lengths (vector<Grob*> const &stems,
 333                           vector<Stem_info> const &stem_infos,
 334                           vector<Real> const &base_stem_ys,
 335                           vector<Real> const &stem_xs,
 336                           Real xl, Real xr,
 337                           bool knee,
 338                           Real yl, Real yr,
 339
 340                           Beam_quant_parameters const *parameters)
 341 {
 342   Real limit_penalty = parameters->STEM_LENGTH_LIMIT_PENALTY;
 343   Drul_array<Real> score (0, 0);
 344   Drul_array<int> count (0, 0);
 345
 346   for (vsize i = 0; i < stems.size (); i++)
 347     {
 348       Grob *s = stems[i];
 349       if (Stem::is_invisible (s))
 350         continue;
 351
 352       Real x = stem_xs[i];
 353       Real dx = xr - xl;
 354       Real beam_y = dx ? yr * (x - xl) / dx + yl * (xr - x) / dx : (yr + yl) / 2;
 355       Real current_y = beam_y + base_stem_ys[i];
 356       Real length_pen = parameters->STEM_LENGTH_DEMERIT_FACTOR;
 357
 358       Stem_info info = stem_infos[i];
 359       Direction d = info.dir_;
 360
 361       score[d] += limit_penalty * max (0.0, (d * (info.shortest_y_ - current_y)));
 362
 363       Real ideal_diff = d * (current_y - info.ideal_y_);
 364       Real ideal_score = shrink_extra_weight (ideal_diff, 1.5);
 365
 366       /* We introduce a power, to make the scoring strictly
 367          convex. Otherwise a symmetric knee beam (up/down/up/down)
 368          does not have an optimum in the middle. */
 369       if (knee)
 370         ideal_score = pow (ideal_score, 1.1);
 371
 372       score[d] += length_pen * ideal_score;
 373
 374       count[d]++;
 375     }
 376
 377   Direction d = DOWN;
 378   do
 379     score[d] /= max (count[d], 1);
 380   while (flip (&d) != DOWN)
 381     ;
 382
 383   return score[LEFT] + score[RIGHT];
 384 }
 385
 386 Real
 387 Beam::score_slopes_dy (Real yl, Real yr,
 388                        Real dy_mus, Real dy_damp,
 389                        Real dx,
 390                        bool xstaff,
 391
 392                        Beam_quant_parameters const *parameters)
 393 {
 394   Real dy = yr - yl;
 395   Real dem = 0.0;
 396
 397   /*
 398     DAMPING_DIRECTION_PENALTY is a very harsh measure, while for
 399     complex beaming patterns, horizontal is often a good choice.
 400
 401     TODO: find a way to incorporate the complexity of the beam in this
 402     penalty.
 403   */
 404   if (fabs (dy / dx) > parameters->ROUND_TO_ZERO_SLOPE
 405       && sign (dy_damp) != sign (dy))
 406     dem += parameters->DAMPING_DIRECTION_PENALTY;
 407
 408   dem += parameters->MUSICAL_DIRECTION_FACTOR
 409     * max (0.0, (fabs (dy) - fabs (dy_mus)));
 410
 411   Real slope_penalty = parameters->IDEAL_SLOPE_FACTOR;
 412
 413   /* Xstaff beams tend to use extreme slopes to get short stems. We
 414      put in a penalty here. */
 415   if (xstaff)
 416     slope_penalty *= 10;
 417
 418   /* Huh, why would a too steep beam be better than a too flat one ? */
 419   dem += shrink_extra_weight (fabs (dy_damp) - fabs (dy), 1.5)
 420     * slope_penalty;
 421
 422   return dem;
 423 }
 424
 425 static Real
 426 my_modf (Real x)
 427 {
 428   return x - floor (x);
 429 }
 430
 431 /*
 432   TODO: The fixed value SECONDARY_BEAM_DEMERIT is probably flawed:
 433   because for 32nd and 64th beams the forbidden quants are relatively
 434   more important than stem lengths.
 435 */
 436 Real
 437 Beam::score_forbidden_quants (Real yl, Real yr,
 438                               Real radius,
 439                               Real slt,
 440                               Real thickness, Real beam_translation,
 441                               Drul_array<int> beam_counts,
 442                               Direction ldir, Direction rdir,
 443
 444                               Beam_quant_parameters const *parameters)
 445 {
 446   Real dy = yr - yl;
 447   Drul_array<Real> y (yl, yr);
 448   Drul_array<Direction> dirs (ldir, rdir);
 449
 450   Real extra_demerit = parameters->SECONDARY_BEAM_DEMERIT / (max (beam_counts[LEFT], beam_counts[RIGHT]));
 451
 452   Direction d = LEFT;
 453   Real dem = 0.0;
 454   Real eps = parameters->BEAM_EPS;
 455
 456   do
 457     {
 458       for (int j = 1; j <= beam_counts[d]; j++)
 459         {
 460           Direction stem_dir = dirs[d];
 461
 462           /*
 463             The 2.2 factor is to provide a little leniency for
 464             borderline cases. If we do 2.0, then the upper outer line
 465             will be in the gap of the (2, sit) quant, leading to a
 466             false demerit.
 467           */
 468           Real gap1 = y[d] - stem_dir * ((j - 1) * beam_translation + thickness / 2 - slt / 2.2);
 469           Real gap2 = y[d] - stem_dir * (j * beam_translation - thickness / 2 + slt / 2.2);
 470
 471           Interval gap;
 472           gap.add_point (gap1);
 473           gap.add_point (gap2);
 474
 475           for (Real k = -radius;
 476                k <= radius + eps; k += 1.0)
 477             if (gap.contains (k))
 478               {
 479                 Real dist = min (fabs (gap[UP] - k), fabs (gap[DOWN] - k));
 480
 481                 /*
 482                   this parameter is tuned to grace-stem-length.ly
 483                 */
 484                 Real fixed_demerit = 0.4;
 485
 486                 dem += extra_demerit
 487                   * (fixed_demerit
 488                      + (1 - fixed_demerit) * (dist / gap.length ()) * 2);
 489               }
 490         }
 491     }
 492   while ((flip (&d)) != LEFT);
 493
 494   if (max (beam_counts[LEFT], beam_counts[RIGHT]) >= 2)
 495     {
 496       Real straddle = 0.0;
 497       Real sit = (thickness - slt) / 2;
 498       Real inter = 0.5;
 499       Real hang = 1.0 - (thickness - slt) / 2;
 500
 501       Direction d = LEFT;
 502       do
 503         {
 504           if (beam_counts[d] >= 2
 505               && fabs (y[d] - dirs[d] * beam_translation) < radius + inter)
 506             {
 507               if (dirs[d] == UP && dy <= eps
 508                   && fabs (my_modf (y[d]) - sit) < eps)
 509                 dem += extra_demerit;
 510
 511               if (dirs[d] == DOWN && dy >= eps
 512                   && fabs (my_modf (y[d]) - hang) < eps)
 513                 dem += extra_demerit;
 514             }
 515
 516           if (beam_counts[d] >= 3
 517               && fabs (y[d] - 2 * dirs[d] * beam_translation) < radius + inter)
 518             {
 519               if (dirs[d] == UP && dy <= eps
 520                   && fabs (my_modf (y[d]) - straddle) < eps)
 521                 dem += extra_demerit;
 522
 523               if (dirs[d] == DOWN && dy >= eps
 524                   && fabs (my_modf (y[d]) - straddle) < eps)
 525                 dem += extra_demerit;
 526             }
 527         }
 528       while (flip (&d) != LEFT);
 529     }
 530
 531   return dem;
 532 }
 533