# How to distribute points evenly on the surface of hyperspheres in higher dimensions?

I am interested in evenly distributing N points on the surface of spheres in dimensions 3 and higher.

To be more specific:

• Given a number of points N and number of dimensions D (where D > 1, N > 1)
• The distance of every point to the origin must be 1
• The minimum distance between any two points should be as large as possible
• The distance of each point to it's closest neighbour doesn't necessarily have to be the same for every point (indeed it's not possible for it to be the same unless the number of points forms the vertices of a platonic solid or if N <= D).

I am not interested in:

• Creating a uniform random distribution on a hypersphere, because I want the minimum distance between any two points to be as large as possible instead of being randomly distributed.
• Particle repulsion simulation type methods, because they are hard to implement and take an extremely long time to run for large N (Ideally the method should be deterministic and in O(n)).

One method that satisfies these criteria is called the fibonacci lattice, but I have only been able to find code implementations for that in 2d and 3d.

The method behind the fibonacci lattice (also known as the fibonacci spiral) is to generate a 1d line that spirals around the surface of the sphere such that the surface area covered by the line is roughly the same at every turn. You can then drop N points equally distributed on the spiral and they will roughly be evenly distributed on the surface of the sphere.

In this answer there is a python implementation for 3 dimensions that generates the following:

I wanted to know whether the fibonacci spiral could be extended to dimensions higher than 3 and posted a question on the maths stack exchange. To my surprise I received two amazing answers which as far as I can tell (because I don't fully understand the maths shown) shows it's indeed possible to extend this method to N dimensions.

Unfortunately I don't understand enough of the maths shown to be able to turn either answer into (pseudo)code. I am an experienced computer programmer, but my maths background only goes so far.

I will copy in what I believe to be the most important part of one of the answers below (unfortunately SO doesn't support mathjax so I had to copy as an image)

Difficulties presented by the above that I struggle with:

• How to resolve the inverse function used for ψn?
• The example given is for d = 3. How do I generate formulae for arbitrary d?

Would anyone here who understands the maths involved be able to make progress towards a pseudo code implementation of either answer to the linked fibonacci lattice question? I understand a full implementation may be quite difficult so I'd be happy with a part implementation that leads me far enough to be able to complete the rest myself.

To make it easier, I've already coded a function that takes spherical coordinates in N dimensions and turns them into cartesian coordinates, so the implementation can output either one as I can easily convert.

Additionally I see that one answer uses the next prime number for each additional dimension. I can easily code a function that outputs each successive prime, so you can assume that's already implemented.

Failing an implementation of the fibonacci lattice in N dimensions, I'd be happy to accept a different method that satisfies the above constraints.

• I understand that the question is essentially "Take the equations from this other answer and turn it into pseudo code". I hope that's an appropriate type of question to ask on here but let me know if it's not. Additionally, let me know if I should copy any information from that answer into this question so that it's less of a "link only" type question. Jul 20, 2019 at 8:59
• Can you edit your question and briefly define the base concepts here? For instance I might be able to implement an n-dimensional Fibonacci lattice if I knew what a Fibonacci lattice is, but not knowing it I unfortunately will skip this question, being low on spare time. Jul 22, 2019 at 10:11
• @LajosArpad I hope I have now added some more detail that will help. Jul 22, 2019 at 10:50
• Thank you for the further information, but I still do not know what a Fibonacci lattice is. You have given some attributes regarding it, but did not define the concept. I will see whether I have time to look into it, but it's inprobable, unfortunately. Jul 22, 2019 at 12:14
• Thanks for the effort. I understand it's quite a complicated concept, and unless you have prior knowledge it probably requires reading the linked question in full at math.stackexchange.com/a/3297830/688579 for a proper understanding. I know that requires quite a bit of effort which is why I've offered all of my rep as a bounty, if I could offer more, then I would. Unfortunately Stack Overflow doesn't support math jax, which limits the amount I can copy from that question into this one without it getting tedious. Jul 22, 2019 at 12:23

Very interesting question. I wanted to implement this into mine 4D rendering engine as I was curious how would it look like but I was too lazy and incompetent to handle ND transcendent problems from the math side.

Instead I come up with different solution to this problem. Its not a Fibonaci Latice !!! Instead I expand the parametrical equation of a hypersphere or n-sphere into hyperspiral and then just fit the spiral parameters so the points are more or less equidistant.

It sounds horrible I know but its not that hard and the results look correct to me (finally :) after solving some silly typos copy/paste bugs)

The main idea is to use the n-dimensional parametrical equations for hypersphere to compute its surface points from angles and radius. Here implementation:

see the [edit2]. Now the problem boils down to 2 main problems:

1. compute number of screws

so if we want that our points are equidistant so they must lye on the spiral path in equidistances (see bullet #2) but also the screws itself should have the same distance between each other. For that we can exploit geometrical properties of the hypersphere. Let start with 2D:

so simply screws = r/d. The number of points can be also inferred as points = area/d^2 = PI*r^2/d^2.

so we can simply write 2D spiral as:

t = <0.0,1.0>
a = 2.0*M_PI*screws*t;
x = r*t*cos(a);
y = r*t*sin(a);


To be more simple we can assume r=1.0 so d=d/r (and just scale the points later). Then the expansions (each dimension just adds angle parameter) look like this:

2D:

screws=1.0/d;          // radius/d
points=M_PI/(d*d);     // surface_area/d^2
a = 2.0*M_PI*t*screws;
x = t*cos(a);
y = t*sin(a);


3D:

screws=M_PI/d;         // half_circumference/d
points=4.0*M_PI/(d*d); // surface_area/d^2
a=    M_PI*t;
b=2.0*M_PI*t*screws;
x=cos(a)       ;
y=sin(a)*cos(b);
z=sin(a)*sin(b);


4D:

screws = M_PI/d;
points = 3.0*M_PI*M_PI*M_PI/(4.0*d*d*d);
a=    M_PI*t;
b=    M_PI*t*screws;
c=2.0*M_PI*t*screws*screws;
x=cos(a)              ;
y=sin(a)*cos(b)       ;
z=sin(a)*sin(b)*cos(c);
w=sin(a)*sin(b)*sin(c);


Now beware points for 4D are just my assumption. I empirically found out that they relate to constant/d^3 but not exactly. The screws are different for each angle. Mine assumption is that there is no other scale than screws^i but it might need some constant tweaking (did not do analysis of the resulting point-cloud as the result look ok to me)

Now we can generate any point on spiral from single parameter t=<0.0,1.0>.

Note if you reverse the equation so d=f(points) you can have points as input value but beware its just approximate number of points not exact !!!

2. generate step on spirals so points are equidistant

This is the part I skip the algebraic mess and use fitting instead. I simply binary search delta t so the resulting point is d distant to previous point. So simply generate point t=0 and then binary search t near estimated position until is d distant to the start point. Then repeat this until t<=1.0 ...

You can use binary search or what ever. I know its not as fast as O(1) algebraic approach but no need to derive the stuff for each dimension... Looks 10 iterations are enough for fitting so its not that slow either.

Here implementation from my 4D engine C++/GL/VCL:

void ND_mesh::set_HyperSpiral(int N,double r,double d)
{
int i,j;
reset(N);
d/=r;   // unit hyper-sphere
double dd=d*d;  // d^2
if (n==2)
{
// r=1,d=!,screws=?
// S = PI*r^2
// screws = r/d
// points = S/d^2
int i0,i;
double a,da,t,dt,dtt;
double x,y,x0,y0;
double screws=1.0/d;
double points=M_PI/(d*d);
dbg=points;
da=2.0*M_PI*screws;
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
x=(t*cos(a))-x0; x*=x;
y=(t*sin(a))-y0; y*=y;
if ((!j)&&(x+y<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y>dd) t-=dtt;
}
if (t>1.0) break;
a=da*t;
as2(i0,i);
}
}
if (n==3)
{
// r=1,d=!,screws=?
// S = 4*PI*r^2
// screws = 2*PI*r/(2*d)
// points = S/d^2
int i0,i;
double a,b,da,db,t,dt,dtt;
double x,y,z,x0,y0,z0;
double screws=M_PI/d;
double points=4.0*M_PI/(d*d);
dbg=points;
da=    M_PI;
db=2.0*M_PI*screws;
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
b=db*t;
x=cos(a)       -x0; x*=x;
y=sin(a)*cos(b)-y0; y*=y;
z=sin(a)*sin(b)-z0; z*=z;
if ((!j)&&(x+y+z<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y+z>dd) t-=dtt;
}
if (t>1.0) break;
a=da*t;
b=db*t;
as2(i0,i);
}
}
if (n==4)
{
// r=1,d=!,screws=?
// S = 2*PI^2*r^3
// screws = 2*PI*r/(2*d)
// points = 3*PI^3/(4*d^3);
int i0,i;
double a,b,c,da,db,dc,t,dt,dtt;
double x,y,z,w,x0,y0,z0,w0;
double screws = M_PI/d;
double points=3.0*M_PI*M_PI*M_PI/(4.0*d*d*d);
dbg=points;
da=    M_PI;
db=    M_PI*screws;
dc=2.0*M_PI*screws*screws;
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
b=db*t;
c=dc*t;
x=cos(a)              -x0; x*=x;
y=sin(a)*cos(b)       -y0; y*=y;
z=sin(a)*sin(b)*cos(c)-z0; z*=z;
w=sin(a)*sin(b)*sin(c)-w0; w*=w;
if ((!j)&&(x+y+z+w<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y+z+w>dd) t-=dtt;
} dt=dtt;
if (t>1.0) break;
a=da*t;
b=db*t;
c=dc*t;
as2(i0,i);
}
}

for (i=0;i<pnt.num;i++) pnt.dat[i]*=r;
for (i=0;i<s1.num;i++) s1.dat[i]*=n;
for (i=0;i<s2.num;i++) s2.dat[i]*=n;
for (i=0;i<s3.num;i++) s3.dat[i]*=n;
for (i=0;i<s4.num;i++) s4.dat[i]*=n;
}


