1 Introduction

The propagation of light through biological tissues is an important research topic in the domain of biomedical optics, serving as a basis for non-destructive in-vivo medical diagnostics e.g. cancer detection, etc [1]. With recent advances in optical sensing and the technology becoming more cost-efficient, these methods draw the attention of food and agricultural industry for quality assessment of fruits and vegetables, determining moisture content in food, fat content in meat, etc [2–5 ]. Moreover, the development of innovative optical sensors for rapid, non-destructive and non-contact food quality monitoring would be promoted once the bulk optical properties (BOP) of such agri-food products can be accurately measured on the production or sorting line.

When light travels through the sample tissue, two phenomena can occur: absorption of light by the tissue’s molecules, and light scattering, caused by refractive index mismatches inside the sample. Light absorption is related to the chemical constituents (e.g., water, protein, fat, carbohydrates), whereas light scattering is associated with the physical characteristics of the tissue (e.g., cellular structure, turgor, density) [1]. As these properties define light propagation in the biological tissue, and accordingly the measured visible (Vis) and near-infrared (NIR) signals, this could open opportunities for simultaneous extraction of chemical (composition) and physical (microstructure) information from optical measurements. However, in a single Vis/NIR spectroscopic measurement, usually reflectance or transmittance, both the effect of absorption and scattering are indissolubly connected and cannot be accurately separated. Consequently, a change in the scattering properties of the measured sample might be misinterpreted as a change in the sample composition and vice versa. On the other hand, multiple spectroscopic measurements in a slightly different measurement configuration are influenced by absorption and scattering in a different way. The combination of these multiple measurement series with an accurate model, which mathematically describes the Vis/NIR light propagation as a function of the sample’s absorption and scattering properties, could allow for separation of the BOP. These multiple spectroscopic measurements can vary in space [6,7 ], time [8–10 ] and/or frequency [11,12 ]. Amongst them, spatially resolved spectroscopy (SRS) is fast, cost-efficient and simple compared to the others. Moreover SRS can be easily adapted for routine laboratory measurements [5].

Two types of SRS methods have already been reported: contact [6,7,13 ] and non-contact measurements [14–16 ]. A contact SRS method uses a fiber probe comprised of one or multiple light illumination fibers and several collection fibers. As the probe is positioned in contact with the sample, this technique is less applicable for on-line industrial applications, mainly due to the risk of cross-contamination and the additional time needed for positioning the probe on the sample. Alternatively, non-contact SRS, which is also known as hyperspectral scatter imaging (HSI), uses a hyperspectral line-scan camera to capture the spatially resolved reflectance profile at the position where the sample is illuminated with a focused white light beam [3,4,17 ]. Thanks to the contactless character of both the illumination and the imaging, this technique has more potential for integration in sorting lines (e.g. fruits, vegetables etc.) in an industrial environment. Therefore, the potential of this SRS technique for contactless food quality control has recently become a relevant research topic [2,3,17 ].

