     cDs[1] += getNormalizedVoxel(volume, volumeParams, samplePos - 2.0*volumeParams.textureSpaceGradientSpacing[1]);
     cDs[2] += getNormalizedVoxel(volume, volumeParams, samplePos - 2.0*volumeParams.textureSpaceGradientSpacing[2]);
     // f' = ( -f(x+2h)+8.f(x+h)-8.f(x-h)+f(x-2h) ) / 12*h
-
-    cDs[0] /= 12.0*volumeParams.worldSpaceGradientSpacing;
-    cDs[1] /= 12.0*volumeParams.worldSpaceGradientSpacing;
-    cDs[2] /= 12.0*volumeParams.worldSpaceGradientSpacing;
-    return transpose(cDs);
+    mat4x3 gradients = transpose(cDs);
+    gradients[0] /= 12.0*volumeParams.worldSpaceGradientSpacing;
+    gradients[1] /= 12.0*volumeParams.worldSpaceGradientSpacing;
+    gradients[2] /= 12.0*volumeParams.worldSpaceGradientSpacing;
+    gradients[3] /= 12.0*volumeParams.worldSpaceGradientSpacing;
+    return gradients;
 }
 
 // Use pre-computed gradients stored in xyz of the current sample. The gradient will be transformed from texture space

    fDs[2] = getNormalizedVoxel(volume, volumeParams, samplePos - volumeParams.textureSpaceGradientSpacing[2]);
    // f' = ( f(x)-f(x-h) ) / h
    fDs = mat3x4(intensity, intensity, intensity) - fDs;
    fDs[0] /= volumeParams.worldSpaceGradientSpacing;
    fDs[1] /= volumeParams.worldSpaceGradientSpacing;
    fDs[2] /= volumeParams.worldSpaceGradientSpacing;
    return transpose(fDs);
}

// Compute world space gradient using higher order central difference: f' = ( -f(x+2h)+8.f(x+h)-8.f(x-h)+f(x-2h) ) / 12*h

    cDs[1] -= getNormalizedVoxel(volume, volumeParams, samplePos - volumeParams.textureSpaceGradientSpacing[1]);
    cDs[2] -= getNormalizedVoxel(volume, volumeParams, samplePos - volumeParams.textureSpaceGradientSpacing[2]);
    // f' = ( f(x+h)-f(x-h) ) / 2*h
    mat4x3 gradients = transpose(cDs);
    gradients[0] /= 2.0*volumeParams.worldSpaceGradientSpacing;
    gradients[1] /= 2.0*volumeParams.worldSpaceGradientSpacing;
    gradients[2] /= 2.0*volumeParams.worldSpaceGradientSpacing;
    gradients[3] /= 2.0*volumeParams.worldSpaceGradientSpacing;
    return gradients;
}

// Compute world space gradient using backward difference: f' = ( f(x)-f(x-h) ) / h
