Unifying radiative transfer models in computer graphics and remote sensing, Part II: A differentiable, polarimetric forward model and validation

The constellation of Earth-observing satellites continuously collects measurements of scattered radiance, which must be transformed into geophysical parameters in order to answer fundamental scientific questions about the Earth. Retrieval of these parameters requires highly flexible, accurate, and fast forward and inverse radiative transfer models. Existing forward models used by the remote sensing community are typically accurate and fast, but sacrifice flexibility by assuming the atmosphere or ocean is composed of plane-parallel layers. Monte Carlo forward models can handle more complex scenarios such as 3D spatial heterogeneity, but are relatively slower. We propose looking to the computer graphics community for inspiration to improve the statistical efficiency of Monte Carlo forward models and explore new approaches to inverse models for remote sensing. In Part 2 of this work, we demonstrate that Monte Carlo forward models in computer graphics are capable of sufficient accuracy for remote sensing by extending Mitsuba 3, a forward and inverse modeling framework recently developed in the computer graphics community, to simulate simple atmosphere-ocean systems and show that our framework is capable of achieving error on par with codes currently used by the remote sensing community on benchmark results.


Introduction
In Part 1 of this work, we examined the history and current state of the art in radiative transfer models in computer graphics and remote sensing and established our motivation for importing knowledge from computer graphics.In this paper, we take first steps toward putting these concepts into action.Our goal in this work is to develop a framework that can match the accuracy and capabilities of multi-angle, spectral, polarimetric instruments on future satellite missions and can leverage the inherent advantages of Monte-Carlo-based solvers (e.g.threedimensional geometry and spatial heterogeneity) as well as other special techniques imported from computer graphics (see Part 1).Although in this work we primarily focus on forward modeling, a medium-term goal of our research is to explore differentiable rendering as an engine for inverse modeling.Lastly, a long-term goal of this project is to release this software for use by the remote sensing community, so flexibility and ease of use are also important considerations.
To this end, we extend a forward and inverse modeling framework recently developed in the computer graphics community, Mitsuba 3 [1], to perform simulations of interest to remote sensing.We describe the pre-existing capabilities of Mitsuba 3 and its unique advantages in Section 2 and our extensions to it in Section 3. Lastly, we validate our framework on recent benchmark tests of simple atmosphere-ocean systems that conform to the plane-parallel layers assumption in Section 4 and show that our framework is capable of the accuracy expected for forward models in remote sensing.We plan to apply our framework to three-dimensional spatially varying media in future work.

Base capabilities of our framework
We chose to build our Monte Carlo radiative transfer framework on top of Mitsuba 3 [1] because it is uniquely suited for tackling this problem: • Physical accuracy.It is a physically based renderer that adheres to radiative transfer principles and solves related problems via Monte Carlo integration.It can optionally simulate polarization via Mueller calculus (see Section 2.1).• Differentiability.Derivatives of arbitrary radiative transfer simulations can be computed via differential radiative transfer (see Part 1, Sec. 4.6).• Modularity.Concepts such as sensors, phase functions, even rendering algorithms themselves, are abstracted into base classes that can be extended to implement specific variants (e.g. a perspective camera or a Rayleigh phase function).These variants are loadable "plug-ins" which can easily be swapped for one another.• Speed.The built-in DR.JIT [2] just-in-time compiler generates platform-optimized megakernels for parallel execution on CPUs or GPUs via LLVM or CUDA/OptiX, respectively (it can also run on CPUs via traditional thread pools).• Accessibility.It is open source, well maintained, and well documented.
Mitsuba 3 consists of a C++ and Python backend containing core classes, functions, and plug-ins and a Python frontend for calling them.Python bindings to the backend are automatically generated at compile time, allowing it to be loaded as a Python package and its classes and functions to be exposed to the user.The system is flexible: a user can either build higher-level functions in the backend and call them in the frontend, or build higher-level functions in the frontend from low-level backend components.For example, in the frontend, a user could load a specific plug-in (e.g. a Lorenz-Mie phase function) and run tests on that alone, load and render a complete scene (as shown in Section 4), or set up a gradient-descent loop for inverse modeling (which we will explore for the remote sensing context in future work).See the Mitsuba 3 tutorial suite for more examples [1].
The base system is an extensible framework with many built-in capabilities including a library of phase functions, bidirectional scattering distribution functions (BSDFs), emitters, sensors, and rendering algorithms that can simulate elastic scattering by surfaces and scattering media.The BSDF library includes a rough dielectric plug-in that can be used to simulate a rough ocean surface.While the base system supports scattering media that varies three-dimensionally in albedo and density (but not phase function), we leave exploration of spatially varying media for future work.

