%=========================================================================
% Author: Martin Smutny, 2023

% For compilation piecewise (see projekt.tex), it is necessary to uncomment it and change
% \documentclass[../projekt.tex]{subfiles}
% \begin{document}

% 	\textbf{Using hyphens instead of dashes} \\
%	A dash is long and there should be whitespace surrounding it. Dash is very often used instead of a dot in sentences: ``This book -- published before the war -- is amazing.'' It is also used for ranges: 	``page 23--26'' or ``success rate of 3--5\,\%.'' More examples can be found in the language handbook \cite{prirucka}.
% Hyphens occur in our IT theses very rarely (at least they should). A good example would be phrases or joining subjects ``Rh-factor'', ``real-time'', ``propane-butane''.
% Print the thesis and have it bound. As soon as you open the bound thesis on a random page, you'll find an error. Honestly, don't worry about it, that's just how it is. Nothing is perfect, so \uv{just ship it}\footnote{\url{http://blog.igor.szoke.cz/2011/08/zkoumejte-bezte-za-bod-odkud-neni.html}} \cite{rady}.

\def\*#1{\mathbf{#1}}

% =============================================================
% Notation
\newcommand{\scatterCoef}{\sigma_{s}}
\newcommand{\absorbCoef}{\sigma_{a}}
\newcommand{\extinctCoef}{\sigma_{t}}
\newcommand{\backScatter}{\sigma_{bs}}
\newcommand{\forwardScatter}{\sigma_{fs}}
\newcommand{\diffuseAttenDown}{K_d}
\newcommand{\phaseFnc}{p}

% =============================================================


\chapter{Introduction}
\label{ch_intro}

%\todo{navodit predstavu o cem bude, zaujmout, popsat vlastnimi slovy reseny problem, proc resen, a jak resen}
%\todo{general introduction to the topic, \\
%only well-known words, no scientific terminology if possible, \\
%broader context of a problem, \\
%why the thesis is important in the given field \\
%what is the relationship with the given field, and other important facts \\
%short summary: past advances in the field, current state, the problem, and the reason why, and how solved
%}
%\todo{no bullshit about faster hardware and whatnot}

%\todo{uvod do problematiky, uvest do sirsiho kontextu} \\


One of the goals of computer graphics is to reproduce natural phenomena in the digital world. Among the natural phenomena is undeniably are oceanic environments. Specifically, the typical green and blue hues of oceanic waters that can be viewed from above and below the waters, as shown in the photos of Figure~\ref{fig:oceans-example}. These typical colors are caused by the complex and highly-variable optical properties in these waters. Sunlight that penetrates into through the ocean surface into the waters causes well-known phenomenon, called \emph{caustics}. 

% TODO figure 
\begin{figure}[H]
\centering
\begin{subfigure}{.40\textwidth}
    \centering
    \includegraphics[width=1.0\linewidth]{figures/ocean-optics/pexels-flo-dahm-1139040.jpg}
\end{subfigure}
\begin{subfigure}{.40\textwidth}
    \centering
    \includegraphics[width=1.0\linewidth]{figures/ocean-optics/pexels-zir-yu-8395909.jpg}
\end{subfigure}
\caption{The left image is a~photo of a~tropic island, by Flo Dahm, pexels.com. The right image is a~photo of a~coral reef in Philippines, by Zir YU, pexels.com.
}
\label{fig:oceans-example}
\end{figure}

Reproducing the visual complexity of underwater environments as close to the perceived reality as possible poses a~challenge in computer graphics. Offline methods in computer graphics deal with simulation of light as it propagates through such volumes in a~straightforward way. On the other hand, real-time methods must look for alternative approaches, that mostly consisted of making simplifications and crude approximations. 

The goal of this thesis is then to, first, investigate methods and techniques in the field of real-time rendering, that try to simulate and render oceanic environments in high quality. Generally, these approaches try to simulate the ocean surface and the interactions of light in the water volume. And second, to design and implement an application that is able to render an underwater scenery in high quality.

A~physically-based bio-optical model is designed to describe optical parameters of a~water volume of oceanic type with non-uniform vertical distribution. A~single scattering and multiple scattering approximation technique is proposed to render the volume through the use of ray marching. Finally, the ocean surface is simulated to incorporate physically-based underwater effects, such as caustics.

The rest of thesis is organized as follows. In Chapter~\ref{ch:light-scatter}, general concepts and principles of light scattering by volumes and surface that bound the water body are explained. Chapter~\ref{ch:optics} studies scattering in natural waters from the point of view of ocean optics. In Chapter~\ref{ch:ocean} the related works for ocean simulation and rendering are summarized, these consist of existing solutions for: ocean surface simulation, modeling of bio-optical parameters, participating media rendering and procedural generation. Lighting effects, such as caustics are subject to Chapter~\ref{cha:caustics}. Finally, Chapters~\ref{ch:design} and \ref{ch:impl} describe the designed and subsequently implemented solution.

\chapter{Fundamentals of Light Scattering Theory}
\label{ch:light-scatter}

This chapter studies the theoretical fundamentals of light scattering in media and from surfaces. The propagation of radiation due to the light being scattered is studied by the field of radiative transfer. It is the radiative transfer equation (RTE) that describes the propagation of light in media. In order to render such media, which is the goal of this thesis, the volume rendering equation (VRE) must be studied. The VRE expresses the RTE in a~form that is suitable for rendering media.

The chapter is then organized as follows: Section~\ref{sec:radiometry} introduces the radiometric quantities, which form the basis of the derivation of the rendering algorithms, that will be used throughout this thesis. Section~\ref{sec:light-media} introduces the concept of participating media and studies the radiative transfer in such media to arrive at the VRE. Section~\ref{sec:light-surface} develops the equations that describe how light is scattered from surfaces that bound the water body, specifically, the planar water surface interface and an opaque sea floor.


% ================================================================================
\section{Radiometric and Photometric Quantities}
\label{sec:radiometry}
%\todo{what is radiometry}

This section provides an overview of radiometric quantities used in rendering based on the work of Pharr et al.~\cite{PJH2016}.

\emph{Radiometry} is the study of measuring electromagnetic radiation. This radiation propagates as waves in an environment. Electromagnetic waves with different \emph{wavelengths} ($\lambda$)
%, that is the distance between two adjacent peaks with the same phase,
have different properties. There is a~wide range of wavelengths from gamma waves with less than 10\,pm in length to extreme low frequency (ELF) radio waves that are tens of thousands of kilometers long. For rendering, only a~small range of these wavelengths is of particular interest, and that is the range that is visible to humans: from approximately 400\,nm to 700\,nm (see Figure~\ref{fig:visible-range}).

% TODO figure visible_range
% TODO missing violet
\begin{figure}[ht]
    \centering
    \includegraphics[width=0.7\textwidth]{figures/theory/spectrum.pdf}
    \caption{The visible spectrum of wavelengths in the context of the full electromagnetic spectrum. (\emph{Image adapted from Real-Time Rendering, 4th edition~\cite{RTR2018:ch8}}.)}
    \label{fig:visible-range}
\end{figure}


Closely related field that is based on radiometric principles is \emph{radiative transfer}. It is the study of the transfer of radiant energy based on empirical observations. Radiative transfer operates at the level of \emph{geometric optics} (also known as \emph{ray optics}), where macroscopic properties of light are sufficient to describe how light interacts with objects much larger than its wavelength. Incorporating phenomena from wave optics models of light is also not uncommon.
%But, for the use of this thesis, geometric optics provides a~sufficient model for the description of light and light scattering.
% TODO explain models of light?



\ifdefined\EXPLAINENERGY
%  radiometry wasn’t originally derived from first principles using the physics of light but was built on an abstraction of light based on particles flowing through space. 
% although connections have since been made between radiometry and Maxwell’s equations, giving radiometry a solid basis in physics
To describe radiometric quantities, it is necessary to start from energy. Radiometry was originally built on an abstraction of light based on particles flowing through space, but also has a~solid basis in physics through Maxwell's equations. There sources of illumination that emit photons. Each photon at a~particular wavelength carries a~particular amount of energy~$Q$ measured in joules (J):
% TODO the equation is needed?
\begin{align*}
    Q = \frac{h c}{\lambda}  ,
\end{align*}
where $c$~is the speed of light, and $h$~is Planck's constant.
\fi
Radiometry was originally built on an abstraction of light based on particles flowing through space. These particles, called photons, are emitted from sources of illumination. Each photon at a~particular wavelength carries a~particular amount of \emph{energy}~$Q$ measured in joules [J]. There are four radiometric quantities that are central to rendering. These fundamental radiometric quantities effectively provide different ways of measuring energy of photons. All of them are generally wavelength-dependent, though that dependence is not explicit in this chapter. 

\ifdefined\SHOWFIGURESUP
\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/theory/flux.pdf}
    \caption{Radiant flux\\$\Phi$}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=0.7\linewidth]{figures/theory/irradiance.pdf}
    \caption{Irradiance\\$E(\*{x})$}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=0.4\linewidth]{figures/theory/radiantIntensity.pdf}
    \vspace{1em}
    \captionsetup{justification=centering}
    \caption{Radiant intensity\\$I(\omega)$}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=0.6\linewidth]{figures/theory/fieldRadiance.pdf}
    \label{fig:radiometric-quantities-radiance}
    \caption{Radiance\\$L(\*{x}, \omega)$}
\end{subfigure}
\label{fig:radiometric-quantities}
\caption{Fundamental radiometric quantities. (\emph{Images from Novák~\cite{Nov2014}.})
}
\end{figure}

\begin{wrapfigure}{r}{0.245\textwidth}
    \includegraphics[width=0.3\textwidth]{figures/theory/flux.pdf}
    \caption{Radiant flux $\Phi$ (\emph{Taken from Novák~\cite{Nov2014}.})
    }
    \label{fig:flux}
\end{wrapfigure}

\begin{wrapfigure}{r}{0.3\textwidth}
    \includegraphics[width=0.3\textwidth]{figures/PJH2016/flux.pdf}
    \caption{Radiant flux emanating from a~point light source measured at spheres surrounding the light.
    (\emph{Image from Pharr et al.~\cite{PJH2016}.})
    }
    \label{fig:flux-sphere}
\end{wrapfigure}
\fi

\vspace{1em}

\noindent\textbf{\emph{Radiant flux}} ($\Phi$) is the amount of energy passing through a~region of space over time, also known as \emph{power}, and is measured in watts [W] or joules per second [J/s].
%\ifdefined\SHOWEQUATIONS
It is defined as the limit of differential energy per differential time:
\begin{align}
    \Phi = \lim_{\Delta t \to 0} \frac{\Delta Q}{\Delta t} = \frac{\text{d}Q}{\text{d}t} .
\end{align}
%\fi
Flux is generally used to describe the total emission from a~light source. Figure~\ref{fig:flux-images} illustrates the flux and shows an example of flux being emitted from a~point light source measured at spheres that surround the light. The total amount of flux measured on the smaller sphere is the same as on the bigger sphere\,---\,although less energy passes through any local region of the larger sphere.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.5\linewidth]{figures/theory/flux.pdf}
    \vspace{1em}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.5\linewidth]{figures/PJH2016/flux.pdf}
\end{subfigure}
\caption{Radiant flux $\Phi$ shown in the left image (\emph{from Novák~\cite{Nov2014}}). The right image (\emph{from Pharr et al.~\cite{PJH2016}}) illustrates radiant flux emanating from a~point light source measured at spheres surrounding the light.
}
\label{fig:flux-images}
\end{figure}

%\begin{wrapfigure}{r}{0.3\textwidth}
%    \centering
%    \includegraphics[width=0.7\linewidth]{figures/theory/irradiance.pdf}
%    \caption{Irradiance\\$E(\*{x})$ (\emph{from Novák~\cite{Nov2014}}). }
%    \label{fig:irradiance}
%\end{wrapfigure}

\noindent\textbf{\emph{Irradiance}} ($E$) can be defined as the average density of power over a~finite area $A$, as illustrated in Figure~\ref{fig:irradiance}. This quantity is referred to as either irradiance, the area density of flux arriving at a~surface, or \emph{radiant exitance} ($M$), the area density of flux leaving a~surface. In the general case, these both terms can be defined as the differential flux per differential area at point $\*{x}$:
\begin{align}
    E(\*{x}) = \frac{\text{d} \Phi (\*{x})}{\text{d}A} .
\end{align}
For the point light source example in Figure~\ref{fig:flux-images}, the area density of radiant flux incident on a~sphere depends on its radius $r$: $E = \Phi/4\pi r^2$. In other words, irradiance at any point on the inner sphere is greater than the irradiance at any point on the outer sphere, since the surface area of the outer sphere is larger. 
%If the point source emits the same amount of power in all directions, then for such a~sphere with radius $r$ the irradiance is: $E = \frac{\Phi}{4\pi r^2 }$ .

In case of a~light source with area~$A$ and flux~$\Phi$ that is illuminating a~surface with normal~$n$, as illustrated in Figure~\ref{fig:lamberts-law}, the irradiance equation shows so-called \emph{Lambert's law}. If the light is shining directly down on the surface then the illuminated area on the surface~$A_1$ is equal to $A$. Irradiance at any point inside $A_1$~is then $E = \Phi / A$. However, when the light is shining at an angle on the surface, the illuminated area on the surface is larger. For small $A$, the area receiving flux, $A_2$, is roughly $A/\cos \theta$. For any point inside $A_2$,~the irradiance is then $E_2 = (\Phi \cos \theta) / A$, in [W/m$^{2}$]. 

\begin{figure}[H]
\centering
\begin{subfigure}{.4\textwidth}
    \centering
    \includegraphics[width=0.55\linewidth]{figures/theory/irradiance.pdf}
    \caption{Irradiance $E(\*{x})$ as the differential flux over a~differential area A. (\emph{Taken from Novák~\cite{Nov2014}}.)}
    %\vspace{0.1em}
    \label{fig:irradiance}
\end{subfigure}
\begin{subfigure}{.59\textwidth}
    \centering
    \includegraphics[width=0.6\textwidth]{figures/PJH2016/lambers_law.pdf}
    \caption{Lambert's law: irradiance incident at a~surface is proportional to the cosine of the angle between the light direction and the surface normal. (\emph{Image from Pharr et al.~\cite{PJH2016}.})
    }
    \label{fig:lamberts-law}
\end{subfigure}
\caption{The left image illustrates the irradiance, the right image explains the Lambert's law.}
\label{fig:irradiance-law}
\end{figure}


%\vspace{1em}
\begin{wrapfigure}{r}{0.3\textwidth}
    \centering
    \includegraphics[width=0.45\linewidth]{figures/theory/radiantIntensity.pdf}
    \caption{Radiant intensity $I(\omega)$ (\emph{from Novák~\cite{Nov2014}})}
    \label{fig:radiant-intensity}
\end{wrapfigure}

\noindent\textbf{\emph{Radiant intensity}} ($I$) is the angular density of power emitted by an infinitesimal light source centered within the unit sphere. This quantity is measured in [W/sr], where $sr$ (\emph{steradians}) is the unit to measure solid angles, for example, the entire unit sphere is a~solid angle of $4\pi$\,sr, and a~hemisphere is $2\pi$\,sr. Intensity over the entire sphere of directions is then:
\begin{align}
    I = \frac{\Phi}{4\pi} .
\end{align}
More generally, it is defined as the differential flux emitted along a~differential cone of directions (shown in Figure~\ref{fig:radiant-intensity}):
\begin{align}
    I = \frac{\text{d}\Phi}{\text{d}\omega} ,
\end{align}
where $\omega$ is a~normalized vector indicating the directions on the unit sphere.

\vspace{1em}

\begin{wrapfigure}{r}{0.3\textwidth}
    \centering
    \includegraphics[width=0.6\linewidth]{figures/theory/fieldRadiance.pdf}
    \caption{Field radiance (\emph{from Novák~\cite{Nov2014}})
    }
    \label{fig:radiance}
\end{wrapfigure}

\noindent\textbf{\emph{Radiance}} ($L$) measures irradiance (or radiant exitance) with respect to solid angles, it is a~measure of electromagnetic radiation in a~single ray. In comparison to irradiance, which gives the differential power per differential area at point $\*{x}$, radiance additionally distinguishes the directional distribution of power. It is defined as:
\begin{align}
    L(\*{x}, \omega) = \frac{\text{d}E_{\omega}(\*{x})}{\text{d}\omega} =
                        \frac{\text{d}\Phi(\*{x},\omega)}{\text{d}A^{\perp} \text{d}\omega } ,
\end{align}
where $E_{\omega}$ denotes irradiance at the surface that is perpendicular to the direction~$\omega$, hence the $\cos \theta$ term from Lambert's law is hence unnecessary. This type of radiance is sometimes more specifically called \emph{field radiance}. In terms of flux, radiance is the density of flux with respect to both area and solid angle, then $dA^{\perp}$ is the projected area $dA$ on a~hypothetical surface perpendicular to $\omega$ as shown in Figure~\ref{fig:radiance}. 

Radiance is the final and most fundamental of all radiometric quantities. If radiance is given, then all other quantities can be computed as integrals of radiance over areas and directions.

%\todo{NEEDED?}
In shading equations, radiance often appears in two forms: \emph{incident radiance} and \emph{outgoing radiance} (as illustrated in Figure~\ref{fig:incident-outgoing-radiance}). The incident radiance, denoted by $L_i (\*{x}, \omega)$, describes radiance arriving at point $\*{x}$~on the surface of an object from direction~$\omega$. The outgoing (or exitant) reflected radiance from the surface at that point is denoted by $L_o (\*{x}, \omega)$. In both cases, by convention, the direction vector $\omega$~points away from $\*{x}$. 

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.45\textwidth]{figures/theory/incident-outgoing-radiance.pdf}
    \caption{The incident radiance function $L_i(\*{x}, \omega)$~describes the distribution of radiance arriving at a~point as a~function of position and direction. The exitant radiance function $L_o(\*{x}, \omega)$~describes the distribution of radiance leaving the point. By convention, $\omega$ points away from the surface. (\emph{Image adapted from Pharr et al.~\cite{PJH2016}.})}
    \label{fig:incident-outgoing-radiance}
\end{figure}

% TODO cite PBRT3:ch5, RTR4:ch8

%\subsection{Photometry}

While radiometry deals with purely physical quantities, a~related field, called \emph{photometry}, takes into account the sensitivity of the human eye. It is the study of visible electromagnetic radiation in relation to its perception by the human visual system. For each radiometric quantity there is a~corresponding photometric quantity. To compute a~photometeric quantity, a~spectral radiometric quantity is integrated against the spectral response curve, which describes the relative sensitivity of the human eye to various wavelengths.

\emph{Luminance} is often used as a~measure of how bright flat surfaces appear to a~human eye. For example, street lighting has a~luminance around 10\,nits, high dynamic range (HDR) screens typically range from 500 to 1000\.nits, while clear sky has~a luminance of about 8000\,nits. Table~\ref{tab:quantities} provides a~summary of radiometric quantities and their photometric equivalents.

% TODO cite PBRT3:ch5, RTR4:ch8

\begin{table}[H]
\begin{tabular}{|lll|lll|}
\hline
\textbf{Radiometric}    &        & \textbf{Unit}         & \textbf{Photometric}        &  & \textbf{Unit}                       \\ \hline
Radiant energy & $Q$    & joule (J)    & Luminous energy    &  & talbot (T)                 \\ \hline
Radiant flux   & $\Phi$ & watts (W)    & Luminous flux      &  & lumen (lm)                 \\ \hline
Intensity      & $I$    & W/sr         & Luminous intensity &  & lm/sr = candela (cd)       \\ \hline
Irradiance     & $E$    & W/$m^2$      & Illuminance        &  & lm/$m^2$ = lux (lx)          \\ \hline
Radiance       & $L$    & W/($m^2$ sr) & Luminance          &  & lm/($m^2$ sr) = cd/$m^2$ = nit \\ \hline
\end{tabular}
\caption{Radiometric quantities and their corresponding photometric quantities.}
\label{tab:quantities}
\end{table}


% ================================================================================

% ================================================================================
% ================================================================================
% ================================================================================
\section{Interaction of Light with Media}
\label{sec:light-media}
%\section{Participating Media}
% Intro - from physics
Computer graphics, for convenience reasons, often assume that scenes are composed of solid objects separated by vacuum, however, this assumption is inaccurate in many real-world situations: fog, smoke, and earth's atmosphere attenuate and scatter light from particles that are distributed throughout the space. Volumes or regions in 3D space with a~large amount of very small particles inside that participate in the overall light transport are called  \emph{participating media} (or simply ``media''). And it is these interactions that give media and physical object the appearance we observe. Figure~\ref{fig:media-examples} shows examples of such media in the real-world.

% TODO Figure: fig:media-examples Examples of various types of participating media
% sky, underwater, fog

\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/brno-sunset.jpg}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/UNDERWATER-LANDSCAPE.jpg}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/forest.jpg}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/brno-mraky.jpg}
\end{subfigure}
\caption{Examples of various types of participating media. Credit: second image from the left by Roban Kramer, flicker.com; the third image from the left by Ferdinand Kozeluh, pixabay.com.}
\label{fig:media-examples}
\end{figure}

% Overview
To goal of this section is to arrive at volume rendering equation, an equation that is used to render participating media. For that, first, the properties of participating media are introduced, then the radiative transfer equation (RTE) is formulated. The last subsections deal with the terms used in VRE and with the various models of scattering.

This section is based on the information found in literature that studies volume rendering in production: Fong et al.~\cite{FWKH2017}, Chapter 11 from Pharr et al.~\cite{PJH2023:ch11}.


%that are based on the assumption that instead of simulating each individual interaction with particle
%makes certain assumptions
%    that allow for treatment
%    probabilisitc process
% TODO cite sources

%
%Scattering in volumes?
\subsection{Properties of Participating Media}

% what are volumes, for our purposes
% volume vs solid
    % other assumptions

% TODO move it up
Volumes or media are a~collection of particles, ranging from atoms and molecules to any particle size. To simulate light transport in these media, interactions between the particles and photons of light are modeled as probabilistic processes. Under certain assumptions the collisions of the photons with the particles can be viewed as statistically independent. 
%The average density of particles must be low enough for the size of the particles to be negligible compared to the average distance between them and the positions of the particles are taken as independent of each other.

% how do we describe volumes:
%   coefficient
When photons travel through such a~volume, they collide with the particles that make up the volume. The results of these collisions define the radiance distribution throughout the volume. Instead of simulating each collision separately, the collisions are modeled as~probability fields. The chance that a~photon collides with a~particle at location~$\*{x}$~is defined by a~coefficient $\sigma(\*{x})$, the probability density of collision per unit distance traveled in the medium. Its units are reciprocal distance ($m^{-1}$) and it is usually a~spectrally varying quantity.

% start with radiance, use from Section of radiance
There are three main physical processes happening when photon collides with a~particle (illustrated in Figure~\ref{fig:processes}) that affect the distribution of radiance inside the volume:
% physical processes
\begin{itemize}
    \item \emph{Absorption}: describes the collision when a~photon is absorbed by the volume, expressed with the absorption coefficient $\absorbCoef(\*{x})$. This causes a~reduction in the amount of radiance along a~ray. Physically, the energy of the photon is converted to another form of energy, such as heat.
    \item \emph{Scattering}: describes when a~photon as a~result of a~collision is scattered into a~different direction defined by the \emph{phase function}. The radiance along a~ray can either decrease due to the photons being deflected out of the ray, also called \emph{out-scattering}, or it can increase when photons are scattered into the path of the ray, called \emph{in-scattering}. Both of these effects are expressed using the scattering coefficient $\scatterCoef(\*{x})$.
    \item \emph{Emission}: causes an increase in radiance due to chemical, thermal, or nuclear processes that convert energy of the photon into visible light. The change in radiance is also determined by how efficiently is light absorbed and emitted by the volume.
\end{itemize}

% TODO above the description of physical processes
% figure wih images of each of the processes 
\begin{figure}[H]
\centering
\begin{subfigure}{1.05\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/processes.pdf}
\end{subfigure}
\caption{Illustration of the different events that can happen when light ray propagating in direction $\omega$~collides with a~particle at location $\*{x}$~resulting in a~change in radiance. 
}
\label{fig:processes}
\end{figure}

% TODO WHERE TO PUT THIS

The properties of a~volume describing how light interacts with the volume are defined by absorption and scattering coefficients, together with the phase function and emission. These properties can be constant, then the medium is called \emph{homogeneous} medium, or they can vary throughout space, then the medium is \emph{inhomogeneous}. The described processes may also behave differently at different wavelengths of light. 
% \emph{homogeneous} or \emph{heterogenous} (or \emph{inhomogeneous})

As a~collective term for the total reduction in radiance due to absorption and out-scattering, so-called \emph{attenuation} (or \emph{extinction}, also informally referred to as the density), is introduced. Physically, it can be interpreted as a~collision where both scattering and absorption are involved. It is defined as sum of the absorption and scattering coefficients:
\begin{align}
    \extinctCoef  = \absorbCoef + \scatterCoef .
    \label{eq:extinction}
\end{align}

