Reference optical phantoms for diffuse optical spectroscopy. Part 1 – Error analysis of a time resolved transmittance characterization method

doi:10.1364/OE.18.011495

Reference optical phantoms for diffuse optical spectroscopy. Part 1 – Error analysis of a time resolved transmittance characterization method

Jean-Pierre Bouchard, Israël Veilleux, Rym Jedidi, Isabelle Noiseux, Michel Fortin, and Ozzy Mermut »View Author Affiliations

Optics Express, Vol. 18, Issue 11, pp. 11495-11507 (2010)
http://dx.doi.org/10.1364/OE.18.011495

View Full Text Article

Acrobat PDF (1187 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Full Text
Article Info
References (25)
Cited By
Figures (7)
Metrics

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⁻¹ (11.3%) and 0.67 cm⁻¹ (6.8%) are obtained for the absorption coefficient and reduced scattering coefficient respectively.

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

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]

]. 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

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]

–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]

], interest in the absolute accuracy of phantoms optical properties is more recent [8

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

,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]

–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]

]. 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 (

μ_{a}

) and 41% for the scattering coefficient (

μ_{s}^{'}

) for a given phantom [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]

]. 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

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]

,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]

]. 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

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]

] 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

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]

]. 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

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]

], 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.

Fig. 1 Sensitivity of a transmittance TPSF trace of a 2 cm thick slab to a perturbation of 1% in the optical properties.

Y_{0}

(blue) is the reference TPSF.

Y_{δ μ_{a}}

and

Y_{δ μ_{s}^{'}}

refers to the perturbed TPSF traces. Difference between the perturbed TPSF and the reference TPSF are plotted in red and green for a

μ_{a}

and a

μ_{s}^{'}

perturbations respectively.

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

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

]: (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

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”.

]. 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

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

]. 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.

Fig. 2 Experimental setup for time resolved transmittance measurements.

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

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

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

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]

], 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

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

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

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

]. 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

\hat{u}

instead of

\hat{μ}

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 = (x, y, z)

, propagation direction

\hat{u}

and current step size s, boundary crossing event are detected when the following condition is met:

\begin{array}{l} C y l i n d r i c a l {(x + s u_{x})}^{2} + {(y + s u_{y})}^{2} > R^{2}, \\ Re c \tan g u l a r | x + s u_{x} | > X or | y + s u_{y} | > Y, \end{array}

(1)

where 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 (μ_{a}, μ_{s}^{'}, g; ρ, t)

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

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

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]

] 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

μ_{s}

and

μ_{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_{μ_{a} = 0} (\hat{r}, \hat{μ}, t)

is the response of absorption-free media to a temporal Dirac delta function

μ_{s}

excitation, then the relation:

L_{μ_{a}} (\hat{r}, \hat{u}, t) = L_{μ_{a} = 0} (\hat{r}, \hat{u}, t) \exp (- μ_{a} v t),

(2)

can be used to compute the case of a uniform non zero absorption coefficient [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]

]. This property can be understood very intuitively. Since all photons travel at the same speed v, at any given time t they all have travelled the same distance

v t

and have therefore experienced the same absorption factor

\exp (- μ_{a} v t)

TPSF were calculated for

μ_{s}

ranging from 9 cm⁻¹ to 74 cm⁻¹ and tabulated into a reference database. The diffusion approximation was used to determine the required step size between successive

