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