1.

Introduction

The biomedical community relies on phantoms mimicking the optical properties of tissues to develop, test, and/or calibrate biomedical optical devices. These tissue-simulating phantoms can be built using different types of scatterers [${\mathrm{TiO}}_2/{\mathrm{Al}}_2{\mathrm{O}}_3$ powder, lipid emulsion (i.e., intralipid), size-calibrated microspheres of polymer, or glass], an absorber (ink, molecular dyes), and supporting materials (aqueous suspension, silicon, polyurethane resin).^1–4 The chosen components are then mixed and sometimes cured (solid phantoms). Ultimately, the optical properties of the resulting phantoms (absorption coefficient ${\mu }_a$; scattering coefficient ${\mu }_s$; scattering anisotropy factor $g$; refractive index $n$) need to be measured.

The base material refractive index, $n$, can be estimated by an independent measurement² or is known from the literature.⁵ Different model-based approaches can be used to infer the remaining optical parameters from experimental data. Photon transport in biological tissues and phantoms can be simulated by Monte Carlo modeling^6,7 that launches and tracks a large number of photons subjected at random to absorption and scattering events. Monte Carlo simulations can accurately model the experiment by considering every significant parameter and the specific geometry of the sample but are computationally time consuming as are inversion calculation algorithms based on Monte Carlo simulations.^8–10 The optical properties of the sample can be obtained by solving the radiative transport equation (RTE), which is difficult but possible under some simplifying hypotheses. Prahl’s adding-doubling (AD) procedure¹¹ solves the RTE for samples with homogeneous optical properties in an infinite parallel-plane slab geometry and does not consider the polarization of light. The AD algorithm first solves the RTE for a single layer in the single-scattering approximation and computes the total reflectance and transmittance. The reflectance and transmission for a slab twice as thick are achieved by successively adding/doubling the thin slab with proper boundary conditions, and this iteration is repeated until the final desired thickness is reached. The inverse adding-doubling (IAD) routine¹² solves for the optical properties of the sample using the measured total reflectance, total transmittance, and additionally the measured unscattered transmittance of ballistic photons that follow Beer–Lambert’s law. Another simplifying approach in solving the RTE is to use the diffusion approximation valid for a medium with high albedo, i.e., ${\mu }_s\gg {\mu }_a$. Under these conditions, the RTE is simplified into the diffusion equation that depends on ${\mu }_a$, and on the reduced scattering coefficient, ${\mu }_s^{\prime }={\mu }_s(1-g)$.^13,14 The diffusion equation can be formally solved for different types of photon sources.¹⁵

In association with these models, many approaches have been used to study the optical properties of biomedical tissues and phantoms. Time-resolved experimental techniques based on the temporal spreading of light pulses while traveling in the sample offer the natural ability to decouple the absorption and scattering effects, since these events differently affect the distribution of photons.¹⁶ Indeed, the solution of the time-dependent form of the RTE for a temporal Dirac delta function photon source is the product of the response for an absorption-free medium and a factor ${\mathrm{e}}^{-{\mu }_{a}vt}$ (where $v$ is the speed of light).¹⁷ This property has been used to precisely measure ${\mu }_a$ of liquid intralipid/India ink phantoms by using the added absorber technique in which calibrated amounts of absorber are added to the scattering solution.¹⁸ The natural logarithm of the ratio of the signals measured for two different concentrations of absorber is then linearly dependent on $\mathrm{\Delta }{\mu }_a$, the difference between the corresponding ${\mu }_a$ values. The absorption coefficient per unit concentration can be accordingly estimated. Some authors also used this property of the RTE to optimize the computation time of look-up-tables (LUT) representing the distribution of photon’s time-of-flight as a function of the optical properties of the sample. They did so by limiting Monte Carlo simulations to absorption-free media with different scattering properties, the absorption contribution being later introduced by multiplying the results by ${\mathrm{e}}^{-{\mu }_{a}vt}$.^19,20 Within the diffusion approximation, measurements of liquid phantoms in the time-domain with different source-detector separation have been carried out in reflection for infinite and semi-infinite scattering media^20,21 and in transmission for liquid phantoms contained in a cuvette assuming an infinite plane-parallel slab geometry.¹⁸ Transmission measurements in the time domain are favored since they are less sensitive to detector-source interactions.¹⁹ Solid phantoms have been analyzed with time-domain reflection measurements of semi-infinite scattering media^22,23 and transmission measurements of turbid slabs.^24–26 Using an inverse Monte Carlo procedure, Bouchard et al.¹⁹ measured ${\mu }_a$ and ${\mu }_s^{\prime }$ of solid tissue-simulating phantom slabs manufactured by the Institut National d’Optique (Biomimic™, INO Inc., Quebec, Canada) with respective uncertainties (coverage factor $k=1$) of about 6% and 3.5%.

