Combination of diffuse reflectance and transmittance spectroscopy to obtain optical properties of liquid phantoms

Abstract. We present a simple fiber-optics probe system that could be used in any lab for convenient determination of optical properties of liquid phantoms based on diffuse reflectance and transmittance measurements in the visible/near-infrared region. We employed Monte Carlo simulations to determine the optimal system setup and to test the inverse algorithm employed to extract the optical properties from measured reflectance and transmittance. The inverse algorithm involved obtaining the fit merit function for values within the optical property range and determining the minimum. The performance of the method was tested by predictive error and validated using similar matrix of milk–ink phantoms on reflectance and transmittance. In the range of optical properties of phantoms with optical properties of 0 to 0.5  cm  −  1 for μa and 20 to 140  cm  −  1 for μs, the median prediction error for the test phantoms at 630 nm was 1.51% for μs and 8.82% for μa. The median difference in predicted values versus expected values was 1.15  cm  −  1 for μs and 0.01  cm  −  1 for μa. In comparison with other techniques, our method was a simple, fast, and convenient way to determine optical properties of liquid phantoms.


Introduction
The determination of optical properties for biological tissue is vital to the development of medical diagnostic devices. The absorption coefficient μ a may provide information about the tissue chromophores, whereas the scattering coefficient μ s and anisotropy factor g may be used to characterize the various scattering components of the tissue. Convenient and accessible methods for measuring optical properties are extremely useful for research in areas not only for diagnostic purposes, (e.g., variation of light absorption in tissue to monitor changes in blood oxygenation and concentrations of oxygenated and deoxygenated hemoglobin) 1 but also in therapeutic applications (e.g., knowing the optical properties of the tissue for calculating the light dosage delivered to tissue by laser therapy). 2 Optical phantoms are often used to assist with the design and calibration of medical optical systems. Such phantoms have defined optical properties and are easy to use within the optical setup. Since optical phantoms are used quite often, the optical properties must be known in order to successfully characterize a system. With the case of liquid phantoms, sometimes the optical properties can change over time or from batch to batch of scatter or absorbers used. Therefore, we propose a simple optical setup for the purpose of determining the optical properties of liquid phantoms quickly, efficiently, and at a low cost. The proposed system is a very simple system that can easily be assembled, calibrated, and can perform measurements within a day. The system does not require a large budget and could be used for any lab, even for those who are not familiar with optics. Although within this paper we present one wavelength, the use of spectrometers allows us to obtain optical properties over a range of wavelengths. In addition to providing the optical properties, this system could also be a valuable educational tool to demonstrate how to create and test the optical properties of simple liquid phantoms due to its easy setup.
A variety of techniques used to measure optical properties have been reviewed. 3,4 A common method is based on integrating sphere set up, which may consist of a single or double integrating spheres. 5,6 This method acquires spectra from three measurement modes, namely collimated transmittance, total transmittance, and total diffuse reflectance, followed by the inversion of these measurements using the adding-doubling method. A technique based solely on diffuse reflectance measurement has significant potential in medical diagnostics and monitoring as it offers the advantages of noninvasive or minimally invasive measurements. Principles and applications of tissue optical reflectance and transmittance have been comprehensively discussed by Wilson and Jacques. 7 Different instrumental techniques based on diffuse reflectance measurements have been devised, using time-resolved, 8 frequency-domain, 9,10 and continuous wave 11,12 methods. However, these methods require bulky and expensive equipment, as well as a large sample volume.
In addition to the equipment, look-up tables have been used in diffuse optical spectroscopy to determine the optical properties of tissue. Some studies used Monte Carlo (MC) simulations, 13,14 whereas others used calibration phantoms 15,16 to create the look-up table. In these studies, only diffuse reflectance measurements were used. Here we use both diffuse reflectance and diffuse transmittance measurements of a standard cuvette to create a look-up table for determining optical properties. The reflectance and transmittance values would have different trends as the scattering or absorption coefficient increases (e.g., increase in reflectance and decrease in transmittance). Therefore, the look-up tables would possibly provide unique information regarding the optical properties based on the measurement values.
In this paper, we present a simple, convenient fiber-optics probe system for the determination of optical properties of liquid phantoms based on diffuse reflectance (R) and transmittance (T) spectroscopy measurements. The system collects reflectance and transmittance measurements simultaneously for wavelengths in the visible and near infrared range for a small volume of liquid without moving the sample. 17 A second system was used to determine the absorption and scattering properties separately before mixing the different components of the liquid phantom to be tested with the first system. MC simulations were used to determine the optimal setup parameters and test the inverse algorithm. Then a calibration model was established with a set of milk-ink phantoms with optical properties within a range typical for biological tissue. Then the inverse algorithm, using the phantom calibration look-up-table, was performed to extract μ a and μ s from the diffuse reflectance and transmittance measurements. The performance of the method was further tested by predictive error and validated using a test set of milk-ink phantoms to obtain reflectance and transmittance. Figure 1 shows the schematic view of fiber-optic system setup used for the diffuse reflectance and transmittance measurements. The system consists of a halogen light source (HL-2000, Ocean Optics, B.V, Netherlands), bifurcated fiber-optic probe (Thorlabs GmbH, Deutschland), a cuvette holder and spectrometers (FLAME-S-VIS-NIR-ES, Ocean Optics, B.V, Netherlands) with a wavelength grating from 350 to 1000 nm, and a laptop computer. Only one wavelength is needed for measurements, so the spectrometer can be replaced with a photodiode. However, with the spectrometer, optical properties of several wavelengths can be obtained all at once. The fiberoptic probe contained two 1000-μm diameter fibers for light delivery and light collection, with an outer diameter of 2.1 mm. The center-to-center distance between the source and detector fibers for reflectance is 1065 μm. The fiber length was <1 m and attenuation within the fiber was considered negligible. There were three cuvette sizes that were tested with 1-mm, 5-mm, and 1-cm width. The 5-mm width cuvette was used for the phantom calibration. The cuvette holder was printed by 3-D printing technology, which considered the difference in cuvette thickness. The turbid medium samples were placed into the cuvette for the measurements. The input light power of the spectrometer for the reflectance measurement was calibrated with the input light power of the spectrometer for the transmittance measurement using a waterfilled cuvette before the turbid medium measurement. To minimize any interference from the background light, the dark measurements were recorded prior to the actual measurements and subtracted from the raw data. We collected the total diffuse reflectance and transmittance data of one sample simultaneously using two spectrometers. Reflectance and transmittance were calculated by the following equations: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 2 8 5

