1 Introduction↩︎

Light transport in volumetric scenes is complex. As light propagates through participating media, it is continuously scattered and absorbed, often undergoing multiple scattering events. These interactions generate long, intertwined transport paths that span the volume and tightly entangle illumination with the spatially varying optical properties of the medium. The radiance observed at any point is therefore not the result of a small number of local interactions, but the cumulative outcome of a dense superposition of indirect light contributions. Computational models of this process typically rely on radiative transfer simulations that use stochastic sampling of light paths [1].

Here, we study the corresponding inverse problem: given image observations of a scene containing participating media, we aim to reconstruct its spatially varying optical properties. This challenging task is commonly formulated as an optimization problem, in which volumetric properties are estimated by minimizing a reconstruction objective using gradient-based methods. Existing approaches can be broadly divided into two categories. First, stochastic light path sampling is made differentiable [2], [3], enabling the optimization to account for full radiative transfer (Fig. 1c). However, these methods demand substantial methodological and engineering effort to achieve acceptable efficiency [4]–[6] and to control bias in the estimated gradients [7]. In a second line of work, simplified models of light transport are employed, often in combination with neural scene representations [8]–[11] (Fig. 1a, b). While these approaches benefit from general-purpose, easy-to-use machine-learning frameworks, they typically fail to capture accurate optical properties under global illumination. To date, these two lines of work have remained largely disjoint.

In this work, we explore a way to reconcile the generality of physically complete light transport formulations with the flexibility and practicality of modern machine learning. Our approach shifts complexity from transport-specific algorithm design to learned representations and optimization, thereby decoupling physical modeling from specialized numerical treatment.

Specifically, we represent both the optical properties of the participating media and the full light field of the scene using neural fields [12], and formulate their joint estimation as a continuous optimization problem (Fig. [fig:teaser]). Crucially, this optimization involves only minimal explicit simulation of light transport (Fig. 1d). Global illumination is incorporated through a residual objective derived from the Radiative Transfer Equation (RTE), following a physics-informed formulation based on local differential operators [13]–[15]. Although the supervision signals induced by the RTE are local, their propagation through optimization implicitly resolves the dense superposition of light paths, leading the joint estimate of the light field and the medium’s optical properties to converge to a globally consistent light transport equilibrium without explicit transport modeling. The input observations and illumination conditions are incorporated through corresponding data constraints. In addition, we include a data term based on the Volume Rendering Equation (VRE), which integrates radiance along primary viewing rays only. This term mitigates the well-known low-frequency bias of physics-informed neural networks [16], [17] by coupling the reconstructed light field to view-ray integrals (we avoid the term "spectral bias", commonly used in the machine learning literature, to prevent confusion with the wavelength dependence of light). Importantly, explicit higher-order light transport simulation is avoided, since the RTE-based formulation accounts for global illumination. This substantially decouples inverse rendering research from renderer-specific algorithm design.

Using our formulation, we demonstrate the feasibility of reconstructing spatially varying, color-resolved scattering, absorption, and phase function parameters from multi-view observations under known illumination conditions, with competitive performance. Moreover, because the formulation relies on general-purpose neural optimization rather than transport-specific algorithmic machinery, it naturally extends beyond reconstruction to generative modeling of participating media with physically meaningful optical properties under full global illumination, learned directly from image observations. In summary, our contributions are:

• A novel formulation for volumetric inverse rendering that accounts for full global illumination with only minimal reliance on transport-specific simulation.

• A joint estimation framework for a scene’s optical properties and full light field, both represented as neural fields and constrained by first principles of global light transport.

• Experimental demonstrations of optical property reconstruction and the application of the formulation to generative modeling of physically meaningful participating media.

All source code and trained models are available at : https://neural-radiative-transfer.mpi-inf.mpg.de/

2 Related Work↩︎

Here, we review related work on inverse light transport (Sec. [sec:inv-rte]), physics-informed learning (Sec. [sec:pinn]), and physical generative models (Sec. [sec:related-generative]).

2.1 Inverse Light Transport↩︎

2.1.1 Discretization-Based Inverse Radiative Transfer↩︎

Inverse radiative transfer problems have traditionally been solved by numerically discretizing the radiative transfer equation [18], [19]. Traditional methods typically represent the light field and the optical properties using linear basis expansions defined on a computational mesh, resulting in a finite-dimensional approximation of the governing equation [20]. The optical properties are then inferred by iteratively solving the resulting system. These approaches have been successfully applied in fields such as atmospheric remote sensing [21], [22], medical imaging [23] and constitute the foundation of many established radiative transfer frameworks [24]. While traditional methods rely on discretization, our proposed method is based on a physics-informed learning framework, where we represent both the light field and the optical properties as neural networks. The governing physics and observational constraints are then enforced through the optimization objective.

2.1.2 Differentiable Monte Carlo Light Transport↩︎

Inverse volume rendering is traditionally addressed by differentiating stochastic light transport simulations, typically based on path tracing, to optimize the optical properties of participating media (Fig. 1c). Gradients of image measurements with respect to scene parameters are estimated by differentiating Monte Carlo estimators of the rendering equation, enabling optimization under physically complete light transport [2], [25], [26]. This framework has been applied to volumetric and translucent materials, allowing inverse estimation of optical properties in participating media [27]–[31]. While conceptually general, differentiable Monte Carlo methods are computationally expensive and often limited by high gradient variance, motivating extensive work on improving efficiency and stability [4]–[6]. Moreover, multiple scattering requires differentiating an expectation over a stochastic branching process. Obtaining gradients which are useful in optimization requires substantial methodological and engineering effort and remains an active area of research [7], [32], [33]. Our approach sidesteps these challenges by avoiding stochastic path sampling altogether, while still accounting for physically complete light transport.

