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

Michael W. Bishop; Samuel D. Butler; Michael A. Marciniak

doi:10.1117/1.OE.60.9.094103

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

Michael W. Bishop, Samuel D. Butler, Michael A. Marciniak

Author Affiliations +

Optical Engineering, Vol. 60, Issue 9, 094103 (September 2021). https://doi.org/10.1117/1.OE.60.9.094103

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%.

1. Background

Light reflected off surfaces can be modeled according to the direction of the incoming light, light direction upon reflection, and surface normal. Such a description is widely employed in scene generation, medical imaging, multi-layer color printing, resource monitoring, weather modeling, and paint development industries.¹^–⁷ For a complete dataset, one must make assumptions about the direction, magnitude, and wavelength dependency of light reflection off the observed material. Such a reflection description is given by a bidirectional reflectance distribution function (BRDF). The BRDF is intended to explain the reflection hemisphere only (the transmission hemisphere is handled separately when applicable but is beyond the scope of this paper), and to be independent of the overall brightness of the illumination. The BRDF, $f_{r}$ , is defined as the ratio of the reflected radiance ( $L_{s}$ ) to the incident irradiance ( $E_{i}$ ) as⁸

Eq. (1)

f_{r} ({\hat{ω}}_{i}, {\hat{ω}}_{s}, λ) = \frac{d L_{s} ({\hat{ω}}_{i}, {\hat{ω}}_{s}, λ)}{d E_{i} ({\hat{ω}}_{i}, λ)},

where

f_{r}

gives the reflectance per solid angle in a given direction in spherical coordinates in units of

s r^{- 1}

. The wavelength,

λ

, is mentioned here for completeness, as the surface parameters 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

{\hat{ω}}_{i}

and

{\hat{ω}}_{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⁹

Eq. (2)

L_{s} = \int_{Ω_{+}} f_{r} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) L_{i} ({\hat{ω}}_{i}) \cos θ_{i} d {\hat{ω}}_{i},

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¹⁰

Eq. (3)

f_{r} = \frac{L_{s}}{E_{i}} = \frac{ρ}{π},

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.⁷

BRDF models may further separate the directional volumetric term from the surface and Lambertian terms leading to a universal BRDF form given as⁷

Eq. (4)

f_{r} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) = ρ_{s} S ({\hat{ω}}_{i}, {\hat{ω}}_{s}) + ρ_{v} V_{d} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) + \frac{ρ_{d}}{π},

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.¹^,⁶^,⁷ 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,⁷

Eq. (5)

f_{r} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) = ρ_{s} P ({\hat{ω}}_{i}, {\hat{ω}}_{s}) D ({\hat{ω}}_{h}) F (θ_{d}) G ({\hat{ω}}_{i}, {\hat{ω}}_{s}) σ (θ_{i}, θ_{s}) + ρ_{v} V_{d} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) + \frac{ρ_{d}}{π},

where the surface reflection function

S

from Eq. (4) is the product of

P ({\hat{ω}}_{i}, {\hat{ω}}_{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 ({\hat{ω}}_{h})

is the microsurface normal distribution,

F (θ_{d})

is the unpolarized form of Fresnel reflection term,

G ({\hat{ω}}_{i}, {\hat{ω}}_{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.³^–⁵^,⁷^,¹¹

The Cook–Torrance model is a common microfacet BRDF model and takes the form,⁷^,¹²

Eq. (6)

f_{r} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) = 4 ρ_{s} D_{b} (θ_{h}) F (θ_{d}) G_{c} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) σ (θ_{i}, θ_{s}) + \frac{ρ_{d}}{π} .

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

Eq. (7)

σ (θ_{i}, θ_{s}) = \frac{1}{4 \cos θ_{i} \cos θ_{s}},

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

Eq. (8)

D_{b} (θ_{h}) = \frac{1}{π m^{2} \cos^{4} θ_{h}} \exp [- {(\frac{\tan θ_{h}}{m})}^{2}],

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¹³).

θ_{h}

is the scattered angle rotated in the microfacet normal geometry.

F (θ_{d})

is the unpolarized form of Fresnel reflection and

G_{c} ({\hat{ω}}_{i}, {\hat{ω}}_{s})

is the Blinn geometric function derived from shadowing and obscuration off V-shaped grooves is given as¹⁴

Eq. (9)

G_{c} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) = \min [1, (\frac{2 \cos θ_{h} \cos θ_{s}}{\cos θ_{d}}), (\frac{2 \cos θ_{h} \cos θ_{i}}{\cos θ_{d}})],

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¹⁵

