# High-quality RegionPlot3D for logical combinations of predicates

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:

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):

• Finally, thanks for the thread. :)
– Kuba
Commented May 25, 2014 at 21:14

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"

• Fantastic! Can you write a little about how it works? The Refine trick is very interesting and a little mysterious.
– user484
Commented May 25, 2014 at 17:33
• @RahulNarain, thanks. I added a few words about how it works. Commented May 25, 2014 at 21:12
• Thanks for the explanation. That's very ingenious! By the way, I would find the line easier to understand if written as Refine[reg, p] && Refine[! reg, ! p], unless there's some subtle difference I'm not seeing.
– user484
Commented May 25, 2014 at 21:30
• @RahulNarain, I think they should be equivalent, and I agree it's clearer that way. Commented May 25, 2014 at 21:34
• @j0equ1nn, the orange coloring is the new default style for 3D graphics (since version 10 I think). You can get the old coloring with the option PlotTheme -> "Classic". To get different colors for each bounding surface you could use something like ContourStyle :> RandomColor[] Commented Oct 5, 2017 at 13:45

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.

• Nice idea. Not really thought this through, but if you use f = BooleanConvert[LogicalExpand[#], "ESOP"] & to convert the predicate to an exclusive sum of products, I think you should split the region into distinct parts which you can use your code on. e.g. Show @@ (myRegionPlot3D[#, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}] & /@ f[1/4 <= x^2 + y^2 + z^2 <= 1 && (x <= 0 || y >= 0)]) Commented May 24, 2014 at 21:55
• @Simon: That's neat! I wasn't aware of BooleanConvert.
– user484
Commented May 24, 2014 at 23:23