Estimating surface orientation from microfacet Mueller matrix bidirectional reflectance distribution function models in outdoor passive imaging polarimetry

Abstract. Representative examples from 3 years of measurements from JPL’s ground-based multiangle spectropolarimetric imager (GroundMSPI) are compared to a Mueller matrix bidirectional reflectance distribution function (mmBRDF). This mmBRDF is used to model polarized light scattering from solar illuminated surfaces. The camera uses a photoelastic-modulator-based polarimetric imaging technique to measure linear Stokes parameters in three wavebands (470, 660, and 865 nm) with a ±0.005 uncertainty in degree of linear polarization. GroundMSPI measurements are made over a range of scattering angles determined from a fixed viewing geometry and varying sun positions over time. This microfacet mmBRDF model predicts an angle of the linear polarization that is consistently perpendicular to the scattering plane and therefore is only appropriate for rough surface types. The model is comprised of a volumetric reflection term plus a specular reflection term of Fresnel-reflecting microfacets. The following modifications to this mmBRDF model are evaluated: an apodizing shadowing function, a Bréon or Gaussian microfacet scattering density function, and treating the surface orientation as an additional model parameter in the specular reflection term. The root-mean-square error (RMSE) between the GroundMSPI measurements and these various forms of the microfacet mmBRDF model is reported. Four example scenes for which a shadowed-Bréon microfacet mmBRDF model yields realistic estimates of surface orientation, and the lowest RMSE among other model options are shown.


Introduction
Ground-based sun photometers do not provide adequate spatial or optical sampling to constrain the inverse problem of global aerosol estimation. 1 A downward-looking polarimeter collects light from both the atmosphere and the Earth's surface and can provide information for spatiotemporal aerosol estimation on a global scale. 2 Many atmospheric imaging studies are conducted over the ocean to avoid optical reflectance from an unknown landscape. 3 The assumption that surface reflection properties are known is not valid for aerosol retrievals over land. If the surface reflection properties are unknown, only multiple-viewing-angle measurements of both intensity and polarization are able to provide the relevant aerosol parameters with sufficient accuracy for climate research. 4,5 Currently, models of atmospheric scattering are more advanced than models of surface reflectance because the polarized light scattering from the Earth's surface is globally diverse. Analysis to date has been semiempirical 6 since entirely physical models of the Earth's polarized light scattering fail to correspond with measurements. 7 Therefore, semiempirical models are derived from extensive measurements and include a small number of parameters that are fit to measurements. 8,9 To extract the optical and microphysical properties of aerosols over urban areas, it is essential to have more accurate models for polarized surface reflectance coupled to atmospheric radiative transfer to allow the simultaneous retrieval of aerosol and surface properties. 10,11 For example, sensitivity studies have reported that an error on the order of 10 −3 for the surface polarized reflectance can lead to an error in aerosol optical thickness (AOT) of 0.04. 7,12 This work employs microfacet mmBRDF modeling methods used in the initial analysis of GroundMSPI measurements. 13 This work builds upon those results by reporting results on a larger diversity of scenes, augmenting the model with a shadowing function, considering two different scattering density functions, and including an estimate of the objects' surface orientation. The model predicts a polarization orientation that is invariant to scattering angle and 90 deg from the scattering plane for unpolarized illumination and materials of real-valued index. GroundMSPI measurements for which this polarization orientation is observed are usually rough and natural materials: e.g., grass, foliage, asphalt, and soil. A root-mean-square error (RMSE) comparison between microfacet mmBRDF models and model variants are reported for these GroundMSPI measurements.