Related to the extinction coefficient is the \emph{single-scattering albedo} or simply \emph{medium's albedo} $\rho$, defined as:
\begin{align}
    \rho = \frac{\scatterCoef}{\scatterCoef + \absorbCoef} = \frac{\scatterCoef}{\extinctCoef} .
    \label{eq:single-scat-albedo}
\end{align}
It is a~measure of the medium's reflectivity. It determines the amount of scattered radiance in the medium. In another words, it is the ratio or probability of scattering versus absorption at a~scattering event. The value of $\rho$~is in range [0, 1]: a~value close to $0$~indicates that most of the light is absorbed, resulting in a~murky medium, such as dark exhaust smoke; a~value close to $1$~indicates that instead of the light being absorbed, it is mostly scattered, resulting in a~brighter medium, such as clouds (as shown in the rightmost image in Figure~\ref{fig:media-examples}), or the earth's atmosphere. Water is another example of a~medium that is highly absorptive, as shown in the second image from the left in Figure~\ref{fig:media-examples}. Figure~\ref{fig:increasing-coefs} illustrates the effects of absorption and scattering on color of a~liquid inside a~container. Highly scattering and low absorptive liquids exhibit milky appearance, while highly absorptive liquids with low scattering are very murky. The combination of both high absorption and scattering gives a~muddy appearance.

% - used to describe volume , while the albedo usualy constant for a~single type of volume
% while the extiction changes - driven by a~density field represented by voxels or a~procedural function

% TODO photos and examples of homogenou vs heterogenous
% TODO examples of media
\begin{figure}[H]
    \centering
\begin{subfigure}{.8\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/RTR4-increasing-coefs.png}
\end{subfigure}
\caption{Photos of containers with liquid with varying absorption and scattering. (\emph{Image from Real-time Rendering, 4\,th edition~\cite{RTR2018}.})}
\label{fig:increasing-coefs}
\end{figure}

\subsection{Light Propagation in Participating Media}

%First, each  - first RTE - necessary terms are explained to arrive at VRE

As mentioned in the previous section, there are three events that affect the radiance along a~light ray propagating through a~medium: absorption, scattering, and emission. Then the change of radiance along the light ray traveling in direction $\omega$ from point $\*{x}$ due to these events is given by the \emph{radiative transfer equation} (RTE)~\cite{Cha1960}:
\begin{align}
    \frac{\text{d}L(\*{x},\omega)}{\text{d}\*{x}} &= 
        \underbrace{
            \absorbCoef(\*{x}) L_e(\*{x},\omega)
        }_\text{emission}
        +
        \underbrace{
            \sigma_s(\*{x}) L_i(\*{x},\omega)
        }_\text{in-scattering}
        -
        \underbrace{
            \sigma_a(\*{x}) L(\*{x},\omega)
        }_\text{absorption}
        -
        \underbrace{
            \sigma_s(\*{x})L(\*{x},\omega)
        }_\text{out-scattering} \\
        &= 
        \underbrace{
            L_v (\*{x}, \omega)
        }_\text{contribution}
        -
        \underbrace{
            \sigma_t(\*{x}) L(\*{x}, \omega)
        }_\text{extinction}
    \label{eq:RTE}
\end{align}
The RTE defines the distribution of radiance in media. In other words, the derivative in radiance along the ray as the amount of radiance lost due to absorption, or out-scattering, denoted by the extinction term, and the amount of radiance gained, $L_v$, due to in-scattering (light coming from any direction scattered towards $\omega$) and emission~\cite{FWKH2017}.

The RTE provides an intuitive way to understand light transport in media. In order to compute the medium contribution for a~ray of given length the RTE needs to be integrated into the \emph{volume rendering equation} (VRE). The integration of RTE to arrive at VRE is studied, for instance, in Fong et al.~\cite{FWKH2017}.

Then the general VRE can be defined in the following form:
\begin{align}
    L(\*{x}, \omega) &=
    \underbrace{
        T(\*{x},\*{x}_s) L_s(\*{x}_s, \omega)
    }_\text{surface contribution}
    +
    \underbrace{
        \int_{t=0}^{S} T(\*{x},\*{x}_t) L_v(\*{x}_t, \omega) \text{ d}t
    }_\text{volume contribution}, \\
    \label{eq:vre}
    L_v(\*{x},\omega) &=
        \absorbCoef(\*{x}) L_e(\*{x}, \omega) +
        \scatterCoef(\*{x}) L_i(\*{x}, \omega) .
\end{align}
It is the measurement of the total radiance $L(\*{x},\omega)$ along a~ray propagating in the direction $\omega$ starting from $\*{x}$. This is illustrated in Figure~\ref{fig:vre}. The integrals are parameterized along the negative part of the ray, with the increasing ray parameter $t$~ending at a~boundary (surface or other volume) with the ray parameter $S$. The surface contribution term expresses the contribution due to the incident radiance at the end of the ray $L_s$ at point $\*{x}_s = \*{x} - S\omega$ that is attenuated by the \emph{transmittance} function $T$ (described in Section~\ref{sec:trans}). The integral on the right side accumulates all the emitted or in-scattered radiance (described in Section~\ref{sec:inscat}) by the volume into $\omega$, attenuated by $T$ at the parameterized location $\*{x}_t = \*{x} - t\omega$.

% TODO fig:vre
\begin{figure}[H]
    \centering
\begin{subfigure}{.9\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/vre.pdf}
\end{subfigure}
\caption{Illustration of the volume rendering equation. A~ray is propagating along a~negative direction $\omega$ starting at point $\*{x}$ with parameter $t$ along the points $\*{x}_t = \*{x} - t\omega$ ending at a~surface point $\*{x}_s = \*{x}-S\omega$. The radiance at point $\*{x}$ is the result of the incident radiance $L_s$ at $\*{x}_s$, the emitted $L_e$ and the in-scattered radiance $L_i$ all attenuated by the medium, depicted as a~cloud-like shape.}
\label{fig:vre}
\end{figure}

The VRE describes the radiance in terms of where it comes from, rather than where it is traveling to, as expressed by the RTE. The VRE ``looks'' backwards accumulating all of the contributions with integrals while the direction $\omega$ points in the direction of radiance flow.

The incident radiance at a~surface point in direction $\omega$~can be expressed as
\begin{align}
    L_s(\*{x},\omega) = \absorbCoef(\*{x}) L_e(\*{x}, \omega) +
                        \int_{\Omega^{+}} f_r(\*{x},\omega,\omega_i) L(\*{x},\omega_i) (\*{n} \cdot \omega_i) \text{ d}\omega_i ,
\end{align}
where $f_r$~is a~hemispherical function called the \emph{bidirectional reflectance distribution function} (BRDF), a~surface equivalent to phase function, and $\*{n}$~is the surface normal. The surface contribution is gathered from the upper hemisphere of directions $\Omega^{+}$, in respect to $\*{n}$. In the absence of any participating medium, i.e., $\absorbCoef = \scatterCoef = 0$, and $L_e = 0$, the volume contributions are zero hence giving the standard rendering equation for surfaces (Kajiya, 1986~\cite{Kaj1986}). Section~\ref{sec:brdf} describes the BRDF in more detail. Transmittance function $T$ and the terms related to scattering, such as in-scattering $L_i$ and phase function are further described in the sections below.

%Since the focus of this thesis is on volume rendering, this part will not be discussed any further (refer the reader to Pharr et al.~\cite{PJ2016} instead). On the other hand, the transmittance function $T$ and the terms related to scattering, such as in-scattering $L_i$ and phase function, are relevant and require further inspection in the sections below.

\subsection{Transmittance}
\label{sec:trans}
The extinction term in the RTE (in Equation~\ref{eq:RTE}) describes the amount of radiance lost along a~ray due to out-scattering or absorption,
\begin{align}
    \frac{\text{d}L(\*{x},\omega)}{\text{d}\*{x}} = -\extinctCoef(\*{x}) L(\*{x}, \omega) ,
\end{align}
can be computed using the \emph{transmittance} $T$, which gives the fraction of radiance that is transmitted (able to get through) between two points in~a medium:
\begin{align}
    T(\*{x}_a, \*{x}_b) = e^{-\tau}
    \text{ , where }
    \tau = \int_{\*{x}=\*{x}_a}^{\*{x}_b} \extinctCoef(\*{x}) || \text{d}\*{x} || .
    \label{eq:trans}
\end{align}
This relationship is known as the \emph{Beer-Lambert law}. The \emph{optical thickness} $\tau$ (or sometimes \emph{optical depth}) is unitless and represents the amount of light attenuation. For a~constant $\extinctCoef$ in a~homogeneous medium, the integral can be evaluated, giving simply \emph{Beer's law}:
\begin{align}
    T(\*{x}_a, \*{x}_b) = e^{-\extinctCoef ||\*{x}_a - \*{x}_b ||} .
    \label{eq:beers-law}
\end{align}
Transmittance has certain observable properties. It is always between 0\,and\,1. Transmittance from and to the same point is 1, $T(\*{x},\*{x}) = 1$, and in vacuum $\extinctCoef = 0$ and so $T(\*{x},\*{x}') = 1$ for all $\*{x}'$. If the attenuation coefficient satisfies the directional property $\extinctCoef(\omega) = \extinctCoef(-\omega)$, i.e., it does not vary with direction but can vary with position, then the transmittance between the same points, $T(\*{x}_a,\*{x}_b) = T(\*{x}_b,\*{x}_a)$, is the same in both directions. Furthermore, transmittance is multiplicative along points on a~ray in any media: 
\begin{align}
     T(\*{x}_a,\*{x}_c) = T(\*{x}_a,\*{x}_b) T(\*{x}_b,\*{x}_c) 
\end{align}
for all points $\*{x}_b$ between $\*{x}_a$~and $\*{x}_c$. This is a~useful property for implementations, since it allows for computation of transmittance between two bounding points along a~ray by taking the product of transmittance at multiple points in between the bounding points along the same ray. Figure~\ref{fig:trans} shows the transmittance function at various depths for different values of the extinction coefficient.

% TODO graph image and example
\begin{figure}[H]
    \centering
\begin{subfigure}{.6\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/trans.pdf}
\end{subfigure}
\caption{Transmittance function as a~function of depth and extinction $\extinctCoef = 0.5$ (blue), $\extinctCoef = 1.0$ (red), and $\extinctCoef = 2.0$ (green). The higher the extinction, the less light travels through medium.}
\label{fig:trans}
\end{figure}


\subsection{Scattering and Phase Functions}
\label{sec:inscat}

Both the radiative transfer equation (RTE) and volume rendering equation (VRE) (in Equations~\ref{eq:RTE},~\ref{eq:vre}) contain the out-scattering and in-scattering terms. While out-scattering accounts for reduction in radiance along a~ray due to scattering in different directions, in-scattering accounts for an increase in radiance due to scattering from other directions (as illustrated in Figure~\ref{fig:vre}):
\begin{align}
    L_i (\*{x}, \omega) = \scatterCoef(\*{x}) \int_{\Omega} p(\*{x}, \omega, \omega_i) L(\*{x}, \omega_i) \text{ d}\omega_i  ,
    \label{eq:inscat}
\end{align}
where $\Omega$ denotes the spherical domain around the position $\*{x}$ and $p$~is the phase function. The $\scatterCoef$ in front of the integral models the incoming radiance of all the directions that scattered into the direction of the ray. 

Phase function models the scattering of light at a~macro level.
%Under the assumption that particles are far enough from each other so that the inter-particle interactions can be ignored, the scattering of light is modeled at a~macro level using phase functions. 
It describes the probability and angular distribution of scattering directions. In other words, it is the probability that a~light's forward direction of travel will scatter in a~given direction. Figure~\ref{fig:phase-fnc} illustrates how the light, traveling in forward direction, is scattered at point~$x$ based on a~distribution of scattering directions characterized by the phase function $p(\theta)$~(in red) with two lobes\,---\,a~smaller backward-scattering and larger forward-scattering lobe. The light is scattered backwards to camera\,A at angle, $\theta_A$, between the light's direction (orange arrow) and the direction to the camera. Analogously, the light is scattered also to camera\,B, but in greater amounts due the larger forward-scattering lobe.

\begin{figure}[H]
    \centering
\begin{subfigure}{.9\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/phaseFnc.pdf}
\end{subfigure}
\caption{Light is scattered according to a~phase function (in red). Camera\,A receives a~small portion of light due to backward-scattering lobe, while camera\,B receives a~larger amount due to stronger forward-scattering lobe.}
\label{fig:phase-fnc}
\end{figure}

For a~phase function to be a~probability distribution it must be energy-conserving and -preserving. A~phase function $p$ between two directions $\omega$~and $\omega'$~is normalized over the sphere domain $\Omega$:
    \begin{align}
        \int_{\Omega} p(\omega, \omega') \text{ d}\omega' = 1 .
    \end{align}
Phase functions have also the property called \emph{reciprocity}, which means, that the interchange of the two directions $\omega$~and $\omega'$~results in the same value of the phase function.


% many phase functions have been developed:
%   - parameterized models, can be used to fit a~function with a~small number of parameters to measured data
%   - analytic models - derived from the scattered radiance distribution from particles with known shape and material (e.g., spherical water droplets)

There are two types of phase functions based on the modeled angular distributions: isotropic and anisotropic. \emph{Isotropic} phase functions have an equal probability of scattering incoming light into any direction independent of the two input directions. The \emph{isotropic phase function} is one of the simplest phase functions:
\begin{align}
    p(\theta) = \frac{1}{4\pi} .
\end{align}

\emph{Anisotropic} phase functions can exhibit complicated angular distribution. These are typically modeled accurately by \emph{Mie scattering} or approximated by \emph{Rayleigh scattering}. Depending on the relative size of a~particle $s_p = \frac{2\pi r}{\lambda}$, with radius $r$ and the considered wavelength of light $\lambda$, different types of scattering are identified:
\begin{itemize}
    \item $s_p \ll 1$: Rayleigh scattering (e.g. air),
    \item $s_p \approx 1$: Mie scattering,
    \item $s_p \gg 1$: Geometric scattering.
\end{itemize}
There are also other components that can influence the scattering result, such as the index of refraction, participating media content, etc.

\vspace{1em}

\noindent\textbf{\emph{Rayleigh Scattering}} was derived to model scattering of light from molecules in the air and describes scattering by spherical particles that are much smaller than the wavelength of the light. It is usually used to describe light scattering in the earth's atmosphere. The function (in Equation~\ref{eq:rayleigh-phase}) has two lobes, as shown in Figure~\ref{fig:rayleigh-phase}, a~backward and forward scattering one, relative to the light direction. 
\begin{align}
    p(\theta) = \frac{3}{16\pi}(1 + \cos^2 \theta) .
    \label{eq:rayleigh-phase}
\end{align}
Rayleigh scattering is also highly wavelength-dependent. Scattering coefficient of Rayleigh scattering shows that short wavelengths, such as blue or violet light, are scattered much more than long wavelengths, such as red light~\cite{Hil2016}.

\begin{figure}[H]
\centering
\begin{subfigure}{.47\textwidth}
    \centering
    \includegraphics[width=.75\linewidth]{figures/theory/Rayleigh.pdf}
    \caption{Polar plot of the Rayleigh phase function as a~function of angle $\theta$~between the incident light and out-scattering directions. It is a~two-lobed function with the same chance of forward and backward scattering.}
    \label{fig:rayleigh-phase}
\end{subfigure}
\begin{subfigure}{.47\textwidth}
    \centering
    \includegraphics[width=.75\linewidth]{figures/theory/HG.pdf}
    \caption{Polar plot of the Henyey-Greenstein phase as a~function of angle $\theta$. The parameter $g$~increases from 0~(green) to 0.3~(orange) and 0.6~(blue), resulting in a~strong forward scattering lobe to the right.}
    \label{fig:HG-plot}
\end{subfigure}
\caption{Polar plots of Rayleigh (on the left) and Henyey-Greenstein (on the right) phase functions, representing Rayleigh scattering and Mie scattering, respectively.}
\end{figure}

% Mie Theory is the solution of Maxwell's Equations for a~plane wave incident onto a~sphere with the same size or larger than the wavelength of light.
\noindent\textbf{\emph{Mie Scattering}} is used to model scattering by particles that are about the same size as the light's wavelength. Phase function for a~specific particle modeled by Mie scattering has typically complex distribution with strong and sharp directional lobes making them expensive to compute. Instead, such phase functions are most commonly approximated by \emph{Henyey-Greenstein} (HG) phase function~\cite{HG1940}. HG was designed to be a~good fit to measured scattering data that exhibit strong backward or forward scattering. For instance, smoke, fog, or dust-like participating media. The HG phase function is evaluated using
\begin{align}
    p_{HG} (\theta, g) = \frac{1}{4\pi} \frac{1-g^2}{ (1 + g^2 - 2g \cos \theta)^{1.5}} ,
    \label{eq:phase-HG}
\end{align}
where the parameter $-1 \leq g \leq 1$~is used to control the asymmetry of the phase function or the distribution of the scattered light in the following directions: backward ($g<0$), isotropic ($g=0$), or forward ($g > 0$). Figure~\ref{fig:HG-plot} shows a~polar plot of the HG with the parameter $g$~increasing from 0~(green) to 0.3~(orange) and 0.6~(blue). To represent more complex phase functions, for example, with forward and backward scattering lobes at the same time, it is possible to blend between multiple HG~\cite{Hil2016}.

\vspace{1em}

\noindent\textbf{\emph{Geometric Scattering}} occurs for very large particles and it is usually simulated on a~macro level using complex scattering phase functions. Rainbow effect is an example of this type of scattering. It is caused by internal reflection of light inside water particles in the air, dispersing the sun light into a~visible spectrum over a~small visual angle ($\approx 3$\,degrees) of the resulting backward scattering~\cite{Hil2016}.


% Wronski: Depending on the medium particles, amount of light that takes part in those processes can be different. One example of scattering models is Rayleigh scattering. It is scattering of very small particles (like air particles), responsible for blue color of the sky. It is very isotropic, uniform, but wavelength dependent – scattering is stronger for shorter wavelengths and absorption is negligible. On the other hand, so called Mie scattering of bigger particles (like aerosols or dust) has very anisotropic shape with strong forward lobe and much higher absorption proportion.

% Phase function is a function that describes how much light is scattered in all direction. It is a function of angle between the light vector and outgoing direction vectors. It has property of energy conservation, so integral over all directions must be equal to 1 (or under 1 if it contains information about absorption baked in it). Some phase functions can be very complex and described from various models or real captured data



% -============================================================================
\section{Interaction of Light with Surfaces Bounding Water Body}
\label{sec:light-surface}
% TODO rename better so it is the closest to what is happening here
% sections based on RTR4:ch8 and https://oceanopticsbook.info/view/surfaces/the-level-sea-surface

This section develops the equations that describe how light is scattered by surfaces that bound the water body, such as water surface or an opaque sea floor. Section~\ref{sec:air-water-interface} discusses the basics of reflection and refraction of light by a~level water surface. Section~\ref{sec:fresnel} then describes the equations that govern the amount of light that is transmitted and reflected. Lastly, Section~\ref{sec:brdf} introduces the bidirectional reflectance distribution function, the fundamental quantity describing how an opaque surface reflects light. Most of the information in this section is based on Chapter \emph{The Level Sea Surface} from Ocean Optics Web Book~\cite{OOWB2021-level}, and on Chapter\,9 from Real-Time Rendering, 4th edition~\cite{RTR2018:ch9}. The rest is cited where appropriate.

The concepts described below are modeled from the point of view of geometric optics (i.e., light is abstracted in form of rays). Any wavelength-dependency of these concepts is taken as implicit. Geometric optics provides sufficient level of abstraction that allows a~real-world wind-blown sea surface to be modeled as a~collection of locally flat patches that can be arbitrarily tilted. 


%- what needs to be defined - organization
    %- how it all fits together

\subsection{Interaction of Light with Planar Air-Water Interface}
\label{sec:air-water-interface}

From an optical perspective, an object surface is represented by a~two-dimensional interface separating two volumes. When a~ray of light strikes a~surface, two aspects of the surface determine the resulting behavior: the substances on either side and the geometry of the surface. In this section an air-water interface, i.e. a~level water surface separating air and water, is considered as an example of a~perfectly plat plane (as illustrated in Figure~\ref{fig:air-water}).

The media on each side are described by an optical property, called \emph{index of refraction} (IOR). The IOR is an abstraction of the interactions of light at molecular level. It defines how the medium affects light of a~given wavelength. In its full form, it is defined as a~complex number $n + i\kappa$, where $\kappa$ is the \emph{attenuation index} and it refers to the absorption of the medium, and the real part of the IOR, denoted by $n$,~determines the direction of the refracted ray. Only the real part of IOR will be considered in any further relations. For an interface separating two media, the incident light can originate on either side of interface. The difference in values of the indices of refraction causes a~discontinuity in density, which causes the incident light to scatter differently based on the different wavelengths that make up the light (this wavelength-dependency is taken as implicit and is not explicitly represented).

Figure~\ref{fig:air-water} illustrates a~planar surface with normal $\*{n}$ separating two media: air, and water, denoted by indices of refraction $n_a$, and $n_w$, respectively. A~ray of light, originating on the air side, strikes the surface from direction $\*{i}$ causing the ray to reflect and refract in directions $\*{r}$ and $\*{t}$, respectively. Because the directions of the incident, reflected and refracted rays lie in the same plane, the amount of incident light is divided into its reflected and refracted directions. Figure~\ref{fig:water-air} illustrates the analogous case but now the light is incident from the water side.

\begin{figure}[H]
\centering
    \raisebox{2.23em}{
\begin{subfigure}{.355\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/air-incident.pdf}
    \caption{Light incident on air-water interface is reflected and refracted.}
    \label{fig:air-water}
\end{subfigure}
    }
\begin{subfigure}{.535\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/water-incident.pdf}
    \caption{Light incident on water-air interface (orange) is reflected and refracted. The red (or darker orange) direction shows total internal reflection when the incident angle $\theta_i$ is greater than the critical angle $\theta_c$.}
    \label{fig:water-air}
\end{subfigure}
\caption{Illustrates a~planar interface between air and water with indices of refraction $n_a$ and $n_w$, respectively. The left image shows the light incident from the air side, while the right image shows light incident from the water side. Incident, reflected, and refracted directions are denoted $\*{i}, \*{r}$, and $\*{t}$, respectively, and $\*{n}$ is the surface normal. 
}
\label{fig:air-water-interface}
\end{figure}

The angle between the incident direction and the normal is called the \emph{incident angle}: 
\begin{align}
    \theta_i = \cos^{-1}(|\*{i}\cdot \*{n}|) .
\end{align}
Based on the \emph{Law of Reflection}, the reflected angle is always equal to the incident angle: $\theta_r = \theta_i$. 

The incident and transmitted angles are related by \emph{Snell's law}:
\begin{align}
    n_1 \sin \theta_1 = n_2 \sin \theta_2 ,
    \label{eq:snells}
\end{align}
where subscripts 1 and 2 refer to the media on the side where the incident light originates, and on the side the light is transmitted to, respectively.

For air-incident light, as illustrated in Figure~\ref{fig:air-water}, the IOR of air is: $n_1 = n_a = 1$ and the IOR of water is $n_2 = n_w$ with $n_a < n_w$. According to Snell's law, now simplified as $\sin \theta_i = n_w \sin \theta_t$, then the angle of transmission $\theta_t$ is given by:
\begin{align}
    \theta_t = \sin^{-1} \left(\frac{1}{n_w} \sin \theta_i \right) .
\end{align}
The reflected and refracted direction are given by the following equations:
\begin{align}
    \*{r} &= \*{i} - 2 (\*{i} \cdot \*{n}) \*{n} , \\
    \*{t} &= \frac{1}{n_w} (\*{i} - c\*{n}) ,
\end{align}
where $c = \*{i}\cdot \*{n} + \sqrt{ (\*{i}\cdot \*{n})^2 + n_w^2 -1 }$.

For the water-incident case, as illustrated in Figure~\ref{fig:water-air}, Snell's law reads $n_w \sin \theta_i = \sin \theta_t$, which gives the angle of transmission as:
\begin{align}
    \theta_t = \sin^{-1} (n_w \sin \theta_i) .
    \label{eq:water-air-trans-angle}
\end{align}
The relations between the reflected and refracted directions are then:
\begin{align}
    \*{r} = \*{i} - 2 (\*{i}\cdot \*{n}) \*{n} , \\
    \*{t} = n_w \*{i} - c\*{n} ,
\end{align}
where $c = n_w\*{i}\cdot \*{n} - \sqrt{(n_w\*{i}\cdot \*{n})^2 - n_w^2 + 1}$.

In the case where $n_1$ < $n_2$, as seen in the air-water interface in Figure~\ref{fig:air-water}, the incident light originates on the side with the lower IOR, this results in a~reflection that is called \emph{external reflection}. In the other case, where the light originates on the denser water side (as shown in Figure~\ref{fig:water-air}), when $n_1 > n_2$ this case of reflection is called \emph{internal reflection}. 

In most cases, when internal reflection occurs the light is transmitted from the more dense medium to the less dense. However, for certain values of the incident angle there is no refraction and all light incident onto the water side is reflected back into the water. This is called \emph{total internal reflection}. This phenomenon occurs for values of incident angle $\theta_i$ that are greater than the \emph{critical angle for total internal reflection}, 
\begin{align}
    \theta_c = \sin^{-1} \left(\frac{1}{n_w} \right) ,
    \label{eq:crit-angle}
\end{align}
in these cases there is no real solution for the inverse sine in Equation~\ref{eq:water-air-trans-angle}. The right part of Figure~\ref{fig:water-air} illustrates this with the dotted line representing the critical angle.

For the air-water case, when $\theta_i = 0^{\circ}$, the incident light is perpendicular to the surface, $\*{i} = \*{n}$, this is called \emph{normal incidence}. When light strikes the surface at grazing incident angle of $\theta_i = 90^{\circ}$ then it is transmitted into the water at the critical angle: $\theta_t = \theta_c$. This causes the light from the entire sky to be transmitted through the surface into a~cone of half angle $\theta_c$, which is known as \emph{Snell's cone} or \emph{Snell's window} (shown in Figure~\ref{fig:snellswindow}). 
% TODO
%Real-world photos of these two phenomena are shown in Figure~\ref{fig:water-phenomena}.

\begin{figure}[H]
    \centering
\begin{subfigure}{.6\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/snellswindow.png}
\end{subfigure}
\caption{Illustrates the incident light being refracted into a~cone of directions causing a~phenomenon called \emph{Snell's window}. (\emph{Image from Cheng et al.~\cite{CCZTG2020}.})
}
\label{fig:snellswindow}
\end{figure}

%Figure two images: total internal reflection and Snel's cone or snells window
% fig:water-phenomena

\subsection{Fresnel's Equations for Unpolarized Light}
\label{sec:fresnel}

Previous section described the relations between the angles and directions of the incident, reflected and transmitted light. However, the amount reflected or transmitted light was not addressed. This information is given by \emph{Fresnel's equations}. 

The same notations and assumptions from the previous section also hold here: a~planar interface, e.g., the flat water surface separating air and water (as shown in Figure~\ref{fig:air-water}, both characterized by indices of refraction $n_a = 1$ and $n_w$, respectively), and a~beam of unpolarized incident light with some irradiance measured on a~surface normal to the direction of propagation. The incident light on the flat surface is split into a~reflected part and a~refracted part. The fraction of this incident irradiance that is reflected by the air-water surface is given by the \emph{Fresnel reflectance $R_F$}:
\begin{align}
    R_F(\theta_i ) = \frac{1}{2} \left\{ 
        \left[ \frac{ \sin(\theta_i - \theta_t) }{ \sin(\theta_i + \theta_t) } \right]^2  
        +
        \left[ \frac{ \tan(\theta_i - \theta_t) }{ \tan(\theta_i + \theta_t) } \right]^2  
    \right\} ,
    \label{eq-fresnel-full}
\end{align}
which is in range 0 to 1, and holds for $\theta_i \neq 0$. For normally incident light, i.e. $\theta_i = 0$, the reflectance is
\begin{align}
    R_{F} (\theta_i = 0) = \left( \frac{n_w - 1}{n_w + 1} \right)^2 .
    \label{eq-fresnel-normal}
\end{align}
The reflectance for normal incidence can be thought of as the characteristic specular color of the substance. Although, the value of $R_F(\theta_i)$ is dependent only on the incident angle, in reality, it varies continuously over the visible spectrum. The general Fresnel equations also use the complex part of IOR, however for water at near-UV to near-IR wavelengths (see Figure~\ref{fig:visible-range} for the whole electromagnetic spectrum) the difference is negligible.

Due to the conservation of energy, the reflected and transmitted energy must equal the incident energy. Therefore the transmitted fraction of the incident light is $R_T = 1 - R_F$. 

%Although, the value of $R_F(\theta_i)$ is dependent only on the incident angle, in reality, it varies continuously over the visible spectrum. The general Fresnel equations use the complex IOR $m = n + i\kappa$, where $\kappa(\lambda) = \lambda \absorbCoef(\lambda) / 2\pi$ is the complex part. Then Equation~\ref{eq-fresnel-normal} is
%\begin{align}
%    R_{F} (\theta_i = 0) = \left| \frac{m - 1}{m + 1} \right|^2
%    = \frac{ (n_w - 1)^2 + \kappa^2 }{ (n_w + 1)^2 + \kappa^2 } .
%    \label{eq-fresnel-normal-full}
%\end{align}

Regarding the air- and water-incident light in the cases described in the previous section, Equations~\ref{eq-fresnel-full} and~\ref{eq-fresnel-normal} hold for both of these cases. For water-incident light that undergoes total internal reflection, i.e., $\theta_i \geq \theta_c$, then $R_F(\theta_i) = 1$. Figure~\ref{fig:fresnel} shows the Fresnel reflectance for the range of indices of refraction of water at visible wavelengths.


\begin{figure}[H]
    \centering
\begin{subfigure}{.7\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/theory/fresnel.pdf}
\end{subfigure}
\caption{Comparison of water-incident (internal) and air-incident (external) reflectance curves at a~air-water interface for various values of the index of refraction of water. The shown value $\theta_i = 47.33$ is the critical angle for $n_w = 1.36$.
\label{fig:fresnel}
}
\end{figure}


\subsection{The BRDF}
\label{sec:brdf}
% https://pbr-book.org/3ed-2018/Color_and_Radiometry/Surface_Reflection#sec:brdf
% https://www.oceanopticsbook.info/view/surfaces/lambertian-brdfs
Generally, surfaces (not just planar interfaces) scatter incident light into various directions and at various wavelengths. This section describes how surfaces reflect the incident light. 
%directional and spectral distribution of the reflected light, i.e., the amount of light reflected, per wavelength, in particular direction.
%- only consider reflection, no ss
%- abstraction 

The \emph{bidirectional reflectance distribution function} (BRDF) is a~function describing reflection from a~surface~\cite{PJH2016}. Figure~\ref{fig:brdf} illustrates the BRDF: it gives the amount of radiance leaving the surface in the outgoing direction $\omega_o$ towards the camera, $L_o(\*{x},\omega_o)$, as a~result of incident radiance along the direction $\omega_i$, $L_i(\*{x},\omega_i)$.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{figures/theory/brdf.pdf}
    \caption{The BRDF is a~function that describes how much incident light along $\omega_i$ is reflected from the surface in the direction $\omega_o$.}
    \label{fig:brdf}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{figures/theory/lambertian-brdf.pdf}
    \caption{An example of a~perfectly diffuse BRDF, so-called \emph{Lambertian}, that reflects the incident light equally in all directions.}
    \label{fig:lambertian}
\end{subfigure}
\caption{The bidirectional reflectance distribution function (BRDF) illustrated on the left image. The right image shows one of the simplest BRDF.}
\end{figure}

Formally, the BRDF is defined as the ratio of the reflected exitant radiance to incident irradiance. The differential irradiance at point $\*{x}$ is the incident radiance on the surface over a~differential cone of directions $\omega_i$:
\begin{align}
    dE(\*{x}, \omega_i) = L_i (\*{x}, \omega_i) (\*{n} \cdot \omega_i) \text{ d}\omega_i .
\end{align}
Then the BRDF of a~surface for a~pair of directions $\omega_i$ and $\omega_o$ is defined as a~fraction of the incident irradiance~\cite{PJH2016}:
\begin{align}
    f_r(\*{x}, \omega_o, \omega_i) =
        \frac{\text{d}L_o(\*{x},\omega_o)}{\text{d}E(\*{x},\omega_i)}
        = \frac{\text{d}L_o(\*{x},\omega_o)}{L_i(\*{x},\omega_i) (\*{n}\cdot \omega_i) \text{ d}\omega_i} .
\end{align}
It is the differential amount of radiance that is reflected in direction $\omega_o$ proportional to the differential irradiance at point $\*{x}$.

To compute the radiance $L_o$ leaving the surface at $\*{x}$ into the direction $\omega_o$ due to the reflected radiance from the incident direction $\omega_i$, the BRDF is integrated over the upper hemisphere of directions $\Omega^{+}$~\cite{RTR2018:ch9}:
\begin{align}
    L_o(\*{x},\omega_o) =
        \int_{\Omega^{+}} f_r(\*{x},\omega_o,\omega_i) L_i(\*{x},\omega_i) (\*{n}\cdot\omega_i) \text{ d}\omega_i .
\end{align}
This is called the \emph{reflectance equation}. 

Physically based BRDFs have the following properties:
\begin{itemize}
    \item \emph{Reciprocity}: For all pairs of directions $\omega_i$, $\omega_o$:
        $f_r(\*{x},\omega_i, \omega_o) = f_r(\*{x},\omega_o,\omega_i)$, i.e., interchange of the directions with each other produces the same value.
    \item \emph{Conservation of energy}: a~surface cannot reflect more energy than the incident energy:
    \begin{align}
        \int_{\Omega^{+}} f_r(\*{x},\omega_o,\omega') \cos \theta' \text{d}\omega' \leq 1, \hspace{1em} \forall\omega_o .
    \end{align}
\end{itemize}

% https://pbr-book.org/3ed-2018/Reflection_Models/Lambertian_Reflection

One of the simplest BRDFs is \emph{Lambertian}~\cite{PJH2016}. It is used to model perfectly diffuse surfaces, i.e., surfaces that reflect incident illumination equally in all directions (as illustrated in Figure~\ref{fig:lambertian}). It is often used as a~reasonable approximation to many real-world surfaces. The Lambertian BRDF has a~constant value that when evaluated gives the following result:
\begin{align}
    f_r(\*{x},\omega_o, \omega_i) = \frac{\rho}{\pi} ,
\end{align}
where $\rho$ is the reflectivity of the surface and it varies from zero (completely absorbing, black surface) to one (completely reflecting, white). 
%Lambertian surface reflects light equally into all directions, such surfaces are called (among other names) \emph{diffuse}.

The planar surface, shown in Figure~\ref{fig:air-water}, from the previous section, is the opposite of diffuse surface\,---\,a~perfectly \emph{specular} surface, or mirror\,---\,where the incident direction $\*{i}$ is always reflected about the surface normal $\*{n}$.




% ================================================================================
% ================================================================================
%\subsection{Participating Media}
%\todo{assumptions - what will be discussed - concerned with underwater rendering}
%- first must understand the differences in these assumptions (should be explained earlier)
%- depends on what is water, what do want to model in the water

% TODO about coefs. https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/inherent-optical-properties 

% - General volumetric rendering for intro RTR p. 623 (605)
%     - approaches for representation of spatially varying participating media

% ================================================================================
% ================================================================================
% interactions of light with the Earth’s oceans
% 
% or strictly just Optical properties of Natrural Waters
\chapter{Optical Properties of Natural Waters}
\label{ch:optics}

This chapter describes optical properties of real-world natural waters. Natural waters is collective term for earth's fresh and saline waters, such as ponds, lakes, rivers, but also seas, oceans. However, for the purpose of this thesis, the focus is mainly on seas and oceans (used interchangeably throughout this text), both of which have complex and highly variable optical properties, as shown in Figure~\ref{fig:oceans-example}, and both are extensively studied by the field of ocean optics. The terms described here provide baseline knowledge and guidelines for studying and designing bio-optical models of water in the next chapters. The information presented in this chapter is based on chapters from an online resource called \emph{Ocean Optics Web Book}~\cite{OOWB2020}, a~continually developing community resource for optical oceanography. Optical oceanography is a~vast field spanning across multiple other fields, of which only a~tiny fraction is mentioned in this text. However, the reader is free to consult the material if necessary.

Optical properties of natural waters are introduced in Section~\ref{sec:optics-props} and put into a~larger context via the simplified ``ocean-optics organization chart'' in Figure~\ref{fig:optics-diagram}. Furthermore, the optical properties of natural waters cannot be simply approximated based on the optical properties of pure water. Natural waters, both fresh and saline, contain various mixtures of dissolved and particulate matter. These solutes and particles themselves are optically significant and directly influence the optical properties of waters. Consequently, to be able to accurately reproduce the colors of natural waters, it is also necessary to consider the optical properties of such constituents. Section~\ref{sec:optics-constituents} describes optically important constituents in oceans. 

Constituents in natural waters are highly variable in their composition and concentration, therefore show large temporal and spatial variations. Because of such great variability in the optical properties of natural waters, it is hard to find easily tabulated data or simple models. Optical oceanography utilizes coupling between constituent properties and optical properties and results from optical measurements to deduce information. These measurements cross to other areas, such as biological oceanography, marine photochemistry, underwater visibility, ocean color remote sensing, etc~\cite{OOWB2021-iaop-intro}.

Then the following Sections~\ref{sec:optics-absorb}, and~\ref{sec:optics-scatter} study how these constituents affect absorption and scattering in natural waters.


%polarization is ignored, only unpolarized light is considered
% TODO see https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/introduction pre last paragraph

\section{Optical Properties of Natural Waters}
\label{sec:optics-props}
% https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/introduction
Optical oceanography divides optical properties of waters into two categories~\cite{OOWB2021-iaop-intro}: \emph{Inherent optical properties} (IOPs) and \emph{Apparent optical properties} (AOPs). The IOPs describe the properties of a~medium, they do not depend on the ambient light field in the medium. They are the various measures of the absorption and scattering properties of a~water body. 

The AOPs are properties of a~medium that depend both on the medium (the IOPs) and on directional structure of the ambient light field, and are regular and stable enough to be useful descriptors of the water body. Commonly used AOPs are various reflectances, average cosines, and diffuse attenuation coefficients. 

The IOPs and AOPs are mutually exclusive descriptors of optical properties, i.e., an optical parameter cannot be both an IOP and an AOP. However, there are also other measures of the optical properties of a~water body that do not fit into either of the categories. Radiometric variables, such as the various irradiances, are neither IOPs nor AOPs. 

The radiative transfer equation (RTE) connects the optical properties of the water body and light within the water. Figure~\ref{fig:optics-diagram} shows how the RTE is used to compute radiometric variables based on the IOPs and boundary conditions. Then the radiometric variables can be used to compute the AOPs. The graph shows the most commonly used IOPs, AOPs, radiometric variables, and how they are related. Boundary conditions specify the radiance incident onto a~water body and how it is transmitted through the water surface and reflected by the surface and the ocean bottom. Once the IOPs, along with the boundary conditions are given and the radiance distribution is known, then the RTE can be solved.

\begin{figure}[H]
\centering
\begin{subfigure}{.58\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/optics-diagram.pdf}
\end{subfigure}
\caption{Overview of relevant ocean optical properties and their relationships. The inherent optical properties and the boundary conditions of the oceanic environment are fed into the radiative transfer equation to obtain radiometric variables, that are in turn used to obtain the apparent optical properties. (Adapted from the full diagram by Mobley~\cite{OOWB2021-iaop-intro}).
}
\label{fig:optics-diagram}
\end{figure}

\subsection{Apparent Optical Properties}
\label{sec:iops-aops}
The absorption and scattering properties of a~medium such as sea water, are described by its IOPs. The fundamental IOPs that describe how water absorbs and scatters light are the \emph{absorption coefficient} and the \emph{volume scattering function}. Once these two IOPs are known, then it is also known how the medium interacts with unpolarized light. All other related coefficients can be derived from them~\cite{OOWB2021-iaop-iop}. 

%- another commonly used IOP is the \emph{single-scattering albedo} $\rho$, TODO ref to VR
%    - $\rho$ is near one in waters at wavelengths where the attenuation is due primary to scattering
%    - $\rho$ is near zero, for waters and wavelengths where the attenuation is due primarily to absorption

%\subsection{Useful AOPs}
On the other hand, there are no fundamental AOPs. The most commonly used AOPs are those that are the least influenced by external environmental conditions (such as sun location or sky conditions) and those that provide the most useful information. The most commonly used AOPs are the various \emph{reflectances} and \emph{diffuse attenuation functions} (or ``$K$\,functions'')~\cite{OOWB2021-iaop-iop}.

For example, the \emph{irradiance reflectance} $R$~can be related to absorption and backscatter coefficient and chlorophyll concentrations. Also called \emph{irradiance ratio} is the ratio of upwelling ($E_u$) to downwelling ($E_d$) plane irradiances at certain depths (illustrated in Figure~\ref{fig:reflectance}): $R = E_u / E_d$. It is a~measure of how much of the radiance traveling downwards is reflected upwards in any direction. Measurement data for this AOP are widely available, because both the $E_u$ and $E_d$ can be measured in a~straightforward way~\cite{OOWB2020-iaop-refl}. 

$K$\,functions are often used as an approximation of the absorption coefficient. For example, the \emph{diffuse attenuation coefficient for downwelling plane irradiance} $\diffuseAttenDown$ is used to approximate downwelling irradiance $E_d$ at certain depths as it gets attenuated. In other words, it is a~measure of the irradiance attenuated in the water volume per unit of depth. In these use cases, the $\diffuseAttenDown$ is typically assumed to be independent of depth, which can be a~good approximation only for typical oceanic conditions in homogeneous water, and when far enough below the surface~\cite{OOWB2021-iaop-kfnc}. Figure~\ref{fig:kfunction} shows how $\diffuseAttenDown$ approximates absorption coefficient.

\begin{figure}[H]
\centering
\begin{subfigure}{.47\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/irradiances.pdf}
    \caption{Illustrates the irradiance reflectance $R$~as a~measure of upwelling irradiance to downwelling irradiance at depth $z$.~The image shows how radiance along the light rays, which travel downwards and are reflected upwards, contributes to the irradiance reflectance.
    }
    \label{fig:reflectance}
\end{subfigure}
% TODO remake graph to pdf
\begin{subfigure}{.45\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/KfncAbs.png}
    \caption{$K$\,functions at 10\,m depth for waters with chlorophyll concentration and only elastic scattering. (\emph{The graph is provided as is by Mobley~\cite{OOWB2021-iaop-kfnc}}.)
    }
    \label{fig:kfunction}
\end{subfigure}
%\label{fig:media-examples}
\caption{Left image illustrates the irradiance reflectance $R$, the right image illustrates the diffuse attenuation coefficient $\diffuseAttenDown$.}
\end{figure}

\ifdefined\SHOWEXTENDED
For example, the \emph{irradiance reflectance} can be related to absorption and backscatter coefficient and chlorophyll concentrations. Also called \emph{irradiance ratio} is the ratio of upwelling ($E_u$) to downwelling ($E_d$) plane irradiances. It is a~measure of how much of the radiance traveling in all downward directions is reflected upward into any direction at any depth within a~water column, or in the air just above the sea surface. Measurement data for this AOP are widely available, because both the $E_u$ and $E_d$ can be measured in a~straightforward way~\cite{OOWB2020-iaop-refl}.

K functions are often used as an approximation of the absorption coefficient. For example, the \emph{diffuse attenuation coefficient for downwelling plane irradiance} $\diffuseAttenDown$ is used to approximate downwelling irradiance $E_d$ at certain depths. In these use cases, the $\diffuseAttenDown$ is assumed to be independent of depth, which can be a~good approximation only for typical oceanic conditions in homogeneous water, and when far enough below the surface~\cite{OOWB2021-iaop-kfnc}.
\fi



\ifdefined\SHOWEXTENDED
Candidate AOPs that provide them most useful information about water bodies and which are the least influenced by external environmental conditions such as sun location or sky conditions.

Various reflectances are one of the most commonly used AOPs. The field of \emph{ocean color remote sensing} developed algorithms that relate irradiance reflectance $R$~to quantities such as absorption and backscatter coefficients and chlorophyll concentrations. The \emph{sepctral irradiance reflectance} or \emph{irradiance ratio}, $R(z, \lambda)$, is defined as the ratio of spectral upwelling to downwelling plane irradiances:
\begin{align}
    R(z, \lambda) = \frac{E_u(z,\lambda)}{E_d(z, \lambda)} .
\end{align}
It is a~measure of how much of the radiance traveling in all downward directions is reflected upward into any direction at any depth $z$~within a~water column, or in the air just above the sea surface. Measurement data for irradiance reflectance are widely available, because both the $E_u$ and $E_d$ can be measured in a~straightforward way.
% TODO cite https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/reflectances

\emph{Diffuse attenuation functions}, or commonly ``K functions'', are often used as an approximation of the absorption coefficient. For example, the \emph{diffuse attenuation coefficient for downwelling plane irradiance} $\diffuseAttenDown$ is used to approximate downwelling irradiance $E_d$ at certain depth $z$:
\begin{align}
    E_d(z, \lambda ) = E_d (0, \lambda) e^{-\diffuseAttenDown (\lambda) z} .
    \label{eq:diffuse-atten-down}
\end{align}
In these use cases the K~functions are assumed to be independent of depth, which is good enough approximation for typical oceanic conditions, for which the incident lighting is provided by the Sun and sky, in homogeneous water, and when far enough below the surface. But, in reality the conditions are more complicated~\cite{OOWB}.
% https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/k-functions
\fi


%coupling between constituent properties and optical properties
%- optical measurements can be used to deduce information
%- connections between optical properties and the biological, chemical and geological constituents of natural waters
%- connections between optics and biology, chemistry, and geology
%Optical oceanography utilizes results from the biological, chemical, geological and physical subdisciplines of limnology and oceanography, to 
%to find precise and easily tabulated data or simple models. Optical properties are tightly coupled with constituent properties, that cross to other areas, such as biological oceanography, marine photochemistry, underwater visibility, ocean color remote sensing, sediment load, etc.


% Intro into absoprtion and scattering

%both absorption and scattering are equally important to the prediction and understanding of light propagation.  add figure 1:
% https://www.oceanopticsbook.info/view/scattering/scattering-overview


\section{Optically Important Constituents of the Ocean}
\label{sec:optics-constituents}

% https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/inherent-optical-properties
% TODO maybe move somewhere else
% The IOPs in the ocean depend on the following physical characteristics of the particulate and dissolved substances:
% \begin{itemize}
%     \item Composition: the absorption of materials, that make up the particle or dissolved substance, and the index of refraction of the materials relative to the surrounding water.
%     \item Morphology: differently sized and shaped particles scatter light differently.
%     \item Concentration: the number of particles in a~given volume of water, described by the particle size distribution, or the amount of a~dissolved substance in the water.
% \end{itemize}

% https://www.oceanopticsbook.info/view/optical-constituents-of-the-ocean/introduction-to-optical-constituents-of-the-ocean
Marine particles are extremely complex and vary in their composition, pigmentation, shape, internal structure, and packaging~\cite{OOWB2021-const-intro}. The IOPs in the ocean depend on physical characteristics of particulate and dissolved substances. Their composition determines absorption (based on the material), scattering (based on index of refraction), and different sizes and shapes scatter light differently~\cite{OOWB2021-iaop-iop}. To understand their interaction with light some important idealizations have to be made. It is impractical to study each individual particle (there are more than a billion bacteria in one ml of seawater), therefore particles are grouped based on similar properties.

In order to model the optical properties of particles, that is, to describe these properties using analytical or numerical means, empirical data is needed (size, shape, index of refraction, etc.). These can be either obtained from microscopy or need to be deduced from other measurements. The specific models for the different types of particles are usually based on observations, either in the lab or in the field, and on numerical calculations with inputs obtained from observations (e.g. using an optical model such as Mie theory to fit empirical data)~\cite{OOWB2021-const-intro}.

The in-water constituents that affect aquatic optical properties are traditionally grouped into the following categories: sea water (water and inorganic dissolved materials), colored dissolved organic material (CDOM), phytoplankton, non-phytoplankton organic particles, inorganic particles, bubbles. The constituents are identified based on the way their optical properties are measured and are often grouped by similar optical properties~\cite{OOWB2021-const-intro}. 
%There are also classifications by mean particle size of the different constituents, biogeochemical, biological (living, non-living, single-celled, taxonomy\,--\,species, etc.), and the operational classification that uses filtration to filter captured matter and identify particles based on pore size.

\begin{figure}[H]
\centering
\begin{subfigure}{.46\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/dropParticles.png}
    \caption{Particles in a~drop of sea water through microscope. Largest is phytoplankton cell, rest is mineral material.}
\end{subfigure}
\begin{subfigure}{.46\textwidth}
    \centering
    \includegraphics[width=0.9\linewidth]{figures/ocean-optics/flasksCDOM.png}
    \caption{A~flask of clear water (left) and CDOM (right) from filtering coastal water.}
\end{subfigure}
\caption{Photos of ocean constituents from Mobley~\cite{OOWB2021-const-seas}.}
\label{fig:constituents-photo}
%\caption{.}
\end{figure}

%The most commonly modeled constituents are first introduced, and  optically classified
The most commonly modeled constituents are introduced in the list below:
\begin{itemize}
    \item Particulate matter: collective term for the all organic and inorganic, animate and inanimate particles. Is classified into organic\,---\,all living or once-living particulate matter (such as phytoplankton and detritus)\,---\,and inorganic, all lithogenic and mineralic matter (such as sand and calcium)~\cite{OOWB2021-const-intro}.
    \item Phytoplankton: microscopic, single-celled organisms (form chains or colonies), comprised of more than 10,000 species, that possess various pigments (such as chlorophyll-A) and photosynthesize. Absorb sun light, and causes scattering~\cite{OOWB2021-const-phyt}.
    \item Colored dissolved organic matter (CDOM): also referred to as ``yellow matter'', dissolved organic materials (such as picoplankton, viruses, and colloids) that absorb and scatter light, hence the yellowish color~\cite{OOWB2021-const-cdom}.
    \item Non-algal particles (NAP): particulate matter that is not a~phytoplankton pigment, includes all living and detrital organic matter such as the non-pigemnted portion of phytoplankton cells, detritus, viruses. May also include inorganic mineral particles (e.g., shells, clay, silt, and sand) and most recently also tiny particles of plastic~\cite{OOWB2021-const-intro}.
    \item Detritus: non-living particulate matter, including dead bacteria, phytoplankton, fragments of cells, shells, marine now aggregates, etc. Both absorbs and scatters light~\cite{OOWB2021-const-nap}.
    \item Minerals: also clay, silt, and sand, are fractions of typically lithogenic or mineralic particles of different sizes. Mineral particles enter ocean from river discharge, erosion of coastal cliffs, sediment resuspension, and deposition of atmospheric dust. Both absorb and scatter light~\cite{OOWB2021-const-nap}.
\end{itemize}

% figure phytoplankton 
% figure sizes at https://www.oceanopticsbook.info/view/optical-constituents-of-the-ocean/introduction-to-optical-constituents-of-the-ocean


% TODO Jerlov classification of types of waters??

The physical characteristics of the dissolved and particulate substances in the ocean vary by orders of magnitude, therefore the IOPs also vary. The absorption and scattering coefficients in turbid coastal waters, that contain high concentrations of phytoplankton, mineral particles, and dissolved organic matter, can be four orders of magnitude larger than those of pure water. The volume scattering function can also vary (at a~given angle and wavelength) by orders of magnitude between open ocean and coastal waters. Optical oceanography tries to understand how the various constituents of the ocean determine IOPs~\cite{OOWB2021-iaop-iop}.



\section{Absorption in Natural Waters}
\label{sec:optics-absorb}
% TODO cite https://www.oceanopticsbook.info/view/absorption/absorption-overview
% TODO cite https://www.oceanopticsbook.info/view/absorption/absorption-by-oceanic-constituents
Absorption is a~physical process caused by the interaction of light with matter, in which case the light can disappear, with its radiant energy being converted into other forms, such as heat. Both absorption and scattering affect visibility in water in significant ways~\cite{OOWB2021-abs-intro}. Figure~\ref{fig:ocean-absorption} illustrates how absorption makes a~distant image darker. Wavelength-dependent absorption also alters the color of a~distant image, as shown in the right image of the figure.

% TODO figure 1
\begin{figure}[H]
\centering
\begin{subfigure}{.7\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/spongeAbsorb.jpg}
\end{subfigure}
\caption{Simulated effect of absorption on a~photo of a~red sponge seen up close and from distance of 5\,m in clear water. The image is darker and the red color has been removed, additionally scattering has blurred fine details (\emph{Photos by Mobley~\cite{OOWB2021-abs-intro}.})}
\label{fig:ocean-absorption}
\end{figure}

Absorption is highly spectrally dependent compared to scattering. Concentration and composition of the particulate and dissolved constituents and water itself determine the magnitude and spectral features. Due to the IOPs being conservative, the magnitude of the absorption coefficient varies linearly with the concentration of the absorbing material. In theory, the absorption coefficient can be expressed as the sum of the absorption coefficient of each oceanic component. In reality, it is not possible to measure the absorption properties of each component individually, and thus the individual components are usually grouped into similarly absorbing constituents based upon similarity in their optical properties:
\begin{align}
    \absorbCoef(\lambda) = \absorbCoef^{w} (\lambda) + 
                           \absorbCoef^{phyt} (\lambda) +
                           \absorbCoef^{NAP} (\lambda) +
                           \absorbCoef^{CDOM} (\lambda) ,
\end{align}
where the superscripts $w$, $phyt$, $NAP$, and $CDOM$ denote water, phytoplankton, non-algal particles, and colored dissolved organic matter, respectively~\cite{OOWB2021-abs-intro}.

% TODO plot and bigger
% TODO show? should be?
\ifdefined\SHOWABSROPTIONSPECTRA
Figures~\ref{fig:abs-spectra} show common absorption spectra for pure water, general phytoplankton, generic NAP, and generic CDOM. Absorption by pure water (Fig.~\ref{fig:abs-spec-water}) is weak in the blue and strong in the red regions and varies with temperature and salinity. Phytoplankton absorption (Fig.~\ref{fig:abs-spec-phyto} ) shows high spectral variability with peaks in blue and red regions due to the ubiquitous presence of pigment chlorophyll-A. Absorption by non-algal particle (Fig.~\ref{fig:abs-spec-nap}) is strongest in the blue, decreasing approximately exponentially to the red. Colored dissolved organic matter absorption (Fig.~\ref{fig:abs-spec-cdom}) is similar to that of NAP (due to the similarity in composition), but generally exhibits a~steeper exponential slope. 

% Figures 4 side by side
\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/absorbWater.png}
    \label{fig:abs-spec-water}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/absorbPhyto.png}
    \label{fig:abs-spec-phyto}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/absorbNAP.png}
    \label{fig:abs-spec-nap}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/absorbCDOM.png}
    \label{fig:abs-spec-cdom}
\end{subfigure}
\label{fig:abs-spectra}
\caption{Absorption spectra of constituents, from left to right: pure water, generic phytoplankton, generic non-algal particle, generic colored dissolved organic matter~\cite{OOWB2021-abs-intro}.}
\end{figure}
\fi


For clear open ocean environments with very low concentrations of suspended and dissolved material, the absorption coeﬃcient is dominated by water. Minimum in the blue region of water absorption spectrum therefore causes the blue color of the seawater. For other oceanic  and/or coastal environments with high concentrations of suspended and dissolved material, the absorption coefficient is dominated by the material and, in general, the minimum absorption is in the green region, hence the green hue of such environments~\cite{OOWB2021-abs-intro},~\cite{OOWB2021-abs-const}.

% want figure? figure 6



\section{Scattering in Natural Waters}
\label{sec:optics-scatter}

Scattering denotes a~physical process when light changes its direction on interaction with matter. There are two types of scattering in water: \emph{elastic scattering} and \emph{inelastic scattering}. The former describes only the change of direction, while the latter additionally considers a~change of wavelength of the light. While absorption is generally more studied phenomenon than scattering, both are equally important to the prediction and understanding of light propagation. Scattering through small angles blurs the edges of an image, and can cause other effects on distant objects, such as twinkling or mirages. 

In oceans, the field of ocean-color remote sensing measures light leaving the oceans due to backscattering (at angles greater than 90$^{\circ}$) and multiple forward scattering (at angles of a~few tens of degrees or greater). Scattering diverts the direction of light from a~geometric straight-line, thus making the effect of absorption more pronounced. Any interaction of light with matter can cause scattering\,---\,a~change in direction or wavelength. Light is scattered by particles of oceanic constituents of all kinds of sizes~\cite{OOWB2021-scat-intro}. 
% https://www.oceanopticsbook.info/view/scattering/scattering-overview

\subsection{Scattering by Oceanic Constituents}

Fundamentally, elastic scattering occurs when there is a~change in the real part of the index of refraction (IOR), as described in Section~\ref{sec:air-water-interface}. This happens when light travels from a~region with one IOR into a~region with a diﬀerent IOR. The term \emph{surface scattering} refers to scattering caused by the change in IOR at the boundary between two media, such as at the air-water interface. On the other hand, the term \emph{volume scattering} refers to scattering caused by a~change in IOR by the presence of discrete particles in the volume, by thermal fluctuations in density, or by turbulent mixing of ﬂuids with diﬀerent physical properties. Often, scattering by particles and by pure water is modeled.

Scattering by homogeneous spherical particles of some IOR within a different IOR is typically modeled using \emph{Mie Theory} (introduced in~Section~\ref{sec:inscat}). Mie Theory is also being applied to approximate scattering by non-spherical and/or non-homogeneous particles. Such particles are common in nature: ice crystals in clouds, atmospheric dust, resuspended sediments in water~\cite{OOWB2021-scat-phys}.

Scattering by pure water is often modeled by \emph{Einstein-Smoluchowski Theory}. Even the purest water has some scattering. Water exhibits fluctuations in density and in the IOR due to the random movement of water molecules. The Einstein and Smoluchowski formula is used as the volume scattering function. Elastic scattering by water has also similarities in angular shape and spectral behavior to \emph{Rayleigh scattering} (introduced in~Section~\ref{sec:inscat}). There is some effect of salinity on the scattering coefficient by water (30$\,\%$ increase for range of salinities normally observed in the ocean), much less so of temperature (4$\,\%$ between 0 and 26$\,^{\circ}\,C$, and even less for pressure~\cite{OOWB2020-const-water}. 

% TODO cite https://www.oceanopticsbook.info/view/scattering/physics-scattering
% TODO cite https://www.oceanopticsbook.info/view/optical-constituents-of-the-ocean/water

\subsection{Index of Refraction of Water}
Index of refraction (IOR), in its full form, defines how a~medium affects light of a~given wavelength. The complex part relates to the absorption of water and varies with wavelength. The real part of IOR of water also varies with wavelength in the visible spectrum. The IOR of water is influenced more or less by the following quantities: temperature, salinity and pressure. Higher temperature and salinity causes higher scattering. On the other hand, pressure has small effect on scattering, because water is very incompressible.

For near-surface waters, IOR is in range of 1.33\,to\,1.35. Including extreme depths, overall, IOR of water varies by less than 3\,\% over the entire range of these parameters~\cite{OOWB2020-const-water}.
%\footnote{ Refer to the source for tables showing various values of IOR for varying pressure, temperature and salinity at different wavelengths. }. 

% TODO cite https://www.oceanopticsbook.info/view/optical-constituents-of-the-ocean/water

%he temperature-dependent ﬂuctuations in the index of refraction appear in Eq. (5) via the various thermodynamic quantities. The equation shows that scattering will be stronger for a higher temperature; this is because the random molecular motions are then greater, and consequently the density ﬂuctuations are larger. 


\subsection{Volume Scattering Function}
% https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/volume-scattering-function-vsf
% https://www.oceanopticsbook.info/view/scattering/scattering-overview

This section is based on Chapters~\cite{OOWB2021-scat-intro},~\cite{OOWB2021-iaop-vsf} from Ocean Optics Web Book~\cite{OOWB2020}.

\emph{Volume scattering function} (VSF) describes the angular distribution of unpolarized light scattered from its incident direction $\omega$ into direction $\omega'$ at a~wavelength $\lambda$. For turbulent aquatic environments, scattering is almost always assumed to be \emph{azimuthally symmetric} about the incident direction. In order words, the scattering is dependent only on the angle $\theta$~between the two directions. This holds true for spherical particles or randomly oriented non-spherical particles for unpolarized incident irradiance, which are present in such environments.

The VSF is denoted by $\beta(\theta, \lambda)$ and its units are [m$^{-1}$ sr$^{-1}$]. In comparison to the \emph{volume scattering phase function} (or simply phase function, as introduced in Section~\ref{sec:inscat}), $p(\theta, \lambda)$ [$sr^{-1}$], the VSF is related to radiant intensity emanating at an angle $\theta$ from a~volume, rather than from a~single particle. It is the scattered intensity per unit incident irradiance per unit volume. Phase function and the VSF are related in the following way:
\begin{align}
    p(\theta, \lambda) \equiv \frac{\beta(\theta, \lambda)}{\scatterCoef (\lambda)}\quad[sr^{-1}].
\end{align}
This can be seen as partitioning the $\beta(\theta, \lambda)$ into the magnitude of the total scattering, expressed by $\scatterCoef(\lambda)$ in [m$^{-1}$], and the angular distribution of the scattered light, which is $p(\theta, \lambda)$ in [sr$^{-1}$]. In other words, a~phase function can be thought of, as the corresponding VSF without the magnitude of the intensity of the scattered light.

Integrating $\beta(\theta, \lambda)$ over all directions (solid angles) gives the \emph{scattering coefficient}, or the total scattered power per unit incident irradiance and unit volume of water:
\begin{align}
    \scatterCoef (\lambda) = 2\pi \int_{0}^{\pi} \beta(\theta, \lambda) \sin \theta \text{ d}\theta .
    \label{eq:scatter-coef}
\end{align}
The term $2\pi$ is based on an assumption, that the scattering process is \emph{azimuthally symmetric} about the incident direction, that is, the scattering depends only on the \emph{scattering angle} $0 \leq \theta \leq \pi$. This assumption holds if the medium is isotropic for unpolarized incident light. This is often taken as a~reasonable assumption for natural waters in which the particles are randomly oriented by turbulence.

This integration can be divided into \emph{forward scattering}, $0 \leq \theta < \pi/2$, and \emph{backward scattering}, $\pi/2 \leq \theta \leq \pi$, part. Then the corresponding forward and backward scattering coefficients are, respectively:
\begin{align}
    \forwardScatter (\lambda) &= 2\pi \int_{0}^{\pi/2} \beta(\theta, \lambda) \sin \theta \text{ d}\theta ,\\
    \backScatter (\lambda) &= 2\pi \int_{\pi/2}^{\pi} \beta(\theta, \lambda) \sin \theta \text{ d}\theta .
\end{align}

% https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/visualizing-vsfs
Figure~\ref{fig:vsfs} illustrates the range of VSFs that can be found in various oceanic waters. VSFs are plotted as functions of the scattering angle, with the ordinate being in logarithmic scale to highlight the wide range of VSFs and to show the smallest scattering angles. Oceanic VSFs typically vary by 5 or 6 orders of magnitude over the range of very small to large scattering angles. The different VSFs can also vary by over two orders of magnitude at a~given scattering angle. The scattering phase functions are plotted in the rightmost graph. The functions appear similar at low angles due to the logarithmic ordinate, but they vary by almost an order of magnitude at large scattering angles~\cite{OOWB2020-iaop-vis}.
The asymmetry parameter $g$ of the phase function is often used as a~convenient measure of the shape of the phase function. Typical ocean waters have $g$ values in range of 0.8 to 0.95. In other words, oceanic phase functions have very large values for small angles $\theta$.
% TODO cite https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/volume-scattering-function-vsf
% TODO figure 1

\begin{figure}[H]
\begin{subfigure}{.495\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/VSFs-cropped.png}
\end{subfigure}
\begin{subfigure}{.495\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/PFs-cropped.png}
\end{subfigure}
\caption{Example of volume scattering functions (on the left) and the corresponding scattering phase functions (on the right). Measured in clear, open ocean water (the blue curve), in a~harbor (the green curve), in very productive coastal waters (the red curve). The corresponding scattering coefficients $b$~and the asymmetry parameter $g$ are also shown. 
(\emph{The graphs are provided only in this form by Mobley~\cite{OOWB2020-iaop-vis}.}).
}
\label{fig:vsfs}
% https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/visualizing-vsfs
\end{figure}

\section{Classification of Natural Waters}
\label{sec:classify-waters}
%https://www.oceanopticsbook.info/view/inherent-and-apparent-optical-properties/classification-schemes
Optical oceanography classifies waters into water types based on their optical properties. For instance, waters can be classified based their colors, then the underlying optical property is primarily absorption (as scattering contributes to intensity but less to color). The most common classifications to water types are based on~\cite{OOWB2021-iaop-clas}:
\begin{itemize}
    \item Visual comparison of the water color to the color of a~standard color sample using human eye.
    \item Analysis of AOP spectra (such as Jerlov water types).
    \item The relative importance of chlorophyll in determining the optical properties of water (Case\,1 vs. Case\,2).
    \item Analysis of IOPs.
\end{itemize}

The most commonly used classifications in bio-optical models, which are used for rendering waters, are the \emph{Jerlov classification} and \emph{Case\,1 versus Case\,2 water types}. These two classifications are described in the sections below.

\subsection{Jerlov Classification}
\label{sec:jerlov-waters}

The Jerlov classification divides water bodies based on the spectral diffuse attenuation coefficient $\diffuseAttenDown$. The $\diffuseAttenDown$ and the absorption coefficient are related and are linked to the ocean color observed. Observations of the $\diffuseAttenDown$ were discretized into five typical open-ocean types, labelled, I, IA, IB, II, and III, and nine typical coastal types, labelled 1 to 9~\cite{OOWB2021-iaop-clas}. 

Today, the Jerlov water types are commonly used as a~convenient one-parameter classification scheme describing water clarity for essentially all natural bodies of water~\cite{SM2015}. Solonenko and Mobley, 2015~\cite{SM2015} present a~method to determine IOPs for all Jerlov water types. They provide tables of sets of absorption and scattering coefficients for the different water types.

% TODO image
% https://www.researchgate.net/publication/356445702_Water_Classification_Based_on_Underwater_Monocular_Image
\begin{figure}[H]
    \centering
\begin{subfigure}{.7\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/jerlov.pdf}
\end{subfigure}
\caption{Illustrates the color of a~completely white object seen in each of the 10\,Jerlov water types at varying depths. (\emph{Image adapted from Toledo et al.~\cite{TVD2021}}.)}
\label{fig:jerlov}
\end{figure}


\subsection{Case\,1 versus Case\,2 Water Types}

The changes in spectral shape of the irradiance reflectance $R = E_u / E_d$ (described in Section~\ref{sec:iops-aops}) separates waters into two types or \emph{cases}. Case\,1 type are waters with concentration of phytoplankton that is high compared to other particles. In contrast, in case\,2 the inorganic particles are dominant. In both cases, the dissolved yellow substance (CDOM) is present in variable amounts. An ideal case\,1 is water that is pure culture of phytoplankton and an ideal case\,2 is a~suspension of nonliving material with a~zero concentration of pigments.

Today, the definitions of the two types are used in relation to the influence of constituents on optical parameters:
\begin{itemize}
    \item  Case\,1: waters whose optical properties are determined primarily by phytoplankton and other constituents\,---\,CDOM and detritus and degradation produces\,---\,then co-vary together with the phytoplankton concentration.
    \item Case\,2: rest of the waters, that is, waters whose optical properties do not covary with the phytoplankton concentration, but instead are significantly influenced by other constituents such as mineral particles, CDOM, or microbubbles.
\end{itemize}

% TODO figure 
% https://www.researchgate.net/publication/216027395_Remote_Sensing_of_Ocean_Colour_in_Coastal_and_Other_Optically-Complex_Waters
\begin{figure}[ht]
    \centering
\begin{subfigure}{.6\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/ocean-optics/cases.pdf}
\end{subfigure}
\caption{Illustration of Case\,1 and Case\,2 waters as top and bottom parts of a~spectrum between varying concentrations\,---\,the representative colors\,---\,of the constituents. (\emph{Image adapted from IOCCG~\cite{Ioc2000}.})
}
\label{fig:cases}
\end{figure}

In reality, natural waters are not simply either Case\,1 or Case\,2, but there is a~continuous spectrum of water optical properties and their dependence on chlorophyll and other water constituents~\cite{OOWB2021-iaop-clas}. Figure~\ref{fig:cases} shows an illustration of such a~spectrum using a~gradient.

% ================================================================================
% ================================================================================
\chapter{Existing Approaches for Real-Time Ocean Simulation and Rendering}
\label{ch:ocean}

This chapter provides an overview of real-time ocean simulation and rendering techniques in computer graphics. These techniques are concerned with a~range of topics: modeling, animation and simulation of ocean's surface or water body particles, rendering of the ocean surface with or without the light transport inside the water body, rendering of effects across the ocean surface, such as foam and spray, and lastly rendering of water volumes underwater with or without the underwater phenomena, such as caustics. The survey presented by Darles et al.~\cite{DCGG2011} is the testament to the sheer volume of these techniques till the year 2011, since then the number of techniques has only increased.

The previous chapter showed the complexity and variability of inherent optical properties of oceans. Showing that the light transport in such media is dependent not just on the optical properties of the media itself, but also on the boundary conditions, that is, the conditions of the sky, the water surface, and the ocean bottom. Also, having in mind the aim of this thesis, the following types of techniques are the focus of this thesis:
\begin{itemize}
    \item water surface simulation and animation for deep ocean waters, in Section~\ref{sec:ocean-water},
    \item bio-optical models used for modeling properties of water volumes, in Section~\ref{sec:ocean-models},
    \item participating media rendering of water volumes, either from above or below the surface, along with more general approaches that might be suitable, in Section~\ref{sec:ocean-media},
    \item procedural generation topics, that are used for modeling or simulation of both water surfaces and surfaces of the ocean bottom, in~Section~\ref{sec:proc-gen},
    \item underwater effects, such as caustics, are the subject of the next chapter.
\end{itemize}

% https://www.oceanopticsbook.info/view/scattering/scattering-overview

% ================================================================================
\section{Ocean Surface Simulation Techniques for Deep Waters}
\label{sec:ocean-water}

Ocean surface simulation refers to methods that model and animate the detailed surfaces of water bodies, or general fluid simulations that can be applied to ocean simulation.

Ocean is complex due to its highly dynamic behavior, ranging from quiet sea, to agitated ocean with small turbulent waves to enormous shorebreaks. This motion is influenced by multiple phenomena occurring at small and large scales. Physically-based methods utilize information obtained by the field of \emph{oceanography}. Oceanography defines the behavior of the ocean surface (or ocean dynamics) depending on its location. It classifies the wave motion of the ocean surface into categories depending its location: deep ocean waters (far from the coast), intermediate areas or shallow water areas (close to the shore). Waves of deep waters exhibit characteristic oscillatory behavior, while waves of shallow waters interact with shallow ocean floor, causing them to break near the shore. Each category is modeled using a~different set of techniques and approaches~\cite{DCGG2011}.

% TODO images of tecniques of shorebreaks and deep wates


Darles et al.~\cite{DCGG2011} provides a~survey of ocean simulation methods. They divide the methods into two categories based on the models built to simulate certain types of ocean dynamics: 
\begin{itemize}
    \item \textbf{Models for deep water simulation}. These approaches use empirical laws from oceanographic research to build parametric, spectral or hybrid models. They, generally focus on simulation of deep water waves, without breaking waves. 
    \item \textbf{Fluid simulation for shallow water}. These are physically-based methods, that rely on \emph{Navier-Stokes Equations} (NSE) to simulate the dynamic behavior of ocean waves as they interact with the bottom. They are mainly used to simulate breaking waves near the shore.
\end{itemize}

This section further investigates methods for simulation of deep ocean waves of the free ocean surface. These methods build models either based on oceanographic observations, such as the parametric models, or based on (primarily) oceanographic measurements, such as spectral models. The goal of these methods is to compute the representation of the surface in the spatial domain, that is, to compute the concrete horizontal and vertical coordinates of each point on the surface. Methods based on parametric models, typically describe the surface directly in the spatial domain. While the methods based on spectral models, describe the surface in spectral domain, that is, the distribution of amplitudes and frequencies according to the modeled wave spectrum, and then they must convert this representation to the spatial domain~\cite{DCGG2011}.

%$\omega$ is used for angular frequency 

% --------------------------------------------------------------------------------
%\subsection{Spatial Domain Approaches}
\subsection{Spatial Domain Approaches}
These methods generate and animate the ocean surface directly in spatial domain. The surface is usually computed as sum of periodical functions that evolve temporarily using a~phase difference. 

\textbf{Sinusoidal Wave Model} serves as a~basis for later techniques. It simplifies the complex motion of the ocean surface as a~two-dimensional ideal surface represented by a~sine wave with constant parameters. It also assumes uniform mean depth and small amplitudes in relation to the size of the water body. The wave travels in the horizontal direction and is vertically displaced for each horizontal coordinate $x$ at time $t$ giving the elevation as:
\begin{align}
    \eta (x, t) = A \cos (kx - \omega t) ,
\end{align}
where $A$~is the amplitude, $k = \frac{2\pi}{\lambda}$ is the \emph{wavenumber}, $\omega = 2\pi f$ is the \emph{angular frequency} in [rad/s], $f = \frac{1}{T}$ is the \emph{wave frequency} in [Hz], and $T$~is the wave period in [sec]. Figure~\ref{fig:sinusoid} shows such a~sinusoidal wave. One may assume a~fixed position on the wave, then the surface elevation for this point evolves through time. Or, one may assume a~fixed point in time, then the surface elevation can be observed for all positions. The highest points of the wave height are called \emph{crests}, and the lowest points are called \emph{troughs}. This model thus provides a~wave where all crests and troughs have the same elevation~\cite{Gam2018}. %The difference in elevation between crest and through is the \emph{wave height}. It is defined as $H = 2A$~\cite{Gam2018}.

\begin{figure}[H]
    \centering
\begin{subfigure}{.8\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/surfaces/sinusoidalWave.pdf}
\end{subfigure}
\caption{A sinusoidal wave model represented by surface elevation $\eta$ at position $x$ or at time $t$ with amplitude A, wavelength $\lambda$, and wave period $T$.}
\label{fig:sinusoid}
\end{figure}


This method can be extended into a~sum of three-dimensional waves represented as a~height map for each horizontal point $(x, z)$:
\begin{align}
    h(x, z, t) = -y_0 + \sum_{i=1}^{N} A_i \cos(k_i^{x} x + k_i^z z - \omega_i t) ,
\end{align}
then the wavenumber becomes \emph{wavevector}, $\*{k} = (k_x, k_z)$, the direction of travel of the wave, the wavenumber becomes the magnitude of the wavevector: $k = ||\*{k}||$, and $y_0$ is the height of the free surface~\cite{DCGG2011}.

Instead of sinusoids, Fournier~\cite{FR1986} propose method called \textbf{Gerstner waves}. This approach can represent a~variety of wave shapes, including waves with a~sharp crest. The method allows for control of the waves to even simulate shape of breaking waves. To yield a~more realistic result, a~sum of Gerstner waves is commonly evaluated~\cite{DCGG2011}. Figure~\ref{fig:gerstner} shows an ocean surface obtained using the method by Fournier et al.~\cite{FR1986}.

\begin{figure}[H]
    \centering
\begin{subfigure}{.4\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/surfaces/fournier.png}
\end{subfigure}
\caption{Ocean surface obtained using the Gerstner waves method by Fournier et al.~\cite{FR1986}.}
\label{fig:gerstner}
\end{figure}

% --------------------------------------------------------------------------------
% TODO should introduce them more closely? no
%- spatial methods compute sum of periodical functions
%- spectral approaches use a wave spectrum

\subsection{Spectral Domain Approaches}
Spectral approaches describe the ocean surface using \emph{wave spectrum}, a~distribution of wave energy among different wave frequencies or wavelengths obtained from theoretical or measured data. The spectral representation is then transformed into the spatial representation using \emph{inverse fast Fourier transform} (IFFT).

The ocean surface is assumed to be a~linear superposition of a~large number of sinusoids with many different wavelengths, frequencies and phases, each traveling in different directions. The wave spectrum defines the distribution of wave energy among the different wave frequencies or wavelengths.

\textbf{Tessendorf}~\cite{Tes2001} computes wave height field at horizontal point $\*{x} = (x,z)$ as the sum of sinusoids with complex, time-dependent amplitudes:
\begin{align}
    h(\*{x}, t) = \sum_{\*{k}} \tilde{h}(\*{k}, t) e^{i\*{k} \cdot \*{x}} ,
    \label{eq:tess}
\end{align}
where $t$ is the time, $\*{k}$ is the wavevector. The height amplitude Fourier components, $\tilde{h}(\*{k}, t)$, determine the structure of the surface and are based on oceanographic data. Furthermore, the slopes of the wave height field can also be computed using IFFT, similarly to finite differences.

This approach produces waves on a~patch with horizontal dimensions $L_x \times L_z$ of certain number of discrete sample points on a~grid $N \times M$. This patch is perfectly periodic, so one such patch can be tiled seamlessly. There are certain restrictions on the dimensions to fully utilize the computational power of FFT: only square patches of sampled points can be produced: $N = M$, ranging from 16 to 2048, and only in powers of two~\cite{Tes2001}. Figure~\ref{fig:tess-output} shows an example of the computed wave height field.

There are various wave spectra to model deep sea waves at various states, e.g. turbulent sea, or less crested waves of the free ocean. For instance, Tessendorf~\cite{Tes2001} presents a~method to generate a~random ocean wave height field using Gaussian pseudo-random generator and a~theoretic wave spectrum, called \emph{Phillips spectrum}. This method is based on the observed Gaussian distribution of waves in deep waters and it was later used in famous computer-generated scenes for \emph{Waterworld} and \emph{Titanic productions}~\cite{DCGG2011}, as shown in Figure~\ref{fig:tess-output}. Gamper~\cite{Gam2018} provides a~comparison of various wave spectra.


% TODO figure of Tessendorf s tess-output

\begin{figure}[H]
    \centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=0.74\linewidth]{figures/surfaces/waveheight-phillips.png}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/surfaces/tessendorf.png}
\end{subfigure}
\caption{The right image shows the wave height field generated with Phillips spectrum. The left image shows a~rendered ocean surface. (\emph{Renderings by Tessendorf~\cite{Tes2001}.})
}
\label{fig:tess-output}
\end{figure}


Typically, the waves produced by the FFT methods have rounded peaks and troughs, but in even in fairly good weather the waves are sharply peaked at their tops, and flattened at the bottoms. Tessendorf~\cite{Tes2001} also describes method called \textbf{Choppy waves} to create more realistic waves. This is achieved by additional displacement of the waves in horizontal dimensions. The displacement vector field is computed using the Fourier amplitudes of the height field. 
%Figure~\ref{fig:choppy} shows the difference between the original and the choppy waves method. 
The method can cause some waves to pass through itself and invert due to large enough amplitudes. %similarly to Gerstner waves as shown in Figure~\ref{fig:gerstner}, 
This can be tested for by computing the Jacobian of the transformation from the original horizontal point to the displaced one.



\subsection{Limitations and Mitigation Strategies for Spectral Domain Approaches}

Although the results obtained look realistic, there are limitations of these methods due to the fixed resolution of the height map: the waves look unrealistic on close-up, visible tiling or repetition artefacts are produces when look from a~distance. This requires a~new surface to be generated with higher resolution, which might not be feasible. Instead, implementations employ various LOD schemes and/or methods to remove the tiling artifact. 

There are hybrid approaches that use combination of spatial and spectral definitions to gain advantages of both. For instance, using a~combination of a~spectrum and a~sum of sinusoids to simulate breaking waves, and Tessendorf's approach~\cite{Tes2001} for deep waters.

Darles et al.~\cite{DCGG2011} provides an overview of LOD approaches to Tessendorf's method~\cite{Tes2001}. Some approaches represent the ocean with multiple maps at various resolutions to account for low frequencies (i.e., the global motion of the waves) and high frequencies (the fine details). Other approaches use adaptive sampling of spectrum depending on the distance to the viewer. 

Rydahl~\cite{Ryd2009} presents an approach to remove tiling by transforming the world coordinates using fractal noise before sampling the displacement maps of the surface.  Nvidia, 2011~\cite{Nvi2011} blends Perlin noise and FFT generated crests to gain detail at near sight (due to FFT) and remove tiling patterns at distance (due to Perlin noise). Dupuy and Bruneton~\cite{DB2012} instead compute a~number of patters, each sampling a~different part of the wave spectrum to remove artifacts viewed at near and far distances. This way the ocean surface has an analytical expression.

%Noise functions are described in Section~\ref{sec:proc-gen} about procedural generation.

% ============================================================================
% ============================================================================
% ============================================================================
\section{Bio-Optical Models}
\label{sec:ocean-models}
% PHYSICAL MODELS: (part of Light transport WITHIN THE WATER)

This section provides an overview of various bio-optical models employed by different techniques for rendering light transport in water volumes, not just real-time techniques. The bio-optical model is the term for the underlying physical model used in these techniques. It is made up of the optical parameters of the ocean, thus being physically-based. These parameters are properties of participating media, e.g., absorption or scattering coefficients, of various oceanic constituents, as described in Chapter~\ref{ch:optics}. 

This section focuses solely on the way the absorption and scattering are modeled: which constituents are chosen, and how are their optical parameters modeled. Typically, these include the optical parameters shown in Table~\ref{tab:models-notation}. Furthermore, only elastic scattering is considered.

\begin{table}[h]
\centering
\begin{tabular}{ll}
\hline
\textbf{Symbol} & \textbf{Description}            \\ \hline
$\sigma_a$      & Absorption coefficient          \\
$\sigma_s$      & Scattering coefficient          \\
$\sigma_t$      & Extinction coefficient          \\
$p(\theta)$     & Phase function                  \\ \hline
$\sigma_{bs}$   & Backscattering coefficient      \\
$K_d$           & Diffuse attenuation coefficient \\
$C$             & Chlorophyll concentration       \\ \hline
\end{tabular}
\caption{The notation used for the various optical parameters.}
\label{tab:models-notation}
\end{table}

\noindent\textbf{Premože and Ashikhmin}~\cite{PA2001} use bio-optical model that supports many different water types, such as deep ocean water, muddy coastal water, and fresh water bodies. They assume homogeneous water bodies and define values of four parameters: $\scatterCoef, \absorbCoef, \backScatter$, and $\diffuseAttenDown$ (see Table~\ref{tab:models-notation} for the used notation). All the optical parameters are wavelength-dependent, with the exception of the scattering coefficient, which is derived from~single-wavelength measurements for certain types of waters. 

% They use the classification by Jerlov~\cite{TODOCITEJERLOV} 
%(TODO? see Section~\ref{sec:JerlovWaterTypes} for detailed description)  TODO should explained in sep. section?
%to divide natural waters into types based on the $\diffuseAttenDown (\lambda)$:
%\begin{itemize}
%\item Types I to III: water types for open ocean waters, where from the clearest (type\,I) to the most turbid (type\,III).
%\item Types 1\,---\,9: correspond to coastal waters, from the clearest (type\,1) to the most turbid (type\,9).
%\end{itemize}

% TODO explain turbidity in ocean section
% TODO needed?
%The authors introduce a~single cumulative \emph{turbidity parameter} to obtain spectral data of $\diffuseAttenDown (\lambda)$ and values of $\scatterCoef$~by interpolation of input data. 
The absorption coefficient is approximated directly from the value of $\diffuseAttenDown$. Any depth dependence on the parameter is ignored.
%, where the used wavelength-dependent values of $\diffuseAttenDown $ are in Table~\ref{tab:PA2001-diffuse-atten-coef}.
Scattering coefficient is computed separately for two cases of waters, each with a~different major cause of scattering: 1) Cloudy or muddy waters due to suspended particles (of silt, plant fibers, sawdust, chemicals, and microorganisms) from single-wavelength measurements, 2) waters with the abundance of phytoplankton and other particulate and dissolved material. The second case, takes into account the molecular scattering coefficient of water, and parameterizes the spectral scattering coefficient of the pigment (chlorophyll\,A and phytoplankton) by their concentrations. The renderings using this model are shown in Figure~\ref{fig:PA2001-outputs}.

\ifdefined\SHOWEQUATIONS
Scattering is modeled separately for two cases of waters, each with a~different major cause of scattering:
\begin{itemize}
    \item Cloudy or muddy waters due to suspended particles (of silt, plant fibers, sawdust, chemicals, and microorganisms): to obtain $\scatterCoef (\lambda)$ and $\backScatter (\lambda)$, the following relations are used:
    \begin{align}
        \scatterCoef(\lambda ) &= \scatterCoef(\lambda_0) \frac{m\lambda + i}{m\lambda_0 + i} , \\
        \backScatter(\lambda ) &= 0.01829 \scatterCoef(\lambda ) + 0.00006 ,
    \end{align}
    where $m = -0.00113$, $i = 1.62517$, and $\scatterCoef(\lambda_0)$ is set using data from Mobley~\cite{Mob1994} at $\lambda_0 = 514\,nm$ depending on the water type (see  Table~\ref{tab:PA2001-scatter-coef-lambda0} for the values of $\scatterCoef(\lambda_0)$).

    \item Waters with the abundance of phytoplankton and other particulate and dissolved material:  \emph{pigment concentration parameter} $C$,~related to concentration of Chlorophyll\,A and Phytoplankton (in $mg/m^{3}$), is introduced to obtain $\backScatter$ instead:
    \begin{align}
        \backScatter(\lambda ) &= 0.5 b_w(\lambda ) + B_b(\lambda ) b_p , \\
        B_b(\lambda ) &= 0.002 + 0.02 (0.5 - 0.25 \log_{10} C) \frac{550}{\lambda} , \\
        b_p(\lambda ) &= 0.3 C^{0.62} ,
    \end{align}
    where $b_w (\lambda)$ is the molecular scattering coefficient of water, $B_b(\lambda)$ is the ratio of backscattering and scattering coefficients of the pigments, and $b_p$ is the scattering coefficient of the pigment. The relations and values of $C$~are provided from Morel (1988)~\cite{Mor1988}. Table~\ref{tab:PA2001-concentration} shows the typical concentration of chlorophyll and phytoplankton only in open ocean waters. Figure~\ref{fig:PA2001-outputs} shows renderings of varying optical parameters.
\end{itemize}
\fi

\ifdefined\SHOWDATATABLES
\begin{table}[h]
\centering
\begin{tabular}{l|cccccccc}
\hline
\multicolumn{1}{c|}{$\lambda$} & \multicolumn{8}{c}{Water type} \\ \hline
$[nm]$ & I & II & III & 1 & 3 & 5 & 7 & 9 \\
310 & 15 & 37 & 65 & 180 & 240 & 350 &  &  \\
350 & 6.2 & 17.5 & 32 & 120 & 170 & 230 & 300 & 390 \\
375 & 3.8 & 12.2 & 22 & 80 & 110 & 160 & 210 & 300 \\
400 & 2.8 & 9.6 & 18.5 & 51 & 78 & 110 & 160 & 240 \\
425 & 2.2 & 8.1 & 16 & 36 & 54 & 78 & 120 & 190 \\
450 & 1.9 & 6.8 & 13.5 & 25 & 39 & 56 & 89 & 160 \\
475 & 1.8 & 6.2 & 11.6 & 17 & 29 & 43 & 71 & 123 \\
500 & 2.7 & 7.0 & 11.5 & 14 & 22 & 36 & 58 & 99 \\
525 & 4.3 & 7.6 & 11.6 & 13 & 20 & 31 & 49 & 78 \\
550 & 6.3 & 8.9 & 12.0 & 12 & 19 & 30 & 46 & 63 \\
575 & 8.9 & 11.5 & 14.8 & 15 & 21 & 33 & 46 & 58 \\
600 & 24 & 26 & 29.5 & 30 & 33 & 40 & 48 & 60 \\
625 & 31 & 33.5 & 37.5 & 37 & 40 & 48 & 54 & 65 \\
650 & 36 & 40 & 44.5 & 45 & 46 & 54 & 63 & 76 \\
675 & 42 & 46.5 & 52 & 51 & 56 & 65 & 78 & 92 \\
700 & 56 & 61 & 66 & 65 & 71 & 80 & 92 & 110 \\ \hline
\end{tabular}
\caption{Diffuse attenuation coefficient for downwelling irradiance, $\diffuseAttenDown(\lambda )$, for different water types, in $10^{-2} m^{-1}$~\cite{PA2001}}
\label{tab:PA2001-diffuse-atten-coef}
\end{table}

\begin{table}[h]
\parbox{.45\linewidth}{
    \centering
    \begin{tabular}{l|l}
    \hline
    Water type             & $\scatterCoef [m^{1}]$ \\ \hline
    Clear ocean (type I)   & 0.037       \\
    Coastal ocean (type 1) & 0.219       \\
    Turbid harbor (type 9) & 1.824       \\ \hline
    \end{tabular}
    \caption{Total scattering coefficient $\scatterCoef$ at $\lambda_0 = 514\,nm$ for different water types~\cite{PA2001}.}
    \label{tab:PA2001-scatter-coef-lambda0}
}
\hfill
\parbox{.45\linewidth}{
    \centering
    \begin{tabular}{l|l}
    \hline
    Water type & $C [\frac{mg}{m^3}]$ \\ \hline
    type I     & 0.0 -- 0.1           \\
    type II    & 0.5                  \\
    type III   & 1.5 -- 2.0           \\ \hline
    \end{tabular}
    \caption{Pigment concentration $C$ for different open water types~\cite{PA2001}.}
    \label{tab:PA2001-concentration}
}
\end{table}
\fi

\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/PA2001/open-deep-1.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/PA2001/open-deep-2.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/PA2001/tropical-3.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/PA2001/muddy-coastal-4.png}
\end{subfigure}
\caption{Renderings of different water types, from left to right: two variations of open ocean deep water, tropical water, and muddy coastal water. (\emph{Renderings from Premože~\cite{PA2001}.)}}
\label{fig:PA2001-outputs}
\end{figure}

\noindent\textbf{Cerezo and Serón~\cite{CS2002}} model waters with high phytoplankton concentration\,---\,case 1 waters\,---\,in which optical parameters can be correlated to chlorophyll concentration.

Absorption coefficient is the sum of the spectral absorption coefficient of pure water, and the chlorophyll-specific absorption coefficient parameterized by chlorophyll concentration~$C$.

Scattering coefficient is modeled using single total scattering coefficient, which is only related to $C$. The phase function takes into account: 1) scattering by pure water, modeled by Rayleigh scattering, that accounts for the anisotropy of the water molecules, 2) scattering by particulate matter, modeled from measured data to account for high peaks in the forward direction of the phase function. Scattering coefficient of water is derived from single-wavelength measurements, and the scattering coefficient of particulate matter is taken simply as the difference between the total scattering coefficient and scattering coefficient of water.
\ifdefined\SHOWEQUATIONS
For the absorption, the absorption by phytoplankton cells due to its pigments is taken into account.
Specifically, the pigment chlorophyll which has strong absorption bands in the blue and in the red part of the visible spectrum of light. The absorption coefficient in the terms of chlorophyll concentration $C$~(in $mg/m^{3}$) is expressed as: 
\begin{align}
    \absorbCoef (\lambda) = (a_w (\lambda) + 0.06 a_c^{*}(\lambda) C^{0.65})
                            (1 + 0.2 e^{-0.014 (\lambda - 440)}) ,