Representing polarization
Mitsuba 3 uses Mueller calculus to represent polarization states, which is also common in remote sensing.Mueller calculus uses a Stokes vector s ∈ R 4 to describe polarization: where I is the total radiance, Q and U describe the linearly polarized radiance and its orientation relative to a reference frame, and V describes circularly polarized radiance; ϵ and μ are the electric permittivity and magnetic permeability of the medium, respectively. 1Graphically, one can think of the Stokes vector as representing the average polarization ellipse, with orientation angle ψ and ellipticity angle χ (see Fig. 1).
For the reference frame of Q and U, Mitsuba 3 uses a right-handed coordinate system with the z-axis pointing along the direction of light propagation (as illustrated in Fig. 1).To be consistent with the results presented in Chowdhary et al. [5], for radiance arriving at the sensor we align the x-axis of the reference plane with the plane spanned by the sensor direction and the global z-axis of the scene (see Fig. 5, this plane is shaded in blue beneath the blue sensor direction arrow), and a positive sign of Q indicates parallel to this plane.See the Mitsuba 3 documentation for more details on its polarization conventions [1].
The degree of polarization (DP) and degree of linear polarization (DLP) of the light are defined as: where DP = 0 is completely unpolarized, 0 < DP < 1 is partially polarized, and DP = 1 is fully polarized light.Any 4-vector that produces DP < 0 or DP > 1 is not a physically valid Stokes vector.In this work, we focus solely on linear polarization since circular/elliptical polarization is generally ignored in the remote sensing context (see Gassó and Knobelspiesse [6] for a recent counterargument) and modern-day polarimeters typically only collect linear polarization data.
A Stokes vector can be transformed into another by multiplication with a Mueller matrix M ∈ R 4×4 such that s o = M s i , which typically represents an interaction of light with a material, such as an optical element (e.g.polarizer), surface, or particle.A Mueller matrix can also rotate the reference frame of a Stokes vector.
Stokes vectors and Mueller matrices are convenient to use in a rendering context because they are natural and physically valid extensions of scalar, unpolarized radiancein a path or light tracer, any scalar radiance quantity can be substituted for a Stokes vector and any scatter event can be substituted for a Mueller matrix.Stokes vectors may be added to represent a superposition of rays of light, and their addition correctly accounts for their cumulative polarization states.However, there is the extra complication of tracking the reference frames of each Stokes vector and rotating them appropriately before addition or multiplication.

Our extensions to the framework
We implemented several extensions to Mitsuba 3 in order to make it useful for a remote sensing context that include, but are not limited to: support for polarized phase functions, Rayleigh phase functions, Lorenz-Mie phase functions with optional distributions of particle size (see to clarify that our Stokes vectors carry radiance, which should be proportional to the flow of electromagnetic energy, which should be proportional to the medium properties ϵ and μ.In practice, this detail becomes relevant upon refraction at a dielectric surface, which changes the solid angle of the diverted beam and its radiance by a factor of the squared relative index of refraction (sometimes called the n 2 law) [3,4].Section 3.2), and a custom path tracer tailored for special handling of atmosphere-ocean systems (see Section 3.1).

