Drawing 3D sphere in C/C++

Question

I am looking for an algorithm which can draw a nice looking 3D sphere on small resolution. I found Bresenham's circle algorithm but it's for 2D drawing. I just need spheres borders (I don't need it filled). I also googled for a solution of the problem but I didn't find anything. This article doesn't help (what is the brute force algorithm?). I can't use any OpenGL libraries, I need vanilla C/C++ solution. Thank you in advance.

Answer 1

if I get it right you want to render all surface Voxels of sphere

The brute force is O(R^3). If you just project rays from plane and compute the 3-th coordinate then you get O(R^2) but to make sure that no Voxels are missing you have to do this projection from all 3 planes which is still O(R^2)

It look like this:

on LED cube 16x16x16 simulation. Now the algorithm:

1. compute visible bounding box

   no need to render whole rendering space just the sphere so center +/- radius...

2. take one plane (XY for example)

   Cast rays from all x,y points inside bounding box which is just 2 for loops and compute the z coordinates where the ray hits via sphere equation:

   (x-x0)^2 + (y-y0)^2 + (z-z0)^2 = R^2
   

   so

   z=z0 +/- sqrt(R^2 - (x-x0)^2 - (y-y0)^2)
   

   and render the two voxels. The int sqrt(int x) for limited size (like LED Cube/Screen or Voxel space) can be done via LUT lookup table to speed things up.

3. do the step #2 for all planes (xy,yz,xz)

The code in C++ looks like this:

//---------------------------------------------------------------------------
//--- LED cube class ver: 1.00 ----------------------------------------------
//---------------------------------------------------------------------------
#ifndef _LED_cube_h
#define _LED_cube_h
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
const int _LED_cube_size=16;
//---------------------------------------------------------------------------
class LED_cube
    {
public:
    int n,map[_LED_cube_size][_LED_cube_size][_LED_cube_size];

    LED_cube()              { n=_LED_cube_size; }
    LED_cube(LED_cube& a)   { *this=a; }
    ~LED_cube()             { }
    LED_cube* operator = (const LED_cube *a) { *this=*a; return this; }
    //LED_cube* operator = (const LED_cube &a) { /*...copy...*/ return this; }
    void cls(int col);                                  // clear cube with col 0x00BBGGRR
    void sphere(int x0,int y0,int z0,int r,int col);    // draws sphere surface with col 0x00BBGGRR
    void glDraw();                                      // render cube by OpenGL as 1x1x1 cube at 0,0,0
    };
//---------------------------------------------------------------------------
void LED_cube::cls(int col)
    {
    int x,y,z;
    for (x=0;x<n;x++)
     for (y=0;y<n;y++)
      for (z=0;z<n;z++)
       map[x][y][z]=col;
    }
//---------------------------------------------------------------------------
void LED_cube::sphere(int x0,int y0,int z0,int r,int col)
    {
    int x,y,z,xa,ya,za,xb,yb,zb,xr,yr,zr,xx,yy,zz,rr=r*r;
    // bounding box
    xa=x0-r; if (xa<0) xa=0; xb=x0+r; if (xb>n) xb=n;
    ya=y0-r; if (ya<0) ya=0; yb=y0+r; if (yb>n) yb=n;
    za=z0-r; if (za<0) za=0; zb=z0+r; if (zb>n) zb=n;
    // project xy plane
    for (x=xa,xr=x-x0,xx=xr*xr;x<xb;x++,xr++,xx=xr*xr)
     for (y=ya,yr=y-y0,yy=yr*yr;y<yb;y++,yr++,yy=yr*yr)
        {
        zz=rr-xx-yy; if (zz<0) continue; zr=sqrt(zz);
        z=z0-zr; if ((z>0)&&(z<n)) map[x][y][z]=col;
        z=z0+zr; if ((z>0)&&(z<n)) map[x][y][z]=col;
        }
    // project xz plane
    for (x=xa,xr=x-x0,xx=xr*xr;x<xb;x++,xr++,xx=xr*xr)
     for (z=za,zr=z-z0,zz=zr*zr;z<zb;z++,zr++,zz=zr*zr)
        {
        yy=rr-xx-zz; if (yy<0) continue; yr=sqrt(yy);
        y=y0-yr; if ((y>0)&&(y<n)) map[x][y][z]=col;
        y=y0+yr; if ((y>0)&&(y<n)) map[x][y][z]=col;
        }
    // project yz plane
    for (y=ya,yr=y-y0,yy=yr*yr;y<yb;y++,yr++,yy=yr*yr)
     for (z=za,zr=z-z0,zz=zr*zr;z<zb;z++,zr++,zz=zr*zr)
        {
        xx=rr-zz-yy; if (xx<0) continue; xr=sqrt(xx);
        x=x0-xr; if ((x>0)&&(x<n)) map[x][y][z]=col;
        x=x0+xr; if ((x>0)&&(x<n)) map[x][y][z]=col;
        }
    }
//---------------------------------------------------------------------------
void LED_cube::glDraw()
    {
    #ifdef __gl_h_
    int x,y,z;
    float p[3],dp=1.0/float(n-1);
    glEnable(GL_BLEND);
    glBlendFunc(GL_ONE,GL_ONE);

    glPointSize(2.0);

    glBegin(GL_POINTS);

    for (p[0]=-0.5,x=0;x<n;x++,p[0]+=dp)
     for (p[1]=-0.5,y=0;y<n;y++,p[1]+=dp)
      for (p[2]=-0.5,z=0;z<n;z++,p[2]+=dp)
        {
        glColor4ubv((BYTE*)(&map[x][y][z]));
        glVertex3fv(p);
        }
    glEnd();
    glDisable(GL_BLEND);
    glPointSize(1.0);
    #endif
    }
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------

class usage:

LED_cube cube;
cube.cls(0x00202020); // clear space to dark gray color
int a=cube.n>>1;      // just place sphere to middle and size almost the whole space
int r=a-3;
cube.sphere(a,a,a,r,0x00FFFFFF);
cube.glDraw();        // just for mine visualization you have to rewrite it to your rendering system

If you want to use C only then decompose class to just global functions and variables and translate C++ operators x++,--,+=,-=,*=,... to C style x=x+1,...

Answer 2

Based on the link, it looks like you're more interested in voxel algorithms for spheres, rather than graphics per se; say something like this page helps with. You don't want a ball, but the surface only.

Midpoint circle algorithm can be used to draw 3D voxel spheres: consider the sphere as a stack of slices, and each slice contains a circle.

In practice, you use two nested midpoint circles, the outer defining the radius for the inner one. (Although a naive algorithm drawing circles on top of each other will likely leave holes in the voxels, the midpoint circle algorithm exploits symmetries, and if correctly implemented, no such holes should occur.)

You build six caps in tandem, like carving a cube into a sphere. Since the surface slopes on each cap are always less than 1, going outwards on a cap will at most change each coordinate by 1, so holes cannot occur.

The problem with this approach is the complexity: each point you calculate may affect up to 48 voxel cells. (At each cap, each point is calculated within an octant, and therefore affects eight cells. There are six caps, and 6*8=48.)

I suggest a different approach.

The equation for the surface of a r-radius sphere centered at x₀, y₀, z₀, is

(x - x₀)² + (y - y₀)² + (z - z₀)² = r²

With integer coordinater, the grid points are rarely exactly on the sphere surface, allow a range of values:

RR_MIN ≤ (x - x₀)² + (y - y₀)² + (z - z₀)² ≤ RR_MAX

where RR_MIN and RR_MAX are constants; specifically, minimum and maximum distance squared to the sphere center.

I recommend using doubled coordinates for general cases. This allows you to select whether the center of the sphere is centered at a coordinate (implying odd diameter), or centered between two adjacent cordinates (implying even diameter).

If you have a SIZE×SIZE×SIZE grid of voxels (lights, building blocks, whatever), then

int sphere_surface(const char x, const char y, const char z,
                   const char size, const int rrmin, const int rrmax)
{
    const int center = size - (size & 1); /* Size rounded down to even integer */
    const int dx = center - x - x,
              dy = center - y - y,
              dz = center - z - z;        /* Doubled coordinates */
    const int rr = dx*dx + dy*dy + dz*dz; /* Distance squared */
    return (rrmin <= rr) && (rr <= rrmax);
}

returns 1 if point (x,y,z) is within the surface region of the sphere centered in the cube. (Technically, it returns if the distance from that point to the center of the size-sized cube is within sqrt(rrmin)/2 and sqrt(rrmax)/2, inclusive.)

The "correct" values of rrmin and rrmax are highly context dependent. rrmax is typically somewhere near size*size (remember, the function uses doubled coordinates), with rrmin somewhat less.

For example, if you have a 3×3×3 grid, you only want the six center cells on each face to fulfill the condition; you can do that with size=3, rrmin=1, rrmax=4:

If you have a 4×4×4 grid, you want the four center cells on each face to fulfill the condition (so a total of 24 of the 64 cells are considered to be on the sphere surface); you can do that with size=4, rrmin=11, rrmax=11:

With a 5×5×5 grid, we get to the interesting side effects of the above algorithm.

size=5, rrmin=8, rrmax=16 yields a very "angular" sphere, almost a cube standing on a corner:

size=5, rrmin=12, rrmax=20 yields my favourite approximate sphere:

size=5, rrmin=16, rrmax=24 yields a rounded cube (each face a 3×3 flat):

Obviously, using rrmin=0 includes all inner cells, too, yielding a ball instead of just the surface of a sphere.

As the grid size increases, the more variants of each size sphere you can represent.

The above function is especially useful on microcontrollers, because you can simply loop through your lattice, updating the state at each point, as you wish. Furthermore, most microcontrollers are tight on memory, but have very fast (single-clock) addition, subtraction, and multiplication instructions. (Although a 16×16-bit multiplication with 32-bit result typically takes two or more instructions.)

A typical microcontroller does not have the ROM/flash capability to store enough interesting voxel patterns, and only a few have DMA capability from an SD card via SPI (so you could just load "frames" of voxel patterns from a microSD card), but functions like the above can produce interesting shapes with few inputs -- and particularly inputs you can interpolate.

The above function can also be adapted for approximated antialising (by comparing rr to rrmin and rrmax), in case your voxels are not just binary, but e.g. PWM-controlled LEDs.

---

Since visualizing the above is typically a bit difficult, here is the quick little hack that I used to generate the above images. It outputs an SVG image to standard output, and ASCII slices to standard error, when given size, rrmin, and rrmax as command-line parameters.

#include <stdlib.h>
#include <string.h>
#include <stdio.h>

#define   BORDER  2
#define   XC(x,y,z) ((x)*16 + (y)*12)
#define   YC(x,y,z) ((x)*6 - (y)*8 - (z)*17)

static int xt = 0;
static int yt = 0;

static void fcube(FILE *out, const int x, const int y, const int z, const int fill)
{
    const int xc = xt + XC(x,y,z);
    const int yc = yt + YC(x,y,z);
    fprintf(out, "<path d=\"M%d,%dl16,6,12,-8,0,-17,-16,-6,-12,8z\" fill=\"#%06x\" stroke=\"#000\" />\n", xc, yc, fill & 0xFFFFFF);
    fprintf(out, "<path d=\"M%d,%dl16,6,12,-8m-12,8l0,17\" fill=\"none\" stroke=\"#000\" />\n", xc, yc-17);
}

static unsigned long rrmin = 0UL;
static unsigned long rrmax = 0UL;
static int           center = 0;

static int surface(const int x, const int y, const int z)
{
    /* Doubled coordinates: */
    const long dx = 2*x - center,
               dy = 2*y - center,
               dz = 2*z - center;
    const unsigned long rr = dx*dx + dy*dy + dz*dz;

    return (rrmin <= rr) && (rr <= rrmax);
}

int main(int argc, char *argv[])
{
    int  width, height;
    int  size, x, y, z;
    char dummy;

    if (argc != 4 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) {
        fprintf(stderr, "\n");
        fprintf(stderr, "Usage: %s SIZE RRMIN RRMAX\n", argv[0]);
        fprintf(stderr, "Where\n");
        fprintf(stderr, "       SIZE    is the size of the voxel cube\n");
        fprintf(stderr, "       RRMIN   is the minimum distance squared, and\n");
        fprintf(stderr, "       RRMAX   is the maximum distance squared,\n");
        fprintf(stderr, "               using doubled coordinates.\n");
        fprintf(stderr, "\n");
        fprintf(stderr, "Examples:\n");
        fprintf(stderr, "       %s 3 1 4\n", argv[0]);
        fprintf(stderr, "       %s 4 11 11\n", argv[0]);
        fprintf(stderr, "       %s 5 8 16\n", argv[0]);
        fprintf(stderr, "       %s 5 12 20\n", argv[0]);
        fprintf(stderr, "       %s 5 16 24\n", argv[0]);
        fprintf(stderr, "\n");
        return EXIT_FAILURE;
    }

    if (sscanf(argv[1], " %d %c", &size, &dummy) != 1 || size < 3) {
        fprintf(stderr, "%s: Invalid size.\n", argv[1]);
        return EXIT_FAILURE;
    }

    if (sscanf(argv[2], " %lu %c", &rrmin, &dummy) != 1) {
        fprintf(stderr, "%s: Invalid rrmin.\n", argv[2]);
        return EXIT_FAILURE;
    }

    if (sscanf(argv[3], " %lu %c", &rrmax, &dummy) != 1 || rrmax < rrmin) {
        fprintf(stderr, "%s: Invalid rrmax.\n", argv[3]);
        return EXIT_FAILURE;
    }

    /* Calculate coordinate range. */
    {   int xmin = XC(0,0,0);
        int ymin = YC(0,0,0);
        int xmax = XC(0,0,0);
        int ymax = YC(0,0,0);

        for (z = 0; z <= size; z++)
            for (y = 0; y <= size; y++)
                for (x = 0; x <= size; x++) {
                    const int xc = XC(x,y,z);
                    const int yc = YC(x,y,z);
                    if (xc < xmin) xmin = xc;
                    if (xc > xmax) xmax = xc;
                    if (yc < ymin) ymin = yc;
                    if (yc > ymax) ymax = yc;
                } 

        xt = BORDER - xmin;
        width = xmax - xmin + 2*BORDER;

        yt = BORDER - ymin;
        height = ymax - ymin + 2*BORDER;
    }

    center = size - 1;

    /* SVG preamble. */
    printf("<?xml version=\"1.0\"?>\n");
    printf("<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 %d %d\">\n", width, height);
    printf("<path d=\"M%d,%dL%d,%d,%d,%d,%d,%d,%d,%d,%d,%dz\" fill=\"#f7f7f7\" stroke=\"#666666\"/>\n",
            xt+XC(   0,   0,   0), yt+YC(   0,   0,   0),
            xt+XC(size,   0,   0), yt+YC(size,   0,   0),
            xt+XC(size,size,   0), yt+YC(size,size,   0),
            xt+XC(size,size,size), yt+YC(size,size,size),
            xt+XC(0,   size,size), yt+YC(   0,size,size),
            xt+XC(0,      0,size), yt+YC(   0,   0,size));
    printf("<path d=\"M%d,%dL%d,%d,%d,%dM%d,%dL%d,%d\" fill=\"none\" stroke=\"#666666\"/>\n",
            xt+XC(   0,   0,   0), yt+YC(   0,   0,   0),
            xt+XC(   0,size,   0), yt+YC(   0,size,   0),
            xt+XC(size,size,   0), yt+YC(size,size,   0),
            xt+XC(   0,size,   0), yt+YC(   0,size,   0),
            xt+XC(   0,size,size), yt+YC(   0,size,size));
    for (z = 0; z < size; z++)
        for (y = size - 1; y >= 0; y--)
            for (x = 0; x < size; x++)
                if (surface(x,y,z))
                    fcube(stdout, x, y, z, 0x00CCFF);
    printf("</svg>\n");

    for (z=0; z < size; z++) {
        for (y = 0; y < size; y++) {
            for (x = 0; x < size; x++)
                fputc(surface(x,y,z) ? 'X' : '.', stderr);
            fputs("   ", stderr);
            for (x = 0; x < size; x++)
                fputc(surface(x,z,y) ? 'X' : '.', stderr);
            fputs("   ", stderr);
            for (x = 0; x < size; x++)
                fputc(surface(y,z,x) ? 'X' : '.', stderr);
            fputc('\n', stderr);
        }
        fputc('\n', stderr);
    }

    return EXIT_SUCCESS;
}

I didn't bother to finesse the output; you can quite easily e.g. choose different colors for each face, maybe add shadows to the background planes, et cetera.

The above images were created with this program, then converted to PNG using GIMP, but I recommend using your browser to view the generated files locally.

Questions?