]> git.donarmstrong.com Git - lilypond.git/blobdiff - lily/beam-quanting.cc
* The grand 2005-2006 replace.
[lilypond.git] / lily / beam-quanting.cc
index 550a0a877c92197c8ce71967618b9e793513584b..90ace1f4012a5dbd2dc2ea71943c21dce1127b2d 100644 (file)
@@ -1,40 +1,52 @@
 /*
   beam-quanting.cc -- implement Beam quanting functions
-  
+
   source file of the GNU LilyPond music typesetter
-  
-  (c)  1997--2003 Han-Wen Nienhuys <hanwen@cs.uu.nl>
+
+  (c) 1997--2006 Han-Wen Nienhuys <hanwen@xs4all.nl>
   Jan Nieuwenhuizen <janneke@gnu.org>
-  
 */
 
+#include "beam.hh"
 
-
-#include <math.h>
+#include <algorithm>
+using namespace std;
 
 #include "warn.hh"
-#include "grob.hh"
 #include "staff-symbol-referencer.hh"
-#include "beam.hh"
 #include "stem.hh"
-#include "paper-def.hh"
-#include "group-interface.hh"
+#include "output-def.hh"
+#include "pointer-group-interface.hh"
 #include "align-interface.hh"
 
-const int INTER_QUANT_PENALTY = 1000;
-const Real SECONDARY_BEAM_DEMERIT  = 10.0;
-const int STEM_LENGTH_DEMERIT_FACTOR = 5;
-
-// possibly ridiculous, but too short stems just won't do
-const int STEM_LENGTH_LIMIT_PENALTY = 5000;
-const int DAMPING_DIRECTION_PENALTY = 800;
-const int MUSICAL_DIRECTION_FACTOR = 400;
-const int IDEAL_SLOPE_FACTOR = 10;
+Real
+get_detail (SCM alist, SCM sym, Real def)
+{
+  SCM entry = scm_assq (sym, alist);
 
-const Real ROUND_TO_ZERO_SLOPE = 0.05;
-const int ROUND_TO_ZERO_POINTS = 4;
+  if (scm_is_pair (entry))
+    return robust_scm2double (scm_cdr (entry), def);
+  return def;
+}
 
-extern bool debug_beam_quanting_flag;
+void
+Beam_quant_parameters::fill (Grob *him)
+{
+  SCM details = him->get_property ("details");
+
+  INTER_QUANT_PENALTY = get_detail (details, ly_symbol2scm ("inter-quant-penalty"), 1000.0);
+  SECONDARY_BEAM_DEMERIT = get_detail (details, ly_symbol2scm ("secondary-beam-demerit"), 10.0);
+  STEM_LENGTH_DEMERIT_FACTOR = get_detail (details, ly_symbol2scm ("stem-length-demerit-factor"), 5);
+  REGION_SIZE = get_detail (details, ly_symbol2scm ("region-size"), 2);
+  BEAM_EPS = get_detail (details, ly_symbol2scm ("beam-eps"), 1e-3);
+
+  // possibly ridiculous, but too short stems just won't do
+  STEM_LENGTH_LIMIT_PENALTY = get_detail (details, ly_symbol2scm ("stem-length-limit-penalty"), 5000);
+  DAMPING_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("damping-direction-penalty"), 800);
+  MUSICAL_DIRECTION_FACTOR = get_detail (details, ly_symbol2scm ("musical-direction-factor"), 400);
+  IDEAL_SLOPE_FACTOR = get_detail (details, ly_symbol2scm ("ideal-slope-factor"), 10);
+  ROUND_TO_ZERO_SLOPE = get_detail (details, ly_symbol2scm ("round-to-zero-slope"), 0.02);
+}
 
 static Real
 shrink_extra_weight (Real x, Real fac)
@@ -42,7 +54,6 @@ shrink_extra_weight (Real x, Real fac)
   return fabs (x) * ((x < 0) ? fac : 1.0);
 }
 
-
 struct Quant_score
 {
   Real yl;
@@ -54,20 +65,19 @@ struct Quant_score
 #endif
 };
 
