Several functions in one RegionPlot with different grayscale

Question

I try to plot the following, but something seems to be wrong in my code and I can't figure out how to fix it:

 RegionPlot3D[{CountRoots[\[Sigma]q33[\[Omega], V1, V2, 
    V3], {\[Omega], 0, Infinity} == 0], 
  CountRoots[\[Sigma]q33[\[Omega], V1, V2, 
    V3], {\[Omega], 0, Infinity} == 1], 
  CountRoots[\[Sigma]q33[\[Omega], V1, V2, 
    V3], {\[Omega], 0, Infinity} == 2], 
  CountRoots[\[Sigma]q33[\[Omega], V1, V2, 
    V3], {\[Omega], 0, Infinity} == 3 ], 
  CountRoots[\[Sigma]q33[\[Omega], V1, V2, 
    V3], {\[Omega], 0, Infinity} == 4]}, {V1, 0, 50}, {V2, 0, 
  50}, {V3, 0, 100}, {GrayLevel[i, 0.5], {i, 0, 0.8, 0.2}}, 
 Mesh -> None, LabelStyle -> (FontSize -> 14)]

The number of roots in omega of the function \[Sigma]q33 depends on the other Variables V1, V2 and V3. I would like to visualize that. The regions where CountRoots gives different values should be in different gray levels, all with an opacity of 0.5. I think the problem I have is that I do not know how to use GrayLevel to define different shades of gray for the different areas.

And do you have an idea how to define the CountRoots-functions all in one, maybe with Map or something like that? Here is the \[Sigma]q33, in case it is necessary to answer the question:

 \[Sigma]q33[\[Omega]_, V1_, V2_,V3_] = (1/48 ( Exp[-V1] + Exp[-V2]) (1 + Exp[-V3])  + 1/48 (Exp[-(V2 - V1)/2] - Exp[(V2 - V1)/2])/(Exp[(V2 - V1)/2] + 
       Exp[-(V2 - V1)/2]) (( Exp[-V1] - Exp[-V2]  Exp[-V3]) - ( 
       Exp[-V2] - Exp[-V1]  Exp[-V3])) + 
   1/8 (( Exp[-V1] - Exp[-V2]  Exp[-V3]) + ( 
       Exp[-V2] - Exp[-V1]  Exp[-V3])) (-2 (  
       Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
        Exp[-V1]  Exp[-V3]) (-Exp[-(V1 + V2)/2 - V3] + 
          Exp[-(V2 + V1)/2])/(Exp[(V2 - V1)/2] + 
          Exp[-(V2 - V1)/2])/(\[Omega]^2 + (  
           Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
            Exp[-V1]  Exp[-V3]) ^2)) + 
   1/2 (( Exp[-V1] - Exp[-V2]  Exp[-V3]) - ( 
       Exp[-V2] - Exp[-V1]  Exp[-V3])) (-( 
        Exp[-V1] + Exp[-V2]  Exp[-V3] - Exp[-V2] - 
         Exp[-V1]  Exp[-V3]) (-Exp[-(V1 + V2)/2 - V3] + 
         Exp[-(V2 + V1)/2])/(Exp[(V2 - V1)/2] + 
         Exp[-(V2 - V1)/2]) (2 \[Omega]^2 - (  
           Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
            Exp[-V1]  Exp[-V3])^2)/((\[Omega]^2 + (  
             Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
              Exp[-V1]  Exp[-V3])^2) (4 \[Omega]^2 + (  
             Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
              Exp[-V1]  Exp[-V3])^2))) + (( 
       Exp[-V1] - Exp[-V2]  Exp[-V3]) + ( 
       Exp[-V2] - Exp[-V1]  Exp[-V3])) ((  
        Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
         Exp[-V1]  Exp[-V3]) (-Exp[-(V1 + V2)/2 - V3] + 
          Exp[-(V2 + V1)/2])/(Exp[(V2 - V1)/2] + 
          Exp[-(V2 - V1)/2]) ((  
            Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
             Exp[-V1]  Exp[-V3])^2 - 
          5 \[Omega]^2)/(6 (9 \[Omega]^2 + (  
              Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
               Exp[-V1]  Exp[-V3]) ^2) (\[Omega]^2 + (  
              Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
               Exp[-V1]  Exp[-V3])^2)) + ( 
         Exp[-V1] + Exp[-V2]  Exp[-V3] - Exp[-V2] - 
          Exp[-V1]  Exp[-V3]) ^2 ((  
          Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
           Exp[-V1]  Exp[-V3]) (-Exp[-(V1 + V2)/2 - V3] + 
            Exp[-(V2 + V1)/2])/(Exp[(V2 - V1)/2] + 
            Exp[-(V2 - V1)/2]) (11 \[Omega]^2 - (  
              Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
               Exp[-V1]  Exp[-V3])^2)/(2 (\[Omega]^2 + (  
                Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
                 Exp[-V1]  Exp[-V3])^2) (4 \[Omega]^2 + (  
                Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
                 Exp[-V1]  Exp[-V3])^2) (9 \[Omega]^2 + (  
                Exp[-V1] +  Exp[-V2]  Exp[-V3] +  Exp[-V2] +  
                 Exp[-V1]  Exp[-V3])^2)))))

Answer 1

ee[v1_?NumericQ,v2_?NumericQ,v3_?NumericQ]:= CountRoots[σq33[ω,v1,v2,v3],{ω,0,Infinity}];

Show[Table[RegionPlot3D[
   ee[v1, v2, v3] == i, {v1, 0, 50}, {v2, 0, 50}, {v3, 0, 100}, 
   Mesh -> None, LabelStyle -> (FontSize -> 14), 
   PlotStyle -> Directive[GrayLevel[ i/10 + 1/10, 3/10]], 
   PlotPoints -> 7, BoundaryStyle -> None], {i, 0, 4}]]

If you want a better plot, increment PlotPoints, but it is too CPU intensive for my poor machine.

Edit

here you have a much better image done by @R.M without specifying PlotPoints

Another one, made by @GustavoBandeira with PlotPoints->12

@cormullion made a fly-around movie of this (too big for uploading here), but here's a still:

Answer 2

I can't get this code to run here. Try this hyperbolic tetrahedron in cage form :

RegionPlot3D[ (x - 1)^2 + (y - 1)^2 + (z - 1)^2 >= 1 && x^2 - y^2 + z^2 >= 0 && 
               x^2 + y^2 - z^2 >= 0 && -x^2 + y^2 + z^2 >= 0,
              {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, 
              Axes -> False, Boxed -> False, PlotRange -> All, PlotPoints -> 90,
              ColorFunction -> "Pastel", MeshStyle -> {Tube[0.25/16]},
              PlotStyle -> {Opacity[0.0]}, ViewPoint -> {0, 0, 5},
              ImageSize -> {674, 501}, Mesh -> 10]