In order to estimate the BOP from the acquired optical measurements, the relation between both needs to be resolved. Light propagation in turbid media can be accurately described by the radiative transfer theory (RTT) which neglects the wave properties and only considers the energy transport [18]. However, due to the geometrical complexity of the microstructure and the phenomenon of multiple light scattering in real biological media, an analytical solution of the RTT is practically impossible. When the radiance in a homogeneous medium is considered at a sufficient distance from the illumination point, where the assumption of isotropic radiance holds, and when scattering is dominant over absorption, the RTT can be approximated with the diffusion theory [19–21 ]. However, in the NIR, this last condition is rarely fulfilled for biological tissues as the absorption by the water molecules generally dominates the scattering. For these samples, the RTT can be solved numerically based on stochastic Monte Carlo (MC) simulations. This method has been shown to provide high accuracy for many different sample geometries and any range of absorption and scattering properties [22,23 ]. However, to reduce stochastic noise, the Monte Carlo method requires to simulate the propagation of many photons, which makes this method very time-consuming. The simulation speed can be further increased by using modern computers with high computational power, employing parallel computing approaches [24–26 ] and/or rescaling originally-simulated photons for each new set of BOP [27–31 ]. Nevertheless, despite of these improvements, the required time to accurately simulate the light propagation for a single set of BOP is still too high. Moreover, depending on the BOP, and with that the number of simulated photons required for an accurate result, the minimum time needed for a single forward simulation is still in the order of seconds to minutes [29,31 ]. The inverse problem, which links a set of BOP to an optical measurement, is typically ill-posed, leading to non-uniqueness and possible instability. Therefore, the inverse problem is generally solved with an iteration over the forward model, independent of the model type (e.g. MC method, diffusion theory, etc.). Generally, the nonlinear optimizer needs to evaluate many sets of BOP before convergence is reached. As a result, the inverse MC method is very time-consuming, taking several minutes for the estimation of each BOP set. Therefore it is not suitable for fast and on-line practical applications [27,28,31 ]. Additionally, following these theoretical light propagation models (e.g. MC methods, diffusion theory, etc.), the geometric relationships between the detection system (e.g. imaging or probe) and the sample still needs to be implemented separately. Generally, this relation is not straightforward and approximations needs to be introduced [17].

To overcome the limitations of the above mentioned theoretical solutions, the use of metamodels has been proposed to directly link a set of measured SRS profiles to their respective BOP [13]. Moreover, the metamodels serve as a high-way bridge between the design space (input parameters: BOP) and the performance space (output parameters: SRS profiles). Following this approach, the estimation procedure is seriously simplified and more robust, as both the light propagation characteristic as well as the measurement geometry is correctly accounted for in a single forward model. Dam et al. [32] used polynomial regression to relate SRS profiles measured with a contact probe to their respective BOP. Nevertheless, due to the simplicity of the models, their methodology was only accurate for a very small range of absorption and scattering properties, which translates into a small wavelength-window for samples with limited variability. Other researchers used neural network, trained on the SRS profiles simulated with the MC method, to avoid the time-intensive MC simulations in future calculations [27,33 ]. The major drawback of these computationally fast neural networks is their incapability of correctly handling the stochastic noise present in Monte Carlo results or the instrumental noise typical for SRS measurement. Moreover, as the noise is incorporated into the neural network, it will influence the simulated output and consequently the inverse estimation process. In contrast, stochastic data-based surrogate models can handle noise which is typical for measurements or data simulated with stochastic techniques. Moreover, they can solve non-linear, intensive design problems with minimal computational effort, resulting into an efficient model which enables tasks such as visualization, design space exploration, sensitivity analysis and optimization [34–38 ]. Therefore, this metamodeling approach can be used as a computationally efficient method to estimate the BOP from SRS profiles [13]. Moreover, in a previous study, this approach has been successfully tested to link the BOP of a designed set of optical phantoms to their respective SRS profiles, measured with a reflectance probe which was in contact with the samples [13].

The main objective of the current study is to evaluate the performance of the surrogate metamodels if trained on hyperspectral scatter images. Due to the fundamentally different measurement geometry (e.g. collection angle, etc.) compared to a contact probe, the SRS profiles generated with this contactless technique are significantly different [17]. The metamodels should take this into account, without resorting to approximations generally used. The novelty of this research can be found in the use of (1) a contactless technique to measure the SRS profiles, in combination with (2) stochastic data-based surrogate metamodels which were trained on (3) a designed set of optical phantoms, to finally derive the BOP of turbid samples. To the best of our knowledge, this research presents the first time a hyperspectral imaging technique is used to estimate the BOP without resorting to the diffusion theory [2–4,17 ].

The developed metamodels were evaluated with respect to the prediction of SRS profiles from provided BOP (forward validation), as well as its performance in combination with an iterative optimization routine (inverse validation). The latter allows for the estimation of the sample’s BOP from contactless HSI measurements. Double integrating sphere (DIS) and unscattered transmittance measurements, in combination with the inverse adding-doubling (IAD) routine, were used to obtain the reference BOP needed for model construction and validation.

