Neural Microfacet Fields for Inverse Rendering

The core idea in emission-absorption volume rendering is that light accumulates along rays, with “particles” along the ray both emitting and absorbing light. The color measured by a camera pixel corresponding to a ray with origin $\mathbf{p}_{c}$ and direction $-\hat{\boldsymbol{\omega}}_{o}$ is:

	$\displaystyle L(\mathbf{p}_{c},\hat{\boldsymbol{\omega}}_{o})$	$\displaystyle=\int_{0}^{\infty}T(t)\sigma(\vec{\mathbf{r}}(t))L_{o}(\vec{\mathbf{r}}(t),\hat{\boldsymbol{\omega}}_{o})dt,$		(1)
	$\displaystyle\text{where}\;\;T(t)$	$\displaystyle=\exp\left(-\int_{0}^{t}\sigma(\vec{\mathbf{r}}(t^{\prime}))dt^{\prime}\right),$		(2)

where $\vec{\mathbf{r}}(t)=\mathbf{p}_{c}-t\hat{\boldsymbol{\omega}}_{o}$ is a camera ray, $\sigma(\mathbf{p})$ is the density at point $\mathbf{p}$ in the volume, $T$ denotes transmittance along the ray, and $L_{o}$ is the outgoing radiance. This formula is often approximated numerically using quadrature, following [26]:

	$\displaystyle L(\mathbf{p}_{c},\hat{\boldsymbol{\omega}}_{o})$	$\displaystyle\approx\sum_{j=0}^{N-1}w_{j}L_{o}(\vec{\mathbf{r}}(t_{j}),\hat{\boldsymbol{\omega}}_{o}),$		(3)
	$\displaystyle\text{where}\;\;w_{j}$	$\displaystyle=T_{j}\big{(}1-\exp(-\sigma(\vec{\mathbf{r}}(t_{j}))(t_{j+1}-t_{j}))\big{)},$		(4)
	$\displaystyle\text{and}\;\;T_{j}$	$\displaystyle=\exp\left(-\sum_{k=0}^{j-1}\sigma(\vec{\mathbf{r}}(t_{k}))(t_{k+1}-t_{k})\right).$		(5)

In this volume rendering paradigm, multiple 3D points can contribute to the color of a ray, with nearer and denser points contributing most.

In surface rendering, and assuming fully-opaque surfaces, the color of a ray is determined solely by the light reflected by the first surface it encounters. Consider that the ray $\vec{\mathbf{r}}$ from camera position $\mathbf{p}_{c}$ in direction $-\hat{\boldsymbol{\omega}}_{o}$ intersects its first surface at a 3D position $\mathbf{p}$. The ray color is then:

$$L(\mathbf{p}_{c},\hat{\boldsymbol{\omega}}_{o})=\int_{\mathbb{S}^{2}}f(\mathbf{p},\hat{\boldsymbol{\omega}}_{o},\hat{\boldsymbol{\omega}}_{i})L_{i}(\mathbf{p},\hat{\boldsymbol{\omega}}_{i})(\hat{\mathbf{n}}(\mathbf{p})\cdot\hat{\boldsymbol{\omega}}_{i})^{+}d\hat{\boldsymbol{\omega}}_{i},$$

(6)

where $\hat{\boldsymbol{\omega}}_{i}$ is the direction of incident light, $\hat{\mathbf{n}}(\mathbf{p})$ is the surface normal at $\mathbf{p}$, $f$ is the BRDF describing the material of the surface at $\mathbf{p}$, $L_{i}$ is the incident radiance, and $(\hat{\mathbf{n}}(\mathbf{p})\cdot\hat{\boldsymbol{\omega}}_{i})^{+}$ is a truncated cosine lobe (i.e. its negative values are clipped to zero) facing outward from the surface. Note that this equation is recursive: $L_{i}$ inside the integral may be the outgoing radiance $L_{o}$ coming from a different scene point. The integral in Equation 6 is also typically approximated by discrete (and often random) sampling, and is the subject of a rich body of work [11, 40].