Methods
A Stokes vector uniquely quantifies the polarization state of a light field using four numbers s ¼ ½I; Q; U; V. 14 In this paper, lowercase bold letters are used to denote vectors and uppercase bold letters are used for matrices. The MSPI camera architecture uses a time-varying retardance in the optical path to modulate the orientation of the linearly polarized component of the incoming light, described by the Stokes components Q (excess of horizontally over vertically oriented polarized light) and U (excess of 45 deg over 135 deg oriented polarized light). 15,16 These orientations associate the Stokes vectors with a coordinate system; usually either the scattering or meridian planes. In this coordinate system, a polarization orientation 90 deg from the scattering plane results in U ¼ 0. Circular polarization V is typically negligible for sunlight scattered by the atmosphere or Earth's surface 17 and therefore not measured by GroundMSPI. The units of the linear Stokes components are radiance but can alternatively be expressed as a dimensionless quantity, the degree of linear polarization (DoLP) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 5 7 6 DoLP ¼ where I is the total radiance. The angle of the linear polarization (AoLP) is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 6 3 ; 5 1 0 AoLP ¼ and is undefined DoLP ¼ 0. Note that DoLP and radiance measurements are invariant to the coordinate system of Q and U. AoLP, Q, and U depend on the coordinate system. For the remainder of the paper, AoLP s ∕Q s ∕U s and AoLP m ∕Q m ∕U m will be used to denote quantities with respect to the scattering (s) or meridian (m) planes. The meridian plane contains the global horizontal surface normal and the view vector. The solar illumination vector (i.e., the direction between the sun and the object) is denoted asr i ¼ ½cos ϕ i sin θ i ; − sin ϕ i sin θ i ; − cos θ i where the hat denotes a unit vector. The scattering plane contains the solar illumination vector, and the view vector that describes the direction that MSPI is pointed and is denoted byr r ¼ ½cos ϕ r sin θ r ; − sin ϕ r sin θ r ; − cos θ r . Our conventions are vectors point in the direction of photon travel, polar angles are measured from nadir, and azimuth angles are measured clockwise from north, see Fig. 1. The variable Δϕ ¼ ϕ i − ϕ r is used when only the difference in azimuth angles has an effect. The scattering angle Ω is defined as the angle between the view and illumination vectors E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 6 3 ; 2 4 5 cos Ω ¼r i ·r r ¼ μ i μ r þ sin θ r sin θ i cosðΔϕÞ; (3) where μ i ¼ cos θ i , μ r ¼ cos θ r , and Δϕ ¼ ϕ r − ϕ i . Given the convention of view and illumination vectors pointing in the direction of photon travel Ω ¼ 180 deg indicates complete backscattering. The results of this paper are reported in Sec. 3, where the scattering angle is used to create two-dimensional graphs comparing measurements and model at various three-dimensional acquisition geometries. GroundMSPI measurements are acquired in time-sequences throughout the day; details in Table 1. Note that the viewing geometry, i.e.,r r , is constant during this acquisition. Therefore, the scattering angles reported in this paper are varied by changes to the illumination direction, i.e.,r i . Figure 2 shows an example GroundMSPI measurement of a cowboy boot sculpture made of concrete. Figure 2 (Video 1) is a static image and a time-sequence from 08:39 to 16:37 MST is available as supplementary material. In Sec. 3, the same measurements are reported in Fig. 3(d) range where the range of scattering angles is 60 deg to 120 deg. Figure 5(c) offers another depiction of the acquisition geometry where the blue vector is r r , and the illumination vector is r i .

