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];
res = Thread[{triangs, params}];
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.