Frequency domain photon migration (FDPM) measurement methods have been developed using an intensity-modulated light source and the detection at some distance from the source of the shift of the mean value of the signal, the amplitude attenuation, and the phase shift relative to the incident light. Solving the diffusion equation for a point source with sinusoidally modulated intensity in the infinite medium approximation or the semi-infinite medium approximation with extrapolated boundary conditions, ${\mu }_a$ and ${\mu }_s^{\prime }$ could be determined by fitting the model to the data. FDPM offers the possibility of measuring at multiple source-detection distances using a single modulation frequency^27,28 [multiple distance (MD) method] or measuring at a fixed source-detection distance with multiple modulation frequencies²⁹ [multiple frequency (MF) method]. The MF method implies a nonlinear fitting procedure of the model to the data, whereas for the MD method, the model is linearly dependent on the distance between the detection points.³⁰ The MD method probes different volumes of the sample, which causes variations of the optical properties in the case of heterogeneous samples. The MD method is favored in the case of homogeneous samples, whereas the MF method can be used for heterogeneous samples.³¹

Spatial-domain measurement techniques of the optical properties of liquid biomedical phantoms have been carried out in the continuous wave (CW) domain, based on source-detector multidistance measurements of fluence rate under the diffusion approximation and for point-like isotropic sources in infinite medium. Under these assumptions, an effective attenuation coefficient, ${\mu }_{\mathrm{eff}}=\sqrt{3{\mu }_a{\mu }_s^{\prime }}$, was estimated. The scattering properties of the samples were measured with the added absorber technique³² at visible wavelengths. The method of water absorption is similar and can be expanded to the near-infrared regime as well.^33,34 For different dilutions of the diffusive medium, i.e., different water concentrations, and assuming that the absorption coefficient of water was known precisely, the optical properties of the diffusive medium were derived from the measurement of ${\mu }_{\mathrm{eff}}$. Both methods were aimed at uncoupling absorption from scattering coefficients in ${\mu }_{\mathrm{eff}}$. Based on this measurement technique, a diffusive liquid reference standard has been proposed using intralipid as a scatterer with well-defined optical properties, ${\mu }_a$ and ${\mu }_s^{\prime }$, to be measured each with variations smaller than 2%.³⁵ One should note that both the added absorber technique and the water absorption method are particularly suited for liquid phantoms. However, liquid phantoms are not ideal for dissemination or round-robin tests; therefore, solid phantoms are preferred as reference standards.

Fiber-based multidistance reflectance measurements of tissue or solid biomedical phantoms in the CW domain require that great care is taken to ensure optimal optical coupling at the sample-fiber interfaces.³⁶ Also, measuring tissue in that experimental situation requires that the contact pressure between the sample and the fibers remains low in order to not modify the optical properties. Fixed distances between the source and detection fibers imply a limited number of measurement points that make prior knowledge of the sample optical properties necessary to optimize fiber separation. To circumvent these problems, noncontact measurement techniques allowing continuous source-probe interdistances have been developed using a laser beam or a projected point source and imaging the sample surface or a single surface point at MDs.^37,38

Structured illumination is a noncontact alternative version of the measurement of the optical properties of tissues or phantoms. Under normal incidence, a sinusoidal pattern is projected on the samples. Assuming a linear medium, the diffusion equation becomes a one-directional second-order Helmholtz equation for the fluence rate as a function of sample depth, yielding an expression of the diffuse reflectance dependent on the optical parameters, ${\mu }_a$ and ${\mu }_s^{\prime }$, and the modulation frequency. The diffuse reflectance shows a different sensitivity to the absorption and scattering properties of the sample for the low- and high-frequency regimes of projected patterns, respectively, allowing the uncoupling of ${\mu }_a$ and ${\mu }_s^{\prime }$.^39–41 Even though this technique is based on the assumption that the diffuse approximation is valid, forward Monte Carlo simulation showed that this limit can be tweaked to allow using structured illumination of a variety of samples with lower albedos.

