geometry - algorithm to rasterize and fill a hypersphere?

Question

I'm trying to rasterize and fill a hypersphere. In essence, I have a d-dimensional grid of fixed size and a sphere (center, radius) and want to find out which cells of the grid overlap with the sphere and store their coordinates.

I am aware of the Midpoint circle algorithm which takes advantage of 8-way mirroring and produces the outer cells (border) of a circle. I have also altered the linked wikipedia code so as to fill the circle (that is, to produce the coordinates of all cells inside the border).

However I am unaware of any algorithms for higher dimension. For example in 4d, I've been thinking of implementing by producing all possible circles like in the following pseudocode. The basic idea is that since a 4d sphere is (x-x0)2 + (y-y0)**2 + (z-z0)**2 + (k-k0)**2 = r2, this is equal to (x-x0)2 + (y-y0)**2 = r2 - (z-z0)**2 - (k-k0)**2. Since I know how to draw a circle, i just need to produce all circles for all possible values of z and k.

assume center=(x0,y0,z0,k0) and radius r

for all dimensions equal or higher than 2://this is z and k
  //make a list of possible values this dimension can take
  //from z0 to z0+radius with a step of 1
  all_lists.append([dim0,dim0+1,...,dim0+radius])

produce the product of all the lists in all_lists
//now i have a list [[z0,k0],[z0,k0+1],....,[z0+1,k0],[z0+1,k0+1],....,[z0+radius,k0],...[z0+radius,k0+radius]]

for every element l of the list, compute the radius of the circular "cut"
  l.append(r**2 - z**2 - k**2)

Now call the Midpoint Circle Algorithm, but for every (x,y) pair that it produces, we need to export 4 points, namely (x,y,±z,±k)

This question seems relevant, but I don't understand the answer.

Answer 1

well no one answers for some time so here is simple and obvious C++ solution of mine:

//---------------------------------------------------------------------------
const int N=10;                     // number of dimensions
struct point { double a[N]; };      // N-dimensional point
#define color DWORD                 // type of your color property
//---------------------------------------------------------------------------
// N x nested for(a=a0;a<=a1;a+=da) return false if ended
// it could be optimized a lot but for clarity I leave it as is
// ix=-1 at start and N = number of dimensions / nested count
bool nested_for(double *a0,double *a,double *a1,double *da,int &ix,int N)
    {
    if (ix<0)
        {
        int i;
        if (N<1) return false;                          // N too low
        for (i=0;i<N;i++) a[i]=a0[i];
        for (i=0;i<N;i++) if (a[i]>a1[i]) return false; // a0>a1 for some i that is wrong !!!
        ix=N-1;
        return true;
        }
    for (;;)                                            // a+=da
        {
        a[ix]+=da[ix];
        if (a[ix]<=a1[ix]) break;
        if (!ix) return false;                          // end of nested for
        a[ix]=a0[ix];
        ix--;
        }
    ix=N-1;

    return true;
    }
//---------------------------------------------------------------------------
void draw_point(point &p,color c)
    {
    // draw your point ...
    }
//---------------------------------------------------------------------------
void fill_hypersphere(point &p0,double r,color c)
    {
    int i,ix;
    bool init;
    point a,b,a0,a1,da;
    double aa,rr=r*r;
    for (i=0;i<N;i++) a0.a[i]=-r;           // loop boundary for all axises <-r,+r>
    for (i=0;i<N;i++) a1.a[i]=+r;
    for (i=0;i<N;i++) da.a[i]=1.0;          // step between pixels
    for (ix=-1;nested_for(a0.a,a.a,a1.a,da.a,ix,N);) // N x nested for
        {
        aa=0.0;                             // aa = distance from zero ^ 2
        for (i=0;i<N;i++) aa+=a.a[i]*a.a[i];
        if (aa<=rr)                         // if it is inside sphere
            {                               // compute the translated point by p0 to b
            for (i=0;i<N;i++) b.a[i]=p0.a[i]+a.a[i];
            draw_point(b,c);                // and draw/fill it or whatever
            }
        }
    }
//---------------------------------------------------------------------------

This is brute force that can be speeded up using ray castig see:

• Drawing 3D sphere in C/C++

which can boost speed considerably especially in higher dimensions.

[Edit1] just some testing done

screenshot is taken from test app output for code above. Viewing XY plane (z=0) for 1D,2D,3D and 4D hyperspheres. I am not sure about 1D but the rest is OK (not sure if hyperspace is defined for 1D or should be only a point instead).

Even the pixel count ~ Volume looks very similar so the algorithm and code should be OK. Be aware that complexity is O(N!) where N is number of dimensions and runtime is c0*(N!)*r where c0 is constant time, r is radius and N is number of dimensions.

[Edit2] Now how to visualize more than 3D? There are 2 common approaches:

1. projection

   either Orthographic (parallel rays as on images above) or Perspective (the moe distant stuff is smaller). The latter reveals the hidden stuff for example 4D axis aligned Tesseract with orthographic projection to 3D is just a cube, but with Perspective its cube inside bigger cube with all 16 corners interconnected connected...

