Primary facility for traceable measurement of the BSSRDF

: The bidirectional scattering-surface reflectance distribution function (BSSRDF) is the function that describes the variation of the radiance of the elementary areas of a surface with respect to the directionally incident radiant flux on that surface. Measurements of the BSSRDF are important for characterizing the translucency of objects and for obtaining those optical parameters affecting volume scattering. However, to date, no traceable measurements of this function are available, since, likely due to its technical complexity, no standard measurement procedure has been established. We have developed a primary facility for measuring the BSSRDF based on a gonio-spectrophotometer with spatial resolution in the collection, spectral resolution in the irradiation, and angular resolution both for irradiation and collection directions. The BSSRDF of twelve homogeneous and translucent samples (with controlled values for the mean diameter of the scattering particles and their concentration) have been measured with relative uncertainties below 3% inside the irradiated area. Some results are shown and commented. This primary facility will allow the BSSRDF scale to be transferred to other instruments.


Introduction
The study and characterization of appearance has become very important in different industrial areas in the last decades.This has led to the International Commission on Illumination (CIE) to include it among the top 10 priorities in its research strategy, and to provide recommendations on the measurement of appearance [1].Color, gloss, translucency and texture must be quantified objectively by measurement scales to provide a complete description of the appearance of an object.These measurement scales rely on the measurement of the reflected and transmitted light by the objects, and its spectral, angular and spatial distributions.
Reflectance (or transmittance) are defined as the ratio of the reflected (or transmitted) radiant flux to the incident radiant flux in given conditions [2].These conditions include the size of the irradiated area on the surface (A i ), the size of the collection area of the surface for which the optical radiation is evaluated (A r ), the irradiation and collection solid angles (ω i and ω r ) or the irradiation and collection directions (defined by the polar, θ i and θ r , and azimuth, ϕ i and ϕ r , spherical coordinates), as well as some non-geometrical conditions such as the polarization [3][4][5][6].Measurements of reflectance or transmittance at specific conditions are only partial measurements and cannot be used for a general description of the reflectance or the transmittance of an object.In order to provide a general geometrical description for reflectance, Nicodemus et al. [7] defined in 1977 the Bidirectional Reflectance Distribution Function (BRDF) as a derivative relating the reflected radiance from a surface in a given direction (L r ) to the irradiance on it from another direction (E i ), by assuming a uniform and directional illumination with a large enough area on an homogeneous and isotropic surface: In this equation, f r denotes the BRDF.This function, when known for a number of well selected bidirectional geometries (pairs of irradiation and collection directions), allows the reflectance of opaque surfaces to be calculated by integration at any geometrical condition [7,8].However, it is insufficient for a complete description of translucent materials, since, in such case, the light scatters within the volume of the object and can emerge from positions outside the irradiated area, and this spatial distribution is not addressed by the BRDF.There is a spatial-dependent BRDF, the Bidirectional Texture Function or BTF [9], but it is defined to describe the texture of a surface, not the translucency.
A more general function is needed for describing these materials, the Bidirectional Scattering-Surface Reflectance Distribution Function (BSSRDF), which was also defined by Nicodemus et al. [7] as it follows: where f ssr denotes the BSSRDF and L ssr is the emergent radiance from an elementary area on the surface, located at x r , and in a specific direction, r r .This radiance is produced by the incident radiant flux Φ i at an elementary area located at x i from the direction r i , as it is shown in Fig. 1.Definitions of radiometric quantities can be found in the literature [10].Although it is known the relation between BSSRDF and reflectance and transmittance [11], the BSSRDF measurement has not been addressed in metrological terms and, to date, no recommendations from international institutions have been given.Some systems have been designed to measure the BSSRDF, mostly in the field of computer graphics [12][13][14][15][16][17][18][19][20], because there is a need of BSSRDF data for evaluating empirical models, for reproducing appearance and for obtaining the phase function and the scattering and absorption coefficients of materials.All these systems are based on camera acquisitions and illumination with a small irradiation area, being the main difference between them their capacity to realize irradiation and collection directions.However, these measuring systems do not provide traceable measurements of the BSSRDF, and do not allow the BSSRDF scale to be transferred to other instruments.In this article, we present a primary facility for measuring the BSSRDF, which should be beneficial for the general improvement of this kind of measurements, providing with reliable data to the scientific-technical community and giving a closer insight into reflectance and transmittance measurements.It is an adaptation of the gonio-spectrophotometer at IO-CSIC [21,22], originally conceived for BRDF measurements.The facility and method are described in this article, and, to show its potential, some results are given and discussed, those from the BSSRDF measurements of 12 homogeneous and translucent samples, with controlled values for the mean diameter of the scattering particles and their concentration, which provide a wide range of translucency levels.