Measurements of the total reflectance and total transmittance of thin homogeneous biomedical samples can be performed with integrating spheres. These measurements are then used as input data to a light propagation model inversion routine to estimate ${\mu }_a$ and ${\mu }_s^{\prime }$. Technically, only a single-integrating sphere is necessary to measure the total reflectance or the total transmittance at the expense of repositioning the sample or the sphere,^42,43 both of which may add to the overall uncertainties. For improved measurement accuracy, measurements with double-integrating spheres (DIS) are favored, since the total reflectance and total transmittance are acquired simultaneously with the sample in the same position, under the same illumination conditions and mechanical constraints. However, relevant application of the light propagating model in a DIS setup requires rigorous control of the experimental conditions (e.g., ensuring good contact between the sample and the sample port of both spheres, mitigated finite-size effect of the sample, etc.), so that the model-based interpretation of the measured optical properties of the sample yields accurate values of ${\mu }_a$ and ${\mu }_s^{\prime }$. The optical properties of the sample are solved using a Monte Carlo inversion routine^10,44,45 or an IAD routine.^46–48 An additional measurement of the unscattered transmittance allows independent determination of ${\mu }_s$ and $g$. However, low signal levels make this measurement difficult for optically thick samples. Also, multiple scattered photons can be collected leading to erroneous values of the unscattered transmittance. Consequently, the anisotropy factor is often an educated guess from size distribution and Mie theory² or is borrowed from literature values.^9,49

At this time, there is no National Institute of Standards and Technology (NIST)-traceable measurement standard for the calibration of biomedical optical devices, where the key measurements are tissue absorption and scattering properties. This report presents progress toward filling this gap. Our choice of the DIS measurement technique in association with the inversion of the AD algorithm for the implementation of this national reference instrument was based on the ability to cover a broad range of optical properties, offers a fast inversion procedure, and allows for the measurement of solid samples as well as liquid phantoms in a cuvette. Moreover, the DIS measurement technique is even qualified as the “golden standard” method by some authors⁴⁸ in reference to the numerous studies it was based on. Hence, we establish a DIS-based measurement platform and a data analysis, graphical user interface (GUI), in which variations of experimental configurations and model-based analysis to compute absorption and scattering properties of the samples are implemented. Analysis of the data (total reflectance and total transmittance) is performed with an inversion procedure similar to Prahl’s IAD,¹² but our procedure allows for the computation of a complete uncertainty budget for each sample measurement. In our previous work,⁵⁰ the development of our setup and the method used to compute an uncertainty budget for the measurement of the optical properties of a sample were explained. Since repeatability and validation of the measurement results are essential qualifications of a standardized measurement platform; in this work, we expand to include the measurement of samples with different thicknesses and compare the results to the optical properties measured by the phantom manufacturer. For repeatability tests, three commercially available INO solid biomedical phantoms having the same optical properties and the same lateral sizes but different thicknesses were measured in a CW regime at different sample orientations and polarization angles of the incident light. For validation tests, the results are compared to those measured by the sample manufacturer on an exemplar sample of the same batch.

2.

Materials and Methods

2.1.

Experimental Setup

The experimental setup is shown in Fig. 1. The polarization of an incident HeNe laser beam ($\lambda =543.5\mathrm{nm}$, beam diameter 0.80 mm, beam divergence 0.86 mrad, minimum output power 0.75 mW, JDS Uniphase Corporation, Milpitas, California; or $\lambda =632.8\mathrm{nm}$, beam diameter 0.80 mm, beam divergence 0.01 mrad, minimum output power 5 mW, Research Electro-Optics, Inc., Boulder, Colorado) is controlled by a linear Glan–Taylor polarizer, P, which is followed by a beam splitter, BS. A reference signal is detected by photodiode ${\mathrm{D}}_1$ to compensate for intensity variations of the laser source. The light is at normal incidence to a sample, S, held between two integrating spheres, R and T (UMBK-190, Gigahertz Optik, Türkenfeld, Germany; internal diameter: 196 mm; entrance port diameter: 25.1 mm; sample port diameter: 63.5 mm reduced to 38.1 mm by a removable port diameter reducer; detector port diameter: 12.7 mm; coated with ODM98 synthetic material). The total reflectance signal and the total transmittance signal are measured by photodiodes ${\mathrm{D}}_2$ and ${\mathrm{D}}_3$, respectively (${\mathrm{D}}_1$, ${\mathrm{D}}_2$, ${\mathrm{D}}_3$; blue enhanced response silicon detector, Edmunds Optics Inc., Barrington, New Jersey). The photocurrent is amplified by three P-9202-6 amplifiers (Gigahertz Optik) connected to a DAQ board (NI PCI-6110, National Instruments Corporation, Austin, Texas) and monitored by a LabVIEW program for data acquisition. A NIST-traceable 99% reflectance standard, $R_{\mathrm{Std}}$, blocks the exit port of the transmittance sphere during the measurements.

