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