3
$\begingroup$

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]

enter image description here

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]

enter image description here

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

$\endgroup$
  • 1
    $\begingroup$ What version are you using?, my 11.3 says that BDR is not a graphics3d primitive. $\endgroup$ – Kuba Dec 11 '18 at 7:31
2
$\begingroup$

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

enter image description here

$\endgroup$
  • $\begingroup$ @ HenrikSchumacher Very helpful! Where could I find further documentation to RegionMeshMeshNearestCellIndex? Thanks. $\endgroup$ – Ulrich Neumann Dec 11 '18 at 15:55
  • $\begingroup$ @Ulrich As always with these undocumented functions such as Region`Mesh`MeshNearestCellIndex: Nowhere. Unfortunately. =/ $\endgroup$ – Henrik Schumacher Dec 11 '18 at 16:03
  • $\begingroup$ It's a pity, thanks $\endgroup$ – Ulrich Neumann Dec 11 '18 at 16:19
  • $\begingroup$ @HenrikSchumacher I'm seeing Part::partw: Part 3 of {Identity,Identity} does not exist. warning in your call to Cases $\endgroup$ – M.R. Dec 11 '18 at 18:45
  • $\begingroup$ @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. $\endgroup$ – Henrik Schumacher Dec 11 '18 at 19:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.