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 "paper-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
  34 const Real BEAM_EPS = 1e-3;
  35
  36 // possibly ridiculous, but too short stems just won't do
  37 const int STEM_LENGTH_LIMIT_PENALTY = 5000;
  38 const int DAMPING_DIRECTION_PENALTY = 800;
  39 const int MUSICAL_DIRECTION_FACTOR = 400;
  40 const int IDEAL_SLOPE_FACTOR = 10;
  41
  42 const Real ROUND_TO_ZERO_SLOPE = 0.05;
  43 const int ROUND_TO_ZERO_POINTS = 4;
  44
  45 extern bool debug_beam_quanting_flag;
  46
  47 static Real
  48 shrink_extra_weight (Real x, Real fac)
  49 {
  50   return fabs (x) * ((x < 0) ? fac : 1.0);
  51 }
  52
  53
  54 struct Quant_score
  55 {
  56   Real yl;
  57   Real yr;
  58   Real demerits;
  59
  60 #if DEBUG_QUANTING
  61   String score_card_;
  62 #endif
  63 };
  64
  65
  66 /*
  67   TODO:
  68
  69    - Make all demerits customisable
  70
  71    - One sensible check per demerit (what's this --hwn)
  72
  73    - Add demerits for quants per se, as to forbid a specific quant
  74      entirely
  75
  76 */
  77
  78 int best_quant_score_idx (Array<Quant_score>  const & qscores)
  79 {
  80   Real best = 1e6;
  81   int best_idx = -1;
  82   for (int i = qscores.size (); i--;)
  83     {
  84       if (qscores[i].demerits < best)
  85         {
  86           best = qscores [i].demerits ;
  87           best_idx = i;
  88         }
  89     }
  90
  91   return best_idx;
  92 }
  93
  94 MAKE_SCHEME_CALLBACK (Beam, quanting, 1);
  95 SCM
  96 Beam::quanting (SCM smob)
  97 {
  98   Grob *me = unsmob_grob (smob);
  99
 100   SCM s = me->get_property ("positions");
 101   Real yl = gh_scm2double (gh_car (s));
 102   Real yr = gh_scm2double (gh_cdr (s));
 103
 104
 105   /*
 106     Calculations are relative to a unit-scaled staff, i.e. the quants are
 107     divided by the current staff_space.
 108
 109    */
 110   Real ss = Staff_symbol_referencer::staff_space (me);
 111   Real thickness = Beam::get_thickness (me) / ss ;
 112   Real slt = Staff_symbol_referencer::line_thickness (me) / ss;
 113
 114   SCM sdy = me->get_property ("least-squares-dy");
 115   Real dy_mus = gh_number_p (sdy) ? gh_scm2double (sdy) : 0.0;
 116
 117   Real straddle = 0.0;
 118   Real sit = (thickness - slt) / 2;
 119   Real inter = 0.5;
 120   Real hang = 1.0 - (thickness - slt) / 2;
 121   Real quants [] = {straddle, sit, inter, hang };
 122
 123   int num_quants = int (sizeof (quants)/sizeof (Real));
 124   Array<Real> quantsl;
 125   Array<Real> quantsr;
 126
 127   /*
 128     going to REGION_SIZE == 2, yields another 0.6 second with
 129     wtk1-fugue2.
 130
 131
 132     (result indexes between 70 and 575)  ? --hwn.
 133
 134   */
 135
 136
 137
 138   /*
 139     Do stem computations.  These depend on YL and YR linearly, so we can
 140     precompute for every stem 2 factors.
 141    */
 142   Link_array<Grob> stems=
 143     Pointer_group_interface__extract_grobs (me, (Grob*)0, "stems");
 144   Array<Stem_info> stem_infos;
 145   Array<Real> base_lengths;
 146   Array<Real> stem_xposns;
 147
 148   Drul_array<bool> dirs_found(0,0);
 149   Grob *common[2];
 150   for (int a = 2; a--;)
 151     common[a] = common_refpoint_of_array (stems, me, Axis(a));
 152
 153   Grob * fvs = first_visible_stem (me);
 154   Grob *lvs = last_visible_stem (me);
 155   Real xl = fvs ? fvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
 156   Real xr = fvs ? lvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
 157
 158   /*
 159     We store some info to quickly interpolate.
 160
 161     Sometimes my head is screwed on backwards.  The stemlength are
 162     AFFINE linear in YL and YR. If YL == YR == 0, then we might have
 163     stem_y != 0.0, when we're cross staff.
 164
 165    */
 166   for (int i= 0; i < stems.size(); i++)
 167     {
 168       Grob*s = stems[i];
 169
 170       Stem_info si (Stem::get_stem_info (s));
 171       si.scale (1 / ss);
 172       stem_infos.push (si);
 173       dirs_found[stem_infos.top ().dir_] = true;
 174
 175       bool f = to_boolean (s->get_property ("french-beaming"))
 176          && s != lvs && s != fvs;
 177
 178       base_lengths.push (calc_stem_y (me, s, common, xl, xr,
 179                                       Interval (0,0), f) / ss);
 180       stem_xposns.push (s->relative_coordinate (common[X_AXIS], X_AXIS));
 181     }
 182
 183   bool xstaff= false;
 184   if (lvs && fvs)
 185     {
 186       Grob *commony = fvs->common_refpoint (lvs, Y_AXIS);
 187       xstaff = Align_interface::has_interface (commony);
 188     }
 189
 190   Direction ldir = Direction (stem_infos[0].dir_);
 191   Direction rdir = Direction (stem_infos.top ().dir_);
 192   bool is_knee = dirs_found[LEFT] && dirs_found[RIGHT];
 193
 194
 195   int region_size = REGION_SIZE;
 196   /*
 197     Knees are harder, lets try some more possibilities for knees.
 198    */
 199   if (is_knee)
 200     region_size += 2;
 201
 202   /*
 203     Asymetry ? should run to <= region_size ?
 204    */
 205   for (int i = -region_size ; i < region_size; i++)
 206     for (int j = 0; j < num_quants; j++)
 207       {
 208         quantsl.push (i + quants[j] + int (yl));
 209         quantsr.push (i + quants[j] + int (yr));
 210       }
 211
 212   Array<Quant_score> qscores;
 213
 214   for (int l =0; l < quantsl.size (); l++)
 215     for (int r =0; r < quantsr.size (); r++)
 216       {
 217         Quant_score qs;
 218         qs.yl = quantsl[l];
 219         qs.yr = quantsr[r];
 220         qs.demerits = 0.0;
 221
 222         qscores.push (qs);
 223       }
 224
 225   /* This is a longish function, but we don't separate this out into
 226      neat modular separate subfunctions, as the subfunctions would be
 227      called for many values of YL, YR. By precomputing various
 228      parameters outside of the loop, we can save a lot of time. */
 229   for (int i = qscores.size (); i--;)
 230     {
 231       Real d =  score_slopes_dy (qscores[i].yl, qscores[i].yr,
 232                                  dy_mus, yr- yl,
 233                                  xr - xl,
 234                                  xstaff);
 235       qscores[i].demerits += d;
 236
 237 #if DEBUG_QUANTING
 238       qscores[i].score_card_ += to_string ("S%.2f",d);
 239 #endif
 240     }
 241
 242   Real rad = Staff_symbol_referencer::staff_radius (me);
 243   int beam_count = get_beam_count (me);
 244   Real beam_translation = get_beam_translation (me) / ss;
 245
 246   Real reasonable_score = (is_knee) ? 200000 : 100;
 247   for (int i = qscores.size (); i--;)
 248     if (qscores[i].demerits < reasonable_score)
 249       {
 250         Real d = score_forbidden_quants (qscores[i].yl, qscores[i].yr,
 251                                      rad, slt, thickness, beam_translation,
 252                                      beam_count, ldir, rdir);
 253         qscores[i].demerits += d;
 254
 255 #if DEBUG_QUANTING
 256         qscores[i].score_card_ += to_string (" F %.2f", d);
 257 #endif
 258       }
 259
 260   for (int i = qscores.size (); i--;)
 261     if (qscores[i].demerits < reasonable_score)
 262       {
 263         Real d=score_stem_lengths (stems, stem_infos,
 264                                  base_lengths, stem_xposns,
 265                                  xl, xr,
 266                                  is_knee,
 267                                  qscores[i].yl, qscores[i].yr);
 268         qscores[i].demerits +=  d;
 269
 270 #if DEBUG_QUANTING
 271         qscores[i].score_card_ += to_string (" L %.2f", d);
 272 #endif
 273       }
 274
 275   int best_idx = best_quant_score_idx (qscores);
 276
 277
 278 #if DEBUG_QUANTING
 279   SCM inspect_quants = me->get_property ("inspect-quants");
 280   if (debug_beam_quanting_flag
 281       && gh_pair_p (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 (debug_beam_quanting_flag)
 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    /*
 410      almost zero slopes look like errors in horizontal beams.
 411     */
 412    if (fabs (dy) > 1e-3
 413        && fabs (dy / dx) < ROUND_TO_ZERO_SLOPE)
 414      dem += ROUND_TO_ZERO_POINTS;
 415
 416    return dem;
 417 }
 418
 419 static Real
 420 my_modf (Real x)
 421 {
 422   return x - floor (x);
 423 }
 424
 425
 426 /*
 427   TODO: The fixed value SECONDARY_BEAM_DEMERIT is probably flawed:
 428    because for 32nd and 64th beams the forbidden quants are relatively
 429    more important than stem lengths.
 430 */
 431 Real
 432 Beam::score_forbidden_quants (Real yl, Real yr,
 433                               Real radius,
 434                               Real slt,
 435                               Real thickness, Real beam_translation,
 436                               int beam_count,
 437                               Direction ldir, Direction rdir)
 438 {
 439   Real dy = yr - yl;
 440   Drul_array<Real> y(yl,yr);
 441   Drul_array<Direction> dirs(ldir,rdir);
 442
 443   Real extra_demerit = SECONDARY_BEAM_DEMERIT / beam_count;
 444
 445   /*
 446     Inside the staff, inter quants are forbidden.
 447    */
 448   Real dem = 0.0;
 449   Direction d = LEFT;
 450   do
 451     {
 452       if (fabs (y[d]) <= (radius + 0.5) && fabs (my_modf (y[d]) - 0.5) < BEAM_EPS)
 453         dem += INTER_QUANT_PENALTY;
 454     }
 455   while ((flip (&d))!= LEFT);
 456
 457
 458   for (int j = 1; j <= beam_count; j++)
 459     {
 460       do
 461         {
 462           /*
 463             see if the outer staffline falls in a beam-gap
 464
 465             This test is too weak; we should really check all lines.
 466            */
 467           Direction stem_dir = dirs[d];
 468           Real gap1 =  y[d] - stem_dir * ((j-1) * beam_translation + thickness / 2 - slt/2 );
 469           Real gap2 = y[d] - stem_dir * (j * beam_translation - thickness / 2 + slt/2);
 470
 471           Interval gap;
 472           gap.add_point (gap1);
 473           gap.add_point (gap2);
 474
 475           for (Real k = - radius ;
 476                k <= radius + BEAM_EPS; k += 1.0)
 477             if (gap.contains (k))
 478               dem += extra_demerit;
 479         }
 480       while ((flip (&d))!= LEFT);
 481     }
 482
 483
 484
 485   // todo: use beam_count of outer stems.
 486   if (beam_count >= 2)
 487     {
 488       Real straddle = 0.0;
 489       Real sit = (thickness - slt) / 2;
 490       Real inter = 0.5;
 491       Real hang = 1.0 - (thickness - slt) / 2;
 492
 493       // hmm, without Interval/Drul_array, you get ~ 4x same code...
 494       if (fabs (y[LEFT] - dirs[LEFT] * beam_translation) < radius + inter)
 495         {
 496           if (dirs[LEFT] == UP && dy <= BEAM_EPS
 497               && fabs (my_modf (y[LEFT]) - sit) < BEAM_EPS)
 498             dem += extra_demerit;
 499
 500           if (dirs[LEFT] == DOWN && dy >= BEAM_EPS
 501               && fabs (my_modf (y[LEFT]) - hang) < BEAM_EPS)
 502             dem += extra_demerit;
 503         }
 504
 505       if (fabs (y[RIGHT] - dirs[RIGHT] * beam_translation) < radius + inter)
 506         {
 507           if (dirs[RIGHT] == UP && dy >= BEAM_EPS
 508               && fabs (my_modf (y[RIGHT]) - sit) < BEAM_EPS)
 509             dem += extra_demerit;
 510
 511           if (dirs[RIGHT] == DOWN && dy <= BEAM_EPS
 512               && fabs (my_modf (y[RIGHT]) - hang) < BEAM_EPS)
 513             dem += extra_demerit;
 514         }
 515
 516       if (beam_count >= 3)
 517         {
 518           if (fabs (y[LEFT] - 2 * dirs[LEFT] * beam_translation) < radius + inter)
 519             {
 520               if (dirs[LEFT] == UP && dy <= BEAM_EPS
 521                   && fabs (my_modf (y[LEFT]) - straddle) < BEAM_EPS)
 522                 dem += extra_demerit;
 523
 524               if (dirs[LEFT] == DOWN && dy >= BEAM_EPS
 525                   && fabs (my_modf (y[LEFT]) - straddle) < BEAM_EPS)
 526                 dem += extra_demerit;
 527             }
 528
 529           if (fabs (y[RIGHT] - 2 * dirs[RIGHT] * beam_translation) < radius + inter)
 530             {
 531               if (dirs[RIGHT] == UP && dy >= BEAM_EPS
 532                   && fabs (my_modf (y[RIGHT]) - straddle) < BEAM_EPS)
 533                 dem += extra_demerit;
 534
 535               if (dirs[RIGHT] == DOWN && dy <= BEAM_EPS
 536                   && fabs (my_modf (y[RIGHT]) - straddle) < BEAM_EPS)
 537                 dem += extra_demerit;
 538             }
 539         }
 540     }
 541
 542   return dem;
 543 }
 544
 545