6

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;

inside or ouside?

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

  • 1
    One ad hoc poposal: take one point of the path (or better a point inside the contour), connect it through a straight line (say) with the point in question, and then use this method to count the intersections of the straight line with the path. If it is even, the point is inside the contour, otherwise it is outside. – user121799 Mar 13 '18 at 0:21
  • If the point is on path, this might not work. For example: imagine a point on the "north" of the "southern peninsula" on the picture above; the number of intersections will be even for 'r' and odd for 'q'. And to find a point inside is a problem on its own (shouldn't be too difficult though). And even if the same thing is done the other way around (as Marcel Krüger have suggested below), intersections in MP do not seem to be a very reliable tool for such a task. – Sergey Slyusarev Mar 13 '18 at 7:49
4

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.

  • This solution as it is doesn't seem to be very reliable. Take, for example, these path and point: p := (0, 0) -- (0, 10cm) -- (10cm, 10cm) -- (10cm, 0cm) -- cycle; z0 = (5cm, 5cm); 'z0' is clearly inside, but the output is 'OUT' (it seems to have something to do with passing through angles). or here: p := (0, 0) -- (0, 10cm) -- (1cm, 10cm) -- (4cm, 6cm) -- (10cm, 6cm) -- (10cm, 0cm) -- cycle; z0 = (5cm, 7cm); 'z0' is outside, but the output is 'IN'. Plus, it feels that there could be some problems if the cutting line is tangent to something. – Sergey Slyusarev Mar 13 '18 at 7:29
  • 1
    @SergeySlyusarev I added another solution which should be more reliable – Marcel Krüger Mar 13 '18 at 10:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.