1. Introduction

Physics-based rendering (PBR) engines attempt to generate photorealistic images by modeling light-matter interaction in a physically plausible way. PBR has become the standard rendering method in animation, gaming, and computer graphic research. More recently, PBR engines have incorporated polarization into rendering algorithms [1,2]. The appearance of specular and diffuse highlights, as well as the shape of an object, is improved by material models based on polarization [3]. These models are polarized bidirectional reflection distribution functions (pBRDFs), which are Mueller- or Jones-matrix-valued functions that depend on material properties and acquisition geometry. Polarization-aware rendering engines traditionally utilize analytic pBRDF Mueller matrix models that are limited in the range of materials that they describe. The first and only publicly available polarized material measurement database for computer graphic applications is the KAIST pBRDF database [4]. The database contains Mueller matrices measured for 25 materials that cover a wide range of scattering geometries. A Mueller matrix is a $4\times 4$ real-valued matrix with 16 independent parameters [5]. Li et. al proposed a pBRDF model which uses a single depolarization parameter to quantify the relative weights between a completely depolarizing and a non-depolarizing Mueller matrix (e.g. Fresnel reflection or transmission) [6]. A Mueller matrix with a triply degenerate (TD) coherency eigenspectrum can be expressed using this single depolarization parameter pBRDF without significant loss of information. Instead of 16 parameters for a general Mueller matrix, the TD form has eight parameters: six for the normalized nondepolarizing Mueller matrix, a single depolarization parameter (e.g the largest eigenvalue), and the unpolarized reflectance. Interpolating these eight TD parameters to unmeasured scattering geometries maintains a physically-realizable Mueller matrix. Using the TD form to represent a Mueller measurement compresses the data by a factor of two, as well as provides a physical constraint for interpolation. Rendering computations from the TD model offer a probabilistic interpretation. The depolarization parameter can be used to compute the probability for a given light-matter interaction. Then, for each ray, a randomly generated selection is made between a non-depolarizing Mueller matrix or a completely depolarizing Mueller matrix. The rays are summed incoherently after the ray trace is completed. In this work, rendered polarimetric images using the TD model, and the TD model computed with polarimetric importance sampling, are compared to the original KAIST database.

The contributions of this work are:

2. Related work

Radiometric BRDF models are scalar-valued and when polarization is required these models are expanded to pBRDF models which are Mueller or Jones matrix valued. Computer vision, computer graphics, atmospheric science, and optical design utilize pBRDF models. Hyde et al. expanded the polarized Cook-Torrance model by including an ideal depolarizer weighted by the albedo of the material. This pBRDF intends to correctly model realistic light-matter interactions by incorporating shadowing/masking functions and a Lambertian model for diffuse reflection. Baek et al. included a polarized diffuse term that models partial polarization from subsurface scattering. Kondo et al. adopted this model, but included an unpolarized albedo component as an additional term in the pBRDF [7,8]. In Kondo et al., a coaxial illumination configuration restricted their evaluation to in-plane scattering. They also demonstrated surface normal retrieval, refractive index, and inverse rendering using their model [8]. The pBRDF model accurately models plastic materials, but does not capture multiple scatterings in microfacets and metallic materials. For analytical pBRDF models, the SCATMECH C++ class library is commonly used for in-plane and out-of-plane light scattering analysis [9]. SCATMECH contains simple and complex pBRDF models, which can be used for multilayer stacked materials.

Representing a Mueller matrix as an ensemble, or set, of Jones matrices, has been appreciated for decades [10]. Mueller matrices that are non-depolarizing can be expressed in a Jones formalism; this class of matrices is referred to as Mueller-Jones matrices. The particular choice of a convex sum of up to four Jones matrices, called the Cloude coherency matrix, first appeared nearly 30 years ago to characterize depolarization in RADAR polarimetry [11]. In that paper, the entropy of the coherency eigenvalues was proposed as a metric for the magnitude of depolarization. The concept of entropy describes the degree of randomness in a process. For example, the entropy is zero if the Mueller matrix is a Mueller-Jones matrix, and the eigenspectrum only has one non-zero value. If all eigenvalues are equal, the entropy equals one, which corresponds to the ideal depolarizer Mueller matrix. Other metrics of depolarization, such as the depolarization index, are related to the eigenspectrum entropy [12,13].

Mathematically, most pBRDF models are a weighted sum of two or three Mueller matrices, where each weight is a free parameter of the model. Li et al. developed a pBRDF model where the relative weights of an ideal depolarizer and a Mueller-Jones matrix are derived from the largest coherency eigenvalue. This pBRDF model demonstrates the utility of the largest coherency eigenvalue in specifying the fractional contribution of specular versus diffuse reflection given a triple-degenerate eigenspectrum. The depolarization parameter varies between 0.25 and 1.0, where a small and large value correspond to interactions where diffuse and specular dominate, respectively. In addition, Li et al. demonstrated the largest coherency eigenvalue dependence on surface texture, albedo, and acquisition geometry. As albedo increases, the polarization entropy increases due to additional subsurface scattering events. The results from Li et al. are consistent with Umov’s effect, which states that the degree of polarization is inversely related to albedo [14]. In Li et al., estimating the largest coherency eigenvalue from as few as two polarimetric measurements is suggested for pBRDF extrapolation from partial polarimetry.