2.1.3 Inverse Rendering with Learned Representations↩︎

Inverse rendering approaches based on learned representations commonly adopt simplified light transport models. A prominent example is the emission–absorption formulation (Fig. 1a), which represents volumes as purely emissive and neglects explicit illumination and scattering [8], [34], [35], making it amenable to automatic differentiation without special treatment. Single-scattering approaches constitute the next level of modeling complexity by accounting for direct illumination and a single light bounce (Fig. 1b). A common paradigm in this setting is to infer surface-based material properties, such as albedo and normals [36]–[41]. Accuracy can be further improved by explicitly modeling incident illumination using specialized components [9], [42]–[44] and by introducing physics-inspired regularization terms that encourage consistency with simplified reflectance models [45]. Recent work has also explored learning-based solutions for volumetric participating media with anisotropic scattering [11], but remains limited to restricted transport models. Multi-bounce illumination in inverse rendering has been considered more rarely. In this setting, surface-based representations can leverage precomputed radiance transfer [46], caching techniques [47], or explicit global-illumination decomposition strategies [48], [49]. For volumetric scenes, however, existing approaches typically rely on strong assumptions, such as homogeneous media [50], or on simplified treatments of higher-order light transport [10]. In contrast, our approach does not rely on such restrictive assumptions and directly targets physically complete volumetric light transport.

2.2 Physics-Informed Learning for Light Transport and Scene Reconstruction↩︎

Physics-informed neural networks (PINNs) integrate physical laws, typically expressed as (integro-)differential equations, into neural network training by penalizing violations of these equations in the loss function [13]. Incorporating the governing equations of volumetric light transport into such frameworks has been explored in both forward and inverse settings [51]–[53]. These works focus on ultra low-frequency optical properties and light fields, and often are further restricted to low-dimensional (1D or 2D) settings, whereas our approach operates in full 3D and aims to reconstruct spatially complex, high-frequency volumetric structure, as required in inverse graphics. This is achieved using an explicit integration scheme along primary rays, which helps alleviate the low-frequency bias of neural networks in general [17] and of PINN-based formulations in particular [16].

Physics-informed learning has also been explored for light transport in surface-based scenes, in both forward [14] and inverse [15] settings. In these surface-based formulations, non-local transport through free space is naturally handled via stochastic ray sampling, whereas our volumetric formulation enforces light transport locally throughout the volume by imposing the RTE as a pointwise constraint. Optimizing hard geometry, i.e., surfaces with sharp discontinuities in PINN-based inverse rendering remains challenging and has not yet been demonstrated, whereas our method naturally supports the optimization of soft geometry (e.g., continuous volumetric media such as clouds or fog).

Another line of work addresses the neural reconstruction of dynamic volumes by incorporating physical soft constraints derived from fluid dynamics [54], [55]. While these approaches model the transport of matter in participating media, our work instead focuses on the transport of light.

2.3 Generative Models of Optical Properties↩︎

Recent work has explored generative modeling of scene appearance and optical properties, typically under simplified image formation models. A first line of work focuses on surface-based representations and learns generative models of material appearance under simple shading assumptions, often limited to direct illumination [56]–[60]. These methods do not model volumetric light transport or participating media. A second line of work considers volumetric scene representations and learns generative models under emission–absorption rendering [61]–[65]. While some approaches target volumetric phenomena such as participating media, the representations remain appearance-driven and do not explicitly model physical optical parameters or enforce global illumination. To our knowledge, no prior work enables generative modeling of physical volumetric optical properties, such as absorption, scattering, and phase functions, under full global illumination.

3 Radiative Transfer Preliminaries↩︎

The full volumetric light transport in scenes with participating media is governed by the Radiative Transfer Equation (RTE) [1], [66], an integro-differential equation that describes how radiance evolves through the medium under absorption, emission, and scattering. We denote by \(\TextOrMath{L}{L}(\TextOrMath{\mathbf{x}}{\mathbf{x}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}) \in \TextOrMath{\mathbb{R}}{\mathbb{R}}^3\) the radiance at spatial location \(\TextOrMath{\mathbf{x}}{\mathbf{x}}\in \TextOrMath{\mathbb{R}}{\mathbb{R}}^3\) traveling in direction \(\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}\in \TextOrMath{\mathbb{S}}{\mathbb{S}}^2\), which characterizes the spatio-angular distribution of light throughout the volume. We refer to this function as light field, avoiding the term radiance field to prevent confusion with emission–absorption-based scene representations [8], in which “radiance” denotes a local appearance term that is conceptually closer to than to the steady-state radiance satisfying the RTE. We model the spectral dependence of the light field and of all other quantities using three (RGB) channels.

