9
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • 3
    $\begingroup$ Why not just apply a linear transformation sending the polygon's plane to the xy plane and proceed from there? $\endgroup$ Commented Sep 7, 2012 at 15:07
  • $\begingroup$ @whuber I was thinking about that, but I'm not that familiar with FindGeometricTransform. $\endgroup$ Commented Sep 7, 2012 at 15:09
  • 2
    $\begingroup$ @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. $\endgroup$ Commented Sep 7, 2012 at 15:22
  • $\begingroup$ @Simon, for polygons with inexact coordinates, it'd be safer to use Newell's algorithm to compute the normal... $\endgroup$ Commented Sep 7, 2012 at 15:31
  • $\begingroup$ @SimonWoods That's a great idea. Want to write it up as an answer? $\endgroup$ Commented Sep 7, 2012 at 15:32

3 Answers 3

Reset to default
7
$\begingroup$

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]

Mathematica graphics

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

Mathematica graphics

Improve this answer
$\endgroup$
11
$\begingroup$

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.

Improve this answer
$\endgroup$
10
  • 1
    $\begingroup$ 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? $\endgroup$ Commented Sep 7, 2012 at 15:50
  • 1
    $\begingroup$ @J.M. Then you might prefer this: basis[vertices_] := First[SingularValueDecomposition[vertices - Mean /@ vertices]][[All, 1 ;; 2]] $\endgroup$ Commented Sep 7, 2012 at 16:16
  • 1
    $\begingroup$ 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 &])? $\endgroup$ Commented Sep 7, 2012 at 16:21
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Jan 16, 2014 at 22:39
  • 1
    $\begingroup$ @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]] $\endgroup$ Commented Jan 18, 2014 at 16:14
4
$\begingroup$

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)
Improve this answer
$\endgroup$
7
  • 2
    $\begingroup$ You know there's a VectorAngle[] function, don't you? $\endgroup$ Commented Sep 7, 2012 at 14:58
  • $\begingroup$ @J.M. Doh! Well, that'll make this simpler-looking... $\endgroup$ Commented Sep 7, 2012 at 15:01
  • 1
    $\begingroup$ 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}])]. $\endgroup$ Commented Sep 7, 2012 at 15:10
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Sep 7, 2012 at 15:25
  • $\begingroup$ @J.M. What are the benefits of your approach? $\endgroup$ Commented Sep 7, 2012 at 15:32

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.