This work uses the KAIST 16-element Mueller matrix pBRDF database to validate a TD pBRDF model. This model compresses Mueller matrices to eight parameters and constrains interpolated values at unmeasured acquisition geometries to physical Mueller matrices. Furthermore, a computational method of polarimetric importance sampling the TD pBRDF model is presented and demonstrated.

3. Methods

3.1 Operational definition of Mueller matrices

For linear light-matter interactions, the polarimetric measurement equation relates a $4\times 4$ Mueller matrix $\mathbf {M}$ to a scalar-valued intensity measurement

Here, $p$ is a noise-free measurement from a polarimeter, $\dagger$ denotes the transpose, $\mathbf {a}$ and $\mathbf {g}$ are $4\times 1$ vectors of Stokes parameters that describe the polarization state analyzer (PSA) and the polarization state generator (PSG), respectively. The PSG is the incident polarization state. The PSA is the polarization-sensitive transmission (e.g. diattenuation) of the optics between the light-matter interaction and the detector. A radiometric (i.e. nonpolarimetric) version of Eq. (1) sets the PSA/PSG to unpolarized states; $\mathbf {g}=\mathbf {a}^{\dagger }=[1, 0, 0, 0]^{\dagger }$. The measurement of this unpolarized light-matter interaction is $p=[\mathbf {M}]_{00}=M_{00}$, which is called the average throughput (i.e. reflectance or transmission). The leftmost column of the Mueller matrix is called the polarizance vector [5]. This is the Stokes vector of exitant light when unpolarized light is incident. The top row of the Mueller matrix is called the diattenuation vector. This is the Stokes vector that forms an inner product with the incident Stokes vector to compute the radiant exitance.

For a Stokes vector of the form $\mathbf {S}=[S_0,S_1,S_2,S_3]^{\intercal }$ the degree of linear polarization (DoLP) is the ratio of the linearly polarized radiance to the total radiance ($S_0$)

The polarizance DoLP of Fresnel transmission is used in pBRDF modeling for partially polarized diffuse scattering and does not depend on illumination direction.

Polarizance DoLP versus angle of incidence (AOI) for the specular and diffuse models are plotted in Fig. 1. Here, the $\mathrm {DoLP_s}$ increases up to Brewster’s angle. For the diffuse model, $\mathrm {DoLP_d}$ monotonically increases versus AOI.

 

Fig. 1. The blue and orange lines are the viewing angle dependence of polarizance DoLP for specular and diffuse reflection, given in Eq. (4) and Eq. (5); respectively. Here ${n_0} = 1.55$ and $\mathrm {DoLP}_s$ peaks at Brewster’s angle $\theta _B = \mathrm {atan}^{-1} (n_0)$ while $\mathrm {DoLP}_d$ increases monotonically with viewing angle.

Download Full Size | PDF

The angle of linear polarization (AoLP) is

The Mueller matrix is a $4\times 4$ real value matrix that generally has 16 independent parameters and relates an input Stokes vector to an output Stokes vector. If fully polarized light remains fully polarized, the light-matter interaction is coherent (e.g. non-depolarizing) and the output Stokes vector can be written as

The annotation $\widehat {\cdot }$ differentiates a non-depolarizing Mueller matrix from a Mueller matrix with depolarization. The non-depolarizing Mueller matrix has 7 independent parameters: 6 polarimetric parameters plus the average throughput. These matrices are also called pure Mueller matrices or Jones-Mueller matrices because the 7 independent parameters are expressed as a Jones matrix $\mathbf {J}$ related to the Mueller matrix by

Note that both $\mathbf {J}$ and $e^{i\zeta }\mathbf {J}$ are solutions to Eq. (8). In other words, the conversion from a Jones matrix to a Mueller-Jones matrix is insensitive to absolute phase.

Eigenanalysis is the de facto standard method for the analysis of linear systems. However, the eigenvectors of a Mueller matrix are not constrained to be physically-realizable Stokes parameters. Instead, the Mueller matrix can be converted to a Hermitian matrix where eigenanalysis is used to express the Mueller matrix as a convex sum of up to four Mueller-Jones matrices. The Cloude coherency matrix is

The inverse relation from the Cloude coherency matrix to the Mueller matrix is

The coefficients of the Mueller matrix are $M_{kl}=\mathrm {tr}\left (\boldsymbol {\Pi }_{kl}\mathbf {H}\right ).$ The coefficients of the coherency matrix are $H_{kl}=\frac {1}{4}\mathrm {tr}\left (\boldsymbol {\Pi }_{kl}\mathbf {M}\right ).$ The eigendecomposition of the Cloude coherency matrix is

Each coherency eigenvector $\mathbf {c}_{n}$ is related to a Jones matrix

A Mueller-Jones matrix is denoted $\widehat {\mathbf {M}}$ because its Cloude coherency matrix is rank one. In other words, the Mueller-Jones matrix can be expressed as a single Jones matrix.

3.2 Triple degeneracy in the coherency eigenspectrum

There are infinite decompositions to express a given Mueller matrix as a weighted sum of other Mueller matrices. The eigendecomposition of the Cloude coherency matrix is a convex sum of up to four Mueller-Jones matrices that can express a given Mueller matrix [11,15]

