Merit functions and measurement schemes for single parameter depolarization models

Mueller polarized bi-directional scattering distribution functions (pBSDFs) are 4x4 matrix-valued functions which depend on acquisition geometry. The most popular pBSDF is a weighted sum between a Fresnel matrix and an ideal depolarizer. This work's main contribution is relating the relative weight between an ideal depolarizer and Fresnel matrix to a single depolarization parameter. Rather than a 16-dimensional matrix norm, this parameter can form a one-dimensional merit function. Then, instead of a full Mueller matrix measurement, a scheme for pBSDF fitting to only two polarimetric measurements is introduced. Depolarization can be mathematically expressed as the incoherent addition of coherent states. This work shows that, for a Mueller matrix to be in the span of a Fresnel matrix and an ideal depolarizer, the weights in the incoherent addition are triply degenerate. This triple degeneracy is observed in five different colored opaque plastics treated with nine different surface textures and measured at varying acquisition geometries and wavebands.


Introduction
Bi-directional scattering distribution functions (BSDFs) describe a material's radiometric response when illuminated and observed from different angles. This work proposes a polarized BSDF (pBSDF) model inspired by the spectral decomposition [2] of Mueller matrix measurements of opaque, diffuse, plastic materials. A Mueller pBSDF is a 16-element matrix that quantifies the polarization and depolarization effects of a light-matter interaction. Mueller matrices can describe non-depolarizing, partially depolarizing, and completely depolarizing light-matter interactions. A Mueller matrix has sixteen degrees of freedom (DoF) from which 7 are associated with nondepolarizing properties: 3 DoF for retardance, 3 DoF for diattenuation, and 1 DoF for throughput (e.g. reflectance, transmittance). The other 9 DoF are associated with depolarization [1,3]. Depolarization can be mathematically expressed by a convex incoherent sum of four or fewer coherent states. The diattenuation and retardance determine the most significant coherent state [1,4]. The weights in this convex sum are 3 of the 9 DoF for depolarization. When these weights are triply degenerate (i.e the smallest three are equal) the depolarization is reduced to 1 DoF. In this work, an approximate triple degeneracy is observed over variations in surface texture and albedo of 5 different colored opaque plastics treated with 9 different surface textures and illuminated by 662, 524, and 451nm wavebands. The albedo, texture, and scattering geometry change the depolarization DoF. Fresnel reflection reasonably approximates the non-depolarizing properties of all measurements.
The most popular pBSDF model is a weighted sum between a non-depolarizing Mueller matrix (e.g. Fresnel relfection) and an ideal depolarizer. Early Mueller pBSDF work by Bickel et al. defined polarizing interactions as arising from perturbations to an idealized smooth surface [5]. Bickel et al. pointed to an idealized Lambertian surface as the "opposite case" to a perfectly smooth surface and posited that all "real-world surfaces" lie between an ideally smooth and ideally Lambertian surface texture. A Fresnel reflection matrix and ideal depolarizer can be used to represent the ideally specular smooth and ideally diffuse Lambertian surface textures, respectively. A polarimetric interpretation of the original unpolarized microfacet BSDF model uses the ideal depolarizer and Fresnel reflection matrix as component Mueller matrices [6][7][8][9].

Background
The microfacet BSDF model was first developed to produce more accurate off-specular scattering in computer graphics to improve the appearance of object renderings [7,17]. The pBSDF models by Baek et al. [10] in 2018 and Kondo et al. [12] in 2020 produce well-matched visual renderings for non-polarized and linearly-polarized illumination sources and observers. Many forms of microfacet distribution functions have since been developed for different object types for applications beyond computer graphics, such as remote sensing [8,15,18,19]. A Mueller pBSDF model is also used in applications such as the realistic rendering of complex scenes, inverse models for material recognition and shape reconstruction, or synthetic training data for neural networks.
The component Mueller matrices and the microfacet distribution functions in a Mueller pBSDF are changed to improve the agreement between a Mueller pBSDF and measurements. Baek et al. proposed a model [10] which replaces the ideal depolarizer in a basic Mueller BSDF model, which is a sum of an ideal depolarizer and Fresnel reflection matrix, with a Mueller matrix that traces a hypothetical ray propagation path through a scattering material. This hypothetical ray path consists of Fresnel transmission into the material followed by bulk scattering represented by an ideal depolarizer and then a final Fresnel transmission out of the material. This model was used to improve surface normal estimation [10] when applied to a 3 × 3 partial Mueller matrix. This 3 × 3 partial Mueller matrix omits the fourth row and fourth column of the full 4 × 4 Mueller matrix. In 2020, Baek et al. tested a full 4 × 4 implementation of this model and observed disagreement in polarimetric accuracy and irradiance when compared to their measurements. Baek et al. 2020 also conclude that rendering a large variety of materials requires the use of look-up tables, which this work does not refute. The Kondo model [12] extends the Baek model [10] as a 4 × 4 model by reintroducing the ideal depolarizer as a third component Mueller matrix. Each component Mueller matrix is normalized, and the weights for each component are constrained to sum to the measured total luminance for an object. The Kondo et al. model is a 3 × 3 partial Mueller matrix model applied only to linear polarization states. Data created using this three-component model is used as part of a training data set for a convolutional neural network to create realistic polarimetric renderings of scenes with elaborate, fine detail. A full 4 × 4 Mueller extension of this model, which assumes sub-surface bulk scattering can be partially isotropically depolarizing, is used as one of the benchmark models in this work.

