如何在更高维度的超球面上均匀分布点？
How to distribute points evenly on the surface of hyperspheres in higher dimensions?
我感兴趣的是在 3 维及更高维的球体表面均匀分布 N 个点。
更具体地说：
- 给定点数 N 和维数 D（其中 D > 1，N > 1）
- 每个点到原点的距离必须为1
- 任意两点之间的最小距离应尽可能大
- 每个点与其最近邻点的距离不一定对每个点都相同（事实上，它不可能相同，除非点数形成柏拉图立体的顶点或如果 N <= D）。
我不感兴趣：
- 在超球体上创建均匀随机分布，因为我希望任意两点之间的最小距离尽可能大，而不是随机分布。
- 粒子排斥模拟类型的方法，因为它们很难实现并且需要很长时间才能 运行 对于大 N（理想情况下，该方法应该是确定性的并且在 O(n) 内）。
满足这些标准的一种方法称为斐波那契点阵，但我只能在 2d 和 3d 中找到它的代码实现。
斐波那契格子（也称为斐波那契螺线）背后的方法是生成一条围绕球体表面盘旋的一维线，使得该线所覆盖的表面积在每一圈都大致相同。然后你可以在螺旋上平均分布N个点，它们将大致均匀分布在球体表面。
在 this answer 中有一个针对 3 个维度的 python 实现，它生成以下内容：
我想知道斐波那契螺旋是否可以扩展到大于 3 的维度，并在数学堆栈交换上发布了一个问题。令我惊讶的是，我收到了 two amazing answers 据我所知（因为我不完全理解所显示的数学）表明确实可以将此方法扩展到 N 维。
不幸的是，我对所显示的数学理解不够，无法将任一答案转换为（伪）代码。我是一名经验丰富的计算机程序员，但我的数学背景仅此而已。
我将复制我认为是以下答案之一的最重要部分（不幸的是，SO 不支持 mathjax，所以我不得不复制为图像）
我遇到的上述困难：
- 如何求解用于ψn的反函数？
- 给出的示例是 d = 3。如何为任意 d 生成公式？
这里是否有人了解所涉及的数学能够在链接斐波那契格问题的任一答案的伪代码实现方面取得进展？我知道完整的实施可能非常困难，所以我很高兴有部分实施可以让我走得足够远，能够自己完成其余部分。
为了方便起见，我已经编写了一个函数，该函数采用 N 维球坐标并将其转换为笛卡尔坐标，因此实现可以输出任何一个，因为我可以轻松转换。
此外，我看到一个答案对每个附加维度使用下一个质数。我可以轻松编写一个函数来输出每个连续的素数，因此您可以假设它已经实现了。
如果无法在 N 维中实现斐波那契格，我很乐意接受满足上述限制的不同方法。
作为部分答案，您可以使用 Newton's method 来计算 f 的倒数。使用 x 作为牛顿迭代的初始点是一个不错的选择，因为 f(x) 与 x 的距离永远不会超过 1 个单位。这是一个 Python 实现：
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
有关牛顿法应用的一个很好的事实是，每当 cos(2*x) = -1（除以 0）时，您会自动得到 sin(2*x) = 0，因此 f(x) = x。在这种情况下，永远不会进入 while 循环并且 f_inv 只是 returns 原始 x.
我们有n个点，分别是P1，...，Pn。我们有一个维数 d。每个(i = 1,n)点可以表示为：
Pi = (pi(x1), ..., pi(xd))
我们知道
D(Pi, 0) = 1 <=>
sqrt((pi(x1) - pj(x1))^2 + ... + (pi(xd) - pj(xd))^2) = 1
任意点之间的最小距离，MD为
MD <= D(Pi, Pj)
当且仅当 MD 不能更高时，解决方案才可以接受。
如果d = 2，那么我们有一个圆并在上面放点。圆是具有以下属性的多边形：
- 它有n个角
- n -> 无穷大
- 每一边的长度相似
因此，n 个角的多边形，其中 n 是一个大于 2 的有限数，而且，每增加一个 n 时，每条边的长度都相似，它更接近一个圆。请注意，d = 2 中的第一个多边形是三角形。我们有一个角度，我们的最小角度单位是360度/n。
现在，如果我们有一个正方形并在其上均匀分布点，那么通过 base transformation 将我们的正方形转换为圆形应该是精确解，或者非常接近它。如果它是精确解，那么对于 d = 2 的情况，这是一个简单的解。如果它 仅 非常接近，那么通过近似的方法我们可以确定解是什么在您选择的给定精度内。
我会在 d = 3 的情况下使用这个想法。我会解决一个立方体的问题，这个问题要简单得多，并使用基础变换将我的立方体点转换为我的球体点。我会在 d > 3 上使用这种方法，解决超立方体的问题并将其转换为超球体。当您将点均匀分布在 d 维超立方体上时，请使用曼哈顿距离。
请注意，我不知道超立方体转化为超球体的解是精确解还是接近于精确解，但如果不是精确解，那么我们可以通过近似来提高精度。
因此，这种方法是解决问题的方法，但就时间复杂度而言，这不一定是最佳方法，因此，如果有人深入研究了斐波那契格子区域并知道如何将其推广到更多维度，那么 his/her 答案可能比我的更适合接受。
如果您定义了 Taylor series of f(x). You will get a polynomial series of x which can be inverted.
，则可以确定 f(x) = x - 0.5sin2x 的倒数
非常有趣的问题。我想将其实现到 mine 4D rendering engine 中，因为我很好奇它会是什么样子，但我太懒惰且无能，无法从数学方面处理 ND 超越问题。
相反，我想出了不同的解决方案来解决这个问题。它不是 Fibonaci Latice!!! 相反，我将 hypersphere or n-sphere 的参数方程扩展为 hyperspiral 然后只拟合螺旋参数所以这些点或多或少是等距的。
我知道这听起来很可怕，但并不难，结果对我来说是正确的（在解决了一些愚蠢的拼写错误 copy/paste 错误之后，结果对我来说是正确的）
主要思想是使用超球体的n-dimensional参数方程从角度和半径计算其表面点。这里实现：
参见 [edit2]。现在问题归结为两个主要问题：
计算螺丝数
所以如果我们希望我们的点是等距的，那么它们必须等距离地位于螺旋路径上（参见项目符号 #2）而且螺钉本身也应该具有相同的距离彼此之间。为此，我们可以利用超球体的几何特性。让我们从 2D 开始：
如此简单screws = r/d。点数也可以推断为points = area/d^2 = PI*r^2/d^2.
所以我们可以简单地将 2D 螺旋写为：
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（稍后缩放点）。然后扩展（每个维度只添加角度参数）如下所示：
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);
注意 4D 点只是我的假设。我根据经验发现它们与 constant/d^3 有关但不完全相关。每个角度的螺钉都不同。我的假设是除了 screws^i 没有其他比例，但它可能需要一些不断的调整（没有对结果 point-cloud 进行分析，因为结果对我来说还不错）
现在我们可以从单个参数 t=<0.0,1.0>.
生成螺旋上的任意点
请注意，如果您反转等式，那么 d=f(points) 您可以将点数作为输入值，但请注意它只是近似点数，并不准确！！！
在螺旋上生成台阶，使点等距
这是我跳过代数混乱并改用拟合的部分。我只是二分搜索增量 t 所以结果点是 d 远离前一点。因此，只需生成点 t=0，然后在估计位置附近进行二进制搜索 t，直到距起点 d 远。然后重复此操作直到 t<=1.0 ...
您可以使用二进制搜索或其他任何方式。我知道它不如 O(1) 代数方法快，但不需要为每个维度推导内容...看起来 10 次迭代就足以拟合，所以它也没有那么慢。
这里是我的 4D 引擎的实现 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;
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;
}
其中 n=N 是设置的维度，r 是半径，d 是点之间的所需距离。我使用了很多未在此处声明的东西，但重要的是 pnt[] 列出了 object 的点列表，并且 as2(i0,i1) 从索引 i0,i1 处的点添加了行到网格。
这里有几张截图...
3D视角：
4D视角：
4D cross-section 超平面 w=0.0:
同理，点数更多，半径更大：
形状随着动画的旋转而变化...
[Edit1] 更多 code/info
这是我的引擎网格 class 的样子：
//---------------------------------------------------------------------------
//--- 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;
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 (!i0) for (j=0;j<n;j++) t0.add(p0[j]);
if (!i1) for (j=0;j<n;j++) t0.add(p1[j]);
if (!i2) for (j=0;j<n;j++) t0.add(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
if (t0.num==n2) for (j=0;j<t0.num;j++) S2.add(t0.dat[j]);
if (t0.num==n3) for (j=0;j<t0.num;j++) S3.add(t0.dat[j]);
}
}
_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
{
for (j=0;j<n;j++) S4.add(p0[j]);
for (j=0;j<n;j++) S4.add(p1[j]);
for (j=0;j<n;j++) S4.add(p2[j]);
for (j=0;j<n;j++) S4.add(p3[j]);
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 (!i0) for (j=0;j<n;j++) t0.add(p0[j]);
if (!i1) for (j=0;j<n;j++) t0.add(p1[j]);
if (!i2) for (j=0;j<n;j++) t0.add(p2[j]);
if (!i3) for (j=0;j<n;j++) t0.add(p3[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
if (t0.num==n2) for (j=0;j<t0.num;j++) S2.add(t0.dat[j]);
if (t0.num==n3) for (j=0;j<t0.num;j++) S3.add(t0.dat[j]);
if (t0.num==n4) for (j=0;j<t0.num;j++) S4.add(t0.dat[j]);
}
}
}
}
glColor4ubv((BYTE*)(&col));
if (render==_render_Wireframe)
{
// add points from higher primitives
for (i=0;i<S2.num;i++) S1.add(S2.dat[i]);
for (i=0;i<S3.num;i++) S1.add(S3.dat[i]);
for (i=0;i<S4.num;i++) S1.add(S4.dat[i]);
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
//---------------------------------------------------------------------------
我使用我的动态列表模板是这样的：
List<double> xxx; 与 double xxx[]; 相同
xxx.add(5); 将 5 添加到列表末尾
xxx[7] 访问数组元素（安全）
xxx.dat[7]访问数组元素（不安全但直接访问速度快）
xxx.num是数组实际使用的大小
xxx.reset()清空数组并设置xxx.num=0
xxx.allocate(100) 为 100 项预分配 space
因此您需要将其移植到您可以使用的任何列表（例如 std:vector<>）。我还使用 5x5 变换矩阵，其中
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
将点转换为全局或局部坐标（通过将直接或逆矩阵乘以点）。你可以忽略它，因为它在渲染中只使用过一次，你可以复制点而不是（不旋转）......在相同的 header 中也有一些常量：
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;
const double rad=180.0/M_PI;
我还在变换矩阵 header 中集成了向量和矩阵数学模板，所以 vector<n> 是 n 维向量，matrix<n> 是 n*n 方阵，但它使用仅用于渲染，所以您可以再次忽略我.如果您对这里感兴趣，请查看所有这些内容的来源链接：
枚举和降维仅用于渲染。 cfg 保存每个维度应该如何减少到 2D。
AnsiString 是来自 VCL 的自重定位字符串，因此请使用您在环境中获得的 char* 或字符串 class。 DWORD 只是无符号的 32 位整数。希望我没有忘记什么...
之前的所有答案都有效，但仍然缺少实际代码。少了两块真块，这个实现一般。
- 我们需要计算
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 个点生成以下图像：
关于如何执行此操作，我有了另一个疯狂的想法。它与我以前的方法完全不同，因此有了新的答案...
好吧，其他答案之一建议在超立方体表面上创建均匀分布的点，然后将点到超立方体中心的距离归一化为超空间半径，并将其用于排斥粒子模拟。我过去在 3D 中这样做并取得了不错的效果，但在更高的维度上，这将非常缓慢或被 BVH 之类的结构复杂化。
但这让我开始思考倒过来做这个怎么样。所以非线性地分布超立方体上的点所以在归一化之后点变成线性分布在超球面上...
让我们从简单的二维开始
所以我们只需在 +/-45 deg 之间步进角度并计算绿点。角度步长 da 必须精确除以 90 deg 并给出点密度。所以所有的 2D 点都将是所有面的 +/-1.0 和 tan(angle) 的组合。
完成所有点后，只需计算每个点到中心的大小并重新缩放，使其等于超球体半径。
这可以很容易地扩展到任何维度
2D以上的每个维度只需加一for cycle ang angle即可迭代
这里是 2D、3D、4D 的 C++ 示例，使用我之前回答的引擎：
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;
pnt.add( x); pnt.add(-y);
pnt.add( x); pnt.add(+y);
pnt.add(-y); pnt.add( x);
pnt.add(+y); pnt.add( x);
}
}
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;
pnt.add( x); pnt.add( y); pnt.add(-z);
pnt.add( x); pnt.add( y); pnt.add(+z);
pnt.add( x); pnt.add(-z); pnt.add( y);
pnt.add( x); pnt.add(+z); pnt.add( y);
pnt.add(-z); pnt.add( x); pnt.add( y);
pnt.add(+z); pnt.add( x); pnt.add( y);
}
}
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;
pnt.add( x); pnt.add( y); pnt.add( z); pnt.add(-w);
pnt.add( x); pnt.add( y); pnt.add( z); pnt.add(+w);
pnt.add( x); pnt.add( y); pnt.add(-w); pnt.add( z);
pnt.add( x); pnt.add( y); pnt.add(+w); pnt.add( z);
pnt.add( x); pnt.add(-w); pnt.add( y); pnt.add( z);
pnt.add( x); pnt.add(+w); pnt.add( y); pnt.add( z);
pnt.add(-w); pnt.add( x); pnt.add( y); pnt.add( z);
pnt.add(+w); pnt.add( x); pnt.add( y); pnt.add( z);
}
}
for (int i=0;i<pnt.num/n;i++) as1(i);
rescale(r,n);
}
//---------------------------------------------------------------------------
n=N 是维度 r 是半径，da 是 [rad] 中的角度步长。
以及透视 2D/3D/4D 预览：
这里有更多点和更好的 3D 尺寸：
立方体图案略微可见，但点距对我来说还可以。很难在GIF上看到它，因为背面点与正面点合并...
这是没有归一化为球体的 2D 正方形和 3D 立方体：
正如您在边缘上看到的那样，点密度要小得多...
预览仅使用透视投影，因为这不会生成网格拓扑，仅生成点，因此无法生成横截面...
另请注意，这会在边缘产生一些重复点（我认为对于某些镜子，将角度循环少一次迭代应该可以解决这个问题，但懒得实施）
强烈建议阅读此文：
我感兴趣的是在 3 维及更高维的球体表面均匀分布 N 个点。
更具体地说：
- 给定点数 N 和维数 D（其中 D > 1，N > 1）
- 每个点到原点的距离必须为1
- 任意两点之间的最小距离应尽可能大
- 每个点与其最近邻点的距离不一定对每个点都相同（事实上，它不可能相同，除非点数形成柏拉图立体的顶点或如果 N <= D）。
我不感兴趣：
- 在超球体上创建均匀随机分布，因为我希望任意两点之间的最小距离尽可能大，而不是随机分布。
- 粒子排斥模拟类型的方法，因为它们很难实现并且需要很长时间才能 运行 对于大 N（理想情况下，该方法应该是确定性的并且在 O(n) 内）。
满足这些标准的一种方法称为斐波那契点阵，但我只能在 2d 和 3d 中找到它的代码实现。
斐波那契格子（也称为斐波那契螺线）背后的方法是生成一条围绕球体表面盘旋的一维线，使得该线所覆盖的表面积在每一圈都大致相同。然后你可以在螺旋上平均分布N个点，它们将大致均匀分布在球体表面。
在 this answer 中有一个针对 3 个维度的 python 实现，它生成以下内容：
我想知道斐波那契螺旋是否可以扩展到大于 3 的维度，并在数学堆栈交换上发布了一个问题。令我惊讶的是，我收到了 two amazing answers 据我所知（因为我不完全理解所显示的数学）表明确实可以将此方法扩展到 N 维。
不幸的是，我对所显示的数学理解不够，无法将任一答案转换为（伪）代码。我是一名经验丰富的计算机程序员，但我的数学背景仅此而已。
我将复制我认为是以下答案之一的最重要部分（不幸的是，SO 不支持 mathjax，所以我不得不复制为图像）
我遇到的上述困难：
- 如何求解用于ψn的反函数？
- 给出的示例是 d = 3。如何为任意 d 生成公式？
这里是否有人了解所涉及的数学能够在链接斐波那契格问题的任一答案的伪代码实现方面取得进展？我知道完整的实施可能非常困难，所以我很高兴有部分实施可以让我走得足够远，能够自己完成其余部分。
为了方便起见，我已经编写了一个函数，该函数采用 N 维球坐标并将其转换为笛卡尔坐标，因此实现可以输出任何一个，因为我可以轻松转换。
此外，我看到一个答案对每个附加维度使用下一个质数。我可以轻松编写一个函数来输出每个连续的素数，因此您可以假设它已经实现了。
如果无法在 N 维中实现斐波那契格，我很乐意接受满足上述限制的不同方法。
作为部分答案，您可以使用 Newton's method 来计算 f 的倒数。使用 x 作为牛顿迭代的初始点是一个不错的选择，因为 f(x) 与 x 的距离永远不会超过 1 个单位。这是一个 Python 实现：
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
有关牛顿法应用的一个很好的事实是，每当 cos(2*x) = -1（除以 0）时，您会自动得到 sin(2*x) = 0，因此 f(x) = x。在这种情况下，永远不会进入 while 循环并且 f_inv 只是 returns 原始 x.
我们有n个点，分别是P1，...，Pn。我们有一个维数 d。每个(i = 1,n)点可以表示为：
Pi = (pi(x1), ..., pi(xd))
我们知道
D(Pi, 0) = 1 <=>
sqrt((pi(x1) - pj(x1))^2 + ... + (pi(xd) - pj(xd))^2) = 1
任意点之间的最小距离，MD为
MD <= D(Pi, Pj)
当且仅当 MD 不能更高时，解决方案才可以接受。
如果d = 2，那么我们有一个圆并在上面放点。圆是具有以下属性的多边形：
- 它有n个角
- n -> 无穷大
- 每一边的长度相似
因此，n 个角的多边形，其中 n 是一个大于 2 的有限数，而且，每增加一个 n 时，每条边的长度都相似，它更接近一个圆。请注意，d = 2 中的第一个多边形是三角形。我们有一个角度，我们的最小角度单位是360度/n。
现在，如果我们有一个正方形并在其上均匀分布点，那么通过 base transformation 将我们的正方形转换为圆形应该是精确解，或者非常接近它。如果它是精确解，那么对于 d = 2 的情况，这是一个简单的解。如果它 仅 非常接近，那么通过近似的方法我们可以确定解是什么在您选择的给定精度内。
我会在 d = 3 的情况下使用这个想法。我会解决一个立方体的问题，这个问题要简单得多，并使用基础变换将我的立方体点转换为我的球体点。我会在 d > 3 上使用这种方法，解决超立方体的问题并将其转换为超球体。当您将点均匀分布在 d 维超立方体上时，请使用曼哈顿距离。
请注意，我不知道超立方体转化为超球体的解是精确解还是接近于精确解，但如果不是精确解，那么我们可以通过近似来提高精度。
因此，这种方法是解决问题的方法，但就时间复杂度而言，这不一定是最佳方法，因此，如果有人深入研究了斐波那契格子区域并知道如何将其推广到更多维度，那么 his/her 答案可能比我的更适合接受。
如果您定义了 Taylor series of f(x). You will get a polynomial series of x which can be inverted.
，则可以确定 f(x) = x - 0.5sin2x 的倒数非常有趣的问题。我想将其实现到 mine 4D rendering engine 中，因为我很好奇它会是什么样子，但我太懒惰且无能，无法从数学方面处理 ND 超越问题。
相反，我想出了不同的解决方案来解决这个问题。它不是 Fibonaci Latice!!! 相反，我将 hypersphere or n-sphere 的参数方程扩展为 hyperspiral 然后只拟合螺旋参数所以这些点或多或少是等距的。
我知道这听起来很可怕，但并不难，结果对我来说是正确的（在解决了一些愚蠢的拼写错误 copy/paste 错误之后，结果对我来说是正确的）
主要思想是使用超球体的n-dimensional参数方程从角度和半径计算其表面点。这里实现：
参见 [edit2]。现在问题归结为两个主要问题：
计算螺丝数
所以如果我们希望我们的点是等距的，那么它们必须等距离地位于螺旋路径上（参见项目符号 #2）而且螺钉本身也应该具有相同的距离彼此之间。为此，我们可以利用超球体的几何特性。让我们从 2D 开始：
如此简单
screws = r/d。点数也可以推断为points = area/d^2 = PI*r^2/d^2.所以我们可以简单地将 2D 螺旋写为：
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（稍后缩放点）。然后扩展（每个维度只添加角度参数）如下所示：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);注意 4D 点只是我的假设。我根据经验发现它们与
constant/d^3有关但不完全相关。每个角度的螺钉都不同。我的假设是除了screws^i没有其他比例，但它可能需要一些不断的调整（没有对结果 point-cloud 进行分析，因为结果对我来说还不错）现在我们可以从单个参数
生成螺旋上的任意点t=<0.0,1.0>.请注意，如果您反转等式，那么
d=f(points)您可以将点数作为输入值，但请注意它只是近似点数，并不准确！！！在螺旋上生成台阶，使点等距
这是我跳过代数混乱并改用拟合的部分。我只是二分搜索增量
t所以结果点是d远离前一点。因此，只需生成点t=0，然后在估计位置附近进行二进制搜索t，直到距起点d远。然后重复此操作直到t<=1.0...您可以使用二进制搜索或其他任何方式。我知道它不如
O(1)代数方法快，但不需要为每个维度推导内容...看起来 10 次迭代就足以拟合，所以它也没有那么慢。
这里是我的 4D 引擎的实现 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;
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;
}
其中 n=N 是设置的维度，r 是半径，d 是点之间的所需距离。我使用了很多未在此处声明的东西，但重要的是 pnt[] 列出了 object 的点列表，并且 as2(i0,i1) 从索引 i0,i1 处的点添加了行到网格。
这里有几张截图...
3D视角：
4D视角：
4D cross-section 超平面 w=0.0:
同理，点数更多，半径更大：
形状随着动画的旋转而变化...
[Edit1] 更多 code/info
这是我的引擎网格 class 的样子：
//---------------------------------------------------------------------------
//--- 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;
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 (!i0) for (j=0;j<n;j++) t0.add(p0[j]);
if (!i1) for (j=0;j<n;j++) t0.add(p1[j]);
if (!i2) for (j=0;j<n;j++) t0.add(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
if (t0.num==n2) for (j=0;j<t0.num;j++) S2.add(t0.dat[j]);
if (t0.num==n3) for (j=0;j<t0.num;j++) S3.add(t0.dat[j]);
}
}
_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
{
for (j=0;j<n;j++) S4.add(p0[j]);
for (j=0;j<n;j++) S4.add(p1[j]);
for (j=0;j<n;j++) S4.add(p2[j]);
for (j=0;j<n;j++) S4.add(p3[j]);
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 (!i0) for (j=0;j<n;j++) t0.add(p0[j]);
if (!i1) for (j=0;j<n;j++) t0.add(p1[j]);
if (!i2) for (j=0;j<n;j++) t0.add(p2[j]);
if (!i3) for (j=0;j<n;j++) t0.add(p3[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
if (t0.num==n2) for (j=0;j<t0.num;j++) S2.add(t0.dat[j]);
if (t0.num==n3) for (j=0;j<t0.num;j++) S3.add(t0.dat[j]);
if (t0.num==n4) for (j=0;j<t0.num;j++) S4.add(t0.dat[j]);
}
}
}
}
glColor4ubv((BYTE*)(&col));
if (render==_render_Wireframe)
{
// add points from higher primitives
for (i=0;i<S2.num;i++) S1.add(S2.dat[i]);
for (i=0;i<S3.num;i++) S1.add(S3.dat[i]);
for (i=0;i<S4.num;i++) S1.add(S4.dat[i]);
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
//---------------------------------------------------------------------------
我使用我的动态列表模板是这样的：
List<double> xxx; 与 double xxx[]; 相同
xxx.add(5); 将 5 添加到列表末尾
xxx[7] 访问数组元素（安全）
xxx.dat[7]访问数组元素（不安全但直接访问速度快）
xxx.num是数组实际使用的大小
xxx.reset()清空数组并设置xxx.num=0
xxx.allocate(100) 为 100 项预分配 space
因此您需要将其移植到您可以使用的任何列表（例如 std:vector<>）。我还使用 5x5 变换矩阵，其中
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
将点转换为全局或局部坐标（通过将直接或逆矩阵乘以点）。你可以忽略它，因为它在渲染中只使用过一次，你可以复制点而不是（不旋转）......在相同的 header 中也有一些常量：
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;
const double rad=180.0/M_PI;
我还在变换矩阵 header 中集成了向量和矩阵数学模板，所以 vector<n> 是 n 维向量，matrix<n> 是 n*n 方阵，但它使用仅用于渲染，所以您可以再次忽略我.如果您对这里感兴趣，请查看所有这些内容的来源链接：
枚举和降维仅用于渲染。 cfg 保存每个维度应该如何减少到 2D。
AnsiString 是来自 VCL 的自重定位字符串，因此请使用您在环境中获得的 char* 或字符串 class。 DWORD 只是无符号的 32 位整数。希望我没有忘记什么...
之前的所有答案都有效，但仍然缺少实际代码。少了两块真块，这个实现一般。
- 我们需要计算
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 个点生成以下图像：
关于如何执行此操作，我有了另一个疯狂的想法。它与我以前的方法完全不同，因此有了新的答案...
好吧，其他答案之一建议在超立方体表面上创建均匀分布的点，然后将点到超立方体中心的距离归一化为超空间半径，并将其用于排斥粒子模拟。我过去在 3D 中这样做并取得了不错的效果，但在更高的维度上，这将非常缓慢或被 BVH 之类的结构复杂化。
但这让我开始思考倒过来做这个怎么样。所以非线性地分布超立方体上的点所以在归一化之后点变成线性分布在超球面上...
让我们从简单的二维开始
所以我们只需在 +/-45 deg 之间步进角度并计算绿点。角度步长 da 必须精确除以 90 deg 并给出点密度。所以所有的 2D 点都将是所有面的 +/-1.0 和 tan(angle) 的组合。
完成所有点后，只需计算每个点到中心的大小并重新缩放，使其等于超球体半径。
这可以很容易地扩展到任何维度
2D以上的每个维度只需加一for cycle ang angle即可迭代
这里是 2D、3D、4D 的 C++ 示例，使用我之前回答的引擎：
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;
pnt.add( x); pnt.add(-y);
pnt.add( x); pnt.add(+y);
pnt.add(-y); pnt.add( x);
pnt.add(+y); pnt.add( x);
}
}
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;
pnt.add( x); pnt.add( y); pnt.add(-z);
pnt.add( x); pnt.add( y); pnt.add(+z);
pnt.add( x); pnt.add(-z); pnt.add( y);
pnt.add( x); pnt.add(+z); pnt.add( y);
pnt.add(-z); pnt.add( x); pnt.add( y);
pnt.add(+z); pnt.add( x); pnt.add( y);
}
}
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;
pnt.add( x); pnt.add( y); pnt.add( z); pnt.add(-w);
pnt.add( x); pnt.add( y); pnt.add( z); pnt.add(+w);
pnt.add( x); pnt.add( y); pnt.add(-w); pnt.add( z);
pnt.add( x); pnt.add( y); pnt.add(+w); pnt.add( z);
pnt.add( x); pnt.add(-w); pnt.add( y); pnt.add( z);
pnt.add( x); pnt.add(+w); pnt.add( y); pnt.add( z);
pnt.add(-w); pnt.add( x); pnt.add( y); pnt.add( z);
pnt.add(+w); pnt.add( x); pnt.add( y); pnt.add( z);
}
}
for (int i=0;i<pnt.num/n;i++) as1(i);
rescale(r,n);
}
//---------------------------------------------------------------------------
n=N 是维度 r 是半径，da 是 [rad] 中的角度步长。
以及透视 2D/3D/4D 预览：
这里有更多点和更好的 3D 尺寸：
立方体图案略微可见，但点距对我来说还可以。很难在GIF上看到它，因为背面点与正面点合并...
这是没有归一化为球体的 2D 正方形和 3D 立方体：
正如您在边缘上看到的那样，点密度要小得多...
预览仅使用透视投影，因为这不会生成网格拓扑，仅生成点，因此无法生成横截面...
另请注意，这会在边缘产生一些重复点（我认为对于某些镜子，将角度循环少一次迭代应该可以解决这个问题，但懒得实施）
强烈建议阅读此文：