Special handling of next-event estimation beneath an ocean surface
We use path tracing to estimate the radiance arriving at the sensor.A path is a set of n positions ("path vertices") throughout a scene where valid scatter events may occur, i.e. on surfaces and within scattering media, beginning at the sensor and ending at a light source.Path tracing estimates the radiance arriving at the sensor by constructing lightcarrying paths vertex by vertexthe next vertex is chosen by sampling a new direction for the path proportional to the current vertex's out-scattering radiance, preferably.
For example, the outgoing radiance L o leaving a surface at x in direction ω2 is given by: (3) which is the sum of emitted radiance L e and reflected radiance L r .The reflected radiance is an integral of incoming radiance L i (x, ω i ) from all directions ω i over the sphere of directions Ω.This light is weighted by the surface's cosine-weighted BSDF f s .In a vacuum, the incoming and outgoing radiances are related as where the ray tracing operator r(x, ω) = x z := x − zω returns the closest surface point along the ray from x in direction − ω.To account for scattering media, Eq. (4) becomes: (5) The outgoing radiance at the end of the ray L o (r(x, ω), ω) is attenuated by the transmittance, and we must additionally integrate the inscattered radiance at each point x t := x − tω through the medium up to the end of the ray.The medium's extinction coefficient is σ, and τ = ∫ d 0 σ(x a +tω) dt is the medium's optical depth along the line segment from x a to x b of length d (we ignore emissive scattering media for brevity).We assume the phase function f p integrates to the singlescattering albedo over the sphere of directions.A simple path tracing algorithm would first sample a random free-flight distance t, and depending on whether this stays within the medium or hits a surface, would then sample a new direction ω i proportional to either the phase function f p or the BSDF f s .
Radiance from sources other than emissive surfaces, e.g. a directional light over some solid angle or an environment map, can be gathered when a path escapes the scene.The efficiency of this simple routine depends on the roughness of materials within the scene and the relative size and placement of light sourcesscenes with a mixture of concentrated (small in size or solid angle) light sources and low-roughness materials are particularly challenging since the integrand becomes more peaky.When this happens, the probability that a path randomly sampled from the BSDF will hit a small light source is low, but when it does, the contribution to the Monte Carlo estimator will be huge, which results in high variance (see Part 1, Fig. 2a).An improvement commonly used in graphics is to include next-event estimation (NEE; or local/ directional estimate), and to combine sampling strategies using multiple importance sampling (MIS), as we discussed in Part 1.
Unfortunately, a straightforward application of NEE and MIS provides little benefit in typical atmosphere-ocean systems since they contain particularly tricky light paths that require extra thought: because the ocean surface is refractive, scattering events that occur beneath the ocean surface cannot leverage traditional NEE.To see why, consider a path that has traveled beneath the ocean surface and is situated somewhere within the scattering ocean body: during NEE, a new ray is generated in the direction of the emitter; since this ray does not directly "see" the emitterthe ocean surface is in the wayit returns zero radiance.Simply ignoring this occlusion is unfortunately also not an option since the ray would be refracted at the ocean-atmosphere interface away from the emission source and hence gather no radiance (Fig. 2b).NEE does not help here because light from above the ocean surface cannot reach a scattering event within the ocean body along a straight line.This type of light path is in fact a caustic path.In graphics there are known general solutions to this problem such as bidirectional path tracing [7], Metropolis light transport [8,9], and photon mapping techniques [10][11][12][13], as well as specialized NEE techniques that account for one intervening reflection or refraction event [14][15][16].While general, these approaches tend to be too heavy-weight for our more simplified setting.
We instead developed an importance sampling scheme to make such Positions and directions match those defined in Eq. ( 8), and the BSDF lobe at x 1 for each incoming direction is shown in purple.
refractive connections which remains light-weight and efficient by taking advantage of the plane-parallel assumptions in the atmosphereocean system (Fig. 2(c,d)).Instead of splitting incident radiance into direct and indirect parts, as is done in traditional NEE, we split incident radiance into caustic subpaths that connect a scattering vertex within the ocean to a light source via two straight path segments, and all other subpaths: We approximate the incident radiance along the "other" subpaths using a standard 1-sample Monte Carlo estimator that samples a new direction ω i proportional to the phase function f p .The subset of length-2 caustic paths expands to: where ω 1 is the direction entering the scattering event at x (from the ocean surface) and ω 2 is the emission direction of the light source; then x 1 = r(x, ω 1 ) is a point on the ocean surface and x 2 = r(x 1 , ω 2 ) is a point on the light source.
In the atmosphere-ocean setting, it is common to model incoming solar radiance using a light source that emits constant radiance L e along a single direction ω ℓ from an infinite distance.Using a Dirac delta distribution, defined as ∫ δ(x)dx = 1 and δ(x) = 0 for all x ∕ = 0, its emitted radiance can be expressed as L e δ(ω − ω ℓ ).Substituting this for L e (x 2 , ω 2 ) in Eq. 6 collapses the inner integral, yielding: ) L e dω 1 . ( Here, x 2 corresponds to the point where the ray from x 1 traveling in direction −ω ℓ exits the atmosphere.We can then define a simple 1-sample Monte Carlo estimator for Eq. ( 8) by sampling a direction ω 1 (〈L〉 := 〈L caustic (x, ω)〉 for brevity): The last question is how to sample a direction ω 1 .Standard NEE can only be successfully applied at the ocean surface point x 1 (Fig. 2b).This would correspond to importance sampling Eq. ( 9) with a PDF proportional to the phase function f p .We found this leads to significant variance due to variation in the BSDF term f s , so we instead developed a strategy that importance samples such two-segment paths proportional to the BSDF: we use the light's emission direction ω ℓ and the BSDF's existing sampling strategy to sample a new direction ω 1 that refracts into the ocean body (Fig. 2c).By the reciprocal nature of light transport, any path starting at x and traveling along the reverse of this BSDF-generated ray will refract toward the emission direction and gather radiance successfully.This sampling step gives us an initial ray to start our subpath: its origin is x and its direction −ω 1 is the reverse of the BSDF-generated ray.When the path reaches the ocean surface, we evaluate the surface BDSF in the emitter direction ω ℓ , and additionally evaluate the transmittance through the ocean and atmosphere (Fig. 2d).This method assumes that the surface is planar and spatially homogeneous, i.e. the location of entry/exit at the surface does not matter.It is also important not to double count the radiance along length-2 caustic subpaths: if a path generated by phase function sampling (to estimate the second part of Eq. 6) happens to refract out of the ocean at the surface, traditional NEE should not be performed at that surface scattering event because its radiance is accounted for in Eqs. ( 7) and ( 8).We must also account for this sampling step by dividing the gathered radiance by the sampling PDF.If f s is a rough microfacet BSDF, a common choice for modeling the ocean surface, a sampling strategy with a PDF perfectly proportional to f s is available [18,19], which causes f s (ω 1 ) and p(ω 1 ) to cancel each other out.Any remaining variance in Eq. ( 9) is due to the phase function and transmittance evaluations.
We validate our strategy on a test case based on the AOS-II atmosphere-ocean model described in Section 4.1 in Fig. 3.The scene is composed of an ocean surface modeled by a microfacet BSDF and an ocean body beneath filled with a homogeneous scattering medium.We estimate the ocean-leaving radiance (ignoring surface reflection) since that is the only component handled by our custom NEE method, and also examine the effects of varying the roughness of the ocean surface.We compare to the standard NEE strategy that generates successful direct connections to the light only at the ocean surface.Our method is relatively invariant to roughness, whereas the basic NEE exhibits higher variance as the surface approaches specular.In every case, our method results in a variance reduction of three to four decimal orders of magnitude.It should be possible to further reduce the remaining variance by performing MIS between traditional NEE (importance sampling proportional to the phase function) and our NEE (importance sampling proportional to the BSDF), though we leave this for future work.
We believe our two-segment NEE strategy is just the starting point for further improved path construction in plane-parallel systems.Our current strategy can be seen as a (long photon volume × point) estimator [20] specialized for a planar refractive surface.This suggests that additional variance reduction may be possible by developing e.g.(photon volume × long beam) strategies, and more.We hope to explore this in the future.

Lorenz-Mie scattering
We follow the implementation of Mishchenko and Yang [17] and Mishchenko et al. [21] for Lorenz-Mie scattering: they describe in detail how to compute the phase function and scattering and extinction cross sections for a given parameter set by summing the necessary spherical functions using stable recurrence relations and integrating those Fig. 3.We demonstrate the variance reduction of our NEE (solid) versus basic NEE (dashed) on the AOS-II atmosphere-ocean model described in Section 4.1 at θ ℓ = 60 ∘ , θ v = 60 ∘ .The scene is composed of a rough ocean surface and a homogeneous ocean body.Here, we have isolated just the ocean-leaving radiance (and ignore surface reflection) to isolate the tricky scenario in question and highlight the difference between the methods.We show the effect of modifying the roughness (α), or mean square slope, of the ocean surface across a range of values: our method performs consistently across all values, while basic NEE becomes more and more noisy as the roughness decreases (the surface approaches specular).In all cases, our method consistently produces 1,000-10,000 × lower variance than traditional NEE.
quantities over a range of particle sizes using Gaussian quadrature.We implemented size distributions as a base class so that arbitrary distributions can easily be added to our framework, such as gamma, lognormal, power-law, etc.
Because Lorenz-Mie scattering is prohibitively expensive to evaluate for every scatter event, we pre-tabulate Lorenz-Mie phase functions over a range of scattering angles to use at runtime.A Lorenz-Mie phase function typically exhibits strong forward scattering, and we would like to create a PDF for sampling directions from this distribution that matches its shape.To that end, we used an existing capability of Mitsuba 3 that computes a piecewise-linear PDF proportional to tabulated data (see e.g.Pharr et al. [22] for details on constructing, sampling, and evaluating such a PDF).We use the top-left Mueller matrix value to construct the tabulated PDF and draw directional samples from this distribution at runtime.We then evaluate the Mueller matrix elements by linearly interpolating each tabulated element in the sampled direction and weight the result by the sampling PDF.While the tabulation resolution can be chosen arbitrarily, for our results in Section 4 we use tabulated data in 0.1 ∘ increments for scattering angles from 0 ∘ to 10 ∘ and 1 ∘ increments from 10 ∘ to 180 ∘ .

Results
To validate the accuracy of our framework, we compare to all of the benchmark results presented in Chowdhary et al. [5].These benchmark results are simulated top-of-atmosphere and just-above-surface polarimetric measurements for atmosphere-ocean systems over several wavelengths and a range of sun and viewing geometries.Because the majority of the radiative transfer codes evaluated in Chowdhary et al. [5] are deterministic and therefore operate under the plane-parallel layers assumption, we made our atmosphere-ocean systems also conform to that assumption for comparison purposes.

Scene descriptions and parameters
The scene composition of the four atmosphere-ocean systems is illustrated in Fig. 5.Each atmosphere-ocean system contains some combination of a homogeneous scattering atmosphere, a wind-ruffled ocean surface, and a homogeneous scattering and absorbing ocean body (see Table 1 for a summary).We list the material parameters in Tables 2 and 3, which match those used in Chowdhary et al. [5]: for the atmosphere we ignore the effects of molecular absorption by gases and borrow optical depths from literature [23]; for the ocean body we assume the ocean is composed of pure seawater and borrow its scattering and absorption properties from literature [24,25] for AOS-I, AOS-II, and AOS-III (the ocean body for AOS-IV is more complex and described in detail below); the ocean floor is modeled by a fully absorbing surface.The ocean surface is modeled by an isotropic Cox-Munk BSDF [26] with a roughness value corresponding to a wind speed of ∼7 m/s at 12 m above the surface, which translates to a mean square slope α 2 = 0.03884 (the default parameter in Mitsuba 3 is the root mean square slope, α).In order to match the results of Chowdhary et al. [5] as closely as possible, we omit the shadowing-masking term for microfacets and do not enforce conservation of energy.
The final AOS-IV model differs from the other models in that it uses a bio-optical phase function [27] for the ocean that is a weighted average of pure seawater, detritus and mineral particles, and phytoplankton particles.In the following, F indicates a polarized phase function that expands into a Mueller matrix.This phase function is given by: The pure seawater scattering matrix F w is modeled by Rayleigh scattering; b w and b p are the scattering coefficients of pure seawater and particles, respectively.The particle scattering matrix F p is further subdivided into detritus & mineral and phytoplankton contributions, given by: The scattering matrices F ph and F dm are modeled by Lorenz-Mie scattering with a power-law (or Junge) distribution of particle sizes; a dm is   the fraction of particles represented by detritus & minerals3 ; σ ph and σ dm are the scattering cross sections of phytoplankton and detritus & mineral particles, respectively.We list the parameters for this phase function in AOS-IV, including those for Lorenz-Mie scattering, in Table 3. Chowdhary et al. [5] helpfully provide their tabulation of this phase function for AOS-IV in their supplemental material; we show that we can also compute this phase function from the given parameters using our own Lorenz-Mie scattering implementation in Fig. 4. For the AOS-IV results presented in Fig. 20,19,13,12, we use the tabulated data provided by Chowdhary et al. [5] to achieve the most faithful comparisonthey provide data in 0.1 ∘ increments for scattering angles from 0 ∘ to 10 ∘ and in 1 ∘ increments from 10 ∘ to 180 ∘ .

Comparison to benchmark results
Chowdhary et al. [5] use eGAP [28,29] to compute their reference results, an adding-doubling radiative transfer code currently used at NASA/GISS, to an accuracy of 1 × 10 −6 for models AOS-I, AOS-II, and AOS-III and 1 × 10 −5 for model AOS-IV.They provide results in units of reflectance, defined as: where s i ∈ {I, Q, U} is the final Stokes vector component (i.e.radiance value) recorded by the sensor, L e is the (unpolarized) incoming extraterrestrial solar flux (set to π in all cases for simplicity), and μ ℓ = cosθ ℓ (see Fig. 5).
We compare reflectance values to the reference results in terms of relative percent error, defined as: and degree of linear polarization in terms of absolute error (since the quantity itself is a percentage): Fig. 4. We compare our Lorenz-Mie scattering computations to two references: results computed using Mishchenko and Yang [17]'s code (top row) and tabulated results provided in Chowdhary et al. [5] (bottom row).We use the Lorenz-Mie scattering parameters in Table 3 to compute the four unique Mueller matrix elements of a collection of spherical particles representing detritus & minerals in the ocean (top-left), phytoplankton in the ocean (top-right), and mixtures of these phase functions at two different ratios (bottom-left and bottom-right, a dm is the fraction of detritus & mineral particles and a ph = 1 − a dm is the fraction of phytoplankton particles).Our results match both references very closely, with some minor artifacts from Gaussian quadrature visible in F 34 .We use 200 Gaussian quadrature points to integrate over the entire domain of particle sizes for our results (Mishchenko and Yang [17]'s code subdivides the integration domain into 100 subintervals and uses 100 Gaussian quadrature points to integrate over each subinterval; the number of quadrature points used in Chowdhary et al. [5] is unknown).
Fig. 5.Each atmosphere-ocean scene in Chowdhary et al. [5] is some combination of a homogeneous atmosphere, a rough ocean surface, and a homogeneous ocean body stacked in parallel layers (see Tables 1 to 3 for the composition and parameters of each).We illustrate the geometric setup of light and viewing angles.
Our results, shown in Figs.7 to 20, were computed using our Mitsuba 3 framework with ∼5.37 × 10 8 samples for all parameter sets.Since our Monte Carlo code would naturally exhibit some random noise, we conducted 20 trials per parameter set using a different random seed for eachthe mean of each set of trials is connected by dark purple line and a kernel density plot ("violin" plot) illustrating the spread of each set of trials is shown in light purple for each parameter set.The mean of all trials is also equivalent to a Monte Carlo estimate using 20 ×5.37 ×10 8 total samples for that parameter set.We also ran convergence tests with increasing sample counts in order to assess the convergence behavior of our Monte Carlo estimates.This provides an indication of whether any remaining error at high sample counts is due to variance (which would vanish at even higher sample counts) or bias (a pernicious error that would not vanish at higher sample counts).Examples shown in Fig. 6 show the convergence of the reflectance estimate for one parameter set of the AOS-I and AOS-II models.The AOS-I estimate appears to be converging to the reference result, whereas the AOS-II model exhibits some slight bias, converging to |ε R | ∼ 0.04%.
These validation tests demonstrate that our frameworkwhile exhibiting a small amount of bias in some casesgenerally reaches an error threshold on par with deterministic codes currently used by the remote sensing community.Therefore, it should be acceptable for use in remote sensing applications.

