Raytracing render на C

В данной статье хочу поделиться собственным подходом и полученными результатами.

Несколько слов про рейтрейсинг

Представим что у нас есть камера (для упрощения будем считать, что наша камера это материальная точка). Также у нас есть плоскость проектирования, которая является набором пикселей. Теперь, от камеры, к каждому пикселю плоскости проектирования будем проводить векторы (первичные лучи) и для каждого луча находить ближайший объект сцены, который он пересекает. Цветом, что соответствует точке пересечения, можно закрасить соответствующий пиксель на плоскости проектирования. Повторив эту процедуру для всех точек плоскости проектирования — мы получим изображение трехмерной сцены. Если ограничиться только этими операциями, то полученный алгоритм будет называться Ray casting.

Рекурсивный Ray casting == Ray tracing

• Из точки пересечения первичного луча с объектом сцены можно выпустить вторичный луч — направленный к источнику света. Если он пересекает какой-либо объект сцены, значит данная точка объекта находится в тени. Этот прием позволяет получать геометрически правильные тени.
• Если объект обладает зеркальной поверхностью, то по законам геометрической оптики можно рассчитать отраженный луч — и запустить для него процедуру рейтрейсинга (рекурсивно). Затем смешать полученный от зеркального отражения цвет с собственным цветом поверхности. Этот прием позволяет моделировать зеркальные поверхности (подобным образом можно моделировать прозрачные поверхности).
• Можно рассчитать расстояние от камеры до точки пересечения луча с объектом сцены. Значение расстояния можно использовать для моделирования тумана: чем дальше объект находится от камеры тем меньше интенсивность его цвета (для расчетов можно использовать, например, экспоненциальную функцию убывания интенсивности).

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

• Цвет каждого пикселя можно рассчитать независимо от других (т.к. лучи, выпускаемые из камеры никак не влияют друг на друга — их можно обрабатывать параллельно)
• Алгоритм не содержит ограничений на геометрию объектов (можно использовать богатый набор примитивов: плоскости, сферы, цилиндры и т.п.)

Архитектура рендера

Все объекты хранятся в куче. Рендер оперирует массивом указателей на 3D объекты (на самом деле все объекты еще упакованы в kd-дерево, но пока будем считать, что дерево отсутствует). На данный момент есть 2 типа примитивов: трианглы и сферы.

Как научить рендер оперировать различными примитивами?

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

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

Язык программирования С не имеет синтаксических конструкций для программирования в парадигме ООП, но в данном случае на помощь приходят структуры и указатели на функции.

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

