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
  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       Real il = gh_scm2double (gh_car (inspect_quants));
 278       Real ir = gh_scm2double (gh_cdr (inspect_quants));
 279
 280       int i = 0;
 281
 282       Real mindist = 1e6;
 283       for (; i < qscores.size(); i ++)
 284         {
 285           Real d =fabs (qscores[i].yl-il) + fabs (qscores[i].yr - ir);
 286           if (d < mindist)
 287             {
 288               best_idx = i;
 289               mindist= d;
 290             }
 291         }
 292       if (mindist > 1e5)
 293         programming_error ("Could not find quant.");
 294     }
 295 #endif
 296
 297   me->set_grob_property ("positions",
 298                          gh_cons (gh_double2scm (qscores[best_idx].yl),
 299                                   gh_double2scm (qscores[best_idx].yr)));
 300 #if DEBUG_QUANTING
 301   if (debug_beam_quanting_flag)
 302     {
 303       qscores[best_idx].score_card_ += to_string ("i%d", best_idx);
 304
 305       // debug quanting
 306       me->set_grob_property ("quant-score",
 307                              scm_makfrom0str (qscores[best_idx].score_card_.to_str0 ()));
 308     }
 309 #endif
 310
 311   return SCM_UNSPECIFIED;
 312 }
 313
 314 Real
 315 Beam::score_stem_lengths (Link_array<Grob> const &stems,
 316                           Array<Stem_info> const &stem_infos,
 317                           Array<Real> const &base_stem_ys,
 318                           Array<Real> const &stem_xs,
 319                           Real xl, Real xr,
 320                           bool knee,
 321                           Real yl, Real yr)
 322 {
 323   Real limit_penalty = STEM_LENGTH_LIMIT_PENALTY;
 324   Drul_array<Real> score (0, 0);
 325   Drul_array<int> count (0, 0);
 326
 327   for (int i=0; i < stems.size (); i++)
 328     {
 329       Grob* s = stems[i];
 330       if (Stem::invisible_b (s))
 331         continue;
 332
 333       Real x = stem_xs[i];
 334       Real dx = xr-xl;
 335       Real beam_y = dx ? yr *(x - xl)/dx + yl * ( xr - x)/dx : (yr + yl)/2;
 336       Real current_y = beam_y + base_stem_ys[i];
 337       Real length_pen = STEM_LENGTH_DEMERIT_FACTOR;
 338
 339       Stem_info info = stem_infos[i];
 340       Direction d = info.dir_;
 341
 342       score[d] += limit_penalty * (0 >? (d * (info.shortest_y_ - current_y)));
 343
 344       Real ideal_diff = d * (current_y - info.ideal_y_);
 345       Real ideal_score = shrink_extra_weight (ideal_diff, 1.5);
 346
 347       /* We introduce a power, to make the scoring strictly
 348          convex. Otherwise a symmetric knee beam (up/down/up/down)
 349          does not have an optimum in the middle. */
 350       if (knee)
 351         ideal_score = pow (ideal_score, 1.1);
 352
 353       score[d] += length_pen * ideal_score;
 354
 355       count[d] ++;
 356     }
 357
 358   Direction d = DOWN;
 359   do
 360     {
 361       score[d] /= (count[d] >? 1);
 362     }
 363   while (flip (&d) != DOWN);
 364
 365   return score[LEFT]+score[RIGHT];
 366 }
 367
 368 Real
 369 Beam::score_slopes_dy (Real yl, Real yr,
 370                        Real dy_mus, Real dy_damp,
 371                        Real dx,
 372                        bool xstaff)
 373 {
 374   Real dy = yr - yl;
 375   Real dem = 0.0;
 376
 377   /*
 378     DAMPING_DIRECTION_PENALTY is a very harsh measure, while for
 379     complex beaming patterns, horizontal is often a good choice.
 380
 381     TODO: find a way to incorporate the complexity of the beam in this
 382     penalty.
 383    */
 384   if (fabs (dy/dx) > ROUND_TO_ZERO_SLOPE
 385       && sign (dy_damp) != sign (dy))
 386     {
 387       dem += DAMPING_DIRECTION_PENALTY;
 388     }
 389
 390    dem += MUSICAL_DIRECTION_FACTOR * (0 >? (fabs (dy) - fabs (dy_mus)));
 391
 392
 393    Real slope_penalty = IDEAL_SLOPE_FACTOR;
 394
 395    /* Xstaff beams tend to use extreme slopes to get short stems. We
 396       put in a penalty here. */
 397    if (xstaff)
 398      slope_penalty *= 10;
 399
 400    /* Huh, why would a too steep beam be better than a too flat one ? */
 401    dem += shrink_extra_weight (fabs (dy_damp) - fabs (dy), 1.5)
 402      * slope_penalty;
 403
 404    /*
 405      almost zero slopes look like errors in horizontal beams.
 406     */
 407    if (fabs (dy) > 1e-3
 408        && fabs (dy / dx) < ROUND_TO_ZERO_SLOPE)
 409      dem += ROUND_TO_ZERO_POINTS;
 410
 411    return dem;
 412 }
 413
 414 static Real
 415 my_modf (Real x)
 416 {
 417   return x - floor (x);
 418 }
 419
 420
 421 /*
 422   TODO: The fixed value SECONDARY_BEAM_DEMERIT is probably flawed:
 423    because for 32nd and 64th beams the forbidden quants are relatively
 424    more important than stem lengths.
 425 */
 426 Real
 427 Beam::score_forbidden_quants (Real yl, Real yr,
 428                               Real radius,
 429                               Real slt,
 430                               Real thickness, Real beam_translation,
 431                               int beam_count,
 432                               Direction ldir, Direction rdir)
 433 {
 434   Real dy = yr - yl;
 435   Drul_array<Real> y(yl,yr);
 436   Drul_array<Direction> dirs(ldir,rdir);
 437
 438   Real extra_demerit = SECONDARY_BEAM_DEMERIT / beam_count;
 439
 440   /*
 441     Inside the staff, inter quants are forbidden.
 442    */
 443   Real dem = 0.0;
 444   Direction d = LEFT;
 445   do
 446     {
 447       if (fabs (y[d]) <= (radius + 0.5) && fabs (my_modf (y[d]) - 0.5) < 1e-3)
 448         dem += INTER_QUANT_PENALTY;
 449     }
 450   while ((flip (&d))!= LEFT);
 451
 452
 453   for (int j = 1; j <= beam_count; j++)
 454     {
 455       do
 456         {
 457           /*
 458             see if the outer staffline falls in a beam-gap
 459
 460             This test is too weak; we should really check all lines.
 461            */
 462           Direction stem_dir = dirs[d];
 463           Real gap1 =  y[d] - stem_dir * ((j-1) * beam_translation + thickness / 2 - slt/2 );
 464           Real gap2 = y[d] - stem_dir * (j * beam_translation - thickness / 2 + slt/2);
 465
 466           Interval gap;
 467           gap.add_point (gap1);
 468           gap.add_point (gap2);
 469
 470           if (gap.contains (radius))
 471             dem += extra_demerit;
 472         }
 473       while ((flip (&d))!= LEFT);
 474     }
 475
 476
 477
 478   // todo: use beam_count of outer stems.
 479   if (beam_count >= 2)
 480     {
 481       Real straddle = 0.0;
 482       Real sit = (thickness - slt) / 2;
 483       Real inter = 0.5;
 484       Real hang = 1.0 - (thickness - slt) / 2;
 485
 486       Real eps = 1e-3;
 487
 488       // hmm, without Interval/Drul_array, you get ~ 4x same code...
 489       if (fabs (y[LEFT] - dirs[LEFT] * beam_translation) < radius + inter)
 490         {
 491           if (dirs[LEFT] == UP && dy <= eps
 492               && fabs (my_modf (y[LEFT]) - sit) < eps)
 493             dem += extra_demerit;
 494
 495           if (dirs[LEFT] == DOWN && dy >= eps
 496               && fabs (my_modf (y[LEFT]) - hang) < eps)
 497             dem += extra_demerit;
 498         }
 499
 500       if (fabs (y[RIGHT] - dirs[RIGHT] * beam_translation) < radius + inter)
 501         {
 502           if (dirs[RIGHT] == UP && dy >= eps
 503               && fabs (my_modf (y[RIGHT]) - sit) < eps)
 504             dem += extra_demerit;
 505
 506           if (dirs[RIGHT] == DOWN && dy <= eps
 507               && fabs (my_modf (y[RIGHT]) - hang) < eps)
 508             dem += extra_demerit;
 509         }
 510
 511       if (beam_count >= 3)
 512         {
 513           if (fabs (y[LEFT] - 2 * dirs[LEFT] * beam_translation) < radius + inter)
 514             {
 515               if (dirs[LEFT] == UP && dy <= eps
 516                   && fabs (my_modf (y[LEFT]) - straddle) < eps)
 517                 dem += extra_demerit;
 518
 519               if (dirs[LEFT] == DOWN && dy >= eps
 520                   && fabs (my_modf (y[LEFT]) - straddle) < eps)
 521                 dem += extra_demerit;
 522             }
 523
 524           if (fabs (y[RIGHT] - 2 * dirs[RIGHT] * beam_translation) < radius + inter)
 525             {
 526               if (dirs[RIGHT] == UP && dy >= eps
 527                   && fabs (my_modf (y[RIGHT]) - straddle) < eps)
 528                 dem += extra_demerit;
 529
 530               if (dirs[RIGHT] == DOWN && dy <= eps
 531                   && fabs (my_modf (y[RIGHT]) - straddle) < eps)
 532                 dem += extra_demerit;
 533             }
 534         }
 535     }
 536
 537   return dem;
 538 }
 539
 540