Technical details and performance
All results were generated on a cluster node with 2.90 GHz Intel Xeon Gold 6226R CPUs.Mitsuba 3 can be compiled to run either on CPU (via traditional thread pools or LLVM infrastructure) or GPU (via CUDA/ OptiX) -we used the CPU variant with LLVM optimizations.This variant, as well the CUDA/OptiX-based GPU variant, rely on the integrated DR.JIT [2] just-in-time compiler to generate platform-optimized megakernels for parallel execution.See Jakob et al. [2] for more technical details.
Due to the independence of traced paths, Monte Carlo codes such as ours are "embarrassingly parallel" -in theory, the minimum runtime would be the cost of tracing one path (plus setup of inputs and processing of outputs), if enough cores were available.Therefore, the "efficiency" of our code is primarily a product of 1) how many samples are taken, and 2) how many cores are available.Running one trial of our results in Section 4 (i.e.calculating one final radiance estimate corresponding to one parameter set) with ∼5.37 × 10 8 samples (i.e. each sample is a full path) and 64 cores generally takes 10 to 40 s total for AOS I-III and 1 to 2.5 min for AOS-IV, depending on the model and parameters.
If comparing performance in terms of wall clock time, most deterministic codes will have the advantagehowever, our framework comes with other advantages, including the inherent advantages of Monte-Carlo-based models discussed in Part 1, as well as the frameworkspecific advantages discussed in Section 2.