2 Materials and methods

2.1 Liquid phantoms

Liquid optical phantoms were prepared as aqueous mixtures of a scattering component: Intralipid 20% (Fresenius Kabi, Sweden), and a water-soluble absorbing dye: methylene blue (Sigma Aldrich, Belgium). Methylene blue (MB) is a non-scattering dye which allows for independent design of the scattering and absorption properties. Moreover, MB has a clear and sharp absorption peak in the considered wavelength range, making it easier to interpret the effects of scattering and absorption on the measured diffuse reflectance spectra.

In total, sixteen liquid phantoms were prepared by mixing four concentrations of the absorber (0, 0.74, 1.48 and 2.22 ml of a 400 µM stock solution of MB) with four concentrations of the scatterer (2, 4, 6 and 8 ml of the Intralipd 20% solution). The mixtures were diluted with distilled water to obtain phantoms of 100 ml. Sixteen prepared phantoms were labeled with a letter and a number, in which an increasing letter (A, B, C and D) indicates an increasing scattering level, while an increasing number (1, 2, 3 and 4) corresponds to an increasing absorption level. The phantoms were thoroughly mixed before taking an aliquot of the sample for measuring in the DIS and unscattered transmittance measurement system. The phantoms were carefully mixed and allowed to rest before scanning by the HSI system.

2.2 Reference bulk optical properties of liquid phantoms

A pre-dispersive double integrating sphere (DIS) and unscattered transmittance setup was used to acquire the diffuse reflectance and total and unscattered transmittance of the sample slabs. These measurements are considered to be the ‘golden standard’ to estimate the BOP for turbid media [39–41 ]. The sample was illuminated with a wavelength-tunable laser setup, especially designed to obtain high signal-to-noise spectra in the 550 – 2250 nm wavelength range [39,42 ]. Diffuse reflectance and total transmittance were measured simultaneously with detectors mounted on the wall of two integrating spheres, positioned respectively in front and behind the cuvette which was filled with the sample. An integrating spheres is a hollow spherical cavity with its interior covered with a diffuse reflective coating to obtain the total power coming from the sample in a single measurement, independent of the original direction and spatial information of the light. The unscattered transmittance, on the other hand, was measured in a separate path with the detector positioned 1.5 m behind the sample and a series of slits placed in between to limit the fraction of scattered photons collected [41,43 ]. For a more extensive description of the measurement setup, the calibration and measurement procedure and a thorough validation, the reader is referred to Aernouts et al. [39]. Moreover, this extensive validation study illustrates the high repeatability and signal-to-noise ratio (SNR) of the system to obtain the Vis/NIR BOP for very turbid samples [39].

The liquid phantom samples were carefully pipetted into a custom made cuvette with borosilicate glass walls of 1.1 mm thickness (Borofloat 33®, Schott, Germany) with 1 mm path length and the cuvette was placed into the sample holder. Measurements were performed only in the 550 – 950 nm spectral range with an interval step of 5 nm, because the HSI setup is limited to this range. Each liquid phantom sample was measured in three random ordered replicates and the cuvette was thoroughly cleaned after each measurement to avoid associated errors.

