1. Introduction

In recent years, studies on near-infrared (NIR) optical imaging as a potential biomedical diagnostic modality have attracted increasing attention [1–7]. Unlike conventional optical imaging that uses endogenous contrast of tissue optical properties between normal and diseased tissues to detect and/or diagnose diseases, optical fluorescence imaging utilizes exogenous fluorescent agents to enhance and/or characterize optical contrast between normal and diseased tissues. The use of fluorescence techniques has appealing features like selective absorption spectrum of contrast agents, distinct fluorescence spectrum and lifetime, cellular and molecular fluorescence probes, etc. These traits can be used to provide tissue functional information, cellular and molecular information, as well as structural information.

Two common diffuse optical imaging methods are the time-domain and frequency-domain photon migration measurements. The most challenging task in diffuse optical imaging is the fact that NIR light is highly scattered in human tissues. Unlike X-ray and PET that depend on the ballistic transmitted component of radiation emission, optical imaging has to deal with highly diffused and back-scattered light, and usually involves an inverse problem in image reconstruction, in which tissue optical properties from a given experimental setup and a given set of measurements are determined through an inverse optimization procedure using predicted boundary measurements in forward modeling of photon migration in the given tissue geometry with known optical tissue parameters.

The majority researchers used iterative methods to solve an inverse problem through a convergent algorithm of appropriate objective functions in order to find distributions of tissue optical parameters. For example, O’Leary et al. [1] formulated an inverse problem for fluorescence lifetime tomography using diffuse photon density waves. Ntziachristos et al. [2] presented quantitative optical imaging of human breast in vivo by using NIR diffuse optical tomography (DOT) after the administration of indocyanine green (ICG) for contrast enhancement and examined the optical imaging concurrently with magnetic resonance imaging (MRI). Recently Godavarty et al. [4] developed three-dimensional (3-D) fluorescence tomographic imaging in the frequency-domain using NIR contrast agents in a large tissue-mimic phantom.

Diffusion model based on a diffuse approximation equation was adopted for predicting photon migration in tissues (the so-called forward modeling) in the majority of optical imaging studies. For instance, Wu et al. [5] investigated fluorescence tomographic imaging using early-arriving photons and Laplace transforms of an analytical diffusion approximation of photon migration. However, the diffusion model is only an approximation to the rigorous radiative transfer equation (RTE) and its applicability is limited to scattering-dominant thick media. Biologic tissues are heterogeneous, and contain some regions that are highly absorbing and/or lowly scattering. For more accurate forward modeling, Chang et al. [6] utilized Fourier transform to convert time-resolved RTEs to the frequency domain for describing photon migration in highly scattering media. Klose and Hielscher [7] also considered fluorescence image reconstruction using the stationary RTE as the forward model. Guo and his co-authors solved the transient RTE in time-domain for ultrafast laser transport in biologic tissues using the discrete-ordinates method (DOM) [8, 9] and the radiation element method [10]. An RTE-based forward modeling may consume several hours in computation using a general PC. Even with the use of a simplified diffusion model, a forward modeling may still take dozens of minutes for 3-D heterogeneous tissue geometries.

In the same time, an inverse problem in image reconstruction needs dozens of to hundreds of times of forward modeling. Given the complexity of both the distributions of tissue optical properties and a forward modeling approach, diffuse optical imaging may require several hours to several days to reconstruct an image. Further, an inverse problem is always ill-posed and may lead to divergent results.

In this treatise, a fast and stable optical imaging method is proposed to detect tumor based on relative fluorescence concentration contrast between cancerous tumor and surrounding healthy tissue, rather than the inverse solution of the distribution of tissue optical parameters. Variation in fluorescence concentration distribution due to preferential uptake or washout in diseased tissues leads to detection of diseases. Here we first simulate temporal fluorescence signals using the discrete-ordinates method for a model fluorescence imaging system. The simulated fluorescence data in detectors located on the tissue boundaries are then used for image reconstruction in which a relative fluorescence emission source distribution is established by allocating the detected fluorescence signals to tissue internal voxels. The area with greater fluorescence emission is considered to be cancerous region because the fluorescent agent has a close affinity with cancer. Finally, a reconstructed 3-D image and tomographic images will be presented and compared with the actual tumor shape, size, and location. In actual experimental study in the future, the simulation of fluorescence signals is not needed. Since we only reconstruct the fluorescence emission distribution, the many difficulties associated with the complexities of tissue optical properties, tissue structures, forward modeling, and inverse optimization, are all avoided in the proposed imaging method. Thus, the image reconstruction is expected to be very fast and stable. We will demonstrate our method in the imaging of a small tumor embedded in a model tissue.

 

