Compute gradient for voxel data efficiently

Question

What is the most efficient way of computing the gradient for fixed sized voxel data, such as the source code below. Note that I need the gradient at any point in space. The gradients will be used for estimating normals in a marching cubes implementation.

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

};

Update 1 A demo project of the problem can be found here:

https://github.com/mortennobel/OpenGLVoxelizer

Currently the output is like the picture below (based on MooseBoys code):

Update 2 The solution that I'm looking for must give fairly accurate gradients, since they are used as normals in a visualisation and visual artefacts like the ones below must be avoided.

enter image description here

Update 2 Solution from the user example is:

enter image description here

I'm interested why you need the gradient of any point. It seems this would be any easy chance for optimisation. Also the regular spacing of points (as opposed to typical shading where triangle size varies) may lend itself to better approaches than get_gradient point by point. Have you looked at such a possibility? — user3125280, Jan 24 '14 at 04:00
There is a potential performance optimization, since I only need gradient at every vertex generated by marching cubes and these vertices are always have two of the axis located at an integer number (in other words this could simplify the calculations to two linear interpolations and one trilinear interpolation). — Mortennobel, Jan 27 '14 at 02:06
if you have the gradient (which is a vector [x,y,z]) at each grid point, then trilinear interpolation (of the gradient) would become linear interpolation along one axis (the non-integral one). On the other hand you could calculate the gradient as the derivative of the interpolated density field. I don't know what pros/cons that has, but it will be slower. Which way were you intending? — user3125280, Jan 27 '14 at 03:07
Where are you getting your data from? I've got amazingly smooth results before by using the original source of the data the voxels were generated with instead of the voxels themselves. Provided that the data is in a form which you can differentiate, such as a sum of smooth fields, this gives beautiful results. — Andy Newman, Jan 28 '14 at 11:37
@AndyNewman: The voxel data I work on is the result of a FEM (finite element method) based computation (more specifically multigrid topology optimization in 3D). AFAIK the gradient the result must be computed based on the voxels using an interpolation schema. My main concern is how to get the best visual result using a reasonable fast approach (I actually don't think performance will be a huge issue). — Mortennobel, Jan 28 '14 at 13:37

score 2 · Accepted Answer · answered Jan 24 '14 at 02:46

2

The following produces a linearly interpolated gradient field:

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

answered Jan 24 '14 at 02:46

MooseBoys

6,641
1
19
43

Isn't there a more accurate way of estimate the gradients? I tried to implement your suggestion, but was quite disappointed with the quality of the gradients found. – Mortennobel Jan 28 '14 at 08:23
1

@Mortennobel What is your definition of accuracy? There are a number of ways to interpolate values between discrete data points, linear being a common one. Another is cubic interpolation, which produces smoother results (both first and second derivatives are continuous), but causes overshoot. – MooseBoys Jan 28 '14 at 09:19
Well I'm just looking for a method that gives a visually pleasing result (a more smooth result that is). – Mortennobel Jan 28 '14 at 09:39
2

@Mortennobel You asked for an efficient way not for an accurate one (or visually smoother as you call it), Probably you can add that as a comment to the original question? – Sergio Ayestarán Jan 28 '14 at 19:40

score 2 · Answer 2 · answered Jan 30 '14 at 14:37

2

One important technique for optimization in many implementations involves time/space trade off. As a suggestion, anywhere you can pre-calc and cache your results may be worth looking at.

answered Jan 30 '14 at 14:37

Bradley Thomas

4,060
6
33
55

score 2 · Answer 3 · answered Aug 28 '14 at 14:52

In general Sobel filters provide slightly nicer results than simple central tendency, but take longer to compute (the Sobel is essentially a smooth filter combined with central tendency). A classic Sobel requires weighting 26 samples, while central tendency only requires 6. However, there is a trick: with GPUs you can get hardware-based trilinear interpolation for free. That means you can compute a Sobel with 8 texture reads, and this can be done in parallel across the GPU. The following page illustrates this technique using GLSL http://www.mccauslandcenter.sc.edu/mricrogl/notes/gradients For your project you would probably want to compute the gradients on the GPU and use GPGPU methods to copy the results back from the GPU to the CPU for further processing.

example · Answer 4 · 2014-01-30T20:09:27.277

MooseBoys already posted a component-wise linear interpolation. It is discontinuous in the y and z component though, whereever (int)x changes from one value to the next (same thing for the other components). This might cause such a rough picture as you are seeing it. If you have enough performance to spare you can improve this by considering not just (int)x but (int)(x+1) aswell. This might look like the following

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

This can probably be written a bit more concise, but maybe this way you can still see what is happening. If this still is not smooth enough for you will have to look into cubic / spline interpolation or similar.

I implemented your suggestion, but the result looks quite bad. (See update to the question for a screenshot). — Mortennobel, Jan 29 '14 at 03:32
@Mortennobel uff. It should be at least as good as mooseboys version. maybe a sign error of some sorts. let me review my code once i'm home from work. — example, Jan 29 '14 at 10:06
@Mortennobel There was a "-1" missing in the calculation of yd0. unfortunately i cannot compile your code without too much effort - but I cannot find any other mistakes in the code. hope this was the mistake that caused the artifacts... — example, Jan 30 '14 at 21:19
That fixed the artifacts. However I don't find the visual result better than the solution suggested by MooseBoys. — Mortennobel, Jan 31 '14 at 01:54

Compute gradient for voxel data efficiently

4 Answers4