Intersection of a plane with surface of ListSurfacePlot3D

If I have a list of points in 3D space that are only roughly located on a surface, this surface can be visualized with ListSurfacePlot3D. How can I find the intersection of this approximate surface with a plane, that spans between two vectors u and v? And to continue this, how would I find the intersection of the resulting line with another plane that spans between two vectors m and n?

The dataset of points is available here: https://www.dropbox.com/s/rlj91jrh1bp4g2c/data.txt?dl=0

• Crossposted here. Sep 24 '19 at 16:23

Being the dataset in the form of a number of points on a surface, it is not a trivial problem to compute the intersection between such a surface and a plane.

One way can be to find the sets of points on the surface that are closer to the given plane more than a preset threshold. Here is a way to do this (note that I am not using ListSurfacePlot3D at all here):

distancePointFromPlane[point_, planeV1_, planeV2_,
planeTranslation_] :=
With[{orthogonalVec = #/Norm@# &@Cross[planeV1, planeV2]},
Norm@Dot[orthogonalVec, point - planeTranslation]
];
distancePointsFromPlane[points_, planeV1_, planeV2_,
planeTranslation_] :=
With[{orthogonalVec = #/Norm@# &@Cross[planeV1, planeV2]},
Norm@Dot[orthogonalVec, # - planeTranslation] & /@ points
];
findPointsCloseToPlane[points_, planeV1_, planeV2_, planeTranslation_,
Range@Length@points,
distancePointsFromPlane[points, planeV1, planeV2,
planeTranslation]
} // Select[#, #[] <= maxDistance &] & // Map@First;

planeV1 = {1, 1, 0};
planeV2 = {1, -1, 1};
planeTranslation = {0, 0, 0.7};
dataPoints =
Import["./data.txt"] // StringReplace[#, "\n" -> ""] & //
ToExpression;
With[{pts = dataPoints},
Graphics3D[{
{Red, PointSize@0.02, Point@pts[[
findPointsCloseToPlane[pts, planeV1, planeV2, planeTranslation,
0.01]
]]},
Point@pts,
Opacity@0.3, InfinitePlane[planeTranslation, {planeV1, planeV2}]
}]
]
Graphics3D[{Point@dataPoints}] The results are not great though. You will get better results if the points of the surface are denser. A quick and dirty way to do this is to use ListSurfacePlot3D to build an interpolating surface, and then retrieve the points from the generated graphics:

planeV1 = {1, 1, 0};
planeV2 = {1, -1, 1};
planeTranslation = {0, 0, 0.7};
dataPoints =
Import["./data.txt"] // StringReplace[#, "\n" -> ""] & //
ToExpression;
With[{pts = ListSurfacePlot3D[dataPoints][[1, 1]]},
With[{goodPoints =
findPointsCloseToPlane[pts, planeV1, planeV2, planeTranslation,
0.01]},
Graphics3D[{
{Red, PointSize@0.02, Point@pts[[goodPoints]]},
Point@pts,
Opacity@0.3, InfinitePlane[planeTranslation, {planeV1, planeV2}]
}]
]] There is still plenty of room to improve, but it already looks much better I think.

• that's not too bad.. what would be a less dirty way to increase the density of the points even more? Or can I iteratively feed the points from ListSurfacePlot3D to a new ListSurfacePlot3D and get an even denser set of points? Sep 24 '19 at 22:06
• @piiipmatz you could try playing with the value of the MaxPlotPoints option of ListSurfacePlot3D.
– glS
Sep 26 '19 at 16:00

Update: Using MeshFunctions to get the intersection of the plane and ListSurfacePlot3D surface:

We can use the points p1, v1 and v2 to get the equation of the plane ip:

Simplify[Cross[v1, v2].({x, y, z} - p1)]

7/5 + x - y - 2 z

and use it as the MeshFunctions option value to show the intersection of ListSurfacePlot3D and ip:

lsp = ListSurfacePlot3D[data, Axes -> False,
MeshFunctions -> {7/5 + # - #2 - 2 #3  &}, Mesh -> {{0.}},
MeshStyle -> Directive[Red, Thick],
MeshShading -> {Opacity[.5, Yellow], Opacity[.5, Blue]},
BoxRatios -> 1];

Show[lsp, Graphics3D[{Lighting -> "Neutral", Point @ data,
Red, Sphere[#, .01] & /@ pts, EdgeForm[], Opacity[.3, Green], ip}]] We can do the same for arbitrary function f that defines a plane through f == 0. For example,

f[x_, y_, z_, w_]:= x - y + 2 z - w; (* the plane x - y + 2 z  == w *)

Show[ListSurfacePlot3D[data,
Axes -> False,
MeshFunctions -> {f[#, #2, #3, 1] &}, Mesh -> {{0}},
MeshStyle -> Directive[Red, Thick],
MeshShading -> {Opacity[.5, Blue], Opacity[.5, Yellow]},
BoxRatios -> 1, BoundaryStyle -> None],
ListPointPlot3D[data],
ContourPlot3D[f[x, y, z, 1] == 0, {x, 0, 1}, {y, 0, 1}, {z, 0, 1} ,
BoundaryStyle -> None,
ContourStyle -> Opacity[.25, Purple],
Mesh -> None]] Using the same InfinitePlane as in glS's answer:

p1 = {0, 0, 7/10};
v1 = {1, 1, 0};
v2 = {1, -1, 1};
dist = .01;

ip = InfinitePlane[p1, {v1, v2}];

Points within dist of ip can be obtained using Select and RegionDistance:

pts = Select[data, RegionDistance[ip, #] <= dist &];

We can also use ip as the value for ClippingPlanes directive to style parts of the surface above and below the plane differently:

lsp = ListSurfacePlot3D[data, Axes -> False,
Mesh -> None, PlotStyle -> Opacity[.7, Yellow],  BoxRatios -> 1];

Graphics3D[{Lighting -> "Neutral", Point @ data,
Red, Sphere[#, .01] & /@ pts, lsp[],
ClipPlanes -> ip, ClipPlanesStyle -> Opacity[.3, Green],
Opacity[.7, Blue], lsp[] /. Directive[__] :> {}}] • How does Graphics3D determine the parts above and below? This must include determining the line at the intersection of the green plane with the surface that is interpolated from the points? Sep 24 '19 at 13:59
• @piiipmatz, ClipPlanes >> Details says "ClipPlanes -> {{a1,b1,c1,d1}, {a2,b2,c2,d2},...} specifies a list of clipping planes defined as ai x + bi y + ci z + di <=0"
– kglr
Sep 24 '19 at 14:30
• my question was more related to the actual algorithm behind the scene! I.e. it is nice if the plot function can do all that stuff, but what if I wanted to use the intersection for further calculations? Using MeshFunctions you get this very nice red curve in the plot but what are the coordinates of the points on this curve? Sep 24 '19 at 21:56
• @piiipmatz, if you have the patience/perserverence you can inspect the output of ClearAttributes[RegionInfinitePlaneCuboidClip, {Protected, ReadProtected}]; Graphics3D[{Sphere[]}, ClipPlanes -> InfinitePlane[{{0, 0, 0}, {0, 1, 1}, {1, 1, 2}}]]; ??RegionInfinitePlaneCuboidClip for an example of what goes on "behind the scene".
– kglr
Sep 24 '19 at 22:19
• @piiipmatz, re "what are the coordinates of the points on this curve", you can try lsp=ListSurfacePlot3D[data, Axes -> False, MeshFunctions -> {f[#, #2, #3, 1] &}, Mesh -> {{0}}, MeshStyle -> Directive[Red, Thick], MeshShading -> {Opacity[.5, Blue], Opacity[.5, Yellow]}, BoxRatios -> 1, BoundaryStyle -> None]; pointsoncurve=Cases[Normal[lsp], Line[x]:>x, All].
– kglr
Sep 24 '19 at 22:21