typedef  struct { 	// указатель на область данных, содержащую параметры конкретного примитива 	// в случае полигона: координаты вершин, цвет (или битмапка с текстурой), свойства материала 	// в случае сферы: координаты центра, радиус, и т.п. 	void * data;      	// указатель на функцию, которая умеет обрабатывать примитив, 	// на который ссылается указатель data 	Boolean (*intersect)(const void * data,                          const Point3d ray_start_point,                          const Vector3d ray,                          // сюда будет сохранятся точка пересечения луча и примитива                          Point3d * const intersection_point);  	// получение цвета в точке пересечения 	// здесь можно возвращать просто цвет объекта 	// или обеспечить процедурную генерацию текстуры 	// или использовать загруженную из файла текстуру :) 	Color (*get_color)(const void * data,                        const Point3d intersection_point);      // получение вектора нормали в точке пересечения (используется для рассчёта освещённости)     // в случае сферы и полигона - вектор нормали можно рассчитать аналитически     // как вариант, вместо аналитечских рассчётов - объект может содержать карту нормалей     Vector3d (*get_normal_vector)(const void * data,                                   const Point3d intersection_point);  	// вызывается рендером перед удалением Object3d 	// тут можно освобождать ресурсы, связанные с объектом 	// например, удалять текстуры 	void (*release_data)(void * data); } Object3d; 

Такой подход позволяет локализировать код относящийся к каждому 3D примитиву в отдельном файле. Следовательно, довольно легко вносить изменения в структуры 3D объектов (например, добавить поддержку текстур или карт нормалей), или даже добавлять новые 3D примитивы. При этом — не нужно вносить изменения в код рендера. Получилось что-то на подобии инкапсуляции.

// ... инклуды  typedef struct {     Point3d center;     Float radius;     Color color; } Sphere;  // ... декларации функций  // "конструктор" сферы Object3d * new_sphere(const Point3d center,            const Float radius,            const Color color) {          Sphere * sphere = malloc(sizeof(Sphere));     sphere->center = center;     sphere->radius = radius;     sphere->color = color;      // упаковываем сферу в обобщённую структуру 3D объекта     Object3d * obj = malloc(sizeof(Object3d));     obj->data = sphere;      // добавляем функции, которые умеют работать со структурой Sphere 	obj->get_color = get_sphere_color;     obj->get_normal_vector = get_sphere_normal_vector; 	obj->intersect = intersect_sphere;     obj->release_data = release_sphere_data;          return obj; }  // цвет сферы всюду одинаковый // но здесь можно реализовать процедурную генерацию текстуры static Color get_sphere_color(const void * data,                  const Point3d intersection_point) { 	const Sphere * sphere = data; 	return sphere->color; }  // вычисляем аналитически вектор нормали в точке на поверхности сферы // сюда можно впилить Bump Mapping static Vector3d get_sphere_normal_vector(const void * data,                          const Point3d intersection_point) {     const Sphere * sphere = data;         // vector3dp - служебная функция, которая создаёт вектор по двум точкам     Vector3d n = vector3dp(sphere->center, intersection_point);     normalize_vector(&n);     return n; }  // пересечение луча и сферы static Boolean intersect_sphere(const void * data,                  const Point3d ray_start_point,                  const Vector3d ray,                  Point3d * const intersection_point) {     const Sphere * sphere = data;         // чтобы найти точку пересечения луча и сферы - нужно решить квадратное уравнение     // полную реализацию метода можно найти в исходниках на GitHub }  // "деструктор" сферы void release_sphere_data(void * data) { 	Sphere * sphere = data; 	// если бы сфера содержала текстуры - их нужно было бы здесь освободить 	free(sphere); } 

void example(void) { 	Object3d * objects[2];  	// красная сфера, с центром в точке (10, 20, 30) и радиусом 100 	objects[0] = new_sphere(point3d(10, 20, 30), 100, rgb(255, 0, 0));  	// зелёный треугольник с вершинами (0, 0, 0), (100, 0, 0) и (0, 100, 0) 	objects[1] = new_triangle(point3d(0, 0, 0), point3d(100, 0, 0), point3d(0, 100, 0), rgb(0, 255, 0));  	Point3d camera = point3d(0, 0, -100);  	Vector3d ray = vector3df(0, -1, 0);  	int i; 	for(i = 0; i < 2; i++) { 		Object3d * obj = objects[i];  		Point3d intersection_point;  		if(obj->intersect(obj->data, camera, ray, &intersection_point)) { 			// вот так можно получить цвет в точке пересечения луча и примитива 			Color c = obj->get_color(obj->data, intersection_point); 			// делаем что-нибудь ещё :-) 			// например, можно искать ближайшую точку пересечения 			// и нам совсем не важно, что именно лежит в массиве objects! 		} 	} } 

Ради справедливости, стоит заметить, что такая архитектура влечёт за собой много жонглирования указателями.

Скорость рендеринга

В качестве тестового примера — я решил визуализировать сцену с фракталом «Пирамида Серпинского», зеркальной сферой и подставкой с текстурой. Пирамиду Серпинского довольно удобно использовать для тестирования скорости рендеринга — т.к. задавая различные уровни рекурсии можно генерировать различное количество трианглов:

Изначально скорость рендеринга была удовлетворительной только для тех сцен, которые содержали меньше тысячи полигонов. Поскольку у меня 4х ядерный процессор — было принято решение ускорить рендеринг путём распараллеливания.

POSIX Threads

Первое впечатление: семантика Java Threads очень близка к pthreads. Поэтому, особых проблем с пониманием модели POSIX потоков не было. Было приянто решение запилить свой Thread Pool 🙂

Thread Pool содержал очередь для тасков. Каждый таск являлся структурой, содержащей ссылку на функцию, которую нужно выполнить, а также ссылку на список аргументов. Доступ к очереди тасков регулировался посредством мютекса и condition-переменной. Для удобства, весь рендеринг инкапсулирован в отдельной функции, одним из аргументов которой — является количество потоков:

void render_scene(const Scene * const scene,              const Camera * const camera,              Canvas * canvas,              const int num_threads);  // структура Scene содержит 3D объекты, источники света, цвет фона и параметры тумана // структура Camera содержит координаты, углы наклона и фокус камеры // структура Canvas эмулирует 2х мерный холст - именно в него и ренедрится изображение 

Эта функция содержала код, связующий цикл рендеринга и пул потоков, в обязанности которого входило:

• разбивать холст на несколько частей
• для каждой части холста формировать отдельный таск на рендеринг
• отправлять все таски в пул и ожидать завершения

Но, к сожалению, на 2х ядерном ноуте рендер периодически зависал или падал с ошибкой Abort trap 6 (причём на 4х ядерном это ни разу не воспроизвелось). Я пытался пофиксить это печальный баг, но вскоре нашёл более привлекательное решение.

OpenMP

OpenMP берёт на себя всю заботу по созданию и распределению тасков, а также организовывает барьер для ожидания завершения рендеринга. Достаточно дописать всего несколько директив чтобы распараллелить однопоточный код, а баги с зависанием исчезли 🙂

void render_scene(const Scene * const scene,              const Camera * const camera,              Canvas * canvas,              const int num_threads) {          const int width = canvas->width;     const int height = canvas->height;     const Float focus = camera->focus;          // устанавливаем количество потоков     omp_set_num_threads(num_threads);          int x;     int y;     // всего две строчки, для того чтобы распараллелить рендеринг :-)     #pragma omp parallel private(x, y)     #pragma omp for collapse(2) schedule(dynamic, CHUNK)     for(x = 0; x < width; x++) {         for(y = 0; y < height; y++) {             const Vector3d ray = vector3df(x, y, focus);             const Color col = trace(scene, camera, ray);             set_pixel(i, j, col, canvas);         }     } } 

Рендеринг немного ускорился, но, увы, сцены содержащие больше тысячи примитивов — по прежнему рендерились очень медленно.

Просчёт пересечения луча с примитивами — относительно ресурсоёмкая процедура. Проблема заключается в том, что для каждого луча проверяются пересечения со всеми примитивами сцены (ищется самый близкий примитив, пересекаемый лучом). Можно ли каким-нибудь образом исключить те объекты, с какими луч точно не пересекается?

Kd-Tree

Давайте разобьём сцену с объектами на две части: «левая» и «правая» (обозначим их как L и R, соответственно):

Поскольку мы делим сцену на части параллельно осям координат — то можно относительно быстро просчитать какую часть сцены пересекает луч. Возможны следующие варианты:

• Луч пересекает только часть сцены L
  Можно не просматривать объекты, содержащиеся в части R
  (на картинке — луч 1)
• Луч пересекает только часть сцены R
  Можно не просматривать объекты, содержащиеся в части L
  (на картинке — луч 3)
• Луч пересекает сначала часть сцены L, а потом часть сцены R
  Сначала просматриваем объекты, принадлежащие части сцены L — если пересечение было найдено, тогда объекты принадлежащие части сцены R можно не просматривать. Если луч не пересекает ни один объект из части L — придётся просмотреть объекты из части R
  (на картинке — луч 2)
• Луч пересекает сначала часть сцены R, а потом часть сцены L
  Такая же логика поиска пересечений, как и в предыдущем варианте (только сначала рассматриваем часть сцены R)

Но, точно таким же образом, чтобы сократить количество просматриваемых полигонов — можно продолжить рекурсивно дробить сцену на всё меньшие участки. Полученная иерархическая структура, содержащая сегменты сцены, с привязанными к ним примитивами — называется kd-дерево.

Давайте, для примера, рассмотрим процесс трассировки луча 2:

1. Луч пересекает сначала левый сегмент сцены (L), потом правый — R
2. Из части L — луч пересекает только сегмент LL
3. Сегмент LL не содержит никаких объектов
4. Переходим к правой половине сцены — R
5. В правой половине сцены луч пересекает сначала сегмент RL, а потом RR
6. В части сцены RL было найдено пересечение луча с объектом
7. Трассировка завершена

Благодаря организации древовидной структуры данных, мы исключили те объекты, которые заведомо не пересекаются лучом. Но есть ещё очень важный нюанс — эффективность kd-дерева зависит от того, каким образом мы разбиваем сцену на части. Как правильно выбирать места разбиения сцены?

Surface Area Heuristic

Вероятность попадания луча в какой-нибудь сегмент сцены — пропорциональна площади этого сегмента. Чем больше примитивов содержится в участке сцены — тем больше пересечений нужно будет просчитывать при попадании луча в этот сегмент. Таким образом, можно сформулировать критерий разбиения: нужно минимизировать произведение количества примитивов и площади сегмента, в котором они содержатся. Данный критерий называется Surface Area Heuristic (SAH):

Давайте, на простом примере, рассмотрим эвристику в действии:

Таким образом, использование SAH позволяет эффективно разделять пустое пространство и 3D объекты — что очень положительно сказывается на производительности рендера.

Наблюдения

В рендере реализована возможность подсчёта среднего количества пересечений на 1 пиксель изображения: каждый раз, когда рассчитывается пересечение луча с каким-нибудь объектом — увеличивается значение счётчика, по окончании рендеринга — значение счётчика делится на количество пикселей в изображении.

Полученные результаты варьируются для различных сцен (т.к. построение kd-дерева зависит от взаимного расположения примитивов). На графиках изображена зависимость среднего количества пересечений на 1 пиксель изображения от количества примитивов:

Можно заметить что эта величина значительно меньше количества примитивов, содержащихся в сцене (если бы не было kd-дерева — то получилась бы линейная зависимость, т.к. для каждого луча нужно было бы искать пересечение со всеми N примитивами). Фактически, используя kd-дерево мы понизили сложность алгоритма рейтрейсинга от O(N) до O(log(N)).

Ради справедливости, стоит заметить что одним из недостатков данного решения является ресурсоёмкость построения kd-дерева. Но для статических сцен — это не критично: загрузили сцену, подождали пока построится дерево — и можно отправляться путешествовать с камерой по сцене 🙂

Сцена, содержащая 16386 треугольников:

Загрузка моделей

Научившись рендерить большое количество примитивов — появилось желание загружать 3D модели. Был выбран довольно простой и распостранённый формат — OBJ: модель хранится в текстовом файле, который содержит список вершин и список полигонов (каждый полигон задаётся точками из списка вершин).

В процессе написания парсера OBJ формата, было замечено, что у многих моделей содержится ещё список нормалей к каждой вершине полигона. Векторы нормалей вершин можно интерполировать для получения вектора нормали в любой точке полигона — этот приём позволяет симулировать гладкие поверхности (см. затенение по Фонгу). Данную фичу было довольно легко реализовать в рамках текущей архитектуры рендера (нужно было просто добавить векторы нормалей в структуру Triangle3d, и написать функцию, которая интерполирует их для любой точки полигона).

typedef struct {      // координаты вершин 	Point3d p1; 	Point3d p2; 	Point3d p3;      // вектор нормали, рассчитанный по трём вершинам     // если используется затенение по Фонгу - тогда нам этот атрибут не важен     Vector3d norm;          Color color;      Material material;          // если текстура отсутствует - закрашиваем триангл цветом, который указан в color     Canvas * texture;      // текстурные координаты     // используем только тогда, когда есть текстура     Point2d t1;     Point2d t2;     Point2d t3;          // векторы нормали в вершинах     // используем только в случае затенения по Фонгу     Vector3d n1;     Vector3d n2;     Vector3d n3;      // есть ещё несколько служебных атрибутов     // которые используются для сокращения вычислений } Triangle3d; 

Пример отрендериной сцены, содержащей порядка 19000 полигонов:

Пример отрендериной сцены, содержащей примерно 22000 полигонов:

Ради интереса решил добавить скайбокс и загрузить модели автомобилей:

Сцена содержит порядка 100000 полигонов (kd-дерево строилось несколько минут)

Заключение

Я рад что выбрал данную задачу во время изучения С — т.к., она позволила узнать различные аспекты экосистемы языка, и получить красивые результаты.

Исходники рендера на гитхабе: github.com/lagodiuk/raytracing-render (в описании проекта — рассказано как запустить демку).

В процессе изучения С мне очень помогли 2 книги:

1. Брайан Керниган, Деннис Ритчи — Язык программирования Си — вначале испытывал некоторые затруднения при чтении данной книги. Но после прочтения Head First C — снова вернулся к этой книге. В данной книге есть много примеров и задачек, которые помогли в изучении
2. David Griffiths, Dawn Griffiths — Head First C — эта книга очень понравилась тем, что в ней даётся общее видение экосистемы языка С (как устроена память, как это работает на уровне ОС, в общих чертах описан make, valgrind, POSIX потоки, юнит тестирование, и есть даже несколько страниц про Arduino)

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

Ещё, хочу порекомендовать очень полезный ресурс: ray-tracing.ru/ — здесь, в доступной форме, описаны различные алгоритмы, применяемые для фотореалистичного рендеринга (в частности — kd-деревья и SAH).