Materials of Varying Albedo and Texture
The object ensemble is a group of red, orange, yellow, green, and blue LEGO Duplo bricks which are roughened using a belt sander equipped with nine different grits of sandpapers. This process creates 45 distinct bricks which differ in color and texture. A white-light interferometer is used to measure each sample's surface profiles and mean surface roughness to quantify surface roughness; see Table 3 in the Appendix.
The 45 bricks pictured in Fig. 1 are measured using a custom large-aperture Mueller matrix imaging polarimeter called the RGB950 [20]. The RGB950 measures objects under narrow-band illumination at 662 ± 11.17, 524 ± 17.31, and 451 ± 9.78 nm wavelengths. A single Mueller matrix measurement is taken at a given source position, camera position, and illumination waveband. In this study, 30 unique combinations of camera position and sample rotation are selected. These 30 geometries are reported in Tab. 2 using angle of incidence onto the macrosurface and angle of exitance from the macrosurface . The and angles are referenced to the central position in the 3 × 3 brick tower in Fig. 1(b). An 11 × 11 pixel region of interest (ROI) from each brick image is selected for analysis. This corresponds to a 3.2 × 3.2mm projected area on each brick. Each ROI on the 3 × 3 brick towers pictured in Fig 1b has a slightly different incident and exitant propagation vectors. Randomizing the texture positions within an image allows verification of trends with texture that are not dependent on acquisition geometry.

Acquisition Geometry
In this work, unit vectors are denoted by a hat · and all vectors are boldface. All vectors point in the direction light travels (i.e. k-vector direction) [10,11]. The object surface normal is n which also referred to as a macronormal to distinguish from microfacet normals. The convention n = z = {0, 0, 1} is adopted. The incident propagation direction of the illumination is = [sin( ) cos( ), sin( ) sin( ), cos( )] where is the source zenith angle and is the source azimuth angle. Therefore backscattering configurations only occur when > 90 • . The exitant propagation direction after a light-matter interaction which reaches the camera is = [sin( ) cos( ), sin( ) sin( ), cos( )]. Therefore if = , the observer is looking into the source and the direction of light travel is unchanged.
Microfacets are planar structures on a larger macrosurface that only reflect in the specular direction. Since the angle of incidence and the angle of specular reflection are equal, the surface orientation produces a specular reflection defined by the incident and exitant propagation vectors. An ensemble of sub-resolution microfacets is used to model a polarized contribution to light scattering in off-specular directions. Surface texture influences the distribution of these microfacet orientations [21]. For example, the ensemble of microfacet orientations for a perfectly smooth mirror would have no deviation from the macronormal. The micronormal m is the surface normal of a microfacet on a larger macrosurface. A microfacet distribution, which can be interpreted as a probability density function (pdf) on m, are designed to describe various surface types, see Sec. C.
The halfway vector h is the surface normal of an object which would produce a specular reflection for a given source and camera position; see Fig. 2. This work adopts the notation for the halfway and difference angles introduced by Rusinkiewicz [22]. The difference angle is (a) Bricks ordered top-to-bottom from smoothest to the roughest texture and labelled T1-T9.
(b) Bricks arranged for measurements and texture location is shuffled. Fig. 1. A set of 45 objects with varying albedo and texture is created from five colors of LEGO DUPLO™ bricks sanded by nine grits of sandpaper. (a) Bricks ordered top-to-bottom from smoothest to roughest with texture labels T1-T9. (b) The texture location is shuffled when bricks are imaged; see Table 3 in Appendix D for details. For the convention of vectors pointing away from the surface, the incident propagation vector becomes − . The halfway vector h is the surface normal of an object which would produce a specular reflection. The surface normal of the measured object is n (which is parallel to z) and called the macrosurface normal. The halfway angle ℎ between h and n is the deviation of the macrosurface from a specular geometry. The difference angle is the angle between h and − . The difference angle is = (180 • − Ω)/2, i.e. half of the angle supplementary to the scattering angle Ω. The scattering angle Ω is the angle between the incident and exitant propagation vectors. the angle of incidence onto a specular microfacet defined as the angle between vectors − and h. The halfway angle is the angle between n and h. The cosine of the difference angle is equivalent to cos( ) = − · h. The difference angle itself is also equivalent to = (180 • − Ω)/2, or half of the supplementary angle for the scattering angle Ω. The scattering angle Ω is the angle between the incident and exitant propagation vectors calculated as cos (Ω) = · .

