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#.

Answer 1

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.

Answer 2

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.