Draw a sphere using 3D pixels (voxels)

Question

Can you suggest an algorithm that can draw a sphere in 3D space using only the basic plot(x,y,z) primitive (which would draw a single voxel)?

I was hoping for something similar to Bresenham's circle algorithm, but for 3D instead of 2D.

FYI, I'm working on a hardware project that is a low-res 3D display using a 3-dimensional matrix of LEDs, so I need to actually draw a sphere, not just a 2D projection (i.e. circle).

The project is very similar to this:

... or see it in action here.

One possibility I have in mind is this:

• calculate the Y coordinates of the poles (given the radius) (for a sphere centered in the origin, these would be -r and +r)
• slice the sphere: for each horizontal plane p_i between these coordinates, calculate the radius of the circle obtained by intersecting said plane with the sphere => r_i.
• draw the actual circle of radius r_i on plane p_i using Bresenham's algorithm.

FWIW, I'm using a .NET micro-framework microprocessor, so programming is C#, but I don't need answers to be in C#.

Solution 1

The simple, brute force method is to loop over every voxel in the grid and calculate its distance from the sphere center. Then color the voxel if its distance is less than the sphere radius. You can save a lot of instructions by eliminating the square root and comparing the dot product to the radius squared.

Pretty far from optimal, sure. But on an 8x8x8 grid as shown, you'll need to do this operation 512 times per sphere. If the sphere center is on the grid, and its radius is an integer, you only need integer math. The dot product is 3 multiplies and 2 adds. Multiplies are slow; let's say they take 4 instructions each. Plus you need a comparison. Add in the loads and stores, let's say it costs 20 instructions per voxel. That's 10240 instructions per sphere.

An Arduino running at 16 MHz could push 1562 spheres per second. Unless you're doing tons of other math and I/O, this algorithm should be good enough.

Solution 2

I don't believe running the midpoint circle algorithm on each layer will give the desired results once you reach the poles, as you will have gaps in the surface where LEDs are not lit. This may give the result you want, however, so that would be up to aesthetics. This post is based on using the midpoint circle algorithm to determine the radius of the layers through the middle two vertical octants, and then when drawing each of those circles also setting the points for the polar octants.

I think based on @Nick Udall's comment and answer here using the circle algorithm to determine radius of your horizontal slice will work with a modification I proposed in a comment on his answer. The circle algorithm should be modified to take as an input an initial error, and also draw the additional points for the polar octants.

• Draw the standard circle algorithm points at y0 + y1 and y0 - y1: x0 +/- x, z0 +/- z, y0 +/- y1, x0 +/- z, z0 +/- x, y0 +/- y1, total 16 points. This forms the bulk of the vertical of the sphere.
• Additionally draw the points x0 +/- y1, z0 +/- x, y0 +/- z and x0 +/- x, z0 +/- y1, y0 +/- z, total 16 points, which will form the polar caps for the sphere.

By passing the outer algorithm's error into the circle algorithm, it will allow for sub-voxel adjustment of each layer's circle. Without passing the error into the inner algorithm, the equator of the circle will be approximated to a cylinder, and each approximated sphere face on the x, y, and z axes will form a square. With the error included, each face given a large enough radius will be approximated as a filled circle.

---

The following code is modified from Wikipedia's Midpoint circle algorithm. The DrawCircle algorithm has the nomenclature changed to operate in the xz-plane, addition of the third initial point y0, the y offset y1, and initial error error0. DrawSphere was modified from the same function to take the third initial point y0 and calls DrawCircle rather than DrawPixel

public static void DrawCircle(int x0, int y0, int z0, int y1, int radius, int error0)
{
  int x = radius, z = 0;
  int radiusError = error0; // Initial error state passed in, NOT 1-x

  while(x >= z)
  {
    // draw the 32 points here.
    z++;
    if(radiusError<0)
    {
      radiusError+=2*z+1;
    }
    else
    {
      x--;
      radiusError+=2*(z-x+1);
    }
  }
}

