1. Introduction

Rough surfaces are all around us. When we generate realistic images, the task of visualizing them arises. Fig. 1 shows examples of such visualizations created by us: frosted glass with objects visible through it, and a rough car interior panel. Fortunately, in many lighting simulations and optical design tasks it is sufficient and more effective to replace real geometry of rough surface by a smooth surface with certain optical characteristics. The definition of scattering properties for a smooth boundary between two media is a simple task and the light scattering can be easily simulated using Snell’s law of refraction and reflection. However, in case of rough boundary the definition of light scattering is more complex and can be expressed via Bidirectional Scattering Distribution Function (BSDF). BSDF determines output angular light transformation (refraction and reflection) in the dependence of input light conditions, angles of the input light.

 

Figure 1: Examples of rough surface visualization

 

In simple cases, when light scattering by whole plate is only important, the direct BSDF measurements may be sufficient. The ordinary way of BSDF measurements for the rough surface is presented on Fig. 2a. The sample, one side of which is rough, is illuminated with a parallel light beam under the specific incident directions, then an angular light distribution of transmitted light (BTDF - Bi-Directional Transmittance Distribution Function) and reflected light (BRDF – Bi-Directional Reflectance Distribution Function) is measured. In other words, such BSDF measured for the whole sample works in cases when we can ignore object thickness. The examples of such ordinary BSDF applications (Fig. 2b) may be various diffuse films, thin plates, layers.

 

Figure 2: BSDF application in the simplified “one-sheet” and more “solid” models

 

However, in plenty of cases the direct usage of measured BSDF is impossible. As an example we can consider a light guiding plate with rough surface (Fig. 2c). Correct simulation of light propagation in this optical system requires BSDF from each side of the rough surface that includes BSDF from the material side. The BSDF measurement from material side is impossible or very expensive because we cannot place the light source and detector inside of the material. Another problem is the significant inaccuracy of BSDF measurements for big illumination angles because of light leakage inside of measured samples, shadowing of sample illumination and some other reasons.

The mentioned problems related to BSDF measurement result in the development of many approaches and methods for BSDF reconstruction. One of the main purposes of this paper is an analysis and verification of the popular methods of BSDF reconstruction. The paper contains an overview of most prominent approaches and their comparison done on the base of real measured samples with rough surfaces.

2. Overview of BSDF reconstruction methods

Generally, methods of BSDF reconstruction can be divided into two main groups:

1. Analytical methods. The analytical methods are based on the theory of physics (optics) or empirical formulas. The methods represent rough surface models published by Ward, Cook-Torrance, Phong, etc. The main advantages of analytical approaches are high efficiency because analytical solutions are fast to calculate. This is important because optimization procedure is typically used to get parameters of analytical functions describing required BSDF shape. The disadvantage of the approaches is their approximation. They use approximate algorithms to describe complex optical effects like a masking or shadowing of the incident light illuminating of the rough surface (Fig. 3a, 3b) and a interreflection of light on rough profile (Fig. 3c). During BSDF reconstruction this can introduce noticeable inaccuracy for surfaces with big roughness.

2. Numerical methods. These approaches are based on simulation of light propagation through models of rough surfaces. In the given paper two main numerical approaches are considered which are based on distribution of microfacet normals or heights. These approaches are more correct than analytical ones from viewpoint of optical theory but require noticeable calculation resources.

 

Figure 3: Masking-shadowing and multiple reflections of light on rough surface

 

The described classification is not the only one. For example, the methods for BSDF reconstruction can be also divided depending on what optics, geometrical (ray) or wave, is applied. In our work we investigated all these groups of methods.

2.1. Analytical methods of BSDF reconstruction

Most of the methods for defining a light scattering (BSDF) through a rough surface are based on the “microfacet” model. In the case of the “microfacet” model rough surface with complex geometry is presented with a set of flat smooth surfaces (micro facets), see Fig. 4. When boundary (microfacet) is smooth the transmission, reflection can be easily simulated using Snell’s and Fresnel laws of refraction and reflection. So it is possible to calculate a general light scattering through the rough surface knowing general density distribution of microfacet slopes or their normals. One of the earlier attempts to model light reflection from a rough surface is described in [1]. It was restricted with reflection light component only but it was a basis for developing one of the more well-known microfacet models introduced by Cook and Torrance [2]. A lot of different modifications of the microfacet model have been developed at that time [3-5]. The next developments are related to extensions of reflection “microfacet” models with the support of anisotropy, sampling with correct weights, application of Backmann distribution [6, 7], and development of alternative sampling methods with fitted separate approximations [8].

 

Figure 4: Microfacet presentation of a rough surface

 

Shlick [9] develops more simple approximation to the Cook-Torrance model with the help of rational approximations with the application of the Fresnel formula widely adopted nowadays. Ashikhmin and Shirley [10] introduced an anisotropic reflection model on the base of Phong microfacet distribution including correct importance sampling. Then an energy-conserving reflection model [11] is introduced. It is derived from arbitrary microfacet distributions, though this formulation involves numerically estimating integrals without closed-form solutions.

