如何在更高维度的超球面上均匀分布点？
Posted
技术标签:
【中文标题】如何在更高维度的超球面上均匀分布点？【英文标题】：How to distribute points evenly on the surface of hyperspheres in higher dimensions? 【发布时间】：2019-11-29 01:58:17 【问题描述】：我有兴趣在 3 维或更高维度的球体表面上均匀分布 N 个点。
更具体一点：
给定点数 N 和维度数 D（其中 D > 1，N > 1） 每个点到原点的距离必须是1 任何两点之间的最小距离应尽可能大 每个点到它最近的邻居的距离不一定对每个点都相同（实际上不可能相同，除非点的数量构成柏拉图立体的顶点，或者如果 N我不感兴趣：
在超球面上创建均匀随机分布，因为我希望任意两点之间的最小距离尽可能大，而不是随机分布。 粒子排斥模拟类型的方法，因为它们难以实现并且需要非常长时间才能运行大 N（理想情况下，该方法应该是确定性的并且在 O(n) 内）。满足这些标准的一种方法称为斐波那契格，但我只能找到 2d 和 3d 中的代码实现。
斐波那契晶格（也称为斐波那契螺旋线）背后的方法是生成一条围绕球体表面盘旋的 1d 线，使得该线所覆盖的表面积在每个转弯处大致相同。然后你可以将 N 个点均匀分布在螺旋上，它们将大致均匀分布在球体表面。
在this answer 中有一个 3 维的 python 实现，它生成以下内容：
我想知道斐波那契螺旋是否可以扩展到高于 3 的维度，并在数学堆栈交换上发布了一个问题。令我惊讶的是，我收到了two amazing answers，据我所知（因为我不完全理解所显示的数学）表明确实可以将此方法扩展到 N 维。
不幸的是，我对所显示的数学了解不够，无法将任一答案转换为（伪）代码。我是一位经验丰富的计算机程序员，但我的数学背景仅此而已。
我将复制以下答案之一中我认为最重要的部分（不幸的是，SO 不支持 mathjax，因此我不得不将其复制为图像）
我遇到的上述困难：
如何求解用于ψn的反函数？ 给出的示例是 d = 3。如何为任意 d 生成公式？了解所涉及数学的任何人都能够在链接斐波那契格子问题的任一答案的伪代码实现方面取得进展吗？我知道完整的实现可能非常困难，所以我很高兴部分实现能够引导我走得足够远，能够自己完成其余部分。
为了方便起见，我已经编写了一个函数，它采用 N 维的球坐标并将它们转换为笛卡尔坐标，因此实现可以输出任何一个，因为我可以轻松转换。
此外，我看到一个答案使用每个附加维度的下一个素数。我可以很容易地编写一个函数来输出每个连续的素数，所以你可以假设它已经实现了。
未能在 N 维中实现斐波那契晶格，我很乐意接受满足上述约束的不同方法。
【问题讨论】：
我知道问题本质上是“从其他答案中提取方程式并将其转换为伪代码”。我希望这是一个合适的问题类型，但如果不是，请告诉我。此外，如果我应该将该答案中的任何信息复制到此问题中，请告诉我，以使其不再是“仅链接”类型的问题。 你能编辑你的问题并在这里简要定义基本概念吗？例如，如果我知道斐波那契格子是什么，我也许可以实现一个 n 维斐波那契格子，但不知道它我很遗憾会跳过这个问题，因为空闲时间很少。 @LajosArpad 我希望我现在添加了一些对您有帮助的细节。 感谢您提供更多信息，但我仍然不知道斐波那契格子是什么。您已经给出了一些关于它的属性，但没有定义这个概念。我会看看我是否有时间研究它，但不幸的是，这不太可能。 感谢您的努力。我知道这是一个相当复杂的概念，除非您有先验知识，否则可能需要在math.stackexchange.com/a/3297830/688579 完整阅读链接的问题以正确理解。我知道这需要相当多的努力，这就是为什么我提供我所有的代表作为赏金，如果我能提供更多，那么我会。不幸的是，Stack Overflow 不支持数学 jax，这限制了我可以从该问题复制到该问题的数量，而不会变得乏味。 【参考方案1】：作为部分答案，您可以使用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 只会返回原始 x。
【讨论】：
很好，这很好地解决了反函数。剩下的唯一问题是如何为任意 d 中的角度生成公式。 不错且简短的实现。【参考方案2】：我们有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
当且仅当 MD 不能更高时，解决方案才是可接受的。
如果 d = 2，那么我们有一个圆并在上面放点。圆是一个多边形，具有以下属性：
它有 n 个角 n -> 无穷大 每一边的长度相似因此，n 个角的多边形，其中 n 是一个有限数且大于 2，而且，每条边的长度相似，每次增加 n 时，它更接近一个圆。请注意，d = 2 中的第一个多边形是三角形。我们只有一个角度，我们的最小角度单位是 360 度/n。
现在，如果我们有一个正方形并将点均匀分布在上面，那么通过base transformation 将我们的正方形转换为圆形应该是精确解，或者非常接近它。如果它是精确解，那么这对于 d = 2 的情况是一个简单的解。如果它仅非常接近，那么通过近似方法，我们可以确定解在一个给定您选择的精度。
我会在 d = 3 的情况下使用这个想法。我会解决一个立方体的问题，这个问题要简单得多，并使用基础变换将我的立方体点转换为我的球体点。我会在 d > 3 上使用这种方法，解决超立方体的问题并将其转换为超球面。当您将点均匀分布在 d 维的超立方体上时，请使用曼哈顿距离。
请注意，我不知道超立方体转换为超球面的解是精确解还是接近它，但如果不是精确解，那么我们可以通过近似来提高精度。
所以，这种方法是解决问题的一种方法，但就时间复杂度而言，这不一定是最好的方法，所以，如果一个人已经深入研究了斐波那契格区域并且知道如何将其推广到更多维度，那么他的/她的回答可能比我的更容易被接受。
如果您定义了 f(x) 的 Taylor series，则可以确定 f(x) = x - 0.5sin2x 的倒数。你会得到一个多项式系列 x which can be inverted。
【讨论】：
所以我们在超立方体表面上等分点，然后把它变成一个超球体，我们将所有点的向量从原点调整为长度为 1。除非我误解了什么你的意思是基础变换，这将导致点在超立方体的角处更加聚集。 @Karl 我同意该解决方案不太可能按原样接受（由于您所说的原因），但也许它可以用来为您提到的粒子排斥方法设置初始状态在 cmets 中。如果初始状态分布良好，那么收敛到最优状态可能会更快。 @JohnColeman 在过去的 4 年里，我一直在研究解决这个问题的粒子排斥方法。我研究的一个领域是使用此答案中描述的技术播种粒子排斥方法（如果我正确理解基本变换）。结果还不错，但我现在想研究确定性方法，这就是为什么我想专注于斐波那契格。 @Karl 我们不使用欧几里得几何计算点之间的距离，而是根据我的想法使用曼哈顿距离（添加不同维度的距离）。这当然只是一个起点。如果这恰好导致根据规范的确定性平均分布，那么这是一个解决方案。如果不是，那么这是一个很好的起点，但在这种情况下，它是不确定的。最好知道是否有人有时间检查初始结果是否符合标准，如果不符合标准，离标准还有多远。 @LajosArpad 这似乎是一个充满希望的开始【参考方案3】：非常有趣的问题。我想在mine 4D rendering engine 中实现它，因为我很好奇它会是什么样子，但我太懒了，也没有能力从数学方面处理 ND 超越问题。
相反，我想出了不同的解决方案来解决这个问题。它不是 Fibonaci Latice !!! 相反，我将 hypersphere or n-sphere 的参数方程扩展为 hyperspiral，然后拟合螺旋参数，因此点或多或少是等距的.
我知道这听起来很可怕，但并不难，而且在解决了一些愚蠢的错别字复制/粘贴错误后，结果对我来说看起来是正确的（最后:)）
主要思想是使用超球面的n维参数方程从角度和半径计算其表面点。这里实现：
algorithm to rasterize and fill a hypersphere?参见[edit2]。现在问题归结为两个主要问题：
计算螺丝数量
因此，如果我们希望我们的点是等距的，那么它们必须以等距离放置在螺旋路径上（参见项目符号 #2），而且螺钉本身也应该彼此之间具有相同的距离。为此，我们可以利用超球面的几何特性。让我们从 2D 开始：
很简单screws = r/d。点数也可以推断为points = area/d^2 = PI*r^2/d^2。
所以我们可以简单地将二维螺旋写成：
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 之外没有其他比例，但它可能需要一些不断的调整（没有对生成的点云进行分析，因为结果对我来说看起来不错）
现在我们可以从单个参数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[] 列出了对象的点列表，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;
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 项目预分配空间
因此您需要将其移植到您可以使用的任何列表中（例如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
将点转换为全局或局部坐标（通过点乘直接或逆矩阵）。您可以忽略它，因为它在渲染中只使用过一次，您可以复制这些点（不旋转）...在同一个标题中还有一些常量：
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;
我还在变换矩阵头中集成了向量和矩阵数学模板，所以vector<n> 是 n 维向量，matrix<n> 是 n*n 方阵，但它仅用于渲染，所以你可以再次忽略它。如果您对这里感兴趣的几个链接，所有这一切都源自：
枚举和降维仅用于渲染。 cfg 包含如何将每个维度缩减为 2D。
AnsiString 是来自 VCL 的自重定位字符串，因此请使用 char* 或您在环境中获得的字符串类。 DWORD 只是无符号的 32 位整数。希望我没有忘记什么...
【讨论】：
我已授予赏金，因为这看起来是最有希望的解决方案，而赏金即将到期。但是，由于未声明的变量，我无法运行此代码，您能解决这个问题吗？ @Karl 我在 [edit1] 中添加了更多（缺失）代码和描述。我没有发布其他网格代码，因为恐怕我接近 30K 的答案限制（网格的完整代码有 23KB 不包括帮助文件）如果您遗漏了其他内容，请评论我。 @Karl 我刚刚更新了一点代码（3D 和 4D 的螺丝和点比更好）【参考方案4】：之前的所有答案都有效，但仍然缺少实际代码。少了两件真实的东西，这是一般实现的。
-
我们需要计算
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 个点的以下图像：
【讨论】：
【参考方案5】：我对如何做到这一点有了另一个疯狂的想法。它与我以前的方法完全不同，因此有了新的答案......
其他答案之一建议在超立方体表面上创建点的均匀分布，然后将点到超立方体中心的距离标准化为超空间半径，并将其用于排斥粒子模拟。我过去在 3D 中这样做并取得了不错的效果，但在更高维度上，这将非常缓慢或被 BVH 之类的结构复杂化。
但这让我开始思考如何倒退。所以在超立方体上非线性分布点，所以在归一化后，点在超球面上线性分布......
让我们从简单的 2D 开始
所以我们只需在+/-45 deg 之间步进角度并计算绿点。角度步长da 必须准确划分90 deg 并给出点密度。所以所有 2D 点将是所有面的 +/-1.0 和 tan(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 立方体：
正如您在边缘看到的那样，点密度要小得多...
预览仅使用透视投影，因为这不会生成网格拓扑，仅生成点，因此无法生成横截面...
还要注意，这会在边缘产生一些重复的点（我认为对于某些镜子来说，将角度循环少一次迭代应该可以解决这个问题，但懒得实现）
我强烈建议阅读以下内容：
uniform points on sphere using Gaussian distribution PRNG【讨论】：
以上是关于如何在更高维度的超球面上均匀分布点？的主要内容，如果未能解决你的问题，请参考以下文章
在Java中如何在更高级别调用函数时显示函数参数？ [复制]