Under the assumption of steady state, the RTE takes the form \[\begin{align} \overbrace{ (\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}\cdot \nabla)\,\TextOrMath{L}{L}(\TextOrMath{\mathbf{x}}{\mathbf{x}}, \TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}) }^{\textrm{Transport}} &= - \overbrace{ \TextOrMath{\sigma_\mathrm{t}}{\sigma_\mathrm{t}}(\TextOrMath{\mathbf{x}}{\mathbf{x}})\,\TextOrMath{L}{L}(\TextOrMath{\mathbf{x}}{\mathbf{x}}, \TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}) }^{\textrm{Extinction}} + \overbrace{ \TextOrMath{\sigma_\mathrm{a}}{\sigma_\mathrm{a}}(\TextOrMath{\mathbf{x}}{\mathbf{x}})\,\TextOrMath{L_\mathrm{e}}{L_\mathrm{e}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}, \TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}) }^{\textrm{Emission}} \\ &\quad + \underbrace{ \TextOrMath{\sigma_\mathrm{s}}{\sigma_\mathrm{s}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}) \int_{\TextOrMath{\mathbb{S}}{\mathbb{S}}^2} \TextOrMath{f_\textrm{p}}{f_\textrm{p}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}, \TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}, \TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}')\,\TextOrMath{L}{L}(\TextOrMath{\mathbf{x}}{\mathbf{x}}, \TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}')\,\mathrm{d}\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}'} _{\textrm{In-scattering}}\!\!. \end{align} \label{eq:rte}\tag{1}\] The medium is parameterized by spatially-varying RGB absorption and scattering coefficients \(\TextOrMath{\sigma_\mathrm{a}}{\sigma_\mathrm{a}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}), \TextOrMath{\sigma_\mathrm{s}}{\sigma_\mathrm{s}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}) \in \TextOrMath{\mathbb{R}}{\mathbb{R}}^3\), with the extinction coefficient given by \(\TextOrMath{\sigma_\mathrm{t}}{\sigma_\mathrm{t}}(\TextOrMath{\mathbf{x}}{\mathbf{x}})=\TextOrMath{\sigma_\mathrm{a}}{\sigma_\mathrm{a}}(\TextOrMath{\mathbf{x}}{\mathbf{x}})+\TextOrMath{\sigma_\mathrm{s}}{\sigma_\mathrm{s}}(\TextOrMath{\mathbf{x}}{\mathbf{x}})\). Volumetric emission is described by the emitted radiance \(\TextOrMath{L_\mathrm{e}}{L_\mathrm{e}}(\TextOrMath{\mathbf{x}}{\mathbf{x}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}) \in \TextOrMath{\mathbb{R}}{\mathbb{R}}^3\). Scattering is governed by the spatially varying phase function \(\TextOrMath{f_\textrm{p}}{f_\textrm{p}}(\TextOrMath{\mathbf{x}}{\mathbf{x}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}') \in \TextOrMath{\mathbb{R}}{\mathbb{R}}_{\geq 0}\), which specifies the normalized angular redistribution of light scattered from incident direction \(\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}'\) into outgoing direction . The left-hand side of Eq. 1 describes infinitesimal directional transport of radiance through the medium, while the right-hand side accounts for local loss due to extinction and local gain due to emission and in-scattering. Although the scattering term integrates radiance over directions, all interactions in the RTE are local in space. Through this angular coupling induced by scattering, radiance at any location and direction arises from a dense superposition of indirect light transport contributions spanning the entire volume, corresponding to full global illumination.

An explicit equation for (,) can be obtained by integrating both sides of Eq. 1 . This integral form has a long history in radiative transfer physics [66], and is known in computer graphics as the Volume Rendering Equation (VRE) [1], [35]: \[\begin{align} \TextOrMath{L}{L}(\TextOrMath{\mathbf{x}}{\mathbf{x}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}) &= \int_{0}^{t_b} \overbrace{ \exp\!\left( - \int_{0}^{t} \TextOrMath{\sigma_\mathrm{t}}{\sigma_\mathrm{t}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}_s)\,\mathrm{d}s \right) }^{\textrm{Transmittance}} \Bigg[ \overbrace{ \TextOrMath{\sigma_\mathrm{a}}{\sigma_\mathrm{a}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}_t)\,\TextOrMath{L_\mathrm{e}}{L_\mathrm{e}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}_t,\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}) \vphantom{ \exp\!\left( - \int_{0}^{t} \TextOrMath{\sigma_\mathrm{t}}{\sigma_\mathrm{t}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}_s)\,\mathrm{d}s \right) } }^{\textrm{Emission}} \\ &\qquad\quad + \underbrace{ \TextOrMath{\sigma_\mathrm{s}}{\sigma_\mathrm{s}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}_t) \int_{\TextOrMath{\mathbb{S}}{\mathbb{S}}^2} \TextOrMath{f_\textrm{p}}{f_\textrm{p}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}_t,\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}')\, \TextOrMath{L}{L}(\TextOrMath{\mathbf{x}}{\mathbf{x}}_t,\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}')\,\mathrm{d}\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}' }_{\textrm{In-scattering}} \Bigg] \,\mathrm{d}t. \end{align} \label{eq:vre}\tag{2}\] Here, points along the ray \(\TextOrMath{\mathbf{x}}{\mathbf{x}}+ \tau\,\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}\) are denoted by \(\TextOrMath{\mathbf{x}}{\mathbf{x}}_\tau\) for \(\tau \ge 0\), and \(t_b\) is the distance to the volume boundary. The volume rendering equation expresses radiance along a ray as an integral over local source terms, with the local extinction coefficient accumulated into a transmittance factor that encodes cumulative attenuation.

Both the differential RTE and the integral VRE are inherently recursive due to scattering-induced self-coupling of radiance; however, this recursion is expressed implicitly in the RTE through spatially local gain–loss constraints, whereas in the VRE it appears explicitly through nested, spatially nonlocal ray integrals.

4 Method↩︎

