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.