The key to our method is a novel combination of the volume rendering and surface rendering paradigms: we model a density field as in volume rendering, and we model outgoing radiance at every point in space using surface-based light transport (approximated using Monte Carlo ray sampling). Volume rendering with a density field lends itself well to optimization: initializing geometry as a semi-transparent cloud creates useful gradients (see Figure 3), and allows for changes in geometry and topology. Using surface-based rendering allows modeling the interaction of light and materials, and enables recovering these materials.

We combine these paradigms by modeling a microfacet field, in which each point in space is endowed with a volume density and a local micro-surface. Light accumulates along rays according to the volume rendering integral of Equation 1, but the outgoing light of each 3D point is determined by surface rendering as in Equation 6, using rays sampled according to its local micro-surface. This combination of volume-based and surface-based representation and rendering, shown in Figure 2, enables us to optimize through a severely underconstrained inverse problem, recovering geometry, materials, and illumination simultaneously.

We represent geometry using a low-rank tensor data structure based on TensoRF [8], with small modifications described in Appendix C. Our model stores both density $\sigma$ and a spatially-varying feature that is decoded into the material’s BRDF at every point in space. We initialize our model at low resolution and gradually upsample it during optimization (see Appendix C for details).

Similar to prior work [37, 41], we use the negative normalized gradient of the density field as a field of “volumetric normals.” However, like [20], we found that numerically computing spatial gradients of the density field using finite differences rather than using analytic gradients leads to normal vectors that we can use directly, without using features predicted by a separate MLP. Additionally, these numerical gradients can be efficiently computed using 2D and 1D convolution using TensorRF’s low-rank density decomposition (see Appendix C). These accurate normals are then used for rendering the appearance at a volumetric microfacet, as will be discussed in Sections 4.3 and 4.5.

Our volumetric normals are regularized using the orientation loss $\mathcal{R}_{o}$ introduced by Ref-NeRF [41]:

$$\mathcal{R}_{o}=\sum_{j}w_{j}\max(0,-\hat{\mathbf{n}}(\mathbf{p}_{j})\cdot\hat{\boldsymbol{\omega}}_{o})^{2},$$

(7)

where $\hat{\boldsymbol{\omega}}_{o}$ is the view direction facing towards the camera, and $\hat{\mathbf{n}}(\mathbf{p}_{j})$ is the normal vector at the $j$th point along the ray. The orientation loss $\mathcal{R}_{o}$ penalizes normals that face away from the camera yet contribute to the color of the ray (as quantified by weights $w_{j}$).

Because our volumetric normals are derived from the density field, this regularizer has a direct effect on the reconstructed geometry: it decreases the weight of backwards-facing normals by decreasing their density or increasing the density between them and the cameras, thereby promoting hard surfaces and improving reconstruction. Note that unlike Ref-NeRF, we do not use “predicted normals” for surface rendering, as our Gaussian-smoothed derivative filter achieves similar effect.

We write our spatially varying BRDF model $f$ as a combination of diffuse and specular components:

$\displaystyle f(\hat{\boldsymbol{\omega}}_{o},\hat{\boldsymbol{\omega}}_{i})=\frac{\rho}{\pi}(1-F_{r}(\hat{\boldsymbol{h}}))+F_{r}(\hat{\boldsymbol{h}})f_{s}$

$\displaystyle(\hat{\boldsymbol{\omega}}_{o},\hat{\boldsymbol{\omega}}_{i}),$

(8)

