A Geometrical Optics Polarimetric Bidirectional Reflectance Distribution Function for Dielectric and Metallic Surfaces

A polarimetric bidirectional reflectance distribution function (pBRDF), based on geometrical optics, is presented. The pBRDF incorporates a visibility (shadowing/masking) function and a Lambertian (diffuse) component which distinguishes it from other geometrical optics pBRDFs in literature. It is shown that these additions keep the pBRDF bounded (and thus a more realistic physical model) as the angle of incidence or observation approaches grazing and better able to model the behavior of light scattered from rough, reflective surfaces. In this paper, the theoretical development of the pBRDF is shown and discussed. Simulation results of a rough, perfect reflecting surface obtained using an exact, electromagnetic solution and experimental Mueller matrix results of two, rough metallic samples are presented to validate the pBRDF. © 2009 Optical Society of America OCIS codes: (290.1483) BSDF, BRDF, and BTDF; (290.5855) Scattering, polarization; (290.5880) Scattering, rough surfaces; (160.3900) Metals References and links 1. D. E. Barrick, “Theory of HF and VHF propagation across the rough sea—parts I and II,” Radio Sci. 6, 517–533 (1971). 2. C. Eckart, “The scattering of sound from the sea surface,” J. Acoust. Soc. Am. 25, 566–570 (1953). 3. E. Y. Harper and F. M. Labianca, “Scattering of sound from a point source by a rough surface progressing over an isovelocity ocean,” J. Acoust. Soc. Am. 58(2), 349–364 (1975). 4. K. Krishen, “Correlation of radar backscattering cross sections with ocean wave height and wind velocity,” J. Geophys. Res. 76, 6528–6539 (1971). 5. B. W. Hapke, “A theoretical photometric function for the lunar surface,” J. Geophys. Res. 68(15), 4571–4586 (1963). 6. D. S. Kimes, “Modeling the directional reflectance from complete homogeneous vegetation canopies with various leaf-orientation distributions,” J. Opt. Soc. Am. A 1(7), 725–737 (1984). 7. R. Hegedüs, A. Barta, B. Bernáth, V. B. Meyer-Rochow, and G. Horváth, “Imaging polarimetry of forest canopies: how the azimuth direction of the sun, occluded by vegetation, can be assessed from the polarization pattern of the sunlit foliage,” Appl. Opt. 46(23), 6019–6032 (2007). 8. G. Zonios, I. Bassukas, and A. Dimou, “Comparative evaluation of two simple diffuse reflectance models for biological tissue applications,” Appl. Opt. 47(27), 4965–4973 (2008). 9. J. Xia and G. Yao, “Angular distribution of diffuse reflectance in biological tissue,” Appl. Opt. 46(26), 6552–6560 (2007). 10. J. F. Blinn, “Models of light reflection for computer synthesized pictures,” in SIGGRAPH 1977 Proceedings, vol. 11, pp. 192–198, Special Interest Group on Graphics and Interactive Techniques (Computer Graphics, 1977). #118213 $15.00 USD Received 6 Oct 2009; accepted 5 Nov 2009; published 18 Nov 2009 (C) 2009 OSA 23 November 2009 / Vol. 17, No. 24 / OPTICS EXPRESS 22138 11. R. L. Cook and K. E. Torrance, “A reflectance model for computer graphics,” in SIGGRAPH 1981 Proceedings, vol. 15, pp. 307–316, Special Interest Group on Graphics and Interactive Techniques (Computer Graphics, 1981). 12. X. D. He, K. E. Torrance, F. X. Sillion, and D. P. Greenberg, “A comprehensive physical model for light reflection,” in SIGGRAPH 1991 Proceedings, vol. 25, pp. 175–186, Special Interest Group on Graphics and Interactive Techniques (Computer Graphics, 1991). 13. M. A. Greiner, B. D. Duncan, and M. P. Dierking, “Bidirectional scattering distribution functions of maple and cottonwood leaves,” Appl. Opt. 46(25), 6485–6494 (2007). 14. P. Y. Ufimtsev, Fundamentals of the Physical Theory of Diffraction (John Wiley & Sons, Inc., Hoboken, NJ, 2007). 15. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge University Press, New York, NY, 1999). 16. P. Beckmann and A. Spizzichino, The Scattering of Electromagnetic Waves from Rough Surfaces (Artech House, Inc., Norwood, MA, 1963). 17. A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, New York, NY, 1997). 18. C.-H. An and K. J. Zeringue, “Polarization scattering from rough surfaces based on the vector Kirchoff diffraction model,” in Proc. SPIE, vol. 5158, pp. 205–216 (The International Society for Optical Engineering (SPIE), 2003). 19. D. A. McNamara, C. W. I. Pistorius, and J. A. G. Malherbe, Introduction to the Uniform Geometrical Theory of Diffraction (Artech House, Inc., Norwood, MA, 1990). 20. K. E. Torrance and E. M. Sparrow, “Theory for off-specular reflection from roughened surfaces,” J. Opt. Soc. Am. 57(9), 1105–1114 (1967). 21. B. P. Sandford and D. C. Robertson, “Infrared reflectance properties of aircraft paints,” in Proceedings of IRIS Targets, Backgrounds and Discrimination (1985). 22. M. P. Fetrow, D. Wellems, S. H. Sposato, K. P. Bishop, T. R. Caudill, M. L. Davis, and E. R. Simrell, “Results of a new polarization simulation,” in Proc. SPIE, vol. 4481, pp. 149–162 (The International Society for Optical Engineering (SPIE), 2002). 23. R. G. Priest and T. A. Germer, “Polarimetric BRDF in the microfacet model: theory and measurements,” in Proceedings of the 2000 Meeting of the Military Sensing Symposia Specialty Group on Passive Sensors, pp. 169–181 (Infrared Information Analysis Center, 2000). 24. R. G. Priest and S. R. Meier, “Polarimetric microfacet scattering theory with applications to absorptive and reflective surfaces,” Opt. Eng. 41(5), 988–993 (2002). 25. D. Wellems, S. Ortega, D. Bowers, J. Boger, and M. Fetrow, “Long wave infrared polarimetric model: theory, measurements and parameters,” J. Opt. A: Pure Appl. Opt. 8(10), 914–925 (2006). 26. D. Wellems, M. Serna, S. H. Sposato, M. P. Fetrow, K. P. Bishop, S. A. Arko, and T. R. Caudill, “Spectral polarimetric BRDF model and comparison to measurements from isotropic roughened glass,” in Workshop on Multi/Hyperspectral Sensors, Measurements, Modeling and Simulation (U.S. Army Aviation and Missile Command, Huntsville, AL, 2000). 27. K. K. Ellis, “Polarimetric bidirectional reflectance distribution function of glossy coatings,” J. Opt. Soc. Am. A 13(8), 1758–1762 (1996). 28. H. G. Tompkins and E. A. Irene, Handbook of Ellipsometry (William Andrew, Inc., Norwich, NY, 2005). 29. R. Anderson, “Matrix description of radiometric quantities,” Appl. Opt. 30(7), 858–867 (1991). 30. D. S. Flynn and C. Alexander, “Polarized surface scattering expressed in terms of a bidirectional reflectance distribution function,” Opt. Eng. 34(6), 1646–1650 (1995). 31. F. E. Nicodemus, “Radiance,” Am. J. Phys. 31, 368–377 (1963). 32. F. E. Nicodemus, “Directional reflectance and emissivity of an opaque surface,” Appl. Opt. 4(7), 368–377 (1965). 33. J. R. Schott, Fundamentals of Polarimetric Remote Sensing (SPIE Press, Bellingham, WA, 2009). 34. J. R. Shell, “Polarimetric Remote Sensing in the Visible to Near Infrared,” Ph.D. dissertation, Chester F. Carslon Center for Imaging Science, Rochester Institute of Technology, Rochester, NY (2005). 35. Y. Sun, “Statistical ray method for deriving reflection models of rough surfaces,” J. Opt. Soc. Am. A 24(3), 724–744 (2007). 36. W. S. Bickel and W. M. Bailey, “Stokes vectors, Mueller matrices, and polarized scattered light,” Am. J. Phys. 53(5), 468–478 (1985). 37. M. G. Gartley, S. D. Brown, and J. R. Schott, “Micro-scale surface and contaminate modeling for polarimetric signature prediction,” in Proc. SPIE, vol. 6972 (The International Society for Optical Engineering (SPIE), 2008). 38. J. R. Maxwell, J. Beard, S. Weiner, D. Ladd, and S. Ladd, “Bidirectional Reflectance Model Validation and Utilization,” Tech. Rep. AFAL-TR-73-303, Air Force Avionics Laboratory, Wright-Patterson Air Force Base, OH (1973). 39. M. G. Gartley, “Polarimetric Modeling of Remotely Sensed Scenes in the Thermal Infrared,” Ph.D. dissertation, Chester F. Carslon Center for Imaging Science, Rochester Institute of Technology, Rochester, NY (2007). 40. A. F. Peterson, S. L. Ray, and R. Mittra, Computational Methods for Electromagnetics (IEEE Press, New York, NY, 1998). 41. R. M. Axline and A. K. Fung, “Numerical computation of scattering from a perfectly conducting random surface,” IEEE Trans. Antennas Propag. AP-26(3), 482–488 (1978). 42. E. I. Thorsos, “The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian rough#118213 $15.00 USD Received 6 Oct 2009; accepted 5 Nov 2009; published 18 Nov 2009 (C) 2009 OSA 23 November 2009 / Vol. 17, No. 24 / OPTICS EXPRESS 22139 ness spectrum,” J. Acoust. Soc. Am. 83(1), 78–92 (1988). 43. A. K. Fung and M. F. Chen, “Numerical simulation of scattering from simple and composite random surfaces,” J. Opt. Soc. Am. A 2(12), 2274–2284 (1985). 44. M. F. Chen and S. Y. Bai, “Computer simulation of wave scattering from a dielectric random surface in two dimensions—cylindrical case,” J. Electromagn. Waves Appl. 4(10), 963–982 (1990). 45. E. Compain, S. Poirier, and B. Drevillon, “General and self-consistent method for the calibration of polarization modulators, polarimeters, and Mueller-matrix ellipsometers,” Appl. Opt. 38(16), 3490–3502 (1999). 46. LabSphere, Inc., “A guide to reflectance coatings and materials,” http://www.labsphere.com/tecdocs.aspx. 47. Luxpop, Inc. http://www.luxpop.com/.


