Симуляция поверхности океана в реальном времени с использованием 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}