Polarization
Spectrally incoherent light can be fully polarized, partially-polarized, or completely unpolarized. The polarization state and irradiance of reflected, transmitted, and scattered light after a lightmatter interaction is, in general, dependent upon the polarization state of the incident light. Mueller calculus describes the polarization transformation of linear light-matter interactions. Mueller calculus is required instead of Jones calculus for describing partial polarization. Both polarizing and depolarizing effects from material interactions are described using the 4 × 4 Mueller matrix, while the polarization state of light is described using the 4 × 1 vector of Stokes parameters. A vector of Stokes parameters S ( ) is defined along a propagation direction, e.g. or . The four Stokes parameters which describe all possible polarization states of light are where the M 00 element is the reflectance for unpolarized incident light. A normalized Mueller matrix is divided by this reflectance m ( , , n) = M ( , , n)/ M ( , , n) 00 .
Here the normalized Mueller matrix is lowercase, similar to the normalized Stokes vector, and the notation [M] 00 = 00 is used. In general, the reflectance, and therefore M 00 , is a function of incident and exitant propagation directions which is related to the scalar-valued BRDF by [23] ( , Here is the BRDF in units of inverse steradians, the reflectance is unitless, and the differential projected solid angle is Ω = cos sin where and are the zenith and azimuth angles of the propagation vector . Here, the subscript indicates that the differential projected solid angle definition applies to both and . A normalized Mueller matrix operating on an incident Stokes vector S ( ) can be scaled by the reflectance to produce the exitant Stokes vector S ( ) = ( , )m( , , n)S ( ). . In (a), 662nm illumination of the red object is high albedo, and in (b) 451nm illumination is low albedo. The white dots' locations are the 30 capture geometries measured, and MATLAB's natural neighbor interpolation is used to generate images. The most pronounced differences between the Mueller matrices in (a) high albedo and (b) low albedo conditions are the diagonal elements at lower scattering angles. High albedo illumination produces more depolarization than low albedo illumination due to an increase in multiple scattering events. Both high and low albedo are anisotropic depolarizers since for isotropic depolarization, see Eq.13, the on-diagonal linear m 11 , m 22 and circular m 33 elements are equal.
The normalized Mueller matrix is an idealized Lambertian reflector since for all measurement geometries, the element [m] 00 equals one. Although the normalization loses absolute radiometry, the polarimetry of the light-matter interaction is preserved. The incident to exitant and the incident to exitant polarization ellipse are identical for both normalized and unnormalized Mueller matrices.
This work analyzes normalized Mueller matrices measurements at 30 geometries; see Fig.3 for examples. In this normalized form, the changes to the relative magnitude of the non-m 00 elements at varying geometries can be assessed independently from the reflectance changes. Measuring the Mueller matrix at a wide range of incident and exitant geometries requires varying the exposure settings with scattering angle. Varying exposure settings were selected to maximize the detector's dynamic range, thus maximizing the polarimetric accuracy. However, the absolute radiometry required for an unpolarized BSDF profile is lost when the exposure is varied. Any Mueller pBSDF model P( , , n) can be factored into a reflectance component ( , , n) and a normalized Mueller matrix model component p( , , n). The normalized pBSDF can be scaled by an unnormalized BSDF to return a polarization-dependent reflectance at every measurement geometry. Figure 3 compares normalized Mueller matrix measurements with 451nm and 662nm illumination for which the smoothest (T1) red brick is low and high albedo, respectively. These low and high albedo measurements differ most in the on-diagonal elements, which correspond to isotropic depolarization, see Eq. 13. Normalized Mueller matrix elements 11 , 22 , and < (1, 1.54) geometries. The same elements for the high albedo Mueller matrix in Fig 3a are ≈ 0. The difference is especially noticeable for measurements < 57 • . The low albedo Mueller matrix elements 01 , 02 , 10 , 20 , 23 , and 32 also have increased magnitude for geometries which are ≈ 0 for the high albedo Mueller matrix.
For the Mueller matrix visualization in Fig.3, measurements are represented in a signed halfway angle versus difference angle space. The halfway angle is specified between 0 • and 90 • , but in signed space, an artificial sign is attached to ℎ . This sign of ℎ is assigned positive when the azimuth angle of the halfway vector, ℎ , is between 90 • and 270 • ; the sign is assigned negative otherwise. This sign convention is further discussed in Appendix A.

Polarimetry
A Mueller matrix is estimated from a series of images acquired at varying Polarization State Analyzer (PSA) and Polarization State Generator (PSG) states [24]. For linear light-matter interactions the relationship between a noise-free scalar-valued irradiance and an object's Mueller matrix is = a Mg.
Here a is the 4 Stokes parameters describing the PSA, g is the Stokes parameters of the PSG, is a noise-free irradiance, and denotes the transpose of a real-valued vector. Consider forming a single 16 × 1 vector from the PSA and PSG Stokes vectors, Then the irradiance in Eq. 8 can be rewritten as the inner-product of vectors where the dependence of the irradiance on the Mueller matrix is written explicitly as (M) and ì M is a 16 × 1 vector of the Mueller elements. A series of irradiance values can be expressed as where W is a 16 × matrix called the polarimetric measurement matrix and each row can be written as a Kronecker product between the ℎ PSA/PSG

Depolarization
Depolarization refers to a reduction in the after a light-matter interaction. A Mueller matrix is depolarizing if the degree of polarization is greater for the incident light than for the exitant light, i.e.
For every Mueller matrix M there is a 4 × 4 complex-valued Hermitian coherency matrix C related by [4] where are Pauli-spin matrices and Here the dagger † indicates a complex conjugate transpose operation and ⊗ indicates a Kronecker product. Since the coherency matrix is Hermitian it can be written as where is the rank, are the real and non-negative eigenvalues in descending order, and c are the orthonormal eigenvectors of the coherency matrix. Each eigenvector of the coherency matrix in Eq. 16 corresponds to a Jones matrix where = 0, 1, 2, 3 indicates the ℎ eigenvector and [·] indicates the ℎ element of the vector. A Jones matrix is related to a non-depolarizing Mueller matrix (i.e. Mueller-Jones matrix), by where the · notation indicates that the coherency matrix of any Mueller-Jones matrix is rank one. The Mueller-Jones matrices computed from the orthogonal eigenvectors of the coherency matrix can be incoherently summed to express the original Mueller matrix This treatment is also called spectral decomposition. Spectral decomposition separates a depolarizing Mueller matrix into non-depolarizing parts. Alternate decompositions of a partially depolarizing Mueller matrix split M into the combination of a fully depolarizing component and either a fully depolarizing D(0) component [27] or an isotropically depolarizing component D( ) [28]. Each m is the Mueller-Jones matrix associated with the eigenvector c in Eq. 16 which is constrained to be normalized. The sum of the eigenvalues, Polarization entropy is a scalar value between 0 and 1 which is related to the depolarization index of a Mueller matrix and is calculated from the coherency matrix eigenspectrum [4,29,30].
where a normalized Mueller matrix m is used to calculate , which constrains the sum of the eigenvalues to one. A polarization entropy of zero indicates the Mueller matrix is a Mueller-Jones matrix because the coherency matrix is rank one, i.e. only one eigenvalue of the coherency matrix is non-zero. A Fresnel reflection matrix is an example of a Mueller-Jones matrix.  Figure 5 shows the entropy and largest normalized eigenvalue 0 which corresponds to the largest basis Mueller matrices shown in Fig. 4. The fractional contribution of m 0 to the measurement is given by the associated eigenvalue 0 . At large on-specular scattering angles, 0 approaches one. For small off-specular scattering geometries, 0 is minimized. The eigenvalue 0 is smaller for high albedo measurements than low albedo measurements.

Triply Degenerate Eigenspectrum
Four measurements are compared in Fig.5 to illustrate how polarization entropy is dependent on measurement geometry, albedo, and surface texture. For both albedos, polarization entropy increases as decreases or as the measurement geometry moves away from specular (i.e. as the halfway angle ℎ increases). As polarization entropy increases, the largest eigenvalue decreases. The approximate triple degeneracy 1 ≈ 2 ≈ 3 applies to both albedos and textures and for all measurement geometries. For the high albedo cases, the eigenvalues approach equal magnitude as the entropy approaches one. At each measurement geometry, a higher albedo yields a higher polarization entropy due to increased bulk scattering. Figure 6 shows that polarization entropy also generally increases as the surface texture becomes rougher. Other authors have established observations and relations between surface texture and polarized light scattering [21].
There is a unique Mueller matrix solution for maximum polarization entropy. The polarization entropy equals one when all coherency eigenvalues are equal. If a flat eigenspectrum of = 1/4 for = 0, 1, 2, 3 is applied to Eq. 16, then the coherency matrix is proportional to the identity matrix and the associated Mueller matrix is proportional to the ideal depolarizer D(0), given in Eq. 13. Substituting the ideal depolarizer's normalized eigenspectrum [1/4,1/4,1/4,1/4] into Eq. 19 yields Here a Mueller matrix complementary to a normalized Mueller-Jones matrix has been defined by The complementary matrix is neither a normalized matrix nor a Mueller-Jones matrix, but is a physically realizable Mueller matrix [31,32]. This separation between a normalized Mueller-Jones matrix and its complementary matrix is motivated by the triply degenerate eigenspectrum, 1 ≈  increases. At a given measurement geometry, the polarization entropy is greater for high albedo objects due to more bulk scattering. The approximate triple degeneracy relationship 1 ≈ 2 ≈ 3 applies to both albedo objects for all measurement geometries.
Here 1/4 < 0 ≤ 1 because it is the largest eigenvalue of the coherency matrix of a normalized Mueller matrix. When a triply degenerate eigenspectrum is substituted into Eq. 19 the normalized Mueller matrix expression is The two terms in these equivalent

Measuring the Largest Normalized Eigenvalue
The triple degeneracy of a Mueller matrix is an important constraint which leads to simplified expressions for measuring the normalized largest eigenvalue. Given this triply degenerate assumption, the Mueller matrix is The dominant normalized Mueller-Jones matrix m 0 , the normalized largest eigenvalue 0 , and the unpolarized reflectance [M] 00 are all dependent on acquisition geometry. The measurement equation, given in Eq. 8, applied to Eq. 26 yields where a D(0)g = 0.5 assumes that the PSA state a and the PSG state g are both fully-polarized.
If the dominant Mueller-Jones process is known or assumed then only 0 and M 00 remain as unknowns. Consider two noise-free measurements 1 = a 1 Mg 1 and 2 = a 2 Mg 2 using fully-polarized PSA/PSG states. The difference over sum of these two measurements is Here Δ will be zero when the measurements are equal, which indicates 0 = 1/4 and the Mueller matrix is an ideal depolarizer. The other extreme is Δ = ±1, which only occurs when one of the measurements is zero and the Mueller matrix M is non-depolarizing. In the general case, the normalized largest eigenvalue is computed from two measurements by 0 ( , , n) = Both the largest normalized eigenvalue 0 and the most significant Mueller-Jones matrix m 0 depend on measurement geometry ( , , n), as shown in Fig. 4 and Fig. 9 (a). Without any assumptions concerning the eigenspectrum, the eigenvalue can be computed from the Mueller matrix using Eq. 16. The capability to formulate a pBSDF model from a smaller quantity of measurements than required to formulate a Mueller matrix is a way to utilize a triple degeneracy assumption by computing the largest normalized eigenvalue from Eq. 29.

Mueller pBSDF Models
Conventionally, a Mueller pBSDF model's components describe the potential ray paths light may experience in a light-matter interaction. The base model and bulk model both include Mueller matrices which are direct descriptions of ray paths. Fresnel reflection from a microfacet is a rotated Fresnel reflection matrix designated as F . The . indicates a rotation to adjust for the frame of reference changes between the incident plane and exitant plane. Appendix B.3 explains frame of reference rotation. The dependence of rotated Fresnel reflection on acquisition geometry could be made explicit by denoting F 0 , 1 ( , , n), but this dependence is assumed in this section. For brevity, the matrix is written as F 0 , 1 . The base model analyzed in this work is a normalized polarimetric interpretation of the original Cook and Torrence model [7]. This polarimetric interpretation separates the Mueller matrix into a non-depolarizing and a fully depolarizing component, which is an alternate decomposition of a Mueller matrix from a spectral decomposition [27]. Another potential ray path is transmission into a material, depolarizing bulk scatter inside the material, and transmission out of a material; denoted F 1 , 0 D( ) F 0 , 1 . This ray path is considered as a third term in the bulk scattering model p (2) , which is a modified form of the Baek [10] and Kondo models [12]. The complementary model p (0) introduced in this work uses a direct ray path description with the Fresnel reflection F term, but the second term is a complementary matrix Q( f 0 , 1 ) which comes from analysis of the coherency matrix; see Sec.3.5. The GGX microfacet distribution [33] is applied to every model, which adds the fit parameter related to surface roughness. The GGX distribution is a state-of-the-art microfacet distribution used in both the Baek et al. model and the Kondo et al. to describe surface texture effects. Appendix C describes the GGX distribution.
Each pBSDF model type is fit over measurements in ℎ vs. space using least-squares fitting. The least-squares fitting routine aims to minimize the mean squared error of simulated irradiances computed from the polarimetric measurement matrix W over all measurement geometries, Δ(m, p|W) (see Eq. 37).

The Complementary Model
The component Mueller matrices of the complementary model p (0) are the normalized rotated Fresnel reflection matrix f 0 , 1 with the m 00 element set to zero and its complementary Mueller and ( , , n; ) = ( , , n; ) ( , , n; ) Here ( , , n; ) is a scalar-valued function which combines the microfacet distribution function ( , , n; ), the associated shadowing-masking function ( , , ), and other established geometrical factors [7,34]. The fit parameter is albedo-dependent, which is denoted by the subscript. Use of the optional fit parameter is dependent on the selected microfacet distribution function. If a microfacet distribution with a fit parameter is applied, then p (0) is a two parameter model.
The Mueller matrix complementary to Fresnel reflection is denoted Q( f 0 , 1 ); see Eq. 22 for the complementary Mueller matrix definition. The complementary model is a simplification of the coherency matrix assuming a triply degenerate eigenspectrum; see Section 3.5. The most significant Mueller-Jones matrix m 0 is assumed to be Fresnel reflection. Figures 7 (a), (b), and (d) show the Mueller matrices D(0), f 0 , 1 ( , , n), and Q( f 0 , 1 ( , , n)) plotted in signed ℎ versus space.
The Fresnel reflection term F 0 , 1 from a specularly oriented microfacet is weighted using a term which depends on measurement geometry and up to two fit parameters. Fit parameter is albedo-dependent, denoted by the subscript. The fit parameter may be included if the microfacet distribution function selected for ( , , n; ) includes a fit parameter.  [16]. However, since the models in this work are evaluated in a normalized form, the and components can be consolidated into one parameter.
The bulk scattering component traces a path of transmission from air into the media, bulk scattering within the media, and propagation from media back to air. The p (2) model includes three fit parameters: , , and . Parameters denoted with a subscript are wavelength-dependent. This bulk scattering component P is an extension of the diffuse component of the partial Mueller matrix pBSDF model proposed by Baek et al. [10] and used in Kondo et al. [12]. Instead of assuming the bulk scattering inside the material's volume is completely depolarizing, partial isotropic depolarization is allowed. The ratio of partial isotropic depolarization is a fit parameter. As approaches 0, more scattering events take place inside the volume of the material. A value of = 1 is equivalent to the 4 × 4 identity matrix and indicates no scattering events occur. In the case that = 0, the p (2) model is a normalized implementation of the Kondo et al. model with a delay factor equal to zero. In the case that no D(0) term is included and = 0, the p (2) model is a normalized implementation of the Baek model [10].

Measurement Agreement
Two Mueller matrices can be compared by the mean squared error (MSE) of simulated irradiance values from a given polarimetric measurement matrix W where [·] is the ℎ irradiance value. In this work, = 40 and the polarimetric measurement matrix is computed from the PSA/PSG pairs of the RGB950 imaging Mueller Matrix polarimeter [20]. For fitting pBSDF Mueller matrix models, the distance metric in Eq. 36 is averaged over measurement geometries where the incident and exitant propogation directions of the ℎ measurement geometry are [ ] and [ ] , respectively. In this work, = 30 and measurement geometries are reported in Table  2 in Appendix A.

pBSDF Model Fitting Results
Models p (0) , p (1) , and p (2) are evaluated through comparing Δ (m, p|W) , which is the MSE defined in Eq. 37 averaged over all entropies (E). The average MSE for each model over all entropies, over only low entropies, and over only high entropies are given in Tab. 1. Measurement agreement over all 4050 measurements is reported in the Δ column. The Δ (m, p|W) are similar within ±0.002 for the p (0) , p (1) , and p (2)  The addition of more fit parameters and a potential ray path in the p (2) bulk model does not increase polarimetric accuracy for measurements where ≥ 0.565. Larger 0 values means that differences in 1 , 2 , 3 are relatively smaller in magnitude, but Tab. 1 shows that Δ <0.565 is greater than Δ ≥0.565 over all tested models. This is attributable to the relative magnitude elements in a normalized Mueller matrix. For the purpose of analysis in this work, half of all measurements are categorized as high albedo and the other half are categorized as low albedo. Polarization entropies ≥ 0.816 are considered high entropy, with a total of 2029 out of 4050 measurements falling into this category. Low entropy measurements have larger magnitude in the normalized Mueller matrix elements (see: Fig. 3) that produce a larger mean simulated irradiance, Δ .
Measurements with low entropy include low albedo measurements taken at on-specular and near-specular geometries ( ℎ < 10 • ) over all , high albedo measurements of smoother textures (e.g. T1, T2, T3, T4) taken at on-specular and near-specular geometries, and high albedo measurements taken of rougher textures (e.g. T5, T6, T7, T8, T9) at on-specular and near-specular geometries when > 65 • . All other measurements are categorized as high entropy measurements. 0.0080 ± 0.0093 0.0093 ± 0.0097 0.0100 ± 0.0104 Table 1. The average MSE (Eq. 37) over all measurements and over measurements divided into quartiles of entropy ranges. The average MSE are given for the complementary p (0) , base p (1) , and bulk p (2) models. The standard errors for MSE over all measurements and over each quartile are indicated by the ± values. Over all entropies, the complementary, base, and bulk models return approximately similar performance. The bulk model performs best in the lower entropy measurements ( < 0.565), while the complementary models perform best for higher entropy samples ( ≥ 0.565).  Fig. 9. In signed halfway angle ℎ and difference angle space the largest normalized eigenvalue 0 in (a) calculated from Mueller matrix measurements, (b) estimated from the p (0) model using Eq. 37, and (c) estimated from Eq. 39. In (a) 0 has a larger magnitude along the ℎ = 5 • axis than both estimates in (b,c). The same bricks, B1 and R9, are also used to show m 0 in Fig. 4.

Estimating 0
If the exact measured value for 0 is plugged into the complementary model (Eq. 30) in place of the existing F 0 , 1 00 / ( , , n; ) term, then the resulting Δ(m, p (0) |W) is 0.0102. Therefore, this measurement agreement is the best possible performance achievable from assuming triple degeneracy and perfect estimation of 0 . Using the exact measured 0 value in Eq.24, which is an alternate parameterization of the base model p (1) , produces the same measurement agreement. The p (2)  = F 0 , 1 ( , , n) 00 4( h · n) (− · n) ( · n) ( , , n; ) ( , , n; ) . Figure 9 compares this estimated 0 to the measured 0 . The measured eigenvalue 0 has a larger magnitude along the ℎ = 0 axis than the 0 value calculated from p (0) ( , , n). By fitting 0 , the overall Δ is increased to 0.0172. The increase in Δ is not surprising because Δ is not the merit function used in fitting 0 . Fitting a distribution function to the normalized largest eigenvalue minimizes Δ 0 as the merit function. The mean square error between the measured and estimated 0 is where [ 0 ] is the normalized largest eigenvalue of a Mueller matrix at the ℎ geometry. This eigenvalue can be computed from the spectral decomposition of a Mueller matrix, see Eq.19, or estimated from a reduced number of measurements as described in Sec.3.9. Compared to the figure of merit for simulated irradiance values at varying measurement geometries in Eq. 37, the figure of merit for the normalized largest eigenvalue in Eq. 39 contains just the single summation over measurement geometries. For an estimated 0 using the p (0) model, Δ 0 ) = 0.0101. A distribution function that more closely describes 0 may be better suited to modeling opaque plastics than the GGX distribution. Figure 10 demonstrates measurement agreement to the 0 value when using Eq.38 as a model and Eq.39 as a merit function. The plots are arranged in a 3 × 3 grid which indicates the bricks' position in the sample plane. The title of each plot is a key that indicates which color and texture brick is positioned in each location during measurement.