Discussion, limitations, & future work
Our framework can theoretically handle most of the complex lightbased phenomena listed in the first paragraph of Part 1, with the caveat that appropriate geometry and material properties would need to be provided.Some exceptions include inelastic scattering (e.g.fluorescence), as wavelengths of propagating light are assumed to be fixed, and spatially varying phase functions in scattering media, which would be necessary to simulate media with a spatially varying distribution of particle sizes.These features could be added in theory, but would require significant modifications to the base framework.
The extensive suite of validation tests in Section 4 establishes that our current framework is an acceptable forward model for simple atmosphere-ocean systems.Although these results conform to the assumption of plane-parallel layers for the sake of comparison, we hope to illustrate the relative advantages of our framework on more complicated cases with spatially varying scattering media and threedimensional geometry in the future.We hope to incorporate several of the relevant threads of research discussed in Part 1, Sec. 4, such as multiple-scattering microfacet models and lower-variance estimators of spectrally and spatially varying scattering media, into our framework in order to push the accuracy and efficiency of our forward model beyond the current state of the art in remote sensing.Concurrent work, a project called Eradiate [33], has also begun to develop a forward modeling framework for remote sensing that uses Mitsuba 3 as its underlying backbone; however, their project aims to build a user-friendly wrapper of remote-sensing-oriented scripts and tools that use Mitsuba 3 "under the hood."Our primary focus is to improve the inverse modeling capabilities of Mitsuba 3 on polarimetric measurements and our forward model is a stepping stone from which to begin that exploration.
Differentiable rendering with Mitsuba 3 has already shown promising results in retrieving the roughness parameter of a microfacet model and spatially varying optical depth of a scattering medium in an unpolarized context [34].Therefore, a natural place for us to begin is to test the retrieval of similar parameters from the simple atmosphere-ocean scenes used in Section 4 in a polarized context.We will also incorporate our Lorenz-Mie scattering model in order to retrieve polarization-sensitive parameters such as size and complex indices of refraction of aerosols.
While differentiable rendering has many desirable qualities, it is not a panacea: there are several practical issues that will need to be examined.For instance, differentiable rendering inherits the drawbacks of whatever optimizer (e.g.stochastic gradient descent, ADAM [35]) it is used with, such as converging to a local optimum instead of a global one, choosing an appropriate learning rate, momentum, regularizer, etc.Since it is a Monte Carlo approach, determining an appropriate sample count for gradient estimates introduces yet another hyperparameter into the optimization.There is also the question of computational "waste": a naive optimization would forget intermediate results computed in prior steps and perform a fresh, costly Monte Carlo simulation at every iteration, even though the parameter changes may be very small.To mitigate this, some works have begun to design optimization methods for differentiable rendering that reuse information from prior steps [36][37][38].Lastly, a guiding principle of this work is accessibility, so our framework will be made available publicly and promptly at the first author's Github page so that it can be used by the scientific community.