Eq. (10)

L_{source} ({\hat{ω}}_{s}, λ) = L_{self} ({\hat{ω}}_{s}, λ) + \int_{Ω_{+}} f_{r} ({\hat{ω}}_{i}, {\hat{ω}}_{s}, λ) L_{i} ({\hat{ω}}_{i}, λ) \cos θ_{i} d {\hat{ω}}_{i},

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.¹⁶

2. Motivation

One standard method for BRDF model fitting to scatterometry measurements uses a non-linear least squares curve fit (lsqcurvefit) function in MATLAB^®.²^,⁶^,¹⁷^,¹⁸ 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.¹⁹ 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.²⁰^–²² 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.²⁰

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.

Fig. 1

(a), (b), 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.

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.⁷

3. 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.⁶^,¹⁸ 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 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

Eq. (11)

f_{v e} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) = ρ_{v} D_{b} (θ_{h e}) F (θ_{d e}),

where

θ_{d e}

and

θ_{h e}

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}

.

θ_{d e}

and

θ_{h e}

are then given as

Eq. (12)

θ_{d e} = \frac{1}{2} \cos^{- 1} [\cos θ_{i} \cos θ_{s} + \sin θ_{i} \sin θ_{s} \cos (ϕ - π)]

and

Eq. (13)

θ_{h e} = \cos^{- 1} [\frac{\cos θ_{i} + \cos θ_{s}}{2 \cos θ_{d e}}]

thereby defining the full model as

Eq. (14)

f_{m} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) = ρ_{s} σ (θ_{i}, θ_{s}) D_{b} (θ_{h}) F (θ_{d}) G_{c} ({\hat{ω}}_{i}, {\hat{ω}}_{s}) + \frac{ρ_{d}}{π} + ρ_{v} D_{b} (θ_{h e}) F (θ_{d e}) .

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.

Table 1

Upper and lower parameter bounds used in the recursive fitting algorithm.

Parameter	Symbol	Lower bound	Upper bound
Diffuse fit	$ρ_{d}$	0	1
Surface fit	$ρ_{s}$	0	100
Volume fit	$ρ_{v}$	0	100
Facet slope dist.	m	0.00001	10
Real index	n	1	100
Imaginary index	$k$	0	100

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.⁶^,⁷^,¹⁸ Fitting performance was based on the square of the mean standard error ( ${MSE}^{2}$ ) over all observation angles. MSE is calculated as

Eq. (15)

MSE = \frac{1}{n} ‖ \ln (\vec{x}) - \ln (\vec{f}) ‖,

where

\vec{x}

is the measured BRDF data, and

\vec{f}

is the BRDF model given in Eq. (14). To specify error within an observation sub-region, the Euclidean norm with

n

elements given as

Eq. (16)

‖ \ln (\vec{x}) - \ln (\vec{f}) ‖ = \sqrt{\sum_{k = 1}^{n} {| \ln ({\vec{x}}_{k}) - \ln ({\vec{f}}_{k}) |}^{2}},

can be squared allowing error sub-regions to be conveniently expanded and described with the

{MSE}^{2}

as

Eq. (17)

{(MSE)}^{2} = \frac{1}{n^{2}} (\sum_{k = 1}^{n} {| \ln ({\vec{x}}_{k}) - \ln ({\vec{f}}_{k}) |}^{2}) .

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.

