pnPoly::usage = " pnPoly[{testx,testy},{{x1,y1},{x2,y2},...}] return true if {testx,testy} in inside polygon. Polygon can be closed or not. A point will be inside exactly one member of a polygonal partitioning. C-code by W. Randolph Franklin (http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html). J. Wesenberg, 2008"; pnPoly[{testx_, testy_}, pts_List] := Xor @@ (( Xor[#[[1, 2]] > testy, #[[2, 2]] > testy] && ((testx - #[[2, 1]]) < (#[[1, 1]] - #[[2, 1]]) (testy - #[[2, 2]])/(#[[1, 2]] - #[[2, 2]])) ) & /@ Partition[pts, 2, 1, {2, 2}])

**A test for the non-believers:**

Needs["ComputationalGeometry`"] pol = (#[[ConvexHull[#]]]) &[RandomReal[{0, 1}, {12, 2}]]; Graphics[{PointSize[Large], {FaceForm[LightGray], EdgeForm[Black], Polygon[pol]}, If[pnPoly[#, pol], {Blue, Point[#]}, {Red, Point[#]}] & /@ RandomReal[{0, 1}, {21, 2}] }] Clear[pol]

**Original.**
The code is ported from this C-version by W.R. Franklin. Note that the Mathematica code is even shorter :)

int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy) { int i, j, c = 0; for (i = 0, j = nvert-1; i < nvert; j = i++) { if ( ((verty[i]>testy) != (verty[j]>testy)) && (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) ) c = !c; } return c; }

## 1 comment:

Thats a really useful one. Is it possible to apply this to 3D shapes like parallelogram ??

Post a Comment