高效计算体素数据的梯度

Question

如何高效计算固定大小的体素数据（例如下面的源代码）的梯度？请注意，我需要在任意空间点处计算梯度。这些梯度将用于在Marching Cubes实现中估计法线。

#import <array>

struct VoxelData {
    VoxelData(float* data, unsigned int xDim, unsigned int yDim, unsigned int zDim)
    :data(data), xDim(xDim), yDim(yDim), zDim(zDim)
    {}

    std::array<float,3> get_gradient(float x, float y, float z){
        std::array<float,3> res;
        // compute gradient efficiently
        return res;
    }

    float get_density(int x, int y, int z){
        if (x<0 || y<0 || z<0 || x >= xDim || y >= yDim || z >= zDim){
            return 0;
        }
        return data[get_element_index(x, y, z)];
    }

    int get_element_index(int x, int y, int z){
        return x * zDim * yDim + y*zDim + z;
    }

    const float* const data;

    const unsigned int xDim;
    const unsigned int yDim;
    const unsigned int zDim;

};

更新1 该问题的演示项目可以在此处找到：

https://github.com/mortennobel/OpenGLVoxelizer

目前的输出结果如下图所示（基于MooseBoys代码）：

更新2 我正在寻找的解决方案必须给出相当精确的梯度，因为它们用作可视化中的法线，并且必须避免像下面这样的视觉伪影。

更新2 用户示例中的解决方案是：

- Mortennobel

---

我很感兴趣为什么你需要任意点的梯度。这似乎是一个很容易优化的机会。此外，点的常规间距（与典型阴影不同，三角形大小会变化）可能比逐点获取梯度更适合采用更好的方法。你有考虑过这种可能性吗？ - user3125280

存在一种潜在的性能优化，因为我只需要在由 Marching Cubes 生成的每个顶点处计算梯度，并且这些顶点始终具有两个坐标轴位于整数位置（换句话说，这可以简化为两个线性插值和一个三线性插值的计算）。 - Mortennobel

如果在每个网格点上都有梯度（即一个向量[x，y，z]），那么三线性插值（梯度）将沿着一个轴（非整数轴）变为线性插值。另一方面，您可以将梯度计算为插值密度场的导数。我不知道这种方法的优缺点，但速度会更慢。您打算采用哪种方式？ - user3125280

你的数据从哪里来？我之前使用数据的原始来源而不是体素本身生成的数据，得到了非常平滑的结果。只要数据以可以区分的形式存在，比如光滑场的总和，这会得到很好的结果。 - Andy Newman

@AndyNewman：我处理的体素数据是基于有限元方法（FEM）计算的结果（更具体地说，是3D多重网格拓扑优化）。据我所知，必须使用插值模式基于体素计算结果来计算梯度。我的主要关注点是如何使用合理快速的方法获得最佳的视觉效果（实际上，我认为性能不会是一个巨大的问题）。 - Mortennobel

Answer 1

以下代码生成线性插值渐变场：

std::array<float,3> get_gradient(float x, float y, float z){
    std::array<float,3> res;
    // x
    int xi = (int)(x + 0.5f);
    float xf = x + 0.5f - xi;
    float xd0 = get_density(xi - 1, (int)y, (int)z);
    float xd1 = get_density(xi, (int)y, (int)z);
    float xd2 = get_density(xi + 1, (int)y, (int)z);
    res[0] = (xd1 - xd0) * (1.0f - xf) + (xd2 - xd1) * xf; // lerp
    // y
    int yi = (int)(y + 0.5f);
    float yf = y + 0.5f - yi;
    float yd0 = get_density((int)x, yi - 1, (int)z);
    float yd1 = get_density((int)x, yi, (int)z);
    float yd2 = get_density((int)x, yi + 1, (int)z);
    res[1] = (yd1 - yd0) * (1.0f - yf) + (yd2 - yd1) * yf; // lerp
    // z
    int zi = (int)(z + 0.5f);
    float zf = z + 0.5f - zi;
    float zd0 = get_density((int)x, (int)y, zi - 1);
    float zd1 = get_density((int)x, (int)y, zi);
    float zd2 = get_density((int)x, (int)y, zi + 1);
    res[2] = (zd1 - zd0) * (1.0f - zf) + (zd2 - zd1) * zf; // lerp
    return res;
}

- MooseBoys

---

有没有更准确的估计梯度的方法？我尝试了你的建议，但是发现找到的梯度质量相当令人失望。 - Mortennobel

@Mortennobel 你对准确性的定义是什么？在离散数据点之间插值有许多方法，其中一种常见的是线性插值。另一种是三次样条插值，它产生更平滑的结果（一阶和二阶导数都是连续的），但会导致过冲。 - MooseBoys

我只是在寻找一种可以产生视觉效果良好的方法（更加平滑的结果）。 - Mortennobel

@Mortennobel，你要求的是一种高效的方式而不是准确的方式（或者像你所说的视觉上更平滑的方式），也许你可以将其作为评论添加到原始问题中？ - Sergio Ayestarán

Answer 2

通常情况下，Sobel滤波器提供的结果比简单的中心趋势要好一些，但计算时间更长（Sobel本质上是一个平滑滤波器和中心趋势的组合）。经典的Sobel需要加权26个样本，而中心趋势只需要6个。然而，有一个技巧：使用GPU可以免费获得基于硬件的三线性插值。这意味着您可以使用8个纹理读取计算Sobel，并且可以在GPU上并行进行。以下页面使用GLSL说明了这种技术：http://www.mccauslandcenter.sc.edu/mricrogl/notes/gradients对于您的项目，您可能希望在GPU上计算梯度，并使用GPGPU方法将结果从GPU复制回CPU以进行进一步处理。

