# How to determine if a point is within a concave 3D polyhedron?

Generally, if you have a 3D polyhedron and wanted to check if a point was within it, you would use something like a ConvexHullMesh to create a region, which you can then use RegionMemberQ to check if a point was within it.

But, this technique will not work for concave polyhedra. I have a programme which generates points to make a surface from. This works well, and I have posted the points & surface in a Pastebin.

points = Import["https://pastebin.com/raw/190HQui1"];
polygon = Import["https://pastebin.com/raw/d3MRBb8K"];

rmesh = Region[polygon];
Show[rmesh, points]


Now, how would I check if a point is within this shape?

I feel it is worth noting that RegionDistance[polygon]works, but only generates a 2-dimensional object - which works as expected - but we want to know if we are in the polyhedron. ConvexHullMesh[polygon] is a poor approximation .

There are these solutions to determine if a point is within a 2D polygon (even a convex one) (1 2). But they don't seem directly applicable to the 3D case.

You can try this:

polygon = Import["https://pastebin.com/raw/d3MRBb8K"];
pts = Union @@ polygon[[1]];
nf = Nearest[pts -> "Index"];
R = BoundaryMeshRegion[pts, Polygon[DeleteDuplicates@*Flatten /@ Map[nf, polygon[[1]], {2}]]];
f = RegionMember[R]

• This results in very strange behaviour with RegionDistance[R], e.g if you run rdf = RegionDistance[R]; DensityPlot[rdf[{x, y, 1.1}], {x, -10, 20}, {y, -10, 20}]  You can see all the values which are smaller than 0 are considered in the volume!
– Tomi
Jun 9, 2020 at 23:01
• That is indeed very odd. Maybe BoundaryMeshRegion does not orient the faces correctly, leading to some face normals that point inward instead of outward. I am not sure whether this is to be considered a bug. It seems that BoundaryMeshRegion @*ToBoundaryMesh is be the more robust way to create the BoundaryMeshRegion. Jun 10, 2020 at 2:50

Here is an alternative approach using SignedRegionDistance that seems pretty fast, but I have not compared it to @Henrik Schumacher's answer. It took about 5 seconds to test 100,000 points on my machine.

