Pull to refresh

Моделирование большого количества взаимодействующих друг с другом частиц

Reading time 6 min
Views 30K
Рассмотрим ситуацию, когда необходимо обрабатывать столкновения между объектами. Как вы в этом случае поступите? Вероятно, самым простым решением будет проверить каждый объект с каждым другим объектом. И это правильное решение, и все будет замечательно до тех пор пока объектов не много. Как только их станет порядка нескольких тысяч, вы заметите, что все стало как-то медленно работать. А если частиц несколько десятков тысяч или сотен? Тогда все замрет. Вот здесь уже интересно, на какие хитрости и оптимизации вы пойдете, чтобы решить такую проблему.

Для простоты, будем рассматривать 2D случай, частицы круглые, радиус частиц у всех одинаковый.

Содержание


1. Обзор алгоритмов
1.1. Полный перебор
1.2. Sweep & Prune
1.3. Регулярная сеть
2. Некоторые оптимизации
2.1. Sweep & Prune
2.2. Регулярная сеть
3. Сравнение скорости выполнения
4. Приложение (программа и исходный код)
5. Заключение



1. Обзор алгоритмов


1.1. Полный перебор

Это самый простой и в то же время самый медленный из всех возможных способов. Идет проверка между всеми возможными парами частиц.

    void Update()
    {
        for (int i = 0; i < NUM_P; ++i) {
            for (int j = i + 1; j < NUM_P; ++j) {
                Collision(&_particles[i], &_particles[j]);
            }
        }
    }

Сложность: O(n ^ 2)

Плюсы:
* Очень прост в понимании и реализации
* Не чувствителен к разным размерам частиц
Минусы:
* Самый медленный из всех

1.2. Sweep & Prune

Суть данного алгоритма заключается в том, что сперва мы сортируем все частицы по их левой границе вдоль оси OX, а затем каждую частицу проверяем только с теми у которых левая граница меньше правой границы текущей частицы.

Попробую описать работу алгоритма в картинках.
Начальное состояние:

Как видим, частицы не упорядочены.
После сортировки по левой границе частицы(позиция по X минус радиус частицы), получаем следующую картину:

Теперь можно приступать к перебору. Проверяем каждую частицу с каждой другой частицей до тех пор пока правая граница первой частицы больше левой границы второй частицы.

Рассмотрим на примере.
0 частица проверит столкновение с 1 и 2:

1 со 2 и 3.
2 только с 3.
3 проверит столкновение с 4, 5 и 6:

Думаю, что суть ясна.

Моя реализация данного алгоритма:

    void Update()
    {
        qsort(_sap_particles, NUM_P, sizeof(Particle *), CompareParticles);

        for (int i = 0; i < NUM_P; ++i) {
            for (int j = i + 1; j < NUM_P && _sap_particles[i]->pos.x + RADIUS_P > _sap_particles[j]->pos.x - RADIUS_P; ++j) {
                Collision(_sap_particles[i], _sap_particles[j]);
            }
        }
    }

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

Так же этот алгоритм можно применить для более сложных объектов, в этом случае следует использовать их AABB(Axis-Aligned Bounding Box, минимальный ограничивающий прямоугольник).

Сложность: O(n ^ 2) — в худшем случае, O(n * log(n)) — в среднем случае

Плюсы:
* Прост в реализации
* Не чувствителен к разным размерам частиц
* Хорошая производительность
Минусы:
* Требуется немного времени, чтобы понять алгоритм

1.3. Регулярная сеть

Как вы, вероятно, уже догадались из названия, суть алгоритма сводится к тому, что все пространство делится на равномерную сеть из маленьких квадратиков, размер которых равен диаметру частицы. Каждый такой квадратик(ячейка) этой сети представляет собой массив.
Выглядеть это может, например, так:
    const int _GRID_WIDTH = (int)(WIDTH / (RADIUS_P * 2.0f));
    const int _GRID_HEIGHT = (int)(HEIGHT / (RADIUS_P * 2.0f));

    std::vector<Particle *> _grid[_GRID_WIDTH][_GRID_HEIGHT];

На каждой итерации главного цикла эта сеть очищается и заново заполняется. Алгоритм заполнения предельно прост: индекс ячейки частицы вычисляется путем деления обеих координат на диаметр частицы и отбрасыванием дробной части. Пример:
    int x = (int)(_particles[i].pos.x / (RADIUS_P * 2.0f));
    int y = (int)(_particles[i].pos.y / (RADIUS_P * 2.0f));

Таким образом, для каждой частицы необходимо вычислить индекс и внести ее в ячейку с этим индексом(добавить элемент в массив).
Теперь остается только пройтись по каждой ячейке и проверить все ее частицы со всеми частицами соседних, окружающих ее 8 ячеек, при этом не забыв проверить с самой собой.

Это можно сделать так:
    void Update()
    {
        for (int i = 0; i < _GRID_WIDTH; ++i) {
            for (int j = 0; j < _GRID_HEIGHT; ++j) {
                _grid[i][j].clear();
            }
        }

        for (int i = 0; i < MAX_P; ++i) {
            int x = (int)(_particles[i].pos.x / (RADIUS_P * 2.0f));
            int y = (int)(_particles[i].pos.y / (RADIUS_P * 2.0f));

            _grid[x][y].push_back(&_particles[i]);
        }
    
        // далее много циклов проверки с соседними ячейками
    }

Сложность: O(n)

Плюсы:
* Самый быстрый из всех
* Прост в реализации
Минусы:
* Требуется дополнительная память
* Чувствителен к разным размерам частиц(требуется модификация)