Where n=N are set dimensionality, r is radius and d is desired distance between points. I am using a lot of stuff not declared here but what isimportant is just that pnt[] list the list of points of the object and as2(i0,i1) adds line from points at indexes i0,i1 to the mesh.

Here few screenshots...

3D perspective:

4D perspective:

4D cross-section with hyperplane w=0.0:

and the same with more points and bigger radius:

the shape changes with rotations in which its animated ...

[Edit1] more code/info

This is how my engine mesh class look like:

//---------------------------------------------------------------------------
//--- ND Mesh: ver 1.001 ----------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _ND_mesh_h
#define _ND_mesh_h
//---------------------------------------------------------------------------
#include "list.h"     // my dynamic list you can use std::vector<> instead
#include "nd_reper.h" // this is just 5x5 transform matrix
//---------------------------------------------------------------------------
enum _render_enum
{
_render_Wireframe=0,
_render_Polygon,
_render_enums
};
const AnsiString _render_txt[]=
{
"Wireframe",
"Polygon"
};
enum _view_enum
{
_view_Orthographic=0,
_view_Perspective,
_view_CrossSection,
_view_enums
};
const AnsiString _view_txt[]=
{
"Orthographic",
"Perspective",
"Cross section"
};
struct dim_reduction
{
int view;               // _view_enum
double coordinate;      // cross section hyperplane coordinate or camera focal point looking in W+ direction
double focal_length;
dim_reduction()     { view=_view_Perspective; coordinate=-3.5; focal_length=2.0; }
dim_reduction(dim_reduction& a) { *this=a; }
~dim_reduction()    {}
dim_reduction* operator = (const dim_reduction *a) { *this=*a; return this; }
//dim_reduction* operator = (const dim_reduction &a) { ...copy... return this; }
};
//---------------------------------------------------------------------------
class ND_mesh
{
public:
int n;              // dimensions

List<double> pnt;   // ND points        (x0,x1,x2,x3,...x(n-1))
List<int>    s1;    // ND points        (i0)
List<int>    s2;    // ND wireframe     (i0,i1)
List<int>    s3;    // ND triangles     (i0,i1,i2,)
List<int>    s4;    // ND tetrahedrons  (i0,i1,i2,i3)
DWORD col;          // object color 0x00BBGGRR
int   dbg;          // debug/test variable

ND_mesh()   { reset(0); }
ND_mesh(ND_mesh& a) { *this=a; }
~ND_mesh()  {}
ND_mesh* operator = (const ND_mesh *a) { *this=*a; return this; }
//ND_mesh* operator = (const ND_mesh &a) { ...copy... return this; }

void as1(int a0)                     { s1.add(a0); }
// init ND mesh
void reset(int N);
void set_HyperTetrahedron(int N,double a);              // dimensions, side
void set_HyperCube       (int N,double a);              // dimensions, side
void set_HyperSphere     (int N,double r,int points);   // dimensions, radius, points per axis
void set_HyperSpiral     (int N,double r,double d);     // dimensions, radius, distance between points
// render
void glDraw(ND_reper &rep,dim_reduction *cfg,int render);   // render mesh
};
//---------------------------------------------------------------------------
#define _cube(a0,a1,a2,a3,a4,a5,a6,a7) { as4(a1,a2,a4,a7); as4(a0,a1,a2,a4); as4(a2,a4,a6,a7); as4(a1,a2,a3,a7); as4(a1,a4,a5,a7); }
//---------------------------------------------------------------------------
void ND_mesh::reset(int N)
{
dbg=0;
if (N>=0) n=N;
pnt.num=0;
s1.num=0;
s2.num=0;
s3.num=0;
s4.num=0;
col=0x00AAAAAA;
}
//---------------------------------------------------------------------------
void ND_mesh::set_HyperSpiral(int N,double r,double d)
{
int i,j;
reset(N);
d/=r;   // unit hyper-sphere
double dd=d*d;  // d^2
if (n==2)
{
// r=1,d=!,screws=?
// S = PI*r^2
// screws = r/d
// points = S/d^2
int i0,i;
double a,da,t,dt,dtt;
double x,y,x0,y0;
double screws=1.0/d;
double points=M_PI/(d*d);
dbg=points;
da=2.0*M_PI*screws;
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
x=(t*cos(a))-x0; x*=x;
y=(t*sin(a))-y0; y*=y;
if ((!j)&&(x+y<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y>dd) t-=dtt;
}
if (t>1.0) break;
a=da*t;
as2(i0,i);
}
}
if (n==3)
{
// r=1,d=!,screws=?
// S = 4*PI*r^2
// screws = 2*PI*r/(2*d)
// points = S/d^2
int i0,i;
double a,b,da,db,t,dt,dtt;
double x,y,z,x0,y0,z0;
double screws=M_PI/d;
double points=4.0*M_PI/(d*d);
dbg=points;
da=    M_PI;
db=2.0*M_PI*screws;
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
b=db*t;
x=cos(a)       -x0; x*=x;
y=sin(a)*cos(b)-y0; y*=y;
z=sin(a)*sin(b)-z0; z*=z;
if ((!j)&&(x+y+z<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y+z>dd) t-=dtt;
}
if (t>1.0) break;
a=da*t;
b=db*t;
as2(i0,i);
}
}
if (n==4)
{
// r=1,d=!,screws=?
// S = 2*PI^2*r^3
// screws = 2*PI*r/(2*d)
// points = 3*PI^3/(4*d^3);
int i0,i;
double a,b,c,da,db,dc,t,dt,dtt;
double x,y,z,w,x0,y0,z0,w0;
double screws = M_PI/d;
double points=3.0*M_PI*M_PI*M_PI/(4.0*d*d*d);
dbg=points;
da=    M_PI;
db=    M_PI*screws;
dc=2.0*M_PI*screws*screws;
dt=0.1*(1.0/points);
for (t=0.0,i0=0,i=1;;i0=i,i++)
{
for (dtt=dt,j=0;j<10;j++,dtt*=0.5)
{
t+=dtt;
a=da*t;
b=db*t;
c=dc*t;
x=cos(a)              -x0; x*=x;
y=sin(a)*cos(b)       -y0; y*=y;
z=sin(a)*sin(b)*cos(c)-z0; z*=z;
w=sin(a)*sin(b)*sin(c)-w0; w*=w;
if ((!j)&&(x+y+z+w<dd)){ j--; t-=dtt; dtt*=4.0; continue; }
if (x+y+z+w>dd) t-=dtt;
} dt=dtt;
if (t>1.0) break;
a=da*t;
b=db*t;
c=dc*t;
as2(i0,i);
}
}

