Principal components analysis on the spectral bidirectional reflectance distribution function of ceramic colour standards

The Bidirectional Reflectance Distribution Function (BRDF) is essential to characterize an object’s reflectance properties. This function depends both on the various illumination-observation geometries as well as on the wavelength. As a result, the comprehensive interpretation of the data becomes rather complex. In this work we assess the use of the multivariable analysis technique of Principal Components Analysis (PCA) applied to the experimental BRDF data of a ceramic colour standard. It will be shown that the result may be linked to the various reflection processes occurring on the surface, assuming that the incoming spectral distribution is affected by each one of these processes in a specific manner. Moreover, this procedure facilitates the task of interpolating a series of BRDF measurements obtained for a particular sample. © 2011 Optical Society of America OCIS codes:(290.1483) BSDF, BRDF, and BTDF; (030.5630) Radiometry; (120.5820) Scattering measurements. References and links 1. F. E. Nicodemus, J. C. Richmond, and J. J. Hsia, “Geometrical considerations and nomenclature for reflectance,” National Bureau of Standards Monograph (National Bureau of Standards, 1977), Vol. 160. 2. C. Bordier, C. Andraud, and J. Lafait, “Model of light scattering that includes polarization effects by multilayered media,” J. Opt. Soc. Am. A25, 1406–1419 (2008). 3. L. Simonot, “Photometric model of diffuse surfaces described as a distribution of interfaced Lambertian facets,” Appl. Opt.48,5793–5801 (2009). 4. R. L. Cook and K. E. Torrance, “A reflectance model for computer graphics,” Technical report, Computer Graphics (ACM, 1981), Vol. 15, No. 3. 5. B. T. Phong, “Illumination for computer generated pictures,” Commun. ACM 18(6), 311–317 (1975). 6. G. J. Ward, “Measuring and modelling anisotropic reflection,” Comput. Graphics 26,265–272 (1992). 7. X. D. He, K. E. Torrance, F. X. Sillion, and D. P. Greenberg, “A comprehensive physical model for light reflection,” Technical report, Computer Graphics (1991), Vol. 25, No. 4. 8. J. F. Blinn, “Models of light reflection for computer synthesized pictures,” Comput. Graphics 11,192–198 (1977). 9. E. P. F. Lafortune, S. C. Foo, K. E. Torrance, and D. P. Greenberg, “Non-linear approximation of reflectance functions,” Technical report (Cornell University, 1997). 10. S. H. Westin, H. Li, and K. E. Torrance, “A comparison of four BRDF models,” Technical report PCG-04-02, Program of Computer Graphics (Cornell University, April 2004). 11. A. Ngan, F. Durand, and W. Matusik, “Experimental analysis of BRDF models,” Eurographics Symposium on Rendering, K. Bala and P. Dutre, eds. (2005). #145188 $15.00 USD Received 27 Jun 2011; revised 19 Aug 2011; accepted 23 Aug 2011; published 19 Sep 2011 (C) 2011 OSA 26 September 2011 / Vol. 19, No. 20 / OPTICS EXPRESS 19199 12. I. G. E. Renhorn, and G. D. Boreman, “Analytical fitting model for rough-surface BRDF,” Opt. Express 16(17), 12892–12898 (2008). 13. J. L. Simonds, “Application of characteristic vector analysis to photographic and optical response data,” J. Opt. Soc. Am.53,968–971 (1963). 14. J. M. Ĺopez-Alonso, J. Alda, and E. Bernab éu, “Principal-component characterization of noise for infrared images,” Appl. Opt.41,320–331 (2002). 15. A. Ferrero, J. Alda, J. Campos, J. M. L ópez-Alonso, and A. Pons, “Principal components analysis of the photoresponse nonuniformity of a matrix detector,” Appl. Opt. 46,9–17 (2007). 16. A. M. Rabal, A. Ferrero, J. L. Fontecha, A. Pons, J. Campos, A. Corr óns, and A. M. Rubio, “Goniospectrophotometer for low-uncertainty measurements of bidirectional scattering distribution function (BSDF),” Proceedings of CIE Expert Symposium on“Spectral and Imaging Methods for Photometry and Radiometry,” Publication CIE x036:2010 (CIE, Vienna, Austria, 2010), pp. 79–84. 17. T. A. Germer and C. C. Asmail, “Goniometric optical scatter instrument for out-of-plane ellipsometry measurements,” Rev. Sci. Instrum. 70,3688–3695 (1999). 18. D. Hünerhoff, U. Grusemann, and A H öpe,“New robot-based gonioreflectometer for measuring spectral diffuse reflection.” Metrologia43,S11–S16 (2006). 19. F. B. Leloup, S. Forment, P. Dutr é, M. R. Pointer, and P. Hanselaer, “Design of an instrument for measuring the spectral bidirectional scatter distribution function,” Appl. Opt. 47(29), 5454–5467 (2008). 20. R. Corey, M. Kissner, and P. Saulnier, “Coherent backscattering of light,” Am. J. Phys. 63,561–564 (1995). 21. T. J. Papetti, W. E. Walker, C. E. Keffer, and B. E. Johnson, “Coherent backscatter: measurement of the retroreflective BRDF peak exhibited by several surfaces relevant to ladar applications,” Proc. SPIE 6682, 66820E (2007).


