Robust method of determining microfacet BRDF parameters in the presence of noise via recursive optimization

Abstract. Accurate bidirectional reflectance distribution function (BRDF) models are essential for computer graphics and remote sensing performance. The popular microfacet class of BRDF models is geometric-optics-based and computationally inexpensive. Fitting microfacet models to scatterometry measurements is a common yet challenging requirement that can result in a model being fit as one of several unique local minima. Final model fit accuracy is therefore largely based on the quality of the initial parameter estimate. This makes for widely varying material parameter estimates and causes inconsistent performance comparisons across microfacet models, as will be shown with synthetic data. We proposed a recursive optimization method for accurate parameter determination. This method establishes an array of local minima best fits by initializing a fixed number of parameter conditions that span the parameter space. The identified solution associated with the best fit quality is extracted from the local array and stored as the relative global best fit. This method is first applied successfully to synthetic data, then it is applied to several materials and several illumination wavelengths. This method proves to reduce manual parameter adjustments, is equally weighted across incident angles, helps define parameter stability within a model, and consistently improves fit quality over the high-error local minimum best fit from lsqcurvefit by an average of 71%.

relative to the wavelength are important, but will not be further considered in this work since it deals with fitting BRDF data that is typically taken with a laser assumed to be at a single wavelength. Incident and scattered vectors are given in spherical coordinates defined with respect to the surface normal byω i andω s , respectively. This convenient description allows the development of reflectance models in terms of illumination and observation geometry. The BRDF is most commonly used in conjunction with the Rendering Equation. The reflection portion of the rendering equation at an intersection point with a material is 9 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 ; 1 1 6 ; 6 5 1 where L i is the incident radiance, L s is the reflected radiance, and Ω þ represents the reflection hemisphere.
For the simplest BRDF, assume a Lambertian surface. The BRDF, f r , is related to the unitless hemispherical reflectance, ρ, as 10 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 ; 1 1 6 ; 5 5 9 where Eq. (3) defines such a material to have a perfectly diffuse BRDF, independent of illumination or observation geometry. For monochromatic light incident on an idealized mirror surface, reflection behaves specularly and the angle of incidence equals the angle of reflection.
In practice, materials are not accurately described by a purely diffuse nor purely specular model, but can be well characterized as having both diffuse and specular components. These contributions can come from surface reflection or volumetric scatter within the material. 7 BRDF models may further separate the directional volumetric term from the surface and Lambertian terms leading to a universal BRDF form given as 7 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 4 1 9 f r ðω i ;ω s Þ ¼ ρ s Sðω i ;ω s Þ þ ρ v V d ðω i ;ω s Þ þ ρ d π ; (4) where ρ s is the surface fitting parameter, ρ v is now the directional volumetric fitting parameter, V d is the directional volumetric scatter function, and ρ d is the diffuse fitting parameter following Lambertian scatter. The last two terms constitute the volumetric scatter. One common model class that takes this approach is the microfacet model. 1,6,7 This class of geometric optics-based models is typically wavelength agnostic and dependent on the surface structure of the material, which is defined in the model as a distribution of infinitesimally small specularly reflecting facets. In general, microfacet models take the form, 7 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 ; 1 1 6 ; 2 9 6 f r ðω i ;ω s Þ ¼ ρ s Pðω i ;ω s ÞDðω h ÞFðθ d ÞGðω i ;ω s Þσðθ i ; θ s Þ þ ρ v V d ðω i ;ω s Þ þ ρ d π ; where the surface reflection function S from Eq. (4) is the product of Pðω i ;ω s Þ, a model specific prefactor accounting for terms not found in other models of similar form, θ d is the angle of incidence in microsurface orientation relative to the overall surface orientation, Dðω h Þ is the microsurface normal distribution, Fðθ d Þ is the unpolarized form of Fresnel reflection term, Gðω i ;ω s Þ is a geometric attenuation term, and σðθ i ; θ s Þ is a cross section conversion term. Wanner and others provide an overview of kernel based BRDF models and Butler provided a robust categorization of BRDF microfacet BRDF models in 2014. [3][4][5]7,11 The Cook-Torrance model is a common microfacet BRDF model and takes the form, 7,12 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 ; 1 1 6 ; 1 6 1 Referencing Eq. (5), prefactor P is equal to four in the Cook-Torrance model and can be ignored by scaling the surface parameter by 1∕4 to remove the prefactor entirely. The cross section conversion term is defined as which converts the surface scatter from spherical scatter, common in physics, to planar scatter more representative of a typical interface. D b ðθ h Þ is the Beckmann (Gaussian) distribution of microfacet surface-normal orientations relative to the macroscopic surface orientation with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 6 6 5 where m is the surface tangent angle parameter of the distribution (technically only true for rough surfaces; for polished surfaces it is more complex but the microfacet model begins to break down anyway in this limit 13 ). θ h is the scattered angle rotated in the microfacet normal geometry. Fðθ d Þ is the unpolarized form of Fresnel reflection and G c ðω i ;ω s Þ is the Blinn geometric function derived from shadowing and obscuration off V-shaped grooves is given as 14 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 5 6 0 where the first term is used when no geometric attenuation occurs; the second term occurs for surfaces displaying obscuration, meaning the microfacet's angle of reflection is large with respect to the macrosurface normal; and the third term describes surfaces displaying shadowing, meaning geometric attenuation occurs due to shadowing from microfacets at large angles with respect to the macrosurface normal. Cook-Torrance, like many other microfacet models, does not include a directional volumetric scatter term (ρ v ¼ 0). One can then apply the full BRDF using the full rendering equation, given as 15 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 4 3 1 to define the observed source radiance (L source ) in terms of both the self-emitted component (L self ) and the surrounding incident radiance (L i ) scattered according to the BRDF model. Accurately defining the behavior of the BRDF for a specific material is commonly dependent on accurate anchoring of model parameters to scatterometry measurements. 16

