Open Access
1 January 2005 Simple peak shift analysis of time-of-flight data with a slow instrumental response function
Goro Nishimura, Mamoru Tamura
Author Affiliations +
Abstract
Analysis of time-of-flight (TOF) data is sometimes limited by the instrumental response function, and optical parameters are extracted from the observed response curve by several mathematical methods, such as deconvolution. In contrast to this, we demonstrate that a method using shifts of the peak time of the response curve with different source-detector separations can yield the average path length of the light traveling in a tissue-like sample without deconvolution. In addition, combining the intensity information allows us to separate the scattering and absorption coefficients. This simple method is more robust in signal-to-noise ratio than the moment analysis, which also does not require the deconvolution procedure, because the peak position is not significantly dependent on the baseline fluctuation and the contamination of the scattering. The analysis is demonstrated by TOF measurements of an Intralipid solution at 800 nm, and is applied to the measurements at 1.29 µm, where the temporal response of photomultiplier tubes is not sufficiently good.

1.

Introduction

Near-infrared (NIR) light enables us to characterize the physiological state of tissue in vivo because of its good penetration in tissues.1 One of the main problems of tissue spectroscopy at the NIR region is a multiple scattering process of light, which randomizes the optical path, resulting in difficulty in the quantification of optical parameters. The multiple scattered light propagates as an energy diffusion because the properties of the wave of light are almost completely lost by the multiple scattering process. This diffusive property of the propagation can be well characterized by macroscopic optical parameters, such as the scattering and absorption coefficients of a medium. Therefore, the measurement of light traveling in a tissue with a macroscopic scale (approximately centimeters), which is essentially this diffusive light, determines such optical parameters. To extract the parameters, the most direct method is time domain measurement of the light transport. The distribution of the path of the detected photons in a sample can be determined by the temporal response function (TRF), and thus the method is called the time-of-flight (TOF) method. The TRF can be described by the solution of a time-dependent optical diffusion equation, or more generally, by that of a time-dependent transport equation of light. Consequently, the parameters are determined by the solutions with an ideal model for the optical geometry and the physical structures of a practical sample. This method is most successful in characterizing the tissue’s optical parameters,2 3 4 5 6 and is also applied to optical computed tomography.7 8

In TOF measurements, the observed TRF is broadened by the instrumental response function (IRF), and thus it generally requires fitting to a model function that takes such broadening effect into account. In contrast to this conventional approach, the mean path length can be directly given by the first moment of the temporal response curve, because the contribution of the IRF is just a constant offset.9 This approach would be of great advantage to someone who is only interested in the mean path length, because with this method, the mean path length can be computed with a certain degree of accuracy. In addition, the analysis can be extended to the second moment of the TRF to separate the scattering and the absorption coefficients.10 However, in some measurements, the quality of the temporal profile is not good enough to calculate these moments because of a noisy background (shown in the latter part of this work), and it seems that the advantages of these approaches are limited. In this work, we focus on the peak of the temporal profile, which is easy to determine from data with less quality, and analyze its shift with different distances from the source position. The theoretical concept and error of the method are discussed using phantom measurements. We have validated an application of a TOF measurement on a human forearm at 1.29 μm, which was presented in our recent work.11

2.

Materials and Method

Two types of TOF measurements were conducted. The first set of measurements at 800 nm was conducted using a time-correlated single photon counting (TCSPC) system with a mode-locked Ti:sapphire laser (Tsunami, Spectra Physics, California). The light pulses were coupled to a 200-μm multimode fiber (FT200UMT, ThorLabs, New Jersey), and the light propagated in the sample was detected by a 62.5-μm GI multimode fiber (GIF625, ThorLabs). In these measurements, two kinds of photomultiplier tubes (PMTs), a multichannel plate (MCP) PMT (R2809-07, Hamamatsu Photonics, Hamamatsu, Japan), and a PMT module (H7141-40, Hamamatsu Photonics), whose IRFs were about 100 and 800 ps at full width at half maximum, respectively, were used. The response function of the PMT modules is significantly slow because it is not designed for time-resolved measurements. The experimental setup for the second set of measurements at 1.29 μm for a practical application is shown in Fig. 1. A mode-locked Ti:sapphire laser (Tsunami, Spectra Physics) pumped an optical parametric oscillator (Opal, Spectra Physics), which generated about 165-mW laser pulses at 1.29 μm with a 76-MHz repetition rate. The full width of half maximum of the pulse was about 110 fs with measurements by an autocorrelator (Spectra Physics). The wavelength of the NIR laser pulse was monitored by the wavelength of the second harmonic generation using a lithium triborate (LBO) crystal. The pulses were sampled by a wedge glass plate beamsplitter and detected for a reference for a time-to-amplitude converter (TAC; Ortec 567, EG&G, Oak Ridge, Tennessee) with a Ge:photodiode (G8198-02, Hamamatsu Photonics) with an amplifier (WP33, Kuranishi, Tokyo, Japan). The main laser beam was coupled to a 50-μm multimode fiber (FG050GLA, ThorLabs) and illuminated a sample. The light traveling in the sample was detected by a NIR photomultiplier tube (NIR-PMT; H9170-75, Hamamtsu Photonics) through a 100-μm multimode fiber (F-MLD-T, Newport, Irvin, California). The PMT signal was amplified (4774D/F, Hewlett Packard, Santa Rosa, California), discriminated (CFD; TC454, Tennelec, Oak Ridge, Tennessee), and fed to the signal input for the TAC. The TAC output was accumulated by a multichannel analyzer (MCA; 35plus, Camberra, Meriden, Connecticut).