Introduction
An object's reflectance properties are fully characterized in each of its surface points by the Bidirectional Reflectance Distribution Function (BRDF).This function describes the way in which a surface reflects the incoming optical radiation; that is, the amount of radiance per unit of irradiance that is reflected in each direction as a function of the illumination direction [1].This way, the determination of a sample's BRDF requires a thorough sampling with different angular configurations.Moreover, the sampling becomes more complex if the object's reflecting properties are such that symmetry and isotropy simplifications are not possible.In this sense, it would be extremely useful to have physical optics-based analytical models [2][3][4][5][6][7][8][9][10][11][12].Unfortunately, taking into account the great diversity of existing surfaces, the quest for a unique model is an unfeasible task.Nonetheless, it is possible to provide a set of guidelines and methods that can be applied to each and every object.In this context, when searching for geometric-based or diffraction-based models, it is a common approach to assume that the BRDF is made up of three independent components: • Specular contribution ( f r,sp ): accounts for mirror-like reflections from the mean plane of the reflecting surface.It can be computed based on Fresnel reflectivity.It is observed when the surface's roughness is small compared with the wavelength.All the photons coming from a given direction, when they interact with the surface, are reflected in the same specular direction.
• Directional diffuse contribution ( f r,dd ): accounts for the field scattered to the hemisphere but with a certain directional character.It occurs when the photons coming from a given direction, when they interact with the surface, change their direction in a different but correlated way.
• Uniform diffuse contribution ( f r,ud ): accounts for the field diffused uniformly over all the hemisphere.In this case the photons coming from a given direction, when they interact with the surface, change their direction in a different and uncorrelated way.
Consequently, each component originates from a different type of photon-surface interaction.A three-component BRDF is, unquestionably, the result of a phenomenological simplification, since in reality there should be an undetermined number of potential interactions.
The Principal Components Analysis (PCA), which is a multivariable analysis technique, has proven to be effective for the analysis or characterization of multivariate response data [13][14][15].This procedure makes it possible to reduce the multidimensionality of the response data to an intrinsically minimum number.The goal of this work is to show how the PCA technique applied to BRDF data may help to the identification of the different interaction processes occurring on the surface, assuming that each one of these processes has a specific impact on the incoming spectral distribution (the probability of a photon undergoing an interaction of a particular type must be a function of its energy; i.e. its wavelength).The results emerging from this type of procedure allow us to gain a deeper understanding of the individual contributions to the BRDF, and simplify to a great extent the search for models.
The present article starts with a brief review of PCA; afterwards it is shown how, based on experimental data, PCA makes it possible to perform a more detailed analysis of BRDF that could be used as starting point to develop specific models for different surface types.

