CC BY
41
ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ
Вестн. Ом. ун-та. 2010. №4. С. 164-167.
УДК 004.92
. . . ЖИДКОСТИ
Предложен параллельный алгоритм построения восьмеричного дерева по значениям скалярного поля и извлечения из него изоповерхности. Рассмотрены особенности реализации алгоритма на графическом процессоре на базе архитектуры CUDA. Приведено сравнение данного алгоритма с алгоритмом на основе регулярной сетки по объёму занимаемой памяти и времени построения изоповерхности. Результаты показали, что при большой детализации поверхности реализация на базе октодерева оказывается эффективнее и является оптимальной в случаях, когда поток жидкости занимает небольшую часть объемлющего пространства.
Ключевые слова: параллельный алгоритм, визуализация поверхности жидкости, архитектура CUD А.
Введение
Существует множество задач, требующих моделирования и визуализации потоков жидкостей в неоднородной среде. В данной статье рассматривается один из способов визуализации жидкости как скалярного поля. Визуализация разделена на два этапа: создание восьмеричного дерева (октодерева) и построение на его основе поверхности методом движущихся кубов [1]. Классическим способом использования восьмеричных деревьев на графическом процессоре (GPU) является их создание на центральном процессоре и копирование его представления на графический процессор. В нашем случае дерево создаётся полностью на GPU с помощью параллельного алгоритма. Для создания изоповерхности вершины узлов-ячеек октодерева помечаются как внешние либо как внутренние, создавая, таким образом, границу, по которой проходит изоповерхность.
Алгоритм построения поверхности основан на статье «Highly Parallel Surface Reconstruction» [2], в которой авторы представляют способ создания октодерева для массива ориентированных точек, т. е. точек и нормалей к ним, полученных с ЗО-сканера [3]. Затем данный алгоритм создаёт сглаженную изоповерхность восьмеричного дерева с помощью техники движущихся кубов [4]. В нашем случае применяется модифицированный алгоритм, предназначенный для создания поверхности потока частиц, не содержащего информации об ориентации точек или границы поверхности.
©А.Ю. Суравикин, В.В. Коробицын, 2010
Построение октодерева с помощью CUDA
Предлагаемый алгоритм обрабатывает массив позиций частиц и создаёт восьмеричное дерево в виде массива узлов, разделённого по уровням. Также для этого дерева создаются массивы вершин и рёбер. Перечислим шаги предлагаемого алгоритма и дадим некоторые комментарии.
Расчёт ограничивающего параллелепипеда. Для расчёта параллелепипеда, содержащего все частицы, используется функция transform_reduce библиотеки thrust [5]. Пользовательские операции реализуются в виде функторов, т. е. структур с переопределенным оператором ( ) [6]. Из массива частиц формируется массив пустых параллелепипедов, который редуцируется в один параллелепипед, содержащий все исходные параллелепипеды [7].
Создание ключей узлов дерева. Для
определения положения узла в восьмеричном дереве используется «перемешанный» ключ, состоящий из набора трёхбитных кодов. Код характеризует область, занимаемую каждым из 8 потомков дерева, в порядке xyz: если значение бита х равно
1, потомок занимает октант, находящийся справа по оси х, иначе - октант, находящийся слева по оси х. Таким образом, смешанный xyz-ключ узла дерева на уровне D определяется как битовая строка, которая показывает путь от корня дерева до листа. Позиция частицы трансформируется в значение ключа следующим образом: сначала позиция частицы Р с помощью минимальных координат ограничивающего куба Дшп и вектора масштабирования S переводится в пространство куба стороной 2D: S = 2° / (Бтах-Дшп), Pr= (Р - Вт т) 1 S. Умножение на вектор S производится покомпонентно. Координаты центра корневого узла изначально равны (2D_1; 2,h 1; 2,h 1). Далее центр узла вычисляется на основе значения ключа для каждого уровня иерархии дерева.
Сортировка ключей. Для реализации параллельного поиска количества частиц в каждом узле необходимо отсортировать ключи, соответствующие позициям частиц. Для сортировки используется функция sorting: :radix_sort() библиотеки thrust, использующая алгоритм поразрядной сортировки (radix sort) [8].
Поиск и подсчёт уникальных ключей. Расчёт числа частиц и удаление повторяющихся кодов узлов происходит в отсортированном массиве в три этапа:
1) расчёт количества уникальных кодов;
2) расчёт адресов уникальных кодов с помощью сканирования массива-счётчика;
3) заполнение массива уникальных кодов и соответствующего массива количеств повторяющихся элементов. Эти шаги заменяются суммированием одинаковых элементов атомарной операцией atomicAdd (элементарная операция архитектуры CUDA, выполняемая на графическом процессоре за один такт) и выделением уникальных элементов.
Расчёт адресов узлов октодерева. Адреса узлов рассчитываются аналогично поиску уникальных ключей и подсчёту их количества, однако в данном случае узлы объединяются в группы по 8 элементов таким образом, чтобы все узлы такой группы принадлежали одному предку. Ключи отсортированы по возрастанию и пронумерованы в соответствии со структурой дерева, поэтому адреса групп узлов рассчитываются как префиксная сумма массива, содержащего результаты сканирования уникальных кодов. При сканировании значения соседних ключей сравниваются не напрямую, а определяется, принадлежат ли они одному предку. Для этого используется логическое выражение ((left | right) & 0xfffffff8), где left - значение текущего кода, right - значение следующего кода. Таким образом, если соседние узлы имеют одинаковые биты, кроме трех последних, то они принадлежат одному предку.
Создание массива узлов на уровне листьев. Массив листьев дерева заполняется данными узлов по вычисленным на предыдущем шаге адресам групп, к которым прибавляется 3 младших бита ключа текущего узла. В узел записывается его ключ, количество содержащихся в нем частиц, а также массив потомков, который заполняется специальным значением, соответствующим несуществующим узлам.
Создание массивов узлов на высоких уровнях. После создания массива листьев мы формируем массив их предков в три этапа: 1) расчёт уникальных ключей для узлов-родителей; 2) заполнение массива адресов узлов-родителей аналогично подсчету адресов уникальных кодов;
166
3) создание массива узлов-родителей по рассчитанным ранее адресам. Таким образом, из массива узлов-листьев получили массив их предков. Тем же способом создаются и более высокие уровни иерархии вплоть до корневого узла.
Поиск соседей для каждого узла. Поиск соседних узлов в отличие от расчета узлов-предков начинается с корневого дерева и далее снижается до листьев. Каждому узлу в массив соседей записывается 27 значений: 26 соседей и сам узел. Соседей корневого узла не существует, поэтому в центральный узел помещается адрес корня, а остальные - специальным значением ОСТ RE E_I N VALI D_KEY.
Заранее рассчитанные таблицы локальных соседей LUTparent и локальных потомков LUTchild показывают отношение соседних узлов для родительских узлов и узлов-потомков. Программа-ядро поиска соседей обрабатывает все узлы и заполняет каждый из элементов массива neighbors каждого узла.
Создание вершин для узлов дерева. Массив вершин заполняется аналогично ранее рассмотренным элементам, т. е. сначала рассчитывается количество вершин в текущем узле, а затем адреса вершин в создаваемом массиве. Однако в данном случае есть своя особенность: одна вершина может принадлежать нескольким соседним узлам, и, следовательно, каждая вершина может быть продублирована несколько раз. Чтобы избежать излишнего выделения памяти, общая вершина хранится только в её узле-собст-веннике, который определяется как узел, имеющий минимальное значение xyz-кода. Чтобы получить множество вершин, которыми владеет рассматриваемый узел, мы используем побитовые операции с инвертированной маской вершины. В масках установленные биты соответствуют локальному номеру вершины.
Создание рёбер для узлов дерева. Массив рёбер создаётся аналогично массиву вершин. Для каждого из 13 предыдущих соседей создаётся 12-битная маска рёбер. Ребро представляет собой адреса двух смежных вершин, поэтому может оказаться так, что одной из смежных вершин нет в данном узле. В таком случае эта вершина ищется в соседних узлах.
____________А. Ю. Суравикин, В. В. Коробицын
Извлечение изоповерхности из октодерева
Для рендера изоповерхности скалярного поля необходимо создать вершинный буфер и буфер индексов. Создание каждого буфера разделяется на два этапа аналогично созданию массивов элементов октодерева: подсчёт адресов элементов, заполнение буферов. Рассмотрим последовательно эти этапы.
Подсчёт вершин и индексов изоповерхности. В соответствии с методом движущихся кубов полигоны создаются в узлах, пересекающих изоповерхность поля, причём вершины этих полигонов находятся на рёбрах кубов. Программа-ядро параллельно подсчитывает количество вершин и отмечает ребра, на которых будут располагаться вершины. Сначала необходимо определить, какие вершины текущего узла внутренние, а какие - внешние по отношению к поверхности. Мы полагаем вершину внутренней, если у нее есть все 7 соседних узла, иначе вершина считается внешней. Количество вершин и ребер в каждом узле сохраняется в соответствующем элементе массива, чтобы далее создать массив их адресов, а также для определения наличия вершины в узле.
Создание буфера вершин. Простейший расчёт позиции вершины проходит для каждого ребра, которое пересекает изоповерхность и заключается в поиске среднего арифметического между позициями вершин узла, находящихся на концах пересекаемого ребра. Линейная интерполяция исходных вершин ребра с коэффициентом, зависящим от значения скалярного поля в этой точке, более сложная в реализации, но создаёт более близкую аппроксимацию изоповерхности.
Создание буфера индексов. Массив индексов заполняется в ядре, которое выполняется для каждого узла октодерева. Для получения локальных индексов рёбер, на которых находятся вершины, формирующие треугольники поверхности, используется таблица, хранящая множество ребер для каждой комбинации внутренних и внешних вершин. Индексы треугольников для узла дерева вычисляются с помощью адреса текущего узла и массивов-адресов.
Используя CUDA API и графическое API, можно скопировать эти буферы в вершинный и индексный буферы OpenGL
или DirectX для рендера, либо заранее установить указатели на буферы графических API как выходные данные для программ-ядер.
Результаты
Согласно предложенному алгоритму была разработана программа на основе архитектуры CUDA. Для вычислительного эксперимента мы использовали графический процессор GeForce GTX 260 с объемом графической памяти 896 МБ. Результаты эксперимента приведены в табл. 1 и
2. Видно, что октодерево становится более эффективным, чем равномерная сетка, когда уровень дерева (логарифм размерности поля по основанию 2) больше 6. Значение, отмеченное знаком * показывает теоретический размер памяти, так как мы не смогли выделить такой объём памяти на нашей экспериментальной базе.
Таблица 1 Сравнение объёмов используемой памяти алгоритмов извлечения изоповерхности
Размеры поля (уровень дерева) Объём памяти
для октодерева, МБ для равномерной сетки, МБ
32 (5) 16,32 2,04
64 (6) 32,03 16,14
128 (7) 48,88 128,57
256 (8) 91,94 1026,26*
Таблица 2 Сравнение времени построения изоповерхности
Размеры поля (уровень дерева) Время расчета
для октодерева, мс для равномерной сетки, мс
32 (5) 26,3 6,4
64 (6) 37,7 12,3
128 (7) 43,4 58,5
256 (8) 48,6 *
В статье представлена реализация алгоритма, который выполняется полностью на СРи (подробнее см. [9]). Он подходит для таких моделей потоков жидкости, в
которых частицы заполняют менее 10 % объёма пространства. Создание октодере-ва в реальном времени с быстрым доступом к данным узла и его соседям может оказаться полезным не только для создания изоповерхности, но и для других целей, например, для определения и отсечения невидимых объектов. Алгоритм извлечения изоповерхности создаёт аппроксимацию поверхности с производительностью, достаточной для интерактивного рендера, но его необходимо доработать, добавив сглаживание скалярного поля и расчёт интерполированных вершин и нормалей.
ЛИТЕРАТУРА
[1] Wilhelms J., Gelder A. V. Octrees for faster isosurface generation // ACM Trans. Graph. 1992. Vol. 11. №. 3. P. 201-227.
[2] Zhou K., Gong М., Huang X., Guo B. Highly Parallel Surface Reconstruction. URL: http://www. kunzhou.net/2008/MSR-TR-2008-53.pdf (дата обращения: 01.04.2010).
[3] Bajaj С. L, Bernardini F., Xu G. Automatic reconstruction of surfaces and scalar fields from 3d scans // Proceedings of SIGGRAPH'95. 1995. P. 109-118.
[4] Lorensen W, Cline H. Marching cubes: a high resolution 3D surface construction algorithm // Computer Graphics. 1987. Vol. 21. №4. P. 163169.
[5] Hoberock J., Bell N. Thrust: C++ Template Library for CUDA. URL: http://code.google.com/ р/thrust/ (дата обращения: 14.05.2010).
[6] Hickey R. Callbacks in C++ Using Template Functors. 1994. URL: http://www.tutok.sk/fastgl/ callback.html (дата обращения: 14.05.2010).
[7] Hoberock J., Bell N. Thrust library usage examples 2010. URL: http://thrust.googlecode.com/ files/examples-v1.2.zip (дата обращения: 14.05.2010).
[8] Hoberock J., Bell N. Sorting Performance Optimizations. 2009. URL: http://www.meganewtons. com/2009/08/sorting-performance-optimizations. html (дата обращения: 14.05.2010).
[9] Суравикин А. Ю., Коробицын В. В. Визуализа-
ция скалярного поля на основе динамического построения изоповерхности // Математические структуры и моделирование. 2009. № 20.
С. 121-133.
