How to Rasterize a Sphere

Question

So, I am trying to create a sphere that has "blocks" on its outside, kinda like as if built in Minecraft. (I don't know what the terminology for the outside of a circle is). The problem is, I can't figure out how to get an equation like the Midpoint Circle Algorithm to work for a sphere. Preferably in lua or java so I can read any answers easier. And I don't want to know how to calculate a point on the sphere with trig, I already know how to do that.

Answer 1

I think the easiest way to do this is something like the Midpoint Circle Algorithm, extended to 3D.

First, lets figure out which blocks we want to fill. Assuming an origin in the middle of block (0,0,0) and radius R:

1. We only want to fill boxes inside the sphere. Those are exactly the boxes (x,y,z) such that x²+y²+z² <= R²; and
2. We only want to fill boxes with a face showing. If a box has a face showing, then at least one of its neighbors is not in the sphere, so: (|x|+1)²+y²+z² > R² OR x²+(|y|+1)²+z² > R² OR x²+y²+(|z|+1)² > R²

It's the 2nd part that makes it tricky, but remember that (|a|+1)² = |a|² + 2|a| + 1. If, say, z is the largest coordinate of a box that is inside the sphere, and if that box has a face showing, then the z face in particular will be showing, because x²+y²+(|z|+1)² = x²+y²+z²+2|z|+1, and that will be at least as big as the analogous values for x and y.

So, it's pretty easy to calculate the boxes that are 1) inside the sphere, 2) have z as their largest coordinate, and 3) have the largest possible z value, i.e., adding 1 to z results in a box outside the sphere. Additionally, 4) have positive values for all x,y,z.

The coordinates of these boxes can then be reflected 24 different ways to generate all the boxes on the surface of the sphere. Those are all 8 combinations of signs of the coordinates times all 3 choices for which axis has the largest coordinate.

Here's how to generate the points with positive x,y,z and z largest:

maxR2 = floor(R*R);
zx = floor(R);
for (x=0; ;++x)
{
    //max z for this x value.
    while(x*x+zx*zx > maxR2 && zx>=x)
        --zx;
    if (zx<x)
        break; //with this x, z can't be largest

    z=zx;
    for(y=0; ;++y)
    {
        while(x*x+y*y+z*z > maxR2 && z>=x && z>=y)
            --z;
        if (z<x||z<y)
            break; //with this x and y, z can't be largest
        FILLBOX(x,y,z); //... and up to 23 reflections of it
    }
}

NOTE: If it matters to you, be careful when you calculate the reflections so that you don't draw, for example, (0,y,z) and (-0,y,z), because that's the same box twice. Also don't swap axes with the same value, because again that would draw the same box twice (e.g., if you have (1,5,5), don't swap y and z and draw again.

NOTE ALSO that R doesn't have to be an integer. It'll look a little nicer if you add 0.5 to it.

Here's an example that takes all of the above into account (you need a browser that supports webgl) https://jsfiddle.net/mtimmerm/ay2adpwb/

Answer 2

You can use the Midpoint Circle Algorithm or Bresenham's Circle Algorithm in nested loops. The outer loop determines the integer valued radius of the circle at different z distances from the origin, while the inner loop calculates the x and y elements of the circle comprising the sphere's cross section perpendicular to the Z axis at location z.

Answer 3

This is written in my own version of MiniBasic.

10  REM Shaded textured sphere 
20 INPUT width 
30 INPUT height 
40 INPUT bwidth 
50 INPUT bheight 
60 REM Sphere radius 
70 LET r = 100 
80 REM light direction 
90 LET lx = 0 
100 LET ly = 0.2 
110 LET lz = -1 
120 LET ambient = 0.1 
130 REM camera z position 
140 LET cz = -256 
150 REM Sphere colour 
160 LET red = 255 
170 LET green =0 
180 LET blue = 100 
190 REM Normalize light  
200 LET len = SQRT(lx*lx+ly*ly+lz*lz) 
210 LET lx = -lx/len 
220 LET ly = -ly/len 
230 LET lz = -lz/len 
240 FOR i = 0 TO height -1 
250 FOR ii = 0 TO width -1 
260 LET x = ii - width /2 
270 LET y = i - height/2 
280 LET dpx = x 
290 LET dpy = y 
300 LET dpz = -cz 
310 LET a = dpx*dpx + dpy*dpy + dpz*dpz 
320 LET b = 2 * ( dpz  * cz)  
330 LET c = cz * cz 
340 LET c = c - r * r 
350 LET bb4ac = b * b - 4 * a * c 
360 IF ABS(a) < 0.001 OR bb4ac < 0 THEN 560 
370 LET mu1 = (-b + SQRT(bb4ac)) / (2 * a) 
380 LET mu2 = (-b - SQRT(bb4ac)) / (2 * a) 
390 LET spx = mu1 * x 400 LET spy = mu1 * y 
410 LET spz = mu1 * 256 - 256 
420 LET nx = spx / r 
430 LET ny = spy / r 
440 LET nz = spz / r 
450 LET costh = nx*lx + ny *ly + nz*lz 
460 IF costh > 0 THEN 480 
470 LET costh = 0 
480 LET lum = ambient + (1 - ambient) * costh 
490 LET v = ACOS(nz)/PI 495 LET nx = nx  * 0.999 
510 LET u = ACOS(nx * SIN(PI * v)) / PI   
535 LET v = -ACOS(ny)/PI + 1.0  
536 LET u = 1 - u 
540 GETPIXELB u * bwidth, v * bheight, red, green, blue 
550 SETPIXEL ii, i,  red * lum,  green *lum,  blue *lum 
560 NEXT ii 
570 NEXT i 

http://www.malcolmmclean.site11.com/www/UserPrograms.html#texsphere

Answer 4

Here is a python example using Matt Timmermans example. This is very fast if using numba jits fuctions and parallel prange if you have multiple spheres. svoxel was added for any starting coordinate. The grid simply places 1s where sphere voxels are in a 3d grid. Also added an optional fill inside step. NOTE: This will only work for non-negative coordinates so the sphere needs to be in positive space for all voxels using svoxels coordinate.

#/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
from numba.experimental import jitclass
from numba import types, njit, prange


spec = [
    ('radius', types.float64),
    ('svoxel', types.int64[:]),
    ('grid', types.b1[:,:,:])
]

@jitclass(spec)
class RasterizeSphere(object):
    def __init__(self, svoxel, radius, grid):
        self.grid = grid
        self.svoxel = svoxel
        self.radius = radius
        R2 = np.floor(self.radius**2)
        zx = np.int64(np.floor(self.radius))
        x = 0
        while True:
            while x**2 + zx**2 > R2 and zx >= x:
                zx -= 1
            if zx < x:
                break
            z = zx
            y = 0
            while True:
                while x**2 + y**2 + z**2 > R2 and z >= x and z >= y:
                    z -= 1
                if z < x or z < y:
                    break
                self.fill_all(x, y, z)
                ###### Optional fills the inside of sphere as well. #######
                for nz in range(z):
                    self.fill_all(x, y, nz)
                y += 1
            x += 1

    def fill_signs(self, x, y, z):
        self.grid[x + self.svoxel[0], y + self.svoxel[1], z + self.svoxel[2]] = True
        while True:
            z = -z
            if z >= 0:
                y = -y
                if y >= 0:
                    x = -x
                    if x >= 0:
                        break
            self.grid[x + self.svoxel[0], y + self.svoxel[1], z + self.svoxel[2]] = True

    def fill_all(self, x, y, z):
        self.fill_signs(x, y, z)
        if z > y:
            self.fill_signs(x, z, y)
        if z > x and z > y:
            self.fill_signs(z, y, x)

@njit(parallel=True, cache=True)
def parallel_spheres(grid):
    # prange just to show the insane speedup for large number of spheres disable visualization below if using large amount of prange.
    for i in prange(2):
        radius = 4.5
        svoxel = np.array([5, 5, 5])
        max = np.int64(np.ceil(radius**2))
        rs = RasterizeSphere(svoxel, radius, grid)
        points = np.where(rs.grid)
        return np.array([*points])


def main():
    # Make max large enough to hold the spheres.
    max = 100
    points = parallel_spheres(np.zeros((max, max, max), dtype=bool))
    fig = plt.figure()
    ax = plt.axes(projection='3d')
    ax.scatter3D(*points)
    plt.show()


if __name__ == '__main__':
    main()