Fig. 2

Algorithm flow chart to determine microfacet BRDF parameters from an array of local minima.

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.

4. 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.

Fig. 3

(a) Idealized test data 1 without error and (b) final test data 1 with random error ranging between 0% and 10% applied to the ideal data.

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.

Fig. 4

Plot (a) shows the fitting algorithm finds multiple local fits to test dataset 1. The relative global best fit is identified from the local minima and extracted as the global solution as shown in (b).

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.

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.

Parameter	Symbol	Min converged	Max converged	Best fit	Actual
Diffuse fit	$ρ_{d}$	$2.29 \times 10^{- 14}$	0.109	0.100	0.100
Surface fit	$ρ_{s}$	$2.498 \times 10^{- 14}$	47.51	11.27	10.00
Volume fit	$ρ_{v}$	0.0052	72.5	0.01	0.01
Facet slope dist.	$m$	0.10	9.995	0.10	0.10
Real index	$n$	1.00	10.5	2.94	3.00
Imaginary index	$k$	0.000118	5.16	0.900	1.00

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.

Fig. 5

(a) Idealized test data 2 without error and (b) final test data 2 with random error ranging between 0% and 10% applied to the ideal data.

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.

Fig. 6

Plot (a) shows the fitting algorithm finds multiple local fits (solid lines) to test dataset 2 (asterisks). The relative global best fit is identified from the local minima and extracted as the global solution shown in (b).

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.

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.

Fig. 8

Plot (a) shows the fitting algorithm finds multiple local fits (solid line) to test dataset 3 (asterisks). The relative global best fit is identified from the local minima and extracted as the global solution as shown in (b).

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.

5. 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.

Fig. 9

Recursive optimization identified four local minima, shown in solid blue and solid black lines, (a) for PNT65VIS data shown in plot and (b) the extracted relative global best fit solution in plot. About 76% improvement potential is observed.

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.

Fig. 10

Recursive optimization identified six local minima, shown in solid blue and solid black lines, (a) for PNT66MWIR and (b) the extracted relative global best fit solution. About 87% improvement potential is observed.

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.

Fig. 11

Recursive optimization identified 49 local minima, shown in solid blue and solid black lines, for (a) PNT66UV and (b) the extracted relative global best fit solution. About 87% improvement potential is observed.

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.

Fig. 12

Recursive optimization identified three unique local minima, shown in solid blue and solid black lines, for (a) STD00698 and (b) the extracted relative global best fit solution. About 18% improvement potential is observed.

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).

Table 3

An average of seven unique local minima exist for a given set of initial conditions using the standard method. Recursive optimization has the potential to improve fit quality an average of 71% over the high error locally minima. Incident angles in the BRDF data obtained from the Air Force Research Laboratory and used for the study are also listed.

Material	θi (Deg.)	Unique local minima	High error	Low error	Improvement potential
65 MWIR	30, 60	9	2.05E−04	5.73E−05	0.72
65 LWIR	30, 60	3	1.80E−03	5.81E−04	0.68
65 NIR	30, 60	2	8.57E−06	3.04E−06	0.64
65 UV	30, 60	2	3.11E−05	2.09E−05	0.33
65 VIS	30, 60	4	1.58E−05	3.75E−06	0.76
66 MWIR	30, 60	6	1.40E−03	1.83E−04	0.87
66 NIR	30, 60	8	6.27E−04	4.47E−05	0.93
66 UV	30, 60	49	3.69E−04	4.73E−05	0.87
36375 NIR	40, 60, 80	4	8.15E−05	1.28E−05	0.84
36495 VIS	30, 60, 85	6	6.20E−05	4.22E−05	0.32
01006 NIR	20, 60, 75	3	4.50E−03	2.07E−04	0.95
01014 NIR	20, 60, 75	6	1.24E−02	3.40E−03	0.73
STD00696	0, 20, 40, 60, 80	2	6.40E−04	4.72E−05	0.93
STD00698	0, 20, 40, 60, 80	3	1.28E−04	1.05E−04	0.18
STD00699	0, 20, 40, 60, 80	3	1.90E−03	1.14E−04	0.94

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.