Principal Components Analysis
The Principal Components Analysis is a powerful statistical technique capable of identifying and quantifying orthogonal contributions to the total variance of a collection of data.The key element is the definition of a multidimensional random variable F = {F 1 , F 2 , . . ., F N }, having M realizations of it.For the BRDF data, this variable is made up of a concatenation of N F i spectra, each one corresponding to the measured BRDF for a given angular configuration.The set of spectral values extracted from the different spectra corresponding to a particular wavelength is taken as a particular realization of the multidimensional variable; the number of realizations, M, coincides with the number of wavelengths.In the case tackled in this paper, the collection of data is arranged as a matrix having N columns and M rows.Each row contains the BRDF for every angular configuration at a specific wavelength.The covariance matrix of the spectra is obtained as an N × N matrix, S F .This covariance matrix is typically non diagonal, thus revealing the inner correlations between the spectra.The PCA method diagonalizes this covariance matrix and produces three types of elements: N eigenvalues, γ j , N eigenvectors, e i, j , and N eigenspectra, A j .The eigenvalues are arranged in decreasing order.They quantify the importance and the contribution to the total variance of the data of their associated eigenspectra.Actually, γ j is the variance of the eigenspectra A j .The eigenvectors can be seen as the coef- ficients of the transformation from the correlated variables given by the spectra to a new set of uncorrelated variables expressed by the eigenspectra.These eigenvectors can be arranged as an N × N orthogonal matrix, E T = E −1 ( T means transposition).Finally, the eigenspectra, A j , describe uncorrelated data that provide spectral insight about the different contributions to the BRDF.A j is calculated as: where Fi is obtained by subtracting the mean from the original i-th spectrum, The importance of a given eigenspectrum, A j , within the data set is quantified by its associated eigenvalue, γ j .By using the eigenvectors it is possible to migrate from the experimental coordinate system, {F 1 , . . ., F N }, to the eigenspectra coordinate system, {A 1 , . . ., A N }, and vice versa.This transformation can be written as follows: From this equation we may obtain a filtered version of the spectra by selecting a customized subset of eigenvectors and their associated eigenspectra.This is easily done by choosing a subset of subscripts within the sum.Then, it is also possible to remove or to select a given collection of contributions characterized by their associated eigenspectra.An interesting calculation can be done to express the variance of the original spectra as a combination of the variances associated with the eigenspectra.Thus, the variance of the spectrum F i is given by: where Fλ i is the value on the wavelength identified as λ , in the zero-mean spectrum Fi .Equation (2) incorporates the principal component decomposition in the calculation of Fλ i 2 as follows: After including the previous equation into the variance calculation, it is possible to go from the sum along the M wavelengths, labelled with λ super-scripts, to the double sum in j and k.This yields the following expression: where we applied the fact that the eigenvectors are uncorrelated, which allows us to cancel all the terms having j = k.When j = k the resulting term is the variance of the eigenvectors A j , which, by definition, is equal to γ j .This previous equation can be written as the following matrix product: where S2 is a column vector containing the values of the variance of the original spectra, ; G is a column vector containing the eigenvalues, G j = γ j ; and E2 is a matrix whose elements are the squares of the components of the eigenvectors, E2 i, j = e 2 i, j .Following the previous derivation we may conclude that the contribution of a given spectrum, A j , to the variance of a selected original spectrum, F i , is given by e 2 i, j γ j .Similarly, the contribution to the standard deviation of the original spectrum i due only to the eigenspectrum j, is given by: Due to the fact that ∑ N j=1 e 2 i, j = 1, the elements of E2 represent the relative contribution to the variance of the original spectra of the variance of the eigenspectra obtained by the PCA method.

