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

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 7:31

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[RegionMeshMeshNearestCellIndex[R, pts]]


{{2, 744}}

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


• @ 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 RegionMeshMeshNearestCellIndex: 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