Input to our method is a set of posed RGB images \(\{\TextOrMath{I}{I}_\TextOrMath{k}{k}\}\) providing multi-view observations of a scene containing heterogeneous participating media. We assume a known and fixed environment illumination. Our goal is to reconstruct the optical properties of the participating media, namely spatially varying, color-resolved absorption and scattering coefficients \(\TextOrMath{\sigma_\mathrm{a}}{\sigma_\mathrm{a}}(\TextOrMath{\mathbf{x}}{\mathbf{x}})\) and \(\TextOrMath{\sigma_\mathrm{s}}{\sigma_\mathrm{s}}(\TextOrMath{\mathbf{x}}{\mathbf{x}})\), as well as spatially varying phase functions \(\TextOrMath{f_\textrm{p}}{f_\textrm{p}}(\TextOrMath{\mathbf{x}}{\mathbf{x}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}')\), under physically complete volumetric light transport.

While our formulation focuses on reconstructing the optical properties of a single scene, it naturally generalizes to a generative setting when observations from multiple scenes are available. Such a generative model can be used to populate virtual scenes with a diverse set of assets, such as clouds. It could also be employed in an inverse setting, where a latent code is optimized to explain an observation (a photo of a participating medium), thereby providing a strong prior on the distribution of plausible 3D volumes. In both settings – reconstruction and generation – the recovered physical scene representations enable physically meaningful analysis, novel-view synthesis, relighting under novel illumination, and other operations that rely on interpretable optical properties.

Physically accurate light transport in participating media involves complex, highly coupled interactions (Sec. [sec:background]) that are typically addressed using specialized simulation techniques. Here, we instead explore how scene properties can be estimated under these interactions through general-purpose neural optimization, with only minimal reliance on explicit transport simulation.

We first introduce our general framework for volumetric inverse rendering (Sec. [sec:rep+opt]). We then describe how this approach can be extended to the generative setting (Sec. [sec:generative]), before providing implementation details (Sec. [sec:implementation]).

4.1 Neural Representation and Optimization↩︎

Central to our approach is the representation of the light field and the optical properties of the participating medium as neural fields (Fig. 2, left). Specifically, we employ a neural field \(\TextOrMath{\TextOrMath{L}{L}_\TextOrMath{\theta}{\theta}}{\TextOrMath{L}{L}_\TextOrMath{\theta}{\theta}}: (\TextOrMath{\mathbf{x}}{\mathbf{x}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}) \rightarrow \TextOrMath{\mathbb{R}}{\mathbb{R}}^3\), which encodes the spatio-angular light field of the entire scene [67], [68], and a neural field \(\TextOrMath{M_\TextOrMath{\phi}{\phi}}{M_\TextOrMath{\phi}{\phi}}: \TextOrMath{\mathbf{x}}{\mathbf{x}}\rightarrow \left( \TextOrMath{\sigma_\mathrm{a}}{\sigma_\mathrm{a}}, \TextOrMath{\sigma_\mathrm{s}}{\sigma_\mathrm{s}}, \TextOrMath{g}{g}\right)\), which encodes the participating medium through spatially varying absorption and scattering coefficients and a scalar parameter that defines a spatially varying Henyey–Greenstein [69] phase function . The subscripts and denote the optimizable parameters of the respective neural fields. To keep the notation concise, we omit explicit parameter subscripts for the medium properties in the remainder of this section and implicitly treat , , and as outputs of .

Disentanglement of the light field and medium properties, as well as the incorporation of data, is achieved through physics-informed optimization objectives combined with data constraints (Fig. 2, right), rather than through specialized architectural mechanisms for separating illumination and material properties. Accordingly, the remainder of this section details the individual constraints that form the core of our approach. While our ultimate objective is the recovery of the output of , most optimization objectives are imposed on to facilitate disentanglement of illumination and medium.

The central objective enforcing physically consistent global light transport is the formulation of Eq. 1 as a residual:

Here, we omit the emission term, since we consider only environment illumination, and \(\TextOrMath{\sigma_\mathrm{a}}{\sigma_\mathrm{a}}(\TextOrMath{\mathbf{x}}{\mathbf{x}})\), \(\TextOrMath{\sigma_\mathrm{s}}{\sigma_\mathrm{s}}(\TextOrMath{\mathbf{x}}{\mathbf{x}})\), and \(\TextOrMath{f_\textrm{p}}{f_\textrm{p}}(\TextOrMath{\mathbf{x}}{\mathbf{x}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}')\) are obtained from the neural medium field . Note that this formulation enforces the radiative transfer constraints at continuously sampled collocation points \((\TextOrMath{\mathbf{x}}{\mathbf{x}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}})\) within the spatio-angular domain (Fig. 1d) and is purely local in space. The first term in the expectation is evaluated using automatic differentiation. We approximate the in-scattering integral using Monte Carlo sampling over the angular domain, evaluated locally at each collocation point’s spatial location \(\TextOrMath{\mathbf{x}}{\mathbf{x}}\), without tracing rays or explicitly tracking light through the volume.

To ensure a unique solution of the radiative transfer equation, the light field is additionally constrained by inflow boundary conditions at the volume boundary. Let \(\TextOrMath{B}{B}\subset \mathbb{R}^3\) denote the spatial domain of the scene, with boundary and outward unit normal \(\TextOrMath{\mathbf{n}_\TextOrMath{B}{B}}{\mathbf{n}_\TextOrMath{B}{B}}(\TextOrMath{\mathbf{x}}{\mathbf{x}})\). The inflow boundary consists of all boundary position–direction pairs for which radiance enters the domain: \[\TextOrMath{\Gamma}{\Gamma} = \left\{ (\TextOrMath{\mathbf{x}}{\mathbf{x}},\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}})\;\middle|\; \TextOrMath{\mathbf{x}}{\mathbf{x}}\in \TextOrMath{\partial\TextOrMath{B}{B}}{\partial\TextOrMath{B}{B}},\; \TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}\in \TextOrMath{\mathbb{S}}{\mathbb{S}}^2,\; \TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}\cdot \TextOrMath{\mathbf{n}_\TextOrMath{B}{B}}{\mathbf{n}_\TextOrMath{B}{B}}(\TextOrMath{\mathbf{x}}{\mathbf{x}}) < 0 \right\}. \label{eq:inflowboundary}\tag{4}\] Boundary conditions are enforced through