Here, the coherency eigenvalues are the weights of the Mueller-Jones matrices and the sum is convex because these weights sum to one, $\sum ^{R}_n\xi _n=1$. The lowercase in $\widehat {\mathbf {m}}$ denotes the unit throughput ($[\widehat {\mathbf {m}}]_{00}=1$), which is a simple normalization of the Mueller matrix by $[\mathbf {M}]_{00}$. Each normalized Mueller-Jones matrix $\widehat {\mathbf {m}}_n$ is calculated from the coherency eigenvectors. The convention of descending order results in $\xi _0\geq \xi _1\geq \xi _2\geq \xi _3$. Therefore, the Mueller-Jones matrix $\widehat {\mathbf {m}}_0$ will be referred to as the dominant coherent process, since the associated eigenvalue is largest.

In the case where the coherency eigenspectrum is triply degenerate (TD) $\xi _1=\xi _2=\xi _3$, the Mueller matrix simplifies to

The TD Mueller matrix in Eq. (17) can be converted to a coherency matrix using Eq. (9) to yield

A Mueller matrix measurement of a pink silicone sphere is shown in Fig. 2(a). To aid in interpretation, the $M_{00}$ element is replaced by an RGB image of the object at an arbitrary scattering geometry. The RGB950 Mueller imaging polarimeter was used to make these measurements at 524 nm [16,17]. Scattering is in-plane and $\theta _d=20^{\circ }$ for the center pixel. The TD compressed Mueller matrix image is shown in Fig. 2(b). The DoLP images of the polarizance vector (i.e. first column of the Mueller matrix) are compared in Fig. 2(c-e). The AoLP images of the polarizance vector are compared in Fig. 2(f-h). For both AoLP images, the orientation of polarization is approximately normal to the surface of the sphere, which is consistent with a diffuse polarized scattering model. The average DoLP and AoLP deviation between these images (see Eq. (24) and Eq. (25)) is 0.5% and 21$^{\circ }$. The $\mathrm {AoLP}^{TD}$ in Fig. 2(f) has smoother spatial variation closer to the on-specular scattering geometry, as compared to the $\mathrm {AoLP}^{RGB950}$ in Fig. 2(g). This indicates the sensitivity of AoLP to changes in scattering geometry at smaller angles in TD compression, which truncates the Mueller matrix to its dominant process, compared to the original. AoLP is often used in computer graphics and computer vision for shape reconstruction. The polarizance AoLP shown in Fig. 2(e) could be measured by a linear Stokes camera with unpolarized illumination. The polarizance AoLP in Fig. 2(f) is computed from the dominant process and thus a full Mueller measurement is required. Therefore, the enhanced sensitivity to surface geometry in Fig. 2(f) compared to Fig. 2(e) reflects the retention of more information from the Mueller measurement.

 

Fig. 2. Normalized Mueller measurement of a pink silicone sphere in a) measured at 524 nm. To aid in interpretation, the $M_{00}$ element is replaced by an RGB image of the sphere at an arbitrary scattering geometry. In b) the normalized measurement is compressed to the TD model, see Eq. (17). The DoLP images of the polarizance vector (i.e. first column of the Mueller matrix) are shown for c) RGB950 measurement ($\mathrm {DoLP^{RGB950}}$) and d) TD compressed ($\mathrm {DoLP^{TD}}$). The DoLP difference between the RGB950 and the TD compressed images is shown in e). The AoLP images of the polarizance vector are shown for f) RGB950 measurement ($\mathrm {AoLP^{RGB950}}$) and g) TD compressed ($\mathrm {AoLP^{TD}}$). The difference in AoLP between the original and TD compressed image is shown in h). The DoLP and AoLP deviation (see Eq. (24) and Eq. (25)) of the sphere is 0.5% and 21$^{\circ }$, respectively. The $\mathrm {AoLP^{TD}}$ has a smooth spatial variation closer to the center of the sphere, compared to the $\mathrm {AoLP^{RGB950}}$. This indicates the sensitivity of AoLP to changes in scattering geometry at smaller angles in TD compression compared to no compression. For both AoLP images, the orientation of polarization is approximately normal to the surface of the sphere, which is consistent with diffuse scattering.

Download Full Size | PDF

3.3 Polarimetric importance sampling

Equation (16) can be restated as the equivalence between a normalized Mueller matrix and the average value of an ensemble of normalized Mueller-Jones matrices

If $\mathbf {L}$ is the ray tracing process, then implementing the RHS of Eq. (20) makes it possible to model depolarization effects using Jones calculus for light-matter interactions since $\widehat {\mathbf {m}}$ is a Mueller-Jones matrix.

Depolarization is often described as the result of physically averaging over disorder in a system, e.g. many scattered ray paths inside an integrating sphere. This work is original by treating the Mueller-Jones matrix (or Jones matrix) itself as a random variable in an otherwise deterministic description of the pBRDF. Depolarization is computed as an expectation, or average value, of the random Mueller-Jones matrix in the ray trace. Performing an average over the ray trace result, given a Mueller-Jones matrix for each light-matter interaction, allows coherent states to be propagated through the ray trace and averaged at the detector plane. Similarly to sampling ray geometry based on the BRDF, this proposed method importance samples the polarization-dependent parameters of a material’s pBRDF.

Normalization by $M_{00}$ in Eq. (16) is equivalent to factoring out the radiometric BRDF from the pBRDF. This factorization ensures that the eigenvalues sum to one and is necessary for a probabilistic interpretation. An advantage of this normalization is that the normalized polarization model can be scaled by an established BRDF model.

