## Abstract

Development, production quality control and calibration of optical tissue-mimicking phantoms require a convenient and robust characterization method with known absolute accuracy. We present a solid phantom characterization technique based on time resolved transmittance measurement of light through a relatively small phantom sample. The small size of the sample enables characterization of every material batch produced in a routine phantoms production. Time resolved transmittance data are pre-processed to correct for dark noise, sample thickness and instrument response function. Pre-processed data are then compared to a forward model based on the radiative transfer equation solved through Monte Carlo simulations accurately taking into account the finite geometry of the sample. The computational burden of the Monte-Carlo technique was alleviated by building a lookup table of pre-computed results and using interpolation to obtain modeled transmittance traces at intermediate values of the optical properties. Near perfect fit residuals are obtained with a fit window using all data above 1% of the maximum value of the time resolved transmittance trace. Absolute accuracy of the method is estimated through a thorough error analysis which takes into account the following contributions: measurement noise, system repeatability, instrument response function stability, sample thickness variation refractive index inaccuracy, time correlated single photon counting system time based inaccuracy and forward model inaccuracy. Two sigma absolute error estimates of 0.01 cm^{−1} (11.3%) and 0.67 cm^{−1} (6.8%) are obtained for the absorption coefficient and reduced scattering coefficient respectively.

© 2010 OSA

## 1. Introduction

Reference optical tissue-mimicking phantoms exhibiting stable and accurately characterized properties are a mandatory tool for the development, validation, calibration and quality control of any biomedical spectroscopic or imaging device [1]. In the development stage, phantoms are repeatedly used to test, debug and optimize the system being built. Multi-center clinical trials require properly referencing and characterization of the instruments among all the sites to ensure data consistency, enable cross-comparison and validation, before adoption of the new technology. Variability in tissue measurements on living subjects is often problematic in most biomedical diagnostic application. This variability is inherent to the population of subjects for which the application is developed and therefore cannot be avoided. On the other hand, instrument to instrument measurements variability, under ideal conditions, can and must be minimized to the highest possible extent. Sharing a “golden” phantom set can effectively reduce this variability but is not an interesting approach in the long term. In this situation, long term consistency of the results produced by a diagnostic instrument would rely on the stability and integrity of the “golden” phantom set used to calibrate every instrument in use. A practical primary characterization method that can determine the optical properties of a phantom with a known degree of accuracy and without comparison to a “golden” reference is preferred. The technique should ideally be rapid and cost effective to be compatible within the process of phantoms production. For example, it would ideally be used as a quality assurance/control step to deliver a calibration certificate for every unit produced, therefore ensuring traceability.

While characterization techniques have been mostly developed in the 1990s [2–7], interest in the absolute accuracy of phantoms optical properties is more recent [8,10–15]. In 2005, results from the application of the MEDPHOT performance assessment protocol on eight different instruments have shown inter-system variation up to 32% for absorption coefficient ( ${\mu}_{a}$ ) and 41% for the scattering coefficient ( ${\mu}_{s}^{\prime}$ ) for a given phantom [13]. These results revealed the need for better assessment of the absolute accuracy of the measurement methods used to evaluate optical properties. Recent studies concerning characterization of liquid phantoms have shown that the error on an increase of absorption and scattering coefficients can be reduced to 2% when the absorbing or scattering ingredients can be incrementally added to the phantom without modifying the measurement setup [12,14]. Unfortunately, liquid phantoms are not stable over a long term, and permanent accessibility to the calibration apparatus is required to work with them, a situation that is clearly not desirable on a production floor or in a clinical setting.

Solid tissue simulating phantoms [16,17] offer the benefits of ease-of-use and long term stability in their optical characteristics. They are therefore the preferred option for instrument response standardization [1]. In this communication, we present a time resolved transmittance characterization technique suitable for solid tissue phantoms that is compatible with a volume production environment. A detailed description of the experimental characterization apparatus and accompanying data analysis is first presented. A thorough evaluation of the absolute accuracy of phantoms is then developed taking into account the random and systematic error sources.

## 2. Time resolved transmittance characterization method

