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

enter image description here

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

enter image description here

then take the top of each stack of points:

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

enter image description here

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

enter image description here

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

  • 1
    $\begingroup$ 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. $\endgroup$ Dec 4, 2022 at 18:16
  • $\begingroup$ @azerbajdzan the example in my question is a simplified toy problem, because my actual problem is too complicated to post $\endgroup$
    – Chris K
    Dec 4, 2022 at 18:40
  • $\begingroup$ 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. $\endgroup$ Dec 4, 2022 at 18:54
  • $\begingroup$ 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. $\endgroup$ Dec 4, 2022 at 18:57
  • 1
    $\begingroup$ @userrandrand All of your answers were very helpful, so it was hard to pick one. thanks! $\endgroup$
    – Chris K
    Dec 5, 2022 at 18:15

3 Answers 3


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

enter image description here

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_] := 
    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]

enter image description here

  • $\begingroup$ 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. $\endgroup$
    – Chris K
    Dec 5, 2022 at 0:57
  • $\begingroup$ @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." $\endgroup$
    – Michael E2
    Dec 5, 2022 at 2:11
  • $\begingroup$ Perfect, also works on my actual problem ;) $\endgroup$
    – Chris K
    Dec 5, 2022 at 6:35


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

enter image description here


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
  • 2
    $\begingroup$ 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}]; $\endgroup$ Dec 5, 2022 at 9:48
  • 1
    $\begingroup$ @userrandrand Thanks your example. $\endgroup$
    – cvgmt
    Dec 5, 2022 at 12:26

Your Answer

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

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