Real-Time Wave Simulation of Large-Scale Open Sea Based on Self-Adaptive Filtering and Screen Space Level of Detail

Abstract

The real-time simulation technology of large-scale open sea surfaces has been of great importance in fields such as computer graphics, ocean engineering, and national security. However, existing technologies typically have performance requirements or platform limitations, and the two are often impossible to balance. Based on the camera-view-based screen space level of detail strategy and virtual camera pose adaptive filtering strategy proposed in this article, we have developed a fast and comprehensive solution for rendering large-scale open sea surfaces. This solution is designed to work without the need for special hardware extensions, making it easy to deploy across various platforms. Additionally, it enhances the degrees of freedom of virtual camera movement. After conducting performance tests under various camera poses, our filtering strategy was found to be effective. Notably, the time cost of simulation using 60 waves at the height of 6 m above sea level was only 0.184 ms. In addition, we conducted comparative experiments with four state-of-the-art algorithms currently in use, and our solution outperformed the others with the best performance and suboptimal visual effects. These results demonstrate the superiority of our approach in terms of both efficiency and effectiveness.

1. Introduction

With the widespread application of ocean engineering (such as in navigation systems [1] or marine simulators [2]), scientific visualization, film visual effects, video games, and other fields, the importance of real-time rendering of large-scale open sea scenes has significantly increased in recent years. The challenge at hand is to simulate and render large-scale ocean surfaces while maintaining a balance between efficiency and quality on power-limited platforms like web and mobile devices. Existing industrial solutions for addressing ocean surface rendering issues primarily involve techniques such as Fast Fourier Transform (FFT) for wave simulation [3], which relies on bandwidth and general-purpose computing, Gerstner wave stacking technology [4], which is constrained by iteration times, and the method of bump mapping [5], which often lacks rendering details.

Eric Bruneton et al. proposed a GPU (Graphic Processing Unit) implementation of ocean rendering in their work [6], which involves the analytical calculation of the normal of the ocean surface height field. They also introduced smooth filtering based on a hierarchical representation of mixing geometry. However, their method lacks dynamic filtering strategies for camera roaming, and the screen space strategy they employed faces challenges in handling changes in camera pose. In response to these limitations, we have implemented several improvements to address these issues. In our approach, considering the requirements for cross-platform compatibility, ease of deployment, and the need to balance performance and effectiveness in large-scale ocean surface rendering, we have designed a forward rendering solution without using the General-Purpose Graphics Processing Unit (GPGPU). This approach relies solely on the basic programmable graphics pipeline to avoid multiple passes and does not depend on any other specialized hardware extensions or compute shaders for non-graphics-related tasks. Additionally, we have utilized screen space level of detail and a self-adaptive filtering strategy to achieve a balance between performance and rendering quality. Our method offers advantages such as extremely low time cost and relatively high quality, making it suitable for applications in the field of ocean engineering.

Our main contributions are:

• In Section 3.2, we present a screen space level of detail strategy based on camera perspective and a novel adaptive filtering strategy based on virtual camera pose;
• We propose a large-scale open sea surface modeling solution that does not rely on special hardware extensions. Our pipeline is lightweight and handy for deploying. It is described in detail in Section 3.3;
• We have increased the degree of freedom of virtual camera motion under the screen space strategy, as discussed in Section 3.2.

A comparative analysis of our proposed method against four existing solutions in terms of efficiency has been conducted, specifically focusing on the time required for ocean surface height field modeling per frame. This comparison serves to highlight the advancements achieved through the implementation of our approach. Additionally, comprehensive tests have been performed to assess the modeling time cost at various camera heights, aiming to substantiate the efficacy of our level of detail and filtering strategy. The resulting renderings from different camera perspectives further underscore the enhanced degrees of freedom in camera movement facilitated by our methodology. To validate the usability of our system, we administered system usability scale questionnaires, the outcomes of which provide empirical support for the user-friendliness of our system. Detailed descriptions of our experimental procedures and findings can be found in Section 4.

2. Related Work

In this section, we will discuss four main topics: waveform simulation techniques for water rendering, level of detail schemes, modeling schemes, and ocean engineering applications. These areas cover important techniques and practical applications in the field of water rendering and simulation, providing a comprehensive perspective to understand and explore the latest developments in the field of water simulation.

2.1. Waveform Simulation Techniques for Water Rendering

The technique of simulating waveforms plays a crucial role in rendering dynamic water surfaces in the field of computer graphics. Various methods have been employed to achieve this, including linear wave superposition, particle-based water simulation, physics-based approaches, pre-rendering, and statistical models. Linear wave superposition involves the accumulation of different linear wave functions to construct the undulating surface, which is particularly effective for deep water. One commonly used method is the Gerstner wave introduced by Gerstner in 1802 [7]. Particle-based water simulation, initially proposed by Yuksel in 2007 [8] and improved by Jeschke et al. [9,10], uses particles to represent individual water waves, enabling wave reflection and interaction with dynamic objects. This approach simplifies the dynamic simulation of three-dimensional water waves to the level of simulating a particle system moving on a plane. Physics-based methods [11] simulate water by solving the Navier–Stokes equation, suitable for offline rendering. Pre-rendering methods [12,13] involve the prior calculation and storage of water surface wave information, commonly employed in applications requiring high realism. Statistical model methods [14,15], based on FFT, find widespread use in the film and video game industries [16] due to their exceptional realism and global controllability.

2.2. Level of Detail Schemes