Introduction
Rough surface scattering has been an active area of research for nearly half a century. One of the early areas of research dealt with radio [1] and acoustic [2,3] wave scatter from the ocean surface. Measurement of this scatter led to methods for sensing ocean wave heights [4]. In another application, in preparation for the NASA Apollo missions, analysis of light scattered from the lunar surface led researchers to conclude that the moon's surface is composed of a particulate material [5]. Study of the scattering from rough surfaces has also been applied to predict reflections from tree canopies [6,7], biological/medical sensing [8,9], and computer/video graphics [10][11][12].
Scattering from a rough surface is typically modeled using a bidirectional reflectance distribution function (BRDF) or its polarimetric counterpart, a polarimetric BRDF. BRDFs are generally classified in two main types-empirical and analytical BRDFs. As the name implies, empirical BRDFs are formulated from measurements. An example of such a BRDF can be found in Ref. [13]. Analytical BRDFs are typically derived using either physical optics or geometrical optics. Physical optics BRDFs rely on the electromagnetic physical optics approximation [14,15] (known as the Kirchoff approximation in optics) to predict the scatter from rough surfaces. The seminal work in this type of BRDF is that of Beckmann and Spizzichino [16]. Another work which presents an excellent description of the Kirchoff approximation of rough surface scattering is that of Ishimaru [17]. More recently, Beckmann and Spizzichino's BRDF has been extended to include shadowing (described below) and polarization [12,18]. Geometrical optics BRDFs rely on the ray approximation of light [15,19]. The seminal paper in this BRDF genre is that of Torrance and Sparrow [20]. Since Torrance and Sparrow, numerous geometrical optics BRDFs have been developed. These include BRDFs specialized to predict IR signatures of aircraft [21], full polarimetric geometrical optics BRDFs [22][23][24][25][26] , and BRDFs derived to predict scatter from multilayer coatings [27]. Both types of analytical BRDFs discussed here require the surface roughness features to be several times larger than the wavelength of the incident light. Since physical optics BRDFs are based on a more sound approximation, they tend to be more accurate than geometrical optics BRDFs. However, geometrical optics BRDFs tend to be simpler in mathematical form and thus more numerically efficient. They also tend to be more physically intuitive than their physical optics counterparts. It is for these reasons that the BRDF introduced in this paper is a geometrical optics BRDF.
The geometrical optics BRDF introduced here is a full polarimetric BRDF (pBRDF). It is very similar in form to the pBRDF introduced by Priest and Germer [23,24]; however, this pBRDF includes a shadowing function (in particular, the Torrance and Sparrow [20] shadowing function) and a Lambertian (diffuse) pBRDF component. It is shown that the pBRDF satisfies reciprocity and conserves energy. Section 2 of this paper introduces the theoretical form of the pBRDF. Section 3 compares the pBRDF prediction of the scatter from a rough, perfectly reflecting surface to that of an exact, electromagnetic solution. In Section 4, Mueller matrix predictions are made using the pBRDF and compared to measurements made using a Mueller matrix ellipsometer [28] in order to validate the model. Lastly, the paper is concluded with a Macroscopic surface scattering geometry. Light subtending solid angle dω i is incident from the (θ i , φ i ) direction on a small area dA of a much larger rough surface with complex index of refraction η = n − jκ. Light is scattered and observed within solid angle dω r in the (θ r , φ r ) direction. summary of the work presented.