where $\rho$ is the RGB albedo, $\hat{\boldsymbol{h}}=\frac{\hat{\boldsymbol{\omega}}_{o}+\hat{\boldsymbol{\omega}}_{i}}{\|\hat{\boldsymbol{\omega}}_{o}+\hat{\boldsymbol{\omega}}_{i}\|}$ is the half vector, $F_{r}(\hat{\boldsymbol{h}})$ is the Fresnel term, $f_{s}(\hat{\boldsymbol{\omega}}_{o},\hat{\boldsymbol{\omega}}_{i})$ is the specular component of the BRDF for outgoing view direction $\hat{\boldsymbol{\omega}}_{o}$ and incident light direction $\hat{\boldsymbol{\omega}}_{i}$. The spatial dependence of these terms on the point $\mathbf{p}$ is omitted for brevity. We use the Schlick approximation [35] for the Fresnel term:

$$F_{r}(\hat{\boldsymbol{h}})=F_{0}(\mathbf{p})+(1-F_{0}(\mathbf{p}))(1-\hat{\boldsymbol{h}}\cdot\hat{\boldsymbol{\omega}}_{o})^{5},$$

(9)

where $F_{0}(\mathbf{p})\in[0,1]^{3}$ is the spatially varying reflectance at the normal incidence at the point $\mathbf{p}$, and we base our specular BRDF on the Cook-Torrance BRDF [38], using:

$$f_{s}(\hat{\boldsymbol{\omega}}_{o},\hat{\boldsymbol{\omega}}_{i})=\frac{D(\hat{\boldsymbol{h}};\alpha,\hat{\mathbf{n}},\hat{\boldsymbol{\omega}}_{o})G_{1}(\hat{\boldsymbol{\omega}}_{o},\hat{\boldsymbol{h}})g(\hat{\boldsymbol{\omega}}_{i},\hat{\boldsymbol{\omega}}_{o})}{4\left(\hat{\mathbf{n}}\cdot\hat{\boldsymbol{\omega}}_{o}\right)\left(\hat{\mathbf{n}}\cdot\hat{\boldsymbol{\omega}}_{i}\right)},$$

(10)

where $D$ is a Trowbridge-Reitz distribution [39] (popularized by the GGX BRDF model [42]), $G_{1}$ is the Smith shadow masking function for the Trowbridge-Reitz distribution, and $g$ is a shallow multilayer perceptron (MLP) with a sigmoid nonlinearity at its output. The distribution $D$ models the roughness $\alpha$ of the material, and it is used for importance sampling, as described in Section 4.5). The MLP $g$ captures other material properties not included in its explicit components.

The parameters for each of the diffuse and specular BRDF components are stored as features $\boldsymbol{x}$ in the TensoRF representation (alongside density $\sigma$), allowing them to vary in space. We compute the roughness $\alpha$, albedo $\rho$ and reflectance at the normal incidence $F_{0}$ by applying a single linear layer with sigmoid activation to the spatially-localized features $\boldsymbol{x}$. Details about the architecture of the MLP $g$ and its input encoding can be found in Appendix B.

We can approximately evaluate the rendering equation integral in Equation 6 more efficiently by assuming all microfacets at a point $\mathbf{p}$ have the same irradiance $E(\mathbf{p})$:

	$\displaystyle\int_{\mathbb{S}^{2}}$	$\displaystyle f(\mathbf{p},\hat{\boldsymbol{\omega}}_{o},\hat{\boldsymbol{\omega}}_{i})L_{i}(\mathbf{p},\hat{\boldsymbol{\omega}}_{i})(\hat{\mathbf{n}}(\mathbf{p})\cdot\hat{\boldsymbol{\omega}}_{i})^{+}d\hat{\boldsymbol{\omega}}_{i}$		(11)
	$\displaystyle\approx\int_{\mathbb{S}^{2}}$	$\displaystyle F_{r}(\hat{\boldsymbol{h}})f_{s}(\mathbf{p},\hat{\boldsymbol{\omega}}_{o},\hat{\boldsymbol{\omega}}_{i})L_{i}(\mathbf{p},\hat{\boldsymbol{\omega}}_{i})(\hat{\mathbf{n}}(\mathbf{p})\cdot\hat{\boldsymbol{\omega}}_{i})^{+}d\hat{\boldsymbol{\omega}}_{i}$
		$\displaystyle+\frac{\rho(\mathbf{p})}{\pi}E(\mathbf{p})\int_{\mathbb{S}^{2}}(1-F_{r}(\hat{\boldsymbol{h}}))d\hat{\boldsymbol{\omega}}_{i}.$		(12)

The irradiance $E(\mathbf{p})$, defined as:

$$E(\mathbf{p})=\int_{\mathbb{S}^{2}}L_{i}(\mathbf{p},\hat{\boldsymbol{\omega}}_{i})(\hat{\mathbf{n}}(\mathbf{p})\cdot\hat{\boldsymbol{\omega}}_{i})^{+}d\hat{\boldsymbol{\omega}}_{i},$$

(13)

can then be easily evaluated using an irradiance environment map approximated by low-degree spherical harmonics, as done by Ramamoorthi and Hanrahan [32]. At every optimization step, we obtain the current irradiance environment map by integrating the environment map with spherical harmonic functions up to degree $2$, and combining the result with the coefficients of a clamped cosine lobe pointing in direction $\hat{\mathbf{n}}(\mathbf{p})$ to obtain the irradiance $E(\mathbf{p})$ [32]. Equation 4.3 can then be importance sampled according to $D$ and integrated using Monte Carlo sampling of incoming light.

We sample half vectors $\hat{\boldsymbol{h}}$ from the distribution of visible normals $D_{\hat{\boldsymbol{\omega}}_{o}}$ [18], which is defined as:

$$D_{\hat{\boldsymbol{\omega}}_{o}}(\hat{\boldsymbol{h}})=\frac{G_{1}(\hat{\boldsymbol{\omega}}_{o},\hat{\boldsymbol{h}})|\hat{\boldsymbol{\omega}}_{o}\cdot\hat{\boldsymbol{h}}|D(\hat{\boldsymbol{h}})}{|\hat{\boldsymbol{\omega}}_{o}\cdot\hat{\mathbf{n}}|}.$$

(14)

However, to perform Monte Carlo integration of the rendering equation, we need to convert from half vector space to the space of incoming light, which requires multiplying by the determinant of the Jacobian of the reflection equation $\hat{\boldsymbol{\omega}}_{i}=2(\hat{\boldsymbol{\omega}}_{o}\cdot\hat{\boldsymbol{h}})\hat{\boldsymbol{h}}-\hat{\boldsymbol{\omega}}_{o}$, which is $4(\hat{\boldsymbol{\omega}}_{o}\cdot\hat{\boldsymbol{h}})$ [42]. Multiplying by the $\hat{\boldsymbol{\omega}}_{i}\cdot\hat{\mathbf{n}}$ term from Equation 6 as well as the Jacobian and Equation 14 results in the following Monte Carlo estimate:

	$\displaystyle L(\mathbf{p},\hat{\boldsymbol{\omega}}_{o})\approx\frac{1}{N}\sum_{n=1}^{N}$	$\displaystyle F_{r}(\hat{\mathbf{h}}^{n})g(\boldsymbol{x},\hat{\boldsymbol{\omega}}_{o},\hat{\boldsymbol{\omega}}_{i}^{n})L_{i}(\mathbf{p},\hat{\boldsymbol{\omega}}_{i}^{n})$
		$\displaystyle+(1-F_{r}(\hat{\mathbf{h}}^{n}))\frac{\rho(\mathbf{p})}{\pi}E(\mathbf{p}),$		(15)
	$\displaystyle\text{where }\hat{\mathbf{h}}^{n}\sim$	$\displaystyle D_{\hat{\boldsymbol{\omega}}_{o}}(\,\,\cdot\,\,;\alpha(\mathbf{p}),\hat{\mathbf{n}}(\mathbf{p}),\hat{\boldsymbol{\omega}}_{o}),$
	$\displaystyle\text{and }\hat{\boldsymbol{\omega}}_{i}^{n}=$	$\displaystyle 2(\hat{\boldsymbol{\omega}}_{o}\cdot\hat{\mathbf{h}}^{n})\hat{\mathbf{h}}^{n}-\hat{\boldsymbol{\omega}}_{o}.$

