Пишем симуляцию океанских волн на языке программирования Ü

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

Сборка

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

По-хорошему её описание требует отдельной статьи, но пока её нету, я опишу вкратце здесь, как она работает. Это отдельный исполняемый файл, который читает описание проекта в особом формате и запускает компилятор, чтобы собрать итоговую программу. Программы описываются не на каком-то отдельном языке, а на самом же Ü, в декларативном виде. Система сборки читает этот файл сборки (именуемый build.u), компилирует его и вызывает функцию для получения информации о проекте. Далее она сама заботится о том, как его собирать, осуществляет поиск зависимостей, инкрементальную сборку, распараллеливание и т. д.

Пример файла простейшего проекта:

import "/build_system.uh"

fn GetPackageInfo( BK::BuildSystemInterface &mut build_system_interface ) : BK::PackageInfo
{
    ust::ignore_unused( build_system_interface );

    var BK::BuildTarget mut target{ .target_type = BK::BuildTargetType::Executable };
    target.source_files.push_back( "main.u" );
    target.name= "hello_world";

    return BK::PackageInfo{ .build_targets= ust::make_array( move(target) ) };
}

Запускается сборка одной строкой:

u.._build_system build

Зависимости

Кроме непосредственно языка/компилятора/системы сборки/стандартной библиотеки нужен ещё код для утилитарных нужд, вроде создания окна, рисования в него, получения пользовательского ввода. Для этих целей я выбрал библиотеку SDL2, которая кроcсплатформенна и позволяет достаточно просто всё это делать.

Внимательный читатель может заметить, что эта библиотека написана на Си. Как же тогда можно её использовать в Ü? Скомпоновать код на Ü с SDL2 не составляет никаких проблем. Система сборки предоставляет возможность указать внешние библиотеки, с которыми нужна компоновка. Делается это примерно так:

target.external_libraries.push_back( "-lSDL2" );

А как вызывать код из такой библиотеки из Ü кода? Тут на помощь приходит утилита, которая позволяет создавать Ü биндинги для заголовочных файлов Си. Называется она u.._cpp_header_converter и входит в поставку Ü. Её название немного вводит в заблуждение, т. к. с C++ кодом она почти не умеет работать, зато с Си кодом работает прекрасно. Она создаёт Ü код объявлений внешних функций, структур, перечислений, некоторых define-констант, которые вычитываются из переданного заголовочного файла Си. Написана данная утилита на основе clang и посему не имеет проблем с парсингом Си кода.

Так вот, чтобы работать с SDL2, нужно всего лишь обработать этой утилитой заголовочный файл SDL.h:

u.._cpp_header_converter /usr/include/SDL2/SDL.h -o SDL.uh -- -std=c11

Делается это конечно не руками. Вместо этого в описании проекта в файле build.u указывается особый шаг сборки, который выполняется системой сборки по-необходимости. Далее полученный таким образом файл просто импортируется и код из него может быть вызван, примерно так:

import "/SDL.uh"

fn GetCurrentTime() : u64
{
    return unsafe( SDL_GetPerformanceCounter() );
}

Вспомогательный математический функционал

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

Наличие перегрузки операторов в языке Ü положительно сказалось на эргономике этих классов. Я перегрузил арифметические операторы для этих классов — там, где это имеет смысл.

Пример перегрузки операторов:

template</type T/>
struct Complex
{
    type ThisType= Complex</T/>;

    fn constructor( T in_re, T in_im )
        ( re= in_re, im= in_im )
    {}

    op+( ThisType& l, ThisType& r ) : ThisType
    {
        return ThisType( l.re + r.re, l.im + r.im );
    }

    op*( ThisType& l, ThisType& r ) : ThisType
    {
        return ThisType( l.re * r.re - l.im * r.im, l.re * r.im + r.re * l.im );
    }

    op*( ThisType& c, T scalar ) : ThisType
    {
        return ThisType( c.re * scalar, c.im * scalar );
    }

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

Интерфейс класса генератора псевдослучайных чисел:

class PseudoRandomGenerator
{
public:
    var u32 c_max= 32767u;
    var u32 c_max_plus_one= 32768u;

public:
    fn constructor() = default;
    fn constructor( u32 seed );