-
 /*
   TODO:
-  
-   - Make all demerits customisable
 
-   - One sensible check per demerit (what's this --hwn)
+  - Make all demerits customisable
 
-   - Add demerits for quants per se, as to forbid a specific quant
-     entirely
+  - One sensible check per demerit (what's this --hwn)
 
+  - Add demerits for quants per se, as to forbid a specific quant
+  entirely
 */
 
-int best_quant_score_idx (Array<Quant_score>  const & qscores)
+int
+best_quant_score_idx (Array<Quant_score> const &qscores)
 {
   Real best = 1e6;
   int best_idx = -1;
@@ -75,39 +85,43 @@ int best_quant_score_idx (Array<Quant_score>  const & qscores)
     {
       if (qscores[i].demerits < best)
        {
-         best = qscores [i].demerits ;
+         best = qscores [i].demerits;
          best_idx = i;
        }
     }
 
   return best_idx;
 }
-  
-MAKE_SCHEME_CALLBACK (Beam, quanting, 1);
+
+MAKE_SCHEME_CALLBACK (Beam, quanting, 2);
 SCM
-Beam::quanting (SCM smob)
+Beam::quanting (SCM smob, SCM posns)
 {
   Grob *me = unsmob_grob (smob);
 
-  SCM s = me->get_grob_property ("positions");
-  Real yl = gh_scm2double (gh_car (s));
-  Real yr = gh_scm2double (gh_cdr (s));
+  Beam_quant_parameters parameters;
+  parameters.fill (me);
 
-  Real ss = Staff_symbol_referencer::staff_space (me);
-  Real thickness = gh_scm2double (me->get_grob_property ("thickness")) / ss;
-  Real slt = me->get_paper ()->get_realvar (ly_symbol2scm ("linethickness")) / ss;
+  Real yl = scm_to_double (scm_car (posns));
+  Real yr = scm_to_double (scm_cdr (posns));
 
+  /*
+    Calculations are relative to a unit-scaled staff, i.e. the quants are
+    divided by the current staff_space.
 
-  SCM sdy = me->get_grob_property ("least-squares-dy");
-  Real dy_mus = gh_number_p (sdy) ? gh_scm2double (sdy) : 0.0;
-  
+  */
+  Real ss = Staff_symbol_referencer::staff_space (me);
+  Real thickness = Beam::get_thickness (me) / ss;
+  Real slt = Staff_symbol_referencer::line_thickness (me) / ss;
+
+  Real dy_mus = robust_scm2double (me->get_property ("least-squares-dy"), 0);
   Real straddle = 0.0;
   Real sit = (thickness - slt) / 2;
   Real inter = 0.5;
   Real hang = 1.0 - (thickness - slt) / 2;
   Real quants [] = {straddle, sit, inter, hang };
-  
-  int num_quants = int (sizeof (quants)/sizeof (Real));
+
+  int num_quants = int (sizeof (quants) / sizeof (Real));
   Array<Real> quantsl;
   Array<Real> quantsr;
 
@@ -116,28 +130,26 @@ Beam::quanting (SCM smob)
     wtk1-fugue2.
 
 
-    (result indexes between 70 and 575)  ? --hwn. 
+    (result indexes between 70 and 575)  ? --hwn.
 
   */
 
-
-  
   /*
     Do stem computations.  These depend on YL and YR linearly, so we can
     precompute for every stem 2 factors.
-   */
-  Link_array<Grob> stems=
-    Pointer_group_interface__extract_grobs (me, (Grob*)0, "stems");
+  */
+  Link_array<Grob> stems
+    = extract_grob_array (me, "stems");
   Array<Stem_info> stem_infos;
   Array<Real> base_lengths;
-  Array<Real> stem_xposns;  
+  Array<Real> stem_xposns;
 
-  Drul_array<bool> dirs_found(0,0);
+  Drul_array<bool> dirs_found (0, 0);
   Grob *common[2];
   for (int a = 2; a--;)
-    common[a] = common_refpoint_of_array (stems, me, Axis(a));
+    common[a] = common_refpoint_of_array (stems, me, Axis (a));
 
-  Grob * fvs = first_visible_stem (me);
+  Grob *fvs = first_visible_stem (me);
   Grob *lvs = last_visible_stem (me);
   Real xl = fvs ? fvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
   Real xr = fvs ? lvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
@@ -148,45 +160,48 @@ Beam::quanting (SCM smob)
     Sometimes my head is screwed on backwards.  The stemlength are
     AFFINE linear in YL and YR. If YL == YR == 0, then we might have
     stem_y != 0.0, when we're cross staff.
-    
-   */
-  for (int i= 0; i < stems.size(); i++)
+
+  */
+  for (int i = 0; i < stems.size (); i++)
     {
-      Grob*s = stems[i];
-      stem_infos.push (Stem::get_stem_info (s));
+      Grob *s = stems[i];
+
+      Stem_info si (Stem::get_stem_info (s));
+      si.scale (1 / ss);
+      stem_infos.push (si);
       dirs_found[stem_infos.top ().dir_] = true;
 
-      bool f = to_boolean (s->get_grob_property ("french-beaming"))
-        && s != lvs && s != fvs;
+      bool f = to_boolean (s->get_property ("french-beaming"))
+       && s != lvs && s != fvs;
 
       base_lengths.push (calc_stem_y (me, s, common, xl, xr,
-                                     Interval (0,0), f));
+                                     Interval (0, 0), f) / ss);
       stem_xposns.push (s->relative_coordinate (common[X_AXIS], X_AXIS));
     }
 