\end{align}
where $a_w(\lambda)$ is the absorption coefficient of pure water, $a_{c}^{*}(\lambda)$ is non-dimensional chlorophyll-specific absorption coefficient. $C$~can range from $0.01\,mg/m^{-3}$ in the clearest open ocean waters, to $10\,mg/m^{-3}$ in productive coastal upwelling regions, to $100\,mg/m^{-3}$ in estuaries or lakes.

For scattering, the following relation for the \emph{total scattering coefficient} is employed:
\begin{align}
    \scatterCoef(\lambda) = \left( \frac{550}{\lambda} \right) 0.30 C^{0.62}  ,
\end{align}
and the phase function is modeled with two terms:
\begin{align}
    p(\theta ) &=
        \underbrace{
            \frac{ b_w(\lambda) }{ \scatterCoef(\lambda) } \tilde{\beta}_w(\theta)
        }_\text{water molecules}
        +
        \underbrace{
            \frac{ b_p(\lambda) }{ \scatterCoef(\lambda) } \tilde{\beta}_p(\theta)
        }_\text{particles} ,
\end{align}
where the first term corresponds to the contribution of water molecules and the second one to the particles contribution. $b_w$ and $b_p$ are the scattering coefficients for the water molecules and for the other particles, respectively:
\begin{align}
    b_w (\lambda) &= 16.06 \left( \frac{\lambda_0}{\lambda} \right)^{4.32}
                    \beta_w (90^{\circ}, \lambda_0) , \\
    b_p (\lambda) &= \scatterCoef( \lambda) - b_w ( \lambda) .