Motivation
One standard method for BRDF model fitting to scatterometry measurements uses a non-linear least squares curve fit (lsqcurvefit) function in MATLAB ® . 2,6,17,18 This function reads in an initial parameter guess and finds final parameters that correspond to the best non-linear fit to the specified model for the provided data. The number of iterations in the optimization depends on the fitting tolerance. Fit performance is partially based on the accuracy of the initial guess, which can be done with parameter tweaking (i.e., human-assisted adjustment of parameters that result in substantially better fits) to manually fit a model to each dataset. However, this method will only converge to one local minimum for that particular set of initial conditions. 19 Two independent initial guesses may converge to two unique local minima that provide differing fit quality. Additionally, manually fitting a model can bias the initial guess to a particular incident angle, making this method largely ineffective for BRDF measurements across multiple incident angles. A similar nonlinear fit function is commonly implemented when fitting BRDF parameters in Mathematica ® software. [20][21][22] The FindFit and NonlinearModelFit functions minimize a norm of the residual function to give estimated model parameters. As is often the case in nonlinear fitting techniques, there may be several local minima and determining the global minimum can be difficult and computationally expensive. 20 Figure 1 shows example scatterometry data of a material illuminated at incident angles of 30 deg and 60 deg. The data are fit with the Cook-Torrance model given in Eq. (6) using three independent initial parameter estimates and the common lsqcurvefit approach. An ideal fitting routine converges to one solution with globally minimized fit error and is independent of the quality of the initial parameter estimate. However, Fig. 1 shows that using the standard fit approach converges to three different local minima "best" fit solutions, each with unique model parameter estimates resulting in dramatically different fits.
Accordingly, it is desirable to reduce manual manipulation and obtain a global minimum for the optimal fit. Doing so can help determine a true best fit for the model and improve upon the fit performance of a model. This paper compares the standard fit routine with a recursive optimization method tailored to interrogate an array of unique local minima for a given model and dataset. A modified Cook-Torrance model with six fit parameters (to be described in the next section) is used to compare fit quality results using the standard fit approach and the recursive fit approach to 15 datasets, although the number of parameters to fit does vary with choice of BRDF model; some common models with number of parameters are given in the referenced source. 7