Measuring system
A sketch of the gonio-spectrophotometer is shown in Fig. 2. It is composed of three subsystems: the irradiating system, which is fixed; the sample-positioning system, a 6-axis robot arm; and the collection system, a CCD camera on a circular rail.This design allows any irradiation and collection direction to be realized by a coordinates system transformation that is well described in the literature [23,24].An incandescence lamp (S) was used with a monochromator in Czerny-Turner configuration (Mc) (Bentham Instruments Ltd, model TMc300), with a 300 mm focal length and two diffraction gratings [1200 g/mm (250 nm -1200 nm) and 830 g/mm (500 nm -1800 nm)] to provide spectral resolution.Originally, the gonio-spectrophotometer was designed with a Köhler optical system to obtain a very uniform irradiation on the sample, needed in BRDF measurements when radiance is evaluated.However, for measuring BSSRDF it is required a small irradiation spot on the sample, because the smaller it is, the better insight into the spatial dependence of the BSSDRF is obtained.The design of the gonio-spectrophotometer allows the irradiated area to be modified by modifying the diaphragm size (P2) at the second lens of the optical system (L2), since this diaphragm is imaged on the sample plane by the third lens (L3).Uniformity is not needed for the irradiated area in BSSRDF measurements, because the relevant quantity is the incident radiant flux, which is integrated in the measurement throughout this irradiated area.Therefore, by neglecting uniformity, it was possible to add a convergent lens (L1) to focus the beam on the diaphragm P2 to increase the available flux.Thus, the modified optical system consists of three lenses (L1, L2 and L3) with focal lengths of 75 mm, 75 mm and 500 mm, respectively.
A filter wheel (FW) is located between L2 and L3, with five positions: clear position (no-filter), dark position (opaque filter) and three neutral-density filters with nominal transmittances of 10%, 1% and 0.01%.Before the filter wheel, an uncoated plate of fused silica (W) reflects about 8% of the incoming radiant flux towards a silicon photodiode (Mon), in order to control the variation of the radiant power of the light source.In addition, a diaphragm (P3) at L3 allows the irradiance on the sample and the solid angle of the incident beam to be modified.It was completely open, providing an irradiation solid angle of 3.5 × 10 −4 sr, approximately.After L3, the incident beam is steered by a first mirror (M45') oriented at -45 o , and a second mirror (M45) oriented at 45 o towards the sample (Sam).This periscopic arrangement allows retroreflectance measurements to be accomplished when the second mirror is replaced by a beam-splitter.
The sample positioning is realized by a 6-axis robot arm (R6) (Stäubli TX-40) that orientates automatically and quickly the sample to any direction.A vacuum sucker at the robot arm holds the samples.The possible effect in the measurement of the reflection of the transmitted light on the robot-arm holder was minimized by leaving 2 cm between the sample and the black surface of the holder.
Regarding the collection system, a scientific CCD camera (CCD) (Qimaging, model Rollera XR) is used as detector to provide with spatial resolution to the collection.This camera has a dynamic range of 12 bit, but it was extended to a higher dynamic range (HDR) by combining images with different integration times through an algorithm.The spatial resolution is 695 pixels × 520 pixels, with a pixel size of 12.9 µm × 12.9 µm.A Navitar Zoom 7000 18:108 mm Macro Lens was used with an f/4 number and a working focal length of 108 mm.The working distance was 360 mm, allowing a resolution on the sample plane of 41.8 µm per pixel.The camera is placed on a platform that can move through a cogwheel ring, which allows an automatic movement around the sample.

Measurement equation
The BSSRDF is defined theoretically in Eq. ( 2) as the derivative of the emergent radiance with respect to the incident radiant flux, in given directions, which can be interpreted in radiometry as the ratio between differential elements of radiance and radiant flux.Since it is not possible to measure differential quantities at the laboratory, this definition is arranged as the ratio between the average emergent radiance, within a narrow enough collection solid angle and a small enough measurement area on the surface of the sample, and the average incident radiant flux, within a narrow enough irradiation solid angle and a small enough irradiation area: In this context, 'narrow or small enough' means that the so measured value of the quantity (BSSRDF in this case) does not change for narrower or smaller sizes.
The emergent radiance, L ssr (x r , r r ), can be expressed in terms of the radiant flux reflected by the surface element located at x r in the r r direction as follows [10]: where A p is the apparent area of the surface element that is under evaluation, which depends on the collection polar angle, θ r .Thus, A r is the actual area of this surface element, which is determined by the area of the field of view of a pixel of the CCD on the sample plane (A fov ) and can be regarded as constant, and ω r is the collection solid angle, also a constant parameter of the measuring system.Thus, Eqs. ( 3) and ( 4) can be rearranged as: Radiant fluxes can be related with the responses of the pixels of the CCD sensor, which compose the images.These responses (N, in units of 'counts'), after dark subtraction, directly depend on the collected radiant flux (Φ), the external quantum efficiency of the pixels (η e ), the transmittance of the objective lens (τ) and the exposure time (t exp ): where the factor λ/hc (with h and c being the Planck and the speed of light in vacuum constants, respectively) represents the energy of a photon with wavelength λ, and F is the conversion factor between counts and photoelectrons (F = 5 photoelectrons/count for our camera).
It must be noticed here that the transmittance of the objective lens and the external quantum efficiency of the pixels are not the same when measuring the direct radiant flux (without sample) than when measuring the reflected radiant flux by the sample, which correspond to the two different radiant fluxes in Eq. ( 5).It is explained by the different ways the light reaches the camera: In the first case [Fig.3 (1)], the measuring solid angle is determined by the irradiation solid angle, which means that the light is reaching the camera only from directions concentrated near the normal direction, while in the second case [Fig.3 (2)], it is limited by the collection solid angle, so that the light is reaching the camera in all the attainable directions and it is transmitted in a different way by the objective lens.In addition, the pixels are irradiated under different beam geometries, and it affects their quantum efficiency.Therefore, these factors do not cancel out in the ratio of radiant fluxes when calculating the BSSRDF, as we will show.The measurement equation can be expressed in terms of dark-subtracted camera response as follows: where "k" refers to the pixel-related position on the surface of the sample, and "i" and "r" refer to the irradiation and collection directions, respectively.The summation of N k,i considers only the pixels that are being irradiated in the incident-radiant-flux configuration.This equation includes the transmittance of the neutral-density filter (τ nd ) used in this configuration.The aim of this filter is to obtain similar levels of counts in both incident-and reflected-radiant-flux configurations (see Fig. 3), avoiding in this way non-linearity issues.Each variable in Eq. ( 7) has an associated standard uncertainty of measurement, which is related to the relative standard uncertainty of the BSSRDF [u r (f ssr,k )].It can be expressed as follows: where 'u r ' denotes relative uncertainty.

Radiometric characterization of the facility
Specific experimental procedures were used for determining the value of each quantity involved in the measurement equation [Eq.( 7)] and its uncertainty.They are described in the following subsections.

Area of the field of view of a pixel, A fov
A vernier caliper was located and aligned as a sample, with the external jaws set at a distance of 20.00 mm ± 0.05 mm.From a camera image at that position, a value of 1.748 × 10 −9 m 2 was obtained for A fov , with a relative standard uncertainty of 0.54%.

Collection polar angle, θ r
The uncertainty of θ r is due to the angular inaccuracy of the system (robot arm and rotary stage of the camera) and the alignment inaccuracy, limited by the spot size of the alignment laser beam.
It was estimated to be 0.1 o (≊ 1.7 × 10 −3 rad).It impacts in the measurement trough the cosine [Eq.( 7)], and its relative standard uncertainty is given as: Note that this expression diverges for a collection polar angle of 90 o because in that case cos θ r equals zero.In practice, no measurement is attainable at that condition.

Dark-subtracted pixel response, N k
There are four main uncertainty sources for the dark-subtracted pixel response: photon noise of the optical radiation to be evaluated, readout noise, thermally-generated electrons and optical radiation strayed in the laboratory.To estimate them, two series of images at nine different integration times (from 1 ms to almost 66 s) were taken: one with the camera completely covered by its lid (covered-camera condition) and the other with the camera uncovered but shuttering the incident beam (uncovered-camera condition).The Fixed-Pattern Noise (FPN) [25], previously evaluated as an average at the minimum integration time (1 ms), was subtracted to all images.
Since photon noise follows a Poisson distribution, its uncertainty is given by the square root of the number of collected photons, unlike the readout noise (σ r ), which is independent of this number.Its contribution to the uncertainty was estimated as the standard deviation of the reading of the pixels at the minimum integration time (first image of the series), in the covered-camera condition.In this case, strayed photons are not present and the noise contribution of thermally-generated electrons is negligible due to the short integration time.A value of 6.87 electrons rms was obtained.The noise from the thermally-generated electrons follows a Poisson distribution too, and it depends on the thermal electrons rate (s t ).This is estimated by subtracting the mean signal of the pixels at the minimum integration time to those at the other integration times, in the covered-camera condition, and dividing by each integration time.At a room temperature of 23 o C, a mean rate of 0.73 electrons per second per pixel was obtained, which was stable for the longest integration times, as expected.Similarly, the strayed-photons noise depends on the stray-light rate (s s ).It is estimated by subtracting, for each integration time, the mean signal of the pixels in the covered-camera condition to that in the uncovered-camera condition, and dividing by the corresponding integration time.A mean rate of 6.69 electrons per second per pixel was obtained for our dark conditions and it was also stable for the longest integration times.These two latter noises contribute to the dark signal, being the difference that the strayed-photons noise does not appear when the lid of the camera is on the macro lens.
These values can be used to estimate the standard uncertainty due to the noise of the dark-subtracted pixel response with the following equation: where F is the conversion factor of the camera.A factor 2 appears in this equation because two acquisitions are needed to subtract the dark response.The relative uncertainty decreases inversely to the square root of the incoming photons.
The photoresponse non-uniformity (PRNU) [26,27] was characterized by using as uniform field an integrating sphere of 150 mm of diameter, with an exit-port diameter of 40 mm, irradiated externally by the irradiating system of the goniospectrophotometer.This PRNU was 4.9% with broadband illumination, and it was corrected for each image.The residual uncertainty is negligible with respect to the other uncertainty sources.Regarding the stray light between the objective lens and the CCD sensor, it is estimated from an image in the incident-radiant-flux configuration [Fig.3 (1)] as being negligible (it is shown in Fig. 4(a) that no significant signal is obtained outside the irradiated area).

Collection solid angle, ω r
A specific procedure has been designed to determine the collection solid angle.A mirror is placed in the robot arm, aligned as a sample, and oriented at 45 o from the incidence beam, so that the specularly-reflected radiant flux enters the camera located at 90 o .The mirror is sequentially orientated in different directions so that the radiant flux enters the objective lens from different directions.In order to sample uniformly the objective lens, we used 11 evenly spaced collection polar angles, θ r , from 40 o to 50 o , and 11 collection azimuthal angles, ϕ r , from 175 o to 185 o .The total collected radiant flux at each direction is assessed by the camera in terms of counts, and the collection solid angle is calculated from this angular distribution (Fig. 4(b)), using the surface integral definition of the solid angle as: where S(θ, ϕ) refers to the sum of the dark-subtracted response of all pixels of the image at each collection geometry.The integration is normalized to the central signal, S(θ 0 , ϕ 0 ), i.e., the signal at the specular direction (θ 0 = 45 o , ϕ 0 = 180 o ).
In practice, a discrete expression of Eq. ( 11) is used: The uncertainty of the collection solid angle is obtained by the uncertainty propagation method described in the literature [28].In this derivation, the uncertainty of the relative angular displacements has been omitted because it is negligible with respect to the uncertainty of the absolute angular position determined by alignment [u(θ n ) = 0.1 o ].

Responsivity ratio, τ
The transmittance of the objective lens (τ) multiplied by the external quantum efficiency of the pixel (η e ) is proportional to the pixel responsivity.As explained before, these two variables can vary with the angular distribution of the incoming radiant flux, and, therefore, they are not canceled out when calculating the ratio between radiant fluxes in Eq. (7).In consequence, the ratio τ i η e,i /τ r η e,r must be accounted in the measurement equation.It can be calculated by using the data represented in Fig. 4(b), where the summation of values within 44 o ⩽ θ r ⩽ 46 o and 179 o ⩽ ϕ r ⩽ 181 o was considered directly related to the responsivity of the detection system in the incident-radiant-flux configuration [Fig.3 (1)], while the summation across all geometries included in the collection solid angle was considered directly related to the responsivity in the reflected-radiant-flux configuration [Fig.3 (2)].It was calculated using the following expression: where S(θ n , ϕ m ) are the signal values and M c and M t refer to the amount of collection directions included for averaging incident-radiant-flux responsivity and reflected-radiant-flux responsivity, respectively.
The relative uncertainty has been estimated as the standard deviation of the M c signal values used for the incident-radiant-flux averaging normalized to their mean.

Exposure time, t exp
Different exposure times were used for each pixel as a function of the response.When averaging several acquisitions (multiexposition), the exposure time is the sum of the integration times of the single acquisitions.In general, a minimum integration time of 10 ms was used to avoid non-linearity issues, which are usually more relevant when shorter integration times are used [29].Also, similar integration times were used for direct-and reflected-flux configurations, by using a neutral-density filter in the former one.The time resolution and its uncertainty are very small with respect to the minimum integration time used, and its contribution to the total relative uncertainty is negligible.

Transmittance of the neutral-density filter, τ nd
The spectral transmittance of the neutral-density filter was calibrated in the measuring system itself, having a spectrally-averaged transmittance of 10.93%, and a relative standard uncertainty of 0.25% [21].

Results and discussion
We have measured the BSSRDF of the 12 translucent samples shown in Fig. 5.They were manufactured by Covestro Deutschland AG, with controlled values for the mean diameter of the scattering particles and their concentration in mass (see Table 1).All samples have a thickness of 6.4 mm and they present a wide range of translucency levels, from almost transparent to almost opaque.The BSSRDF was measured at 40 different bidirectional measurement geometries, with two polar angles of incidence θ i = {15 o , 45 o } and only one azimuthal angle of incidence ϕ i = 0 o , and some collection directions obtained by combination of polar angles with values These angles are defined at the sample reference system as shown in Fig. 6(a).In the following representations, reflectance symmetry with respect to the incidence plane was assumed to show the complete reflection hemisphere.The facility allows for spectral assessment, although discussion on the spectral dependence of the studied samples is not given here.A coordinate system transformation from the camera reference system (CRS) to the sample reference system (SRS) is used to relate the position of every pixel in the image (x ′ i , y ′ i ) with the position on the sample surface (x i , y i ), which can also be defined by polar coordinates (a and α) as shown in Fig. 6(b).The origin of the SRS is related to the centroid of the irradiance image (Fig. 4(a)) (x ′ 0 , y ′ 0 ) for each irradiation wavelength.The observed projection due to the collection angle θ r must be compensated, and the data need to be rotated according to the value of the azimuth angle ϕ r , as it follows: Some quantities describe exclusively the measuring system, and are independent of the sample to be measured.They were experimentally characterized according to the methods described in section 4. Their values and relative standard uncertainties are given in Table 2. Apart from the response of the pixels, the error source that predominate in the relative uncertainty of the BSSRDF is the determination of the responsivity ratio (τ i η e,i /τ r η e,r ), while the determination of the collection solid angle or the cosine of the collection polar angle contribute with a comparatively small relative standard uncertainty.Some results are given here to convey the measurement capability of the system.The dependence of the BSSRDF on the distance a to the center of the irradiated area (distance between A i and A r in Fig. 1) is shown in Fig. 7(a), for each sample described in Table 1 and shown in Fig. 5, at the measurement geometry with θ i = 15 o , ϕ i = 0 o , θ r = 10 o and ϕ r = 0 o , irradiated with λ = 550 nm.Their relative standard uncertainty is shown in Fig. 7(b).The BSSRDF values represented in Fig. 7(a) provide information about the translucency level of each sample.Those samples most opaque (S1, S2 and S3) and most transparent (S7, S8 and S9) in appearance (see Fig. 5) present steeper slopes, while the slopes of the most translucent samples (S5, S6, S11 and S12) are much smoother.These slopes characterizing the samples must be examined outside the irradiated area, because the measured value of the BSSRDF within it might depend on the size of this area [11], and, in that case, it does not characterize exclusively the optical properties of the sample anymore.In order to address this dependence, measurements must be carried out with different irradiation areas, what is allowed by our measuring system by using different diameters for P2 (Fig. 2).
The relative standard uncertainty of the represented measurement geometry is minimum within the irradiated area (between 2% and 3%, depending on the sample), and it continuously increases toward outer positions due to the decreasing emergent radiance, i.e., the lower response of the pixels.The more translucent the sample is, the larger the radiance for a given outer position and the lower the relative uncertainty.
To show in more detail the light interaction with the sample, the angular distribution of the BSSRDF (which can be understood by looking at the red dome of Fig. 6(b)) of five different positions on the surface of the sample S6 (see Fig. 5) is represented in Fig. 8.This representation corresponds to an irradiation at θ i = 15 o with λ = 550 nm.The color code represents the value of the decimal logarithm of the BSSRDF at each collection direction.In this kind of images, the Cartesian coordinates (X d , Y d ) are related to θ r and ϕ r as it follows [30]: where θ r must be expressed in radians.These five positions represented in Fig. 8 correspond to surface elements placed at the same distance from the irradiation position (a = 2.5 mm), but at different polar angles, α, on the plane of the sample, so that the observed variations depend only on the direction of the light propagation inside the sample.
In general, according to the presented measurements, the BSSRDF tends to increase for higher collection polar angles.In addition, the maximum value uses to appear when the azimuthal collection angle coincides with α, which agrees with the conclusions from Donner et al. [30] on simulated data, showing that maximum values of the angular distribution have to be distributed around the single scattering plane when low-order scattering events dominate.These angular distributions are related to the phase function (angular distribution of light intensity scattered by a particle at a given wavelength) and the scattering coefficient.The availability of reliable BSSRDF measurements will be very useful to determine them with the support of the Radiative Transfer Equation (RTE) [16,31].

Conclusions
The objective of this work is to make available to the scientific-technical community a primary facility to provide BSSRDF measurements with traceability to the International System of Units (SI).This quantity is very important for studying the scattering in materials, and for understanding and reproducing translucency.This primary facility has been developed, characterized and tested at IO-CSIC.The error sources limiting the relative uncertainty of the BSSRDF measurement are the response of the pixels and the determination of the ratio of the responsivities of the camera in direct-and reflected-flux conditions.The BSSRDF of twelve homogeneous and translucent samples, with a wide range of translucency levels, has been measured in the facility.The BSSRDF angular distributions, obtained for different positions on the surface with respect to the irradiated area, can be related to the phase function and to the scattering coefficient.Therefore, the availability of reliable BSSRDF measurements must be very useful to determine these quantities with support of the Radiative Transfer Equation.

Fig. 1 .
Fig. 1.Geometrical variables involved in the description of the BSSRDF measurement.

Fig. 4 .
Fig. 4. Spatial distribution of the irradiated area on the plane of the sample.

Fig. 5 .
Fig. 5. Translucent samples from Covestro Deutschland AG used in this work.

Fig. 8 .
Fig. 8. Angular distribution of the BSSRDF of five different positions on the surface of the sample S6 when it is being irradiated at θ r = 15 o with λ = 550 nm.