Microfacet models are widely used in computer graphics; experimental data appear for verification of scattering models. For example, different models of BRDF reconstruction (“Ward”, “Ward-Duer”, “Blin-Phong”, “Cook-Torrance”, “Lafortune et al”, “He et al”, “Ashikhmin-Shirley”) are compared with real measurements in [12]. A set of development is related to the derivation of the refraction part of scattered light [13]. There are investigations that take into account thin effects such as the shadowing masking, multi-inter-reflections on elements of rough surface, using of importance sampling [2, 14, 15].

The reflection models based on wave optics are proposed in [16]. The method can simulate a wider range of surface effects than microfacet models. However, wave approaches are much more expensive to calculate and as a rule very approximate in support of thin effects as a multi reflection on profile with rough surfaces.

The numerical simulations of transmissions models are performed in [17-20, 37]. The “GGX” microfacet model was introduced in [21]. It is an improved variant of the Cook-Torrance microfacet model supporting reflection as well as refraction and shadowing-masking. The [21] work contains numerical data comparison of different analytical models and demonstrates a lot of advantages relative to other analytical methods of BSDF reconstruction of a rough surface. The “GGX” model is considered one of the most accurate, flexible and wide-used analytical approaches. It supports both reflection and refraction components, masking-shadowing, and importance sampling and shows more accurate output in relation to the Cook-Torrance model [21].

So, the “GGX” model is selected for examination in our paper as representative of the analytical group of methods. Typically, an analytical model is represented with two base functions. The first function, denoted as D(m), is a microfacet distribution function. It describes the statistical distribution of surface normal m over microsurfaces. The second bi-directional function, denoted as G(i, o, m), describes what fraction of the microsurface with normal m is visible in both directions i and o (Fig. 5). Typically, the shadowing-masking function has relatively little influence on the shape of the BSDF except for near grazing angles or for very rough surfaces but is needed to maintain energy conservation.

 

Figure 5: Micro vs. macro surface

 

The bi-directional function G(i, o, m ) can be approximated with two omnidirectional shadowing terms:

(1)

where G1 is derived from the microfacet distribution D as described in [14, 15].

We used the “GGX” model with the following microfacet distribution and masking shadowing function D(m), parameter specifies surface roughness:

(2)

and omnidirectional masking-shadowing function:

(3)

where is the angle between m and n, between v and n, and is the positive characteristic function (which equals one if > 0 and zero if ≤ 0). v equals either to i or o vectors (Fig. 5). Note the function is rather similar to the well-known Beckmann distribution used in the Cook-Torrance model.

So the process of BSDF reconstruction consists of the definition of the parameter - degree of surface microroughness for which generated BSDF gives more close results to measurement data. It will be considered in the next chapters in more detail.

2.2. Numerical methods of BSDF reconstruction

Nowadays with increasing of computer’s power new approaches for BSDF reconstruction have been developed in [22, 23, 25, 26]. Part of them is based on pure numerical methods in which a BSDF is calculated by ray tracing simulation through an explicit geometry model of rough surface. The method based on the normal density distribution of rough surfaces is proposed in [26, 27]. In this method the micro-relief is simulated with the help of distribution of normals represented with an analytical function having a set of parameters defined with the help of the optimization process. The process of BSDF calculation is presented in Fig. 6.

 

Figure 6. “Normals” method of BSDF reconstruction

 

The approach is maximally natural and transparent. To calculate BSDF the flat boundary presenting a rough surface is illuminated with parallel light from both sides of the boundary. Typically, stochastic (Monte Carlo) ray tracing is used. Each time ray hits the boundary, normal is defined with a probability according to analytical function – normal density distribution. Then ray reflection, refraction is defined according to Snell’s law. The transformed light is registered with detectors which finally form resultant BSDF. The main problem here is the definition of analytical function specifying normal density distribution and its parameters. In [27] it was proposed to use two analytical functions like Gauss and Cauchy (Fig. 7):

	(4)
	(5)

where Θ is an angular variable specifying angle of surface normal, is zero angle specifying the position of function maximum and it should be equal to zero for most of cases for rough surfaces with roughness distribution close to normal. So, these functions are used in our work. The two main parameters δ and n specify shape of function of normal density distribution and can be defined with the help of the optimization process, which is presented in the scheme in Fig. 7 and consists of several main steps:

1. The first step includes an input of measured BSDF and other sample parameters affecting light propagation like refractive index and thickness.

2. An objective function for optimization and parameters of illumination and observation to be used in light simulation are defined at the second step. The measured BSDF of whole sample can be used directly as objective function, sometimes it is recalculated to ordinary angular intensity distribution for simplification. The detector parameters (angular, spatial resolution, distance to measured surface) during simulation are chosen maximally close to parameters of real detectors used in measurements. As a rule only small illumination angles close to normal of measured sample are used during optimization. It is done because the accuracy of measurements decreases significantly for incident angles far from normal direction.