    // Produces value in range [ 0, c_max ]
    fn Next( mut this ) : u32;

private:
    u32 state_= 0u;
}

Кроме того я реализовал на его основе код получения пар псевдослучайных чисел с нормальным распределением — на основе marsaglia polar method. Позже будет пояснено, зачем это нужно.

Волны

Итак, у нас стоит задача просимулировать и нарисовать поверхность океана. А что самое заметное на водной поверхности? Конечно же волны!

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

_{Простая синусоида на плоскости. Так океан конечно не выглядит, но это только начало.}

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

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

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

_{Визуализация спектра. Цвета обозначают фазу, чем больше яркость, тем больше амплитуда.}

Итак, у нас есть спектр, что делать с ним дальше? А дальше мы можем получить карту высот волн из него. Наивный подход заключался бы в том, чтобы посчитать синусоиды для всех направлений и длин волн и сложить их. Но поскольку их сильно много, это было бы сильно медленно, нужно было бы вычислять множество тригонометрических функций для каждого текселя карты высот. Поэтому вместо этого используется более быстрый подход — двухмерное быстрое обратное преобразование Фурье.

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

_{Итоговая карта высот, составленная из множества элементарных волн по спектру из предыдущего изображения.}

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

Анимация волн во времени

Генерация карты высот поверхности воды это конечно интересно, но не достаточно. Эта поверхность должна со временем изменяться, а не абы как, а чтобы волны двигались, причём более-менее естественным образом.

Движение любой волны можно выразить как изменение её фазы со временем. В моём случае такое изменение выразимо как умножение комплексной амплитуды для каждой волны на комплексное число, получаемое возведением 𝑒 в степень 𝑖, помноженной на угол поворота. Звучит это мудрёно, но на самом деле это возведение в степень есть ничто иное, как вычисление косинуса для действительной части и синуса для мнимой.

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

Модификация формы волн

Пока что я считал все волны синусоидальными. Но в действительности форма волн на воде несколько иная, синусоидальное приближение справедливо только для волн, у которых амплитуда сильно меньше их длины. Для более высоких волн их форма это не синусоида, а трохоида (волна Герстнера). Данная форма волны получается при движении каждой точки поверхности воды по круговой траектории, как вверх-вниз, так и в стороны.

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

_{Карта высот одиночной волны Герстнера.}

_{Она же в трёх измерениях — для большей наглядности.}

В принципе можно было бы рисовать поверхность океана, используя регулярную сетку полигонов со сдвигом вершин на посчитанные вышеизложенным образом величины. Но у меня отрисовка работает с картой высот (см. ниже), а значит мне нужно получить модифицированную вычисленными ранее сдвигами карту высот. Делается это отдельным шагом, который осуществляет примерно следующее: для каждой позиции X и Y в исходной карте высот значение высоты записывается в позицию X + dx, Y + dy итоговой карты высот. Делается это в два прохода — сначала со сдвигом по оси X а затем по оси Y. Кроме этого значения линейно интерполируются и гарантируется, что итоговая карта высот не содержит прорех.

_{Итоговая карта высот с учётом боковых сдвигов волн.}

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

Раскраска

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

Базовый цвет воды это оттенок синего. В местах, где наклон поверхности значительный (определяется нормалью), подмешивается немного зелёного цвета — дабы имитировать просвечивание воды солнцем. Кроме того на гребнях волн (в местах, где значение карты высот в точке выше среднего значения соседей) к цвету воды подмешивается белый цвет — чтобы имитировать морскую пену.

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

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

_{Итоговая раскраска поверхности.}

Рисование

Итак, у меня есть карта высот с цветами. Как же её отобразить на экране? Можно было бы это всё нарисовать полигональной сеткой с использованием растеризации на GPU, используя какой-нибудь OpenGL или Vulkan. Но мне это кажется слишком банальным, посему я рисую карту высот на CPU и не через растеризацию полигонов. Производительности CPU не хватило бы для растеризации количества полигонов, достаточного для детального отображения океанских волн.

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

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

Он также даёт искажения формы отображаемой карты высот, если камера не смотрит точно в направлении горизонта, что особо заметно при сильных наклонах. Но с этим можно смириться, многие старые игры имели такую неправильную перспективу (Duke Nukem 3D, Chasm: the Rift).

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

_{Итоговая поверхность воды с синусоидальной формой волн.}

_{Итоговая поверхность воды со скорректированной формой волн.}

Ещё снимки экрана на разных масштабах:

Результаты и выводы

Карта высот поверхности океана вычисляется в разрешении 512x512 и рисуется в окно разрешением 1024x768 пикселей. С данными настройками удаётся добиться частоты кадров около 30 Гц. Большую часть времени кадра (более 3/4) занимает вычисление карты высот и её раскраска, из этого времени около половины уходит на вычисление двухмерного быстрого обратного преобразования Фурье для трёх карт — смещению волн по высоте и по осям X и Y.

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

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

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

Благодарности

Благодарю пользователя Youtube под псевдонимом Acerola за видео, повествующее о его реализации симуляции океана. Оно меня весьма вдохновило на создание моей версии такой симуляции.

Отдельно хотел бы похвалить использованную в ходе процесса разработки IDE ecode. Она имеет встроенную поддержку языка программирования Ü — подсветку синтаксиса, поддержку языкового сервера, отладки. Мои благодарности автору этой IDE за принятие запросов на слияние, реализующих эту поддержку.

Ссылки и использованные ресурсы

• Исходный код. Там же доступно демо для скачивания.
• Репозиторий языка Ü. Там же доступны сборки.
• Видео с демонстрацией
• Damien Hinsinger, Fabrice Neyret, Marie-Paule Can, Interactive Animation of Ocean Waves
• Ocean surface simulation