Methodology
Fifteen measurements of nine samples were made at incident angles ranging from θ i ¼ 0 deg to θ i ¼ 85 deg. Five different illumination wavelengths, one at a time, provided spectral sampling at 0.325 μm (UV), 0.6328 μm (VIS), 1.06 μm (NIR), 3.39 μm (MWIR), and 10.6 μm (LWIR); see Table 3. Not all samples were illuminated at all wavelengths. Measurements for this work were conducted at the optical measurements facility (OMF) at the Air Force Research Laboratory, Materials and Manufacturing Directorate. Samples are primarily paint coatings with glossy (specular) or matte (diffuse) visual characteristics. Samples were chosen to exhibit scatter with a variety of volumetric scatter components. 6,18 The incident angles at which data were taken are shown in Table 3; for each incident angle, several scattered angles that are in the plane of incidence were , and (c) each show BRDF model fit quality using the standard approach on identical data but three independent initial parameter estimates. Standard curve fitting converged to three unique local minima. The set shows that only local optimization occurs and final fit quality is based on the quality of the initial parameter estimate.
sampled to obtain the BRDF value. The data obtained was polarimetric, but in this work we are using an unpolarized BRDF by averaging the s and p polarization BRDF components.
The recursive fit method is anchored to data using a modified Cook-Torrance model. The Cook-Torrance model is commonly given without a directional volumetric scatter term; accordingly, this study includes a semi-empirical directional volume term empirically based on subsurface scatter given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 6 6 3 where θ de and θ he are modified incident and scattered angles rotated for microsurface coordinates to aid in forming a volumetric backscatter lobe. These variables are related to their more commonly used forward scatter variables θ h and θ d but are at ϕ ¼ ϕ i instead of ϕ ¼ ϕ i þ π, resulting in a secondary backscatter peak centered about ϕ s ¼ ϕ i that is observed in the data obtained from the Air Force Research Laboratory. Given Δϕ is the difference between the incident, ϕ i , and scattered, ϕ s , azimuthal angles, one may conveniently allow ϕ i to be zero for all cases so that Δϕ is equivalently expressed as ϕ s . θ de and θ he are then given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 5 4 9 and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 4 9 9 thereby defining the full model as For each dataset, an array of unique best fit local minima is determined by initializing 250 sets of randomized parameter values that span the defined parameter bounds given in Table 1.
A non-linear least squares fit is applied for each of the initial conditions, and each converges to its respective locally optimized fit. Each incident angle is fit simultaneously and the error is calculated for each local minimum.
The hardware used for the nonlinear curve fit was a standard off-the-shelf MacBook Pro laptop purchased in 2018. Run time varied with dataset being analyzed, but did not exceed 10 minutes. No attempts were made to optimize the code for speed, and the code was not written to take advantage of parallel processing, so this execution time could be greatly enhanced with some effort, as the solution presented is trivially parallel.
Fits are calculated by fitting to the natural log of the measured BRDF data. BRDF often varies several orders of magnitude and drops off very quickly with deviations from the specular peak; therefore, the logarithmic method allows for an emphasized fit over all observation angles, not just the model's forward specular peak. 6,7,18 Fitting performance was based on the square of the mean standard error (MSE 2 ) over all observation angles. MSE is calculated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 6 ; 7 3 5 MSE ¼ 1 n k lnðxÞ − lnðfÞk; (15) wherex is the measured BRDF data, andf is the BRDF model given in Eq. (14). To specify error within an observation sub-region, the Euclidean norm with n elements given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 6 7 6 can be squared allowing error sub-regions to be conveniently expanded and described with the MSE 2 as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 6 ; 6 0 0 This study evaluates MSE 2 over the full observation region.
From the array of local minima, the total number of unique local minima is determined. Randomly initializing parameters provides unique unbiased sampling of the parameter space allowing insight to the various local minima that could be observed if one were to use the standard single fit method. The single fit providing the lowest fit error is taken as the relative global best fit. The relative global best fit parameter set is stored and acts to define the final solution for comparison across materials. A summarized graphical outline of the described parameter fitting algorithm is shown below in Fig. 2.
To compare effectiveness of the recursive method across datasets, multiple metrics are considered. First, the total number of unique local minima is identified, indicating the confidence of any one solution if the standard single fit method were utilized. Second, the fit quality of the local fit with the highest MSE 2 is compared to the fit quality of the local fit with the lowest MSE 2 . This identifies the fit variance for each material. Finally, an improvement potential is determined given as the percentage of reduced error between the best and the worst fitting local minima. This is the realized improvement if one were to fit data using the standard single fit method with convergence to the worst local minimum, compared to using the recursive optimization method, which identifies the relative global best fit solution.