for (i=0;i<pnt.num;i++) pnt.dat[i]*=r;
for (i=0;i<s1.num;i++) s1.dat[i]*=n;
for (i=0;i<s2.num;i++) s2.dat[i]*=n;
for (i=0;i<s3.num;i++) s3.dat[i]*=n;
for (i=0;i<s4.num;i++) s4.dat[i]*=n;
}
//---------------------------------------------------------------------------
void ND_mesh::glDraw(ND_reper &rep,dim_reduction *cfg,int render)
{
int N,i,j,i0,i1,i2,i3;
const int n0=0,n1=n,n2=n+n,n3=n2+n,n4=n3+n;
double a,b,w,F,*p0,*p1,*p2,*p3,_zero=1e-6;
vector<4> v;
List<double> tmp,t0;        // temp
List<double> S1,S2,S3,S4;   // reduced simplexes
#define _swap(aa,bb) { double *p=aa.dat; aa.dat=bb.dat; bb.dat=p; int q=aa.siz; aa.siz=bb.siz; bb.siz=q; q=aa.num; aa.num=bb.num; bb.num=q; }

// apply transform matrix pnt -> tmp
tmp.allocate(pnt.num); tmp.num=pnt.num;
for (i=0;i<pnt.num;i+=n)
{
v.ld(0.0,0.0,0.0,0.0);
for (j=0;j<n;j++) v.a[j]=pnt.dat[i+j];
rep.l2g(v,v);
for (j=0;j<n;j++) tmp.dat[i+j]=v.a[j];
}
// copy simplexes and convert point indexes to points (only due to cross section)
S1.allocate(s1.num*n); S1.num=0; for (i=0;i<s1.num;i++) for (j=0;j<n;j++) S1.add(tmp.dat[s1.dat[i]+j]);
S2.allocate(s2.num*n); S2.num=0; for (i=0;i<s2.num;i++) for (j=0;j<n;j++) S2.add(tmp.dat[s2.dat[i]+j]);
S3.allocate(s3.num*n); S3.num=0; for (i=0;i<s3.num;i++) for (j=0;j<n;j++) S3.add(tmp.dat[s3.dat[i]+j]);
S4.allocate(s4.num*n); S4.num=0; for (i=0;i<s4.num;i++) for (j=0;j<n;j++) S4.add(tmp.dat[s4.dat[i]+j]);

// reduce dimensions
for (N=n;N>2;)
{
N--;
if (cfg[N].view==_view_Orthographic){}  // no change
if (cfg[N].view==_view_Perspective)
{
w=cfg[N].coordinate;
F=cfg[N].focal_length;
for (i=0;i<S1.num;i+=n)
{
a=S1.dat[i+N]-w;
if (a>=F) a=F/a; else a=0.0;
for (j=0;j<n;j++) S1.dat[i+j]*=a;
}
for (i=0;i<S2.num;i+=n)
{
a=S2.dat[i+N]-w;
if (a>=F) a=F/a; else a=0.0;
for (j=0;j<n;j++) S2.dat[i+j]*=a;
}
for (i=0;i<S3.num;i+=n)
{
a=S3.dat[i+N]-w;
if (a>=F) a=F/a; else a=0.0;
for (j=0;j<n;j++) S3.dat[i+j]*=a;
}
for (i=0;i<S4.num;i+=n)
{
a=S4.dat[i+N]-w;
if (a>=F) a=F/a; else a=0.0;
for (j=0;j<n;j++) S4.dat[i+j]*=a;
}
}
if (cfg[N].view==_view_CrossSection)
{
w=cfg[N].coordinate;
_swap(S1,tmp); for (S1.num=0,i=0;i<tmp.num;i+=n1)               // points
{
p0=tmp.dat+i+n0;
if (fabs(p0[N]-w)<=_zero)
{
}
}
_swap(S2,tmp); for (S2.num=0,i=0;i<tmp.num;i+=n2)               // lines
{
p0=tmp.dat+i+n0;              a=p0[N];              b=p0[N];// a=min,b=max
p1=tmp.dat+i+n1; if (a>p1[N]) a=p1[N]; if (b<p1[N]) b=p1[N];
if (fabs(a-w)+fabs(b-w)<=_zero)                             // fully inside
{
continue;
}
if ((a<=w)&&(b>=w))                                         // intersection -> points
{
a=(w-p0[N])/(p1[N]-p0[N]);
}
}
_swap(S3,tmp); for (S3.num=0,i=0;i<tmp.num;i+=n3)               // triangles
{
p0=tmp.dat+i+n0;              a=p0[N];              b=p0[N];// a=min,b=max
p1=tmp.dat+i+n1; if (a>p1[N]) a=p1[N]; if (b<p1[N]) b=p1[N];
p2=tmp.dat+i+n2; if (a>p2[N]) a=p2[N]; if (b<p2[N]) b=p2[N];
if (fabs(a-w)+fabs(b-w)<=_zero)                             // fully inside
{
continue;
}
if ((a<=w)&&(b>=w))                                         // cross section -> t0
{
t0.num=0;
i0=0; if (p0[N]<w-_zero) i0=1; if (p0[N]>w+_zero) i0=2;
i1=0; if (p1[N]<w-_zero) i1=1; if (p1[N]>w+_zero) i1=2;
i2=0; if (p2[N]<w-_zero) i2=1; if (p2[N]>w+_zero) i2=2;
if (i0+i1==3){ a=(w-p0[N])/(p1[N]-p0[N]); for (j=0;j<n;j++) t0.add(p0[j]+a*(p1[j]-p0[j])); }
if (i1+i2==3){ a=(w-p1[N])/(p2[N]-p1[N]); for (j=0;j<n;j++) t0.add(p1[j]+a*(p2[j]-p1[j])); }
if (i2+i0==3){ a=(w-p2[N])/(p0[N]-p2[N]); for (j=0;j<n;j++) t0.add(p2[j]+a*(p0[j]-p2[j])); }
if (t0.num==n1) for (j=0;j<t0.num;j++) S1.add(t0.dat[j]);// copy t0 to target simplex based on points count
}
}
_swap(S4,tmp); for (S4.num=0,i=0;i<tmp.num;i+=n4)               // tetrahedrons
{
p0=tmp.dat+i+n0;              a=p0[N];              b=p0[N];// a=min,b=max
p1=tmp.dat+i+n1; if (a>p1[N]) a=p1[N]; if (b<p1[N]) b=p1[N];
p2=tmp.dat+i+n2; if (a>p2[N]) a=p2[N]; if (b<p2[N]) b=p2[N];
p3=tmp.dat+i+n3; if (a>p3[N]) a=p3[N]; if (b<p3[N]) b=p3[N];
if (fabs(a-w)+fabs(b-w)<=_zero)                             // fully inside
{
continue;
}
if ((a<=w)&&(b>=w))                                         // cross section -> t0
{
t0.num=0;
i0=0; if (p0[N]<w-_zero) i0=1; if (p0[N]>w+_zero) i0=2;
i1=0; if (p1[N]<w-_zero) i1=1; if (p1[N]>w+_zero) i1=2;
i2=0; if (p2[N]<w-_zero) i2=1; if (p2[N]>w+_zero) i2=2;
i3=0; if (p3[N]<w-_zero) i3=1; if (p3[N]>w+_zero) i3=2;
if (i0+i1==3){ a=(w-p0[N])/(p1[N]-p0[N]); for (j=0;j<n;j++) t0.add(p0[j]+a*(p1[j]-p0[j])); }
if (i1+i2==3){ a=(w-p1[N])/(p2[N]-p1[N]); for (j=0;j<n;j++) t0.add(p1[j]+a*(p2[j]-p1[j])); }
if (i2+i0==3){ a=(w-p2[N])/(p0[N]-p2[N]); for (j=0;j<n;j++) t0.add(p2[j]+a*(p0[j]-p2[j])); }
if (i0+i3==3){ a=(w-p0[N])/(p3[N]-p0[N]); for (j=0;j<n;j++) t0.add(p0[j]+a*(p3[j]-p0[j])); }
if (i1+i3==3){ a=(w-p1[N])/(p3[N]-p1[N]); for (j=0;j<n;j++) t0.add(p1[j]+a*(p3[j]-p1[j])); }
if (i2+i3==3){ a=(w-p2[N])/(p3[N]-p2[N]); for (j=0;j<n;j++) t0.add(p2[j]+a*(p3[j]-p2[j])); }
if (t0.num==n1) for (j=0;j<t0.num;j++) S1.add(t0.dat[j]);// copy t0 to target simplex based on points count
}
}
}
}
glColor4ubv((BYTE*)(&col));
if (render==_render_Wireframe)
{
// add points from higher primitives
glPointSize(5.0);
glBegin(GL_POINTS);
glNormal3d(0.0,0.0,1.0);
if (n==2) for (i=0;i<S1.num;i+=n1) glVertex2dv(S1.dat+i);
if (n>=3) for (i=0;i<S1.num;i+=n1) glVertex3dv(S1.dat+i);
glEnd();
glPointSize(1.0);
glBegin(GL_LINES);
glNormal3d(0.0,0.0,1.0);
if (n==2)
{
for (i=0;i<S2.num;i+=n1) glVertex2dv(S2.dat+i);
for (i=0;i<S3.num;i+=n3)
{
glVertex2dv(S3.dat+i+n0); glVertex2dv(S3.dat+i+n1);
glVertex2dv(S3.dat+i+n1); glVertex2dv(S3.dat+i+n2);
glVertex2dv(S3.dat+i+n2); glVertex2dv(S3.dat+i+n0);
}
for (i=0;i<S4.num;i+=n4)
{
glVertex2dv(S4.dat+i+n0); glVertex2dv(S4.dat+i+n1);
glVertex2dv(S4.dat+i+n1); glVertex2dv(S4.dat+i+n2);
glVertex2dv(S4.dat+i+n2); glVertex2dv(S4.dat+i+n0);
glVertex2dv(S4.dat+i+n0); glVertex2dv(S4.dat+i+n3);
glVertex2dv(S4.dat+i+n1); glVertex2dv(S4.dat+i+n3);
glVertex2dv(S4.dat+i+n2); glVertex2dv(S4.dat+i+n3);
}
}
if (n>=3)
{
for (i=0;i<S2.num;i+=n1) glVertex3dv(S2.dat+i);
for (i=0;i<S3.num;i+=n3)
{
glVertex3dv(S3.dat+i+n0); glVertex3dv(S3.dat+i+n1);
glVertex3dv(S3.dat+i+n1); glVertex3dv(S3.dat+i+n2);
glVertex3dv(S3.dat+i+n2); glVertex3dv(S3.dat+i+n0);
}
for (i=0;i<S4.num;i+=n4)
{
glVertex3dv(S4.dat+i+n0); glVertex3dv(S4.dat+i+n1);
glVertex3dv(S4.dat+i+n1); glVertex3dv(S4.dat+i+n2);
glVertex3dv(S4.dat+i+n2); glVertex3dv(S4.dat+i+n0);
glVertex3dv(S4.dat+i+n0); glVertex3dv(S4.dat+i+n3);
glVertex3dv(S4.dat+i+n1); glVertex3dv(S4.dat+i+n3);
glVertex3dv(S4.dat+i+n2); glVertex3dv(S4.dat+i+n3);
}
}
glEnd();
}
if (render==_render_Polygon)
{
double nor[3],a[3],b[3],q;
#define _triangle2(ss,p0,p1,p2)                 \
{                                           \
glVertex2dv(ss.dat+i+p0);                   \
glVertex2dv(ss.dat+i+p1);                   \
glVertex2dv(ss.dat+i+p2);                   \
}
#define _triangle3(ss,p0,p1,p2)                 \
{                                           \
for(j=0;(j<3)&&(j<n);j++)                   \
{                                       \
a[j]=ss.dat[i+p1+j]-ss.dat[i+p0+j];     \
b[j]=ss.dat[i+p2+j]-ss.dat[i+p1+j];     \
}                                       \
for(;j<3;j++){ a[j]=0.0; b[j]=0.0; }        \
nor[0]=(a[1]*b[2])-(a[2]*b[1]);             \
nor[1]=(a[2]*b[0])-(a[0]*b[2]);             \
nor[2]=(a[0]*b[1])-(a[1]*b[0]);             \
q=sqrt((nor[0]*nor[0])+(nor[1]*nor[1])+(nor[2]*nor[2]));    \
if (q>1e-10) q=1.0/q; else q-0.0;           \
for (j=0;j<3;j++) nor[j]*=q;                \
glNormal3dv(nor);                           \
glVertex3dv(ss.dat+i+p0);                   \
glVertex3dv(ss.dat+i+p1);                   \
glVertex3dv(ss.dat+i+p2);                   \
}
#define _triangle3b(ss,p0,p1,p2)                \
{                                           \
glNormal3dv(nor3.dat+(i/n));                \
glVertex3dv(ss.dat+i+p0);                   \
glVertex3dv(ss.dat+i+p1);                   \
glVertex3dv(ss.dat+i+p2);                   \
}
glBegin(GL_TRIANGLES);
if (n==2)
{
glNormal3d(0.0,0.0,1.0);
for (i=0;i<S3.num;i+=n3) _triangle2(S3,n0,n1,n2);
for (i=0;i<S4.num;i+=n4)
{
_triangle2(S4,n0,n1,n2);
_triangle2(S4,n3,n0,n1);
_triangle2(S4,n3,n1,n2);
_triangle2(S4,n3,n2,n0);
}
}
if (n>=3)
{
for (i=0;i<S3.num;i+=n3) _triangle3 (S3,n0,n1,n2);
for (i=0;i<S4.num;i+=n4)
{
_triangle3(S4,n0,n1,n2);
_triangle3(S4,n3,n0,n1);
_triangle3(S4,n3,n1,n2);
_triangle3(S4,n3,n2,n0);
}
glNormal3d(0.0,0.0,1.0);
}
glEnd();
#undef _triangle2
#undef _triangle3
}
#undef _swap
}
//---------------------------------------------------------------------------
#undef _cube
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------


