25

我有兴趣在 3 维或更高维的球体表面上均匀分布 N 个点。

更具体:

  • 给定点数 N 和维度数 D(其中 D > 1,N > 1)
  • 每个点到原点的距离必须为 1
  • 任意两点之间的最小距离应尽可能大
  • 每个点到它最近的邻居的距离不一定对每个点都相同(实际上它不可能相同,除非点的数量形成柏拉图立体的顶点或如果 N <= D )。

我对以下内容不感兴趣:

  • 在超球面上创建均匀随机分布,因为我希望任意两点之间的最小距离尽可能大,而不是随机分布。
  • 粒子排斥模拟类型的方法,因为它们很难实现并且需要很长时间才能运行大 N(理想情况下,该方法应该是确定性的并且在 O(n) 中)。

满足这些标准的一种方法称为斐波那契格,但我只能找到 2d 和 3d 中的代码实现。

斐波那契格子(也称为斐波那契螺旋线)背后的方法是生成一条围绕球体表面盘旋的 1d 线,使得该线所覆盖的表面积在每一圈都大致相同。然后,您可以将 N 个点均匀分布在螺旋上,它们将大致均匀分布在球体表面。

这个答案中,有一个 3 维的 python 实现,它生成以下内容:

在此处输入图像描述

我想知道斐波那契螺旋是否可以扩展到高于 3 的维度,并在数学堆栈交换中发布了一个问题。令我惊讶的是,我收到了两个惊人的答案,据我所知(因为我不完全理解所显示的数学)表明确实可以将此方法扩展到 N 维。

不幸的是,我对能够将任一答案转换为(伪)代码的数学了解不够。我是一位经验丰富的计算机程序员,但我的数学背景仅此而已。

我将复制我认为是以下答案之一的最重要部分(不幸的是,SO不支持mathjax,所以我不得不复制为图像)

在此处输入图像描述

我遇到的上述困难:

  • 如何求解用于 ψn 的反函数?
  • 给出的示例是 d = 3。如何为任意 d 生成公式?

了解所涉及数学的任何人都能够在链接斐波那契格子问题的任一答案的伪代码实现方面取得进展吗?我知道完整的实现可能非常困难,所以我很高兴部分实现能够引导我走得足够远,能够自己完成其余部分。

为了更容易,我已经编写了一个函数,该函数采用 N 维的球坐标并将它们转换为笛卡尔坐标,因此实现可以输出任何一个,因为我可以轻松转换。

此外,我看到一个答案对每个附加维度使用下一个素数。我可以很容易地编写一个函数来输出每个连续的素数,所以你可以假设它已经实现了。

如果未能在 N 维中实现斐波那契晶格,我很乐意接受满足上述约束的不同方法。

4

5 回答 5

9

非常有趣的问题。我想在我的 4D 渲染引擎中实现它,因为我很好奇它会是什么样子,但是我太懒了,也没有能力从数学方面处理 ND 超越问题。

相反,我想出了不同的解决方案来解决这个问题。它不是斐波那契胶乳!!!相反,我将超球面或 n 球面的参数方程扩展为超螺旋,然后仅拟合螺旋参数,因此这些点或多或少是等距的。

我知道这听起来很可怕,但它并不难,结果对我来说看起来是正确的(最后:) 在解决了一些愚蠢的错别字复制/粘贴错误之后)

主要思想是使用超球面的 n 维参数方程从角度和半径计算其表面点。这里实现:

[edit2]。现在问题归结为两个主要问题:

  1. 计算螺丝数量

    因此,如果我们希望我们的点是等距的,那么它们必须以等距离位于螺旋路径上(参见项目符号#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.0d=d/r稍后再缩放点)。然后扩展(每个维度只是添加角度参数)看起来像这样:

    二维:

    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)您可以将点作为输入值,但要注意它只是近似的点数不准确!!!

  2. 在螺旋上生成台阶,所以点是等距的

    这是我跳过代数混乱并改用拟合的部分。我只是二进制搜索增量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透视:

3D透视

4D透视:

4D透视

具有超平面的 4D 横截面w=0.0

4D 横截面 w=0.0

与更多的点和更大的半径相同:

4D 横截面 w=0.0 HQ

形状随着动画的旋转而变化......