Fig. 1. Sketch of a model imaging system.

Download Full Size | PDF

2. Simulation and imaging analysis

2.1 Simulation of fluorescence signals

The tissue model for 3-D fluorescence imaging is shown in Fig. 1 where a cubical tissue phantom of side length of 20mm is embedded with a tumor–like object of 4×4×4mm³ in size at the tissue center. After administration of fluorescent agents, exciting laser pulses are incident from the center of one surface of the tissue phantom (e.g., from the plane at X=0) and the consequent fluorescence signals are received by detector arrays located on the four side surfaces of the tissue.

The fluorescence signals are simulated using the discrete-ordinates method for the coupled radiation transfer equations that govern photon transport of both the excitation laser and consequent fluorescence in the model tissue:

where c is the speed of light in the tissue, I^l the light intensity in a discrete direction sl with three direction cosines ξ, η and µ, and σ _e the extinction coefficient that is the sum of absorption coefficient σ_a and scattering coefficient σ_s. The subscripts 1 and 2 in the symbols denote the excitation laser and fluorescence light, respectively. The radiation source terms due to scattering can be formulated in the discrete-ordinates format as

in which Φ ^il represents the scattering phase function, and ω the scattering albedo defined as ω=σ _s /σ _e . A quadrature set of n discrete ordinates with appropriate angular weights wⁱ (i=1, 2, …, n) is used. The choice of quadrature scheme is arbitrary, and n=N(N+2) in a S _N DOM scheme.

$S_c^l$ is the source contribution of the collimated laser irradiation and can be written as

where ξ^c represents the collimated laser incident direction, δ the Dirac delta function, and I ₀ the laser beam intensity irradiated to the surface of the tissue.

$S_F^l$ is the excited fluorescence emission [11] and can be formulated as

where C_v, α_F , and τ are the molarity, quantum yield, and lifetime of fluorescence agent, respectively, and ν is the frequency of radiation.

 

Fig. 2. A representative fluorescence signal.

Download Full Size | PDF

2.2 Image reconstruction

A typical time-resolved fluorescence signal from the present simulation results is shown in the solid curve in Fig. 2. This signal can be divided into two parts. The rising part to the left of the maximum fluorescence intensity I _max depends on early arriving photons of ballistic and snake components, and provides timing information, while the falling part after the I _max is affected by fluorescence lifetime. Of course, the whole curve is also strongly influenced by the medium property through which the fluorescent light passed. Let T_’i,j,k be the time measured at half I _max position in the rising part and T_’mi,j,k be the time measured at the peak intensity point. The indices of ‘i’, ‘j’ and ‘k’ indicate the sequence numbers of a detector in the X, Y and Z directions, respectively. Both T’ _i,j,k and T’_mi,j,k include photon flight time of the excitation laser from incident position to tumor position and the consequent fluorescent photon flight time from the tumor (the fluorescent agent has a close affinity with cancer) to detector (i,j,k). Therefore, the fluorescence flight times from the tumor to detector (i,j,k) at half and full I _max, respectively, are corrected as follows:

in which the subscript 0 represents the detector located at the laser incident position. The times T’₀ and T’_m0 equal to twice of the photon flight time from the laser incident position to the tumor position, respectively measured at half and full I _max.

The majority of fluorescent photons cannot travel straightly to the tissue boundary because of multiple scattering events associated with highly scattering turbid media [5–10]. It is nearly impossible to directly measure the photon flight time along a straight line between the tumor center and a detector except for the case that the tissue is extremely thin which might have no real meaning for noninvasive tumor detection. Nevertheless, we need to estimate the distance between the tumor center and a detector. We hypothesize that the flight time T_i,j,k of early arriving fluorescence photons (the rising part in the fluorescence signal curve) is governed by a Gaussian probability distribution:

where P is the probability, and σ² is the variance of the distribution that can be expressed by

In order to see whether the rising curve in Fig. 2 can be assumed as a Gaussian profile, we have drawn a dotted curve that is symmetric to the rising curve against the I _max line. Clearly these two curves form a distribution that is very close to a Gaussian profile. We have examined this assumption using the numerical data of this study and previous experimental data of laser transport in highly scattering media [13]. The outcome was satisfactory.

Now the fluorescence flight time along a straight line between the tumor center and a detector is estimated as

The average distance between the tumor and a detector located at (i, j, k) is then calculated by

The equation for a spherical surface centered at detector (i, j, k) with a radius of R_i,j,k is

in which (D_i , D_j , D_k ) are the coordinates of the detector location. We define a voxel space line as X=V_i ’ and Y=V_j ’, and solve for the intersection points between the spherical surface and the space line.

When ${R_{i\mathit{,}j\mathit{,}k}}^2$ -(V_i′ ’-D_i )²-(V_j′ ’-D_j )²≥0, the intersection points are

Otherwise, there is no intersection point. (V_i’,’,V_j’ ’,V_k ’) is the voxel and/or pixel grid position for image reconstruction.

If the intersections are located in the tissue geometry, the value of I _max detected by this detector is then allocated to the voxel (V_i’ ’,V_j’ ’,V_k ’) and its adjacent eight voxels. Repeating the same procedure for every detector, we can get a relative intensity of fluorescence emission from each voxel contributed from all detectors:

Here, E(V_i ’,V_j ’,V_k ’) is the relative fluorescence intensity of voxel (V_i ’,V_j ’,V_k ’), I_D (D_i , D_j , D_k ; V_i ’,V_j ’,V_k ’) is the intensity contributed from detector (D_i , D_j , D_k ) to voxel (V_i ’,V_j ’,V_k ’), and I _max(D_i , D_j , D_k ) is the maximum fluorescence intensity at detector (i, j, k). The computational grid is N×M×K.

3. Results

In the simulation model, we considered pulsed excitation light of pulse width 1 ps and wavelength at 780 nm. ICG dye was used as the fluorescent agent. Fluorescence signals were received for wavelength greater than 800 nm. The absorption coefficient and reduced scattering coefficient of the model tissue were 0.0023mm^-1 and 1.018mm^-1 against the excitation laser, and 0.0031mm^-1 and 0.864mm^-1 against the fluorescence, respectively. The absorption coefficient of the ICG solution was 0.03mm^-1 against the excitation laser (assuming C_v0=1µM), and 0.0025mm^-1 against the fluorescence. Those property values were used by Godavarty et al. [4]. We considered an uptake distribution of the fluorescent dye in the model tissue. The dimensionless molarity C_v/C_v0 of the ICG solution was between 0.9 and 1.0 in the tumor region, while it was randomly distributed in the surrounding healthy tissue with values between 0.01 and 0.06. The uptake ratio varied between 15:1 and 100:1. The coupled radiation transfer equations were solved with the S ₁₀ DOM method [9]. A computational grid of 25×25×25 was used in the simulation.

 

Fig. 3. Reconstructed 3-D tumor image.

Download Full Size | PDF

Figure 3 exhibits the reconstructed 3-D tumor image for the imaging model shown in Fig. 1. A filter technique based on the gray-level histogram concept [12] was used to enhance the image contrast and quality. The blue block in Fig. 3 has greater fluorescence emission value and stands for the embedded cubical tumor. In order to see clearly the location of the reconstructed tumor image, the 3-D image is projected onto three planes normal to the X, Y, and Z directions, respectively. The three projections are shown in Figs. 4a, 4b, and 4c. The squares enclosed with black dash-dot lines show the actual tumor location, size, and shape. The red areas enclosed with bright green lines are the location, size, and shape of the detected tumor. It is seen that the images closely match the actual embedded tumor. Slight displacements of tumor image locations are observed. This is caused by the estimation of the fluorescence flight time along a straight path in order to get the distance information. It is very difficult to precisely measure this timing information because of the highly scattering phenomenon of NIR light in biologic tissues.