The BOP (µ_a and µ_s’) were derived from the diffuse reflectance and total and unscattered transmittance measurements, following the inverse adding-doubling (IAD) routine [41,44 ]. In the adding-doubling (AD) method, the RTT is used to calculate the total reflectance and transmittance for a single ‘infinitesimally’ thin sample layer for which the single scattering assumption is fulfilled. This plane-parallel layer is then ‘doubled’ and the reflectance and transmittance of the doubled layer are calculated. This process of doubling is repeated until the desired thickness of the homogeneous sample is reached. Different layers (e.g., glass slides, etc.) can be ‘added’, together with their contributions, taking into account internal reflections at the boundaries [42]. The AD method allows to calculate the diffuse reflectance and total and unscattered transmittance very accurately, while maintaining a high degree of flexibility and time efficiency [44]. However, as the forward method calculates the reflectance and transmittance for a tissue layer with known bulk optical properties, it has to be inverted to allow extraction of the bulk absorption (μ_a) and reduced scattering coefficients (μ_s’) from the measured spectra. Moreover, these BOP, which are the inputs for the AD method, are iteratively changed until the simulated reflectance and transmittance values match with the measured ones [42,44 ]. As the method needs to compensate for the specular reflectance at the sample-cuvette and cuvette-air boundaries, the wavelength-dependent refractive indices of the cuvette glass and sample are required. Refractive index information for the cuvette windows was provided by the manufacturer (Schott, Germany), while the refractive index of the intralipid phantoms was calculated with the formula for lipid particles in water [45]. The reference BOP for the three replicates of each phantom, acquired with IAD from DIS and unscattered transmittance measurements, were averaged and used as input for the metamodel and for the validation.

2.3 Hyperspectral scatter imaging

A HSI setup was developed for contactless acquisition of spatially resolved reflectance profiles in the wavelength range from 400 to 950 nm. The system has four main units: a hyperspectral camera system, a sample presentation unit, a light source and a computer.

The hyperspectral camera system consists of a high performance 12-bit panchromatic CCD camera (TXG14NIR, Baumer, Germany), a prism-grating-prism based imaging spectrograph (ImSpector V10, Spectral Imaging Ltd., Oulu, Finland) and a 12.5 – 75 mm focusing lens with f/1.8 (Cosmicar, Pentax, Japan). The camera has an optical sensitivity from 400 to 950 nm, 1392 by 1040 effective pixel image resolution and is based on GigE Vision® fast Ethernet data transfer. The imaging spectrograph covers the spectral range from 400 to 1000 nm with a resolution of 9 nm, which is well-compatible with the sensitivity of the camera.

The light source unit consists of a quartz tungsten halogen lamp (Alphabright, East Sussex, UK) generating broadband light (300 – 1550 nm) with 20 W electrical power and 280 lumen flux. The light is coupled into a 200 μm optical fiber (FC-IR200-1.5, Avantes, Eerbeek, The Netherlands) ending on a 74-VIS collimating lens (Ocean Optics, Duiven, The Netherlands) with 350 – 2000 nm range. The lens focuses the light onto the sample surface as a spot of 1.0 mm diameter at 7.5 cm distance, under an angle of 15° with respect to the vertical axis to avoid measuring the specular reflection of light. This angle was chosen as small as possible, based on the design of the HSI system. Moreover, it was confirmed that the influence of this angle on the symmetry of the captured images is negligible [46].

The sample presentation consists of an X-Y-Z axis metric stage (Edmund Industrial Optics, NJ, USA) mounted on a breadboard (Thorlabs, NJ, USA). This unit allows to precisely position the sample during each measurement to ensure an accurate height (13.5 cm) distance between sample surface and the camera lens of the measurement system. A schematic illustration of the entire HSI setup is presented in Fig. 1 . The hyperspectral imaging system was placed in a dark chamber to exclude the effects of ambient light during the measurements.

 

Fig. 1 Schematic illustration of the hyperspectral scatter imaging system.

Download Full Size | PDF

