Real‐Time Underwater Spectral Rendering

The light field in an underwater environment is characterized by complex multiple scattering interactions and wavelength‐dependent attenuation, requiring significant computational resources for the simulation of underwater scenes. We present a novel approach that makes it possible to simulate multi‐spectral underwater scenes, in a physically‐based manner, in real time. Our key observation is the following: In the vertical direction, the steady decay in irradiance as a function of depth is characterized by the diffuse downwelling attenuation coefficient, which oceanographers routinely measure for different types of waters. We rely on a database of such real‐world measurements to obtain an analytical approximation to the Radiative Transfer Equation, allowing for real‐time spectral rendering with results comparable to Monte Carlo ground‐truth references, in a fraction of the time. We show results simulating underwater appearance for the different optical water types, including volumetric shadows and dynamic, spatially varying lighting near the water surface.


Introduction
The ability to simulate the appearance of underwater scenes is crucial for a wide number of applications, in particular as an analysisby-synthesis approach in domains such as underwater imaging and remote sensing.For such applications, beyond a plausible and appealing appearance, the simulation needs to provide physical correctness, considering both the optical properties of water as well as the spectral response curve of the sensors.Furthermore, water absorption and scattering have very characteristic spectral profiles with strong wavelength dependencies, so spectral rendering is a must for accurate color reproduction of underwater scenes.Spectral light transport simulation of participating media is, on its own, a daunting and time-consuming task, particularly because interactions between light and matter can happen at every differential point of its trajectory, and because of its recursive nature, where every differential scattering interaction generates new light trajectories.The rendering time of offline approaches, on the order of hours, is not practical for many applications for which interactivity is a requirement.
In this paper, we propose a spectral, physically-based, real-time rendering algorithm to simulate underwater light transport at practical frame rates (more than 60 fps in an NVIDIA GeForce GTX 1660 Super).Our method aims to bridge the gap between common real-time solutions in game engines, that are neither spectral nor physically-based enough, and offline simulation techniques, that are prohibitively expensive for interactive applications.Such a speedup in simulations has benefits for a number of applications, including underwater computer vision, visual ecology applications, submarine or remotely operated vehicle training tasks using AR/VR (Augmented/Virtual Reality) setups, or the entertainment industry.
To achieve real-time frame rates, we separate single scattering from multiple scattering, and provide physically-based approximations for each while preserving physical accuracy as much as possible.For single scattering, we rely on froxel lighting [Vos14,Hil15,Kov20] for spatially-varying illumination while accounting for the physically-based behavior of light transport (Section 4.1).For multiple scattering we rely on the apparent optical property of diffuse downwelling attenuation to characterize underwater downwelling irradiance, as usual in ocean optics (Section 3.2).Through some derivations on the Radiative Transfer Equation [Cha50], this enables us to devise an expression that can be evaluated in constant time per wavelength (Section 4.2).
Our results work at interactive frame rates, for a wide number of real-world water types.This simulation of different water types is made possible thanks to their routine characterization in oceanography, and the availability of the real world measurements of the coefficients [SM15, WBF * 03].We analyze the importance of spectral rendering on underwater color reproduction.Our results also provide a rather accurate match compared to a full light transport simulation using path tracing in a minimal fraction of the time.
We also show how our method can be easily integrated with a number of real-time techniques in physically based pipelines, which allows to simulate a variety of underwater phenomena such as moving caustics, the dielectric water surface or the Snell window, as we can see in Figure 1.