Fiber-Optic System Setup to Measure Optical Properties of Liquid Phantoms
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 2 3 0 where I r is the measured diffuse reflectance intensity, I t is the measured transmittance intensity, I 0 is the input light power that subtracted the dark measurement, and I 0 0 is calibrated input light power for reflectance spectroscopy measurement.

Collimated Transmittance Setup
With the collimated transmittance setup, the absorption coefficient (μ a ) of India ink (Winsor & Newton ® black India ink) and the scattering coefficient (μ s ) of milk (Dawn ® low-fat milk, produced in Ireland) can be measured one by one with a broad wavelength range. The purpose of this system (Fig. 2) was to determine μ s and μ a separately to determine the optical properties of the ink and milk before mixing the phantom together. Figure 2 shows a sketch of collimated transmittance setup to measure the total attenuation coefficient (μ t ¼ μ a þ μ s ). A halogen lamp (HL-2000, Ocean Optics, B.V., Netherlands), spectrometer (FLAME-S-VIS-NIR-ES, Ocean Optics, B.V, Netherlands), and collimation lens (LA4647, Thorlabs GmbH, Deutschland) were used in this system. 18 The distance between the sample and detector was large enough to allow only a small acceptance angel (0.46 deg) in order to minimize the collection of scattered light.
The total attenuation coefficient μ t is calculated from the Beer-Lambert law: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 5 4 5 IðλÞ where d is the sample path length, I 0 is the incident intensity, and IðλÞ is the measured intensity. For solely absorbing samples, μ a ≫ μ s , the measured attenuation coefficient equals the absorption coefficient, μ t ¼ μ a . For solely scattering sample, μ s ≫ μ a , therefore, μ t ¼ μ s , like it is the case for milk dilutions. To avoid the influence of multiple scattered light on the measurement, either the sample can be diluted or the thickness of the cuvette can be decreased. In this study, the milk was diluted with water to a concentration of <5% (volume, milk to water), which is low enough to avoid the influence of multiple scattered light on the measurement. This was tested by measuring the linearity of the attenuation of light as a function of concentration. We compared our experimental results using set up of Fig. 2 with Mie theory simulation results and they match well (data not shown).
We have performed several transmittance measurements on different concentrations of ink and milk, and a linear fit was made for each. The absorption coefficient of India ink and the scattering coefficient of low-fat milk have been determined by extrapolating absorption/scattering coefficient values obtained in an added solution series experiment at low ink/milk concentrations. We assume that the absorption and scattering coefficients are not altered by the mixing.

