# Extracting surface from RegionPlot3D

As in this question, I'd like to extract a surface from RegionPlot3D. Unfortunately, the suggestion to use ContourPlot3D in that question's comments doesn't work for my actual problem, which is rather involved and has a very steep gradient in the objective function.

Here's a simpler example that demonstrates the same issues as my actual problem:

reg = RegionPlot3D[z < (1 + 0.1 y)/(x + 10^-5), {x, 0, 1}, {y, 0, 20}, {z, 0, 10}]


I'm 99% of the way there, but my solution has some artifacts. My idea was to extract the points as here:

pts = Flatten[Cases[Normal@reg, Line[pts_] :> pts, Infinity], 1];
ListPointPlot3D[pts]


then take the top of each stack of points:

top = GroupBy[pts, Part[#, 1 ;; 2] & -> Last, Max];
ListPointPlot3D[top]


The points look right, but when I try to get the surface with ListPlot3D there are weird blades I need to get rid of:

ListPlot3D[top, Mesh -> False, PlotStyle -> {Gray, Opacity[0.5]}]


Any ideas on how to either remove those artifacts in ListPlot3D, or otherwise extract the top surface from RegionPlot3D?

• How did you find that: "Unfortunately, the suggestion to use ContourPlot3D in that question's comments doesn't work for my actual problem"? I think it works perfectly well. Dec 4, 2022 at 18:16
• @azerbajdzan the example in my question is a simplified toy problem, because my actual problem is too complicated to post Dec 4, 2022 at 18:40
• ImplicitRegion[And[z < (1 + 0.1 y)/(x + 10^-5), 0 <= x <= 1, 0 <= y <= 20, 0 <= z <= 10 ], {x, y, z}] // BoundaryDiscretizeRegion // MeshPrimitives[#, 2] &//RegionUnion has dimension 2 as checked by RegionDimension. It does not look exactly like the plot in the image but maybe that is because of the scale. Dec 4, 2022 at 18:54
• I also tried to find an exact region given that you have equations using ImplicitRegion[ And[z < (1 + 0.1 y)/(x + 10^-5), 0 <= x <= 1, 0 <= y <= 20 , 0 <= z <= 10] // Rationalize[#, 0] & // CylindricalDecomposition[#, {x, y, z}, "Boundary"] &, {x, y, z}] but the region is missing pieces of the original region's surface. Dec 4, 2022 at 18:57
• @userrandrand All of your answers were very helpful, so it was hard to pick one. thanks! Dec 5, 2022 at 18:15

### Simple case (works here)

If no mesh of the surface is ever parallel to the x,y or z plane then one can utilize the VertexNormals option for light shading to obtain the normals of the meshes in the x,y or z plane. These normals look like {0,0,1}, {0,-1,0}, {0,1,0} etc but as floating point real numbers. Hence, we may remove them like this :

(Note the usage of Normal below to convert the GraphicsComplex structure into an ordinary list of graphics primitives and directives to facilitate pattern matching)

(p[[2]] below has the form VertexNormals-> {values__} )

Normal@reg /.
p_Polygon :>
Nothing /;
AllTrue[p[[2, 2]],
MatchQ[Abs@Rationalize[#], {OrderlessPatternSequence[0, 1, 0]}] &]


That works except maybe if one explicitly changes the VertexNormals option via NormalsFunction or something.

### More complicated scenario with piecewise constant surfaces (not the case here)

If the surface does have meshes that are parallel to the x, y or z plane then the direction of normals is not enough and one has to use the coordinates of these planes. Such a scenario is likely rare unless one has a piece wise function or a bad mesh but for completeness a code for such a scenario is included below.

The code below takes as argument a polygon poly and a couple {n,val} where n=1,2 or 3 for x,y or z and val is the constant value taken on the plane:

withinPlane[couple_][poly_] :=
AllTrue[Rationalize[poly[[1]],
0], #[[couple[[1]] ]] == couple[[2]] &];


Then one may remove the meshes belonging to planes at the boundaries of the plot:

Normal@reg /.
p_Polygon :>
Nothing /;
Or @@ Through@{withinPlane[{1, 0}], withinPlane[{1, 1}], ,
withinPlane[{2, 0}], withinPlane[{2, 20}], withinPlane[{3, 0}],
withinPlane[{3, 10}]}@p


Which leads to the same plot above.

Perhaps this:

With[{pts = First@Cases[reg, GraphicsComplex[p_, ___] :> p, Infinity]},
reg /. p_Polygon :>
Nothing /; !
FreeQ[Length /@ DeleteDuplicates /@ Transpose@pts[[p[[1, 1]]]], 1]
]


• Awesome, thanks! I figured that you used Mesh -> False to eliminate the "jungle gym" effect. Also, adding PlotRangePadding -> 0 seems to be an easy way to get rid of the rest of the frame, leaving just the surface floating there. Dec 5, 2022 at 0:57
• @ChrisK Actually I was in a hurry -- had to go somewhere -- and getting rid of Mesh seemed likely to simplify things. You can get rid of the boundary edges with BoundaryStyle -> None, if that's what you meant by "frame." Dec 5, 2022 at 2:11
• Perfect, also works on my actual problem ;) Dec 5, 2022 at 6:35

Edit

• We locate the level of the lastGraphicsGroup and use Delete to delete it's upper level by Most.
• And DeleteCases the Line.
indexs = Position[reg, GraphicsGroup];
reg1 = Delete[reg, indexs // Last // Most];
DeleteCases[reg1, _Line, -1]


Original

Another possible way.

• We add FaceForm[] to the last GraphicsGroup to erase the surface and replace Line to Nothing to erase all the lines.
Insert[reg, FaceForm[], Position[reg, GraphicsGroup[_]] // Last] /.
Line -> Nothing

• Or ReplacePart the last GraphicsGroup to Nothing.
ReplacePart[reg, GraphicsGroup -> Nothing,
Position[reg, GraphicsGroup[_]] // Last] /. Line -> Nothing

• Very nice this method even works with piecewise constant meshes like RegionPlot3D[If[0.2 < x < 0.7 && 8 < y < 16, z < 4, z < (1 + 0.1 y)/(x + 10^-5)], {x, 0, 1}, {y, 0, 20}, {z, 0, 10}]; Dec 5, 2022 at 9:48
• @userrandrand Thanks your example. Dec 5, 2022 at 12:26