Conclusion
This work applies the spectral decomposition of normalized Mueller matrix measurements to pBSDF modeling. The object ensemble consists of red, orange, yellow, green, and blue bricks roughened using nine textures of sandpaper listed in Tab Spectral decomposition [2] analysis on these 4050 Mueller matrix measurements reveals an approximate triple degeneracy in the smallest three eigenvalues 1 , 2 , and 3 ; see Sec.3.6. When the spectral decomposition is triply degenerate, the pBSDF model's depolarization is completely parameterized by the largest normalized eigenvalue, denoted 0 , which represents the fractional weight of the most significant Mueller-Jones matrix. The most significant Mueller-Jones matrix for these back-scattering measurements is well-approximated by the Fresnel reflection matrix. Figure 4 exemplifies how this Fresnel reflection approximation is applicable over varying albedo and texture. Polarization entropy (Eq. 20), which has an inverse relationship with the largest coherency eigenvalue 0 , increases as surface texture becomes rougher; see Fig. 6. The magnitude of 0 decreases as acquisition geometry moves away from specular orientations, scattering angle approaches normal incidence, and/or albedo increases; see Fig. 5. Triple degeneracy in the eigenspectrum allows the most popular implementation of a back-scattering pBSDF model, the weighted addition of an ideal depolarizer D(0) and Fresnel reflection F 0 , 1 to be represented using only a single, measurable parameter: 0 . The capability to formulate a pBSDF model from a smaller quantity of measurements than required to formulate a Mueller matrix is a way to utilize a triply degenerate assumption by computing the largest normalized eigenvalue from Eq. 29.
This work assesses three normalized models: the p (0) complementary model (Eq.30), the p (1) base model (Eq.32), and the p (2) bulk model (Eq.33). Normalized Mueller matrix measurements m, normalized Mueller pBSDF models p, and the respective polarimetric measurement matrices W are used to evaluate the polarimetric accuracy of a model. Normalized Mueller matrices set the relative irradiance equal to one by dividing each element by the M 00 element. Comparisons using normalized Mueller matrices assess polarimetric accuracy in the non-M 00 elements independently from relative irradiance.
The proposed p (0) complementary model produces similar measurement agreementΔ(m, p|W) as the existing pBSDF models p (1) and p (2) when assessed using the measurement agreement of the simulated normalized irradiance; see Eq.37 forΔ(m, p|W) in Sec.3.9.Δ(m, p (0) |W) is similar toΔ(m, p (1) |W) andΔ(m, p (2) |W) withinΔ(m, p|W) ± 0.0002. Table 1 lists the average measurement agreementΔ over different quadrants of data separated by polarization entropy. For measurements where entropy is ≥ 0.0565, the p (0) proposed model produces the closest measurement agreement by aΔ improvement between 0.0003 and 0.0020. For measurements where entropy is < 0.0565, the p (2) bulk model produces the closest measurement agreement by an improvement toΔ of 0.0045.
One advantage of the proposed complementary model p (0) is the relationship to the measurable quantity 0 . Section 3.9 describes a method for measuring 0 . Another advantage of the complementary model p (0) is a normalized Mueller matrix output which does not require a division by M 00 step. Any normalized Mueller pBSDF model in this work can be easily adapted for applications sensitive to relative irradiance using Eq.6, Eq.7, and any of the many existing non-polarized, scalar BSDF models.
The spectral decomposition analysis that produced the p (0) model also leads to a novel method of fitting a Mueller pBSDF model, described in Sec. 4.2. The second merit function introduced in this work is Δ 0 , the measurement agreement for the largest coherency eigenvalue, 0 ; see Sec.4.2. Equation 31 is one potential distribution function to describe the profile of 0 over varying , , and n. The 0 parameter is useful for Mueller pBSDF modeling because it can either be computed from a full Mueller matrix or directly measured. A full Mueller matrix requires at least 16 polarized measurements, but a direct measurement of 0 requires at least two polarized measurements. Section 3.7 describes a novel method of measuring 0 , assuming the coherency matrix for the light-matter interaction under observation has a triply degenerate eigenspectrum.
For Mueller matrix measurements that are not triply degenerate more than two component Mueller matrices are appropriate. If the most significant Mueller-Jones matrix is not known a priori, then more fit parameters need to be added to the pBSDF model. In this work, these two conditions have been used to demonstrate simplified merit functions for pSBDF fitting and simplified measurements schemes for pBSDF characterization.

