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