Phantom Samples
The scattering and absorption coefficients of the milk and the India ink were determined by the collimated transmittance setup described in Sec. 2.2 before mixing. The obtained scattering/ absorption coefficient of the phantoms were derived as a combination of the milk/ink with pure water, each of which the concentration by volume in the phantoms. The absorption of pure ink is much higher than that of typical biological substances, therefore, the inks were prediluted into water solution and we used it as a stock ink-water solution. The absorption coefficients μ a of the ink dilutions were calculated as 85% of the attenuation coefficients μ t because 15% of the attenuation was attributed to scattering by the carbon particles in the original ink solution. 19 In order to create the look-up table measurements to form the basis for the calibration with the liquid phantoms prepared as aqueous mixtures of milk and Indian ink were divided, there were two groups: a calibration set of a matrix of 7 × 7 liquid phantoms (49 samples) were prepared for building the look-up table and a validation set of similar matrix phantom samples to perform the validation experiments. The liquid phantom set was prepared to cover the range of scattering coefficients at 630 nm with 20 to 140 cm −1 and absorption coefficient values at 630 nm with 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, and 0.5 cm −1 .

Monte Carlo Simulations
MC simulations were run using CUDAMCML, 20 which is a GPU-based accelerated form of the original MCML code. 21 The simulation was set up so that light travelled from a point source into a medium of finite length (infinite radius) from air (n ¼ 1) to water (n ¼ 1.33). For the look-up table, 1 billion photon packages were launched. To test the optical properties algorithm, only 1 million were used. The simulation with 1 billion photon packages took about an hour, whereas the 1 million photon packages took about <1 min. The anisotropic factor was assumed to be 0.9, similar to biological tissue. The absorption coefficient was varied from 0.1 to 1.0 cm −1 in steps of 0.1 cm −1 . The scattering coefficient was varied from 10 to 100 cm −1 in steps of 10 cm −1 . The diffuse reflectance and diffuse transmittance measurements from the MC simulations were gathered from the detector locations and integrated over the detector area. There were two measurements: one in reflectance and the other in transmittance. This formed the basis of the look-up table. Then there were >600 MC simulations run with random optical properties to test the accuracy of the inverse algorithm using the MC look-up table.
To determine the optimal system parameters for the system, CUDAMCML simulations were used to model three different measurement geometries. Cuvettes of 1-mm, 5-mm, and 1-cm width were considered. The transmission fiber was either considered to be directly across from the source fiber [Fig 1(c); detector 1] on the other side of the cuvette or slightly off centered at the same distance as the source and reflectance fibers were [Fig 1(c); detector 2]. Either detector 1 or detector 2 was chosen at one time for testing to see which detector yielded more accurate results. A look-up-table was created for each detector geometry (reflectance and detector 1, reflectance and detector 2). To evaluate the ability to predict the optical properties for the different geometries, a training set was used to generate a look-up table that was used to assess the accuracy of the predictions. Since there was difficulty matching the experimental data with the MC simulations using a single calibration factor (see Fig. 3), experimental data with known optical properties were used to generate the look-up table to determine sequential phantom optical properties. The MC simulations were used to create the look-up-table that would help determine the optimal geometry of the setup and to test the algorithm to determine optical properties. Therefore, a new look-up-table was created for the experiment with phantoms of known optical properties.

