11
$\begingroup$

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?

my 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 convex hull.

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.

Improve this question
$\endgroup$
0

3 Answers 3

Reset to default
13
$\begingroup$

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]
Improve this answer
$\endgroup$
2
  • $\begingroup$ 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! $\endgroup$ Commented Jun 9, 2020 at 23:01
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Jun 10, 2020 at 2:50
11
$\begingroup$

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

Points in concave mesh

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

Timing Comparison

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

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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). $\endgroup$ Commented Jun 6, 2020 at 6:45
  • 1
    $\begingroup$ @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. $\endgroup$ Commented Jun 6, 2020 at 14:01
  • $\begingroup$ "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. $\endgroup$ Commented Jun 6, 2020 at 14:16
6
$\begingroup$

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 Compiled version of the above.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.