2. Некоторые оптимизации


2.1. Sweep & Prune

Данный алгоритм можно немного улучшить, сделав выбор по какой оси сортировать частицы. Оптимальная ось та, вдоль которой находится наибольшее количество частиц.
Оптимальная ось OY:

Оптимальная ось OX:



2.2. Регулярная сеть

Оптимизация №1:
Используем более эффективное представление сети:
    const int _MAX_CELL_SIZE = 4;
    const int _GRID_WIDTH = (int)(WIDTH / (RADIUS_P * 2.0f));
    const int _GRID_HEIGHT = (int)(HEIGHT / (RADIUS_P * 2.0f));

    int _cell_size[_GRID_WIDTH * _GRID_HEIGHT];
    Particle *_grid[_GRID_WIDTH * _GRID_HEIGHT * _MAX_CELL_SIZE];

Думаю, не стоит говорить о том, что одномерный массив быстрее двумерного, а так же, что статический массив константной длины быстрее динамического std::vector, не говоря уже о том, что все это мы объединили. Так же пришлось завести еще один массив, который будет говорить о том сколько частиц находится в каждой из ячеек.
И да, мы ввели ограничение на количество частиц в ячейке, и это не совсем хорошо, т.к. грозит тем, что при сильном сдавливании, в ячейку потребуется записать значительно большее количество частиц, а т.к. это невозможно, некоторые столкновения не обработаются. Однако, при хорошем разрешении столкновений данной проблемы можно избежать.

Оптимизация №2
Вместо того, чтобы на каждой итерации главного цикла полностью очищать массив сетки, мы будем лишь вносить в нее те изменения, которые произошли.
Пробегаем по всем частицам каждой ячейки сети, вычисляем новый индекс частицы, и если он не равен текущему, то удаляем частицу из текущей ячейки и добавляем на свободное место в новую ячейку.
Удаление происходит путем записи на ее место последней частицы из этой же ячейки.
Вот мой код, который это реализует:
    for (int i = 0; i < _GRID_WIDTH * _GRID_HEIGHT; ++i) {
        for (int j = 0; j < _cell_size[i]; ++j) {
            int x = (int)(_grid[i * _MAX_CELL_SIZE + j]->pos.x / (RADIUS_P * 2.0));
            int y = (int)(_grid[i * _MAX_CELL_SIZE + j]->pos.y / (RADIUS_P * 2.0));
            
            if (x < 0) {
                x = 0;
            }
            if (y < 0) {
                y = 0;
            }
            if (x >= _GRID_WIDTH) {
                x = _GRID_WIDTH - 1;
            }
            if (y >= _GRID_HEIGHT) {
                y = _GRID_HEIGHT - 1;
            }

            // если индекс частицы изменился и в новой ячейке есть свободные места
            if (x * _GRID_HEIGHT + y != i && _cell_size[x * _GRID_HEIGHT + y] < _MAX_CELL_SIZE) {
                _grid[(x * _GRID_HEIGHT + y) * _MAX_CELL_SIZE + _cell_size[x * _GRID_HEIGHT + y]++] = _grid[i * _MAX_CELL_SIZE + j];
                // на место старой частицы записываем последнюю из этой же ячейки
                _grid[i * _MAX_CELL_SIZE + j] = _grid[i * _MAX_CELL_SIZE + --_cell_size[i]];
            }
        }
    }

Оптимизация №3
Мы используем проверку с соседними 8 ячейками. Это избыточно, достаточно проверять лишь 4 соседних ячейки.
Например, так:

С остальными ячейками проверка пройдет либо раньше, либо позже:


Оптимизация №4
Можно немного ускорить работу алгоритма увеличив локальность данных в памяти. Это позволит чаще читать данные из кеша процессора, и меньше обращаться к оперативной памяти. Но не стоит это делать на каждой итерации главного цикла, т.к. это убьет всю производительность, однако, можно проводить данную операцию, например, раз в секунду, или раз в несколько секунд, если сцена изменяется незначительно.
Увеличить локальность можно изменив порядок частиц в главном массиве, в соответствии с их расположением в сети.
Т.е. первые 4 частицы будут из 0 ячейки, следующие 4 из 1 и так далее.

3. Сравнение скорости выполнения


Количество частиц Bruteforce (мс) Sweep & Prune (мс) Регулярная сеть (мс)
1000 4 1 1
2000 14 1 1
5000 89 4 2
10000 355 10 4
20000 1438 28 7
30000 3200 55 11
100000 12737 140 23

Как видно из таблицы время работы алгоритмов соответствует их сложности, регулярная сеть растет почти линейно :).

4. Приложение (программа и исходный код)


http://rghost.ru/35826275

В архиве находится проект написанный в CodeBlocks, а так же бинарники под Linux, но т.к. для создания окна и вывода графики я использовал библиотеку ClanLib, то проект должен без проблем скомпилироваться и под Windows.

Управление в программе:
цифра 1 — Bruteforce
цифра 2 — Sweep & Prune
цифра 3 — Регуляная сеть
левая кнопка мышки — перемещение «расталкивалки».

5. Заключение


В данной статье приведен краткий обзор наиболее распространенных алгоритмов, использующихся для обработки столкновений. Однако, это не ограничивает область применения этих алгоритмов, ничто не мешает использовать их для любого другого взаимодействия. Основная задача, которую они на себя берут — это более удобное представление данных, для возможности быстрого отсеивания заранее неинтересных нам пар.
Tags:
Hubs:
+143
Comments 45
Comments Comments 45

Articles