Results
The experimental measurements upon which PCA was applied were obtained by means of a gonio-spectrophotometer GEFE [16].This instrument consists of a fixed, uniform and collimated light source (Xenon lamp) and a robot-arm that places the sample, making it possible to change automatically and simultaneously the two illumination spherical coordinates (θ i , φ i ) and the two observation ones (θ s , φ s ).These spherical coordinates are defined relative to the Sample Coordinate System, whose z axis is parallel to the sample's normal direction (Fig. 1).Moreover, thanks to its periscopic system with beamsplitter, it is possible to measure under retro-reflection conditions.This configuration yields important information about the BRDF, as will be shown further on.For the data acquisition a CS-2000 Konica Minolta spectroradiometer is used.This device operates within the visible range [380 nm -780 nm], performing spectral measurements with a 1 nm spectral sampling interval and a 4 nm bandwidth.The relative uncertainty of the BRDF measurement with k=2 is less than 0.02 (limited by the linearity correction when the measurement ranges within more than 5 decades).This uncertainty is comparable to other instruments with similar features [17][18][19].Figure 2 shows the BRDF spectra depicted in a semilogarithmic plot.The spectra in the specular region are relatively flatter than those in the diffuse region, and their respective average values differ in more than two orders of magnitude.Moreover, each set of spectra in the specular region corresponds to a particular angle of incidence θ i .Figures 3 and 4 show that the bigger the angle θ i , the higher the resulting BRDF (with the exception of the 0 o -configuration).Furthermore, these figures clearly reveal the effect of θ i upon the spectral distribution: Whereas the spectral distributions at 40 o , 50 o , 60 o and 70 o are very similar, those recorded at 10 o and 20 o are slightly different from the rest, and also different from each other.On the other hand, within the diffuse region all spectra are rather similar (Fig. 2).The high frequency noise in the spectra (e.g. between 450 nm and 500 nm and around 770 nm) is an artefact produced by the division of the measured radiance and the measured irradiance on the sample to obtain the BRDF, and it corresponds to the narrow lines of the Xe lamp used in GEFE.
By dispensing with the spectral information (i.e., by representing only the value for a particular wavelength or the average value), the BRDF's angular distribution can be plotted in greater detail.For this purpose, diagrams such as those shown in Fig. 5 are used.This figure represents the BRDF data for λ = 500 nm measured within the incidence plane.This plane comprises two half-planes, each of which is characterized by a particular angle φ s .We will say that, by defini-  tion, the half-plane containing the incoming light is φ s = 0 o , whereas the half-plane containing the specularly reflected light is φ s = 180 o .
A very interesting effect can be observed in Figs. 3 and 5.In Fig. 3, the specular BRDF spectrum at θ i = 0 o has a higher value that those spectra at θ i = 10 o , 20 o and 30 o .On the other hand, an increase of the BRDF at retro-reflection conditions is observe in Fig. 5.Both effects might be explained by the coherent backscattering of light (CBS) [20,21].
Since Fig. 5 represents the incidence plane, we will call this figure 0:180, according to the notation (first half-plane: second half-plane).Each plot in this figure represents the measured BRDF as a function of the observation polar angle θ s for a specific illumination polar angle θ i .A negative θ s indicates, by definition, that the observation is carried out in the first half-plane.
Figure 6 shows the resulting BRDF curve (λ = 500 nm) when observing from outside the incidence plane (half-planes 30:150), for θ i = 20 o , 30 o , 40 o , 50 o , 60 o and 70 o .This figure reveals more clearly the presence of the BRDF's diffuse component.
As the reader will observe, the diffuse BRDF decreases considerably when either the illumination (θ i ) or the observation (θ s ) angle is large (grazing angles), whereas for smaller angles it remains relatively constant.There is also a certain rise in BRDF with θ s within the 150-halfplane (close to the half-plane that contains the specular direction); in fact, the larger θ i is, the steeper the increase is.Therefore, it can be concluded that the diffuse BRDF is more uniform within the small-angle range (applicable both for illumination (θ i ) as well as for observation (θ s ) angles).The weak point of this type of representation is that it is not possible to assess the spectral dependency with the different angular configurations.We will see that this type of representations can be used to represent and highlight precisely this spectral variation, resorting for this purpose to PCA.