I use mine dynamic list template so:

List<double> xxx; is the same as double xxx[];
xxx.add(5); adds 5 to end of the list
xxx[7] access array element (safe)
xxx.dat[7] access array element (unsafe but fast direct access)
xxx.num is the actual used size of the array
xxx.reset() clears the array and set xxx.num=0
xxx.allocate(100) preallocate space for 100 items

so you need to port it to any list you have at disposal (like std:vector<>). I also use 5x5 transform matrix where

void ND_reper::g2l    (vector<4> &l,vector<4> &g);  // global xyzw -> local xyzw
void ND_reper::l2g    (vector<4> &g,vector<4> &l);  // global xyzw <- local xyzw


convert point either to global or local coordinates (by multiplying direct or inverse matrix by point). You can ignore it as its used just once in the rendering and you can copy the points instead (no rotation)... In the same header are also some constants:

const double pi   =    M_PI;
const double pi2  =2.0*M_PI;
const double pipol=0.5*M_PI;
const double deg=M_PI/180.0;


I got also vector and matrix math template integrated in the transform matrix header so vector<n> is n dimensional vector and matrix<n> is n*n square matrix but its used only for rendering so again you can ignore it. If youre interested here few links from whic all this was derived:

The enums and dimension reductions are used only for rendering. The cfg holds how should be each dimension reduced down to 2D.

