python - 球体表面点的均匀分布?

标签 python c++ algorithm math geometry

<分区>

我正在寻找一种算法来优化球体表面上的 N 个点分布(不是随机的,我说的是 Tammes problem 意义上的最优分布)。 我想在 C++ 中这样做,但问题不是语言问题,而是算法问题。这是我目前使用的算法(在 python 中):

# Function to distribute N points on the surface of a sphere 
# (source: http://www.softimageblog.com/archives/115)
def uniform_spherical_distribution(N): 
    pts = []   
    inc = math.pi * (3 - math.sqrt(5)) 
    off = 2 / float(N) 
    for k in range(0, int(N)): 
        y = k * off - 1 + (off / 2) 
        r = math.sqrt(1 - y*y) 
        phi = k * inc 
        pts.append([math.cos(phi)*r, y, math.sin(phi)*r])   
    return pts

对于大量的点,视觉结果似乎相当不错,但对于少量的点(2 到 10 之间),结果似乎一点也不好(半球中有更多点)。有没有一种方法可以增强算法,使其在点数较少的情况下也可以...(答案可以是 C、C++ 或 Python)。

最佳答案

如果您只需要较少的点数,蛮力可能会奏效。我所说的蛮力是指预先计算 2..10 的解决方案。如果 N<=10,则在数组中查找它们.

对于小N就像那样,应该可以解决数学优化问题,以最佳方式分配您的点数。你需要 multistart 但对于小 N这应该不是问题。这是我在 AMPL 中的尝试:

param N;

var x{1..N};
var y{1..N};
var z{1..N};

var d;

param xbest{1..N};
param ybest{1..N};
param zbest{1..N};
param dbest;

maximize obj: d;

all_points_on_the_sphere{i in 1..N}:
  x[i]^2 + y[i]^2 + z[i]^2 = 1;

all_pairs{i in 1..N, j in 1..N : i<j}:
  (x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2 >= d;

fix_first_x:  x[1] = 1;
fix_first_y:  y[1] = 0;
fix_first_z:  z[1] = 0;  
fix_second_z: z[2] = 0;

#############################################

option show_stats 1;
option presolve 10;
option substout 1;
option var_bounds 2;
#option nl_comments 0;
#option nl_permute 0;
option display_precision 0;

option solver "/home/ali/ampl/ipopt";

for {k in 2..10} {
    let N := k;  
    solexpand _con;
    let dbest := -1.0;
    # Multistart
    for {1..2000} {
        for {j in 1.._snvars}
            let _svar[j] := Uniform(-1, 1);
        let d := N;     

        solve;

        if (solve_result_num < 200 and d > dbest) then {
            let dbest := d;
            for {i in 1..N} {
                let xbest[i] := x[i];
                let ybest[i] := y[i];
                let zbest[i] := z[i];
            }
        }
    }

    print "@@@";
    printf "N=%d, min distance %6f\n", N, sqrt(dbest);

    for {i in 1..N}
        printf "(%9.6f, %9.6f, %9.6f)\n", xbest[i], ybest[i], zbest[i];
}

在我的机器上运行这个脚本花了 5 分钟,解决方案是:

N=2, min distance 2.000000
( 1.000000,  0.000000,  0.000000)
(-1.000000,  0.000000,  0.000000)

N=3, min distance 1.732051
( 1.000000,  0.000000,  0.000000)
(-0.500000,  0.866025,  0.000000)
(-0.500000, -0.866025,  0.000000)

N=4, min distance 1.632993
( 1.000000,  0.000000,  0.000000)
(-0.333333, -0.942809,  0.000000)
(-0.333333,  0.471405, -0.816497)
(-0.333333,  0.471405,  0.816497)

N=5, min distance 1.414214
( 1.000000,  0.000000,  0.000000)
(-0.208840,  0.977950,  0.000000)
(-0.000000,  0.000000,  1.000000)
(-0.212683, -0.977121,  0.000000)
( 0.000000,  0.000000, -1.000000)

N=6, min distance 1.414214
( 1.000000,  0.000000,  0.000000)
(-1.000000,  0.000000,  0.000000)
( 0.000000, -0.752754, -0.658302)
( 0.000000,  0.752754,  0.658302)
( 0.000000,  0.658302, -0.752754)
( 0.000000, -0.658302,  0.752754)

N=7, min distance 1.256870
( 1.000000,  0.000000,  0.000000)
(-0.688059, -0.725655,  0.000000)
( 0.210138, -0.488836, -0.846689)
( 0.210138, -0.488836,  0.846688)
(-0.688059,  0.362827,  0.628435)
(-0.688059,  0.362827, -0.628435)
( 0.210138,  0.977672,  0.000000)

N=8, min distance 1.215563
( 1.000000,  0.000000,  0.000000)
( 0.261204, -0.965284,  0.000000)
( 0.261204,  0.565450,  0.782329)
(-0.783612, -0.482642, -0.391165)
( 0.261204, -0.199917, -0.944355)
( 0.261204,  0.882475, -0.391165)
(-0.783612,  0.599750,  0.162026)
(-0.477592, -0.399834,  0.782329)

N=9, min distance 1.154701
( 1.000000,  0.000000,  0.000000)
(-0.500000, -0.866025,  0.000000)
( 0.333333, -0.577350, -0.745356)
(-0.500000,  0.866025,  0.000000)
(-0.666667, -0.000000,  0.745356)
(-0.666667,  0.000000, -0.745356)
( 0.333333, -0.577350,  0.745356)
( 0.333333,  0.577350, -0.745356)
( 0.333333,  0.577350,  0.745356)

N=10, min distance 1.091426
( 1.000000,  0.000000,  0.000000)
(-0.605995,  0.795469,  0.000000)
( 0.404394,  0.816443,  0.412172)
(-0.664045, -0.746251, -0.046407)
( 0.404394, -0.363508, -0.839242)
(-0.664045,  0.002497,  0.747688)
(-0.605995,  0.046721, -0.794096)
( 0.404394, -0.908368,  0.106452)
( 0.255933,  0.703344, -0.663179)
( 0.404394, -0.159620,  0.900548)

通过查看数字,很明显可以通过分析计算某些解决方案(我认识 sqrt(2) 和 sqrt(3) 等)。我相信 N=2、4 和 6 的解是直线 ([-1,0,0], [1,0,0])、四面体、八面体。

没有strong保证以上是点的最佳分布。非线性求解器可能陷入局部最优;作为N增长,局部最优的数量也增长。

您可以将上述解决方案放在一个数组中,然后在 Python、C++ 或您使用的任何语言中使用它们。

关于python - 球体表面点的均匀分布?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21977328/

相关文章:

python - Django - 以一种形式重复n次表单字段

c++ - 间接转发引用

c++ - 查找 x 轴上最大非重叠线数的算法

algorithm - 时间表在有限时间内的最低成本?

python - 在 macOS Big Sur 上将 Python 3.7 升级到 3.9

python - 为什么将此函数应用于未作为参数调用的变量?

python - Django 2.0 使用 views.py 呈现 html 页面 (TemplateDoesNotExist)

c++ - 将字符串拆分为字符串 vector 时出现无效转换错误

c++ - 在 Visual C++ 和 clang 中使用 C++11 unordered_set

algorithm - 从包含在集合 D 中的集合 C 的集合中有效地找到所有集合 S