Consider a light-matter interaction for a probabilistic depolarization model with the triple-degeneracy assumption. Here, a nondepolarizing Mueller-Jones matrix and the ideal depolarizer Mueller matrix are weighted by two associated probabilities. When this ray hits a material, a Mueller-Jones matrix $\widehat {\mathbf {m}}$ is selected according to the probabilistic distribution.

Mueller calculus is required to represent partially-polarized states. An advantage offered by polarimetric importance sampling is that light-matter interactions are either completely polarized or unpolarized. Rays can be combined at the detector plane after the ray trace terminates. This would make it possible to compute depolarization from a radiometric ray trace combined with a polarimetric ray trace using Jones calculus. Jones versus Mueller calculus is computationally advantageous due to a few parameters in the matrices (7 versus 16). Also, quaternions are applicable to Jones (or non-depolarizing Mueller) matrices but not to partially polarized Mueller matrices [18]. Quaternion rotations are a common tool that is used to speed up computer graphics computations.

3.4 KAIST pBRDF dataset and model geometry

The KAIST dataset stores the wavelength-dependent pBRDFs for 25 materials as a 6-dimensional tensor where the scattering geometry is specified in a Rusinkiewicz coordinate system [19]. The pBRDF is defined in the local tbn-space in the direction of light travel, as shown in Fig. 3, where $\widehat {\boldsymbol {\omega }}_i$ and $\widehat {\boldsymbol {\omega }}_o$ are the incident and outgoing directions, respectively. The microfacet model postulates that for every $\widehat {\boldsymbol {\omega }}_i$ to $\widehat {\boldsymbol {\omega }}_o$ there is a subresolution facet that creates a specular reflection [9]. The surface normal of a microfacet is called the halfway vector $\mathbf {h}$ and is calculated by

For in-plane scattering, the zenith angle of the halfway vector $\theta _h=0$. Since the materials are isotropic, they are invariant to the azimuth angle of the halfway vector $\phi _h$. Changing this halfway azimuth angle would be equivalent to rotating a sample about its surface normal. The angle of incidence on the facet is called the difference angle $\theta _d$ and is related to the incident direction and facet surface normal by

The azimuth difference angle $\phi _d$ is the azimuth angle of $-\widehat {\boldsymbol {\omega }}_i$ after a rotation to make $\mathbf {h}$ and $\mathbf {n}$ parallel. Isotropic BRDF databases often assume a bilateral symmetry that truncates the range of $\phi _d$ to $[0,\pi ]$. Polarimetry breaks bilateral symmetry. Scattering above, or below, the plane of incidence are distinct geometrical transforms for linear polarization. Therefore, the pBRDF database parameterizes $\phi _d$ in the range $[-\pi,\pi ]$.

 

Fig. 3. The incident and outgoing directions, $\widehat {\boldsymbol {\omega }}_i$ and $\widehat {\boldsymbol {\omega }}_o$, are in the direction of light travel. Using these zenith and azimuth difference angles ($\theta _d$ and $\phi _d$) as well as the zenith halfway angle ($\theta _h$) to specify acquisition geometry is a widely used convention called Rusinkiewicz coordinates [19].

Download Full Size | PDF

Not all entries in the KAIST pBRDF database correspond to physical Mueller matrices. A necessary and sufficient condition for the Mueller matrix to be non-physical is if any of the eigenvalues of the coherency matrix are negative [20]. The condition for Mueller matrices to be physical in this work differs from the method discussed by Baek et al., which instead used the percentage of valid Stokes matrices [20]. Table 1 reports the percentage of physical Mueller matrices for each material in the KAIST pBRDF database. The diattenuation is the polarization dependence of the reflectance. Objects with lower diattenuation, such as silicone and spectralon, have a higher percentage of physical Mueller matrices. Diffuse objects are easier to measure with a polarimeter than specular objects, such as gold or brass, which have a larger diattenuation. There are fewer physical Mueller measurements for materials with higher diattenuation (e.g. metals) because a larger dynamic range is required in the detection process. Non-depolarizing materials are also more likely to create nonphysical Mueller reconstructions because even slight perturbations by measurement noise can create measurements that are inconsistent with physical Mueller matrices. [5,21,22].

Table 1. Each material in the KAIST Mueller database is measured at 2,989,441 scattering geometries and four wavelengths. If the eigenvalues of the coherency matrix corresponding to each Mueller matrix are non-negative than the Mueller matrix is physical; see Eq. (13). There are fewer physical Mueller measurements for materials of higher diattenuation (e.g. metals) since a larger dynamic range in the detection process is required. Non-depolarizing materials are also more likely to create nonphysical Mueller reconstructions because even slight perturbations by measurement noise can create measurements that are inconsistent with physical Mueller matrices [5].

View Table | View all tables in this article

3.5 Proposed triply degenerate database

For each material, the entries of the KAIST Mueller database are stored as single precision elements that result in a 912Mb shape tensor ($361 \times 91 \times 91 \times 5 \times 4 \times 4$). Here, the dimensions correspond to ($\phi _d$, $\theta _d$, $\theta _h$, $\lambda$, and $\mathbf {M}$). The eight TD parameters replace the 16 elements of the Mueller matrix that compress the database to $361 \times 91 \times 91 \times 5 \times 4 \times 2$, resulting in a 456 Mb file per material. No assumptions are made regarding the dominant coherent process for each acquisition geometry. Here, the $4 \times 2$ matrix contains six real-valued numbers to describe the complex-valued $4 \times 1$ coherency eigenvector $\mathbf {c}_0$ (see Eq. (13)), the largest coherency eigenvalue $\xi _0$, and the average throughput $M_{00}$. The coherency eigenvector is unit length, and the absolute phase is arbitrary (see Eq. (18)), therefore six parameters specify $\mathbf {c}_0$. Therefore, the first element of $\mathbf {c}_0$ can be real-valued and computed from the magnitudes of the other three elements.