Related work
Commonly, full light transport simulations of underwater appearance numerically approximate the physics of the real world through Monte Carlo methods, at the cost of extremely long rendering times.These Monte Carlo methods are often improved by bidirectional techniques such as photon mapping [Jen96], which accelerate the convergence of the projected caustic patterns and their volumetric beams and shadows.Frisvad and colleagues [FCJ07] generalized the Lorenz-Mie theory to rendering natural waters, while Gutierrez et al. [GSMA08] explored the effects of inelastic scattering (common in oceans due to phytoplankton) in the appearance of water.We refer the reader to a complete review of the most advanced offline techniques recently presented by Droske  Real-time techniques.Few works tackle the particularities of physically based rendering of underwater scenes in real time.In game engines, artistic control is often more important than physical accuracy.For example, in Subnautica (2018) the rules of physics are broken to allow for more compelling colors and improved visibility [Dev16].The most common approach in real-time engines is just a tinted fog approach, where radiance is attenuated exponentially.Ambient light reaching both the surface and the medium is usually made depth-dependant, but without following physical interactions [LKA * 17].These approaches are controllable and computationally inexpensive, but do not allow to render different water types directly based on their coefficients.Moreover, real-time approaches usually adopt fog rendering techniques, which are only suitable for the atmosphere since they assume monochromatic extinction.In water bodies, however, extinction has a strong wavelength dependency [ATS * 17, AT18].
Chen et al. [CZSZ22] propose a controllable model for real-time underwater lighting, by making use of a depth-dependant irradiance approximation used in oceanography for ambient illumination.However, they approximate lighting using wideband attenuation coefficients in XYZ space for efficiency, which can lead to large errors in the final images.Furthermore, their approach models ambient illumination as a constant for the whole scene (not accounting for differences in vertical depth), and does not take into account the angular dependence of the veiling light (the underwater saturated color, visible when the viewing path is long enough).
Our proposed solution addresses these limitations, increasing accuracy by computing radiance for more wavelength bands, including both single scattering and caustics, and also allowing to simulate the different water types and depths (see Section 3).Furthermore, our method is flexible, capable of reproducing the captured color by different sensors, making it useful for oceanographic studies that employ consumer cameras.
Underwater light diffusion.Diffusion has been a popular approximation for complex multiple scattering problems in fields such as neutron transport [Cas60] and subsurface scattering rendering [JMLH01].These methods usually rely on reducing these interactions to a one-dimensional problem.In this context, Premoze and Ashikhmin [PA01] developed a model of the diffusion of light in the water volume based on the ocean optics literature and measurements.Although their method was offline and focused on an outside viewer, we propose a similar approach for multiple scattering that allows to simulate the underwater perspective, accounting for the seafloor reflectance, and model different sensors, in real time.We validate our model with respect to full light transport simulations.

Radiative transfer
The Radiative Transfer Equation (RTE) [Lom87,Cha50] models the scalar radiance transfer (i.e.without polarization or inelastic scattering) in a differential manner.In computer graphics, we usually formulate the RTE in its integral form, also known as the volume rendering equation (VRE) [FWKH17].This expression models the final radiance L(o, ω) reaching the camera (at point o looking towards ω, see Table 1 for an overview of the symbols), traveling a distance P: where T (P) = e −σt P is the transmittance of the water for distance P. The extinction coefficient, σt = σs + σa, accounts for outscattering and absorption events.Ls is the radiance at o + Pω, and L i represents the in-scattering: all incident irradiance at any point in the medium scattered towards the sensor.Finally, the integral sums the contribution of every point o + vω along the path.Lm corresponds to the integral below, which accounts for all the medium in-scattering along each differential point in the full path.
The fundamental difficulty of this problem is the computational cost of approximating this integral, as the in-scattering itself at a given point p is yet another nested integral over all directions Ω: where fs is the phase function, which represents the ratio of light being scattered from ω i towards ω, and L(p, ω i ) the incoming radiance from ω i .We aim to find an approximation to Equation (1) that does not require the numerical evaluation of it and Equation (2) while respecting the physical behaviour of radiance fields entering the water bodies.

Ocean optics background
Inherent and apparent optical properties.In the ocean optics literature, the optical properties of water, which determine how it interacts with light, are often classified into Inherent (IOPs) and Apparent (AOPs) Optical Properties.IOPs, such as the absorption and scattering coefficients and the phase function, only depend on the water constituents, and fully describe its light transport properties through Equation (1).AOPs, on the other hand, also depend on the radiance field.For example, the irradiance profile with respect to depth is an AOP, as it depends on the incident irradiance at the surface.However, AOPs are studied because they are easier to measure than IOPs, and many works have been devoted to find relationships between the two.
The so called K functions [Mob94] are the key AOPs in ocean optics, modeling the rate of change of irradiance with respect to depth.Although they do depend on sun angle and other external factors, they are considered quasi-inherent optical properties due to their robustness to external factors [Mob94, LDC * 05, Mob01, Gor89].
Moreover, K functions are widely used in oceanography due to practical aspects, since they are easy to measure in different waters by measuring irradiance at different depths.In many cases, they also show a relationship with biological constituents, and are useful for remote sensing [SB78b,SB78a].
In this work we rely on the diffuse downwelling attenuation Kfunction K d [Mob22], which is defined for the downward irradiance differential as where y is the vertical distance from the surface and E(y) the irradiance at that point.
Asymptotic radiance distribution.All homogeneous, optically deep waters show a characteristic depth-dependent asymptotic radiance distribution.This was first conjectured by Whitney [Whi41b,Whi41a], then given empirical support from actual measurements [Tyl58, Tyl60, JF60, SWOO58]), and finally given mathematical proof through Preisendorfer's radiative transfer model [Pre59].This asymptotic radiance distribution implies that, in homogeneous waters, AOPs such as the K functions become IOPs given enough depth, independent from the boundary conditions at the surface and bottom [Mob94].While the Jerlov water types offer a compact generalization of most commonly encountered optical conditions, they do not capture the large variety seen across the world's natural bodies of water.Thousands of modern-day measurements of K d are available in the SeaWiFS Bio-optical Archive and Storage System (SeaBASS) [WBF * 03].Here, we demonstrate our results using the Jerlov water types for simplicity and also because the corresponding IOPs for them are known, but our method can be used with any real-world K d dataset.

Our method
Near the surface of the ocean, single scattering (Figure 3, green) produces sharp, spatially varying caustic patterns, which are important cues for underwater appearance [DHV * 23].On the other hand, lighting due to multiple scattering in the ocean (see Figure 3, purple) is mainly homogeneous horizontally, and slowly varying vertically.This is related to the asymptotic behaviour of the K functions we mentioned in the previous section [Mob94,Mob22].
From the path radiance Lm(o, ω) in Equation (1) we can obtain the following expression, separating the full contributions of single (L SS (o, ω)) and multiple (L MS (o, ω)) scattering: where L iSS (o + vω, ω) and L iMS (o + vω, ω) account for the single and multiple scattering at each point, respectively.In the follow-ing, we detail how we approximate each of these terms.We omit the spectral dependencies of some terms for simplicity, and explain how they are taken into account at the end of the section.

Single scattering
Single-scattering in-scattering L iSS (y, ω) results from direct sunlight, refracted by the water surface, and attenuated exponentially by the transmission of the medium e −ycos(θsun)σt , where θsun is the vertical angle from the refracted sun direction.We emulate this moving light field by the usual approach [DHV * 23] of procedurally shadowing a directional light source, due to the great efficiency of this approach.
Then, to approximate the resulting spatially varying illumination, we rely on the froxel ligthing technique, widely used in game engines [Vos14, Hil15, Kov20], and implemented in Unity's High Definition Rendering Pipeline (HDRP) [Lag18].Froxel lighting is a two-pass algorithm which first computes a discrete representation of the in-scattering L iSS in the scene, and then traverses this representation to approximate the final single-scattering integral L SS .A summary of this discretization is illustrated in Figure 4, where the view frustum is subdivided in view-space voxels (i.e., froxels).
The second pass approximates the single-scattering integral of each ray by adding the contributions of the corresponding froxels: where o + vω ∈ XV (ω), XV (ω) is the set of froxels for direction ω, and d the length of each froxel along the camera z axis.The L iSS inscattering term is able to consider any phase function, as both the viewing and the light direction direction are known when evaluating the expression.We use a forward-scattering Henyey-Greenstein phase function (with anisotropy g = 0.8).
To handle the air-water refraction of the sunlight, we use two different lights to represent the sun.The above-surface one has the actual direction and affects the sky rendering, while the second directional light, which affects Equation (5) and the underwater objects, is oriented according to Snell's law.