where represents the incoming radiance at the boundary originating from environment illumination.

So far, we have only enforced data-agnostic physical and boundary constraints within the optimization. To incorporate the available image data, we introduce additional constraints on the light field derived from the observations. For an image \(\TextOrMath{I}{I}_{\TextOrMath{k}{k}}\), each pixel \(\TextOrMath{\mathbf{q}}{\mathbf{q}}\in \TextOrMath{\mathbb{R}}{\mathbb{R}}^2\) corresponds to a viewing direction \(\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}_{\TextOrMath{k}{k}}(\TextOrMath{\mathbf{q}}{\mathbf{q}}) \in \TextOrMath{\mathbb{S}}{\mathbb{S}}^2\) originating at the camera center \(\TextOrMath{\mathbf{o}}{\mathbf{o}}_{\TextOrMath{k}{k}} \in \TextOrMath{\mathbb{R}}{\mathbb{R}}^3\). Under known camera calibration, the observed pixel value \(\TextOrMath{I}{I}_{\TextOrMath{k}{k}}(\TextOrMath{\mathbf{q}}{\mathbf{q}}) \in \TextOrMath{\mathbb{R}}{\mathbb{R}}^3\) specifies the radiance arriving at the camera location from the corresponding viewing direction, i.e., traveling along direction \(-\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}_{\TextOrMath{k}{k}}(\TextOrMath{\mathbf{q}}{\mathbf{q}})\) (green arrows in Fig. 1d). We therefore impose this as a pointwise constraint through

where the expectation is taken over all observed views and pixels. This term constrains the light field only at the camera locations and for the set of observed incoming directions.

While the three objectives introduced so far are theoretically sufficient to guide the joint optimization of the light field and the material field , in practice we observe that the resulting solutions are overly smooth and lack high-frequency detail. This aligns with previous observations that neural network–based representations exhibit a bias toward low-frequency solutions even under direct supervision [17], and that related optimization pathologies have been reported for physics-informed formulations relying on differential constraints [16]. We alleviate this effect by introducing an additional objective based on the VRE (Eq. 2 ), which explicitly integrates radiance along primary viewing rays (ray in Fig. 1d) and thereby introduces nonlocal constraints that counteract the oversmoothing induced by purely local supervision:

where the locations \(\TextOrMath{\mathbf{x}}{\mathbf{x}}_\tau\) lie on the primary viewing ray \(\TextOrMath{\mathbf{o}}{\mathbf{o}}_{\TextOrMath{k}{k}} + \tau\,\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}_{\TextOrMath{k}{k}}(\TextOrMath{\mathbf{q}}{\mathbf{q}})\), while radiance is evaluated in the incoming direction \(-\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}_{\TextOrMath{k}{k}}(\TextOrMath{\mathbf{q}}{\mathbf{q}})\).

Crucially, the in-scattering radiance in Eq. 7 is obtained by querying the neural light field , which is jointly optimized under all objectives, with global light transport consistency enforced by the RTE residual Eq. 3 . As a result, although Eq. 7 aggregates radiance only along primary viewing rays, the radiance values integrated during volume rendering are not restricted to single-scattering or direct illumination, but instead evolve to encode the cumulative effect of arbitrarily many scattering interactions throughout the volume. Intuitively, the RTE governs the accumulation of indirect illumination throughout the volume, while the VRE term provides complementary non-local data coupling that anchors the directly observed component along each viewing direction.

In practice, the integrals in Eq. 7 are approximated using stochastic sampling. The line integral along each viewing ray is evaluated using jittered stratified sampling, following standard practice [8], [70]. The directional in-scattering integral appearing in Eq. 7 is handled identically to the corresponding term in the RTE residual Eq. 3 , using Monte Carlo sampling over the angular domain and evaluated locally at each spatial location.

Our complete optimization objective is given by \[\begin{align} \TextOrMath{\mathcal{L}}{\mathcal{L}}(\TextOrMath{\theta}{\theta},\TextOrMath{\phi}{\phi}) &= \TextOrMath{\mathcal{L}}{\mathcal{L}}_{\textrm{RTE}}(\TextOrMath{\theta}{\theta},\TextOrMath{\phi}{\phi}) + \TextOrMath{\lambda}{\lambda}_{\textrm{BC}} \TextOrMath{\mathcal{L}}{\mathcal{L}}_{\textrm{BC}}(\TextOrMath{\theta}{\theta}) \\ &\,+ \TextOrMath{\lambda}{\lambda}_{\textrm{obs}} \TextOrMath{\mathcal{L}}{\mathcal{L}}_{\textrm{obs}}(\TextOrMath{\theta}{\theta}) + \TextOrMath{\lambda}{\lambda}_{\textrm{VRE}} \TextOrMath{\mathcal{L}}{\mathcal{L}}_{\textrm{VRE}}(\TextOrMath{\theta}{\theta},\TextOrMath{\phi}{\phi}), \end{align} \label{eq:full-loss}\tag{8}\] where the scalar weights \(\TextOrMath{\lambda}{\lambda}_{\textrm{BC}}\), \(\TextOrMath{\lambda}{\lambda}_{\textrm{obs}}\), and \(\TextOrMath{\lambda}{\lambda}_{\textrm{VRE}}\) balance the relative contributions of the individual objectives.

