Tell me more ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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 Graphics`Mesh`PointWindingNumber 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?

share|improve this question
1  
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
1  
@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... – J. M. 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
show 1 more comment

2 Answers

up vote 7 down vote accepted

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]}]]

Two plots


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.

share|improve this answer
Can you give more detail about checking if the point is in the plane? – Eli Lansey Sep 7 '12 at 15:34
@Eli: equation 18 here might give you a clue on how to do the checking. (I used it in the coplanarQ[] function I mentioned in another comment.) – J. M. Sep 7 '12 at 15:38
1  
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? – J. M. Sep 7 '12 at 15:50
@J.M. Right; Chop works well. I considered using the centroid, but it potentially has similar numerical problems for many polygons. If you're worried about this, perhaps the most efficient approach is to choose, say, five or six vertices at random and for the basis calculation use three that form a sufficiently large angle between them. A less efficient but very stable and trustworthy method is to take the first two principal components in an SVD of the vertices(considered as a $3$ by $n$ matrix) relative to their centroid: that is, a centered PCA. – whuber Sep 7 '12 at 16:02
Hah, nice! I forgot that PCA can do the job here, even though it is not the most efficient way to go about things... but safe nevertheless. – J. M. Sep 7 '12 at 16:12
show 2 more comments

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)
share|improve this answer
2  
You know there's a VectorAngle[] function, don't you? – J. M. Sep 7 '12 at 14:58
@J.M. Doh! Well, that'll make this simpler-looking... – Eli Lansey Sep 7 '12 at 15:01
1  
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}])]. – J. M. Sep 7 '12 at 15:10
1  
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. – J. M. Sep 7 '12 at 15:25
@J.M. What are the benefits of your approach? – Eli Lansey Sep 7 '12 at 15:32
show 2 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.