Development of Methodology
Prior to application on real BRDF measurement data, the fitting algorithm was validated with three different test datasets that were created from idealized BRDF's. Each test simulated observation of −85 deg < θ s < 85 deg. The BRDF used for test data 1, shown on the left plot in Fig. 3, has incident light at θ i ¼ 45°and known parameters of ρ d ¼ 0.10, ρ s ¼ 10.00, ρ v ¼ 0.01, m ¼ 0.1, n ¼ 3.00, and k ¼ 1.00. From the ideal model, random error ranging between 0% and 10% was placed on the ideal model at every 1 deg of observation. The resulting test data 1 is shown on the right in Fig. 3. Test data 1 was fit to the modified Cook-Torrance model using the recursive optimization algorithm where 13 local minima were identified, 12 of which are very closely spaced and difficult to distinguish visually but are as a whole, the set of clearly poor fits shown on the left plot in Fig. 4. Interestingly, this shows that the naive technique of using the previous best fit as the initial guess for a subsequent curve fit is not likely to succeed, as there is an apparent cluster of local minima within the vicinity of these 12 poor quality fits, whereas the best fit contains radically different fit parameters not likely to be converged on if such an approach is used. The parameters associated with the best local fits were stored and used for the initial parameters in a final fit. The right plot in Fig. 4 shows the fit solution extracted from the set of local minima. The final fit converged to parameter values of ρ d ¼ 0.10, ρ s ¼ 11.27, ρ v ¼ 0.01, m ¼ 0.1, n ¼ 2.94, and k ¼ 0.90. Large changes in ρ s result in small changes in the BRDF for this particular test BRDF configuration; therefore, parameter ρ s varied more significantly from the known value with little impact to the final solution. All but the complex component of the index of refraction were fit to within 10% of the known parameter value, whereas ρ d , ρ v , and m converged to the exact parameter values in the known BRDF.
When examining the entire table of fits for this set of simulated data, the minimum and maximum parameters are summarized in Table 2. In each case, at least 1 local minimum converged to the Min Converged and Max Converged parameters. The Best Fit parameters are from a single local minimum, which contained the minimum residual error overall out of all 250 curve fits and was quite different in value from the other 12 local minima. Of substantial note is the extreme range in parameters that resulted in a best fit. It is also worth noting that the m parameter (Facet slope dist.) is used in a Gaussian distribution, and thus its variation has substantial effect on the BRDF fit quality.
The fitting algorithm was similarly validated using test data 2 shown in Fig. 5. The BRDF used for test data 2, shown on the left plot in Fig. 5, is configured for incident light at θ i ¼ 80 deg to examine the algorithm's resilience with near grazing illumination conditions. The known parameters are ρ d ¼ 0.20, ρ s ¼ 1.50, ρ v ¼ 0.02, m ¼ 0.15, n ¼ 1.15, and k ¼ 0.9. Random error was again placed on the ideal model resulting in test dataset 2 shown on the right in Fig. 5.
Test data 2 was fit to the modified Cook-Torrance model using the recursive optimization and found six unique local minima. Each is shown in blue in the left plot of Fig. 6. The right plot in Fig. 6 shows the relative global minimum solution extracted from the set of local minima.
The dataset 2 solution converged to parameter values of ρ d ¼ 0.20, ρ s ¼ 3.37, ρ v ¼ 0.02, m ¼ 0.15, n ¼ 1.01, and k ¼ 0.87. Again, large changes in ρ s result in small changes in the BRDF, therefore parameter ρ s varied more significantly from the known value with little impact to the fit quality of the final solution. All other parameters converged to within 15% of the known parameter value, whereas ρ d , ρ v , and m converged to the exact known parameter values.
Test data 3 used BRDF shown on the left plot in Fig. 7. This set is configured for incident light at θ i ¼ 30°to examine the algorithm's resilience with near normal illumination conditions. The parameters are also adjusted to give a low volumetric component and are known to be ρ d ¼ 0.15, ρ s ¼ 1.00, ρ v ¼ 0.001, m ¼ 0.005, n ¼ 1.10, and k ¼ 0.86. Random error was again placed on the ideal model resulting in test dataset 3 shown on the right in Fig. 7. Table 2 This table shows max, min, best fit, and actual values of parameters used in simulated data for test data set 1. The max and min values were all returned as a supposed "best fit" from a single iteration, whereas the best fit values are from a single iteration using parameters from the overall minimum error.     Fig. 7 Idealized test data 3 without error (a) and final test data 3 with random error ranging between 0% and 10% applied to the ideal data (b). Test data 3 was fit to the modified Cook-Torrance model using the recursive optimization and found four unique local minima. This dataset has few datapoints at the forward and backscatter specular peaks causing three of the four local minima to ignore these as outliers. Accordingly, these three local minima converge to a Lambertian BRDF, which is clearly incorrect. Recursive optimization however identified one specular description that far outperforms the fit quality of the other local minima. Each unique local minimum is shown as a solid blue line in the left plot of Fig. 8. The right plot in Fig. 8 shows the final solution extracted from the set of local minima.
The dataset 3 solution converged to ρ d ¼ 0.15, ρ s ¼ 2.37, ρ v ¼ 0.001, m ¼ 0.005, n ¼ 1.02, and κ ¼ 0.92. Consistent agreement with the known parameter values is again observed with this third dataset. All but ρ s converged to within 9% of the known values, whereas ρ d , ρ v , and m again converged to the exact parameter values. Consistently accurate m parameter determination is significant for material to material identification. With validation of the methodology using test datasets, the recursive optimization approach was applied to fit the modified Cook-Torrance model to each of the 15 real datasets described at the beginning of this section.