AnsiString is a self relocating string from VCL so either use char* or string class you got in your environment. DWORD is just unsigned 32 bit int. Hope I did not forget something ...

• I've awarded the bounty because this looks like the most promising solution, and the bounty was about to expire. However I haven't been able to get this code to run due to the undeclared variables, would you be able to resolve that issue? Jul 28, 2019 at 21:21
• @Karl I added [edit1] with more (missing) code and descriptions. I did not post the other meshes code as I am afraid I am near 30K limit of answer (full code of the mesh has 23KB not counting helper files) if you are missing something else comment me. Jul 28, 2019 at 21:48
• @Karl I just updated the code a little (better screw and point ratios for 3D and 4D) Jul 28, 2019 at 22:22

All of the previous answer worked, but still lacked actual code. There were two real pieces missing, which this implements generally.

1. We need to compute the integral of sin^(d-2)(x). This has a closed form if you do recursive integration by parts. Here I implement it in the recursive fashion, although for dimension ~> 100 I found numeric integration of sin^d to be faster
2. We need to compute the inverse function of that integral, which for sin^d, d > 1 doesn't have a close form. Here I compute it using binary search, although there are likely better ways as stated in other answers.

These two combined with a way to generate primes results in the full algorithm:

from itertools import count, islice
from math import cos, gamma, pi, sin, sqrt
from typing import Callable, Iterator, List

