В данной работе решается задача моделирования поверхности открытого океана в реальном времени с использованием GPU для вычисления поля высот.
\documentclass[a4paper,12pt]{article}
\usepackage{cmap}
\usepackage[utf8]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{framed}
\usepackage{hyperref}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage[colorinlistoftodos]{todonotes}
\usepackage{wrapfig}
\usepackage{lipsum}
\usepackage{listings}
\usepackage{color}
\usepackage{indentfirst}
\usepackage{times}
\usepackage{textcomp}
\definecolor{mygray}{rgb}{0.4,0.4,0.4}
\definecolor{mygreen}{rgb}{0,0.8,0.6}
\definecolor{myorange}{rgb}{1.0,0.4,0}
\lstdefinestyle{customc}{
belowcaptionskip=1\baselineskip,
breaklines=true,
frame=L,
xleftmargin=\parindent,
language=C,
showstringspaces=false,
basicstyle=\footnotesize\ttfamily,
keywordstyle=\bfseries\color{green!40!black},
commentstyle=\itshape\color{purple!40!black},
identifierstyle=\color{blue},
stringstyle=\color{orange},
numbers=left,
numbersep=12pt,
numberstyle=\small\color{mygray},
}
\lstset{escapechar=@,style=customc}
\newcommand{\HRule}{\rule{\linewidth}{0.5mm}}
\begin{document}
\input{title.tex}
\newpage
\tableofcontents
\newpage
\section{Введение}
\subsection{Исследуемая область}
Исследуемой областью является одно из перспективных направлений в компьютерном моделировании -- симуляция поверхности океана. Кажущееся интуитивным предположение, что поставленная задача эквивалентна задаче моделирования жидкости, на самом деле является неверным.
Моделирование жидкости -- область компьютерной графики, использующая средства вычислительной гидродинамики для реалистичного моделирования, анимации и визуализации жидкостей, газов, взрывов и других связанных с этим явлений, в то время как задача симуляции поверхности жидкости является задачей моделирования волны, что сильно упрощает используемый математический аппарат.
\subsection{Перспективность исследований}
Активному развитию данной области способствуют кинематограф и игровая индустрия, которые используют разработанные методы в своих продуктах. Так например в данной работе частично реализована модель океанской поверхности, которая была использована в фильме "Титаник".
Хотя CPU с каждым годом увеличивают свою мощность, ее все равно еще не достаточно для использование реалистичных моделей поверхности жидкости в компьютерных играх. Решением данной проблемы может стать одно из перспективных направлений исследований в последнее время -- распределенные вычисления на GPU. Именно этот подход используется в данной работе для увеличения скорости работы программы.
Если для кинематографа скорость не является ключевым фактором, то в игровой индустрии -- наоборот. Поэтому, пока у пользователей нет возможности использовать необходимую мощность для просчета реалистичной модели воды, данная область исследований останется перспективной, так как с одной стороны необходимо улучшать внешний вид моделей и их физическое поведение, так с другой стороны -- нельзя выходить за рамки возможного количества вычислений, т.е. необходимо разрабатывать новые модели и более оптимизированные алгоритмы.
\subsection{Постановка задачи}
В данной работе решается задача моделирования поверхности открытого океана в реальном времени с использованием GPU для вычисления поля высот.
\subsection{Возможные методы решения}
Существует много подходов к моделированию и анимации поверхности воды. Большая часть из них основана на аналитической модели суперпозиции волн, и решение здесь либо задается сразу в виде линейной комбинации тригонометрических функций со специально подобранными коэффициентами, либо получается в результате обратного преобразования Фурье со специально заданным спектром. Выбор спектра и определяет сложность, реалистичность и детализацию модели. Такие модели могут быть как очень сложными и дорогими с вычислительной точки зрения, так и довольно простыми.
Моделированием водной поверхности занималось достаточное количество разработчиков, наиболее успешным среди которых был Тессендорф (Tessendorf). Именно его модель частично реализована в данной работе.
\subsection{Используемые технологии}
\begin{itemize}
\item Программа написана на OpenGL 4.1 с использованием GLFW3 в качестве GUI.
\item Для вычисления поля высот используется технология CUDA и библиотека CUFFT для вычисления обратного преобразования Фурье.
\item Библиотека SDL2\_image и Devil необходимы для загрузки текстур для шейдера.
\item В качестве библиотеки линейной алгебры используется библиотека GLM, специально разработанная для использования с OpenGL.
\end{itemize}
\newpage
\section{Используемые математические модели} \label{sec:math}
\subsection{Модель волны}
В данной работе рассморен статистический метод создания волны. Модели такого типа основаны на возможности разложения поля высоты волны в сумму синусоид и косинусоид с использование случайных величин для генерации амплитуд частот. Вычислительно более выгодно использовать БПФ для вычисления этих сумм.
Поле высот с использованием БПФ можно представить в виде сумм синусоид с сложными, изменяющимися со временем амплитудами:
\begin{equation} \label{eq:fft}
h(\mathbf{x},t) = \sum_{\mathbf{k}} \tilde{h}(\mathbf{k},t) \exp{(ik \cdot \mathbf{x})}
\end{equation}
где $\mathbf{x} = (x_1, x_2)$ -- положение точки на двумерной сетке, $t$ -- время, а $\mathbf{k} = (k_1, k_2), k_1 = \dfrac{2 \pi n}{L_1}, k_2 = \dfrac{2 \pi m}{L_2}, -\dfrac{N}{2} \le n < \dfrac{N}{2}, -\dfrac{M}{2} \le m < \dfrac{M}{2}.$
Величины амплитуды однозначно задают все поле высот. Идея статистического метода заключается в создании случайных наборов амплитуд, удовлетворяющих эмпирическим законам океанографии.
Океанографические исследования показали, что уравнение \ref{eq:fft} является достаточно точным представлением ветряных волн, возникающих в открытом океане.
$\tilde{h}(\mathbf{k},t)$ можно представить в виде:
\begin{equation} \label{eq:ht}
\tilde{h}(\mathbf{k},t) = \tilde{h_0}(\mathbf{k}) \exp{(i\omega(k)t)} + \tilde{h_0^{*}}(-\mathbf{k}) \exp{(-i\omega(k)t)}
\end{equation}
где $\tilde{h_0}(\mathbf{k})$ - амплитуды поля высоты в момент времени $t=0$, которые задаются по формуле:
\begin{equation} \label{eq:h0}
\tilde{h_0}(\mathbf{k}) = \dfrac{1}{\sqrt{2}} (\xi_r + i \xi_i) \sqrt{P_h(\mathbf{k})}
\end{equation}
где $ \xi_r,\xi_i \text{\texttildelow} N(0,1)$, а $P_h(\mathbf{k})$ - спектр Филлипса, который задается эмпирической формулой:
\begin{equation} \label{eq:phill}
P_h(\mathbf{k}) = A \dfrac{\exp{(-1/(kL)^2)}}{k^4}|\hat{\mathbf{k}} \cdot \hat{\omega}|^2
\end{equation}
в которой $L = \dfrac{V^2}{g}.$
Именно спектр Филлипса является необходимым эмпирическим законом океанографии, благодаря которому волны становятся похожи на настоящие. Данная формула была получена с помощью эксперементальных данных разного вида, собиравшихся в течение продолжительного количества времени.
\subsection{Модель освещения}
В качестве модели освещения используется модель Блинна-Фонга. Несмотря на то, что существую более точные модели освещения, она является стандартом в компьютерной графике.
Основная идея модели Блинна-Фонга заключается в предположении, что освещенность каждой точки тела разлагается на 3 компоненты:
\begin{enumerate}
\item фоновое освещение (ambient),
\item рассеянный свет (diffuse),
\item бликовая составляющая (specular).
\end{enumerate}
Свойства источника определяют мощность излучения для каждой из этих компонент, а свойства материала поверхности определяют её способность воспринимать каждый вид освещения.
\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{img/blinn-phong.png}
\caption{Модель Блинна-Фонга}
\label{fig:blinn-phong-model}
\end{figure}
Фоновое освещение это постоянная в каждой точке величина надбавки к освещению. Вычисляется фоновая составляющая освещения как:
\begin{equation}
I_a = k_a i_a, \text{ где}
\end{equation}
\begin{itemize}
\item $I_a$ -- фоновая составляющая освещенности в точке;
\item $k_a$ -- свойство материала воспринимать фоновое освещение;
\item $i_a$ -- мощность фонового освещения.
\end{itemize}
Рассеянный свет при попадании на поверхность рассеивается равномерно во все стороны. При расчете такого освещения учитывается только ориентация поверхности (нормаль) и направление на источник света. Рассеянная составляющая рассчитывается по закону косинусов (закон Ламберта):
\begin{equation}
I_d = k_d (\vec{L} \cdot \vec{N}) i_d, \text{ где}
\end{equation}
\begin{itemize}
\item $I_d$ -- рассеянная составляющая освещенности в точке,
\item $k_d$ -- свойство материала воспринимать рассеянное освещение,
\item $\vec{L}$ -- направление из точки на источник,
\item $\vec{N}$ -- вектор нормали в точке,
\item $i_d$ -- мощность рассеянного освещения.
\end{itemize}
\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{img/bpvectors.png}
\caption{Необходимые векторы для модели Блинна-Фонга}
\label{fig:bpvectors}
\end{figure}
Зеркальный свет при попадании на поверхность подчиняется следующему закону: “Падающий и отраженный лучи лежат в одной плоскости с нормалью к отражающей поверхности в точке падения, и эта нормаль делит угол между лучами на две равные части”. Т.о. отраженная составляющая освещенности в точке зависит от того, насколько близки направления на наблюдателя и отраженного луча. Это можно выразить следующей формулой:
\begin{equation}
I_s = k_s (\vec{H} \cdot \vec{N})^{\beta} i_s, \text{ где}
\end{equation}
\begin{itemize}
\item $I_s$ -- зеркальная составляющая освещенности в точке,
\item $k_s$ -- коэффициент зеркального отражения,
\item $\vec{H} = \dfrac{\vec{L} + \vec{V}}{|\vec{L} + \vec{V}|}$ -- ориентация площадки, на которой будет максимальное отражение,
\item $\vec{N}$ -- вектор нормали в точке,
\item $i_s$ -- мощность зеркального освещения,
\item $\beta$ -- коэффициент блеска, свойство материала.
\end{itemize}
Именно зеркальное отражение представляет наибольший интерес, но в то же время его расчет требует больших вычислительных затрат. При фиксированном положении поверхности относительно источников света фоновая и рассеянные составляющие освещения могут быть просчитаны единожды для всей сцены, т.к. их значение не зависит от направления взгляда. С зеркальной составляющей этот фокус не сработает и придется пересчитывать её каждый раз, когда взгляд меняет свое направление.
Во всех вычислениях выше, для рассеянной и зеркальной компонен, если скалярное произведение в правой части меньше нуля, то соответствующая компонента освещенности полагается равной нулю.
\newpage
\section{Детали реализации} \label{sec:code}
\subsection{Скайбокс}
Скайбокс -- объект в трёхмерной графике, играющий роль неба и горизонта. Представляет собой несложную трёхмерную модель (как правило, куб), с внутренней стороны которой натянута текстура неба -- "кубическая текстура".
\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{img/skybox.png}
\caption{Реализованный скайбокс}
\label{fig:skybox}
\end{figure}
Обработка трёхмерной графики требует много вычислительной работы, поэтому "честно" просчитывать объекты, находящиеся на горизонте, было бы расточительством. К тому же, трёхмерное аппаратное обеспечение имеет Z-буферы, которые из-за ограниченной разрядности отбрасывают всё, что находится далеко от камеры.
Поэтому удалённые объекты изображаются крайне примитивно: в виде куба, шесть граней которого — текстуры неба и горизонта. Если отобразить этот куб так, чтобы камера находилась точно в центре, будет казаться, что через камеру действительно видны небо и горизонт.
В данной компьютерной симуляции скайбокс реализован в виде куба с 6-ю различными текстурами, которые накладываются на каждую из граней.
Фрагментный шейдер, используемый для текстурирования имеет тривиальный вид:
\begin{lstlisting}
#version 410 core
in vec3 texCoord;
out vec4 fColor;
uniform samplerCube cubemap;
void main (void) {
fColor = texture(cubemap, texCoord);
}
\end{lstlisting}
\subsection{Поле высот}
Для того, чтобы эффективно реализовать симуляцию, используя распределенные вычисления на GPU, необходимо принимать во внимание следующие 2 факта:
\begin{enumerate}
\item Копирование данных со стороны CPU на сторону GPU является очень затратной операцией,
\item GPU имеет ограниченный объем памяти.
\end{enumerate}
Вторая проблема решиется на уровне постановки задачи -- предполагается, что вся необходимая для симуляции информация, полностью помещается в памяти GPU.
Для того, чтобы решить первую проблему, необходимо написать программу таким образом, чтобы значения, вычисленные на GPU не передавались обратно на сторону CPU, а сразу копировались в OpenGL буффер. Для того, чтобы реализовать такое поведение, CUDA предоставляет 4 функции:
\begin{itemize}
\item cudaGraphicsGLRegisterBuffer -- регистрирует вспомогательную CUDA структуру, которая может обращаться к OpenGL буферу,
\item cudaGraphicsMapResources -- соединяет вспомогательную структуру с OpenGL буфером,
\item cudaGraphicsResourceGetMappedPointer -- возвращает указатель, с помощью которого можно скопировать данные в OpenGL буфер напрямую,
\item cudaGraphicsUnmapResources -- закрывает соединение вспомогательной структуры с OpenGL буфером.
\end{itemize}
\begin{figure}[h]
\centering
\includegraphics[width=\textwidth]{img/screenshot2.png}
\caption{Пример сгенерированного поля высот}
\label{fig:heightfield}
\end{figure}
Вычисление поля высот происходит каждый раз, когда рендерится изображение.
\begin{lstlisting}
// generate wave spectrum in frequency domain
cudaGenerateSpectrumKernel(d_h0, d_ht, spectrum, meshSize, meshSize, curTime, patchSize);
// execute inverse FFT to convert to spatial domain
checkCudaErrors(cufftExecC2C(fftPlan, d_ht, d_ht, CUFFT_INVERSE));
// update heightmap values in vertex buffer
checkCudaErrors(cudaGraphicsMapResources(1, &cuda_heightVB_resource, 0));
checkCudaErrors(cudaGraphicsResourceGetMappedPointer((void **)&g_hptr, &num_bytes, cuda_heightVB_resource));
cudaUpdateHeightmapKernel(g_hptr, d_ht, meshSize, meshSize);
checkCudaErrors(cudaGraphicsUnmapResources(1, &cuda_heightVB_resource, 0));
\end{lstlisting}
В данном коде функции cudaGenerateSpectrumKernel, cufftExecC2C, cudaUpdateHeightmapKernel выполняются на GPU. Функция cufftExecC2C является библиотечной функцией, которая производит в данном случае обратное преобразование Фурье, а оставшиеся функции написаны вручную и имеют следующую реализацию:
\begin{lstlisting}
// generate wave heightfield at time t based on initial heightfield and dispersion relationship
__global__ void generateSpectrumKernel(float2 *h0, float2 *ht,
unsigned int in_width, unsigned int out_width,
unsigned int out_height, float t, float patchSize)
{
unsigned int x = blockIdx.x*blockDim.x + threadIdx.x;
unsigned int y = blockIdx.y*blockDim.y + threadIdx.y;
unsigned int in_index = y*in_width+x;
unsigned int in_mindex = (out_height - y)*in_width + (out_width - x); // mirrored
unsigned int out_index = y*out_width+x;
float2 k;
k.x = (-(int)out_width / 2.0f + x) * (2.0f * CUDART_PI_F / patchSize);
k.y = (-(int)out_width / 2.0f + y) * (2.0f * CUDART_PI_F / patchSize);
float k_len = sqrtf(k.x*k.x + k.y*k.y);
float w = sqrtf(9.81f * k_len);
if((x < out_width) && (y < out_height)) {
float2 h0_k = h0[in_index];
float2 h0_mk = h0[in_mindex];
ht[out_index] = complex_add(complex_mult(h0_k, complex_exp(w * t)),
complex_mult(conjugate(h0_mk), complex_exp(-w * t)));
}
}
// update height map values based on output of FFT
__global__ void updateHeightmapKernel(float *heightMap,
float2 *ht, unsigned int width)
{
unsigned int x = blockIdx.x * blockDim.x + threadIdx.x;
unsigned int y = blockIdx.y * blockDim.y + threadIdx.y;
unsigned int i = y * width +x;
float sign_correction = ((x + y) & 0x01) ? -1.0f : 1.0f;
heightMap[i] = ht[i].x * sign_correction;
}
\end{lstlisting}
Функция generateSpectrumKernel генерирует поле высот из начального поля высот и пройденного времени, а updateHeightmapKernel -- вспомогательная функция, которая реализует смещение точек после обратного преобразования Фурье.
Не менее интересной является реализация функции, которая генерирует начальное поле высот:
\begin{lstlisting}
float Waves::phillips(float Kx, float Ky)
{
float k_squared = Kx * Kx + Ky * Ky;
if (k_squared == 0.0f) {
return 0.0f;
}
float L = windSpeed * windSpeed / g;
float k_x = Kx / sqrtf(k_squared);
float k_y = Ky / sqrtf(k_squared);
float w_dot_k = k_x * windDir.x + k_y * windDir.y;
float phillips = A * expf(-1.0f / (k_squared * L * L))
/ (k_squared * k_squared) * w_dot_k * w_dot_k;
// filter out waves moving opposite to wind
if (w_dot_k < 0.0f) {
phillips *= dirDepend; // dir_depend;
}
return phillips;
}
void Waves::generateH0()
{
for (unsigned int y = 0; y < spectrum; ++y) {
for (unsigned int x = 0; x < spectrum; ++x) {
float kx = (-(int)meshSize / 2.0f + x) * (2.0f * CUDART_PI_F / patchSize);
float ky = (-(int)meshSize / 2.0f + y) * (2.0f * CUDART_PI_F / patchSize);
float P = sqrtf(phillips(kx, ky));
float Er = gauss();
float Ei = gauss();
float h0_re = Er * P * CUDART_SQRT_HALF_F;
float h0_im = Ei * P * CUDART_SQRT_HALF_F;
int i = y * spectrum + x;
h_h0[i].x = h0_re;
h_h0[i].y = h0_im;
}
}
}
\end{lstlisting}
\subsection{Освещение}
Реализация модели освещения Блинна-Фонга имеет стандартный вид:
\begin{lstlisting}
/* vertex shader */
#version 410 core
layout(location = 0) in vec4 meshPos;
layout(location = 1) in float height;
layout(location = 2) in vec2 slope;
uniform mat4 PVM;
uniform vec3 lightPos;
uniform vec3 eyePos;
out vec3 l;
out vec3 h;
out vec3 n;
out vec3 r;
out vec4 pos;
void main() {
vec3 lp = abs(lightPos);
vec3 p = vec3(meshPos.x, 1e+2 * height, meshPos.z);
gl_Position = PVM * vec4(p, 1.0);
p.x = p.x - 1000; p.z = p.z - 1000;
p.y = p.y - 500;
pos = vec4(p, 1.0);
l = normalize(lp - p);
vec3 v = normalize(eyePos - p);
h = normalize((v + l) / length(v + l));
n = normalize(cross( vec3(0.0, slope.y, 1.0 / 256), vec3(1.0 / 256, slope.x, 0.0)));
r = reflect(-l, n);
}
\end{lstlisting}
\begin{lstlisting}
/* fragment shader */
#version 410 core
in vec4 pos;
out vec4 fColor;
uniform vec3 sourceColor;
uniform vec3 diffColor;
uniform vec3 specColor;
uniform vec3 lightPos;
uniform vec3 eyePos;
in vec3 l;
in vec3 h;
in vec3 n;
in vec3 r;
uniform vec3 Ka;
uniform vec3 Kd;
uniform vec3 Ks;
uniform float alpha;
vec3 BlinnPhongModel()
{
return Ka * sourceColor +
Kd * max(dot(n, -l), 0.0) * diffColor +
Ks * max(pow(dot(n, h), alpha), 0.0) * specColor;
}
void main (void) {
vec3 BlinnPhong = exp(-0.8 + 1.2*abs(pos.x/3000+pos.z/3000)) * BlinnPhongModel();
fColor = vec4(BlinnPhong, 0.9);
}
\end{lstlisting}
В фрагментном шейдере используется $\alpha$-канал не равный единице, в данной реализации $\alpha = 0.9$, чтобы была возможность сквозь воду просматривать дно.
\newpage
\section{Подведение итогов}
\subsection{Полученные результаты}
Результат работы программы виден на следующем скриншоте:
\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{img/screenshot.png}
\caption{Скриншот работы программы}
\label{fig:screenshot}
\end{figure}
Полученный результат напоминает жидкость, но океан имеет более сложную текстуру. Можно было бы продолжить исследовать проблему симуляции поверхности океана и добавить такие эффекты, как порывистые волны, брызги, пену, более подходящую для океана модель освещения, интерференцию волн, отражение мира на поверхности воды, каустический эффект и многое другое, которые бы улучшили внешний вид воды, но тема слишком сложная для любительского ознакомления.
\subsection{Вывод}
Симуляция поверхности океана - очень интересный, важный и активно развивающийся раздел моделирования. С помощью распределенных вычислений на GPU можно добиться вычисления очень большой площади поверхности в режиме реального времени с неплохой точностью. В данной работе была реализована самая простая статистическая модель волны, однако даже эта модель позволяет просчитать поведение волн с хорошей точностью.
\newpage
\begin{thebibliography}{9}
\addcontentsline{toc}{section}{\refname}
\bibitem{tessendorf} Tessendorf, J. 2001. Simulating Ocean Water. ACM SIGGRAPH.
\bibitem{nvidia} NVIDIA. 2011. Ocean Surface Simulation. NVIDIA Graphics SDK 11 Direct3D.
\bibitem{gpuOceanSim} Chin-Chih Wang, Jia-Xiang Wu, Chao-En Yen, Pangfeng Liu, Chuen-Liang Chen. Ocean Wave Simulation in Real-time using GPU.
\bibitem {farber2011cuda} Farber, R. 2011, CUDA Application Design and Development, Applications of GPU computing, Elsevier Science
\bibitem{wen-mei} David B. Kirk, Wen-mei W. Hwu. 2012, Programming Massively Parallel Processors: A Hands-on Approach, Newnes
\end{thebibliography}
\end{document}
\usepackage[english,russian]{babel}