LOD (Level of detail), an essential technique in real-time computer graphics, adjusts detail levels at different distances and viewpoints to enhance performance and reduce computational load. Common LOD implementation methods include Geometry Clipmap [17], CDLOD (Continuous Distance-Dependent Level of Detail) [18], and Quadtree [6,19]. Geometry Clipmap utilizes nested regular grids to maintain a seamless and continuous terrain surface, allowing dynamic clipping and the addition of geometric details for improved real-time performance. CDLOD is a real-time continuous detail-level rendering method, avoiding discontinuities between different LODs. Quadtree, a common LOD implementation, dynamically adjusts detail levels at different distances through a hierarchical data structure, providing flexibility and simplicity in LOD schemes.

2.3. Modeling Schemes

Modeling schemes encompass manual modeling, screen space mesh projection, and multi-resolution patch systems. Manual modeling involves adjusting water shape, ripples, and textures using professional tools, and is suitable for static scenes and cases requiring artistic control. Screen space mesh projection [20] employs screen space techniques for real-time water simulation, which is particularly suitable for real-time applications but has limited flexibility. Multi-resolution patch systems [6,21] represent the water surface using patches of varying resolutions, optimizing performance by using low-resolution patches in the distance and high-resolution patches up close to ensure finer details and realism.

2.4. Ocean Engineering Applications

As information technology advances rapidly, water rendering technology finds significant applications in engineering across various domains. In educational platforms, the technology enhances the learning experience for engineering students, providing realistic scenarios for water resource management and hydraulics experiments [22]. In Geographic Information Systems (GIS) [23], water rendering supports map-making and spatial analysis by providing clear visualizations of geographical information such as water bodies and flows. The applications in the field of marine and ocean engineering are focused on the physical modeling of the waves, while rendering technology is essential for visualizing simulation results. In ship engineering [24], water rendering aids in visualizing ship motion and interaction effects in different water bodies, facilitating a deeper understanding of ship performance for design optimization. In fluid dynamics [25], water rendering technology is applied to visualize fluid flow data, aiding engineers in analyzing water flow patterns and optimizing hydraulic structures. In deep-sea petroleum extraction engineering [26], water rendering technology is crucial for simulating water environments, providing real-time ocean data visualization for platform design, operation, and safety analysis. These applications underscore the diverse value of water rendering technology in engineering, offering powerful tools and methods to drive development and innovation in the field.

3. Methodology

The schematic diagram in Figure 1 illustrates the fundamental workflow of our methodology. The diagram is divided into two main sections. The left side delineates the tasks executed on the CPU host, encompassing the generation of wave data utilizing a spectrum formula [27], computation of transformation matrices for subsequent GPU processing, incorporation of camera specifics, and initialization of essential global variables crucial for modeling and shading, such as lighting parameters and temporal settings. On the right side, we elucidate the modeling pipeline executed on the GPU device, which encompasses four primary components. Firstly, in Section 3.1, we present the Gerstner wave model derived from the simplified Euler equation, serving as the theoretical foundation for our wave simulation approach. Secondly, we delve into the self-adaptive screen space mesh generation detailed in Section 3.2, functioning within the vertex shader stage and tessellation control shader stage. Furthermore, Section 3.3 represents the crux of our methodology, introducing an innovative screen space adaptive filtering strategy to render animated waves at varying levels of detail. Lastly, we offer a concise overview in Section 3.4 outlining our shading strategy, which, while not the primary focus of our research, plays a pivotal role in ensuring the seamless operation of the entire system.

3.1. Wave Model

The wave model begins with the Euler Equations of an incompressible fluid:

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{v}}{\partial t}+\left(\mathbf{v}·\nabla \right)\mathbf{v}& =-\nabla p+\mathbf{F}\hfill \\ \hfill \nabla ·\mathbf{v}& =0\hfill \end{array}$$

(1)

where $\mathbf{v}(\mathbf{x},t)$ represents the velocity field, $\mathbf{x}$ represents position, t is time, p represents hydrostatic pressure, and $\mathbf{F}$ represents external force. $\mathbf{F}$ can be represented as a gradient of some scalar field if it is conservative. In addition, for waves that have movements caused by a conservative force, we have ${\oint }_L\mathbf{v}·\mathrm{d}\mathbf{l}=0$, and the velocity field $\mathbf{v}$ can be also represented as the gradient of some scalar function $\phi (\mathbf{x},t)$ where $\mathbf{x}=(x,y,z)$ represent field position. Furthermore, after linearization, if we assume $\mathbf{F}$ is just gravity, and consider only the changing shape of the water surface, this simplification allows us to neglect the complex ocean appearance such as surf breaking with bubbles; underwater motion can also be ignored. This focused approach enables us to derive the following simplified model:

$$\begin{array}{c}\hfill \frac{\partial \phi }{\partial t}=-p+p_a-gy\\ \hfill {\nabla }^2\phi =0\end{array}$$

(2)

Since we constrain wave movements on the ocean surface, we have $p=p_a$, where $p_a$ is atmospheric pressure. After differentiating the first equation by t, we obtain $\frac{{\partial }^2\phi }{\partial t^2}=-g\frac{\partial h}{\partial t}$. and by setting the kinematic term of surface $\frac{\partial y}{\partial t}={\mathrm{v}}_{\mathrm{y}}$, we obtain:

$$\begin{array}{cc}\hfill \frac{{\partial }^2\phi }{\partial t^2}+g\frac{\partial \phi }{\partial y}& =0\hfill \\ \hfill {\nabla }^2\phi & =0\hfill \end{array}$$

(3)

By using bottom boundary conditions $\frac{\partial \phi }{\partial y}$ in $y=-\eta$ and a free-surface boundary condition $\frac{{\partial }^2\phi }{\partial t^2}+g\frac{\partial \phi }{\partial y}=0$ in $y=0$ to solve ${\nabla }^2\phi =0$ [28], we obtain:

$$\begin{array}{c}\hfill \phi =\frac{gA}{\omega }sin(\mathbf{k}·\mathbf{x}-\omega t)\frac{cosh(\parallel \mathbf{k}\parallel (y+\eta \left)\right)}{cosh(\parallel \mathbf{k}\parallel \eta )}\end{array}$$

(4)

In Equation (4), A denotes the amplitude, $\omega$ represents the angular frequency, and $\mathbf{k}$ signifies the wave vector. As depth $\eta \to \infty$, the velocity potential would be:

$$\begin{array}{c}\hfill \phi =\frac{gA}{\omega }sin(\mathbf{k}·\mathbf{x}-\omega t)e^{\parallel \mathbf{k}\parallel y}\end{array}$$

(5)

At the ocean surface, where $y=0$, the surface height field can be determined as follows:

$$\begin{array}{cc}\hfill \phi & =\frac{gA}{\omega }sin(\mathbf{k}·\mathbf{x}-\omega t)\hfill \\ \hfill h& =-\frac{1}{g}\frac{\partial \phi }{\partial t}\hfill \end{array}$$

(6)

To achieve the appearance of choppy waves in practical applications, it is advisable to introduce a slight offset along the gradient direction of the height field $h({\mathbf{x}}_{\perp },t)$, where ${\mathbf{x}}_{\perp }=(x,z)$, as suggested by [3]. Consequently, we have:

$$\begin{array}{cc}\hfill \Delta {\mathbf{x}}_{\perp }& =\frac{1}{\parallel \mathbf{k}\parallel }\nabla h({\mathbf{x}}_{\perp },t)=\frac{1}{\parallel \mathbf{k}\parallel }\left[\begin{array}{c}-{\mathrm{k}}_{x}Asin(\mathbf{k}·\mathbf{x}-\omega t)\\ -{\mathrm{k}}_{z}Asin(\mathbf{k}·\mathbf{x}-\omega t)\end{array}\right]\hfill \end{array}$$

(7)

where $\Delta {\mathbf{x}}_{\perp }$ represents the offset of the horizontal projection of the ocean surface position $\mathbf{p}$ concerning the horizontal projection of field position ${\mathbf{x}}_{\perp }$. Combining $\Delta {\mathbf{x}}_{\perp }$ and $h({\mathbf{x}}_{\perp },t)$, we obtain $\Delta \mathbf{x}={(\Delta {\mathrm{x}}_x,h({\mathbf{x}}_{\perp },t),\Delta {\mathrm{x}}_z)}^{\mathrm{T}}$. Therefore, the parametric equation for the ocean surface position $\mathbf{p}$ can be expressed as follows:

$$\begin{array}{cc}\hfill \mathbf{p}& {=\mathbf{x}|}_{y=0}+\Delta \mathbf{x}=\left[\begin{array}{c}x\\ 0\\ z\end{array}\right]+\left[\begin{array}{c}-\frac{{\mathrm{k}}_x}{\parallel \mathbf{k}\parallel }Asin(\mathbf{k}·\mathbf{x}-\omega t)\\ Acos(\mathbf{k}·\mathbf{x}-\omega t)\\ -\frac{{\mathrm{k}}_z}{\parallel \mathbf{k}\parallel }Asin(\mathbf{k}·\mathbf{x}-\omega t)\end{array}\right]\hfill \end{array}$$

(8)

Equation (8) demonstrates that the shape of waves can be effectively modeled using trochoids. Utilizing the Gerstner wave model $\mathbf{p}(\mathbf{x},t)$, it becomes straightforward to determine the height offset of each surface position relative to the horizontal plane $y=0$ in world space. Additionally, this model facilitates the accurate calculation of the normal vector representation.

$$\begin{array}{c}\hfill \mathbf{n}=\mathrm{normalize}(\frac{\partial \mathbf{p}}{\partial x}\wedge \frac{\partial \mathbf{p}}{\partial z})\end{array}$$

(9)

3.2. Self-Adapting Screen Space Meshing

In this section, our objective is to transform the wave model into a mesh representation. The desired outcome is to achieve the appearance of an infinite ocean surface from any viewpoint, provided that $y>0$ and there is no roll component in the camera’s orientation (the treatment of the camera’s roll component is specifically addressed in Section 3.2). Before proceeding, let us define some fundamental symbols: $\mathbf{P},\mathbf{V}$ denote the projection matrix and the camera view matrix, respectively; homogeneous coordinates $\mathbf{v}$, with $z=0$, are used to represent a vertex point in screen space; and $\mathbf{p}$ is used to denote a mesh vertex in world space. Taking perspective division into account, the point at $z=\infty$ along the z-axis in world space can be depicted as follows:

$$\begin{array}{c}\hfill {\mathbf{x}}_{\infty }={\left(\begin{array}{cccc}0& 0& -1& 0\end{array}\right)}^{\mathrm{T}}\end{array}$$

(10)

Assuming $f=\infty$, namely the far clip plane distance of $\mathbf{P}$ is infinity, the projected position in NDC (Normalized Device Coordinate) space can be expressed as follows:

$$\begin{array}{c}\hfill {\mathbf{v}}_{\infty }=\mathbf{P}\mathbf{V}{\mathbf{x}}_{\infty }\end{array}$$

(11)

After perspective division, the y component of ${\mathbf{v}}_{\infty }$, denoted as $H\in [-1,1]$, represents the maximum height of the projected grid in screen space, as illustrated in Figure 2b. To accommodate adaptive resolution, we leverage the hardware tessellation feature to tessellate the area of the screen that is occupied. Initially, we start with two patches comprising eight vertices $\left\{{\mathbf{v}}_i\right\}$, which cover the screen’s NDC (Normalized Device Coordinates) area of ${[-1,1]}^2$. Given that the camera motion lacks a roll component, the projected grid area must form a rectangle. Consequently, we can straightforwardly adjust the y component of the four upper vertices of the two rectangles to H.

$$\begin{array}{c}\hfill {\mathbf{v}}_i={\left(\begin{array}{cccc}{\mathrm{v}}_{i,x}& H& 0& 1\end{array}\right)}^{\mathrm{T}},i\in \{0,1,4,5\}\end{array}$$

(12)

We have prepared two rectangular patches for tessellation. To implement tessellation, we utilize OpenGL’s control shader and evaluation shader. A self-adaptive tessellation strategy is crucial to ensure that we do not generate an overly fine-grained mesh when the projected area is small. This situation often arises when the angle between the view direction and the $y=0$ plane approaches zero. Let us define N (with $supN$ dependent on the hardware device, typically $supN\le 64$) as the subdivision count, which represents the maximum number of cells that can be subdivided along the width of the rectangle. The self-adaptive subdivision count can then be calculated as follows:

$$\begin{array}{c}\hfill ⌊\frac{(H+1)N}{2}⌋\end{array}$$

(13)

The OpenGL control shader plays a pivotal role in adjusting the tessellation level of detail within a quad domain by modifying its inner and outer levels. Following this, the evaluation stage of the pipeline subdivides the two patches into the desired triangle mesh. After tessellation, we obtain the positions of vertices $\left\{{\mathbf{x}}_i^{\prime }\right\}$, each with a z component of 0. The subsequent step involves implementing inverse projection to transform these vertices from NDC (Normalized Device Coordinates) space back to world space. Initially, it is important to note that the projected vertices $\left\{{\mathbf{b}}_i^{\prime }\right\}$ located on the far clip plane have a z value of 1.

$$\begin{array}{c}\hfill {\mathbf{b}}_i^{\prime }={\left(\begin{array}{cccc}{\mathrm{x}}_x^{\prime }& {\mathrm{x}}_y^{\prime }& 1& 1\end{array}\right)}^{\mathrm{T}},\{i\in 0,1,4,5\}\end{array}$$

(14)

The direction vector from the camera position $\mathbf{o}$ to ${\mathbf{x}}_i^{\prime }$ should be:

$$\begin{array}{c}\hfill {\mathbf{d}}_i=\frac{{\mathbf{V}}^{-1}{\mathbf{P}}^{-1}{\mathbf{b}}_i^{\prime }-\mathbf{o}}{\parallel {\mathbf{V}}^{-1}{\mathbf{P}}^{-1}{\mathbf{b}}_i^{\prime }-\mathbf{o}\parallel }\end{array}$$

(15)

The distance from $\mathbf{o}$ to the horizontal plane at $y=0$ along ${\mathbf{d}}_i$ is represented as $t_i$.

$$\begin{array}{c}\hfill t_i=\frac{{\mathrm{o}}_y}{|{\mathrm{d}}_{i,y}|}\end{array}$$

(16)

So the intersection, namely inverse projected vertex position ${\mathbf{x}}_i$ can be represented as:

$$\begin{array}{c}\hfill {\mathbf{x}}_i=\mathbf{o}+t_i{\mathbf{d}}_i\end{array}$$

(17)

Following the application of inverse projection to all tessellated screen space triangles, they can be manipulated just like any conventional mesh.

Roll Component Processing

Euler angles are commonly used to represent rotations for rigid bodies in engineering practice due to their direct nature. Given a local coordinate basis $\{\mathbf{r},\mathbf{t},\mathbf{f}\}$, it is straightforward to decompose a rotation into three components: pitch, roll, and yaw. These components describe the counterclockwise rotation, in radians, along each axis $(\psi ,\theta ,\phi )$. However, it is worth noting that Euler angles do not align with the convention adopted in ocean engineering and the marine sector in general [29]. To address this, we can use Equation (18) to facilitate the mutual conversion between Euler angles $(\psi ,\theta ,\phi )$ and quaternions $(q_w,q_x,q_y,q_z)$. Additionally, in engineering practice, when $\theta$ approaches $\frac{\pi }{2}$, we can adjust $\theta$ by $\pm \delta \theta$ to avoid the issue of gimbal lock.

$$\begin{array}{cc}\hfill & \left[\begin{array}{c}\phi \\ \theta \\ \psi \end{array}\right]=\left[\begin{array}{c}\mathrm{atan}2(2(q_w{q}_x+q_y{q}_z),1-2(q_x^2+q_y^2))\\ arcsin\left(2(q_w{q}_y-q_z{q}_1)\right)\\ \mathrm{atan}2(2(q_w{q}_z+q_x{q}_y),1-2(q_y^2+q_z^2))\end{array}\right]\hfill \\ \hfill & \left[\begin{array}{c}{q}_w\\ q_x\\ q_y\\ q_z\end{array}\right]=\left[\begin{array}{c}cos(\phi /2)cos(\theta /2)cos(\psi /2)+sin(\phi /2)sin(\theta /2)sin(\psi /2)\\ sin(\phi /2)cos(\theta /2)cos(\psi /2)-cos(\phi /2)sin(\theta /2)sin(\psi /2)\\ cos(\phi /2)sin(\theta /2)cos(\psi /2)+sin(\phi /2)cos(\theta /2)sin(\psi /2)\\ cos(\phi /2)cos(\theta /2)sin(\psi /2)-sin(\phi /2)sin(\theta /2)cos(\psi /2)\end{array}\right]\hfill \end{array}$$

(18)

