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