def int_sin_m(x: float, m: int) -> float:
"""Computes the integral of sin^m(t) dt from 0 to x recursively"""
if m == 0:
return x
elif m == 1:
return 1 - cos(x)
else:
return (m - 1) / m * int_sin_m(x, m - 2) - cos(x) * sin(x) ** (
m - 1
) / m

def primes() -> Iterator[int]:
"""Returns an infinite generator of prime numbers"""
yield from (2, 3, 5, 7)
composites = {}
ps = primes()
next(ps)
p = next(ps)
assert p == 3
psq = p * p
for i in count(9, 2):
if i in composites:  # composite
step = composites.pop(i)
elif i < psq:  # prime
yield i
continue
else:  # composite, = p*p
assert i == psq
step = 2 * p
p = next(ps)
psq = p * p
i += step
while i in composites:
i += step
composites[i] = step

def inverse_increasing(
func: Callable[[float], float],
target: float,
lower: float,
upper: float,
atol: float = 1e-10,
) -> float:
"""Returns func inverse of target between lower and upper

inverse is accurate to an absolute tolerance of atol, and
must be monotonically increasing over the interval lower
to upper
"""
mid = (lower + upper) / 2
approx = func(mid)
while abs(approx - target) > atol:
if approx > target:
upper = mid
else:
lower = mid
mid = (upper + lower) / 2
approx = func(mid)
return mid