3. During the third step an explicit model of the sample with rough surface is generated for normal density distribution for some initial parameters δ and n.

4. Most of modern light simulation software can simulate light propagation through boundary between two dielectric media specified with normal density distribution. So there is no problem to calculate angular light distribution for sample model defined in the previous optimization step.

5. An optimization criterion is defined as root mean square deviation (RMSD) between measured and simulated angular intensity distributions.

6. An optimization criterion (RMSD) calculated on the previous step together with current δ, n parameters transfer to optimizer. An external optimizer of SCIPY library with “Simplex” algorithm was used in our work.

7. The optimizer makes decision to continue optimization process (7.1 in Fig. 7) or to interrupt it (7.2 in Fig. 7) in case the optimization goal is achieved or due to another reason (for example, maximal number of optimization steps is achieved). If the goal is achieved a final model of rough surface based on optimized normal density distribution is generated. In case of simulation there is no problem to place detectors and light sources anywhere including inside of sample material and calculate light scattering from both sides of rough surface, i.e. to calculate BSDF of rough surface.

 

Figure 7. Optimization procedure for “Normals” method of BSDF reconstruction

 

The optimization procedure (Fig. 7) was used to reconstruct BSDF with “Normal” numerical method. More details are presented in [27]. The investigations show the “Normal” method is very effective from viewpoint of calculation speed and fast convergence in optimization procedure during BSDF reconstruction. However, it has evident drawbacks too, namely, it does not support interreflections and masking-shadowing.

Another numerical approach is based on height density distribution and described in [28]. There is some similarity of the ”Heights” and the “Normals” methods. However, an analytical function is used here for another goal: to define 2D height distribution H(x, y).

(6)

 

Figure 8. Definition of “Height” distribution

 

It shows regular grid of points with uniform steps along x and y axes (Fig. 8). Each point in the grid presents a node of microprofile. To define profile height in each node with (x, y) coordinates analytical probability function of one or several parameters can be used. In other words, height in each node is defined according to a probability defined according to normal (Gauss) or some another analytical function specifying height density distribution. In our work two analytical functions were used for “Heights” approach – Gauss and Cauchy:

	(7)
	(8)

Note that formulas (7), (8) are similar to (4), (5) used for “Normal” approach but use z coordinate instead of angular variable, which is defined in the range [0, H_max] and specifies heights distribution.

Both functions depend on four parameters ( σ, H_max , n and z₀). It is rather substantial number of parameters which can complicate process of optimization convergence. However experiments show that in most of the cases the only σ (sigma) is sufficient, “n” (degree) can slightly improve convergence in some cases. H_max can be set to 1 in most of the cases if to set step between nodes of profile grid around the same unit value. z₀ is supposed to be zero (density of heights is symmetrical relatively to H_max). z is in range [0, H_max]. So defines distribution of height density.

According to formula (6-8) height distribution of microprofile can be defined and used for profile geometry generation (Fig. 9).

 

Figure 9. A schematic appearance of microprofile based on analytical heights distribution: (a) – perspective view; (b) – top view

 

The definition of optimal parameters specifying height distribution is fulfilled with optimization procedure similar to “Normals” methods (Fig. 7). The only difference is an explicit model of sample where rough surface is simulated with a geometry based on heights distribution instead of simplified normal density distribution. At present, most of the modern optical software, such as SPEOS, LightTools, Lumicept [35], for example, allow calculating such microgeometry, so BSDF can be calculated without a problem.

Note that the main advantage of the “Height” method is support of all effects: interreflections and shadowing-masking. The process of BSDF reconstruction with this method is more complex relatively to “Normals” because the number of parameters used in reconstruction is increased. In the “Heights” method the H_max parameter defining the maximal scale of micro profile is added to parameters specifying the shape of the height distribution function, namely parameters δ and n (formulas (4) and (5)).

2.3. Methods of BSDF reconstruction based on wave optics

Most of the methods described in the previous chapters use ray optics for simulation of light propagation. However, an application of geometrical optics can be inaccurate. A rough surface is considered as a combination of microelements and their size varies from great to small values, up to sizes comparable with wavelength. Application of geometrical optics theory can result in the noticeable inaccuracy of reconstructed BSDF. Another problem of the geometrical approach is the parasitic influence of measurement noise in case of measured height distribution.

The main problem of wave methods is their extreme complexity. The precise wave methods cannot be applied practically due to the complexity of micro-surface geometry. So an approximate wave solution should be used for BSDF reconstruction. As an example, one of them is described in [16] but it is related to the reflection model only. The more well-known and more usable method to reconstruct BSDF for rough surfaces is based on the Kirchhoff approximation. The method is built on a simple FFT (Fast Fourier transform) based procedure. A more detailed description can be found in [29-34]. The BSDF reconstruction based on the Kirchhoff approximation is developed for both reflection and transmission components and was examined in our work.