The HSI system was spectrally and spatially calibrated as described in [4]. Each phantom was measured ones. For each pixel in the scanned line, the light diffusively reflected by the sample is focused by the zoom lens onto the slit (80 µm width) of the spectrograph which disperses it into its different wavelength components. This light is further projected onto a pixel line of the CCD detector. The hyperspectral image of 696 × 520 pixels has a spectral resolution of 1.15 nm/pixel and a spatial resolution of 0.052 mm/pixel. The scanning line was taken 1.5 mm from the center [Fig. 1 top view] of the illumination spot to avoid saturation of the CCD detector pixels due to high signal intensity close to the illumination point. In order to limit the calculation time and to improve the signal-to-noise ratio, the images were binned by reducing the number of wavelength bands to 81 (from 550 – 950 nm) with a step of 5 nm, corresponding to the DIS and unscattered transmittance measurements. The region of interest (ROI) of the image was then selected based on pre-defined area in the image order to reduce the data size and speed up the computation time. The symmetric two sides in positive and negative direction of the scanning line [Fig. 1 top view] were averaged and converted to relative reflectance using the white and dark reference images, measured on a calibrated 99% reflectance Spectralon® disk (Labsphere Inc., North Sutton, USA) and a black cap, respectively [4]. At each single wavelength, the final SRS profile describes the diffuse reflectance in function of the source-detector distance. The latter varied from 0 to 0.8 cm with a 0.04 cm step, resulting in 21 source-detector distances. These SRS profiles were used as input for training the metamodel (calibration phantoms) and for validation purposes.

2.4 Metamodeling

Metamodeling, also known as surrogate modeling, is defined as the construction of an input/output function based on experimental measurements or simulations. A metamodel based on Kriging regression uses an interpolation method which computes each new value through a Gaussian process. Although the interpolation property of this method has certain advantages, it is not desirable when dealing with measurements subject to noise, which can be characterized by performing several (independent) measurements. In stochastic Kriging, on the other hand, this variance is taken into account to avoid overfitting the model [13]. Therefore, a stochastic Kriging surrogate model was constructed in this study to link the measured SRS values and the respective BOP. This was done by using the ooDACE toolbox, which is a versatile Matlab toolbox that implements the Gaussian Process based Kriging surrogate models [47–49 ]. First, the metamodels were trained on an accurate data set of SRS profiles, covering the relevant range of BOP. Such a data set can be acquired by measuring reference materials with pre-defined BOP, known as optical phantoms [1,39,45 ].

As the metamodels describe the SRS profile as a function of the BOP, the procedure should be inverted in order to estimate the BOP from an SRS profile. Once the metamodels were trained, the inverse problem is solved through an iterative optimization algorithm which uses the forward model as a basis [13]. Moreover, it uses a Nelder-Mead simplex optimization procedure, which adjusts the BOP estimates at a certain wavelength until the simulated SRS profile matches with the measured one. To avoid that the high diffuse reflectance values acquired close to the point of illumination would dominate the BOP estimation, the sum of the squared relative errors was used as the cost function in the optimization procedure [13]. As this optimization problem may not be convex, several combinations of BOP, evenly distributed in the design space were used as starting points for the simplex algorithm. This ensures that at least one of the starting points would lead to the global minimum corresponding to the actual BOP.

2.5 Forward and inverse validation of metamodeling approach using liquid phantoms

The sixteen liquid phantoms were divided into two groups: a calibration set of 12 phantoms and a validation set of 4 phantoms. The metamodels were built on the calibration set, while their performance was evaluated on the phantoms which were not used for calibration. The selection of phantoms for the calibration and validation sets is represented in Table 1 .

Table 1. Selection of calibration set (normal) and validation set (bold) of liquid phantoms.

View Table

In total, 21 metamodels were built, each one corresponding with a single source-detector distance. Every single metamodel relates the diffuse reflectance values at a specific source-detector distance with the respective set of BOP. The metamodels were constructed by linking the measured SRS profiles (output) of the calibration phantoms with the corresponding BOP (input) in the 550 – 950 nm wavelength range. First, the metamodels themselves were evaluated by comparing the diffuse reflectance values predicted by the metamodels with those measured by the HSI setup. Moreover, the root mean squared error was calculated for the calibration samples (RMSE_C) and the validation samples (RMSE_V).

After the forward validation, the performance of the inverse model for estimating the BOP of the validation samples was also investigated. This evaluation step was done by comparing the estimated BOP at all wavelengths for the four validation phantoms (A3, B1, C2 and D3) with the corresponding values obtained from DIS measurements. The prediction performance of the inverse model was then quantified in terms of the root mean squared error of prediction (RMSE_P), which was calculated as:

where:

3 Results and discussion