Multiple scattering
In general, multiple scattering is prohibitively expensive to simulate.We leverage ocean optics knowledge (and measurements), introduced in Section 3.2, to devise a faster-to-evaluate analytical approximation.We follow a similar approach to Preisendorfer [Pre76,Mob94] and make the approximation that multiplescattered radiance at any underwater position depends only on vertical depth according to the diffuse downwelling attenuation coefficient K d as shown in Equation (8).This approximation neglects shadows, that are less prominent in multiple scattering than in single scattering, where we actually account for them (Section 4.1).
We assume that the medium is spatially homogeneous.Therefore, the full downwelling irradiance at any depth y results only Figure 4: The froxel lighting technique approximates volumetric single scattering in two passes.First: the 3D view frustum is discretized into voxels (blue grid), which store the incident illumination at each point, accounting for occlusions (grey objects in the image).Then the RTE integral of each ray is approximated by adding the in-scattering of each of its voxels.For Lss(ω 1 ), in orange, it will thus account for the inscattering at xa, while Lss(ω 2 ), in blue, will be zero due to occlusions.from the irradiance at the surface E D (0).We can compute it by integrating the diffuse downwelling attenuation coefficient: We additionally propose to also account for the radiance being reflected from the seafloor, as illustrated in Figure 5, using a similar approach.We assume a Lambertian horizontal floor at depth y f , with reflectance coefficient r.Then, the upwelling irradiance E U (y) being reflected from it can be obtained through Equation ( 6) as By adding both, upwelling and downwelling, we get an irradiance profile solely based on the vertical axis: The water surface radiance L W S accounts for the dielectric nature of the water surface.Snell's window is the circle in the water surface, above the viewer, outside of which there is total internal reflection.Inside it, the whole sky is refracted.Our method naturally simulates this phenomenon, as we can see above the shark.
By assuming an isotropic phase function and substituting the diffuse in-scattering from Equation (4) for this expression we analytically solve the integral over the whole path L MS , as the expression becomes a first-order ODE: where yv = yo + yωv is the vertical depth v units away of the origin in the direction of the ray.After substituting and integrating, this gets simplified into the final expression.
Note that this approximation has turned the complex computation of multiple scattering into a single O(1) expression (per wavelength, per pixel), accounting for both downwelling irradiance from the sun, as well as the upwelling reflectance from the seafloor.
On top of that expression, we also consider the multiplescattering radiance coming from the surface Ls(o, ω) as: which accounts for the diffuse downwelling attenuation until the surface and transmittance attenuation from the surface to the sensor.Water surface.We also consider the special case of the radiance from the water surface L W S (p WS , ω), at a given point p W S , often visible when looking up.We can compute the dielectric color by combining the refracted radiance from the sky L sky and the reflected one (again using L MS Equation ( 9)), according to the normal n and the well-known Fresnel equations: where F(ω, n, η) is the ratio of refraction with respect to reflection for the given vectors and the index of refraction of water (η = 1.33).L sky (ω f r ), the radiance from the sky, can be a texture or a physically based approximation.ω f r is the refracted vector, computed according to Snell's law, and ω f l the specular reflection.
This term allows to visualize Snell's window, the circle in the water surface directly above the viewer, approximately 96º in diameter, in which the whole sky is refracted underwater [Mob94], as we can see in Figure 6.Outside this window, delimited by the critical angle, there is total internal reflection.

