# Efficiently determining if 3D points are within a surface composed of polygons

This is the 2nd part of a previous question which I edited to make into 2 separate questions: Extracting polygons from 3D contour plot surface

As an extension of my earlier question involving simple convex polyhedra (Checking if a point is in a convex 3D polyhedron), I now want to identify points which are contained in an arbitrary polygonal surface which was generated with a ListContourPlot3D.

Here is the surface defined by polygons:

data = Table[x^2 + y^2 + z^2 + RandomReal[0.1], {x, -2, 2, 0.2},
{y, -2, 2,  0.2}, {z, -2, 2, 0.2}];
plot = ListContourPlot3D[data, Contours -> {1}, Mesh -> None,  ImageSize -> 400];
polygonCoords =
Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}]
Graphics3D[{Opacity[.9], EdgeForm[Opacity[.3]],
Polygon[#,  VertexColors -> Table[Hue[RandomReal[]], {Length[#]}]] & /@
Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}]},
Lighting -> "Neutral", ImageSize -> 400]


Here are some random points

points = RandomReal[{-2,2},{100,3}];
ListPointPlot3D[points]


How do I determine which points are within the surface? It doesn't seem that you can generate a RegionFunction for the surface. Does TetGenLink offer functionality like this?

@JosephO'Rourke provided a hint in his answer (http://mathematica.stackexchange.com/a/15198/65) to my simpler question that some code accompanying his book might provide the solution. The code is available in C and Java. I'm not familiar with MathLink or LibraryLink to integrate with the C code and my Java is very rusty but maybe using @Leonid's Java Reloader (http://mathematica.stackexchange.com/a/6377/65) it might be possible to call the java implementation (source code here: http://cs.smith.edu/~orourke/books/ftp.html).

-

Note that ListContourPlot3D takes the coordinates to be the position indices by default. If you want to keep the coordinates used in generating the data, then you have to include it.

data = Flatten[
Table[{x, y, z, x^2 + y^2 + z^2 + RandomReal[0.1]},
{x, -2, 2, 0.2}, {y, -2, 2, 0.2}, {z, -2, 2, 0.2}], 2];

plot = ListContourPlot3D[data, Contours -> {1}, Mesh -> None];

polygonCoords =
Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}];