- user1677899

Answer 3

在许多实现中，优化的一个重要技术涉及时间/空间权衡。作为建议，您可以预先计算和缓存结果的任何地方都值得考虑。

- Bradley Thomas

Answer 4

MooseBoys已经发布了一个分量线性插值。然而，在y和z分量中，无论(int)x何时从一个值变为另一个值（其他分量也是如此），它都是不连续的。这可能会导致您看到的图片很粗糙。如果您有足够的性能可以改进它，考虑不仅考虑(int)x，还要考虑(int)(x+1)。这可能看起来像下面这样：

std::array<float,3> get_gradient(float x, float y, float z){
    std::array<float,3> res;

    int xim = (int)(x + 0.5f);
    float xfm = x + 0.5f - xi;
    int yim = (int)(y + 0.5f);
    float yfm = y + 0.5f - yi;
    int zim = (int)(z + 0.5f);
    float zfm = z + 0.5f - zi;
    int xi = (int)x;
    float xf = x - xi;
    int yi = (int)y;
    float yf = y - yi;
    int zi = (int)z;
    float zf = z - zi;


    float xd0 = yf*(          zf *get_density(xim - 1, yi+1, zi+1) 
                    + (1.0f - zf)*get_density(xim - 1, yi+1, zi))
                +(1.0f - yf)*(zf *get_density(xim - 1, yi  , zi+1) 
                    + (1.0f - zf)*get_density(xim - 1, yi  , zi));
    float xd1 = yf*(          zf *get_density(xim,     yi+1, zi+1) 
                    + (1.0f - zf)*get_density(xim,     yi+1, zi))
                +(1.0f - yf)*(zf *get_density(xim,     yi  , zi+1) 
                    + (1.0f - zf)*get_density(xim,     yi  , zi));
    float xd2 = yf*(          zf *get_density(xim + 1, yi+1, zi+1) 
                    + (1.0f - zf)*get_density(xim + 1, yi+1, zi))
                +(1.0f - yf)*(zf *get_density(xim + 1, yi  , zi+1) 
                    + (1.0f - zf)*get_density(xim + 1, yi  , zi));
    res[0] = (xd1 - xd0) * (1.0f - xfm) + (xd2 - xd1) * xfm;

    float yd0 = xf*(          zf *get_density(xi+1, yim-1, zi+1) 
                    + (1.0f - zf)*get_density(xi+1, yim-1, zi))
                +(1.0f - xf)*(zf *get_density(xi  , yim-1, zi+1) 
                    + (1.0f - zf)*get_density(xi  , yim-1, zi));
    float yd1 = xf*(          zf *get_density(xi+1, yim  , zi+1) 
                    + (1.0f - zf)*get_density(xi+1, yim  , zi))
                +(1.0f - xf)*(zf *get_density(xi  , yim  , zi+1) 
                    + (1.0f - zf)*get_density(xi  , yim  , zi));
    float yd2 = xf*(          zf *get_density(xi+1, yim+1, zi+1) 
                    + (1.0f - zf)*get_density(xi+1, yim+1, zi))
                +(1.0f - xf)*(zf *get_density(xi  , yim+1, zi+1) 
                    + (1.0f - zf)*get_density(xi  , yim+1, zi));
    res[1] = (yd1 - yd0) * (1.0f - yfm) + (yd2 - yd1) * yfm;

    float zd0 = xf*(          yf *get_density(xi+1, yi+1, zim-1) 
                    + (1.0f - yf)*get_density(xi+1, yi  , zim-1))
                +(1.0f - xf)*(yf *get_density(xi,   yi+1, zim-1) 
                    + (1.0f - yf)*get_density(xi,   yi  , zim-1));
    float zd1 = xf*(          yf *get_density(xi+1, yi+1, zim) 
                    + (1.0f - yf)*get_density(xi+1, yi  , zim))
                +(1.0f - xf)*(yf *get_density(xi,   yi+1, zim) 
                    + (1.0f - yf)*get_density(xi,   yi  , zim));
    float zd2 = xf*(          yf *get_density(xi+1, yi+1, zim+1) 
                    + (1.0f - yf)*get_density(xi+1, yi  , zim+1))
                +(1.0f - xf)*(yf *get_density(xi,   yi+1, zim+1) 
                    + (1.0f - yf)*get_density(xi,   yi  , zim+1));
    res[2] = (zd1 - zd0) * (1.0f - zfm) + (zd2 - zd1) * zfm;
    return res;
}

这段内容可以更简洁地表达，但也许这样你仍然可以看出正在发生的事情。如果这对你来说仍然不够流畅，那么你需要研究三次/样条插值或类似方法。

- example

---

我实现了你的建议，但结果看起来相当糟糕。（请参见问题更新中的屏幕截图。） - Mortennobel

@Mortennobel 呜。它应该至少和mooseboys的版本一样好。也许是某种符号错误。等我下班回家后再检查我的代码。 - example

@Mortennobel 在 yd0 的计算中缺少了一个“-1”。不幸的是，我不能轻易地编译您的代码，但我在代码中找不到其他错误。希望这就是导致伪影的错误。 - example

那样修复了瑕疵。然而，我并没有发现视觉效果比MooseBoys建议的解决方案更好。 - Mortennobel