\end{align}
The phase function that characterizes scattering in pure water is defined as:
\begin{align}
    \tilde{\beta}_w (\theta) = 0.06225 (1 + 0.835 \cos^{2} \theta ) ,
    \label{eq:CS2002-phase-water}
\end{align}
The phase function for pure water is similar to Rayleigh scattering, which is almost isotropic, except for the $0.835$~term which is attributed to the anisotropy of the water molecules. The authors then present data that shows how even a~slight amount of particulate matter, which is the case even for the clearest water, causes high peaks in the forward direction of the phase function and an increase in the scattering coefficient by at least a~factor of ten. To account for this, the particle phase function $\tilde{\beta}_p$ is modeled in a~form given in a~graph for each discrete scattering angle.
\fi
Figure~\ref{fig:CS2002-outputs} shows renderings of varying values of the optical parameters.

\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/CS2002/pure-sea-water-1.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/CS2002/sea-low-phyto-2.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/CS2002/sea-some-phyto-3.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/CS2002/sea-high-phyto-4.png}
\end{subfigure}
\caption{Renderings of seawater with different concentrations of phytoplankton, from left to right: pure seawater only, seawater with low concentration, seawater some concentration, seawater with high concentration. The selected values of parameters are also shown in the source.
(\emph{Renderings from Cerezo et al.~\cite{CS2002}.}).
}
\label{fig:CS2002-outputs}
\end{figure}

\noindent\textbf{Iwasaki et al.}~\cite{IDN2003} and \textbf{Aranha}~\cite{Ara2005} model the phase function using Henyey-Greenstein phase function to account mainly for forward scattering in water. The scattering coefficient is then computed in a~standard way by integrating the volume scattering phase function (in Equation~\ref{eq:scatter-coef}).

\vspace{1em}

\noindent\textbf{Gutierrez et al.}~\cite{GMAS2005} (2005) simulates elastic and inelastic scattering in deep ocean waters by taking into account concentrations of phytoplankton (pigment chlorophyll) and CDOM (yellow matter). Absorption coefficient is modeled as the sum of absorption coefficients of the following constituents from spectral measurements: pure seawater, phytoplankton, and CDOM. Only the absorption coefficient of CDOM is computed from single-wavelength empirical absorption. Varying concentrations of chlorophyll and CDOM, both wavelength-independent, are taken account for both absorption and scattering. The used phase function accounts for: 1) scattering by pure seawater, modeled by Rayleigh scattering, in the same manner as Cerezo et al.~\cite{CS2002}, 2) scattering by suspended particles, modeled using Henyey-Greenstein phase function, with a~high value of the mean cosine parameter $g$ to account for high forward scattering.

\ifdefined\SHOWEQUATIONS
To model absorption, the absorption coefficients of the following constituents are taken into account:
\begin{itemize}
    \item pure seawater: peaks at 760\,nm and reaches minimum at 430\,nm,
    \item phytoplankton: reaches maximum at 350\,nm, rapidly decreasing to almost zero beyond 500\,nm,
    \item yellow matter: $a_y(\lambda) = a_y(440)^{-0.014 (\lambda - 440)}$, where $a_y(440)$ is the empirical absorption at 440\,nm in Mobley~\cite{Mob1994} (1994).
\end{itemize}
For both absorption and scattering, the model takes into account chlorophyll concentration $\Gamma^{c}$, ranging from $0.01$\,to $0.1$, and CDOM concentration $\Gamma^y$, ranging from 0.005 to 0.025, both being wavelength-independent. 

The used phase function is composed of phase functions for pure seawater and for scattering by suspended particles: $p = p_w + p_p$. Where $p_w$ is the phase function used in Cerezo et al.~\cite{CS2002}, defined in Equation~\ref{eq:CS2002-phase-water}, $p_p$ is modelled using Henyey-Greenstein phase function (in Equation~\ref{eq:phase-HG}) with $g = 0.924$. The actual values of the parameters are provided in the graphs of the source material.
\fi

\vspace{1em}

\noindent\textbf{Gutierrez et al.}~\cite{GSMA2008} (2008) builds on the work of Gutierrez et al.~\cite{GMAS2005} (2005) and extends the bio-optical model to also account for complex inelastic scattering events that occur underwater, such as \emph{fluorescence} and \emph{Raman scattering}. Both absorption and scattering coefficients are sum of the respective coefficients of the following constituents: pure water, phytoplankton, CDOM (no elastic scattering), and detrius and minerals. 

The spectral absorption coefficient of pure water is taken from measurements. Both spectral absorption coefficients of CDOM and of detritus and minerals are approximated by fitting an exponential curve to the corresponding single-wavelength measurements. Absorption by phytoplankton is approximated by modulating the spectral absorption coefficient of pigment chlorophyll-A with a~varying chlorophyll concentration parameter.

Spectral scattering coefficients of: a) pure water is based on single-wavelength measurements of the scattering by pure water, b) phytoplankton is only related to its concentration in water, c) detritus and minerals are provided in form of measured data. 

Spectral phase function is modeled as sum of the phase functions of each of the constituents: a) pure water, modeled by Einstein-Smoluchowski theory, b) phytoplankton, modeled by isotropic phase function, c) detritus and minerals, modeled based on the Mie theory with weighted two-terms Henyey Greenstein phase function to reproduce backscattering more accurately.
\ifdefined\SHOWEQUATIONS
The model is based on the published data from oceanography studies. Because only elastic scattering is considered, the inelastic part in the  following equations will be omitted.

The values of optical parameters are obtained as the sum of the contributions of pure water and the dissolved particles and particulate matter present in the water, in the following manner:
\begin{align}
    \absorbCoef (\lambda) &= \alpha_w (\lambda) + \sum_i \alpha_i (\lambda) , \\
    \scatterCoef (\lambda) &= \sigma_w (\lambda) + \sum_i \sigma_i (\lambda) , \\
    p(\lambda, \theta) &= \frac{\sigma_w (\lambda)}{\scatterCoef(\lambda)}
                        p_w (\lambda, \theta) +
                        \sum_i \frac{\sigma_i (\lambda)}{\scatterCoef (\lambda)} 
                       p_i (\lambda, \theta) ,
\end{align}
where superscript $w$~stands for the contribution of the pure water (fresh or salty), the superscript $i$~stands for the constituents in the water body, namely CDOM (yellow matter), phytoplankton, and minerals and organic detritus. 

Spectral absorption coefficients of the constituents are modeled in the following way:
\begin{itemize}
    \item $\alpha_w (\lambda)$: is based on tabulated values from Smith and Baker~\cite{SB1981}, shown in Table~\ref{tab:GSMA2008:absorp}.
    \item $\absorbCoef^y (\lambda)$ (CDOM): by fitting an exponential curve of the form (based on data by Bricaud, Morel and Prieur~\cite{BMP1981}):
    \begin{align}
        \alpha_y (\lambda ) = \alpha_y (\lambda_0) e^{-S_y (\lambda - \lambda_0 )} ,
    \end{align}
    where the reference wavelength $\lambda_0 = 440\,nm$, $S_y$~is the slope of the semilogarithmic absorption curve, usually taken as a~constant with a~value of $0.014\,nm^{-1}$, but can also vary geographically and temporally and be wavelength-dependent. The values of absorption at reference wavelength $\lambda_0$ vary in range between $0.01\,m^{-1}$ to $20\,m^{-1}$, as a~function of turbidity.
    \item $\alpha_p$ (phytoplankton): is expressed using the concentration of the chlorophyll $C$~in the water, which usually ranges from $0.01\,mg/m^{3}$ for open water to $100\,mg/m^{3}$:
    \begin{align}
        \alpha_p (\lambda) = C \alpha_{p}^{*} (\lambda) ,
    \end{align}
    where $C$~can be defined either as the concentration of the main pigment \emph{chlorophyll-a} or as the sum of its concentrations and its degradation products, the \emph{pheopigments}. $\alpha_{p}^{*}$ is the \emph{specific spectral absorption coefficient} (the absorption per unit of concentration) for a~particular species of phytoplankton, given in $m^2/mg$. Typical concentrations were provided (based on work of Sathyendranath, Lazzara and Prieur~\cite{SLP87} ), in Table~\ref{tab:GSMA2008:absorp}, or as roughly correspondence between chlorophyll concentrations and several oceanic water types (given by Morel~\cite{Mor1988} (1988)). 
    \item $\absorbCoef^d$ (detritus and minerals): is approximated by an exponential function ( according to Roesler, Perry and Carder~\cite{RPC1989}):
    \begin{align}
        \alpha_d (\lambda) = \alpha_d (\lambda_0) e^{-S_d (\lambda - \lambda_0)} ,
    \end{align}
    where $\lambda_0 = 400\,nm$ is the reference wavelength, $S_d$ is the exponent coefficient with typical values in range from $0.006\,nm^{-1}$ to $0.014\,nm^{-1}$, with $0.011\,nm^{-1}$ being the most common value. Based on~\cite{BSF2003} average slope $S_d$ of $0.0123\,nm^{-1}$ is used.
\end{itemize}

\begin{table}[h]
\begin{tabular}{l|cccccccc}
\hline
$\lambda [nm]$ & 380 & 440 & 500 & 550 & 610 & 670 & 720 & 780    \\ \hline
$\alpha_w$ & 0.00022 & 0.000145 & 0.000257 & 0.000638 & 0.00289 & 0.0043 & 0.01169 & 0.0236 \\
$\alpha_{p}^{*}$ & 0.025 & 0.035 & 0.02 & 0.01 & 0.007 & 0.015 & 0.001 & 0.0001 \\
\hline
\end{tabular}
\caption{ Values of absorption coefficient for clear water body~$\alpha_w$ in $[cm^{-1}]$~and specific absorption coefficient for phytoplankton~$\alpha_{p}^{*}$ in $[m^2 \cdot mg^{-1}]$~\cite{GSMA2008}. }
\label{tab:GSMA2008:absorp}
\end{table}

Scattering coefficients and phase functions of each considered constituent are given by:
\begin{itemize}
    \item pure water: uses phase function defined by \emph{Einstein-Smoluchowski theory} instead of Rayleigh's scattering, due to its higher accuracy:
    \begin{align}
        \sigma_w (\lambda) &= 16.06 \beta_w (\lambda_0, 90^{\circ})
                             \left( \frac{\lambda_0}{\lambda} \right)^{4.32} , \\
        p_w(\theta) &= 0.06225 (1 + 0.835 \cos^2 \theta) ,
    \end{align}
    where typical values for $\beta_w (\lambda_0, 90^{\circ})$ for both fresh and saline pure water are in range from $14.1\cdot 10^{-4}\,m^{-1}$ to $134.5\cdot 10^{-4}\,m^{-1}$ (given in Morel~\cite{Mor1974}, 1974).
    \item CDOM: no elastic scattering is modeled. 
    \item phytoplankton: uses the following relations based on Gordon and Morel~\cite{GM1983}:
        \begin{align}
            \sigma_p (\lambda ) &= \left( \frac{550}{\lambda} \right) 0.30 C^{0.62} , \\
            p_p &= 1/\pi ,
        \end{align}
        where the constant $0.30$~is selected to fit the collected data with the upper bound being at $0.45$ and the phase function $p_p$ is an isotropic function.
    \item organic detritus and minerals: is modeled based on \emph{Mie theory}. To reproduce backscattering more accurately, instead of one term Henyey-Greenstein phase function (HG),  a~\emph{Two-Terms Henyey-Greenstein phase function} (TTHG) is used:
    \begin{align}
        p_d(\theta , a, g_f, g_b) = a p_{HG}(\theta , g_f) + (1-a) p_{HG}(\theta , g_b) ,
    \end{align}
    where $a$~is a~weighting function between zero and one. The first term defines a~forward scattering lobe, with $g_f \in [0 .. 1]$, and the second term defines a~backscattering lobe, with $g_b \in [-1 .. 0]$, and $p_{HG}$ is the Henyey-Greenstein phase function (defined in Equation~\ref{eq:phase-HG}). The typical values for scattering coefficient $\sigma_d$ were provided based on measured data after Stramski et al.~\cite{SBM2001} due to a~wide variety of particulate matter (shown in Table~\ref{tab:GSMA2008:scatter}).
\end{itemize}
\fi
Figure~\ref{fig:GSMA2008-outputs} shows selected renderings with the colors being mainly caused by absorption and elastic scattering.

\ifdefined\SHOWDATATABLES
\begin{table}[h]
\begin{tabular}{ll|cccccccc}
\hline
& $\lambda [nm]$ & 380    & 440    & 500    & 550   & 610    & 670    & 720   & 780   \\ \hline
detritus & $\sigma_{dt} [m^{-1}]$ & 0.045  & 0.0375 & 0.0325 & 0.03  & 0.0285 & 0.0275 & 0.027 & 0.027 \\
minerals & $\sigma_{m} [m^{-1}]$ & 0.0675 & 0.0525 & 0.05   & 0.045 & 0.04   & 0.036  & 0.034 & 0.032 \\
total & $\sigma_{d} [m^{-1}]$ & 0.1125 & 0.09   & 0.0825 & 0.075 & 0.0685 & 0.0635 & 0.061 & 0.059 \\ \hline
\end{tabular}
\caption{Values of scattering coefficients for detritus $\sigma_{dt}$ and minerals $\sigma_{m}$ (after Stramski et al.~\cite{SBM2001})~\cite{GSMA2008}.}
\label{tab:GSMA2008:scatter}
\end{table}
\fi

\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/GSMA2008/no-concent-1.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/GSMA2008/low-phyto-2.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/GSMA2008/low-turbidity-3.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/GSMA2008/low-phyto-turbidity-4.png}
\end{subfigure}
\caption{Renderings with varying concentrations of chlorophyll and the minerals and detritus turbidity, from left to right: no concentration, low concentration of chlorophyll, low turbidity and no concentration, low concentration of both. The source provides values of the simulated parameters.
(\emph{Renderings from Gutierrez et al.~\cite{GSMA2008}}).
}
\label{fig:GSMA2008-outputs}
\end{figure}


\noindent\textbf{Sedlazeck and Koch}~\cite{SK2011} simulate deep sea underwater images based on physical models. The authors focus on deep sea scenarios, where no natural light exists, but the medium is illuminated in a~few meters radius by point lights before the light gets attenuated. The model uses absorption coefficients from measurements and parameterized phase function. The scattering coefficient is computed by integrating the phase function. All the optic parameters are wavelength-dependent. 

The absorption coefficient is computed as the sum of each considered group of matter: water molecules, dissolved salt, CDOM, phytoplankton, and anorganic particles caused by erosion and similar processes. 

They model forward scattering and backward scattering separately. The forward scattering is approximated by a~linear filter (in essence, a~low pass filter). On the other hand, backward scattering involves an explicit use of phase function. The phase function used is a~combination of the Einstein-Smoluchowski model and the \emph{Kopelevich model}. The former is used to describe scattering by pure water (from single-wavelength measurements), and the latter is used to add scattering by large and small particles in the water (also from single-wavelength measurements of their respective phase functions and varying concentrations).
\ifdefined\SHOWEQUATIONS
\textbf{Sedlazeck and Koch~\cite{SK2011}} simulate deep sea underwater images based on physical models. The authors focus on deep sea scenarios, where no natural light exists, but the medium is illuminated in a~few meters radius by point lights before the light gets attenuated. The model uses parameterized absorption coefficient and phase function with all the optic parameters being wavelength-dependent. 

The absorption coefficient is computed as the sum of each considered group of matter: water molecules, dissolved salt, yellow matter (dissolved organic particles), phytoplankton, and anorganic particles caused by erosion and similar processes. 

Scattering is partitioned into forward ($\theta \in [0, \pi/2]$), and backward scattering ($\theta \in [\pi/2, \pi]$). The forward scattering, which happens mostly in the very small angles at the beginning of the interval, is approximated by a~linear filter (a~low pass filter, in essence). For the backward scattering, the phase function is used explicitly and is modeled as a~combination of the \emph{Einstein-Smoluchowski model} and the \emph{Kopelevich model}. The Einstein-Smoluchowski model is used to describe scattering by pure water:
\begin{align}
    \beta_w(\theta, \lambda) = \beta_w(90^{\circ}, \lambda_0)
                                \left( \frac{\lambda_0}{\lambda} \right)^{4.32} ,
\end{align}
where values for $\lambda_0, \beta_w(90^{\circ}, \lambda_0)$ and for the final $\scatterCoef (\lambda)$ can be found in Mobley~\cite{Mob1994} (1994). The Kopelevich model is used to add scattering at large and small particles in the water. Then the final phase function is given by:
\begin{align}
    p(\theta, \lambda) = \beta_w(\theta, \lambda) + 
                         \nu_s \beta_s^* (\theta) \left( \frac{\lambda_0}{\lambda} \right)^{1.7} +
                         \nu_l \beta_l^* (\theta) \left( \frac{\lambda_0}{\lambda} \right)^{0.3} ,
\end{align}
where $\nu_s$ and $\nu_l$ are the volume concentrations of small and large particles, respectively, and $\beta_s^* (\theta)$ and $\beta_l^* (\theta)$ are the phase functions for small and large particles per unit volume concentration, respectively. The concrete values for these terms and for the  parameterization of absorption coefficient are not specified, but the works of Gutierrez et al.~\cite{GSMA2008} and Mobley~\cite{Mob1994} are cited. 
\fi
Figure~\ref{fig:SK2011-outputs} shows renderings of varying turbidity.

\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/SK2011/1-no-turbid.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/SK2011/2-low-turbid.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/SK2011/3-turbid.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/SK2011/4-high-turbid.png}
\end{subfigure}
\caption{Renderings with increasing turbidity from left to right, caused by the increased concentration of scattering particles.
(\emph{Renderings from Sedlazeck et al.~\cite{SK2011}}.)
}
\label{fig:SK2011-outputs}
\end{figure}


% TODO Shi et al.~\cite{SZZW2012} focus on realistic rendering of polluted water

\noindent\textbf{Ma et al.}~\cite{MXWTC2016} employ a~bio-optical model to render ocean color in the \emph{euphotic zone} of the ocean, the uppermost layer that receives the most sunlight. Their model includes constituents: pure water, phytoplankton, CDOM and inorganic suspended particles. They also consider non-uniform vertical distribution of all bio-optical properties based on typical chlorophyll concentration in such waters. The vertical distribution of chlorophyll concentration $C$~is then approximated for each depth according to a~Gaussian function for some initial concentration and depth range.

The spectral absorption coefficient is the sum of absorption coefficient of: a) pure water, taken from wavelength-dependent measurements, b) phytoplankton, computed from spectral chlorophyll-specific absorption coefficient of phytoplankton based on $C$,~c) CDOM, which is derived from $C$~by fitting an exponential curve to the measured input data, similarly to Gutierrez et al.~\cite{GSMA2008}.

The scattering coefficient is the sum of spectral scattering coefficients of pure water, which is obtained from measurements, and of particulate matter, which is computed using $C$. They model two phase functions to account for: a) scattering by pure seawater, which is modeled based on Rayleigh scattering, b) for scattering by particulate matter, which is derived from measured data. The final phase function is then the sum of these two phase functions weighted by the ratio of their corresponding scattering coefficients to the (total) scattering coefficient.
\ifdefined\SHOWEQUATIONS
\textbf{Ma et al.~\cite{MXWTC2016}} employ a~bio-optical model to render ocean color in the \emph{euphotic zone} of the ocean, the uppermost layer that receives the most sunlight (described in Section~\ref{sec:ocean-zones}). Their model includes three constituents: phytoplankton, CDOM and inorganic suspended particles. Consider non-uniform vertical distribution of bio-optical properties linked to chlorophyll concentration for case\,I waters (as per Jerlov classification~\cite{TODOJERLOV}). The vertical distribution of chlorophyll concentration is then approximated according to a~Gaussian function:
\begin{align}
    C(z) = C_0 + \frac{h}{\mu \sqrt{2\pi}} exp \left[
        - \frac{1}{2} \left( \frac{z - z_{max}}{\mu} \right)^{2}
    \right] ,
\end{align}
where $C_0$ is the background chlorophyll concentration, $z_{max}$ is the depth of the maximal chlorophyll concentration, $\mu$~is the standard deviation, and $h$~determines the amplitude of the maximal chlorophyll concentration. 

The spectral absorption coefficient for a~depth $z$~is then computed as the sum of the absorption coefficients of the considered constituents:
\begin{align}
    \absorbCoef (\lambda, z) &= \alpha_w (\lambda) + \alpha_p (\lambda, z) + \alpha_c (\lambda, z) , \\
    \alpha_p (\lambda, z) &= 0.06 \alpha_{c}^{*} (\lambda) C(z)^{0.05} , \\
    \alpha_c (\lambda, z) &= 0.012 C(z)^{0.65} e^{-0.014 (\lambda - 440)} ,
\end{align}
where $w$, $p$, $c$ stand for water (by Pope and Fry~\cite{PF1997}, 1997), phytoplankton and CDOM, respectively, and $\alpha_{c}^{*}$~is chlorophyll-specific absorption coefficient of phytoplankton (by Morel~\cite{Mor1991}, 1991).

The spectral scattering coefficient that varies with depth $z$~is then decomposed in the form:
\begin{align}
    \scatterCoef (\lambda, z) &= \sigma_w (\lambda) + \sigma_p (\lambda, z) , \\
    \sigma_p (\lambda, z) &= \frac{550}{\lambda} 0.3 C(z)^{0.62} , \\
\end{align}
where $\sigma_w$ is the scattering coefficient for pure water, and $\sigma_p$~is the scattering coefficient of particulate matter, both from Smith and Baker~\cite{SB1981}, 1981. The authors plot the coefficients of the constituents for a~choice of parameters ($C_0 = 0.2\,mg\,m^{-3}, h = 144\,mg\,m^{-3}, \mu = 9\,m, z_{max} = 17\,m$) with varying depth.

The used phase function is composed of two terms:
\begin{align}
    p(\theta, \lambda, z) &=
        \frac{\sigma_w (\lambda)}{\scatterCoef(\lambda, z)} \beta_w(\theta) +
        \frac{\sigma_p (\lambda)}{\scatterCoef(\lambda, z)} \beta_p(\theta) , \\
    \beta_w (\theta) &= \frac{3}{4\pi (3 + p)} (1 + p \cos^2 \theta) ,
\end{align}
where $\beta_w$ is the \emph{Rayleigh scattering phase function} for pure seawater (from Morel~\cite{Mor74}, 1974), with $p = 0.835$~as the \emph{polarization factor}, and $\beta_p$ is the \emph{Petzold average-particle function} for particulate matter approximated derived from the measured data by Petzold~\cite{Pet1972}, 1972. 
\fi
Renderings with this model are shown in Figure~\ref{fig:MXWTC2016-outputs}.

\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/MXWTC2016/1-pure-seawater.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/MXWTC2016/2-high-chlorophyll.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/MXWTC2016/3-high-CDOM.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/MXWTC2016/4-high-particulates.png}
\end{subfigure}
\caption{Renderings of ocean color with varying optical parameters, from left to right: pure seawater, seawater with high chlorophyll concentration, seawater with high CDOM concentration, and seawater with high concentration of suspended particulate matter.
(\emph{Renderings from Ma et al.~\cite{MXWTC2016}.})
}
\label{fig:MXWTC2016-outputs}
\end{figure}


% TODO table of comparison between bio optical models, but in what way?
% cite the sources or write the equations?




% ================================================================================
% ================================================================================
\section{Participating Media Rendering}
\label{sec:ocean-media}
%\section{Rendering Participating Media}

This section describes existing approaches that compute light transport in participating media. First, the scope of the techniques is defined. Then, the related approaches are divided into the following groups:
\begin{enumerate}
    \item Approaches that focus on light transport inside the water volume, computed either from above or below the water surface. 
    \item Approaches based on computer vision algorithm.
    \item Approaches that approximate single scattering in any media using volumetric textures.
    \item Approaches that approximate multiple scattering.
\end{enumerate}
This section uses the notation presented in Table~\ref{tab:models-notation}.

There is a~vast amount of techniques concerning participating media rendering in computer graphics, even those techniques that are concerned with rendering of waters~\cite{DCGG2011}. Light transport in volumes is governed by the radiative transfer equation (RTE, described in Section~\ref{sec:light-media}). As the light enters the water volume, it is scattered multiple times inside the water volume by particulates and molecules and quickly attenuated. Some of the scattered light is propagated back through the surface into the sky due to backscattering, contributing to the radiance field in the environment.

Tessendorf~\cite{Tes2001} in his notes discusses this radiosity problem and provides a~schematic form to compute radiance at any point below the water surface:
\begin{align}
    L_{below} = R_T L_D + R_T L_I + L_{SS} + L_M ,
\end{align}
where $R_T$ is the \emph{Fresnel transmissivity} of the incident light at each point and angle on the water surface, $L_D$ and $L_I$ is the direct light from the sun and the indirect light from the atmosphere that penetrates into the water, respectively, $L_{SS}$ is the single-scattered light, from the sun and the atmosphere, that is scattered once in the volume before arriving at any point, and $L_M$ is the multiply-scattered light, i.e., the single-scattered light that undergoes multiple scattering events in the volume. All the terms depend on each other in complex ways: single-scattered light depends on the direct and indirect light, the multiply-scattered light depends on the single scattered.

Darles et al.~\cite{DCGG2011} divides ocean rendering techniques into two categories based on the order of approximation:
\begin{itemize}
    \item First order approximations. These techniques apply geometrical optics laws to render light-water interactions. The water is assumed to be homogeneous and isotropic.
    \item Multiple order approximations. Techniques that simulate the phenomena that occur inside the water volume by simulating complex light-water interactions.
\end{itemize}
The multiple order approximations usually solve the RTE in some form. This field is full of offline approaches that are able to directly simulate global illumination phenomena, such as caustics. Such techniques are, e.g., photon mapping, or path tracing. These techniques are out of scope of this work, the focus is solely on volume scattering approaches for real-time.

\subsection{First Order Approximations}
The first order approximations, as introduced above, often target real-time. The common workflow of these techniques is following: Ray travels from a~light source to an intersection point on a~water surface, where it is reflected and refracted (as described in Section~\ref{sec:air-water-interface}). The refracted ray then travels inside the water volume, either intersecting the bottom, or more commonly, an environment map is sampled using the direction of travel. The attenuation along the ray is computed using Beer's law. Finally the resulting radiance is weighted using Fresnel transmissivity. There might other simplifications or tweaks, such as parallel incident rays, or water surface is represented as a~horizontal plane~\cite{DCGG2011}.

\vspace{1em}

\noindent\textbf{Tessendorf}~\cite{Tes2001} provides a~simplified approach to render ocean surface based on the optical parameters of deep ocean waters. Applies Beer's law to compute the radiance along a~ray as it travels inside the volume accounting for absorption and out-scattering. The out-scattering part is viewed as an approximation of the intensity loss due to multiple scattering: 
\begin{align}
    L_{uw} = \scatterCoef I e^{-\extinctCoef d} e^{-\diffuseAttenDown z}
           = \scatterCoef I e^{-\extinctCoef d - \diffuseAttenDown z} ,
\end{align}
where $I$ is the intensity of the ray, $d$~is the distance traveled in the volume, and $z$ is depth below the surface. The volume is treated as Lambertian reflector with a~single value for its reflectivity.

\vspace{1em}

\noindent\textbf{Jensen and Goliáš}~\cite{JG2001} introduce real-time rendering approach of oceans. Specifically, they focus on rendering surfaces of deep ocean waters. After a~ray of light hits the water surface, a~reflected and refracted directions are computed. The reflected direction is computed along the surface normal and the refracted direction is computed based on Snell's law, with a~fixed index of refraction of water. To render the refraction effect on the water surface, a~cube map of the underwater environment is used. For local reflections and refractions, i.e. when the ray should intersect objects underwater or above the water surface, a~scene is rendered into a~texture as if it were reflected or refracted from a~flat plane placed below or above the surface. The reflected or refracted ray is then intersected with this plane. The intensity along the reflected and refracted ray is weighted based a~precomputed Fresnel term.

\vspace{1em}

\noindent\textbf{Baboud and Décoret}~\cite{BD2006} use the lighting (and bio-optical) model given by Premože et al.~\cite{PA2001} to render water volumes in real-time. They assume a~water volume that is enclosed by terrain and water surface, both represented as heightfields, in a~small aquarium-like bounding box (see Figure~\ref{fig:BD2006-renderings} for the final renderings). Simplified ray-tracing approach is then used to find, first, the intersection of view ray with the water surface, second, the intersection of the refracted view ray with the terrain (as illustrated in Figure~\ref{fig:BD2006-algorithm}). The water surface is intersected only once, multiple intersection, even if they happen (e.g., due to the presence of high waves), are ignored. The terrain is assumed to be diffuse, and its color underwater is computed according to \emph{Phong shading model}. The color of the sky along the reflected ray is sampled from an environment map.
%in the following fashion:
%~view ray intersects the water surface, and is both refracted and reflected. The reflected ray, of the view ray along the surface normal, can hit either the terrain above the surface or the outer bounding-box. The final color along the reflected ray is then the color of the terrain at the corresponding point or the color of the sky sampled using an \emph{environment map}. 
%- The refracted ray, computed using a~fixed index of refraction of water, is intersected with the terrain. If there is no intersection found i.e., the outer bounding-box is hit, then the environment map is sampled (illustrated in Figure~\ref{fig:BD2006-algorithm}).
%- The final color is then computed as the weighted sum of the computed color along the refracted and reflected rays, based on the Fresnel terms. 

% TODO add figure 4 by side by side each
\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/BD2006/alg1.pdf}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/BD2006/alg2.pdf}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/BD2006/alg3.pdf}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/BD2006/alg4.pdf}
\end{subfigure}
\caption{Four possible ray intersections with scene handled by Baboud et al.~\cite{BD2006} (from left to right): hits the scene bounding box and samples the environment map representing the sky; the diffuse terrain is hit; the water surface is hit, while the reflected ray samples the environment map, the refracted ray intersects ocean floor; the reflected ray samples the terrain color. (\emph{Images taken from Baboud et al.~\cite{BD2006}}.)
}
\label{fig:BD2006-algorithm}
\end{figure}


