Translator Disclaimer
1 November 2007 Time-resolved diffuse optical tomography and its application to in vitro and in vivo imaging
Author Affiliations +
This work reviews our research during the past several years on time-resolved (TR) near-infrared diffuse optical tomography (DOT). Following an introduction of the measuring modes, two proposed schemes of image reconstruction in TR-DOT are described: one utilizes the full TR data, and the other, referred to as the modified generalized pulse spectrum technique (GPST), uses the featured data extracted from the TR measurement. The performances of the two algorithms in quantitativeness and spatial resolution are comparatively investigated with 2-D simulated data. TR-DOT images are then presented for phantom experiments, which are obtained by using a 16-channel time-correlated single photon counting system, and the factors affecting the quantification of the reconstruction are discussed. Finally, in vitro and in vivo imaging examples are illustrated for validating the capibility of TR-DOT to provide not only the anatomical but also the physiological information of the objects.



Diffuse optical tomography (DOT) is an emerging biomedical imaging technology. Compared to the established imaging modalities, e.g., x-radiography, ultrasound, and magnetic resonance imaging, DOT offers distinct advantages in terms of sensitivity to functional changes, safety, and bedside use. The DOT mechanism in providing both pathological and physiological information of biological tissue lies in the spectroscopic reconstruction of the absorption coefficient (μa) and reduced scattering coefficient (μs) , e.g., a tumor can be imaged as a volume of increased absorption coefficient,1, 2 and the metabolic status of tissue can be accessed by monitoring the blood oxygenation, which is obtained from reconstructing absorption coefficients at multiple (at least two) wavelengths.3, 4, 5

The absorption and scattering coefficients express the probabilities of photons being absorbed and scattered per unit length of travel, respectively. The absorption of biological tissue in the near-infrared (NIR) wavelength range comes mainly from water, hemoglobin in blood, myoglobin in muscle, and melanin in skin. Although some of these chromophores, such as water and melanin, exist with a great extent, their contents in tissue can be treated constant during some period of time of physical and physiological activities. Therefore, the constituent contributing to a change in absorption is primarily hemoglobin, which exists within the red blood cell in oxygenated (oxy-) and deoxygenated (deoxy-) forms. Since they have distinct characteristic absorption spectra and weak absorption in the NIR wavelength range, it enables the detection of light propagated through several centimeters of tissue and monitoring of the change in the pathologically or physiologically related hemoglobin concentration.6 However, the scattering of NIR light by biological tissue is usually tens to hundreds of times stronger than the absorption. This fact brings difficulties in collecting and interpreting the measured light, since the signal measured at several centimeters apart from the source is dominated by diffused light.7, 8, 9

In DOT, light from an array of sources is orderly injected into tissue, a part of the light re-emitted from the tissue after multiple scattering is collected by an array of detectors in parallel, and then a model quantifying light propagation in the tissue is employed for reconstructing the distributions of the optical properties in the tissue through best fitting the model prediction and the measured dataset.

To detect the weak light re-emitted from biological tissue, three kinds of measurement systems have been developed so far, i.e., a continues-wave (cw) system, a frequency-domain (FD) system, and a time domain or time-resolved (TR) system. In a cw system, the light intensity is constant or modulated at a frequency not higher than a few tens of kilohertz, and the attenuation of the incident light after passing through the tissue is the measurable parameter.10, 11, 12 As an example of its application, cyclic hemodynamic changes revealed by a cw tomographical system have been demonstrated;13 commercial cw systems for breast imaging have been established by Philips Research Laboratories.14 In a FD system, light sources shine on tissues continuously, but their intensities are modulated at frequencies in the order of tens to hundreds of megahertz, and the measurable parameters are the phase delay and the amplitude attenuation of both the DC and AC components. Many researchers and groups have made great progress in the frequency domain system and its applications.15, 16, 17

In a TR system, pulsed lasers with ultrashort duration irradiate biological tissues, and the temporal distributions of re-emitted photons known as the temporal profiles are measured.5, 18 The up-to-data status of TR-DOT has been reported by the researchers at University College London, who obtained 3-D tomography images of neonatal heads using their purpose-built TR imaging system.19 Comparing the three kinds of measurement methodologies, cw systems are simple in instrumentation and fast in data acquisition, but a multispectral approach with correct selection of the measurement wavelengths is required for the separation of the absorption information from the scattering one.11 Although the total intensity and the mean photon flight time calculated from a temporal profile measured by a TR system are equivalent to the amplitude and phase measured by a FD system, a TR system can provide additional information from making full use of the temporal profiles.1, 20 Despite its high cost of instrumentation and long times for data collection, a TR system using the time-correlated single photon counting technique has proved to be more sensitive than a FD system when applied to a bulky tissue of several centimeters.1

For imaging of optically thick tissues, as in the case of NIR-DOT, the most appropriate model of light propagation (forward mode) is the diffusion equation (DE), which is a P1 approximation of the more physically rigorous radiative transfer equation. Reconstruction of the optical properties with the known source and detector locations and the measured dataset can be achieved by solving the inverse problem of DE, which is intrinsically nonlinear and ill-posed. For a geometry with an arbitrary shape and arbitrary distributions of optical properties, the DE can only be solved numerically, e.g., using a finite element method (FEM).21 The linearization of the inverse problem based on a FEM spatial discretization generally leads to an iteration scheme with a core task of solving a severly underdetermined and ill-posed matrix equation, where the number of measurable quantities (measurement set) is much less than that of unknowns, i.e., optical properties in the FEM nodes, and the inversion is noise sensitive. Efforts have been made to find a reasonably regularized solution to the ill-posed inverse problem. 22, 23, 24, 25, 26, 27, 28, 29, 30, 31 But, challenges still remain in the image reconstruction of DOT, mainly on the improvements of the spatial resolution, the quantitative characterization of heterogeneities, and separation between the absorption and reduced scattering coefficients.

