Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

enter image description here

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

enter image description here

share|improve this question
    
Finally, thanks for the thread. :) –  Kuba May 25 at 21:14

2 Answers 2

up vote 17 down vote accepted

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]

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"

share|improve this answer
    
Fantastic! Can you write a little about how it works? The Refine trick is very interesting and a little mysterious. –  Rahul Narain May 25 at 17:33
    
@RahulNarain, thanks. I added a few words about how it works. –  Simon Woods May 25 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 Narain May 25 at 21:30
    
@RahulNarain, I think they should be equivalent, and I agree it's clearer that way. –  Simon Woods May 25 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}]

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.

share|improve this answer
2  
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 at 21:55
    
@Simon: That's neat! I wasn't aware of BooleanConvert. –  Rahul Narain May 24 at 23:23

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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