The procedure to form the TD pBRDF database from the original Mueller matrices starts by storing $M_{00}$ in the TD database and then the Mueller matrix is normalized by this number. Normalized Mueller matrices are converted to Cloude coherency matrices using Eq. (9) and an eigenanalysis is performed. This eigenanalysis reveals non-physical Mueller matrices as discussed in Sec. 3.1. The largest eigenvalue $\xi _0$ is rounded to 0.25 or 1.0 if the eigenvalue is out of this range. The largest eigenvalue is stored in the TD database. Any negative eigenvalues are set to zero, which is a standard physical constraint [15]. The vector $\mathbf {c}_0$ is the eigenvector of the Cloude coherency matrix associated with the largest eigenvalue. This vector is unit-length and is scaled by a phasor to make the first element real. The real and imaginary parts of the other three elements are stored in the TD database. To convert the TD database into a Mueller matrix $\mathbf {c}_0$ forms a Mueller-Jones $\widehat {\mathbf {m}}_0$ which is substituted in Eq. (17). The TD pBRDF databases, as well as modified Mitsuba 2 plug-ins, are available as supplementary materials, Dataset 1, [23]. A flow diagram illustrating the prior and proposed methods in this work is shown in Fig. 4.

 

Fig. 4. pBRDF evaluation flow diagram. Here, the top path is the existing evaluation sequence introduced by Baek et al. and currently implemented in Mitsuba 2. This method interpolates the Mueller matrix elements for non-measured geometries and does not apply a physicality constraint to the Mueller matrix. The path of the proposed method utilizes the TD pBRDF database. Instead of interpolating the elements of the Mueller matrix, the proposed method interpolates over the coherency eigenspectrum. After interpolation, either the coherency eigenspectrum is directly converted to a Mueller matrix or polarimetric importance sampling can occur (See Sec. 3.3) to select between the dominant process or an ideal depolarizer. Once converted, the Mueller matrix is then evaluated.

Download Full Size | PDF

4. Results

4.1 On-specular renderings

Renderings using unpolarized illumination and specular acquisition geometries were computed for chrome, ceramic alumina, pink silicone, and white silicone. The DoLP of the polarizance vector measures the degree of unpolarized light becoming linearly polarized after reflection; see Sec. 3.1 for more information. Depolarization describes a reduction in the degree of polarization by a light-matter interaction; see Sec. 3.2. The rendered polarizance DoLP versus AOI is shown in Fig. 5 using the pBRDF of the 16 original Mueller matrix elements, the pBRDF of 8 TD elements, and the pBRDF of 8 TD elements with polarimetric sampling; see Eq. (17) for the description of the TD elements. The depolarization parameter $\xi _0$ is also included. Both DoLP and $\xi _0$ increase with AOI for diffuse materials such as pink Fig. 5(c) and white Fig. 5(d) silicone. This monotonic increase is consistent with the Fresnel transmission model, see Fig. 1. Under red illumination, the DoLP for pink silicone is 2% less than white silicone for all measured angles, which is consistent with Umov’s effect, which is the inverse relation between albedo and polarizance.

 

Fig. 5. On-specular polarizance DoLP and the depolarization parameter $\xi _0$ versus angle of incidence (AOI) renderings from Mitsuba 2 at the wavelength 622nm for: a) chrome, b) ceramic alumina, c) pink silicone and d) white silicone materials. The results for the original 16-element Mueller matrix pBRDF (red), TD pBRDF (green), and the mean and standard deviation for the importance sampled TD pBRDF (blue) are compared. The black diamonds are the depolarization parameter $\xi _0$ from Eq. (17). The materials in a) and b) are close to nondepolarizing, and therefore $\xi _0$ is above 0.9 at all AOIs. The silicone materials in c) and d) have a depolarization trend that increases with AOI. The polarizance DoLP also usually increases with AOI with the exception of the peak at $57^{\circ }$ in b) which resembles Brewster’s angle as shown in Fig. 1. The maximum deviation between the TD and original DoLP is 0.08 and usually is much lower except for the chrome material shown in a). Only 45% of the chrome measurements are physical (see Table 1), creating unavoidable deviations since TD is physically constrained. For a relatively small number of ray samples per pixel (spp), spp = 64, the polarimetric importance sampled TD converges to the partially depolarizing Mueller matrix TD pBRDF.

Download Full Size | PDF

For chrome, Fig. 5(a) and ceramic alumina, Fig. 5(b) $\xi _0 > 0.8$ for all measured AOI. Although both materials are highly nondepolarizing, the polarizance of the dominant process varies. Material properties such as refractive index and absorption affect the dominant process. In addition, nondepolarizing materials in the KAIST database have a lower percentage of physical Mueller matrices than depolarizing materials (see Table 1). When constructing the TD database, the coherency eigenvalue $\xi _0$ for the non-physical Mueller matrices is made physical following the methods detailed in Sec. 3.4. A larger percentage of coherency eigenvalues are replaced for materials with low depolarization, resulting in more discrepancies between the original and TD DoLP. The polarizance DoLP of the ceramic material under specular illumination tends toward the DoLP of the polarizance vector of the Fresnel reflection Mueller matrix with $n_0 = 1.55$. In all cases, the probabilistic implementation of the 8-parameter TD pBRDF converges to the deterministic results using 64 spp.