This work presents our achievements over the past several years in the development of image reconstruction algorithms, as well as in phantom, in vitro, and in vivo experimental results.


Image Reconstruction Algorithms


General Framework

As described before, given a forward model F with a vector of the optical properties p , including both the absorption coefficient μa and diffusion coefficient κ=1(3μs) (or scattering coefficient μs ), the forward problem can be expressed as y=F(p) , where y denotes a vector of the measured quantities, and the inverse problem is accordingly written as p=F1(y) . If the actual optical properties p are close to its initial estimation p0 , then the nonlinear inverse problem posed by DOT can be linearized as

Eq. 1

where F is the first-order Frechet derivative of F , Δy=yF(p0) , and Δp=pp0 . When the forward model is solved by a numerical technique, the discrete form of F , referred to as the Jacobian matrix J , can be efficiently obtained using a perturbation method of the DE.

The widely used reconstruction algorithm is mostly based on a nonlinear scheme, where an iterative procedure is utilized to update both the optical properties and Frechet derivative until the best fitting between the measured and the computed quantities is achieved. At the i ’th updating stage, the iteration procedure can be expressed as

Eq. 2

where χ is a column vector of the computed measurements in length of S×D ( S and D are the numbers of source and detector positions used in the measurement, respectively), F(pi) represents the solution of DE by FEM at the i ’th updating stage with the distribution of the optical properties pi ,
is the Jacobian matrix of F with respect to both μa and κ , and δpi=[δμa,δκ]T is the perturbation of the optical properties at the i ’th updating stage. Thus, the main task of image reconstruction is to update the Jacobian matrix for a known value of pi and determine δpi at the iteration stage. For solving the underdetermined and ill-posed matrix equation, the calculation of δpi is generally performed using regularization techniques, e.g., the algebraic reconstruction technique (ART).26


Featured-Data Algorithm

In TR-DOT algorithms, some featured data types are employed for characterizing the temporal profile and saving the computation time in the image reconstruction. The data types prevalently employed in TR-DOT are derived from the Mellin transform of the temporal profile, such as the intensity E , mean time of flight (TOF) t , and high-order central moments such as the variance c2 , etc., which represent the cw component, the mean time position, and the time width of the temporal profile, respectively.27, 28, 29 However, many practices with the Mellin-transform-based data types have shown that the high-order central moments are rather noise sensitive. As a results, only E and t have been put into use.

To overcome these limitations and improve the image quality, we have proposed a new featured-data algorithm that uses the Laplace transform of the measured temporal profile as the data type. In addition, we have also developed an algorithm using the full temporal profile to demonstrate the advantages of the TR scheme over the other two, as well as offer a potential “gold standard” for the development of the featured-data algorithms.


Generalized pulse spectrum technique

The generalized pulse spectrum technique (GPST) is based on converting the time-dependent signal by Laplace transform to the complex frequency domain, and then using the ratio of the Laplace-transformed TR measurements at two distinct real domain frequencies.30, 32 For the time-dependent re-emission χ(ξd,ζs,t) in a tissue volume Ω with a boundary Ω , the Laplace transform of χ(ξd,ζs,t) is given by

Eq. 3

where ξd and ζs denote the detector and source positions, respectively, and f is referred to as the frequency. It can be readily seen that the Jacobian matrix Jp(ξd,ζs,f) (p[μa,κ]) for χ̂(ξd,ζs,f) is in similar form to data-type E .32

The ratio of the Laplace-transformed signals at two distinct frequencies (denoted by f1 and f2 ) is employed as a data type:

Eq. 4


Thus, the Jacobian matrix for the data-type R can be reformularized according to those for χ̂(ξd,ζs,f1) and χ̂(ξd,ζs,f2) :25

Eq. 5



Modified generalized pulse spectrum technique

For the simultaneous reconstruction of μa and κ , the Jacobian matrix in the previous standard GPST needs the computation of the gradients of both the Green’s function of the forward model and the re-emitted flux χ , which is much more difficult and time consuming than for the case of μa -only reconstruction. Therefore, the modified GPST algorithm has been further developed, aiming at substantially reducing the computation time. The key step of the modified GPST is to derive an approximate μs -regarding Jacobian matrix Jμs , that is free of gradient computation of the flux Green’s function and photon density.25, 32

The Jacobian matrix in the modified GPST for R is calculated as below:32

Eq. 6

where Jμa(ξd,ζs,f) is the μa -regarding Jacobian matrix for χ̂(ξd,ζs,f) , and

Eq. 7


It can be seen that, with the modified GPST, the simultaneous reconstruction of μa and μs is computationally as efficient as the μa -only reconstruction from the point of view of the Jacobian matrix computation.


Scaling of the Jacobian matrix

In general, the entries of Jμa(R) are several orders of magnitudes larger than those of Jμs(R) , resulting in the μa -reconstruction more pronounced and the simultaneous reconstruction of both the μa and μs being unbalanced. To enhance the efficiency of the modified GPST algorithm in μs -reconstruction, a scaling strategy is used to equalize the matrices by normalizing the matrix along a row, and the i ’th row of the Jacobian matrix can be written:

Eq. 8

where the scaling factors wμa(i) and wμs(i) are the maxima of the elements at the i ’th row of the μa - and μs -regarding Jacobian matrices for R , Ĵμa , and Ĵμs , respectively. It is noted that the reconstructed optical properties need to be de-equalized after the reconstruction.


Full Time-Resolved-Data Algorithm

According to the general definition of DOT in Eq. 1, the Jacobian matrix in TR-DOT using full TR data can similarly be derived using a perturbation procedure. By adding small perturbations δμa and δκ to μa and κ , respectively, one can calculate the perturbation of the photon density δΦ from the DE and obtain δχ by taking into account the relation between Φ and χ . Then the Jacobian matrix at the time tm can be formalized as:

Eq. 9

where ui(r) is the shape function at the i ’th node of the FEM mesh and Ω is the domain of interest. The following two considerations have to be taken for using the full TR data in DOT image reconstruction.
  • Since the intensity information is used, a mathematical method has to be developed to eliminate the requirement for the absolute intensity scaling, which is difficult in actual measurement. One of the possible methods is to normalize the original TR signals by their respective cw components, and generate a self-normalized temporal profile χ¯(ξd,ζs,t)=χ(ξd,ζs,t)E(ξd,ζs) .

  • The original TR data contains significant redundancy that is dependent on the temporal resolution of the measuring system, and might be easily influenced by measurement noise. To take account of the temporal resolution of the system and improve the robustness of the algorithm to measurement noise, one effective way is to divide a temporal profile into many time slices with a time width of Δt , and integrate the data in this time slice. Then one can obtain a data for the m ’th time slice,

    as illustrated in Fig. 1 .31

Fig. 1

The time slice technique employed in the reconstruction algorithm using the full time-resolved data. A temporal profile was divided into about ten time slices with a width of Δt .


With the prior considerations, the Jacobian matrix can be obtained accordingly,

Eq. 10

where Jp(E) are the Jacobian matrices for E . Jp(Δt) is in similar form to Eq. 9, but replacing χ(ξd,r,tmτ) with χ(Δt)(ξd,r,tmτ) . Similar to the GPST algorithm, a scaling strategy is also needed here for the equalization of the Jacobian matrix.


Comparisions of the Algorithms Using the Featured Data with that Using the Full Time-Resolved Data

The performances of the algorithms based on the featured data types E , t , E and t , t and c2 , R , and on the full TR data, are evaluated on the aspects of the quantitativeness in μa reconstruction, the execution time, and the cross talk between μa and μs ( μs=3κ is used instead of κ , hereafter) in the simultaneous two-parameter reconstruction. All the reconstructions in this section were obtained from a numerically simulated 2D domain with a diameter of 80mm and the background optical properties of μa=0.01mm1 and μs=1.0mm1 . 16 sources and detectors were assumed to be located with an equal spacing around the boundary. The TR re-emissions or their featured data types were calculated from the 2-D diffusion model using a FEM with an appropriate mesh density (2400 triangle elements connected at 1261 nodes) and time interval (20ps) to ensure the numerical accuracy.31


Quantitativeness and Execution Time in μa Reconstruction

There are many methods for evaluating medical images; those in common use are contrast-detail analysis, spatial resolution testing, and receiver operating characteristic analysis. The contrast-detail analysis, which indicates whether targets with different contrasts and with different sizes are observable, is more suitable than spatial resolution testing when the targets have finite contrasts to the background, as in the case of DOT.33 For evaluating images in revealing targets with different contrasts but with identical size, we define the ratio of the reconstructed peak value of the optical property to the true one, referred to as quantitativeness, as a measure. By the definitions, the quantitativeness provides information about not only the visibility of a target but also how correctly the target is reconstructed. In a NIR-DOT image, tissues are usually distinguished according to the contrast of the optical properties, and the pathological changes of tissue can therefore be deduced in terms of the changes in the optical properties, so the quantitativeness will be a good measure for evaluating the reconstructed images.

The aforementioned 2-D domain containing two inclusions with a diameter of 16mm was employed. The two inclusions are positioned symmetrically to the center with a separation of 40mm . Figure 2 shows the reconstructed μa when the two inclusions have the absorption contrast of 2, 3, 5, or 7 to the background while they have the same scattering as the background. We can see that when the contrast is 2, the algorithm using the prior featured data types and the full TR data have almost identical performances in reconstructing μa , and the quantitativeness ratio reaches as high as 85%. As the contrast increases, the peak value of reconstructed μa keeps almost constant for the featured-data algorithms, and thus the quantitativeness ratio decreases. As a whole, the algorithm using the full TR data provides better quality in quantitative reconstruction, e.g., when the contrast is 7, the quantitativeness ratio is only about 0.3 from the GPST algorithm but 0.5 from the full TR data algorithm.

Fig. 2

Influence of the absorption contrast on the reconstructed absorption coefficient by the algorithms using the featured data types and the full TR data. The true absorption coefficients are indicated by the blank rectangles, and the reconstructed ones by the bars of different patterns. The absorption contrasts of the inclusions to the background are 2, 3, 5, and 7, respectively.


Table 1 illustrates the relative computation cost per iteration of the earlier algorithms, assuming that the execution time of the algorithm using t and c2 to be 100%. For μa -only reconstruction, it can be seen that the execution times of the featured-data algorithms using t , E and t , and R are almost comparable to that of the algorithm using E , while that using the full TR data is 50 times that of the algorithm using E . For simultaneous reconstruction of μa and μs , the advantage of the modified GPST algorithm using R is more evident, and its execution time is only 110 that of the featured-data algorithm using t and c2 . As for the full TR-data algorithm, although higher spatial resolution and quantitativeness ratio can be obtained,25 its application to 3-D reconstruction may be limited because of the tremendous time cost for simultaneous reconstruction of μa and μs .

Table 1

A comparison of the relative computation time per iteration among the algorithms using the featured data types and full TR data.

E ⟨t⟩ ⟨t⟩ and c2 E and ⟨t⟩ Rs Full TR data
μa -reconstruction only223010040251200
μa and μs simultaneous reconstruction263010066103000


Cross Talks between μa and μs Reconstruction