The used lighting model (Premože et al.~\cite{PA2001}) is a~crude simplification of the RTE, which introduces two separate underwater light fields: the \emph{diffuse field radiance} $L_{df}$,~to account for light scattered throughout the media, and the \emph{directional radiance}\,---\,the computed radiance along a~ray. The \emph{irradiance ratio} $S = E_u / E_d$, a~ratio of \emph{upwelling irradiance} to \emph{downwelling irradiance} (an IOP characteristic), is used, along with~$E_d$, to compute the value of $L_{df}$ just below the surface. The ratio $S$~is approximated using optical parameters of the medium, and the downwelling irradiance just below the surface is approximated as a~sum of sun and sky contributions. The uniform water body allows for analytical integration of the attenuation along a~ray inside the water volume, while accounting for absorption with the increasing distance traveled and for out-scattering with the increasing depth.

%Baboud et al~\cite{BD2006} also discuss various techniques and optimizations concerning the intersection algorithm, accessing of the heightfields, precomputation of the attenuation integral to a~2D texture parameterized by depth and distance traveled underwater, or even a~1D texture by ignoring the depth dependence, this simplifies the exponentiation to linear blending operation. Also provide a~way to determine the maximum distance along the ray, after which the contribution of the terrain is negligible, an analogy to ``early-out'' for ray marching.

% TODO figure 11: side by side

\begin{figure}[H]
\begin{subfigure}{.325\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/BD2006/BD2006-output-1.png}
\end{subfigure}
\begin{subfigure}{.325\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/BD2006/BD2006-output-2.png}
\end{subfigure}
\begin{subfigure}{.325\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/BD2006/BD2006-output-3.png}
\end{subfigure}
\caption{Renderings of water volumes as proposed by Baboud et al.~\cite{BD2006}, using the lighting model from Premože et al.~\cite{PA2001}. }
\label{fig:BD2006-renderings}
\end{figure}


\noindent\textbf{Wenzel}~\cite{Wen2006} in notes about real-time atmospheric effects elaborates on ocean shading in CryEngine~\cite{cryengine}. In addition to the homogeneous volume, the water surface is represented as a~horizontal water plane. The final color along a~ray underwater is sum of the attenuated scene color and in-scattering contribution. They assume that the light entering the water always propagates straight down towards the seabed (see Figure~\ref{fig:Wen2007-renderings} for illustration). Then the contribution due to in-scattering along a~ray is analytically computed. The incident radiance from the sky is sampled from sky environment map.

\begin{figure}[H]
\centering
\begin{subfigure}{.325\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/Wen2007/schema.png}
\end{subfigure}
\begin{subfigure}{.325\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/Wen2007/under1.png}
\end{subfigure}
\begin{subfigure}{.325\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/Wen2007/under2.png}
\end{subfigure}
\caption{From left to right: observer underwater with rays propagating straight down from a~surface plane; underwater rendering; the same rendering with higher density. (\emph{Images taken from Wenzel~\cite{Wen2006}}.)
}
\label{fig:Wen2007-renderings}
\end{figure}

\subsection{Computer Vision-based Approaches}

\textbf{Sedlazeck et al.}~\cite{SK2011} simulate both color and geometrical effects of deep sea underwater images using physical models. They focus on deep sea scenarios, where there is no natural light, but only point lights that illuminate only few meters before the light is attenuated. They account for light attenuation, scattering, and refraction.

They simulate underwater housing of a~camera, with a~flat glass port. This way the incident light is being refracted twice as it travels underwater: at the water-glass interface, and at the glass-air interface. They also account for camera transmittance and other effects such as vignetting.

To compute the color underwater they extend \emph{Jaffe-McGlamery model}, a color model used in computer vision algorithms for the support of multiple light sources, rendering shadows, and for the use of a~parameterized volume scattering function. Irradiance incident upon the pixel sensor of the camera, using the color model, is expressed as the sum of three irradiance components (shown in Figure~\ref{fig:SK2011-renderings}):
\begin{align}
    E_{total} = E_{direct} + E_{forwardScatter} + E_{backScatter} .
\end{align}
Direct irradiance $E_{direct}$, accounts for point light sources immersed in the media illuminate nearby objects and for the reflected light from visible objects, that are not shadowed by any other objects. The incident ray enters the camera while being attenuated due to out-scattering and absorption. Light being scattered forward at small angles is approximated with a~distance dependent Gaussian filter.
%, where - filter mask depends on the distance he light traveled - and the amount of forward scatter added is weighted by an empirical factor.

The Forward scattered irradiance is again approximated using the distance dependent Gaussian filter of $E_{direct}$.

The backscatter (or \emph{veilling light}) irradiance approximates the light being scattered in all directions multiple times, causing ambient or veiling light close to a~light source. This is approximated by slicing the volume into a~number of planes of certain thickness that are parallel to the camera's near plane. Each such plane, has the irradiance incident upon it computed similar to $E_{direction}$, now for each position on a~plane with the volume scattering function applied, and forward scattered part is also taken into account. Then the contribution due to backscatter is computed as the superposition of all the slices.

% TODO image of camera from fig 1. , side by side, with three outputs from fig 6.

\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/SK2011/alg1-scene.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/SK2011/alg2-direct.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/SK2011/alg3-backscatter.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/SK2011/alg4-combined.png}
\end{subfigure}
\caption{Renderings, form left to right: scene, direct signal, backscatter contribution, combined result. (\emph{Renderings from Sedlazeck et al.~\cite{SK2011}}.)
}
\label{fig:SK2011-renderings}
\end{figure}



\subsection{Single Scattering Approximations}
The following approaches render participating media using voxelized volumes. Though these approaches are used to compute atmospheric scattering (as shown in Figure~\ref{fig:wronski-output}), they provide a~general outline of how to compute single scattering in any media.

\begin{figure}[H]
    \centering
\begin{subfigure}{.6\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/partic/wronski.png}
\end{subfigure}
\caption{Rendered scene with atmospheric scattering using the participating media rendering approach by Wronski~\cite{Wro2014}.
}
\label{fig:wronski-output}
\end{figure}


% \textbf{Elek et al.~\cite{ERDS2014}} rather not ...

%Improves over on Wenzel's~\cite{Wen2007} approach to atmospheric effects in real-time
%- based on Kaplanyan Light Propagation volumes,

\noindent\textbf{Wronski}~\cite{Wro2014} presents unified compute shader-based approach for rendering atmospheric scattering as participating media, called \emph{volumetric fog}. Media have varying density and shadowing term is added, which causes light shafts or godrays. Ray marching is used to render voxelized participating media into a~volumetric texture that is positioned in the view frustum, so-called \emph{frustum aligned 3D texture}. 

The approach allows for parallelization of ray marching steps by storing intermediate results into volumetric textures between multiple passes:
\begin{enumerate}
    \item Density estimation: generates density (attenuation) and scattering coefficient, stores the former in RGB channels and the latter in alpha channel.
    \item Lighting calculation: stores the constant ambient terms from the main light and from multiple point lights that intersect the view frustum. Both are modulated by phase function and stored in RGB.
    \item Ray marching: solves the scattering by linearly marching through the intermediate 3D textures. Applies Beer-Lambert law to compute transmittance which is applied on the in-scattered light and accumulated. Transmittance is stored as a~single value (wavelength-independent) in the alpha channel.
\end{enumerate}
The computed lighting can be applied in forward or deferred shading pipeline by a~simple 3D texture bilinear look-up. The 3D texture is mapped to camera frustum in NDC space with exponential distribution of depth slices. The volume used is $160 \times 90$ with 64 or 128 depth slices, this can produce the effect in distances in world space between 5 and 128\,meters.

% TODO figure

This approach also has limitations and produces artefacts. It lacks high frequency geometry details and overall produces a~rather soft effect, which is suitable for fog. The use of a~discretized volume inherently introduces temporal aliasing and flickering. The author tries to mitigate this by storing only low frequency information for every depth along the ray, and uses uses quadrilinear filtering on the volumetric data to minimize edge artifacts on depth discontinuities.

To mitigate the temporal aliasing or flickering problem of high frequency source information, low-pass filter is used. Also lower resolution and smoother shadow maps are preferred, which shadow leaking being acceptable in participating media.
 %   - and exponential shadow maps algorithm for low detail, smooth shadow maps
 %   - shadow leaking is acceptable in participating media
 %   - downsample and blurred

To mitigate the aliasing even further, especially with animated participating media, various temporal anti-aliasing techniques (e.g., termporal super-sampling) can be used. Instead of regular grid sampling, which produces aliasing, jittered grid sampling can be used. It reduces aliasing but introduces higher frequency noise. This higher frequency signal can then be filtered and blurred using lower resolution kernels (in both temporal and spatial domains).

\vspace{1em}

\noindent\textbf{Hillaire}, 2015~\cite{Hil2015} presents physically based and unified approach to volumetric rendering. Focuses on atmospheric scattering and fog with varying density and scattering occlusion to simulate light shafts. %Allows for volumetric regions of space to be added to scene, with the volumetric lighting applied on to the objects it interacts with, this causes volumetric shadows, light shafts.
%while allowing author tailored materials for the participating media

%Compute single scattering inside 
Uses frustum aligned 3D textures in clip space, similarly to Wronski~\cite{Wro2014}, with the same resolution as screen and with 64\,depth slices. Assumes a~tiled-based deferred pipeline which uses 16x16 tiles with culled light lists. The approach brings the following improvements to the each corresponding step of Wronski's algorithm:
\begin{enumerate}
    \item Material properties: participating media properties are stored into so-called \emph{volume buffer} or V-Buffer, which allows for heterogeneous media: stores wavelength-dependent scattering coefficient and single value for extinction (wavelength-independent). Also stores $g$~parameter for phase functions.
    \item Froxel light scattering: into a single 3D texture per voxel evaluates in-scattered light and extinction, considers indirect light, sun light and local lights, Takes into account shadows, and volumetric shadows for local lights. Introduces \emph{temporal volumetric integration} to mitigate aliasing and flickering. This is accomplished by jittering of samples by the same offset for both the material and scattering samples along the view ray. Then, 
    %similarly to the temporal anti-aliasing technique [Kar2014],
    the current frame is blended with the previous frames, based on exponential moving average.
    \item Final integration: for each voxel along view ray, solves the single scattering by ray marching and integrating. Introduces energy conservative integration, by analytically integrating the scattered light as a~function of transmittance, to mitigate light leaking due to highly scattering participating media thus improving volumetric shadows.
\end{enumerate}
Animated participating media (e.g., fluid simulation) or moving lights can cause artifacts by leaving trails behind.

% TODO images of it being applied in a~scene

% Should be a subsection? or just put it here
%\subsection{Multiple-scattering approximations}

% Iwasaki et al. second order scattering
% Wrenninge et al. and Hil2016
% Monzon et al.

\subsection{Multiple Scattering Approximations}

Computation of multiple scattering is beyond what is feasible in real-time. Most techniques that focus on real-time thus must assume crude simplifications. Only few relevant and simplest are investigated.

\vspace{1em}

\noindent\textbf{Hillaire}, 2016~\cite{Hil2016} approximates multiple scattering in clouds by summing multiple octaves of single scattering, while increasing the properties of participating media with each octave. It can be easily incorporated to ray marching, but it does not represent well complex behavior caused by multiple scattering: side or backward scattering, cone spread, etc.

\vspace{1em}

\noindent\textbf{Monzon et al.}~\cite{MAGM2022} approximate multiple scattering in underwater scenarios based on the commonly used AOP (apparent optical property)\,---\,\emph{diffuse attenuation coefficient for downwelling irradiance}~$\diffuseAttenDown$ (described in more detail in Section~\ref{sec:iops-aops} about IOPs and AOPs). They build on the assumption that multiple scattering underwater is of very low frequency and can therefore be modeled as a~function of depth and wavelength. They further assume non-emissive isotropic medium. This simplification, along with the assumption of constant incident irradiance (radiance from every direction) at every distance from water surface, allows for analytical evaluation of the volume rendering equation within the medium. Figure~\ref{fig:magm2022-output} shows renderings for different values of optical parameters in different waters based on the Jerlov classification (described in Section~\ref{sec:jerlov-waters}).

\begin{figure}[H]
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/MAGM2022/Canon_10D_JerlovIB.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/MAGM2022/Canon_10D_2_Jerlov1C.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/MAGM2022/Canon_10D_2_JerlovIII.png}
\end{subfigure}
\begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/MAGM2022/Canon_10D_Jerlov3C.png}
\end{subfigure}
\caption{Rendered scene with different set of spectral optical parameters for different Jerlov water types. (\emph{Images from Monzon et al.~\cite{MAGM2022}}.)}
\label{fig:magm2022-output}
\end{figure}


% ==============================================================================
%\subsection{Ray Marching}
%-> cites Perlin Hypertexture 1989

%- 1450: Quı́lez, Íñigo, “Rendering Worlds with Two Triangles with Ray Tracing on the GPU in 4096 bytes,” NVScene, Aug. 2008. Cited on p. 454, 594, 752


% ================================================================================
%\section{Approaches for Rendering Underwater Caustics}

% ================================================================================
%\section{Sky Models}

% ================================================================================
\section{Procedural Noise Functions and Implicit Surfaces}
\label{sec:proc-gen}

Procedural generation is an alternative way of generating data or simply assets to the typical human or manual generation of assets. The term \emph{procedural} refers to the description of entities by algorithms and program code rather than by data structures. Its goal is to approximate the seemingly random and unpredictable variety observed in real world by automatically generating various kinds of assets based on parameters and through the use of randomness or pseudo-randomness. There is a~wide range of techniques used for generation of textures, models, and animation~\cite{TMPA2002}.

The focus of this section is to provide an overview of noise functions, their combination, and can they be used to generate or model terrains represented as implicit surfaces. The last part provides a~brief discussion about visualization of such representations. Noise functions can be used to generate, animate or aid in modeling of water surfaces, or generally any surface or terrain, such as the ocean floor~\cite{TMPA2002}.

\subsection{Noise Functions}

The building block of procedural generation is noise, which is usually computed by interpolating between pseudo-random values. There are certain requirements or properties for a~function to be noise that is useful for generation, it should not be an incoherent, completely random noise, such as white noise, because true randomness is an unusual and undesirable phenomenon. The properties of commonly used noise functions usually are~\cite{TMPA2002}:
\begin{itemize}
    \item it produces ``apparently'' random, but not entirely random values like white noise,
    \item it is repeatable, i.e., it gives the same value for the same inputs (e.g., texture coordinates),
    \item it is smooth enough for it to be stable and not show aliasing artifacts,
    \item its values are in range from -1 to 1,
    \item it does not exhibit obvious periodicities or regular patterns.
\end{itemize}

Commonly used noise functions range from one-dimensional to three-dimensional. There are many categories of noise functions based on the way they are computed. One such category is \emph{Lattice noises}. These functions use pseudo-random numbers that are distributed uniformly in a~lattice (a~line for 1D, a~grid for 2D, a~3D grid). \emph{Value noise} is a~type of noise that interpolates between the corner values of each cell to obtain noise. This is done in 1D, by interpolation between two points on a~line, or in 2D, by interpolation between the four corners of a~cell of the grid, and also in 3D, by interpolation between eight corners of a~cube~\cite{GVL2015}. 

Interpolation functions are usually polynomials of a~certain degree. The higher the degree, the smoother values are obtained, but also the computational cost rises. Most commonly used functions  are: linear, and cubic~\cite{GVL2015}. 

Because value noise tends to look ``blocky'', another type, called \emph{Gradient noise}, is oftentimes used instead. Typical example of gradient noise is \emph{Perlin noise}~\cite{Per1985}. In this case, noise is determined at a~point by computing a~pseudo-random gradient at each of the eight nearest vertices on the integer cubic lattice and then smoothly interpolating between them. Perlin noise is used widely in graphics because of its simplicity and performance cost.

% TODO cite https://thebookofshaders.com/11/
% TODO images - value nosie vs gradient noise

% 2D value noise by https://thebookofshaders.com/edit.php#11/2d-vnoise.frag
% 2D gradient noise by https://thebookofshaders.com/edit.php#11/2d-gnoise.frag
% by Quilez www.iquilezles.org

%\subsection{Fractal Brownian Motion}
% https://iquilezles.org/articles/fbm/

A~technique called \emph{fractal Brownian Motion} (fBM), or simply ``fractal noise'' is used to produce noise that contains many incrementally finer details, as observed in nature. This is achieved by adding different iterations of noise, called \emph{octaves}, with successively higher frequencies (\emph{lacunarity}) and smaller amplitudes (\emph{gain}). The simplest variant of fBM of a~dimensional noise is:
\begin{align}
    fBM(x) = \sum_i^{octaves} A_i\cdot noise(x\cdot f_i) ,
\end{align}
where $A_i = A_{i-1} \cdot gain$ is the amplitude and $f_i = f_{i-1} \cdot lacunarity$ is the frequency for $i$-th octave, respectively~\cite{GVL2015}.
% TODO cite https://thebookofshaders.com/13/

\subsection{Implicit Surfaces}
Generally, points on a~surface are described explicitly as a~list of concrete parameters, such as $x$, $y$, and $z$. \emph{Implicit surfaces} provide another way of representing surfaces. Instead of describing a~point on the surface using concrete parameters, an implicit form is used instead:
\begin{align}
    f(x, y, z) = f(\*{p}) = 0 .
    \label{eq:implicit}
\end{align}
A point lies on the implicit surface if the function evaluates to zero. In other cases, it gives the signed distance to the implicit surface: negative distance means that $\*{p}$ is inside the surface, and positive outside. Therefore, these functions are sometimes also called \emph{signed distance functions} (SDFs).

Surfaces in their implicit form are often used for intersection testing with rays, as it is often simpler to find the intersection in this form than in the the corresponding (if any) explicit form. Implicit surfaces can be stored in  a~compact way and evaluated on the fly without increasing memory footprint. Objects described in this form can be easily combined in various ways, such as blending, subtracting, adding, etc~\cite{RTR2018:ch17}.

Quilez~\cite{Qui2013-sdf2d}~\cite{Qui2013-sdf} provides an extensive list of both 2D and 3D SDFs for all kinds of primitives. Some simple examples of implicit surfaces are:
\begin{align}
    f_{sphere}(\*{p}, r) &= ||\*{p}|| - r, \\
    f_{plane}(\*{p})     &= p_y, 
\end{align}
where $f_{sphere}$ is an implicit function for a~sphere at origin with radius $r$, and $f_{plane}$ returns the vertical distance (on $y$-axis) to a~horizontal plane. 
% https://iquilezles.org/articles/distfunctions2d/
% https://iquilezles.org/articles/distfunctions/

% TODO image illustrating signed distance fields
% https://iquilezles.org/articles/distfunctions2d/

The normal of an implicit surface is described by the partial derivatives of the function:
\begin{align}
    \*{n} = \nabla f (x, y, z) = \left(
        \frac{\text{d}f}{\text{d}x},
        \frac{\text{d}f}{\text{d}y},
        \frac{\text{d}f}{\text{d}z}
    \right) ,
\end{align}
where $\nabla f$~denotes the gradient of the function. To compute the normal using this equation, i.e., analytically, the function must be differentiable, and thus also continuous. More commonly, a~numerical method called \emph{central differences} is used instead. To approximate the normal, the implicit function is sampled in each direction in small steps:
\begin{align}
    n_x = \nabla f_x \approx 
        f(\*{p} + \epsilon \*{e}_x) - f(\*{p} - \epsilon \*{e}_x) , \\
    n_y = \nabla f_y \approx
        f(\*{p} + \epsilon \*{e}_y) - f(\*{p} - \epsilon \*{e}_y) , \\
    n_z = \nabla f_z \approx
        f(\*{p} + \epsilon \*{e}_z) - f(\*{p} - \epsilon \*{e}_z) ,
\end{align}
where $\*{e}_x = (1,0,0)$, $\*{e}_y = (0,1,0)$, $\*{e}_z = (0,0,1)$, and $\epsilon$ is a~small number representing the step size~\cite{RTR2018:ch17}.

Primitives can be easily combined together to create whole scenes. Quilez~\cite{Qui2013-sdf} defines all kinds of operators for combining two primitives together into one primitive. For instance, union, subtraction, intersection and their smooth variants, as shown in Figure~\ref{fig:sdf-operators}. Primitives can also be repeated over a~region of space using so-called \emph{Domain repetition operator} that uses modulo operator. Furthermore, primitives can be easily distorted, deformed or blended together.

\begin{figure}[H]
\centering
\begin{subfigure}{.65\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/procGen/sdf-operators.png}
\end{subfigure}
\begin{subfigure}{.34\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/procGen/sdf-operators-smooth.png}
\end{subfigure}
\caption{Images illustrate how primitives represented as SDFs can be combined together using various operators. The left image shows two primitives (box and sphere) and the results of operators (from left to right): union, subtraction, intersection. The right image shows the same primitives repeated with the smooth variants of the same operators applied. (\emph{Images from Quilez~\cite{Qui2013-sdf}.})
}
\label{fig:sdf-operators}
\end{figure}

% https://iquilezles.org/articles/distfunctions/
% TODO figure fig:sdf-operators: union subtract, intersection and smooth variants 
% domain repetation image of grass blades

To visualize a~set of implicit functions, usually one of the two methods are employed: \emph{ray marching} or \emph{polygonalization}. Polygonalization converts implicit surfaces into polygon meshes using algorithms, such as \emph{marching cubes} or \emph{marching tetrahedra}. These meshes can then be rasterized in a~traditional way~\cite{RTR2018:ch17}. Ray marching, on other hand, directly  samples the implicit functions along rays, computing the final color based on ``what is seen''. Ray marching is described in more detail in the section below.

\subsection{Ray Marching}

Ray marching is used to find an intersection of ray with surface. It provides an alternative to analytic ray-intersection methods and rasterization with a~depth buffer. In comparison to ray-intersection methods, which can be used only for primitives with analytic solutions to ray-intersection, ray marching can be thought of as an algorithm for numerical estimation of arbitrary ray-surface intersections. 

Generally, ray marching works by ``marching'' from the ray origin $\*{o}$ along the ray by sampling points $\*{p} = \*{o} + t \omega$ using the ray parameter $t$~as a~distance from the origin. Ray marching terminates when the sampled point intersects or crosses a~surface. It is usually employed to find intersections with implicit surfaces. In this case, the marching terminates when the implicit function evaluates to zero or less.

Marching can be done with fixed steps from the ray origin, which might be too slow or miss an intersection all together, or it can use a~signed distance estimator to find a~safe step size it can march to avoid overstepping the intersection. The safe step size is usually determined by finding a~sphere originating at point $\*{p}$ that touches (i.e., is tangent to) at most one primitive. The safe step size is then the radius of the sphere, since the closest intersection is found at the surface of the sphere. This approach is called \emph{sphere-tracing}~\cite{McG2021}. Figure~\ref{fig:sphere-tracing} illustrates how an intersection is found in a~scene with implicit surfaces using this approach. 

% Apart from the ability to find the closest intersection with the ray, it can be used to efficiently approximate effects such as soft shadows, ambient occlusion or depth of field, because of the way it tracks distances to surfaces in a~scene.

% https://graphicscodex.courses.nvidia.com/app.html?page=about
% TODO https://graphicscodex.courses.nvidia.com/app.html?page=_rn_rayMrch

\begin{figure}[H]
    \centering
\begin{subfigure}{.8\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/procGen/sphere-tracing.pdf}
\end{subfigure}
\caption{Finding an intersection with signed distance fields along a~ray using sphere tracing. The dashed circles indicate a~distance to the closest surface from their centers that is  that is safe to advance along a ray. (\emph{Image from Real-time Rendering, 4th edition~\cite{RTR2018:ch17}}.)
}
\label{fig:sphere-tracing}
\end{figure}


%Quilez~\cite{Qui2008}
% https://iquilezles.org/articles/raymarchingdf/


%\subsection{Rendering Procedurally Generated Terrains}

% https://iquilezles.org/articles/terrainmarching/
% https://iquilezles.org/articles/fbmsdf/

% \subsection{SDFs with Fractal Brownian Noise}
% https://iquilezles.org/articles/fbmsdf/



% ================================================================================
% ================================================================================
% ================================================================================
% ================================================================================


\chapter{Real-Time Caustics Approaches}
\label{cha:caustics}
The interaction of light with curved reflective or refractive surfaces causes the light to defocus in some areas and focus in others. While the defocused areas create shadows, the focused or more dense areas of light cause a~visual phenomenon called \emph{caustics}. For instance, the reflection of light by a~coffee mug causes caustics on the inside and on the outside. Or in a~more relevant case, light incident onto a~curved water surface is refracted and attenuated as it travels through the water volume finally converging on opaque surfaces causing caustics.

This is shown in Figure~\ref{fig:caustics-examples}. In computer graphics, real-time especially, this phenomena is divided into two separate effects: \emph{surface caustics} (image on the left) and \emph{volumetric caustics} (image on the right). Furthermore caustics that are caused by reflection are called \emph{reflected} or \emph{reflective} and caustics caused by refraction are called \emph{refracted} or \emph{refractive}.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.8\linewidth]{figures/caustics/pexels-photo-MIRTO-KON-5250849.jpeg}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{figures/caustics/pexels-jeremy-bishop-2397651.jpg}
\end{subfigure}
\caption{The left image shows surface caustics formed on the shallow floor (by MIRTO KON, pexels.com). The right image shows rays of light underwater, called volumetric caustics (by Jeremy Bishop, pexels.com). 
}
\label{fig:caustics-examples}
\end{figure}

Rendering caustics, especially in participating media, is a~complex task that requires global information, and is therefore generally constrained to offline global illumination techniques, i.e., techniques that solve the volume rendering equation in some form (described in Section~\ref{sec:light-media}). Exact solutions are feasible only in the simplest environments. Renderers, generally, approximate this using numerical methods, such as bidirectional path tracing~\cite{PJH2016} or photon mapping~\cite{Jen2001}.

Real-time rendering methods either fake the caustics altogether, or they try to produce physically-based caustics that correspond to the conditions of the scene (i.e., state of lighting, water surface, etc.). Techniques of the latter type are commonly based on photon mapping. The common workflow of photon mapping is following~\cite{Kim2019}:
\begin{enumerate}
    \item \emph{photon emission}: photons are emitted from light sources into the scene,
    \item \emph{photon tracing}: photons are traced in the scene and scattered by reflective or refractive surfaces (also called reflectors and refractors), until they hit surfaces that are opaque enough, so-called caustics receivers.
    \item \emph{photon gathering}: millions of photons are gathered as they converge at certain places to obtain smooth caustics in \emph{photon map}.
\end{enumerate}
The use of photon mapping is not feasible for real-time because it implies heavy use of ray tracing and it requires a~large number of photons (even tens of million) to obtain smooth images, otherwise the photon map is very noisy. Therefore, the photon mapping-based methods address these issues in various ways:
\begin{itemize}
    \item intersection testing: from crude simplifications of surfaces by simple shapes (e.g., ocean bottom represented by a~plane), to image-space tracing (e.g., ray marching depth buffers), to recent utilization of specialized ray tracing pipelines,
    \item photon gathering: additive blending, usually accompanied by filtering, splatting photons of various sizes (not just point primitives), beam tracing (i.e., tracing triangular beams instead of points), and recently denoising (i.e., blurring).
\end{itemize}

Section~\ref{sec:caustics-surf} provides an overview of certain techniques used for rendering surface caustics, while Section~\ref{sec:caustics-vol} describes certain techniques that render volumetric caustics in media. %Both categories of real-time techniques often employ, in some form, a~simplified approach approach to photon-mapping called \emph{caustics mapping}. 


% top: mention photon mapping and cite{} 
    % gather and scatter
    %- the goal is to gather the photons recorded in the texture
% caustics mapping: compute caustics using a~two-pass photon-mapping like algorithm
% TODO to top
% based on caustics mapping
%   1) construction of the caustics map texture, involves intersection
%   2) application of the caustics map on (diffuse) receivers
%   caustic-forming object
% add the figure 3 of showing how the focus and defocus
%   a specular bounce, or just a~bounce
% reflective and refractive caustics
% must find an intersection with the receiving geometry: ray traced vs. estimate vs. spec.
% point splatting vs. beam tracing - projecting -areas
%   suffers from aliasng 

%\subsection{Caustics Mapping Approaches}
%caustics mapping 
%- similar to shadow mapping
%- photons points are plat into a~light-space texture 
%- sampled by repojection
%- same limitations
%- caustics are rendered only if they are visible to the light source
%
%- first render a \emph{Caustics map}, an image with the refractive/reflective geometry
%- then from it are rays generated as opposed to photon mapping, which directly generates refracted rays on the triangles of the surface
%
%% 4th set of techniques
%- voxelize the receiving geometry and march the voxelized 3D grid
%    - can be an expensive operation, especially when done each frame
%    - poses high implementation complexity
%
%% TODO figure from Gru2019 1.4 remade in draw io

% ------------------------------------------------------------------
\section{Surface Caustics}
\label{sec:caustics-surf}

The techniques presented are divided into the following categories:
\begin{itemize}
    \item Fake caustics: approaches that produce physically incorrect, yet convincingly looking caustics,
    \item Image-space caustics: techniques that trace rays in image space or use other simplifications to estimate intersections with scene geometry, e.g., by marching through depth buffer from the light's point of view,
    \item Ray-traced caustics: recent techniques that intersect rays with  the dynamic scene in world space using DirectX Raytracing (DXR) pipeline, that is supported on certain hardware.
\end{itemize}

\subsection{Fake Caustics}
\label{sec:caustics-fake}

These approaches are usually early attempts at real-time caustics simulation. Caustics are simulated by faking or approximating the observed look. These techniques are purely aesthetics-driven and do not consider realism, i.e., the caustics do not correspond to the real light transport in a~given water volume. Caustic patterns, possible generated offline or created by artists, are stored in textures, which are then projected on surfaces and applied as light maps. These caustics textures might be even animated by using a~set of such textures. This approach has been used in games before~\cite{RTR2018:ch14}.

\vspace{1em}

% Jos Stam
% Stam, Jos. 1996. "Random Caustics: Wave Theory and Natural Textures Revisited." Technical sketch. In Proceedings of SIGGRAPH 1996, p. 151. Available online at http://www.dgp.toronto.edu/people/stam/INRIA/caustics.html
\noindent\textbf{Stam}~\cite{Sta1996} computes a~caustic-like patterns using wave theory (see Figure~\ref{fig:Sta1996-output}). The resulting patterns are periodic both in time and space and are generated for a~ground plane at various depths. Such patterns can be used in form of a~animated texture that is then additively blended with an object's base textures.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.5\linewidth]{figures/Sta1996/caustic_blur.png}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.5\linewidth]{figures/Sta1996/output.jpg}
\end{subfigure}
\caption{Caustics patterns created using Stam's approach~\cite{Sta1996} and stored in a~texture (on the left, \emph{image from Stam~\cite{Sta1996}}), which is then applied onto objects as a~light map (on the right, \emph{image from Guardado et al.~\cite{GSC2004}}).}
\label{fig:Sta1996-output}
\end{figure}

\noindent\textbf{Guardado et al.}~\cite{GSC2004} present purely aesthetics-drive approach for the computation of underwater caustics. They consider only a~subset of arriving rays by making the following assumptions: the sun is directly above, rays hit the ocean floor from directly above, light attenuation is fixed per unit distance, and the ocean floor is represented by a~plane. The approach is based on simplification of the backward ray tracing, i.e., each visible point on the ground represents a~ray that is vertically projected to a~point on the water surface mesh above it (see Figure~\ref{fig:GSC2004-output}). The ray is then refracted by the surface and the intensity along the ray is derived based on how close the its direction is to the vertical direction (sun direction) and how far it has traveled to the surface. The resulting caustics mesh (per-vertex) or texture (per-pixel) is rendered using additive blending on top of the ocean floor (shown in the rightmost image in Figure~\ref{fig:GSC2004-output}).

\begin{figure}[H]
\centering
\begin{subfigure}{.45\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/GSC2004/schematic.jpg}
\end{subfigure}
\begin{subfigure}{.25\textwidth}
    \centering
    \includegraphics[width=0.955\linewidth]{figures/GSC2004/lightMap.jpg}
\end{subfigure}
\begin{subfigure}{.25\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/GSC2004/output-comparison.png}
\end{subfigure}
\caption{The leftmost image illustrates how points on the ground plane are vertically projected to the points on the water surface represented by a~wave equation, are refracted and sample a~precomputed light map (the image in the middle) that corresponds to sun shining from directly above. The rightmost image then shows the rendered caustics. (\emph{Images from Guardado et al.~\cite{GSC2004}}.)
\label{fig:GSC2004-output}
}
\end{figure}

% explain at the top:
%   first order, one bounce
%   beam-tracing approach Watt^Watt
% specular surface -> specular triangle
% caustic mesh approaches

\subsection{Image-space Techniques}

These techniques usually work by tracing photons in image-space or tracing in world space but only with simplified scene, such as a~horizontal plane or other assumptions. These approaches work similarly to shadow mapping: a~scene is rendered from the light's point of view into a~texture called \emph{caustics map}. This texture consists of the visible refractive or reflective geometry. The rays are then generated by sampling the textures and the approximate intersections of the rays with the scene receivers are found. The intersection is usually found by marching the shadow map depth buffer (from the light) or the camera depth buffer (from the camera). 

The techniques suffer from the same visibility limitations as shadow mapping (ignoring issues with resolution): only caustics that are visible from the light are drawn, furthermore it cannot account for reflected/refracted rays that are occluded in both the camera view and the light view.

\vspace{1em}

\noindent\textbf{Jensen et al.}~\cite{JG2001} compute caustics for a~water surface tile represented as a~mesh. They consider only first order rays and a~diffuse surface at a~constant depth. They use beam-tracing approach, illustrated in Figure~\ref{fig:JG2001-output}. For each specular triangle of the water surface a~caustic triangle at the bottom is computed by intersecting refracted rays with the ground plane. These caustics triangles then form a~\emph{caustic mesh} that is drawn into a~texture . The final intensity of a~caustic triangle is based on the ratio of areas of the specular triangle and the caustic triangle. They also make the caustics texture tileable, which also results in an applied anti-aliasing to an effect of a~2\,x\,2 super~sampling. The texture's UV coordinates are computed by projecting the texture from the height of the water surface in the direction of a~ray from the sun, similarly to a~parallel projection. Finally, the sampled intensity is multiplied by dot product between the surface's normal and direction to the sun.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{figures/JG2001/causticMesh.png}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{figures/JG2001/causticsTexture.png}
\end{subfigure}
\caption{The image on the left shows caustics triangles (i.e., triangles at the bottom plane) created by intersecting the refracted rays of each vertex of surface triangles (i.e., on the water surface). The image on right shows the final seamless 256x256 caustics texture with applied 2x2\,anti-aliasing. (\emph{Images from Jensen et al.~\cite{JG2001}}.) 
}
\label{fig:JG2001-output}
\end{figure}

\noindent\textbf{Baboud et al.}~\cite{BD2006} compute caustics using a~two-pass photon-mapping like algorithm that considers only one specular bounce and involves both CPU and GPU.
%They assume a~water volume that is enclosed by a~ground and a~water surface, both represented as height fields, in a~aquarium like manner. 
The first pass renders the water surface from the light source into a~photon texture (illustrated in the leftmost image of Figure~\ref{fig:BD2006-caustics}). Each texel of the texture records the coordinates of the intersection of a~light ray (after being refracted by the water surface) with the ground 
%(uses a~specialized ray tracing approach), 
, and a~photon contribution. The contribution is computed based on the Fresnel transmittance, the traveled distance, the incident angle and a~ray sampling height. The photon texture is then transferred to the CPU-side to construct an \emph{illumination texture} by accumulating at each sample the stored intensity to the position that is stored at that sample. The texture is then transferred back to the GPU where it is used for lighting in the final render pass. 
%The approach can be combined with shadow mapping to also account for shadows of other objects in the scene onto the ground.

Based on the resolution of the two textures (the photon texture and the final caustics texture) the caustics can be very noisy and exhibit aliasing. To reduce aliasing a~random jittering to the light rays can be applied (as illustrated in the Figure~\ref{fig:BD2006-caustics}). The caustics texture is also blurred using an \emph{edge preserving filter} (a~GPU-based \emph{bilateral filter}) to preserve the natural high-frequencies (as shown in the two rightmost images in Figure~\ref{fig:BD2006-caustics}).

\begin{figure}[H]
\centering
\begin{subfigure}{.45\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/BD2006/caustics-alg.pdf}
\end{subfigure}
\begin{subfigure}{.25\textwidth}
    \centering
    \includegraphics[width=0.955\linewidth]{figures/BD2006/caustics-raw.png}
\end{subfigure}
\begin{subfigure}{.25\textwidth}
    \centering
    \includegraphics[width=0.955\linewidth]{figures/BD2006/caustics-denoised.png}
\end{subfigure}
\caption{The leftmost image shows photons emitted from jittered pixel positions of the light frustum through the water surface and onto the ground. The last two images show, from left to right: raw caustics texture, its filtered variant.
(\emph{Images from Baboud et al.~\cite{BD2006}}.)
}
\label{fig:BD2006-caustics}
\end{figure}

\noindent\textbf{Shah et al.}~\cite{SKP2007} present a~real-time technique for rendering caustics from reflective and refractive objects, called \emph{caustics mapping}. The algorithm allows for fully dynamic non-planar geometries. First, from light's view, the receiver geometry is rendered into a~positions texture, then the caustic-forming objects are rendered to generate refracted rays at their vertices. Each such ray finds an estimate intersection with the receiving geometry by interatively sampling the positions texture. Points are splatted into the caustics map by accumulating the computed light intensity using additive alpha blending. The light intensity is estimated from the total flux contribution for each vertex of the caustic-forming object's mesh based on the area of triangles around the vertex to the total area of the mesh. They also account for the absorption as the ray travels through a~refractive object.

The splatting of point primitives leaves gaps in between. Either the geometry of the caustics-forming objects must be further tessellated and/or the caustics map must be low pass filtered (see the Figure~\ref{fig:SP2007-caustics}). Also suffers from problems similar to shadow mapping, such as: fixed resolution of the caustics map, visibility issues, temporal incoherency, etc.

\begin{figure}[H]
\centering
\begin{subfigure}{.325\textwidth}
    \centering
    \includegraphics[width=.85\linewidth]{figures/SP2007/causticsMap-raw.png}
\end{subfigure}
\begin{subfigure}{.325\textwidth}
    \centering
    \includegraphics[width=.85\linewidth]{figures/SP2007/causticsMap-filtered.png}
\end{subfigure}
\begin{subfigure}{.325\textwidth}
    \centering
    \includegraphics[width=0.88\linewidth]{figures/SP2007/cubeCaustics.png}
\end{subfigure}
\caption{Images from left to right show: raw, unfiltered caustics map of a~water ripple; the same filtered caustics map; refractive caustics from animated water surface splatted onto the receiver's geometry. (\emph{Images from Shah et al.~\cite{SKP2007}}.)
}
\label{fig:SP2007-caustics}
\end{figure}

% TODO Hierarchical caustic maps
%\noindent\textbf{Wyman et al.~\cite{WN2009}}
% TODO adaptive caustic maps

% resolution: caustics blurring, multiresolution caustics map, 
% noise caustics triangles, splats

\noindent\textbf{Wyman et al.}~\cite{WN2009} propose adaptive caustics maps with deferred shading. They try to mitigate aliasing issues of caustics mapping, of the similar kind to shadow mapping, by adaptively sampling photons instead of the uniform sampling approach. Refractive objects are sorted and then rendered from back to front, after the background opaque geometry, into buffers. Then the photons are generated only on the refactors by sampling the buffers. 

A~coarse grid of photons is adaptively sampled in 2$\times$2 clusters of photons based on two metrics. The metrics work with the clusters, by either discarding them, or by dividing each photon in the cluster into its own cluster. The first metric tests whether all the photons in a~given cluster miss the refractor, if yes, then none are output, else the photons are divided until a~certain depth, in which case they are splatted into a~caustics map. The second metric works as the first, but continues dividing the cluster until it converges to a~single texel in the caustic map.

Furthermore, the photons are generated on the fly using the geometry shader, so there is no need for an explicit photon buffer. This allows for higher amounts of photons to be splatted, while maintaining interactive frame rates.

\vspace{1em}

\noindent\textbf{Meenrattanakorn et al.}~\cite{ML2022} present approach for rendering both refractive and reflective caustics caused by single specular interface, such as water surface. The approach is able to render caustics outside of the light frustum by considering \emph{cascaded} image spaces, i.e.,
the rays are traced from both the light space and the view space. They assume  deferred shading pipeline used with Reflective Shadow Maps (RSMs)~\cite{DS2005} of the light view and G-buffer of the camera view, and additionally a~separated depth buffer for opaque geometries behind the refractive object. 

RSMs are sampled randomly using blue noise to generate rays and scale the sampled flux by the ratio of the RSMs resolution to the total number of samples. Irrelevant points are removed in order to trace rays only from specular surfaces (using roughness cutoff) and with luminance, which is calculated using the flux, that is high enough to influence the caustics. Then, a~Fresnel factor is computed to probabilistically determine whether the ray is refracted or reflected.

The ray is traced by ray marching the depth buffers: first, from the light view, then the end position is projected to the camera view where it continues marching until the hit is found. The hit point is discarded if the dot product between the ray's direction and a~sampled normal is negative. The final hit point is found by interpolation between information from view-space depth and the scene geometry.

The resulting hit points are then saved into a~buffer as a~tuple of clip-space coordinates and the irradiance is computed. The irradiance is calculated using the world-space area of the pixel, similarly to Kim~\cite{Kim2019} (described in a~section below). The buffer is then used to draw the hit points as point primitives into a~2D photon map.

Finally, the photon map is denoised spatially and temporally. It uses the G-buffer's normal, depth, and history buffers along with a~motion vector. Radiance towards the camera is then evaluted by considering \emph{Lambertian reflectance} for diffuse surfaces.

Even though it can account for caustics outside the light frustum, the light path must be contiguously visible in the light and camera views, otherwise the caustics are lost. Furthermore, panning the camera might show noise in certain situations due to very sparse photon area. And finally, because the irradiance is computed in image space, other than diffuse surfaces are not physically plausible, since the incident vector is not stored in the photon map the radiance due to specular reflectance cannot be evaluated.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.7\linewidth]{figures/ML2022/output-cornell.png}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.7\linewidth]{figures/ML2022/output-sponza.png}
\end{subfigure}
\caption{Rendering of two scenes with a~specular interface producing refractive (emphasized by orange rectangle) and reflective (emphasized by green rectangle) caustics. The reflective shadow maps resolution is 2048$\times$2048 with up-to 2\,million photons. (\emph{Renderings by Meenrattanakorn et al.~\cite{ML2022}}.)
}
\label{fig:ML2022-output}
\end{figure}