Figure 1

Schematic diagram of the experiments in the 1.3-μm region.

526501j.1.jpg

Intralipid solutions of 1 and 2%, which were diluted from a 10% solution of Intralipid (Fresenius Kabi, Uppsala, Sweden), were used for the model samples in the measurements. The edges of source and detector fibers were embedded at about 20 to 26 mm below the surface of the Intralipid solution in a 300-ml beaker. An absorber (Acid Violet 6B, Chugai Kasei, Tokyo, Japan) was dissolved into the Intralipid solution in the measurements at 800 nm. The fiber separation was adjustable with about 1 mm as the accuracy of the separation.

3.

Theory

Temporal response of light traveling in a tissue or tissue-like sample is expressed by a solution of the optical diffusion equation when the distance between the source and detector is sufficiently larger than the mean free transport length 1/μs , where μs is the reduced scattering coefficient. In a homogeneous medium, its response at a distance ρ from the source point, which satisfies the condition ρ≫1/μs , can be approximately expressed by,

Eq. (1)

S(t,ρ)=At(3/2+ν) exp(μact)exp(ρ24Dct),
with ν=0 and ν=1 for an infinite and a semi-infinite configuration, respectively. D=1/(3μs ), μa, and c are the optical diffusion constant, the absorption coefficient, and the light velocity in the medium, respectively. We note that the reflectance from a semi-infinite sample should be calculated under a correct boundary condition. For instance, the extrapolated or partial current boundary condition is much better than the zero boundary one, and detailed comparisons of the solutions can be found in previous research.12 13 However, these solutions are reduced to a solution with the zero boundary condition under ρ2≫μs ′−2. Then, the solution is reduced to the earlier form, which was originally proposed by Patterson, Chance, and Wilson.2

The peak time of the TRF, t peak , satisfies dS(t,ρ)/dt=0, and is consequently derived from Eq. (1),

Eq. (2)

tpeak=12 (3μsμa)1/2 ρc [1+(3+2ν)212μsμaρ2]1/23+2ν4μac.
When the distance and the absorption are large enough, i.e., s μaρ2≫1, the expression can be simplified to a proportionality,

Eq. (3)

tpeak=12 (3μsμa)1/2 ρc+C,
where C is a constant. In principle, the fitting of the experimental data to Eq. (2) yields both the scattering and absorption coefficients. However, the offset of the time origin is dependent on the model. In addition, the IRF causes a time shift of the experimental data. Therefore, the simple proportionality, Eq. (3), is easier and better for use in the analysis.

When the TRF is exactly expressed by the response S(t,ρ), the intensity of the light should be equal to the infinite integral in the time domain of S(t,ρ). Then,

Eq. (4)

S(ρ)=A exp[(3μsμa)1/2ρ]/ρ  with ν=0,A(3μsμa)1/2[1+(3μsμaρ2)1/2] ×exp[(3μsμa)1/2ρ]/ρ2 with ν=1.
Using both the peak time and the intensity with different distances, the scattering and the absorption coefficients are separately obtained by fitting to Eqs. (3) and (4).

The mean transit time 〈t〉=∫ tS(t,ρ)dt/∫ S(t,ρ)dt can be analytically calculated from Eq. (1), and

Eq. (5)

t=1/2(3μs/μa)1/2ρ/c  with ν=0,
1/2(3μs/μa)1/2ρ/c[1+(3μsμaρ2)1/2]1  with ν=1.
Using s μaρ2≫1, both expressions are reduced to the same form,

Eq. (6)

t=12 (3μsμa)1/2 ρc.
Therefore, using Eqs. (3) and (6), the mean path length c〈t〉 can be evaluated from the slope of the peak time with different distances. The error of the estimation of mean path length from the slope is plotted in Fig. 2. The error at distance ρ can be estimated from the mean transit time given by ρdt peak /dρ, and the exact transit time given by Eq. (5). The estimated value from the slope is always small with the infinite configuration in Fig. 2(a). In contrast to this, it is slightly larger, about 2%, at a large (3μs μa)1/2 and at a longer distance of about 50 mm with the semi-infinite one in Fig. 2(b). With practical values, μs=1 mm −1 and μa=0.02∼0.05 mm −1, for tissue in the 750 nm region, (3μs μa)1/2 gives 0.24∼0.39 mm−1, and the error is less than 5% at a distance larger than 15 mm in both cases, infinite and semi-infinite. Therefore, the slope of the shift of the peak time with different distances will approximately give the mean transit time and thus the mean path length in the tissue sample.