Spectral rendering
Until now, we have omitted the wavelength dependency of all coefficients and radiance terms for simplicity.In practice, Equation (9) and Equation (10) are computed per wavelength band, using the corresponding spectral values of the coefficients from the given water type, leading to wavelength-dependent radiance L(λ).
In addition, different sensors have different spectral sensitivities, depending on their spectral response function f (λ).This function returns the RGB triplet that the sensor outputs when excited by wavelength λ.At the top of Figure 7, we can see two different response functions, for the Canon A2200 and the Arriflex D21, re-spectively.As we can see in the bottom row, as the green channel of the first camera shows a stronger response, the same scene looks greener.In our work, we simulate a selection of eighty nine response functions from the most widely used sensors, collected from the literature [GN03, JLGS13, SA23].
The final RGB triplet L RGB corresponding to our L(λ) captured by response f (λ) is obtained as the inner product of both functions: where w is the width of the bands: w = (max(λ) − min(λ))/n, with n the number of wavelength bands evaluated.In practice, we have found n = 8 bands to be a good trade-off between accuracy and computational cost, with higher values bringing progressively diminishing returns in color accuracy.In our implementation, Equation ( 12) and the multiple scattering equations are evaluated per pixel in a single deferred shader.

Results
Our method can be integrated into different rendering pipelines.For this work, we have implemented it as a deferred shader in Unity's High Definition Render Pipeline (HDRP) [Lag18], adding support for spectral data (coefficients and sensor response curves), and combining it with the volumetric buffer computed by HDRP's implementation of froxel lighting.In the following results, we used 256-deep slices for the froxels and set the sun spectrum E D (0) to D65.
Our results show high quality renderings at interactive frame rates even with relatively low-end GPUs: we achieve over 30 frames per second on an NVIDIA GeForce GTX 1050 GPU, and more than 60 frames per second on a more modern NVIDIA GTX 1660 Super.Our physically-based pipeline is able to simulate the appearance of varied Jerlov water types at any depth, captured by a wide variety of sensors.Beyond these results, our method would work with any other water dataset, and any other sensor.
Our video in the supplemental material shows a real-time capture of our interactive application to render different types of ocean waters with different characteristics (see also Figure 9 for a screenshot).First, we start ten meters deep in a Jerlov 1C water (coastal, moderately turbid), and simulate different sensors such as the Nikon D1x or the D220.We then lower depth to 5.6 meters, and change to a Jerlov 3C medium (still coastal, moderately turbid), which presents stronger scattering.We then switch to a Jerlov IB water, and show how in this case light reaches depths of up to 50 meters.Looking up, Snell's window is visible (as in Figure 6).Last, we modify the water's IOPs, showing how our method can be used to render waters beyond Jerlov types: we first double the scattering coefficient, which increases turbidity, then we increase the absorption coefficient by a factor of three.beams from the water surface.Below, at greater depths, the direct sunlight beams get rapidly attenuated, and the diffuse multiple scattering term becomes dominant.
In Figure 10, we show the equivalent table for different coastal waters, where we can see a similar effect.These coastal waters, with higher scattering due to greater concentrations of suspended particles, rapidly extinguish all light at around ten meters below the surface.

Comparison with state of the art
The recent USSim method [CZSZ22] is the closest in the literature both in terms of goal and approach: efficient rendering of controllable, physically accurate underwater imagery.The authors also focus on rendering different Jerlov water types based on their IOPs.However, USSim performs its lighting computations in the XYZ space to reduce computation time, and then converts the resulting colors to sRGB.This introduces large errors in the apperance of the final scenes, as shown in Figure 11 for six different Jerlov waters and depths.Moreover, their depth value is constant for each scene.In contrast, our solution is fully spectral and accounts for varying depths within the scene, leading to much more accurate results with respect to the ground-truth, path-traced images in all six cases.For this comparison, in the path tracer and our spectral method we have applied the CIE XYZ spectral sensitivity curves, then converted to sRGB to bring all methods to the same color space [RKAJ08].
In Table 2, the average per-pixel Root Mean Square Error (RMSE) of each of the scenes is listed for both methods, compared to the references.On average, we achieve half the RMSE.

