High-quality RegionPlot3D for logical combinations of predicates

Question

Unlike RegionPlot, RegionPlot3D copes poorly with logical combinations of predicates (&&, ||), which should result in sharp edges in the region to be plotted. Instead, these edges are usually drawn rounded and sometimes with severe aliasing artifacts. This has been observed in many posts on this site:

• Problem in RegionPlot3D
• on RegionPlot3D and ContourPlot3D
• Tailoring RegionPlot3D with PlotPoints?
• "RegionPlot3D plots always seem to have that 'home-made' look about them..."
• "...the limitations of RegionPlot3D make the edges appear jagged."

One solution, as noted by Silvia, by halirutan, and most recently by Jens, is to use a ContourPlot3D instead with an appropriate RegionFunction, as this produces much higher-quality results. I think it would be useful to have a general-purpose solution along these lines. That is, we want a single function that can be used as a drop-in replacement for RegionPlot3D and will automatically produce high-quality results by setting up the appropriate instances of ContourPlot3D.

Here is a test example, inspired by this post:

RegionPlot3D[1/4 <= x^2 + y^2 + z^2 <= 1 && (x <= 0 || y >= 0),
 {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]

It should look more like this (created by increasing PlotPoints, and even then the edges are not perfectly sharp):

Answer 1

This is based on Rahul's ideas, but a different implementation:

contourRegionPlot3D[region_, {x_, x0_, x1_}, {y_, y0_, y1_}, {z_, z0_, z1_}, 
  opts : OptionsPattern[]] := Module[{reg, preds},
  reg = LogicalExpand[region && x0 <= x <= x1 && y0 <= y <= y1 && z0 <= z <= z1];
  preds = Union@Cases[reg, _Greater | _GreaterEqual | _Less | _LessEqual, -1];
  Show @ Table[ContourPlot3D[
     Evaluate[Equal @@ p], {x, x0, x1}, {y, y0, y1}, {z, z0, z1}, 
     RegionFunction -> Function @@ {{x, y, z}, Refine[reg, p] && Refine[! reg, ! p]},
     RegionBoundaryStyle -> None, opts], {p, preds}]]
 

(update - added RegionBoundaryStyle -> None which is required in v12.2)

Examples:

contourRegionPlot3D[
 (x < 0 || y > 0) && 0.5 <= x^2 + y^2 + z^2 <= 0.99,
 {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, Mesh -> None]

contourRegionPlot3D[
 x^2 + y^2 + z^2 <= 2 && x^2 + y^2 <= (z - 1)^2 && Abs@x >= 1/4,
 {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, Mesh -> None]

contourRegionPlot3D[
 x^2 + y^2 + z^2 <= 0.4 || 0.01 <= x^2 + y^2 <= 0.05,
 {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, Mesh -> None, PlotPoints -> 50]

How it works

Firstly LogicalExpand is used to split up multiple inequalities into a combination of single inequalities, for example to convert 0 < x < 1 into 0 < x && x < 1.

Like Rahul's code, an inequality like x < 1 is converted to the equality x == 1 to define a part of the surface enclosing the region. We do not generally want the entire x == 1 plane though, only that part for which the true/false value of the region function is determined solely by the true/false value of x < 1.

This is done by plotting the surface with a RegionFunction like this:

Refine[reg, p] && Refine[! reg, ! p]

which is equivalant to the predicate "reg is true when p is true, and reg is false when p is false"

Answer 2

This is not a complete answer. It only handles cases where each predicate is a single inequality, and they are combined only by And. So, for example, 0 <= x <= 1 will have to be rewritten as 0 <= x && x <= 1, while for x <= 0 || y >= 0 I don't have a solution yet. A complete general solution is still desired.

Suppose the combined predicate is $p = p_1 \wedge p_2 \wedge \cdots \wedge p_n$. Now $p$ can change from true to false if and only if one of the $p_i$ does while all the other $p_j$ are true. So the boundary of the region defined by $p$ consists of several pieces, namely the boundaries of each $p_i$ restricted to the region $\bigwedge_{j\ne i}p_j$. If $p_i$ is, say, $f\le g$, then said boundary is the contour $f=g$ in the specified region.

And so:

Options[myRegionPlot3D] = Options[ContourPlot3D];
myPatchPlot3D[eq_, preds_, {x_, x0_, x1_}, {y_, y0_, y1_}, {z_, z0_, z1_}, 
  opts : OptionsPattern[]] := 
 ContourPlot3D[eq, {x, x0, x1}, {y, y0, y1}, {z, z0, z1}, 
  RegionFunction -> Function[{x, y, z}, And @@ preds], opts]
myRegionPlot3D[pred_, {x_, x0_, x1_}, {y_, y0_, y1_}, {z_, z0_, z1_}, 
  opts : OptionsPattern[]] := 
 With[{preds = List @@ (pred && x0 <= x && x <= x1 && y0 <= y && y <= y1 && z0 <= z && z <= z1)}, 
  Show @@ Table[
    myPatchPlot3D[Equal @@ preds[[i]], Delete[preds, i], 
     {x, x0, x1}, {y, y0, y1}, {z, z0, z1}, opts], {i, Length@preds}]]

It can't do the example in the question, but it can do this:

myRegionPlot3D[x^2 + y^2 + z^2 <= 2 && x^2 + y^2 <= (z - 1)^2 && Abs@x >= 1/4,
 {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]

which the default RegionPlot3D is bad at:

Another limitation is that the mesh lines don't line up.