μ_{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.

Fig. 3 TPSF database for a 20 mm thick cylindrical sample (D = 55 mm).

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:

M (t) = \int_{ω} \int_{A} L [\hat{r} = (ρ, θ, z = d), \hat{u}, t] d A d Ω,

(3)

where ω 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

{IRF}_{δ} (t)

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 (t) * I R F_{δ} (t)

. 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

τ_{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:

{m |}_{μ_{a}, μ_{s}^{'}} = G \cdot M ([0, Δ t, 2 Δ t, \dots]) * I R F .

(4)

Convolution with the IRF introduces the effect of the finite response time of the TCSPC system and also translates the modeled TPSF to the correct position on the TCSPC system time axis. Convolution with the IRF also has the added benefit of smoothing out the statistical variations of the Monte-Carlo model. A gain factor 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

(μ_{a}, μ_{s}^{'}, G)

to optimize the goodness-of-fit, criteria

χ^{2} = \sum_{n} {(y_{c n} - m_{n})}^{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.

Fig. 4 Fit results for cylindrical samples coming from 4 separate batches B0052, B0053, B0054 and B0055 of our phantom production samples collection [19

]. Note that the residuals have been magnified by a factor of 10.

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

Δ μ_{a} = 0.0006 {cm}^{- 1} (0.7 %)

and

Δ μ_{s}^{'} = 0.027 {cm}^{- 1} (0.3 %)

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

Δ μ_{a} = 0.0017 {cm}^{- 1} (1.8 %)

for the absorption coefficient and

Δ μ_{s}^{'} = 0.12 {cm}^{- 1} (1.2 %)

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

Δ μ_{a} = 0.0005 {cm}^{- 1} (0.62 %)

for the absorption coefficient and

Δ μ_{s}^{'} = 0.04 {cm}^{- 1} (0.4 %)

for the reduced scattering coefficient. These error estimates are conservative since multiple samples can be measured within a 10 minutes time frame.

Fig. 5 Instrument response functions acquired over a 4 hr time period

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

Δ d = 300 μ 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 μ 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

Δ μ_{a} = 0.001 {cm}^{- 1} (1.1 %)

and

Δ μ_{s}^{'} = 0.25 {cm}^{- 1} (2.5 %)

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

Δ μ_{a} = 0.001 {cm}^{- 1} (1.1 %)

and

Δ μ_{s}^{'} = 0.035 {cm}^{- 1} (0.35 %)

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₂ 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 - μ_{s}^{'} / μ_{s} \approx 1 - μ_{s}^{'} / μ_{t}

which neglects the small contribution of absorption to the total attenuation

μ_{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.

μ_{s}^{'}

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₂ concentration was

g = 0.62 \pm 0.015

. The uncertainty on g was calculated using the

μ_{s}^{'}

inaccuracy value obtained in this paper (see section 5) and the

μ_{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

Δ μ_{a} = 0.0001 {cm}^{- 1} (0.1 %)

and

Δ μ_{s}^{'} = 0.003 {cm}^{- 1} (0.03 %)

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

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

]. 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^{'} = 1.01 \cdot t

, for computing the theoretical TPSF from the Monte-Carlo model. The stretching of the axis resulted in offsets of

Δ μ_{s}^{'} = 0.07 {cm}^{- 1} (0.7 %)

and

Δ μ_{a} = 0.0015 {cm}^{- 1} (1.7 %)

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⁻¹) but different absorption coefficients (approximately 0.07 cm⁻¹ and 0.16 cm⁻¹) 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.

Fig. 6 Phantom set fabricated from a single batch to evaluate the model limitations

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

μ_{s}^{'} = 10 {cm}^{-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 .

Table 1 Characterization results for all phantoms

			Low absorption batch		High absorption batch
Geometry	Size	Thickness	$μ_{a}$	$μ_{s}^{'}$	$μ_{a}$	$μ_{s}^{'}$
	mm	mm	cm⁻¹	cm⁻¹	cm⁻¹	cm⁻¹
Cyl	55	10	0.060	9.08	0.146	8.82
Cyl	55	20	0.069	9.28	0.161	9.15
Cyl	55	30	0.071	9.49	0.164	9.32
Rec	30	20	0.072	9.50	0.162	9.33
Rec	50	20	0.070	9.45	0.162	9.35
Rec	80	20	0.070	9.46	0.160	9.20
		Mean	0.069	9.38	0.159	9.20
		Std. dev.	0.0042	0.17	0.0064	0.20
		% variation	6.2%	1.8%	4.0%	2.2%

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

μ_{a}

results of the low absorption set.

Fig. 7 μ _a and μ _s’ dependence on the sample’s lateral dimension and thickness

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

{Δ μ_{a} |}_{2 σ} = 0.01 {cm}^{-1} (11.3 %)

and

{Δ μ_{s}^{'} |}_{2 σ} = 0.67 {cm}^{-1} (6.8 %)

. 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

{Δ μ_{a} |}_{2 σ rnd
                            
                           only} = 0.0017 {cm}^{-1} (2 %)

and

{Δ μ_{s}^{'} |}_{2 σ rnd
                            
                           only} = 0.13 {cm}^{-1} (1.3 %)

. 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.

Table 2 Contribution of errors

Error source	Absolute error (cm⁻¹)		Relative error (%)*
Error source	$μ_{a}$	$μ_{s}^{'}$	$μ_{a}$	$μ_{s}^{'}$
Random error sources
Measurement noise	0.0006	0.027	0.71%	0.27%
System repeatability	0.0017	0.12	1.8%	1.2%
IRF instabilities	0.0005	0.04	0.59%	0.4%
Systematic error sources
Sample thickness	0.001	0.25	1.18%	2.5%
Refractive index	0.001	0.035	1.18%	0.35%
Anisotropy factor	0.0001	0.003	0.1%	0.03%
Time base inaccuracies	0.0015	0.07	1.76%	0.7%
Monte-Carlo model	0.004**	0.20	4.71%	2.2%
Total error (RSS, 1σ)	0.005	0.34	5.7%	3.4%
Total error (RSS, 2σ)	0.010	0.67	11.3%	6.8%

* relative to µ_a = 0.085 cm⁻¹ and µ’_s = 10 cm⁻¹ ** The standard deviation of µ_a obtained from the low absorption series of the sample geometry experiment (see Table 1) was selected because the mean absorption coefficient (0.069 cm⁻¹) was closer to the mean absorption coefficient of the other samples that were used for estimating the other error contributions (0.085 cm⁻¹).

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).

OCIS Codes
(120.3890) Instrumentation, measurement, and metrology : Medical optics instrumentation
(170.6510) Medical optics and biotechnology : Spectroscopy, tissue diagnostics

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: March 10, 2010
Revised Manuscript: April 22, 2010
Manuscript Accepted: May 6, 2010
Published: May 14, 2010

Virtual Issues
Vol. 5, Iss. 10 Virtual Journal for Biomedical Optics

Citation
Jean-Pierre Bouchard, Israël Veilleux, Rym Jedidi, Isabelle Noiseux, Michel Fortin, and Ozzy Mermut, "Reference optical phantoms for diffuse optical spectroscopy. Part 1 – Error analysis of a time resolved transmittance characterization method," Opt. Express 18, 11495-11507 (2010)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-18-11-11495

Sort: Author | Year | Journal | Reset

References

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]
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]
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]
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]
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]
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]
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]
E. Alerstam, S. Andersson-Engels, and T. Svensson, “Improved accuracy in time-resolved diffuse reflectance spectroscopy,” Opt. Express 15, 10434–10448 (2007).
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]
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]
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).
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]
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]
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]
E. Alerstam, S. Andersson-Engels, and T. Svensson, “Improved accuracy in time-resolved reflectance spectroscopy,” Opt. Express 16(14), 10440–10448 (2008). [CrossRef] [PubMed]
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]
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]
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]
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”.
W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in FORTRAN 77, (Cambridge, 1992), Chap. 15.
W. Becker, The bh TCSPC Handbook, Third Edition (Becker & Hickl GmbH, 2008)
L. V. Wang, Biomedical Optocs, Principles and Imaging (Wiley, 2007), Chap. 5.
“(MCML) Monte Carlo for Multi-Layered media, ” http://omlc.ogi.edu/software/mc/
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]
W. Becker, Becker & Hickl, Nahmitzer Damm 30, 12277 Berlin, (personal communication, 2008).

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper