c++ - 如果我有一个3d的立方体网格，如何有效地找到所有与半径r和位置p的球面相交的立方体？

Question

想象一下，您在x，y和z方向上的0D和 d 米之间有一个3D空间(这是一个长度为 d 的立方体)。现在假设我想将其分解为较小的多维数据集的3D数组。我们想要每个方向上的和立方体。这意味着每个多维数据集的长度，宽度和高度为 d / n (您可以假设d / n完全除法)。
可以将其保存在以下数组中:Cube* cubes[n][n][n];现在让我们假设我们可以有一个半径为 r 的球体，并且位置 p ，其中 p 可以在3d空间内的任何位置，而 r 没有任何约束。查找和返回与该球面相交的所有多维数据集的最有效方法是什么？
当前，我正在遍历数组中的所有多维数据集，并使用毕达哥拉斯定理检查每个多维数据集与 p 之间的距离，我在下面为其编写了伪代码。

Cube* getIntersectingCubes(r,p) { 
    Cube* cubes[n][n][n];
    Cube* result[];
    for(int i=0;i<n;i++) {
        for(int j=0;j<n;j++) {
            for(int k=0;k<n;k++) {
                if(distanceBetween(cubes[i][j][k],p) < r) {
                    result.add(cubes[i][j][k]);
                }
            }
        }
    }
    return result;
}

但是，绝对有更好的方法可以做到这一点。如何返回相同的输出，但要检查尽可能少的多维数据集？

Answer 1

如果仍然难以确定坐标系和算法，请从图表开始:

本质上，您拥有的是3D多维数据集网格，您可以按照自己喜欢的任何方式选择点p，只要它们完全描述了多维数据集即可。您要确定的是半径r与哪个多维数据集相交的多维数据集的任何部分。您可以使用三个 vector 来处理。您只需将多维数据集原点视为 vector 。您知道每个立方体的边长(d / n)可用于形成第二个 vector (称它为p或p_extent的范围)，可以将其添加到原点 vector 中以确定立方体中的最远点。因此p_extent等于d/n i + d/n j + d/n k 。
您可以创建一个结构以帮助将坐标作为 vector 处理，例如:

typedef struct {
    double x, y, z;
} vect_pos;

