How to determine if a point is within a concave 3D polyhedron?
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]
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["NDSolve`FEM`"]
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["NDSolve`FEM`"]
(* 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.
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];
res = Thread[{triangs, params}];
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 Compile
d version of the above.