def inside(expr j, p) =
  begingroup
  save n;
  hide(
  save L, P, t, tt,bottom, top;
  boolean abort;
  path L,P;
  pair bottom, top;
  bottom = llcorner p;
  top = urcorner p;
  P := p;
  n := 0;
  if (xpart j <= xpart bottom) or (xpart j >= xpart top) 
    or (ypart j <= ypart bottom) or (ypart j >= ypart top):
  else:
    P := subpath (epsilon,(length P)) of P;
    L := j--(xpart j, -infinity);
    forever:
      t := xpart (P intersectiontimes L);
      exitif t<0;
      tt := t;
      forever:
        exitif t<0;
        exitif not (xpart direction t of P = 0);
 	L := L shifted (-epsilon,0); 
        t := xpart (P intersectiontimes L);
        exitif t<0;
        if t<tt: 
   	  L := L shifted (3epsilon,0); 
          t := xpart (P intersectiontimes L);
        fi
      endfor
      exitif t < 0;
      n := n + 1;
      P := subpath (t+epsilon,(length P)) of P;
    endfor; 
  fi )
if (n>0) and (odd n): true else: false fi
endgroup
enddef;

beginfig(1);
path P;
 P:= (fullcircle scaled 20) shifted (20,20);
 P:= (0,10){down}..(10,10){up} .. (30,10){down} .. (20,10){up}..cycle;
 P := ((0,0)--(10,0)--(10,10)--(0,10)--cycle) rotated 130 shifted (20,10);
% P := (10,10){down}..(30,10){up} .. (0,10){down} .. (20,10){up}..cycle;
% P := ((25,10){down}..(10,0){right}..(30,5){up}..
%   (20,10){up}..(30,15){up}..(20,20){left}..cycle) shifted (0,5);

nx := 50;
ny := 40;
pair endgrid;
endgrid := (40,30);
pickup pencircle xscaled (xpart endgrid/nx) yscaled (ypart endgrid/ny);
pair J;

for k = 1 upto nx:
 for l = 1 upto ny:
  J := (k/nx*xpart endgrid,l/ny*ypart endgrid);
  if inside(J,P):
     drawdot J withcolor green;
  else:
     drawdot J withcolor blue;
  fi
 endfor;
endfor;
draw P withcolor red;
currentpicture := currentpicture scaled 10;

endfig;
end.