Эффективно вычисляйте градиент для воксельных данных

Каков наиболее эффективный способ вычисления градиента для данных вокселей фиксированного размера, таких как исходный код ниже. Обратите внимание, что мне нужен градиент в любой точке пространства. Градиенты будут использоваться для оценки нормалей в реализации марширующих кубов.

#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. Решение из пользовательского примера:

введите здесь описание изображения


person Mortennobel    schedule 22.01.2014    source источник
comment
Мне интересно, зачем вам нужен градиент любой точки. Кажется, это был бы легкий шанс для оптимизации. Кроме того, регулярное расстояние между точками (в отличие от типичного затенения, когда размер треугольника варьируется) может оказаться лучшим подходом, чем get_gradient по точкам. Вы рассматривали такую ​​возможность?   -  person user3125280    schedule 24.01.2014
comment
Существует потенциальная оптимизация производительности, так как мне нужен только градиент в каждой вершине, генерируемой марширующими кубами, и эти вершины всегда имеют две оси, расположенные на целых числах (другими словами, это может упростить вычисления до двух линейных интерполяций и одной трилинейной). интерполяция).   -  person Mortennobel    schedule 27.01.2014
comment
если у вас есть градиент (который является вектором [x, y, z]) в каждой точке сетки, тогда трилинейная интерполяция (градиента) станет линейной интерполяцией вдоль одной оси (неинтегральной). С другой стороны, вы можете рассчитать градиент как производную интерполированного поля плотности. Я не знаю, какие плюсы/минусы у этого есть, но это будет медленнее. Какой путь вы намеревались?   -  person user3125280    schedule 27.01.2014
comment
Откуда вы берете свои данные? Раньше я получал удивительно гладкие результаты, используя исходный источник данных, с помощью которых были сгенерированы воксели, а не сами воксели. При условии, что данные представлены в форме, которую можно дифференцировать, например, в виде суммы сглаженных полей, это дает прекрасные результаты.   -  person Andy Newman    schedule 28.01.2014
comment
@AndyNewman: Воксельные данные, над которыми я работаю, являются результатом вычислений на основе FEM (метод конечных элементов) (точнее, многосеточной оптимизации топологии в 3D). Насколько мне известно, градиент результат должен быть рассчитан на основе вокселей с использованием схемы интерполяции. Меня больше всего беспокоит, как получить наилучший визуальный результат, используя разумный быстрый подход (на самом деле я не думаю, что производительность будет большой проблемой).   -  person Mortennobel    schedule 28.01.2014


Ответы (4)


Следующее создает линейно интерполированное поле градиента:

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;
}
person MooseBoys    schedule 24.01.2014
comment
Нет ли более точного способа оценки градиентов? Я попытался реализовать ваше предложение, но был весьма разочарован качеством найденных градиентов. - person Mortennobel; 28.01.2014
comment
@Mortennobel Каково ваше определение точности? Существует несколько способов интерполяции значений между дискретными точками данных, наиболее распространенным из которых является линейный. Другим является кубическая интерполяция, которая дает более гладкие результаты (и первая, и вторая производные непрерывны), но вызывает перерегулирование. - person MooseBoys; 28.01.2014
comment
Ну, я просто ищу метод, который дает визуально приятный результат (то есть более гладкий результат). - person Mortennobel; 28.01.2014
comment
@Mortennobel Вы просили эффективный способ, а не точный (или визуально более плавный, как вы это называете). Возможно, вы можете добавить это в качестве комментария к исходному вопросу? - person Sergio Ayestarán; 28.01.2014

Одним из важных методов оптимизации во многих реализациях является компромисс между временем и пространством. В качестве предложения, стоит посмотреть везде, где вы можете предварительно вычислить и кэшировать свои результаты.

person Brad Thomas    schedule 30.01.2014

В общем, фильтры Собеля дают немного более хорошие результаты, чем простая центральная тенденция, но требуют больше времени для вычисления (Собеля, по сути, является гладким фильтром в сочетании с центральной тенденцией). Классический Sobel требует взвешивания 26 отсчетов, а центральная тенденция требует только 6. Однако есть хитрость: с GPU вы можете бесплатно получить аппаратную трилинейную интерполяцию. Это означает, что вы можете вычислить Sobel с 8 чтениями текстур, и это можно сделать параллельно через GPU. На следующей странице показан этот метод с использованием GLSL http://www.mccauslandcenter.sc.edu/mricrogl/notes/gradients Для вашего проекта вы, вероятно, захотите вычислить градиенты на графическом процессоре и использовать методы GPGPU для копирования результатов обратно из графического процессора в ЦП для дальнейшей обработки.

person user1677899    schedule 28.08.2014

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

Вероятно, это можно написать немного более лаконично, но, возможно, так вы все же сможете увидеть, что происходит. Если это все еще недостаточно гладко, вам придется изучить кубическую / сплайновую интерполяцию или что-то подобное.

person example    schedule 28.01.2014
comment
Я реализовал ваше предложение, но результат выглядит довольно плохо. (См. обновление вопроса для скриншота). - person Mortennobel; 29.01.2014
comment
@Mortennobel уф. Она должна быть не хуже версии mooseboys. возможно, ошибка знака какого-то рода. позвольте мне просмотреть мой код, как только я вернусь с работы. - person example; 29.01.2014
comment
@Mortennobel В расчете yd0 не хватало -1. к сожалению, я не могу скомпилировать ваш код без особых усилий, но я не могу найти других ошибок в коде. надеюсь, это была ошибка, которая вызвала артефакты... - person example; 31.01.2014
comment
Это исправило артефакты. Однако я не нахожу визуальный результат лучше, чем решение, предложенное MooseBoys. - person Mortennobel; 31.01.2014