[metapost] glyph operator and contours order

[Laurent Méhats]

Le 07/02/2011 21:14, Dan Luecking a écrit :
> 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.

Indeed, I was a bit too much enthusiast ...

> 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;

Let ngl stand for angle((D-pnt) rotated -angle(A-pnt)). I understand that
ngl may be incorrect if pnt lies on D--A, and is incorrect if pnt lies
strictly between D--A and seg: in the first case ngl would be 180 while we
may expect -180, and in the second case we would expect ngl+360 if ngl is
negative or ngl-360 if ngl is positive. Is that right ?

[Thinking aloud] The first case is easily detected, due to the value of
ngl. Would it be interesting to try and detect the second one ? (We may
assume that seg has no inflexion point: if B--C intersects seg then it
would be enough to consider seg on both sides of the intersection point.
Thus we may assume that 'seg&(D--A)--cycle' is convex; that may help.)

>   % 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.

Nice.

>> 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.

I'm learning a lot of things here, many thanks.

Regards,
Laurent Méhats