The spatial resolution and contrast of the reconstructed images of NIR-DOT can be affected by many factors, e.g., the number of the measured data, and the mesh density in the forward and inverse models.34 The absorption-scattering cross talk is another reason of the degradation of the reconstructed images. The absorption-scattering cross talk, which means the localized variations in absorption appear as the localized variations in scattering in the reconstructed image or vice versa, arises from both the intrinsic nonuniqueness, the information incompleteness in the dataset, and inefficiency of the numerical methods used for the reconstruction.35 Because the ranges of the unknown optical properties in DOT, e.g., absorption and the reduced scattering coefficients, and the refractive index if necessary, are generally different, a unique solution for the distributions of the optical properties may not be obtained from the measured dataset.29, 36 Even if only the absorption and reduced scattering coefficients are assigned to be the unknowns, absorption-scattering cross talk may also happen due to the imbalance between the μa - and μs (or κ )-related entries in the Jacobian matrices.36, 37

For evaluating the algorithms on suppressing the cross talk, we considered the aforementioned 2-D domain containing three inclusions.38 The inclusions have the identical diameters of 16mm but different optical properties. Inclusion 1 has a contrast of 3 to the background in both the absorption and scattering. Inclusions 2 and 3 have contrasts of 3 either in the scattering or in the absorption, respectively, as illustrated in the first row of Fig. 3 . The reconstructed images from the algorithms using the featured data types and the full TR data are shown in the other rows of Fig. 3. The following observations can be obtained from the reconstructed images. 1. The scattering-only inclusion (inclusion 2) generates false images in the absorption images for all the algorithms except those using t and c2 , i.e., those algorithms are prone to excessively enhance the absorption reconstruction because of the cross talk between μa and μs reconstruction. 2. Coexisting absorption and scattering contrasts such as these in inclusion 1 degrade the reconstruction of scattering. 3. As a whole, the algorithms using the intensity-related featured data types such as E , E , and t , and the full TR data exhibit an interparameter cross talk, while those using time-related data types, such as t , t , and c2 , and R , are less suffering from the cross talks. Thus, it is strongly recommended that one pay attention to the pseudo-image in reconstruction and the pronounced quantitativeness caused by cross talks when the intensity-related featured data types are employed.

Fig. 3

Reconstructed images from the numerically simulated data for evaluating the cross talk between the absorption and scattering reconstructions when the different featured data types and the full TR data were used in the algorithms.



Tomographic Imaging of In Vitro Tissues

The measurement system used in this work was a typical multichannel time-resolved system based on the time-correlated single photon counting (TCSPC) technique.19 The system contained three pulsed laser diodes emitting ultrashort light pulses with the pulse width of about 100ps at wavelengths 759, 799, and 834nm , respectively. The temporal resolution of the whole system was about 180ps , which is decided mainly by the transit temporal spread of the photon multiplier tubes. The fiber bundles had a coaxial structure with the inner core for the irradiation and the outer annulus for the detection. In the DOT experiment, 16 fibers were fixed around the object using a metal fiber holder at a cross-sectional plane with equal spacing, as shown in Fig. 4 . As one of the 16 fiber bundles was selected to irradiate the object, and the attenuators connected at the other end of the detection bundles were automatically first adjusted to make efficient use of the TCSPC dynamic range. Then the re-emissions from the surface of the object were simultaneously collected by all the detectors at the three wavelengths. As a result, a set of 16×16 temporal profiles for every wavelength were finally recorded.

Fig. 4

Fibers, fiber holder, and the arrangement of 16 optodes for the TR-DOT experiments.


The in vitro imaging experiment was carried out on a chicken leg. Considering the irregular shape of the chicken leg, it was immersed into an Intralipid solution contained in a cylinder, as shown in Fig. 5a , and a differential imaging scheme was employed, i.e., measurements were conducted before and after the immersion of the chicken leg into the solution. For simplicity, no refractive index mismatch was considered in the image reconstruction.

Fig. 5

In vitro imaging of a chicken leg: (a) schematic of the experimental setup showing the measuring plane, and (b) the reconstructed absorption and reduced scattering coefficient images at the measuring plane based on 2-D algorithms. The values of the optical properties are denoted by colors (upper row) and by colored heights (lower row).


The μa and μs images reconstructed with the modified GPST algorithm using R are shown in the left and right columns of Fig. 5b, respectively. The maximum and minimum in the μa and μs images are considered corresponding to the chicken bone. In the reconstructed μa image, the size of the chicken leg is about 32mm in the x direction and 28mm in the direction perpendicular to the x axis. The separation between the chicken leg centroid and the container center is estimated to be about 15mm . All these values conform to the experimental conditions. At 759nm , the μa of the chicken leg muscle and bone reconstructed in this work are estimated to be 0.022 and 0.027mm1 , respectively, and the μs are estimated to be 0.59 and 0.55mm1 , respectively. The larger μa and a smaller μs of the bone than those of the chicken muscle may be mainly caused by the fact that the imaged intact chicken bone has only a thin wall of compact bone and is filled with marrow inside the hollow space.38

The results show that the combination of the TR-DOT system with the simultaneous absorption and scattering reconstruction algorithm based on the modified GPST can qualitatively reveal the inner structure of the tissue and reasonably quantify the optical properties. Differential imaging was successful in eliminating the influences from the instrumental responses and background heterogeneities, as well as in discriminating the hard tissues from the soft ones.


Tomographic Imaging of In Vivo Tissues

The in vivo TR-DOT experiments were carried out on the human forearm and lower leg.


Forearm Diffuse Optical Tomography for Imaging the Anatomical Structure

A hand-gripping exercise was selected for stimulating blood oxygen saturation change in the forearm. 16 fibers were fixed on the forearm with a diameter of 62mm on the imaging plane, and the orientation of the fibers to the forearm is shown in Fig. 6a . During the experiment, an initial set of the temporal profiles was collected while the forearm was at rest. Soon after the initial measurement (also called the reference measurement), the subject was asked to start the hand-gripping exercise. 30s later, the task measurement started. After the optical measurement, a coregistered transverse MRI image was taken, as shown in Fig. 6b.