Figure 2

Estimated error of the mean path length from the slope of the peak time shift with different source-detector distances with (a) infinite and (b) semi-infinite configurations. Each line is calculated with different values of (3μs μa)1/2 from 0.05 (bottom) to 0.35 mm−1 (top) with 0.05 mm−1 steps. The error is estimated from the ratio of the mean path length of the value estimated from the slope to the value calculated analytically.

526501j.2.jpg

In the previous discussion, the IRF is assumed to be ideally narrow and negligible. However, when the IRF is not negligible, which is in most practical cases, the peak time of the observed TRF is not trivially equal to that of the true TRF. The observed TRF is generally given by a convolution with an IRF R(t), and I(t,ρ)=∫ S(t−τ,ρ)R(τ)dτ. In the simplest case, where the true TRF just shifts in the time axis and the shape is invariant with different distances, the peak position of the convolution integral also shifts with the same value as the shift of the true peak position. Therefore, if the true response function is approximately invariant, the shift of the peak position of the observed response function must determine the mean transit time as described before. This intuitive description can be discussed in a general formulation of a convolution integral as indicated next.

When a TRF S(t) has a peak at t peak , it can be expanded as a Taylor series,

Eq. (7)

S(t)=n=0 1n!dnS(T)dTn|T=tpeak(ttpeak)n.
The observed TRF I(t) is given by the convolution integral with response functions S(t) and R(t), and the peak position of I(t) can be calculated from the derivative dI(t)/dt=0. The derivative of I(t) is given by

Eq. (8)

dI(t)dt= dS(T)dT|T=tτR(τ)dτ,
and, replacing S(t) by Eq. (7),

Eq. (9)

dI(t)dt=n=0 1n!dnS(T)dTn|T=tpeak(ttpeakτ)nR(τ)dτ.

Defining the solution of dI(t)/dt=0 as t peak +〈τ〉+ξ, where 〈τ〉≡∫τR(τ)dτ/∫ R(τ)dτ, Eq. (9) can be reduced to an equation of ξ. After a simple modification, the right term of the integral can be expressed by

m=0n(nm)(ttpeakτξ)mk=0nm()k(nmk)×τ˜kξnmk,
where 〈τ˜k〉=∫(τ−〈τ〉)kR(τ)dτ/∫ R(τ)dτ is the k’th moment of the distribution R(τ) around the average 〈τ〉, and eventually, ξ can be obtained as a solution of

Eq. (10)

n=1 1n!dn+1S(T)dTn+1|T=tpeakk=0n()k(nk)τkξnk=0.
When both S(t) and R(t) are even functions, dnS/dtn=〈τ˜n〉=0 with n=(odd). Eventually, ξ=0 is the unique solution of Eq. (10), and t peak obs =t peak +〈τ〉. On the other hand, when the shape of the response function is invariant around t peak with changing parameters, dnS/dtn is also invariant, resulting in a constant solution ξ, which is independent of these parameters, and t peak obs =t peak +〈τ〉+ξ. In general, the response curve is an approximately even function near the peak position, and thus its higher order derivatives dnS/dtn|t=t peak are small, although these are the indices of degree of asymmetry of the function. When the absorption becomes large, the shape change of the TRF is very small with different distances. Consequently, it is expected that the IRF contributes less to the change of the peak shift.

Using some practical values, the previous discussion can be confirmed. Since the IRF is asymmetric around a single peak and has a tail at later time points, the IRF is approximated by a log-normal distribution function R(t)=a exp(−{log[(t−θ)/m]}2/2/σ2)/(t−θ)/σ, where θ, σ, and m are the location, shape, and scale parameters, respectively. The shape parameter is fixed to the value 0.138 obtained by fitting to the IRF of the H9170-75 PMT. The scale parameter is selected to give a width of the response function, which is from 0.1 to 0.5 ns of the full width at e−2 of the maximum. These simulated IRFs and TRFs of the PMT are shown in Fig. 3(a). The simulated response with a 0.4-ns width is close to the PMT response.

Figure 3

The IRFs and TRFs calculated from the convolution of the instrumental function and the solution of the diffusion equation with a semi-infinite configuration. (a) shows the simulated IRFs by a log-normal function with different widths, 0.1, 0.2, 0.3, 0.4, and 0.5 ns, and that of the PMT measured by the laser scattering (thick solid line). (b) and (c) show the TRFs with μs =0.5 mm −1, μa=0.05 mm −1 and μs =0.5 mm −1, and μa=0.01 mm −1, respectively. The curves are calculated by the convolution with two different distances, 20 mm (upper curves) and 40 mm (lower curves), and different widths of the instrumental functions shown by (a). The thick solid curves in (b) and (c) denote the true response function.