Computational cost
The theoretical per-frame rendering time t f (n) of our algorithm for n wavelengths is t f (n) = to + t SS + nt MS (1), where to is the overhead time mainly due to other pipeline passes (GBuffer, postprocessing, antialiasing, etc.), t SS is the cost of single scattering, and t MS the cost to evaluate a single wavelength for multiple scattering.Since only the multiple scattering term depends on the number of wavelengths, the order of the algorithm is O(n) with respect to this parameter.
Figure 12 shows the full frame rendering time (in milliseconds) for the reef scene, varying n between 3 and 250, and rendered in Full HD with our full single and multiple scattering pipeline, measured in an NVIDIA GeForce GTX 1660 Super GPU.In blue, we show the aggregated cost of the standard rendering passes (to), which do not depend on n and add up to 11.06 ms on average.This time is mainly due to the GBuffer pass (50 %), the cascaded shadow maps (10 %), temporal antialiasing and other post-processing effects (such as bloom).In orange, t SS is also independent of n and averages 3.57 ms.In green, we show the cost of the deferred shader evaluating multiple scattering for different wavelength bands.For low values of n, its cost is still lower than single scattering (e.g., at t MS (32) = 2.5 < 3.57 = t SS ms).These low values are the interesting ones in practice, as the spectral resolution of our IOP data is limited to 17 bands.Asymptotically, t MS (n) ≈ 0.056n + 1.17 ms.

Qualitative validation
We have compared our multiple scattering approximation using n = 8 wavelength bands with a path-traced Monte Carlo reference, implemented in Mitsuba2 [NDVZJ19], which not only does not rely on the downwelling approximation of irradiance, but also considers a forward-scattering Henyey-Greenstein phase function (with mean cosine g = 0.8), which is more accurate than our isotropic phase function.We also compare our approximation with an RGB version of our method, evaluating Equations ( 9) and (10) for the three most representative wavelength bands given a specific sensor response curve.Figure 12: Full rendering time (ms) with respect to the number of wavelengths evaluated, using our single (orange, t SS ≈ 3.57 ms) and multiple scattering (green, t MS ≈ 0.056n + 1.17 ms) pipeline.The blue time (to ≈ 11.07 ms) corresponds to standard passes, mainly GBuffer, shadow mapping, antialiasing and postprocessing.Using n = 8 as in the rest of the paper, we measured the total time as t f (8) = t MS (8) + t SS + to = 1.8 + 3.57 + 11.07 = 16.44 ms.
Figure 13 shows two examples of this experiment for two scenarios: Jerlov IB water captured by the Canon 600D sensor (left), and Jerlov II water captured by the Canon 400D.We show the corresponding renders of these scenes with three methods: a groundtruth volumetric spectral path tracer (top), our real-time spectral method (second row), and a real-time RGB version (bottom).As we can see, the spectral version of our method achieves a very similar appearance to the path-traced references in both cases.Please refer to the supplemental material for more comparisons.
Note that these Monte Carlo spectral renders used as reference take approximately four hours to converge, while our spectral method renders each frame in around 16 ms.Even if the hardware running these methods is not the same (Mitsuba is a multithreaded CPU program, while the real-time solutions run in the NVIDIA GeForce 1660 Super GPU), our method is around 6 orders of magnitude faster.For further analysis of the computational cost, see Section 5.2.

Quantitative evaluation
We further compare our multiple scattering spectral solution with the RGB implementation using the path tracer simulation as ground truth.We study six combinations of medium coefficients and underwater depths (3, 5 and 9 meters).We include two scatteringdominated waters (Jerlov 3C and Jerlov II), and three absorptiondominated waters (I, IA, IB).Given a scene defined by its medium and depth, we render it using the spectral response curves of 28 different cameras.
We take the pixel-wise difference with respect to the ground truth, using two different metrics: RMSE and the more perceptual Delta E 2000 over CIELAB [Bra03].The results are summarized in Figure 14 where, for instance, Jerlov3C_d9 means a Jerlov 3C medium at 9 meters depth.For both metrics, our method achieves a significantly higher accuracy.for RGB (blue) and our spectral rendering method (orange), compared to ground-truth path tracing.For both metrics the error is significantly lower with our method.

