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