526501j.3.jpg

The true TRF is calculated by the semi-infinite solution and numerically convolved with R(t). The simulated curves are shown in Figs. 3(b) and 3(c). The scattering coefficient μs =0.5 mm −1 is used because it is thought to be a typical value of tissue at the 1.29-μm region,11 though the value is arbitrarily determined with different distances ρ because μs ρ2 is the scaling factor in Eq. (1). The results with different absorption coefficients μa=0.05 mm −1 and 0.01 mm−1 are shown in Figs. 3(b) and 3(c), respectively. The different curves in the figure show the different distances, 20 mm (upper) and 40 mm (lower), with different widths of the IRF. The IRF gives a shift of the simulated response toward later time points because of the asymmetric shape of the instrumental response. A larger value of μa results in a much narrower and relatively symmetric shape of the response curve. Eventually, the shift of the simulated response from the ideal curve becomes extremely small.

The peak position of the simulated curve is estimated by a Gaussian function with an appropriate fitting range. The peak times are shown in Fig. 4. The bars in Fig. 4(a) denote the half width at 95% of the peak intensity of the fitted Gaussian function. The peak shift of data with smaller absorption is much larger than that with larger absorption. In addition, the nonlinearity of the shift with different ρ at lower absorption becomes significant as expected from the derivation of the proportionality, using s μaρ2≫1. The peak times of the broadened TRFs are larger than the peak times of the ideal TRFs, though the change of the shift is almost the same as the ideal TRFs as indicated by lines. This is clearly shown by the subtraction of the ideal peak times from the peak times of the broadened TRFs in Fig. 4(b). The difference is almost constant with different distances. The data at the lower absorption or μs ′1/2ρ<10 mm 1/2 with the slower instrument slightly deviate from the constant lines. However, the change of the shift of the peak time almost follows the ideal change, and the difference is almost within 5% deviation from the ideal one at μs ′1/2ρ>10 mm 1/2.

Figure 4

Peak times of the simulated TRF and the difference of broadened TRFs from the ideal ones. (a) shows the peak times of the broadened TRFs and ideal ones. Thick and thin lines are the results with μa=0.01 and 0.05 mm−1, respectively. The open and filled symbols denote 0.1 and 0.5 ns in the width of the instrumental response function, respectively. Solid lines denote the peak of the ideal response function. The dashed and broken lines fit the peaks of TRFs to the ideal peaks with a constant offset. The data with μa=0.01 mm −1 are plotted with an offset of 0.2 ns to avoid overlap of the curves in the figure. (b) shows the difference from the ideal curve. Each symbol shows the results of the simulation with the same parameters as in (a). The lines show the constant offsets.

526501j.4.jpg

4.

Experimental Results

First, the measurements at 800 nm were conducted because the measurements of a slow instrument could be compared with those of another, faster one. In the slow instrument measurements, the MCP-PMT was replaced by the PMT module equipped with a level converter from a transistor-transistor logic (TTL) signal to a nuclear instrumentation methods (NIM) signal in our TCSPC setup. The temporal response curves of the light traveling in a 1% Intralipid solution with 2.7 g/l acid violet 6B in a 300-ml beaker with the slow instrument are shown in Fig. 5. The source and detector fiber edges were positioned at the middle of the beaker at a depth 26 mm from the surface. The profile of the TRF with the different distances was almost the same as that of the IRF shown by the thin solid line. In contrast, the peak time of the TRF showed a clear shift toward later time points with an increase in the distance.

Figure 5

Temporal response curves from a 1% Intralipid/Acid Violet solution with different distances from the source position at 800 nm. (a) shows temporal response curves normalized by the acquisition time with ρ=10, 20, 30, and 40 mm, from bottom to top. (b) shows temporal response curves normalized by their peak counts with different distances ρ=10, 20, 30, and 40 mm, from left to right. Thin solid lines in both figures show the IRF.

526501j.5.jpg

The peak positions, which were estimated by a curve fitting of Gaussian function around the peak position, are plotted in Fig. 6(a). The results with the fast instrument are also plotted in this figure. The slopes, calculated at a distance of 15 to 40 mm, are 20.8±0.6 and 19.8±0.3 ps/mm with the slow and fast instruments, respectively. These slopes are in good agreement within the margin of error. Using the light velocity 2.29×1011mm/s, the estimated mean path lengths are (4.7±0.1)ρ and (4.53±0.06)ρ, respectively. The absorbance of 2.7 g/l acid violet 6B in water was 0.137 cm−1 at 800 nm and thus μa=0.0315 mm −1. Using the scattering coefficient of 1% Intralipid at 800 nm, μs =1.0 mm −1, 14 and Eq. (6), the mean path length is expected to be 4.87ρ. This value is in very good agreement with the one estimated from the slope.