We model far field illumination using an environment map, represented using an equirectangular image with dimensions $H\times W\times 3$, with $H=512$ and $W=1024$. We map the optimizable parameters in our environment map parameterization into high dynamic range RGB values by applying an elementwise exponential function.

We use the term primary to denote a ray originating at the camera, and secondary to denote a ray bounced from a surface to evaluate its reflected light (whether that light arrives directly from the environment map, or from another scene point).

To minimize sampling noise, instead of using a single environment map element per secondary ray, we use the mean value over an axis-aligned rectangle in spherical coordinates, with the solid angle covered by the rectangle adjusted to match the sampling distribution $D$ at that point. Concretely, to query the environment map at a given incident light direction, we first compute its corresponding spherical coordinates $(\theta,\phi)$, where $\theta$ and $\phi$ are the polar and azimuthal angles respectively. We then compute the mean value of the environment map over a (spherical) rectangle centered at $(\theta,\phi)$, whose size $\Delta\theta\times\Delta\phi$ we constrain to have aspect ratio $\frac{\Delta\theta}{\Delta\phi}=\sin\theta$. To choose the solid angle of the rectangle, $\Delta\theta\cdot\Delta\phi$, we modify the method from [10], which is based on Nyquist’s sampling theorem (see Appendix A). We compute these mean values efficiently using integral images, also known as summed-area tables [13].

In this section, we describe how the model components representing geometry (Section 4.2), materials (Section 4.3), and illumination (Section 4.4) are combined to render the color of a pixel. For each primary ray, we choose a set of sample points following the rejection-sampling strategy of [21, 28] to prioritize points near object surfaces. We query our geometry representation [8] for the density $\sigma_{j}$ at each point, and use Equation 4 to estimate the contribution weight $w_{j}$ of each point to the final ray color.

To compute the color for each of these points, we compute the irradiance from the environment map and apply Equation 4.3 to obtain $L(\mathbf{p}_{j},\hat{\boldsymbol{\omega}}_{o})$, as described in Section 4.3.

For each primary ray sample with weight $w_{j}$, we allocate $N=\left\lfloor w_{j}M\right\rfloor$ secondary rays, where $M=128$ is an upper bound on the total number of secondary rays for each primary ray (since $\sum_{j}w_{j}\leq 1$). The secondary rays are sampled according to the Trowbridge-Reitz distribution using the normal vector $\hat{\mathbf{n}}(\mathbf{p}_{j})$ and roughness value $\alpha_{j}$ at the current sample.

When computing the incoming light $L_{i}(\mathbf{p},\hat{\boldsymbol{\omega}}_{i}^{n})$, we save memory by randomly selecting a fixed number $R$ of secondary rays to interreflect through the scene while others index straight into the environment map. We importance sample these $R$ secondary rays according to the largest channel of the weighted RGB multiplier $w_{j}\cdot g(\boldsymbol{x}_{j},\hat{\boldsymbol{\omega}}_{i},\hat{\boldsymbol{\omega}}_{o})$. In practice, we find that adding a small amount of random noise to the weighted RGB multiplier before choosing the $R$ largest values improves performance by slightly increasing the variation of the selected rays.

The remaining $N-R$ secondary rays with lower contribution are rendered more cheaply by evaluating the environment map directly rather than considering further interactions with the scene, as described in Section 4.4.

This combined sample color value $L(\mathbf{p}_{j},\hat{\boldsymbol{\omega}}_{o})$ is then weighted by $w_{j}$, and the resulting colors are summed along the primary ray samples to produce pixel values. Finally, the resulting pixel values are tonemapped to sRGB color and clipped to $[0,1]$.