High-quality RegionPlot3D for logical combinations of predicates

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]

enter image description here

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]

enter image description here

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]

enter image description here

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"


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

enter image description here

which the default RegionPlot3D is bad at:

enter image description here

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

Tags:

Regions

Plotting

Recent Posts