Figure 6

Peak time and the first moment of the temporal response curves at 800 nm. (a) peak time with distances from 10 to 40 mm, (b) the intensities, and (c) the first moment referenced by that of the IRF. The right vertical axis shows the mean path length estimated with assumption of the light velocity 2.29×1011mm/s. The circles and boxes denote data with a PMT module and a MCP-PMT, respectively. The solid and broken lines denote best-fit curves.

526501j.6.jpg

The intensities of the light traveling in the sample are shown in Fig. 6(b). The curves are the best fit curves using Eq. (4). The fitting to data at a distance of 15 to 40 mm yields (3μs μa)1/2=0.31±0.01 and 0.288±0.001 mm−1 with the slow and fast instruments, respectively. Although the small difference between the instruments might be due to the uncertainty of neutral density filter attenuation that was critically dependent on the position of the incident light or the relatively low sensitivity of the slow instrument, these results are consistent within about a 5% difference. Combined with the peak shift data, μs =0.98 mm −1 and μa=0.032 mm −1 with the slow instrument were yielded. Similarly, μs =0.87 mm −1 and μa=0.032 mm −1 with the fast instrument were yielded, and μs was slightly inaccurate compared with the expected value of 1 mm−1 in our measurements.

For a comparison with other methods, the first moment of the temporal profile 〈t〉 relative to that of the IRF 〈t0 is shown in Fig. 6(c). The symbols show the results with an estimation using data points from −0.5 to 2.5 ns. The slopes of the first moment with different source-detector distances with the slow and fast instruments are 18.2±0.2 and 19.7±0.2 ps/mm, respectively. The first moment with the slow instrument is slightly smaller than that with the fast one. The first moment at a low count region with the slow instrument slightly deviates from the linear line because of the bad signal-to-noise ratio. It can also be noted that the first moment is significantly dependent on the region of the calculation, and this is especially critical in a low intensity region because the background count level is not flat in all channels and effects the estimation. A simple subtraction, 〈t〉−〈t0〉, directly removes the broadening due to the IRF and gives a mean transit time. Using data at ρ=20 mm, the estimated mean path lengths are (4.5±0.3)ρ and (4.7±0.1)ρ with the slow and fast instruments, respectively. Thus, the simple subtraction gives the mean path length with similar accuracy. However, this value is very critical to the evaluation of the first moment of the instrument, and the relative measurements with different distances reduce the uncertainty of the results. Thus, the peak shift analysis is as compatible a method in terms of accuracy as the first moment analysis.

The scattering and absorption coefficients can be evaluated by a conventional fitting using Eq. (1). Using the data with the fast instrument, the values as μs =0.92±0.01 mm −1 and μa=0.0360±0.0004 mm −1 at ρ=20 mm were yielded. Then, the mean transit time was calculated as 4.4, and this is in good agreement with the prior two methods. The fitting, using an improved solution, will yield much more accurate values, though the values with the previous simple solution are accurate within a 10% margin of error. On the other hand, using the data with the slow instrument, the fitting yielded μs =1.3±0.1 mm −1 and μa=0.048±0.001 mm −1. These values systematically deviate from the expected values and are not yielded with different distances reproducibly. This is because the IRF is too slow to extract the true TRF from the measured TRF. The discussion on the accuracy of the fitting to the diffusion model can be found in the literature.15 Nevertheless, the peak shift analysis accurately estimates the mean path length with the slow instrument.

At a region above 1.1 μm, the speed and sensitivity of the detector are very limited. The dark noise is also very large. The peak shift analysis is currently the only applicable method under this condition. Figure 7 shows TRFs of the scattering light at different positions from the source position. The raw data had a significantly large constant offset due to the very large dark count of the PMT, as shown by Fig. 7(a). In addition, the longer accumulation time enhanced an oscillatory noise, which probably originated from the radio frequency noise of the electronics. In particular, the raw data at longer distances shown in Fig. 7(a) shows that there existed large noise over the whole time range, and the true response could not be identified without the primary knowledge of the timing of the incident light. In addition, it was very difficult to accurately determine the dark count level because of the large oscillatory background noise. Therefore, we simply assumed a constant dark count, which was estimated by the average count just before the response curve began to rise, and was subtracted from the data. These results are shown in Fig. 7(b). The accumulation time of data was also calibrated, so that the vertical axis corresponded to the photon count per second in each bin (6.4 ps/bin). The signal intensity was significantly attenuated with an increase in the distance. The line shape of the data is very similar to that of the instrumental one measured by the laser scattering from a piece of white paper (thin solid line). The small peaks at 1 ns might be due to scattering from the optics because it was dependent on the size of the aperture of the PMT housing. The distance between the input fiber receptor and the PMT is about 10 to 15 cm, and this value suggests that this is a possible origin of these peaks. However, the origin could not be identified because the optics inside the PMT housing could neither be removed nor modified.