As depicted in Figure 2b, it is evident that the camera movement involving yaw and pitch components does not impact the horizontal orientation of the projection plane. Utilizing Euler angles, any arbitrary rotation matrix $\mathbf{R}$ in ${\mathbb{R}}^3$ can be expressed as:

$$\begin{array}{cc}\hfill \mathbf{R}& ={\mathbf{R}}_{\mathbf{f}}\left(\phi \right){\mathbf{R}}_{\mathbf{t}}\left(\theta \right){\mathbf{R}}_{\mathbf{r}}\left(\psi \right)\hfill \end{array}$$

(19)

Since matrix multiplication adheres to the associative law, we have the flexibility to postpone the roll component ${\mathbf{R}}_{\mathbf{f}}\left(\phi \right)$ until after the projected vertices’ offset, denoted as ${\mathbf{p}}_i$, has been performed. This implies that we should separate ${\mathbf{R}}_{\mathbf{f}}$ from the composite rotation matrix $\mathbf{R}$.

$$\begin{array}{c}\hfill \phi =\mathrm{atan}2(\frac{{\mathrm{R}}_{21}}{cos\theta },\frac{{\mathrm{R}}_{11}}{cos\theta }),{\mathrm{R}}_{31}\ne \pm 1\end{array}$$

(20)

Fortunately, in practice, our raw rotation data for camera movement is often represented by Euler angles or a quaternion rather than a rotation matrix. This means that $\phi$ is readily available to us, and Equation (20) is primarily for situations where Euler angles are not accessible to developers, such as in shader redevelopment. The final issue pertains to the local coordinate basis, which involves ${\mathbf{R}}_{\mathbf{f}}\left(\phi \right)$ describing the rotation by the base vector $\mathbf{f}$, namely the forward vector of the camera. Consequently, the roll matrix should be computed using Rodrigues’ rotation formula:

$$\begin{array}{c}\hfill {\mathbf{R}}_{\mathbf{f}}\left(\phi \right)=cos\phi \mathbf{I}+(1-cos\phi )\mathbf{f}{\mathbf{f}}^{\mathrm{T}}+sin\phi {\mathbf{f}}^{\wedge }\end{array}$$

(21)

where ${\mathbf{f}}^{\wedge }$ represents the skew-symmetric matrix constructed by $\mathbf{f}$. ${\mathbf{R}}_{\mathbf{f}}\left(\phi \right)$ should be calculated in the CPU host and then be transmitted into the GPU device as global data for proficiency consideration. We use Equations (19) and (20) to decompose roll component out of camera matrix $V^{\prime }$ used in Algorithm 1 and the shading part, and we use Equation (21) to construct V used in Section 3.2, Algorithms 1 and 2.

Algorithm 1: Calculate Filter Threshold.

Input: View direction $\mathbf{d}$, view matrix $\mathbf{V}$ without roll component, camera position $\mathbf{o}$, distance n from camera to near clip plane, half fov $\theta$, subdivision count N Output: Vertical length L for current grid cell 1  calculate filter threshold; 2  $\mu \leftarrow ntan\theta /N$; 3  $\mathbf{t}\leftarrow {(0,1,0,0)}^{\mathrm{T}}$; 4  $\mathbf{f}\leftarrow {(0,0,-1,0)}^{\mathrm{T}}$; 5  $\mathrm{d}\mathbf{u}\leftarrow \mu \mathbf{t}$; 6  ${\mathbf{d}}^{\mathrm{view}}\leftarrow \mathbf{V}\mathbf{d}$;             // $\phantom{(}$transform $\mathbf{d}$ into camera space 7  ${\mathbf{d}}_{\perp }\leftarrow (0,{\mathrm{d}}_y^{\mathrm{view}},{\mathrm{d}}_z^{\mathrm{view}},0)$;    // $\phantom{(}$${\mathbf{d}}_{\perp }$ and ${\mathbf{d}}_i$ in Figure 3 are equivalent 8  ${\mathbf{q}}_{\perp }\leftarrow n{\mathbf{d}}_{\perp }/({\mathbf{d}}_{\perp }·\mathbf{f})+\mathrm{d}\mathbf{u}$;   // $\phantom{(}$${\mathbf{q}}_{\perp }$ and ${\mathbf{q}}_i$ in Figure 3 are equivalent 9  ${\mathbf{d}}^{\prime }\leftarrow {\mathbf{V}}^{-1}{\mathbf{d}}_{\perp }/\parallel {\mathbf{V}}^{-1}{\mathbf{d}}_{\perp }\parallel$;      // $\phantom{(}$transform ${\mathbf{d}}_{\perp }$ back to world space$\phantom{(}$ 10  ${\mathbf{q}}^{\prime }\leftarrow {\mathbf{V}}^{-1}{\mathbf{q}}_{\perp }/\parallel {\mathbf{V}}^{-1}{\mathbf{q}}_{\perp }\parallel$;     // $\phantom{(}$transform ${\mathbf{q}}_{\perp }$ back to world space$\phantom{(}$ 11  $\alpha \leftarrow |{\mathrm{o}}_y/{{\mathrm{d}}^{\prime }}_y|$; 12  $\beta \leftarrow |{\mathrm{o}}_y/{{\mathrm{q}}^{\prime }}_y|$; 13  $L\leftarrow \parallel \beta {\mathbf{q}}^{\prime }-\alpha {\mathbf{d}}^{\prime }\parallel$; 14  return L

Figure 3. Algorithm for calculating vertical grid length of each screen space vertex.

Figure 3. Algorithm for calculating vertical grid length of each screen space vertex.

Algorithm 2: Meshing and Inverse Projection.

Algorithm 3: Wave.

3.3. Screen Space Level of Detail

