Extracting surface from RegionPlot3D

Question

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?

Answer 1

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.

Answer 2

Perhaps this:

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

Answer 3

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