PCA on BRDF data
PCA was applied to BRDF experimental data, as described in section 2 (N = 401 wavelengths and M = 448 BRDF spectra, with a 448 × 448 covariant matrix).The resulting eigenspectra are shown in Fig. 7.These are the spectra associated with the four highest eigenvalues (in descending order) which are, consequently, the most relevant ones when it comes to computing the total variance of the data (their joint contribution accounts for 99.6% of this variance: eigenvalues 0.9864, 0.0066, 0.0018 and 0.0013 for these four eigenspectra and 0.0003 for the first neglected eigenspectrum).We can see that the eigenspectrum PC1 (Fig. 7(a)) looks very similar to the specular spectra in Figs. 3 and 4. Similarly, there seems to be a large similarity between PC2 (Fig. 7(b)) and diffuse spectra in Fig. 2.
It is important a proper choice of the scale of the data, since different scales return different eigenvalues.In our specific case, choosing absolute spectra as input variables to the PCA returns a better result than choosing relative ones (for instance, normalizing the spectra to their variances) .The reason is that specular spectra have a smaller relative variance than diffuse spectra and, in addition, there are much more diffuse spectra than specular spectra in the data.In consequence, applying PCA with a relative scaling would give a smaller weight to the specular components.This consideration should be done before applying PCA to a similar kind of data.
The values on the Y-axis have a relative importance, since it must be recalled that, according to Eq. ( 2), when reconstructing the original variance the eigenspectrum is multiplied by the coordinate of the corresponding eigenvector, which ranges between -1 and 1.Even though the eigenspectra that are represented in the figures are all oriented in a natural way for the sake of clarity (to make interpretation easier), the eigenspectra obtained following PCA are often inverted (multiplied by -1).In order to determine the true orientation of a given eigenspectrum A j , the simplest way is to locate an spectrum i where this eigenspectrum has a substantial contribution to the total variance; then one only needs to look at the sign of e i, j .If the sign is negative then the eigenspectrum A j has to be inverted in the representation, but not in subsequent calculations.How the contribution of an eigenspectrum to the variance of a given spectrum can be computed is explained in the paragraph that follows Eq. (7).Figures 8 and 9 show the relative contribution to the spectral variance of the first four principal components in a semilogarithmic plot.These plots make use of an angular representation, and they correspond to the half-planes 0:180 and 30:150, respectively.No systematic tendency outside the specular region could be extracted from the data of the contributions.It is worth pointing out that, in each configuration, the relative contributions of the four components must add up to nearly one.If this was not the case for a particular scenario, it would imply that a non-negligible component was being left out; in that situation the corresponding eigenspectrum would have to be added.The impact of leaving out these excluded eigenspectra will be rigorously analyzed and assessed in Section 4.

Physical interpretation of the eigenspectra of the BRDF
In the light of Figs. 8 and 9, one can made the following analysis about the physical meaning of each eigenspectrum: • Eigenspectrum PC1.This component always makes a contribution in specular-reflection scenarios (Fig. 8).It never contributes significantly in other configuration scenarios (Fig. 9).It is the only component that always contributes significantly in specular conditions; that is why it can be stated that eigenspectrum PC1 is influenced by Fresnel reflection.spectral distributions associated to eigenspectra PC2 and PC3 originate from different physical phenomena.
• Eigenspectrum PC4.Even though this component's contribution clearly exceeds that of eigenspectrum PC1 under non-specular observation conditions, its contribution is nonetheless much smaller than that of eigenspectra PC2 and PC3.It is under specular conditions with θ i = 10 o and θ i = 20 o where this component's contribution becomes really significant (see also Fig. 3).Perhaps this could be explained by the presence in the material of a second flat layer, where small incidence angles could be reflected and transmitted, but where large incidence angles would be trapped due to total reflection.
The idea behind this analysis is that it is possible to use the resulting information for identifying the uncorrelated components of the BRDF to interactions in different layers of the surface.A model of layers can be proposed from the PCA information as starting point for interpolation with physical models.It is necessary to notice that uncorrelated variables are not always independent.To make sure that the number of layers is not overestimated, it is possible to apply an algorithm of Independent Component Analysis to the PCA result.