\subsection{Ray Traced Approaches}

These techniques usually trace rays against the whole dynamic scene in world space using DXR. In this case, DXR bounding volume acceleration structures are assumed to be in place. Furthermore, deferred rendering system is employed, with G-buffer for rendering including the depth buffer. Due to rays being traced in world space, these techniques do not suffer from the same visibility limitations as image space techniques do.

Finally, to reach eliminate noise, denoising is applied. The objective of denoiser is to estimate the density of the photons in a~particular pixel by gathering nearby photons and accumulating their contribution at each pixel. This is comparable to finding a~full caustics light path ~\cite{ML2022}.

\vspace{1em}

\noindent\textbf{Kim}~\cite{Kim2019} presents an approach called \emph{screen-space photon mapping} (SSPM) to render reflected and refracted caustics. Instead of storing the photons as volumes in the world space, as it is with image-space techniques, the photons are evaluated and stored directly in screen-space, this is advantageous for deferred shading.

Photons are emitted only to areas bounded by so-called \emph{projection volume}\,---\,a~bounding box that contains objects generating caustics. The bounding box is then projected to the negative direction of directional light to obtain a~rectangular area that is used to compute emission flux for each photon proportionally to the area it takes in the photon grid.

Photons are traced until certain conditions are met or until a~surface hit is opaque enough. Any photon that hits a~location out of screen or that has negligibly small intensity, or is beyond screen depth is discarded. Finally, photons that were not discarded are compressed into a~single pixel no the screen. The final radiance is computed by normalizing the intensity by area of the pixel in world-space units. Due to the radiance being computed in screen-space, the information stored in G-buffer can be used to evaluate specular and diffuse colors based on BRDFs of the surfaces.

The rendered result is shown in Figure~\ref{fig:kim2019-output}. Reflection denoiser from NVIDIA GameWorks was used to gather photons based on camera data, depth buffer, roughness buffer and normal buffer. This technique produces artifacts on close-up due to the high resolution of photons near the screen frustum's near depth.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.7\linewidth]{figures/Kim2019/output-raw.png}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.7\linewidth]{figures/Kim2019/output-denoised.png}
\end{subfigure}
\caption{Images show refracted caustics by water surface in Cornell Box scene rendered with 2k\,$\times$\,2k photons. Left image shows the raw output, while the right shows the denoised variant. (\emph{Renderings by Kim~\cite{Kim2019}}.)
}
\label{fig:kim2019-output}
\end{figure}

\noindent\textbf{Gruen}~\cite{Gru2019s} renders one-bounce caustics in a~scene with a~refractive interface, e.g., water surface. The interface is assumed to be represented by a~triangular mesh in order to compute the triangular ratio similarly to Jensen et al.~\cite{JG2001}. The approach is summarized in five steps:
\begin{enumerate}
    \item A~refracted mesh of the interface mesh is computed by intersecting rays, which originate from each interface's vertex and point in the refracted direction of the incident light, with the scene's receivers. For each interface triangle there is a~resulting refracted triangle making up the refracted mesh.
    \item A~compression ratio is computed for each of the refracted triangle as the ratio of areas of an interface triangle and its corresponding refracted triangle. This, so-called \emph{compression ratio}, is thought of as a~metric of how much the light is focused or defocused.
    \item A~number of refracted light rays is adaptively generated at random locations within each triangle of the interface. The number of samples is based on the compression ratio: the higher the ratio, the more samples are generated; and on the on-screen size of the refracted triangle: the closer the triangle is to the viewer, the more samples are generated. %A~Ray does not get generated if it is shadowed from above, or if it is guaranteed that there will be no on-screen intersections.
    \item Each generated rays is then again intersected with the receivers in the scene. The intersection is then projected to a~screen space pixel where it is accumulated only if it corresponds to the front-most pixel at that position. A~brightness value gets accumulated into a~caustics buffer. The value can be scaled by the compression ratio and/or the amount of light that has been absorbed due to a~medium, similarly to Baboud et al.~\cite{BD2006}.
    \item Finally, the accumulated caustics buffer is denoised at half-resolution using a~set of iterated cross-bilateral blurring steps, that account for differences in view-space, depth, normal and view-space positions. Figure~\ref{fig:Gru2019-output} shows both the raw caustics buffer and the denoised buffer.
\end{enumerate}
Some refracted triangles might be off-screen while some of its generated rays might still produce on-screen caustics. This can be solved by also considering (in Step 3) three additional triangles at different distances along the rays that form the beam. Furthermore, the quality of the compression ratio depends on the geometric complexity of the surface mesh.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.8\linewidth]{figures/Gru2019s/output-raw.png}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.8\linewidth]{figures/Gru2019s/output-denoised.png}
\end{subfigure}
\caption{The left image shows the raw caustics buffer, and the right images shows its denoised variant. (\emph{Renderings by Gruen~\cite{Gru2019s}}.)
}
\label{fig:Gru2019-output}
\end{figure}

% TODO figure with the specular and refracted triangle that is sampled inside


\noindent\textbf{Yang et al.}~\cite{YO2021} present two specialized approaches for rendering caustics cast by water surfaces after one-bounce reflection or refraction. They build on the work of Gruen~\cite{Gru2019v}, \emph{Ray-Guided Water Caustics}, the continuation of his previous work~\cite{Gru2019s} on surface caustics, (explained in Section~\ref{sec:caustics-vol} about volumetric caustics). They intensify the optical details in water caustics by applying any of the two approaches: \emph{Photon Difference Scattering} (PDS) or \emph{Procedural Caustic Mesh} (PCM). They support most light types (including area lights) and also introduce \emph{cascaded caustics maps} (CCMs), i.e., caustics map at various resolutions to cover vast scene regions. 

The PDS technique treats ray hit points as photon footprints. The initial size of the footprint is determined by the resolution of the caustics map and the desired precision set by the user, then its scaling factor, an analogy to compression ratio in Gruen's method~\cite{Gru2019v}, is computed. The adjacent hit points are computed using photon differentials, i.e., using the finite difference of nearby rays (as illustrated in Figure~\ref{fig:YO2021-footprint}). Each photon footprint is then rendered as quad sprite and blended with the scene. The quad is placed perpendicular to the ray direction and then it is projected to the opaque surface. The intensity is cosine weighted by the incident angle and lighting is evaluated using the scene materials. Using elliptic shape instead of quad forms more continuous patterns even without any denoising (as shown in Figure~\ref{fig:YO2021-footprint}). Lower resolutions of caustics map produce blurry caustics, therefore this technique is coupled with the CCMs.

\begin{figure}[H]
\centering
\begin{subfigure}{.415\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/YO2021/ray-diffs.pdf}
\end{subfigure}
\begin{subfigure}{.285\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/YO2021/rectangular.png}
\end{subfigure}
\begin{subfigure}{.285\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/YO2021/elliptic.png}
\end{subfigure}
\caption{The Leftmost image illustrates the size of a~photon footprint as the combination of ray's origin, direction and hit point with the adjacent rays' origins and directions. The hit points of adjacent rays (red) are estimated using photon differentials from the current ray's hit (green). The other two images show the footprints' shape: rectangular (left) and elliptic (right). (\emph{Images from Yang et al.~\cite{YO2021}}.)
}
\label{fig:YO2021-footprint}
\end{figure}

The other approach, PCM, converts the hit points of the regular grid in the caustics map to an intermediate mesh with triangular topology. Any invalid primitive is discarded and contribution and intensity is computed according to its world-space area. Finally, the mesh is rendered into the caustics buffer and blended with the scene. It always produces sharp caustic patterns, however, it suffers from ``black edge'' artifacts at objects corners where the triangles are culled.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.8\linewidth]{figures/YO2021/PCM-2k.png}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.8\linewidth]{figures/YO2021/PDS-2k.png}
\end{subfigure}
\caption{Comparison between caustics maps, with resolution 2048\,$\times$\,2048, generated using procedural caustic mesh (left) and photon difference scattering (right) methods. (\emph{Renderings from Yang et al.~\cite{YO2021}}.)
}
\label{fig:YO2021-output}
\end{figure}
 

% ================================================================================
% ================================================================================
\section{Volumetric Caustics in Single-Scattering Media Methods}
\label{sec:caustics-vol}

Volumetric caustics are crepuscular rays or godrays underwater. It is the phenomenon when light rays are being visible as they propagate through a~water volume. It is caused by photons being out-scattered from their path due to the presence of particles in participating medium.

This section focuses on approaches that render volumetric caustics in single-scattering water volumes. The focus is mainly on refractive caustics, while most of the algorithms support also reflective caustics.

The methods investigated here are approximations of the photon mapping algorithm and most of them compute both surface and volumetric caustics, but only  volume caustics are of the interest here. The techniques described below are divided into two categories based on the level of abstraction, or more specifically, based on the primitive that is rendered: line-based approaches and beam tracing approaches.



% ------------------------------------------------------------------

\subsection{Line-Based Approaches}
These approaches work analogously to point splatting methods of surface caustics\,---\,by splatting line primitives instead of point primitives. Lines start the origin of the refracted ray and terminate at the intersection, while visible pixels on the rasterized line approximate in-scattering contribution. The visibility problem is solved by using a~depth buffer and additive blending. Analogously to splatting point primitives, these methods also suffer from noisy results needing a~greater amount of line primitives. This implies fill-rate intensive applications are instead oftentimes rendered at lower resolutions than the target screen, or must be blurred.

\vspace{1em}

\noindent\textbf{Papadopoulos et al.}~\cite{PP2009} utilize the image-space ray-scene intersection method by Shah et al.~\cite{SKP2007} by tracing a~uniformly distributed grid of photons form the light source. Scattering is simulated per photon by rasterization of lines between the water surface and scene intersection points. The light source tracks the camera to rasterize lines only in visible parts of the camera frustum. 

The final intensity along a~line is computed by applying Beer's law and Mie scattering phase function, accounting for attenuation due to the distance traveled underwater and towards the viewer. The intensity from all the lines are additively blended and low pass filtered (i.e., blurred). Figure~\ref{fig:ppout} shows the output of this method.

% TODO figure difure, photon grid 200 x 200, 100 x 100, accumulated into 128 x 128 render target, 
% also diagram
\begin{figure}[H]
    \centering
\begin{subfigure}{.6\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/caustics/pp-output.png}
\end{subfigure}
\caption{Rendering of method by Papadopoulos et al.~\cite{PP2009}, uses photon grid of 200$\times$200 accumulated into 128$\times$128 render target.}
\label{fig:ppout}
\end{figure}


\noindent\textbf{Hu et al.}~\cite{HDIGY2010} present a~line rendering method similar to Papadopoulos et al.~\cite{PP2009} but also account for multiple refractions and heterogeneous media.

To compute intensity along a~line, the line is sampled around an intersection point $\*{x}_l$ (between the viewing ray the light ray, shown in Figure~\ref{fig:hu-output}) at few regular steps. The final radiance is then the sum of radiances of each sample weighted by Gaussian distribution. That is, the closer the sample is to $\*{x}_l$, the greater is its influence on the final radiance. For heterogeneous media, the transmittance must be sampled and accumulated by marching along the line. The step size can vary based on the frequency of changes in the medium.


% Figure of the line with the point
% Lines is drawn from exit point $\*{p}_1$, on the underside of the specular object, to surface point $\*{p}_2$. $\*{x}_l$ is the intersection point between the view ray and the light ray.
% volume caustics buffer resolution: 1024 x 1024, 

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.7\linewidth]{figures/caustics/hu-drawing.pdf}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.6\linewidth]{figures/caustics/hu-output.png}
\end{subfigure}
\caption{Left image shows a~line drawn from exit point $\*{p}_1$, on the underside of the specular object, to surface point $\*{p}_2$. $\*{x}_l$ is the intersection point between the view ray and the light ray. Right image shows the output with volume caustics buffer of 1024$\times$1024.
(\emph{Images from Hu et al.~\cite{HDIGY2010}}.)
}
\label{fig:hu-output}
\end{figure}


\subsection{Beam Tracing methods}

This category of techniques traces the intensity changes inside the caustic volumes using triangular beams or more generally, prisms. Triangle beams are composed of the top part, the specular triangle of the water surface, and the bottom triangle, the caustic triangle, which is found by intersecting the refracted directions of incident light from each vertex of the specular triangle with the ocean bottom (shown in Figure~\ref{fig:warped}). The vertices of these two triangles are connected in a~fashion that a~triangular beam is formed. As the light focuses and defocuses due to the curvature of the water surface, the refracted directions at the vertices of the specular triangles point in different directions causing the beam to warp, so-called \emph{warped caustic volumes}, forming bilinear patches on sides (shown in Figure~\ref{fig:warped}). Such volumes cannot be easily rasterized, therefore bounding volumes are used to approximate the caustic volumes.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.6\linewidth]{figures/caustics/ernst-warped.pdf}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/caustics/caustic-beams.pdf}
\end{subfigure}
\caption{Left Image by Ernst et al.~\cite{EAMJ2005}: A~warped triangular beam is illustrated and cut through, the top triangle is called specular triangle and bottom triangle is the caustic triangle. It is formed by intersecting the bottom with each of the refracted directions at each vertex of the specular triangle. Right image by Liktor et al.~\cite{LD2011}: Shows the different types of caustic beams due to the divergence and convergence of light.
}
\label{fig:warped}
\end{figure}


\noindent\textbf{Jensen et al.}~\cite{JG2001} uses beam tracing to surface caustics. Triangular beams are formed with intensities according to the compression ratio: a~ratio of the area of the specular triangle to the caustic triangle. The beams are intersected with horizontal planes to form volumetric caustics by computing the ratio for each new caustic triangle that is the result of the beam intersecting a~plane. This results in a~blocky appearance of volumes, due to the ratio being fixed throughout the volume.

\vspace{1em}

\noindent\textbf{Ernst et al.}~\cite{EAMJ2005} presents a~technique that uses so-called \emph{warped volumes}, an optimization that allows for interpolation of intensities inside the volume, thus reducing the blocky appearance resulting in smoother caustics. Homogeneous media are assumed along with generators (the specular objects) represented only as triangular meshes. 

The receivers are not intersected directly, rather, the specular triangles of the generators are extruded and their axis-aligned bounding boxes (so-called \emph{bounding prisms} method, on CPU) are rasterized. Then for each pixel the bounding volume is intersected using a~scan plane along the ray from the camera (illustrated in left image of Figure~\ref{fig:warped}). Intersections of the refracted directions of each vertex of a~specular triangle forms a~caustic triangle. This triangle is then intersected with the camera ray to obtain a~point $\*{p}$ inside the triangle and the distance traveled through the beam $\Delta$. The radiance seen from the camera is then computed analytically and accumulated using additive blending, accounting for the attenuation due to the distance $d$ from the camera to $\*{p}$, and from $\*{p}$ to the specular triangle:
\begin{align}
    L = \Delta e^{-\extinctCoef (d)} \scatterCoef p(\theta) L_i(\*{r}) ,
\end{align}
where $\*{r}$ is the direction of the refracted light inside the beam, $L_i$~is the incident radiance of the caustic beam, and $\theta$ is the angle between the camera ray and $\*{r}$ directions. The $L_i$ is computed for point $\*{p}$ by, first, transforming the point to the space of them (called \emph{beam space}), and second, by trilinearly interpolating the radiance values at the vertices of the caustic triangles using barycentric coordinates. The beam space allows for the interpolation of the compression ratio inside the beam. %The output of this method is shown in Figure~\ref{fig:ernst}.



%Then each visible world space position is tested whether it is inside a~volume, using so-called \emph{point-in-volume} test. This is done by transforming the point to the local space of the beam (so-called \emph{beam-space}) obtaining thus a~caustic triangle, the point could possibly lie in. Point-triangle test is done to confirm this using barycentric coordinates. If the point lies in the triangles then its intensity is computed as interpolation between the intensities of the triangle's vertices reusing the barycentric coordinates. 

\noindent\textbf{Liktor et al.}, 2010~\cite{LD2010} extend the method proposed by Ernst et al.~\cite{EAMJ2005} in the following ways:
\begin{itemize}
    \item A~separate triangular mesh is used for the beam generation, so-called \emph{projected grid}, instead of the mesh geometry of generators. The grid is projected onto the specular objects by rendering them from the light source. This way, the algorithm is independent of the geometric complexity of the generator meshes.
    \item Provides full GPU implementation by extruding the beams in geometry shader. 
\end{itemize}
The scattered radiance towards the camera is computed the same way (using the scan plane), i.e., single-scattering in homogeneous participating media, that is additively blended.

\noindent\textbf{Liktor et al.}, 2011~\cite{LD2011} extend the previously mentioned technique by the authors by introducing:
\begin{itemize}
    \item hierarchical projected grid: a~projected grid that is adaptively refined in geometry shader according to heuristics, such as on the discontinuities of the specular surface, based on the distance to the viewer, and whether the beams are convergent or divergent.
    \item bounding volume extrusion: even more tightly fitting bounding box around the bounding volumes is found by modification of the original bounding prism method, drastically improving for the worst case when the beam is warped, leading to less fill-rate.
\end{itemize}
The radiance scattered towards the viewer is computed the same way as in the previous technique\,---\,by intersecting a~beam with a~scan plane from the camera. The rendering of this method is shown in Figure~\ref{fig:LD2011-output}.

\begin{figure}[H]
    \centering
\begin{subfigure}{.5\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/caustics/ld2011-output.png}
\end{subfigure}
\caption{Rendering of method by Liktor et al.~\cite{LD2011}.}
\label{fig:LD2011-output}
\end{figure}
% TODO figure


\noindent\textbf{Gruen}~\cite{Gru2019v} extends his method for surface caustics~\cite{Gru2019s} (described in Section~\ref{sec:caustics-surf}) to render both surface and volumetric caustics in single scattering participating media using ray tracing and rasterization pipelines. The bounding volumes problem is approximated by adaptively tessellating the triangular beams into smaller beam volumes. 

This in-scattering integral is approximated by tessellating the caustics volumes into small enough volumes that allow for the accumulation of the in-scattering term only at the frontfacing triangles of each volume. This way, there is no need for the camera to intersect the beams directly. Rather, the rastarization pipeline is utilized, along with the ray tracing pipeline to find intersections with a~dynamic scene. In comparison to the surface caustics method by the author, this method is more concerned with refining the beams rather than the surface caustics. This method reuses ideas from caustics mapping, in the following way:
\begin{enumerate}
    \item The compression ratios are computed for the refracted water mesh as the ratio of the areas of the water triangle and the refracted triangle. The water surface mesh must be detailed enough to facilitate high-enough-quality compression ratio.
    \item The visible portion of the water surface is rendered from the light's point of view into \emph{caustics map}\,---\,a~texture containing positions and normals of the surface.
    \item Each pixel of the caustics map forms a~refracted ray that is intersected with the scene and the intersection is stored in the \emph{refracted caustics map} at the same pixel position.
    \item Surface caustics are rendered by projecting the intersection points from the previous step to screen space. Then testing, whether they are the frontmost pixel on the screen. The radiance value is computed, accounting for absorption by the distance traveled, similarly to Baboud et al.~\cite{BD2006}, scaled by the compression ratio, and accumulated into photon map using atomic add operations (analogous to additive blending, because many intersection points can end up on the same screen location).
    \item The water triangles are adaptively tessellated by drawing the water mesh. The compression ratio is thought of as the likelihood that a~triangular beam forms a~non-convex volumes, therefore, the higher the ratio the higher the tessellation factor. Furthermore, it is scaled to account for any light leaks doe to small scene features, or to approximate the in-scattering integral enough. This is rather artistically controlled.
    \item The geometry shader is used to build the triangular beam volumes by projecting the tessellated points to the caustics map and the refracted caustics map to obtain positions of the vertices of the top and bottom triangle, respectively. The vertices of the triangles are then connected to form a~triangular beam.
    \item Finally, the volumes are rasterized and the in-scattering is computed on the pixels of the frontfacing sides of each volume, evaluating phase function for the interpolated ray direction inside the volume, multiplied by the interpolated thickness of the volume. The radiance is additively blended into an offscreen buffer that can be directly used in deferred shading pipelines.
\end{enumerate}
The resulting volumetric caustics buffer can be further denoised similarly, to the surface caustics method. 

%The author provides additional guidelies for the tessellation of the beams. GDC19 presentation of this technique~\cite{Gru2019prez}, the author explains and mentions additional guidelines for the tesellation of the beams:


\ifdefined\VOLUMECAUSTICSNOISE
there the beam ratio (as computed by Jensen et al.~\cite{JG2001}) is interpolated 

volume caustics intensity is computed and accumulated using additive blending.

The so-called point-in-volume test is done 

- assumes that generators are triangular meshes, 
- homogeneous media
- linear interpolation of point inside caustics triangle using barycentric coordinates, 
- intensities are computed at each vertex of the triangle
- each point on a~receiver is tested wether it is in caustic volume
    - point-in-volume test
    - uses beam space to 
- analyiticall evaluation of the scattering contribution:
    attenuation to the point in the volume, to the viewer,
    times thickness of the volume - the distance the eye traveled inside the caustic beam,
    phase function
    
Volumetric Caustics
- considers only homogeneous media,
    allows for analytical evaluation of the scattering contribution due to a~caustic beam:
\begin{align}
    L_e (\vec{\omega}) = \Delta x  e^{-\sigma_t (d_e + d_s)}
                         \sigma_s p(\vec{\omega} \cdot \vec{\omega}') L_i(\vec{\omega}')
    ,
\end{align}
where $L_e$ is the radiance seen from the camera, $\sigma_t$~is the extinction coefficient, $\sigma_s$~is the scattering coefficient, $d_e$~and $d_s$~are the travel distances of the ray from the camera and the caustic beam, $p$~is the phase function of the medium, $L_i$~is the incident radiance of the caustic beam, and $\Delta x$ is the distance traveled of the camera ray inside the caustic beam.

% TODO figure 2 in Ernst



% TODO more from previous work of Ernst

\subsection{Ernst, Akenine-Möller, Jensen. Interactive Rendering of Caustics using Interpolated Warped Volumes}

\todo{what does it do}
- interactive rendering of caustics from refractive and reflective surfaces as well as volumetric caustics
\todo{what does it bring to the table, what does it solves better?}
- using programmable graphics hardware
- previous real-time methods used simple prisms for the caustics volumes and a constant intensity approximation at the receiver
    - here the approach uses interpolated caustics volumes in result renders smooth caustics
    % how does it do it?
    - uses a~formula for evaluating the density of wave-fronts along a caustic ray
        - derived a~formula for evaluating the 
    - developed a~method for rendering caustic volumes bounded by bilinear patches
- both optimizations are suitable for programmable graphics hardware

- mainly builds on techniques: Backward beam tracing, Nishita. Method of Displaying Optical effects within water using accumulation buffer.
- the rendering of volumetric caustics is similar to the apporach for real-time rendering of shadows using the shadow volumes algorithm~\cite{TODO_shadow_volumesis5}.
% Franklin C. Crow. Shadow Algorithms for Computer Graphics. In Computer Graphics (SIGGRAPH 77), pages 242–248. ACM Press, 1977.

\todo{How does it do it, specifically, but not in too much detail, every one can read it}

- removes the blocky appearance of caustics by
    interpolating the light intensity over each caustic triangle
- uses more accurate method, the \emph{warped caustics volumes} for finding bounding volumes around the caustic volumes than the prisms used in previous works.
- uses analytic computations for volume caustics in homogeneous participating media
- the algorithm developed is suited for the standard programmable rasterization pipeline.

Algoritm:

divides geometry in the scene into (caustics) \emph{generators} and \emph{receivers}
where generators, are objects with a~specular component in their BRDF
while receivers always have a~diffuse component

Below is the pseudo code of rendering caustics using a~caustic volume approach:
\begin{verbatim}
    foreach visible point P on a receiver
        foreach caustic volume V
            if P is inside V
                compute & accumulate caustic intensity
\end{verbatim}
- this algorithm is later implemented in a~form that is suitable for GPUs

The description of the algorithm then for simplicity assumes a~scene with one refractive generator (e.g., a~water surface) and one light source.
- it can handle fully dynamic scenes

- caustic volumes are computed on the CPU-side:
    - for each vertex $\mathbf{v}^i$ of the generator, a~refracted vector $\mathbf{r}^i$ is computed using the direction to the light source and the vertex normal $\mathbf{n}^i$
        - then for each triangle $\triangle^{jkl}$ consisting of vertices $\mathbf{v}^j$, $\mathbf{v}^k$, and $\mathbf{v}^l$ a~caustic volume is generated
        - the volume is bounded by the rays $\mathbf{v}^j + t\mathbf{r}^j$, $\mathbf{v}^k + t\mathbf{r}^k$, for $t \geq 0$.
Then the caustics are rendered in two passes:
First pass:
    renders the world space positions of the receivers into a~texture
    - receiver geometry can be rendered in arbitrary form, not just as triangular meshes
Second pass:
    - a~bounding volume is drawn for each caustic volume
    - using the positions from the texture, each final pixel undergoes a~point-in-volume test
    - if point is inside the volume, then the caustic intensity is computed and accumulated in the frame buffer using additive blending.

Limitations:
- accounts only for a~single specular bounce
- no shadowing for caustics
- generators must be triangular meshes
- receivers on the other hand, must be possible to render to into the depth buffer.
- reflected/refracted vector are linearly interpolated over specular triangles
-> depends complexity of the generators to some degree

Warped Caustic Volumes
- improves on the assumption that the volumes are prisms, 
    that is, their sides are planar surfaces
    which is generally not true, due to the different normals at the vertices of a~specular triangle,
    see Figure TODO

% Figure 3. with description

- instead makes certain assumptions and modifications to optimize the point-in-volume test

Instead of planar side surface makes a~more precise approximation for refraction:
    New reflected/refracted directions at vertex $\mathbf{v}^i$ to another $\mathbf{v}^j$ can by created by interpolation of $\mathbf{r}^i$ to $\mathbf{r}^j$

The coordinate system

Point-in-Volume Test
- is executed for every pixel covered by the bounding volume of each caustic volume
% refers to figure 3

- tests if a~point $\mathbf{p}$ in world space is inside a~volume~$V$
\begin{enumerate}
    \item Point $\mathbf{p}$ is transformed into the local coordinate system of $V$, resulting in point $\mathbf{p}'$. All subsequent computations are carried out in this coordinate system.
    \item The vertices $\mathbf{c}^0$, $\mathbf{c}^1$, and $\mathbf{c}^2$ of the caustic triangle $\triangle_c$ are computed using a~virtual plane. The plane with normal $\mathbf{n} = (0,1,0)$ containing $\mathbf{p}'$ intersects the caustic rays, yielding the caustic vertices (as depicted in Figure~\ref{TODO} to the right).
    \item Point-in-triangle test is employed to check whether $\mathbf{p}'$ is inside $\triangle_c$.  The test applied is based on computing three values that correspond to barycentric coordinates, because these values are later needed for linear interpolation. The test is further simplified to two dimensions by projecting the involved points along the local y-axis onto a~plane. If all three computed values have the same sign, then the point is inside the triangle.
\end{enumerate}

Caustics Interpolation
- the same idea: based on the ratio of areas of the specular and caustic triangles
- the intensities inside the caustic triangles are interpolated by computing intensity values per vertex rather than per an entire caustic triangle
- Derives a~function to calculate the intensity at a~given depth inside a~caustic volume for the point $\mathbf{p}'$.

- Defines the double area of a~caustic triangle as:
\begin{align}
    A(\triangle_c) = | (\mathbf{c}^1 - \mathbf{c}^0) \times (\mathbf{c}^2 - \mathbf{c}^0)|
\end{align}
By utilizing the same coordinate system from Figure~\ref{TODO-coord-sys3}, the area function is rewritten as a~function of depth along the $y$-axis using the so-called \emph{triangle area coefficients}. These coefficients are precomputed for each vertex in the following way.

With the assumption that all caustic volumes are described in the same coordinate system, the interpolation is done as follows:
\begin{enumerate}
    \item The area coefficients for a~vertex $\mathbf{v}^i$ of a~generator are computed by summing up of the triangle area coefficients of each triangle around the vertex (basically, a~triangle fan around the vertex is drawn). 
    \item The first area coefficient $\overline{a}_i$ is the sum of the areas of the specular triangles around the vertex. The other two are used to compute the area of the caustic triangle at certain depth $\overline{A}_i (\alpha)$. Then the intensity of the vertex is the ratio of $\overline{a}_i$ and $\overline{A}_i (\alpha)$.
    \item An interpolated intensity at the point $\mathbf{p}'$ in a~caustic volume, defined by the vertices $\mathbf{v}^0$, $\mathbf{v}^1$, and $\mathbf{v}^2$, is then computed using the area coefficients of each vertex and the barycentric coordinates from the point-in-triangle test.
\end{enumerate}
For the interpolation to work, all caustic volumes must be described in the same coordinate system. This is done by flattening the fan around vertex~$\mathbf{v}^i$ for the computation of the coefficients as depicted in Figure~\ref{fig:EAMJ05:flattening}. The local $y$-axis is chosen as the vertex normal~$\mathbf{n}^i$. All vertices in the fan are projected along their reflected/refracted vectors into a~plane with normal $\mathbf{n}^i$ containing $\mathbf{v}_i$.
Also, the $\mathbf{p}$ does not have to be transformed into three different coordinate systems to compute the three vertex intensities, instead the $y$-coordinate of $\mathbf{p}'$ can be used. This is due to the fact, that the area coefficients are computed for every vertex of every triangle of the flattened fan. This way, the coordinate system for the fan is the same as the coordinate system of caustic volume.


% TODO fig:EAMJ05:flattening

Volumetric Caustics
- considers only homogeneous media,
    allows for analytical evaluation of the scattering contribution due to a~caustic beam:
\begin{align}
    L_e (\vec{\omega}) = \Delta x  e^{-\sigma_t (d_e + d_s)}
                         \sigma_s p(\vec{\omega} \cdot \vec{\omega}') L_i(\vec{\omega}')
    ,
\end{align}
where $L_e$ is the radiance seen from the camera, $\sigma_t$~is the extinction coefficient, $\sigma_s$~is the scattering coefficient, $d_e$~and $d_s$~are the travel distances of the ray from the camera and the caustic beam, $p$~is the phase function of the medium, $L_i$~is the incident radiance of the caustic beam, and $\Delta x$ is the distance traveled of the camera ray inside the caustic beam.

To evaluate $L_e$ the values of $d_e$~, $d_s$~and $\Delta x$ are computed in the following manner: 
\begin{enumerate}
    \item A~bounding volume is drawn for each caustic volume. The volume is found using the technique below, called \emph{Bounding prisms}.
    \item For each pixel of the rasterized bounding volume, a~scan plane is intersected with the three caustic rays. The resulting triangle, while having curved edges, is assumed to be linear.
    \item A~ray from the camera in the scan plane is intersected with the triangle edges yielding the entry point $\mathbf{p}_n$~and the exit point $\mathbf{p}_f$ of the volume.
    \item $\Delta x$ is the distance between $\mathbf{p}_n$ and $\mathbf{p}_f$, $d_e$ and $d_s$ are the distances from $(\mathbf{p}_n + \mathbf{p}_f) / 2$ to the camera and to the specular triangle, respectively.
\end{enumerate}

Bounding Prisms
- instead of drawing a~screen-size quad for each caustic volume to point-in-volume test for each pixel
- an axis-aligned bounding box is found and rasterized, the goal is to find a~tighter bounding volume to have less point-in-volume tests
- a~caustic beam is a~triangular prism:
    has specular triangle at the top with vertices $\mathbf{v}^i$,
    it has finite length, there is a~bottom triangle composed of vertices $\mathbf{b}^i$