Fig. 1

Experimental setup (50): P: polarizer; BS: beam splitter; ${\mathrm{D}}_1$, ${\mathrm{D}}_2$, and ${\mathrm{D}}_3$: photodiodes; S: sample; $R_{\mathrm{Std}}$: reflectance standard; R and T spheres: reflectance and transmittance spheres.

3.

Discussion and Conclusion

Figures 2 and 3 present the optical properties ${\mu }_a$ and ${\mu }_s^{\prime }$ measured at $\lambda =543.5\mathrm{nm}$ and $\lambda =632.8\mathrm{nm}$ for three 100-mm-square INO phantoms of thicknesses $d=5.08\mathrm{mm}$, 7.09 mm, and 9.92 mm. To account for reproducibility of the results, each sample was measured eight times at each wavelength by flipping the face hit by the incident laser beam and rotating the sample 90 deg around the $Y$-axis between measurements. The time-domain transmittance technique used by INO requires optically thick samples to allow measurable pulse spreading. This thickness is not applicable in our DIS-based setup, so our results from thinner ($<10\mathrm{mm}$) samples are compared to the results provided by INO at the same wavelengths on a 20-mm-thick phantom with identical optical properties. The uncertainties in the INO results were estimated from the measurements at $\lambda =660\mathrm{nm}$.¹⁹ Figures 2 and 3 present the results with $1\sigma$ error bars (i.e., coverage factor $k=1$). At $\lambda =543.5\mathrm{nm}$, our results are consistent with the values from INO for both ${\mu }_a$ and ${\mu }_s^{\prime }$. At $\lambda =632.8\mathrm{nm}$, the agreement is only for $d=9.92\mathrm{mm}$ for both ${\mu }_a$ and ${\mu }_s$. However, with a coverage factor $k=2$, preferred by NIST,⁵³ the optical properties measured here agree with the results provided by INO. Table 1 summarizes the results and the uncertainty with $k=2$.

Fig. 2

Results at $\lambda =543.5\mathrm{nm}$ of (a) ${\mu }_a$, the absorption coefficient, and (b) ${\mu }_s^{\prime }$, the reduced scattering coefficient, versus the thickness of the sample for 100-mm-square INO phantoms. The error bars (coverage factor $k=1$) are obtained for each experiment by propagation of the type A and type B uncertainties. The NIST results are compared to the values from INO (mean value and error, dotted lines, and dashed lines, respectively).

Fig. 3

Results at $\lambda =632.8\mathrm{nm}$ of (a) ${\mu }_a$, the absorption coefficient, and (b) ${\mu }_s^{\prime }$, the reduced scattering coefficient, versus the thickness of the sample for 100-mm-square INO phantoms. The error bars (coverage factor $k=1$) are obtained for each experiment by propagation of the type A and type B uncertainties. The NIST results are compared to the values from INO (mean value and error, dotted lines, and dashed lines, respectively).

Table 1

Results and uncertainties (k=2) of μa, the absorption coefficient of the sample, and μs′, the reduced scattering coefficient of the sample. The uncertainties on the INO results were estimated from measurements at λ=660.0  nm.19

One major advantage of the time-domain transmittance measurement technique used by INO is that it ensures that the detected photons experience multiple scattering events before detection. Hence, it is less sensitive than DIS to the measurement of the anisotropy factor, $g$. Moreover, the INO procedure uses an inverse Monte Carlo model to take into consideration the photons being reflected or transmitted at the sample boundaries, which is not the case with the AD routine that considers samples of homogeneous optical properties in an infinite plane-parallel slab geometry. However, the geometry of the samples measured here allows for this assumption to be valid.

