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.

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

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
11
  • $\begingroup$ @ HenrikSchumacher Very helpful! Where could I find further documentation to RegionMeshMeshNearestCellIndex? Thanks. $\endgroup$
    – Ulrich Neumann
    Commented Dec 11, 2018 at 15:55
  • $\begingroup$ @Ulrich As always with these undocumented functions such as Region`Mesh`MeshNearestCellIndex: Nowhere. Unfortunately. =/ $\endgroup$
    – Henrik Schumacher
    Commented Dec 11, 2018 at 16:03
  • $\begingroup$ It's a pity, thanks $\endgroup$
    – Ulrich Neumann
    Commented Dec 11, 2018 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.
    Commented Dec 11, 2018 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
    Commented Dec 11, 2018 at 19:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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