Declaration of Competing Interest
None.

Fig. 1 .
Fig. 1.Light travels along the z axis with linear polarization aligned with the y axis (left) and elliptical polarization with orientation angle ψ and ellipticity angle χ (right).

Fig. 2 .
Fig. 2. Traditional next-event estimation (NEE) vs. our modified version to gather radiance within the scattering ocean body in a simple atmosphere-ocean system.Positions and directions match those defined in Eq. (8), and the BSDF lobe at x 1 for each incoming direction is shown in purple.
Our target error regions of | ε R | < 0.01% and |ε DLP | < 0.1% are shown in light blue on each plot, which also gives a visual cue of instances where the y-axis has been scaled to encompass outliers.For reflectance, the mean errors for AOS-I and AOS-III remain within our target error threshold of |ε R | < 0.01% in most cases.Individual trials for AOS-I and AOS-III mostly remain within |ε R | < 0.1% for top of atmosphere and |ε R | < 0.5% for just above surface.AOS-II exhibits more consistent error, with error means typically in the range |ε R | < 0.1%, but sometimes straying into the |ε R | < 0.5 −1% range as the viewing angle θ v increases.Error means for AOS-IV mostly remain within |ε R | < 0.5% for top of atmosphere and |ε R | < 1% for just above surface, with some outliers in |ε R | < 5% for just above surface.This may be due to a difference in implementation in how the tabulated phase function is sampled and evaluated.For degree of linear polarization, the majority of all trials remain comfortably within our target error threshold of |ε DLP | < 0.1% for AOS-I and AOS-III, except for some outliers generally within |ε DLP | < 0.5% for AOS-I and AOS-III just above surface.AOS-II trials meet our target threshold of |ε DLP | < 0.1% for the ϕ v = 0 ∘ case; for other geometries the error is generally within |ε DLP | < 0.5% and in a few cases within | ε DLP | < 1% (at θ ℓ = 60 ∘ , ϕ v = 60 ∘ and θ ℓ = 60 ∘ , ϕ v = 240 ∘ ).AOS-IV consistently meets our target error threshold of |ε DLP | < 0.1% for top of atmosphere except for some outliers within |ε DLP | < 0.5%; for just above surface, the error means remain consistently within |ε DLP | < 0.5% and some outliers within |ε DLP | < 1%.

