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