How to find the cells of a region that intersect a line?

Question

I'm looking for the fastest way to access the mesh cell (i.e. polygon) in a Region that intersects a Line in 3D. For example,

reg = BoundaryDiscretizeRegion[Ball[], PrecisionGoal -> 1, 
   MaxCellMeasure -> 0.1];
Graphics3D[{{Blue, Thick, Line[{{-1, -1, -1}, {0, 0, 0}}]}, 
  Opacity[0.8], reg}, Axes -> True]

Here's a slow naive way:

ps = MeshPrimitives[reg, 2];
intersections = 
  If[EmptyRegion[3] =!= 
      RegionIntersection[#, Line[{{0, 0, 0}, {-1, -1, -1}}]], #, 
     Nothing] & /@ ps;
Graphics3D[{{Red, Line[{{0, 0, 0}, {-1, -1, -1}}]}, intersections}, 
 Axes -> True]

There must be an easy efficient way to find the cell that intersects the line, I'm just not seeing it.

What version are you using?, my 11.3 says that BDR is not a graphics3d primitive. — Kuba, Dec 11 '18 at 07:31

Henrik Schumacher · Accepted Answer · 2018-12-11T10:12:09.463

2

Before I begin, I'd like to say that this is indeed something that I really miss in Mathematica 11.3.

Well, we can try to employ RegionDistance to find the intersection point as follows:

R = DiscretizeRegion[Sphere[], MaxCellMeasure -> {1 -> 0.5}];
L = Line[{{0, 0, 0}, {-1, -1, -1}}];

d = RegionDistance[R];
curve = (1 - t) L[[1, 1]] + t L[[1, 2]];
pt = curve /. FindRoot[d[curve] == 0, {t, 0, 1}, Method -> "Secant"]

{-0.576338, -0.576338, -0.576338}

For more than one intersection, one may use the good ol' MeshFunction approach:

gc = Cases[
   Plot[d[curve], {t, 0, 1}, MeshFunctions -> {#2 &}, 
    Mesh -> {{0.01}}, MeshStyle -> Red],
   _GraphicsComplex,
   ∞
   ];
candidates = Join[{0.}, Sort[gc[[1, 1, Cases[gc, _Point, ∞][[1, 1]], 1]]], {1.}];
pts = curve /. FindRoot[
 d[curve] == 0, 
 {t, #[[1]], #[[2]]}, 
 Method -> "Secant"] & /@ Partition[candidates, 2, 1];

This returns only the intersection points, not the intersected triangles. But the cell indices of the nearest mesh cells of highest dimension(?) can be obtained by

intersectedtriangles = DeleteDuplicates[Region`Mesh`MeshNearestCellIndex[R, pts]]

{{2, 744}}

Graphics3D[{Thick, L, MeshPrimitives[R, intersectedtriangles]}]

edited Dec 11 '18 at 10:12

answered Dec 11 '18 at 08:58

Henrik Schumacher

106,770
7
179
309

@ HenrikSchumacher Very helpful! Where could I find further documentation to RegionMeshMeshNearestCellIndex? Thanks. – Ulrich Neumann Dec 11 '18 at 15:55
@Ulrich As always with these undocumented functions such as Region`Mesh`MeshNearestCellIndex: Nowhere. Unfortunately. =/ – Henrik Schumacher Dec 11 '18 at 16:03
It's a pity, thanks – Ulrich Neumann Dec 11 '18 at 16:19
@HenrikSchumacher I'm seeing Part::partw: Part 3 of {Identity,Identity} does not exist. warning in your call to Cases – M.R. Dec 11 '18 at 18:45
@M.R. Sorry, I cannot reproduce this. =/ But admittedly, the way I used MeshFunctions is not very robust. If Plot finds no intersections, the rest of the code will throw errors. – Henrik Schumacher Dec 11 '18 at 19:54
@HenrikSchumacher what if I only needed the normal of the intersection point, not the intercepting cell - is there a faster way? – M.R. Dec 11 '18 at 21:36
@M.R. Hm. I don't think so... – Henrik Schumacher Dec 11 '18 at 22:08
Why did you choose not to use RegionIntersection by the way? @HenrikSchumacher – M.R. Dec 11 '18 at 22:24
Because RegionIntersection can currently intersect only regions of the same embedding dimension. And I tried to use only methods that have complexity $O(\log(N))$ where $N$ is the number of triangles. – Henrik Schumacher Dec 11 '18 at 22:25
Nice, why not replace the line with a very thin polygon then? – M.R. Dec 11 '18 at 22:29
I tried. But RegionIntersection wants MeshRegions as input and nearly gegenerate triangles are eliminated when initializing a MeshRegion... – Henrik Schumacher Dec 11 '18 at 22:33

How to find the cells of a region that intersect a line?

1 Answers1

Linked