public static void DrawSphere(int x0, int y0, int z0, int radius)
{
  int x = radius, y = 0;
  int radiusError = 1-x;

  while(x >= y)
  {
    // pass in base point (x0,y0,z0), this algorithm's y as y1,
    // this algorithm's x as the radius, and pass along radius error.
    DrawCircle(x0, y0, z0, y, x, radiusError);
    y++;
    if(radiusError<0)
    {
      radiusError+=2*y+1;
    }
    else
    {
      x--;
      radiusError+=2*(y-x+1);
    }
  }
}

---

For a sphere of radius 4 (which actually requires 9x9x9), this would run three iterations of the DrawCircle routine, with the first drawing a typical radius 4 circle (three steps), the second drawing a radius 4 circle with initial error of 0 (also three steps), and then the third drawing a radius 3 circle with initial error 0 (also three steps). That ends up being nine calculated points, drawing 32 pixels each. That makes 32 (points per circle) x 3 (add or subtract operations per point) + 6 (add, subtract, shift operations per iteration) = 102 add, subtract, or shift operations per calculated point. In this example, that's 3 points for each circle = 306 operations per layer. The radius algorithm also adds 6 operations per layer and iterates 3 times, so 306 + 6 * 3 = 936 basic arithmetic operations for the example radius of 4. The cost here is that you will repeatedly set some pixels without additional condition checks (i.e. x = 0, y = 0, or z = 0), so if your I/O is slow you may be better off adding the condition checks. Assuming all LEDs were cleared at the start, the example circle would set 288 LEDs, while there are many fewer LEDs that would actually be lit due to repeat sets.

It looks like this would perform better than the bruteforce method for all spheres that would fit in the 8x8x8 grid, but the bruteforce method would have consistent timing regardless of radius, while this method will slow down when drawing large radius spheres where only part will be displayed. As the display cube increases in resolution, however, this algorithm timing will stay consistent while bruteforce will increase.

Solution 3

Assuming that you already have a plot function like you said:

    public static void DrawSphere(double r, int lats, int longs)
    {
        int i, j;
        for (i = 0; i <= lats; i++)
        {
            double lat0 = Math.PI * (-0.5 + (double)(i - 1) / lats);
            double z0 = Math.Sin(lat0) * r;
            double zr0 = Math.Cos(lat0) * r;

            double lat1 = Math.PI * (-0.5 + (double)i / lats);
            double z1 = Math.Sin(lat1) * r;
            double zr1 = Math.Cos(lat1) * r;

            for (j = 0; j <= longs; j++)
            {
                double lng = 2 * Math.PI * (double)(j - 1) / longs;
                double x = Math.Cos(lng);
                double y = Math.Sin(lng);

                plot(x * zr0, y * zr0, z0);
                plot(x * zr1, y * zr1, z1);
            }
        }
    }

That function should plot a sphere at the origin with specified latitude and longitude resolution (judging by your cube you probably want something around 40 or 50 as a rough guess). This algorithm doesn't "fill" the sphere though, so it will only provide an outline, but playing with the radius should let you fill the interior, probably with decreasing resolution of the lats and longs along the way.

Solution 4

Just found an old q&a about generating a Sphere Mesh, but the top answer actually gives you a short piece of pseudo-code to generate your X, Y and Z :

(x, y, z) = (sin(Pi * m/M) cos(2Pi * n/N), sin(Pi * m/M) sin(2Pi * n/N), cos(Pi * m/M))

Check this Q&A for details :) procedurally generate a sphere mesh

View more solutions

Comments

[Christian Rau]

Well, the problem is that he doesn't seem to want a solid sphere. In this case you got the classical rasterization problem of preventing holes and swellings on the surface of the sphere. Though he could still use your solution and make a second pass removing all interior voxels by checking their 6-neighbours.

[Christian Rau]

Oh right, you don't need a second pass for the removal of interior voxels, as the object is defined implicitly. But in this case you have to do the distance computation more than 1 time per voxel on average.

[Christian Rau]

Isn't this just NominSim's solution?

[japreiss]

Yeah, at first I thought it was filled in the demo video, but on closer inspection it looks like the shell only. I agree with your first comment, doing a morphological operation would be simplest.

[Cedric Mamo]

don't you need to multiply the coordinates by r in the plot calls? as it is r is unused, meaning it will always draw a sphere of radius 1.

[John T]

I don't understand where is the actual part to draw plot

[John T]

Why 2 plots? Whats the difference.