Fig. 6

NIR-DOT measurements on the right forearms with a diameter of 62mm on the imaging plane. (a) 16 fibers were arranged orderly in an anticlockwise direction when viewed from the body to the right hand. (b) The MRI image of the right forearm at the DOT imaging plane in the experiment.


Figure 7 shows the reconstructed absolute optical images measured at 759, 799, and 834nm for the subject in the rest and hand-gripping states, respectively, using the modified GPST algorithm. In the μa images, some higher μa regions are shown to correspond to some of the arteries, veins, and muscle. By comparing the μs image with the MRI image of the subject’s forearm, the spots with high μs are believed to correspond to the bones (the upper one is the radius). The observations indicate that the modified GPST algorithm with the featured data-type R has practical significance, since the simultaneous reconstruction of both μa and μs are necessary for effectively discovering the inner anatomical structures.

Fig. 7

Reconstructed absolute images of the absorption and reduced scattering coefficients ( μa and μs ) of a human forearm at states of rest (two left columns) and hand-gripping (two right columns). The three rows correspond to the images at wavelengths of 759, 799, and 834nm , respectively.


Table 2 compares the optical properties obtained by our NIR-DOT images with those reported in the literature,39, 40, 41, 42, 43, 44 where the arm tissue refers to the muscle and skin. The table shows that the reconstructed optical properties of the bone and arm tissue by our TR-DOT methodology are close to or within the ranges of those of bone and muscle/skin reported in the literature. However, the reconstructed μa and μs of blood vessels by our methodology are smaller and larger than those reported in the literature, respectively. These differences may be caused by the extremely high absorption contrast of blood vessels to the background, as discussed in Sec. 3.1.


Forearm Diffuse Optical Tomography for Imaging the Changes in the Blood Oxygenation

Theoretically, the hemodynamic (or optical) changes can be inferred from the direct subtraction of the individual reconstructions between the hand-gripping and reference states. This usually leads to large errors due to the tiny difference in the optical properties between the two states. In practice, Δμa between the two states is directly reconstructed using the differential image reconstruction scheme. In the reconstruction, a homogeneous background (with optical properties of μaB and μsB ) is employed as the initial values. Supposing that there is only the variation in absorption existing between the two states, we have45

Eq. 11

where Δμa and Δχ are the changes in the absorption coefficient and the measurement set, respectively.

Assuming that oxy-hemoglobin (HbO) and deoxy-hemoglobin (Hb) are the only chromophors in tissue, the physiological information or the changes in the concentrations of HbO and Hb between the two states can be obtained as follows:

Eq. 12

where λ1 and λ2 are the two wavelengths used in the experiment.

The images of the hemoglobin concentration changes in the forearm are shown in Fig. 8 . It can be seen that [ΔHb] and [ΔHbO] are slightly positive and negative in the most regions, respectively, meaning that the whole muscles were deoxygenated during the hand-gripping exercise. Judging from the increases in [Hb] (about 70to110μM ) and decreases in [HbO] (less than 50μM ), region A may correspond to some thin veins or the muscles responsible for the hand-gripping exercise since it has been reported that, by an intense exercise from a rest state, the oxygen saturation of venous blood may fall from 75 to 25%.46 One can also find that region B possesses an increased blood volume ([Hb]+[HbO]) (about 30μM ) and a decreased [HbO]. Thus, region B may represent some thick veins. Compared with the MRI image, region C is believed to be the image of some arteries that provide more blood with high oxygen saturation to compensate for oxygen consumption in muscle.

Fig. 8

Changes in the concentrations (in micromolars) of deoxy-hemoglobin [ΔHb] (left), oxy-hemoglobin [ΔHbO] (middle), and total hemoglobin [ΔHb]+[ΔHbO] (right) between the hand-gripping and rest states of the forearm.



Lower-Leg Diffuse Optical Tomography for Investigating the Size Effect of Targets

As shown in the images in Sec. 3.1, the quantitativeness of the reconstructed images decreases as the absorption contrast increases. It has also been shown that the size of the inclusions considerably influence the quantitativeness.44 Here, we define another measure for the evaluation of the reconstructed images, i.e., size ratio, which is the ratio of the representative size of an inclusion to that of the imaging area, to investigate the size effect of targets on the quantitativeness of the reconstruction. Because the size ratio of the blood vessels in the leg is relatively small, a lower leg was selected as a proper object for the evaluation.

The experiment was conducted on the right lower leg of a female subject. The imaging plane was located at about 13 of the height from her malleolus to the knee, where the size of the two bones and their separation present little variation.46 The diameters of the subject’s lower legs confined in the fiber holder were about 56mm .

Figure 9 shows the reconstructed absorption (μa) and scattering (μs) images at 759 and 834nm using the modified GPST algorithm, respectively. The two large areas with high μs in the scattering images can be assigned to the images of two bones, tibia and fibula, by comparing with the MRI images (not shown here) and the optical fiber orientation in the measurement. The separation between the two bones in the μs images is about 28mm , and the diameters of the tibia and fibula are estimated on average to be 19 and 12mm , respectively. These values coincide well with those measured from the relevant MRI image. The observations showed that the positions of the bones are properly revealed but the shapes are somewhat distorted. The μs value of the bones estimated from the images is in the range of 1.6to2.2mm1 . The other areas with high μs are assigned to the muscle that has approximately a μs value of 1.0to1.4mm1 . The fat and blood is estimated to have a μs value of 0.8to1.0mm1 in terms of the background in the μs images and could not be distinguished very clearly.

Fig. 9

The reconstructed absorption (μa) and scattering (μs) images of a subject’s lower legs at 759 and 834nm using the modified GPST algorithm.


By coregistering the anatomical structure from the MRI image to the μa images, the spots with the μa values as high as 0.04to0.043mm1 near the boundary are most probably the images of the short saphenous vein (S-vein) and the great saphenous vein (G-vein). However, the meaning of other high μa spots near the boundary is uncertain, and the arteries or veins deep inside the legs are unobservable in the μa images. The μa value of the leg tissues including blood estimated from the background of the μa image is around 0.035mm1 , which is much less than that observed in Fig. 7.

The degradation of Fig. 9 compared to Fig. 7 in the spatial resolution and quantitativeness of the μa reconstruction is partly attributed to the small size ratio of the blood vessels to the leg. The previous observation is qualitatively justified using the modified GPST and the full TR algorithms for a numerical 2-D phantom mimicking the lower leg, as shown in Fig. 10 . The circular domain has a diameter of 64mm and the background optical properties of μa=0.025mm1 and μs=1.0mm1 . Two large inner circular regions represent tibia and fibula bones in the lower leg. The bones are assumed to have μa=0.025mm1 (the same as the background) and μs=1.6mm1 . Some small circular regions denote the blood vessels, i.e., posterior tibial artery, short saphenous vein, and great saphenous vein in a leg. The contrast in μa of the blood vessels to the background is assumed to be 16, while the contrast in μs is 1.6.

Fig. 10

Numerical phantom for simulating the imaging plane of the lower leg. Two big circles imitate the bones and four small circles denote the blood vessels.


Figure 11 shows the reconstructed images for this numerical phantom. In the μs images, the reconstructed contrasts in the μs of the bones are 1.65 and 1.42 with the full TR-data algorithm and the modified GPST algorithm, respectively. In μa images, the absorbing inhomogeneities near the boundary (inhomogeneity 1, 3, and 4) are observable in each image using the full TR-data algorithm or the GPST algorithm, while those deep inside the domain are fairly disclosed only with the full TR-data algorithm. It is also noted even with the full TR algorithm that the quantitativeness of the μa reconstruction is merely 0.11, far below the true value. From these results, it is briefly deduced that the small size ratios of the targets might be the primary reason that causes the low spatial resolution and quantitativeness in the absorption images in Fig. 9. This means that the influences from the adverse performances of the algorithm, e.g., intrinsic ill-posedness, information incompleteness, and nonuniqueness, etc., on the image quality are overwhelmingly more pronounced than the experimental uncertainties, and it is more difficult to quantify a small-sized inclusion than a large-sized one.47 Further work has to be carried out for the improvement of the spatial resolution for reduction of the cross talk between the absorption and reduced scattering coefficients, and for quantitative characterization of heterogeneities. The effective ways verified so far include use of multigrid or adaptive meshing strategies, adoption of region-based reconstruction techniques,48 and incorporation of a-priori information from the established imaging modalities, such as MRI,29, 49 ultrasound imaging,50 and x-ray CT.51

Fig. 11

Reconstructed images of the simulated lower leg by the modified GPST algorithm and that using the full TR data, together with the target images.


Table 2

Optical properties of some arm tissues from the reconstructed images in this study and from the literature.

Tissue types λ (nm) μa (mm−1) μs′ (mm−1)
ArmReconstructed images759 to 8340.03 to 0.061.0 to 1.8
TissuesLiterature30, 31, 32, 33, 34700 to 9000.02 to 0.030.45 to 1.2
BloodReconstructed images759 to 8340.06 to 0.081.0 to 1.2
Literature32633 to 6850.13 to 0.40.4 to 0.6
BoneReconstructed images759 to 8340.03 to 0.041.8 to 2.8
Literature35800,8490.027, 0.02151.8, 0.91



This work summarizes part of the theoretical and experimental studies of our group on NIR-TR-DOT in the past several years. The major achievements include the development and demonstration of the two image reconstruction algorithms—the modified GPST algorithm and the full TR-data one. Their performances are quantitatively investigated by simulations and further validated with in vitro and in vivo DOT experiments, combined with a laboratory-equipped TCSPC-based multichannel TR system. It can be concluded that after several years of laboratory study on NIR-TR-DOT technology, high-quality DOT images can be obtained by developing a reconstruction algorithm that can effectively incorporate the time-resolved information. Studies toward the clinical applications of NIR-TR-DOT, such as optical mammography and preterm infant brain imaging, are now being carried out and will be reported in the near future.


Huijuan Zhao and Feng Gao thank the support from the National Natural Science Foundation of China (60578008, 60678049), and the National Basic Research Program of China (2006CB705700). The authors also acknowledge the fellowship from New Energy Development Organization (NEDO), Japan, for fulfilling part of the work.



D. Grosenick, K. T. Moesta, H. Wabnitz, J. Mucke, C. Striszczynski, R. Macdonald, P. M. Schlag, and H. Rinneberg, “Time-domain optical mammography: initial clinical results on detection and characterization of breast tumors,” Appl. Opt., 42 3170 –3148 (2003). 0003-6935 Google Scholar


A. Rice and C. M. Quinn, “Angiogenesis, thrombospodin, and ductal carcinoma in situ of the breast,” J. Clin. Pathol., 55 569 –574 (2002). 0021-9746 Google Scholar


D. A. Boas, D. H. Brooks, C. A. Dimarzio, M. Kilmer, R. J. Gaudette, Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag., 11 57 –1175 (2001). 1053-5888 Google Scholar


G. Strangman, D. A Boas, and J. Sutton, “Non-invasive neuroimaging using near infrared light,” Biol. Psychiatry, 52 676 –693 (2002). 0006-3223 Google Scholar


J. C. Hebden, A. P. Gibson, R. M. Yusof, N. Everdell, E. M. C. Hillman, D. T. Delpy, S. R. Arridge, T. Austin, J. H. Meek, and J. S. Wyatt, “Three-dimensional optical tomography of the premature infant brain,” Phys. Med. Biol., 47 4155 –4166 (2002). 0031-9155 Google Scholar


M. Ferrari, L. Mottola, and V. Quaresima, “Principles, techniques, and limitations of near infrared spectroscopy,” Can. J. Appl. Physiol., 29 463 –487 (2004). 1066-7814 Google Scholar


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 (7), 1117 –1130 (2001). 0031-9155 Google Scholar


F. Gao, H. Zhao, Y. Onodera, A. Sassaroli, Y. Tanikawa, and Y. Yamada, “Image reconstruction from experimental measurements of an multichannel time resolved optical tomographic imaging system,” Proc. SPIE, 4250 351 –361 (2001). 0277-786X Google Scholar


A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol., 50 R1 –R41 (2005). 0031-9155 Google Scholar


R. L. Barbour, H. L. Graber, Y. Pei, S. Zhong, and C. H. Schmitz, “Optical tongraphic imaging of dynamic features of dense-scattering media,” J. Opt. Soc. Am. A, 14 3018 –3036 (2001). 0740-3232 Google Scholar


A. Corlu, R. Choe, T. Durduran, K. Lee, M. Schweiger, S. R. Arridge, E. M. C. Hillman, and A. G. Yodh, “Diffuse optical tomography with spectral constraints and wavelength optimization,” Appl. Opt., 44 (11), 2082 –2093 (2005). 0003-6935 Google Scholar


A. Li, Q. Zhang, J. P. Culver, E. L. Miller, and D. A. Boas, “Reconstructing chromosphere concentration images directly by continuous-wave diffuse optical tomography,” Opt. Lett., 29 (3), 256 –258 (2004). 0146-9592 Google Scholar


R. L. Barbour, C. H. Schmitz, D. P. Klemer, Y. Pei, and H. L. Graber, “Design and initial testing of system for simultaneous bilateral dynamic optical tomographic mammography,” WD4 (2004). Google Scholar


R. J. Grable, D. P. Rohler, and K. L. A. Sastry, “Optical tomography breast imaging,” Proc. SPIE, 2979 197 –210 (2004). 0277-786X Google Scholar


S. Srinivasan, B. W. Pogue, S. Jiang, H. Dehghani, and K. D. Paulsen, “Spectrally constrained chromophore and scattering near-infrared tomography provides quantitative and robust reconstruction,” Appl. Opt., 44 (10), 1858 –1869 (2005). 0003-6935 Google Scholar


G. Yu, T. Durduran, D. Furuya, J. H. Greenberg, and A. G. Yodh, “Frequency-domain multiplexing system for in vivo diffuse light measurements of rapid cerebral hemodynamics,” Appl. Opt., 42 (16), 2931 –2939 (2003). 0003-6935 Google Scholar


T. H. Pham, T. Spott, L. O. Svaasand, and B. J. Tromberg, “Quantifying the properties of two-layer turbid media with frequency-domain diffuse reflectance,” Appl. Opt., 39 (25), 4733 –4745 (2000). 0003-6935 Google Scholar


H. Eda, I. Oda, Y. Ito, Y. Wada, Y. Oikawa, Y. Tsunazawa, M. Takada, Y. Tsuchiya, Y. Yamashita, M. Oda, A. Sassaroli, Y. Yamada, and M. Tamura, “Multichannel time-resolved optical tomographic imaging system,” Rev. Sci. Instrum., 70 (9), 3595 –3601 (1999). 0034-6748 Google Scholar


A. P. Gibson, T. Austin, N. L. Everdell, M. Schweiger, S. R. Arridge, J. H. Meek, J. S. Wyatt, D. T. Delpy, and J. C. Hebden, “Three-dimensional whole-head optical tomography of passive motor evoked responses in the neonate,” Neuroimage, 30 521 –528 (2006). 1053-8119 Google Scholar


S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol., 37 1531 –1559 (1992). 0031-9155 Google Scholar


S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl., 15 R41 –R93 (1999). 0266-5611 Google Scholar


R. J. Gaudette, D. H. Broools, C. A. DiMarzio, M. Kilmer, E. L. Miller, T. Gaudette, and D. A. Boas, “A comparison study of linear reconstruction techniques for diffuse optical tomographic imaging of absorption coefficient,” Phys. Med. Biol., 45 1051 –1070 (2000). 0031-9155 Google Scholar


M. Schweiger, A. P. Gibson, and S. R. Arridge, “Computational aspects of diffuse optical tomography,” IEEE Comput. Sci. Eng., 5 33 –41 (2003). 1070-9924 Google Scholar


H. Jiang, N. Iftimia, J. Eggert, L. Fajardo, and K. Klove, “Near-infrared optical imaging of the breast with model-based reconstruction,” Acad. Radiol., 9 (2), 186 –194 (2002). 1076-6332 Google Scholar


F. Gao, Y. Tanikawa, H. Zhao, and Y. Yamada, “A semi three-dimensional algorithm for time-resolved diffuse optical tomography by use of generalized pulse spectrum technique,” Appl. Opt., 41 (34), 7346 –7358 (2002). 0003-6935 Google Scholar


F. Gao, P. Poulet, and Y. Yamada, “Simultaneous mapping of absorption and scattering coefficients from a three-dimensional model of time-resolved optical tomography,” Appl. Opt., 39 (31), 5898 –5910 (2000). 0003-6935 Google Scholar


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). 0031-9155 Google Scholar


E. M. C. Hillman, “Experimental and theoretical investigations of near infrared tomographic imaging methods and clinical applications,” Univ. College London, (2002). Google Scholar