Our setup compensates for the fluctuations of the laser beam intensity by means of a reference channel that requires a beam splitter and control of the polarized light input by a linear polarizer set at a fixed angular position. Within the range of optical properties of tissues, multiple-scattering events are expected to totally depolarize light.⁵⁴ However, the INO samples have a smooth nonspecular surface and it was a concern that surface reflection captured by the reflectance integrating sphere might show enough of a polarization dependence that the measured optical properties of the sample might be affected. To ensure that the measurements of the reflectance are not sensitive to the chosen polarization angle of the input beam, additional measurements of the optical properties were carried out for different angles of the linear input polarization. For that purpose, a rotating achromatic half-wave plate was set in the incident beam path between the beam splitter and the reflectance sphere. The fast axis of the half-wave plate was initially aligned along the polarization direction set by polarizer P, and the angular orientation of the half-wave plate was incremented by a 15° rotation step. The corresponding polarization orientation, ${\vartheta }_P$, of the incident beam is ${\vartheta }_P=0\mathrm{deg}$,   30 deg, 60 deg, and 90 deg. Figure 4 presents the optical parameters, ${\mu }_a$ and ${\mu }_s^{\prime }$, versus ${\vartheta }_P$ measured for one of the INO phantoms (thickness $d=7.09\mathrm{mm}$) at $\lambda =632.8\mathrm{nm}$. The values present fluctuations well within the error bars at $k=2$, so we can rule out any measurable influence of the linear input polarization on the optical parameters measured with our DIS setup.

Fig. 4

Results at $\lambda =632.8\mathrm{nm}$ of (a) ${\mu }_a$, the absorption coefficient, and (b) ${\mu }_s^{\prime }$, the reduced scattering coefficient, versus $\vartheta$, the angle of the input linear polarization ($\vartheta =0\mathrm{deg}$ corresponds to a vertical linear polarization) for the 100-mm-square INO phantoms of thickness $d=7.09\mathrm{mm}$. The error bars (coverage factor $k=1$) are obtained for each experiment by propagation of the type A and type B uncertainties. The NIST results are compared to the values from INO (mean value and error, dotted lines, and dashed lines, respectively).

In this paper, we presented measurements of the optical properties of solid biomedical phantoms at $\lambda =543.5\mathrm{nm}$ and $\lambda =632.8\mathrm{nm}$ and compared the results to the values obtained from their manufacturer, INO. The samples were measured using our double-integrating-sphere setup in the steady-state domain. The measured data consist of the total transmittance and the total reflectance of the sample, and analysis of the data was based on a custom-made inversion routine of the adding-doubling procedure, which solves the radiative transfer equation iteratively for the data. Our inversion routine computed the optical parameters and the uncertainty budget for each measurement. A GUI was developed to facilitate the data analysis. The samples were three square solid phantoms of thicknesses $d=5.08\mathrm{mm}$, 7.09 mm, and 9.92 mm. There is an agreement within the error bars at $k=2$ between our results and the ones from INO at both $\lambda =543.5\mathrm{nm}$ and $\lambda =632.8\mathrm{nm}$ and for all thicknesses. The uncertainty in the optical parameters measured was estimated to be between 4.0% and 5.0% for ${\mu }_a$ and between 11% and 12% for ${\mu }_s^{\prime }$, compared to 12% and 7%, respectively, obtained by INO.¹⁹ With the DIS measurements of the optical properties, the uncertainties on ${\mu }_s^{\prime }$ were found to be mainly sensitive to the uncertainty on the scattering anisotropy factor, $g$, and the uncertainty on ${\mu }_a$ was mainly dependent on the uncertainty on the incident angle.⁵⁰ As previously mentioned, the time-domain technique is less sensitive to measurements of $g$, and consequently, the uncertainty on ${\mu }_s^{\prime }$ is smaller. Figures 2(a) and 3(a) show ${\mu }_a$ values that are thickness independent, whereas in Figs. 2(b) and 3(b), ${\mu }_s^{\prime }$ values are also thickness independent within the error bars. However, the mean values of the reduced scattering coefficient have a trend of increasing as a function of the sample thickness. A close consideration of this finite sample-size effect needs to be addressed for more careful analysis, which is the scope of future work. The well-collimated laser source used in our measurements helps reduce the uncertainty in the absorption coefficient.⁵⁰ Future work in developing this facility will include extending the capabilities to cover a broad spectral range, analyzing liquid samples, and implementing an in-depth analysis of the associated uncertainties. Interlaboratory comparisons with other national metrology institutes will be investigated.