Figure 7

Temporal response curves from a 2% Intralipid solution with different distances from the source position at 1.29 μm. (a) shows the raw data with 15, 20, 25, and 30 mm, from bottom to top. t acc indicates the accumulation time of the measurement. (b) shows temporal response curves normalized by the accumulation time with different distances, 15, 20, 25, and 30 mm, from top to bottom. The thin solid line denotes the IRF.

526501j.7.jpg

Normalized temporal response curves are shown in Fig. 8. Figure 8(a) shows the normalized curves by the peak count. The peak of the response curve showed a clear shift toward later time points with an increase in the distances. The curve with the longest distance was affected by the large noise around the 1-ns region. The peak time of the curve was estimated by fitting to a Gaussian function, and the curves normalized by the peak time of the curve with a distance of 15 mm are shown in Fig. 8(b). The curves overlap each other and are indistinguishable, except for that with a 30-mm distance. The difference between these curves and the IRF is also very small, as shown in Fig. 8(c), which shows the ratio between the measured TRFs and IRF with the normalization of the peak time. Therefore, the IRF was too slow to resolve the difference of the shape of sample response within the margin of experimental errors. A very small profile difference between the IRF and the measured TRF causes the large error of the conventional fitting method to the theoretical response function, such as Eq. (1), and a large uncertainty of the parameters. For instance, a fitting to the results at ρ=20 mm yields μs =0.077±0.037 and μa=0.060±0.017 mm −1, the values of which are significantly different from the previously obtained values16 and are unreliable because of the large uncertainties.

Figure 8

TRF normalized by the peak intensity. (a) shows the normalized TRFs with distances 15, 20, 25, and 30 mm, from left to right. (b) shows the normalized TRFs, aligned at the peak position. (c) shows the difference between TRF and IRF.

526501j.8.jpg

The intensity of the detected light and the peak time are plotted in Figs. 9(a) and 9(b), respectively. The intensity was estimated by the total count of the TRF, subtracting a constant dark count estimated by the total count of data with no input signal (8.46±0.09 kcps). The estimated intensity at 30 mm was negative because this might be caused by the small intensity of TRF and overestimation of dark count level of data due to the large amount of fluctuating noise. The peak time of the response curve was estimated by Gaussian fitting. The intensity was exponentially attenuated with an increase in the distance. The fitting by Eq. (4) is shown by the solid line in the figure. This line is in good agreement with the data points, and eventually (3μs μa)1/2=0.47±0.01 mm −1 is yielded using c=2.29×1011mm/s. Thus, the experiment satisfied the assumption s μaρ2≫1. On the other hand, the peak time was linearly shifted with an increase in the distance, as shown in Fig. 9(b). Again, the peak time at 30 mm deviates from the other data points because of the distortion of the temporal response curve. The linear line agrees well with the data points except for the point at 30 mm. The slope of the line yielded is 9.7±0.6 ps/mm. This indicates that the mean path length is extended by a factor of 2.2 of the geometric distance between the source and the detector. Eventually, we are able to obtain μs =0.69 mm −1 and μa=0.11 mm −1 with the 2% Intralipid solution at 1.29 μm. These values are almost in agreement with those obtained by the integrated sphere experiment.16 It can be concluded that the absorption of the sample is most likely determined by the water.17

Figure 9

(a) Intensity and (b) peak time of the TRF for the 2% Intralipid solution with different distances from the source position. The line is a fitting to the model.

526501j.9.jpg

The trial of the first moment analysis of the previous data is shown in Fig. 10. Two calculation regions of the first moment yield completely different slopes, −0.6±10 and 75±22 ps/mm for −0.5 to 1.0 and −0.5 to 3.0 ns, respectively, and thus the results are strongly dependent on the region. It is very difficult to determine the appropriate region of the first moment calculation. This might be because of the low sensitivity and large noise of the PMT.

Figure 10

The first moment of TRF with different distances. The calculation range from −0.5 ns to 3 and 1 ns are plotted by circles and squares, respectively. The lines show the best-fit lines.

526501j.10.jpg

5.

Discussion

The peak shift analysis presented here is very simple and has a great advantage when the IRF of PMT is not good enough to use conventional fitting methods. In the visible region, a PMT with good temporal response can be selected in many cases, as the first set of measurements in this work demonstrated. However, it is currently impossible in the 1.3-μm region. The PMT used in this experiment is the fastest and most sensitive one to the authors’ knowledge. Even using this PMT, the response time is about 350 ps, as shown in Fig. 7, and too slow to resolve the temporal broadening by the scattering in tissue.11 Moreover, the highly sensitive PMT has the disadvantage of having a high dark level, which causes the noisy background in the TCSPC data. This is also why it is difficult to apply conventional techniques such as deconvolution and the first moment analysis. The first moment analysis in the visible region with the slow PMT module yields accurate results, because the noise level of the PMT is very small relative to the NIR-PMT. Nevertheless, in such a NIR-PMT case, the peak time and the intensity are the only analyzable parameters from the temporal response curve in this region.