The absorption coefficient and the reduced scattering coefficient cannot independently be measured in a turbid medium. A model-based parameter extraction must be employed which consist of iteratively adjusting unknown optical properties of the light propagation model until it matches the measurements. Measurement of time resolved transmittance of light pulses through the phantoms sample has been selected as our preferred characterization method. Variations in absorption and scattering coefficient have distinct effects on the measured Time Point Spread Function (TPSF) traces (see Fig. 1 ) and are therefore easily decoupled. The availability of broad spectrum supercontinuum sources now makes it possible to use the technique over a very large continuous wavelength range with high optical power. The technique produces an information rich measurement vector. Although a perfect TPSF trace has been shown to contain limited information content that can be summarized by 2 parameters [18], an experimental measurement trace can reveal symptoms of problems in the characterization setup (e.g. flaws such as ghost reflection and lack of proper light isolation between the source and the detection system). The choice of the transmittance geometry is also motivated by accuracy. Transmittance geometry makes it easier to isolate unwanted light to reach the detector. In diffuse reflectance geometry, light can escape the sample close to the source fiber, reflect off surfaces, re-enter the sample close to detection, and contribute to the measured trace in a substantial way. Light baffling can eliminate this interference but is not a trivial solution when the highest level of accuracy is desired. Finally, the technique’s only calibration step is a measurement of the instrument response function.

Optical properties are obtained by adjusting a theoretical model to the experimental data until an optimal fit is achieved. A least squares fit between the experimental data and the theoretical model gives the most likely values for the sought after optical properties. The previous statement is true only if the following two working assumptions are valid [20]: (1) the random noise of the experiment follows a normal distribution and (2) the model accurately represents the underlying physics of the measurements. Time Correlated Single Photon Counting (TCSPC) measurements are affected by noise that follows a Poisson distribution. However, the Poisson distribution can be approximated by a normal distribution with negligible error for photon counts over 100, which is always the case in our measurements; therefore the first assumption for unbiased retrieval of optical properties is valid. The accurate model assumption is much more challenging and to meet it requires careful experiment design, an accurate forward model, and appropriate data pre-treatment. The remainder of this section describes the care taken in these three aspects to ensure bias-free determination of optical properties.

#### 2.1 Sample size

Each individual biomedical application and/or instrument requires its own phantom size and geometry. A mouse or head-shaped phantom will be harder to accurately characterize than a large uniform slab. This is further complicated by the cumbersome need to adapt a characterization setup for every imaginable phantoms geometry that may be produced. Phantoms material is fabricated by batches that are poured in a mold of defined size and desired geometry [19]. The volume that can be poured in a single mold is limited by the exothermic reaction of polymerization that occurs in the mold once the material has been poured. Hence, the fabrication of large phantoms (e.g. head-shaped phantom) usually requires multiple batches. From the above statement, a batch characterization method that samples a small portion of each batch for characterization referencing and quality control purposes is highly desirable. We selected a cylindrical sample shape with a 20 mm thickness and 55 mm diameter. The 20 mm thickness allows for sufficient spreading of the light pulse even for low scattering coefficient values. The phantom diameter is chosen such that the lateral boundaries are further away than the source detection separation therefore limiting the boundary effect, although the model we used for optical properties takes into account the presence of the lateral boundary (see section 2.3).

#### 2.2 Experimental setup

Figure 2 . shows the experimental setup used to measure the time resolved transmittance of the phantoms. A super-continuum laser (SC400, Fianium, UK) generates ~90 ps light pulses at a repetition frequency of 40 MHz. The laser light is filtered to obtain a 660 nm ± 5 nm beam. A small fraction of the light is sent to a synchronization detector (PHD-400, Becker & Hickl, Germany) by means of a microscope slide acting as a beam splitter. The remaining narrow collimated beam is normally incident on the phantom. Light exiting the phantom on the opposite side is collected with a photon counting micro channel plate photomultiplier tube (R3809U, Hamamatsu, Japan) located 8 cm from the exit surface. The signal from the PMT and the synchronization detector are sent to the TCSPC computer board (SPC-730, Becker & Hickl, Germany). The TCSPC system returns 4096 measurement points with a temporal resolution of 6.1 ps. The optical signal is attenuated to maintain a count rate of approximately 200 kHz. This value has been adjusted to avoid the broadening effect that a high count rate has on MCP-PMT Instrument Response Function (IRF) [21]. The IRF is measured with no sample and a piece of thin (< 50 µm) translucent adhesive tape to diffuse the light over the total area of the PMT. Uniform illumination of the detector surface is crucial as the MCP-PMT has a position dependent IRF.