Fig. 13

On average recursive optimization provides 71% fit quality improvement over the high error local minima.

6. 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 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!

Acknowledgments

This material is based upon work supported by the Air Force Office of Scientific Research under award number F4FGA09014J002. The authors would also like to acknowledge and thank Joe Costantino of the Air Force Research Laboratory, Materials and Manufacturing Directorate for selecting and providing the BRDF data used in this work.

References

1.

C. Schlick, “An inexpensive BRDF model for physically-based rendering,” Comput. Graphics Forum, 13 (3), 233 –246 (1994). https://doi.org/10.1111/1467-8659.1330233 CGFODY 0167-7055 Google Scholar

2.

C. B. Madsen, B. K. Mortensen and J. R. Andersen, “Real-time view-dependent visualization of real world glossy surfaces,” in GRAPP, 231 –240 (2008). Google Scholar

3.

W. C. Snyder and Z. Wan, “BRDF models to predict spectral reflectance and emissivity in the thermal infrared,” IEEE Trans. Geosci. Remote Sens., 36 (1), 214 –225 (1998). https://doi.org/10.1109/36.655331 IGRSD2 0196-2892 Google Scholar

4.

R. Latifovic, J. Cihlar and J. Chen, “A comparison of BRDF models for the normalization of satellite optical data to a standard sun-target-sensor geometry,” IEEE Trans. Geosci. Remote Sens., 41 (8), 1889 –1898 (2003). https://doi.org/10.1109/TGRS.2003.811557 IGRSD2 0196-2892 Google Scholar

5.

S. H. Westin, H. Li and K. E. Torrance, “A comparison of four BRDF models,” in Eurographics Symp. Rendering, 1 –10 (2004). Google Scholar

6.

B. E. Ewing, S. D. Butler and M. A. Marciniak, “Improved grazing angle bidirectional reflectance distribution function model using Rayleigh–Rice polarization factor and adaptive microfacet distribution function,” Opt. Eng., 57 (10), 105102 (2018). https://doi.org/10.1117/1.OE.57.10.105102 Google Scholar

7.

S. D. Butler and M. A. Marciniak, “Robust categorization of microfacet BRDF models to enable flexible application-specific BRDF adaptation,” Proc. SPIE, 9205 920506 (2014). https://doi.org/10.1117/12.2061134 Google Scholar

8.

F. E. Nicodemus et al., “Geometrical considerations and nomenclature for reflectance,” (1977). Google Scholar

9.

J. T. Kajiya, “The rendering equation,” ACM SIGGRAPH, 20 (4), 143 –150 (1986). https://doi.org/10.1145/15886.15902 Google Scholar

10.

R. W. Boyd, Radiometry and the Detection of Optical Radiation, John Wiley and Sons, New York (1983). Google Scholar

11.

W. Wanner, X. Li and A. Strahler, “On the derivation of kernels for kernel-driven models of bidirectional reflectance,” J. Geophys. Res.: Atmos., 100 (D10), 21077 –21089 (1995). https://doi.org/10.1029/95JD02371 Google Scholar

12.

R. L. Cook and K. E. Torrance, “A reflectance model for computer graphics,” ACM Trans. Graphics, 1 (1), 7 –24 (1982). https://doi.org/10.1145/357290.357293 ATGRDF 0730-0301 Google Scholar

13.

S. D. Butler, S. E. Nauyoks and M. A. Marciniak, “Comparison of microfacet BRDF model to modified Beckmann–Kirchhoff BRDF model for rough and smooth surfaces,” Opt. Express, 23 (2015). https://doi.org/10.1364/OE.23.029100 OPEXFF 1094-4087 Google Scholar

14.