-  bool xstaff= false;
+  bool xstaff = false;
   if (lvs && fvs)
     {
       Grob *commony = fvs->common_refpoint (lvs, Y_AXIS);
       xstaff = Align_interface::has_interface (commony);
     }
-  
+
   Direction ldir = Direction (stem_infos[0].dir_);
   Direction rdir = Direction (stem_infos.top ().dir_);
-  bool knee_b = dirs_found[LEFT] && dirs_found[RIGHT];
+  bool is_knee = dirs_found[LEFT] && dirs_found[RIGHT];
 
+  int region_size = (int) parameters.REGION_SIZE;
 
-  int region_size = REGION_SIZE;
   /*
-    Knees are harder, lets try some more possibilities for knees. 
-   */
-  if (knee_b)
+    Knees are harder, lets try some more possibilities for knees.
+  */
+  if (is_knee)
     region_size += 2;
 
   /*
     Asymetry ? should run to <= region_size ?
-   */
-  for (int i = -region_size ; i < region_size; i++)
+  */
+  for (int i = -region_size; i < region_size; i++)
     for (int j = 0; j < num_quants; j++)
       {
        quantsl.push (i + quants[j] + int (yl));
@@ -194,15 +209,15 @@ Beam::quanting (SCM smob)
       }
 
   Array<Quant_score> qscores;
-  
-  for (int l =0; l < quantsl.size (); l++)  
-    for (int r =0; r < quantsr.size (); r++)
+
+  for (int l = 0; l < quantsl.size (); l++)
+    for (int r = 0; r < quantsr.size (); r++)
       {
        Quant_score qs;
        qs.yl = quantsl[l];
        qs.yr = quantsr[r];
        qs.demerits = 0.0;
-       
+
        qscores.push (qs);
       }
 
@@ -212,28 +227,31 @@ Beam::quanting (SCM smob)
      parameters outside of the loop, we can save a lot of time. */
   for (int i = qscores.size (); i--;)
     {
-      Real d =  score_slopes_dy (qscores[i].yl, qscores[i].yr,
-                                dy_mus, yr- yl, 
-                                xr - xl,
-                                xstaff);
+      Real d = score_slopes_dy (qscores[i].yl, qscores[i].yr,
+                               dy_mus, yr- yl,
+                               xr - xl,
+                               xstaff, &parameters);
       qscores[i].demerits += d;
 
 #if DEBUG_QUANTING
-      qscores[i].score_card_ += to_string ("S%.2f",d);
+      qscores[i].score_card_ += to_string ("S%.2f", d);
 #endif
     }
 
   Real rad = Staff_symbol_referencer::staff_radius (me);
