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 [ 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 ; scattering coefficient ; scattering anisotropy factor ; refractive index ) need to be measured.
The base material refractive index, , can be estimated by an independent measurement2 or is known from the literature.5 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 modeling6,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) procedure11 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) routine12 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., . Under these conditions, the RTE is simplified into the diffusion equation that depends on , and on the reduced scattering coefficient, .13,14 The diffusion equation can be formally solved for different types of photon sources.15
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.16 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 (where is the speed of light).17 This property has been used to precisely measure of liquid intralipid/India ink phantoms by using the added absorber technique in which calibrated amounts of absorber are added to the scattering solution.18 The natural logarithm of the ratio of the signals measured for two different concentrations of absorber is then linearly dependent on , the difference between the corresponding 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 .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 media20,21 and in transmission for liquid phantoms contained in a cuvette assuming an infinite plane-parallel slab geometry.18 Transmission measurements in the time domain are favored since they are less sensitive to detector-source interactions.19 Solid phantoms have been analyzed with time-domain reflection measurements of semi-infinite scattering media22,23 and transmission measurements of turbid slabs.24–26 Using an inverse Monte Carlo procedure, Bouchard et al.19 measured and of solid tissue-simulating phantom slabs manufactured by the Institut National d’Optique (Biomimic™, INO Inc., Quebec, Canada) with respective uncertainties (coverage factor ) 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, and 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 frequency27,28 [multiple distance (MD) method] or measuring at a fixed source-detection distance with multiple modulation frequencies29 [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.30 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.31
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, , was estimated. The scattering properties of the samples were measured with the added absorber technique32 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 . Both methods were aimed at uncoupling absorption from scattering coefficients in . Based on this measurement technique, a diffusive liquid reference standard has been proposed using intralipid as a scatterer with well-defined optical properties, and , to be measured each with variations smaller than 2%.35 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.36 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, and , 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 and .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 and . 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 and . The optical properties of the sample are solved using a Monte Carlo inversion routine10,44,45 or an IAD routine.46–48 An additional measurement of the unscattered transmittance allows independent determination of and . 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 theory2 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 authors48 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,12 but our procedure allows for the computation of a complete uncertainty budget for each sample measurement. In our previous work,50 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.
Materials and Methods
The experimental setup is shown in Fig. 1. The polarization of an incident HeNe laser beam (, beam diameter 0.80 mm, beam divergence 0.86 mrad, minimum output power 0.75 mW, JDS Uniphase Corporation, Milpitas, California; or , 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 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 and , respectively (, , ; 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, , blocks the exit port of the transmittance sphere during the measurements.
The samples are manufactured by INO, Inc. and are commercially available. They consist of a base material, polyurethane, with titanium dioxide () as a scatterer and carbon black as an absorbing component. The samples are designed to match nominal values of optical parameters (, at ). They have a 100-mm-square base and come in three different thicknesses measured with a dial thickness gauge, , 7.09 mm, and 9.92 mm (uncertainty 0.0058 mm). A custom-made holder was used to hold the sample by opposite corners to minimize the contact area. It allows the sample to freely rotate around the - and -axes to ensure flat contact between the sample port of the reflectance sphere and the sample. The lateral position and tilt angle of the sample ports are adjustable to align the incident laser beam normal to the sample surface. The spheres and sample holder are installed on a linear translator aligned along the laser beam to control and minimize mechanical pressure from sample-port contact while ensuring light-tight contact.
The experimental procedure requires measuring the total reflectance and the total transmittance of the sample. The measured total reflectance is the ratio of the light intensity reflected from the sample set between both spheres’ sample ports to the intensity reflected from a NIST-certified nominal 99% reflectance standard, , set at the sample port of the reflectance sphere. The measured voltages for the signal channel and the reference channel are and , respectively, with , Standard. Background values of both the signal channel and the reference channel, and , respectively, are measured by blocking the incident beam right after the laser. Hence,
The measured total transmittance is the ratio of the light intensity transmitted through the sample set between both spheres’ sample ports to the intensity with no sample. The measured voltages for both the signal channel and the reference channel are and , respectively, with , Empty. Background values of both the signal channel and the reference channel, and , respectively, are measured by blocking the incident beam right after the laser. Hence,
The details on how to measure the optical properties can be found elsewhere.2,12,51 The experimental procedure is based on the substitution method, i.e., the sample port is successively blocked by the sample, the reflectance standard, or left empty. One should note that the measurements in integrating sphere throughput efficiency are sensitive to the reflectance of the internal coating, the internal geometry of the sphere, and the reflectance of the sample placed at the sample port. Hence, corrections reflecting those variations in the sphere are included in our model. The model we used is slightly different from the one used in the IAD algorithm. IAD considers the diameter of the exit port of the transmittance sphere to be set to zero, i.e., it replaces the reflectance of the 99% reflectance standard by the reflectance value of the sphere and, therefore, neglects any effect of the exit port of the transmittance sphere. This is inconsistent with the way the transmittance sphere is calibrated as explained in the IAD manual.12 The thick coating material present in the integrating spheres used in our setup makes it necessary to consider the presence of the exit port of the transmittance sphere. Moreover, the internal geometry of each of our integrating spheres is more complex than the one from the sphere model in the IAD routine. In our case, the additional reflectance surface coming from the baffle and a sample port diameter reducer has to be taken into account. However, for its simplicity, we still use the sphere model from the IAD routine in our computations only after scaling the diameter of the real integrating spheres to the one of a model sphere presenting the same value of the reflectance coefficient of the spheres’ walls. The details are described elsewhere.50
The AD method was used with a modified inversion routine to compute the optical properties of the samples. This modified routine is similar to the IAD procedure from Prahl2,11 but allows for some flexibility in the choice of the model for the integrating sphere. As is the case with IAD, this model takes the normal incident light as well as the light back-scattered by the spheres into account to compute the total hemispherical reflectance and total hemispherical transmittance of the sample.50 A comparison study of the optical properties obtained by analyzing the same set of simulated data by IAD and our inversion routine showed an agreement better than 2% (better than 5% for high-scattering or high-absorbing samples). Using the Guide to the Expression of Uncertainty in Measurement,52 our method also enables the computation of the uncertainty budget for the measurements of the optical properties of the sample by considering the statistical uncertainty (type A uncertainty) on the measurements and the nonstatistical uncertainty (type B uncertainty) on the sample’s parameters, on the integrating spheres’ parameters, on the detectors’ reflectance, and on the calibration standard’s reflectance. The uncertainty propagation is based on the numerical estimation of the sensitivity of the optical parameters of the sample to each experimental parameter and measured signal in order to compute a covariance matrix for the input quantities of the problem. Details on this procedure can be found elsewhere.50 A GUI was developed to facilitate the data analysis. The refractive index of the base material of the phantoms and the scattering anisotropy factor used in the data analysis comes from the values measured by INO at 660 nm, i.e., , and .19
Discussion and Conclusion
Figures 2 and 3 present the optical properties and measured at and for three 100-mm-square INO phantoms of thicknesses , 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 -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 () 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 .19 Figures 2 and 3 present the results with error bars (i.e., coverage factor ). At , our results are consistent with the values from INO for both and . At , the agreement is only for for both and . However, with a coverage factor , preferred by NIST,53 the optical properties measured here agree with the results provided by INO. Table 1 summarizes the results and the uncertainty with .
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
|Thickness (mm)||Parameter||Value (mm−1)||Absolute uncertainty (mm−1)||Relative uncertainty (%)|
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, . 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.54 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, , of the incident beam is , 30 deg, 60 deg, and 90 deg. Figure 4 presents the optical parameters, and , versus measured for one of the INO phantoms (thickness ) at . The values present fluctuations well within the error bars at , so we can rule out any measurable influence of the linear input polarization on the optical parameters measured with our DIS setup.
In this paper, we presented measurements of the optical properties of solid biomedical phantoms at and 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 , 7.09 mm, and 9.92 mm. There is an agreement within the error bars at between our results and the ones from INO at both and and for all thicknesses. The uncertainty in the optical parameters measured was estimated to be between 4.0% and 5.0% for and between 11% and 12% for , compared to 12% and 7%, respectively, obtained by INO.19 With the DIS measurements of the optical properties, the uncertainties on were found to be mainly sensitive to the uncertainty on the scattering anisotropy factor, , and the uncertainty on was mainly dependent on the uncertainty on the incident angle.50 As previously mentioned, the time-domain technique is less sensitive to measurements of , and consequently, the uncertainty on is smaller. Figures 2(a) and 3(a) show values that are thickness independent, whereas in Figs. 2(b) and 3(b), 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.50 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.
Certain commercial materials and equipment are identified in order to adequately specify the experimental procedure. Such identification does not imply recommendation by the National Institute of Standards and Technology. This work was supported by the NIST Innovation in Measurement Science (IMS) program on “Optical Medical Imaging for Clinical Applications” and the NIST intramural research program in the Statistical Engineering Division of the Information Technology Laboratory of NIST.