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.

Update 2 Solution from the user example is:

Answer 1

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