In the previous step, we obtained a flat triangle mesh plane, and the subsequent task is to displace each horizontal vertex by our wave model $\mathbf{p}(x,y,z)$. We must begin by generating wave data. A Gerstner wave is primarily composed of the following components: wavelength $\lambda$, which represents the distance over which the wave completes one full cycle; amplitude A, which denotes the maximum height of the wave; angular frequency $\omega$, which refers to the phase addition per unit time; and wave vector $\mathbf{k}$, which describes the propagation direction of the wave. For optimization purposes, we generate these data [30] using random number generators on the CPU host only once at the beginning and deliver it, along with any other shading-related data, to the GPU device memory through an OpenGL uniform buffer. This eliminates the need for updating to be performed during every render loop.

The Pierson–Moskowitz spectrum, based on measurements of waves taken by accelerometers on British weather ships in the North Atlantic, is utilized [27]. The spectra are characterized by the following form:

$$\begin{array}{c}\hfill S\left(\omega \right)=\frac{\alpha g^2}{{\omega }^5}exp(-\beta {\left(\frac{{\omega }_0}{\omega }\right)}^4)\end{array}$$

(22)

which shows the energy distribution of gravity waves as a function of their frequency, where $\alpha =8.1\times {10}^{-3}$, $\beta =0.74$, ${\omega }_0=g/U_{19.5}$ and $U_{19.5}$ is the wind speed at a height of 19.5 m above the sea surface, the height of the anemometers on the weather ships used by Pierson and Moskowitz in 1964 [27].

Suppose we have N waves with pre-calculated parameters, each animated vertex $\mathbf{x}$ will be offset by ${\sum }_{i=1}^N\Delta {\mathbf{x}}_i$. Consequently, we obtain the animated vertex as follows:

$$\begin{array}{c}\hfill {\mathbf{p}}_i={\mathbf{x}}_i+\sum _{i=1}^N\Delta {\mathbf{x}}_i\end{array}$$

(23)

However, implementing this concept directly on our grid may not be ideal. This is because our screen-space projected grid is regular and uniform, resulting in significant variation in the resolution of the inversely projected horizontal grid for different $t_i$. In particular, the difference in the stride length that a grid represents between the camera’s near clip plane and far clip plane will be quite substantial. This discrepancy leads to a severe aliasing problem, as we attempt to sample waves across all frequencies with an irregular spatial sampling frequency, namely the inverse projected grid resolution. Low-pass filtering becomes essential. Given that we are using a free camera, the filtering strategy should be adaptive. We need to take into account the Nyquist sampling theorem:

$$\begin{array}{c}\hfill f_s\ge 2f_{max}\end{array}$$

(24)

in Equation (24), $f_s$ represents the sample frequency, and $f_{\max }$ denotes the maximum frequency of signal s. If the sampling frequency is too low, it fails to accurately reconstruct the true signal; conversely, if the sampling frequency is excessively high, it results in a large frequency resolution. In practical scenarios, the sampling frequency is typically chosen to be higher than twice the Nyquist frequency (often by a factor of k, where $k\in [3,5]$). In our application, the sampling frequency, equivalent to the resolution of the inverse projected grid, remains fixed. To adhere to the Nyquist limit for each grid, we employ a low-pass filter for adaptation.

Back to physics, the mechanical wave parameter $\lambda$ represents the propagation length of the wave in a full phase, $\omega$ represents the addition of phase per unit of time, and the wave speed v can be represented as:

$$\begin{array}{c}\hfill v=\lambda \omega \end{array}$$

(25)

v indicates how far a wave propagates per unit of time, and when we fix time, we have a low-pass filter:

$$\begin{array}{c}\hfill v\le \frac{L}{k}\end{array}$$

(26)

The illustration of our adaptive algorithm is shown as Figure 3. Assuming we have a 2D camera atlas $\{\mathbf{o},\{\mathbf{t},\mathbf{f}\}\}$, L represent inverse projected grid length projected to $\mathbf{f}$, the viewing direction unit vector. The intersection point ${\mathbf{u}}_i$ between ${\mathbf{d}}_i$ (which is equivalent to ${\mathbf{d}}_{\perp }$ in Algorithm 1) and screen space projected grid can be calculated by solving equation $(\alpha {\mathbf{d}}_i-n\mathbf{f})·\mathbf{f}=0$ for $\alpha$.

$$\begin{array}{c}\hfill {\mathbf{u}}_i=\mathbf{o}+\alpha {\mathbf{d}}_i\end{array}$$

(27)

Since the screen space projected grid is uniform, using $\mu$ to represent the height of each screen space projected grid in world space, the offset of ${\mathbf{u}}_i$ along the top vector $\mathbf{t}$ can be represented as:

$$\begin{array}{c}\hfill \mathrm{d}{\mathbf{u}}_i=\mu \mathbf{t}=\frac{ntan\theta }{N}\mathbf{t}\end{array}$$

(28)

where $\theta =\frac{\mathrm{fov}}{2}$, and we can easily obtain an auxiliary vector ${\mathbf{q}}_i$ (which is equivalent to ${\mathbf{q}}_{\perp }$ in Algorithm 1).

$$\begin{array}{c}\hfill {\mathbf{q}}_i=\frac{{\mathbf{u}}_i+\mathrm{d}{\mathbf{u}}_i-\mathbf{o}}{\parallel {\mathbf{u}}_i+\mathrm{d}{\mathbf{u}}_i-\mathbf{o}\parallel }\end{array}$$

(29)

In Algorithm 1, $\beta$ shares the same principle as t, and the threshold length L can be calculated by:

$$\begin{array}{c}\hfill L=\parallel \beta {\mathbf{q}}_i-t{\mathbf{d}}_i\parallel \end{array}$$