Diffuse downwelling irradiance
The previous results suggest that our choice of an isotropic phase function for the multiple scattering contribution does not introduce as much error as it would in a Monte Carlo simulation, as our references use a more accurate forward-scattering phase function.This can be explained by the fact that the irradiance profile obtained by the downwelling analytical approximation (e −K d y ) closely matches the actual irradiance.
In Figure 15, we show a simplified scenario simulating the (monochromatic) irradiance profiles with depth in an infinite, homogeneous water volume.The lines are the Monte Carlo results for different Jerlov water types (each in a color), with their corresponding IOPs.Crosses are the analytical approximation, E D (y) = e −K d y .Similar experiments can be found in the literature [SM15].This match also means that K d can be recovered by a least-squares op-  timization for other media where it is not routinely measured, as long as their absorption and scattering is known.

Conclusions
We have presented a method that achieves real-time performance for rendering a wide variety of spectral underwater scenes, where it naturally reacts to the varying physical coefficients, including the IOPs, the underwater depth, the sun intensity and the response curve of the sensor.We have integrated our method into the Unity engine HDRP to combine it with some of its physically-based rendering features.
Our method accounts for both the downwelling irradiance from the water surface, as well as the diffuse reflectance of the seafloor, improving its accuracy in shallow scenes.Inspired by underwater optics literature, we have proposed an expression for multiple scattering that models downwelling irradiance as a function of the distance to the surface.Our method also models the dielectric nature of the water surface, and reproduces both single scattering volumetric shadows and caustics.By working in the spectral domain, our results consistently achieve a closer match to ground-truth, pathtraced references.
Our method assumes homogeneous coefficients and an isotropic phase function.The actual phase function of oceanic water is a combination of its water molecules (isotropic), and suspended particles (strongly forward scattering, with anisotropy g = 0.924).In oceanography, it is common to model this with the Fournier-Forand phase function, which for efficiency in Monte Carlo simulations is often approximated with a weighted sum of Henyey-Greenstein phase functions [Mob94, DHV * 23].An interesting avenue for future exploration would be to devise constant-time expressions for both single and multiple scattering considering other analytical phase functions (such as Henyey-Greenstein).
Our method could be improved by tackling its other main limitations: multiple scattering does not consider occlusions or local light sources.Although this limits its applicability, this works well for the usually open geometry of underwater scenes.
Our model has been formally derived from radiative transfer theory and oceanography.Our methodology could be inspiring for other complex media, in which radiative transfer can be combined with heuristic knowledge from a particular field of knowledge (oceanography in our case) to achieve real-time performance.Moreover, another promising research direction would be to apply our model for inverse rendering problems, in order to estimate the physical parameters of water bodies from underwater photographs.
-world K d datasets.The earliest datasets of K d originate from the Swedish Deep Sea Expedition during which oceanographer Nils Gunnar Jerlov made measurements of the downwelling irradiance spectrum for different water bodies as a function of depth, which he later used to calculate K d values for various locations around the globe.Typically referred to as Jerlov Water Types [JER51,Jer76,Jer77], these data categorize the optical conditions of the world's oceans into 10 classes.The waters are divided in five open waters (I-III) (Figure 2, top), and five coastal (1-9C) (Figure 2, bottom).The figure shows their respective wavelength dependent coefficients of absorption σa, scattering σs, extinction σt , and diffuse downwelling attenuation K d .Note that the absorption and scattering coefficients corresponding to the K d values of Jerlov water types were obtained by reverse engineering using the Hydrolight software by Solonenko & Mobley [SM15].