[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 位整数。希望我没有忘记什么...

于 2019-07-28T10:50:17.050 回答
2

作为部分答案,您可以使用牛顿法计算 f 的倒数。在牛顿迭代中使用x作为初始点是一个不错的选择,因为与 .f(x)的距离永远不会超过 1 个单位x。这是一个 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。

于 2019-07-22T11:29:35.123 回答
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 <= D(Pi, Pj)

当且仅当 MD 不能更高时,解决方案才是可接受的。

如果 d = 2,那么我们有一个圆并在其上放置点。圆是具有以下属性的多边形:

  • 它有 n 个角
  • n -> 无穷大
  • 每边的长度相似

所以,一个有 n 个角的多边形,其中 n 是一个有限数并且大于 2,而且,每条边的长度相似,每次我们增加 n 时都更接近一个圆。请注意,d = 2 中的第一个多边形是三角形。我们有一个角度,我们的最小角度单位是 360 度/n。

现在,如果我们有一个正方形并将点均匀分布在上面,那么通过基本变换将正方形转换为圆形应该是精确解,或者非常接近它。如果它是精确解,那么这对于 d = 2 的情况是一个简单的解。如果它只是非常接近,那么通过近似方法,我们可以确定在您选择的给定精度内的解是什么。

我会在 d = 3 的情况下使用这个想法。我会解决一个立方体的问题,这个问题要简单得多,并使用基本变换将我的立方体点转换为我的球体点。我会在 d > 3 上使用这种方法,解决超立方体的问题并将其转换为超球面。当您将点均匀分布在 d 维的超立方体上时,请使用曼哈顿距离。

请注意,我不知道将超立方体转换为超球体的解是精确解还是接近它,但如果它不是精确解,那么我们可以通过近似来提高精度。

因此,这种方法是解决问题的一种方法,就时间复杂度而言,这不一定是最好的方法,因此,如果有人深入研究了斐波那契格区域并知道如何将其推广到更多维度,那么他/她的答案可能是比我更好的接受选择。

如果您定义了 f(x) 的泰勒级数,则可以确定 f(x) = x - 0.5sin2x 的倒数。你会得到一个可以反转的多项式系列 x 。

于 2019-07-22T13:57:56.407 回答
2

之前的所有答案都有效,但仍然缺少实际代码。缺少两个真正的部分,这通常是实现的。

  1. 我们需要计算 的积分sin^(d-2)(x)。如果您按部分进行递归集成,则它具有封闭形式。在这里,我以递归方式实现它,尽管对于维度 ~> 100 我发现数字积分sin^d更快
  2. 我们需要计算该积分的反函数,对于sin^dd > 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 个点生成以下图像:

球体分布

于 2019-12-11T06:03:47.457 回答
2

我对如何做到这一点有了另一个疯狂的想法。它与我以前的方法完全不同,因此有了新的答案......

那么其他答案之一建议在超立方体表面上创建点的均匀分布,然后将点到超立方体中心的距离标准化为超空间半径,并将其用于排斥粒子模拟。我过去在 3D 中这样做,效果很好,但在更高维度上,这将非常缓慢或因 BVH 之类的结构而变得复杂。

但这让我开始思考如何倒退。因此,在超立方体上非线性分布点,因此在归一化后,点在超球面上呈线性分布......

让我们从简单的 2D 开始

二维

因此,我们只需在绿点之间步进角度+/-45 deg并计算绿点。角度步长da必须精确划分90 deg并给出点密度。所以所有的 2D 点都是所有面的+/-1.0组合tan(angle)

当所有点都完成后,只需计算每个点到中心的大小并重新缩放它,使其等于超球面半径。

这可以很容易地扩展到任何维度

2D 以上的每个维度只需为循环角度添加一个即可迭代。

这里使用我之前回答的引擎的 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 4D

这里有更多的点和更好的 3D 尺寸:

3D

立方体图案略微可见,但点距离对我来说看起来不错。很难在 GIF 上看到它,因为后面的点正在与前面的点合并……

这是没有归一化为球体的 2D 正方形和 3D 立方体:

二维正方形 3D立方体

正如您在边缘看到的那样,点密度要小得多......

预览仅使用透视投影,因为这不会生成网格拓扑,仅使用点,因此无法生成横截面...

还要注意这会在边缘产生一些重复的点(我认为对于某些镜子来说,将角度循环少一次迭代应该可以解决这个问题,但懒得实现)

我强烈建议阅读以下内容:

于 2020-07-06T10:57:55.837 回答