I finally got around to cleaning my up my code a bit so I could share it. The examples are all in C#, but with the exception of the Animated GIF library I used it should be extremely simple to port into any other language.

I'm pretty sure that providing any code for the Complex class is overkill, but here's an abridged version anyway:

`public class Complex`

{

public double real;

public double imag;

public Complex()

{

this.real = 0;

this.imag = 0;

}

public Complex(double real, double imag)

{

this.real = real;

this.imag = imag;

}

public void SetValue(double real, double imag)

{

this.real = real;

this.imag = imag;

}

}

The Voxel struct is extremely simple. It's basically just a float that acts as our bucket and a bool value that flags whether or not the bucket counted when rendering.

`public struct Voxel`

{

public bool active;

public float value;

public Voxel(bool active)

{

this.active = active;

this.value = 0;

}

}

Here we init all of the Voxels and carve the sphere out of a 3D array. The diameter should be equal to the width and height of the image you want to render. Setting the "active" bool to true tells the program whether or not the Voxel is part of the sphere. This serves two functions:

1. We can save some time during rendering by ignoring the "false" voxels. (in retrospect we could also check for "null" if we don't initialize them. I'm not sure which is faster in C#.)

2. The algorithm that fills the buckets is still a little buggy and will draw outside the sphere. We can clean it up by disabling rendering on all points outside of the sphere.

`public Voxel[,,] CreateVoxelSphere(int diameter)`

{

Voxel[,,] voxelSphere = new Voxel[diameter, diameter, diameter];

//Diameter = 300

// 0 -> -150

//299 -> 150

//c = im - 150

//300 = 299m - 150

//(300 - 150) / (300 - 1)

float radius = (float)diameter / 2;

float m = diameter / (diameter - 1);

//Calculate the coordinates of the voxel.

for (int x = 0; x < diameter; x++)

{

float xm = (x * m) - radius;

for (int y = 0; y < diameter; y++)

{

float ym = (y * m) - radius;

for (int z = 0; z < diameter; z++)

{

float zm = (z * m) - radius;

float distance = (float)Math.Sqrt((xm * xm) + (ym * ym) + (zm * zm));

if (distance <= radius)

{

voxelSphere[x, y, z] = new Voxel(true);

}

else

{

voxelSphere[x, y, z] = new Voxel(false);

}

//Console.WriteLine(xm + ", " + ym + ", " + zm);

if (z == 0 && y == 0)

{

Console.WriteLine(x + ", " + y + ", " + z);

}

}

}

}

return voxelSphere;

}

Point3D is a struct that is part of System.Windows.Media.Media3D, but all I use it for is to store three doubles that represent 3D space.

`public struct Point3D`

{

public double X;

public double Y;

public double Z;

public Point3D()

{

this.X = 0;

this.Y = 0;

this.Z = 0;

}

public Point3D(double x, double y, double z)

{

this.X = x;

this.Y = y;

this.Z = z;

}

}

This converts a Complex number into a 3D point on the sphere.

`static Point3D ConvertComplexToRiemannCoordinates(Complex input)`

{

double realSquare = input.real * input.real;

double imagSquare = input.imag * input.imag;

double x = (2 * input.real) / (1 + realSquare + imagSquare);

double y = (2 * input.imag) / (1 + realSquare + imagSquare);

double z = (realSquare + imagSquare - 1) / (1 + realSquare + imagSquare);

return new Point3D(x, y, z);

}

This is the slightly buggy algorithm that will fill the sphere. At the end of each iteration in your fractal generator, convert \( Z_n \) and \( Z_{n-1} \) into Point3D's using the above function and send them to this function along with the Voxel Sphere.

`static void RiemannVoxelTravelAlgorithm(Point3D start, Point3D end, Voxel[,,] voxelSphere)`

{

//Rise Over Run

double xyRateOfChange = (end.Y - start.Y) / (end.X - start.X);

//b = y - mx

double xy_yIntercept = start.Y - (xyRateOfChange * start.X);

double yzRateOfChange = (end.Z - start.Z) / (end.Y - start.Y);

double yz_zIntercept = start.Z - (yzRateOfChange * start.Y);

int directionX;

int directionY;

int directionZ;

int nextIntX;

int nextIntY;

int nextIntZ;

if (start.X < end.X)

{

directionX = 1;

nextIntX = (int)Math.Ceiling(start.X);

}

else

{

directionX = -1;

nextIntX = (int)Math.Floor(start.X);

}

if (start.Y < end.Y)

{

directionY = 1;

}

else

{

directionY = -1;

}

if (start.Z < end.Z)

{

directionZ = 1;

}

else

{

directionZ = -1;

}

Point3D currentSubX = new Point3D(start.X, start.Y, start.Z);

Point3D nextSubX = new Point3D(end.X, end.Y, end.Z);

Point3D currentSubY = new Point3D();

Point3D nextSubY = new Point3D();

Point3D currentSubZ = new Point3D();

Point3D nextSubZ = new Point3D();

int numberOfXSteps = (int)Math.Abs(Math.Floor(start.X) - Math.Floor(end.X)) + 1;

Console.Write(numberOfXSteps + " ");

for (int x = 1; x <= numberOfXSteps && x < voxelSphere.GetLength(0); x++)

{

//Is this the last step?

if (x != numberOfXSteps)

{

//If not, calculate the next point.

//Solve for y at the given x.

//y = mx + b

double subX_y = (xyRateOfChange * nextIntX) + xy_yIntercept;

//Solve for z at the given y.

//z = my + b

double subX_z = (yzRateOfChange * subX_y) + yz_zIntercept;

nextSubX.X = nextIntX;

nextSubX.Y = subX_y;

nextSubX.Z = subX_z;

}

else

{

//Otherwise, set the next point to the value of the end point.

nextSubY.X = end.X;

nextSubY.Y = end.Y;

nextSubY.Z = end.Z;

}

//Calculate the number of Y sub-steps.

//A sub-step occurs when the line passes an integer threshold on the y-axis.

int numberOfYSubSteps = (int)Math.Abs(Math.Floor(currentSubX.Y) - Math.Floor(nextSubX.Y)) + 1;

currentSubY = new Point3D(currentSubX.X, currentSubX.Y, currentSubX.Z);

int currentIntY;

if (start.Y < end.Y)

{

currentIntY = (int)Math.Ceiling(currentSubY.Y);

}

else

{

currentIntY = (int)Math.Floor(currentSubY.Y);

}

nextIntY = currentIntY - directionY;

for (int y = 1; y <= numberOfYSubSteps; y++)

{

if (y != numberOfYSubSteps)

{

nextIntY += directionY;

nextSubY.Y = nextIntY;

//Calculate the next sub-point.

//Solve for x at the given y.

//y = mx + b

//y - b = mx

//(y - b) / m = x

//x = (y - b) / m

nextSubY.X = (nextIntY - xy_yIntercept) / xyRateOfChange;

// Solve for z at the given y.

//z = my + b

nextSubY.Z = (yzRateOfChange * nextSubY.Y) + yz_zIntercept;

}

else

{

//Otherwise, set the next point to the value of the end point.

nextSubY.X = nextSubX.X;

nextSubY.Y = nextSubX.Y;

nextSubY.Z = nextSubX.Z;

}

//Calculate the number of Y sub-steps.

//A sub-step occurs when the line passes an integer threshold on the y-axis.

int numberOfZSubSteps = (int)Math.Abs(Math.Floor(currentSubY.Z) - Math.Floor(nextSubY.Z)) + 1;

currentSubZ = new Point3D(currentSubY.X, currentSubY.Y, currentSubY.Z);

int currentIntZ;

if (start.Z < end.Z)

{

currentIntZ = (int)Math.Ceiling(currentSubZ.Z);

}

else

{

currentIntZ = (int)Math.Floor(currentSubZ.Z);

}

nextIntZ = currentIntZ - directionZ;

for (int z = 1; z <= numberOfZSubSteps; z++)

{

if (z != numberOfZSubSteps)

{

nextIntZ += directionZ;

nextSubZ.Z = nextIntZ;

//Calculate the next sub point.

//yz

//Solve for y at the given z.

//z = my + b

//z - b = my

//(z - b) / m = y

//y = (z - b) / m

nextSubZ.Y = (nextIntZ - yz_zIntercept) / yzRateOfChange;

//Solve for x at the current y.

//y = mx + b

//y - b = mx

//(y - b) / m = x

//x = (y - b) / m

nextSubZ.X = (nextSubZ.Y - xy_yIntercept) / xyRateOfChange;

}

else

{

//Otherwise, set the next point to the value of the end point.

nextSubZ.X = nextSubY.X;

nextSubZ.Y = nextSubY.Y;

nextSubZ.Z = nextSubY.Z;

}

int bucketX = (int)Math.Floor(nextSubZ.X);

int bucketY = (int)Math.Floor(nextSubZ.Y);

int bucketZ = (int)Math.Floor(nextSubZ.Z);

if (bucketX >= 0 && bucketX < voxelSphere.GetLength(0) && bucketY >= 0 && bucketY < voxelSphere.GetLength(1) && bucketZ >= 0 && bucketZ < voxelSphere.GetLength(2))

{

//Perform heat map operations here.

//Calculate the distance between currentSub and nextSub.

double deltaX = currentSubZ.X - nextSubZ.X;

double deltaY = currentSubZ.Y - nextSubZ.Y;

double deltaZ = currentSubZ.Z - nextSubZ.Z;

double delta = Math.Sqrt((deltaX * deltaX) + (deltaY * deltaY) + (deltaZ * deltaZ));

voxelSphere[bucketX, bucketY, bucketZ].value += (float)(delta / root3);

}

else

{

//For Debugging

//Console.WriteLine(bucketX + ", " + bucketY + ", " + bucketZ);

}

currentSubZ.X = nextSubZ.X;

currentSubZ.Y = nextSubZ.Y;

currentSubZ.Z = nextSubZ.Z;

}

currentSubY.X = nextSubY.X;

currentSubY.Y = nextSubY.Y;

currentSubY.Z = nextSubY.Z;

}

currentSubX.X = nextSubX.X;

currentSubX.Y = nextSubX.Y;

currentSubX.Z = nextSubX.Z;

nextIntX += directionX;

}

}

This is a really simple function that just finds the maximum value in the given array. It's used in the rendering method that's coming up next.

`private double FindMaxDoubleIn2DArray(double[,] doubleArray)`

{

double maxDouble = double.MinValue;

for (int x = 0; x < doubleArray.GetLength(0); x++)

{

for (int y = 0; y < doubleArray.GetLength(1); y++)

{

if (doubleArray[x, y] > maxDouble)

{

maxDouble = doubleArray[x, y];

}

}

}

return maxDouble;

}

This is the code for the RGB version of animation renderer. Send it three spheres, each generated from the same function with different maximum iterations.

In this code Bitmap is a C# image class and AnimatedGIF is an external library I used. If you're porting this to another language you'll need to find replacements.

`static void VoxelSphereAnimationRGB(string folder, string fractalName, Voxel[,,] voxelSphereR, Voxel[,,] voxelSphereG, Voxel[,,] voxelSphereB, int numberOfFrames, double xAngleMultiplier = 1, double yAngleMultiplier = 0, double xAngleOffset = 0, double yAngleOffset = 0)`

{

double angleRateOfChange = (2 * Math.PI) / numberOfFrames;

double currentAngle = 0;

double currentAngleX = currentAngle * xAngleMultiplier + xAngleOffset;

double currentAngleY = currentAngle * xAngleMultiplier + xAngleOffset;

double sinAngleX = Math.Sin(currentAngleX);

double cosAngleX = Math.Cos(currentAngleX);

double sinAngleY = Math.Sin(currentAngleY);

double cosAngleY = Math.Cos(currentAngleY);

string folderRotation = folder + @"\Rotation";

Directory.CreateDirectory(folderRotation);

string gifFileName = folder + @"\Transparent Rotation - " + fractalName + ".gif";

int diameter = voxelSphereR.GetLength(0);

int radius = diameter / 2;

Bitmap bitmap = new Bitmap(diameter, diameter, PixelFormat.Format24bppRgb);

int valueR = 0;

int valueG = 0;

int valueB = 0;

double maxValueR = 0;

double maxValueG = 0;

double maxValueB = 0;

double logBaseR = 0;

double logBaseG = 0;

double logBaseB = 0;

double[,] combinedValuesR = new double[diameter, diameter];

double[,] combinedValuesG = new double[diameter, diameter];

double[,] combinedValuesB = new double[diameter, diameter];

Complex center = new Complex((diameter - 1) / 2, (diameter - 1) / 2);

Complex point = new Complex(); //Used to store some (x, y) coordinates.

Complex newPoint = new Complex(); //Used to store some (x, y) coordinates.

int newX = 0;

int newY = 0;

//30 fps

using (var gifWriter = AnimatedGif.AnimatedGif.Create(gifFileName, 30))

{

for (int i = 0; i < numberOfFrames; i++)

{

//Reset the output maps.

for (int y = 0; y < diameter; y++)

{

for (int x = 0; x < diameter; x++)

{

combinedValuesR[x, y] = 0;

combinedValuesG[x, y] = 0;

combinedValuesB[x, y] = 0;

}

}

//Rotate the sphere points and update the flattened output map.

for (int z = 0; z < diameter; z++)

{

for (int y = 0; y < diameter; y++)

{

for (int x = 0; x < diameter; x++)

{

if (voxelSphereR[x, y, z].active)

{

//Calculate the new x coordinate by rotating the x,z coordinates about the y-axis

point.setValue(x, z);

newPoint = point.rotateAboutPoint(center, sinAngleX, cosAngleX);

newX = (int)Math.Floor(newPoint.real);

point.setValue(newPoint.imag, y);

newPoint = point.rotateAboutPoint(center, sinAngleY, cosAngleY);

newY = (int)Math.Floor(newPoint.imag);

if (newX >= 0 && newX < diameter && newY >= 0 && newY < diameter)

{

combinedValuesR[newX, newY] += voxelSphereR[x, y, z].value;

combinedValuesG[newX, newY] += voxelSphereG[x, y, z].value;

combinedValuesB[newX, newY] += voxelSphereB[x, y, z].value;

}

}

}

}

if (z % 5 == 0)

{

Console.WriteLine(i + " rotation - " + z);

}

}

//Calculate the log value for the Red output map.

maxValueR = FindMaxDoubleIn2DArray(combinedValuesR);

logBaseR = Math.Pow(Math.E, Math.Log(maxValueR) / 255);

//Calculate the log value for the Red output map.

maxValueG = FindMaxDoubleIn2DArray(combinedValuesG);

logBaseG = Math.Pow(Math.E, Math.Log(maxValueG) / 255);

//Calculate the log value for the Red output map.

maxValueB = FindMaxDoubleIn2DArray(combinedValuesB);

logBaseB = Math.Pow(Math.E, Math.Log(maxValueB) / 255);

//XY

//Create the image from the 2D Arrays.

for (int y = 0; y < diameter; y++)

{

for (int x = 0; x < diameter; x++)

{

if (combinedValuesR[x, y] != 0)

{

valueR = (int)(Math.Log(combinedValuesR[x, y]) / Math.Log(logBaseR));

if (valueR > 255)

{

valueR = 255;

}

else if (valueR < 0)

{

valueR = 0;

}

}

else

{

valueR = 0;

}

if (combinedValuesG[x, y] != 0)

{

valueG = (int)(Math.Log(combinedValuesG[x, y]) / Math.Log(logBaseG));

if (valueG > 255)

{

valueG = 255;

}

else if (valueG < 0)

{

valueG = 0;

}

}

else

{

valueG = 0;

}

if (combinedValuesB[x, y] != 0)

{

valueB = (int)(Math.Log(combinedValuesB[x, y]) / Math.Log(logBaseB));

if (valueB > 255)

{

valueB = 255;

}

else if (valueB < 0)

{

valueB = 0;

}

}

else

{

valueB = 0;

}

bitmap.SetPixel(x, y, Color.FromArgb(valueR, valueG, valueB));

}

}

bitmap.Save(folder + @"\Rotation\" + i + ".bmp");

gifWriter.AddFrame(bitmap, delay: -1, quality: GifQuality.Bit8);

currentAngle += angleRateOfChange;

currentAngleX = currentAngle * xAngleMultiplier + xAngleOffset;

currentAngleY = currentAngle * xAngleMultiplier + xAngleOffset;

sinAngleX = Math.Sin(currentAngleX);

cosAngleX = Math.Cos(currentAngleX);

sinAngleY = Math.Sin(currentAngleY);

cosAngleY = Math.Cos(currentAngleY);

}

}

}

That's pretty much everything. After work today I want to draw up a diagram of how the algorithm is *supposed* to work. I've written it down a couple of times in this thread but I think a visual diagram would be more helpful.