def uniform_hypersphere(d: int, n: int) -> List[List[float]]:
"""Generate n points over the d dimensional hypersphere"""
assert d > 1
assert n > 0
points = [[1 for _ in range(d)] for _ in range(n)]
for i in range(n):
t = 2 * pi * i / n
points[i][0] *= sin(t)
points[i][1] *= cos(t)
for dim, prime in zip(range(2, d), primes()):
offset = sqrt(prime)
mult = gamma(dim / 2 + 0.5) / gamma(dim / 2) / sqrt(pi)

def dim_func(y):
return mult * int_sin_m(y, dim - 1)

for i in range(n):
deg = inverse_increasing(dim_func, i * offset % 1, 0, pi)
for j in range(dim):
points[i][j] *= sin(deg)
points[i][dim] *= cos(deg)
return points


Which produces the following image for 200 points on a sphere:

• Based on your answer, I optimized the code such that it can (approximately) generated a lot of points in high dimensional space (30000 points on 300 dims took around 23s). See pastebin.com/YPpVJMkc Dec 1, 2022 at 21:02

I got another insane idea on how to do this. Its entirely different than my previous approach hence new Answer...

Well one of the other answers suggest creating uniform distribution of points on hypercube surface and then normalize the points distance to center of hypercube to radius of hyperspace and use that for repulsion particle simulation. I did that in past for 3D with good results but in higher dimensions that would be insanely slow or complicated by BVH like structures.

But it got me thinking what about doing this backwards. So distribute the points on hypercube non linearly so after the normalization the points became linearly distributed on hypersphere...

So we simply step angle between +/-45 deg and compute the green points. The angle step da must divide 90 deg exactly and gives the point density. So all the 2D points will be a combination of +/-1.0 and tan(angle) for all faces.

When all points are done simply compute the size of each point to center and rescale it so it will be equal to hypersphere radius.

This can be easily expanded to any dimensionality

Each dimension above 2D just add one for cycle ang angle to iterate.

Here C++ example for 2D,3D,4D using my engine from previous answer of mine:

void ND_mesh::set_HyperSpherePCL(int N,double r,double da)
{
reset(N);

int na=floor(90.0*deg/da);
if (na<1) return;
da=90.0*deg/double(na-1);

if (n==2)
{
int i;
double a,x,y,l;
for (a=-45.0*deg,i=0;i<na;i++,a+=da)
{
x=tan(a); y=1.0;
l=sqrt((x*x)+(y*y));
x/=l; y/=l;
}
}
if (n==3)
{
int i,j;
double a,b,x,y,z,l;
for (a=-45.0*deg,i=0;i<na;i++,a+=da)
for (b=-45.0*deg,j=0;j<na;j++,b+=da)
{
x=tan(a); y=tan(b); z=1.0;
l=sqrt((x*x)+(y*y)+(z*z));
x/=l; y/=l; z/=l;
}
}
if (n==4)
{
int i,j,k;
double a,b,c,x,y,z,w,l;
for (a=-45.0*deg,i=0;i<na;i++,a+=da)
for (b=-45.0*deg,j=0;j<na;j++,b+=da)
for (c=-45.0*deg,k=0;k<na;k++,c+=da)
{
x=tan(a); y=tan(b); z=tan(c); w=1.0;
l=sqrt((x*x)+(y*y)+(z*z)+(w*w));
x/=l; y/=l; z/=l; w/=l;
}
}

for (int i=0;i<pnt.num/n;i++) as1(i);
rescale(r,n);
}
//---------------------------------------------------------------------------