Microfacet Model for the mmBRDF
When light is incident upon a surface, the light-matter interaction can be quantified as a change to the Stokes vector. Linear effects are conventionally described by a 4 × 4 matrix, which relates an incident Stokes vector to a reflected Stokes vector where F λ is called the bidirectional reflectance distribution matrix (mmBRDF), which is in general dependent upon the wavelength λ of the incident s 0 and reflected s r fields; the transpose is denoted by †. The scalar-valued counterpart of the mmBRDF describes nonpolarimetric light scattering [i.e., the first component of the vectors in Eq. (4)] and is  called the bidirectional reflectance distribution function (BRDF). Some authors relate I pol of the incident to reflected fields using a bidirectional polarized reflectance distribution function (BPDF). In general, the mmBRDF, BPDF, and BRDF are all dependent upon the geometric configuration between the camera, the illumination, and the object. The microfacet mmBRDF model used in the work is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 6 3 ; 6 7 5 where F λ is a 4 × 4 matrix and a λ , k, b, and ξ are four scalar-valued dependent variables. The independent variables describing the acquisition geometry are μ i , μ r , and Δϕ. The first term in Eq. (5) describes the unpolarized portion of the total reflected field and the second term describes the polarized portion. Therefore, the variables a λ , k, b describe the unpolarized field and ξ describes the relative intensity of the polarized radiation. In other words, the mmBRDF in Eq. (5) specifies the form of the polarized light scattering as a function of the acquisition geometry. Agreement between this model and GroundMSPI measurements is based only on adjustments to the relative intensity of the polarized radiance through the parameter ξ. The polarized light scattering from specularly reflecting microfacets is described by the function prðβ; ΨÞ. Variations of this mmBRDF model introduce more dependent variables to the polarized term by parameterizing prðβ; ΨÞ. For example, the surface normal of the object is estimated by relating prðβ; ΨÞ to the acquisition geometry.
Microfacets are perfect Fresnel reflecting elements that are small compared to the imaging system resolution. 19 The magnitude of the total reflected polarized light is proportional to the quantity of microfacets that are oriented toward specular reflection at a given acquisition geometry. At a given acquisition geometry, the illumination vector and view vector are bisected by a surface normaln, which reflects specularly. The angle β is betweenn and the vector normal to the surface of the objectr s . If at a given acquisition geometry and surface orientation, there are many microfacets in the β orientation, then the specular reflection will be strong. Conversely, if at that value of β the surface does not have many microfacets the specular reflection is weak. The microfacet structure of a surface is described by a scattering density function prðβ; ΨÞ, which parameterizes the density of microfacets at a polar and azimuth angle; denoted β and Ψ, respectively. The orientation of microfacets is usually assumed to be azimuthally uniform so that prðβ; ΨÞ ¼ prðβÞ. The assumption of azimuthal symmetry imposes an invariance in the polarized light scattering when rotating the source and the view direction together about the surface normal. For top of atmosphere measurements, the surface normal to the land surface can be reasonably approximated by the zenithẑ. The surface normal of each object in a GroundMSPI image is not known a prior. In this work, two approaches to defining β are compared: (1) with respect to the zenithẑ (as in Fig. 3) and (2) with respect to the object's surface normal vectorr s , which is not known a prior and is estimated from the measurements (as in Figs. 4 and 5). The variabler s is inserted into Eq. (5) using the substitutions E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 3 2 6 ; 7 5 2 cosðβÞ ¼r s ·n; μ i ¼ −r s ·r i ; μ r ¼r s ·r r ; wherer s is a vector normal to the surface. The first term in Eq. (5) models the unpolarized portion of the reflected field. Here, D is a 4 × 4 matrix that is zero except for the (1,1) element. Without the matrix D, the first term in Eq. (5) reduces to a scalar-valued BRDF model referred to as a modified Rahman-Pinty-Verstraete (mRPV) that is currently used in aerosol retrievals by NASA's MISR instrument. 20 The mRPV model uses three parameters to describe a nonpolarimetric intensity measurement: a λ to characterize the amount of reflectance at each wavelength, k to characterize the anisotropy of reflectance, and b to characterize the forward/backscattering contribution. 21,22 The second term in Eq. (5) models polarization. Pðn; γÞ is a 4 × 4 matrix of Fresnel reflection coefficients of a material where n is the index of refraction and γ is the angle of incidence. A wavelength subscript on Pðn; γÞ is not used due to the assumption that the refractive index of the surface is spectrally neutral within the measured waveband. Since specular reflection from each microfacet is assumed the scattering angle is related to the angle of incidence by cos Ω ¼ − cos 2γ. Details about the structure and elements of the matrix Pðn; γÞ, for applications in GroundMSPI analysis, have been described in Ref. 1. This model assumes all reflected polarized radiance is generated by a single reflection of light from the top surface. The matrix Mðα 0 Þ rotates an incident Stokes vector from the solar meridian plane to the scattering plane where α 0 is the angle between the two planes. If results are desired in meridian, rather than scattering, coordinates then another rotation matrix can be inserted on the left-hand side of Pðn; γÞ. The second term in Eq. (5) models the polarimetric contribution to the microfacet mmBRDF using only a single linear parameter: ξ. The model describes the functional relationship between the linear Stokes parameters and the acquisition geometry and fitting ξ to the polarimetric measurements simply scales this functional form. This simplification is useful since polarimetric measurements show a similar trend with respect to scattering angle for rough materials and the Earth's natural landscapes. 7 Shadowing effects manifest themselves both in photometric and polarimetric characteristics of scattered radiation. 23 Without a shadowing function, microfacet mmBRDF models assume that every element of the surface contributes to the reflected field. This assumption neglects the shadowing of points on the surface caused by neighboring points, an effect that is increasingly important at large angles of incidence. A shadowing theory described in Refs. 24 and 25 derives the probability that a backscattering specular point on the surface will not be shadowed for a Gaussian slope density function. Analytic results for a non-Gaussian example are presented in Ref. 26. The slope density function is related to the scattering density function prðβÞ from which a shadowing function Sðμ i ; μ r Þ can be calculated so that the shadow-modified scattering density is prðβÞSðμ i ; μ r Þ. Both the Bréon 27,28 and Gaussian scattering density functions are used in this work and reported in Table 2. Measurement comparisons to these models will be computed both with and without a shadowing function.
Note that the Bréon shadowing function only depends on the acquisition geometry. The Gaussian shadowing function depends on both the acquisition geometry and the width of prðβÞ, denoted as σ. This width parameter of the Gaussian model is an additional dependent parameter as compared to the Bréon model. For empirical analysis, an optimal model imposes the strictest constraints possible while characterizing a realistic population of measurements. In this application, mmBRDF constraints are valuable for untangling land surface and atmospheric polarized light scattering for the task of aerosol retrieval. Since the Bréon model imposes a more strict constraint on the surface polarized light scattering, it would be preferable to the Gaussian model in the case where equivalent RMSE values were observed.
A Gaussian scattering density function is also called the Cox-Munk distribution in reference to the first researchers to use this model to describe ocean roughness. 29 Where, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 3 2 6 ; 6 4 2 Derivations of Eqs. (9) and (11) can be found in Refs. 24, 26, and 28.
Recent work on microfacet modeling for polarized light scattering has demonstrated that shadowing functions introduce errors at large grazing angles, i.e., μ i , μ r → 0. 30 In this paper, only GroundMSPI measurements that are directly solar illuminated are considered. Therefore, the range of reported acquisition geometries does not include grazing angles because of the long shadows physical objects cast on each other when the sun is near the horizon.