A. Geometries
In computer graphics literature, the use of a halfway vector h is common in the implementation of microfacet BSDF models [8,10,11,22,35]. For backscattering events, the halfway vector bisects the incident and reflected ray [22]. An example of the halfway vector is depicted in Figure  2. The vectors and point in the direction of light travel; the halfway vector for reflection is defined as In a system where the surface of the subject under observation is a defined microfacet profile and not a statistical model, a reflection only occurs if the micronormal is parallel to the halfway vector m h [8]. In rendered images, the microfaceted surface profile would be created using a statistical model. The halfway angle ℎ describes the deviation from a specular microfacet from the surface normal calculated as the angle between h and n. The cosine of the halfway angle is therefore equivalent to as cos( ℎ ) = h · n. In the special case where = , the halfway vector is equivalent to .

B.1. Jones Matrix Form
The Fresnel matrices are Muller-Jones matrices which describe the polarization effects at a material interface derived from the Fresnel equations. The Fresnel Mueller matrices can be calculated from the Fresnel Jones matrices using The matrix U is the unitary matrix from Eq. 15. The Jones matrix for Fresnel interactions is typically written as where and are Fresnel amplitude coefficients perpendicular and parallel to the plane of incidence, respectively. These coefficients are expressed as or for reflection and or for transmission. The plane of incidence contains the incident propagation vector and the surface normal n of the material. A rotation into the plane of incidence which contains and the halfway vector h is needed to described the Fresnel reflection from a microfacet.
In systems where the x-axis correspond to s-polarized light, Eq. 42 is used. In systems where the y-axis corresponds to s-polarized light, the positions of the s and p polarization in Eq. 42 are swapped. The measurements in this work result in a global coordinate definition where the vertical y-axis corresponds to the orientation of s-polarized light; this geometrical difference is accounted for in model implementation to ensure a geometrical match with the measurement setup.