#### 2.3 Numerical modeling of light transport through the sample

The optical properties are defined in the framework of the widely accepted radiative transfer equations (RTE), which makes it the only possible choice for our light propagation model. A derivation of the radiative transfer equation based modeling of light transport in turbid medium can be found in previous literature [22].

Diffusion approximation (DA) of the RTE is often used to extract optical properties from characterization measurements. Use of the DA should be avoided if possible when accuracy is a concern. Even though the validity of the DA has been extensively studied and proven to be a very good approximation of the RTE [2,24], it is still an approximation and thus can impair the accuracy of the retrieved results by introducing an unnecessary bias. Alerstam et al. have recently demonstrated that using the DA can lead to errors on the order of 20% in the obtained optical properties when compared to those determined using a RTE based model [8].

#### 2.4 Monte Carlo solution of the RTE

Solution of the RTE in the finite geometry of our phantom test samples can only be achieved through the use of the Monte Carlo method. The Monte Carlo simulation used in this study is based on Wang and Jacques’s well known MCML program [23]. This MCML code is well-documented and only modifications made to it will be described here. We implemented two additions to the original MCML code in order to represent the geometry of our samples and the time-dependence of our measurement approach. Notation used herein is consistent with MCML documentation with the exception of the direction vector, which is noted $\widehat{u}$ instead of $\widehat{\mu}$ to avoid confusion with the optical properties.

The first modification is the addition of the appropriate boundary conditions to model the finite size of our phantoms. Two cases have been implemented, cylindrical- and rectangular- shaped phantoms. Given the current photon packet position $r=\left(x,y,z\right)$ , propagation direction $\widehat{u}$ and current step size *s*, boundary crossing event are detected when the following condition is met:

*R*,

*X*and

*Y*represents the position of the boundary in cylindrical or rectangular coordinates (radius or half-length of the phantom). If the conditions are such that a photon packet would cross the boundary, the intersection position and an updated direction vector are computed. The photon packet is then reduced in weight according to Fresnel formulas and propagated with the updated direction.

The second modification consists of adding time resolved capability through registration of the photon packet total time of flight. In the original code, photons are collected into annular bins centered at $r=0$ as they exit the sample. Our modification divides each annular bin into temporal bins in order to retrieve the time point spread function (TPSF) for each radial position. The output of the modified MCML code is noted $T\left({\mu}_{a},{\mu}_{s}^{\prime},g;\rho ,t\right)$ .

#### 2.5 Speeding up Monte Carlo simulations

Time resolved Monte-Carlo simulations of light transport with small statistical variation are very computationally intensive and thus impede the use of full RTE modeling in the calibration phantoms characterization process. Acceleration strategies that exploit scaling properties of the RTE [8,9] cannot be exploited with a finite geometry phantom. Two strategies have been used to overcome this difficulty in routine phantoms characterization. The first one is to use a standardized test sample size [19]. Always using the same phantom geometry facilitates efficiency by running the calculation only once for a particular parameter set and saving the results for further use. Thereafter, a database of calculated responses as a function of optical properties can be constructed, and tabulated (a so-called “look-up” table). In this study, *g* was fixed at the constant value of 0.62 [19] leaving ${\mu}_{s}$ and ${\mu}_{a}$ as the free parameters over which to tabulate the results. Even computing a two-dimensional grid of TPSF is still prohibitively long. This last difficulty can be solved by exploiting a very interesting property of the RTE. If ${L}_{{\mu}_{a}=0}\left(\widehat{r},\widehat{\mu},t\right)$ is the response of absorption-free media to a temporal Dirac delta function ${\mu}_{s}$ excitation, then the relation:

*v*, at any given time

*t*they all have travelled the same distance $vt$ and have therefore experienced the same absorption factor $\mathrm{exp}\left(-{\mu}_{a}vt\right)$ .