The peak shift analysis uses only the information of the change of the peak position in time with different spatial distances (greater than centimeters), and yields the mean path length of the light traveling in the sample. It does not require the deconvolution procedure. In addition, the temporal origin of the observed temporal response can be arbitrary to that of the instrumental one because the measurement is a relative method. Other analyses using moments of the TRF have been performed.9 10 18 Since the first and second moments are obtained as simple sums of those of the true TRF and the IRF, deconvolution is not required in the analysis. For instance, the mean path length can be estimated by the time average of the observed TRF and the IRF without deconvolution. From the measurements at 800 nm in Fig. 5, the peak shift and the first moment analyses provide the same order of accuracy in the estimation of mean path length.

This first moment method, however, is dependent on the accuracy of the time average of the TRF. For instance, it is very difficult for the evaluation to obtain an accurate enough average, because it is critically dependent on the background subtraction in our measurements at 1.29 μm. We assumed a constant background level in the measured TRF. As shown in Fig. 7, the background noise has an oscillatory behavior. Thus, the assumption of the constant background might not be appropriate, and there still exists a difficulty in the determination of the background level. This oscillatory behavior also causes errors in the determination of the first moment. Figure 10 shows that the estimation of the first moment is quite dependent on the region of the calculation and completely unreliable due to background noise. In a similar case, the IRF was quite noisy due to the radio frequency noise for the time-resolved photometry using a digital oscilloscope.19 This obscured the details of the response curve in TOF measurements of the tissue. In both cases, however, the peak is still clearly observed and the peak shift is analyzable. The main point is that the peak position has the highest count of data, which is the reason why it is robust in the analysis with broadened data under a low signal-to-noise ratio. Therefore, the peak shift analysis is another option that is widely applicable.

The disadvantage of this approach is the scanning of the detection position. Consequently, the results are only accurate for light traveling in a homogeneous medium. Since the tissue sample is heterogeneous, this analysis provides only a coarse graining of it to characterize the sample, and conversely the deviation from the linear relationship, which is the key to our analysis, indicates the heterogeneity of the sample. Although we emphasize that the peak shift analysis is less sensitive to the response time of the instrument, the slower response will cause an error in the peak estimation. The broader shape of the temporal response does not resolve the peak shift. For instance, in the case of a 5% uncertainty of the intensity at the peak, the peak count is 400 counts/bin, assuming Poisson distribution of count data. In this case, the peak position can be determined with a change of TRF more than 5%, otherwise the peak position is indistinguishable from the noise. The width of the 5% change of the Gaussian distribution function is given by 0.32σ, where σ is a half width at e−1/2 intensity. This width, using TRFs, is plotted as error bars in Fig. 4. In our measurements at the NIR region, the peak shift was about 10 ps/mm. Optimistically assuming a 100-ps accuracy of the position, the TRF should decay with σ=312 ps. This value is about twofold slower than the measured TRFs, and thus the measurements with a twofold slower IRF are analyzable. In many measurements, the peak time is determined by fitting to data with a region wider than a 5% change, and thus the width of the TRF can be wider than that in the earlier discussion. On the other hand, in practice, the required accuracy of the peak time would be 20 to 50 ps to distinguish the change corresponding to a 10-mm step. Further investigation is necessary into the limit of the width of the IRF.

6.

Conclusion

The peak time of the TRF of light traveling in a tissue-like phantom linearly shifts with different distances between the source and detector positions. This shift is less dependent on the IRF, and therefore, the peak shift can be analyzed by the observed temporal response function without a deconvolution procedure. The slope of the peak time gives the mean path length, and combining the intensity information yields μs and μa of the sample. This analysis can be applied to measurements at the 1.3-μm region, where instrumental response is limited by the available PMT. This simple analysis is widely applicable to characterize a homogeneous sample.

Acknowledgment

The authors acknowledge Hamamatsu Photonics K.K. for the support of the NIR-PMT module.

REFERENCES

1. 

H. M. Heise, “Applications of near infrared spectroscopy in medical sciences,” in Near-Infrared Spectroscopy—Theory, Instruments, and Applications, H. W. Siesler, Y. Ozaki, S. Kawata, and H. M. Heise, Eds., p. 289, Wiley and Sons, New York (2002).

2. 

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 2331 –2336 (1989). Google Scholar

3. 

S. L. Jacques , “Time-resolved reflectance spectroscopy in turbid tissues,” IEEE Trans. Biomed. Eng. , 36 1155 –1161 (1989). Google Scholar

