Open Access
Add to BibSonomy
Abstract
Non-invasive determination of the optical properties is essential for understanding the light propagation in biological tissues and developing optical techniques for quality detection. Simulation-based models provide flexibility in designing the search space, while measurement-based models can incorporate the unknown system responses. However, the interoperability between these two types of models is typically poor. In this research, the mismatches between measurements and simulations were explored by studying the influences from light source and the incident and detection angle on the diffuse reflectance profiles. After reducing the mismatches caused by the factors mentioned above, the simulated diffuse reflectance profiles matched well with the measurements, with R2 values above 0.99. Successively, metamodels linking the optical properties with the diffuse reflectance profiles were respectively built based on the measured and simulated profiles. The prediction performance of these metamodels was comparable, both obtaining R2 values above 0.96. Proper correction for these sources of mismatches between measurements and simulations thus allows to build a simulation-based metamodel with a wide range of desired optical properties that is applicable to different measurement configurations.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Spectroscopic techniques have been widely combined with chemometrics methods to build quality prediction models for various agricultural and food products [1–4]. These models estimate the link between the obtained spectra and the corresponding quality parameters with the help of mathematical techniques. As the acquired spectra are influenced by factors such as light source condition, detection mode and type of instrument [5], transfer of the prediction model from one setup to another is not straightforward. Therefore, a new prediction model has to be built for every new setup, which greatly increases the workload and costs. In response to these limitations, the bulk optical properties (BOP), which are independent of the factors mentioned above, have been proposed as an alternative option [5]. The BOP are mainly composed of the bulk absorption (μa) and reduced scattering coefficient (μs’), which respectively reflect the chemical and physical properties of the tissue [6]. Therefore, separating the BOP is expected to improve the performance of models to predict the quality properties of interest. In general, the field of building measurement configurations and quality detection models is in urgent need of robust and reliable techniques to determine the BOP of complex (biological) materials [7].
The optical methods to detect the BOP can be mainly classified into four categories: integrating sphere [8], time-resolved [9], frequency-domain [10] and spatially resolved [11] measurements. Among these methods, the integrating sphere measurement is the only one that requires destructive sample preparation. However, thanks to the rigorous control of the experimental conditions [12], accurate BOP can be obtained. Therefore, the integrating sphere measurement, which is usually regarded as the “gold standard” for the determination of BOP [13,14], is often used as a reference measurement. From the other methods, the time-resolved and frequency-domain methods are less popular than the spatially resolved measurement, due to their relatively expensive instrument and complex signal conversion [15]. The spatially resolved measurement can be further divided into two categories according to the sensing configuration: optical fiber arrays and reflectance imagery [16]. The limited number of detection points and the requirement for a good contact between the tissue and the fiber [15], which is actually not easy to control, make the optical fiber arrays less suitable for measurements on solid samples compared to the reflectance imagery [7].
After measuring the spatially resolved diffuse reflectance profiles, light transport models are required to link the profiles to the BOP. The most frequently used model is the diffusion equation, an analytical approximation of the radiative transfer theory [17]. However, the requirements of the diffusion approximation are typically not met in the region near the illumination source, where the strongest reflectance signals can be acquired [18]. Moreover, the equation is only applicable for scattering dominated light propagation (μs’ >> μa) [7], which is not always true for biological products. As a consequence, these conditions restrict the applicability of the diffusion equation. Recently, a library concept, also referred to as lookup table, has been proposed to link the diffuse reflectance profiles with the BOP [5]. By interpolating between the values in these lookup tables, metamodels can be created, which can then be used to estimate the BOP [19–22]. The data for training these metamodels can either be generated with experimental measurements [11,21] or with Monte Carlo (MC) simulations [20,23,24]. The experimental measurements-based metamodels are typically trained on data acquired for samples with known BOP (e.g., optical phantoms). This approach allows to estimate the BOP of new samples without having to separately consider the inherent noise and response of the setup system [20,25]. Nevertheless, the prediction capabilities are limited to the same setup for which the model was built, which limits the flexibility of such models. In addition, different combinations of BOP are typically designed by adjusting the concentrations of the absorbing and scattering components. As the BOP of the absorber and scatterer are wavelength-dependent, a large set of samples has to be measured to create a metamodel covering a wide range of BOP values. On the other hand, the MC simulations-based metamodels are more flexible, as simulations can be performed for any desired combination of BOP independent of the wavelength [5]. In the past, the long computation times were a major shortcoming of the MC simulations, but thanks to the development of the Graphics Processing Units (GPU)-based MC, the required time has been greatly reduced.
Several researchers have reported on the estimation of BOP from spatially resolved diffuse reflectance data. Palmer et al. [25] evaluated a Monte Carlo simulation-based inverse model and validated its performance on two sets of optical phantoms, obtaining average fitting errors of 3% and 12%. Similarly, good agreement between the predicted BOP obtained from the Monte Carlo simulation-based inverse model and the expected values was also reported by Rajaram et al. [19] and Hennessy et al. [20]. Both of these researches set the parameters of the fiber geometry fully in the simulation during the forward model generation. The good matching performance indicated that Monte Carlo simulations are a powerful tool to accurately simulate the light propagation in biological tissues. However, as mentioned above, fiber arrays require good contact between the fiber and the sample, which is not easy to control in practice. For the noninvasive and noncontact estimation of the optical properties, other researchers also used CCD cameras to capture the diffuse reflectance profiles. Kienle et al. [15] built a neural network based on the Monte Carlo simulation results and then used the network to predict the optical properties from the diffuse reflectance profiles, obtaining an rms error of 14% for μa and 2.6% for μs’ over the range of 0.02 - 1 cm-1 for μa and 5 - 20 cm-1 for μs’. Zhang et al. [23] trained two neural networks respectively on Monte Carlo simulations and experimental results. The simulation-based neural network applied on simulated data resulted into a mean relative error below 2.5%, while the measurement-based network applied on measured data provided a relative error of 9% for μa and 3% for μs’. However, the relative error of the bulk optical properties of measured samples that were estimated with the simulation-based network could be up to 50% for μa and 25% for μs’, which indicated that the well-trained network based on the simulations could not estimate the optical properties of the measurements properly.
One reason for the large deviation could be caused by the oversimplification of the measurement configuration in the simulations. For instance, the light beam in the simulation is presumed as an infinite narrow beam, while in the real measurement, the light beam has its power, size and shape. Moreover, the illumination angle in the measurement configuration is typically neglected in the simulations. Although, several researchers have reported that an illumination angle within 5 ° - 10 ° [15] will not cause a significant difference in the diffuse reflectance profiles, the detailed influences from such factors haven’t been published. In addition, the detection angle is usually restricted by the settings of the measurement configuration, while in the simulation, all the photons exited at all directions were considered by default. Therefore, in this research, the influences from the key elements (light source, incident and detection angle) of a measurement configuration on the diffuse reflectance mismatch between measurements and simulations will be first explored.
Another possible reason for the difference in the prediction performance could be attributed to the noise in the measurements that was not considered in training the simulation-based neural network. The differences between the calculated and desired output values were generally used to update the connection weights in the neural network during the training process [23]. For the purely simulation-based neural network, the stochastic noise in the simulations was the major source of the difference which can be reduced by increasing the number of photons and was far smaller than the noise in the measurements. Therefore, the larger noise in the measurements being ignored in the simulations would make a difference in the input data, resulting in a poor prediction performance on the measured profiles. To incorporate the measurement noise into the model, Aernouts et al. [11] built a stochastic Kriging (SK) metamodel based on the diffuse reflectance profiles obtained from the imaging measurements to inversely predict the optical parameters. They considered both the intrinsic and extrinsic uncertainty in their SK metamodel. The extrinsic uncertainty refers to the inherent stochastic noise produced in the Gaussian process and the intrinsic uncertainty is represented by the standard deviation of the repetitions in the measurements. Such a measurement-based SK metamodel can closely link the BOP to the measured diffuse reflectance profiles and thereby obtain a good prediction performance (R2 above 0.98), which is especially suitable for predicting the BOP from diffuse reflectance profiles acquired with a fixed measurement configuration. Given the good fitting performance reported for the SK metamodel, its potential for building measurement-based and simulation-based metamodels will be explored.
Although the constraints mentioned above could in theory be lifted by changing the settings (e.g., the incident and detection angle) in the original C code of the simulations, this would mean that the simulation-based metamodel would have to be rebuilt for every small change in the detection configuration. It is hypothesized here that a basic metamodel could be made applicable to different detection configurations with a few modifications. The main objectives of this research were to: (1) explore the influences from the light source and the incident and detection angle on the mismatch between measurements and simulations; (2) build metamodels to link the spatially resolved diffuse reflectance profiles with their respective BOP both for the measurements and the simulations and compare their performance; (3) evaluate the performance of the metamodels to estimate the BOP from measured diffuse reflectance profiles.
2. Materials and methods
2.1 Flowchart
The steps followed to investigate the main objectives of this research are illustrated in Fig. 1. For the 1st objective, the influences from the light source and the incident and detection angle on the diffuse reflectance profiles were explored by changing the settings in the Monte Carlo simulations based on the BOP measured for a few representative calibration samples. Meanwhile, the diffuse reflectance profiles of these samples were measured and then compared with the values obtained in the simulations. The 2nd objective consisted of building the measurement-based and simulation-based metamodel and comparing the fitting performances. The measurement-based metamodel was built based on the BOP and diffuse reflectance profiles measured for the calibration samples, while the simulation-based metamodel was trained based on the designed BOP and their corresponding simulated diffuse reflectance profiles. The 3rd objective was to evaluate the prediction performance of these inverse metamodels. The measured diffuse reflectance profiles of the validation samples were fed to the measurement-based and simulation-based inverse metamodel, and the corresponding BOP obtained from the inverse metamodels were then compared with the measured BOP. The detailed information for each of these steps is discussed in the following sections.
Fig. 1. Flowchart of the steps that were followed to reach the (1) 1st, (2) 2nd and (3) 3rd objective.
2.2 Sample preparation
Naphthol Blue Black (195243, Sigma Aldrich, Missouri, USA), Intralipid 20% (Fresenius Kabi, Sweden) and distilled water were respectively used as the absorber, scatterer and diluent to prepare a series of liquid optical phantoms. To obtain a wide range of BOP combinations, 7 concentration levels of absorber (15, 35, 55, 75, 115, 155 and 195 μM, marked as 1 to 7) and 7 volume concentration levels of scatterer (1.135, 2.270, 3.405, 4.540, 5.675, 6.810 and 7.945%, marked as A to G) were mixed to produce 49 liquid phantoms (black rectangles in Fig. 2) for the measurement-based metamodel building. To evaluate the prediction performance of the metamodels, an extra group of 12 validation phantoms (red dots in Fig. 2) with 6 concentration levels of absorber (20, 45, 70, 95, 135 and 175 μM) and 6 volume concentration levels of scatterer (1.816, 2.724, 3.632, 4.994, 6.356 and 7.718%) were prepared.
Fig. 2. Distribution of the absorption and scattering level of the liquid phantoms (absorption level increases from 1 to 7, scattering level increases from A to G). Samples marked with black rectangles were used to build the measurement-based metamodel; Samples marked with red dots were used to validate the measurement-based and simulation-based metamodels.
2.3 Measurement configuration
A supercontinuum laser (SC450-4, Fianium Ltd., Southampton, UK) with 4 W output power and a spectral coverage from 440 nm to 2400 nm, was coupled into a monochromator (Oriel Cornerstone 260 ¼ m, Newport, Irvine, USA) to produce monochromatic light (530 nm to 970 nm in steps of 5 nm) for both the double integrating spheres (DIS) (RT-060-IG, Labsphere Inc., North Sutton, USA) and the diffuse reflectance imaging measurements. As shown in Fig. 3, after leaving the monochromator, the light at the target wavelength was split by a beam splitter. A small fraction of the light was directed to a Si detector (PDA100A, Thorlabs Inc., New Jersey, USA) to monitor the stability of the laser system, while the rest of the light was controlled by two mirrors to direct to the sample path with three different directions. If the first mirror was up, then the light (marked as solid lines) was guided to the DIS to obtain the total diffuse reflectance and total transmittance of the sample in a cuvette (sample thickness of 0.58 mm) sandwiched between the two spheres. These signals were recorded by Si detectors installed on the spheres. If both the mirrors were down, then the light (marked as dotted lines) was sent to illuminate the sample in a cuvette (sample thickness of 0.16 mm) in the unscattered transmittance measurement path. A Si detector was installed 1.5 m behind the sample to measure mainly the unscattered photons. If only the second mirror was up, then the light (marked as dashed lines) was focused into an optical fiber (multimode, NA 0.22, 200 μm core diameter) for the diffuse reflectance imaging measurements. To exclude the specular reflectance from the diffuse reflectance measurements, the end of the fiber was fixed at an angle of 20 ° with respect to the normal of the sample surface. A CCD camera (TXG-14NIR, Baumer, Frauenfeld, Switzerland) in combination with an extender (2-EX, Pentax, Tokyo, Japan) and a 35 mm fixed focal length lens (67-716, Edmund Optics, New Jersey, USA) with a f/4 aperture, installed right above the sample stage, was used to capture the diffuse reflected light. The detectable angle of the measurement configuration was calculated to be around 1.5 °. For more details about the configuration, the reader is referred to Aernouts et al. [26] and Van Beers et al. [27].
2.4 Bulk optical properties’ estimation
The inverse adding-doubling (IAD) routine developed by Prahl et al. [28] was used to estimate the BOP from the DIS and unscattered transmittance measurements. The parameters used in the routine are detailed in Aernouts et al. [26]. The software code combining these parameters and the measurement results with the IAD routine was executed in Matlab (version 9.1, The Mathworks Inc., Massachusetts, USA). Three repetitions were performed for each phantom and the average values were used in the further analyses.
2.5 Spatially resolved diffuse reflectance measurement
Unlike the phantoms being filled in a thin glass cuvette for the DIS and unscattered transmittance measurements, the phantoms for the diffuse reflectance imaging measurements were loaded in a cylindrical container with a diameter of 3 cm and a depth of more than 6 cm to mimic the light propagation in a semi-infinite object and reduce the influence from the boundaries of the container. As a slight deviation in the distance between the sample surface and the optics (camera and illumination) would lead to a difference in the output intensity as well as the position of the incident beam and the observed reflectance profiles, this distance was adjusted for each sample until the center of the light spot was at the target position in the image before starting the measurement. During a measurement, images were successively acquired at different wavelengths by changing the monochromator setting without moving the sample. Thanks to this procedure, the center of the light spot could be assumed to be the same for all the wavelengths. As the power of the monochromatic light differed from wavelength to wavelength and the samples also absorb light at specific wavelengths, this resulted in variation in the illumination intensity for different wavelengths. The region of the illuminating beam in the image showed significantly higher reflectance signals. Accordingly, the region with reflectance intensities higher than 1/3 of the maximum value was first selected to determine the center of the light spot for each image. The average position of the spot in the high signal wavelength range 650 nm - 800 nm was then used to calculate the center for all the wavelengths (530 nm - 970 nm). After obtaining the position of the center, the pixels with the same distance from the center in the region of interest were averaged and the average value was regarded as the diffuse reflectance at that distance. The definition of this region of interest will be further discussed in the following section. Thus, by this step, the 2-dimensional image has been converted into a diffuse reflectance profile as a function of distance. The average diffuse reflectance profile of 3 repetitions for each phantom was retained for the further data analyses.
2.6 Mismatch between measurements and simulations
Among the prepared samples, 4 phantoms with the max/min concentration levels (A1, A7, G1 and G7) were selected for exploring the diffuse reflectance mismatch between measurements and simulations. This exploration was mainly focused on 3 aspects: light source, incident angle and detection angle.
The influence from light source mainly consists of two aspects: the power and the size of the light beam. The light beam here refers to the light incident on the sample surface. As shown in section 2.3, the light has to go through a series of optical elements when moving from the laser to the sample. So, a small change in any of the elements can alter the light beam incident on the sample surface. In the MC simulations, the power, size and shape (fixed in this research) are used to describe the light source. To make the simulations interchangeable among different configurations, it is necessary to know how these factors affect the diffuse reflectance profiles.
The incident angle refers to the angle that the illuminating beam makes with the normal of the sample surface, as shown in Fig. 4(a). The default setting in the simulations is an illuminating beam that is perpendicular to the sample surface. As in the real measurements, the relevant diffuse reflectance profiles with sufficient signal are relatively close to the illumination spot and the size of the camera is nonnegligible, an angle between the incident light and the normal of the sample surface is unavoidable. Together with this angle, the shape of the light spot changes, causing the central symmetry to be lost. Therefore, the influence from the region of interest (the selected part of the image spreading outward from the center of the light spot in Fig. 4(b)) and the incident angle were studied together.
The detection angle refers to the angle between the line that links the light spot to the edge of the camera aperture and the sample surface normal, as shown in Fig. 4(a). The escaped photons, being counted into the diffuse reflectance profiles in the simulations, can exit the sample in a wide range of directions. However, in the real measurement, the camera is generally installed at a fixed position. Therefore, it is necessary to explore the spatial distribution of the exited photons. In addition, the size of the aperture and the distance between the sample and the camera can differ from configuration to configuration, both determining the detection angles. Therefore, the influence from the detection angle on the diffuse reflectance profile was further explored.
The effect of light source and incident angle were studied in simulation using the Monte Carlo in voxelized media (MCVM) code [29], while the influence of detection angle was studied using CUDAMCML [30].
2.6.1 MCVM
MCVM is a Monte Carlo code that allows to simulate the light propagation in a 3-dimensional heterogeneous medium [29,31]. In MCVM, a sample is split into small voxels in the 3-dimensional structure and a combination of optical properties (e.g., μa, μs’ and refractive index), which corresponds to a tissue type, is assigned to each voxel. This technique allows to simulate the light propagation in tissues with irregular boundaries between the tissue types with different BOP and refractive indices. However, this property was not explored in this study. Two other advantages of this, namely the flexibility in defining the incident light and detailed recording of the state of each escaped photon, were fully exploited.
MCVM simulations with 10,000,000 photons were used to explore the influences from light source and incident angle for the BOP measured for optical phantoms A1, A7, G1 and G7 at 650 nm. The BOP of these phantoms were obtained from the DIS and unscattered transmittance measurements. In exploring the influence from the light source, the illumination beam was assumed to be perpendicular to the sample surface corresponding to an incident angle of 0 °. The region of interest initially covered the full 360 ° of the sample plane, while the detection angle corresponded to a cone of photon exit angles from 0 to 90 ° with the sample surface normal. Additionally, incident angles of 10 °, 20 ° and 30 ° were respectively used in the simulation to explore the influence from the incident angle. During this process, the detection angle still covered 0 to 90 ° of the photon exit angle, while the region of interest, as shown in Fig. 4(b), was divided into 4 parts according to the relative position to the illumination spot: left, right (perpendicular to the plane formed by the illuminating beam and the normal of the sample surface), up (towards the fiber) and down (away from the fiber). As the light spot was still symmetrical in the direction of the illumination fiber, the left and right part were the same. Therefore, only the left, up and down part with specific extraction angles of 30 °, 60 ° and 90 ° were selected for the analyses. The minimal extraction angle of 30 ° was selected to maximize the differences in the different directions and ensure sufficient signal in the selected region of interest. To make the comparison clearer, the profiles were limited to a source-detector distance of 0.15 cm.
2.6.2 CUDAMCML
In studying the influence of the detection angle, the cone of photon exit angles (0 to 90 ° with sample surface normal) was equally divided into 60 groups to include the 1.5 ° detectable angle of the measurement configuration. However, the number of photons captured in such a small detection angle segment of 1.5 ° was quite limited. Therefore, to enhance the stability of the simulation results, CUDAMCML simulations (see below) with more photons were used for this exploration. During this exploration, the incident angle was set at 0 ° and the region of interest was selected as the one that gave the best match between simulated and measured values, which was selected in studying the influence from incident angle. In studying the influence from the angular detection range, 0-1.5 °, 0-6 °, 0-15 °, 0-30 °, 0-60 ° and 0-90 ° were selected for the simulations, the corresponding simulated diffuse reflectance profiles were then compared with the measured values to further explore the influence from the detection angle on the mismatch between measurements and simulations.
Since the MCVM simulations output the detailed state of all the photons, which took a large amount of time and space, CUDAMCML was also used to simulate the inputs for the metamodel. Considering the possible fluctuation in the results of the statistical method, 3 runs were performed for each simulation and the average value was used in further analyses. To make the results more accurate, the number of photons used in the CUDAMCML simulation was increased 20 times compared to the number used in the MCVM simulations, and much higher than the values reported by other researchers [20, 32–34]. The resolution in the depth and radius direction were both set at 0.001 cm, which were the same as the resolutions used in MCVM. Based on this resolution, the number of grids in the radius direction was 150 considering a maximal radius of 0.15 cm. The detectable angle of the measurement configuration was calculated to be around 1.5 °. Hence, the number of grids in the detection angle was set at 60, resulting in detection angle segments of 1.5 ° each. The refractive index of the layer was set at the value reported for water at 650 nm [35]. The layer thickness was set at 6 cm, similar to the value used in the measurements mentioned in section 2.5. To allow a comparison of the metamodels, the BOP in the simulations were set to be similar to the values used in the measurements. Therefore, the anisotropy factor (g) value was set to 0.7, which is close to the values obtained from the measurements at 650 nm. The range for μs’ and μa was first determined close to the BOP of the optical phantoms used in this research, and then the range for μs’ was converted into a range for μs according to μs= μs’/(1-g). The parameters used in the CUDAMCML simulations are summarized in Table 1.
Table 1. Inputs for CUDAMCML simulations
2.6.3 Convolution
The light beam in the simulation mentioned above was regarded as an infinite narrow beam, which is different from the one in the measurements. Therefore, the convolution program (CONV) developed by Wang et al. [36] was used to simulate the response of the light beam with finite size based on the output of MCVM and CUDAMCML. The power, shape and size of the light beam are three important factors in CONV. The signals obtained in the measurements are relative reflectance values, which involve the conversion between the optical and electrical signals. As the signals obtained in the simulation are absolute values in numbers of reflected photons, a scaling factor between these two has to be introduced (discussed in the following section). Therefore, a fixed input power of 1 J was used in the CONV program. A Gaussian beam was considered to simulate the monochromatic illuminating light beam of the measurements. The size of the light beam (1/e2 radius) was determined by the method detailed in Kienle et al. [15]. A mirror was placed at the sample position for the diffuse reflectance measurement. The mirror was positioned perpendicular to the incident light beam. The distance between the mirror and the camera was controlled by adjusting the center of the light spot to the fixed position in the image, similar to the sample measurements. The radius of the monochromatic light beam at 650 nm was calculated to be 0.03 cm. The other parameters (e.g., diffuse reflectance profiles, resolution of radius and number of grid elements) used in the CONV were kept the same as the ones in the simulation output.
2.7 Measurement-based and simulation-based metamodels
Stochastic Kriging (SK) surrogate modeling in the ooDACE toolbox was applied to link the BOP and the diffuse reflectance profiles with a metamodel, as proposed by Watté et al. [21]. SK is commonly based on a Gaussian process and uses the distances between the unknown point and its surrounding points to calculate their correlations and thereby determine the value of the unknown point [37]. During this process, the intrinsic and extrinsic uncertainty are both incorporated in the model [38]. To construct such a model, a set of Kriging models was first built based on a set of initial features, which were selected by a Bayesian approach. After that, new features were successively incorporated into the model and the model performance was evaluated by the maximum likelihood method for each incorporation. The model kept updating until it obtained a high accuracy and flexibility at the same time. Detailed information about the SK model can found in Watté et al. [21] and Chen et al. [39].
In this study, metamodels were respectively built based on the measured and simulated data. In building the measurement-based metamodel, the diffuse reflectance profiles of the 49 phantoms (black rectangles in Fig. 2) with their corresponding BOP in the wavelength range from 530 nm to 970 nm were selected as the inputs. Although the laser power varied with the wavelength, the wavelength was not considered as a variable in the metamodel. To minimize the effect of this wavelength-dependent intensity variation, the measured diffuse reflectance values at a certain wavelength were scaled by the values acquired in an integrating sphere for the same wavelength. These normalized values were then used in the metamodeling. In building the simulation-based metamodel, the designed BOP in Table 1 and their corresponding simulated diffuse reflectance profiles were used. 27 source-detector distances in steps of 0.00115 cm and 31 source-detector distances in steps of 0.001 cm located between 0.04 cm and 0.07 cm were respectively selected to build the measurement-based and simulation-based metamodel. In addition, the standard deviation on the diffuse reflectance profiles for the 3 measured repetitions was used as the intrinsic uncertainty for the measurement-based metamodel. In the case of the simulation-based metamodel, the deviation between measurements and simulations was also included in the intrinsic uncertainty.
The evaluation of the metamodels was quantified for three sample sets: calibration, validation and prediction samples. To make the figures clear, only parts of the samples were selected for the evaluation. The phantoms that combine 4 absorption levels (1, 3, 5 and 7) and 4 scattering levels (A, C, E and G) on the whole wavelength range were selected for the calibration of the measurement-based metamodel, while 19 absorption levels (from 1.61 cm-1 to 16.01 cm-1, in steps of 0.8 cm-1) and 26 reduced scattering levels (from 10.8 cm-1 to 85.8 cm-1, in steps of 3 cm-1) were used to calibrate the simulation-based metamodel. The 12 validation phantoms (red dots in Fig. 2) and 49 calibration phantoms (black rectangles in Fig. 2) at 650 nm were respectively used for the validation and prediction of the simulation-based metamodel. During this procedure, the measured BOP of these phantoms were selected as the inputs for the simulation-based metamodel and the predicted diffuse reflectance profiles, after dividing by a scaling factor (introduced in the following section), were compared with the measured values. To make a comparison between these two metamodels, the 12 validation phantoms at 650 nm were also used for validating the measurement-based metamodel, meanwhile, the 49 calibration phantoms at 650 nm were also selected as a separate prediction set for the measurement-based metamodel. To clearly present the comparison results, the intensities of the diffuse reflectance profiles were plotted on a logarithmic scale.
When the simulation-based metamodel was applied on the measurements, a scaling factor was introduced to eliminate the scale difference between measurements and simulations. Considering the saturated values at source-detector distances below 0.04 cm for samples with high scattering levels (e.g., G1 and G7) and the low signal to noise ratio beyond 0.07 cm for samples with high absorption levels (e.g., A7 and G7), the scaling factor was calculated as an average value of the ratios between simulations and measurements for phantoms A1, A7, G1 and G7 in the 0.04-0.07 cm source-detector range.
2.8 Inverse metamodels
With the metamodel, the diffuse reflectance profiles can be predicted based on the BOP. However, the objective of this research was to predict the BOP from the acquired diffuse reflectance profiles, which requires inversion of the metamodel. Therefore, the metamodel was implemented in an iterative optimization algorithm. This iterative process starts from an initial guess of the BOP and keeps updating the BOP until the corresponding diffuse reflectance values predicted by the metamodel match the acquired diffuse reflectance values. The following cost function Eq. (1) is iteratively minimized to select the set of BOP that gives the best fit to the acquired diffuse reflectance profile.
Where N represents the number of source-detector distances; ${I_{i,meas}}$ and ${I_{i,pre}}$ respectively represent the measured (the value obtained from measurements or MC simulations) and predicted (the value obtained from the metamodels) diffuse reflectance intensity at the ith source-detector distance.
The evaluation of the inverse metamodels was performed on the same sample sets as mentioned in section 2.7. To compare the validation/prediction performance of both inverse metamodels, the measured diffuse reflectance profiles of the 12 validation and 49 calibration phantoms at 650 nm, with the same format as the data in each metamodel training, were selected as the inputs for both the inverse measurement-based and simulation-based metamodels and the predicted BOP were compared with the reference values.
3. Results
3.1 BOP of the optical phantoms
The average BOP for each calibration (black rectangles in Fig. 2) and validation (red dots in Fig. 2) phantom in the wavelength range from 530 nm to 970 nm are presented in Fig. 5. As shown in Fig. 5(a), the μa values of the prepared samples varied between 0 cm-1 and 17.5 cm -1. The main absorption region for all the optical phantoms is located between 530 nm and 700 nm, where the maximal μa value for each sample can be observed around 620 nm. On the other hand, the μa values in the wavelength range from 700 nm to 970 nm are close to zero. The μs values range from 20 cm-1 to 400 cm-1 and the μs’ values range from 10 cm-1 to 100 cm-1, both exponentially decreasing with increasing wavelength. The g value also decreases with increasing wavelength and varies between 0.45 and 0.80. Moreover, the decreasing trend seems to be steeper for the samples with higher scattering level than for those with lower scattering level. Detailed information on the BOP for A1, A7, G1 and G7 at 650 nm is listed in Table 2.
Fig. 5. Measured values for (a) μa, (b) μs, (a) g and (b) μs’ of the liquid phantoms. 1 to 7 represent the absorption levels in increasing order and A to G represent the scattering levels in increasing order, Val refers to the validation phantoms.
Table 2. Bulk optical properties of the phantoms used for comparison
3.2 Mismatch between measurements and simulations
3.2.1 Influence from light source
The diffuse reflectance profiles for optical phantoms A1, A7, G1 and G7 at 650 nm were simulated with different input powers for the same light beam radius (results not shown). All the results indicated that the spectral intensity at the same source-detector distance is in proportion to the input power and that the input power does not change the shape of the diffuse reflectance profiles. Generally, the input power is explained by the number of photons in the measurement, while in the simulation, once the distribution of the diffusely reflected photons is stable, the input power can be conceptually transferred to the weight of the input photon packets, which is directly related to the weight of the exit photon packets, namely the intensity of the obtained diffuse reflectance profiles.
The influence from the size of the light beam (represented in the simulations by its radius) was also explored for the 4 selected phantoms. Light beams with a fixed power of 1 J and radii of respectively 0.01, 0.03 and 0.05 cm were used in the simulation. As shown in Fig. 6, samples illuminated by a light beam with a larger radius tend to have diffuse reflectance profiles with a lower intensity in the short source-detector distance and a flatter decreasing trend with the distance relative to the illumination center. Moreover, the influence from the size of the light beam diminishes with increasing source-detector distance, which is quite obvious for A1 in Fig. 6(a). For phantoms with high absorption and/or scattering levels, the diffuse reflectance quickly reduces over a short source-detector distance, as illustrated in Fig. 6(b), (c) and (d). In these cases, the effect of the light beam size remains visible until the diffuse reflectance values become very low. This indicates that the size of the light beam has an important impact on the diffuse reflectance profiles.
Fig. 6. Comparison of the measured diffuse reflectance profiles at 650 nm for optical phantoms (a) A1, (b) A7, (c) G1 and (d) G7 with the corresponding profiles simulated with radii of the illuminating beam of 0.01, 0.03 and 0.05 cm.
After applying the conditions of the light source used in the measurement (power of 1 J and beam radius of 0.03 cm at 650 nm) into the simulations, the simulated diffuse reflectance profiles (yellow line in Fig. 6) were scaled and compared with the ones extracted from the laser scattering images (blue line in Fig. 6). As shown in Fig. 6, the simulated profiles (yellow, beam radius 0.03 cm) for A1 and G7 generally match with the measured values, except for the zone close to the illumination point. Since all the samples were illuminated by the same light beam, the scaling factor was assumed to be suitable for all the samples. However, in Fig. 6(b) and (c), the simulations either overestimate or underestimate the measured values, which might be linked to a high absorption and scattering, respectively. These two effects seem to partly cancel each other out in the case of G7 with high absorption and high scattering in Fig. 6(d). Nevertheless, there is no perfect match between the simulations and measurements, thus further exploration of the other factors seems also crucial.
3.2.2 Influence from incident angle
Next, the effect of the incident angle on the mismatch between simulations and measurements was explored. As shown in Fig. 7, the diffuse reflectance profiles extracted from the left part with incident angles of 10 °, 20 ° and 30 ° fit well with the ones with 0 ° incident angle (called “standard state/value” in the following section). This indicates that the incident angle had little influence on the profiles at the left side. On the other hand, the deviation between the profiles extracted at the up/down-side and the standard values (for an illumination angle of 0 °) increased with increasing incident angle. In addition, it can be observed that under the same nonzero incident angle, the profiles extracted at the up and down sides are nearly equidistant from the ones extracted for the “standard state” with the same reflectance intensity for the high scattering levels (G1 and G7). For the low scattering levels (A1 and A7), the reflectance profiles extracted at the upside were further away from the “standard value” than those at the downside. This indicates that the illumination angle influences the diffuse reflectance profiles in different ways, depending on the BOP of the sample.
Fig. 7. Comparison of the simulated diffuse reflectance profiles for different combinations of the direction of the region of interest (left, up and down) and incident angle (0 °, 10 °, 20 ° and 30 °) for (a) A1, (b) A7, (c) G1 and (d) G7 at 650 nm. The extraction angle was fixed at 30 °.
As observed above, the influence from the incident angle on the profiles extracted from the left side was quite limited. Therefore, this region of interest was used in the further analyses. Although the solid lines in Fig. 7(c) and (d) matched quite well, small deviations can be noticed in Fig. 7(a) and (b). Thus, the data for phantoms A1 and A7 were used to further explore the influence from the incident angle. As already mentioned above, an incident angle larger than 0 ° changed the shape of the illumination spot, and thus changed the spatial distribution of the diffusely reflected light. Therefore, the combined effect of the direction of the region of interest (left, up and down) and the incident angle were explored, following to that, the area of the region of interest (30 °, 60 ° and 90 °) and the incident angle were combined to further study their influence. To make the comparison clearer, only the diffuse reflectance values within the source-detector distance of 0.05 cm were used. As shown in Fig. 8, the diffuse reflectance profiles start to deviate from the “standard values” (blue lines) for increasing incident angles. However, by increasing the area of the region of interest, the average profiles obtained from that area partially compensate the mismatch caused by the non-zero incident angle. For example, under the incident angle of 20 ° (yellow lines), which was used in the measurement configuration in this study, the diffuse reflectance profiles extracted from the left 90 ° area (yellow dotted lines) matched better with the “standard values” (blue lines) than the ones extracted from the left 30 ° area (yellow solid lines).
Fig. 8. Comparison of the simulated diffuse reflectance profiles with different combinations of the area of the region of interest (30 °, 60 ° and 90 °) and incident angle (0 °, 10 °, 20 ° and 30 °) for optical phantoms (a) A1 and (b) A7 at 650 nm.
3.2.3 Influence from detection angle
To investigate the effect of the detection angle, the diffuse reflectance profiles were collected in every detection angle segment of 1.5 °, 6 of which, with an interval of 15 ° were used for the comparison. As shown in Fig. 9, all the selected diffuse reflectance profiles, no matter at which detection angle, decreased exponentially with increasing source-detector distance. At a specific distance, the diffuse reflectance first increased with the increasing detection angle until reaching a maximum at a detection angle around 45 °. Afterwards, the intensity started to decrease again. It should be noted that the differences in the diffuse reflectance values are not linearly proportional to the difference in the detection angle. As the detection angle deviates more from the 45 ° direction, the difference in the diffuse reflectance values gradually increases. However, when the diffuse reflectance profiles were plotted on a logarithmic scale (results not shown), they were more parallel.
Fig. 9. Effect of the detection angle on the simulated diffuse reflectance profiles at 650 nm for optical phantoms (a) A1, (b) A7, (c) G1 and (d) G7.
In spatially resolved spectroscopy based on reflectance imagery (e.g., hyperspectral laser scatter imaging), the camera used to capture the diffuse reflectance images is typically installed right above the sample. This allows to easily use the obtained image without correcting for the spatial distortion. However, the range of detection angles can vary among different configurations. Therefore, different ranges of detection angles defined from the sample surface normal were used to further explore the influence from the detection angle on the mismatches. As the profiles extracted from the left 90 ° region of interest were least influenced by the incident angle in section 3.2.2, this region of interest of the measured images was used for this comparison. Similar to the comparison in section 3.2.1, the simulated results were first convoluted by applying the conditions of the light source used in the measurements, and then the convoluted profiles were divided by a scaling factor. Since the spectral intensity differed with the detection angle range, a scaling factor was calculated for each of them. In Fig. 10, the simulated diffuse reflectance profiles obtained for different detection angles are compared with the measured ones for the 4 selected optical phantoms.
Fig. 10. Relative deviation of the diffuse reflectance values simulated for different detection angle ranges (0 - 1.5 °, 0 - 6 °, 0 - 15 °, 0 - 30 °, 0 - 60 ° and 0 - 90 °) from the measured values at 650 nm for optical phantoms (a) A1, (b) A7 (c) G1 and (d) G7.
For detection angles within the cone of 30 ° relative to the sample surface normal, the relative deviations are basically within (+/-) 6%, except for the values at the distance within 0.05 cm for A7 and beyond 0.06 cm for G7. The deviations at 0.04 cm distance in A7 may be attributed to an error in locating the center of the illumination spot in the acquired diffuse reflectance image. A high absorption level reduced the diffuse reflectance profiles significantly with increasing source-detector distance, while a low scattering level reduced the overall diffuse intensity, both highlighting the distortion of the light spot caused by the incident angle in the measurement. The slightly larger deviations for G7 may be attributed to the low signal to noise ratio in the measurements, which can also be observed in Fig. 6(d). Furthermore, the relative deviations increase for larger detection angle ranges. This indicates that the detection angle range should be correctly taken into account in the simulations to avoid large mismatches between the simulated and measured diffuse reflectance profiles.
3.2.4 Improvement of the matching performance
The effect of minimizing the impact of the identified sources of mismatch between measurements and simulations is illustrated in Fig. 11. Compared to the traditional approach of radial averaging and ignoring the effect of the detection angle (Fig. 11(a)), the simulation results match quite well with the measurements after implementing these corrections (Fig. 11(b)), with R2 values above 0.99 for all the 4 samples and mean relative errors of 1.83%, 5.18%, 0.96% and 3.57% for A1, A7, G1 and G7, respectively.
Fig. 11. Comparison between measurements and the simulated diffuse reflectance profiles at 650 nm based on (a) the region of interest of the full 360 ° and detection angle range of 0 - 90 ° and (b) the region of interest of left 90 ° and detection angle range of 0 - 1.5 °. The solid lines and symbols respectively represent the measured and simulated profiles.
3.3 Metamodeling
3.3.1 Inputs of the metamodel
After obtaining the good matching performance between measurements and simulations, the measured and simulated diffuse reflectance values were used for metamodel building. For each selected source-detector distance, a metamodel was built to predict the diffuse reflectance value at that distance from the BOP µa and µs’. As mentioned above, Monte Carlo simulations provide the flexibility to use a uniform distribution of the BOP over the design space, as illustrated in Fig. 12(b). In the case of the measurements, the BOP combinations are limited to what can be designed by changing the concentrations of the absorber and scatterer in the liquid phantoms. As shown in Fig. 12(a), many data points are available in the region with low μa values, corresponding to wavelengths outside the absorption peak of the absorber. As shown in Fig. 5(a), the main absorption region of the liquid phantoms is located between 530 nm and 700 nm, while the μa values are close to zero for the remaining wavelengths. Since the μa values are combined with the μs’ values, excluding the μa values in the non-absorption region would result in removing important parts of the μs’ values. As a result, more samples would have to be prepared to obtain a wide range of BOP. In this research, all the measured BOP in the 530 - 970 nm wavelength range were used to build the measurement-based metamodel.
Fig. 12. Visualization of the data used for building the (a) measurement-based and (b) simulation-based metamodels for a source-detector distance of 0.04 cm.
3.3.2 Evaluation of the metamodel
The performance of the metamodels was evaluated by comparing the diffuse reflectance profiles obtained from the metamodel with the values measured/simulated based on the same BOP. As shown in Fig. 13, most of the data points are on the target line for both the measurement-based and simulation-based metamodel, which indicates that the diffuse reflectance values predicted by the metamodels are in good agreement with the corresponding values obtained from the measurements/simulations. The R2-values above 0.98 for all the samples indicate that the SK metamodel is well able to predict the diffuse reflectance values from the BOP. Meanwhile, small deviations can be observed for parts of the calibration samples with higher diffuse reflectance intensities (black rectangles in Fig. 13(a)) in the measurement-based metamodel and parts of the prediction samples with lower diffuse reflectance intensities (blue triangles in Fig. 13(b)) in the simulation-based metamodel. The deviated samples in Fig. 13(a) actually correspond to the data points from the 700 nm to 970 nm wavelength range, where the μa values are low, as shown in Fig. 5(a). The deviation of these data points could be mainly attributed to the variation in the diffuse reflectance values for the different data points corresponding to nearly the same BOP combinations. As the BOP values used for training the simulation-based metamodel were designed to be evenly distributed over a smooth search space, the fitting performance outperformed the one for the measurement-based metamodel. The slightly deviated prediction samples in Fig. 13(b) correspond to the values for phantom A7 at the selected source-detector distances. As can be viewed in Fig. 11(b), after minimizing the impact of the identified sources of mismatch on the diffuse reflectance profiles, the simulated values are still slightly higher than the measured values, especially at the short source-detector distances. Therefore, the slight deviation in Fig. 13(b) can be attributed to the remaining mismatch between simulations and measurements for some specific samples. In general, the diffuse reflectance profiles obtained from the metamodels matched well with the measured/simulated values.
Fig. 13. Scatterplots of the diffuse reflectance values predicted by the metamodels against the (a) measured and (b) simulated values. The red line is the 1:1 line, the dark rectangles, red dots and blue triangles respectively represent the calibration, validation and prediction samples.
3.3.3 Evaluation of the inverse metamodel
The performance of the inverse metamodels was evaluated by comparing the BOP obtained from the inverse metamodel with the corresponding reference values. As shown in Fig. 14(a1) and (b1), the calibration and validation performances of the measurement-based and simulation-based metamodels for the absorption coefficient μa were comparable, with R2 values above 0.98. However, the prediction performance of the simulation-based metamodel was slightly worse than the corresponding performance of the measurement-based metamodel. For the prediction performance (blue triangles in Fig. 14(b1)), one sample out of each absorption level obviously deviated from the target line, with predicted μa smaller than the measured values. These data points corresponded to the samples with low scattering levels (A1 to A7). If these samples were excluded, then the prediction performance of the simulation-based metamodel would be comparable with that of the measurement-based metamodel. As already explained in section 3.3.2, these deviations can be attributed to the remaining mismatch between simulations and measurements for those samples. As can be observed in Fig. 13(b), the predicted diffuse reflectance profiles for these samples were slightly higher than the measured values in the forward metamodel. In other words, the measured diffuse reflectance profiles, as the input of the inverse metamodel, were lower than the desired values. As a consequence, the predicted BOP combination shifted towards lower μa and lower μs’. Moreover, as the resolution of μa values was higher than that of μs’ values in the metamodel, the deviation in the μa values was more obvious. For the scattering coefficient μs’, the calibration performance of the measurement-based metamodel was slightly worse than for the simulation-based metamodel, which can be attributed to the deviating data points in Fig. 13(a). The validation and prediction performances both for the measurement-based and the simulation-based inverse metamodel were quite good, with R2 values above 0.98. In general, the simulation-based metamodel obtained comparable prediction performance on measurements as the measurement-based metamodel, with mean relative errors of 8.19% (validation) and 10.96% (prediction) for μa and 5.12% (validation) and 6.42% (prediction) for μs’.
Fig. 14. Scatterplot of (1) μa and (2) μs’ predicted by the inverse metamodel against the reference values for the (a) measurement-based and (b) simulated-based inverse metamodels. The red line is the 1:1 line, the dark rectangles, red dots and blue triangles respectively represent the calibration, validation and prediction samples.
4. Discussion
In this research, the mismatch between simulated and measured diffuse reflectance profiles was explored for the three key elements of a measurement configuration: light source, incident angle and detection angle. For a configuration composed of optical fiber arrays, the incident angle is perpendicular to the sample surface and the detection angle covers photon exit angles that bare much larger than an imagery-based configuration. Accordingly, the fiber-based systems are more consistent with the default settings in the MC simulations considering photon exit angles from 0 to 90 °. For such a configuration, once the conditions of the light source were applied in the simulations, the simulated diffuse reflectance profiles fitted relatively well with the measured values. This can explain the good prediction performance of the simulation-based model on the measurements in the research of Hennessy et al. [20]. In the case of a reflectance imagery-based configuration which is more suitable for contactless, fast and non-destructive BOP measurement, the incident and detection angle have a large impact on the acquired signals. However, these effects are ignored in most studies. For instance, in the research of Zhang et al. [23], the CCD camera was installed at an angle of 8 ° with respect to the sample surface normal, while the illumination beam coincided with the sample surface normal. This angle of 8 ° was neglected in the simulations, given the reason that the image was symmetric to the incident point as this angle would not cause a significant difference in the diffuse reflectance profiles. As the positions of the camera and the illumination light were exchanged compared to the ones in the present research, this angle can be regarded as the incident angle. According to the results obtained in the present research, such an incident angle would indeed change the radial symmetry of the diffusely reflected light and its influence depends on the BOP of the object as well as the selected source-detector distances. For the samples with high scattering levels, the incident angle made the spectral intensity at the sides towards/away from the illumination fiber deviating from the expected value with a similar extent in opposite directions. Radial averaging thus balanced the deviations. However, for the samples with low scattering levels, the spectral intensity at those two sides still deviated in opposite directions but with different extents, resulting in a mismatch between the radial averaged spectral intensity and the desired values. This may partly explain the poor prediction performance of their simulation-based model on the measurements (relative error increased from 2.5% and 2% to 50% and 25%).
Compared to the BOP prediction performances obtained by other researchers, some improvements were made in this research. Qin et al. [7] used a trust-region nonlinear least-squares fitting algorithm to fit the measured diffuse reflectance profiles for a set of liquid phantoms with the diffusion equation to estimate the optical properties. The average fitting errors in μa and μs’ were respectively 12% and 7%, while for the samples with relatively high μa values (0.4 - 0.8 cm-1), the fitting errors could be up to 17% and 11%. In the present research, not only the μa values increased 20 times, which could be up to around 16 cm-1, but at the same time, the fitting errors were slightly reduced to 10.98% and 6.42%. Similarly, the ranges in μa and μs’ values were both extended compared to the work executed by Aernouts et al. [11] to better match the relevant range for biological products characterized by high absorption levels. Moreover, the simulation-based metamodel presented here provides a more flexible option for the determination of BOP.
In the research of Palmer et al. [25], a model with wide ranges of BOP was built based on Monte Carlo simulations. The validation performance of that model was evaluated on optical phantoms, obtaining average fitting errors within 12%, which were comparable with the results obtained in the present research. However, they measured the diffuse reflectance profiles of the phantoms with a fiber optic probe.
When comparing the measurement-based and simulation-based metamodel, several advantages can be observed for the latter one. The prediction performance of the SK metamodels largely depends on the design of the inputs during the training process. For instance, the value of an unknown point in a SK model is determined by its distances to the surrounding points, while the prediction performance would be worse for the points at the edge of the search space. To avoid this problem, a large number of samples with a wide range of optical properties need to be prepared to cover the desired search space for the training of the measurement-based metamodel. Although the simulation-based metamodel faces the same limitation, the search space can be freely extended to be sufficient to cover the desired space. Moreover, as shown in Fig. 12(a) and Fig. 13(a), the overrepresentation of the low μa values had a negative impact on the quality of the measurement-based metamodel. Therefore, the data points in that region need to be carefully selected to get a more even distribution and without losing regions of the search space. In the case of the simulation-based metamodel, the distribution of the data points can be arbitrarily designed, which is far more convenient.
Another advantage of the simulation-based metamodel is the universality. The inputs (a wide range of BOP) and the outputs (diffuse reflectance profiles as a function of source-detector distance and detection angle) obtained from the simulations in the present research can be regarded as the original database. For the future application on measurements, the light source condition, the resolution in the source-detector distance direction and the detection angle can be easily updated by using the convolution program (CONV) developed by Wang et al. [36]. The difference in the diffuse reflectance intensity caused by the power difference among wavelengths and the signal transformation between measurements and simulations, can be corrected with a scaling factor that is obtained on a small set of samples (e.g., 4 samples used in this research). In this way, the built metamodel can be easily transferred and applied to different measurement configurations.
However, several recommendations can be made towards the further application of the simulation-based metamodel. First, the use of an incident angle of 0° in the simulations is inconsistent with the actual value in the measurements. As a consequence, a diffuse reflectance mismatch between measurements and simulations will be present for some specific samples (e.g., A7 in this research). Decreasing the incident angle in the measurement configuration will help to reduce this mismatch. Second, the detection angle used in this research was calculated based on the central point right below the camera. Since the distance between the camera and the sample was far greater than the size of the camera and because the latter is larger than the source-detector distance of interest, the asymmetric exit angle that can be detected for points away from the center could be neglected. However, for a configuration with opposite conditions, the change in the captured exit angle with the source-detector distance should be considered as well. Otherwise, it may introduce an extra mismatch between the simulations and the measurements. Third, the SK model allows to set the intrinsic uncertainty manually. For the simulation-based metamodel, especially for application on the measurements, the intrinsic uncertainty should be set to a similar level as the measurement noise. Moreover, the remaining relative deviations between measurements and simulations after applying the correction procedures (e.g., 6% in this research) should also be considered in the metamodel training.
5. Conclusions
The effects of the power, size and incident angle of the illumination beam and the detection angle on the diffuse reflectance profiles were investigated with Monte Carlo simulations. The power of the light source was found to determine the intensity level of the diffuse reflectance profiles, while the size of the light beam influenced the decreasing trend of the profiles with increasing source-detector distance, which was most obvious in the region with high signals near the illumination center. Deviation of the incident angle from the normal of the sample surface changed the central symmetry of the light spot, thereby influencing the spatial distribution of the diffusely reflected light. The diffuse reflectance profiles extracted from the sides perpendicular to the illuminating beam were close to the profiles extracted with 0 ° incident angle. Through optimization of the area of the region of interest, the mismatch between simulations and measurements could be further reduced, especially for the samples with low scattering level. At any detection angle within 0 ° - 30 °, the relative deviations between measured and simulated diffuse reflectance values remained below (+/-) 6%. A good match (R2 above 0.99) between the simulated and measured diffuse reflectance profiles over source-detector distances from 0.04 cm to 0.07 cm was achieved after accounting for the detection angle of 0 ° - 1.5 ° and with the profiles extracted from the measured images at the left 90 °. After minimizing the impact of the considered sources of mismatch on the diffuse reflectance profiles, the stochastic Kriging metamodels linking the optical properties with the diffuse reflectance profiles were respectively built for the measurements and the simulations. Small deviations in the predicted diffuse reflectance profiles were found in the wavelength range with extremely low absorption levels in the measurement-based metamodel. Moreover, the fitting performance of the simulations-based metamodel outperformed the one for the measurement-based metamodel thanks to its even distribution of the BOP over the search space. Except for a few specific samples with remaining mismatch between measurements and simulations, the simulation-based inverse metamodel obtained an equally good BOP prediction performance as the one based on measurements, obtaining R2 values above 0.98 both for the μa and μs’.
Funding
KU Leuven (C1 project -3E160393); China Scholarship Council (201606320240).
Disclosures
The authors declare no conflicts of interest.
References
1. S. H. E. J. Gabriëls, P. Mishra, M. G. J. Mensink, P. Spoelstra, and E. J. Woltering, “Non-destructive measurement of internal browning in mangoes using visible and near-infrared spectroscopy supported by artificial neural network analysis,” Postharvest Biol. Technol. 166, 111206 (2020). [CrossRef]
2. W. Lan, B. Jaillais, A. Leca, C. M. G. C. Renard, and S. Bureau, “A new application of NIR spectroscopy to describe and predict purees quality from the non-destructive apple measurements,” Food Chem. 310, 125944 (2020). [CrossRef]
3. P. Theanjumpol, K. Wongzeewasakun, N. Muenmanee, S. Wongsaipun, C. Krongchai, V. Changrue, D. Boonyakiat, and S. Kittiwachana, “Non-destructive identification and estimation of granulation in ‘Sai Num Pung’ tangerine fruit using near infrared spectroscopy and chemometrics,” Postharvest Biol. Technol. 153, 13–20 (2019). [CrossRef]
4. K. Włodarska, J. Szulc, I. Khmelinskii, and E. Sikorska, “Non-destructive determination of strawberry fruit and juice quality parameters using ultraviolet, visible, and near-infrared spectroscopy,” J. Sci. Food Agric. 99(13), 5953–5961 (2019). [CrossRef]
5. R. Lu, R. Van Beers, W. Saeys, C. Li, and H. Cen, “Measurement of optical properties of fruits and vegetables: A review,” Postharvest Biol. Technol. 159, 111003 (2020). [CrossRef]
6. C. Sun, R. Van Beers, B. Aernouts, and W. Saeys, “Bulk optical properties of citrus tissues and the relationship with quality properties,” Postharvest Biol. Technol. 163, 111127 (2020). [CrossRef]
7. J. Qin and R. Lu, “Hyperspectral diffuse reflectance imaging for rapid, noncontact measurement of the optical properties of turbid materials,” Appl. Opt. 45(32), 8366–8373 (2006). [CrossRef]
8. R. Van Beers, B. Aernouts, R. Watté, A. Schenk, B. Nicolaï, and W. Saeys, “Effect of maturation on the bulk optical properties of apple skin and cortex in the 500-1850nm wavelength range,” J. Food Eng. 214, 79–89 (2017). [CrossRef]
9. D. Milej, A. Abdalmalak, D. Janusek, M. Diop, A. Liebert, and K. S. Lawrence, “Time-resolved subtraction method for measuring optical properties of turbid media,” Appl. Opt. 55(7), 1507–1513 (2016). [CrossRef]
10. D. Hu, R. Lu, and Y. Ying, “Spatial-frequency domain imaging coupled with frequency optimization for estimating optical properties of two-layered food and agricultural products,” J. Food Eng. 277, 109909 (2020). [CrossRef]
11. B. Aernouts, C. Erkinbaev, R. Watté, R. Van Beers, N. Do Trong, B. Nicolai, and W. Saeys, “Estimation of bulk optical properties of turbid media from hyperspectral scatter imaging measurements: metamodeling approach,” Opt. Express 23(20), 26049–26063 (2015). [CrossRef]
12. P. Lemaillet, J. P. Bouchard, J. Hwang, and D. W. Allen, “Double-integrating-sphere system at the National Institute of Standards and Technology in support of measurement standards for the determination of optical properties of tissue-mimicking phantoms,” J. Biomed. Opt. 20(12), 121310 (2015). [CrossRef]
13. A. López-Maestresalas, B. Aernouts, R. Van Beers, S. Arazuri, C. Jarén, J. De Baerdemaeker, and W. Saeys, “Bulk optical properties of potato flesh in the 500-1900nm range,” Food Bioprocess Technol. 9(3), 463–470 (2016). [CrossRef]
14. B. Aernouts, R. Van Beers, R. Watté, T. Huybrechts, J. Lammertyn, and W. Saeys, “Visible and near-infrared bulk optical properties of raw milk,” J. Dairy Sci. 98(10), 6727–6738 (2015). [CrossRef]
15. A. Kienle, L. Lilge, M. S. Patterson, R. Hibst, R. Steiner, and B. C. Wilson, “Spatially resolved absolute diffuse reflectance measurements for non-invasive determination of the optical scattering and absorption coefficients of biological tissue,” Appl. Opt. 35(13), 2304–2314 (1996). [CrossRef]
16. M. Pilz, S. Honold, and A. Kienle, “Determination of the optical properties of turbid media by measurements of the spatially resolved reflectance considering the point-spread function of the camera system,” J. Biomed. Opt. 13(5), 054047 (2008). [CrossRef]
17. T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties invivo,” Med. Phys. 19(4), 879–888 (1992). [CrossRef]
18. A. H. Hielscher, S. L. Jacques, L. Wang, and F. K. Tittel, “The influence of boundary conditions on the accuracy of diffusion theory in time-resolved reflectance spectroscopy of biological tissue,” Phys. Med. Biol. 40(11), 1957–1975 (1995). [CrossRef]
19. N. Rajaram, T. H. Nguyen, and J. W. Tunnella, “Lookup table–based inverse model for determining optical properties of turbid media,” J. Biomed. Opt. 13(5), 050501 (2008). [CrossRef]
20. R. Hennessy, S. L. Lim, M. K. Markey, and J. W. Tunnell, “Monte Carlo lookup table-based inverse model for extracting optical properties from tissue-simulating phantoms using diffuse reflectance spectroscopy,” J. Biomed. Opt. 18(3), 037003 (2013). [CrossRef]
21. R. Watté, N. Nguyen, B. Aernouts, C. Erkinbaev, J. De Baerdemaeker, B. Nicolaï, and W. Saeys, “Metamodeling approach for efficient estimation of optical properties of turbid media from spatially resolved diffuse reflectance measurements,” Opt. Express 21(26), 32630–32642 (2013). [CrossRef]
22. X. Zhong, X. Wen, and D. Zhu, “Lookup-table-based inverse model for human skin reflectance spectroscopy: two-layered Monte Carlo simulations and experiments,” Opt. Express 22(2), 1852–1864 (2014). [CrossRef]
23. L. Zhang, Z. Wang, and M. Zhou, “Determination of the optical coefficients of biological tissue by neural network,” J. Mod. Opt. 57(13), 1163–1170 (2010). [CrossRef]
24. M. O. Gökkan and M. Engin, “Artificial neural networks-based estimation of optical parameters by diffuse reflectance imaging under in vitro conditions,” J. Innovative Opt. Health Sci. 10(01), 1650027 (2017). [CrossRef]
25. G. M. Palmer and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms,” Appl. Opt. 45(5), 1062–1071 (2006). [CrossRef]
26. B. Aernouts, E. Zamora-Rojas, R. Van Beers, R. Watté, L. Wang, M. Tsuta, J. Lammertyn, and W. Saeys, “Supercontinuum laser based optical characterization of Intralipid® phantoms in the 500-2250 nm range,” Opt. Express 21(26), 32450–32467 (2013). [CrossRef]
27. R. Van Beers, B. Aernouts, L. León Gutiérrez, C. Erkinbaev, K. Rutten, A. Schenk, B. Nicolaï, and W. Saeys, “Optimal illumination-detection distance and detector size for predicting braeburn apple maturity from Vis/NIR laser reflectance measurements,” Food Bioprocess Technol. 8(10), 2123–2136 (2015). [CrossRef]
28. S. A. Prahl, M. J. C. Van Gemert, and A. J. Welch, “Determining the optical properties of turbid media by using the adding-doubling method,” Appl. Opt. 32(4), 559–568 (1993). [CrossRef]
29. T. Li, H. Gong, and Q. Luo, “MCVM: Monte carlo modeling of photon migration in voxelized media,” J. Innovative Opt. Health Sci. 03(02), 91–102 (2010). [CrossRef]
30. E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt. 13(6), 060504 (2008). [CrossRef]
31. R. Watté, B. Aernouts, R. Van Beers, E. Herremans, Q. T. Ho, P. Verboven, B. Nicolai, and W. Saeys, “Modeling the propagation of light in realistic tissue structures with MMC-fpf: a Meshed Monte Carlo method with free phase function,” Opt. Express 23(13), 17467–17486 (2015). [CrossRef]
32. S. Chugunov and C. Li, “Monte carlo simulation of light propagation in healthy and diseased onion bulbs with multiple layers,” Comput. Electron. Agric. 117, 91–101 (2015). [CrossRef]
33. F. Vaudelle, “Light scattering in thin turbid tissue including macroscopic porosities: A study based on a Monte Carlo model,” Opt. Commun. 425, 91–100 (2018). [CrossRef]
34. C. Sun, B. Aernouts, R. Van Beers, and W. Saeys, “Simulation of light propagation in citrus fruit using monte carlo multi-layered (MCML) method,” J. Food Eng. 291, 110225 (2021). [CrossRef]
35. S. Kedenburg, M. Vieweg, T. Gissibl, and H. Giessen, “Linear refractive index and absorption measurements of nonlinear optical liquids in the visible and near-infrared spectral region,” Opt. Mater. Express 2(11), 1588–1611 (2012). [CrossRef]
36. L. Wang, S. Jacques, and L. Zheng, “MCML-Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef]
37. X. Ruan, Q. Zhou, L. Shu, J. Hu, and L. Cao, “Accurate prediction of the weld bead characteristic in laser keyhole welding based on the stochastic kriging model,” Metals 8(7), 486 (2018). [CrossRef]
38. B. Wang and J. Hu, “Some monotonicity results for stochastic kriging metamodels in sequential settings,” INFORMS J. Comput. 30(2), 278–294 (2018). [CrossRef]
39. X. Chen and Q. Zhou, “Sequential design strategies for mean response surface metamodeling via stochastic kriging with adaptive exploration and exploitation,” Eur. J. Oper. Res. 262(2), 575–585 (2017). [CrossRef]