3.1 Reference bulk optical properties of the liquid phantoms

For each of the 16 liquid phantoms, the BOP in the wavelength range from 550 to 950 nm were estimated with the inverse adding-doubling algorithm from the DIS and unscattered transmittance measurements. The mean results for 3 replicates are presented in Fig. 2 , where the absorption coefficient (µ_a) spectra are grouped (same color) according to the four absorption levels (1 – 4), while the reduced scattering coefficient spectra (µ_s’) are grouped according to the four scattering levels (A – D).

 

Fig. 2 Reference bulk optical properties estimated from DIS measurements for 16 liquid phantoms: (a) absorption coefficients µ_a spectra grouped according to the absorption level (1 – 4) and (b) reduced scattering coefficients µ_s’ spectra grouped according to the scattering level (A – D).

Download Full Size | PDF

In Fig. 2, the absorption coefficient and reduced scattering coefficient spectra are clearly grouped according to the designed levels of scattering and absorption. This indicates that the inverse adding-doubling algorithm was well able to separate the effects of scattering and absorption from the reflectance and transmittance spectra acquired in the DIS and unscattered transmittance setup.

The absorption spectra corresponding to a given absorption level (1, 2, 3 and 4) are very similar with clear absorption peaks at 615 and 670 nm for respectively the methylene blue dimers and monomers [39], where the dimers are more pronounced for the higher concentrations. The minimum and maximum values at the methylene blue peak were about 0.001 cm⁻¹ (no methylene blue) and 1.2 cm⁻¹ (2.2 ml of 400 µM methylene blue), respectively. In Fig. 2(b), a clear decrease in the reduced scattering coefficient values with increasing wavelength can be observed. This can be associated with the scattering properties of the fat globules of the intralipid emulsion which have diameters ranging from 0.1 to 0.65 µm [50]. At each designed scattering level, the reduced scattering coefficient spectra of the four liquid phantoms with different absorption levels overlap very well. The reduced scattering coefficient values range from about 5.0 cm⁻¹ (2 ml of intralipid) until 17 cm⁻¹ (8 ml of intralipid) at 670 nm. Overall, the reference BOP were accurately estimated from the DIS and unscattered transmittance measurements and were in good agreement with values found in literature [39,45 ]. The use of methylene blue as a molecular absorbing agent (non-scattering) resulted in a better (visual) separation of the absorption and scattering effects compared to studies where indian ink was used. Moreover, the absorption spectrum of indian ink is rather independent of the Vis/NIR spectral wavelength and not easily distinguishable from the effect of scattering [13].

3.2 Hyperspectral scatter imaging SRS profiles of the liquid phantoms

A typical 2D hyperspectral image acquired for liquid phantom B3 is illustrated in Fig. 3 . In this figure, each horizontal line represents the reflectance spectrum for a particular source-detector distance on the scanning line [Fig. 1 top view], while the X-axis represents the wavelength dimension. The intensity of the camera is expressed in digital counts (DC), which ranges from 0 until maximum 4096 (12-bit), is given as a color bar.

 

Fig. 3 Illustration of the 2D hyperspectral image of liquid phantom B3 (µ_a = 0.8 cm⁻¹ and µ_s’ = 9 cm⁻¹ at 670 nm). The X-axis represents the wavelength dimension, while each horizontal line represents the reflectance spectrum for a particular source-detector distance on the scanning line.

Download Full Size | PDF

The extracted SRS profiles in the 550 – 960 nm range for the B* phantoms, with the same scattering property, but different absorption levels, are illustrated in Fig. 4 . A clear decrease of the diffuse reflectance values with increasing distance can be observed at each wavelength. This can be explained by the facts that the radial area increases with increasing distance from the point of illumination. Furthermore, the photons exiting the sample at a further distance from the point of illumination have travelled a longer path through the sample and thus have had more chance to be scattered and/or absorbed. Also, for phantoms B2, B3 and B4 a dip in reflectance around 670 nm can be observed. Moreover, the reflectance dip at 670 nm is more significant if the dye concentration is increased as it relates to absorption by the methylene blue dye. The reflectance decreased for wavelengths above 900 nm, indicating the presence of water with an absorption peak at around 980 nm.

 