使用两个 vector ，您可以简单地将它们加在一起(例如 p + p_extent)以形成一个 vector ，该 vector 描述立方体中距原点最远的点。 (您希望多维数据集内的 vector 的方向与原点 vector 具有相同的符号方向(因此，您可以将其与该方向上的单位 vector 相乘)。例如，您的点是 -3 i -1 j 和 2 k ，您可以将 p_extent vector 乘以 -1 i -1 j 1 k 。 。
要将两个 x, y, z加在一起，您可以简单地使用 vector 加法，例如

/* vector addition */
vect_pos vect_add (const vect_pos *a, const vect_pos *b)
{
    vect_pos c = { .x = 0. };
    
    c.x = a->x + b->x;
    c.y = a->y + b->y;
    c.z = a->z + b->z;
    
    return c;
}

在多维数据集的某个角选择坐标的最后一个皱纹需要您解决的问题是，原点 vector 不会总是指向最接近原点的角(上图中的图片，如果您旋转原点计数器， vector_pos将指向该位置，绕轴顺时针方向旋转90度(您可以选择以立方体的中心作为原点，但是计算的次数是原来的两倍)。
您可以通过简单的检查来处理最接近的角点问题，并将原点的角点移动一个立方体的边长，以计算 p是否在最接近的角和最远的角之间(该立方体仍然是相同的立方体)。原点和描述的最远角，您可以简单地使用vector-distance公式计算每个点到原点的距离，如果 r在两者之间，则球体与该立方体相交。
要在计算距离之前选择最接近原点的点，可以使用:

/* point in cube closest to origin */
vect_pos closest_pt (const vect_pos *p, const int len)
{
    vect_pos c = {  p->x >= 0 ? p->x : p->x + len,
                    p->y >= 0 ? p->y : p->y + len,
                    p->z >= 0 ? p->z : p->z + len  };
    
    return c;
}

(在坐标为正的原点 r处选择多维数据集，可以在坐标为负的情况下简单地添加多维数据集的边长，以选择多维数据集上最接近的点作为原点)
对于距离计算，平方和的简单平方根为:

/* vector distance */
double vect_dist (const vect_pos *p)
{
    return sqrt (p->x * p->x + p->y * p->y + p->z * p->z);
}

要检查 0, 0, 0是否落在每个多维数据集的最近距离与最远距离之间，可以将其与一个简短函数一起使用，该简短函数可以在与原点相同的方向上选择范围 vector ，将 vector 相加，然后将每个距离与 r比较，例如

/* does r intersect cube at postition p */
int intersect (vect_pos *p, const double r, const int len)
{
    vect_pos c = closest_pt (p, len),                   /* closest pt to origin */
            p_extent = { .x = c.x >= 0 ? len : -len,    /* length, unit vector */
                         .y = c.y >= 0 ? len : -len,    /* pointing outward */
                         .z = c.z >= 0 ? len : -len },
             p_farthest = vect_add (&c, &p_extent);     /* farthest point in cube */
    double   nearest  = vect_dist (&c),                 /* distance to nearest */
             farthest = vect_dist (&p_farthest);        /* distance to farthest */
    
    return nearest <= r && r <= farthest;               /* return radius in between */
}

基本上就是这样。您可以添加简短功能以方便地输出点，例如

void prn_pos (const vect_pos *p)
{
    printf ("(% d, % d, % d)\n",
            (int)p->x, (int)p->y, (int)p->z);
}

然后添加一个允许输入 r， main()和 r的 d并输出结果:

int main (int argc, char **argv) {
    
    double  r = argc > 1 ? strtod (argv[1], NULL) : .9;     /* r */
    int     d = argc > 2 ? strtol (argv[2], NULL, 0) : 1,   /* d */
            n = argc > 3 ? strtol (argv[3], NULL, 0) : d,   /* n */
            dn = 0,             /* d/n */
            nc = 0,             /* number of cubes in each direction */
            total = 0,          /* total cubes in search grid */
            intersecting = 0;   /* number of cubes that intersect sphere */
    
    if (r < 0 || d < 0) {   /* validate length inputs */
        fputs ("error: negtive length in input of r or d.\n", stderr);
        return 1;
    }
    if (!n || n < 0) {      /* validate n not zero or negative */
        fputs ("error: n - divide by zero or negative.\n", stderr);
        return 1;
    }
    if (d < r) {            /* ensure d >= r */
        int min_d = ceil (r);
        n *= min_d / d;
        d = min_d;
    }
    if (d % n) {            /* validate d equally divisible by n */
        fputs ("error: d not equally divisible by n.\n", stderr);
        return 1;
    }
    
    dn = d / n;                     /* cube side length */
    nc = ceil (r / dn) + 1;         /* limit of cubes needed per-direction */
    total = nc * nc * nc * 8;       /* total number of cubes in search grid */
    
    printf ("\nOut of %d cubes, the following intersect sphere of radius %.2f\n\n",
            total, r);
    
    for (int i = -d; i <= d; i += dn)               /* loop -d to d */
        for (int j = -d; j <= d; j += dn)
            for (int k = -d; k <= d; k += dn) {
                vect_pos p = { i, j, k };           /* vector to cube origin */
                if (intersect (&p, r, dn)) {        /* does it intersect sphere? */
                    prn_pos (&p);                   /* output cube origin */
                    intersecting++;                 /* increment count */
                }
            }
    
    printf ("\nintersecting cubes: %d\n", intersecting);
}

添加标题:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

然后针对数学库 n编译链接，您就有一个可以测试的程序。如果未提供 -lm， r或 d的参数，则该代码将 n， .90和 1作为默认值。请注意，在每个维度上，循环限制从 1到 -(next int above r)至少是最小的，或者由 (next int above r)给出，以大于且不小于 d的限制。
示例使用/输出
带有 r和 r = .9和 d = 1的默认情况:

$ ./bin/cube-sphere-intersect

Out of 64 cubes, the following intersect sphere of radius 0.90

(-1, -1, -1)
(-1, -1,  0)
(-1,  0, -1)
(-1,  0,  0)
( 0, -1, -1)
( 0, -1,  0)
( 0,  0, -1)
( 0,  0,  0)

intersecting cubes: 8

如果在网格上绘制多维数据集，则这些是围绕原点的8个多维数据集。
选择 n = 1， r = 5.5和 d = 6会得到相同的结果，除了立方体的边长不同:

$ ./bin/cube-sphere-intersect 5.5 6 1

Out of 64 cubes, the following intersect sphere of radius 5.50

(-6, -6, -6)
(-6, -6,  0)
(-6,  0, -6)
(-6,  0,  0)
( 0, -6, -6)
( 0, -6,  0)
( 0,  0, -6)
( 0,  0,  0)

intersecting cubes: 8

如果您输入的 n = 1长度实际上是 r，那么您将在此处选择是否将 1视为等于最小或最大长度。有论点说“接触”一个多维数据集并不一定意味着它相交，因为该点在技术上将位于相距一定距离的切线上，但是同样有很好的论据表明该相交(选择是否要使用 r或比较中的 <=留给您)
例如，如果您在 <( 1.0 < sqrt(2))之间输入任何内容，则会发现32个立方体相交，例如

$ ./bin/cube-sphere-intersect 1.1

Out of 216 cubes, the following intersect sphere of radius 1.10

(-2, -1, -1)
(-2, -1,  0)
(-2,  0, -1)
(-2,  0,  0)
(-1, -2, -1)
(-1, -2,  0)
(-1, -1, -2)
(-1, -1, -1)
(-1, -1,  0)
(-1, -1,  1)
(-1,  0, -2)
(-1,  0, -1)
(-1,  0,  0)
(-1,  0,  1)
(-1,  1, -1)
(-1,  1,  0)
( 0, -2, -1)
( 0, -2,  0)
( 0, -1, -2)
( 0, -1, -1)
( 0, -1,  0)
( 0, -1,  1)
( 0,  0, -2)
( 0,  0, -1)
( 0,  0,  0)
( 0,  0,  1)
( 0,  1, -1)
( 0,  1,  0)
( 1, -1, -1)
( 1, -1,  0)
( 1,  0, -1)
( 1,  0,  0)

intersecting cubes: 32

然后对于 1.414...之间的半径，您会发现56个立方体相交，依此类推。
这只是实现它的一种方法。如果您好奇，请重写它以将每个多维数据集的中心用作多维数据集的原点。您不必使用结构，您可以直接对每个 sqrt(2) <= sqrt(3)进行操作，但是该结构有助于使逻辑保持清晰。仔细检查一下，如果您还有其他问题，请告诉我。