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