# Imaging device to measure the reflective and transmissive part of isotropic BSSRDF

Morgane Gerardin, Mathias Paulin, Romain Pacanowski

► To cite this version:

Morgane Gerardin, Mathias Paulin, Romain Pacanowski. Imaging device to measure the reflective and transmissive part of isotropic BSSRDF. Optics Express, 2024, 32 (22), pp.39267-39292. 10.1364/OE.538429. hal-04736931

HAL Id: hal-04736931

https://inria.hal.science/hal-04736931v1

Submitted on 15 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Imaging device to measure the reflective and transmissive part of isotropic BSSRDF

MORGANE GERARDIN,<fnref n="1" />,<fnref n="2" />,* MATHIAS PAULIN,<fnref n="2" /> ID
AND ROMAIN PACANOWSKI<fnref n="1" /> ID

**Abstract:** The bidirectional scattering-surface reflectance distribution function (BSSRDF) is a fundamental scattering quantity that characterizes the appearance of translucent materials. The few existing BSSRDF measurement facilities that ensure traceable results have shown the high complexity of such measurement, while only accounting for light scattered in the reflection hemisphere. We have developed an imaging device for measuring both the reflective and the transmissive part of the spectral BSSRDF. The complete analysis of the uncertainties associated with the facility allows relative measurement errors of around 10%. The measurements performed on two translucent samples and primary comparisons with metrological measurements are provided and discussed.

Published by Optica Publishing Group under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

## 1. Introduction

Appearance of translucent materials is described by the Bidirectional Scattering-Surface Reflectance Distribution Function (BSSRDF) [1], which is defined as the ratio of the spectral radiance $L_{s,\lambda}$ observed, from a specific direction at a given location on the material surface and the incident spectral flux $\Phi_{i,\lambda}$ at location $\mathbf{x}_i$ from the direction $(\theta_i, \phi_i)$. Since observation can be performed either on the front surface of the material or on its back surface, we make a distinction between the reflective and the transmissive part of the BSSRDF.

We refer to BSSRDF ($f_{ssr}$) when scattered light is observed in reflection (upper hemisphere) at the location $\mathbf{x}_r$ on the front surface of the material from the direction $(\theta_r, \phi_r)$:

$$f_{ssr}(\mathbf{x}_i, \mathbf{x}_r, \theta_i, \phi_i, \theta_r, \phi_r, \lambda) = \frac{L_{s,\lambda}(\mathbf{x}_r, \theta_r, \phi_r)}{\Phi_{i,\lambda}(\mathbf{x}_i, \theta_i, \phi_i)}. \quad (1)$$

We define the observation distance $d_r$ from the illumination location as $d_r = |\mathbf{x}_i - \mathbf{x}_r|$. When observed in transmission (lower hemisphere) from the direction $(\theta_t, \phi_t)$, we refer to the BSSTDF ($f_{sst}$, bidirectional scattering-surface transmittance distribution function). In this case, the observation distance $d_t = |\mathbf{x}_{i,\perp} - \mathbf{x}_t|$ is defined as the distance between the observation point $\mathbf{x}_t$ located on the back surface of the material and the orthogonal projection of the illumination point onto the back surface $\mathbf{x}_{i,\perp} = \mathbf{x}_i - h\mathbf{n}$, where h is the material thickness and n the surface normal at $\mathbf{x}_i$.

$$f_{sst}(\mathbf{x}_i, \mathbf{x}_t, \theta_i, \phi_i, \theta_t, \phi_t, \lambda) = \frac{L_{s,\lambda}(\mathbf{x}_t, \theta_t, \phi_t)}{\Phi_{i,\lambda}(\mathbf{x}_i, \theta_i, \phi_i)}. \quad (2)$$

Therefore we define the BSSxDF (Bidirectional Scattering-Surface Distribution Function) as the function that accounts for reflection and transmission. We use such notation in the following to refer to both BSSRDF and BSSTDF. The spatio-angular parameters involved in the description of the BSSxDF are illustrated in Fig. 1.

<footnote n="1">¹Inria, France</footnote>

<footnote n="2">²IRIT, University of Toulouse, CNRS, UT3, CIMI, Toulouse, France</footnote>

<footnote>*gerardin.morgane@gmail.com</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Geometrical variables involved in the BSSDF measurements. The light is incoming on the parallelepipedic sample at point **x**ᵢ with direction (θᵢ, φᵢ). The observation direction (θ_o, φ_o) is defined so the central pixel (cₓ, cᵧ) of the camera sensor images the incident point **x**ᵢ. Each pixel (i, j) of the CCD sensor collects the emerging radiance either in reflection at a point **x**<sub>r</sub> (on the front surface depicted in green and defined by points **B**<sub>1..4</sub>) in the direction (θ<sub>r</sub>, φ<sub>r</sub>), or in transmission at a point **x**<sub>t</sub> (on the back surface depicted in blue) in the direction (θ<sub>t</sub>, φ<sub>t</sub>). θ angles are defined with respect to the sample surface normal (**n** in reflection, **n**⊥ in transmission) and φ are azimuthal angles defined with respect to the sample surface tangent. The observation distance to the illumination point is defined as d<sub>r</sub> = |**x**ᵢ - **x**<sub>r</sub>| in reflection and d<sub>t</sub> = |**x**ᵢ⊥ - **x**<sub>t</sub>| in transmission, with **x**ᵢ⊥ the orthogonal projection of **x**ᵢ onto the back surface.</figcaption></figure>

In this paper, we introduce a measurement device designed for traceable spectral and isotropic BSSDF measurements that accounts for scattering in both reflection and transmission (see Sec. 3). The calibration procedure is fully described in Section 4, along with the related uncertainties in Section 5. In Section 6, we detail the measurement procedure as well as the processing of the acquired data and show the potential of our setup by providing complete spectral BSSRDF and BSSTDF measurements performed on two samples. Discussions about the results and the facility are provided in Section 7.

**2. Previous work**

Capturing the appearance of translucent materials has been of much interest in computer graphics and optics [2–7], either to develop convincing analytical BSSRDF models or to retrieve the intrinsic parameters of the materials (phase function, absorption and scattering coefficient). The
acquisition of BSSxDF is usually performed by illuminating a sample with a highly focused beam while a camera captures the scattering of light inside the material. Such a setup was proposed by Jensen et al. [2] where the light source, the camera and the sample are at fixed positions. Therefore this optical design allows only for one measurement geometry ($\theta_i, \theta_o$) = (0°, 50°). The simplicity of this setup is sufficient to characterize the material translucency under the assumption that the multiple scattering inside the material does not depend on the direction. However, in a more general case, this is too restrictive. Gkioulekas et al. [7] developed a more complete experimental setup to account for possible dependency on the incident and observation directions of homogeneous materials. With a goniometric approach, two motorized rotation stages allow for the variation of both the incident and the observation directions. However, the reachable measurement geometries are limited to the incident plane with azimuth angles $\phi_i$ and $\phi_o$ equal to either 0° or 180°. Since the measurement geometries (i.e., angular configurations of the sample) are scanned one after the other, the acquisition time to capture a dense BSSRDF with a goniometric approach may be very long. Another method that improves the measurement speed, introduced by Mukaigawa and Inoshita [5,6], is to acquire all desired geometries at once. They use a polyhedral mirror system to create virtual point light sources and cameras on the reflection hemisphere. It allows a full sampling of the BSSRDF, and the characterization of heterogeneous materials in a reasonable time by capturing the BSSRDF at different positions on the sample. However, the measurement geometries are discretized due to the limited number of mirrors used in the polyhedral system and cannot be made denser for a finer study on a specific sample. Moreover, this system requires a complex calibration process to correctly retrieve the position and direction of the cameras and light sources. Contrary to our setup, their system is not metrologically qualified and is not able to measure the transmissiive part of the BSSRDF.

Most of the measurement setups found in the literature are only developed to characterize the light scattered in the reflection hemisphere of a material and therefore capture the BSSRDF. Among them, a few allow for measurements of the BSSTDF, with limited geometries: While the setup from Hong et al. [8] only accounts for the normal transmission, the one from Gkioulekas et al. [7] permits to scan different geometries in the incident plane.

As for the spectral dependency of the BSSxDF, it is often neglected by using a white light source [2,5,6]. There is some interest in this regard with the setup from [7]. The measurements can be performed under three different highly monochromatic illuminations using three laser. However, this reduction to three wavelengths is not representative enough of the whole visible spectrum therefore fine variations between wavelengths may be missed.

One of the biggest issue regarding the measurement of BSSxDF is the traceability of the measurements. Indeed, the previously cited acquisition setups are not metrologically calibrated, therefore the results obtained with the different setups cannot be compared to one another or fitted with BSSxDF models unless they are normalized. The lack of physical units induced by such normalization then make all of these measurement protocols unsuitable to be used as a reference.

A primary facility providing a standard measurement procedure has been developed by Santafe et al. [9]. This goniometric design uses of a 6-axes robot arm allowing for any measurement geometry in the reflection hemisphere to be reached. Even the particular case of back-scattering (where light source and detector are at the same position) is made available by constructing the setup on two levels. It also allows for spectral measurements with a high resolution by filtering the light source with a monochromator. However, because the studied sample is fixed on top of the robot arm, the measurement of the BSSTDF is impossible.

Since both the incident flux and the observed radiance are metrologically calibrated, traceable BSSRDF measurements are obtained. Nevertheless, the number of measurements performed is still too limited to be representative of the full subsurface scattering occurring in a material. Although the setup theoretically allows very dense spectral measurements, only a few geometries in reflection in the incident plane at 550 nm have been acquired. One possible reason for this
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Optical design of our setup for BSSxDF measurement. (a) The white laser beam is spectrally filtered by a continuously variable bandpass filter $F_\lambda$, and focused onto the sample's surface. (b) A camera collects the light scattered by the sample in multiple incidence/observation geometries by rotating the three rotation stages $M_0$, $M_1$, and $M_2$. The camera is rotated using $M_0$ to reach various observation directions. $M_1$ and $M_2$ control the rotation of the sample around two different axes depicted by the red dotted lines, allowing the illumination direction to be varied as well. (c) Schematic representation of the device. Both BSSRDF and BSSTDF can be measured by positioning the camera either in front of or behind the sample.</figcaption></figure>

choice would be to reduce the acquisition time that may be very long given the complexity of the design. Moreover, as the aim of this paper [9] was to provide a standard BSSRDF measurement procedure, only synthetic homogeneous samples with controlled particle size were studied. A more complete study would require a larger variety of samples, including natural ones, as well as an increase of the spatial and spectral measurement resolution. Finally, the optical design only allows to analyze the scattering at distance $d \le 6$ mm from the impact point that can be insufficient for the full study of quite translucent materials.

### 3. Acquisition setup

Our measurement device permits spectral measurements of both the BSSRDF and the BSSTDF. In this section we give details about our instrument design. An outline of our setup is shown in Fig. 2.

The instrument is composed of three main parts: the light source, the sample holder, and the detection system.

The sample is illuminated by a focused light beam of known spectral features (cf. Figure 2). We chose a supercontinuum laser (NKT Photonics, model SuperK COMPACT) emitting from 450 to 2400 nm. According to the manufacturer, this light source shows a stability with less than 1% fluctuations after 1 hour. Therefore, we do not account for this parameter as a source of
uncertainty in the following sections. The beam is first focused using a 30 mm achromatic doublet $L_1$ on a continuously variable bandpass filter $F_\lambda$ (Delta Optical Thin Film, model LF102499). This filter is mounted on a motorized translation stage which permits, by moving the spot position on the filter, to select spectral intervals with central wavelengths $\lambda_0$ in the range 400-700 nm and bandwidths $\Delta_\lambda$ ranging from 13 to 19 nm depending on the central wavelengths. The filter is tilted so the reflected part of the laser beam does not enter back in the light source which may damage it. As a consequence and because the filter consists of a stack of interferential thinfilms, the incident beam is dispersed into multiple distinct beams with different spectral signatures. The desired spectral interval is then selected by a diaphragm $D$, and focused onto the sample surface using two achromatic doublets $L_2$ and $L_3$ of focal lengths 75 and 400 mm respectively. This optical system is aligned so that the spot reaching the sample is reduced to less than 1 mm in diameter. For safety purposes we added a beam dump to trap the light that is reflected at the air/spectral filter interface. As the spectral filter may transmit some infrared radiations above 1200 nm, we also add an infrared filter $F_{IR}$ to ensure only visible light is reaching the sample.

In order to vary both the incident direction ($\theta_i$, $\phi_i$) and the observation direction ($\theta_o$, $\phi_o$) the sample is mounted on three motorized rotation stages from Zaber Technologies Inc. The first one (X-RST120AK-E03) denoted as $M_0$ and supporting the camera rail allows for the camera to rotate around the sample, i.e. around the y-axis (cf. Fig. 2), of an angle $m_0$ so that it scans both the reflection and the transmission hemisphere. It only affects the observation direction. The two other ones $M_1$ and $M_2$ (X-RSW60A-E03) also influence the incident direction by rotating the sample's plane of angles $m_1$ and $m_2$ around the vertical y-axis and the horizontal x-axis respectively. Since only three rotations are available out of the four that are necessary to describe the full BSSxDF (cf. Eqs. (1), (2)), we assume the reflective and transmissive properties of the measured samples are isotropic. The sample holder that is fixed on top of the three rotation stages allows for the study of planar samples with various dimensions and thickness (up to 17 cm in width, 15 cm in height and 4 cm in thickness).

We used a Canon EOS 5D Mark IV camera with pixel size of 5.35 µm along with a Sigma Art 35 mm F1.4 DG HSM lens for detection. The camera is fixed on an horizontal rail at a distance of 47 cm from the sample. Such lens and configuration is chosen in order to optimize both the resolution and the depth of field. It will be used at a fixed aperture of $f_\# = 5.6$ resulting in a depth of field of about 15 mm.

# 4. Calibration

To obtain accurate values of the BSSRDF and BSSTDF (Eqs. (1), (2)), our device is fully calibrated. This includes both radiometric calibrations (cf. Sec. 4.1) and geometric calibrations (cf. Sec. 4.2). All symbols used in the following sections are summarized in Table 1.

## 4.1. Radiometric calibration

The radiometric calibration is performed for both the light source (cf. Sec. 4.1.1) and the camera (cf. Sec 4.1.2) in order to respectively estimate, for each spectral interval of interest, reliable values of the incident flux $\Phi_{i,\Delta\lambda}$ and of the per pixel observed radiance $L_{s,\Delta\lambda}$.

### 4.1.1. Light source

We are interested in the radiant flux that is reaching the sample, therefore the light source calibration is performed at the sample surface position using an integrating sphere, at the end of the illumination optical bench.

We used the absolute integrating sphere method [10] to absolutely calibrate the total spectral radiant flux of the light source, a method that is rather fast and easy to put in place [11,12]. We converted the focused beam into a uniform light source using a Newport 819C-IS-5.3 integrating sphere of radius $R$ and coating albedo $\rho = 0.99$ for any wavelengths in the visible spectral range.
<table><caption>Table 1. Mathematical symbols and their physical description used in this article. N.U. corresponds to "No Units".</caption><thead><tr><th>Symbol</th><th>Description</th><th>Units</th></tr></thead><tbody><tr><td>L<sub>s,Δλ</sub></td><td>Radiance per spectral interval</td><td>W.m<sup>-2</sup>.sr<sup>-1</sup></td></tr><tr><td>Φ<sub>i,Δλ</sub></td><td>Incident flux per spectral interval</td><td>W</td></tr><tr><td>E<sub>Δλ</sub></td><td>Irradiance per spectral interval</td><td>W.m<sup>-2</sup></td></tr><tr><td>x<sub>i</sub></td><td>Illumination point on the object</td><td>m</td></tr><tr><td>x<sub>o</sub></td><td>Observation point on the object</td><td>m</td></tr><tr><td>d = |x<sub>i</sub> - x<sub>o</sub>|</td><td>Observation distance</td><td>m</td></tr><tr><td>(θ<sub>i</sub>, φ<sub>i</sub>)</td><td>Light incident direction</td><td>(rad, rad)</td></tr><tr><td>(θ<sub>o</sub>, φ<sub>o</sub>)</td><td>Camera observing direction (i.e., optical axis direction)</td><td>(rad, rad)</td></tr><tr><td>(θ<sub>r</sub>, φ<sub>r</sub>)</td><td>Pixel observation direction</td><td>(rad, rad)</td></tr><tr><td>m</td><td>Camera magnification</td><td>N.U</td></tr><tr><td>f<sub>‡</sub></td><td>Camera aperture number</td><td>N.U</td></tr><tr><td>(i, j)</td><td>pixel coordinates in image plane</td><td>N.U.</td></tr><tr><td>A<sub>p</sub></td><td>Pixel area</td><td>m<sup>2</sup></td></tr><tr><td>g(i, j)</td><td>Flat field factor</td><td>N.U</td></tr><tr><td>V(i, j)</td><td>Pixel value (i.e., Digital Number)</td><td>DN</td></tr><tr><td>f</td><td>Camera sensor response</td><td>(W.s)<sup>-1</sup></td></tr><tr><td>t<sub>e</sub></td><td>Exposure time</td><td>s</td></tr><tr><td>ρ</td><td>Integrating sphere albedo</td><td>N.U</td></tr><tr><td>α</td><td>Fraction of the integrating sphere composed of open ports</td><td>N.U</td></tr><tr><td>R</td><td>Radius of the integrating sphere</td><td>m</td></tr></tbody></table>

The sphere is fixed on the optical bench so that its entrance port is in the focusing plane of the incident beam. The light beam is sufficiently collimated as it enters the integrating sphere (divergence smaller than 1°). We therefore use the sphere in the collimated beam configuration. For each spectral interval of interest ($λ_0$, $Δλ$), the spectral filter is positioned accordingly. We use a spectrometer (Sekonic, model C-7000) to measure the spectral irradiance $E(λ)$ at the exit port of the sphere (see Fig. 3). For reducing the uncertainties, the spectral irradiance is averaged over 5 measurements. Then, the irradiance per spectral interval defined as:

$$E_{\Delta\lambda} = \int_{\lambda_{min}=380\,nm}^{\lambda_{max}=780\,nm} E(\lambda)d\lambda \quad (3)$$

is deduced by numerical integration from 380 nm to 780 nm with a resolution of 1 nm corresponding to the spectrometer spectral resolution.

Thanks to the integrating sphere properties, the incoming light is scattered multiple times so the irradiance is uniform all over the sphere wall. The irradiance $E_{Δλ}$ is then proportional (cf. [13]) to the total radiant flux that is entering the sphere $Φ_{i,Δλ}$ such as:

$$\Phi_{i,\Delta\lambda} = 4\pi R^2 \frac{(1 - \rho(1 - \alpha))}{\rho} E_{\Delta\lambda} \quad (4)$$

where $α$ is the fraction of the sphere area that is occupied by the opened ports (i.e., incoming and exit ports).

**Results.** As shown in Fig. 4, we characterize the absolute power of 13 spectral intervals of interest. The fall off of the absolute power reaching the sample at smaller wavelengths is due
<figure><img src="image_5.png" /><figcaption>Fig. 3. Experimental setup for the per bandwidth light source calibration. The laser beam is filtered while going through the optical system of Fig. 2 so that a specific central wavelengths $\lambda_0$ and associated spectral bandwidth $\Delta_{\lambda}$ is selected. The light enters an integrating sphere which entrance port is in the focus plane of the beam. The spectral irradiance of light scattered by the sphere is collected using a spectrometer $S_p$ (Sekonic, model C-7000).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Light source calibration. (a) The calibration is done for 13 selected central wavelengths $\lambda_0$ with their associated spectral bandwidths $\Delta_{\lambda}$ in (465, 13.8), (480, 14.0), (495, 14.2), (515, 14.3), (530, 14.4), (550, 15.2), (570, 15.9), (590, 16.4), (610, 16.8), (630, 17.1), (655, 17.8), (675, 17.8) and (700, 18.5) nm), chosen for a well covering of the reachable spectral domain. The spectral irradiance of each bandwidth is averaged over 5 acquisitions at the integrating sphere exit port (cf. protocol on Fig. 3). (b) The absolute power reaching the sample is computed using Eq. (4) for each chosen spectral interval. Dots represent the central wavelengths, horizontal bars illustrate the bandwidths. The uncertainty on the absolute power (cf. Sec. 5) are shown by the vertical error bars.</figcaption></figure>

