[metapost] new is_clockwise routine

[Werner LEMBERG]

Dear hobby drawers,


a longer time ago we discussed how to circumvent the turning number
bug in metapost to (re)implement the functions `clockwise' and
`counterclockwise' in mf2pt1.  We finally came up with the following
solution:


----------------------------------------------------------------------


vardef Angle primary d =
  if d <> (0, 0):
    angle d
  else:
    0
  fi
enddef;

vardef postdir expr t of p =    % t is expected to be an integer.
  save a, b, c, d, s;
  pair a, b, c, d;
  path s;

  s := subpath (t, t + 1) of p;

  a = point 0 of s;
  b = postcontrol 0 of s;
  c = precontrol 1 of s;
  d = point 1 of s;

  if a <> b:
    b - a
  elseif a <> c:
    c - a
  elseif a <> d:
    d - a
  else:
    (0, 0)
  fi
enddef;

vardef is_clockwise primary p =
  save res;

  res := 0;

  for t = 0 upto length p - 1:
    res := res
	   + (-Angle (postdir t of p)
	      + Angle (postdir t + 1 of p) + 180) mod 360
	   - 180;
  endfor;

  res <= 0
enddef;

vardef counterclockwise primary c =
  (if is_clockwise c: (reverse c) else: c fi)
enddef;

vardef clockwise primary c =
  (if is_clockwise c: c else: (reverse c) fi)
enddef;


----------------------------------------------------------------------


Unfortunately, this doesn't work, as I've found out meanwhile.  As an
example, it fails for this curve:

  (0,0)
  .. controls (10,-30) and (10,40).. (0,10)
  -- cycle;

Reason is that the above algorithm can't recognize properly that a
Bézier curve element in the path turns more than 180 degrees.  I tried
hard, but I couldn't find a simple solution to overcome this flaw
reliably.

Below I present a completely different solution which should work for
all curves which don't have a self-intersection.  It simply searches
the outermost intersection point of the path with an infinitely long
horizontal line and determines its direction.

Please comment.  It seems to run reasonably fast (at least I couldn't
find a greate difference to the old algorithm), but maybe it can be
streamlined.


    Werner


======================================================================


vardef crossproduct (expr u, v) =
  save u_, v_;
  pair u_, v_;

  u_ := unitvector u;
  v_ := unitvector v;

  abs (xpart u_ * ypart v_ - ypart u_ * xpart v_)
enddef;

vardef makeline primary p =
  save start, loop, d, n;
  pair start, d;

  loop := 0;
  for i = 0.5 step 1 until length p - 0.5:
    % add some randomness to get different lines for each function call
    n = uniformdeviate 0.9 - 0.45 + i;
    start := point n of p;
    d := direction n of p;
    if d <> (0, 0):
      loop := 1;
    fi;
    exitif loop = 1;
  endfor;

  if loop = 0:
    errmessage ("Cannot find a starting point on path for orientation test");
  fi;

  d := unitvector (d rotated 90);

  % construct line orthogonal to tangent which intersects path at least once
  start - eps * d
  -- infinity * d
enddef;

vardef is_clockwise primary p =
  save line, cut, cut_new, res, line_dir, tangent_dir;
  path line;
  pair cut, cut_new, line_dir, tangent_dir;

  res := 1;
  line := makeline p;
  line_dir := direction 0 of line;

  % find outermost intersection
  cut := (0, 0);
  forever:
    cut_new := line intersectiontimes p;
    exitif cut_new = (-1, -1);

    % compute new line if we have a strange intersection
    tangent_dir := direction (ypart cut_new) of p;

    if abs tangent_dir < eps:
      % vector zero or too small
      line := makeline p;
      line_dir := direction 0 of line;
    elseif crossproduct (tangent_dir, line_dir) < 0.02:
      % grazing intersection
      line := makeline p;
      line_dir := direction 0 of line;
    else:
      cut := cut_new;
      line := subpath (xpart cut + eps, infinity) of line;
    fi;
  endfor;

  tangent_dir := direction (ypart cut) of p;
  res := (angle tangent_dir - angle line_dir + 180) mod 360 - 180;

  res < 0
enddef;

vardef counterclockwise primary c =
  (if is_clockwise c: (reverse c) else: c fi)
enddef;

vardef clockwise primary c =
  (if is_clockwise c: c else: (reverse c) fi)
enddef;