The n=N is dimensionality r is radius and da is angualr step in [rad].

And perspective 2D/3D/4D previews:

And here more points and better size for 3D:

The cube pattern is slightly visible but the point distance looks OK to me. Its hard to see it on GIF as the backside points are merging with the front ones...

And this is the 2D square and 3D cube without normalization to sphere:

as you can see on edges is much smaller point density...

Preview is only using perspective projection as this does not generate mesh topology, just the points so cross section is not possible to do...

Also beware this produces some duplicate points on the edges (I think looping the angles one iteration less for some of the mirrors should remedy that but too lazy to implement that)

I strongly suggest to read this:

As a partial answer, you can use Newton's method to compute the inverse of f. Using x as the initial point in the Newton iteration is a good choice since f(x) is never more than 1 unit away from x. Here is a Python implementation:

import math

def f(x):
return x + 0.5*math.sin(2*x)

def f_inv(x,tol = 1e-8):
xn = x
y = f(xn)
while abs(y-x) > tol:
xn -= (y-x)/(1+math.cos(2*xn))
y = f(xn)
return xn


A nice fact about this application of Newton's method is that whenever cos(2*x) = -1 (where you would have division by 0) you automatically have sin(2*x) = 0 so that f(x) = x. In this case, the while loop is never entered and f_inv simply returns the original x.

• Good, this resolves the inverse function quite nicely. The only problem left is how to generate formulae for the angles in arbitrary d. Jul 22, 2019 at 11:42
• Nice and short implementation. Jul 23, 2019 at 19:35

We have n points, which are P1, ..., Pn. We have a dimension number d. Each (i = 1,n) point can be represented as:

Pi = (pi(x1), ..., pi(xd))

We know that

D(Pi, 0) = 1 <=>

sqrt((pi(x1) - pj(x1))^2 + ... + (pi(xd) - pj(xd))^2) = 1

and the minimal distance between any points, MD is

MD <= D(Pi, Pj)

A solution is acceptable if and only if MD could not be higher.

If d = 2, then we have a circle and put points on it. The circle is a polygon with the following properties:

• it has n angles
• n -> infinity
• each side is of similar length

So, a polygon of n angles, where n is a finite number and higher than 2, also, each side is of similar length is closer to a circle each time we increment n. Note that the firs polygon in d = 2 is the triangle. We have a single angle and our minimal angle unit is 360degrees / n.

Now, if we have a square and distribute the points evenly on it, then converting our square into circle via base transformation should be either the exact solution, or very close to it. If it is the exact solution, then this is a simple solution for the case when d = 2. If it is only very close, then with an approach of approximation we can determine what the solution is within a given precision of your choice.

I would use this idea for the case when d = 3. I would solve the problem for a cube, where the problem is much simpler and use base transformation to convert my cube points to my sphere points. And I would use this approach on d > 3, solving the problem for a hypercube and transform it to a hypersphere. Use the Manhattan distance when you evenly distribute your points on a hypercube of d dimensions.

Note that I do not know whether the solution for a hypercube transformed into a hypersphere is the exact solution or just close to it, but if it's not the exact solution, then we can increase precision with approximation.

So, this approach is a solution for the problem, which is not necessarily the best approach in terms of time complexity, so, if one has delved into the Fibonacci lattice area and knows how to generalize it for more dimensions, then his/her answer might be a better choice for acceptance than mine.

The invert of f(x) = x - 0.5sin2x can be determined if you defined the Taylor series of f(x). You will get a polynomial series of x which can be inverted.

• So we equidistribute points on the surface of a hypercube, and then to turn that into a hyphersphere we resize all of the point's vectors from the origin to have a length of 1. Unless I've misunderstood what you mean by the base transformation, that will have the result of points being more bunched up where the corners of the hypercube were. Jul 22, 2019 at 14:12
• @Karl I agree that the solution is unlikely to be acceptable as is (for the reason that you state), but perhaps it could be used to set up the initial state for the particle repulsion approach you allude to in the comments. If the initial state has a good distribution then convergence to the optimal might be faster. Jul 22, 2019 at 16:39
• @JohnColeman I have been researching particle repulsion methods for this problem for the last 4 years. One of the areas I researched was seeding a particle repulsion method using the technique described in this answer (if I understand the base transformation correctly). The results were not bad, but I wish to study deterministic methods now, which is why I'd like to focus on the fibonacci lattice. Jul 22, 2019 at 17:01
• @Karl we do not calculate the distance between the points using Euclidean geometry, but using Manhattan distance (adding the distances of different dimensions) according to my idea. This is of course only a starting point. If this happens to result in a deterministic equidistribution according to the specs, then this is a solution. If not, then this is a good starting point, but in that case it is undeterministic. It would be good to know if someone has the time to check whether the initial result matches the criteria and if not, how far is it from that. Jul 23, 2019 at 10:35
• @LajosArpad It does seem like a promising start Jul 23, 2019 at 10:41
{"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 1, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0}