Tomographic images at any plane of any orientation are easily obtainable by slicing the 3-D image. Figures 5a, 5b, and 5c display the tomographic images along the X-, Y-, and Z-axis, respectively. The eight slices in each figure are selected in the respective axis locations from 7.5mm to 12.5mm, where the detected tumor shape, size, and location at different axis positions in the model tissue are illustrated clearly. The tomographic images sliced around the tumor center capture the actual tumor very well. However, deformations are found for tomographic images sliced at the edge of the tumor. Because the cubical tumor has an abrupt or discontinuous change at its boundaries, it is a formidable task to reconstruct a sharp image for any imaging method. Also such deformations are partially attributable to the varied fluorescence solution uptake. We have assumed a steady change of fluorescent agent uptake from 0.90 at the edge of the tumor to 1.0 at the tumor center and a random distribution of agent uptake in the background. Considering the facts that the tumor image was retrieved from a certain depth inside a highly scattering dense tissue and the embedded model tumor has a very sharp shape, however, the deformations and slight displacements are reasonably acceptable. As a matter of fact, the present images have excellent quality and high contrast as compared with most published optical imaging results in the literature.

In order to validate our method, we further imaged a small tumor embedded at an off-center position of the tissue. Figure 6 displays the reconstructed 3-D tumor image. The tumor was centered at the position of X=15 mm, Y=15 mm, and Z=15 mm. The size of the tumor was reduced to 2.4×2.4×2.4mm³. Other conditions remain the same like the preceding model. The detected tumor in Fig. 6 closely matches the preassigned tumor location and size conditions.

A unique feature of the present method is that it is extremely computationally efficient in the image reconstruction. When a DELL PC with one CPU of 1.7 GHz was used, it took less than five minutes to reconstruct the 3-D image for the demonstrated model imaging problems. This makes the method very attractive to be deployed for clinical application for early detection of cancer. One limitation is that the use of the fluorescence flight time in the calculation of the distance between a tumor and a detector may restrict the estimation of tumor center location to presence of a single tumor. We are developing new means for estimating centers of multiple tumors.

It is worth noticing that the fluorescence signals used in the present article were obtained through numerical simulation based on forward modeling of radiation transfer. In actual experimental studies, the forward modeling is not required. We only need to measure the time-resolved fluorescence signals in the detector arrays surrounding the four side surfaces of the tissue. In this modeling study, an array of 8×8 detectors was placed on each of the side surfaces. In actual measurements, we may reduce the number of detectors without sacrificing the accuracy. For instance, an array of 4×8 detectors may be deployed (each side of the tissue has 8 detectors and all the detectors are connected with a 32-channel streak camera). The array is then moved vertically (e.g., along the Z-direction in Fig. 1) to get fluorescence signals on the whole side surfaces. When an ultrafast laser with repetition rate of 76 MHz is used for excitation, the time gap between two successive short pulses is only 13 ns. The whole measurement can be completed within a very short time period.

 

Fig. 4. Projections of the 3-D tumor image: (a) on the Y-Z plane, (b) on the X-Z plane, (c) on the X-Y plane.

Download Full Size | PDF

 

Fig 5. Tomographic images: (a) along the X-direction, (b) along the Y-direction, (c) along the Z-direction.

Download Full Size | PDF

 

Fig. 6. Reconstructed 3-D image for a small tumor embedded at off-center of the tissue.

Download Full Size | PDF

4. Conclusions

We developed a novel 3-D optical fluorescence imaging method for early detection of small tumor. The image reconstruction in this method does not involve any inverse optimization. It is very efficient and can be completed in minutes by a PC. We demonstrated the feasibility of this method in the imaging of a small tumor embedded in a cubical model tissue. The fluorescence signals were obtained via simulation of the coupled transient RTEs that govern the photon migration of both the excitation laser and consequent fluorescence in the tissue. The simulation was set to mimic actual experiment. The uptake of fluorescent agent was considered to be non-uniformly distributed in the tissue and tumor. The reconstructed 3-D image closely matches the actual tumor in terms of shape, size, and location. Numerous tomographic images are obtainable by slicing the 3-D image. The detected images are of high contrast and high accuracy.