- a~tight plane is found for each edge of the specular triangle containing the edge.
    - describes how to compute such a~plane for each edge that contains the remaining vertex of the specular triangle in the negative half space of the plane
    - the resulting planes contain the entire caustic volume in their negative half spaces
    - that is, the resulting prism must also contain the vertices of the bottom triangle
- once all the planes were found for each edge of the specular triangle, then a~ray is computed for each vertex as the intersection of the adjacent edges of the planes
- finally, a~bottom triangle is found by making sure that all $\mathbf{b}^i$ are inside the the volume.

Figure~\ref{fig:EAMJ2005:} shows the volumetric caustics rendering in a~homogeneous participating media
% TODO Figure 9
Dispersion is also visible by having slightly different refraction indices for R, G, and B.



\subsection{Liktor. Dachsbacher. Real-time Volumetric caustics with Projected Light Beams}

Mainly builds (TODO better word??) on the method of Ernst et al.~\cite{EAMJ2005} outlined in the previous Section.
- method for visualization of volumetric caustics in single-scattering participating media
- can visualize both reflective and refractive one-bounce caustics
Instead of using the geometry of generators a~separate \emph{projected grid} in light's image space is used. Thus making the solution independent from the geometric complexity of generators.
Also, uses the geometry shader to extrude and accumulate the caustic volumes into the final volumetric effect.
- proposes fully GPU-based implementation of the volumetric caustics algorithm
- the projected grid is used for the generation of beams
    - a~vertex grid  projected to the surface in light's image space
    and the beams are extruded  in the geometry shader

Algorithm
- Works similarly to Ernst et al.~\cite{EAMJ2005}, but now mainly the GPU is utilized

Beam generation
\begin{enumerate}
    \item The caustics generator surfaces from the light source are rendered into geometry buffers. The surface world space positions, normals and material properties are stored. Surfaces of arbitrary number and complexity can be used as generators.
    \item A~regular grid in light's image space is projected to the specular surfaces. Each vertex of the grid samples values from the geometry buffers. Any primitives of the grid that are not covered by the generator surfaces are discarded in the geometry shader. This is done by invalidating the values in the buffers prior to rendering.
    \item The triangles that are bases of caustics beams are extruded in the geometry processing stage. Thus new primitives that form the convex boundaries of the caustic volumes are emitted. The lengths along the refracted vectors of each vertex of the specular triangle are set to a~constant value, but can be determined by a~depth-texture lookup or another approximate ray-tracing methods. The triangular mesh bounding the bilinear patches is found using the \emph{Bounding prisms} method of Ernst et al.~\cite{EAJM2005} described in the previous Section. 
    \item The fragment shader then uses the information passed from the geometry processing stage to solve the volume rendering equation by sampling the beams.
    
\end{enumerate}

Fragment shader
- finds the radiance at each sample point in the beam
- uses the method described by Ernst et al.~\cite{EAMJ2005}:
-   the sampled radiance depends on the initial caustic radiance on the surface and on the area ratio of the specular triangle and the parallel cross section of the beam at the sample point

- assumes that the cross section is parallel to the specular triangle and computes the radiance at each point using the area ratio

% TODO figure fig:LD2010:beam-space

% TODO nomenclature vector vs. point: \bf vs \vec
To compute the area of the cross section of the warped volume 
- uses the same special coordinate system as proposed by Ernst et al.~\cite{EAMJ2005} now called \emph{Beam space} in a~manner illustrated in Figure~\ref{fig:LD2010:beam-space}.
- where $\mathbf{v}_i$ are the vertices of the specular triangle, $\mathbf{v}_0$ is the origin, $y$-axis is the triangle normal, $\mathbf{v}_0 - \mathbf{v}_1$ edge is the $x$-axis and the cross product of the latter two normalized vectors is the $z$-axis, thus making an orthogonal basis in beam space. The ray direction vectors $\mathbf{d}_i$ are also scaled so that their $y$-coordinate equals to one in beam space.
- the area function is then derived as in Ernst et al.~\cite{EAMJ2005} and computed using the area coefficients in the geometry shader. 
- the fragment shader then evaluates the volumetric rendering equation by ray-marching in beam space

Geometry shader is also used to remove triangles that would result in incorrect caustics
- due to the projected gird being independent of the scene, it is possible that some triangles are extremely distorted, for example, because they do not lie on the same surface.
- to discard such triangles the following heuristics is used:
    - dot product of the triangle normal $\mathbf{n}_{\triangle}$ and the light direction vector $\mathbf{d}_l$ is taken
    - the caustics effect is said to be most accurate if $\mathbf{n}_{\triangle}$ is almost parallel to $\mathbf{d}_l$
    - if they are almost perpendicular to each other, the the triangle is discarded
the result of the application of this heuristics is shown in Figure~\ref{fig:LD2010:heuristics}

% TODO figure 4: heuristics

Ray-Marching Caustics Volumes
- the volume rendering equation is solved in a~similar manner to Ernst et al.~\cite{EAMJ2005}
- the camera ray is intersected with the beams, 
- then the VRE is solved for a~set of sample points inside them
- assumes VRE for homogenous participating media in the following form:
\begin{align}
    L (\vec{\omega}) = \int_{\Delta x} e^{-\sigma_t (s_1 + s_2)} \sigma_s P(\vec{\omega} \cdot \vec{\omega}') L_{in} (\vec{\omega}') ds ,
\end{align}
% TODO need to describe omega_t and s ?? use YOUR NOTATIONS
where $L{\vec{\omega}}$ is the radiance seen from the camera, $\Delta x$ is the distance inside the beam along the camera ray, $\sigma_t$ and $\sigma_s$ are the extinction and the scattering coefficients, respectively, $P$~is the phase function, and $s_1 +s_2$ is the length traveled by the caustics photons from the specular triangle into the camera. The terms are illustrated in Figure~\ref{fig:LD2010:notations}.

% TODO  Figure~\ref{fig:LD2010:notations}

This equation is then approximated and used in the same form as in Ernst et. al.~\cite{EAMJ2005} TODO ref to eq., here for convience:?? TODO?
The fragment shader then perfoms the sampling analogously to Ernst et. al.~\cite{EAMJ2005}:
\begin{enumerate}
    \item The view ray is transformed into beam space. Using the vertex attributes from the geometry processing stage.
    \item A~scan plane containing the view ray is intersected with the beam's caustic photon directions to form a~cross section triangle. This is illustrated in Figure~\ref{fig:LD2010:notations}. The normal of such a~plane is computed as:
    \begin{align}
        \mathbf{n}_p = \mathbf{r}_{view} \times (\mathbf{r}_{view} \times \mathbf{n}_{\triangle}) ,
    \end{align}
    where $\mathbf{r}_{view}$ is the direction of the view ray, $\mathbf{n}_{\triangle}$ is the triangle normal, which is the $y$-basis in beam space.
    \item The cross-section triangle is intersected with the ray and the radiance at sample point $\mathbf{p}$ is computed by trilinearly interpolating the radiance values at the corners. Since the triangle and ray lie in the same plane, the intersection is done by finding the intersections with the edges (2D line-line intersection). If the view ray does not intersect the triangle, the fragment is discarded.
\end{enumerate}

% TODO figure result??

This approach suffers from:
- must find the right compromise between shading quality and fill-rate, that is determined by the grid resolution, 
- the larger the grid, then the more accurate the caustic effect, but more beams get rasterized to the same pixel
- smaller grid may also introduce sampling artifacts, mainly at geometric discontinuities



\subsection{Liktor, Dachsbacher. Real-Time Volume Caustics with Adaptive Beam Tracing}

method for rendering surface and volume caustics using approximate beam tracing
- supports two-sided refractions, inhomogeneous participating media
- multi-bounce light-surface interactions at real-time frame rates

% TODO figure 1 with te ocean and then the skull, because it is refed later

- single-scattering participating media
- based on caustic mapping, triangle-based volumetric methods

Improves on the previously described method by the same authors in the following aspects:
- the previous variant suffers from high fill-rate due to rasterizing a~large number of beams,
    here the method presented is based on \emph{adaptive beam tracing}
    - renders the scene from the light source's view to capture the first light-surface interaction,
    - then again uses this projected grid to adaptively generate beams
        - it adapts to geometric discontinuities of the specular objects, to the light's intensity, in the case of spot or textured light sources, and to the beam's contribution in the view space 
        - the grid is iteratively refined according to discontinuities in the geometry and photon distribution
    - to handle two-sided refraction and multiple bounces, approximate ray tracing algorithms are used
    - handles full dynamic scenes
    - tighter bounding volumes are found to reduce fill-rate when accumulating the inscattered light from the caustic beams

- combines ideas of caustic mapping approachs with triangle-based caustic rendering
    - adaptive sampling int he light's image space

Overview of the algorithm is shown in Figure~\ref{fig:LD2011:algorithm}. First, all directly visible surfaces from the light source are rasterized into geometry buffers. The geometry buffer is then resampled using a~triangular 2D grid. The grid starts with a~coarse resolution independent of all the surface details and is adaptively refined until certain conditions are met. Then each triangle in the sampling grid forms a~base of a~caustic beam, which is traced across surfaces to accumulate focused light in the volume. For every beam its bounding volume is found and is rasterized to determine the radiance scattered towards the camera depending on the phase function and density of the participating media. 

% TODO Figure~\ref{fig:LD2011:algorithm}

In addition, occlusion and higher order caustics are considered. For this case the primary beams after generation are stored for later passes. To find intersections with the surrounding scene geometry approximate ray tracing is used. 
% TODO should cite Szirmay-Kalos et al. 2005?
The higher order beam tracing suffers from ??
    the beams must be split where a~surface is partially intersected. This is handled similarly to the adaptive grid refinement while using scene geometry impostors rather than the light space geometry buffer.

TODO The algorithm is ... outline, summarized, related .. .
- how each step of the algorithm is addressed in comparison to the previous method
- mainly focus on the improvements
\begin{enumerate}
    \item Beam generation
    - rasterize the specular surfaces visible from the light source to off-screen geometry buffers
    - store: positions, normals, material properties (including transparency, index of refraction, etc.)
    - a~planar grid is projected to the lit surfaces, by
        -  each vertex sampling the corresponding G-buffer location.
    - initially, the grid is regular and coarse
    - each cell of the regular grid is treated as two triangles that later form triangular bases of caustic beams
    - discards invalid triangles:
        - crossing the silhouettes of specular objects, or
        - lying outside

    TODO WHY?
    - also takes into account light's emission characteristics, and the beam's contribution to the final rendered image

    Adaptive refinement
    - first undersampled regions in the grid are detected
    - then the corresponding grid triangles are split into smaller ones
        depends whether a~beam should be kept, subdivided or discarded.
        this is decided based on an empirically set-up heuristic composed of three features??:
        beams are refined if:
        \begin{itemize}
            \item span across geometric discontinuities. 
            \item Beams' energy and their contribution to the final image
                - divergent beams
                - convergent beams
        \end{itemize}
    - the refinement process is implemented using a~geometry shader, with stream back,
        is usually run in 3 to 4 iterations.
        
    \item Bounding volume extrusion
        bounding volume is generated for every beam

        to determine the extrusion length of the beams,
            - the maximum length is assumed to be computed by intersecting the beams with the view frustum
            - occlusion from other objects in the scene is also handled using approximate ray tracing technique described by Szirmay-Kalos et al.~\cite{SKALP2005}.
            this is done by geometry shader
            - also only beams that have a~significant contribution to the final image are considered,
                - that is, if a~beam is divergent, then from a~certain distance the contribution if considered negligible
                - the maximum distance is computed using energy approximation

        to compute the area of the cross-section of a~beam at any distance along the rays, the special coordinate system, here called \emph{beam space}, by Ernst et al.~\cite{EAMJ2005} is applied
        
        bounding volumes
        - generally, the sides of the caustic beams are bilinear patches that cannot be represented by a~finite set of triangles
        - proposes method that improves on the one of Ernst~\cite{EAMJ2005},
            - the original method constructs three touching planes for each edge of the base triangle and one of the top vertices, in a~way that the entire beam liea in the negative half space of the plane (shown in Figure~\ref{fig:LD2011:bounding-volume}. The intersection of the planes yields the edges of the bounding volume.
        - the original method works well for slightly warped beams, but in the worst case, when the beam is strongly warped due to the beam focusing mid way so that the top triangle ``flips over'' (as depicted in Figure~\ref{fig:LD2011:bounding-volume}), it generates overly large volumes
        - introduces \emph{improved bounding volumes}. which results in tighter volumes 
            - works analogously, but uses the top edges as base edges of the bounding prism and then intersects the two prisms with each other
            - simplifies this by:
                - each bounding plane is defined by one triangle edge and one vertex of the opposite triangle
                - so the formed triangle is taken directly as one side of the bounding volume
                - then the final bounding volume is defined by six triangles, selected for each edge independently
                - selects six triangles using the edges of the two triangles and one vertex for each of the opposite triangle
            - whenever there are multiple triangle candidates for any edge, the original bounding prism method is used instead of performing another test.
        
    % TODO Figure~\ref{fig:LD2011:bounding-volume}
    \item Inscattering and surface caustics
    - the fragment shader then computes the amount of radiance scattered towards the viewer
    - this is done by solving the same approximation of the integral in Equation~\ref{eq:LD2010:approx}
    - a~single sample is taken as a~sufficient approximation of the integral (in Equation~\ref{eq:LD2010:integral}, assuming that the beams are narrow enough in image space
    - for the evaluation, the following steps are taken for every fragment:
    \begin{itemize}
        \item Camera ray corresponding to each fragment intersects the beam. If the intersection is found, then the location is computed. This is done using the scan plane-based method by Ernst et al.~\cite{EAMJ2005}. The intersections of the scan plane and the edges of the beam form a~caustic triangle containing the camera ray (shown in Figure~\ref{fig:LD2011:intersection}, with the notations described in the previous Section~\ref{sec:LD2010:intersection}).
        \item To minimize the error caused by the curved edges of the caustic triangle, the normal vector of the plane is chosen to be close to the normal of the caustic triangle, as in the Equation~\ref{eq:LD2010:scan-plane-normal} from the previous method by the same authors.
        \item The length of the camera ray segment inside the caustic triangle is obtained by a~2D ray-triangle edge intersection.
    \end{itemize}

% Figure~\ref{fig:LD2011:intersection}
% todo describe the notation

    - the radiance scattered towards the viewer is computed:
        - 
        - the inscattered radiance is computed as described by Ernst et al.~\cite{EAMJ2005}
        - 
    \item Surface caustics
        - to render surface caustics the very same beams as for volume caustics are used.
        - additionally, caustics receivers are rendered into a~geometry buffer from the camera point of view:
            diffuse surface colors and depth values are stored
        - the previously mentioned fragment shader is extended in the following way:
            during rasterization, the camera-view geometry buffer is sampled and point-in-volume tests are performed on each surface point covered by the beam, similarly to Ernst et al.~\cite{EAMJ2005}
        - 
    \item Multi-bounce caustics
\end{enumerate}

Limitations
- 



\subsection{Gruen. Ray-Guided Volumetric Water Caustics in Single Scattering Media with DXR}

Water surface triangles are adaptively tessellated:
    - the beam compression ratio is thought of as a~value describing the likelihood that a~triangular beam forms a~non-convex volume
    - therefore it used as the initial tessellation factor to subdivide coarse triangular beams into smaller beams
        - sort of a~``guide'' for the tessellation of the beam
    - then this tessellation is factor is scaled appropriately to conform to the following guidelines:
    - sort of guildelines
    \begin{enumerate}
        \item Provides enough slices for well enough approximation of the in-scattering integral.
        \item Prevents the triangular beams from turning convex.
        \item No volumetric light leaks through small scene features.
    \end{enumerate}

    GDC19 presentation of this technique~\cite{Gru2019prez}, the author explains and mentions additional guidelines for the tesellation of the beams:
    \begin{enumerate}
        \item The triangle/beam is culled if the beam compression ratio is too low. This is set by the artist.
        \item The triangle is culled if the beam does not intersect the camera frustum.
            % TODO in clip space - simple , 
        \item Tessellation factor is computed to the size of the caustic map.
        \item Tessellation factor is adapted to the size of small features in the scene.
    \end{enumerate}
    The last two guidelines are related: the caustic map must be high resolution enough to account for small features in the scene, if it is not, then the higher tessellation factor is required, this is usually also set by the artist.

Triangular beam volumes are built

    Geometry shader is run for every tessellated water triangles to build the triangulated hull of the corresponding triangular beam

    - the 3D vertices of the incoming triangle are projected to the (refracted) caustics map space
    - vertices of the specular triangle are sampled from the caustic map
    - vertices of the base triangle are sampled from the refracted caustics map

    - the triangles are generated in such a~way that their normals point out of the bounding volume

    thickness
        - the area of triangle

\fi% ================================================================================





% ================================================================================

% ---------------------------------------------------------------------------
% ===========================================================================

% ================================================================================
% ================================================================================
% ================================================================================
\chapter{Design of Technique for Rendering Underwater Scenes}
\label{ch:design}
% something like light transport underwater
% solving the light transport inside the water

The goal of this thesis is to visualize underwater scenery in high quality using a~selection of graphics techniques. Theory behind rendering participating media was given in Section~\ref{sec:light-media}. The underlying field of ocean optics was studied and is described in Chapter~\ref{ch:optics}. 

Then, various techniques for simulation and rendering ocean environments were investigated in Chapter~\ref{ch:ocean}. Specifically: the simulation of water surfaces, existing bio-optical models that abstract the optical properties of waters, participating media approaches that render such media, and procedural generation for generation and representation of various surfaces that bound the water body. Chapter~\ref{cha:caustics} then focused solely on the techniques that simulate lighting effects inside water bodies, i.e., caustics.

First the scene representation is proposed. Then the technique to render such scene is designed. A~subset of the designed techniques was implemented and the results are shown in Chapter~\ref{ch:impl}.

% ================================================================================
\section{Scene Representation}
\label{sec:design-scene-repr}

As described in the previous chapters, light transport in underwater environments is influenced by many complex processes and components. The radiative transfer equation connects the optical properties of the water body with the light within the water. Before the light can be simulated in the water body, the boundary conditions and the optical properties of the water body must be defined. The scene is composed of the following components that try to reproduce an oceanic environment:
\begin{enumerate}
    \item Sky model: the Preetham's precomputed sky model~\cite{PSB1999} was used.
    \item Water surface model: based on the Tessendorf's choppy model~\cite{Tes2001} for simulation of waves of deep ocean waters.
    \item Terrain: the ocean floor, simply terrain, is modeled as a~signed distance field/function in combination with fractal Brownian Motion approach.
    \item Bio-Optical Model: a~physically-based model that tries to reproduce optical properties of case\,1 waters. A~vertical distribution of phytoplankton is defined and based on that, the other optical parameters are deduced, this results in vertically non-uniform parameters.
\end{enumerate}
The following sections describe these models in more detail. The description and survey of sky models are out of scope of this thesis. The choice of the sky model was driven purely by performance and ease of use.

\subsection{Water Surface Model}
\label{sec:design-surface}

Water surface is modeled using the Choppy waves method by Tessendorf~\cite{Tes2001} for deep ocean waves. Therefore the surface is represented as a~heightfield, specifically, as a~wave height field computed for each horizontal point (in Equation~\ref{eq:tess}) on a~synthesized patch of ocean using IFFT. The Phillips spectrum, as proposed by Tessendorf~\cite{Tes2001}, is also used. It is a~semi-empirical model for wind-driven waves in a~fully developed sea based on statistical analysis of observations of the ocean surface. The spectrum sampled in the fourier domain using gaussian random number generator, with mean 0 and standard deviation 1. 

Apart from the wave height field, $h(\*{x},t)$, the slope vectors $\*{s}(\*{x},t)$ and the displacement vectors, $\*{D}(\*{x},t)$, are also computed using IFFT for each horizontal position $\*{x} = (x, z)$ and time point $t$. Then the three-dimensional position $\*{p}_w$ on a~grid on the water surface at $t$ after displacement is:
\begin{align}
    \*{p}_w = \left[\begin{array}{c}
        x + u D_x (\*{x}, t) \\
        h(\*{x}, t) \\
        z + u D_z (\*{x}, t)
    \end{array}
    \right]^T ,
\end{align}
where $u \in \mathbb{R}$ is a~user controlled parameter to scale the displacement vector, usually kept at value 1. 

The normals $\*{n}_w$ are computed using the slopes $\*{s}(\*{x},t)$. The slopes must also account for the displacement~\cite{Gam2018}:
\begin{align}
    \*{s}(\*{x},t) &= \left[
        \frac{
    \displaystyle 
            \frac{\text{d}h(\*{x},t)}{\text{d}x}
        }{
    \displaystyle 
        1 + u\frac{\text{d}D_x(\*{x},t)}{\text{d}x}
        },
        \frac{
    \displaystyle 
            \frac{\text{d}h(\*{x},t)}{\text{d}z}
        }{
    \displaystyle 
        1 + u\frac{\text{d}D_z(\*{x},t)}{\text{d}z}
        }
    \right] , \\
    \*{n}_w &= \frac{ (-s_x, 1, -s_z)}{||(-s_x,1,-s_z)||} ,
\end{align}
where $D_x$ and $D_z$ are partial derivatives of $\*{D}$ in $x$ and $z$ dimensions, respectively, these must also be computed separately.

Concerning optics, the water surface is treated as~collection of small planar facts, each tilted according to its normal. Then, the intersection of light ray with the surface at point $\*{p}_w$ is treated according to the equations and assumptions presented in Section~\ref{sec:air-water-interface}.



\subsection{Terrain Representation and Generation}

Ocean floor or seabed is modeled as a~terrain. The terrain must be represented in a~way that facilitates simple ray-surface intersection. This is accomplished through the combination of signed distance functions (SDFs), therefore, the terrain is a~signed distance field. Then the intersection of a~ray with the terrain can be found through the use of ray marching with distance estimation\,---\,sphere tracing\,---\,as described in Section~\ref{sec:proc-gen}. 

Quilez~\cite{Qui2019} discusses how the detail obtained in a~manner described by fractal Brownian Motion (fBM) can be used with SDFs. Traditional 3D noise is not a~distance function, it is a~signed field/function that does not measure distances. A~3D noise that is also SDF is obtained by summing multiple octaves/layers of sphere SDFs, in the following way:
\begin{enumerate}
    \item the base noise/function for fBM must be a~smooth SDF, e.g., big sphere,
    \item an infinite grid of spheres of random sizes is computed using domain repetition operator, each sphere has a~random radius that is smaller than half the edge-length of the grid,
    \item for a~given point in space, the SDF is evaluated for the 8\,spheres at each corner of a~grid cell,
    \item instead of summing, each additional layer is combined with the base SDF using the SDF operators, such as \emph{min, smoothmin}.
\end{enumerate}
Combination of multiple layers of spheres with increasingly smaller radii can produce highly detailed 3D terrain that is also SDF.


\subsection{Bio-Optical Model of Water Volume}
\label{sec:design-model}
The bio-optical model of the water volume is represented by the following inherent optical properties (IOPs): absorption coefficient, scattering coefficient, and phase function. Optical properties of the following constituents are considered: water, phytoplankton, and colored dissolved matter (CDOM, commonly known as yellow matter). Only elastic scattering is modeled. While both water and phytoplankton show elastic scattering, the elastic scattering produced by CDOM is shown to be negligible~\cite{OOWB2021-const-cdom} and is therefore ignored. 
% TODO cite https://www.oceanopticsbook.info/view/optical-constituents-of-the-ocean/colored-dissolved-organic-matter

The bio-optical properties are modeled to represent waters with high phytoplankton concentration, in which optical parameters can be correlated to chlorophyll concentration. It is assumed that the vertical distribution of the properties is non-uniform and that it correlates to chlorophyll concentration. The vertical distribution of \emph{chlorophyll concentration} $C$~at depth $z$~is approximated using a~Gaussian function, as employed by Ma et al.~\cite{MXWTC2016}:
\begin{align}
    C(z) = C_0 + \frac{h}{\sigma \sqrt{2\pi}} exp \left[
        - \frac{1}{2} \left( \frac{z - z_{max}}{\sigma} \right)^{2}
    \right] ,
\end{align}
where $C_0$ is the background chlorophyll concentration, $z_{max}$ is the depth of the maximal chlorophyll concentration, $\sigma$~is the standard deviation, and $h$~determines the amplitude of the maximal chlorophyll concentration. This vertical distribution of chlorophyll concentration is then used to determine the vertical distribution of other IOPs. 

\subsubsection{Absorption}
The \emph{total absorption coefficient} $\absorbCoef(\lambda, z)$ is composed as:
\begin{align}
    \absorbCoef(\lambda, z) = \absorbCoef^w (\lambda) + \absorbCoef^p(\lambda, z) + \absorbCoef^c (\lambda, z) ,
\end{align}
where superscripts $w$, $p$, and $c$ stand for water, phytoplankton and CDOM, respectively. 

The spectral absorption coefficient of water is modeled according to measurements of the absorption spectrum of pure water by Pope and Fry, 1997~\cite{PF1997}. Figure~\ref{fig:design-water-curves} shows the wavelength-dependent measurements of $\absorbCoef^w$ for the range wavelengths according to visible spectrum. The absorption in the blue (i.e., the blue part of the visible spectrum as shown in Figure~\ref{fig:visible-range}) is very low while increasingly rising to the maximum in red. 

The absorption coefficient of phytoplankton is parameterized as follows~\cite{MXWTC2016}:
\begin{align}
    \absorbCoef^p(\lambda, z) = 0.06 \alpha_c^{*}(\lambda) C(z)^{0.05} ,
\end{align}
where $\alpha_c^{*}$ is \emph{chlorophyll-specific absorption coefficient} of phytoplankton describing the absorption per unit of concentration for a~particular species of phytoplankton, given in m$^2$/mg. The wavelength-dependent measurements of $\alpha_c^{*}$ are based on the data provided by Bricaud et al., 1995~\cite{BBMC1995}. They provide data for absorption spectra of phytoplankton normalized by chlorophyll, shown in Figure~\ref{fig:design-phyto-abs}. The absorption by phytoplankton demonstrates high spectral variations with peak in blue, decreasing approximately exponentially to the red. 

Absorption coefficient of CDOM is parameterized as~\cite{MXWTC2016}:
\begin{align}
    \absorbCoef^c (\lambda, z) = 0.012 C(z)^{0.65} e^{-S_y (\lambda - 440)} .    
\end{align}
The coefficient is computed by fitting an exponential absorption curve where $S_y$~is the slope taken as a~constant, with a~value of $0.014\,$nm$^{-1}$.

% figures of absorption graphs
\subsubsection{Scattering}
The \emph{total scattering coefficient} $\scatterCoef(\lambda, z)$ is decomposed into:
\begin{align}
    \scatterCoef(\lambda, z) = \scatterCoef^{w}(\lambda) + \scatterCoef^p(\lambda, z) ,
\end{align}
where the scattering coefficient for water $\scatterCoef^{w}$ is modeled using measurements for pure seawater by Smith and Baker, 1981~\cite{SB1981}. The scattering coefficient phytoplankton $\scatterCoef^{p}$ is given by Gordon and Morel, 1983~\cite{GM1983}, which has shown that even small concentrations of phytoplankton contribute to the total elastic scattering in the water:
\begin{align}
    \scatterCoef^{p}(\lambda, z) = \left( \frac{550}{\lambda} \right) 0.30 C(z)^{0.62} ,
\end{align}
where the constant $0.30$\,was selected to fit the data collected from many water types, while the upper bound has a~value of $0.45$~\cite{GSMA2008}.

The phase function is composed of two scattering phase functions:
\begin{align}
    p(\theta, \lambda, z) &=
        \frac{\scatterCoef^w (\lambda)}{\scatterCoef(\lambda, z)} p_w(\theta) +
        \frac{\scatterCoef^p (\lambda)}{\scatterCoef(\lambda, z)} p_p(\theta) ,
\end{align}
where $p_w$ and $p_p$ are the scattering phase functions for water and for phytoplankton, respectively. The phase function for pure water is modeled using Einstein-Smoluchowski theory, as given by Gutierrez et al.~\cite{GSMA2008}:
\begin{align}
    p_w (\theta) = 0.06225 (1 + 0.835 \cos^2 \theta) .
\end{align}
The scattering phase function for phytoplankton is modeled using the \emph{Petzold average-particle function} based on measurements of VSFs by Petzold, 1972~\cite{Pet1972}. The average-particle phase function, which satisfies the normalization condition, was made as the average of the three VSFs measured in different waters (clear ocean, coastal ocean, and turbid water) and is used as a~reference particle phase function~\cite{OOWB2020-scat-petz}.
% TODO cite https://www.oceanopticsbook.info/view/scattering/petzolds-measurements
Figure~\ref{fig:petzold-average} shows the average phase function, in general, forward scattering dominates over scattering in backward direction.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.0\linewidth]{figures/design/PopeFry97.pdf}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.0\linewidth]{figures/design/SmithBaker81-scatter.pdf}
\end{subfigure}
\caption{Left image shows the absorption coefficient of pure water measured by Pope and Fry, 1997~\cite{PF1997}. The right image shows the scattering coefficient of seawater provided by Smith and Baker, 1981~\cite{SB1981}.}
\label{fig:design-water-curves}
\end{figure}

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.0\linewidth]{figures/design/BricaudEtAl95.pdf}
    \caption{The chlorophyll-specific absorption coefficient of phytoplankton describing the absorption per unit of concentration for a~particular species of phytoplankton, measured by Bricaud et al., 1995~\cite{BBMC1995}.}
    \label{fig:design-phyto-abs}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.0\linewidth]{figures/design/PetzoldAvg.pdf}
    \caption{Log-linear plot of Petzold average-particle phase function, used as phase functions for phytoplankton and particles. The function is the average of phase functions measured in different types of waters by Petzold, 1972~\cite{Pet1972}.}
    \label{fig:petzold-average}
\end{subfigure}
\label{fig:design-curves}
\end{figure}

% TODO final table of equations, similar to GSMA2008 or Hil2016

% The concentration of the chlorophyll in the water usually ranges from 0.01 mg/m 3 for open waters to 100 mg/m 3 . The spec- tral absorption coefficient of the phytoplankton is usually ex- pressed as a function of this concentration C as:~\cite{GSMA2008}

%Scattering by water
%  R. C. Smith, K. S. Baker, "Optical properties of the clearest natural waters (200-800 nm)," Appl. Opt.,20, 177-184 (1981). 

%\todo{Implementation: Sampling of the spectrum}
% https://analytik.co.uk/wp-content/uploads/2019/06/light-visible-spectrum-wavelength.png
% rgb
% 645 nm, 526 nm, 444 nm
% unreal takes centimeters
% spectrum -> RGB https://physics.stackexchange.com/questions/196747/converting-spectra-to-rgb-colour


% ================================================================================
% ================================================================================
\section{Scene Rendering}

%The scene is composed of 
%
%- what must be addressed
%    - how does the light flow, repeat if you must
%- how it is addressed
%- refs to subsections

% ---------------------------------------------------------------------------------
% TODO remove this subsection just make this the real into
%\subsection{Rendering Water Volume}
% Participating Media
%\section{Real-Time Scattering in Water Volume}

Using the designed scene representation and the designed bio-optical model, a~technique is proposed to render the scene.

A~bio-optical model, proposed in Section~\ref{sec:design-model}, tries to reproduce the optical parameters of Case\,1 waters. This is accomplished by defining a~vertical distribution of chlorophyll concentration that is used to determine all the other parameters. The volume rendering equation (VRE) (in Equation~\ref{eq:vre}) must be modified to reflect this depth dependence on the optical parameters. Moreover, it is not feasible to compute the VRE directly especially if interactivity is to be reached. Therefore the VRE must be approximated appropriately for this types of participating media.

Starting from the VRE for non-emissive media (henceforth just VRE),
\begin{align}
    L(\*{x}, \omega) &=
    \underbrace{
        T(\*{x},\*{x}_s) L_s(\*{x}_s, \omega)
    }_\text{surface term}
    +
    \int_{t=0}^{S} T(\*{x},\*{x}_t)
    \underbrace{
        \scatterCoef(\*{x}) L_i(\*{x}_t, \omega)
    }_\text{in-scattering term}
    \text{ d}t ,
    \label{eq:vre-noemis}
    %L_i (\*{x}, \omega) = \scatterCoef(\*{x}) \int_{\Omega} p(\*{x}, \omega, \omega_i) L(\*{x}, \omega_i) \text{ d}\omega_i  ,
\end{align}
by utilizing the properties of scattering in waters the VRE is approximated specifically for rendering underwater. It is approximated by two independent contributions: single scattering and multiple scattering. Then the radiance entering the camera at position $\*{c}$ from direction $\omega$, as illustrated in Figure~\ref{fig:design-final}, is then given by:
\begin{align}
    L(\*{c}, \omega) &= L_{sscat}(\*{c},\omega) +
                        L_{medium}(\*{c},\omega) .
\end{align}
The single scattering part $L_{sscat}$ models the contributions of the surface term and single bounce of the in-scattering term. Specifically, it accounts for: 1) $L_s$ the radiance at the surface point $\*{x}_s$, as a~result of the exitant radiance $L_o$ and one-bounce surface caustics at that point; 2) the in-scattering integral due to volume scattering is approximated by volumetric caustics method. The single scattering part is further described in Section~\ref{sec:design-ss}.

The multiple scattering part $L_{medium}$ models the contributions from the media by approximating the in-scattering integral based on the general scattering properties of waters. It also accounts for the non-uniform vertical distribution of bio-optical properties. This part is described in Section~\ref{sec:design-ms}.

% TODO figure how it is computed
\begin{figure}[H]
    \centering
\begin{subfigure}{.8\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/design/singlescatteringDesign.pdf}
\end{subfigure}
\caption{Radiance entering the camera at position $\*{c}$ is computed by sampling along the ray from the camera to the closest intersection point $\*{p}$ on a~surface. Each sample accumulates transmittance and evaluates contribution due to multiple scattering $L_{medium}$. The dark blue arrows show contributions due to single scattering at each sample, this is explicitly computed using techniques that simulate caustics.}
\label{fig:design-final}
\end{figure}

% single scattering
\subsection{Single Scattering Approximation}
\label{sec:design-ss}

%- problems with traditional single scattering, as used in the manner by Wronski~\cite{Wro2014} and Hillaire~\cite{Hil2015}
%    - ref to figure 
%- cannot be used due to the refraction - water surface interface as illustrated in figure
%    - the goal is to approximate in-scattering by taking a~single bounce
%    - Wronski and Hillaire do that by ray marching along the direction casting rays towards the light source - a backward like ray tracing manner
%- assume only one directional light in the scene - the sun, this approach cannot be done, because it is not possible to know which refracted directions contribute to the point using this backward like ray tracing approach
%    - there can be multiple incident light rays that get refracted to the same point
%    - would have to try all the directions or make crude approximations as certain methods
%
%- instead of this we can take another route, approximating this part not as part of single scattering as a~whole but by specialized techniques for rendering caustics that work in a~forward ray tracing like manner
%
%- physics: due to the curvature of the water surface, change in IOR - focusing and defocusing of light, produces light streaks, godrays that intersect with the terrain or seabed to form caustics
%    - ref to the final figure
%    - must be addressed separatelly -- ref to sections
%    - the godrays, called volumetric caustics, and the intersection of the beams with the terrain --- surface caustics
%
%Final construciton:
%- once these two are addressed, they solve the single bounce, what is left is to account for the exitant radiance from a point on the terrain,
%    - simplification - assumed to be diffuse, Lambertian surface
%    %- described in section (todo needed?) 
%    - ray marching to account for the bio-optical model - compute transmittance
%        - to each point on the terrain
%
%- two issues with the banding:
%    - ray marching, or sphere tracing step of the underlying terrain
%    - transmittance being accumulated at these discrete steps
%    -> accumulate transmittance in a~texture and blur it
%
%- light scattered from a~punctual light source, i.e., a~light source represented by a~single infinitesimal point
%- ref to image showing how VRE works
%
%From below:
%This section 
%Builds on the related works for participating media rendering 
%problems that must be addressed:

The single scattering approximation, as computed by Wronski~\cite{Wro2014} and Hillaire~\cite{Hil2015} cannot be computed in this scenario directly. The refractive  interface, represented by water surface, causes the incident light rays to be refracted and additionally, due to the surface being curved, the refracted rays point in different directions. This causes the rays to focus and defocus in various ways, thus forming caustics (as described in Chapter~\ref{cha:caustics}).

In the single scattering approximation, the in-scattering integral is computed by ray marching towards the light source for each sample along the view ray. However, for this scenario, there might be multiple light rays contributing to a~sample from multiple refracted directions of the incident light, or there might not be any direction that contributes directly (a~defocused region). Therefore, the in-scattering integral, as it is accumulated in the single scattering approximation, cannot be computed in this backward ray tracing manner, i.e., by marching the rays from the camera.

Offline approaches can compute this by tracing a~large number of rays for each view sample to the water surface, effectively accumulating contributions from all possible directions. There are also real-time approaches that can approximate the in-scattering integral by making aggressive assumptions about the incident light and the interface. All of the aforementioned approaches are not feasible for the proposed representation of the scene. Instead of computing the in-scattering as part of the single scattering approximation, there are specialized techniques that approximate both surface and volumetric caustics\,---\,the result of the in-scattering integral.