Fig. 4 Hyperspectral SRS profiles for the liquid phantoms of group B with the same scattering coefficient level µ_s’ = 9 cm⁻¹ and different levels of absorption: (a) B1: µ_a = 0 cm⁻¹, (b) B2: µ_a = 0.4 cm⁻¹, (c) B3: µ_a = 0.8 cm⁻¹, (d) B4: µ_a = 1.2 cm⁻¹ at 670 nm.

Download Full Size | PDF

3.3 Bulk optical property estimation of the liquid phantoms

The performance of the forward metamodels is presented in Fig. 5 as scatter plots that show the predicted versus measured diffuse reflectance values at four selected source-detection distances. The validation is done for all wavelengths of both the calibration (blue) and validation phantoms (red). The reflectance values predicted by the metamodels are very close to the measured reflectance values. This good performance can also be noticed from the high R ² and low RMSE values for calibration ($R_C^2$ = 0.969 – 0.999 and RMSE_C = 0.089 – 0.596%) and validation ($R_V^2$ = 0.944 – 0.999 and RMSE_V = 0.123 – 0.581%). Moreover, for the different source-detector distances, both the R ² and RMSE for calibration and validation set are nearly similar, indicating the high robustness of the metamodels. The further the distance from the illumination spot, the lower the prediction capability of the metamodel. This is probably due to the decrease in signal-to-noise ratio (SNR) with increasing source-detection distance. As a result, some deviation of the points was observed for the furthest distance [Fig. 5(d)], resulting in a slightly lower R ² value of 0.969 (calibration) and 0.944 (validation). However, as the signal at that distance was relatively low and noisy, this is acceptable. Therefore, it can be concluded that the forward metamodels was able to predict well the measured SRS profiles for the liquid phantoms based on the corresponding BOP.

 

Fig. 5 Scatter plots of predicted versus measured reflectance values for 12 liquid calibration phantoms (blue) and four validation phantoms (red) at four source-detector distances: (a) 0.1 cm, (b) 0.2 cm, (c) 0.4 cm, and (d) 0.8 cm.

Download Full Size | PDF

The inverse model was validated by comparing the estimated BOP of the 4 validation phantoms with the reference values obtained from DIS and unscattered transmittance measurements. The SRS profiles were obtained from the measured hyperspectral scatter images, while the BOP were estimated from the SRS profiles through the iterative inversion of the forward metamodels. In Fig. 6 , the predicted BOP for these validation phantoms are plotted against the corresponding reference values. The obtained results (red and blue dots) show that both the predicted absorption and reduced scattering coefficient values for the four validation phantoms closely match the reference values [section 2.2] with $R_P^2$ values of 0.912 and 0.997, and RMSE_P of 0.068 cm⁻¹ and 0.226 cm⁻¹ for µ_a and µ_s’, respectively. Moreover, the relations between the measured and reference BOP are indicated by the blue lines in Fig. 6.

 

Fig. 6 Scatter plots of predicted versus measured bulk optical properties for the four validation liquid phantoms: (a) absorption coefficients µ_a ; (b) reduced scattering coefficients µ_s’. The red dots represent the data of the 570 – 900 nm range, while the blue dots represent the 550 – 570 nm and 900 – 950 nm range. The blue line is a linear fit to all the data (blue and red dots), while the red line is the linear fit to the data of the reduced wavelength range (red dots).

Download Full Size | PDF

To get a better view on the prediction performance in function of the spectral dimension, the predicted BOP spectra for the four validation phantoms (A3, B1, C2, and D4) are illustrated in Fig. 7 as a function of the wavelength, together with the reference BOP spectra.

 

Fig. 7 Predicted and measured bulk optical properties spectra for the four liquid validation phantoms: (a) absorption coefficients µ_a; (b) reduced scattering coefficients µ_s’.

