def inside(expr j, p) =
  begingroup
  save n;
  hide(
  save done, L, P, t, bottom, top;
  boolean done;
  path L,P;
  numeric t[];
  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:
    L := j--(xpart j, ((ypart bottom)-100));
    done := false;
    forever:
     (t[n],whatever) = L intersectiontimes P;
      if (t[n]<0) and (n>0) :  % find upward match
        P := reverse P;
        numeric t[];
        L := j--(xpart j, ((ypart bottom)-100));
        forever:
          (t[n],whatever) = L intersectiontimes P;
          exitif t[n]<0;
          L := subpath (t[n]+epsilon, length L) of L;
          n := n + 1;
        endfor;
        done := true;
      fi
      exitif done;
      L := subpath (t[n]+epsilon,length L) of L;
      n := n + 1;
    endfor; 
  fi )
if (n>1) and (not 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 := unitsquare scaled 10 rotated 30 shifted (20,10);
% P := (10,10){down}..(30,10){up} .. (0,10){down} .. (20,10){up}..cycle;
P := (0,10){down}..(10,0){right}..(30,5){up}..
    (20,10){up}..(30,15){up}..(20,20){left}..cycle;
draw P withcolor red;

nx := 100;
ny := 100;
pair endgrid;
endgrid := (30,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;
currentpicture := currentpicture scaled 10;

endfig;
end.