[metapost] point "inside" cycle

[Boguslaw Jackowski]

> Does MetaPost have a clever way of determining whether a point z is
> "inside" a closed path (cycle) p?

I believe that the method based on computing a winding number for a given 
point and curve (idea suggested by Larry Siebenmann; described in details
in "Computing the area and winding number for a B\´ezier curve" by B.J.,
http://tug.org/TUGboat/tb33-1/tb103jackowski.pdf ) is sufficiently
efficient and robust for this purpose. Of course, because of general 
problems with discrete geometry, there is an unavoidable problem with 
points nearly touching the curve.

Here you have my implementation that I use in for making fonts:

vardef windingangle(expr p,q) = % |p| -- point, |q| -- B\'ezier segment
  save a,b,v;
  a=length(p-point 0 of q); b=length(p-point 1 of q);
  if min(a,b)<2 eps: % MF and MP are not the masters of precision, we'd 
better stop now
   errhelp "It is rather not advisable to continue. Will return 0.";
   errmessage "windingangle: point unsafely near upon B\'ezier segment 
(dist="
      & decimal(min(a,b)) & ")";
   0
  else:
   v:=mock_arclength(q); % |v| denotes both length and angle
   if (v>=a) and (v>=b): % possibly too long B\'ezier arc
    windingangle(p, subpath (0, 1/2) of q) + windingangle(p, subpath (1/2, 1) of q)
   else:
    v:=angle((point 1 of q)-p)-angle((point 0 of q)-p);
    if v>180: v:=v-360; fi
    if v<-180: v:=v+360; fi
    v
   fi
  fi
enddef;

vardef mock_arclength(expr p) = % |p| -- B\'ezier segment
  % |mock_arclength(p)>=arclength(p)|
  length((postcontrol 0 of p)-(point 0 of p)) +
  length((precontrol 1 of p)-(postcontrol 0 of p)) +
  length((point 1 of p)-(precontrol 1 of p))
enddef;

vardef windingangles (expr p,q) = % |p|, |q| -- cyclic curves
  save a, b, c; a:=b:=0;
  for t:=1 upto length(q):
   a:=a+windingangle(point 0 of p, subpath(t-1,t) of q);
  endfor;
  for t:=1 upto length(p):
   b:=b+windingangle(point 0 of q, subpath(t-1,t) of p);
  endfor;
  (a,b)
enddef;

In particular, I use it to check whether two disjoint curves are
placed inside each other:

tertiarydef a inside b =
  if path a: % |and path b|
   begingroup
    if (a intersectiontimes b)<>(-1,-1):
     false
    else:
     save a_,b_; (a_,b_)=windingangles(a,b); % |a_|, |b_| can be either 0 or $\pm$360
     (abs(abs(a_)-360)<eps) and (abs(b_)<eps)
    fi
   endgroup
  else: % |numeric a and pair b|
   begingroup
    (a>=xpart b) and (a<=ypart b)
   endgroup
  fi
enddef;

Cheers -- Jacko

-- 
BOP s. c.
ul. Bora-Komorowskiego 24, 80-377 Gdansk, Poland
tel. (+48 58) 553 46 59,  fax (+48 58) 511 03 81
bop at bop.com.pl, http://www.bop.com.pl