Assumptions for Passive Outdoor Scenes
Only GroundMSPI measurements at 660 nm are included in this work given that the spectral invariance hypothesis is well supported for polarized light scattering from the Earth's surface. GroundMSPI measurements have been analyzed with respect to the spectral invariance hypothesis in other work. 1,31 The spectral invariance hypothesis assumes that both the magnitude and the angular shape of the polarimetric term in the mmBRDF are spectrally neutral. Researchers have noted only small exceptions, such as wavelength dependencies in the shortwave infrared (SWIR) at large scattering angles, that could affect aerosol retrievals over land when the AOT is small. 32 A conventional approach is to fit the parameters of a mmBRDF model for land surface reflectance in the SWIR where, in general, the atmospheric contribution is small. Since the polarized surface reflectance depends weakly on the wavelength, the fitted values for the parameters are used at other wavelengths to predict the contribution of polarized light scattering from the land surface.  We will restrict the modeling to unpolarized solar illumination [i.e., s 0 ¼ ½1;0; 0;0 in Eq. (4)]; this precludes scattered skylight that is polarized. In the first GroundMSPI paper, 1 modeling solar illumination as partially polarized did not significantly improve the model and measurement agreement for data collected on a clear day. When the index of refraction is real-valued, the Fresnel reflection coefficients are also real-valued. In this work, n ¼ 1.5 is assumed. Putting together the two assumptions of a real index of refraction and unpolarized illumination leads to an important constraint on the model for the reflected Stokes vector: the third component is zero (i.e., U s ¼ 0) and therefore AoLP s ¼ 90 deg at all acquisition geometries. Other researchers also use this assumption. 7 Others have also noted AoLP s ≠ 90 deg in their measurements of land surface scattering with errors up to 15%. 33 Following methods from Ref. 1, the polarimetric and nonpolarimetric parameters are estimated sequentially. This is because specular reflection from the microfacets contributes to the total radiance field, as well as the polarized field, but the polarized field is unaffected by the remaining parameters. First ξ (and σ if using a Gaussian facet distribution model) are estimated from a least-squares (LS) fit to Q s measured at a series of different scattering angles. These estimates parameterize an estimate of the first term in Eq. (5), which is then subtracted then from the measured radiance and an LS fit to the first term in Eq. (5) (i.e., a λ , k, and b) is performed.
Although the surface orientation affects both the nonpolarimetric and the polarimetric term in this work, we choose to estimater s from an LS fit of the polarimetric term (i.e., simultaneous with estimating ξ) and then use that estimate when fitting the other parameters in the nonpolarimetric term. Figure 3 shows four GroundMSPI intensity images alongside double y-axis graphs of linear Stokes measurements and the predictions from four models: Bréon and Gaussian; each with and without shadowing [see Eq. (7)]. Here, the surface orientation is assumed to be horizontal and therefore the surface normal vector equals the zenith. The four RMSE values of the models for each scene are in Table 3.