4.2 Convergence test of polarimetric importance sampling

The convergence properties of the polarimetric importance sampled TD pBRDF (see Eq. (21)) with respect to spp is shown Fig. 6. Here, the percent error of the TD pBRDF DoLP with and without importance sampling is computed. The error bars are one standard deviation from ten independent trials of the polarimetric importance sampling. The geometry of the sampled ray is fixed in a specular configuration ($AOI = 30^{\circ }$), hence the TD pBRDF DoLP does not change with spp. The rendered DoLP using polarimetric importance sampling for ceramic alumina, Fig. 6(a), reaches an asymptotic value for a lower spp compared to pink silicone, Fig. 6(b). In general, a lower $\xi _0$ indicates more depolarization and a higher entropy of the coherency eigenspectrum. In the extreme cases that $\xi _0=1$ or $\xi _0=0.25$ the probabilities for a non-depolarizing or completely depolarizing interaction become certainties (see Eq. (21)) and there is no probabilistic variance. The probabilities of a nondepolarizing or completely depolarizing interaction are equal when $\xi _0=0.375$. Convergence to the same solution as a partially-depolarizing Mueller matrix can be expected to require a higher spp as these probabilities approach equality.

 

Fig. 6. The rendered probabilistic importance sampled polarizance DoLP ($\mathrm {DoLP^{Prob}}$ error percentage as a function of spp. $\mathrm {DoLP^{Prob}}$ is compared to the rendering without importance sampling. Here, the sampled ray geometry is fixed in a specular configuration ($AOI = 30^{\circ }{}$), therefore, the non-importance sampled DoLP does not change with spp. The error between the two renderings is shown for difference materials in a) and b) where the error bars are one standard deviation from ten independent trials of the TD pBRDF computed with importance sampling. The percentage error approaches zero as the number of spp increases. The TD pBRDF rendering with importance sampling converges at a lower spp for ceramic alumina compared to the pink silicone.

Download Full Size | PDF

4.3 Methods and figures of merit to compare rendered images

Renderings are generated using Mitsuba 2, a polarization-aware PBR engine [1]. The number of spp varies for each scene to reach visually assessed convergence. Stokes images are rendered from an unpolarized source. For $s_0$, a composite $S_0$ of the red, green, and blue wavebands is formed and then normalized by the largest pixel value in each color band. The AoLP and DoLP images rendered at a wavelength of 633 nm are shown in Figs. 7, 8, and 9. Each scene is rendered using the pBRDF of the original Mueller database, the TD compressed database, and TD with polarimetric importance sampling, as described in Sec. 3.3.

 

Fig. 7. Renderings of RGB, DoLP, and AoLP images (left to right columns) of a pink silicone sphere for: a) the original pBRDF database, b) the compressed TD pBRDF database, and c) the TD pBRDF database computed with polarimetric importance sampling. The rendered AoLP is approximately normal to the sphere, which matches the diffuse light scattering model (see Eq. (5)) and the measured AoLP shown in Fig. 2. The DoLP deviation (see Eq. (24)) as compared to the original pBRDF is 1.10% for both the TD pBRDF and the TD pBRDF with polarimetric importance sampling. The AoLP deviation (see Eq. (25)) as compared to the original pBRDF is 16.90$^{\circ }$ for both the TD pBRDF and the TD pBRDF with polarimetric importance sampling.

Download Full Size | PDF

 

Fig. 8. Rendering of the Baek et al. 2020 robot scene using a) the original KAIST, b) the compressed TD, and c) the TD pBRDF database with polarimetric importance sampling. Umov’s effect is evident with the low DoLP for the red pieces of the robot and a high DoLP for the blue pieces [7,14]. The DoLP deviation (see Eq. (24)) as compared to the original pBRDF database is 5.89% for the TD pBRDF and 5.94% for the TD pBRDF with polarimetric importance sampling. The AoLP deviation (see Eq. (25)) as compared to the original pBRDF database is 19.22$^{\circ }$ for the TD pBRDF and 20.41$^{\circ }$ for the TD pBRDF with polarimetric importance sampling.

Download Full Size | PDF

 

Fig. 9. Rendering of RGB, DoLP, and AoLP (left to right columns) of a room with a light source behind an open door developed by Benedikt Bitterli. In a) original KAIST, b) the compressed TD, and c) TD pBRDF with polarimetric importance sampling. The teapot material (left-right) isis chrome, ceramic alumina, and pink silicone. The DoLP deviation (see Eq. (24)) compared to the original pBRDF database is 2.2% for the TD pBRDF and 3.0% for the TD pBRDF with polarimetric importance sampling. The AoLP deviation (see Eq. (25)) compared to the original pBRDF database is 30.0$^{\circ }$ for the TD pBRDF and 34.1$^{\circ }$ for the TD pBRDF with polarimetric importance sampling.

Download Full Size | PDF

