1 Introduction

Materials like non-diffuse metallic paints, varnish coatings, and effect paints have complex optical properties that produce fascinating appearance in manufactured products. Metallic paints contain metal flakes, causing the incident light to be specularly reflected. Effect paints are made using thin metal oxide layers on transparent mica platelets [28] and contain pearlescent pigments. The multi-layered structure of pearlescent pigments helps increasing changes in visual appearance of a material, with respect to the incident and viewing directions [28], including angle-dependent spectral reflectance [22]. These materials, often referred as “goniochromatic” [32], are commonly used in the print and packaging industry.

Fig. 1
figure 1

Materials used in our study, rendered inside a Cornell box using the Löw et al. BRDF model [26], under a constant (flat) illumination spectrum. The materials have been assigned, in turn, to different shapes, in order to better display the quality of our results. In a, a blue-green goniochromatic material (measured with our setup) has been assigned to the “Bunny”, the “Blue Metallic Paint” (from the MERL dataset [30]) to the “Happy Buddha” and a gold sample (measured) to the “Asian Dragon” (all models are from the Stanford 3D Scanning Repository). b and c are obtained from a, by rotating anticlockwise the assignment order of the materials to the shapes

Full size image

Such materials are produced using different printing techniques (e.g. offset, gravure, screen printing [22]) and contribute to some of the main challenges in the printing line, including:

  1. 1.

    performing fast and easy process control measurements;

  2. 2.

    synthetically reproduce and match visual properties of packaging materials for customer approval and quality control;

  3. 3.

    communicate material appearance across production and quality control departments in a print production line.

To characterise and communicate the visual appearance properties of goniochromatic print and packaging materials during print production, typically bidirectional reflectance distribution function (BRDF) measurements are required [10, 32]. Commercially available devices, such as multi-angle spectrophotometers and goniospectrophotometers, could be used to perform such bidirectional measurements [22]. However, such devices prove to be slow and relatively expensive, therefore not suited for inline measurements in the production process.

Image-based measurements [14] could represent an efficient, fast, and a practical method to accurately estimate bidirectional reflectance, also in the context of the print and packaging industries, able to satisfy the needs of control and inline quality evaluation. Furthermore, back-scattering measurements often provide enough information to analytically characterise the reflectance properties of a given material [1, 16, 17], whereas additional measurements can be used to improve the initial estimates [9].

Building upon the above, in this paper we estimate the BRDF of print samples using a small set of salient measurements taken with a simple, low-cost image-based setup, easy to integrate in inline quality evaluation for print and packaging industries. These measurements are used to estimate the optimal set of parameters of a commonly used BRDF model [26], by means of a genetic algorithm (GA)-based BRDF fitting method. We demonstrate that even for print samples showing goniochromatic optical properties, typically challenging to capture, we are able to obtain visually satisfactory renderings.

The main contributions of this paper are:

  1. 1.

    the use of common BRDF models to faithfully represent the appearance of non-diffuse, goniochromatic print samples described in this paper;

  2. 2.

    the use of retro-reflective in-plane measurements as a key to successfully represent the appearance of non-diffuse packaging print samples;

  3. 3.

    a GA BRDF fitting method, rather than the commonly used Nelder–Mead down-hill simplex algorithm, to obtain the optimal set of BRDF model parameters.

Along with the print samples, we use an additional non-diffuse sample, “Blue Metallic Paint” (BMP), from the MERL dataset [30]. The BMP material is used to assess the performance of in-plane BRDF measurements to faithfully represent the appearance of materials with visual characteristics comparable to the print samples measured, as opposite to a full measurement dataset. Figure 1 demonstrates the results achievable with our setup and approach, showing some renderings of the materials used in our study, fitted using our GA method. Finally, we compare measurements taken using our setup with the ones of a commercially available goniospectrophotometer.

2 Background and related work

The reflectance properties of a homogeneous, opaque material can be described using the BRDF, defined by Nicodemus et al. [37] as

$$\begin{aligned} f_\mathrm{r}\left( \mathbf {l},\mathbf {v}\right) =\frac{{\mathrm{d}L}_\mathrm{r}\left( \mathbf {v} \right) }{{\mathrm{d}E}_\mathrm{i}\left( \mathbf {l}\right) }= \frac{{\mathrm{d}L}_\mathrm{r}\left( \mathbf {v} \right) }{{L}_\mathrm{i}\left( \mathbf {l}\right) \cos {\theta _\mathrm{i}}d\omega _\mathrm{i}}. \end{aligned}$$
(1)

In Eq. (1), \(\mathbf {l}\) and \(\mathbf {v}\) are incident and viewing direction unit vectors, \(E_\mathrm{i}\) is incident spectral irradiance, \(L_\mathrm{i}\) is incident spectral radiance (flux per unit area, per unit solid angle (\(\omega _\mathrm{i}\))), \(L_\mathrm{r}\) is the reflected spectral radiance and d is the differential. The unit of a BRDF is inverse steradian [1/sr]. There exist a variety of possible designs for BRDF measurement setups [14]; measured data can be left in tabular form, or represented in a compact way by means of BRDF models, either phenomenological or physically based, depending upon the specific needs of the application. In the remainder of this section, we discuss related work on BRDF measurement and models. For a comprehensive survey on the taxonomy of BRDF measurement setups and reflectance models, we refer the reader to the work by Guarnera et al. [14].