B.2. Mueller Matrix Form
The polarization effects of a specular reflection from a microfacet are described by the nondepolarizing Fresnel reflection matrix (43) An asterisk indicates a complex conjugate and Eq. 42 has been converted from a Jones to a Mueller matrix. The Fresnel reflection coefficients and are and ( 0 , 1 , ) = 1 cos( ) − 0 cos( ) where 0 is the starting refractive index, 1 is the ending refractive index, is the incident angle upon the material microfacet, and is calculated using Snell's law: 0 sin ( ) = 1 sin ( ). For a material which is illuminated by light which propagates through air first, 0 = 1 while 1 > 1. The polarization effects of light passing from inside the media back to air toward the camera are described by the non-depolarizing Fresnel transmission matrix The Fresnel transmission coefficients and are calculated using Eq. 47 and Eq. 48.
In Eq. 47 and 48, 0 is the ambient material's refractive index (e.g. air) and 1 is the observed material refractive index. The angle is the incident angle upon the boundary between media from the 0 side, and = arcsin( 0 sin( )/ 1 ) as calculated from Snell's Law. For transmission from air into a non-air material, 0 = 1 and 1 > 1. For transmission from a non-air material out to air, 0 > 1 and 1 = 1.

B.3. Mueller Matrix Rotation
A polarization ellipse is specified by two orthogonal directions in the transverse plane [24]. A rotation of this coordinate system, which preserves the transverse plane but changes the local coordinate system is described by Eq. 49.
where is the angle between two coordinate systems. This matrix rotates positive angles in the counterclockwise direction about the propagation direction. The local coordinate system of a Mueller matrix can be specified using the unit vectors , and , which are adapted from Priest & Germer [23]. , corresponds to the macrosurface while , corresponds to the microsurface. The subscripts or indicate whether the vectors correspond to incident or exitant directions of travel for a light-matter interaction.
The unit vector = ( × h)/| × h| is perpendicular to the photon travel direction and the halfway vector h; which is to say is the unit vector perpendicular to the micro-incident plane formed by and h. The unit vector = ( × n)/| × n| is perpendicular to a photon travel direction and the macronormal n. The vector is perpendicular to the macro-incident plane formed by and n. The angle from to rotates the coordinate system from the macro-incident plane to the plane of incidence for a Fresnel microfacet, i.e. the micro-incident plane [23]. The four-quadrant tangent function is used to calculate the incident rotation angle [9] tan ( ) = ( × n) · h ( h · n) − ( · h) ( · n) where identities (a × b) · (c × d) = (a · c) (b · d) − (a · d) (b · c) and (a × b) × (a × c) = ((a × b) · c)a have been used. In principle, Eq. 50 should be equal to | × | · = sin ( ) cos ( ) , but in practice this has led to rotation angles constrained between 0 • and 90 • .
A similar rotation is needed for the exitant light. The unit vector = −(| × n)/| × n| is perpendicular to a photon travel direction and the macronormal n; the vector is perpendicular to the macro-exitant plane. The unit vector = −( × h)/| × h| is perpendicular to the photon travel direction and the halfway vector h; the vector is perpendicular to the macro-exitant plane [23]. The angle is . (51)

B.4. Rotated Fresnel Matrices
The Fresnel matrices in Eq. 43 and Eq. 46 are applied to light-material interactions in the plane of incidence defined by the incident propagation vector and the micronormal, here called the micro-incident plane. The Fresnel matrix F is rotated from the macro-incident plane by to the micro-incident plane, where the Fresnel matrix is applied, and then to the macro-exitant plane by where is the angle from the macro-incident plane to the micro-incident plane and is the angle from micro-incident plane to the macro-exitant plane. Eq. 50 and Eq. 51 in Appendix A are used to calculate the transverse rotation angles.

C. GGX Microfacet Distribution Function
Microfacet distribution functions ( , , n; ) are probability distribution functions used to describe how much light scatters after a material interaction from a given source position , exitant propagation direction , and (optionally) surface roughness . Shadowing-masking functions ( , , n; ) are optionally applied to BSDFs to adjust for the effects of steep microfacet orientations that shadow neighboring microfacets. This work applies the GGX microfacet distribution function to the Fresnel reflection matrix.
Microfacet distribution functions ( , , n; ) are distinct from microfacet response functions ( , , n; ) which are sometimes used in other works to describe a statistical distribution of surface normals over a surface [8,10]. Microfacet distribution functions ( , , n; ) are preferrentially used in this work for being probability distribution functions which integrate to 1. They are generally related through the relationship ( , , n; ) = ( h · n) ( , , n; ).
The GGX microfacet function [8,33] is designed to be applied to both reflection and transmission events with a consideration for conservation of radiance included for the case of transmission. The GGX distribution was developed by fitting reflection measurements from and transmission measurements through different finishes of roughened glass [8]. Eq. 53 and Eq. 54 describe the GGX distribution function and its associated shadowing-masking function, respectively. The GGX distribution includes an fit parameter which tunes the shape of the microfacet distribution curve. Not all microfacet distribution functions include a parameter.  Table 2 reports the position of the illumination source and camera as and , respectively. The center of the illumination source and camera's optical axis are the same height as the center of the sample plane. The angles of incidence and angles of observation quoted in Tab. 2 are in reference the the center brick in the 3x3 LEGO tower arrangement (Fig. 1). Exact regions of interest in the other eight off-center bricks will have and values which vary up to ±8 • depending on the measurement geometry under observation.