Fig. 6 .
Fig. 6. (Absolute value of) Relative percent error in reflectance, |ε R |, for AOS-I and AOS-II models at top of atmosphere, for θ ℓ = 30 ∘ , θ v = 30 ∘ , ϕ v = 60 ∘ .The mean of each set of 20 random trials is connected by dark purple line and the standard deviation of the trials is shown as error bars; both axes are a log scale.The AOS-I model shows a pattern of converging toward the reference result, while AOS-II shows a slight bias that approaches |ε R | 0.04%.

Fig. 7 .
Fig. 7. Relative percent error in reflectance (ε R ) for Chowdhary et al. [5] AOS-I model at top of atmosphere, for all parameter sets.The mean of each set of 20 random trials is connected by dark purple line and a kernel density plot ("violin plot") illustrating the spread of each set of trials is shown in light purple for each parameter set.Our target error threshold of |ε R | < 0.01% is shown in light blue.

Fig. 8 .
Fig.8.Same as Fig.7, but with sensor at just above surface level instead of top of atmosphere.

Fig. 9 .
Fig.9.Same as Fig.7, but for AOS-II model (note that there is no distinction between just above surface and top of atmosphere for AOS-II because there is no atmosphere).

Fig. 12 .
Fig. 12. Same as Fig. 7, but for AOS-IV model at top of atmosphere.

