## 1.

## Introduction

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
$\left({\mu}_{a}\right)$
and reduced scattering coefficient
$\left({\mu}_{s}\right)$
, 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
${P}_{1}$
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.

## 2.

## Image Reconstruction Algorithms

## 2.1.

### General Framework

As described before, given a forward model $\mathbf{F}$ with a vector of the optical properties $\mathbf{p}$ , including both the absorption coefficient ${\mu}_{a}$ and diffusion coefficient $\kappa =1\u2215\left(3{\mu}_{s}^{\prime}\right)$ (or scattering coefficient ${\mu}_{s}^{\prime}$ ), the forward problem can be expressed as $\mathbf{y}=\mathbf{F}\left(\mathbf{p}\right)$ , where $\mathbf{y}$ denotes a vector of the measured quantities, and the inverse problem is accordingly written as $\mathbf{p}={\mathbf{F}}^{-1}\left(\mathbf{y}\right)$ . If the actual optical properties $\mathbf{p}$ are close to its initial estimation ${\mathbf{p}}_{0}$ , then the nonlinear inverse problem posed by DOT can be linearized as

## Eq. 1

$$\Delta \mathbf{y}\approx {\mathbf{F}}^{\prime}\left({\mathbf{p}}_{0}\right)\Delta \mathbf{p},$$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

$$[\mathbf{\chi}-\mathbf{F}\left({\mathbf{p}}_{i}\right)]=\mathbf{J}\left({\mathbf{p}}_{i}\right)\delta {\mathbf{p}}_{i},\phantom{\rule{1em}{0ex}}{\mathbf{p}}_{i+1}={\mathbf{p}}_{i}+\delta {\mathbf{p}}_{i},$$^{26}

## 2.2.

### 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)
$\u27e8t\u27e9$
, and high-order central moments such as the variance
${c}_{2}$
, 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
$\u27e8t\u27e9$
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.

## 2.2.1.

#### 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
$\chi ({\xi}_{d},{\zeta}_{s},t)$
in a tissue volume
$\Omega $
with a boundary
$\partial \Omega $
, the Laplace transform of
$\chi ({\xi}_{d},{\zeta}_{s},t)$
is given by

## Eq. 3

$$\widehat{\chi}({\xi}_{d},{\zeta}_{s},f)={\int}_{0}^{+\propto}\chi ({\xi}_{d},{\zeta}_{s},t)\mathrm{exp}(-ft)dt,$$^{32}

The ratio of the Laplace-transformed signals at two distinct frequencies (denoted by ${f}_{1}$ and ${f}_{2}$ ) is employed as a data type:

## Eq. 4

$$R({\xi}_{d},{\zeta}_{s},{f}_{1},{f}_{2})=\frac{\widehat{\chi}({\xi}_{d},{\zeta}_{s},{f}_{2})}{\widehat{\chi}({\xi}_{d},{\zeta}_{s},{f}_{1})}.$$Thus, the Jacobian matrix for the data-type
$R$
can be reformularized according to those for
$\widehat{\chi}({\xi}_{d},{\zeta}_{s},{f}_{1})$
and
$\widehat{\chi}({\xi}_{d},{\zeta}_{s},{f}_{2})$
:^{25}

## Eq. 5

$${\widehat{\mathbf{J}}}_{p}({\xi}_{d},{\zeta}_{s},{f}_{1},f)=\frac{{\mathbf{J}}_{p}({\xi}_{d},{\zeta}_{s},{f}_{2})}{\widehat{\chi}({\xi}_{d},{\zeta}_{s},{f}_{1})}-\mathbf{R}({\xi}_{d},{\zeta}_{s},{f}_{1},{f}_{2})\frac{{\mathbf{J}}_{p}({\xi}_{d},{\zeta}_{s},{f}_{1})}{\widehat{\chi}({\xi}_{d},{\zeta}_{s},{f}_{1})}.$$## 2.2.2.

#### Modified generalized pulse spectrum technique

