POV Ray 3.7 code
//===================================================================
// 3D quaternion Mandelbrot-Mandelbrot-Mandelbrot set generation
// by a cellular automaton (POV-Ray 3.7)
//===================================================================

#include "colors.inc"  // Standard Color definitions
#include "glass.inc"   //  Glass pigment effect
#include "metals.inc"   // Metal pigment effect
global_settings {assumed_gamma 1.0}
background { color rgb <00,0.0,0.0> }

camera {
    location <-20, 20, -300>
    look_at  <-5, 0,  0>
}

light_source {
	< -20, 20, -10>
	rgb <1.000000, 1.000000, 1.000000> * 2.0
}

light_source {
	< 20, 20, -10>
	rgb <1.000000, 1.000000, 1.000000> * 3.0
}

#declare R =  0.1;     // set radius value
#declare L = 2;  // set square area side
#declare n = 200;     // set number of pixels per area side
#declare st = 1;   // set increment step
#declare N = 40;    // set number of cycles
#declare Th = -45;   // set rotation angles
#declare Ph = 30;

#declare K = array[2*n+1][2*n+1][2*n+1];   // escape rate 3D matrix
  #for (p, -n, n, st)
  #declare IncX = p*L/n;   // Mand1 X increment
  #for (q, -n, n, st)
  #declare IncY = q*L/n;    // Mand1 Y1 increment
    #for (r, -n, n, st)
      #declare IncZ = r*L/n;  // Mand2 Y2 increment
      #declare X = 0;   // start MandX
      #declare Y = 0;   // start MandY
      #declare Z = 0;   // start MandZ
      #declare K[p+n][q+n][r+n] = 0;

     #for (k,0,N)
       #declare XX = X*X - Y*Y - Z*Z + IncX;  // cycle MandX
       #declare YY = 2*X*Y + IncY;  // cycle MandY
       #declare ZZ = 2*X*Z + IncZ;   // cycle MandZ
       #declare X = XX;
       #declare Y = YY;
       #declare Z = ZZ;
       #declare W = X*X +Y*Y + Z*Z;

caption(W < R)     // escape if   0.0015
     	 #declare K[p+n][q+n][r+n] = k;
#end  // end if
#end  // end for k
#end  // end for q
#end  // end for p
#end  // end for r

#declare Rnd_1 = seed (1153);
#declare Rnd_2 = seed (553);
#declare Rnd_3 = seed (876);
#declare pp = 0;
#declare qq = 0;
#declare rr = 0;

union{
#for(j,0,pow(n,3))    // partial values n*n/4, n*n, n*n*n/64, n*n*n/16, n*n*n
   #declare p = pp + st*pow(-1,int(n*rand(Rnd_1)));
   #declare q = qq + st*pow(-1,int(n*rand(Rnd_2)));
   #declare r = rr + st*pow(-1,int(n*rand(Rnd_3)));

   #if ( abs(p) < n)
   #if( abs(q) < n)
   #if (abs(r) < n)
   #if(K[p+n][q+n][r+n] > 0)
   #declare pp = p;
   #declare qq = q;
   #declare rr = r;
 	 	sphere {
  		  < p, q, r >, 1   //  adding 3d axis
     			 texture {
    			  pigment { color Col_Glass_Yellow }
      		 }
      		 finish { ambient rgb <0.3,0.1,0.1>
     			  diffuse .3
     			 reflection .3
     			 specular 1 }   // plot sphere

translate < 30, 15, 0 >
rotate < 0, Th, Ph > }
#end // end if
#end // end if
#end // end if
#end // end if