Download Full Size | PDF

Visual inspection of the curves in Fig. 7(a) shows a very good match between the absorption coefficient spectra predicted from the hyperspectral scatter images and those obtained from the DIS measurements. Especially, the absorption peak of methylene blue is predicted very well for the different phantoms, without crosstalk between the reduced scattering and bulk absorption coefficient spectra. However, the prediction of the µ_a values in the 550 – 570 nm and 900 – 960 nm range was less reliable [Fig. 7(a)], while the µ_s’ of the highest scattering level (phantom D) was slightly underestimated for the wavelengths below 570 nm [Fig. 7(b)]. This can be explained by the low SNR of the measured reflectance values, mainly due to the low sensitivity of the camera at these wavelengths. Additionally, above 900 nm, the SNR is even further reduced due to high absorption by the sample’s water molecules. To illustrate the effect of the wavelength range on the model performance, the data points corresponding to the 570 – 900 nm range (higher SNR) were indicated in red [Fig. 6] and the performance parameters were recalculated. Moreover, the red lines [Fig. 6] represent the linear relation between the predicted and measured BOP, resulting in an $R_P^2$ of 0.980 and 0.998, and RMSE_P of 0.032 cm⁻¹ and 0.197 cm⁻¹ for µ_a and µ_s’, respectively. This clearly indicates that the model performs even better in the 570 – 900 nm, were the SNR of the SRS profiles is sufficient.

Apart from the good prediction performance of the inverse modeling approach, also the required computation time is important for practical usage. Training the metamodels is the most computational intensive step, taking about 52 minutes on a standard desktop computer (Intel® Core i5-4570 CPU @ 3.2 GHz) for a single source-detector distance or single metamodel. In total, as 21 metamodels were constructed, the entire calibration process took approximately 18 hours to complete. Once the metamodels were constructed, it only required 9.4 milliseconds on the same computer to estimate the SRS profile (21 source-detector distances) for a single set of BOP. Moreover, the calculation time is independent of the BOP itself. This is substantially faster relative to the seconds or minutes which is required by the optimized MC algorithms [29,31 ]. Nevertheless, as the used Nelder-Mead simplex optimization is relatively simple, many sets of BOP (1305 on average) were evaluated to find the BOP which correspond to a measured SRS profile. Therefore, the inverse estimation required 12.5 seconds on average to complete, which is too still high for predicting the sample’s BOP in an on-line application. Therefore, future research should focus on the optimization of this inverse estimation. For instance, the BOP estimate at a neighboring wavelength can be used as a starting point for the optimizer since the BOP of neighboring wavelengths are highly correlated. This could reduce the number of required starting points to find the global minimum corresponding to the actual BOP.

Finally, it should also be noted that the applicability of the approach followed in this study is limited to turbid media which are homogeneous on a scale which can be visualized with the imaging system ( ± 50 µm) [51]. Nevertheless, this approach could be further extended to cases of multi-layered sample structures. However, to train the extended model, a significantly larger set of calibration phantoms would be required to cover all the possible variability (BOP and thickness of each layer).

4 Conclusions

In this study, a metamodeling approach based on the stochastic Kriging concept has been elaborated and combined with an iterative inversion to estimate the bulk optical properties from spatially resolved reflectance profiles obtained with contactless hyperspectral scatter imaging. Forward validation of these models showed that the predicted reflectance values matched well with the measured ones (R ² ≥ 0.944). Next, these metamodels were incorporated in an to estimate bulk optical properties from the measured reflectance profiles. Validation of the iterative inversion scheme showed accurate and fast prediction of the absorption coefficient and reduced scattering coefficient for an independent set of validation phantoms, with $R_P^2$ values of respectively 0.980 and 0.998, and RMSE_P values of 0.032 cm⁻¹ and 0.197 cm⁻¹. This confirms the potential of the inverse metamodeling approach for rapid and accurate bulk optical property estimation of turbid media. Future research should be focused on validating this approach for actual biological tissues and samples.