Assumed Surface Orientation
In Fig. 3(a), a large region of gravel is selected and the average linear Stokes values in this region are fit to the microfacet mmBRDF models. For this example, the assumed horizontal surface orientation is correct. The Bréon models, both with and without shadowing, produce lower RMSE than the Gaussian models. The shadowed-Bréon model RMSE is 17%, and the unshadowed-Bréon model RMSE is 22%. In Fig. 3(b), a region on the wall of a brick building is selected. The true surface orientation is vertical yet the model assumes horizontal. Even with this modeling mismatch, all four microfacet mmBRDF models trend correctly with the measurement dependence on scattering angle. RMSE of unshadowed-Bréon and unshadowed-Gaussian models are both 12%, adding a shadowing function increases RMSE in both cases. The region used for analysis in Fig. 3(c) is the leafy part of a tree and RMSE of three of the models is 22 AE 1%; adding shadowing to the Bréon model reduces the RMSE to 9%. Measurements in Fig. 3(d) are of smooth cement and the Bréon model actually trends nearly opposite of the measurements; the Gaussian RMSE is more than four times lower. As shown in Table 3, shadowing does not improve the RMSE for the Gaussian distribution. In summary, the model with the lowest RMSE is the shadowed-Bréon for the tree and gravel measurements, the unshadowed-Gaussian for the smooth cement measurement, and a tie between these two models for the brick measurement.

Estimated Surface Orientation
Here two angular parameters to describe the object's surface normal are included in the LS parametric fit, see Figs. 4 and 5. This surface normal isr s defined in Eq. (6). As reported in Table 4, all RMSE values are below 8% and most RMSE values are below 5% when the object's surface normal is estimated as part of the microfacet mmBRDF model. These RMSE differences are barely noticeable in the Stokes plots in Figs. 4 and 5. The difference in the estimated surface normals is very noticeable in Figs. 4 and 5. The four GroundMSPI measurements are split into the two figures: one for which the true value of the surface normal is known in Fig. 4, and the other for which the surface normal is unknown or undefined in Fig. 5.
In Fig. 5(a), the true-value of the surface orientation is undefined since the selected pixels in the tree contain numerous leaves, shadows, and even the background seen through the tree branches. The purpose of including this example is to demonstrate the improvement in RMSE even when the surface orientation parameter is undefined. In Fig. 5(c), the smooth cement is on a south-facing vertical surface that is angled toward the sky. Both the shadowed-Gaussian and shadowed-Bréon are south-facing. The shadowed-Bréon is close to perfectly vertical. The shadowed-Gaussian is angled Table 3 RMSE values between GroundMSPI measurements and models in Fig. 3.   close to a 45 deg angle, which appears steeper than the truevalue. The true-value of the surface orientation is somewhere between the shadowed-Gaussian and shadowed-Bréon surface orientation estimates. As reported in Table 4, the addition of a shadowing function perturbs the RMSE values by a few percent, both increasing and decreasing depending upon the sample. The Gaussian and Bréon RMSE values are very similar and the rank-ordering of model RMSE is different for each of the four examples. For gravel shadowed-Bréon and shadowed-Gaussian are tied at 1%. For the brick and tree examples, the unshadowed-Bréon and unshadowed-Gaussian are tied at 3% and 7%, respectively. The Gaussian model RMSE is 2% for the smooth cement example, both shadowed and unshadowed.

Conclusions
We applied modeling methods used in the initial analysis of GroundMSPI measurements 1 to four example datasets collected at the University of Arizona. The microfacet mmBRDF models for this paper were modified and RMSE comparisons were made for shadowed-and unshadowed-Gaussian and shadowed-and unshadowed-Bréon scattering density functions. The RMSE for both an assumed value and an estimated value of the objects' surface normal were compared. To estimate the surface normal, two additional parameters were added to the model and this improved the RMSE by a factor of 5 to 6. For the two examples where the true value of the surface normal is well-defined and known the estimated surface normal is most accurate for the shadowed-Bréon model. This analysis is preliminary. Future work will include a greater variety of materials and acquisition geometries from more University of Arizona GroundMSPI datasets to characterize the utility of microfacet mmBRDF models in aerosol retrieval. Also, other popular mmBRDF models and shadowing functions will be included in future RMSE comparisons. Incorporating prior knowledge and constraints to estimates of surface orientation will be used to improve the microfacet scattering density functions. The statistical significance of an additional parameter in the Gaussian scattering density function, as compared to the one parameter Bréon model, will be accessed using stepwise regression and F-tests on the error.