For the simultaneous reconstruction of
${\mu}_{a}$
and
$\kappa $
, 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
$\chi $
, which is much more difficult and time consuming than for the case of
${\mu}_{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
${\mu}_{s}^{\prime}$
-regarding Jacobian matrix
${\mathbf{J}}_{{\mu}_{s}^{\prime}}$
, 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

$${\widehat{\mathbf{J}}}_{p}({\xi}_{d},{\zeta}_{s},{f}_{1},{f}_{2})=[\frac{{\mathbf{J}}_{{\mu}_{a}}({\xi}_{d},{\zeta}_{s},{f}_{2})}{\widehat{\chi}({\xi}_{d},{\zeta}_{s},{f}_{1})}{C}_{1}-R({\xi}_{d},{\zeta}_{s},{f}_{1},{f}_{2})\frac{{\mathbf{J}}_{{\mu}_{a}}({\xi}_{d},{\zeta}_{s},{f}_{1})}{\widehat{\chi}({\xi}_{d},{\zeta}_{s},{f}_{1})}{C}_{2}],$$## Eq. 7

$$[{C}_{1},{C}_{2}]=\{\begin{array}{ll}[1,1]& p\u220a{\mu}_{a}\\ \{3\kappa [{\mu}_{a}\left(\mathbf{r}\right)+{f}_{2}\u2215c],3\kappa [{\mu}_{a}\left(\mathbf{r}\right)+{f}_{1}\u2215c]\}& p\u220a{\mu}_{s}^{\prime}\end{array}\phantom{\}}.$$It can be seen that, with the modified GPST, the simultaneous reconstruction of ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ is computationally as efficient as the ${\mu}_{a}$ -only reconstruction from the point of view of the Jacobian matrix computation.

## 2.2.3.

#### Scaling of the Jacobian matrix

In general, the entries of ${\mathbf{J}}_{{\mu}_{a}}^{\left(R\right)}$ are several orders of magnitudes larger than those of ${\mathbf{J}}_{{\mu}_{s}^{\prime}}^{\left(R\right)}$ , resulting in the ${\mu}_{a}$ -reconstruction more pronounced and the simultaneous reconstruction of both the ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ being unbalanced. To enhance the efficiency of the modified GPST algorithm in ${\mu}_{s}^{\prime}$ -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

$$\widehat{\mathbf{J}}\left(i\right)=\left[\frac{{\widehat{\mathbf{J}}}_{{\mu}_{a}}\left(i\right)}{{w}_{{\mu}_{a}}\left(i\right)}\phantom{\rule{0.3em}{0ex}}\frac{{\widehat{\mathbf{J}}}_{{\mu}_{s}^{\prime}}\left(i\right)}{{w}_{{\mu}_{s}^{\prime}}\left(i\right)}\right],$$## 2.3.

### 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 $\delta {\mu}_{a}$ and $\delta \kappa $ to ${\mu}_{a}$ and $\kappa $ , respectively, one can calculate the perturbation of the photon density $\delta \Phi $ from the DE and obtain $\delta \chi $ by taking into account the relation between $\Phi $ and $\chi $ . Then the Jacobian matrix at the time ${t}_{m}$ can be formalized as:

## Eq. 9

$$\{\begin{array}{l}{\mathbf{J}}_{{\mu}_{a}}^{\left(\chi \right)}({\xi}_{d},{\zeta}_{s},{t}_{m})={\left[\begin{array}{c}-{\int}_{0}^{{t}_{m}}{\int}_{\Omega}\chi ({\xi}_{d},\mathbf{r},{t}_{m}-\tau )\Phi (\mathbf{r},{\zeta}_{s},\tau )c{u}_{1}\left(\mathbf{r}\right)d\mathbf{r}d\tau \\ \vdots \\ -{\int}_{0}^{{t}_{m}}{\int}_{\Omega}\chi ({\xi}_{d},\mathbf{r},{t}_{m}-\tau )\Phi (\mathbf{r},{\zeta}_{s},\tau )c{u}_{N}\left(\mathbf{r}\right)d\mathbf{r}d\tau \end{array}\right]}^{T}\\ {\mathbf{J}}_{\kappa}^{\left(\chi \right)}({\xi}_{d},{\zeta}_{s},{t}_{m})={\left[\begin{array}{c}-{\int}_{0}^{{t}_{m}}{\int}_{\Omega}{\nabla}_{\mathbf{r}}\chi ({\xi}_{d},\mathbf{r},{t}_{m}-\tau )\bullet {\nabla}_{\mathbf{r}}\Phi (\mathbf{r},{\zeta}_{s},\tau ){u}_{1}\left(\mathbf{r}\right)d\mathbf{r}d\tau \\ \vdots \\ -{\int}_{0}^{{t}_{m}}{\int}_{\Omega}{\nabla}_{\mathbf{r}}\chi ({\xi}_{d},\mathbf{r},{t}_{m}-\tau )\bullet {\nabla}_{\mathbf{r}}\Phi (\mathbf{r},{\zeta}_{s},\tau ){u}_{N}\left(\mathbf{r}\right)d\mathbf{r}d\tau \end{array}\right]}^{T}\end{array}\phantom{\}},$$• 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 $\overline{\chi}({\xi}_{d},{\zeta}_{s},t)=\chi ({\xi}_{d},{\zeta}_{s},t)\u2215E({\xi}_{d},{\zeta}_{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 $\Delta t$ , and integrate the data in this time slice. Then one can obtain a data for the $m$ ’th time slice,

$${\overline{\chi}}^{\left(\Delta t\right)}({\xi}_{d},{\zeta}_{s},{t}_{m})={\int}_{0}^{\Delta t}\overline{\chi}({\xi}_{d},{\zeta}_{s},{t}_{m}+\tau )d\tau ,$$as illustrated in Fig. 1 .^{31}

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

## Eq. 10

$${\widehat{\mathbf{J}}}_{p}^{\left(\Delta t\right)}({\xi}_{d},{\zeta}_{s},t)=\frac{{\mathbf{J}}_{p}^{\left(\Delta t\right)}({\xi}_{d},{\zeta}_{s},t)}{E({\xi}_{d},{\zeta}_{s})}-{\overline{\chi}}^{\left(\Delta t\right)}({\xi}_{d},{\zeta}_{s},t)\frac{{\mathbf{J}}_{p}^{\left(E\right)}(\xi ,{\zeta}_{s})}{E({\xi}_{d},{\zeta}_{s})}\phantom{\rule{1em}{0ex}}p\u220a[{\mu}_{a},{\mu}_{s}^{\prime}],$$## 3.

## 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$
,
$\u27e8t\u27e9$
,
$E$
and
$\u27e8t\u27e9$
,
$\u27e8t\u27e9$
and
${c}_{2}$
,
$R$
, and on the full TR data, are evaluated on the aspects of the quantitativeness in
${\mu}_{a}$
reconstruction, the execution time, and the cross talk between
${\mu}_{a}$
and
${\mu}_{s}^{\prime}$
(
${\mu}_{s}^{\prime}=3\u2215\kappa $
is used instead of
$\kappa $
, 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
$80\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
and the background optical properties of
${\mu}_{a}=0.01\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
and
${\mu}_{s}^{\prime}=1.0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
. 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
$\left(20\phantom{\rule{0.3em}{0ex}}\mathrm{ps}\right)$
to ensure the numerical accuracy.^{31}

## 3.1.

### Quantitativeness and Execution Time in ${\mu}_{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 $16\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ was employed. The two inclusions are positioned symmetrically to the center with a separation of $40\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . Figure 2 shows the reconstructed ${\mu}_{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 ${\mu}_{a}$ , and the quantitativeness ratio reaches as high as 85%. As the contrast increases, the peak value of reconstructed ${\mu}_{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.

Table 1
illustrates the relative computation cost per iteration of the earlier algorithms, assuming that the execution time of the algorithm using
$\u27e8t\u27e9$
and
${c}_{2}$
to be 100%. For
${\mu}_{a}$
-only reconstruction, it can be seen that the execution times of the featured-data algorithms using
$\u27e8t\u27e9$
,
$E$
and
$\u27e8t\u27e9$
, 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
${\mu}_{a}$
and
${\mu}_{s}^{\prime}$
, the advantage of the modified GPST algorithm using
$R$
is more evident, and its execution time is only
$1\u221510$
that of the featured-data algorithm using
$\u27e8t\u27e9$
and
${c}_{2}$
. 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
${\mu}_{a}$
and
${\mu}_{s}^{\prime}$
.

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

${\mu}_{a}$ -reconstruction only | 22 | 30 | 100 | 40 | 25 | 1200 |

${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ simultaneous reconstruction | 26 | 30 | 100 | 66 | 10 | 3000 |

## 3.2.

### Cross Talks between ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ 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
${\mu}_{a}$
- and
${\mu}_{s}^{\prime}$
(or
$\kappa $
)-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
$16\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
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
$\u27e8t\u27e9$
and
${c}_{2}$
, i.e., those algorithms are prone to excessively enhance the absorption reconstruction because of the cross talk between
${\mu}_{a}$
and
${\mu}_{s}^{\prime}$
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
$\u27e8t\u27e9$
, and the full TR data exhibit an interparameter cross talk, while those using time-related data types, such as
$\u27e8t\u27e9$
,
$\u27e8t\u27e9$
, 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.

## 4.

## Tomographic Imaging of *In Vitro* Tissues

*In Vitro*

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
$100\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
at wavelengths 759, 799, and
$834\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, respectively. The temporal resolution of the whole system was about
$180\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
, 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\times 16$
temporal profiles for every wavelength were finally recorded.

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.

The
${\mu}_{a}$
and
${\mu}_{s}^{\prime}$
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
${\mu}_{a}$
and
${\mu}_{s}^{\prime}$
images are considered corresponding to the chicken bone. In the reconstructed
${\mu}_{a}$
image, the size of the chicken leg is about
$32\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
in the
$x$
direction and
$28\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
in the direction perpendicular to the
$x$
axis. The separation between the chicken leg centroid and the container center is estimated to be about
$15\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
. All these values conform to the experimental conditions. At
$759\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, the
${\mu}_{a}$
of the chicken leg muscle and bone reconstructed in this work are estimated to be 0.022 and
$0.027\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
, respectively, and the
${\mu}_{s}^{\prime}$
are estimated to be 0.59 and
$0.55\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
, respectively. The larger
${\mu}_{a}$
and a smaller
${\mu}_{s}^{\prime}$
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.

## 5.

## Tomographic Imaging of *In Vivo* Tissues

*In Vivo*

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

## 5.1.

### 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 $62\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ 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. $30\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ later, the task measurement started. After the optical measurement, a coregistered transverse MRI image was taken, as shown in Fig. 6b.

Figure 7 shows the reconstructed absolute optical images measured at 759, 799, and $834\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ for the subject in the rest and hand-gripping states, respectively, using the modified GPST algorithm. In the ${\mu}_{a}$ images, some higher ${\mu}_{a}$ regions are shown to correspond to some of the arteries, veins, and muscle. By comparing the ${\mu}_{s}^{\prime}$ image with the MRI image of the subject’s forearm, the spots with high ${\mu}_{s}^{\prime}$ 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 ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ are necessary for effectively discovering the inner anatomical structures.

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
${\mu}_{a}$
and
${\mu}_{s}^{\prime}$
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.

## 5.2.

### 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,
$\Delta {\mu}_{a}$
between the two states is directly reconstructed using the differential image reconstruction scheme. In the reconstruction, a homogeneous background (with optical properties of
${\mu}_{aB}$
and
${\mu}_{sB}^{\prime}$
) is employed as the initial values. Supposing that there is only the variation in absorption existing between the two states, we have^{45}

## Eq. 11

$$\Delta \chi \bullet \chi ({\mu}_{aB,},{\mu}_{sB}^{\prime})=\mathbf{F}\left({\mu}_{a}\right)+{\mathbf{J}}_{{\mu}_{a}}\left({\mu}_{a}\right)\Delta {\mu}_{a},$$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

$$\{\begin{array}{l}\Delta \left[\mathrm{Hb}\right]=-\frac{\Delta {\mu}_{a}^{{\lambda}_{2}}{\epsilon}_{\mathrm{Hb}\mathrm{O}}^{{\lambda}_{1}}-\Delta {\mu}_{a}^{{\lambda}_{1}}{\epsilon}_{\mathrm{Hb}\mathrm{O}}^{{\lambda}_{2}}}{{\epsilon}_{\mathrm{Hb}}^{{\lambda}_{1}}{\epsilon}_{\mathrm{Hb}\mathrm{O}}^{{\lambda}_{2}}-{\epsilon}_{\mathrm{Hb}}^{{\lambda}_{2}}{\epsilon}_{\mathrm{Hb}\mathrm{O}}^{{\lambda}_{1}}}\\ \Delta \left[\mathrm{Hb}\mathrm{O}\right]=\frac{\Delta {\mu}_{a}^{{\lambda}_{2}}{\epsilon}_{\mathrm{Hb}}^{{\lambda}_{1}}-\Delta {\mu}_{a}^{{\lambda}_{1}}{\epsilon}_{\mathrm{Hb}}^{{\lambda}_{2}}}{{\epsilon}_{\mathrm{Hb}}^{{\lambda}_{1}}{\epsilon}_{\mathrm{Hb}\mathrm{O}}^{{\lambda}_{2}}-{\epsilon}_{\mathrm{Hb}}^{{\lambda}_{2}}{\epsilon}_{\mathrm{Hb}\mathrm{O}}^{{\lambda}_{1}}}\end{array}\phantom{\}},$$The images of the hemoglobin concentration changes in the forearm are shown in Fig. 8
. It can be seen that
$\left[\Delta \mathrm{Hb}\right]$
and
$\left[\Delta \mathrm{Hb}\mathrm{O}\right]$
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
$70\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}110\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$
) and decreases in [HbO] (less than
$-50\phantom{\rule{0.3em}{0ex}}\mu \mathrm{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
$(\left[\mathrm{Hb}\right]+\left[\mathrm{Hb}\mathrm{O}\right])$
(about
$30\phantom{\rule{0.3em}{0ex}}\mu \mathrm{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.

## 5.3.

### 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
$1\u22153$
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
$56\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
.

Figure 9 shows the reconstructed absorption $\left({\mu}_{a}\right)$ and scattering $\left({\mu}_{s}^{\prime}\right)$ images at 759 and $834\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ using the modified GPST algorithm, respectively. The two large areas with high ${\mu}_{s}^{\prime}$ 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 ${\mu}_{s}^{\prime}$ images is about $28\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , and the diameters of the tibia and fibula are estimated on average to be 19 and $12\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , 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 ${\mu}_{s}^{\prime}$ value of the bones estimated from the images is in the range of $1.6\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}2.2\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ . The other areas with high ${\mu}_{s}^{\prime}$ are assigned to the muscle that has approximately a ${\mu}_{s}^{\prime}$ value of $1.0\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1.4\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{\u20131}$ . The fat and blood is estimated to have a ${\mu}_{s}^{\prime}$ value of $0.8\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1.0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{\u20131}$ in terms of the background in the ${\mu}_{s}^{\prime}$ images and could not be distinguished very clearly.

By coregistering the anatomical structure from the MRI image to the ${\mu}_{a}$ images, the spots with the ${\mu}_{a}$ values as high as $0.04\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}0.043\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{\u20131}$ 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 ${\mu}_{a}$ spots near the boundary is uncertain, and the arteries or veins deep inside the legs are unobservable in the ${\mu}_{a}$ images. The ${\mu}_{a}$ value of the leg tissues including blood estimated from the background of the ${\mu}_{a}$ image is around $0.035\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{\u20131}$ , 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 ${\mu}_{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 $64\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ and the background optical properties of ${\mu}_{a}=0.025\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=1.0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ . Two large inner circular regions represent tibia and fibula bones in the lower leg. The bones are assumed to have ${\mu}_{a}=0.025\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ (the same as the background) and ${\mu}_{s}^{\prime}=1.6\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ . 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 ${\mu}_{a}$ of the blood vessels to the background is assumed to be 16, while the contrast in ${\mu}_{s}^{\prime}$ is 1.6.

Figure 11
shows the reconstructed images for this numerical phantom. In the
${\mu}_{s}^{\prime}$
images, the reconstructed contrasts in the
${\mu}_{s}^{\prime}$
of the bones are 1.65 and 1.42 with the full TR-data algorithm and the modified GPST algorithm, respectively. In
${\mu}_{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
${\mu}_{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}

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

Arm | Reconstructed images | 759 to 834 | 0.03 to 0.06 | 1.0 to 1.8 |

Tissues | Literature^{30, 31, 32, 33, 34} | 700 to 900 | 0.02 to 0.03 | 0.45 to 1.2 |

Blood | Reconstructed images | 759 to 834 | 0.06 to 0.08 | 1.0 to 1.2 |

Literature^{32} | 633 to 685 | 0.13 to 0.4 | 0.4 to 0.6 | |

Bone | Reconstructed images | 759 to 834 | 0.03 to 0.04 | 1.8 to 2.8 |

Literature^{35} | 800,849 | 0.027, 0.0215 | 1.8, 0.91 |

## 6.

## Summary

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.

## Acknowledgments

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.

## References

*in vivo*diffuse light measurements of rapid cerebral hemodynamics,” Appl. Opt., 42 (16), 2931 –2939 (2003). https://doi.org/10.1364/AO.42.002931 0003-6935 Google Scholar

*in vitro*chicken leg by use of time-resolved near infrared optical tomography,” Phys. Med. Biol., 47 (11), 1979 –1993 (2002). https://doi.org/10.1088/0031-9155/47/11/310 0031-9155 Google Scholar

*In vivo*measurements of the wavelength dependence of tissue-transport scattering coefficients between 760 and $900\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ measured with time-resolved spectroscopy,” Appl. Opt., 36 386 –396 (1997). 0003-6935 Google Scholar

*In vivo*optical characterization of human tissues from $610\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1010\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ by time-resolved reflectance spectroscopy,” Phys. Med. Biol., 46 2227 –2237 (2001). https://doi.org/10.1088/0031-9155/46/8/313 0031-9155 Google Scholar

*in vivo*breast tissue: implementation of a Laplacian-type regularization to incorporate magnetic resonance structure,” J. Biomed. Opt., 10 (5), 051504 (2005). https://doi.org/10.1117/1.2098627 1083-3668 Google Scholar