如何在更高维度的超球面上均匀分布点？
Kar*_*arl 18 math geometry pseudocode
我对在尺寸为3或更高的球面上均匀分布N个点感兴趣。
更加具体：
- 给定多个点N和多个维D（其中D> 1，N> 1）
- 每个点到原点的距离必须为1
- 两点之间的最小距离应尽可能大
- 每个点到它最近的邻居的距离不必对于每个点都相同（实际上，除非点的数量形成柏拉图式实体的顶点，或者如果N <= D，则不可能相同。 ）。
我对以下内容不感兴趣：
- 在超球面上创建均匀的随机分布，因为我希望任意两点之间的最小距离尽可能大而不是随机分布。
- 粒子排斥模拟类型的方法，因为它们难以实现并且需要花费很长的时间才能运行较大的N（理想情况下，该方法应该是确定性的，并且以O（n）为单位）。
满足这些条件的一种方法称为斐波那契晶格，但我只能在2d和3d中找到该方法的代码实现。
斐波纳契晶格（也称为斐波纳契螺旋）背后的方法是生成绕球体表面成螺旋形的一维线，以使该线所覆盖的表面积每转大致相同。然后，您可以丢掉均匀分布在螺旋上的N个点，它们将大致均匀地分布在球体的表面上。
在此答案中，有一个针对3个维度的python实现，可生成以下内容：
我想知道斐波那契螺旋是否可以扩展到大于3的尺寸，并在数学堆栈交换中发布了一个问题。令我惊讶的是，我收到了两个令人惊讶的答案，据我所知（因为我不完全理解所显示的数学）表明确实有可能将该方法扩展到N维。
不幸的是，我对所显示的数学知识还不够了解，无法将任何一个答案都转换成（伪）代码。我是一位经验丰富的计算机程序员，但是我的数学背景仅此而已。
我将复制我认为是以下答案之一最重要的部分（不幸的是，SO不支持mathjax，因此我必须复制为图像）
我遇到的上述困难：
- 如何解析用于？n的反函数？
- 给出的示例是d = 3的。如何为任意d生成公式？
在座的任何人都可以理解所涉及的数学知识，从而能够朝着链接斐波那契晶格问题的任一答案的伪代码实现取得进展？我知道完整的实施可能很困难，因此我对部分实施感到满意，该实施可以使我足够自己完成其余的工作。
为简化起见，我已经编写了一个函数，该函数将N个维度的球面坐标转换为笛卡尔坐标，因此该实现可以输出任意一个，因为我可以轻松进行转换。
另外，我看到一个答案为每个附加维使用下一个质数。我可以轻松地编写一个输出每个连续素数的函数，因此可以假定已经实现了。
如果未能在N个维度上实现斐波那契晶格，我很乐意接受满足上述约束的另一种方法。
很有趣的问题。我想将它实现到我的 4D 渲染引擎中，因为我很好奇它会是什么样子，但我太懒惰和无能，无法从数学方面处理 ND 超越问题。
相反，我想出了不同的解决方案来解决这个问题。它不是斐波那契格子！！！相反，我一个的参数方程扩展超球面或正球状成hyperspiral，然后只适合螺旋参数，使这些点都或多或少等距。
我知道这听起来很可怕，但它并不难，在解决了一些愚蠢的错别字复制/粘贴错误后，结果对我来说看起来是正确的（终于:)）
主要思想是使用超球面的 n 维参数方程从角度和半径计算其表面点。这里实现：
见[edit2]。现在问题归结为两个主要问题：
计算螺钉数
因此，如果我们希望我们的点是等距的，那么它们必须以等距分布在螺旋路径上（参见第 2 项），而且螺钉本身之间也应具有相同的距离。为此，我们可以利用超球面的几何特性。让我们从 2D 开始：
就这么简单
screws = r/d。点数也可以推断为points = area/d^2 = PI*r^2/d^2。所以我们可以简单地将二维螺旋写为：
Run Code Online (Sandbox Code Playgroud)t = <0.0,1.0> a = 2.0*M_PI*screws*t; x = r*t*cos(a); y = r*t*sin(a);为了更简单，我们可以
r=1.0这样假设d=d/r（稍后只需缩放点）。然后扩展（每个维度只是添加角度参数）看起来像这样：二维：
Run Code Online (Sandbox Code Playgroud)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：
Run Code Online (Sandbox Code Playgroud)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：
Run Code Online (Sandbox Code Playgroud)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);现在当心 4D 点只是我的假设。我凭经验发现它们与
constant/d^3但不完全相关。每个角度的螺丝都不一样。我的假设是没有其他比例，screws^i但它可能需要一些不断的调整（没有对结果点云进行分析，因为结果对我来说还不错）现在我们可以从单个参数生成螺旋上的任何点
t=<0.0,1.0>。请注意，如果您反转等式以便
d=f(points)您可以将点作为输入值，但要注意它只是近似的点数而不是精确的！！！在螺旋上生成台阶，因此点是等距的
这是我跳过代数混乱并使用拟合的部分。我只是二分搜索增量，
t所以结果点d与前一点相距很远。因此，只需生成点t=0，然后t在估计位置附近进行二分搜索，直到d距离起点较远。然后重复这个直到t<=1.0...您可以使用二分搜索或其他任何方式。我知道它不像
O(1)代数方法那么快，但不需要为每个维度推导出东西......看起来 10 次迭代足以拟合，所以它也不是那么慢。
这是我的 4D 引擎C++/GL/VCL 的实现：
t = <0.0,1.0>
a = 2.0*M_PI*screws*t;
x = r*t*cos(a);
y = r*t*sin(a);
在哪里n=N设置维度，r是半径，d是点之间的所需距离。我使用了很多这里没有声明的东西，但重要的是pnt[]列出对象的点列表，并as2(i0,i1)添加从索引点i0,i1到网格的线。
这里有几个截图...
3D视角：
4D视角：
具有超平面的 4D 横截面w=0.0：
与更多的点和更大的半径相同：
形状随着其动画的旋转而变化......
[Edit1] 更多代码/信息
这是我的引擎网格类的样子：
//---------------------------------------------------------------------------
//--- 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; }
// add simplex
void as1(int a0) { s1.add(a0); }
void as2(int a0,int a1) { s2.add(a0); s2.add(a1); }
void as3(int a0,int a1,int a2) { s3.add(a0); s3.add(a1); s3.add(a2); }
void as4(int a0,int a1,int a2,int a3){ s4.add(a0); s4.add(a1); s4.add(a2); s4.add(a3); }
// 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;
x0=0.0; pnt.add(x0);
y0=0.0; pnt.add(y0);
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;
x0=t*cos(a); pnt.add(x0);
y0=t*sin(a); pnt.add(y0);
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;
x0=1.0; pnt.add(x0);
y0=0.0; pnt.add(y0);
z0=0.0; pnt.add(z0);
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;
x0=cos(a) ; pnt.add(x0);
y0=sin(a)*cos(b); pnt.add(y0);
z0=sin(a)*sin(b); pnt.add(z0);
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;
x0=1.0; pnt.add(x0);
y0=0.0; pnt.add(y0);
z0=0.0; pnt.add(z0);
w0=0.0; pnt.add(w0);
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;
x0=cos(a) ; pnt.add(x0);
y0=sin(a)*cos(b) ; pnt.add(y0);
z0=sin(a)*sin(b)*cos(c); pnt.add(z0);
w0=sin(a)*sin(b)*sin(c); pnt.add(w0);
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)
{
for (j=0;j<n;j++) S1.add(p0[j]);
}
}
_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
{
for (j=0;j<n;j++) S2.add(p0[j]);
for (j=0;j<n;j++) S2.add(p1[j]);
continue;
}
if ((a<=w)&&(b>=w)) // intersection -> points
{
a=(w-p0[N])/(p1[N]-p0[N]);
for (j=0;j<n;j++) S1.add(p0[j]+a*(p1[j]-p0[j]));
}
}
_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
{
for (j=0;j<n;j++) S3.add(p0[j]);
for (j=0;j<n;j++) S3.add(p1[j]);
for (j=0;j<n;j++) S3.add(p2[j]);
continue;
}
if ((a<=w)&&(b>=w)) // cross section -> t0
{
t0.num=0;
之前的所有答案都有效，但仍然缺乏实际的代码。缺少两个真实的部分，这通常是实现的。
- 我们需要计算 的积分
sin^(d-2)(x)。如果您按部分进行递归积分，则它具有封闭形式。在这里，我以递归方式实现它，尽管对于维度 ~> 100 我发现了数值积分sin^d更快 - 我们需要计算该积分的反函数，对于
sin^d，d > 1没有闭合形式。在这里，我使用二分搜索来计算它，尽管其他答案中可能有更好的方法。
这两者与生成素数的方法相结合，得到完整的算法：
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
对于球体上的 200 个点，生成以下图像：