-  int beam_count = get_beam_count (me);
-  Real beam_translation = get_beam_translation (me);
+  Drul_array<int> edge_beam_counts
+    (Stem::beam_multiplicity (stems[0]).length () + 1,
+     Stem::beam_multiplicity (stems.top ()).length () + 1);
+
+  Real beam_translation = get_beam_translation (me) / ss;
 
-  Real reasonable_score = (knee_b) ? 200000 : 100;
+  Real reasonable_score = (is_knee) ? 200000 : 100;
   for (int i = qscores.size (); i--;)
     if (qscores[i].demerits < reasonable_score)
       {
        Real d = score_forbidden_quants (qscores[i].yl, qscores[i].yr,
-                                    rad, slt, thickness, beam_translation,
-                                    beam_count, ldir, rdir); 
+                                        rad, slt, thickness, beam_translation,
+                                        edge_beam_counts, ldir, rdir, &parameters);
        qscores[i].demerits += d;
 
 #if DEBUG_QUANTING
@@ -244,12 +262,12 @@ Beam::quanting (SCM smob)
   for (int i = qscores.size (); i--;)
     if (qscores[i].demerits < reasonable_score)
       {
-       Real d=score_stem_lengths (stems, stem_infos,
-                                base_lengths, stem_xposns,
-                                xl, xr,
-                                knee_b,
-                                qscores[i].yl, qscores[i].yr);
-       qscores[i].demerits +=  d;
+       Real d = score_stem_lengths (stems, stem_infos,
+                                    base_lengths, stem_xposns,
+                                    xl, xr,
+                                    is_knee,
+                                    qscores[i].yl, qscores[i].yr, &parameters);
+       qscores[i].demerits += d;
 
 #if DEBUG_QUANTING
        qscores[i].score_card_ += to_string (" L %.2f", d);
@@ -258,47 +276,55 @@ Beam::quanting (SCM smob)
 
   int best_idx = best_quant_score_idx (qscores);
 
-
 #if DEBUG_QUANTING
-  SCM inspect_quants = me->get_grob_property ("inspect-quants");
-  if (debug_beam_quanting_flag
-      && gh_pair_p (inspect_quants))
+  SCM inspect_quants = me->get_property ("inspect-quants");
+  if (to_boolean (me->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-quanting")))
+      && scm_is_pair (inspect_quants))
     {
-      Real il = gh_scm2double (gh_car (inspect_quants));
-      Real ir = gh_scm2double (gh_cdr (inspect_quants));
+      Drul_array<Real> ins = ly_scm2interval (inspect_quants);
 
       int i = 0;
 
       Real mindist = 1e6;
-      for (; i < qscores.size(); i ++)
+      for (; i < qscores.size (); i++)
        {
-         Real d =fabs (qscores[i].yl-il) + fabs (qscores[i].yr - ir);
+         Real d = fabs (qscores[i].yl- ins[LEFT]) + fabs (qscores[i].yr - ins[RIGHT]);
          if (d < mindist)
            {
              best_idx = i;
-             mindist= d;
+             mindist = d;
            }
        }
       if (mindist > 1e5)
-       programming_error ("Could not find quant.");
+       programming_error ("can't find quant");
     }
 #endif
+
+  Interval final_positions;
+  if (best_idx < 0)
+    {
+      warning (_ ("no feasible beam position"));
+      final_positions = Interval (0, 0);
+    }
+  else
+    {
+      final_positions = Drul_array<Real> (qscores[best_idx].yl,
+                                         qscores[best_idx].yr);
+    }
   
-  me->set_grob_property ("positions",
-                        gh_cons (gh_double2scm (qscores[best_idx].yl),
-                                 gh_double2scm (qscores[best_idx].yr)));
 #if DEBUG_QUANTING
-  if (debug_beam_quanting_flag)
+  if (best_idx >= 0
+      && to_boolean (me->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-quanting"))))
     {
       qscores[best_idx].score_card_ += to_string ("i%d", best_idx);
-      
+
       // debug quanting
-      me->set_grob_property ("quant-score",
-                            scm_makfrom0str (qscores[best_idx].score_card_.to_str0 ()));
+      me->set_property ("quant-score",
+                       scm_makfrom0str (qscores[best_idx].score_card_.to_str0 ()));
     }
 #endif
 
-  return SCM_UNSPECIFIED;
+  return ly_interval2scm (final_positions);
 }
 
 Real
@@ -306,60 +332,63 @@ Beam::score_stem_lengths (Link_array<Grob> const &stems,
                          Array<Stem_info> const &stem_infos,
                          Array<Real> const &base_stem_ys,
                          Array<Real> const &stem_xs,
-                         Real xl, Real xr, 
-                         bool knee, 
-                         Real yl, Real yr)
+                         Real xl, Real xr,
+                         bool knee,
+                         Real yl, Real yr,
+
+                         Beam_quant_parameters const *parameters)
 {
-  Real limit_penalty = STEM_LENGTH_LIMIT_PENALTY;
+  Real limit_penalty = parameters->STEM_LENGTH_LIMIT_PENALTY;
   Drul_array<Real> score (0, 0);
   Drul_array<int> count (0, 0);
-  
-  for (int i=0; i < stems.size (); i++)
+
+  for (int i = 0; i < stems.size (); i++)
     {
-      Grobs = stems[i];
-      if (Stem::invisible_b (s))
+      Grob *s = stems[i];
+      if (Stem::is_invisible (s))
        continue;
 
       Real x = stem_xs[i];
-      Real dx = xr-xl;
-      Real beam_y = dx ? yr *(x - xl)/dx + yl * ( xr - x)/dx : (yr + yl)/2;
+      Real dx = xr - xl;
+      Real beam_y = dx ? yr * (x - xl) / dx + yl * (xr - x) / dx : (yr + yl) / 2;
       Real current_y = beam_y + base_stem_ys[i];
-      Real length_pen = STEM_LENGTH_DEMERIT_FACTOR;
-      
+      Real length_pen = parameters->STEM_LENGTH_DEMERIT_FACTOR;
+
       Stem_info info = stem_infos[i];
       Direction d = info.dir_;
 
-      score[d] += limit_penalty * (0 >? (d * (info.shortest_y_ - current_y)));
-      
+      score[d] += limit_penalty * max (0.0, (d * (info.shortest_y_ - current_y)));
+
       Real ideal_diff = d * (current_y - info.ideal_y_);
       Real ideal_score = shrink_extra_weight (ideal_diff, 1.5);
-      
+
       /* We introduce a power, to make the scoring strictly
          convex. Otherwise a symmetric knee beam (up/down/up/down)
          does not have an optimum in the middle. */
       if (knee)
        ideal_score = pow (ideal_score, 1.1);
-      
+
       score[d] += length_pen * ideal_score;
 
-      count[d] ++;
+      count[d]++;
     }
 
   Direction d = DOWN;
   do
-    { 
-      score[d] /= (count[d] >? 1);
-    }
-  while (flip (&d) != DOWN);
+    score[d] /= max (count[d], 1);
+  while (flip (&d) != DOWN)
+    ;
 
-  return score[LEFT]+score[RIGHT];
+  return score[LEFT] + score[RIGHT];
 }
 
 Real
 Beam::score_slopes_dy (Real yl, Real yr,
                       Real dy_mus, Real dy_damp,
                       Real dx,
-                      bool xstaff)
+                      bool xstaff,
+
+                      Beam_quant_parameters const *parameters)
 {
   Real dy = yr - yl;
   Real dem = 0.0;
@@ -370,35 +399,26 @@ Beam::score_slopes_dy (Real yl, Real yr,
 
     TODO: find a way to incorporate the complexity of the beam in this
     penalty.
-   */
-  if (fabs (dy/dx) > ROUND_TO_ZERO_SLOPE
+  */
+  if (fabs (dy / dx) > parameters->ROUND_TO_ZERO_SLOPE
       && sign (dy_damp) != sign (dy))
-    {
-      dem += DAMPING_DIRECTION_PENALTY;
-    }
-
-   dem += MUSICAL_DIRECTION_FACTOR * (0 >? (fabs (dy) - fabs (dy_mus)));
+    dem += parameters->DAMPING_DIRECTION_PENALTY;
 
+  dem += parameters->MUSICAL_DIRECTION_FACTOR
+    * max (0.0, (fabs (dy) - fabs (dy_mus)));
 
-   Real slope_penalty = IDEAL_SLOPE_FACTOR;
+  Real slope_penalty = parameters->IDEAL_SLOPE_FACTOR;
 
-   /* Xstaff beams tend to use extreme slopes to get short stems. We
-      put in a penalty here. */
-   if (xstaff)
-     slope_penalty *= 10;
+  /* Xstaff beams tend to use extreme slopes to get short stems. We
+     put in a penalty here. */
+  if (xstaff)
+    slope_penalty *= 10;
 
-   /* Huh, why would a too steep beam be better than a too flat one ? */
-   dem += shrink_extra_weight (fabs (dy_damp) - fabs (dy), 1.5)
-     * slope_penalty;
+  /* Huh, why would a too steep beam be better than a too flat one ? */
+  dem += shrink_extra_weight (fabs (dy_damp) - fabs (dy), 1.5)
+    * slope_penalty;
 
-   /*
-     almost zero slopes look like errors in horizontal beams. 
-    */
-   if (fabs (dy) > 1e-3
-       && fabs (dy / dx) < ROUND_TO_ZERO_SLOPE)
-     dem += ROUND_TO_ZERO_POINTS;
-   
-   return dem;
+  return dem;
 }
 
 static Real
@@ -407,124 +427,106 @@ my_modf (Real x)
   return x - floor (x);
 }
 
-
 /*
   TODO: The fixed value SECONDARY_BEAM_DEMERIT is probably flawed:
-   because for 32nd and 64th beams the forbidden quants are relatively
-   more important than stem lengths.
+  because for 32nd and 64th beams the forbidden quants are relatively
+  more important than stem lengths.
 */
 Real
 Beam::score_forbidden_quants (Real yl, Real yr,
                              Real radius,
                              Real slt,
                              Real thickness, Real beam_translation,
-                             int beam_count,
-                             Direction ldir, Direction rdir)
+                             Drul_array<int> beam_counts,
+                             Direction ldir, Direction rdir,
+
+                             Beam_quant_parameters const *parameters)
 {
   Real dy = yr - yl;
-  Drul_array<Real> y(yl,yr);
-  Drul_array<Direction> dirs(ldir,rdir);
-  
-  Real extra_demerit = SECONDARY_BEAM_DEMERIT / beam_count;
+  Drul_array<Real> y (yl, yr);
+  Drul_array<Direction> dirs (ldir, rdir);
 
-  /*
-    Inside the staff, inter quants are forbidden.
-   */
-  Real dem = 0.0;
-  Direction d = LEFT;
-  do
-    {
-      if (fabs (y[d]) <= (radius + 0.5) && fabs (my_modf (y[d]) - 0.5) < 1e-3)
-       dem += INTER_QUANT_PENALTY;
-    }
-  while ((flip (&d))!= LEFT); 
+  Real extra_demerit = parameters->SECONDARY_BEAM_DEMERIT / (max (beam_counts[LEFT], beam_counts[RIGHT]));
 
+  Direction d = LEFT;
+  Real dem = 0.0;
+  Real eps = parameters->BEAM_EPS;
 
-  for (int j = 1; j <= beam_count; j++)
+  do
     {
-      do
+      for (int j = 1; j <= beam_counts[d]; j++)
        {
-         /*
-           see if the outer staffline falls in a beam-gap
-           
-           This test is too weak; we should really check all lines.
-          */
          Direction stem_dir = dirs[d];
-         Real gap1 =  y[d] - stem_dir * ((j-1) * beam_translation + thickness / 2 - slt/2 );
-         Real gap2 = y[d] - stem_dir * (j * beam_translation - thickness / 2 + slt/2);
+
+         /*
+           The 2.2 factor is to provide a little leniency for
+           borderline cases. If we do 2.0, then the upper outer line
+           will be in the gap of the (2, sit) quant, leading to a
+           false demerit.
+         */
+         Real gap1 = y[d] - stem_dir * ((j - 1) * beam_translation + thickness / 2 - slt / 2.2);
+         Real gap2 = y[d] - stem_dir * (j * beam_translation - thickness / 2 + slt / 2.2);
 
          Interval gap;
          gap.add_point (gap1);
          gap.add_point (gap2);
 
-         if (gap.elem_b (radius))
-           dem += extra_demerit;
+         for (Real k = -radius;
+              k <= radius + eps; k += 1.0)
+           if (gap.contains (k))
+             {
+               Real dist = min (fabs (gap[UP] - k), fabs (gap[DOWN] - k));
+
+               /*
+                 this parameter is tuned to grace-stem-length.ly
+               */
+               Real fixed_demerit = 0.4;
+
+               dem += extra_demerit
+                 * (fixed_demerit
+                    + (1 - fixed_demerit) * (dist / gap.length ()) * 2);
+             }
        }
-      while ((flip (&d))!= LEFT); 
     }
+  while ((flip (&d)) != LEFT);
 
-
-  
-  // todo: use beam_count of outer stems.
-  if (beam_count >= 2)
+  if (max (beam_counts[LEFT], beam_counts[RIGHT]) >= 2)
     {
       Real straddle = 0.0;
       Real sit = (thickness - slt) / 2;
       Real inter = 0.5;
       Real hang = 1.0 - (thickness - slt) / 2;
 
-      Real eps = 1e-3;
-
-      // hmm, without Interval/Drul_array, you get ~ 4x same code...
-      if (fabs (y[LEFT] - dirs[LEFT] * beam_translation) < radius + inter)
-       {
-         if (dirs[LEFT] == UP && dy <= eps
-             && fabs (my_modf (y[LEFT]) - sit) < eps)
-           dem += extra_demerit;
-         
-         if (dirs[LEFT] == DOWN && dy >= eps
-             && fabs (my_modf (y[LEFT]) - hang) < eps)
-           dem += extra_demerit;
-       }
-
-      if (fabs (y[RIGHT] - dirs[RIGHT] * beam_translation) < radius + inter)
-       {
-         if (dirs[RIGHT] == UP && dy >= eps
-             && fabs (my_modf (y[RIGHT]) - sit) < eps)
-           dem += extra_demerit;
-         
-         if (dirs[RIGHT] == DOWN && dy <= eps
-             && fabs (my_modf (y[RIGHT]) - hang) < eps)
-           dem += extra_demerit;
-       }
-      
-      if (beam_count >= 3)
+      Direction d = LEFT;
+      do
        {
-         if (fabs (y[LEFT] - 2 * dirs[LEFT] * beam_translation) < radius + inter)
+         if (beam_counts[d] >= 2
+             && fabs (y[d] - dirs[d] * beam_translation) < radius + inter)
            {
-             if (dirs[LEFT] == UP && dy <= eps
-                 && fabs (my_modf (y[LEFT]) - straddle) < eps)
+             if (dirs[d] == UP && dy <= eps
+                 && fabs (my_modf (y[d]) - sit) < eps)
                dem += extra_demerit;
-             
-             if (dirs[LEFT] == DOWN && dy >= eps
-                 && fabs (my_modf (y[LEFT]) - straddle) < eps)
+
+             if (dirs[d] == DOWN && dy >= eps
+                 && fabs (my_modf (y[d]) - hang) < eps)
                dem += extra_demerit;
            }
-         
-         if (fabs (y[RIGHT] - 2 * dirs[RIGHT] * beam_translation) < radius + inter)
+
+         if (beam_counts[d] >= 3
+             && fabs (y[d] - 2 * dirs[d] * beam_translation) < radius + inter)
            {
-             if (dirs[RIGHT] == UP && dy >= eps
-                 && fabs (my_modf (y[RIGHT]) - straddle) < eps)
+             if (dirs[d] == UP && dy <= eps
+                 && fabs (my_modf (y[d]) - straddle) < eps)
                dem += extra_demerit;
-             
-             if (dirs[RIGHT] == DOWN && dy <= eps
-                 && fabs (my_modf (y[RIGHT]) - straddle) < eps)
+
+             if (dirs[d] == DOWN && dy >= eps
+                 && fabs (my_modf (y[d]) - straddle) < eps)
                dem += extra_demerit;
            }
        }
+      while (flip (&d) != LEFT);
     }
-  
+
   return dem;
 }
 
-