# 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 May 25 '14 at 21:14

## 2 Answers

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]},
opts], {p, preds}]]


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. – Rahul May 25 '14 at 17:33
@RahulNarain, thanks. I added a few words about how it works. – Simon Woods May 25 '14 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. – Rahul May 25 '14 at 21:30
@RahulNarain, I think they should be equivalent, and I agree it's clearer that way. – Simon Woods May 25 '14 at 21:34

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)]) – Simon Woods May 24 '14 at 21:55
@Simon: That's neat! I wasn't aware of BooleanConvert. – Rahul May 24 '14 at 23:23