2. cross section

   You just cut N-dimensional space by a N-dimensional hyperplane and any intersection of your object and hyper plane will give you N-1 dimensional object. This can be applied recursively until hit 3D and render with standard methods.

Both approaches can be combined (some dimensions are reduced by cross sections others with projection ...)

Here some more examples of 4D hypersphere (centered (0,0,0,0) and low poly count so its not a wireframe mess):

Here higher poly count hyper sphere cross-section (W=0.3):

As you can see there are more "grid" like features than the standard parametrical spherical coordinates generated mesh.

However cross section requires that that rendered objects are defined by simplexes covering the objects volume (even surface ones need some kind of volume) otherwise the implementation will get vary nasty filled with edge cases handling.

For more info about 4D rendering see:

• how should i handle (morphing) 4D objects in opengl?

To get back to hypersphere:

according to wiki n-sphere we can describe the surface points of n-sphere by parametric equations:

Where all the angles except last are in interval <0,PI> and the last one is <0,2*PI>. This allows to construct n-sphere manifold/mesh directly. In my engine it looks like this:

//---------------------------------------------------------------------------
void ND_mesh::set_HyperSphere     (int N,double r) // dimensions, radius
    {
    // ToDo
    const int d1=12;
    const int d2=d1+d1;
    const double da1=    M_PI/double(d1-1);
    const double da2=2.0*M_PI/double(d2);
    int i;
    reset(N);
    if (n==2)
        {
        int i0,i1;
        int ia,ja;
        double a;
        pnt.add(0.0);
        pnt.add(0.0);
        for (ja=d2-1,ia=0,a=0.0;ia<d2;ja=ia,ia++,a+=da2)
            {
            pnt.add(cos(a));
            pnt.add(sin(a));
            i0=ja+1;
            i1=ia+1;
            s3(i0,i1,0);
            }
        }
    if (n==3)
        {
        int i00,i01,i10,i11;
        int ia,ib,ja,jb;
        double a,b;
        pnt.add(0.0);
        pnt.add(0.0);
        pnt.add(0.0);
        for (ja=d1-1,ia=0,a=0.0;ia<d1;ja=ia,ia++,a+=da1)
         for (jb=d2-1,ib=0,b=0.0;ib<d2;jb=ib,ib++,b+=da2)
            {
            pnt.add(cos(a));
            pnt.add(sin(a)*cos(b));
            pnt.add(sin(a)*sin(b));
            i00=(ja*d2)+jb+1;
            i01=(ja*d2)+ib+1;
            i10=(ia*d2)+jb+1;
            i11=(ia*d2)+ib+1;
            s4(i00,i01,i11,0);
            s4(i00,i11,i10,0);
            }
        }
    if (n==4)
        {
        int i000,i001,i010,i011,i100,i101,i110,i111;
        int ia,ib,ic,ja,jb,jc;
        double a,b,c;
        pnt.add(0.0);
        pnt.add(0.0);
        pnt.add(0.0);
        pnt.add(0.0);
        for (ja=d1-1,ia=0,a=0.0;ia<d1;ja=ia,ia++,a+=da1)
         for (jb=d1-1,ib=0,b=0.0;ib<d1;jb=ib,ib++,b+=da1)
          for (jc=d2-1,ic=0,c=0.0;ic<d2;jc=ic,ic++,c+=da2)
            {
            pnt.add(cos(a));
            pnt.add(sin(a)*cos(b));
            pnt.add(sin(a)*sin(b)*cos(c));
            pnt.add(sin(a)*sin(b)*sin(c));
            i000=(ja*d1*d2)+(jb*d2)+jc+1;
            i001=(ja*d1*d2)+(jb*d2)+ic+1;
            i010=(ja*d1*d2)+(ib*d2)+jc+1;
            i011=(ja*d1*d2)+(ib*d2)+ic+1;
            i100=(ia*d1*d2)+(jb*d2)+jc+1;
            i101=(ia*d1*d2)+(jb*d2)+ic+1;
            i110=(ia*d1*d2)+(ib*d2)+jc+1;
            i111=(ia*d1*d2)+(ib*d2)+ic+1;
            _cube(i000,i001,i010,i011,i100,i101,i110,i111);
            }
        }

    for (i=0;i<pnt.num;i++) pnt.dat[i]*=r; // radius
    }
//---------------------------------------------------------------------------

Where pnt[] is list of points and _cube adds a cube constructed from 5 tetrahedrons covering its volume. As you can see this creates 2D disc, 3D ball, and 4D sphere (not full volume just layer between manifold neighboring nodes) "surface".

The code just creates n-sphere grid points (with constant angular increase) and then using the 2,4,or 8 neighboring points (and center of sphere for the 2D/3D) to add rendering primitive to object (triangle, tetrahedron, cube).

The renderer then just reduce dimension to 3D and render.

There is also one other approach for this and that is ray tracing. The above stand but with ray tracing we can use algebraic representation of object so we do not need any mesh nor manifolds nor topology. Simply compute closest intersection of a hypersphere and ray (which is just a point) and render accordingly ...