Inverse Method Using Interpolation of MC Simulations
The information from the MC simulations was used to create a look-up table for determining the design parameters for the system and to test the inverse algorithm. The diffuse reflectance (R) and diffuse transmittance (T) outputs were formatted into matrices, in which the columns and rows corresponded to particular absorption and scattering coefficient inputs (Fig. 3). These reflectance and transmittance matrices can be used to calculate the error function. However, since the final answer is unlikely to be at the same exact values that were used for the MC simulations, the interpolation (MATLAB function: interp2) of optical properties in between the simulated values was calculated. 22 Therefore, the look-up table becomes the diffuse reflectance MC measurements (R) and the diffuse transmittance MC measurements (T) with interpolated points (see Fig. 3). The entire fit merit function χ2 can be calculated using E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 7 2 3 where the value R M was the reflectance measurement and T M was the transmittance measurement taken by the probe. Finally, σ R and σ T are the standard deviation of R and T, correspondingly. The standard deviation was taken for the whole matrix to output a singular value. The objective function can be quickly and easily calculated for all interpolated points. Then the minimum of the fit merit function can be located, which should correspond closely to the optical properties of the input measurement. To increase the accuracy of the solution, these steps are repeated with a finer interpolation around the points that correspond to solution (in order to decrease the amount of memory used). The optical properties that correspond to the new minimum of the objective function were the final output (MATLAB code available in CodeOcean https://codeocean.com/capsule/8373142/tree/v1) (Fig. 4). The range of optical properties was divided into 81 bins (i.e., μ a ¼ 0.1 to 0.2 cm −1 and μ s ¼ 10 to 20 cm −1 ). There were eight simulations run per bin, in which the combination of optical properties was randomly chosen within the range of the bin (i.e., μ a ¼ 0.14 cm −1 and μ s ¼ 18 cm −1 ). A total of 648 random simulations were run using a million photon packages. The same random optical properties were run for each path length that was examined. Fig. 4 Flowchart for the inverse model to obtain optical properties. First, the values between the data points of the MC look-up table are interpolated. The new interpolated matrix is used to determine the fit merit function (χ2). The optical properties that correspond to the minimum of χ2 is determined and the grid surrounding the optical properties is interpolated into a finer mesh and χ2 is determined again for a more accurate solution. Finally, the corresponding optical properties of the second fit merit function minimum is determined.
To test the performance of the algorithm, a prediction test was performed using the simulated diffuse reflectance and transmittance measurement data for the calibration model as input to extract μ a and μ s . The relative predictive errors were calculated as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 6 9 9 relative error ð%Þ ¼ where μ pred is the predicted value and μ ref is the true value of either μ s or μ a .

Monte Carlo Simulation and Geometry Consideration
The diffuse reflectance and transmittance values for the random MC simulations were used to determine the best source-detector distance and cuvette width. The percent error was calculated from the output of the inverse algorithm compared with the input of the MC simulation. Table 1 shows the sum of percent errors of μ a and μ s output values of the algorithm. The lowest error was a cuvette length of 5 mm for detector 2 at a source-detector distance of 1065 μm. With the exception of detector 2 at a cuvette path length of 1 mm, as the sourcedetector distance increased the percent error decreased. Therefore, we chose to use the fibers with 1065-μm source-detector distance and a fiber radius of 500 μm. Figure 5 shows the absolute percent error binned based on the optical properties. For μ a , the highest error was when μ s was 10 to 20 cm −1 or when the scattering coefficient was at the lower range. The absolute error for μ s did not seem to show a pattern, but error is <1%. The mean absolute errors for μ a and μ s were 0.28% and 1.17%, respectively. These results are based on ideal conditions with no noise.
The expected optical coefficient was plotted against the predicted values. The correlation between the expected and predicted values was very good with R 2 ≈ 1 for both μ a and μ s .

Phantom Measurement, Calibration Model, and Forward Validation
As previously demonstrated, 12,23 MC light-propagation simulation models used in conjunction with multivariate analysis can accurately extract optical properties from reflectance measurements. We chose to calibrate the system on a set of liquid phantoms mixed with milk and India ink. Plots of diffuse reflectance and transmittance at 630 nm as a function of the absorption coefficient μ a and the scattering coefficient μ s are shown in Fig. 6. The dots represent the experiment data of reflectance and transmittance measurements of 49 liquid phantoms combining the high and low levels of scattering and absorption coefficient.

Comparison of Predicted Versus True Values Based on Phantom Calibration
The calibration phantoms were used to train the algorithm and then samples were newly created to obtain a test set. The logarithm of the reflectance and transmittance measurements from the test set was fed back into the algorithm and then compared with the known values. Figure 7 shows the fitted correlation plots for both μ a and μ s for the test data. The absorption coefficient has a high R 2 value of 0.93 and a slope of 0.90. The scattering coefficient has an R 2 value of 1.00 and a slope of 0.97. Both predicted coefficients have very high agreement with expected values.   Figure 8 shows the box plots for the percent error for μ a and μ s . The median percent error for μ a was 8.82% while for scattering the percent error was 1.51%. The scattering coefficient had a lower range of percent error compared with the absorption coefficient. The difference between the predicted and expected values for the absorption coefficient is < ∼ 0.05 cm −1 and for the scattering coefficient is <5 cm −1 . The median difference value for μ a and μ s was 0.01 and 1.15 cm −1 , respectively.

Discussion and Conclusions
The inverse algorithm can predict the absorption and scattering properties from diffuse reflectance spectroscopy data of a random sample with very high accuracy with the optimal system parameters. The random samples were generated from the MC simulations, so this was an ideal test case scenario with no noise added into the system. With the use of the phantom samples, the same algorithm (calibrated to the training phantoms) was also able to acquire a relatively accurate prediction of the optical properties. The algorithm was able to predict the scattering coefficient much better than the absorption coefficient. There was higher error for the lower absorption samples. Therefore, the current system may not be able to resolve absorption <0.05 cm −1 or the look-up table needs to have a finer resolution in the lower absorption optical property range. The reason behind this phenomenon was that the reflectance and transmittance measurements changed more rapidly as the scattering coefficient increased compared with the absorption coefficient that does not change very much with increasing absorption (see Figs. 3 and 6). However, the median difference between the predicted and expected μ a was 0.01 cm −1 Fig. 8 Box plots of (a) percent error and (b) difference in the predicted versus expected value for μ a and μ s . and could closely predict the higher absorption values. The phantoms were used to calibrate the system over the MC simulations because there were factors in the optical setup that could not be accounted for in the simulations.
The minimum technique used in the current study to determine the optical properties was an efficient optimization algorithm that does not require an initial guess (as with Newton methods) that can adversely affect the final output. However, the range of absorption and scattering coefficients is limited and determined by the extensiveness of the calibration phantom properties. Therefore, optical properties outside these initial calibration phantoms will most likely not be accurate. The system could be modified or used in a few different ways. First, the accuracy of the absorption might increase if more measurements are added to the system. The additional measurement could be more than 1 mm away from the central line of the incoming light, so as to allow for a longer path length that could lead to more absorption. The increase in absorption events may lead to a greater change in reflectance or transmittance as μ a increases. Also only liquid phantoms can be used in the current setup, but the measurement only requires a small amount to determine the optical properties. Measurements with solid phantoms may be possible, but care would have to be taken in order to create or cut test samples to be put into the holder since a different thickness in the media would give different and incorrect results. Finally, since the system uses spectrometers, the reflectance and transmittance of multiple wavelengths can be determined all at once and can also be placed into the algorithm. Here we present one wavelength for simplicity. This paper demonstrated a simple, convenient fiber-optic probe system for the determination optical properties of liquid phantoms based on diffuse reflectance and transmittance spectroscopy measurement applicable in small sample volume. Results obtained indicated the general concept is achieved. Further study is required to improve the accuracy of μ a , especially for low absorbing phantoms. Nevertheless, in comparison with other methods, our proposed method is simple, fast, and convenient. Additionally, the system setup and liquid phantoms can be a valuable educational tool for the purpose of introducing biophotonics to students who have not had any experience in the field. The current setup could even be modified to have a single-laser light source with photodiode detectors to simplify the system (and reduce cost) to obtain optical properties at one wavelength. Since the code is freely available, the optical properties could be obtained by entering a few values and the optical properties would quickly be outputted. In a similar example, the Beer-Lambert law was explained to students using different concentrations of instant coffee and obtaining the absorption coefficient using spectroscopy. 24 The hands-on experienced helped with the understanding of absorption coefficients and its relation to optical imaging.
The developed system could be used for future real-time determination of optical properties of biological tissue experiment. Furthermore, the measurement data obtained by this approach could be used as an empirical calibration method with a look-up table-based inverse model for accurate extraction of optical properties from a wide range absorption and scattering coefficients of biological tissue. Jacqueline Gunther received her PhD in biomedical engineering from Columbia University in 2016. Her PhD work involved using diffuse optical tomography to monitor breast cancer patients undergoing neoadjuvant chemotherapy and develop methods to predict response to treatment. Also in 2016, she joined the Biophotonics@Tyndall group at Tyndall National Institute, University College Cork, Ireland, as a postdoctoral researcher. Her work involves fundamental research of light propagation in biological media, pulse oximetry, and acoustooptical tomography.
Huihui Lu received her PhD in chemistry from the National University of Ireland. Her prior experience includes medical device industry of nanoparticle-based immunoassay development, point-of-care diagnostics device instrumentation, and in Tyndall National Institute, integrating silicon photonic platform for biosensing applications and multimodality optical techniques for cancer guided surgery. Her research interests include biosensor and assay solutions, silicon photonic integration for on-chip sensing, and multimodality spectroscopy techniques for point-ofcare diagnosis devices and healthcare applications.
Stefan Andersson-Engels is a head of the Biophotonics@Tyndall group at Tyndall National Institute, University College Cork, Ireland, and was appointed there in 2016. He has set up a new Biophotonics Lab at Tyndall, with a number of translational projects in collaboration with clinicians and industrial partners. The more fundamental aspects of the research are toward developing biophotonics tools for deep tissue imaging, including acousto-optical tomography and how to best employ upconverting nanoparticles as a contrast agent. Other areas of general interest are optical spectroscopy and photodynamic therapy.