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

Как всегда, в нашем распоряжении только голый компилятор, никаких сторонних библитек использовать нельзя. Я не хочу заморачиваться с оконными менеджерами, обработкой мыши/клавиатуры и тому подобным. Результатом работы нашей программы будет простая картинка, сохранённая на диск. Я совершенно не гонюсь за скоростью/оптимизацией, моя цель — показать основные принципы.

Итого, как в таких условиях нарисовать вот такую картинку за 180 строчек кода?



Давайте я даже анимированную гифку вставлю (шесть метров):



А теперь разобьём всю задачу на несколько этапов:

Этап первый: прочитать предыдущую статью

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

Этап второй: рисуем одну сферу

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



Код смотреть здесь, но давайте я приведу основной прямо в тексте статьи:

#define _USE_MATH_DEFINES
#include <cmath>
#include <algorithm>
#include <limits>
#include <iostream>
#include <fstream>
#include <vector>
#include "geometry.h"

const float sphere_radius   = 1.5;

float signed_distance(const Vec3f &p) {
    return p.norm() - sphere_radius;
}

bool sphere_trace(const Vec3f &orig, const Vec3f &dir, Vec3f &pos) {
    pos = orig;
    for (size_t i=0; i<128; i++) {
        float d = signed_distance(pos);
        if (d < 0) return true;
        pos = pos + dir*std::max(d*0.1f, .01f);
    }
    return false;
}

int main() {
    const int   width    = 640;
    const int   height   = 480;
    const float fov      = M_PI/3.;
    std::vector<Vec3f> framebuffer(width*height);

#pragma omp parallel for
    for (size_t j = 0; j<height; j++) { // actual rendering loop
        for (size_t i = 0; i<width; i++) {
            float dir_x =  (i + 0.5) -  width/2.;
            float dir_y = -(j + 0.5) + height/2.;    // this flips the image at the same time
            float dir_z = -height/(2.*tan(fov/2.));
            Vec3f hit;
            if (sphere_trace(Vec3f(0, 0, 3), Vec3f(dir_x, dir_y, dir_z).normalize(), hit)) { // the camera is placed to (0,0,3) and it looks along the -z axis
                framebuffer[i+j*width] = Vec3f(1, 1, 1);
            } else {
                framebuffer[i+j*width] = Vec3f(0.2, 0.7, 0.8); // background color
            }
        }
    }

    std::ofstream ofs("./out.ppm", std::ios::binary); // save the framebuffer to file
    ofs << "P6\n" << width << " " << height << "\n255\n";
    for (size_t i = 0; i < height*width; ++i) {
        for (size_t j = 0; j<3; j++) {
            ofs << (char)(std::max(0, std::min(255, static_cast<int>(255*framebuffer[i][j]))));
        }
    }
    ofs.close();

    return 0;
}

Класс векторов живёт в файле geometry.h, описывать я его здесь не буду: во-первых, там всё тривиально, простое манипулирование двух и трёхмерными векторами (сложение, вычитание, присваивание, умножение на скаляр, скалярное произвдение), а во-вторых, gbg его уже подробно описал в рамках курса лекций по компьютерной графике.

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

Итак, в функции main() у меня два цикла: второй цикл просто сохраняет картинку на диск, а первый цикл — проходит по всем пикселям картинки, испускает луч из камеры через этот пиксель, и смотрит, не пересекается ли этот луч с нашей сферой.

Внимание, основная идея статьи: если в прошлой статье мы аналитически считали пересечение луча и сферы, то сейчас я его считаю численно. Идея простая: сфера имеет уравнение вида x^2 + y^2 + z^2 — r^2 = 0; но вообще функция f(x,y,z) = x^2 + y^2 + z^2 — r^2 определена во всём пространстве. Внутри сферы функция f(x,y,z) будет иметь отрицательные значения, а снаружи сферы положительные. То есть, функция f(x,y,z) задаёт для точки (x,y,z) расстояние (со знаком!) до нашей сферы. Поэтому мы просто будем скользить вдоль луча до тех пор, пока либо нам не надоест, либо функция f(x,y,z) станет отрицательной. Функция sphere_trace() именно это и делает.

Этап третий: примитивное освещение

Давайте закодим простейшее диффузуное освещение, вот такую картинку я хочу получить на выходе:



Как и в прошлой статье, для простоты чтения я сделал один этап = один коммит. Изменения можно смотреть тут.

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

Vec3f distance_field_normal(const Vec3f &pos) {
    const float eps = 0.1;
    float d = signed_distance(pos);
    float nx = signed_distance(pos + Vec3f(eps, 0, 0)) - d;
    float ny = signed_distance(pos + Vec3f(0, eps, 0)) - d;
    float nz = signed_distance(pos + Vec3f(0, 0, eps)) - d;
    return Vec3f(nx, ny, nz).normalize();
}

В принципе, конечно, поскольку мы рисуем сферу, то нормаль можно получить гораздо проще, но я сделал так с заделом на будущее.

Этап четвёртый: давайте нарисуем паттерн на нашей сфере

А давайте нарисум какой-нибудь паттерн на нашей сфере, например, вот такой:



Для этого в предыдущем коде я изменил всего две строчки!

Как я это сделал? Разумеется, у меня нет никаких текстур. Я просто взял функцию g(x,y,z) = sin(x) * sin(y) * sin(z); она опять же определена во всём пространстве. Когда мой луч пересекает сферу в какой-то точке, то значение функции g(x,y,z) в этой точке мне задаёт цвет пикселя.

Кстати, обратите внимание на концентрические круги по сфере — это артефакты моего численного подсчёта пересечения.

Этап пятый: displacement mapping

Для чего я захотел нарисовать этот паттерн? А он мне поможет нарисовать вот такого ёжика:



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

Чтобы это сделать, достаточно изменить три строчки в нашем коде:

float signed_distance(const Vec3f &p) {
    Vec3f s = Vec3f(p).normalize(sphere_radius);
    float displacement = sin(16*s.x)*sin(16*s.y)*sin(16*s.z)*noise_amplitude;
    return p.norm() - (sphere_radius + displacement);
}

То есть, я изменил расчёт расстояния до нашей поверхности, определив его как x^2+y^2+z^2 — r^2 — sin(x)*sin(y)*sin(z). По факту, мы определили неявную функцию.

Этап шестой: очередная неявная функция

А почему это я оцениваю произведение синусов только для точек, лежащих на поверхности нашей сферы? Давайте переопределим нашу неявную функцию вот так:

float signed_distance(const Vec3f &p) {
    float displacement = sin(16*p.x)*sin(16*p.y)*sin(16*p.z)*noise_amplitude;
    return p.norm() - (sphere_radius + displacement);
}

Разница с предыдущим кодом совсем маленькая, лучше посмотреть дифф. Вот что получится в итоге:



Таким образом мы можем определять несвязные компоненты в нашем объекте!

Этап седьмой: псевдослучайный шум

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



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



А потом просто их просуммировать:



Подробнее прочитать можно здесь и здесь.

Давайте добавим немного кода, генерирующего этот шум, и получим такую картинку:



Обратите внимание, что в коде рендеринга я не изменил вообще ничего, изменилась только функция, которая «мнёт» нашу сферу.

Этап восьмой, финальный: добавляем цвет

Единственное, что я изменил в этом коммите, это вместо равномерного белого цвета я наложил цвет, который линейно зависит от величины приложенного шума:

Vec3f palette_fire(const float d) {
    const Vec3f   yellow(1.7, 1.3, 1.0); // note that the color is "hot", i.e. has components >1
    const Vec3f   orange(1.0, 0.6, 0.0);
    const Vec3f      red(1.0, 0.0, 0.0);
    const Vec3f darkgray(0.2, 0.2, 0.2);
    const Vec3f     gray(0.4, 0.4, 0.4);

    float x = std::max(0.f, std::min(1.f, d));
    if (x<.25f)
        return lerp(gray, darkgray, x*4.f);
    else if (x<.5f)
        return lerp(darkgray, red, x*4.f-1.f);
    else if (x<.75f)
        return lerp(red, orange, x*4.f-2.f);
    return lerp(orange, yellow, x*4.f-3.f);
}

Это простой линейный градиент между пятью ключевыми цветами. Ну а вот картинка!

Заключение

Эта техника трассировки лучей называется ray marching. Домашнее задание простое: скрестить предыдущий рейтрейсер с блэкджеком и отражениями с нашим взрывом, да так, чтобы взрыв ещё и освещал всё вокруг! Кстати, этому взрыву сильно не хватает полупрозрачности.