J. F. Blinn, “Models of light reflection for computer synthesized pictures,” ACM SIGGRAPH, 11 (2), 192 –198 (1977). https://doi.org/10.1145/965141.563893 Google Scholar

15.

D. S. Immel, M. F. Cohen and D. P. Greenberg, “A radiosity method for non-diffuse environments,” ACM SIGGRAPH, 20 (4), 133 –142 (1986). https://doi.org/10.1145/15886.15901 Google Scholar

16.

F. M. Cady et al., “BRDF error analysis,” Proc. SPIE, 1165 154 –164 (1990). https://doi.org/10.1117/12.962845 PSISDG 0277-786X Google Scholar

17.

MATLAB, “Version 9.4.0.813654 (R2018a),” (2020). Google Scholar

18.

B. E. Ewing and S. D. Butler, “Grazing angle experimental analysis of modification to microfacet BRDF model for improved accuracy,” Proc. SPIE, 10644 106441I (2018). https://doi.org/10.1117/12.2303744 PSISDG 0277-786X Google Scholar

19.

Jr. J. E. Dennis, Nonlinear Least Squares, Academic Press(1977). Google Scholar

20.

Mathematica, “Version 12.1,” Champaign, Illinois (2020). Google Scholar

21.

I. G. Renhorn and G. D. Boreman, “Analytical fitting model for rough-surface BRDF,” Opt. Express, 16 (17), 12892 –12898 (2008). https://doi.org/10.1364/OE.16.012892 OPEXFF 1094-4087 Google Scholar

22.

D. Chan, S. Games and A. Blizzard, “Material advances in Call of Duty: WWII,” (2018). Google Scholar

Biography

Michael W. Bishop is a research physicist for the United States Space Force’s Space and Missile Systems Center in Los Angeles, California. He received his BS degree in physics from Northern Arizona University and MS degree in applied physics from the Air Force Institute of Technology. His research interests include bidirectional reflectance distribution functions, scatterometry, and space systems survivability. He is an SPIE member, has published one refereed paper, two additional technical papers, and presented at multiple scientific conferences.

Samuel D. Butler is a lieutenant colonel in the United States Air Force and an assistant professor at the Air Force Institute of Technology (AFIT). He received his BS degree in applied physics with a computer science emphasis from Brigham Young University, MS degree in applied physics, and PhD from AFIT. His research interests include bidirectional scatter. He has advised 7 MS and PhD students, and been published 19 technical papers. He is a member of SPIE.

Michael A. Marciniak is a professor at the Air Force Institute of Technology (AFIT). He received his BS degree in mathematics from St. Joseph’s College, India, BSEE from the University of Missouri-Columbia, and MSEE and PhD degrees from AFIT. His research interests include bidirectional scatter distribution and scatterometry of optical metasurfaces. He has advised 12 PhD dissertations and 55 MS theses, and published 40 refereed and 80 other technical papers. He is a senior member of SPIE.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Michael W. Bishop, Samuel D. Butler, and Michael A. Marciniak "Robust method of determining microfacet BRDF parameters in the presence of noise via recursive optimization," Optical Engineering 60(9), 094103 (6 September 2021). https://doi.org/10.1117/1.OE.60.9.094103

Received: 6 May 2021; Accepted: 20 August 2021; Published: 6 September 2021

Access the abstract

JOURNAL ARTICLE
15 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

KEYWORDS

Bidirectional reflectance transmission function

Data modeling

Optimization (mathematics)

Optical engineering

Reflection

Solids

Performance modeling

1.

Background

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

2.

Motivation

Fig. 1

3.

Methodology

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Table 1

Eq. (15)

Eq. (16)

Eq. (17)

Fig. 2

4.

Development of Methodology

Fig. 3

Fig. 4

Table 2

Fig. 5

Fig. 6

Fig. 7

Fig. 8

5.

Results with Measured Data

Fig. 9

Fig. 10

Fig. 11

Fig. 12

Table 3

Fig. 13

6.

Conclusion

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years