How to check if a given point is inside some closed path in MetaPost

Question

If there is some arbitrary closed path (without self intersections) in MetaPost, is there a simple way to check if a given point falls inside or outside of it?

For example:

path p;
pair q, r;
p := (-3cm, 0) .. (-4cm, 1cm) .. (2cm, 3cm) .. (3cm, 1cm) .. (-1cm, -1cm) .. (1cm, -2cm) .. (-2cm, -1cm) .. cycle;
q := (0, 0);
r := (0, 1cm);
draw p;
fill (fullcircle scaled 1mm) shifted q withcolor red;
fill (fullcircle scaled 1mm) shifted r withcolor green;

In this case I want to get 'false' for 'q', the red one, and 'true' for 'r', the green one regarding their relation to 'p'.

Answer 1

Based on Get all intersection points in MetaPost and implementing the even-odd rule:

  vardef is_inside(expr p,q) =
    begingroup;
    save cut, inside;
    path cut;boolean inside;
    cut := q -- (ulcorner p + (-2,2));
    inside := false;
    forever:
      cut := cut cutbefore p;
      exitif length cuttings = 0;
      cut := subpath(4epsilon, length cut) of cut;
      inside := not inside;
    endfor;
    inside
    endgroup
  enddef;

  beginfig(1) ;
  path p;
  p := (-0.5cm,0cm) .. (5cm,5cm) .. (0cm,3cm) .. (-7cm,9cm) .. cycle;
  z0 = (2cm,3.5cm);
  draw p;
  draw z0 withpen pencircle scaled 2;
  draw origin withpen pencircle scaled 1;

  if is_inside(p,z0):
    label.lrt("IN", origin);
  else:
    label.lrt("OUT", origin);
  fi
  endfig ;

An alternative rule is the non-zero winding rule. This is an implementation of windingnumber for MetaPost by Bogusław Jackowski. I extended the inside operator to allow testing if a point is inside a path:

vardef mock_arclength(expr p) = % |p| -- Bézier 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 windingangle(expr p,q) = % |p| -- point, |q| -- Bézier segment
  save a,b,v;
  a=length(p-point 0 of q); b=length(p-point 1 of q);
  if min(a,b)<2eps: % MP is not the master of precision, we’d better stop now
    errhelp "It is rather not advisable to continue. Will return 0.";
    errmessage "windingangle: point unsafely near Bézier 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ézier 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 windingnumber (expr p,q) = % |p| -- point, |q| -- Bézier spline
    save a; a:=0;
    for t:=1 upto length(q):
        a:=a+windingangle(p, subpath(t-1,t) of q);
    endfor
    a/360
enddef;
tertiarydef a inside b =
    if path a: % |and path b|; |a| and |b| must not touch each other
        begingroup
            save a_,b_; (a_,b_)=
                (windingnumber(point 0 of a,b), windingnumber(point 0 of b,a));
            (abs(a_-1)<eps) and (abs(b_)<eps)
        endgroup
    elseif pair a: % |and path b|; |a| must not lie on |b|
        begingroup
            abs(windingnumber(a,b))>eps
        endgroup
    else: % |numeric a and pair b|
        begingroup
            (a>=xpart b) and (a<=ypart b)
        endgroup
    fi
enddef;

This is based on a topological invariant, so it should be more reliable. At least it should fail if it can't compute the right result.