The Kirchhoff method should be applied to the surfaces containing smooth roughness (without breaks) or consisting of enough great facets. The local condition of applicability looks like this:

(9)

where λ is a wavelength (in the medium where scattered light propagates), R is a "typical" curvature radius of roughness, is a local angle of incidence.

The wave-based approach does not consider the multiple reflections. The limitation can be expressed in form:

(10)

where is a characteristic roughness length, is RMS of height deviation from a flat surface.

The method also does not take into account shadowing and masking (shadowing is for occlusion of illumination direction, masking is for occlusion of observation of one). This limitation can be expressed as:

(11)

where is an illumination/observation angle that is counted from a normal to a flat surface. In the wave model scattering is calculated for the infinite periodic surface. If there is no seamless conjugation between opposite sample edges, an artifact scattering by periodic conjugation can arise. It is negligible for a large relief sample but can be quite serious for a small one. The calculations are done for non-polarized illumination.

3. Verification of the BSDF reconstruction methods

3.1. Set of samples for verification

Before describing the samples to be used in the investigation let’s consider a profile of rough surface, Fig. 10.

 

Figure 10. Parameters of microroughness

 

Several wide-used parameters describe a profile rough surface. These parameters will be used for the description of measured samples of rough surfaces, so shortly consider them. The first parameter Ra is the most common one and is calculated using the formula (12):

(12)

is an arithmetical mean deviation of the assessed profile. The next parameter is a root mean square

(13)

The next two parameters specify average values of valleys (heights below the mean line) and ledges (heights above the mean line) over the assessed profile, Rv and Rp correspondently:

;

(14)

The next parameter is the most trivial. Rz is maximal profile height and is calculated using parameters from (14):

(15)

One more well-known parameter is RzJIS or Rz5. It is related to the Japanese industrial format. It is based on the five highest peaks and lowest valleys over the entire sampling length (l in Fig. 10).

(16)

And the last two advanced parameters present Rsk – skewness and Rku – Kurtosis:

;

(17)

The eight samples made of acryl with refractive index = 1.49 are selected for investigation. One surface in each sample has roughness and another is smooth. Two types of measurements are fulfilled for all samples:

1. A height distribution was measured with the precise Taylor Hobson’s profilometer.

2. Light transmission distribution was measured with goniophotometer GP-200 [35, 36]  by Murakami Color Research Laboratory, see Fig. 11.

 

Figure 11. Input measured data of investigated samples

 

The parameters of all eight profiles have been calculated based on measured height distributions and using formulas (12)-(17). These parameters are combined in Table 1. The #1-#8 in the first line are sample identifiers. Additionally, the second and third rows of the  table 1 present size of the measured fragment on the sample and resolution of measurements – number of measured profile points along with both x and y directions. The step between measurement points was constant.

 

Table 1 . Profile parameters of measured samples