(30)

The overall process of the modeling is shown in Algorithm 2, and the vertex attributes in the output $\mathcal{V}$ of Algorithm 2 will be interpolated by the GPU for shading in the pixel shader.

3.4. Rendering

Since our focus is on modeling, we employ only fundamental shading techniques for rendering. This includes adopting a constant value for sun radiance and utilizing environmental mapping instead of atmospheric scattering models for global illumination.

Although the Algorithm 2 outputs interpolated vertex attributes to the pixel shader, including $\mathbf{p}$, $\partial \mathbf{p}/\partial u$, $\partial \mathbf{p}/\partial v,\mathit{\sigma }$, we can still follow the approach of [6] to execute Algorithm 3 per pixel in the pixel shader, and we also follow the approach of [6] to combine three different ways of dealing with the effects of the waves on scene appearance: displacement maps, normal perturbation, and BRDF (Bi-directional Reflection Distribution Function) modification with the help of the GLSL (OpenGL Shading Language) function smoothstep and the threshold L. Although not as accurate as analytical calculations, in fact using attribute interpolation for coloring is still feasible. If interpolation is used for shading, more patches and vertices are needed to improve the grid resolution of screen space. At the same time, due to the issue of uniform grid resolution of screen space, certain hacks (such as volumetric fog) are needed to alleviate artifacts near the far clip plane.

Given the conversion of the filtering scale from grid to pixel, effectively resulting in an increased sampling rate, it is necessary for the filtering scale, denoted as P, to be reduced:

$$\begin{array}{c}\hfill P=\frac{L}{\mu h}\end{array}$$

(31)

The variable h signifies the height resolution of the screen and $\mu h$ is the sampling frequency. Concurrently, given that the control flow on a GPU significantly depletes register resources, as highlighted in [31], an optimization strategy is adopted from [6,31]. This approach involves pre-sorting the waves on the CPU host and employing a threshold $T=Pk$ as the criterion for filtering levels. Subsequently, a mapping is utilized to directly ascertain the starting index, where k represents the Nyquist sampling frequency.

$$\begin{array}{c}\hfill ⌊{log}_2\left(\frac{T}{inf\left\{{\lambda }_i{\omega }_i\right\}}\right)\frac{C}{{log}_2\left(\frac{sup\left\{{\lambda }_i{\omega }_i\right\}}{inf\left\{{\lambda }_i{\omega }_i\right\}}\right)}⌋\end{array}$$

(32)

where C represents the size of the pre-sorted wave array. The visualization result of mapped indices is shown in Figure 4b. It can be seen that as the z value increases, the filtering interval continues to increase, which means that we can directly calculate the array interval instead of judging if $\left|\lambda \omega \right|\ge kL$ in Algorithm 3.

4. Results

We executed our system, developed in C++ and GLSL 4.5, on an Intel Core i5-12400F processor (Intel, Santa Clara, CA, USA) with 16 GB of memory (six cores) and an NVIDIA GeForce RTX 2060 GPU (NVIDIA, Santa Clara, CA, USA). Within this system, our core algorithms are integrated, and we will present our key findings in the subsequent sections.

To ensure a fair comparison, we exclusively focus on the time required for wave modeling and analytical normal calculation, both essential for shading purposes. Complex shading processes like atmospheric scattering or subsurface scattering are not taken into account in our evaluation.

With a $1280\times 800$ resolution and a camera height of 6 m, two quad patches, and maximum tessellation level $N=64$, using a NVIDIA Nsight Graphics tool to capture the frame, we obtain ≥1500 fps (disable v-sync) with 60 waves with $\lambda \in (0.01,30)$ (1500 fps contains the time cost brought by the context of GLFW window form library and imgui GUI library). As described in Section 3.4, we analytically calculate normals per pixel for shading. Since a large number of high-frequency waves are filtered out by low-pass filters in the modeling process, this gives a $0.184$ ms time cost of the superposition of waves totally per frame, and when the camera height rises to 60 m, the total time cost of modeling and analytic normal calculating has been shrunk to $0.152$ ms, and if we remove our self-adaptive filtering model, the time cost rises to close to $0.412$ ms, no matter what camera height is selected, which proves our LOD strategy is effective. To verify the effectiveness of our self-adaptive filtering method and show its balance between performance and speed, we compare our algorithm with three current solutions: the built-in ocean water system of the Unreal 5 Engine based on Gerstner waves and a tile-based LOD system; NVIDIA WaveWorks 2.0 based on FFT and a Continuous Distance-Dependent LOD system; and a Quadtree tessellation LOD system.

To ensure the utmost fairness, after several tests, we used 100 waves in the Gerstner-wave-based method (ours, Unreal 5, Quadtree tessellation). As for the FFT-based method, known for its remarkable capabilities in simulating large-scale waveforms, WaveWorks 2.0, we adopted the configuration from NVIDIA’s official samples. The comparative outcomes are presented in

Table 1. Average time (ms) comparison among our self-adaptive filtering work, Unreal 5 ocean system, NVIDIA WaveWorks 2.0, and Quadtree tessellation. Resolution is set to $1280\times 800$ and all data come from NVIDIA Nsight Graphics and Unreal 5 built-in performance profiler.

Method	Time Cost (ms)	FPS	FPS (Disable v-Sync)
Ours	0.553	165	≥1328
WaveWorks 2.0	3.872	165	≥563
Unreal 5 ocean system	9.542	89	≥85
Quadtree tessellation	6.147	165	≥328

.