Results with Measured Data
Following the recursive optimization approach, each of the 15 datasets were fit with the modified Cook-Torrance model. Results show that the standard single fit approach may converge to one of multiple possible local minima, that is, there were not any instances of consistent convergence to only one fit solution. This suggests that the standard single fit method should not be used for comparing multiple models as the quality of the local fit will be unclear. Figure 9 shows the BRDF model fit to PNT65VIS data and the four unique local minima identified. Recursive optimization extracts the relative global solution with the highest fit quality and is shown on the right in Fig. 9. Depending on the initial parameter estimates, the standard non-linear least squares fit will converge to one of four locally optimized minima, three of which fail to represent the specular backscatter peak. Accordingly, recursive optimization improves fit quality 76% over the local fit with the poorest fit quality, defined as the high error fit. Figure 10 shows the BRDF model fit to PNT66MWIR data and the six unique local minima identified. Five minima agree at near-normal observations but provide unique fits at near-grazing observations. Additionally, five of the six local minima underestimate the forward scattered BRDF and underestimate the backscatter BRDF. Recursive optimization extracts the relative global fit with the highest fit quality and is shown on the right in Fig. 10. Recursive optimization shows improvement potential of 87% over the high error local fit. Figure 11 shows the BRDF model fit to PNT66UV data and 49 unique local minima identified. Each local minimum is determined by a minimum of 1% difference in MSE 2 and confirmation of unique parameter values found in that fit. This example again demonstrates that fit  quality is likely to vary for this material when the standard single fit method is used. Additionally, establishing parameter values to associate with this material will prove difficult without identification of possible local minima. The majority of local minima at θ i ¼ 30°underestimate the BRDF at near normal observations and overestimate the BRDF at grazing observations. Recursive optimization extracts the relative global solution with the highest fit quality and is shown on the right in Fig. 11. Recursive optimization shows improvement potential of 87% over the high error local fit.   Figure 12 shows the BRDF model fit to STD00698 data and three unique local minima identified. All local minima agreed well at near normal observations but differed at grazing observations. The recursive optimization solution shown on the right in Fig. 11 provides improvement potential of 18% over the high error local fit.
Results for the remaining 11 datasets are summarized in Table 3. The total number of unique local minima identified is provided along with the MSE 2 of the best fit (low error) solution and the local fit with the poorest fit (high error).
In all datasets multiple unique local minima were identified. PNT65 MWIR was not shown here but provided the second highest number of unique local minima at nine. Additionally, PNT01006 NIR showed the largest improvement potential of 95%. A visual comparison of improvement potential across all datasets is provided in Fig. 13.

Conclusion
An unsupervised, yet reliable recursive optimization routine was proposed as an improved approach for consistent determination of BRDF parameter values. This method improves upon the standard single fit approach that is susceptible to sub-optimal fit conclusions due to local minimization constraints that may result with a naive unsupervised routine, as was demonstrated in this paper using synthetic data. A modified Cook-Torrance model with a semi-empirical directional volume term was used with the recursive optimization method and validated against three synthetic test datasets, as well as implemented on scatterometry measurements of several materials. On average, seven unique local minima were identified in the datasets and the relative global best fit solution was determined. The identified best fit using the algorithm proposed here provided 71% improvement potential over the standard fit approach.
This method proves to be stable as improvements were observed over all provided material types and illumination wavelengths. Recursive optimization does not require manual tweaking of the model parameters prior to fitting nor after a local minimum is found. The method is equally weighted across incident angles since manual fitting for initial parameter estimates is not needed. Additionally, recursive optimization helps define fit stability for any one particular model as it provides a range of locally optimized minima within the defined parameter space. Finally, confidence in fit quality is improved allowing for accurate model to model comparison.
Although we determined 250 iterations were best for this particular BRDF model employed, the reader should be aware that we reached this number by adding iterations until new unique local minima were no longer obtained, and that result depends upon the BRDF model being used. In our case, we found the BRDF model in this paper had as many as 49 unique local Fit Improvement Potential (%) Fig. 13 On average recursive optimization provides 71% fit quality improvement over the high error local minima.
minima (see Table 3). The same procedure will otherwise work for other models, but to obtain the number of iterations for another BRDF model, one should iterate repeatedly until no new unique local minima are found. However, this work has demonstrated that one should not use a single iteration of lsqcurvefit to determine BRDF parameters, as it may fail even with synthetic data when a small amount of noise is added!