4.2 Extension to Generative Modeling↩︎

Beyond inverse rendering of individual scenes, our approach naturally extends to learning a generative model of volumetric scene properties from a distribution of observations across multiple scenes. Specifically, we assume access to an extended set of posed observations \(\{\TextOrMath{I}{I}_{\TextOrMath{l}{l},\TextOrMath{k}{k}}\}\), where indexes scenes and indexes images within each scene, as before. In our proof-of-concept application, we assume the same illumination conditions across scenes.

We employ an auto-decoder framework [71], where each scene is associated with an optimizable latent code \(\TextOrMath{\mathbf{z}}{\mathbf{z}}_{\TextOrMath{l}{l}} \in \TextOrMath{\mathbb{R}}{\mathbb{R}}^{\TextOrMath{{d_{\mathbf{z}}}}{{d_{\mathbf{z}}}}}\). This latent code is provided as an additional input to both and , conditioning them on the corresponding scene. Training proceeds as in the scene-specific case described in Sec. [sec:rep+opt], with only two changes. First, the scene-specific latent codes \(\{\TextOrMath{\mathbf{z}}{\mathbf{z}}_{\TextOrMath{l}{l}}\}\) are jointly optimized with the neural field parameters and , and are additionally regularized with an \(L_2\) penalty [71]. Second, the observation-related objectives \(\TextOrMath{\mathcal{L}}{\mathcal{L}}_\mathrm{obs}\) (Eq. 6 ) and \(\TextOrMath{\mathcal{L}}{\mathcal{L}}_\mathrm{VRE}\) (Eq. 7 ) are extended to account for the distribution over scenes. Eq. 6 now reads \[\TextOrMath{\mathcal{L}}{\mathcal{L}}_{\mathrm{obs}}(\TextOrMath{\theta}{\theta}, \{\TextOrMath{\mathbf{z}}{\mathbf{z}}_{\TextOrMath{l}{l}}\}) = \mathbb{E}_{\TextOrMath{l}{l}, \TextOrMath{k}{k},\TextOrMath{\mathbf{q}}{\mathbf{q}}} \Big[ \big\| \TextOrMath{I}{I}_{\TextOrMath{l}{l}, \TextOrMath{k}{k}}(\TextOrMath{\mathbf{q}}{\mathbf{q}}) - \TextOrMath{\TextOrMath{L}{L}_\TextOrMath{\theta}{\theta}}{\TextOrMath{L}{L}_\TextOrMath{\theta}{\theta}}\big(\TextOrMath{\mathbf{o}}{\mathbf{o}}_{\TextOrMath{k}{k}},\,-\TextOrMath{\boldsymbol{\omega}}{\boldsymbol{\omega}}_{\TextOrMath{k}{k}}(\TextOrMath{\mathbf{q}}{\mathbf{q}}), \TextOrMath{\mathbf{z}}{\mathbf{z}}_{\TextOrMath{l}{l}}\big) \big\|_2^2 \Big], \label{eq:obs-loss-generative}\tag{9}\] and Eq. 7 is modified analogously. The resulting model defines a generative distribution over volumetric scenes, jointly capturing light fields and material property fields, which can be sampled by randomly drawing a latent code and evaluating the corresponding neural fields to obtain a complete scene.

It is worth noting that alternative inverse rendering formulations, including differentiable stochastic light transport, could in principle also be extended to learn distributions over medium properties, provided that the underlying scene representation supports such variability. However, doing so would require an independent light transport simulation for each sampled scene instance and for each training iteration, incurring the full cost of global illumination evaluation repeatedly during optimization. In contrast, our formulation jointly models a distribution over both medium properties and light fields, allowing the optimization to exploit shared structure across scenes and to amortize the resolution of global light transport through the learned light field representation. While the light field remains an auxiliary quantity, this joint modeling enables efficient learning of scene-level variability under physically complete volumetric illumination.

4.3 Implementation Details↩︎

Our implementation relies exclusively on general-purpose machine-learning components provided by PyTorch [72], without requiring any custom operations or rendering-specific kernels. We minimize Eq. 8 using the Adam [73] optimizer with default settings, and \(\TextOrMath{\lambda}{\lambda}_{\textrm{BC}} = \TextOrMath{\lambda}{\lambda}_{\textrm{obs}}=1\), and \(\TextOrMath{\lambda}{\lambda}_{\textrm{VRE}}=70\).

As observed in prior work on both differentiable sampling-based and learning-based inverse rendering [7], [36], scene initialization plays a crucial role in obtaining high-quality solutions. Accordingly, following common practice, we adopt a two-stage optimization procedure. In the first stage, we train an emission–absorption-based reconstruction model  [8] (Fig. 1a) using a neural hash grid architecture [74]. The color output of is discarded, as it represents a non-physical emission term. In the second stage, we jointly optimize the neural medium field and the light field , both parameterized by SIREN networks [75]. The medium field multiplicatively modulates the density output of , which provides a strong high-frequency initialization for the volumetric structure. We enforce positivity of and via an exponential mapping, and bound using a hyperbolic tangent. Experiments using an albedo–extinction parameterization yielded slightly inferior results. In the generative setting, latent codes are incorporated differently, with using multiplicative modulation of its hash grid features via an MLP-based encoding of prior to decoding, and and concatenating the latent code with the positional and directional inputs. The described architectural design choices yielded the best results in a reasonably comprehensive pilot search.

