Ray Tracing to find a subset of 2 D set of points

In computational geometry, there is a classical problem of partitioning a set of points in a 2D into those enclosed by a polygon and those outside. Currently, the best way is supposed to be a ray tracing method where a ray is drawn from a candidate point in the positive X-direction and number of intersections with the polygon is counted. If the number is odd, the point is inside, if it is even, it is outside.

In Mathematica 8.0 (or probably earlier), there is a new facility to get coordinates of points inside a 2D plot. With a right-mouse click inside the plot, choose “Get Coordinates”, then draw the boundary of the area desired. Points are chosen which form the polygon to divide. Then Ctrl+C copies the coordinates of the points. Assign it to the tlist in the code below. Then run the whole code.

You get a partition of the points inside and outside the polygon. The code is optimized as much as possible so that time taken is a minimum. But more optimization for speed is not ruled out.

alist=RandomReal[{-5,5},{10000,2}];
tlist={{-3.567,3.357},{-3.772,3.357},{-4.079,2.678},{-3.977,2.451},{-3.67,2.282},{-3.567,1.999},{-3.704,1.433},{-3.943,0.6971},{-4.045,-0.03853},{-3.738,-0.2083},{-3.67,-0.5478},{-3.84,-1.17},{-3.84,-1.906},{-3.533,-2.359},{-3.123,-2.019},{-2.953,-1.793},{-2.748,-1.906},{-2.44,-2.698},{-1.962,-2.925},{-1.518,-2.019},{-1.279,-1.057},{-1.04,-0.4912},{-0.8697,0.1312},{-1.279,1.15},{-1.143,1.999},{-1.075,2.847},{-1.348,3.074},{-1.689,2.565},{-2.133,1.489},{-2.372,1.546},{-2.133,2.678},{-2.167,3.13},{-2.474,2.282},{-2.611,2.168},{-2.782,2.395},{-2.645,2.904},{-2.577,3.357},{-2.85,3.583},{-3.021,3.583}};
tl=Append[tlist,tlist[[1]]];
plcomp=Show[ListPlot[tl,PlotStyle->Black,Joined->True],ListPlot[alist,PlotStyle->{PointSize[0.001],Hue[0.01]}],Frame->True,PlotRange->All,Axes->False];
Ptin[tp_,tl_:tl]:=Module[{mlx,mly,upy,crp},
{mlx,mly}=Transpose@((#-tp)&/@tl);
upy=Ratios[mly];
crp=Flatten@Position[upy,_?Negative];
Times@@((mlx[[#]]-mly[[#]](((mlx[[#+1]]-mlx[[#]])/(mly[[#+1]]-mly[[#]]))))&/@crp)<0
];
Timing@Show[plcomp,ListPlot[Select[alist,Ptin],PlotStyle->{PointSize[0.001],Green}]]