M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol., 44 1699 –1717 (1999). 0031-9155 Google Scholar


J. C. Ye, K. J. Webb, R. P. Millane, and T. J. Downar, “Modified distorted Born iterative method with an approximate Frechet derivative for optical diffusion tomography,” J. Opt. Soc. Am. A, 16 1814 –1826 (1999). 0740-3232 Google Scholar


F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, “Time-resolved diffuse optical tomography using a modified generalized pulse spectrum technique,” IEICE Trans. Inf. Syst., E85-D (1), 133 –142 (2002). 0916-8532 Google Scholar


F. Gao, H. Zhao, and Y. Yamada, “Improvement of image quality in diffuse optical tomography by use of full time-resolved data,” Appl. Opt., 41 (4), 778 –791 (2002). 0003-6935 Google Scholar


B. W. Pogue, S. C. Davis, X. Song, B. Brooksby, H. Dehghani, and K. D. Paulsen, “Image analysis methods for diffuse optical tomography,” J. Biomed. Opt., 11 (3), 033001 (2006). 1083-3668 Google Scholar


P. K. Yalavarthy, H. Dehghani, B. W. Pogue, and K. D. Paulsen, “Critical computational aspects of near infrared circular tomographic imaging: analysis of measurement number, mesh resolution and reconstruction basis,” Opt. Express, 14 6113 –6127 (2006). 1094-4087 Google Scholar


Y. Pei, H. L. Graber, and R. L. Barbour, “Normalized-constraint algorithm for minimizing inter-parameter crosstalk in DC optical tomography,” Opt. Express, 11 (2), 97 –109 (2001). 1094-4087 Google Scholar


P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol., 50 (4), R1 –R43 (2005). 0031-9155 Google Scholar


Y. Xu, X. Gu, T. Khan, and H. Jiang, “Absorption and scattering images of heterogeneous scattering media can be simultaneously reconstructed by use of dc data,” Appl. Opt., 41 (25), 5427 –5437 (2002). 0003-6935 Google Scholar


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


C. R. Simpson, M. Kohl, M. Essenpreis, and M. Cope, “Near infrared optical properties of ex vivo human skin and subcutaneous tissues measured using the Monte Carlo inversion technique,” Phys. Med. Biol., 43 2465 –2478 (1998). 0031-9155 Google Scholar


W. Cheong, S. A. Prahl, and A. J. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron., 26 (12), 21662184 (1990). 0018-9197 Google Scholar


G. Zaccanti, A. Taddeucci, M. Barilli, P. Bruscaglioni, and F. Martilli, “Optical properties of biological tissues,” Proc. SPIE, 2389 513 –521 (1995). 0277-786X Google Scholar


S. J. Matcher, M. Cope, and D. T. Delpy, “In vivo measurements of the wavelength dependence of tissue-transport scattering coefficients between 760 and 900nm measured with time-resolved spectroscopy,” Appl. Opt., 36 386 –396 (1997). 0003-6935 Google Scholar


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


M. Firbank, M. Hiraoka, M. Essenpreis, and D. T. Delpy, “Measurement of the optical properties of the skull in the wavelength range 650950nm,” Phys. Med. Biol., 38 503 –510 (1993). 0031-9155 Google Scholar


H. Zhao, F. Gao, Y. Tanikawa, K. Homma, and Y. Yamada, “Time-resolved diffuse optical tomographic imaging for the provision of both anatomical and functional information about biological tissue,” Appl. Opt., 44 (10), 1906 –1916 (2005). 0003-6935 Google Scholar


L. A. Paunescu, “Tissue blood flow and oxygen consumption measured with near-infrared frequency-domain spectroscopy,” (2001). Google Scholar


F. Gao, H. Zhao, Y. Tanikawa, K. Homma, and Y. Yamada, “Influences of target size and contrast on near infrared diffuse optical tomography—a comparison between featured-data and full time-resolved schemes,” Opt. Quantum Electron., 37 (13–15), 1287 –1304 (2005). 0306-8919 Google Scholar


S. Srinivasan, B. W. Pogue, H. Dehghani, S. Jiang, X. Song, and K. D. Paulsen, “Improved quantification of small objects in near-infrared diffuse optical tomography,” J. Biomed. Opt., 9 (6), 1161 –1171 (2004). 1083-3668 Google Scholar


B. Brooksby, S. Jiang, H. Dehghani, B. W. Pogue, K. D. Paulsen, J. Weaver, C. Kogel, and S. P. Poplack, “Combining near-infrared tomography and magnetic resonance imaging to study in vivo breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” J. Biomed. Opt., 10 (5), 051504 (2005). 1083-3668 Google Scholar


Q. Zhu, S. H. Kurtzma, P. Hegde, S. Tannenbaum, M. Kane, M. Huang, N. G. Chen, B. Jagjivan, and K. Zarfos, “Utilizing optical tomography with ultrasound localization to image heterogeneous hemoglobin distribution in large breast cancers,” Neoplasia, 7 (3), 263 –270 (2005). 1522-8002 Google Scholar


Q. Zhang, T. J. Brukilacchio, A. Li, J. J. Stott, T. Chaves, E. Hillman, T. Wu, M. Chorlton, E. Rafferty, R. H. Moore, D. B. Kopans, and D. A. Boas, “Coregistered tomographic x-ray and optical breast imaging: initial results,” J. Biomed. Opt., 10 024033 (2005). 1083-3668 Google Scholar
©(2007) Society of Photo-Optical Instrumentation Engineers (SPIE)
Huijuan Zhao, Feng Gao, Yukari Tanikawa, and Yukio Yamada "Time-resolved diffuse optical tomography and its application to in vitro and in vivo imaging," Journal of Biomedical Optics 12(6), 062107 (1 November 2007).
Published: 1 November 2007

Back to Top