Needs["NDSolveFEM"]
points = Import["https://pastebin.com/raw/190HQui1"];
polygon = Import["https://pastebin.com/raw/d3MRBb8K"];
(* Convert into BoundaryMeshRegion *)
bmr = BoundaryMeshRegion[ToBoundaryMesh[Region[polygon]]];
(* create a SignedRegionDistance function *)
srdf = SignedRegionDistance[bmr];
(* create some random coodinates *)
crd = RandomReal[10, {100000, 3}];
(* If srdf is <0, then point is in region *)
inRegQ = PositionIndex[srdf[#] < 0 & /@ crd];
(* Show outside Points in Red and inside in Green *)
Show[Graphics3D[{{Red, Point[crd[[inRegQ[False]]]]}, {Green,
Point[crd[[inRegQ[True]]]]}}]]
(* Show points in region only *)
Show[RegionPlot3D[bmr, PlotStyle -> Directive[Yellow, Opacity[0.25]],
Mesh -> None], Graphics3D[{{Green, Point[crd[[inRegQ[True]]]]}}]]


# Timing Comparison

Since Henrik was so kind to speed up my code, I replicated some repeated timings on the various permutations.

(* Henrik's Answer *)
polygon = Import["https://pastebin.com/raw/d3MRBb8K"];
pts = Union @@ polygon[[1]];
nf = Nearest[pts -> "Index"];
R = BoundaryMeshRegion[pts,
Polygon[DeleteDuplicates@*Flatten /@ Map[nf, polygon[[1]], {2}]]];
f = RegionMember[R];
Needs["NDSolveFEM"]
(* Convert into BoundaryMeshRegion *)
bmr = BoundaryMeshRegion[ToBoundaryMesh[Region[polygon]]];
(* create SignedRegionDistance function based on bmr *)
srdfbmr = SignedRegionDistance[bmr];
(* create SignedRegionDistance function based on R*)
srdfr = SignedRegionDistance[R];
(* create some random coodinates *)
crd = RandomReal[10, {100000, 3}];
(* Henrik's Solution *)
{timeHS, inRegQ} = RepeatedTiming@PositionIndex[f[crd]];
(* Tim Laska's Original Solution *)
{timeTL, inRegQ} =
RepeatedTiming@PositionIndex[srdfbmr[#] < 0 & /@ crd];
(* Tim Laska's With Henrik's UnitStep Suggestion *)
{timeHSSug, inRegQ} =
RepeatedTiming@
PositionIndex[{True, False}[[UnitStep[srdfbmr[crd]] + 1]]];
(* Tim Laska's With Henrik's Polygon *)
{timeTLR, inRegQ} =
RepeatedTiming@PositionIndex[srdfr[#] < 0 & /@ crd];
(* Tim Laska's With Henrik's UnitStep Suggestion and His Polygon *)
{timeHSSugPoly, inRegQ} =
RepeatedTiming@
PositionIndex[{True, False}[[UnitStep[srdfr[crd]] + 1]]];
data = {{"Henrik's Answer", timeHS}, {"Tim's Original",
timeTL}, {"Tim's with Henrik's UnitStep",
timeHSSug}, {"Tim's with Henrik's Poly",
timeTLR}, {"Tim's with Henrik's Poly and UnitStep",
timeHSSugPoly}};
data = SortBy[data, Last];
Text@Grid[Prepend[data, {"Method", "Time(s)"}],
Background -> {None, {Lighter[Yellow, .9], {White,
Lighter[Blend[{Blue, Green}], .8]}}},
Dividers -> {{Darker[Gray, .6], {Lighter[Gray, .5]},
Darker[Gray, .6]}, {Darker[Gray, .6], Darker[Gray, .6], {False},
Darker[Gray, .6]}}, Alignment -> {{Left, Right, {Left}}},
ItemSize -> {{20, 5}}, Frame -> Darker[Gray, .6], ItemStyle -> 14,
Spacings -> {Automatic, .8}]


On my machine, Henrik's UnitStep suggestion boosted performance about 3x. The performance of RegionMember and SignedRegionDistance are similar with Henrik's suggestion.

• Mapping srdf makes it quite slow. {True, False}[[UnitStep[srdf[crd]] + 1]] is almost 4 times faster; and so is RegionMember. I don't know what ToBoundaryMesh does exactly, but for some reason, using the region Ras created by me instead of bmr seems to be slightly faster (and leads to the same results). Jun 6, 2020 at 6:45
• @HenrikSchumacher Thank you very much for your insight. I will try to incorporate your suggestions soon. I suspect the triangle count changes as ToBoundaryMesh will try to make the polygons isotropic. Jun 6, 2020 at 14:01
• "will try to make the polygons isotropic."Ah right, I get it. This is a good thing for many tasks, but for this one, it is just not necessary. Jun 6, 2020 at 14:16

Here is a method that takes around 2-2.5 times longer than the one from @TimLaska. It has the advantage that it can perhaps be made considerably faster using Compile. It is code from here that I adjusted slightly for the problem at hand.

The main idea is to find boundary triangles that a ray from the outside to the given point can intersect. We count these; odd means point is inside. I used a random transformation to avoid zero denominators that can arise with data that is too well "aligned" with one or more coordinate axes.

points0 = Import["https://pastebin.com/raw/190HQui1"];
pgon0 = Import["https://pastebin.com/raw/d3MRBb8K"];

SeedRandom[1234];
randpt = RandomReal[1, 3];
translate = TranslationTransform[randpt];
randdir = RandomReal[1, 3];
theta = RandomReal[Pi];
rotate = RotationTransform[theta, randdir];
transform = Composition[rotate, translate];

rmesh0 = Region[pgon0];

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

triangles =
Map[transform,
Flatten[Map[makeTriangles, rmesh0[[1, 1]]], 1], {2}];
verts = Map[transform, points0[[All, 1, 1]]];

flats = Map[Most, triangles, {2}];
pts = verts;
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];

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

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


Creating the main structure takes a bit of up-front time.

AbsoluteTiming[
gbox = 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] &, gbox[[j]]], {j,
Length[gbox]}];]

(* Out[2893]= {1.47625, Null} *)

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

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


Running it takes 8.8 seconds.

SeedRandom[12345];
crd = Map[transform, RandomReal[10, {100000, 3}]];
AbsoluteTiming[inRegQ = Map[isInside, crd];]

(* Out[2906]= {8.83544, Null} *)


The code from Tim Laska took around 4.3 seconds on this machine for the same point set. I suspect that could be attained by a Compiled version of the above.