Figure 5 and Table 1 present the visualization outcomes of four distinct methods. Within Figure 5b, it becomes apparent that NVIDIA WaveWorks 2.0 stands out in terms of visual fidelity. Nonetheless, this superiority is offset by significant drawbacks, including elevated computation times on general-purpose hardware and intricate challenges in engineering execution, which collectively hinder its practical deployment. Furthermore, it necessitates additional multi-resolution FFT computations to fulfill anti-aliasing prerequisites, thereby exacerbating the IO load. As a result, FFT-based approaches find their niche primarily in high-end console or PC gaming environments, as well as in the realm of cinematic visual effects. Meanwhile, their application might not extend seamlessly to more resource-constrained platforms such as web applications and mobile games, or to specific use cases within ocean engineering, including navigation systems [1] and marine simulation software [2]. These considerations do not align with our present objectives.

In the Quadtree experiment, we utilized the traditional camera-position-based tessellation strategy for Level of Detail (LOD) purposes and calculated the normals analytically without applying any filtering to the model. The results are depicted in Figure 5c. It is evident that as $z\to \infty$, pixels near the horizontal line display significant aliasing issues. Unlike screen space tessellation, the irregularity and lack of uniformity in Quadtree nodes pose challenges in designing an effective low-pass filter. Additionally, rendering the ocean surface controlled by Quadtree nodes necessitates additional draw calls, resulting in increased time costs.

To corroborate the effectiveness of our method, we performed a series of tests aimed at assessing the rendering time required for a single frame across various camera heights $y\in (6,66)$ within a designated scene. In the interest of ensuring fairness, we contrasted these results with the time required to activate wind waves in WaveWorks 2.0, since the shading algorithm utilized by WaveWorks 2.0 demands a substantially higher computational effort in contrast to our approach. The results, depicted in Figure 6, unequivocally demonstrate a decrease in rendering time cost as the camera height y increases. This observation unequivocally underscores the efficiency of our Level of Detail (LOD) technique.

The results from the Unreal 5 experiment, illustrated in Figure 5d, demonstrate a traditional ocean rendering pipeline that opts for interpolating geometry information within the fragment shader instead of utilizing analytic calculations. The grid near the horizon results in a discontinuous patchwork effect. Furthermore, this method capitalizes on high-quality textures crafted by artists to enhance visual expression. As a result, the specular component of the rendered ocean might not achieve the same level of realism when compared to alternative methodologies.

To substantiate the enhancement in the degrees of freedom for camera movement, we further selected camera viewpoints from varied perspectives for rendering. The outcomes derived are depicted in Figure 7. Upon substituting V with $V^{\prime }$ as described in Algorithm 2, the issue arising from improper handling of roll components, as illustrated in Figure 8 and further confirmed by the depiction in Figure 2b, validates that our approach indeed amplifies the camera movement’s degrees of freedom without detriment to the screen space algorithm.

To evaluate the effectiveness of our method, we designed a questionnaire based on the System Usability Scale (SUS), encompassing the 10 questions presented in

Table 2. System usability scale questionnaire design and the average score.

	Qusetion	Average Score
A	I think that I would like to use this system frequently.	4.3
B	I found the system unnecessarily complex.	1.2
C	I thought the system was easy to use.	4.5
D	I think that I would need the support of a technical person to be able to use this system.	1.6
E	I found the various functions in this system were well integrated.	4.3
F	I thought there was too much inconsistency in this system.	0.6
G	I would imagine that most people would learn to use this system very quickly.	4.1
H	I found the system very cumbersome to use.	0.8
I	I felt very confident using the system.	4
J	I needed to learn a lot of things before I could get going with this system.	1.3
	Final Score	89.25

. The ratings for each questionnaire range from 1 to 5. We invited 10 experts from relevant fields to interact with our system and complete the questionnaire. The scoring procedure for each question involves a transformation of the participant’s responses: for questions A, C, E, and G, the original score is reduced by 1, while for questions B, D, F, H, and J, it is subtracted from 5. These adjusted scores are then aggregated and multiplied by 2.5, converting the initial range of 0–40 to a standardized scale of 0–100. Based on research, a SUS score above 68 would be considered above average and anything below 68 is below average; however, the best way to interpret the scores is to produce a percentile ranking. The outcomes derived from the questionnaire are depicted in Figure 9.

Although there are some disadvantages in appearance due to the use of a relatively simple shading algorithm, the results presented in Table 1, Figure 5, and Table 2 do demonstrate that our approach outperforms the existing methods in terms of modeling efficiency. This superior performance can be attributed to the inherent benefits of utilizing screen space level of detail, where the granularity of primitives is determined by the screen resolution rather than the complexity of the scene itself. Concurrently, our novel self-adaptive filtering strategy adeptly addresses the challenge of non-uniform sampling across the projected grid. Moreover, in comparison to the FFT-based method, our approach of employing analytic calculations enables us to achieve competitive shading outcomes in simulations with a significantly reduced wave count.

5. Conclusions

This article introduces an enhanced screen space level of detail strategy, which is innovatively based on the camera view, coupled with an adaptive filtering strategy based on virtual camera poses. The combination of these strategies results in a detailed open sea surface rendering solution that is easy to deploy without the necessity for specialized hardware extensions. By enhancing the motion degrees of freedom of virtual cameras while also improving performance and appearance, this approach represents a significant advancement in the field of ocean engineering. For example, in aircraft [32,33] or marine simulators [2], our method can effectively overcome the lack of support for camera rolling in conventional screen space strategies while simultaneously achieving the efficiency of screen space strategy, which provides a new solution for naval applications. Drawing inspiration from NVIDIA’s Variable Rate Shading (VRS) technique [34]—which intelligently applies varying shading rates across different regions—we envisage integrating our method with VRS in future works. We aspire to leverage our method within the virtual reality domain, aiming to propel the frontier of technological advancements further.