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 "align-interface.hh"
  16 #include "international.hh"
  17 #include "output-def.hh"
  18 #include "pointer-group-interface.hh"
  19 #include "staff-symbol-referencer.hh"
  20 #include "stem.hh"
  21 #include "warn.hh"
  22
  23 Real
  24 get_detail (SCM alist, SCM sym, Real def)
  25 {
  26   SCM entry = scm_assq (sym, alist);
  27
  28   if (scm_is_pair (entry))
  29     return robust_scm2double (scm_cdr (entry), def);
  30   return def;
  31 }
  32
  33 void
  34 Beam_quant_parameters::fill (Grob *him)
  35 {
  36   SCM details = him->get_property ("details");
  37
  38   INTER_QUANT_PENALTY = get_detail (details, ly_symbol2scm ("inter-quant-penalty"), 1000.0);
  39   SECONDARY_BEAM_DEMERIT = get_detail (details, ly_symbol2scm ("secondary-beam-demerit"), 10.0);
  40   STEM_LENGTH_DEMERIT_FACTOR = get_detail (details, ly_symbol2scm ("stem-length-demerit-factor"), 5);
  41   REGION_SIZE = get_detail (details, ly_symbol2scm ("region-size"), 2);
  42   BEAM_EPS = get_detail (details, ly_symbol2scm ("beam-eps"), 1e-3);
  43
  44   // possibly ridiculous, but too short stems just won't do
  45   STEM_LENGTH_LIMIT_PENALTY = get_detail (details, ly_symbol2scm ("stem-length-limit-penalty"), 5000);
  46   DAMPING_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("damping-direction-penalty"), 800);
  47   MUSICAL_DIRECTION_FACTOR = get_detail (details, ly_symbol2scm ("musical-direction-factor"), 400);
  48   IDEAL_SLOPE_FACTOR = get_detail (details, ly_symbol2scm ("ideal-slope-factor"), 10);
  49   ROUND_TO_ZERO_SLOPE = get_detail (details, ly_symbol2scm ("round-to-zero-slope"), 0.02);
  50 }
  51
  52 static Real
  53 shrink_extra_weight (Real x, Real fac)
  54 {
  55   return fabs (x) * ((x < 0) ? fac : 1.0);
  56 }
  57
  58 struct Quant_score
  59 {
  60   Real yl;
  61   Real yr;
  62   Real demerits;
  63
  64 #if DEBUG_QUANTING
  65   std::string score_card_;
  66 #endif
  67 };
  68
  69 /*
  70   TODO:
  71
  72   - Make all demerits customisable
  73
  74   - One sensible check per demerit (what's this --hwn)
  75
  76   - Add demerits for quants per se, as to forbid a specific quant
  77   entirely
  78 */
  79
  80 int
  81 best_quant_score_idx (Array<Quant_score> const &qscores)
  82 {
  83   Real best = 1e6;
  84   int best_idx = -1;
  85   for (int i = qscores.size (); i--;)
  86     {
  87       if (qscores[i].demerits < best)
  88         {
  89           best = qscores [i].demerits;
  90           best_idx = i;
  91         }
  92     }
  93
  94   return best_idx;
  95 }
  96
  97 MAKE_SCHEME_CALLBACK (Beam, quanting, 2);
  98 SCM
  99 Beam::quanting (SCM smob, SCM posns)
 100 {
 101   Grob *me = unsmob_grob (smob);
 102
 103   Beam_quant_parameters parameters;
 104   parameters.fill (me);
 105
 106   Real yl = scm_to_double (scm_car (posns));
 107   Real yr = scm_to_double (scm_cdr (posns));
 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->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
 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_QUANTING
 317   if (best_idx >= 0
 318       && to_boolean (me->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-quanting"))))
 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 (Link_array<Grob> const &stems,
 333                           Array<Stem_info> const &stem_infos,
 334                           Array<Real> const &base_stem_ys,
 335                           Array<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 (int 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