5

对于固定大小的体素数据，计算梯度的最有效方法是什么，例如下面的源代码。请注意，我需要空间中任何点的渐变。梯度将用于估计行进立方体实现中的法线。

#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 解决方案是：

在此处输入图像描述

4

4 回答 4

2

下面产生一个线性插值梯度场：

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;
}
于 2014-01-24T02:46:18.300 回答
2

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

于 2014-01-30T14:37:15.837 回答
2

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

于 2014-08-28T14:52:23.780 回答
1

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;
}

这可能可以写得更简洁一些，但也许这样你仍然可以看到正在发生的事情。如果这仍然不够平滑，您将不得不研究三次/样条插值或类似的。

于 2014-01-28T18:47:55.957 回答