In all experiments, we use 10k collocation points to evaluate Eq. 3 , 32 directional samples to approximate the in-scattering integrals in Eq. 3 and Eq. 7 , 10k boundary samples for Eq. 5 , and 100 samples per viewing ray for evaluating Eq. 7 . All samples are drawn randomly from the corresponding continuous domains at each optimization iteration. To avoid optimization gradient bias arising from stochastic sampling of in-scattering directions, we draw independent samples for the forward and backward passes [76]. The runtime of our method is almost entirely attributable to evaluating and optimizing the four objective terms \(\TextOrMath{\mathcal{L}}{\mathcal{L}}_{\mathrm{RTE}}\) (38%), \(\TextOrMath{\mathcal{L}}{\mathcal{L}}_{\mathrm{BC}}\) (7%), \(\TextOrMath{\mathcal{L}}{\mathcal{L}}_{\mathrm{obs}}\) (2%), and \(\TextOrMath{\mathcal{L}}{\mathcal{L}}_{\mathrm{VRE}}\) (53%).

5 Evaluation↩︎

We first describe our evaluation setup, followed by an evaluation of our approach on volumetric scene reconstruction (Sec. [sec:eval-recon]) and generation (Sec. [sec:eval-gen]). We conclude with additional analysis (Sec. [sec:eval-analysis]).

Data Our work explores a novel formulation for modeling global illumination in volumetric inverse rendering, rather than aiming to develop a robust end-to-end pipeline for imperfect real-world data. To our knowledge, real multi-view datasets of volumetric scenes with known illumination do not yet exist. Accordingly, we evaluate our approach on synthetic scenes with known camera parameters and illumination conditions, with access to ground-truth volumetric properties for evaluation. To cover a broad range of volumetric configurations, we construct a custom dataset of 50 scenes. Each scene is assembled from a random set of volumetric elements drawn from the Disney Clouds dataset [77], which spans a continuum of structures from diffuse, low-frequency regions to sharply defined boundaries. For each selected element, we define RGB absorption and scattering coefficients, and , by scaling its base volumetric density with independently sampled random colors. This yields a high dynamic range of and (coefficients take values up to 5, corresponding to an expected 20 light–volume interactions per path) and produces diverse combinations of absorption and scattering, resulting in a wide range of optical behaviors spanning weakly absorbing, strongly scattering, and highly opaque elements. We additionally assign each element a Henyey-Greenstein phase function parameter sampled uniformly from \([-0.5, 0.5]\). A scene is formed by randomly selecting between one and four augmented elements, applying independent random spatial shifts, rotations, and isotropic scalings, and combining them into a single volume. Illumination is provided by ten environment maps from Polyhaven [78]. For each scene, we randomly choose two different environment maps and render the scene under both conditions to enable evaluation of relighting. Observations are obtained by volumetric path tracing from 500 camera viewpoints placed on a sphere around the volume and oriented toward its center; we reserve 50 views for testing. In addition to this multi-illumination dataset, we construct a single-illumination variant in which the same environment map is shared across all scenes for training and evaluating the generative setting. Finally, to enable direct comparison with baseline implementations that do not support anisotropic scattering, we also provide a variant with isotropic scattering in which all phase functions are fixed to \(\TextOrMath{g}{g}=0\). All datasets and the scripts used to generate them will be released upon publication.

Baselines We consider two baseline approaches. First, we compare against Differential Ratio Tracking (DRT) [7], which accounts for global illumination by differentiating through stochastic light path sampling, enabling the optimization of a 3D grid of optical parameters. Like our approach, the method relies on an emission–absorption initialization. To ensure a meaningful comparison, we initialize both methods using our emission–absorption solution, which we found yields improved performance for DRT compared to its default setup. For quality evaluation, we perform an equal-time comparison by running DRT beyond its recommended iteration count until its total runtime matches that of our method. The resulting slightly improved DRT results are those reported. For runtime comparisons, however, we use DRT’s recommended hyperparameters. DRT’s publicly available implementation supports only isotropic phase functions for inverse rendering.

Our second baseline is an approach based on a learned representation. Since we are not aware of a learning-based method that publicly provides a complete implementation for volumetric inverse rendering under global illumination, we compare against TensorIR [49] under fixed illumination as a representative baseline. TensorIR approximates second-order light transport using ray marching through an emissive volume and is not designed to recover physical volumetric optical properties.

5.1 Reconstruction↩︎

We evaluate all methods on the large-scale dataset with isotropic scattering using two evaluation protocols. First, we compute the mean squared error (MSE) of the reconstructed absorption, scattering, and extinction coefficients, as well as the phase function parameters, by densely sampling the corresponding quantities throughout the volume; the error of is weighted by the local scattering coefficient prior to aggregation. Second, we evaluate novel-view synthesis on test views under a novel illumination condition. Images for Differential Ratio Tracking and our approach are synthesized by re-rendering the reconstructed volumes using volumetric path tracing. We evaluate image quality using PSNR, SSIM [79], and LPIPS [80] on logarithmically tonemapped images to account for high dynamic range.