g = Graphics3D[{Opacity[.09], EdgeForm[Opacity[.3]],
Polygon[#,
VertexColors -> Table[Hue[RandomReal[]], {Length[#]}]] & /@
Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}]},
Lighting -> "Neutral", ImageSize -> 400, Axes -> True];


The side of an outwardly oriented (by the order of the vertices) convex polygon can be determined by the determinant of the differences of the point X and three consecutive vertices P, Q, R. In this case, the signs should all be positive to indicate the point is inside the polyhedron.

side[{P_, Q_, R_, ___}, X_] := Det@Differences[{X, P, Q, R}];
insideQ[polyhedron_, point_] := And @@ Positive[side[#, point] & /@ polyhedron];


Example:

points = RandomReal[{-1.5, 1.5}, {1000, 3}];

pg = Graphics3D[
Point[points,
VertexColors -> (insideQ[polygonCoords, #] & /@
points /. {True -> Red, False -> Directive[Opacity[0.5], Blue]})],
Axes -> True
];

Show[g, pg]


-

Here is a more general approach. It is based on the 2D method from here. It assumes the polyhedron is not self-intersecting but imposes no requirement of convexity or even connectedness, other than that it be closed and bounded. Strictly speaking, I think this will work for an unbounded polyhedron provided it contains no vertical ray.

For ease of exposition consider that the object lies above the x-y plane. The idea is to make a grid in the x-y plane, and for each box figure out all triangles of the surface that at least partly lie over that box. This is done in a preprocessing step. It is a bit crude in that the code makes no effort to avoid overcounting, or to be "smart" about iterations that find cases where an entire grid box might be inside some triangle. We simple rule out duplicates at the end. It is possible (though I think I got this much correct) that some candidates actually do not intersect a claimed grid box. In practice if this is happening, it is not at a to an extent sufficien to cause bad speed issues.

Once given a point we use the grid box it lies over to find candidate triangles that a vertical line from x-y plane to infinity, through that point, intersects. This may be a subset of the set of triangles we associated to that grid box in the preprocessing. Next figure out how many of these surface triangles the given point lies above. If even, it is outside, else inside.

Once again I made no effort to encapsulate the code, and there are global variables used inside several functions. This could of course be improved if the end goal is to have generally robust production code. Also I took no great pain to speed the preprocessing. For a polyhedron of a few thousand triangles this part takes maybe a couple of seconds. Significantly larger examples will require either more patience or perhaps adept use of Compile.

I will say that I found a need to rotate the polyhedron ever so slightly to avoid a degenerate situation wherein a triangle in 3-space projects to a segment in the x-y plane. That is to say, the preprocessing works with the rotated version. The strategy is to rotate the given points in the same way, and use them for the in-or-out testing, but plot the original polyhedron and test points. About all I can say is I hope I got this part correct. (The lack of 1/0 error messages, and a decent-looking plot, give me some confidence in this respect. But that warm feeling hardly a proof of correctness and for all I know might actually indicate incontinence.)

Here is the code. Some is taken from my previous response. Also I use a variant of the running example, this time modified to have more faces.

First the data:

data = Flatten[
Table[{x, y, z, x^2 + y^2 + z^2 + RandomReal[{-.1, .1}]}, {x, -2,
2, 0.2}, {y, -2, 2, 0.2}, {z, -2, 2, 0.2}], 2];

plot = ListContourPlot3D[data, Contours -> {3}, Mesh -> None,
ImageSize -> 400, MaxPlotPoints -> 30]
makeTriangles[tri : {aa_, bb_, cc_}] := {tri}
makeTriangles[{aa_, bb_, cc_, dd__}] :=
Join[{{aa, bb, cc}}, makeTriangles[{aa, cc, dd}]]


Now rotate the data a small amount around some random axis through the origin. Might be better to do this with one centered at the barycenter, but I didn't go that far.

{theta, axis} = {RandomReal[{.1}], RandomReal[1, {3}]}

(* Out[164]= {0.0581784648287, {0.606522492066, 0.357476057555,
0.623359098046}} *)


We will base the polygons on this rotated data.

plotr = plot /. plot[[1]] :> Rotate[plot[[1]], theta, axis];

polygonCoords =
Cases[Normal[plotr[[1]]], Polygon[x_, ___] :> x, {0, Infinity}];
polygonCoords // Length

(* Out[120]= 2966 *)


Here is the opaque version of the unrotated form.

g =
Graphics3D[{Opacity[.09], EdgeForm[Opacity[.3]],
Polygon[#,
VertexColors -> Table[Hue[RandomReal[]], {Length[#]}]] & /@
Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}]},
Lighting -> "Neutral", ImageSize -> 400, Axes -> True]


makeTriangles[tri : {aa_, bb_, cc_}] := {tri}
makeTriangles[{aa_, bb_, cc_, dd__}] :=
Join[{{aa, bb, cc}}, makeTriangles[{aa, cc, dd}]]

Clear[triangles, pts];
triangles = Flatten[Map[makeTriangles, polygonCoords], 1];
triangles = triangles;
pts = Union[Flatten[triangles, 1]];
Length[triangles]


(* Out[247]= 5924 *)

So we are working with around 6000 triangles.

flats = Map[Most, triangles, {2}];
xcoords = pts[[All, 1]];
ycoords = pts[[All, 2]];
zcoords = pts[[All, 3]];
xmin = Min[xcoords];
ymin = Min[ycoords];
xmax = Max[xcoords];
ymax = Max[ycoords];
zmin = Min[zcoords];
zmax = Max[zcoords];


I'll use a 30x30 grid.

n = 30;
xspan = xmax - xmin;
yspan = ymax - ymin;
dx = 1.05*xspan/n;
dy = 1.05*yspan/n;
midx = (xmax + xmin)/2;
midy = (ymax + ymin)/2;
xlo = midx - 1.05*xspan/2;
ylo = midy - 1.05*yspan/2;


The code below finds triangles that lie over grid boxes. A triangle does so if a vertex is over the grid box, or an edge projects to a segment intersecting the grid box (remember to assign the triangle to both neighbors) or a vertex of the grid box lies inside the projected triangle, in which case there are now four neighbors that get assigned that triangle.

edges[{a_, b_, c_}] := {{a, b}, {b, c}, {c, a}}

vertexBox[{x1_, y1_}, {xb_, yb_, dx_, dy_}] := {Ceiling[(x1 - xb)/dx],
Ceiling[(y1 - yb)/dy]}

segmentBoxes[{{x1_, y1_}, {x2_, y2_}}, {xb_, yb_, dx_, dy_}] :=

Module[
{xmin, xmax, ymin, ymax, xlo, xhi, ylo, yhi, xtable, ytable, xval,
yval, index},
xmin = Min[x1, x2];
xmax = Max[x1, x2];
ymin = Min[y1, y2];
ymax = Max[y1, y2];
xlo = Ceiling[(xmin - xb)/dx];
ylo = Ceiling[(ymin - yb)/dy];
xhi = Ceiling[(xmax - xb)/dx];
yhi = Ceiling[(ymax - yb)/dy];
xtable = Flatten[Table[
xval = xb + j*dx;
yval = (((-x2)*y1 + xval*y1 + x1*y2 - xval*y2))/(x1 - x2);
index = Ceiling[(yval - yb)/dy];
{{j, index}, {j + 1, index}}
, {j, xlo, xhi - 1}], 1];
ytable = Flatten[Table[
yval = yb + j*dy;
xval = (((-y2)*x1 + yval*x1 + y1*x2 - yval*x2))/(y1 - y2);
index = Ceiling[(xval - xb)/dx];
{{index, j}, {index, j + 1}}
, {j, ylo, yhi - 1}], 1];
Union[Join[xtable, ytable]]
]

pointInsideTriangle[
p : {x_, y_}, {{x1_, y1_}, {x2_, y2_}, {x3_, y3_}}] := With[
{l1 = -((x1*y - x3*y - x*y1 + x3*y1 + x*y3 - x1*y3)/
(x2*y1 - x3*y1 - x1*y2 + x3*y2 + x1*y3 - x2*y3)),
l2 = -(((-x1)*y + x2*y + x*y1 - x2*y1 - x*y2 + x1*y2)/
(x2*y1 - x3*y1 - x1*y2 + x3*y2 + x1*y3 - x2*y3))},
Min[x1, x2, x3] <= x <= Max[x1, x2, x3] &&
Min[y1, y2, y3] <= y <= Max[y1, y2, y3] && 0 <= l1 <= 1 &&
0 <= l2 <= 1 && l1 + l2 <= 1]

faceBoxes[
t : {{x1_, y1_}, {x2_, y2_}, {x3_, y3_}}, {xb_, yb_, dx_, dy_}] :=
Catch[Module[
{xmin, xmax, ymin, ymax, xlo, xhi, ylo, yhi, xval, yval, res},
xmin = Min[x1, x2, x3];
xmax = Max[x1, x2, x3];
ymin = Min[y1, y2, y3];
ymax = Max[y1, y2, y3];
If[xmax - xmin < dx || ymax - ymin < dy, Throw[{}]];
xlo = Ceiling[(xmin - xb)/dx];
ylo = Ceiling[(ymin - yb)/dy];
xhi = Ceiling[(xmax - xb)/dx];
yhi = Ceiling[(ymax - yb)/dy];
res = Table[
xval = xb + j*dx;
yval = yb + k*dy;
If[pointInsideTriangle[{xval, yval},
t], {{j, k}, {j + 1, k}, {j, k + 1}, {j + 1, k + 1}}, {}]
, {j, xlo, xhi - 1}, {k, ylo, yhi - 1}];
res = res /. {} :> Sequence[];
Flatten[res, 2]
]]

gridBoxes[pts : {a_, b_, c_}, {xb_, yb_, dx_, dy_}] := Union[Join[
Map[vertexBox[#, {xb, yb, dx, dy}] &, pts],
Flatten[Map[segmentBoxes[#, {xb, yb, dx, dy}] &, edges[pts]], 1],
faceBoxes[pts, {xb, yb, dx, dy}]
]]


So here is the main preprocessing step.

Timing[
gb = DeleteCases[
Map[gridBoxes[#, {xlo, ylo, dx, dy}] &,
flats], {a_, b_} /; (a > n || b > n), 2];
grid = ConstantArray[{}, {n, n}];
Do[Map[AppendTo[grid[[Sequence @@ #]], j] &, gb[[j]]], {j,
Length[gb]}];]

(* Out[267]= {1.850000, Null} *)


Now we need code to decide when, given a point and a candidate triangle from the grid box for that point, decides whether the triangle really lies over the points projection into the grid box.

planeTriangleParams[
p : {x_, y_}, {p1 : {x1_, y1_}, p2 : {x2_, y2_}, p3 : {x3_, y3_}}] :=
With[
{den = x2*y1 - x3*y1 - x1*y2 + x3*y2 + x1*y3 - x2*y3},
{-((x1*y - x3*y - x*y1 + x3*y1 + x*y3 - x1*y3)/
den), -(((-x1)*y + x2*y + x*y1 - x2*y1 - x*y2 + x1*y2)/den)}]

getTriangles[p : {x_, y_}] := Module[
{ix, iy, triangs, params, res},
{ix, iy} = vertexBox[p, {xlo, ylo, dx, dy}];
triangs = grid[[ix, iy]];
params = Map[planeTriangleParams[p, flats[[#]]] &, triangs];
Select[res,
0 <= #[[2, 1]] <= 1 &&
0 <= #[[2, 2]] <= 1 && #[[2, 1]] + #[[2, 2]] <= 1.0000001 &]
]


This last code counts the triangle on the polyhedron surface that a given point lies over, and declares inside iff that number is odd. For purposes of catching bad cases, I have a Print that fires when the number of faces a vertical through the point intersects is odd. It caught several bugs for me along the way.

countAbove[p : {x_, y_, z_}] := Module[
{triangs = getTriangles[Most[p]], threeDtriangs, lambdas, zcoords,
zvals},
threeDtriangs = triangles[[triangs[[All, 1]]]];
lambdas = triangs[[All, 2]];
zcoords = threeDtriangs[[All, All, 3]];
zvals =
Table[zcoords[[j, 1]] +
lambdas[[j, 1]]*(zcoords[[j, 2]] - zcoords[[j, 1]]) +
lambdas[[j, 2]]*(zcoords[[j, 3]] - zcoords[[j, 1]]), {j,
Length[zcoords]}];
If[OddQ[Length[triangs]] && OddQ[Length[Select[zvals, z > # &]]],
Print[{p, triangs, Length[Select[zvals, z > # &]]}]];
Length[Select[zvals, z > # &]]
]

isInside[{x_, y_,
z_}] /; ! ((xmin <= x <= xmax) && (ymin <= y <= ymax) && (zmin <=
z <= zmax)) := False
isInside[p : {x_, y_, z_}] := OddQ[countAbove[p]]


Here is a test using the figure above. It was constructed in such a way that min and max coordinate values are all around -1.7 and 1.7 respectively, so I'll use random points in the cube from -2 to 2 in all coordinates.

points = RandomReal[{-2., 2.}, {10000, 3}];
pointsr = Map[RotationTransform[theta, axis], points];

Timing[
pg = Graphics3D[
Point[points,
VertexColors -> (isInside /@ pointsr /. {True -> Red,
False -> Directive[Opacity[0.15], Blue]})], Axes -> True,
AxesLabel -> {"x", "y", "z"}];]

(* Out[272]= {2.860000, Null} *)


Now the picture.

Show[g, pg, AxesLabel -> {"x", "y", "z"}]


Given the number of issues I ran into along the way in coding this, I do not expect ever to be upping it to four dimensions.

-

Here is a fast method that will "often" work. Roughly, it requires that the convex polygon have no sharp angles between faces. Preprocessing goes as follows.

Create triangles from the polygons. So a 5-gon with vertices {a,b,c,d,e} would become the set of triangles {{a,b,c},{a,c,d},{a,d,e}}.

For each vertex we average it's star (set of points connected by one edge). The idea is to provide what we hope will be an inward-pointing vector for that vertex.

For each triangle we need an inward pointing normal. We attempt to get this quickly by "aligning with the stars" (no, I didn't plan my notation to use that line, it just...sorta...fell...from the heavens). That is to say, we take a normal for a given triangle, and, if it is more than 90 degrees from the star for the first vertex of that triangle, then we negate it. We now check this against the the stars for the other two vertices. With luck we get a consistent alignment there. If not, we go with the direction that agrees with two stars and hope for the best.

Use Nearest to create a NearestFunction that, for a given point, will find the closest vertex and also the list of inward-pointing vectors of triangles containing that vertex.

That ends the preprocessing. The code I use is not properly modularized and uses AppendTo, so it might not be optimal if there are many vertices/triangles. That could easily be improved, were it to become a bottleneck.

The inside[] predicate finds the closest vertex to a given point, and the normals of its containing triangles. If computes the vector from nearest vertex to the given point. If that is less than 90 degrees from all the surrounding normal vectors then the point is deemed to be inside, else outside.

I'll show this for the running example, using the setup from @Michael E2.

data = Flatten[
Table[{x, y, z, x^2 + y^2 + z^2 + RandomReal[0.1]}, {x, -2, 2,
0.2}, {y, -2, 2, 0.2}, {z, -2, 2, 0.2}], 2];
plot = ListContourPlot3D[data, Contours -> {1}, Mesh -> None,
ImageSize -> 400];
polygonCoords =
Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}];

g = Graphics3D[{Opacity[.1], EdgeForm[Opacity[.3]],
Polygon[#,
VertexColors -> Table[Hue[RandomReal[]], {Length[#]}]] & /@
Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}]},
Lighting -> "Neutral", ImageSize -> 400];

makeTriangles[tri : {aa_, bb_, cc_}] := {tri}
makeTriangles[{aa_, bb_, cc_, dd__}] :=
Join[{{aa, bb, cc}}, makeTriangles[{aa, cc, dd}]]

Clear[triangles, pts];
triangles = Flatten[Map[makeTriangles, polygonCoords], 1];
pts = Union[Flatten[triangles, 1]];

Clear[triangleList];
Map[(triangleList[#] = {}) &, pts];
MapIndexed[AppendTo[triangleList[#], #2[[1]]] &, triangles, {2}];

Clear[ptStars, means];
Map[(ptStars[#] = {}) &, pts];
Map[AppendTo[ptStars[#], triangles[[triangleList[#]]]] &, pts];
Map[(ptStars[#] = Complement[Union[Flatten[ptStars[#], 2]], {#}]) &,
pts];
Map[(means[#] = Total[ptStars[#]]/Length[ptStars[#]]) &, pts];

Clear[triangleStars];
triangleStars =
Table[Join @@ Table[ptStars[triangles[[j, k]]], {k, 3}], {j,
Length[triangles]}];

grad[tri_] := Cross[tri[[1]] - tri[[2]], tri[[1]] - tri[[3]]]

tri = triangles[[j]];
pt = tri[[1]];
, {j, Length[triangles]}];


In this case nothing is printed, ergo the normals all aligned normally, so to speak.

nf = Nearest[Map[# -> {#, triangleList[#]} &, pts]];

inside[pt_] := Module[
{near = nf[pt][[1]], diff, triangs, ptgrads},
diff = near[[1]] - pt;
triangs = near[[2]];
And @@ (Map[diff.# > 0 &, ptgrads])
]


This example will give some indication of overall speed.

points = RandomReal[{-1, 1}, {10000, 3}];
Timing[pg =
Graphics3D[
Point[points,
VertexColors -> (inside /@ points /. {True -> Red,
False -> Directive[Opacity[0.5], Blue]})], Axes -> True];]

(* Out[355]= {0.340000, Null} *)


Here is the picture. It looks basically correct but it is certainly possible that some small percentage of points are incorrectly marked.

Show[g, pg]


-
Twinkle twinkle... stellar answer! –  Yves Klett Feb 25 '13 at 20:34
Wow, thanks (to both answerers). Its going to take me a bit to benchmark on real data and work on compiling your answer. I promise I will accept an answer soon! –  s0rce Feb 26 '13 at 3:34
I actually have a "better" answer (does not rely on convexity or other niceness). But it is buggy and marks a small percentage of points inside that actually are not. Will post if I get it sorted out. –  Daniel Lichtblau Feb 26 '13 at 23:44

In Version 10 we have BoundaryDiscretizeGraphics and RegionMember to the rescue. So here we go:

Reproducing the graphics:

data = Flatten[Table[{x, y, z, x^2 + y^2 + z^2 + RandomReal[0.1]}, {x, -2, 2,
0.2}, {y, -2, 2, 0.2}, {z, -2, 2, 0.2}], 2]; (* Thanks to Taliesin Beynon for the tip *)

plot = ListContourPlot3D[data, Contours -> {1}, Mesh -> None, ImageSize -> 400];

polygonCoords = Cases[Normal[plot[[1]]], Polygon[x_, ___] :> x, {0, Infinity}]

gr = Graphics3D[{Opacity[.9], EdgeForm[Opacity[.3]],
Polygon[#, VertexColors -> Table[Hue[RandomReal[]], {Length[#]}]] & /@
polygonCoords}, Lighting -> "Neutral", ImageSize -> 400]


We create a Region using BoundaryDiscretizeGraphics:

reg = BoundaryDiscretizeGraphics[gr]


We create the points:

pts = RandomReal[{-2, 2}, {1000, 3}];


Now we ask which points are inside the Region:

inreg = RegionMember[reg, #] & /@ pts;


Visualization time:

Show[HighlightMesh[reg, Style[2, Opacity[0.5]]],
Graphics3D[{Red, Point[Pick[pts, inreg]], Black, Point[Pick[pts, inreg, False]]}],
Boxed -> True]


OR with the original graphics:

Show[gr /. Opacity[_] -> Opacity[0.3],
Graphics3D[{Red, Point[Pick[pts, inreg]], Black,
Point[Pick[pts, inreg, False]]}], Boxed -> True]


-
+1 These now region-functions are so exciting. One can even control the polygons with MaxCellMeasure. –  eldo Jul 20 '14 at 19:26
@eldo Thanks. I agree that these new Region-functions are exciting. They were one of the main reasons I quickly upgraded to V10. –  RunnyKine Jul 20 '14 at 19:29
I had just finished typing up the same solution when you posted this! Nice. Note that it would be good to generate the data via data = Flatten[Table[{x, y, z, x^2 + y^2 + z^2 + RandomReal[0.1]}, {x, -2, 2, 0.2}, {y, -2, 2, 0.2}, {z, -2, 2, 0.2}], 2]; so that the coordinate system of the resulting region is the same as what was implicitly used by the Table, instead of the coordinates being the indices of that original array. –  Taliesin Beynon Jul 20 '14 at 20:22
@TaliesinBeynon, Interesting that we were typing up the same method at the same time :). Thanks for the tip on generating the data, I was wondering why the points were off initially. I'll update with your suggestion. –  RunnyKine Jul 20 '14 at 20:37