[metapost] glyph operator and contours order

[Dan Luecking]

At 09:37 AM 2/6/2011, Laurent Méhats wrote:

>The 'even-odd rule' wikipedia page led me to the 'point in polygon
>problem' and then to the 'winding number algorithm': "the winding number
>of a closed curve in the plane around a given point is an integer
>representing the total number of times that curve travels counterclockwise
>around the point"; if the winding number is zero then the point lies
>outside the curve, else it lies inside.
>
>vardef wnum (expr ccl, pnt)=
>   save sum;
>   sum:=0;
>   for t=1 upto length ccl:
>     sum:=sum+
>       angle((point t of ccl-pnt) rotated -angle(point t-1 of ccl-pnt));
>   endfor
>   sum/360
>enddef;

This doesn't work except for polygons. The above
calculation only deals with the endpoints of a segment,
while the segment itself could go around the point on
either side.

One modification can make it work: instead of evaluating
only at (point t of ccl) for integer t, step through the
loop with fractional step. How small a fraction? Small enough
so that the length of arc traveled between successive values
of t is less than the the distance from the path ccl to pnt.
This would either require a calculation to get that distance,
or one could vary the step size depending on the current
distance:

% seg is a single Bezier segment (length=1)
vardef delta_theta (expr seg, pnt)=
   save A,B,C,D, M;
   A := point 0 of seg;
   B := postcontrol 0 of seg;
   C := precontrol 1 of seg;
   D := point 1 of seg;
   % M = upper bound on derivative of seg.
   M := 3*max(abs(B-A),abs(C-B),abs(D-C));
   if M=0: 0 % trivial path, return 0
   else:
     save sum; sum := 0;
     save curr_t,next_t; curr_t := 0;
     M := M + 1; % to avoid overflow when dividing.
     forever:
       % precision problem: next_t = curr_t is possible
       % if numerator is very small and M rather large.
       next_t := curr_t + abs((point curr_t of seg) - pnt)/M;
       next_t := max(next_t,curr_t + epsilon);
       next_t := min(1,next_t); % don't overshoot end of seg.
       sum := sum +
         angle(((point next_t of seg) - pnt)
                rotated -angle((point curr_t of seg) - pnt));
     exitif next_t = 1;
     endfor
     sum % return sum
   fi
enddef;

Then get the winding angle by summing
     delta_theta (subpath (j-1,j) of ccl, pnt).
(Warning: I did not test this code.)

Another approach is to simply keep track of the
quadrants the cycle moves through: count +/-0.25
for each such transition. This can be detected, without
angle computations, by noting a change in sign of
xpart or ypart. The size of M ensures (at least
mathematically) that only one of these can change sign
at once. One subtle point: one has to decide which quadrant
each semi-axis belongs to.

>Except for precision problems, this should do the trick

One only needs accuracy within +/-180 degrees to
distinguish inside from outside. Unfortunately, angle
and rotation precision are inversely proportional to the
length of the vectors, so a complete analysis of the
precision would have to take into account the location of
a point in relation to the path.

The axis crossing algorithm has no unusual precision
problems. Apart, that is, from such problems inherent in
metapost. For example, if pnt is close enough to the cycle
(say within epsilon), then a change in t of epsilon (smallest
possible) could jump multiple quadrants.


Dan


Daniel H. Luecking
Department of Mathematical Sciences
Fayetteville, Arkansas
http://www-cs-faculty.stanford.edu/~knuth/iaq.html