TPSF were calculated for ${\mu}_{s}$ ranging from 9 cm^{−1} to 74 cm^{−1} and tabulated into a reference database. The diffusion approximation was used to determine the required step size between successive ${\mu}_{s}$ values. The criterion was that linear interpolation between the two computed TPSF shall not introduce an error greater than 1%. An average of $120\cdot {10}^{6}$ photon packets were launched for each case. Figure 3
shows the entire TPSF database for a cylindrical phantom having a 55 mm diameter and a thickness of 20mm.

#### 2.6 Data pre-treatment and analysis

When measuring with a sample, we are in fact measuring the following function defined in the model time reference:

*ω*denote the acceptance solid angle of the detector and

*A*is the exposed output surface of the sample.

A first data pre-treatment step consists of correcting for the DC dark count level that is not taken into account in the model. The dark count level is first estimated by taking the average of 100 data points taken in a portion of the TPSF trace where only noise is present. This dark count level is then subtracted from the raw measurement vector.

Another very important correction is needed to reflect the temporal axis differences between the IRF measurement and the sample measurement. We first illustrate this correction with an idealized experiment where the IRF is a Dirac delta function. Time zero is defined in the model as the time photons are incident on the surface of the sample. To determine the position of this origin on the time axis of the TCSPC system, the ideal experiment would be to position the detector in the plane of the first surface of the sample and measure the ${\text{IRF}}_{\delta}\left(t\right)$ trace. The detector would then have to be moved back by a distance equal to the sample thickness *d* to be located on the output surface of the sample and to measure the experimental trace $M\left(t\right)\ast IR{F}_{\delta}\left(t\right)$ . This method would yield a correct referencing of the time axis for the sample measurement. In practice, we prefer to leave the detector fixed to minimize the manipulation steps. The required sample thickness dependant retardation ${\tau}_{s}=d/c$ was thus applied by shifting the measured TPSF in time through data interpolation. The resulting pre-conditioned measurement vector corrected for the baseline noise and sample thickness is defined as ${y}_{c}$ .

The measurement vector has to be compared to a modeled vector **m** to extract optical properties. The model vector was obtained by the following expression:

*G*is introduced to account for the measured arbitrary amplitude output by the system. Fitting of

**m**to the measurement vector ${y}_{c}$ was achieved by varying the three floating parameters $\left({\mu}_{a},{\mu}_{s}^{\prime},G\right)$ to optimize the goodness-of-fit, criteria ${\chi}^{2}={\displaystyle {\sum}_{n}{\left({y}_{cn}-{m}_{n}\right)}^{2}}$ , using the Levenberg-Marquardt algorithm. A very stringent fit window, defined as 1% of signal peak for both the rising edge and the falling edge of the trace, was selected. Figure 4 . shows typical fit results obtained with our characterization procedure. The RMS value of the residuals (shown magnified by a factor of 10) is less than 0.3% of the maximum amplitude of the signal.

## 3. Error analysis: Sources of random errors

#### 3.1 Measurement noise

Short term random fluctuations in a measured TCSPC trace come from the shot noise associated with the optical signal itself and the dark counts. Measurement noise can be modeled but an experimental determination is more straightforward and convincing. The effect of measurement noise was therefore estimated by measuring 4 TPSF traces in sequence for six different samples. The averaged standard deviation of the fitted optical properties for each sample was $\Delta {\mu}_{a}=0.0006\text{\hspace{0.17em}}{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(0.7\%\right)$ and $\Delta {\mu}_{s}^{\prime}=0.027\text{\hspace{0.17em}}{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(0.3\%\right)$ .

#### 3.2 System repeatability

A general system repeatability experiment was conducted to quantify variation that can occur when the complete measurement sequence is repeated from powering up to final sample TPSF measurement. The following sequence has been repeated five times for a single sample:

- a) power up the system and wait five minutes for warm-up,
- b) measure the IRF,
- c) insert the sample in the sample holder,
- d) measure the sample TPSF,
- e) repeat step c) and d) three times randomly rotating the sample each time,
- f) shut down the system.

