Raytracing render на C

в 9:57, , рубрики: 3d графика, openmp, ray tracing, Анимация и 3D графика, Программирование, метки: , ,

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

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

Raytracing render на C

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

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

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

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

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

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

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

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

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

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

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

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

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

Описание структуры Object3d из исходников рендера

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;

Raytracing render на C

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

В качестве примера: код сферы из исходников рендера

sphere.c

// ... инклуды

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

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

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

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

Raytracing render на C

Изначально скорость рендеринга была удовлетворительной только для тех сцен, которые содержали меньше тысячи полигонов. Поскольку у меня 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, соответственно):
Raytracing render на C
Поскольку мы делим сцену на части параллельно осям координат — то можно относительно быстро просчитать какую часть сцены пересекает луч. Возможны следующие варианты:

  • Луч пересекает только часть сцены L
    Можно не просматривать объекты, содержащиеся в части R
    (на картинке — луч 1)

  • Луч пересекает только часть сцены R
    Можно не просматривать объекты, содержащиеся в части L
    (на картинке — луч 3)

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

  • Луч пересекает сначала часть сцены R, а потом часть сцены L
    Такая же логика поиска пересечений, как и в предыдущем варианте (только сначала рассматриваем часть сцены R)

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

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

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

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

Surface Area Heuristic

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

Raytracing render на C

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

Raytracing render на C

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

Наблюдения

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

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

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

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

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

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

Пример очень простой модели: tetrahedron.obj

# tetrahedron

# vertexes:
v 1.00 1.00 1.00
v 2.00 1.00 1.00
v 1.00 2.00 1.00
v 1.00 1.00 2.00

# triangles:
f 1 3 2
f 1 4 3
f 1 2 4
f 2 3 4

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

Структура 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 полигонов:
Raytracing render на C

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

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

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

Raytracing render на C

Заключение

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

Исходники рендера на гитхабе: 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).

Автор: stemm

Источник

Поделиться

* - обязательные к заполнению поля