4. 

M. Oda , Y. Yamashita , G. Nishimura , and M. Tamura , “A simple and novel algorithm for time-resolved multiwavelength oximetry,” Phys. Med. Biol. , 41 551 –562 (1996). Google Scholar

5. 

R. Cubeddu , A. Pifferi , P. Taroni , A. Torricelle , and G. Valentini , “Noninvasive absorption and scattering spectroscopy of bulk diffusive media: An application to the optical characterization of human breast,” Appl. Phys. Lett. , 74 874 –876 (1999). Google Scholar

6. 

A. Torricelli , A. Pifferi , P. Taroni , E. Giambattistelli , and R. Cubeddu , “In vivo optical characterization of human tissues from 610 to 1010 nm by time-resolved reflectance spectroscopy,” Phys. Med. Biol. , 46 2227 –2237 (2001). Google Scholar

7. 

E. M. C. Hillman , J. C. Hebden , M. Schweiger , H. Dehghani , F. E. W. Schmidt , D. T. Delpy , and S. R. Arridge , “Time resolved optical tomography of the human forearm,” Phys. Med. Biol. , 46 1117 –1130 (2001). Google Scholar

8. 

H. Zhao , F. Gao , Y. Tanikawa , Y. Onodera , M. Ohmi , M. Haruna , and Y. Yamada , “Imaging of in vitro chicken leg using time-resolved near-infrared optical tomography,” Phys. Med. Biol. , 47 1979 –1993 (2002). Google Scholar

9. 

H. Zhang , M. Miwa , T. Urakami , Y. Yamashita , and Y. Tsuchiya , “Simple subtraction method for determining the mean path length traveled by photons in turbid media,” Jpn. J. Appl. Phys. , 37 700 –704 (1998). Google Scholar

10. 

A. Libert , H. Wabnitz , D. Grosenick , M. Mo¨ller , R. Macdonald , and H. Rinneberg , “Evaluation of optical properties of highly scattering media by moments of distributions of times of flight of photons,” Appl. Opt. , 42 5785 –5792 (2003). Google Scholar

11. 

G. Nishimura and M. Tamura, “Peak time analysis of TOF data with limitation of the temporal resolution and its application for measurements on a human forearm at 1.29 μm,” preseted at Biomed. Topl. Mtg. on CD-ROM, 2004, Washington, D.C., Paper WF40, OSA.

12. 

A. H. Hielscher , S. L. Jacques , L. Wang , and F. K. Tittel , “The influence of boundary conditions on the accuracy of diffusion theory in time-resolved reflectance spectroscopy of biological tissues,” Phys. Med. Biol. , 40 1957 –1975 (1995). Google Scholar

13. 

A. Kienle and M. S. Patterson , “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A , 14 246 –254 (1997). Google Scholar

14. 

H. J. van Staveren , C. J. M. Moes , J. van Marle , S. A. Prahl , and M. J. C. van Gemert , “Light scattering in Intralipid-10% in the wavelength range of 400–1100 nm,” Appl. Opt. , 30 4507 –4514 (1991). Google Scholar

15. 

V. Ntziachristos and B. Chance , “Accuracy limits in the determination of absolute optical properties using time-resolved NIR spectroscopy,” Med. Phys. , 28 1115 –1124 (2004). Google Scholar

16. 

T. L. Troy and S. N. Thennadil , “Optical froperties of human skin in the NIR wavelength range of 1000–2200 nm,” J. Biomed. Opt. , 6 (2), 167 –176 (2001). Google Scholar

17. 

G. M. Hale and M. R. Querry , “Optical constants of water in the 200-nm to 200-μm wavelength region,” Appl. Opt. , 12 555 –563 (1973). Google Scholar

18. 

E. M. C. Hillman , J. C. Hebden , F. E. W. Schmidt , S. R. Arridge , M. Schweiger , H. Dehghani , and D. T. Delpy , “Calibration techniques and datatype extraction for time-resolved optical tomography,” Rev. Sci. Instrum. , 71 3415 –3427 (2000). Google Scholar

19. 

G. Nishimura and M. Tamura , “Simple setup for nano-second time-resolved spectroscopic measurements by a digital storage oscilloscope,” Phys. Med. Biol. , 48 N283 –N290 (2003). Google Scholar
©(2005) Society of Photo-Optical Instrumentation Engineers (SPIE)
Goro Nishimura and Mamoru Tamura "Simple peak shift analysis of time-of-flight data with a slow instrumental response function," Journal of Biomedical Optics 10(1), 014016 (1 January 2005). https://doi.org/10.1117/1.1854684
Published: 1 January 2005
Lens.org Logo
CITATIONS
Cited by 8 scholarly publications.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Tissue optics

Absorption

Scattering

Deconvolution

Error analysis

Geometrical optics

Light scattering

Back to Top