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