to the laser source properties. Indeed, the laser does not emit far under 450 nm, so its power density is weaker at these wavelengths leading to small flux values. The decrease above 650 nm is caused by a smaller transmission efficiency of the spectral filter at larger wavelengths.

### 4.1.2. Camera radiometric calibration

The pixels $(i, j)$ of the CCD camera sensor only outputs digital values $V(i, j)$ (in count DN). A proper photometric calibration allows to convert such value into the radiance that is actually reaching each pixel [14–16]. For any pixel $(i, j)$ of the sensor, the corresponding radiance can be modeled (cf. [14]) as:

$$L_{s,\Delta_{\lambda}}(i,j) = \frac{4(m+1)^2 f_{\#}^2 f^{-1}(V_{\Delta_{\lambda}}(i,j))}{A_p \pi g(i,j) t_e} \quad (5)$$
Equation (5) involves camera fixed parameters (magnification $m$ and pixel area $A_p$), and controllable acquisition parameters (aperture $f_\#$ and exposure time $t_e$). It also accounts for the flat field $g(i,j)$ and the inverse camera response function $f^{-1}$ that both require to be calibrated.

**Camera noise:** The camera is subject to various noises that affect the pixel values. These should be corrected in order to account for accurate digital values $V(i,j)$.

For Canon cameras, in order to account for both positive and negative noise fluctuations caused by the analog to digital units conversion, a constant offset called bias is added to the pixel values by the manufacturer. This noise is measured by taking multiple images with no light entering the pixels. It is done by taking pictures in the dark (without removing the lens cap) and with the smallest exposure time available on the camera (1/8000 s). The bias value is then computed as the mean value of the pixels.

A second type of noise is the thermal noise, also known as dark current, that depends on the temperature of the sensor. The longer the sensor is exposed to light, the higher is the sensor temperature and the higher is the thermal noise. For each exposure time available on the camera, ten dark frames similar to the previous bias frames are taken and averaged. The mean pixel value corresponds to the sum of both the bias noise and the thermal noise at this exposure time.

Finally, we only consider as accurate the value of the well exposed pixels, corresponding to the pixels that are neither saturated nor underexposed. We empirically define a threshold value of 5% of the maximum DN value the sensor can reach above which a pixel is considered as saturated, and a threshold of 2% of the bias value under which a pixel is underexposed.

**Results.** Since we are using a RGB camera, the different noises has been studied for each channel independently. An identical bias value of 512 DN is obtained for the three channels. However, depending on the channel, this value is determined with a different statistical uncertainty (±2.46 DN, ±5.54 DN, ±7.71 DN, for R, G and B channel respectively). The fluctuations of the thermal noise appears to be very small compared to the bias value (±1 DN). To account for this negligible noise, we defined the extended bias value of 513 DN as a constant offset that has to be subtracted to every pixels of an image to obtain corrected pixel values.

According to the pixel bit-depth of the camera (14 bit), we chose a threshold value of 15563 DN. Above this maximum value, pixels are considered as saturated. There is no correction in this case, the corresponding pixels are then ignored during the image processing (cf. Sec. 6.3).

**Flat field:** Flat field is mainly the result of the camera lens vignetting and is aperture dependent. It is characterized by imaging a uniform light source as shown on Fig. 5. As described by [15–17], we used an integrating sphere lit by a Thorlabs OSL2 halogen light source. To remove possible internal defects of the integrating sphere, the flat field images were captured through a white diffuse screen, as shown in Fig. 5. Ten pictures were acquired at an ISO setting of 100, an aperture of $f_\# = 5.6$ and an exposure time of 1/15 s. This exposure time was chosen to maximize the signal on the sensor without reaching saturation. A mean image is computed by stacking the 10 individual pictures. For each RGB channel of the sensor, the longitudinal profile is fitted with a third degree polynomial, which is extended to the whole flat field by assuming a symmetry around the optical axis. The use of other models to fit the flat field profile or raise the symmetry assumption are discussed in the Supplement 1.

**Results.** By denoting $r(i,j)$ the euclidian distance between the pixel of coordinates $(i,j)$ from the optical center $(c_x, c_y)$, and $r = r(i,j) - c_x$, we obtain the following flat field fittings:

$$g(i,j) = \begin{cases} 6.0e^{-13}r^3 - 3.2e^{-8}r^2 + 1.9e^{-4}r + 0.70, & \text{for Red channel.} \\ 5.6e^{-13}r^3 - 3.1e^{-8}r^2 + 1.8e^{-4}r + 0.71, & \text{for Green channel.} \\ 4.9e^{-13}r^3 - 2.7e^{-8}r^2 + 1.6e^{-4}r + 0.75, & \text{for Blue channel.} \end{cases} \quad (6)$$

Figure 6 shows a comparison between the flat field measurement and its correction fitted by Eq. (6). To better attest the robustness of our method, we also conducted the experiment with exposures of 1/8 s and 1/30 s, leading to very similar flat field profiles.
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Flat field experimental set-up. The sphere is illuminated using a divergent halogen light source $S_{OSL2}$. At the exit port, the light source is scattered by a white diffuser $D_{iff}$ to neglect integrating sphere internal defects, and then collected by the Canon camera C.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Flat field correction. (a) Flat field images are taken at exposure 1/15 s. Each RGB channel are treated independently, averaged over multiple acquisitions and normalized. The lens vignetting results in lower pixel values on the sides of the pictures. (b) The flat field of each channel is fitted using a third degree polynomial (cf. Eq. (6)). (c) Since flat field is a function of the pixels coordinates, the absolute uncertainty associated with the fit results in a 2D map. The highest error are reached on the very sensor's sides. For a smaller region of interest in the pictures (delimited by the white rectangle) that is effectively imaging the sample, $\sigma_{g(i,j)} \le 0.1$.</figcaption></figure>

**Camera spectral response:** The camera spectral response $f$ is the quantification of how the sensor responds to a given illumination, as a function of the wavelength. When the sensor is illuminated at a given bandwidth $\Delta_{\lambda}$, it establishes a correlation between a pixel value $V_{\Delta_{\lambda}}(i,j)$ and the energy it receives that is the product of the flux reaching the pixel $\Phi_{\Delta_{\lambda}}(i,j)$ and the camera exposure time $t_e$.

$$V_{\Delta_{\lambda}}(i,j) = f(\Phi_{\Delta_{\lambda}}(i,j)t_e) \qquad (7)$$

In this work, we made the assumption that all the pixels of the sensor react to light in the same way, so the camera response defined in Eq. (7) is the same for every pixel. The spatial non-uniformity are corrected by the flat field calibrated in the previous paragraph.

The camera spectral response is measured by imaging the SuperK Compact supercontinuum laser calibrated in Sec. 4.1.1. Similarly to Burggraaff et al. [17], we used an integrating sphere to produce a uniform illumination, as shown in Fig. 7. An accurate but too expensive experimental setup would involve a monochromator to filter the light source [18]. Instead we use our optical bench to lit the sphere with the laser spectrally filtered using the continuously variable bandpass
filter. In this case, radiance at the exit port of the sphere is uniform and related to the measured irradiance (Eq. (3)) by the relationship:

$$L_{\Delta\lambda} = \frac{\rho E_{\Delta\lambda}}{\pi} \qquad (8)$$

with ρ the albedo of the sphere coating. The flux reaching a pixel is then related to the irradiance of the light source by the combination of Eqs. (5) and (8). The energy is changed by varying the exposure time of the camera from the lowest available on the camera to the highest before saturation (from 1/8000 s to 4 s).

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Experimental setup to measure the camera response curve. The integrating sphere is placed such that its entrance port is in the focusing plane of the optical system of Fig. 2. It is lit using the filtered incident light beam. The camera C images the exit port of the sphere that acts as a uniform light source.</figcaption></figure>

For each spectral interval (λ₀, Δλ) of interest, and for each exposure, 5 pictures of the sphere exit port are acquired and stacked into an average picture. We used an ISO setting of 100.

In order to neglect flat field effects, the camera response is determined for center pixels, on a window of 50 × 50 pixels. Each channel of the camera are treated independently. The inverse camera response is then obtained by interpolating the camera response for every possible DN value a pixel can reach (DN ∈ [[0, 2¹⁴ − 1 = 16383]] for a 14-bit sensor), and tabulated for further use.

**Results.** Figure 8 shows the response and the inverse response of the camera when illuminated by a spectral interval centered on λ₀ = 530 nm. For plotting the camera response, different levels of energy reaching the sensor are obtained by changing the exposure time.

Similar response are computed for every spectral interval of interest to characterize the complete spectral response of the camera over the visible range.

These results cannot be verified in absolute since to our knowledge there is no other measurement of the EOS 5D Mark IV spectral response. However, the spectral response has the same tendencies as the one obtained on other Canon camera [19,20], which tends to validate our protocol.

The full radiometric calibration of both the incident flux and observed radiance ensures a reliable absolute BSSxDF value as defined in Eqs. (1), (2).
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Camera response calibration results. (a) Camera response at $\lambda_0$ = 530 nm along with the measurement uncertainties (color bars). Scale of the error bars has been extended (×500) for visibility. Each curve represents the sensor response for a given RGB channel. The amount of energy reaching the sensor is increased by increasing the exposure time. (b) Camera inverse response curves at $\lambda_0$ = 530 nm, tabulated for each RGB channel. The energy values are interpolated so that every possible DN value of the pixels is related to an energy level. The inverse fitted response are plotted (dot curves), with their uncertainties (color areas). Scale of the error has been extended (×20) for visibility.</figcaption></figure>

## 4.2. Geometric calibration

The geometric calibration includes the calibration of the camera intrinsic and extrinsic parameters (cf. Sec. 4.2.1), the calibration of the sample’s position (cf. Sec. 4.2.2) and the motor axes calibration (cf. Sec. 4.2.3).

### 4.2.1. Camera geometric calibration

We need to calibrate both the camera intrinsic parameters (focal length $f'$, optical center ($c_x, c_y$) in pixels, and distortion parameters) as well as the position of the camera relative to the sample surface.

This calibration is performed using a classic 3D vision calibration process using a ChArUco board, that is a combination of Aruco markers [21] and a checkerboard, as calibration target. We used a 11x11 board with chessboard squares of size 12 mm and ArUco markers of size 9 mm. The board is mounted on the sample holder and photographed in multiple different point of view by moving the rotation stages of known motor angles. For later calibration purposes (cf. Sec. 4.2.3), each point of view is reached by moving only one rotation stage. We use the OpenCV [22] library to retrieve the intrinsic and extrinsic camera parameters.

**Results.** We used a distorted pinhole camera model. It is a generalized pinhole model considering that the camera sensor is a surface rather than a plane in order to account for the correction of radial, tangential and thin prism distortion [23]. A 3D point $P$ in world coordinate system is projected onto the pixel $p$ following:

$$p = A[R|t]P \qquad (9)$$

with $A$ the camera intrinsic matrix and $[R|t]$ a rotation-translation matrix.

The matrix $A$ accounts for the parameters that are intrinsic to the camera and therefore do not depend on the scene viewed. With this model, we obtained:
$$A = \begin{bmatrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 6598.42 & 0 & 3438.62 \\ 0 & 6584.61 & 2215.78 \\ 0 & 0 & 1 \end{bmatrix} \qquad (10)$$

Since we are using a circular lens, the x-axis and y-axis focal lengths, respectively $f_x$ and $f_y$ are pretty similar (6598.42 vs. 6584.61 pixels, which are equivalent to 35.30 vs 35.23 mm). As expected, the optical center of the model ($c_x, c_y$) = (3438.62, 2215.78) in pixels is close to the center of the image (3372, 2251). These small differences may be the result of slight mechanical misalignment between the camera body and lens.

For each point of view of the camera, the transformation matrix $[R|t]$ projects a point in world coordinates system towards a point in the camera coordinates system. One such matrix is then obtained for each viewpoint used in the calibration.

Usually, camera lenses show mostly radial distortion and slight tangential distortion. The results from our modeling indeed confirm such hypothesis. By denoting ($x_d, y_d$) the location of the distorted pixels, ($x, y$) the location of the undistorted pixels in the image plane, and $r^2 = x^2 + y^2$, the distortion is corrected using the following model:

$$\begin{cases} x_d = x(1 + k_1 r^2 + k_2 r^4 + k_3 r^6) + 2p_1 xy + p_2(r^2 + 2x^2) \\ y_d = y(1 + k_1 r^2 + k_2 r^4 + k_3 r^6) + p_1(r^2 + 2y^2) + 2p_2 xy \end{cases} \qquad (11)$$

where ($k_1, k_2, k_3$) = (-0.269, 10.0, -162) and ($p_1, p_2$) = (-0.00331, 0.00315) are the parameters describing respectively the radial and the tangential distortion. According to the modeling results, the distortion parameters relative to the thin prism distortion are null therefore this kind of distortion is negligible. Otherwise it would be required to complete the distortion model in Eq. (11).

For the establishment of this camera model, we used 60 positions of the camera around the ChArUco board moving the rotation stages of the optical bench. The re-projection error describing the accuracy of the model in estimating the position of the board across these different points of view is about 0.61 pixels in image space, which corresponds to 44.7 µm in object space.

### 4.2.2. Material sample geometric calibration

The design of the measurement device allows the study of parallelepiped material samples with various dimensions. There is no specific setting to define the sample position on the optical bench (see Fig. 2), so it is adjusted to select a particular point on the sample's surface to be illuminated. Therefore, in general, the sample coordinate system defined by the normal $n$, the tangent and the bi-tangent (see Fig. 1), and whose origin is the barycenter of the front face (i.e., the face receiving the illumination) is not aligned with the one from the geometric calibration of the camera (cf. Sec. 4.2.1). In order to estimate correctly the position of the sample with respect to the camera, we need to compute the rigid transformation between the sample position and the origin of the world coordinate system.

The location of the sample is defined by the position of the four vertices of its front surface $B_{1..4}$. In the ideal case where all coordinate systems are perfectly aligned, these positions are easy to retrieve knowing the origin and the dimensions of the sample. In practice, these are triangulated using the direct linear transformation (DLT) [24] that we adapted to account for 9 different camera viewpoints. Since the repeatability of the rotation stages is pretty high (<0.087 mrad for $M_0$, <0.349 mrad for $M_1$ and $M_2$), we used a subset of the camera positions that were used for the establishment of the camera model. The projection matrices from world space to camera space in these configurations have then already been estimated. For each viewpoint, the pixel coordinates of the vertices $B_i$ are manually selected on the picture. The DLT provides
for each vertex the corresponding 3D position that minimizes the reprojection error across the multiple viewpoints. The displacement of the sample from its ideal location is finally computed as the rigid transformation between the initial positions and the triangulated set of points [25].

### 4.2.3. Motor axes calibration

Since the rotation stages manufacturing implies inaccuracies as well as their fixing on the optical table, the actual motor axes of the rotation stages may be slightly different from the theoretical ones. We calibrate them to ensure a precise sample position prediction for a set of given motor angles.

Our calibration objective for each rotation stage is to establish its rotation axis (its direction and its 3D position) depicted on Fig. 2(b). More precisely, for each rotation stage, we use the pictures of the ChArUco board taken for the establishment of the camera model (cf. Sec. 4.2.1). These pictures are taken for different motor angles of the rotation stage. For a given set of pictures, the movement of a given corner of the ChArUco board must be contained in a plane that can be computed by solving a least-square problem. The normal of that plane corresponds to the rotation axis of the motor stage. To compute the center of the rotation, the 3D positions of the corners are projected onto the plane and their 2D coordinates are fitted with an ellipse [26]. Since they are many corners on a ChArUco board, we obtain a set of axes (resp. a set of centers of rotation) that we average directionally (resp. spatially) to obtain the final rotation axis and its 3D position. Performing such calibration for each rotation stage independently leads to the complete calibration of the rotation axes of the three motors.

**Results.** Each position of the camera being reached by a given combination of the motor angles $m_0, m_1$, and $m_2$, the rotation of the camera relative to its original position (along the sample surface normal, pointing towards the sample) is defined by the matrix $R_c$ such as:

$$R_c = R_{x_2}(m_2)R_{y_1}(m_1)R_{y_0}(m_0) \quad (12)$$

where $R_X$ is the rotation matrix around the axis X, and $x'_k$ and $y'_k$ are the calibrated rotation axes of the rotation stage $M_k$ mainly rotating around the x-axis or y-axis respectively.

Our motor axes calibration process is assessed by comparing at a given geometry ($m_0, m_1, m_2$) the triangulated position of the sample's front face from a photograph and the position of the vertices of a mesh with same dimensions as the object projected into an image using the camera model from Section 4.2.1 rotated according to Eq. (12). Such comparison is illustrated by Fig. 9. This calibration allows to reliably predict the position of any point of the object, for any rotation of the camera, within 1.6 mm (cf. Sec. 5.2).

<figure><img src="image_11.png" /><figcaption>(a) Material sample photograph.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Before calibration.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) After calibration.</figcaption></figure>

**Fig. 9.** Calibration of $m_0$ rotation axis. (a) Picture of the sample fixed on the sample holder, captured by the measurement camera (cf. Fig. 2(b)). The camera is rotated using the motor $m_0$. (b) Estimation of the sample position by projecting a virtual mesh of same dimensions as the sample onto the image plane assuming that the motor axis $m_0$ turns around the y-axis=(0, 1, 0). The photograph shows a clear positioning offset between the virtual mesh on the real-sample. (c) After the motor axis calibration, the camera is supposed to turn around the axis (-0.01868, -0.9994, 0.02756). Consequently, the projected virtual mesh superimposes much better to the real sample.
The different geometric calibrations described above allow to accurately retrieve the spatio- angular parameters of the BSSSDF that are the positions **x**<sub>i</sub> and **x**<sub>o</sub> of the incoming light point and the observation point on the sample respectively, as well as the angles defining the incident light direction and pixel dependent observation directions.

## 5. Uncertainty analysis

<table><caption>Table 2. Uncertainties associated with the parameters involved in the radiometric BSSSDF measurements. The uncertainties from manufacturer are highlighted in cyan.</caption><thead><tr><th>Symbol</th><th>Description</th><th>Absolute value</th><th>Relative value</th></tr></thead><tbody><tr><td>σ<sub>ρ</sub></td><td>Integrating sphere albedo</td><td>0.00265</td><td>2.68 × 10<sup>-3</sup></td></tr><tr><td>σ<sub>D<sub>s</sub></sub></td><td>Integrating sphere diameter</td><td>0.001 cm</td><td>10<sup>-6</sup></td></tr><tr><td>σ<sub>D<sub>i</sub></sub></td><td>Entrance port diameter</td><td>0.004 cm</td><td>2 × 10<sup>-3</sup></td></tr><tr><td>σ<sub>D<sub>o</sub></sub></td><td>Exit port diameter</td><td>0.003 cm</td><td>10<sup>-3</sup></td></tr><tr><td>σ<sub>E<sub>Δλ</sub></sub></td><td>Irradiance per bandwidth</td><td>∈ [0.392, 4.46] × 10<sup>-2</sup> mW.m<sup>-2</sup></td><td>∈ [0.677, 1.68] × 10<sup>-4</sup></td></tr><tr><td>σ<sub>Φ<sub>i,Δλ</sub></sub></td><td>Incident flux</td><td>∈ [0.0044, 0.037] mW</td><td>0.098</td></tr><tr><td>σ<sub>f<sub>‡</sub></sub></td><td>Aperture</td><td>Negligible</td><td>Negligible</td></tr><tr><td>σ<sub>t<sub>exp</sub></sub></td><td>Exposure time</td><td>Negligible</td><td>Negligible</td></tr><tr><td>σ<sub>A<sub>p</sub></sub></td><td>Pixel area</td><td>0.005 µm<sup>2</sup></td><td>2 × 10<sup>-3</sup></td></tr><tr><td>σ<sub>m</sub></td><td>Magnification</td><td>1.8 × 10<sup>-5</sup></td><td>2.5 × 10<sup>-4</sup></td></tr><tr><td>σ<sub>g(i,j)</sub></td><td>Flat field</td><td>≤ 0.1</td><td>≤ 0.14</td></tr><tr><td>σ<sub>f<sup>-1</sup></sub></td><td>Camera inverse response</td><td>R channel: ≤ 7.24 × 10<sup>-11</sup> mJ<br>G channel: ≤ 3.10 × 10<sup>-10</sup> mJ<br>B channel: ≤ 1.48 × 10<sup>-9</sup> mJ</td><td>≤ 8.63 × 10<sup>-3</sup><br>≤ 4.34 × 10<sup>-3</sup><br>≤ 8.93 × 10<sup>-3</sup></td></tr><tr><td>σ<sub>L<sub>s,Δλ</sub>(i,j)</sub></td><td>Observed radiance</td><td>Sample dependent<br>Alabaster: ≤ 548 W.m<sup>-2</sup>.sr<sup>-1</sup><br>Fimo clay: ≤ 548 W.m<sup>-2</sup>.sr<sup>-1</sup></td><td>Sample dependent<br>≤ 0.136<br>≤ 0.136</td></tr></tbody></table>

<table><tr><td></td><td>Uncertainty from manufacturer</td></tr><tr><td></td><td>Computed Uncertainty</td></tr></table>

In order to obtain more reliable results, we establish all uncertainties associated with the BSSSDF measurement. In this section we study both the absolute radiometric uncertainty related to the BSSSDF value (σ<sub>f<sub>xxx</sub></sub>) and the geometric uncertainties related to the angular and spatial parameters of the measurements (σ<sub>θ<sub>i</sub></sub>, σ<sub>φ<sub>i</sub></sub>, σ<sub>θ<sub>r</sub></sub>, σ<sub>φ<sub>r</sub></sub>, σ<sub>x<sub>i</sub></sub>, σ<sub>x<sub>r</sub></sub>). For our uncertainties characterization, we rely on Farrance et al. review of rules for calculating uncertainty [27]. Detailed computations are given in the Supplement 1. The uncertainties of every parameter are summarized in Table 2.

### 5.1. Radiometric uncertainty

By definition of the BSSSDF (cf. Eqs. (1), 2), the relative uncertainty of the BSSSDF value is related to both the uncertainty of the incident flux and of the observed radiance:

$$u^2(f_{ssx}) = u^2(\Phi_{i,\Delta\lambda}) + u^2(L_{s,\Delta\lambda}) \quad (13)$$

where $u(x) = \frac{\sigma_x}{x}$ is the relative uncertainty of variable $x$.

Each parameter involved in Eqs. (4) and (5) has an uncertainty of measurement that propagates to the BSSSDF value [27,28]. Some uncertainties are given by the manufacturers (integrating
sphere characteristics, or camera pixel size for examples). The others are derived from the light source and camera radiometric calibration procedures (cf. Sec. 4.1.1 and 4.1.2).

**Central wavelength:** The uncertainty of the central wavelength is directly related to the uncertainty of the translation stage that moves the spectral filter. According to the manufacturer and our measurements the translation stage positional uncertainty is 2 µm. Such value implies an uncertainty on the position of the central wavelength of 0.015 nm.

**Incident flux:** From Eq. (4), the uncertainty associated with the incident flux depends on multiple uncorrelated physical parameters from the light source and the calibration setup. It therefore depends on the spectral interval and can be expressed as follows:

$$ \sigma^2(\Phi_{i,\Delta\lambda}) = \left(\frac{\partial\Phi_{i,\Delta\lambda}}{\partial D_s}\right)^2 \sigma_R^2 + \left(\frac{\partial\Phi_{i,\Delta\lambda}}{\partial \rho}\right)^2 \sigma_\rho^2 + \left(\frac{\partial\Phi_{i,\Delta\lambda}}{\partial \alpha}\right)^2 \sigma_\alpha^2 + \left(\frac{\partial\Phi_{i,\Delta\lambda}}{\partial E_{\Delta\lambda}}\right)^2 \sigma_{E_{\Delta\lambda}}^2 \quad (14) $$

First of all, parameters from the integrating sphere are involved. According to the integrating sphere specifications, and the coating calibration certificate [29], the coating albedo $\rho$ is known up to $\sigma_\rho = 0.00265$ and the sphere diameter up to $\sigma_R = 10^{-3}$ cm. Through the term $\alpha$, the incident flux also depends on the entrance and exit port diameters that we measured ten times using a caliper leading to respectively $D_i = 2.509 \pm 0.004$ cm and $D_o = 2.526 \pm 0.003$ cm.

As for the uncertainty associated with the irradiance per bandwith $E_{\Delta\lambda}$, it is given by:

$$ \sigma^2(E_{\Delta\lambda}) = \sum_{w=\mu}^{M-1} \left( \frac{\lambda_{w+1} - \lambda_w}{2} \right)^2 [\sigma^2(E(\lambda_w)) + \sigma^2(E(\lambda_{w+1}))] \quad (15) $$

where the sum on $w$ is the result of the integration over the spectral range of the spectrometer [$\lambda_\mu = 380$ nm, $\lambda_M = 780$ nm] using the trapezoidal rule, and $\sigma^2(E(\lambda)) = \frac{1}{5 \times 4} \sum_{m=1}^5 (E_m(\lambda) - E(\lambda))^2$ is the uncertainty of the spectral irradiance averaged over 5 measurements. Each measured spectral irradiance $E_m(\lambda)$ is affected by the spectrometer accuracy ($\sigma(E_m(\lambda)) = \frac{1}{\sqrt{3}}(0.05E_m(\lambda) + 10^{-12})$). Note that according to the spectrometer spatial resolution, $\lambda_{w+1} - \lambda_w = 1$ nm in our case. See Supplement 1 for more details.

Finally a relative uncertainty of $u(\Phi_{i,\Delta\lambda}) = 0.098$ was obtained for every spectral interval of interest in the visible spectral range.

**Observed radiance:** The observed radiance (Eq. (5)) is a function of the intrinsic parameters of the camera and more specifically of the sensor's response to light. Even though the response is assumed to be identical for every pixels, the vignetting induced by the camera lens involves that the uncertainty associated with the radiance is a spectral 2D map accounting for the pixel coordinates.

$$ \sigma^2(L_{s,\Delta\lambda}(i,j)) = \left(\frac{\partial L_{s,\Delta\lambda}}{\partial A_p}\right)^2 \sigma_{A_p}^2 + \left(\frac{\partial L_{s,\Delta\lambda}}{\partial g(i,j)}\right)^2 \sigma_{g(i,j)}^2 + \left(\frac{\partial L_{s,\Delta\lambda}}{\partial m}\right)^2 \sigma_m^2 + \left(\frac{\partial L_{s,\Delta\lambda}}{\partial f^{-1}}\right)^2 \sigma_{f^{-1}}^2 \quad (16) $$

The aperture and the exposure time are discrete variables of the camera from the manufacturer. We considered both uncertainties related to these parameters as negligible.

Since the sensor is composed of squared pixels with known size $s_p = 5.35 \pm 0.005$ µm, the uncertainty related to the area of one pixel is computed to be $\sigma_{A_p} = 0.05$ µm².

Theoretically, the magnification is defined on the optical axis. Therefore we considered a constant value for every pixels of the sensor. It corresponds to the average of the magnification factors computed for every sides of every 9 × 9 mm² ArUco marker detected on an image of the ChArUco board taken at normal observation. Knowing the pixel size and that, from the camera calibration, the markers are detected with an uncertainty of $\sigma_p = 0.61$ pixels, we obtain $m = 0.073 \pm 1.8 \times 10^{-5}$.

Since the flat field accounts for the non uniformity over the camera sensor, the associated uncertainty is computed on each pixel as the root mean square error between the flat field
measurements and the prediction of Eq. (6) at this pixel. This results in a 2D uncertainty map (see Fig. 6(c)) where the maximum values lie on the sides where the vignetting effect has more influence. For a sample that is imaged in a smaller region of interest centered in the pictures, the uncertainty would not exceed an absolute value of $\sigma_{g(i,j)} = 0.1$. Considering the smallest value along a longitudinal profile across the flat field image, it corresponds to a maximum relative uncertainty of $u(g(i,j)) = 0.14$.

The non saturated data from the inverse response are approximated using an adequate analytical model (cf. Fig. 8(b)). If saturation is reached by none of the RGB channels, they all show a linear response to light and so is the inverse response. However, when one of the RGB channels with a nearly linear response reaches saturation, the other channels show a non-linear response (cf. Fig. 8). The inverse response is then fitted using a third degree polynomial. By definition, the camera response and so the inverse response is null when no light is incoming on the sensor, and strictly positive otherwise. The fits are then performed under the constrained $f^{-1}(0) = 0$. The standard deviations of the fitted parameters are given by the diagonal elements of the fit covariance matrix (see Supplement 1 Sec. 1.C for more details). We obtain the absolute inverse camera response uncertainty $\sigma_{f^{-1}}$ by propagating the parameters errors through the fitted model. It depends on the pixel value, as for the actual inverse response. The results of the fits and uncertainties evaluation performed at $\lambda_0 = 530$ nm are given on Fig. 8. For other wavelengths, we refer the reader to the Supplement 1.

The radiance uncertainty is sample dependent and therefore will be computed for every studied sample by evaluating Eq. (16) (cf. Sec. 6).

## 5.2. Geometric uncertainties

In this section we detail the geometric uncertainties of our BSSxDF measurement. The geometric uncertainties regroup spatial (i.e., $\sigma_{x_i x_r}$) and angular ($\sigma_{\theta_r, \phi_i, \theta_r, \phi_r}$) uncertainties.

At the moment, our device assumes that the normal (see vector **n** in Fig. 1) of the material sample as well as the tangent and bi-tangent are constant over its surface.

The spatial positions $\mathbf{x}_i$ and $\mathbf{x}_r$ are computed numerically by using our camera model (cf. Sec. 4.2.1) and the virtual parallelepiped mesh representing our material sample (cf. Section 4.2.2). Indeed, for each image pixel, we find the intersection between the ray, which starts at the camera position ($P$) and crosses the pixel, and the virtual plane defined by the four vertices $\mathbf{B}_{1..4}$ (cf. Fig. 1).

The uncertainty $\sigma_P$ of the camera is related to the calibration procedure of the extrinsic parameters (cf. Sec. 4.2.1) for which the reprojection error in terms of distance is very low (less than a pixel which corresponds to a physical distance of 44, $7 \times 10^{-6}$ meters). Furthermore, we rely on the very high repeatability ($<0.087$ mrad for $M_0$ and $<0.349$ mrad for $M_1$ and $M_2$) of the three rotation stage motors, which allows us to assume that, once the motor axes have been calibrated, the position of the camera and the material sample are always correctly retrieved.

Consequently, the uncertainties $\sigma_{x_i}$ and $\sigma_{x_r}$ of the spatial parameters depend mostly on the calibrated positions $\mathbf{B}_{1..4}$ of the virtual plane. These positions are found thanks to a triangulation procedure (cf. Sec. 4.2.2), which yields the uncertainty results for $\mathbf{B}_{1..4}$ of 1.6 mm (see Section 1.D of Supplement 1 for more details).

Finally, the angular uncertainties $\sigma_{\theta_i}$, $\sigma_{\phi_i}$ (resp. $\sigma_{\theta_r}$, $\sigma_{\phi_r}$) are directly related to the 3D coordinates uncertainties $\sigma_{x_i}$ (resp. $\sigma_{x_r}$) and the position uncertainty of the camera $\sigma_P$, which yields a smaller than 0.05 degree uncertainty for $\sigma_{\theta_i}$ (resp. $\sigma_{\theta_r}$ and smaller than 2 degree for $\sigma_{\phi_i}$ (resp. $\sigma_{\phi_r}$) (see Supplement 1 for a full derivation).
# 6. Measurements results

## 6.1. Material samples

We have performed both BSSRDF and BSSTDF measurements on two different parallelepiped samples with known dimensions and showing different levels of translucency (see Fig. 10):

* A white Fimo clay cut to a thickness of 0.78±0.0029 cm. This sample is very homogeneous.

* An alabaster block with thickness 2.6±0.0029 cm. Although it appears quite heterogeneous, the sample is more translucent so lateral light transport in the material is more important.

<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Top: Photograph of the measured samples under ambient illumination and bottom: under a 655 nm normal illumination $\theta_i = 0°$ during BSSRDF measurement at $(\theta_o, \phi_o) = (20°, 0°)$ (b,d) and BSSTDF measurement at $(\theta_o, \phi_o) = (200°, 180°)$ (c,e). For alabaster (b,c), short exposure times (0.067s in reflection, 0.0167s in transmission) are used. Since Fimo clay (d,e) is less translucent, longer exposure times are required to visualize the lateral transport (0.25s in reflection, 2s in transmission).</figcaption></figure>

## 6.2. Measurement parameters

### 6.2.1. Geometries

We chose a set of 720 measurement geometries $(\theta_i, \phi_i, \theta_o, \phi_o)$ to acquire the BSSDF in both the reflection and the transmission hemisphere. Since our device allows for three degrees of freedom, we chose the incident direction such that $\theta_i \in \{0°, 10°, 20°, 30°, 40°, 50°, 60°\}$ and $\phi_i \in \{0°, 60°, 120°, 180°, 240°, 300°\}$, as well as the observation direction $\theta_o \in \{0°, 10°, 20°, 30°, 40°, 50°, 60°\}$. For each combination of $(\theta_i, \theta_o, \phi_i)$, it exists two possible values for the azimuthal observation angle $\phi_o$. The same geometries are chosen for both the reflective and the transmissive
hemisphere. Among the defined geometries, some are duplicated, some cannot be reached by the system, and other are obstructing the incident light beam. These are removed leading to the final set of measurement geometries.

### 6.2.2. Exposure times and wavelengths selection

For every position of the sample and the camera, the 13 spectral intervals centered in $\lambda_0 \in 465, 480, 495, 515, 530, 550, 570, 590, 610, 630, 655, 675, 700$ nm and calibrated in Sec. 4.1.1 are used. Under each illumination, we acquire multiple pictures at different exposures $t_e$. Such high dynamic range acquisition is necessary to well characterize the lateral transport of light inside the materials. Indeed, the higher energy around the point $\mathbf{x}_i$ implies quick saturation of the pixels in $\mathbf{x}_i$ neighborhood. A contrario, for points at long distance from $\mathbf{x}_i$ low energy values require longer exposure time.

Different exposures are chosen depending on the translucency of the sample:

* For the Fimo clay sample, we use 10 fixed exposure times per measurement geometry. For each exposure, the measurement is averaged over 5 pictures. For the BSSRDF we use 0.001s, 0.002s, 0.004s, 0.008s, 0.016s, 0.033s, 0.067s, 0.125s, 0.25s, 0.5s, as for the BSSTDF we use 0.004s, 0.008s, 0.016s, 0.033s, 0.067s, 0.125s, 0.25s, 0.5s, 1s, 2s.

* For the alabaster sample, we use only three exposure times per measurement geometry. These were chosen on the fly during the experiment to ensure at least 70% of the pixels in a region of interest around incident spot $\mathbf{x}_i$ were well exposed. Only one picture is taken for each exposure. We use exposure times between 0.002s and 0.125s for the BSSRDF, and between 0.008s and 2.0s for the BSSTDF.

Finally, the aperture of the camera lens is fixed at $f_H = 5.6$ for all measurements. This choice allows a depth of field of about 15 mm, while taking advantage of the sensor resolution. Indeed, in this configuration, the resolution is limited by the diffraction characterized by an Airy disk of 6.1 to 9.5 µm depending on the wavelength, which remains quite close to the pixel size of 5.35 µm. Improving the resolution would require to increase the aperture at the expense of the depth of field. This aperture is then a good compromise to increase both depth of field and picture resolution.

### 6.3. Data processing

One complete BSSXDF measurement is performed automatically in about 60 hours. Considering that the camera needs a delay of one second between the acquisition of two pictures for the sensor to cool down, about 12 min (resp. 15 min) are necessary to fully characterize one measurement geometry of the Fimo clay in reflection (resp. in transmission), which corresponds to 650 pictures (10 exposures x 5 acquisitions x 13 spectral intervals). For alabaster, this reduces to 39 s to 1 min 24 s depending on the exposure times that are selected as only 39 pictures (3x1x13) are taken. The full measurement data require from 0.8 TB (alabaster sample) to 2.5 TB (Fimo clay sample) to be stored depending on the number of exposures that are used.

Each acquired picture is processed to retrieve the per pixel BSSXDF value and spatio-angular parameters. An overview of the complete measurement process is summarized in Fig. 11. More details about acquisition and data processing schemes are given in the Supplement 1 (cf. Sec. 2., Fig. S9.).

First, the pictures are corrected to get rid of the camera noise (cf. Sec. 4.1.2). Since the thermal noise ($N_D = \pm 1$ DN) is negligible compared to the bias noise ($N_B = 512$ DN), this correction is performed by subtracting a slightly overestimated bias value that accounts for the thermal noise fluctuations ($N_B = 513$ DN).

We then proceed to the selection of well exposed pixels. Saturated pixels are removed according to an empirical threshold value of 5% of the maximum DN value the sensor can reach, while
<figure><img src="image_15.png" /><figcaption>**Fig. 11.** Overview of the measurement process. A set of parameters describing the measurement configurations and the camera parameters are provided by the user. Our optical device automatically proceeds to the acquisition. For each measurement configuration, the data processing scheme converts acquired pictures into a BSSxDF data file requiring much less storage capacity. More details about acquisition and data processing schemes are given in Figure S9 of Supplement 1.</figcaption></figure>

under exposed pixels are removed according to an empirical threshold value of 2% of the bias value. For the acquisition in the reflection hemisphere, apart from hot pixels of the camera sensor, the saturated pixels are the one imaging the incident spot on the front face of the sample. This allows to estimate the spot position on the images. We chose to ignore these pixels in the processing as they would provide erroneous values, which explains the lack of BSSRDF data near $d_r = 0$ on the results from Figs. 12 and 13.

<figure><img src="image_16.png" /><figcaption>**Fig. 12.** BSSxDF profiles of alabaster sample in the plane of incidence and illuminated with different spectral intervals. BSSxDF data (dots) are plotted against the distance $d_r = x_i - x_r$ between the observation position $x_r$ on the sample and the illumination spot position $x_i$. (a) Results in reflection are obtained for the measurement geometry ($\theta_i, \phi_i, \theta_o, \phi_o$) = (0°, 0°, 20°, 0°), with relative uncertainties lying in the range $u(f_{ssr}) \in [0.1, 0.16]$. (b) Results in transmission are obtained for the measurement geometry ($\theta_i, \phi_i, \theta_o, \phi_o$) = (0°, 0°, 200°, 180°), with relative uncertainties lying in the range $u(f_{sst}) \in [0.1, 0.12]$. The sample is quite heterogeneous resulting in the small variations across the profiles, and shows an important dependency on the wavelength.</figcaption></figure>

Thanks to the geometric calibration (cf. Sec. 4.2), the acquired scene is virtually reproduced using a mesh of same dimensions as the measured sample whose original position is determined through triangulation and rigid transformation. The position of the camera with respect to the sample in any measurement geometry is computed knowing the angles and motor axes of the rotation stages. The spatio-angular parameters are retrieved using ray-tracing from the camera to the first intersection with the virtual mesh. The incident parameters ($\theta_i, \phi_i, \mathbf{x_i}$) are
<figure><img src="image_17.png" /><figcaption>**Fig. 13.** BSSxDF profiles of Fimo clay sample in the plane of incidence and illuminated with different spectral intervals. BSSxDF data (dots) are plotted against the distance $d = x_i - x_r$ between the observation position $x_r$ on the sample and the illumination spot position $x_i$. (a) Results in reflection are obtained for the measurement geometry ($\theta_i, \phi_i, \theta_o, \phi_o$) = (0°, 0°, 20°, 0°), with uncertainties lying in the range $u(f_{ssr}) \in [0.12, 0.21]$. (b) Results in transmission are obtained for the measurement geometry ($\theta_i, \phi_i, \theta_o, \phi_o$) = (0°, 0°, 200°, 180°), with relative uncertainties lying in the range $u(f_{sst}) \in [0.1, 0.11]$. As the sample is very homogeneous, the profiles are very smooth. The dropout on the BSSTDF profile is due to the vertical heterogeneity of the sample on Fig. 10(e). The sample shows a slight spectral dependency in transmission. The BSSTDF measured signal below 550 nm and above 655 nm was too low to be detected with our system.</figcaption></figure>

computed considering a ray passing through the pixel imaging the center of the light spot, while the observation parameters ($\theta_r, \phi_r, x_r$) are computed considering the previously selected well exposed pixels. The BSSxDF values of each pixels of a captured image are estimated using Eq. (5). Only the pixels whose distance to the mesh lies in the focus region of the camera lens are actually converted into radiometric values and stored. The first and last focus plane (respectively at a distance of 46.6 cm and 48.1 cm from the camera) of the camera are computed using the intrinsic parameters from camera calibration, the acquisition aperture, as well as a circle of confusion with a diameter that we defined at 3 pixels to be of same order of magnitude than the diffraction limit. Outside this range, pixels are considered to be out of focus. Their values result from an average over a neighborhood, and are therefore not reliable to characterize accurately the BSSxDF at position $x_r$.

According to this process, all the acquired images are processed independently from one another, and corrected by the exposure time at which they were taken through Eq. (5). For one measurement geometry and one spectral illumination, a pixel of coordinates $(i, j)$ can be defined as valid (i.e., well exposed and in focus) for different exposures. Therefore multiple radiance values are computed for the same observation direction and location on the sample. The measured radiance at this pixel is then defined as the average of the multiple radiance values over the number of exposures for which the pixel is valid, and over the number of measurement repetitions at these exposures. Since pixels are processed independently, no debayering algorithm is applied. For each RGB channel, we use the corresponding calibrated camera response curve (cf. Sec. 4.1.2).

After processing, for each measurement geometry, the per pixel spatio-angular parameters and spectral BSSxDF values are stored in a Python numpy binary table, which is processed for BSSxDF visualization. This reduces the required storage to 24 GB for alabaster sample and 1.6 GB for the clay sample.
## 6.4. Visualization

We demonstrate the capabilities of our system by visualizing some results obtained on the samples described in Sec. 6.1. Results in both reflection and transmission are discussed here.

The BSSRDF and BSSTDF profiles in Figs. 12 and 13 show the dependency on the distance $d_t = |\mathbf{x}_i - \mathbf{x}_r|$ in the incident plane, where $\mathbf{x}_i$ is the position of the illumination spot on the front face in reflection, and its projection on the sample's back face in transmission. Because of the oblique observation and light propagation into the samples, the maximum of the BSSTDF profiles are not obtained at $d = 0$ mm. In order to reduce the measurement noise, the data are spatially averaged over 1 mm wide areas along the profile.

With our system and data processing procedure, we are able to characterize the subsurface scattering at distances up to 3 cm from the spot. While in reflection the profiles are highly peaked around the spot, these are much smoother in transmission. Since the Fimo clay is more opaque, it shows a sharper and narrower profile in reflection. It also shows a much lower BSSTDF than alabaster despite being thinner (maximum BSSTDF value of about $60 \text{ m}^{-2}\text{sr}^{-1}$ for Fimo clay and $10^4 \text{ m}^{-2}\text{sr}^{-1}$ for alabaster).

These results also demonstrates how important the spectral dimension of the measurement is. Although the alabaster sample appears whitish under ambient light, its measurement under different illuminations show great variations. Therefore it is essential to proceed with the measurements at different spectral intervals to fully characterize this material. However in the case of Fimo clay, the spectral dependency is much more negligible. This observation allows to limit the spectral resolution of the measurement and therefore to drastically reduce the acquisition time. This can also be done for the benefit of a greater spatial angular resolution for example.

The smoothness of the profile provides information about the homogeneity of the sample. The heterogeneity of the alabaster sample results in slight variations on the profile values that are systematic whatever the wavelength of interest (at $d = -7$ mm in reflection or at $d = 3$ mm in transmission for example). It also highly affects the symmetry of the transmission profiles. In comparison, the Fimo clay appears very homogeneous except for a dropout at $d = -1$ mm that is the signature of a thickness discontinuity of the sample visible on Fig. 10(e).

For both samples and for every measurement configuration, the BSSSDF relative uncertainty lies between 10% and 20%. Higher uncertainties are obtained in reflection for observation point close to the illumination spot, due to the possible saturation of the sensor. Otherwise, it is mainly dominated by the uncertainty of the incident flux, and more precisely by the contribution of the integrating sphere albedo $\rho$. A first solution for the reducing the BSSSDF uncertainty would require to minimize the uncertainty $\sigma_{\rho}$. Since this parameter is being characterized by metrological institutes, its value is already measured with a high accuracy. Another solution would be to establish a more precise absolute incident flux calibration protocol.

The angular distribution ($\theta_r$, $\Delta\phi = |\phi_i - \phi_r|$) of the measurements in reflection for a given incidence $\theta_i$ and a given distance $d_r$ from the illumination spot $\mathbf{x}_i$ is illustrated by the polar representations on Figure 14. The missing geometries in the incident plane are explained by the fact that our design does not allow for back scattering to be measured. Except for normal incidence that implies an ambiguity on the definition of $\phi_i$, for small incidences, the measurements are quite uniformly distributed across the reflective hemisphere. However, for more grazing incidences the acquired data tend to be denser around specific geometric configurations, leading to larger unconsidered areas on the reflective dome. This distribution may be improved with a better choice of the measurement geometries depending on which incident direction is at stake.
<figure><img src="image_18.png" /><figcaption>**Fig. 14.** Angular distribution ($\theta_r$, $\Delta\phi = |\phi_i - \phi_r|$) of the measurements performed at 550 nm in reflection on alabaster for a given incidence $\theta_i$ and for a given range of distances from the illumination spot $d_r = |\mathbf{x}_i - \mathbf{x}_r| \in [2, 2.1 \text{ cm}]$. Each plotted point represents one BSSRDF measurement. The total number of valid measurements after data processing (cf. Sec. 6.3) is displayed as $N_o$. Because of the ambiguity on $\phi_i$ at normal incidence, only nearly in plane measurements can be performed. Out of plane measurements are obtained while increasing the incidence. For grazier incidences, some domain of the reflective hemisphere remains completely ignored (gray areas).</figcaption></figure>

## 7. Discussion and future work

Our measurement device has been developed for the characterization of the appearance of translucent materials. The filtered white light source allows for a spectral illumination by selecting specific spectral intervals covering the visible spectral range. Our goniometer design offers a high angular resolution, while the detection part of the system offers a better spatial resolution and a wider range of studyable observation positions on the sample than previous measurement device. Thanks to the full radiometric and geometric calibration, we obtain traceable BSSRDF results that are comparable to the results from a metrological institute [9] (cf. Fig. 15). As far as we know, our device is the first one that also allows for the measurement in transmission, offering new insights on the relationship between the translucency and the thickness of materials. The device is fully automated and only requires a set of measurement geometries and a set of illumination central wavelengths to proceed to the full BSSRDF acquisition.

<figure><img src="image_19.png" /><figcaption>**Fig. 15.** Validation of the BSSRDF measurements performed on sample S9 from [9], at 550 nm, for $\theta_i = 45^\circ$, $\phi_i = 0^\circ$, $\theta_o = 0$, $\phi_o = 180^\circ$. The profile obtained with our device (green round dots) are compared with the profiles from [9] obtained with the device from CSIC (black triangle dots). Our device allows to measure the BSSRDF up to ±4 cm away from the incident spot. To the width of the CSIC profiles, our results matches the one from CSIC with a nearly constant offset.</figcaption></figure>
The system still has some limitations regarding the accuracy and the efficiency of the BSSDF measurements, that we envision to improve. First, the BSSRDF accuracy may benefit from a new definition of some parameters involved in the radiometric conversion from pixel value to radiance (cf. Eq. (5)), in particular the camera aperture. It is defined as the focal length $f'$ and external pupil diameter $D$ ratio $f_\# = f'/D$ and is assumed to be perfectly known from the camera manufacturer. This definition is however an approximation since the lens diaphragm is not perfectly circular. Imaging the camera diaphragm using a mirror [30] would allow the measurement of the aperture along with an associated uncertainty. The uncertainty on other parameters may also be minimized. For the camera inverse response it is possible, for a given illumination spectral interval, to only consider the RGB channel with the smallest uncertainty. This would however be at the expense of the spatial resolution: due to the Bayer matrix of the sensor, only one pixel over four (respectively two over four) would be processed in case the R or B channel is selected (respectively G channel). Another solution to reduce the uncertainties would be to modify the calibration protocol. Indeed, we used a third-party spectrometer for both the flux and the camera calibration. Since the BSSDF is defined as a ratio, a relative calibration by imaging the light source directly with the camera will be sufficient, eliminating the use of the spectrometer and therefore its associated uncertainty. However, such calibration requires the use of optical density to reduce the signal and avoid damaging the camera sensor. As for the uncertainties on the spatio-angular parameters, these are limited by the user ability to manually select the sample vertices used in the triangulation process (cf. Sec. 4.2.2). A more reliable procedure would be to automatically detect the sample position using calibrated target, such as ArUco markers, fixed on the sample surface. Finally, the camera lens may involve coatings that introduce spectral distortion. Its correction would improve the spectral accuracy of the measured BSSDF.

Improving the measurement efficiency relies on the shortening of the acquisition time. To do so, several options are envisioned. The more naive one is to consider less measurement geometries or less illumination spectral intervals, but this comes with a loss in spatial or spectral resolution of the acquired BSSDF. A first improvement that does not degrade the resolution of the measurement is to optimize the rotation of the motorized stages in order to reduce the delay that is necessary for the device to go from one measurement geometry configuration to the next one. The choice of exposure times could also be improved, in order to maximize the high dynamic range imaging using a limited number of exposures selected on the fly for each sample. Finally, another solution would be to use a curved homogeneous material sample and to take advantage of the known topography to acquire a larger range of observation geometries in a single image. Since fewer images would be taken, this reduces the need for storage capacity as well.

Our system also suffers from a limited depth of field we are aiming to improve. So far, a fixed camera focus setting, chosen so that the front face of the sample is at focus, is used for the whole acquisition of the BSSDF. The position of the back face of thick samples may therefore lies before the first focusing plane of the camera, resulting in an averaging of the BSSTDF. A more reliable approach would be to define two focus settings depending on whether the reflective or transmissive part of the BSSDF is measured. An actual increasing of the depth of field would require to stack acquisitions performed for a range of focus settings. The resulting image can be reconstruct by selecting the sharpest pixel of each view. This would however multiply the number of images taken and therefore increase the acquisition time.

The full BSSDFs of two samples have already been measured. We are aiming to perform measurements in both reflection and transmission on a wide variety of translucent materials with various thicknesses. All measurements will be made publicly available for further use in the community.

In order to attest the validity of our BSSRDF measurements, we would like to compare our results with the one obtained by other existing traceable setups [31]. We obtained primary
comparisons with the measurements performed with the device from Santafe et al. [9], developped at the Instituto de Óptica 'Daza de Valdés' of the Spanish National Research Council (IO-CSIC). We used the sample denoted as S9 in their work and performed acquisitions for a set of geometries in the reflective hemisphere at 550 nm. A comparison of the in plane BSSRDF profiles is given on Fig. 15. To the width of the IO-CSIC profile, our results match the previous work with an offset that appears to be nearly constant across the profiles (cf. Fig. 15). We are considering the establishment of a finer comparison factor between the two measurements as a future work. This difference may come from the uncertainty on the incident flux calibration procedure or from a constant noise that has been misidentified in one or the other of the devices design. It is noticeable that our design allows to study observation points that are further away from the impact point $x_i$ (up to ±4cm against ±6mm for CSIC).

We are aiming to extend such study, which implies comparing out of plane measurements as well as the spectral dependency of the BSSRDF on various reference samples. Such validation would be a first step towards the establishment of a calibration standard for BSSRDF measurements that would be transposable to any acquisition device. However the validation of BSSTDF remains unsolved, since our instrument is the first to offer a full transmission characterization.

Another future work would be to make use of the measured data for the inverse problem solving. The propagation of light inside translucent materials is a function of intrinsic parameters that remain unknown and difficult to estimate independently from one another from experiment. These parameters are involved in the radiative transfer (RT) theory, therefore the use of differentiable rendering techniques [32,33] should allow to retrieve a solution set of values from the BSSDF measurements. In computer graphics, these intrinsic parameters are used to evaluate the RT equation and render images of corresponding translucent objects with complex geometries, which would be beneficial either for the reproduction of existing objects such as cultural heritage artefacts for example, or for the virtual prototyping of new objects made out of characterized translucent materials.

## 8. Conclusion

We developed a BSSDF measurement device for the spectral characterization of volumetric light scattering in translucent materials for parallelepiped samples, in both reflection and transmission. The full calibration process associated with the facility allows for traceable measurements. By providing the measurements performed on two different translucent samples, we demonstrate the device ability to fully characterize the BSSDF, spectrally, spatially, and angularly. The discussion about the results opens to multiple future works to either reduce the measurement uncertainties or improve the acquisition time of the device.
**Funding.** Institut National de Recherche en Informatique et en Automatique (INRIA); Agence Nationale de la Recherche (ANR-17-CE23-0017, Project VIDA); Agence Nationale de la Recherche (ANR-11-LABX-0040, Labex CIMI).

**Acknowledgments.** The authors would like to thank Pablo Santafé-Gabarda and Alejandro Ferrero from the Instituto de Óptica 'Daza de Valdés' of the Spanish National Research Council (IO-CSIC) for sharing their BSSRDF measurements, allowing for a metrological validation of the results.

**Disclosures.** The authors declare no conflicts of interest.

**Data availability.** Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

**Supplemental document.** See Supplement 1 for supporting content.

**References**

1. F. E. Nicodemus, J. C. Richmond, J. J. Hsia, et al., *Geometrical considerations and nomenclature for reflectance*, (National Institute of Standards and Technology, Gaithersburg, MD, 1977).

2. H. Wann Jensen, S. R. Marschner, M. Levoy, et al., "A practical model for subsurface light transport," in *Seminal Graphics Papers: Pushing the Boundaries*, Volume 2, (2023), pp. 319–326.

3. M. Goesele, H. P. Lensch, J. Lang, et al., "Disco: acquisition of translucent objects," in *ACM SIGGRAPH 2004 Papers*, (2004), pp. 835–844.

4. S. G. Narasimhan, M. Gupta, C. Donner, et al., "Acquiring scattering properties of participating media by dilution," in *ACM SIGGRAPH 2006 Papers*, (2006), pp. 1003–1012.

5. Y. Mukaigawa, S. Tagawa, J. Kim, et al., "Hemispherical confocal imaging using turtleback reflector," in *Computer Vision–ACCV 2010: 10th Asian Conference on Computer Vision*, Queenstown, New Zealand, November 8-12, 2010, Revised Selected Papers, Part I 10. (Springer, 2011), pp. 336–349.

6. C. Inoshita, S. Tagawa, M. A. Mannan, et al., "Full-dimensional sampling and analysis of bssrdf," *Information and Media Technologies* **8**, 1105–1109 (2013).

7. I. Gkioulekas, S. Zhao, K. Bala, et al., "Inverse volume rendering with material dictionaries," *ACM Trans. Graph.* **32**, 1–13 (2013).

8. S. Hong, C. Lee, and S. Chin, "Physically based optical parameter database obtained from real materials for real-time material rendering," *Journal of Visual Languages Computing* **48**, 29–39 (2018).

9. P. Santafé-Gabarda, A. Ferrero, N. Tejedor-Sierra, et al., "Primary facility for traceable measurement of the bssrdf," *Opt. Express* **29**, 34175–34188 (2021).

10. Y. Ohno and Y. Zong, "Nist facility for total spectral radiant flux calibration," in *Proc., Symposium of Metrology*, (2004).

11. Y. Ohno, "New method for realizing a luminous flux scale using an integrating sphere with an external source," *J. Illum. Eng. Soc.* **24**(1), 106–115 (1995).

12. Y. Ohno, "Detector-based luminous-flux calibration using the absolute integrating-sphere method," *Metrologia* **35**(4), 473–478 (1998).

13. R. Boyd, *Radiometry and the Detection of Optical Radiation* (Wiley, 1983).

14. R. Fiete and Society of Photo-optical Instrumentation Engineers, *Modeling the Imaging Chain of Digital Cameras*, SPIE Digital Library (SPIE Press, 2010).

15. M. Pagnutti, R. E. Ryan, G. Cazenavette, et al., "Laying the foundation to use raspberry pi 3 v2 camera module imagery for scientific and engineering purposes," *J. Electron. Imaging* **26**(1), 013014 (2017).

16. O. Burggraaff, N. Schmidt, J. Zamorano, et al., "Standardized spectral and radiometric calibration of consumer cameras," *Opt. Express* **27**(14), 19075–19101 (2019).

17. C. A. Coburn, A. M. Smith, G. S. Logie, et al., "Radiometric and spectral comparison of inexpensive camera systems used for remote sensing," *International Journal of Remote Sensing* **39**(15-16), 4869–4890 (2018).

18. J. Nakamura, *Image sensors and signal processing for digital still cameras* (CRC press, 2017).

19. S.-H. Baek, I. Kim, D. Gutierrez, et al., "Compact single-shot hyperspectral imaging using a prism," *ACM Trans. Graph.* **36**(6), 1–12 (2017).

20. S. Garrido-Jurado, R. Muñoz-Salinas, F. J. Madrid-Cuevas, et al., "Automatic generation and detection of highly reliable fiducial markers under occlusion," *Pattern Recognit.* **47**, 2280–2292 (2014).

21. J. Jiang, D. Liu, J. Gu, et al., "What is the space of spectral sensitivity functions for digital color cameras?" in *2013 IEEE Workshop on Applications of Computer Vision (WACV)*, (IEEE, 2013), pp. 168–179.

22. OpenCV Documentation, "Camera Calibration," https://docs.opencv.org/4.x/d9/d0c/group__calib3d.html.

23. J. Weng, P. Cohen, M. Herniou, et al., "Camera calibration with distortion models and accuracy evaluation," *IEEE Trans. Pattern Anal. Machine Intell.* **14**(10), 965–980 (1992).

24. R. Hartley and A. Zisserman, *Multiple view geometry in computer vision* (Cambridge university press, 2003).

25. O. Sorkine-Hornung and M. Rabinovich, "Least-squares rigid motion using svd," *Computing* **1**, 1–5 (2017).

26. R. Halir and J. Flusser, "Numerically stable direct least squares fitting of ellipses," in *Proc. 6th International Conference in Central Europe on Computer Graphics and Visualization. WSCG*, vol. 98 (Citeseer, 1998), pp. 125–132.
27. I. Farrance and R. Frenkel, "Uncertainty of measurement: a review of the rules for calculating uncertainty components through functional relationships," *The Clinical Biochemist Reviews* **33**, 49 (2012).

28. I. BIPM, I. IFCC, I. ISO, and O. IUPAP, "Guide to the expression of uncertainty in measurement, jcgm 100: 2008, gum 1995 with minor corrections," *International Organization for Standardization Geneva ISBN (2008)*.

29. Labsphere, "8° Hemispherical Reflectance Calibration Certificate," https://terraref.ncsa.illinois.edu/clowder/files/5de554fb4f0cc808be61b522/blob/ (2019).

30. A. Lucat, R. Hegedus, and R. Pacanowski, "Diffraction effects detection for hdr image-based measurements," *Opt. Express* **25**(22), 27146-27164 (2017).

31. P. Santafé-Gabarda, L. Gevaux, A. Ferrero, et al., "Development and comparison of primary facilities for traceable bssrdf measurements," *J. Opt.* **26**(3), 035601 (2024).

32. X. Deng, F. Luan, B. Walter, et al., "Reconstructing translucent objects using differentiable rendering," in *ACM SIGGRAPH 2022 Conference Proceedings*, (2022), pp. 1–10.

33. C. Zhang, Z. Yu, and S. Zhao, "Path-space differentiable rendering of participating media," *ACM Trans. Graph.* **40**(6), 1–16 (2021).