Param./Samp..	#1	#2	#3	#4	#5	#6	#7	#8
Size (mm x mm)	0.37x0.37	0.95x0.95	0.95x0.95	0.95x0.95	0.37x0.37	0.95x0.95	0.95x0.95	0.95x0.95
Resolution	1024x1024	1333x1333	1330x1330	1332x1332	1024x1024	1332x1332	1`3x1333	667x677
Ra (μm)	0.178	0.456	0.668	0.738	1.170	2.038	2.596	10.724
Rq (μ m)	0.232	0.581	0.866	0.956	1.466	2.669	3.308	13.456
Rv	1.754	2.908	5.870	5.400	5.726	13.731	16.889	33.913
Rp	0.594	1.613	3.380	2.369	3.588	7.817	7.847	40.7834
Rz	2.340	4.521	9.251	7.769	9.314	21.548	24.736	74.697
Rz5	2.329	4.519	9.235	7.762	9.313	21.542	24.727	74.582
Rsk	-0.786	-0.700	-0.839	-0.928	-0.560	-0.643	-0.658	0.012
Rku	5.156	3.654	4.482	4.492	3.275	4.279	3.920	2.813

 

The images of investigated profiles are presented in Fig. 12. For convenience the profiles in Fig. 12 are placed in order of their root mean square (Rq) increasing from the left to the right and from up to down.

 

Figure 12. The appearance of measured profiles

 

GP-200 goniophotometer [35, 36] was selected for measurements because it has very advanced characteristics, such as angular resolution = 0.6°, very small angular step = 0.1° and wide range of observation directions = ± 90°. The high angular resolution is very important in case of our investigation because part of measured samples has very small roughness, comparable with wavelength, so angular transmission is supposed to have a very narrow shape. The measurements of transmission are done for five angles of incident light direction = 0°, 15°, 30°, 45° and 60° in the single plane of light incidence. The goniometer GP-200 outputs data in the relative shape calibrated to measurements without sample. So, for the correlation of measured and simulated data, the same calibration process is fulfilled in simulation, as it is explained in [31].

3.2. Set of methods to be verified

The first two methods selected for verification are based on the measured heights distribution and fulfilled in the Lumicept lighting simulation software [35]. The software has special instruments for direct simulation of rough microgeometry on the base of numerical height distribution, apart from it has physically accurate Monte-Carlo ray tracing and BSDF generator allowing to calculate BSDF based on ray as well as wave optics (Kirchhoff approximation).

The first method is based on ray optics. It is denoted as “Measured_Heights (ray)”. In the given method explicit geometry is created as the boundary between air medium with refractive index = 1 and dielectric medium with refractive index = 1.49 (the refractive index of sample material). The boundary is illuminated under different incident directions from 0° to 85° with parallel light from both sides: from the air and dielectric. Ray propagation through the rough surface is based on Fresnel, Snell laws. The detectors are placed above and below the boundary of the rough surface and detect transmitted and reflected light. Then BSDF is generated based on calculated data. It supports all complex effects such as interreflections on microrelief, masking and shadowing explained in the second chapter. The main restriction of the method is an applicability of ray optics. It can be inaccurate for the sample with small roughness (with sizes close to the wavelength). Another possible drawback of the method is also related to ray optics. It is the high sensitivity of generated BSDF to the quality of measurements. The different steps between measured nodes or noise can result in a noticeable difference in BSDF shape.

The second method is denoted as “Measured_Heights (wave)”. It uses measured sample profile, i.e. height distribution, too. However, light propagation is realized here analytically based on Kirchhoff approximation. The disadvantages of the approach are listed in chapter 2.3 above.

It should be pointed out that measured profiles are not used directly in this investigation. To minimize possible errors related to the quality of height distribution measurements, application of ray or wave optics, and other possible reasons optimization procedure is run for each profile. It is explained in [26] and similar to optimization procedure presented for “Normals” approach in Fig. 7. The purpose of the optimization is to obtain the transmitted light distribution maximally close to the measured one. And parameters of optimization are scaling and filtration of microrelief. The scaling allows to increase/reduce microroughness and filtration allows reducing measurement noise. These ways of profile modification are presented in Figure 13.

 

Figure 13: Microrelief modification

 

The verification of the next three methods is the main goal of this work. They do not require measurements of height distribution which can be expensive or simply not available. The method denoted as “GGX” was selected as the best representative of analytical approach. A utility was created to generate BSDF based on the analytical formulas (2) and (3). To define an optimal parameter of roughness ( according to formulas (2) and (3)) an optimization procedure, similar to the one presented in Fig. 7 was executed. The optimization goal was to obtain light transmission distribution maximally close to measurements.

The generation of BSDF with two numerical approaches denoted as “Normals” and “Heights” is very similar to “GGX” and is explained in chapter 2.1 and in [31, 32] in more details. The parameters to reconstruct normal and heights distribution are defined with an optimization procedure with an objective function to obtain maximal closeness to measured transmission.

So, finally, we have five methods to be verified: two based on measured profile and angular sample “Measured_ Heights (ray)”, “Measured_ Heights (wave)” and three methods based only on angular sample transmission “GGX”, “Normals” and “Heights”.

3.3. Visual verification of methods

Comparison of measured versus simulated light transmission through a plate with rough surface is used for verification of different BSDF reconstruction methods. However, such comparison can be not sufficient. The BSDF of rough surface can have complex shape and even small inaccuracy in its generation can result in defects visible in the image, appearance of some artifacts. Especially it can be noticeable if BSDF is attached to complex curved objects which are illuminated under grazing angles. So, it is also preferable to verify how BSDF samples are visualized under some realistic conditions.

A special model aimed at visualization was prepared, see Fig. 14. The scene presents a virtual model of a special measuring box JUDGE-II by X-Rite [39]. It has surfaces close to diffuse and several luminescent tube lamps emulating daylight. The several objects: a plate, a sphere and a torus are placed into the measuring box. The reconstructed BSDF is attached to the external surface of the test objects. Internal surfaces are simulated as ideally smooth and have perfect Fresnel properties. The medium of all objects has the refractive index = 1.49, which corresponds to measured samples.

 

Figure 14: Scheme for visualization

 

The scene is observed at a finite distance with a special sensor emulating the human eye or camera. The image is generated with the help of simple forward Monte Carlo ray tracing technique in Lumicept [35]. Although it is not the most effective tool nowadays from viewpoint of efficiency and calculation speed, and generated images, as a rule, contain noise, however it is a more reliable and safe tool because of its simplicity.

4. Results

The results of the simulation are presented in two variants:

1. As graphs with angular distribution of transmitted light intensity. A special scene to simulate the characteristic as precisely as possible has been prepared, which is maximally close to the measurement scheme of GP-200 goniophotometer [31]. The simulation was done for normal incident direction of parallel light in one plane corresponding to the plane of light incidence (“sigma” = 0deg). All six graphs (one measured with GP-200 + all five reconstructed methods) are combined into the single graph picture.

2. As images generated as it is specified in section 3.3. The images are generated with the help of simple forward Monte Carlo ray tracing renderer in Lumicept simulation system [35]. The simulation is fulfilled for all five methods of BSDF reconstruction explained in section 3.2.

 

Figures 15 and 16 present graphs of angular intensity distribution of transmitted light for normal incident direction (sigma = 0deg). More results for other incident directions are published in [40].

 

Figure 15: Angular intensity distribution of transmitted light for samples #1, #2, #3, and #4.

 

Figure 16: Angular intensity distribution of transmitted light for samples #5, #6, #7, and #8

 

Fig. 17 presents images generated for samples #1-#4. Each row in the figure presents the same simulated sample and rows presents different approached of BSDF reconstruction. Fig. 18 presents samples #5-#8.

 

Figure 17: Images for samples #1-#4.

 

Figure 18: Images for samples #5-#8.

 

The general numerical difference (error) between measured and simulated angular intensity distribution of transmitted light is estimated as root mean square deviation (RMSD) reduced to the maximal value of measured intensity in relative shape (*100%):

(18)

where 𝐼_𝑚 is the measured intensity and 𝐼_𝑠 is calculated. Index i means the value of the intensity defined for a specific direction of illumination and observation. All observation directions (in the range of ±90deg with step = 0.1deg) and all illumination directions (sigma = 0, 15, 30, 45 and 60deg) are used in the calculation of the difference.

The value of “error” for all samples and all BSDF reconstruction methods is combined in Table 2. The best result (the lowest error) is highlighted (bolded) for each sample.

 

Table 2. Error for different BSDF reconstruction methods

Sample name	Measured_Heights (ray)	Measured_heights (wave)	Analytical (“GGX”)	Numerical (“Normals”)	Numerical (“Heights”)
Sample #1 (Rq = 0.23μ m)	1.85%	1.83%	3.07%	2.94%	3.15%
Sample #2 (Rq = 0.58μm)	1.45%	1.64%	1.91%	1.81%	1.90%
Sample #3 (Rq = 0.87μm)	1.29%	1.33%	2.02%	1.57%	2.37%
Sample #4 (Rq = 0.96μm)	0.85%	1.07%	2.07%	1.79%	1.96%
Sample #5 (Rq = 1.47μm)	1.51%	3.75%	2.06%	1.69%	2.04%
Sample #6 (Rq = 2.67μm)	2.80%	4.60%	3.03%	1.27%	2.68%
Sample #7 (Rq = 3.31μm)	2.93%	4.14%	3.19%	1.68%	2.75%
Sample #8 (Rq = 13.46μm)	5.56%	27.22%	3.86%	5.23%	3.82%

5. Conclusions

As we can see from the results in Table 2 the most of the investigated methods work well. The exception is “Measured_Heights (wave)” based on Kirchhoff approximation, where we see a noticeable error for samples #5-#8 with Rq > 1μm. This is also clearly seen on graphs (Fig. 16). It is quite predictable analyzing the restrictions of the Kirchhoff approximation based method listed in chapter 2.3. So, the “Measured_Heights (wave)” method cannot be recommended for samples with roughness Rq > 1μm. From the other side, analyzing graphs for samples #1-#3 (Fig. 15), the wave optics approach gives more close results in the shape of angular distribution of transmitted light. Moreover the wave approach almost does not require optimization of measured height distribution unlike the “Measured_Heights (ray)” based on ray optics. It can be explained with the big sensitivity of the ray approach to the quality of microrelief measurements (noise, step between measured nodes).

The analytical “GGX” method (the improved Cook-Torrance model) works reasonably for all samples. In the case of “GGX” the noticeable inaccuracy appears only for samples with big roughness. So, considering its simplicity because only one parameter manages BSDF shape and analytical type of calculation, the method can be recommended for modeling rough surfaces with average microroughness.

Comparison of methods based on measured height distribution (“Measured_Heights (ray)” and “Measured_Heights (wave)”) for samples with small microroughness versus all other methods demonstrates that agreement between measured and simulated intensity is better for methods which use measured geometry of microroughness, especially for big illumination angles. Both numerical methods show good agreement with measurements practically for all examined samples. The numerical “Normals” method is slightly better than numerical “Heights” in the area of general transmission estimation, has better convergence during optimization, and is simpler in calculation. Surprisingly, “Measured_Heights (ray)” works not well for sample #8 (Fig. 16, Table 2). One of the possible reasons is too small measurement area of height distribution, so it is just not representative.

In general the methods based on measured height distribution are supposed to be more precise because the real profile geometry is used during ray transformation. In the case of the “Measured_Heights (ray)” method the interreflection, shading, and masking effects are supported in the whole volume because it uses the Monte Carlo ray tracing. This method can suffer from the restrictions of ray optics, inaccuracy of measurements of height distribution or measured fragments are not representative. However, these drawbacks are overcome with modification of microrelief during optimization procedure, at least partially. Thus the “Measured_Heights (ray)” method can be considered as reference (“etalon”) one in visual comparison with other methods, so all other methods are compared to it.

During visual comparison (Fig. 17 and 18) we see the images of the flat plane with average illumination and observation inclination are similar to each other. The situation with the curved object is more complicated. The images generated with the analytical “GGX” method are similar to the reference (etalon) image in the case of small and average microroughness. The effect of “a dark ring in the sphere” is absent for any images created with the “GGX” BSDF, however the curved objects look darker for samples with big roughness. Likely the energy conversation works not so well in approximations of analytical methods for samples with noticeable microroughness. The numerical “Normals” approach generates quite good images for samples #1-#4 (Fig. 17) but images for rougher samples #5-8 (Fig. 18) have noticeable artifacts like bright edge ring in the sphere. The reason for the effect is evident: the approach does not support interreflection, masking, and shadowing. From the viewpoint of visual appearance the numerical “Height” method shows the best results (most close to the etalon images) for all samples.

Summarizing all simulated data we can recommend the numerical “Heights” distribution method as more accurate in case the precise simulation is required and there are no measurements of microrelief geometry. In case of rather small roughness the analytical “GGX” or the numerical “Normals” can be sufficient.

6. References

[1] K. E. Torrance, E. M Sparrow, Theory for off-specular reflection from roughened surfaces. Journal of Optical Society of America 57, 9 (1967), 1105–1114.

[2] R. L. Cook, K. E. Torrance, A reflectance model for computer graphics. ACM Transactions on Graphics 1, 1 (Jan. 1982), 7–24.

[3] C. Kelemen, L. Szirmay-Kalos, A microfacet based coupled specular-matte BRDF model with importance sampling. Eurographics Short Presentations (2001)

[4] B. Van Ginneken, M. Stavridi, J. J. KoenDerink, Diffuse and Specular Reflectance from Rough Surfaces. Applied Optics 37 (Jan. 1998), 130–139.

[5] S. C. Pont, J. J. KoenDerink, Bidirectional reflectance distribution function of specular surfaces with hemispherical pits. Journal of the Optical Society of America A 19 (Dec. 2002), 2456–2466.

[6] G. J. W. Larson, Measuring and modeling anisotropic reflection. In Computer Graphics (Proceedings of SIGGRAPH 92), pp. 265–272, (July 1992)

[7] B. Walter, Notes on the Ward BRDF. Technical Report, PCG-05-06, Cornell Program of Computer Graphics, (Apr. 2005)

[8] J. Lawrence, S. Rusinkiewicz, R. Ramamoorthi, Efficient BRDF importance sampling using a factored representation. ACM Transactions on Graphics 23, 3, 496–505, (Aug. 2004)

[9] C. Schlick, An inexpensive BRDF model for physically-based rendering. Computer Graphics Forum 13, 3, pp. 233–246, (1994).

[10] M. Ashikhmin, P.S. Shirley, An anisotropic phong BRDF model. Journal of Graphics Tools 5, 2, pp. 25–32, (2000)

[11] M. Ashikhmin, S. Premoze, P.S. Shirley, A microfacet-based BRDF generator. In Proceedings of ACM SIGGRAPH 2000, pp. 65–74, (July 2000).

[12] J.Stamp, An illumination model for a skin layer bounded by rough surfaces. In Rendering Techniques 2001: 12th Eurographics Workshop on Rendering (June 2001), pp. 39–52.

[13] M. I. Sancer, Shadow Corrected Electromagnetic Scattering from Randomly Rough Surfaces. IEEE Trans. on Antennas and Propagation 17, 577–585, (1969).

[14] B. G. Smith, Geometrical shadowing of a random rough surface. IEEE Trans. on Antennas and Propagation, pp. 668–671, (1967).

[15] C. Bourlier, G. Berging, J. Saillard, One- and two-dimensional shadowing functions for any height and slope stationary uncorrelated surface in the monostatic and bistatic configurations. IEEE Trans. on Antennas and Propagation 50, pp. 312–324., (Mar. 2002)

[16] H.D. He, K. E. Torrance, F.X. Sillon, D. P. Greenberg, A comprehensive physical model for light reflection. In Computer Graphics (Proceedings of SIGGRAPH 91), pp. 175–186, (July 1991).

[17] J. E. Rogers, D. K. Edwards, Bidirectional reflectance and transmittance of a scattering-absorbing medium with a rough surface. In Thermophysics Conference (May 1975).

[18] T. A. Germer, Polarized light diffusely scattered under smooth and rough interfaces. In Polarization Science and Remote Sensing, vol. 5158 of Proceedings of the SPIE, pp. 193–204., (Dec. 2003)

[19] J. A. Sanchez-Gil, M. Nieto-Vesperinas, Light scattering from random rough dielectric surfaces. Journal of the Optical Society of America A 8 (Aug. 1991), 1270–1286

[20] M. Nieto-Vesperinas, J. A. Sanchez-Gil, A. J. Sant, J. C. Dainty, Light transmission from a randomly rough dielectric diffuser: theoretical and experimental results. Optics Letters 15, pp. 1261–1263, (Nov. 1990)

[21] B. Walter, S. R. Marschner, H. Li, K. E. Torrance, Microfacet Models for Refraction through Rough Surfaces. Rendering techniques, p.18 ^th , (2007).

[22] Ghosh, A., “Cook-Torrance Model,” Computer Vision, https://doi.org/10.1007/978-0-387-31439-6_531 (2016)

[23] N. Seylan, S. Ergun, A. Öztürk “BRDF Reconstruction Using Compressive Sensing”// 21st International Conference on Computer Graphics, Visualization and Computer Vision 2013. – pp. 88-94. ISBN 978-80-86943-74-9 (2013).

[24] J. B. Nielsen, H. W. Jensen, R/ Ramamoorthi “On Optimal, Minimal BRDF Sampling for Reflectance Acquisition”// ACM TOG 34(6). pp. 1-11 (2015)

[25] J. Filip, M. Havlí, R. Vávra “Adaptive highlights stencils for modeling of multi-axial BRDF anisotropy”// The Visual Computer, Volume 33 (2017).

[26] Sokolov, V.G., Zhdanov D.D., Potemin, I.S., Garbul, A.A., Voloboy, A.G., Galaktionov, V.A., Kirilov N., Reconstruction of scattering properties of rough air-dielectric boundary, Optical Review. 23(5), 834-841, DOI: 10.1007/s10043-016-0250-6, (2016)

[27] Bogdanov N., Zhdanov A.D., Zhdanov D.D., Potyomin I.S., et al, Bidirectional Scattering Function Reconstruction Method Based on Optimization of the Distribution of Microrelief Normals // Light & Engineering, Volume 27, ¹1 (2019)

[28] V. Sokolov, D. Zhdanov, I. Potemin, A. Zhdanov, N. Deryabin. A Bidirectional Scattering Function Reconstruction Method Based on Optimization of Microrelief Heights Distribution, CEUR Workshop Proceedings, 2020, volume 2744, doi: 10.51130/graphicon-2020-2-3-6

[29] E. I. Thorsos, The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian roughness spectrum, J. Acoust. Soc. Am. 83, pp. 78–92, (1988)

[30] P. Beckmann, A. Spizzichino, The Scattering of Electromagnetic Waves from Rough Surfaces, Artech, Norwood, Mass., (1987)

[31] A. Ishimaru and J. S., Chen, Numerical simulation of the second-order Kirchhoff approximation from very rough surfaces and a study of backscattering enhancement, J. Acoust. Soc. Am. 88, 1846–1850 (1990)

[32] J. Richard, Wombell and John A. DeSanto, "Reconstruction of rough-surface profiles with the Kirchhoff approximation," J. Opt. Soc. Am. A 8 , pp. 1892-1897 (1991)

[33] J.Z. Buchwald & C.-P. Yeang, "Kirchhoff's theory for optical diffraction, its predecessor and subsequent development: the resilience of an inconsistent theory", Archive for History of Exact Sciences, vol. 70, no. 5, pp. 463–511; doi:10.1007/s00407-016-0176-1 (2016)

[34] E. Hecht  (2017). "10.4 Kirchhoff's Scalar Diffraction Theory". Optics (5th (Global) ed.). Pearson. pp. 532–535. ISBN 978-1-292-09693-3 (2017)

[35] Lumicept – Hybrid Light Simulation Software, URL: http://www.integra.jp/en

[36] Muracami Color Research Laboratory, https://www.mcrl.co.jp/english/products/p_color_sp/detail/GP200.html .

[37] V.Sokolov, I. Potemin, D.Zhdanov,Boris Barladian, Simulation of the BSDF measurements for scattering materials with GP-200 goniophotometer for light guiding plates, Proc. SPIE 11783, Modeling Aspects in Optical Metrology VIII, 1178305, https://doi.org/10.1117/12.2592378 (2021)

[38] D. Antensteiner,  S. Stolc “Full BRDF Reconstruction Using CNNs from Partial Photometric Stereo-Light Field Data” // Workshop on Light Fields for Computer Vision at ECCV-2017, – pp. 13-21 (2017)

[39] X-Rite, Judge II, URL https://www.bhphotovideo.com/c/product/597881-REG/X_Rite_JUS75U30A_Judge_II_S.html

[40] V. Sokolov, I. Potemin, A. Voloboy. Comparison of BSDF reconstruction methods for rough surfaces. Proceedings of the 31st International Conference on Computer Graphics and Vision (GraphiCon 2022), Ryazan, Russia, September 19-22, 2022 [In print].