PCA's usefulness for the development of BRDF models
The task of finding an analytical model that fits the highly complex BRDF data can be simplified to a large extent if the contributions of the different reflections occurring on the sample can be identified, isolated and analyzed separately.Generally this is not the case, since samples have a complex surface structure that is not known down to the detail.Since PCA yields orthogonal components, it is worth considering whether or not it would be possible to reconstruct the BRDF from a relatively small number of PCA components and a few fitting parameters.This way, an empirical model would be available, allowing us to interpolate the BRDF for any given illumination/observation geometry and to compute the sample's reflectance or reflectance factor with greater accuracy.
In the framework of the principal components, the reconstruction of the BRDF data is carried out by using only those components ( Fi, j ) that interest us: Fi, j = e i, j A j (8) where Fi, j is the resulting reconstruction of Fi when only A j is used.If we recall that, as stated in a previous section, Fi = F i − F i , then it can be inferred that the F i corresponding to the principal component j can be computed as F i = Fi + κ i, j F i , where κ i, j represents the relative contribution associated with A j to the i spectrum of the BRDF.If we assume that this fraction is proportional to this component's contribution, then under this approximation it can be said that: κ i, j = E2 i, j .This approximation can be used to perform the model's fitting, taking κ i, j as free parameter with an initial value E2 i, j .Strictly speaking, the number of principal components that is required to explain a BRDF spectrum must be the one that is able to reproduce the original data (within a given uncertainty).Figure 10 represents, for the (0:180) half-planes, the relative average error ε r,i resulting from the reconstruction of the angular configurations i that relies, respectively, on the four and five eigenspectra associated with the highest eigenvalues.The average error is computed as follows: where F r,i represents the reconstruction of Fi .The figure reveals, as was to be expected, that the more eigenspectra are used for the reconstruction, the smaller the resulting average error.For the four-component case considered in this work, the relative average error is around 1%, and clearly falls below this value when one additional eigenspectra are used (around 0.5%).

Conclusions
The PCA method has been applied to the BRDF that was experimentally measured on a sample glossy green ceramic colour standard.The analysis leads to the conclusion that only 4 principal components are responsible for most of the variance of the BRDF.The relative weight of each one of these components is linked to the specific geometry of the reflection.Thus, component 1 is more significant (i.e., plays a more significant role) when specular geometry is considered, whereas 2, 3 and 4 are more significant in those geometries where there is diffuse reflection.What makes these components 2, 3 and 4 different is the relative weight that they have in the calculation of the variance of the spectral distribution observed in the specular direction.Consequently, there is a relationship between the principal components and the physical phenomena that govern reflection in this sample.In this context, a model having a physical meaning can be set out based on these principal components.
It has been demonstrated that this principal-component-based model can explain very well the BRDF spectra even when only 4 components are included.With PCA it is possible to represent in a simple way and simultaneously both the spectral and the angular nature of the BRDF.As a result, PCA becomes a powerful tool to undertake the complex task of interpreting this kind of data.Moreover, PCA facilitates the task of interpolating the various BRDF measurements corresponding to a particular sample.The reason lies in its underlying physical model: in PCA, each of the resulting principal components corresponds to a particular component of the BRDF of different nature; also, each principal component can be fitted independently from the rest.

3. 1 .
BRDF data The BRDF measurements shown below are from a glossy Ceramic Colour Standard, CCS (Glossy Green 5 cm × 5 cm).The BRDF was sampled at 448 different angular configurations; each one corresponds to a specific combination of θ i and θ s (which take the values 0 o , 10 o , 20 o , 30 o , 40 o , 50 o , 60 o and 70 o ), φ i (which take only the value 0 o ) and φ s (which take the values 0 o , 30 o , 60 o , 90 o , 120 o , 150 o , 180 o ).

Fig. 5 .Fig. 6 .
Fig. 5. Angular representation of the BRDF measurements at 500 nm in the incidence plane (0:180).Every plot shows BRDF values versus observation angles and corresponds to a different polar incidence angle (θ i ).

Fig. 8 .
Fig. 8. Relative contributions to the spectral variance of the first four principal components (half-plane 0:180) in a semilogarithmic plot.

Fig. 9 .
Fig. 9. Relative contributions to the spectral variance of the first four principal components (half-plane 30:150) in a semilogarithmic plot.
• Eigenspectrum PC2.It predominates over the other components, except for the case where specular reflection is involved.It represents the BRDF's diffuse component.