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