Fig. 13 .
Fig. 13.Same as Fig. 7, but for AOS-IV model at just above surface.

Fig. 14 .
Fig. 14.Absolute error in degree of linear polarization (ε DLP ) for Chowdhary et al. [5] AOS-I model at top of atmosphere, for all parameter sets.The mean of each set of 20 random trials is connected by dark purple line and a kernel density plot ("violin plot") illustrating the spread of each set of trials is shown in light purple for each parameter set.Our target error threshold of |ε DLP | < 0.1% is shown in light blue.

Fig. 15 .
Fig. 15.Same as fig.14, but with sensor at just above surface level instead of top of atmosphere.

Fig. 16 .
Fig.16.Same as Fig.14, but for AOS-II model (note that there is no distinction between just above surface and top of atmosphere for AOS-II because there is no atmosphere).

Fig. 17 .
Fig. 17.Same as Fig. 14, but for AOS-III model at top of atmosphere.

Fig. 18 .
Fig. 18.Same as Fig. 14, but for AOS-III model at just above surface.

Table 1
Scene composition for atmosphere-ocean models.

Table 2
Scene parameters for atmosphere, ocean surface, and ocean body in all atmosphere-ocean models, where applicable, per wavelength (λ).

Table 3
[5]ameters for AOS-IV model's bio-optical phase function.Both Lorenz-Mie phase functions are calculated at λ = 550 nm and are assumed to be spectrally invariant.The low chlorophyll mix ratio is used for λ = 350 nm and λ = 450 nm and the high chlorophyll mix ratio is used for λ = 550 nm and λ = 650 nm (see Chowdhary et al.[5]for discussion).Particle sizes are integrated over the domain 0.01 − 100μm.Abbreviations: scattering coefficient (b), scattering cross section (σ), index of refraction (IOR), mix ratio of detritus & minerals to phytoplankton (a dm ).