The standard deviation obtained from all 15 measurements (five repetitions times three random orientations) is $\Delta {\mu}_{a}=0.0017{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(1.8\%\right)$ for the absorption coefficient and $\Delta {\mu}_{s}^{\prime}=0.12{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(1.2\%\right)$ for the reduced scattering coefficient.

#### 3.3Instrument response function (IRF) instability

The system response can drift between the IRF measurement and the sample measurement. To minimize the error introduced by this instability, a measurement of the IRF is performed at the beginning of the phantom characterization session. An upper bound to the contribution of the IRF instability to the total uncertainty has been determined by analyzing a single TPSF trace with 20 instrument response functions acquired over a time period of 4 hrs (see Fig. 5 ). The IRF instability contribution to the total error has been determined by taking the standard deviation of the 20 retrieved optical properties. This procedure was repeated for TPSF traces measured on 4 different phantoms samples. The average of the standards deviations obtained from the 4 repetition gives an IRF instability contribution of $\Delta {\mu}_{a}=0.0005{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(0.62\%\right)$ for the absorption coefficient and $\Delta {\mu}_{s}^{\prime}=0.04{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(0.4\%\right)$ for the reduced scattering coefficient. These error estimates are conservative since multiple samples can be measured within a 10 minutes time frame.

## 4. Error analysis: Sources of systematic errors

#### 4.1 Sample thickness inaccuracy

Our batch characterization samples are machined to a cylindrical shape with a certain degree of dimensional accuracy. By taking the standard deviation of the measured thickness of 80 standard samples (see section 2.1), we have determined the accuracy in thickness to be $\Delta d=300\mu m$ . Our data analysis is based on pre-computed TPSF traces from a Monte-Carlo model that assumes a nominal sample thickness of 20 mm. The sample thickness variation, even though it can be characterized to a high level of accuracy, is not taken into account in our analysis and induces a bias on the optical properties. The magnitude of this bias was estimated by fitting an experimental trace using both a correct (case 1) and an erroneous (case 2) sample thickness, ${d}_{0}$ and ${d}_{0}+300\text{\hspace{0.05em}}\mu m$ respectively, and determining the difference in the obtained optical properties. A diffusion approximation, DA, slab model [7] was used to allow free variation of the sample thickness. The error estimate therefore assumes that the bias induced by the use of the DA is common to both cases and thus subtracts out when taking the difference in retrieved optical properties. The process was repeated with the 4 experimental traces shown in Fig. 4. The average of the bias values obtained for the 4 traces gives sample thickness inaccuracy bias estimates of $\Delta {\mu}_{a}=0.001{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(1.1\%\right)$ and $\Delta {\mu}_{s}^{\prime}=0.25{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(2.5\%\right)$ .

#### 4.2 Refractive index inaccuracy

Error in the evaluation of the refractive index of the bulk sample has a direct impact on the recovered optical properties. As described in [19], the refractive index of the polyurethane used for phantoms fabrication was determined by a time-of-flight experiment. The value obtained for the refractive index was $n=1.521\pm 0.006$ .

The impact of this refractive index uncertainty on the retrieved optical properties uncertainty was evaluated in a similar fashion for the sample thickness inaccuracy. 4 experimental traces were analyzed with both the correct and the erroneous values of the refractive index. Retrieved optical properties for each index value were then subtracted to obtain the biases induced by the refractive index inaccuracy. The average of the bias values obtained for the 4 traces gives a refractive index inaccuracy bias estimate of $\Delta {\mu}_{a}=0.001{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(1.1\%\right)$ and $\Delta {\mu}_{s}^{\prime}=0.035{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(0.35\%\right)$ .

#### 4.3 Anisotropy factor inaccuracy

Error in the evaluation of the anisotropy factor *g* may also impact the recovered optical properties. The *g* factor used in the Monte-Carlo model was determined experimentally as described in [19]. In brief, phantom batches with TiO_{2} particles but no absorber were prepared and machined into thin wedges in addition to our standard characterization cylinders. The thickness of the wedged samples was selected to insure single scattering regime in transmission. The anisotropy factor was calculated using $g=1-{\mu}_{s}^{\prime}/{\mu}_{s}\approx 1-{\mu}_{s}^{\prime}/{\mu}_{t}$ which neglects the small contribution of absorption to the total attenuation ${\mu}_{t}$ because no absorber was used. The total attenuation coefficient of a given batch was determined by measuring the coherent (non-scattered) transmission of the thin wedges. Samples were wedged to allow measurements at differential thicknesses on the same sample to experimentally factor out the contribution of Fresnel reflection at the interfaces. ${\mu}_{s}^{\prime}$ was determined by characterizing the standard cylinders coming from the same batch using the technique described in section 2 assuming a *g* value of 0.59. The mean value of the anisotropy factor obtained for the various TiO_{2} concentration was $g=0.62\pm 0.015$ . The uncertainty on *g* was calculated using the ${\mu}_{s}^{\prime}$ inaccuracy value obtained in this paper (see section 5) and the ${\mu}_{t}$ standard deviation observed experimentally.

To evaluate the possible impact of this uncertainty on the retrieved optical properties, Monte-Carlo simulations were used to generate traces for *g* values excursion of 0.015. These traces were then treated like experimental input vectors to recover the optical properties using our reference database (which assumes a *g* value of 0.62). The dependence on the *g* value was found to be relatively weak. The average of the bias values are of $\Delta {\mu}_{a}=0.0001{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(0.1\%\right)$ and $\Delta {\mu}_{s}^{\prime}=0.003{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(0.03\%\right)$ .

#### 4.3 Time base inaccuracy

The time axis of a TCSPC trace is determined by the system’s time-to-analog (TAC) converter. This electronics component (internal to the TCSPC system) measures the time between a photon detection event and a reference synchronization pulse by integrating a current source in a capacitor, therefore converting a time delay into a voltage that can be digitized into a numerical value. TCSPC system time bases are calibrated by the manufacturer by sending pulses with known delays to the CFD (detector) input and the synchronizing input using a delay generator. The calibration error in this time base is estimated to be 1% according to the manufacturer [25]. The implication of this unknown time stretching on the uncertainty in the retrieved optical properties has been estimated by using a stretched version of the time axis vector, ${t}^{\prime}=1.01\cdot t$ , for computing the theoretical TPSF from the Monte-Carlo model. The stretching of the axis resulted in offsets of $\Delta {\mu}_{s}^{\prime}=0.07{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(0.7\%\right)$ and $\Delta {\mu}_{a}=0.0015{\text{cm}}^{-1}\text{\hspace{0.17em}}\left(1.7\%\right)$ in the optical properties.

#### 4.4 Forward model inaccuracy

The chosen approach to evaluate the limitations of the model is to compare the recovered optical properties from phantoms of various shapes and sizes made from a single batch of polyurethane mix. Two complete phantom sets of equal reduced scattering coefficients (~10 cm^{−1}) but different absorption coefficients (approximately 0.07 cm^{−1} and 0.16 cm^{−1}) were casted into molds and machined into cylinder and rectangular blocks for a total of six different geometries (Fig. 6
.). More details about our phantom fabrication process including scatterer and absorber calibrations can be found in [19,17]. This phantom set allows evaluation of the model dependence with respect to two geometrical parameters: lateral size and thickness.

Each phantom was characterized using a Monte-Carlo model taking into account its specific geometry, as described in section 2.3. Proper modeling of the light escaping from the sides of the phantom gets particularly important as the lateral size is reduced. For example, results from our simulations are showing that given ${\mu}_{s}^{\prime}=10\text{\hspace{0.17em}}{\text{cm}}^{\text{-1}}$ , approximately 20% of the incoming light is lost in this manner for the 30 mm wide rectangular phantom. For each sample, three TPSF were acquired with the light beam normally incident on its center. Averages of the recovered optical properties are compiled in Table 1 .

As evidenced in Fig. 7 . (left graphs), minimal dependence on the lateral dimension of the rectangular samples is observed. This suggests that boundary conditions are properly modeled in the modified MCML code. However, the model appears not as robust with respect to phantom thickness. When plotted against sample thickness, the recovered absorption and scattering coefficients (Fig. 7., right graphs) show a clear trend that cannot be attributed to measurement noise. The highest relative variability, expressed as the standard deviation of the values over the mean, is observed for the ${\mu}_{a}$ results of the low absorption set.

## 5. Error analysis budget

Summarized in Table 2 are the estimated results of the random and systematic error sources evaluated in this work. Worst case values were selected for each contribution. A 2 sigma root mean square sum of all contributor gives uncertainties of ${\Delta {\mu}_{a}|}_{2\sigma}=0.01{\text{cm}}^{\text{-1}}\text{\hspace{0.17em}}\left(11.3\%\right)$ and ${\Delta {\mu}_{s}^{\prime}|}_{2\sigma}=0.67{\text{cm}}^{\text{-1}}\text{\hspace{0.17em}}\left(6.8\%\right)$ . It is to be noted that the precision of the technique (its ability to detect small relative changes in the optical properties), is only affected by the random fluctuations for which an RSS sum gives ${\Delta {\mu}_{a}|}_{2\sigma \text{\hspace{0.17em} rnd \hspace{0.17em} only}}=0.0017{\text{cm}}^{\text{-1}}\text{\hspace{0.17em}}\left(2\%\right)$ and ${\Delta {\mu}_{s}^{\prime}|}_{2\sigma \text{\hspace{0.17em} rnd \hspace{0.17em} only}}=0.13{\text{cm}}^{\text{-1}}\text{\hspace{0.17em}}\left(1.3\%\right)$ . The technique is therefore very sensitive to small changes in the optical properties. Highly accurate values of the ratio of optical properties could therefore be obtained since the systematic portion of the errors would cancel out. Ratios of absorption or reduces scattering coefficient could be calculated between two phantoms or between to values obtained at different wavelengths for the same phantom.

## 6. Conclusion

The primary objective here was to evaluate the absolute accuracy of determined optical properties from a solid phantom when using the time resolved transmittance measurement method. The random contribution of measurement noise and IRF instability could be reduced by averaging multiple measurements. The sample thickness contribution could be reduced by tighter manufacturing tolerances on the dimensions of the test sample. Generation of a 2-parameter TPSF database that can be interpolated for a measured thickness and scattering coefficient is also an effective strategy for reducing this contribution. Refinement to the refractive index measurements and time base accuracy evaluations can also be explored, but the largest error contributor remains the model inaccuracy. However, improvement on the robustness of the model can be non trivial. Bias imposed by boundary effects, the sample thickness correction, and the RTE modeling have already been taken into account. Investigation of the root cause of the remaining biases will be the focus of future efforts. Once those last issues are resolved, the proposed technique and data analysis presented herein could serve as a standard method to determine the optical properties of turbid tissue-mimicking media.

## References and links

**1. **B. W. Pogue and M. S. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt. **11**(4), 041102 (2006). [CrossRef] [PubMed]

**2. **F. Martelli, D. Contini, A. Taddeucci, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. II. Comparison with Monte Carlo results,” Appl. Opt. **36**(19), 4600–4612 (1997). [CrossRef] [PubMed]

**3. **S. A. Prahl, M. J. C. van Gemert, and A. J. Welch, “Determining the optical properties of turbid media by using the adding-doubling method,” Appl. Opt. **32**(4), 559–568 (1993). [CrossRef] [PubMed]

**4. **J. W. Pickering, S. A. Prahl, N. van Wieringen, J. F. Beek, H. J. C. M. Sterenborg, and M. J. C. van Gemert, “Double-integrating sphere system for measuring the optical properties of tissue,” Appl. Opt. **32**(4), 399–410 (1993). [CrossRef] [PubMed]

**5. **M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. **28**(12), 2331–2336 (1989). [CrossRef] [PubMed]

**6. **J. B. Fishkin, P. T. C. So, A. E. Cerissi, S. Fantini, M. A. Franceschini, and E. Gratton, “Frequency-domain method for measuring spectral properties in multiple-scattering media: methemoglobin absorption spectrum in a tissuelike phantom,” Appl. Opt. **34**(7), 1143–1155 (1995). [CrossRef] [PubMed]

**7. **D. Contini, F. Martelli, and G. Zaccanti, “Photon migration through a turbid slab described by a model based on diffusion approximation. I. Theory,” Appl. Opt. **36**(19), 4587–4599 (1997). [CrossRef] [PubMed]

**8. **E. Alerstam, S. Andersson-Engels, and T. Svensson, “Improved accuracy in time-resolved diffuse reflectance spectroscopy,” Opt. Express **15**, 10434–10448 (2007).

**9. **A. Kienle and M. S. Patterson, “Determination of the optical properties of turbid media from a single Monte Carlo simulation,” Phys. Med. Biol. **41**(10), 2221–2227 (1996). [CrossRef] [PubMed]

**10. **C. Chen, J. Q. Lu, H. Ding, K. M. Jacobs, Y. Du, and X.-H. Hu, “A primary method for determination of optical parameters of turbid samples and application to intralipid between 550 and 1630 nm,” Opt. Express **14**(16), 7420–7435 (2006). [CrossRef] [PubMed]

**11. **L. Spinelli, F. Martelli, A. Farina, A. Pifferi, A. Torricelli, R. Cubeddu, and G. Zaccanti, “Accuracy of the nonlinear fitting procedure for time-resolved measurements on diffusive phantoms at NIR wavelength,” Proc. SPIE **717424**, 1–10 (2009).

**12. **L. Spinelli, F. Martelli, A. Farina, A. Pifferi, A. Torricelli, R. Cubeddu, and G. Zaccanti, “Calibration of scattering and absorption properties of a liquid diffusive medium at NIR wavelengths. Time-resolved method,” Opt. Express **15**(11), 6589–6604 (2007). [CrossRef] [PubMed]

**13. **A. Pifferi, A. Torricelli, A. Bassi, P. Taroni, R. Cubeddu, H. Wabnitz, D. Grosenick, M. Möller, R. Macdonald, J. Swartling, T. Svensson, S. Andersson-Engels, R. L. P. van Veen, H. J. C. M. Sterenborg, J.-M. Tualle, H. L. Nghiem, S. Avrillier, M. Whelan, and H. Stamm, “Performance assessment of photon migration instruments: the MEDPHOT protocol,” Appl. Opt. **44**(11), 2104–2114 (2005). [CrossRef] [PubMed]

**14. **F. Martelli and G. Zaccanti, “Calibration of scattering and absorption properties of a liquid diffusive medium at NIR wavelengths. CW method,” Opt. Express **15**(2), 486–500 (2007). [CrossRef] [PubMed]

**15. **E. Alerstam, S. Andersson-Engels, and T. Svensson, “Improved accuracy in time-resolved reflectance spectroscopy,” Opt. Express **16**(14), 10440–10448 (2008). [CrossRef] [PubMed]

**16. **T. Moffitt, Y.-C. Chen, and S. A. Prahl, “Preparation and characterization of polyurethane optical phantoms,” J. Biomed. Opt. **11**(4), 041103 (2006). [CrossRef] [PubMed]

**17. **M. L. Vernon, J. Freàchette, Y. Painchaud, S. Caron, and P. Beaudry, “Fabrication and characterization of a solid polyurethane phantom for optical imaging through scattering media,” Appl. Opt. **38**(19), 4247–4251 (1999). [CrossRef]

**18. **A. R. Pineda, M. Schweiger, S. R. Arridge, and H. Barrett, “Information content of data types in time-domain optical tomography,” J. Opt. Soc. Am. A **23**(12), 2989–2996 (2006). [CrossRef]

**19. **J.-P. Bouchard, National Optics Institute, 2740 Einstein, Québec, Qc, G1P 4S4 are preparing a manuscript to be called “Reference optical phantoms for diffuse optical spectroscopy. Part 2 - Fabrication”.

**20. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in FORTRAN 77, (Cambridge, 1992), Chap. 15.

**21. **W. Becker, *The bh TCSPC Handbook, Third Edition* (Becker & Hickl GmbH, 2008)

**22. **L. V. Wang, Biomedical Optocs, Principles and Imaging (Wiley, 2007), Chap. 5.

**23. **“(MCML) Monte Carlo for Multi-Layered media, ” http://omlc.ogi.edu/software/mc/

**24. **F. Martelli, M. Bassani, L. Alianelli, L. Zangheri, and G. Zaccanti, “Accuracy of the diffusion equation to describe photon migration through an infinite medium: numerical and experimental investigation,” Phys. Med. Biol. **45**(5), 1359–1373 (2000). [CrossRef] [PubMed]

**25. **W. Becker, Becker & Hickl, Nahmitzer Damm 30, 12277 Berlin, (personal communication, 2008).