2.1 Image-based BRDF measurement setups

A number of image-based measurement setups, making use of one or more cameras as sensors, have been proposed [27, 29, 31, 35, 38, 47, 51], since they allow to perform bidirectional measurements in a fast and relatively inexpensive way. Ngan et al. [35] presented an image-based measurements to measure anisotropic materials like velvet by wrapping it around a cylinder in different orientations. A conceptually similar setup was used in [44] to investigate the suitability of image-based measurements to estimate the reflectance of isotropic packaging materials, represented using two well-known reflectance models (Cook–Torrance [7] and Ward [50]). Three flexible packaging materials, with different optical properties ranging from fairly diffuse to goniochromatic, were measured using both their setup and a goniospectrophotometer; both their setup and the goniospectrophotometer cannot measure retro-reflectance. For BRDF fitting, the Nelder–Mead down-hill simplex algorithm [34] was used, driven by a RMS-based error cost function. Their results show a large relative error, in particular, for the goniochromatic print sample, due to the optimisation algorithm often converging to local minima. In [44], both BRDF models used have 3 free parameters (enforced in the Cook–Torrance model by replacing the Fresnel term with a constant). Phenomenological and physically based models with a higher, yet reasonable amount of free parameters, could provide better generalisation properties for packaging materials.

An important line of research is related to defining the number of measurements needed, that is, finding an optimal and minimal number of acquisitions to represent a material [38]. Xu et al. [51] demonstrated that two image-based measurements can suffice to estimate a BRDF, under the assumption that the reflectance lies in the subspace spanned by the MERL dataset. Their measurement setup used a single near-field fixed camera with multiple lighting directions enabling a simple and fast acquisition method. The need for reducing the number of acquisitions is not limited to BRDF acquisition, with more complex reflectance representations, such as BSSRDFs, presenting additional challenges [13].

2.2 BRDF models

The Lafortune model [24] is a generalisation to multiple steerable lobes of cosine lobe-based models, such as Phong [40]. The generalisation is achieved using a \(3\times 3\) matrix, in which the direction vectors are defined to a fixed local coordinate system with respect to the surface normal. For a single specular lobe, the Lafortune model can be written as:

$$\begin{aligned} f_\mathrm{r}\left( \mathbf {l},\mathbf {v}\right) = \frac{\rho _\mathrm{d}}{\pi } + \rho _\mathrm{s}\left[ C_{x}l_ {x}v_{x} + C_{y}l_{y}v_{y} + C_{z}l_{z}v_z\right] ^{\alpha }. \end{aligned}$$
(2)

In the above, \(\rho _\mathrm{d}\) and \(\rho _\mathrm{s}\) are, respectively, the diffuse and specular albedo, while \(C_x\), \(C_y\), \(C_z\), and \(\alpha \) controls the shape and orientation of the specular lobe, retro-reflection (with \(C_x\), \(C_y\), \(C_z\) as positive), and anisotropy (with \(C_x\ne C_y\)). \(l_{x,y,z}\) and \(v_{x,y,z}\) are direction components of the incident (\(\mathbf {l}\)) and viewing (\(\mathbf {v}\)) direction vectors. Given the six free parameters per lobe, the Lafortune model is potentially more versatile than the Ward model, thanks also to the possibility of emulating the Fresnel effect by using an additional lobe with increasing intensity towards grazing angles. Therefore, it might be expressive enough to fit complex reflectances, while still being efficient and following both reciprocity and energy conversation principles of a BRDF.

Bagher et al. [2] introduced the Shifted Gamma micro-facet distribution with the Cook–Torrance model. Such a distribution results in a more accurate reflectance representation than the Beckmann distribution. Löw et al. [26] introduced two isotropic models for accurate and efficient rendering of glossy surfaces, either based on the Rayleigh-Rice light scattering theory or on the micro-facet theory; both models makes use of a modified version of the ABC model [5, 6]. In particular, the micro-facet model introduced in [26] Eq. (3) is based on the Cook–Torrance model [7]:

$$\begin{aligned} f_\mathrm{r}(\mathbf {l},\mathbf {v})=\frac{k_\mathrm{d}}{\pi }+\frac{S( \sqrt{1-(\mathbf {n}\cdot \mathbf {h}})F(\theta _\mathrm{h})G(\mathbf {n}\cdot \mathbf {l},\mathbf {n}\cdot \mathbf {v})}{(\mathbf {n}\cdot \mathbf {l})(\mathbf {n}\cdot \mathbf {v})}. \end{aligned}$$
(3)

In Eq. (3), \(k_\mathrm{d}\) controls the diffuse component, G and F are, respectively, the geometrical attenuation and Fresnel factors as defined in [7] and given in Eqs. (5) and (6). \(\theta _\mathrm{h}\) is the half angle between the normal and the halfway vector \(\mathbf {h}\), \(\mathbf {l}\) and \(\mathbf {v}\) are the incident and viewing direction vectors, \(\mathbf {n}\) is a normal at a point on the surface. Finally, S is the ABC-based micro-facet distribution, reported in Eq. (4):

$$\begin{aligned} S\left( f\right) =\frac{A}{\left( 1+Bf^2\right) ^C}. \end{aligned}$$
(4)

In the above equation, B and C, respectively, control the width of the specular peaks and the fall-off rate of wide-angle scattering, while A is a scaling factor for the specular component. Therefore, Eq. (4) represents a non-normalised distribution. The term f is defined as \(\sqrt{f_{x}^{2}+f_{y}^{2}}\) where \(f_{x}=\left( \sin \theta _\mathrm{r} \cos \phi _\mathrm{r}-\sin \theta _\mathrm{i} \right) /\lambda \) and \(f_{y}=\left( \sin \theta _\mathrm{r} \sin \phi _\mathrm{r} \right) /\lambda \). \(\lambda \) is wavelength of the incident light.

$$\begin{aligned}&G=\min \left\{ 1,\frac{2\left( \mathbf {n}\cdot \mathbf {h} \right) \left( \mathbf {n}\cdot \mathbf {v} \right) }{\left( \mathbf {v}\cdot \mathbf {h}\right) },\frac{2\left( \mathbf {n}\cdot \mathbf {h} \right) \left( \mathbf {n}\cdot \mathbf {l} \right) }{\left( \mathbf {v}\cdot \mathbf {h}\right) } \right\} \end{aligned}$$
(5)
$$\begin{aligned}&F=\frac{\left( g-c \right) ^{2}}{2\left( g+c \right) ^{2}}\left\{ 1+\frac{\left[ c\left( g+c \right) -1 \right] ^{2}}{\left[ c\left( g-c \right) +1 \right] ^{2}} \right\} . \end{aligned}$$
(6)

In Eq. (6), \(c=\mathbf {v}\cdot \mathbf {h}\), \(g=\eta ^{2}+c{^2}-1\) and \(\eta \) is the index of refraction. In the following, we will refer to the model described by Eq. (3) as “ABC model”.

While the models used in our work are all analytical BRDF models, it is worth to mention that a large number of surface reflectance representations fall in the data-driven class [14], that aims at representing measured reflectance data in a suitable function space, for instance, using factored representations [48]. Soler et al. [45] presented a method for learning a nonlinear manifold of measured BRDFs, starting from a set of reflectance measurements. The measurements are mapped into a 2d latent space, in which novel points can be computed by interpolation and mapped back to the 4d BRDF measurement space.

2.3 BRDF fitting metrics

Estimating the optimal set of the parameters for a reflectance model, given an optimisation algorithm, cost function, and the measured reflectance data, is a common task required to extrapolate material BRDF data, for instance, to be used in rendering. A wide range of different fitting metrics has been used in the previous work. Lafortune et al. [24] minimise the mean-square error of the reflectance multiplied by the cosine of both the incident and outgoing direction. The cost function defined by Löw et al. [26] makes use of a logarithmic function, with the understanding that it yields a better visual reproduction of wide-angle scattering compared to the previously used metrics [35]. Fores et al. [12] used psychometric experiments to demonstrate that a cube root-based fitting metric is perceptually more uniform compared to a RMS error-based metric, and does not depend on the analytical model used. In recent years, perceptually motivated metrics have been further explored [41], and they proved to be useful also in gamut mapping tasks [46]. Recently, Lagunas et al. [25] used a deep learning architecture with a novel loss function to learn a feature space that is well correlated with visual appearance similarity of different materials. Guarnera et al. [15] proposed the use of a perceptually based image similarity metric, which accounts for both colour differences and gradient distribution. However, their approach requires renderings of the input material in a specific setting. In general, the choice of the cost function for fitting is not obvious, and depends also on the sample to be measured and the reflectance model used [26].

Fig. 2
figure 2

The two flexible samples used in our paper, wrapped around a cylinder, used in our measurement setup

Full size image
Fig. 3
figure 3

Colour shift obtained from the spectral radiance factor measurements of the blue-green sample surface when measured at different viewing directions [44]

Full size image

3 Method

3.1 Measurement samples

In our paper, we focus our attention on the two flexible packaging samples reported in Fig. 2. The gold sample is a metallic gold thin cardboard commonly used for decorative purposes in print and packaging industry, while the blue-green sample is a packaging paper printed using effect pigments and varnish coatings. Both samples are non-diffuse with the blue-green sample also being goniochromatic. Figure 3 shows the spectral shift in the blue-green sample with the change in viewing direction. A Munsell white N9/ sheet (MW), produced according to the ANSI standards, was measured along with the gold and the blue-green samples and used as reference white for bidirectional reflectance calculations of samples measured with our setup. Along with the print samples, a metallic paint sample (“Blue Metallic Paint”—BMP), from the MERL dataset [31], is used to assess the performance of in-plane BRDF measurements in representing the appearance of print samples, as opposite to full BRDF measurements.

3.2 BRDF measurement setup

Our measurement setup for flexible samples is schematically represented in Fig. 4. In order to measure a sample, this is wrapped around a cylinder of known radius. Each point on the curved sample surface corresponds to an incident (\(\theta _\mathrm{i}\)) and reflection (\(\theta _\mathrm{r}\)) angle with respect to the surface normal and incident direction (\(\theta _\mathrm{L}\)) of the light source in the setup. Our setup performs in-plane measurements (azimuthal angles \(\phi _\mathrm{i} = \phi _\mathrm{r}=0\)) and the captured image records the radiance (\(L_\mathrm{r}\left( \theta _\mathrm{i},\theta _\mathrm{r}\right) )\) exiting from the sample surface, expressed in terms of digital pixel values (R, G, B); the (R, G, B) values clearly depend also on the per-channel camera sensor spectral sensitivity (\(\bar{r}\), \(\bar{g}\), \(\bar{b}\)), other than on the material properties. Figure 4 (“Top View”) shows a schematic diagram of our measurement setup. A 16 bit Nikon D200 DSLR camera was used for the measurements. A film projector, consisting of a halogen tungsten lamp, was used as the light source. The radius of the cylinder used is 56 mm, while the distance between the cylinder and the light source is 1 m, as well as the distance between the cylinder and the camera. Please refer to [42] for details about the estimation of the incident (\(\theta _\mathrm{i}\)) and reflection (\(\theta _\mathrm{r}\)) angles; details about the accuracy of our measurement setup are reported in [43].

Fig. 4
figure 4

Sample measurements: with our setup, we take measurements at four different illumination directions (\(\theta _\mathrm{L}\))

Full size image

The samples were measured at four different illumination directions (\(\theta _\mathrm{L}= 0^{\circ }, -20^{\circ }, -30^{\circ }\), and \(-40^{\circ }\)) (Fig. 4). In order to capture retro-reflected light from the sample surface, \(\theta _\mathrm{L}= 0^{\circ }\) incident light direction was used during the measurements. Due to practical constraints, it was not possible to have incident light direction (\(\theta _\mathrm{L} = 0^{\circ }\)) in-plane with the camera as it blocks the camera view. In order to overcome this limitation, the samples were illuminated at approximately \(\phi _\mathrm{L}= 4.6^{\circ }\) (see “Side View” in Fig. 4). Since the azimuthal angle (\(\phi _\mathrm{L}\)) is fairly small, we consider these measurements as approximately in-plane, i.e. with \(\phi _\mathrm{i}=\phi _\mathrm{r}=0^{\circ }\).

In order to compare our measurements with the ones from a professional device, the samples were measured using both our setup and a goniospectrophotometer, the Murakami’s GCMS-3B [33] (GCMS in the following). The GCMS records the spectral radiance factor (390–730 nm at 10 nm intervals) at anormal incident (\(\theta _\mathrm{i}\)) and reflection (\(\theta _\mathrm{r}\)) angles in the range of \(+80^{\circ }\) to \(-80^{\circ }\) at \(5^{\circ }\) intervals. GCMS uses a tungsten halogen light bulb as a light source and a silicon photo-diode array as a detector. The sample lays flat on a plate, which rotates between anormal angles \(\pm 80^{\circ }\) with respect to the incident light source, the latter normal to the sample surface; the instrument performs automatic correction for the change in illumination and viewing area due to sample rotation. The reference white plate used in the instrument is assumed to be a perfect reflecting diffuser. Therefore, we calculate its BRDF as \(\beta =\pi f_\mathrm{r}\).

Following the definition of radiance factor [39], the discussions in [20], and using the Munsell White N9 reflectivity (78.66%), we calculate the bidirectional reflectance for the sample as follows:

$$\begin{aligned} f_\mathrm{r}\left( \theta _\mathrm{i},\phi _\mathrm{i},\theta _\mathrm{r},\phi _\mathrm{r},\lambda \right) \approx 0.79*\frac{\pi L_\mathrm{r}\left( \theta _\mathrm{i},\phi _\mathrm{i},\theta _\mathrm{r},\phi _\mathrm{r},\lambda \right) }{L_\mathrm{r}^\mathrm{PRD}\left( \theta _\mathrm{i},\phi _\mathrm{i},\theta _\mathrm{r},\phi _\mathrm{r},\lambda \right) }. \end{aligned}$$
(7)

In Eq. (7), \(L_\mathrm{r}\) and \(L_\mathrm{r}^\mathrm{PRD}\) are radiance at the sample and the perfect reflecting diffuser (PRD) surface. \(\theta \) and \(\phi \) are the polar and azimuth angles, respectively. Indexes i and r are incident and reflected radiation and \(\lambda \) is the wavelength. PRD not being real, in practice, reference white materials like the spectralon tile that can be traceable to a metrological reference or a transfer standard is commonly used as a PRD [19]. We use the MW, which is wrapped around the cylinder along with the gold and blue-green sample, as a PRD in Eq. (7).

3.3 BRDF fitting

3.3.1 Choice of the fitting metric

Print samples, such as the ones used in [44], show non-diffuse and goniochromatic properties which are challenging to visualise. Due to these complex optical properties, we test two different fitting metrics. The metric \(M_{1}\), commonly used in the previous work, is given in Eq. (8). It uses a \(\cos {\theta _\mathrm{i}}\) assuming a uniform incoming radiance at the sample surface thus giving more weight to error in the specular region [35].

$$\begin{aligned} M_{1} =\mathop {\sum }\nolimits _\mathrm{RGB}\sqrt{\frac{\sum _{P}\left[ \left( f_{\mathrm{r}_\mathrm{m}}\left( \mathbf {l}, \mathbf {v}\right) \cos \theta _\mathrm{i} \right) -\left( f_{\mathrm{r}_\mathrm{e}}\left( \mathbf {l}, \mathbf {v} \right) \cos \theta _\mathrm{i} \right) \right] ^{2}}{P}}.\nonumber \\ \end{aligned}$$
(8)

In Eq. (8), P is the measurement at each pixel and \(f_{\mathrm{r}_\mathrm{m}}\) and \(f_{\mathrm{r}_\mathrm{e}}\) are the bidirectional reflectance measured and estimated using the reflectance model, respectively. \(\theta _\mathrm{i}\) is the anormal incident angle. The cost function \(M_{2}\) as defined by Löw et al. [26] is able to produce more visually accurate results, and it is reported in Eq. (9).

$$\begin{aligned}&M_{2}=\mathop {\sum }\nolimits _\mathrm{RGB}\nonumber \\&\quad \sqrt{\frac{\sum _{P}\left[ \text {ln}\left( 1+f_{\mathrm{r}_\mathrm{m}}\left( \mathbf {l},\mathbf {v} \right) \cos \theta _\mathrm{i} \right) - \text {ln}\left( 1+f_{\mathrm{r}_\mathrm{e}}\left( \mathbf {l},\mathbf {v} \right) \cos \theta _\mathrm{i} \right) \right] ^{2}}{P}}. \end{aligned}$$
(9)

In Eq. (9), similar to the \(M_{1}\) cost function, \(f_{\mathrm{r}_\mathrm{m}}\) and \(f_{\mathrm{r}_\mathrm{e}}\) are the bidirectional reflectance measured and estimated using the reflectance model, respectively, and \(\theta _\mathrm{i}\) is the anormal incident angle.

The BRDFs of the print samples, as well as the BMP sample (from the MERL dataset), were estimated using an optimal set of BRDF parameters for the two reflectance models described in Sect. 2, Lafortune [24] and micro-facet model by Löw et al. [26]. The print samples were measured both using our setup and the GCMS instrument.

Fig. 5
figure 5

Directional vectors in the coordinate system used in our setup, with respect to the surface normal at point (P) on the curved sample surface

Full size image

3.3.2 Lafortune model

Figure 5 shows the directional vectors of the Lafortune model, in our setup coordinate system. Since both our setup and the GCMS instrument perform in-plane measurements, and the samples used are isotropic (\(C_{x}=C_{y}=C_{xy}\)), Eq. (2) can be rewritten as in Eq. (10), where we report also the normalisation factor used.

$$\begin{aligned}&f_{\mathrm{r}_\mathrm{RGB}}(\mathbf {l},\mathbf {v}) \nonumber \\&\quad =\frac{\rho _{\mathrm{d}_\mathrm{RGB}}}{\pi }\nonumber \\&\qquad + \frac{\alpha +2}{2\pi \left[ \max \left( \mid C_{xy}\mid ,\mid C_{z}\mid \right) \right] ^\alpha }\left[ -C_{xy}\sin \theta _\mathrm{i}\sin \theta _\mathrm{r}\right. \nonumber \\&\qquad \left. +C_{z}\cos \theta _\mathrm{i}\cos \theta _\mathrm{r} \right] ^\alpha \end{aligned}$$
(10)

where diffuse (\(\rho _\mathrm{d}\)) and specular (\(\rho _\mathrm{s}\)) albedo are optimised per channel. The Lafortune model parameters \(\rho _{\mathrm{d}_\mathrm{RGB}}\), \(\rho _{\mathrm{s}_\mathrm{RGB}}\), \(C_{xy}\), \(C_{z}\), and \(\alpha \) were optimised using \(M_{1}\) cost function, Nelder–Mead down-hill simplex algorithm [34] as the optimisation tool, and the measured data from our setup (all \(\theta _\mathrm{L}\) directions); additionally, measurements from the GCMS were used for comparisons.

3.3.3 ABC model

Using individual diffuse (\(\rho _\mathrm{d}\)) and specular (\(\rho _\mathrm{s}\)) component albedo per channel, the micro-facet ABC model from Eq. (3) can be rewritten as given in Eq. (11) to estimate the sample BRDF.

$$\begin{aligned}&f_{\mathrm{r}_\mathrm{RGB}}\left( \mathbf {l}, \mathbf {v} \right) \nonumber \\&\quad = \frac{k_{\mathrm{d}_\mathrm{RGB}}}{\pi } \nonumber \\&\qquad +\frac{S_\mathrm{RGB}\left( \sqrt{1-\left( \mathbf {n}\cdot \mathbf {h} \right) } \right) F\left( \theta _\mathrm{h} \right) G\left( \mathbf {n}\cdot \mathbf {l},\mathbf {n}\cdot \mathbf {v} \right) }{\left( \mathbf {n}\cdot \mathbf {l} \right) \left( \mathbf {n}\cdot \mathbf {v} \right) }. \end{aligned}$$
(11)

\(S_\mathrm{RGB}\) is the modified ABC distribution with parameter A (in \(S_\mathrm{RGB}\)) being used as a scaling parameter per channel for the specular component albedo and \(k_{\mathrm{d}_\mathrm{RGB}}\) is the diffuse component albedo.

To find a salient measurement dataset for analytically estimating material BRDF using the micro-facet ABC model, we performed in total eight optimisations consisting of two cost functions (\(M_{1}\) and \(M_{2}\)), and four different sets of measurements. Three of these sets of measurements represent different subsets of the measurements made using our setup, as detailed in the following:

  1. 1.

    Illumination direction \(\theta _\mathrm{L}= 0^{\circ }\) (which includes retro-reflective measurements);

  2. 2.

    Illumination direction \(\theta _\mathrm{L}= \{0^{\circ }, -40^{\circ }\}\) (as in the previous point, plus one additional direction to further improve the estimates),

  3. 3.

    All illumination directions: \(\theta _\mathrm{L}= \{0^{\circ }, -20^{\circ },-30^{\circ },-40^{\circ }\}\).

In addition to the above, for each sample, we fitted the GCMS measurements. As for the BMP sample, in-plane measurements from the MERL dataset were used, testing both the \(M_1\) and \(M_2\) metrics.

Estimating an optimal set of ABC model parameters using Nelder–Mead down-hill simplex algorithm proved to be difficult, as A in (\(S_\mathrm{RGB}\)) is not normalised. The model parameters, \(k_{\mathrm{d}_\mathrm{RGB}}\), \(A_\mathrm{RGB}\) (in S), B, C, and \(\eta \) (in F), were therefore optimised for the three samples using the \(M_{1}\) and \(M_{2}\) cost functions and the GA method instead, as detailed in the next subsection.

3.3.4 BRDF fitting algorithm

Fitting measured BRDF data to an analytical BRDF model typically implies optimising for the set of parameters that minimises a given cost function. The cost function, defined over the range of possible parameter values of the BRDF model, measures the difference between the acquired reflectance data and its representation using the selected model. Optimisation algorithms, such as the Nelder–Mead down-hill simplex or Powells, may converge to a local minimum or a saddle point, rather than finding global minima, in particular, in the challenging cases involving non-convex objective functions.

To address this issue, an evolutionary algorithm such as a GA-based method could be used, in particular, in situations where a large number of BRDF model parameters need to be optimised. Due to their applicability both in constrained or unconstrained nonlinear systems, GA methods have been successfully used in computer graphics, for instance, to derive new BRDF models [4], to represent measured subsurface scattering data in a compact way [23], to derive and appearance-preserving mapping between the parameter space of any two arbitrary analytical BRDF models [15], for application-specific tone mappings [8] and to estimate unknown illumination spectra in facial appearance acquisition setups [13].

In the context of BRDF fitting, among the potential advantages, there is an increased probability to have in output a set of model parameters derived from a global minimum of the objective function. In fact, GA test a number of different solutions (represented by the population) at any given step of the optimisation (i.e. generation). Thus, by controlling the population size and the number of stall generations (i.e. number of consecutive generations that do not lead to an improved solution), it is possible to converge to more accurate results, while there is no theoretical guarantee to reach a global minimum. Furthermore, GA does not require the user to specify an initial guess of the parameters, a challenging task due to the complex effect of the parameters on material appearance [15, 36], in particular, in the presence of goniochromatic materials and model parameters with no clear bounds.

The parameters range has a significant impact both on the quality of the solution and on the fitting time. This is particularly true for the ABC model, in which the micro-facet distribution is not normalised and the parameters controlling it (A, B, and C) do not have a clear upper bound, as well as the parameter to control the Fresnel term. To address this issue, we rely on the fitting results available in the supplemental material of [26], under the assumption that the reflectances of our materials lie in the subspace spanned by the MERL dataset. Similar assumptions about the gamut of the MERL dataset have been used in the previous work [30, 41, 51]. Indeed, a wider range of parameters could be used. However, to the best of our knowledge, this has been done so far only for the purpose of conducting detailed experiments on surface appearance perception [49].

In our implementation, with a single panmictic population of 100 individuals, parents for the next generation are selected using the stochastic universal sampling algorithm [3], children are given by the weighted arithmetic mean of two parents, where the weight depends on the fitness values of the parents, and small random mutations are obtained by enforcing a direction in the change which is consistent with the last successful generation, with a step length accounting for the boundaries derived from the MERL dataset, as described in the above.

4 Results

Mitsuba 0.6 [21] was used to render the estimated materials BRDF. To display our results, we used the geometry and lighting described in Havran et al. [18]. Figures 7, 9, and 11 show the renderings obtained using the optimised reflectance models, discussed in Sect. 5.

Fig. 6
figure 6

G-channel gold (a, b) and blue-green (c, d) sample measurements using GCMS instrument and estimation using the Lafortune model. In a and c, we used data measured with our setup, using all the available \(\theta _\mathrm{L}\) directions, whereas in b and d reflectance data was acquired with the GCMS instrument measurements. In all cases, the \(M_{1}\) fitting metric was used

Full size image
Fig. 7
figure 7

Comparison of the images rendered using the Lafortune model, from measurements obtained with our setup (a, c) and with the GCMS instrument (b, d). In all cases, the \(M_{1}\) fitting function was used, along with the Nelder–Mead algorithm

Full size image
Fig. 8
figure 8

G-channel BMP sample measurements from the MERL dataset and estimated using Lafortune model (a) and ABC model (b, c) derived from the subset of in-plane MERL data. For the ABC model, we report both \(M_{1}\) and \(M_{2}\) metrics

Full size image
Fig. 9
figure 9

“Blue Metallic Paint” (BMP) sample renderings. In a, b, and c, in-plane data from the MERL dataset are used, while in d the full BRDF data is used. The limitations deriving from using only in-plane measurements a–c are noticeable, when compared to the renderings obtained by fitting the BRDF with a complete set of measurements (i.e. including out-of-plane data). Depending on the model, the lack of out-plane-data might result in visible blur (a), or reduced visual contrast between diffuse areas and specular peaks

Full size image
Fig. 10
figure 10

Gold sample (metric \(M_1\) first column, metric \(M_2\) second column) and blue-green sample (metric \(M_1\) third column, metric \(M_2\) fourth column), fitted using the ABC model. The first row refers to measurements derived from the GCMS instrument; the second, third, and fourth row refer to measurements from our setup, respectively, using all the available \(\theta _\mathrm{L}\) measurements, only \(\theta _\mathrm{L}=0^{\circ }\) and \(\theta _\mathrm{L} = \{0^{\circ }\), \(40^{\circ }\}\). In all cases, the plots refer to the G-channel

Full size image
Fig. 11
figure 11

Gold sample (first and second column) and blue-green sample (third and fourth column), rendered using the micro-facet ABC model. The reflectance data is measured using the GCMS instrument (first row) and using our setup (second, third, and last row). The rather limited resolution of the GMSC measurements (\(5^{\circ }\)) fails to capture the appearance of the material samples, while our setup allows to derive more faithful results

Full size image
Fig. 12
figure 12

Relative error calculated using Eq. (12) between the measured and estimated measurement data using ABC model optimised with different measurement datasets

Full size image

Figure 6 compares the Lafortune model fits parameters obtained from measurements using our setup (all available \(\theta _\mathrm{L}\) directions) and the GCMS instrument. In all cases, the Nelder–Mead down-hill simplex algorithm converged to local minima, thus preventing to reach satisfactory fits. Since the samples are isotropic, it follows that \(C_{x} = C_{y}\), which causes the BRDF, and hence the cost function, to assume the same value when \(\{ C_{xy}=\xi , C_{z}=\chi \}\) and \(\{C_{xy}=\chi , C_{z}=\xi \}\), with \(\{\chi ,\xi \} \in \mathbb {R}\) (see Eq. 10). Therefore, fitting only in-plane measurements to the Lafortune model leads to additional issues and sub-optimal fits. The same consideration applies regardless of the device used to acquire the in-plane measurements (i.e. GCMS, our setup or the in-plane only data extracted from the MERL dataset) and the fitting algorithm. Therefore, in our experiments we did not further explore the use of the Lafortune model. Figure 7 shows renderings obtained using the Lafortune model.

To assess the quality of the fits achievable using in-plane BRDF measurements, rather than the full BRDF, we used the subset of in-plane reflectance data for BMP sample in the MERL dataset; Fig. 8 reports the results of the experiment, while Fig. 9 shows the corresponding renderings, including as a reference the reconstructed appearance of the material using the whole BRDF data.

Figure 10 compares the fits obtained for the micro-facet ABC model using the GA algorithm, relying on measurements from the GCMS (first row) and from our setup (second to last row); Fig. 11 reports the corresponding renderings.

To objectively compare the effectiveness of the measurements used to estimate the BRDF parameters, we used the relative error (Err), computed using Eq. (12), which accounts for the maximum value in the measurements. In Eq. (12), \(f_{\mathrm{r}_\mathrm{m}}\) represents the measurements obtained using our setup, \(f_{\mathrm{r}_\mathrm{e}}\) represents the data estimated using the optimised micro-facet ABC model, and N is the total number of data points (P).

$$\begin{aligned} \mathrm{Err}=\frac{1}{3}\sum _\mathrm{RGB}\frac{1}{N}\frac{\sum _{P}\mid f_{\mathrm{r}_\mathrm{m}}-f_{\mathrm{r}_\mathrm{e}} \mid }{\max (f_{\mathrm{r}_\mathrm{m}})}. \end{aligned}$$
(12)

Figure 12 shows the relative error using different measurement datasets, for both material samples and metrics.

5 Discussion

With reference to Figs. 6 and 7, in-plane measurements were not sufficient to achieve satisfactory fits for the Lafortune model, regardless the measurement device used. In fact, the \(C_{x} = C_{y}\) condition for isotropic materials, along with the use of in-plane only measurements, resulted into optimisation converging to sub-optimal local minima. Even though the estimated BRDF shows a good fit with the measurements (Fig. 6), the renderings obtained fail to display the goniochromatic properties of the sample. Therefore, for some combinations of surface reflectance and analytical BRDF model, this represents a limitation of our measurement setup, since using out-of-plane measurements could result in more robust estimates. Gold sample renderings (Fig. 7a) show a greenish colour cast, which we believe is due to the spectral sensitivity functions (\(\bar{r}\), \(\bar{g}\), and \(\bar{b}\)) of the camera that was used as a detector.

As for the micro-facet ABC model, the parameter A of the micro-facet distribution acts as a scaling factor for the specular term, thus resulting into a non-normalised distribution. This consideration, along with the lack of a clear upper bound for all the model parameters, suggested the use of a GA method to fit the acquired reflectance data, instead of commonly used optimisation tools. For the print samples measured in this paper, the micro-facet ABC model, fitted using our GA method, allows to obtain visually satisfactory renderings that correctly display non-diffuse and goniochromatic properties. However, our setup is able to acquire few measurements at grazing angles, due to design limitations. Therefore, the lack of information about surface reflectance at grazing angle may affect the estimation of the model parameter related to the material refractive index (\(\eta \)), thus affecting the quality of the estimated Fresnel effect. A possible solution would be replacing the cylinder in our setup with an elliptical surface.

The micro-facet ABC model parameters were derived from different sets of measurements with our setup, comparing the results. In-plane measurements that include the retro-reflective slice of the BRDF (\(\theta _\mathrm{L}=0^{\circ }\)), allowed us to obtain visually satisfactory renderings for the measured samples (Fig. 11). The inclusion of additional measurements (\(\theta _\mathrm{L} = -40^{\circ }\)) increases the quality of the renderings, in particular, for the goniochromatic blue-green sample. The limited resolution of the GCMS instrument (at \(5^{\circ }\) intervals) fails to capture the specular and goniochromatic properties. In comparison, our setup provides a sparser set of measurements, locally more dense. In practice, the density of our measurements depends on the cylinder radius on which the sample is wrapped around, the distance between detector and the sample, and the resolution of the camera used as a detector. Performing measurements using our setup is also faster compared to measuring using the GCMS instrument, as expected for an image-based acquisition setup.

With respect to the relative error calculated between the measurements and estimated data (Fig. 12), the performance of both \(M_{1}\) and \(M_{2}\) metrics is numerically rather similar. Fits obtained using the logarithmic cost function (\(M_{2}\)) led to more realistic renderings, able to faithfully convey the goniochromatic properties of the blue-green sample. The difference between BMP sample rendered using the two different fits achieved for the ABC model, using the \(M_{1}\) and \(M_{2}\) metrics, is noticeable (Fig. 9).

6 Conclusion

We characterise the surface reflectance of two print samples displaying complex optical properties by fitting their BRDF to commonly used reflectance models. Goniochromatic and non-diffuse optical properties are rendered using the estimated BRDF.

In-plane retro-reflective measurements taken with our setup, along with the GA method as a BRDF fitting tool, allowed to estimate an optimal set of the reflectance model parameters. Renderings obtained using the micro-facet ABC BRDF model show that using just in-plane retro-reflective measurements is salient enough to render the reflectance properties of the print samples measured. However, more measurements led to more accurate renderings, in particular, for the goniochromatic sample. In-plane measurements obtained from our setup, as well as from the goniospectrophotometer, did not allow to derive satisfactory fits for the Lafortune model, given its analytical definition for isotropic materials. We believe out-of-plane measurements would be required, in order to solve the resulting ambiguities.

Our measurement setup represents a simple and fast BRDF measurement tool and, along with the GA-based BRDF fitting, could be used to build upon the existing methods [11] for acquiring and rendering discrete sparkles for both isotropic and anisotropic packaging materials, along with goniochromatism and specularity.