改进 C++ 算法以查找半径为 r 的球体内的所有点

Question

语言/编译器：C++ (Visual Studio 2013) 经验：~2 个月

我在 3D 空间中的矩形网格中工作（大小：xdim by ydim by zdim）其中，“xgrid、ygrid 和 zgrid”分别是 x、y 和 z 坐标的 3D 数组。现在，我有兴趣找到以点“(vi,vj,vk)”为中心的半径为“r”的球体内的所有点。我想将这些点的索引位置存储在向量“xidx,yidx,zidx”中。对于单个点，该算法可以工作并且足够快，但是当我希望迭代 3D 空间内的许多点时，我会遇到很长的运行时间。

有没有人对我如何改进这个算法在 C++ 中的实现有任何建议？运行一些分析软件后，我在网上找到了（非常困，Luke stackwalker），似乎“std::vector::size”和“std::vector::operator[]”成员函数让我的代码陷入困境。非常感谢任何帮助。

注意：由于我不知道球体内有多少体素，因此我将向量 xidx,yidx,zidx 的长度设置为大于必要的长度，然后在函数末尾擦除所有多余的元素。

void find_nv(int vi, int vj, int  vk, vector<double> &xidx, vector<double> &yidx, vector<double> &zidx, double*** &xgrid, double*** &ygrid, double*** &zgrid, int r, double xdim,double ydim,double zdim, double pdim)

{
double xcor, ycor, zcor,xval,yval,zval;
vector<double>xyz(3);
xyz[0] = xgrid[vi][vj][vk];
xyz[1] = ygrid[vi][vj][vk];
xyz[2] = zgrid[vi][vj][vk];
int counter = 0;

// Confine loop to be within boundaries of sphere
int istart = vi - r;
int iend = vi + r;
int jstart = vj - r;
int jend = vj + r;
int kstart = vk - r;
int kend = vk + r;

if (istart < 0) {
    istart = 0;
}
if (iend > xdim-1) {
    iend = xdim-1;
}
if (jstart < 0) {
    jstart = 0;
}
if (jend > ydim - 1) {
    jend = ydim-1;
}
if (kstart < 0) {
    kstart = 0;
}
if (kend > zdim - 1)
    kend = zdim - 1;

//-----------------------------------------------------------
 // Begin iterating through all points
//-----------------------------------------------------------
for (int k = 0; k < kend+1; ++k)
{
    for (int j = 0; j < jend+1; ++j)
    {
        for (int i = 0; i < iend+1; ++i)
        {
            if (i == vi && j == vj && k == vk)
                continue;
            else
            {
            xcor = pow((xgrid[i][j][k] - xyz[0]), 2);
            ycor = pow((ygrid[i][j][k] - xyz[1]), 2);
            zcor = pow((zgrid[i][j][k] - xyz[2]), 2);
            double rsqr = pow(r, 2);
            double sphere = xcor + ycor + zcor;
            if (sphere <= rsqr)
            {
                xidx[counter]=i;
                yidx[counter]=j;
                zidx[counter] = k;
                counter = counter + 1;
            }

            else
            {
            }
            //cout << "counter = " << counter - 1;
            }
        }
    }
}

// erase all appending zeros that are not voxels within sphere
xidx.erase(xidx.begin() + (counter), xidx.end()); 
yidx.erase(yidx.begin() + (counter), yidx.end()); 
zidx.erase(zidx.begin() + (counter), zidx.end()); 

return 0;

Answer 1

您似乎已经在这类事情上使用了我最喜欢的技巧，摆脱了相对昂贵的平方根函数，只使用半径和中心到点距离的平方值。

另一种可能加速^(a)的可能性是替换所有：

xyzzy = pow (plugh, 2)

调用更简单：

xyzzy = plugh * plugh

您可能会发现删除函数调用可以加快处理速度，但幅度很小。

如果您可以确定目标数组的最大大小，另一种可能性是使用实数数组而不是向量。我知道他们使矢量代码尽可能地优化，但它仍然无法匹配固定大小的数组以提高性能（因为它必须完成固定大小数组所做的所有事情 plus 处理可能的扩展） .

同样，这可能只会以更多内存使用为代价提供非常微小的改进，但以空间换时间是一种经典的优化策略。

除此之外，请确保您明智地使用编译器优化。在大多数情况下，默认构建具有低级别的优化，以使调试更容易。将其升级为生产代码。

---

^(a) 与所有优化一样，您应该衡量，而不是猜测！这些建议就是：建议。它们可能会也可能不会改善这种情况，因此由您来测试它们。

Answer 2

您最大的问题之一，并且可能阻止编译器进行大量优化的问题是您没有使用网格的常规特性。

如果你真的使用常规网格，那么

xgrid[i][j][k] = x_0 + i * dxi + j * dxj + k * dxk
ygrid[i][j][k] = y_0 + i * dyi + j * dyj + k * dyk
zgrid[i][j][k] = z_0 + i * dzi + j * dzj + k * dzk

如果你的网格是轴对齐的，那么

xgrid[i][j][k] = x_0 + i * dxi
ygrid[i][j][k] = y_0 + j * dyj
zgrid[i][j][k] = z_0 + k * dzk

在核心循环中替换这些应该会显着提高速度。

Answer 3

你可以做两件事。减少您正在测试包含的点数，并将问题简化为多个 2d 测试。

如果你沿着 z 轴看球体，你会在球体中拥有 y+r 到 yr 的所有点，使用这些点中的每一个，你可以将球体切割成包含所有点的圆圈x/z 平面限制在您正在测试的特定 y 处的圆半径。计算圆的半径是解决直角三角形底边长度问题的简单方法。

现在您正在测试一个立方体中的所有点，但球体的上部范围不包括大多数点。上述算法背后的想法是，您可以将在球体的每一层测试的点限制在包含该高度处圆半径的正方形内。

这是一个简单的手绘图形，从侧面显示球体。

在这里，我们正在查看半径为 ab 的球体切片。由于您知道直角三角形的长度 ac 和 bc，您可以使用毕达哥拉斯定理计算 ab。现在你有一个简单的圆圈，你可以测试其中的点，然后向下移动，它减少长度 ac 并重新计算 ab 并重复。

现在，一旦你有了它，你实际上可以做更多的优化。首先，您不需要针对圆圈测试每个点，您只需要测试四分之一的点。如果您测试圆左上象限（球体的切片）中的点，那么其他三个点中的点只是同一点的镜像，该点从确定的点向右、底部或对角偏移位于第一象限。

最后，您只需要对球体的上半部分进行圆形切片，因为下半部分只是上半部分的镜像。最后，您只测试了球体中四分之一的遏制点。这应该是一个巨大的性能提升。

我希望这是有道理的，我现在不在机器旁，我可以提供样品。

Answer 4

这里的简单事情是从球体中心进行 3D 泛洪填充，而不是在需要访问较小点时迭代封闭的正方形。此外，您应该实现洪水填充的迭代版本以提高效率。

Flood Fill