# How to check if a 3D point is in a planar polygon?

Following up on ndroock1's question, I naively tried to apply the solutions to a 3D point and polygon and they didn't work. For example, functions involving ArcTan that are used in kguler's answer don't work with three arguments, Mac's answer doesn't consider the 3rd dimension at all, the undocumented function GraphicsMeshPointWindingNumber doesn't work in 3D, and complex numbers only map to 2D planes. And so on.

So, is there any way to check if a 3D point is in a 3D planar polygon?

-
Why not just apply a linear transformation sending the polygon's plane to the xy plane and proceed from there? –  whuber Sep 7 '12 at 15:07
@whuber I was thinking about that, but I'm not that familiar with FindGeometricTransform. –  Eli Lansey Sep 7 '12 at 15:09
@EliLansey, you could compute the normal n to the plane by taking the Cross product of two vectors in the plane, and then use RotationMatrix[{n,{0,0,1}}] to get a rotation matrix which will rotate everything into the xy plane. –  Simon Woods Sep 7 '12 at 15:22
@Simon, for polygons with inexact coordinates, it'd be safer to use Newell's algorithm to compute the normal... –  Guess who it is. Sep 7 '12 at 15:31
@SimonWoods That's a great idea. Want to write it up as an answer? –  Eli Lansey Sep 7 '12 at 15:32

In Version 10, we can use RegionMember to select points that are within a region; whether we're dealing with points in a 3D polygon or a 2D polygon embedded in 3D. Let's take a look at the latter case which is what the OP asks:

We create a triangle embedded in 3D, discretize it and collect points we know for sure are in the triangle:

tri = Triangle[{{0, 0, 0}, {1, 0, 0}, {0, 1, 1}}];
dt = DiscretizeRegion[tri, MaxCellMeasure -> 0.000001];
intri = MeshCoordinates[dt]; (* points in triangle *)


Now we generate random points and mix them with the points created above:

pts = Join[RandomReal[1, {5000, 3}], intri];


We create a RegionMemberFunction and apply it to all points to determine which points are inside/outside the polygon:

rm = RegionMember[dt]


Notice how this tells us we're dealing with a 2D region embedded in 3D. We apply it to all points:

inout = rm[pts];
colors = inout /. {True -> Red, False -> Blue};


Visualize:

Graphics3D[{Transpose[{colors, Point /@ pts}], EdgeForm[Black], Opacity[0], tri}]


-

Just work in local coordinates for the polygon's plane. This requires finding a change of basis matrix and then applying it (which is just matrix multiplication).

One way to obtain a basis is to use any three points on the polygon, assuming they are not collinear. Thus:

basis[{x_, y_, z_, ___}] := Orthogonalize[{y - x, z - x}] // Transpose;


(Orthogonalize is not really necessary but it assures the transformation preserves distances and angles for applications that might need that.)

Right-multiplication by basis[...] converts 3D coordinates to 2D coordinates (ignoring any component of the 3D coordinate orthogonal to the polygon's plane). Then apply whatever 2D algorithm you prefer.

### Example

Let's generate a random 3D polygon and another random point:

{vertices, p} = Through[{Most, Last}[RandomReal[NormalDistribution[0, 1], {4, 3}]]];


To illustrate the use of basis, here are the 3D and local 2D renderings of this configuration:

Graphics3D[{Polygon[vertices], PointSize[0.02], Darker[Red], Point[p]}]
With[{a = basis[vertices]},
Graphics[{Lighter[Gray], Polygon[ vertices . a], Darker[Red], PointSize[0.02],  Point[p . a]}]]


One check that the point actually is in the plane of the polygon is achieved this way:

inPlane[p_, basis_, origin_] := Abs[Det[Append[Transpose[basis], p - origin]] ] < 10.^(-12)


The origin must be a point known to be in the plane (such as one of the polygon's vertices). E.g.,

inPlane[#, basis[vertices], First[vertices]] &  /@ Append[vertices, p]


{True, True, True, False}

verifies that the triangle does lie in its plane and, as it happens, this particular point does not.

-
One way to safen the implementation of basis[] would be to take the Mean[] of the polygon's points and use that as the x... on another note, for your inPlane[], might it be better to do the test as Chop[Det[(* stuff *)]] == 0? –  Guess who it is. Sep 7 '12 at 15:50
@J.M. Then you might prefer this: basis[vertices_] := First[SingularValueDecomposition[vertices - Mean /@ vertices]][[All, 1 ;; 2]] –  whuber Sep 7 '12 at 16:16
I'd do First[SingularValueDecomposition[vertices - Mean /@ vertices, 2]] myself... but wait, shouldn't that be # - Mean[vertices] & /@ vertices (i.e., Standardize[vertices, Mean, 1 &])? –  Guess who it is. Sep 7 '12 at 16:21
@Rainer Assuming the points are very close to coplanar, the smallest eigenvalue will be close to zero and much smaller than the other two. The eigenvectors associated with the other eigenvalues--that is, the first two principal components--thereby span a plane passing close to most of the points. This gives two solutions: (1) The cross product of the first two principal components will be a normal vector for that plane. (2) The third principal component will be orthogonal to the other two, whence it too is a normal vector for the plane. –  whuber Jan 16 '14 at 22:39
@whuber Thanks for the hints I now succeeded in defining a nice, stable function to get the normal vector of a 3D polygon even if this is not exactly planar: GetPolygonNormal[vertices_] := Last[Last[ SingularValueDecomposition[# - Mean[vertices] & /@ vertices]]\[Transpose]] –  Rainer Jan 18 '14 at 16:14

I use the algorithm described here for convex 3D polygons. Basically, if a point is inside a polygon, the sum of the angles between the point and each pair of vertices should be $2\pi$, otherwise it's outside the polygon. The angle between two vectors is given by $$\theta=\arccos\left[{\vec a\cdot\vec b\over \|\vec a\|\|\vec b\|}\right].$$ So,

inPolyQ[poly_,pt_]:=2.π==Total[ArcCos[
Dot@@@#/Times@@@Map[Norm,#,{2}]&@Transpose@{#,RotateRight[#]}
]]&@(#-pt&/@poly)


This only works for a convex polygon, however. Additionally, this doesn't check if the polygon is, in fact, planar.

Edit As per J.M.'s comment, this can be simplified using the VectorAngle function:

inPolyQ[poly_,pt_]:=2.π==Total[
VectorAngle@@@Transpose@{#,RotateRight[#]}
]&@(#-pt&/@poly)

-
You know there's a VectorAngle[] function, don't you? –  Guess who it is. Sep 7 '12 at 14:58
@J.M. Doh! Well, that'll make this simpler-looking... –  Eli Lansey Sep 7 '12 at 15:01
Here is a not-very-efficient way to check if a set of given points are all coplanar: coplanarQ[pts_?MatrixQ] := If[Length[pts] < 4, True, And @@ ((Chop[Det[PadRight[#, {4, 4}, 1]]] == 0) & /@ Subsets[pts, {4}])]. –  Guess who it is. Sep 7 '12 at 15:10
An alternative implementation for inPolyQ[]: inPolyQ[pt_?VectorQ, poly_?MatrixQ] := Chop[Total[VectorAngle @@@ Partition[(# - pt) & /@ poly, 2, 1, {1, 1 - 2 Boole[TrueQ[First[poly] == Last[poly]]]}]] - 2 π] == 0. –  Guess who it is. Sep 7 '12 at 15:25
@J.M. What are the benefits of your approach? –  Eli Lansey Sep 7 '12 at 15:32