Here's a potentially better algorithm, as well as a faster implementation of the one you're using: http://geomalgorithms.com/a03-_inclusion.html
I think the optimization is that they check for intersection between a horizontal line and a line segment, which is a LOT faster than checking whether two segments intersect.
EDIT: I ported the code to SB, it's not optimal but should be pretty good
OPTION DEFINT 'crossing number DEF CN_PNPOLY(PX,PY,VX[],VY[]) VAR CN,I FOR I=0 TO LEN(VX)-2 IF (VY[I]<=PY && VY[I+1]>PY) || (VY[I]>PY && VY[I+1]<=PY) THEN VAR VT#=(PY-VY[I])/(VY[I+1]-VY[I]) IF PX<VX[I]+VT#*(VX[I+1]-VX[I]) THEN INC CN ENDIF NEXT RETURN CN AND 1 END 'winding number DEF WN_PNPOLY(PX,PY,VX[],VY[]) VAR WN,I FOR I=0 TO LEN(VX)-2 IF VY[I]<=PY THEN IF VY[I+1]>PY && ISLEFT(VX[I],VY[I],VX[I+1],VY[I+1],PX,PY)>0 THEN INC WN ELSEIF VY[I+1]<=PY && ISLEFT(VX[I],VY[I],VX[I+1],VY[I+1],PX,PY)<0 THEN DEC WN ENDIF NEXT RETURN WN END DEF ISLEFT(X0,Y0,X1,Y1,X2,Y2) RETURN (X1-X0)*(Y2-Y0)-(X2-X0)*(Y1-Y0) ENDRemember that the first point needs to be repeated at the end (so the array for a triangle will actually have a length of 4) The winding number test seems to be slightly faster (as expected) but I only did a few tests. But if you're trying to fill a polygon, there are much more efficient ways than checking whether every single point is inside the shape. Here's one: http://alienryderflex.com/polygon_fill/