Methodology
Consider the scattering geometry shown in Fig. 1. Light, subtending solid angle dω i and centered on polar angle θ i and azimuth angle φ i , is incident on a small area dA of a larger rough surface with complex index of refraction η = n − jκ. Light is scattered from the surface and observed within solid angle dω r at polar angle θ r and azimuth angle φ r . The BRDF is defined as the ratio of the scattered radiance to the incident irradiance: where dL r is the scattered radiance, dE i is the incident irradiance, and φ = |φ r − φ i | (i.e., surface is isotropic and homogeneous). Note that the incident irradiance is equivalent to the product of the incident radiance and the projected solid angle cos θ i dω i [20,[29][30][31][32][33][34][35]. In the general polarimetric case, the scattered radiance and incident irradiance are replaced by Stokes vectors and the BRDF (now a pBRDF) by a Mueller matrix, i.e. [23,24,26,29,30,33,34,36,37], Note that f of Eq. (1) equals F 00 where the subscript 00 is the element in the first row and first column of the pBRDF Mueller matrix [33,34]. It is common in literature to express the scattered radiance as the sum of the radiance which leaves the surface after one reflection and the radiance which leaves the surface after multiple reflections [11, 20, 25-27, 33-35, 37, 38]: implying that Fig. 2. Scattering geometry of a single microfacet. The angle α is the polar angle from the mean surface normal to the microfacet normal n. The angle β is the incident angle onto and reflected angle from a microfacet as measured from the microfacet normal. The angle γ i is the angle between the macroscopic plane of incidence and the scattering plane of the microfacet (depicted in the figure as the plane containing the vectors n and t). Likewise, the angle γ r is the angle between the macroscopic plane of reflection and the scattering plane of the microfacet.
It follows that light leaving the surface after a single scattering event (for convenience termed the first reflection) models the specular component of reflection. This component carries with it all polarimetric information gained from interaction with the material surface. The multiple scattering term models the diffuse, or Lambertian component of reflection. One can glean this by considering the effect on polarization after multiple, random surface reflections. In general, the first reflection is partially polarized and directed in the specular direction relative to the local surface normal of the illuminated area. If a portion of the first reflection is incident on another part of the material surface, the second reflection is a partially polarized version of the first, directed, once again, in the specular direction relative to the local surface normal of that illuminated area. This being repeated numerous times results in scatter which is unpolarized and of uniform radiance throughout the scattering hemisphere. Therefore, a pBRDF can be expressed as the sum of a polarized, specular component and of an unpolarized, diffuse component [11, 20, 25, 26, 33-35, 37, 38]: As noted by Sun [35], Maxwell and Beard [38], and Ellis [27], this assumption is not always valid; however, it is an assumption which is very common in literature and, for the purpose of keeping the mathematical form of the pBRDF simple, is utilized here. In the next subsection, the form of the specular pBRDF component is shown and discussed.

Specular pBRDF component
The pBRDF in this paper makes use of the microfacet surface model introduced by Torrance and Sparrow [20]. The model assumes that a rough surface is composed of a collection of randomly (according to some distribution) oriented facets each scattering light in the manner stipulated by Fresnel's equations (see Fig. 2). For Fresnel's equations to be an accurate model for reflection from the surface of a microfacet, the size of the facet must be large compared to the incident wavelength λ . This implies that the "roughness" of the macroscopic surface should be large compared to λ . As discussed by Sun [35], the microfacet surface model can be where P is the distribution function modeling the orientation of the facets around the mean surface normal (z direction in Figs. 1 and 2), M is the Mueller matrix modeling the polarimetric scattering from the material surface, and G is the visibility function (shadowing/masking factor) [23][24][25][26][33][34][35]. Note that the angles α and β are derived using spherical trigonometry: cos α = (cos θ i + cos θ r ) / (2 cos β ) cos 2β = cos θ i cos θ r + sin θ i sin θ r cos φ .
For the sake of brevity, the derivation of Eq. (6) is not shown. A full derivation of the expression can be found in Ref. [35]. It should be noted that Eq. (6) differs from that given in Refs. [33,34] by a factor of cos α in the denominator. Shell [34] notes that the impact of this factor is minimal since the BRDF magnitude rapidly decreases with increasing α such that division by a decreasing cos α is negligible.
Since the height distributions of most natural surfaces are Gaussian [35], the facet distribution function utilized in this pBRDF is that of Beckmann [16]: Note that while this facet distribution function is very common, utilized in Refs. [11,23,24,26,35], other facet distributions do exist. For example, Torrance and Sparrow [20] utilize a simpler, unnormalized distribution function of Gaussian shape, Shell [34] and Gartley [39] discuss where E s i and E s r are the s-pol (perpendicular polarization), incident and reflected, complex electric field components, E p i and E p r are the p-pol (parallel polarization), incident and reflected, complex electric field components, and r s and r p are the complex Fresnel field reflection coefficients for the s-and p-pol, respectively [23][24][25][26]33,34,39]. Note that s-and p-polarization for E s i , E s r , E p i , and E p r are defined with respect to the macroscopic coordinate system (see Fig. 1); however, the complex Fresnel field reflection coefficients are defined with respect to the microfacet coordinate system (see Fig. 2). Thus, it is necessary to perform coordinate system rotations to align the macroscopic planes of incidence and reflection with the microfacet planes of incidence and reflection. This fact explains the rotation matrices in Eq. (9). Relating the angles γ i and γ r to the macroscopic angles is, once again, accomplished using trigonometry [23-26, 33, 34, 39]: Before converting the Jones matrix into a Mueller matrix, it is important to discuss briefly the physical interpretation of the Jones scattering matrix in Eq. (9). The Jones matrix elements, T ss , T ps , T sp , and T pp , can be interpreted as modeling how s-or p-pol incident light couples into s-or p-pol reflected light. For instance, the T ps element models how incident s-pol light couples into p-pol reflected light. The other elements can be interpreted in a similar manner. For scattering in the specular plane (i.e., φ = π), intuition leads one to conclude that the Jones scattering matrix in Eq. (9) becomes diagonal, i.e., incident s-pol and incident p-pol couple into reflected s-pol and reflected p-pol, respectively. This conclusion is easily confirmed by substituting φ = π into Eq. (10). Converting the Jones matrix in Eq. (9) to a Mueller matrix is performed using the analysis in Ref. [30]: where * is the complex conjugate operation. For the sake of brevity, the expressions for all 16 elements are not shown. The Mueller matrix represented in Eq. (11) occurs when observation is confined to the specular plane as is done for the measurement results presented in Section 4. The full expressions can be found in Refs. [23-26, 33, 34, 39]. Note that the Mueller matrix in Eq. (11) has a similar physical interpretation as the Jones scattering matrix discussed above. For instance, the M 23 element of the above Mueller matrix models how the fourth Stokes parameter (circular polarization) couples into the third Stokes parameter (linear polarization) upon reflection.
This expression is derived assuming that each microfacet comprises one side of a symmetric vshaped groove (see Fig. 3). Mathematically, G determines the fraction of an illuminated microfacet which contributes to the scattered radiance. Physically, as is evident from Fig. 3, G models the incident and reflected light blocked by adjacent microfacets. Most of the BRDFs/pBRDFs in literature include a shadowing/masking function of some form [11,20,25,26,35,38]. The function is instrumental in keeping the BRDF bounded and thus satisfying the conservation of energy. A notable exception to this is the pBRDF of Priest and Germer [23,24]. The lack of G causes their pBRDF to asymptotically approach infinity as the angle of incidence or observation approaches grazing [25,26]. The desired specular component of the pBRDF can now be formed by substituting the facet distribution function P [Eq. (8) Fig. 3. Scattering geometry of a v-shaped groove. The top subfigure depicts shadowing while the bottom subfigure depicts masking. Shadowing occurs when the angle of incidence approaches grazing. Similarly, masking occurs when the angle of observation nears grazing. stated above, G plays a critical role in keeping the pBRDF bounded and thus a realistic physical model. In order to demonstrate the function's importance, consider the pBRDF predictions shown in Fig. 4. The figure shows traces comparing the F 00 elements of the Priest and Germer pBRDF [23,24] and the pBRDF in Eq. (13) for θ i = 45 • , 60 • , 75 • , and 85 • with 2 1/2 σ h / = 0.3. The pBRDFs are evaluated in the specular plane (φ = π) and using a perfect reflecting surface, i.e., a perfect electric conductor (PEC). The figure clearly shows that the pBRDF in Eq. (13) remains bounded while the Priest and Germer pBRDF diverges as θ r approaches 90 • . Having developed and discussed the specular pBRDF component, attention can now be turned to the diffuse component.

Diffuse pBRDF component
Before the diffuse component of the pBRDF can be developed, the concept of directional hemispherical reflectance (DHR) must be introduced. The DHR is defined as the ratio of the total energy reflected into the entire hemisphere above a material surface to the total energy incident from a particular direction [23-26, 30, 33-35, 39]: Note that ρ DHR 1, otherwise the pBRDF violates the conservation of energy (assuming a passive material). The stated condition becomes an equality when the surface is a PEC. Substituting Eq. (5) into the DHR expression and applying the equality condition (a PEC surface) .
Note that Eq. (16) represents the fraction of scattered light not comprising the specular lobe, or equivalently, the fraction of light which is scattered multiple times. Therefore, generalizing Eq. (16) to a surface other than a PEC is simply a matter of multiplying F d, PEC 00 times M 00 [25,26]. Eq. (16) possesses two notable characteristics which make it an attractive model for the diffuse pBRDF component. First, it depends only on the angle of incidence and the statistical properties of the rough surface. One's intuition dictates that for a "smooth" surface light leaves the surface after a single reflection and therefore most of the scattered radiance is contained within the specular lobe. This can be verified mathematically by noting that as σ h → 0, the facet distribution function in Eq. (13) becomes a Dirac delta function. Substituting this expression into Eq. (16) results in F d 00 = 0. Likewise, as surface roughness increases, one would expect light to be scattered multiple times before leaving the material surface. Mathematically this can be verified by observing that F s 00 → 0 as σ h → ∞. Substituting F s 00 = 0 into Eq. (16) produces the trivial result F d 00 = 1 (i.e., pure diffuse scattering). The second characteristic of note in favor of modeling the diffuse pBRDF component in the manner outlined above is that no fitted coefficients are required to model the strength of the diffuse pBRDF component. The use of coefficients, whose values are determined by fitting the BRDF to measured data, is a common feature in other BRDFs [11,20,33,34,[37][38][39].

Summary of theory
Combining the specular pBRDF term [Eq. (13)] and the diffuse pBRDF term [Eq. (16)] produces the desired result [25,26] Note that since it is assumed to be unpolarized, the diffuse component only contributes to the F 00 element of the pBRDF Mueller matrix.
In order to show that the above expression satisfies electromagnetic reciprocity, θ i and φ i must be switched with θ r and φ r . It is easy to show that doing so produces the same expression as that in Eq. (17); thus, the pBRDF satisfies the reciprocity condition. Proving that Eq. (17) conserves energy requires one to show that ρ DHR 1. Note that the conservation of energy is enforced when finding the value of the diffuse pBRDF component (detailed above). Therefore, Eq. (17) conserves energy as well.
Summarizing the theory, Eq. (17) possesses two characteristics which distinguishes it from existing geometrical optics pBRDFs in literature. The first is the addition of the shadowing/masking function G. As discussed above, G keeps Eq. (17) bounded and thus a realistic physical model. The second is the development of a diffuse pBRDF component. As previously stated, this component depends only on physical parameters and does not need to be fit to measured data. In the next section, predictions made using Eq. (17) of a rough, PEC surface are compared to Method of Moments [40] (MoM) solutions for the purpose of validating the model.

Simulation
Before analyzing the simulation results, a brief background on the MoM is warranted. The MoM is a technique to solve integral equations which arise frequently in electromagnetics. The problem of interest here is a 15, 000λ long, random (surface height is Gaussian distributed) PEC surface. The surface is invariant in the z direction (see Fig. 5) significantly reducing the number of unknowns in the problem. Also, as is shown in Fig. 5, only s-pol is considered here. The electric field integral equation (EFIE) for the scattering problem depicted in Fig. 5 (assuming plane wave excitation) is formulated by applying the transverse electric field boundary condition at the random, PEC surface, i.e., E z where Z 0 is the intrinsic impedance of free-space (approximately 377 Ω), ρ ρ ρ = xx + yy is the observation vector, ρ ρ ρ = xx + yy is the source vector, k i = x sin θ i − y cos θ i is the propagation vector of the incident field, J z is the current induced on the PEC surface by the field, and H (2) 0 is a zeroth order Hankel function of the second kind. Note that the integral in Eq. (18) is over the parameterized surface contour denoted by C . The unknown in Eq. (18) is the surface current J z . Note that assuming J = 2n × H i forms the basis of the physical optics, or Kirchoff approximation [14][15][16][17]. In the MoM, J z is expanded in a set of basis functions (in this case, fixed width pulses): After simplification, the resulting expression is then tested using another set of functions (in this case, Dirac delta functions located at the center of each pulse) to produce an N × N matrix equation where N is the number of unknowns, i.e., ⎡ Note, for example, that the 2N element of the MoM matrix shown above models how the N th source current segment contributes to the scattered field at the 2 nd observation segment. The other elements of the MoM matrix can be interpreted in a similar manner. Solving Eq. (20) yields the unknown current. Once J z is computed, the scattered field can be found at any observation point by α n C n exp j 2π λ x sin θ r + y cos θ r dC n (21) where C n is the segment of the parameterized contour represented by the n th pulse and ρ = x 2 + y 2 1/2 is the Euclidean distance from the origin. The second line of the above expression assumes that the observation point ρ is in the far-field as defined by Fraunhofer [15]. Note that the MoM solution shown above is a coherent field solution. Since the incoherent solution is the one desired, the 15, 000λ surface is divided up into M = 100 partitions and the scattered field from each partition is summed incoherently [41][42][43][44]. Also, in order to minimize the effect of edge diffraction from the surface partitions, i.e., approximate an infinite surface, J z is windowed using a Gaussian taper: where the index m represents the m th surface partition, x m is the center of the m th partition, and w is the taper width [41][42][43][44]. The incoherent, far-field reflectance distribution can now be found from where E z r,m is the scattered field from the m th partition [41][42][43][44]. Note that the above expression is the average incoherent radar cross section (RCS) of the random, PEC surface normalized by the effective illumination length. Detailed analysis of these steps can be found in Refs. [41][42][43][44].
Having provided the necessary background on the MoM, attention can now be turned to the simulation results.
The simulation results are shown in Fig. 6. As mentioned above, the simulation surface is a 15, 000λ long, random (surface height is Gaussian distributed) PEC surface. The Gaussian surface is generated as shown in Ref. [43] with roughness equal to 2 1/2 σ h / = 0.3. The traces on the figure are far-field reflectance distributions for θ i = 10 • , 30 • , 45 • , 60 • , and 75 • . Note that the reflectance distributions in the figure are normalized with respect to their values at the specular angles (θ i = θ r ), and observation for both the MoM and pBRDF predictions is in the specular plane (φ = π). Overall, the pBRDF predictions match very well with the exact, MoM solutions. At some incident angles, the pBRDF predictions deviate from the MoM solutions; however, there is almost unanimous agreement between the pBRDF and MoM solutions on the locations and magnitudes of reflectance maxima. Note that ripples are visible in the MoM solution traces for θ i = 45 • , 60 • , and 75 • . These ripples are caused by constructive and destructive interference in the scattered field. Incoherent scatter should not interfere; however, as discussed above, the MoM solution is coherent. The incoherent scatter is being approximated by summing the scattered field incoherently over 100 partitions of the entire 15, 000λ surface. Although this should be sufficient [41], some interference is still occurring. Summing the scattered field incoherently over more partitions should lessen the interference ripples; however, the cost is a longer simulation run time. The results shown here are sufficient to demonstrate the validity of the pBRDF. In the next section, Mueller matrix element predictions made using Eq. (17) are compared to experimental measurement results in order to further validate the pBRDF.

Mueller matrix measurement results
The instrument used to collect the Mueller matrix data presented here is an ellipsometer [28] at the Air Force Research Laboratory (AFRL), Wright-Patterson Air Force Base, Ohio. A photograph of the ellipsometer is shown in Fig. 7. It, like all ellipsometers, consists of two arms-the polarization state generator (PSG) and the polarization state analyzer (PSA). The material under test (MUT) is placed in between the PSG and PSA in a sample holder which is able to rotate. The PSG of the instrument shown in Fig. 7 consists of a 1064 nm laser followed by polarization optics mounted on a stationary optical rail. Note that the polarization state of the laser is set using a polarizing beam splitter. The PSG polarization optics consist of two half-wave plates and a quarter-wave plate. They are contained within a mechanical housing which allows them to be moved precisely into and out of the source beam. This allows the PSG to generate four independent polarization states to interrogate the MUT. As is common, the PSA of the ellipsometer shown in Fig. 7 is a mirror image of the PSG. It consists of a set of polarization optics followed by a horizontal linear polarizer and a detector mounted on a rotating base. The PSA is able to rotate independently of the MUT sample holder, thus allowing any (θ i , θ r ) to be measured. The polarization optics used in the PSA are identical to and are contained within the same type of mechanical housing as that of the PSG. This allows the PSA to analyze four independent polarization states. Overall, the instrument is able to make 16 independent, polarimetric measurements of a MUT, thereby providing all the necessary information to deduce the MUT's Mueller matrix.
Before the MUT is measured, the instrument is calibrated using the Eigenvalue Calibration Method (ECM) [45]. The ECM is a calibration technique developed by Compain et al. [45] in which the Mueller matrices of a set of known standards are measured in order to compute the experimental Mueller matrices for the PSG and PSA. Once these matrices have been determined (see Ref. [45] for details), the desired Mueller matrix of the MUT can be found by where S is a 4×4 matrix of measured (ellipsometric) irradiances, W is the experimental Mueller matrix of the PSG, and A is the experimental Mueller matrix of the PSA. In the measurement results presented here, the standards used to calibrate the ellipsometer are a no sample measurement, a linear polarizer, and a quarter-wave plate. The linear polarizer is measured at 0 • , 60 • , and 120 • ; the quarter-wave plate is measured at 45 • . The Mueller matrix measurement results are shown in Figs. 8 and 9. Fig. 8 shows the Mueller matrix measurement results for LabSphere Infragold [46]. The complex index of refraction for gold (η = 0.285 − j7.3523) is obtained from Ref. [47]. The figure shows the experimental M 01 , M 11 , M 22 , and M 23 elements compared to predictions made using the pBRDF in Eq. (17). Note that all measurements are made in the specular plane (φ = π), and the Mueller matrix elements are normalized with respect to the M 00 element as annotated on the figure. The plotted values for the measured Mueller matrix elements of LabSphere Infragold are the means of 256 irradiance measurements. The bars on the figure represent ±1σ , i.e., one standard deviation of those 256 measurements. Note that the LabSphere Infragold results are consistent with those published by Priest and Meier [24]. Although there are discrepancies between the measured Mueller matrix elements and the pBRDF (especially at 10 • and 20 • ), the pBRDF predictions agree well with the measurements. The most important aspect of the results is that the pBRDF captures the trend of the data (i.e., the physics of the material surface interaction) as the observation angle (or equivalently incident angle) approaches grazing. Figure 9 shows the Mueller matrix measurement results for flame sprayed aluminum (FSA). The complex index of refraction for aluminum is η = 1.226 − j10.413 [47]. Shown in the figure are the same Mueller matrix elements as Fig. 8  specular plane (φ = π), and the Mueller matrix elements are normalized with respect to the M 00 element. Once again, the plotted values for the measured Mueller matrix elements of FSA are the means of 256 irradiance measurements. The bars on the figure represent ±1σ of those 256 measurements. As is the case in Fig. 8, discrepancies do exist between the measured data and the pBRDF predictions; however, the pBRDF, once again, captures the trend of the measured data.
Before concluding, it is worth discussing a possible cause of the discrepancies observed in Figs. 8 and 9. It is assumed that the published values of η for gold and aluminum are accurate for the LabSphere Infragold and FSA samples measured in this experiment. As noted in Ref. [47], the value of η can vary greatly depending on sample quality, sample preparation, or measurement technique. If, instead of using the published values for η, the best values for η are found via nonlinear least squares, one obtains for LabSphere Infragold η = 0.4364 − j5.2526 and for FSA η = 0.8886 − j3.4602. While it is possible that the best-fit index for LabSphere Infragold could be more representative of the true η value for the specimen (considering how LabSphere Infragold is manufactured [46]), the best-fit FSA index is more difficult to explain. One possibility is that the FSA specimen used in this experiment is slightly oxidized (i.e., a thin coating of Al 2 O 3 ). This hypothesis would explain the sharper than predicted rise in the measured M 22 /M 00 element in Fig. 9.

Conclusion
In this paper, a geometrical optics pBRDF is presented. As discussed, the pBRDF is composed of a specular (single reflection) component and a diffuse (multiple reflection) component. The specular component, derived using the microfacet surface model [20], is shown to consist of a facet distribution function, a Mueller matrix modeling the polarimetric scattering from the material surface, and a visibility (shadowing/masking) function. Each one of these constituent functions is discussed in detail. The diffuse component is derived using the DHR and the conservation of energy. It is shown that a diffuse pBRDF component derived in this fashion de- pends only on geometrical parameters (angle of incidence, surface height standard deviation, and surface correlation length) and does not require coefficients fit to measured data. Taken as a whole (i.e., adding the specular and diffuse components), the pBRDF presented in this paper is shown to satisfy electromagnetic reciprocity and the conservation of energy. Lastly, in order to validate the pBRDF, predictions made using the pBRDF are compared to MoM solutions of a rough, PEC surface and experimental Mueller matrix data for two rough, metallic samples. It is shown, via these results, that the pBRDF accurately models the physics of the light/material surface interaction. The pBRDF presented in this paper possesses two characteristics which distinguishes it from existing geometrical optics pBRDFs in literature. The first is the addition of the visibility (shadowing/masking) function. As shown and discussed, the visibility function keeps the pBRDF bounded and thus a realistic physical model. The second is the development of a diffuse pBRDF component. This component allows for better modeling of rough, reflective surfaces which tend to depolarize light via multiple surface reflections.