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   if (best_idx < 0)
  81     {
  82       programming_error ("no best beam quant score");
  83       best_idx = 0;
  84     }
  85
  86   return best_idx;
  87 }
  88
  89 MAKE_SCHEME_CALLBACK (Beam, quanting, 1);
  90 SCM
  91 Beam::quanting (SCM smob)
  92 {
  93   Grob *me = unsmob_grob (smob);
  94
  95   SCM s = me->get_property ("positions");
  96   Real yl = scm_to_double (scm_car (s));
  97   Real yr = scm_to_double (scm_cdr (s));
  98
  99   /*
 100     Calculations are relative to a unit-scaled staff, i.e. the quants are
 101     divided by the current staff_space.
 102
 103   */
 104   Real ss = Staff_symbol_referencer::staff_space (me);
 105   Real thickness = Beam::get_thickness (me) / ss;
 106   Real slt = Staff_symbol_referencer::line_thickness (me) / ss;
 107
 108   Real dy_mus = robust_scm2double (me->get_property ("least-squares-dy"), 0);
 109   Real straddle = 0.0;
 110   Real sit = (thickness - slt) / 2;
 111   Real inter = 0.5;
 112   Real hang = 1.0 - (thickness - slt) / 2;
 113   Real quants [] = {straddle, sit, inter, hang };
 114
 115   int num_quants = int (sizeof (quants) / sizeof (Real));
 116   Array<Real> quantsl;
 117   Array<Real> quantsr;
 118
 119   /*
 120     going to REGION_SIZE == 2, yields another 0.6 second with
 121     wtk1-fugue2.
 122
 123
 124     (result indexes between 70 and 575)  ? --hwn.
 125
 126   */
 127
 128   /*
 129     Do stem computations.  These depend on YL and YR linearly, so we can
 130     precompute for every stem 2 factors.
 131   */
 132   Link_array<Grob> stems
 133     = extract_grob_array (me, ly_symbol2scm ("stems"));
 134   Array<Stem_info> stem_infos;
 135   Array<Real> base_lengths;
 136   Array<Real> stem_xposns;
 137
 138   Drul_array<bool> dirs_found (0, 0);
 139   Grob *common[2];
 140   for (int a = 2; a--;)
 141     common[a] = common_refpoint_of_array (stems, me, Axis (a));
 142
 143   Grob *fvs = first_visible_stem (me);
 144   Grob *lvs = last_visible_stem (me);
 145   Real xl = fvs ? fvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
 146   Real xr = fvs ? lvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
 147
 148   /*
 149     We store some info to quickly interpolate.
 150
 151     Sometimes my head is screwed on backwards.  The stemlength are
 152     AFFINE linear in YL and YR. If YL == YR == 0, then we might have
 153     stem_y != 0.0, when we're cross staff.
 154
 155   */
 156   for (int i = 0; i < stems.size (); i++)
 157     {
 158       Grob *s = stems[i];
 159
 160       Stem_info si (Stem::get_stem_info (s));
 161       si.scale (1 / ss);
 162       stem_infos.push (si);
 163       dirs_found[stem_infos.top ().dir_] = true;
 164
 165       bool f = to_boolean (s->get_property ("french-beaming"))
 166         && s != lvs && s != fvs;
 167
 168       base_lengths.push (calc_stem_y (me, s, common, xl, xr,
 169                                       Interval (0, 0), f) / ss);
 170       stem_xposns.push (s->relative_coordinate (common[X_AXIS], X_AXIS));
 171     }
 172
 173   bool xstaff = false;
 174   if (lvs && fvs)
 175     {
 176       Grob *commony = fvs->common_refpoint (lvs, Y_AXIS);
 177       xstaff = Align_interface::has_interface (commony);
 178     }
 179
 180   Direction ldir = Direction (stem_infos[0].dir_);
 181   Direction rdir = Direction (stem_infos.top ().dir_);
 182   bool is_knee = dirs_found[LEFT] && dirs_found[RIGHT];
 183
 184   int region_size = REGION_SIZE;
 185   /*
 186     Knees are harder, lets try some more possibilities for knees.
 187   */
 188   if (is_knee)
 189     region_size += 2;
 190
 191   /*
 192     Asymetry ? should run to <= region_size ?
 193   */
 194   for (int i = -region_size; i < region_size; i++)
 195     for (int j = 0; j < num_quants; j++)
 196       {
 197         quantsl.push (i + quants[j] + int (yl));
 198         quantsr.push (i + quants[j] + int (yr));
 199       }
 200
 201   Array<Quant_score> qscores;
 202
 203   for (int l = 0; l < quantsl.size (); l++)
 204     for (int r = 0; r < quantsr.size (); r++)
 205       {
 206         Quant_score qs;
 207         qs.yl = quantsl[l];
 208         qs.yr = quantsr[r];
 209         qs.demerits = 0.0;
 210
 211         qscores.push (qs);
 212       }
 213
 214   /* This is a longish function, but we don't separate this out into
 215      neat modular separate subfunctions, as the subfunctions would be
 216      called for many values of YL, YR. By precomputing various
 217      parameters outside of the loop, we can save a lot of time. */
 218   for (int i = qscores.size (); i--;)
 219     {
 220       Real d = score_slopes_dy (qscores[i].yl, qscores[i].yr,
 221                                 dy_mus, yr- yl,
 222                                 xr - xl,
 223                                 xstaff);
 224       qscores[i].demerits += d;
 225
 226 #if DEBUG_QUANTING
 227       qscores[i].score_card_ += to_string ("S%.2f", d);
 228 #endif
 229     }
 230
 231   Real rad = Staff_symbol_referencer::staff_radius (me);
 232   Drul_array<int> edge_beam_counts
 233     (Stem::beam_multiplicity (stems[0]).length () + 1,
 234      Stem::beam_multiplicity (stems.top ()).length () + 1);
 235
 236   Real beam_translation = get_beam_translation (me) / ss;
 237
 238   Real reasonable_score = (is_knee) ? 200000 : 100;
 239   for (int i = qscores.size (); i--;)
 240     if (qscores[i].demerits < reasonable_score)
 241       {
 242         Real d = score_forbidden_quants (qscores[i].yl, qscores[i].yr,
 243                                          rad, slt, thickness, beam_translation,
 244                                          edge_beam_counts, ldir, rdir);
 245         qscores[i].demerits += d;
 246
 247 #if DEBUG_QUANTING
 248         qscores[i].score_card_ += to_string (" F %.2f", d);
 249 #endif
 250       }
 251
 252   for (int i = qscores.size (); i--;)
 253     if (qscores[i].demerits < reasonable_score)
 254       {
 255         Real d = score_stem_lengths (stems, stem_infos,
 256                                      base_lengths, stem_xposns,
 257                                      xl, xr,
 258                                      is_knee,
 259                                      qscores[i].yl, qscores[i].yr);
 260         qscores[i].demerits += d;
 261
 262 #if DEBUG_QUANTING
 263         qscores[i].score_card_ += to_string (" L %.2f", d);
 264 #endif
 265       }
 266
 267   int best_idx = best_quant_score_idx (qscores);
 268
 269 #if DEBUG_QUANTING
 270   SCM inspect_quants = me->get_property ("inspect-quants");
 271   if (to_boolean (me->get_layout ()->lookup_variable (ly_symbol2scm ("debug-beam-quanting")))
 272       && scm_is_pair (inspect_quants))
 273     {
 274       Drul_array<Real> ins = ly_scm2interval (inspect_quants);
 275
 276       int i = 0;
 277
 278       Real mindist = 1e6;
 279       for (; i < qscores.size (); i++)
 280         {
 281           Real d = fabs (qscores[i].yl- ins[LEFT]) + fabs (qscores[i].yr - ins[RIGHT]);
 282           if (d < mindist)
 283             {
 284               best_idx = i;
 285               mindist = d;
 286             }
 287         }
 288       if (mindist > 1e5)
 289         programming_error ("can't find quant");
 290     }
 291 #endif
 292
 293   me->set_property ("positions",
 294                     ly_interval2scm (Drul_array<Real> (qscores[best_idx].yl,
 295                                                        qscores[best_idx].yr)));
 296 #if DEBUG_QUANTING
 297   if (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