When user-defined pBRDF databases are input to Mitsuba 2, continuous random variables are sampled to create independent ray directions. This process requires interpolation of the pBRDF database to unmeasured acquisition geometries. There are various methods to interpolate complex numbers and each offers different advantages [24,25]. Interpolating the real and imaginary parts of two complex numbers forms a line between them in the complex plane. Interpolation of the magnitude and phase of two complex numbers forms an arc between them in the complex plane. In this work, linear interpolation of the real and imaginary parts of the six parameters $\mathbf {c}_0$ is required, as phase unwrapping is not a current functionality of Mitsuba 2 pBRDF interpolation. Unfortunately, this constraint yields lower interpolated magnitudes as compared to polar coordinate interpolation. These lower interpolated magnitudes created higher interpolated magnitudes when the first element of $\mathbf {c}_0$ was calculated from the magnitudes of the other three elements. This bias in the interpolated magnitudes created a distinct geometric region where the polarizance DoLP was incorrectly zero. Overparameterizing the TD model by including an additional parameter for $\mathbf {c}_{0}$ minimizes this artifact. Phase unwrapping for complex-valued parameters in the pBRDF interpolation would be a more elegant solution, which is not available at this time. Therefore, the TD model used in this work contains nine parameters with eight degrees of freedom. Using nine parameters, the dimensions for each material in the database is $361 \times 91 \times 91 \times 5 \times 3 \times 3$ where the file size is 513.2Mb. The pBRDF plug-ins and material models can be found in the supplementary materials, Dataset 1, [23].

The rendered images from various pBRDF databases are quantitatively compared using two figures of merit: the average DoLP and AoLP standard deviation. The average DoLP variance between two images: DoLP$^{a}$ and DoLP$^{b}$ over all pixels is

The tabulated results of $\sigma _{\mathrm {DoLP}}$ and $\sigma _{\mathrm {AoLP}}$ for the scenes rendered in this work are shown in Table 2. The number of rays required for the probabilistic importance sampling to converge to the partially-depolarizing Mueller matrix implementation is a function of the scene geometry and ray path depth. The sphere uses direct illumination, where there is only a single surface reflection. For visual convergence, all rendered spheres use an spp = 512. For the robot and room scene, spp = 2,048 and spp = 4,096 for convergence, respectively, for all rendered pBRDFs excluding the probabilistic pBRDF. Here, the probabilistic pBRDF requires additional rays and converges when spp = 4,096 and spp = 16,384 for the robot and room scene, respectively. These scenes include indirect illumination configurations, which require additional rays for the polarimetric importance sampled pBRDF to converge.

Table 2. Figures of merit to compare renderings from three different pBRDF databases. The average deviation of DoLP and AoLP between the images $a$ and $b$ is denoted by $\sigma (\mathrm {DoLP}^{a},\mathrm {DoLP}^{b})$ and $\sigma (\mathrm {AoLP}^{a},\mathrm {AoLP}^{b})$ in Eq. (24) and Eq. (25), respectively. In the table superscripts denote the original pBRDF, TD pBRDF, and TD pBRDF computed using polarimetric importance sampling.

View Table | View all tables in this article

4.4 Sphere renderings

Polarimetric renderings of polarizance AoLP and DoLP of a pink silicone sphere are shown in Fig. 7 using the original KAIST, the TD, and the TD computed with importance sampling pBRDF. The angle between the light source and camera is $25^{\circ }$. Table 2 reports summary statistics of the comparisons. The DoLP deviation (see Eq. (24)) as compared to the original pBRDF is 1.1% for both the TD pBRDF and the TD pBRDF with polarimetric importance sampling. The AoLP deviation (see Eq. (25)) as compared to the original is 16.9$^{\circ }$ for both the TD pBRDF and the TD pBRDF with polarimetric importance sampling.

In Fig. 7(c) the DoLP and AoLP polarizance renderings are computed using polarimetric importance sampling with 512 spp visually converging to the partially-depolarizing Mueller matrix rendering in Fig. 7(b). The rendered polarizance AoLP is normal to the surface and the DoLP increases radially for this material due to diffuse reflection. Measurements of a pink silicone sphere reported in Fig. 2 have a consistent polarizance AoLP and DoLP dependency. This results suggest that the dominant Mueller-Jones matrix is similar to Fresnel transmission. Both the TD compression of the measurement and the rendering from the KAIST database maintains most of this pattern in the polarizance AoLP and DoLP. An area of disagreement in the AoLP is close to the center of the sphere. The variation in AoLP is larger, compared to the measurement and the rendering from the original KAIST database, for the TD pBRDF. At these acquisition geometries, the viewing angle is relatively small and the TD pBRDF polarizance AoLP in Fig. 7(b) is more sensitive to changes in the surface geometry, as compared to the renderings from the original KAIST database in Fig. 7(a). A similar effect was noted and discussed in Fig. 2 when measurements were compressed to the TD model. Another area of disagreement is along the edges of the sphere. Here, the TD assumption does not capture the variations of the DoLP for larger acquisition angles which leads to higher values compared to the original pBRDF.

4.5 Scene renderings