Figure 2 :
Figure 2: Jerlov water types, five open (top) and five coastal (bottom).They are ordered from lowest (left) to highest (right) scattering to absorption ratio (single-scattering albedo).

Figure 3 :
Figure 3: The final radiance L results from the contributions of the surface of the object L S (black), and path radiance (dotted line).The path radiance is the integral of the medium in-scattering (orange) at each point.Each of the particles in the medium can contribute in-scattering, either directly (single scattering, green), or through several interactions (multiple scattering, purple).

Figure 5 :
Figure 5: At depth y, the total irradiance (E(y), red) results from adding the descending irradiance (E D (y), purple) and the upwelling irradiance, previously reflected on the ocean floor (E U (y), gold).

Figure 6 :
Figure6: The water surface radiance L W S accounts for the dielectric nature of the water surface.Snell's window is the circle in the water surface, above the viewer, outside of which there is total internal reflection.Inside it, the whole sky is refracted.Our method naturally simulates this phenomenon, as we can see above the shark.

Figure 7 :
Figure 7: The response functions of two sensors (top), Canon A2200 and Arriflex D21, are used to render the same underwater scene (bottom) of Jerlov IA, a very clear ocean water, typically of blue appearance to the human eye.The first sensor shows a stronger green channel response than the second one, which results in a clear difference in the RGB values they capture.

Figure 8
Figure8shows a reef scene for different open water types, where the 3D model was created from the images in Sea-thru dataset D3[AT19].Each row shows increasing depth.In the top row (y f = 4 m), we can see how the increasing single-scattering albedo (σs/σt ) from left to right correlates with stronger single-scattering

Figure 8 :
Figure 8: Real-time results of the different Jerlov open waters (I to II) in each column (increasing scattering to the right).Jerlov III is omitted due to low visibility.Each row shows increasing depth.In shallower depths (top), single scattering is much more prevalent.

Figure 9 :
Figure 9: Screenshot of our interactive application, where the user can change multiple parameters defining the scene and its depth, the type of water and its IOPs, the sensor characteristics, and the incoming light.Please refer to the supplemental video.

Figure 10 :
Figure 10: Real-time results of different Jerlov coastal waters (7C and 9C omitted due to low visibility); the rows show increasing depth.In these turbid coastal waters, radiance gets extinguished at around ten meters deep.On the other hand, in shallow depths, single scattering volumetric effects are more clearly visible.Scene Chen et al.Ours Jerlov 3C d3 0.138 0.059 Jerlov 3C d9 0.057 0.041 Jerlov I d5 0.333 0.117 Jerlov IA d9 0.249 0.107 Jerlov IB d5 0.239 0.107 Jerlov II d3 0.212 0.107 Average 0.205 0.090 Table 2: Average per-pixel RMSE of each method (left: Chen et al. [CZSZ22], right: ours) with respect to the references.These values correspond to the renders shown in Figure 11.Our method consistently produces closer matches to the ground truth, showing on average half the RMSE.

Figure 11 :
Figure 11: Comparison with the state-of-the-art method of Chen and colleagues (USSim) [CZSZ22], for six different scenes varying water type and depth.PT refers to the path-traced, ground-truth solution.Our method the ground-truth references much more closely.

Figure 14 :
Figure 14: Summary of the average RMSE and Delta E 2000for RGB (blue) and our spectral rendering method (orange), compared to ground-truth path tracing.For both metrics the error is significantly lower with our method.

©
2024 The Authors.Computer Graphics Forum published by Eurographics and John Wiley & Sons Ltd.

Figure 15 :
Figure 15: Irradiance profiles with depth of an infinite homogeneous body of various Jerlov water types.Lines: Monte Carlo simulation using the spectral IOPs from the literature [SM15].Crosses: irradiance E(y) approximation using only K d , E(y) = e −K d y (Equation (6)).

Table 1 :
Main symbols used in the paper.