Quantitative results are listed in Tab. ¿tbl:tab:quant-comparison?, left, while qualitative comparisons are shown in Fig. 3. We observe that our method achieves significantly more accurate recovery of the scattering and extinction coefficients than differentiable stochastic path sampling. At the same time, image-level novel-view synthesis quality is highest for the stochastic path-sampling approach, reflecting the benefit of explicit Monte Carlo simulation for optimizing pixel-level appearance. The learning-based baseline did not produce competitive results on our challenging dataset, despite extensive efforts to ensure a fair evaluation. Specifically, we tested configurations using LDR input observations, fixed the environment illumination to the ground truth, and explored a range of training hyperparameters. As we were able to reproduce the authors’ reported results on the datasets used in their evaluation, we hypothesize that the method may not be well suited to the intricate volumetric structures present in our dataset.

A quantitative evaluation of reconstruction quality for our method with anisotropic scattering is provided in Tab. ¿tbl:tab:quant-anisotropic?. We observe that incorporating anisotropic scattering as an additional degree of freedom leads to only a moderate reduction in reconstruction accuracy compared to the isotropic case, while none of the baseline implementations support recovery of anisotropic scattering coefficients. Fig. [fig:teaser] further demonstrates the accuracy of the recovered scene properties for an anisotropic medium under complex relighting, including strong local light sources.

5.2 Generation↩︎

To evaluate the generative capabilities of our approach, we train a generative model using our large-scale dataset with a single shared illumination condition and a latent dimensionality of \(\TextOrMath{{d_{\mathbf{z}}}}{{d_{\mathbf{z}}}}= 4\). We then synthesize 30 scenes by sampling latent codes and render the corresponding scenes under the training illumination. To assess the faithfulness of the learned distribution over optical properties, we additionally render the same generated scenes under randomly sampled novel illuminations. In both cases, the distribution of the resulting renderings is compared to that of the corresponding real images – either from the training set or from the multi-illumination dataset – using the FID score [81]. For both evaluations, we extract random crops from each image to improve the robustness of the distributional comparison. The corresponding scores are reported in Tab. ¿tbl:tab:quant-comparison?, right, while qualitative samples are shown in Fig. 4. Our results show that the model learns a generative distribution over physical volumetric optical properties – rather than appearance – yielding scene samples that support consistent global illumination and relighting, a capability that has not been demonstrated by prior volumetric generative approaches.

5.3 Ablations↩︎

Our formulation is intentionally minimal, relying on a small set of physically motivated objectives (Sec. [sec:abl-objectives]) and two generic neural fields for scene representation (Sec. [sec:abl-rep]), with solutions obtained via optimization (Sec. [sec:abl-optimization]).

5.3.1 Objectives↩︎

Most of our optimization objectives are essential: removing any of the RTE residual, boundary conditions, or observation constraints renders the problem underconstrained and prevents convergence to a meaningful solution. Specifically, this leads to physically inconsistent solutions, neglect of illumination, and disregard of observations, respectively. We therefore focus our ablation study on the only component of our formulation that appears optional, namely the VRE-based objective (Eq. 7 ). As shown in Fig. 5, omitting the VRE term leads to overly smooth reconstructions that lack high-frequency detail. This behavior is consistent with the known low-frequency bias of physics-informed neural networks relying solely on local differential constraints. Importantly, this effect is largely orthogonal to architectural choices such as positional encoding [8], [82], which only provide the capacity to represent high-frequency functions but still require sufficiently direct supervision for this capacity to be utilized. Incorporating the nonlocal VRE objective effectively provides such supervision and yields substantially sharper and more faithful reconstructions.

5.3.2 Representation↩︎

We use neural fields to represent both the light field and the volumetric properties. For reconstruction, the former is inherently a 5D function and thus naturally represented using a neural field, whereas the latter is 3D and therefore, in principle, amenable to straightforward discretization. However, as shown in Tab. ¿tbl:tab:ablation?, a grid-based representation yields inferior results compared to our neural representation in this setting. Moreover, grid-based representations are fundamentally incompatible with our generative modeling formulation.

5.3.3 Optimization↩︎

Our optimization proceeds in two stages. First, an emission–absorption model is obtained, whose density is subsequently modulated by a function optimized in the second stage to yield the final medium properties. As shown in Tab. ¿tbl:tab:ablation?, a direct, single-stage optimization yields inferior results compared to our two-stage approach.

6 Conclusion↩︎

We have demonstrated that, for global-illumination-aware inverse rendering, a global-illumination renderer is not essential. By representing both the optical properties of participating media and the full light field as neural fields, and by constraining them through predominantly local formulations of radiative transfer, our optimization-based approach enables the reconstruction of high-quality, physically meaningful volumetric parameters. Beyond inverse rendering, the proposed formulation naturally supports generative modeling of participating media with physical optical properties, learned directly from image observations.

Because both the scene representation and the optimization objective are expressed in terms of principled, physics-based constraints, the framework is agnostic to the dimensionality of the underlying signal and compositional with respect to constraints. This opens promising directions for extending the approach to spectral and transient light transport, as well as dynamic scenes, while preserving physical consistency. Additionally, while our experiments focus on environment illumination, the proposed framework can be extended to other lighting representations, such as local light sources, by incorporating the emission term into the radiative transfer equation. Further, exploring the application of Eq. 7 to estimate incident illumination within the volume could lead to hybrid approaches that combine our neural formulation with explicit simulation. We believe this perspective provides a flexible foundation for future work on volumetric scene understanding and synthesis.

Acknowledgments↩︎

We thank Felix Mujkanovic and Stephanie Wenxin Liu for their assistance with the evaluation of this work.
Open Access funding enabled and organized by Projekt DEAL.

References↩︎