Эффективно вычислить градиент для данных вокселей
Что является наиболее эффективным способом вычисления градиента для данных вокселей фиксированного размера, таких как исходный код ниже. Обратите внимание, что мне нужен градиент в любой точке пространства. Градиенты будут использоваться для оценки нормалей в реализации марширующих кубов.
#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
Решение из пользовательского примера:
Решение
Следующее создает линейно интерполированное градиентное поле:
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;
}
Другие решения
В общем, фильтры Собела дают немного более хорошие результаты, чем простая центральная тенденция, но требуют больше времени для вычислений (фильтр Собела по сути является гладким фильтром в сочетании с центральной тенденцией). Классический Sobel требует взвешивания 26 выборок, в то время как центральная тенденция требует только 6. Однако есть хитрость: с помощью графических процессоров вы можете получить аппаратную трилинейную интерполяцию бесплатно. Это означает, что вы можете вычислить Sobel с 8 чтениями текстур, и это можно сделать параллельно через GPU. Следующая страница иллюстрирует эту технику с использованием GLSL
http://www.mccauslandcenter.sc.edu/mricrogl/notes/gradients
Для вашего проекта вы, вероятно, захотите вычислить градиенты на GPU и использовать методы GPGPU, чтобы скопировать результаты обратно из GPU в CPU для дальнейшей обработки.
Одним из важных методов оптимизации во многих реализациях является компромисс между временем и пространством. В качестве рекомендации, возможно, стоит посмотреть в любом месте, где вы можете предварительно рассчитать и кэшировать свои результаты.
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;
}
Это может быть написано более кратко, но, возможно, таким образом вы все еще можете видеть, что происходит. Если это все еще не достаточно гладко, вам придется изучить кубическую / сплайн-интерполяцию или подобное.