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 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 extern bool debug_beam_quanting_flag;
  44
  45 static Real
  46 shrink_extra_weight (Real x, Real fac)
  47 {
  48   return fabs (x) * ((x < 0) ? fac : 1.0);
  49 }
  50
  51
  52 struct Quant_score
  53 {
  54   Real yl;
  55   Real yr;
  56   Real demerits;
  57
  58 #if DEBUG_QUANTING
  59   String score_card_;
  60 #endif
  61 };
  62
  63
  64 /*
  65   TODO:
  66
  67    - Make all demerits customisable
  68
  69    - One sensible check per demerit (what's this --hwn)
  70
  71    - Add demerits for quants per se, as to forbid a specific quant
  72      entirely
  73
  74 */
  75
  76 int best_quant_score_idx (Array<Quant_score>  const & qscores)
  77 {
  78   Real best = 1e6;
  79   int best_idx = -1;
  80   for (int i = qscores.size (); i--;)
  81     {
  82       if (qscores[i].demerits < best)
  83         {
  84           best = qscores [i].demerits ;
  85           best_idx = i;
  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 = gh_scm2double (gh_car (s));
 100   Real yr = gh_scm2double (gh_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   SCM sdy = me->get_property ("least-squares-dy");
 113   Real dy_mus = gh_number_p (sdy) ? gh_scm2double (sdy) : 0.0;
 114
 115   Real straddle = 0.0;
 116   Real sit = (thickness - slt) / 2;
 117   Real inter = 0.5;
 118   Real hang = 1.0 - (thickness - slt) / 2;
 119   Real quants [] = {straddle, sit, inter, hang };
 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   int beam_count = get_beam_count (me);
 242   Real beam_translation = get_beam_translation (me) / ss;
 243
 244   Real reasonable_score = (is_knee) ? 200000 : 100;
 245   for (int i = qscores.size (); i--;)
 246     if (qscores[i].demerits < reasonable_score)
 247       {
 248         Real d = score_forbidden_quants (qscores[i].yl, qscores[i].yr,
 249                                      rad, slt, thickness, beam_translation,
 250                                      beam_count, ldir, rdir);
 251         qscores[i].demerits += d;
 252
 253 #if DEBUG_QUANTING
 254         qscores[i].score_card_ += to_string (" F %.2f", d);
 255 #endif
 256       }
 257
 258   for (int i = qscores.size (); i--;)
 259     if (qscores[i].demerits < reasonable_score)
 260       {
 261         Real d=score_stem_lengths (stems, stem_infos,
 262                                  base_lengths, stem_xposns,
 263                                  xl, xr,
 264                                  is_knee,
 265                                  qscores[i].yl, qscores[i].yr);
 266         qscores[i].demerits +=  d;
 267
 268 #if DEBUG_QUANTING
 269         qscores[i].score_card_ += to_string (" L %.2f", d);
 270 #endif
 271       }
 272
 273   int best_idx = best_quant_score_idx (qscores);
 274
 275
 276 #if DEBUG_QUANTING
 277   SCM inspect_quants = me->get_property ("inspect-quants");
 278   if (debug_beam_quanting_flag
 279       && gh_pair_p (inspect_quants))
 280     {
 281       Drul_array<Real> ins = ly_scm2interval (inspect_quants);
 282
 283       int i = 0;
 284
 285       Real mindist = 1e6;
 286       for (; i < qscores.size(); i ++)
 287         {
 288           Real d =fabs (qscores[i].yl- ins[LEFT]) + fabs (qscores[i].yr - ins[RIGHT]);
 289           if (d < mindist)
 290             {
 291               best_idx = i;
 292               mindist= d;
 293             }
 294         }
 295       if (mindist > 1e5)
 296         programming_error ("Could not find quant.");
 297     }
 298 #endif
 299
 300   me->set_property ("positions",
 301                          ly_interval2scm (Drul_array<Real> (qscores[best_idx].yl,
 302                                           qscores[best_idx].yr)));
 303 #if DEBUG_QUANTING
 304   if (debug_beam_quanting_flag)
 305     {
 306       qscores[best_idx].score_card_ += to_string ("i%d", best_idx);
 307
 308       // debug quanting
 309       me->set_property ("quant-score",
 310                              scm_makfrom0str (qscores[best_idx].score_card_.to_str0 ()));
 311     }
 312 #endif
 313
 314   return SCM_UNSPECIFIED;
 315 }
 316
 317 Real
 318 Beam::score_stem_lengths (Link_array<Grob> const &stems,
 319                           Array<Stem_info> const &stem_infos,
 320                           Array<Real> const &base_stem_ys,
 321                           Array<Real> const &stem_xs,
 322                           Real xl, Real xr,
 323                           bool knee,
 324                           Real yl, Real yr)
 325 {
 326   Real limit_penalty = STEM_LENGTH_LIMIT_PENALTY;
 327   Drul_array<Real> score (0, 0);
 328   Drul_array<int> count (0, 0);
 329
 330   for (int i=0; i < stems.size (); i++)
 331     {
 332       Grob* s = stems[i];
 333       if (Stem::is_invisible (s))
 334         continue;
 335
 336       Real x = stem_xs[i];
 337       Real dx = xr-xl;
 338       Real beam_y = dx ? yr *(x - xl)/dx + yl * ( xr - x)/dx : (yr + yl)/2;
 339       Real current_y = beam_y + base_stem_ys[i];
 340       Real length_pen = STEM_LENGTH_DEMERIT_FACTOR;
 341
 342       Stem_info info = stem_infos[i];
 343       Direction d = info.dir_;
 344
 345       score[d] += limit_penalty * (0 >? (d * (info.shortest_y_ - current_y)));
 346
 347       Real ideal_diff = d * (current_y - info.ideal_y_);
 348       Real ideal_score = shrink_extra_weight (ideal_diff, 1.5);
 349
 350       /* We introduce a power, to make the scoring strictly
 351          convex. Otherwise a symmetric knee beam (up/down/up/down)
 352          does not have an optimum in the middle. */
 353       if (knee)
 354         ideal_score = pow (ideal_score, 1.1);
 355
 356       score[d] += length_pen * ideal_score;
 357
 358       count[d] ++;
 359     }
 360
 361   Direction d = DOWN;
 362   do
 363     {
 364       score[d] /= (count[d] >? 1);
 365     }
 366   while (flip (&d) != DOWN);
 367
 368   return score[LEFT]+score[RIGHT];
 369 }
 370
 371 Real
 372 Beam::score_slopes_dy (Real yl, Real yr,
 373                        Real dy_mus, Real dy_damp,
 374                        Real dx,
 375                        bool xstaff)
 376 {
 377   Real dy = yr - yl;
 378   Real dem = 0.0;
 379
 380   /*
 381     DAMPING_DIRECTION_PENALTY is a very harsh measure, while for
 382     complex beaming patterns, horizontal is often a good choice.
 383
 384     TODO: find a way to incorporate the complexity of the beam in this
 385     penalty.
 386    */
 387   if (fabs (dy/dx) > ROUND_TO_ZERO_SLOPE
 388       && sign (dy_damp) != sign (dy))
 389     {
 390       dem += DAMPING_DIRECTION_PENALTY;
 391     }
 392
 393    dem += MUSICAL_DIRECTION_FACTOR * (0 >? (fabs (dy) - fabs (dy_mus)));
 394
 395
 396    Real slope_penalty = IDEAL_SLOPE_FACTOR;
 397
 398    /* Xstaff beams tend to use extreme slopes to get short stems. We
 399       put in a penalty here. */
 400    if (xstaff)
 401      slope_penalty *= 10;
 402
 403    /* Huh, why would a too steep beam be better than a too flat one ? */
 404    dem += shrink_extra_weight (fabs (dy_damp) - fabs (dy), 1.5)
 405      * slope_penalty;
 406
 407 #if 0
 408    /*
 409      almost zero slopes look like errors in horizontal beams.
 410     */
 411    /*
 412      This causes too much problems, because horizontal depends on
 413      horizontal spacing details.  These errors should be dealt with
 414      through concaveness. --hwn.
 415     */
 416    if (fabs (dy) > 1e-3
 417        && fabs (dy / dx) < ROUND_TO_ZERO_SLOPE)
 418      dem += ROUND_TO_ZERO_POINTS;
 419 #endif
 420
 421    return dem;
 422 }
 423
 424 static Real
 425 my_modf (Real x)
 426 {
 427   return x - floor (x);
 428 }
 429
 430
 431 /*
 432   TODO: The fixed value SECONDARY_BEAM_DEMERIT is probably flawed:
 433    because for 32nd and 64th beams the forbidden quants are relatively
 434    more important than stem lengths.
 435 */
 436 Real
 437 Beam::score_forbidden_quants (Real yl, Real yr,
 438                               Real radius,
 439                               Real slt,
 440                               Real thickness, Real beam_translation,
 441                               int beam_count,
 442                               Direction ldir, Direction rdir)
 443 {
 444   Real dy = yr - yl;
 445   Drul_array<Real> y(yl,yr);
 446   Drul_array<Direction> dirs(ldir,rdir);
 447
 448   Real extra_demerit = SECONDARY_BEAM_DEMERIT / beam_count;
 449
 450   /*
 451     Inside the staff, inter quants are forbidden.
 452    */
 453   Real dem = 0.0;
 454   Direction d = LEFT;
 455   do
 456     {
 457       if (fabs (y[d]) <= (radius + 0.5) && fabs (my_modf (y[d]) - 0.5) < BEAM_EPS)
 458         dem += INTER_QUANT_PENALTY;
 459     }
 460   while ((flip (&d))!= LEFT);
 461
 462
 463   for (int j = 1; j <= beam_count; j++)
 464     {
 465       do
 466         {
 467           /*
 468             see if the outer staffline falls in a beam-gap
 469
 470             This test is too weak; we should really check all lines.
 471            */
 472           Direction stem_dir = dirs[d];
 473           Real gap1 =  y[d] - stem_dir * ((j-1) * beam_translation + thickness / 2 - slt/2 );
 474           Real gap2 = y[d] - stem_dir * (j * beam_translation - thickness / 2 + slt/2);
 475
 476           Interval gap;
 477           gap.add_point (gap1);
 478           gap.add_point (gap2);
 479
 480           for (Real k = - radius ;
 481                k <= radius + BEAM_EPS; k += 1.0)
 482             if (gap.contains (k))
 483               dem += extra_demerit;
 484         }
 485       while ((flip (&d))!= LEFT);
 486     }
 487
 488
 489
 490   // todo: use beam_count of outer stems.
 491   if (beam_count >= 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       // hmm, without Interval/Drul_array, you get ~ 4x same code...
 499       if (fabs (y[LEFT] - dirs[LEFT] * beam_translation) < radius + inter)
 500         {
 501           if (dirs[LEFT] == UP && dy <= BEAM_EPS
 502               && fabs (my_modf (y[LEFT]) - sit) < BEAM_EPS)
 503             dem += extra_demerit;
 504
 505           if (dirs[LEFT] == DOWN && dy >= BEAM_EPS
 506               && fabs (my_modf (y[LEFT]) - hang) < BEAM_EPS)
 507             dem += extra_demerit;
 508         }
 509
 510       if (fabs (y[RIGHT] - dirs[RIGHT] * beam_translation) < radius + inter)
 511         {
 512           if (dirs[RIGHT] == UP && dy >= BEAM_EPS
 513               && fabs (my_modf (y[RIGHT]) - sit) < BEAM_EPS)
 514             dem += extra_demerit;
 515
 516           if (dirs[RIGHT] == DOWN && dy <= BEAM_EPS
 517               && fabs (my_modf (y[RIGHT]) - hang) < BEAM_EPS)
 518             dem += extra_demerit;
 519         }
 520
 521       if (beam_count >= 3)
 522         {
 523           if (fabs (y[LEFT] - 2 * dirs[LEFT] * beam_translation) < radius + inter)
 524             {
 525               if (dirs[LEFT] == UP && dy <= BEAM_EPS
 526                   && fabs (my_modf (y[LEFT]) - straddle) < BEAM_EPS)
 527                 dem += extra_demerit;
 528
 529               if (dirs[LEFT] == DOWN && dy >= BEAM_EPS
 530                   && fabs (my_modf (y[LEFT]) - straddle) < BEAM_EPS)
 531                 dem += extra_demerit;
 532             }
 533
 534           if (fabs (y[RIGHT] - 2 * dirs[RIGHT] * beam_translation) < radius + inter)
 535             {
 536               if (dirs[RIGHT] == UP && dy >= BEAM_EPS
 537                   && fabs (my_modf (y[RIGHT]) - straddle) < BEAM_EPS)
 538                 dem += extra_demerit;
 539
 540               if (dirs[RIGHT] == DOWN && dy <= BEAM_EPS
 541                   && fabs (my_modf (y[RIGHT]) - straddle) < BEAM_EPS)
 542                 dem += extra_demerit;
 543             }
 544         }
 545     }
 546
 547   return dem;
 548 }
 549
 550