Two scenes containing multiple surface interactions are rendered in Fig. 8 and Fig. 9. All objects use materials from the KAIST pBRDF database. Figure 8 replicates the robot scene in Baek 2020 et al. Figure 9 is a scene from Bitterli [26] where three teapots sit on a table and light enters the room from an ajar door. The teapot material (left to right) is: chrome, ceramic alumina, and pink silicone. Both scenes include indirect scattering events, so the path depth is unlimited. For the probabilistic implementation, Fig. 8 and Fig. 9, require spp = 4,096 and spp = 16,384, respectively, to reach the visually assessed convergence. The variation in spp requirement is due to ray depth. The deviation between the renderings using the TD pBRDF and the original KAIST pBRDF is reported in Table 2. In addition to only retaining the dominant process of the KAIST database, the TD database physically constrains the interpolated and non-interpolated Mueller matrices. Both the compression and interpolation techniques used for the TD renderings result in differences between the KAIST and TD rendered images.

In Fig. 8, the DoLP deviation (see Eq. (24)) as compared to the original pBRDF database is 5.9% for the TD pBRDF and the TD pBRDF with polarimetric importance sampling. Umov’s effect is evident with the low DoLP for the red pieces of the robot and a high DoLP for the blue pieces [7,14]. The AoLP deviation (see Eq. (25)) as compared to the original pBRDF database is 19.2$^{\circ }$ for the TD pBRDF and 20.4$^{\circ }$ for the TD pBRDF with polarimetric importance sampling. The most notable area of AoLP disagreement between the original and TD renderings is in the upper right corner on the shaded wall. The green triangle and white stripes on the walls are made of mint silicone and ceramic alumina, respectively. The AoLP renderings of mint silicone are in closer agreement with those of ceramic alumina. Differences between original and TD pBRDFs for ceramic alumina were observed at low AOI in Fig. 5(b). In the shadowed region, more noise is observed in the probabilistic TD rendering than the partially depolarizing Mueller matrix rendering.

In Fig. 9 the DoLP deviation (see Eq. (24)) as compared to the original pBRDF database is 2.2% for the TD pBRDF and 3.0% for the TD pBRDF with polarimetric importance sampling. The AoLP deviation (see Eq. (25)) as compared to the original pBRDF database is 30.0$^{\circ }$ for the TD pBRDF and 34.1$^{\circ }$ for the TD pBRDF with polarimetric importance sampling. The rendered polarizance AoLP of the TD pBRDF is more sensitive to surface geometry compared to that of the original KAIST pBRDF. The same effect was observed in Fig. 2 and Fig. 7 for the measured and rendered spheres; respectively. The differences between TD AoLP in Fig. 9(b) and original pBRDF in Fig. 9(a) can be observed on the door where AoLP = $90^{\circ }$ at the center of the door. For the TD renderings, the AoLP rotates counterclockwise toward 135$^{\circ }$ and clockwise toward 45$^{\circ }$ for the top and bottom of the door, respectively. The rendering using the original KAIST pBRDF is insensitive to this geometric change in the view direction. The illumination on the left side of the scene is caused by indirect scattering events. In this region of the scene, the ray depth is highest, and the polarimetric importance sampled rendering requires more rays to achieve the same visual appearance as the non-importance sampled TD rendering.

5. Conclusion

The contributions of this work are: compression which retains dominant polarimetric properties, interpolation to non-measured geometries which maintains physical constraints, and convergence demonstrations of polarimetric importance sampling.

The materials in the KAIST database include highly diffuse depolarizing materials (e.g. pink silicone), as well as nondepolarizing metals (e.g. chrome); see Table 1. An eight-parameter physically constrained model for $4\times 4$ Mueller matrices for the 25 materials in the KAIST pBRDF database has been tested. Comparisons between the original and compressed TD databases were done by generating polarimetric renderings of scenes containing multiple surface interactions, non-trivial ray depth, and varying surface geometries. For a rendered pink silicone sphere with a single surface interaction, for a rendered pink silicone sphere from the KAIST pBRDF database with a single surface interaction, the average DoLP deviation was 1.1% for both the TD compressed and probabilistic TD compressed pBRDF model as compared to the uncompressed model. The deviation in AoLP between the TD and the original pBRDF is 16.9$^{\circ }$. Other rendering comparisons are reported in Table 2.

The retention of physically constrained Mueller matrices when interpolating to nonmeasured geometries is an advantage offered by the TD database. However, the implementation in this work was limited by the need for phase unwrapping to interpolate complex numbers. A 9-parameter version of the TD database was created that includes a dependent variable to avoid interpolation artifacts. A more complete and satisfactory solution would require an interpolation scheme designed for complex-numbers; see Sec 4.3.

The expression of the Mueller matrix for the TD model can be restated as the average value of a normalized Mueller-Jones matrix and a completely depolarizing Mueller matrix. A new method of polarimetric importance sampling was demonstrated by relating the depolarization parameter to the probability of a coherent/fully polarized or unpolarized light-matter interaction. Mueller calculus is required to represent partially polarized states. An advantage offered by polarimetric importance sampling is that light-matter interactions are either completely polarized or unpolarized. Rays can be combined at the detector plane after the ray trace terminates. This would make it possible to compute depolarization from a radiometric ray trace combined with a polarimetric ray trace using Jones calculus. Although the probabilistic importance sampling in this work was computed using a Mueller polarimetric rendering engine, a potential computational advantage would be to use a Jones formalism. Jones calculus requires fewer parameters to describe material polarization properties and to track the polarization state upon light-matter interactions. In addition, quaternion mathematics that already exist for most PBR rendering engines may simplify and decrease the computation time for operations in a Jones formalism compared to Mueller calculus [18]. In this work, the pBRDF and ray direction were sampled independently. Convergence and realism could potentially be improved by coupling pBRDF and ray-direction importance sampling.