% TODO figure design-ss

The single scattering approximation, in Equation~\ref{eq:design-ss}, is then split into parts, where each part is computed separately, in the following way: the surface term, and the in-scattering integral. The surface term is computed the same way as in the single scattering approximation. The evaluation of the in-scattering integral in the volume results in light shafts or godrays, called \emph{volumetric caustics} in underwater environments. The intersection of volumetric caustics with the scene produces surface caustics. These both effects, though part of the same phenomenon, are then simulated explicitly using specialized techniques that target real-time.
\begin{align}
    L_{sscat}(\*{c},\omega) &= T(\*{c}, \*{p})
                L_o(\*{p},\omega) + L_{surfCaustics}(\*{p},\omega) + L_{volCaustics}(\*{c},\*{p},\omega) ,
    \label{eq:design-ss}
\end{align}


\subsubsection{Surface Term}
The surface term is computed by ray marching along the view ray accumulating transmittance until a~point $\*{p}$ on a~surface is reached (as illustrated in Figure~\ref{fig:design-final}). Due to the proposed bio-optical model with non-uniform vertical distribution of media properties, transmittance cannot be computed analytically and must be iteratively accumulated by sampling along the ray.
\begin{align}
    L_{sscat}^{surf}(\*{c},\omega) =
        e^{-\int_{\*{x}=\*{c}}^{\*{p}} \extinctCoef(z) ||\text{d}\*{x}||}
        L_o(\*{p}, \omega) ,
\end{align}
where $z$ is the depth below the water surface of the sample. $L_o$ is the exitant radiance at surface point $\*{p}$. Underwater surfaces are modeled as diffuse. Since there is a~lack of measurements and models of the BRDF for actual ocean bottom materials, the bottom being a~Lambertian reflecting surface is a~typical and reasonable assumption~\cite{OOWB2021-lambert}. 
% TODO cite https://oceanopticsbook.info/view/surfaces/lambertian-brdfs
On the other hand, the water surface is modeled as diffuse with specular highlights using \emph{Blinn-Phong reflection model}~\cite{Bli1997}. The exitant radiance from point $\*{p}_{w}$ on the water surface and from point $\*{p}_t$ on the terrain towards the view direction $\omega$ is evaluated, similarly to Baboud et al.~\cite{BD2006}:%(as illustrated in Figure~\ref{fig:design-L_ws}:
% TODO cite James F. Blinn (1977). "Models of light reflection for computer synthesized pictures". Proc. 4th Annual Conference on Computer Graphics and Interactive Techniques: 192–198. CiteSeerX 10.1.1.131.7741. doi:10.1145/563858.563893. S2CID 8043767
\begin{align}
    L_o^{w}(\*{p}_{w}, \omega) &= (1 - R_F) L_t(\*{p}_t, \*{t}) + R_F L_r(\*{p}_{w},\*{r}) , \\
    %L_r(\*{p}_{w},\*{r}) &= (\*{n}_{w} \cdot \*{H})^{32} (\*{n}_w \cdot \*{l}) L_{sun} L_{sky}(\*{r}), \\
    L_t(\*{p}_t,\*{t}) &= T(\*{p}_w,\*{p}_t) L_o^{t}(\*{p}_t, \*{-t}) L_{surfCaustics}(\*{p}_t) 
        % TODO + \backScatter L_
    , \\
    L_o^t(\*{p},\omega) &= A_t(\*{p}) (\*{n}_t \cdot \omega) ,
    %\*{H} &= \frac{\*{l} + \omega}{||\*{l} + \omega||} ,
\end{align}
where $R_F$ is Fresnel reflectance, $L_t$ is the transmitted radiance from point on the terrain $\*{p}_t$ underwater along the refracted direction $\*{t}$, and $L_r$ is the reflected radiance from the reflected direction $\*{r}$ at the point on the water surface, and
%$\*{H} = \frac{\*{l} + \omega}{||\*{l} + \omega||}$ is the half-vector
$A_t$ is the albedo or reflectance of the terrain.

 %TODO figure: fig:design-L_ws


%- single scattering approximation 
%From above:
%- apply do model by Baboud
%- the corresponding water surface position and normal is sampled
%- the refracted and reflected ray is computed
%Blinn-Phong model:
%- the amount of radiance along the reflected ray is by sampling the sky model
%- the amount of radiance along the refracted ray is computed as if underwater
%- combined using Fresnel coefficients


%\subsubsection{Terrain Rendering}

\subsubsection{Surface and Volumetric Caustics}
\label{sec:design-caustics}

Both surface and volumetric caustics are computed as part of one algorithm, thus reusing already computed information, such as intersections. The technique reuses ideas in ray-guided volumetric caustics by Gruen~\cite{Gru2019v}, and in \emph{procedural caustics mesh} (PCM) method by Yang et al.~\cite{YO2021}. The approach are photon mapping-based and in-scattering integral is approximated using beam tracing. It works in the following steps:
\begin{enumerate}
    \item A~Photon regular grid is created, a~triangular mesh independent from the geometry of the water surface.
    \item Caustics map is created by rendering the water surface from the point of view of the light, storing positions and normals.
    \item Refracted caustics map is then computed: each pixel [x,y] in the caustics map is treated as a~ray, and is refracted along the surface. Each ray is then traced in the scene until an intersection is found, which is guaranteed. The intersection is found by evaluating the terrain SDF using sphere tracing. The world space position of the intersection is stored at the same pixel location [x,y] in the refracted caustics map. Additionally, the radiance contribution of the photon is stored. % TODO analytically??
    \item Compression ratios are computed for each triangle of the photon grid by projecting each photon to the caustics map, to retrieve the specular triangle, and to the refracted caustics map, to retrieve the caustic triangle. The ratios are then stored in a~buffer.
    \item Surface caustics are accumulated in surface caustics map (SCM) in the similar manner as Meenrattanakorn et al.~\cite{ML2022}: by storing the photon grid in a~vertex buffer and drawing it either as point primitives, or as triangles, addively blending the contributions. 
    \item Volume caustics are computed by drawing the photon grid as a~triangular mesh and projecting each photon to the caustics map, thus sampling the positions and normals of the water surface. Each triangle is tessellated based on the compression ratio and the on-screen size of its caustic triangle (obtained by projecting each photon to the refracted caustics map). From the tessellated vertices a~ray is generated to intersect the terrain again thus forming new specular and caustic triangles. These are then built into the triangular beam, which has its frontface rasterized, additively blending computed contribution per pixel.
\end{enumerate}

Both the surface caustics and volumetric caustics produce noise, and therefore must be denoised.



\ifdefined\NOTSUREYET
- directional light source representing the sun
    - assumption - case 1 waters - incident angle near normal incidence
    - sun is high above

- regular photon grid
    - decoupled from the geometric complexity of the water surface
- first hit - emission
    - specular and refractive surface, water surface
    - no occlusion from light view
    - water surface textures are projected to the photon grid

    - set up refracted rays, compute flux:
    
radiant intensity emitted into direction = flux * max(0, dot(n,l))

according to Kim~\cite{Kim2019} emission flux of a~infinitesimal directional light is
\begin{align}
    \Phi_e = l_e \frac{Resolution(Tex_{ws})}{N} ,
\end{align}
where $l_e = \text{light color} \cdot \text{light intensity}$ is the radiance from the light source.
where $N$ number of samples
% ML2022


- goal: trace only the photons that influence caustics, as Meenrattanakorn et al.~\cite{ML2022}:
    compute luminance using the flux:
\begin{align}
    lum = R_T \Phi_e \text{max}(0, \*{n}_{ws} \cdot \*{l}) ,
\end{align}
where $\*{n}_{ws}$ is the normal of the water surface at that point,
$\*{l}$ is the direction to the light source.
$R_T$ is the Fresnel transmittance, which takes the  IOR of both the air and water
- any photons with luminance less than 0.04 are discarded
    \todo{what does that mean?}

2. Photon Tracing

3. Photon 
- compute irradiance stored in a~pixel (Meenrattanakorn et al.~\cite{ML2022}):
\begin{align}
    E_p = \frac{\Phi_e}{A_p} (-\*{d} \cdot \*{n}_{t})
\end{align}
where $A_p$ is the area of pixel in world space

saved into a photon buffer: as tuple of clip-space coordinates, and the calculated irradiance

the buffer is used as a~vertex buffer to render each photon as a~point primitive, resulting in a~2D photon map

- additive blending with the scene
- make sure that the pixel is the front-most using camera depth buffer
\fi

\subsection{Multiple Scattering Approximation}
\label{sec:design-ms}

Multiple scattering is approximated based on the rationale by Monzon et al.~\cite{MAGM2022}. Multiple scattering is caused by the light being scattered a~number of times between particles before reaching the eye and it plays a~crucial role in underwater illumination. It is assumed to be very low frequency and so it can be modeled as a~function of depth and wavelength. Due to the designed bio-optical model (in Section~\ref{sec:design-model}), which assumes non-uniform vertical distribution of bio-optical properties, the proposed analytical solution by the authors cannot be used directly therefore a~solution that is more suited must be derived.

The contribution from the medium (based on VRE for non-emissive media in Equation~\ref{eq:vre-noemis}) due to multiple scattering is defined as:
\begin{align}
    L_{medium}(\*{c},\omega) &=
        \int_{t=0}^{S} T(\*{c}, \*{x}_t) L_i(\*{x}_t,\omega) \text{ d}t ,
        \label{eq:medium-vre} \\
    L_i (\*{x}, \omega) &= \scatterCoef(\*{x}) \int_{\Omega} p(\*{x}, \omega, \omega_i) L(\*{x}, \omega_i) \text{ d}\omega_i  ,
\end{align}
where $\*{c}$ is the view camera position, $\omega$ is the direction toward the view, and $S$~is the distance to the closest intersection along $\omega$ (as illustrated by Figure~\ref{fig:vre}).  

Assuming an isotropic media, where light is distributed uniformly in the sphere, the phase function becomes constant $p(\*{x},\omega,\omega_i) = \frac{1}{4\pi}$. Furthermore, the in-scattering integral describing radiance from every direction\,---\,the incident irradiance\,---\,is approximated by the vertical varying incident downwelling irradiance at the surface of the water volume $L_{ws}(\lambda)$ that relates to the diffuse attenuation coefficient for downwelling irradiance $\diffuseAttenDown(\lambda)$. This is a~fair assumption for general use case, since most of the downwelling irradiance underwater comes from a~light cone directly above due to so-called Snell's window (shown in Figure~\ref{fig:snellswindow}), even under different zenith angles~\cite{MAGM2022}. The incident irradiance due to multiple scattering will be constant at every depth $z$~below the surface as:
\begin{align}
    L_i(\*{x},\omega) = L_{ws} e^{-\diffuseAttenDown(\lambda)} z .
\end{align}

As described in Section~\ref{sec:iops-aops}, $\diffuseAttenDown(\lambda)$ is a~physically measurable coefficient that measures the irradiance attenuated in the water volume per meter per wavelength. Though it could be obtained from tabulated data (e.g., by Solonenko and Mobley, 2015~\cite{SM2015}) for case\,1 waters, it would not account for the non-uniform vertical distribution of bio-optical properties posed by the designed bio-optical model, therefore $\diffuseAttenDown$ must be computed. 

There are various relations between absorption coefficient and backscattering coefficient $\backScatter$ that approximate $\diffuseAttenDown$. According to Lee et al., 2005~\cite{LDCDSR2005} the value of $\diffuseAttenDown$ can be estimated using the following relation that is based on analysis and numerical simulations of the radiative transfer equation:
\begin{align}
    \diffuseAttenDown(\lambda, z) &= m_0 \absorbCoef(\lambda, z) +
                                    \nu \backScatter(\lambda, z),
\end{align}
where the parameters $m_0$~and $\nu$~describe the contributions of $\absorbCoef$ and $\backScatter$~to $\diffuseAttenDown$, respectively. It was found that $m_0 approx (1 + 0.005 \theta_z)$ where $\theta_z$~is the solar zenith angle above the water surface, and $\nu$~can be expressed as:
\begin{align}
    \nu = m_1 (1 - m_2 e^{-m_3 \absorbCoef(\lambda, z)}),
\end{align}
where the values of $m_1$, $m_2$, and $m_3$ are 4.18, 0.52, and 10.8, respectively. To obtain the value of backscattering coefficient in waters where phytoplankton and other particulate are dominant, Morel, 1988~\cite{Mor1988} introduced the following relations:
\begin{align}
    \backScatter(\lambda, z) &= 0.5 \scatterCoef^w(\lambda,z) +
                                B_b \scatterCoef^{p}(\lambda,z) , \\
    B_b &= 0.002 + 0.02 (0.5 - 0.25 \log_{10} C(z))
\end{align}

Now, having all the necessary terms, the contribution from the medium, from Equation~\ref{eq:medium-vre}, is expanded step by step (omitting any wavelength-dependency for clarity):
\begin{align}
        \nonumber
    L_{medium}(\*{c},\omega) &=
        \int_{t=0}^{S} T(\*{c}, \*{x}_t)
            \scatterCoef(\*{x}) \int_{\Omega} p(\*{x}, \omega, \omega_i) 
                L(\*{x}, \omega_i)
            \text{ d}\omega_i
        \text{ d}t  \\
        \nonumber
     &= \int_{t=0}^{S} T(\*{c}, \*{x}_t)
            \scatterCoef(z) \frac{1}{4\pi} \int_{\Omega}
                L(\*{x}, \omega_i)
            \text{ d}\omega_i
        \text{ d}t  \\
        \nonumber
     &= \int_{t=0}^{S} T(\*{c}, \*{x}_t)
            \scatterCoef(z) \frac{1}{4\pi} L_{ws} e^{-K_d(z) z}
        \text{ d}t  \\
        %\nonumber
     &= \frac{1}{4\pi} L_{ws} 
            \int_{t=0}^{S} T(\*{c}, \*{x}_t) \scatterCoef(z) e^{-K_d(z) z}
        \text{ d}t ,
        \label{eq:medium-exp}
\end{align}
where $z$~is the depth, or distance to the water surface along the vertical axis, of each sample.

Furthermore, the proposed energy-conserving analytical integration by Hillaire, 2015~\cite{Hil2015} is applied. It analytically integrates the scattered light over a~range according to both $\extinctCoef$ and a~scattered light sample $I$, as well as an integration depth. Applying this proposed analytical solution to the integration in Equation~\ref{eq:medium-exp}, for a~light sample $I(z) = \scatterCoef(z) e^{-K_d(z) z}$, gives the final form of the contribution from the medium:
\begin{align}
    L_{medium}(\*{c},\omega) &= \frac{1}{4\pi} L_{ws} 
            \int_{t=0}^{S} \frac{I(z) - I(z) e^{-\extinctCoef(z) ||\*{c} - \*{x}_t||}}{\extinctCoef(z)}
        \text{ d}t .
    \label{eq:medium}
\end{align}
This integration requires less number of samples compared to the trivial sampling of the integral to make participating media look correct with increasing density (i.e., extinction).


% perhaps before the single scattering?
%\subsection{Terrain Rendering}


\ifdefined\DESIGNCAUSTICS
\section{Caustics}
\subsection{Surface Caustics}

1. Computing the intensity

Emission flux of a~directional light:
\begin{align}
    l_e = light color \times light intensity \\
    \Phi_e^D = l_e \frac{a_l}{p_w p_h}
\end{align}
each photon has an area $1 / (p_w p_h)$ for light area $a_l$

The flux is saved as irradiance on the pixel:

where $A_p$ is the world-space area of a~pixel $p$:
\begin{align}
    A_p = 
\end{align}
\cite{Kim2019}, \cite{ML2022}

the brightness value is scaled by:
    - compression ratio
    - the amount of light that has been absorbed by the length of the way the ray travels through the water volume \cite{BD2006}
\cite{Gru2019s}

Problem with big open scenes: pixels in the distance capture many more rays from the caustics map
, hence they get a~lot brighter:
    1) estimate the world-space size of each pixel
    2) divide the incoming energy by this area to normalize brightness
for volumetrics:
    bring down the slicing complexity in the distance
\cite{Gru2019v}

\subsection{Volumetric Caustics}
\fi



% ================================================================================
% ================================================================================
% ================================================================================
\chapter{Implementation and Evaluation of Technique for Rendering Underwater Scenes}
\label{ch:impl}

The designed technique from Chapter~\ref{ch:design} was not implemented in its entirety. Only a~subset of techniques was implemented and experimented with. A~standalone C++ application with the use of Vulkan API~\cite{vulkan} was developed for this purpose. In Figure~\ref{fig:gallery} below is a~gallery of screenshots of experiments with those techniques. The application implements deferred shading pipeline and the following components:
\begin{itemize}
    \item Water Surface. Was implemented as designed. It is represented as a~triangular mesh\,---\,a~regular grid of 256$\times$256 vertices that is displaced. There are two textures computed each frame: displacement map and normal map (with the same resolution as the grid). The textures are computed  on CPU-side using the FFTW library~\cite{fftw}. The execution was parallelized and runs on all available threads. Uses both, the parallel execution of FFTW plans, and OpenMP library to parallelize for loops.
    \item Camera. An FPS camera was implemented. It uses Euler angles.
    \item Sky Model. Implements the Preetham's sky model~\cite{PSB1999}.
    \item Camera G-Buffer. Has the following attachments: world space positions, world space normals, albedo, and depth. 
    \item Terrain. Simple combination of multiple SDFs, possible to use multiple octaves of fBM SDF.
    \item Surface caustics. Implemented as described.
    \item Volumetric caustics. Implemented partially, without the tessellation part.
    \item The multiple scattering approach was implemented fully and tested.
\end{itemize}

Furthermore, the application implements its own Vulkan abstraction. And additionally uses the following libraries:
\begin{itemize}
    \item Dear ImGuI~\cite{imgui}: to expose user-controllable parameters in GUI form,
    \item CPU side profiler,
    \item Shaderc library~\cite{shaderc}: allows for compilation of GLSL shaders at runtime.
\end{itemize}

The application can be compiled using CMake. The instructions are given in a~readme text file in the root folder.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/results/water-above.png}
    \caption{Water surface rendered with camera above the surface.}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/results/terr-above.png}
    \caption{Terrain refracted from above the water surface.}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/results/terr.png}
    \caption{Terrain below the water surface, rendered using the computer vision algorithm model by Sedlazeck et al.~\cite{SK2011}.}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/results/oceanic-terr.png}
    \caption{Terrain below the water surface, rendered using the computer vision algorithm model by Sedlazeck et al.~\cite{SK2011}.}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/results/snellswindow.png}
    \caption{Snell's window below the water surface, rendered using the computer vision algorithm model by Sedlazeck et al.~\cite{SK2011}.}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/results/multiple.png}
    \caption{SDF terrain underwater rendered using the designed multiple scattering approach and the proposed bio-optical model.}
\end{subfigure}
\caption{Screenshots of experiments with various techniques and methods.}
\label{fig:gallery}
\end{figure}

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/results/caustics.png}
    \caption{Surface caustics rendered on a~terrain represented as a~plane. Caustics show blocky appearance and noise.}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/results/multiple-scattering-w-caust.png}
    \caption{Surface caustics rendered on the SDF terrain combined with the designed multiple scattering approach.}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/results/beams.png}
    \caption{Volumetric caustics using the beam tracing method intersecting terrain represented by a~simple horizontal plane.}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=1.\linewidth]{figures/results/combined.png}
    \caption{The resulting combined scene: water surface, surface caustics on SDF terrain and volumetric beams.}
\end{subfigure}
\caption{Continuation of screenshots of experiments with various techniques and methods.}
\label{fig:gallery1}
\end{figure}





%\chapter{Evaluation}

\chapter{Conclusion}

The goal of this work is was study, design and implement an application that is able to render underwater environments in real time. Existing solutions for real-time rendering and simulation of oceans were extensively investigated. Furthermore, techniques that simulate underwater caustics were investigated. 

A~technique was designed to render underwater scenes in high-quality and in real-time. Single scattering and multiple scattering approximations were designed that use ray marching to render water volumes. A~physically-based bio-optical model for oceanic waters was designed. This model defines non-uniform vertical distribution of bio-optical parameters. The multiple scattering approximation was specifically designed by exploiting the properties of waters and accounting for the proposed bio-optical model. 

The designed technique was not implemented in its entirety. Only a~subset of techniques were implemented and experimented with. The implemented techniques run in real-time on relatively modern hardware. Research and investigation of existing techniques and the subsequent designed was prioritized instead.


% results

% implementaion

% future



%=========================================================================

\ifdefined\SHOWCAUSTICSMETHODS
%\subsection{Gruen. Real-Time Ray-Traced One-Bounce Caustics}
goal: investigation how real-time ray tracing can simplify and innovate on current methods for rendering real-time caustics
- assumes scenes with DXR bounding volume acceleration structures already set
    - see Section about implementation details in the article.
- this allows to find the specular intersections of refracted and reflected light rays and a~dynamic 3D scene
    - in many games: avoid finding intersections with the full game scene, only compute intersections with a~ground plane or the inside of a~3D box
    - for shadows, that are the result of many incoherent refracted/reflected light rays, shadow map is used, - that is not physically correct - it does not reflect how in reality the water surface moves

- conceptually the approach is more similar to photon mapping~\cite{Jen2001}
does not suffer from the limitations of image space methods:
    - they find approximate intersections
    - the refracted/reflected light rays may be occluded in both the primary view (camera view) and the view from light \todo{figure 1.4, rather should be put into overview}
- there are some mitigation strategies but those lead to more implementation complexity and may have higher runtime cost and memory consumption over this DXR based approach

- also there are approaches that voxelize the receiving geometry and march the voxelized 3D grid
    - while, these can yield good results, depending on the resolution of the grid, voxelization is usually not cheap, also raymarching the 3D grid also has a~high performance cost.
    - the implementation complexity of such techniques is above this approach

- can compute refracted and reflected caustics, 
- assumes that the refractive/reflective interface (from now on, just interface) is represented as a~mesh with certain complexity in number of vertices
Algorithm overview
divided into x steps:
\begin{enumerate}
    \item \textbf{Compute a~refracted/reflected mesh.}
        - construct a~refracted/reflected ray $R$~for each vertex of the interface
        with the origin at the vertex's position and with the direction pointing along the refracted/reflected direction of incident light.
        - Intersect each ray~$R$ with the receiving geometry. The intersections are the positions of vertices of the resulting refracted/reflected mesh. 
        - The resulting mesh has the same number of vertices and the same triangle count as the interface. Figure~\ref{fig:Gru2019-refracted-mesh} illustrates this step for refracted rays.

    \begin{figure}[ht]
        \centering
        \includegraphics[width=0.7\textwidth]{figures/Gru2019s/refracted-mesh.png}
        \caption{Computed refracted mesh from refracted rays starting at the vertices of interface's triangles and intersecting the scene~\cite{Gru2019s}.}
        \label{fig:Gru2019-refracted-mesh}
    \end{figure}
        
    \item \textbf{Compute a~compression ratio for each refracted/reflected triangle.}
        For each triangle of the interface there is a~corresponding triangle in the refracted mesh. The refracted triangles can have either bigger or smaller area than their respective interface triangle, corresponding to the focusing and defocusing of the refracted light within a~triangular beam, as shown in Figure~\ref{fig:Gru2019-beams}. The same applies for reflection. The ratio of such areas for each refracted/reflected triangle is computed as:
        \begin{align}
            compression ratio = \frac{ area(\triangle_{interface}) }
                                     { area(\triangle_{refracted}) }
        \end{align}
        
    \begin{figure}[ht]
        \centering
        \includegraphics[width=0.7\textwidth]{figures/Gru2019s/focus-light.png}
        \caption{Refracted incident light forms beams that focus or defocus light~\cite{Gru2019s}.}
        \label{fig:Gru2019-beams}
    \end{figure}
        
    \item \textbf{Generate refracted/reflect light rays.}
        A~number of samples is generated for each interface's triangle. The higher the refractive/reflective compression ratio, the more samples are generated. Also the on-screen size of the refracted/reflected triangle is considered, so that more samples get generated for triangles closer to the viewer.
        Each sample is represented as a~random world-space position in the interface's triangle that can be seen from the light source. 
        A~ray, or a~photon, is emitted from this position along the the interpolated refracted/reflected direction. The visibility of the sample can be determined by a~shadow map or by shooting a~ray towards the light source.

\begin{figure}[H]
\centering
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{figures/Gru2019s/sample-generation.png}
    %\caption{Environment map of the sun directly above.}
    %\label{fig:GS2004-sun-angle}
\end{subfigure}
\begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width=.5\linewidth]{figures/Gru2019s/three-triangles.png}
    %\caption{Intensity of the sun disk based on the angle of incidence}
    %\label{fig:GS2004-graph-intensity}
\end{subfigure}
\label{fig:GS2004-environment-map}
\caption{Environment map of the sun directly above and the way it is sampled using the normal of the wave equation~\cite{GS2004}.}
\end{figure}
        
    \item \textbf{Intersect and accumulate brightness.} 
        An intersection of each such ray with the receiving scene geometry is computed. Each intersection is then projected to screen space. If the position is on screen and corresponds to the front-most pixel in that position, then a~brightness value is accumulated in that screen location. The brightness value can be scaled by the compression ratio, and/or by the amount of light that has been absorbed as the ray traveled through a medium like water (based on Baboud and Décoret~\cite{BD2006}).
    \item \textbf{Denoise the resulting buffer from the previous step.} This results in denoised and smooth caustics buffer. A~set of iterated cross-bilateral blurring steps is used. It accounts for differences in view-space depth, normal and view-space positions.
\end{enumerate}
The caustics buffer is then used to shed light on the receiving scene.

To generate the positions from Step ?? a~simple seed based random number generator is used. The generator produces ``random'' barycentric positions inside an interface's triangle based on its ID/index as the seed to make sure that the generated noise is table across a~triangle. Any change in the number of samples for a~triangle between frames then only results in addition or removal of samples while the positions in the triangle remain the same.
\todo{keep? then cite}
\lstset{style=c++}
\begin{lstlisting}[language=c++]
float3 RandBarycentrics(inout uint randSeed)
{
    float s = RandFloat(randSeed);
    float t = RandFloat(randSeed);

    [flatten]if (s + t > 1.0f)
        s = 1.0f - s, t = 1.0f - t;
        
    return float3(1.0f - s - t, s, t);
}
\end{lstlisting}

Limitations, what to consider:
- The refracted/reflected mesh from Step~1 needs to be detailed enough to allow for the computation of a~high enough quality compression ratio, it does not need to be as fine as possible to follow every detail of the receiving geometry. 
- In Step ?? the rays are formed for each interface's triangle within an enlarged view frustum. That is done to also account for refracted/reflected rays that enter the view frustum from outside the normal view frustum. 
- Some refracted/reflected triangles may still end up off-screen even if some of the generated rays cause on-screen caustics. Then the maximum screen-space size of the refracted triangle and three additional triangles is taken. These three triangles are formed by the rays of a~refracted/reflected triangular beam at the intersection points at different distances, as depicted in Figure~\ref{fig:Gru2019-three-triangles}.
- To speed up the denoising, the brightness can be accumulated into a~half-resolution buffer.
- It is possible to extend the algorithm to incorporate multi-bounce caustics. Secondary rays can be cast from intersections of the refracted rays with the scene. Then for each new generated intersection of a~secondary ray, the brightness is accumulated at the projected screen space position of the intersection.

The results 
- raw caustics buffer without denoising
- after denoising has been applied

\todo{figure 1.9, 1.10}

\begin{figure}[H]
\centering
\begin{subfigure}{.7\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{figures/Gru2019s/result-raw.png}
    %\caption{Environment map of the sun directly above.}
    %\label{fig:GS2004-sun-angle}
\end{subfigure}
\begin{subfigure}{.7\textwidth}
    \centering
    \includegraphics[width=.9\linewidth]{figures/Gru2019s/result-denoised.png}
    %\caption{Intensity of the sun disk based on the angle of incidence}
    %\label{fig:GS2004-graph-intensity}
\end{subfigure}
\label{fig:GS2004-environment-map}
\caption{Environment map of the sun directly above and the way it is sampled using the normal of the wave equation~\cite{GS2004}.}
\end{figure}


\cite{Gru2019s}

%\subsection{Kim. Caustics Using Screen-Space Photon Mapping}

Technique based on photon mapping~\cite{Jen2001}, which is normally used as a~global illumination technique for rendering caustics and indirect lighting, is here specialized to render only caustics 
- technique to render caustics with photon mapping in screen space with hardware ray tracing and a~screen-space denoiser in real-time.

The use of traditional photon mapping in real-time games requires many photons to obtain smooth images. There were other approaches such as \emph{image space photon mapping} (ISPM) that store a~photon as a~volume in the world. Because there are much 
Traditional photon mapping

This approach is called \emph{Screens-space photon mapping} (SSPM), and builds on \emph{image space photon mapping} (ISPM)~\cite{McGuireLuebkeISPM} real-time photon mapping technique. The ISPM stores photons as volumes in the world. Because of the low amount of such photons, they must be spread (scattered) is some way. 

This screen-space photon mapping technique is used purely for caustics, that are caused by both refracted photons and reflected photons.
- also works with transparent objects - in RT manner

The algorithm is divided into following stages described further in more detail, here in summary: First, a~photon is emitted and traced from a~light source to an opaque surface, after which it gets stored in screen space. The emission and tracing stages use DXR to perform ray tracing in world space instead of screen space. The resulting render target texture used to store these photons then becomes a~\emph{screen-space photon map}. Due to photons being stored as texels in the screen-space texture, the screen-space photon map contains noisy caustics. In order to obtain a~smooth photon map, \emph{photon gathering} or \emph{photon denoising} is performed. Finally, the photon map is used for lighting. Additionally, the approach assumes a~deferred rendering system: a G-buffer, including a~depth buffer or roughness buffer, for storing and gathering photons in screen space.

\begin{enumerate}
    \item \textbf{Photon emission}: Emitted photons carry color, intensity and direction. Once a~photon stops in the world, it can be stored in a~single pixel containing only \emph{photon flux} (color and intensity). Since more than one photon can be stored in a~single pixel, the direction is not preserved.
    For a~point light, photons are emitted from the position of the light. Then the \emph{emission flux} $\Phi_e$ of a~photon has the intensity of a~light ray,
\begin{align}
    \Phi_e = \frac{l_e}{4\pi} \frac{1}{p_w p_h} ,
    \label{eq:Kim2019-emission-flux-p}
\end{align}
where $p_w$ and $p_h$ are the sizes of photon rays (??) and $l_e$ is the radiance of the light from the light source:
\begin{align}
    l_e = light color \times light intensity
\end{align}
For a~directional light, photons spread to the whole scene without attenuation. To consider only photons emitted to a~certain area, a~bounding box called a~\emph{projection volume} is used. Its is a~bounding box that encompasses objects that generate caustics as shown in figure~\ref{fig:Kim2019-projection-volume} by volume $V$. The box is then projected to the negative direction of the directional light to obtain a~rectangular light area shown as area $A$. Photons are them emitted only to the projection volume, thus concentrating on objects generating caustics. Additionally, this allows a~control between consistent quality or consistent performance, by having constant photon emission area or constant number of rays.

\todo{fig:Kim2019-projection-volume}

Then the emission flux of a~directional light is
\begin{align}
    \Phi_e^D = l_e \frac{a_l}{p_w p_h} ,
    \label{eq:Kim2019-emission-flux-d}
\end{align}
where $a_l$ is the area of light area $A$~with a~certain resolution in world-space units, and $1/(p_w p_h)$ is the area of each photon for the area $a_l$, since each photon corresponds to a~portion of the light area.

    \item \textbf{Photon tracing}:
        Once a~photon is emitted, then it is ray traced in the world until the maximum number of ray tracing steps is reached or a~surface hit by the photon is opaque enough. 
        When an object is hit, its material is evaluated. If the surface is opaque enough, the photon will stop, and the following conditions are evaluated to find out whether it will be discarded or stored in screen space:
        \begin{enumerate}
            \item The photon's intensity is negligibly small, i.e., the intensity was reduced by passing through transparent objects.
            \item The location to be stored is out of the screen.
            \item Photon's scene depth from the camera is beyond the values stored in the depth buffer.
            \item The photon is part of direct lighting depicted by the red rays in Figure~\ref{fig:Kim2019-projection-volume}. Such photons are redundant since direct light is assumed to be processed in a~separate pass.
        \end{enumerate}
    
    \item \textbf{Storing photons}:
    A~photon is compressed and stored into a~single pixel. The following denoise stage scatters the pixel's photon onto its neighbors, i.e., it gets gathered from its neighbors.

    First, the area $a_p$~of a~pixel in world-space units is computed as
    \begin{align}
        a_p = \left( \frac{2d \tan (\theta_x / 2)}{w} \right)
              \left( \frac{2d \tan (\theta_y / 2)}{h} \right) ,
    \end{align}
    where $w$~and $h$ are the width and height of the screen, respectively; $\theta_x$~and $theta_y$~are the $x$~and $y$~angles of the field of view, respectively; and $d$~is the distance from the eye to the pixel in world space.
    
    Then the flux of a~photon is converted to radiance $l_p$ to be stored in a~pixel as
    \begin{align}
        l_p = \frac{\Phi_e}{a_p} = \left(
            4 \cdot 4\pi \tan (\frac{\theta_x}{2}) \tan (\frac{\theta_y}{2})
        \right)^{-1}
        \frac{l_e}{d^2} \left( \frac{w h}{p_w p_h} \right) .
    \end{align}
    Additionally, due to it being a~screen-space technique, the photon can use the eye vector when being stored. That is, the photon can use the BRDFs of the surfaces to have the specular color as well as diffuse color evaluated.

    \item \textbf{Photon gathering.}
To get an accurate and smooth caustics in real-time, the photons are gathered in screen space using the reflection denoiser from NVIDIA GamesWorks Ray Tracing~\cite{nvidia-denoiser}. Photons are treated as reflections with some tricks. The inputs to the denoiser are: camera data, depth buffer, roughness buffer and normal buffer. A~world is constructed from the camera matrix and depth, the neighbor's reflections are gathered based on normal and roughness. Parameters of the denoiser were then set by experimentation. The comparison between before and after denoised rendering is shown in Figure~\ref{fig:Kim2019-results}.

    \item \textbf{Lighting}
    The photon map is applied in a~standard manner as screen-space lighting on deferred rendering system. That is, a~photon is treated as an additional light for the corresponding pixel.
    
\end{enumerate}

\todo{fig:Kim2019-results}
% Cornell box scene rendered with 2k \times 2k photons from point lights



Limitations:
artifacts and limitations
- When the camera is close to the caustic surfaces photons do not gather well. That is, because pixels that are near the screen frustum's near depth have high resolutions for photons.
- Transparent objects are drawn to the depth buffer as opaque objects are. This way photons can be stored and therefore caustics are visible on transparent objects but they cannot be seen beyond transparent objects.
- The denoising process might cause some of the photons to be blurred out. Additional intensity can be added to caustics for various artistic purposes, however making the result ``physically incorrect''.
A~bilateral photon denoiser can be used instead. Which downsamples a~photon map beforehand. A~downsampled photon map is then used for near-depth surfaces.

\cite{Kim2019}

%\subsection{Meenrattanakorn, Lambers. Real-Time Caustics Using Cascaded Image-Space Photon Tracing}


Algorithm:

- goal/aim: rendering caustics induced by single specular surface
- assumes that the receiver surface is diffuse, with certain supported light paths
- prerequisite data:
    Refractive shadow maps of the light view, 
    G-Buffer of the camera view
to identify the first ray hit from the source

\begin{enumerate}
    \item \textbf{Ray generation}:
        - from RSMs texels as ray origins
        - sampled randomly using blue noise, a low-discrepancy sequence of pseudorandom numbers
        - to preserve the radiance intensity of the final caustics regardless of the number of samples, the sampled flux is scaled by the ratio of the RSMs resolution to the total number of samples
    \item \textbf{Removal of irrelevant sampling points}
        from all the RSMs samples only certain photons meaningfully influence caustics
        takes into account:
            roughness of the material,
            current luminance
        Roughness of the material is used to determinate whether the surface is specular or not.
        The cutoff was set to 0.3\,out of 1.0. Surface with the roughness value greater than this cutoff is assumed to be diffuse, thus is unlikely to bounce the light ray to a~specific direction as a~specular surface would.
        Based on the flux, the current luminance is computed. If its value is too low, then the photon is unlikely to contribute to the final caustics. The cutoff value was set to 0.04.

    Then the ray is set up using the information from the RSMs at the corresponding texture coordinate in a~standard manner.
    The origin is a~view-space position which is reconstructed using the depth value and the projection matrix.
    direction:
        - determine whether the light ray is reflected or refracted
        - based on \emph{Fresnel factor}
        - takes in IOR of both source and destination media
        - the factor is interpreted as a~probability that the light ray is refracted
            if a~random noise value exceeds this factor then the ray is refracted, else is reflected
    \item \textbf{Tracing rays}:
        
    \item \textbf{Irradiance estimation}:
\end{enumerate}


\fi


% For compilation piecewise (see projekt.tex), it is necessary to uncomment it
% \end{document}