1. Introduction

Fluorescence/ bioluminescence diffuse optical tomography (F/BDOT) has been studied for some medical applications [1,2]. Near-infrared light with the wavelength range from about 700 nm to 900 nm propagates diffusely in biological tissues which strongly scatter and weakly absorb near-infrared light. By measuring the light intensities at the surface of the medium and by solving an inverse problem, we can reconstruct the tomographic images of the parameters such as the absorption and scattering coefficient inside the measured objects [3–5]. Based on this principle, F/BDOT primarily reconstructs the tomographic images of the light source strengths in the medium where the light sources are fluorescence/bioluminescenc molecular probes.

One of the promising applications of F/BDOT is molecular imaging of small animals. In drug developments, the process of the drug delivery in small animals must be monitored by some methods. By combining F/BDOT technology with fluorescence/bioluminescence-labele drugs, in-vivo molecular imaging of the drug delivery process can be realized. F/BDOT also allows us to study disease and to evaluate the efficien y of a treatment over time for an identical animal [1]. Another potential application of FDOT is the screening of breast cancer. Since the optical properties of lesions are different from healthy tissues due to angiogenesis etc., the breast cancer can be diagnosed by reconstructing the scattering and absorption coeffi cients [6–11]. Now it is expected that some fluorescen contrast agents targeting tumors can enhance the image contrast of the cancer without relying on the endogenous components such as hemoglobin, water or fat [12–14]. FDOT can provide the information of the fluorescen properties, e.g. quantum yield and lifetime, which are affected by pH and oxygenation of diseased tissues, and these information must be useful for diagnoses. Schemes of F/BDOT including the acquisition and analysis of measurement data are still at the stage of development. The reconstruction techniques of the fluorescen properties are actively discussed [15].

The methods of image reconstruction by solving the inverse problem for F/BDOT are roughly categorized into two major approaches. One approach recovers the fluorescenc properties in the equations describing light propagation by solving a nonlinear optimization problem based on gradient methods, e.g. Newton-Raphson method [16–21]. The photon diffusion equation (PDE) is often used to describe the diffuse propagation of light in tissues. The other inverts the integral equation relating the fluorescenc source parameters to the measurement data with the Green’s functions for PDE [22–28]. The former approach can reconstruct the image of the distribution of the fluorescenc properties, such as the fluorophor concentration, the quantum yield, the lifetime, and various parameters in PDE such as the absorption and scattering coefficient of the medium. The residual errors between the measurement data and the forward solution are minimized in the least squares sense. In every iteration of the optimization process, the gradient to update the parameters is calculated. This approach requires high computational costs generally, because the finit element method (FEM) calculation is required iteratively in many realistic cases where the analytical solution of PDE cannot be obtained. And the reconstructed images depend on some parameters such as initial values for iterative updating process.

The latter approach solves the linear vector-matrix system equations describing the relation between the measurement data and the source term. The latter approach using the conventional methods, such as pseudo-inversion, the algebraic reconstruction technique (ART) and the least-squares QR method (LSQR), is computationally less intensive than the former approach [20–22, 25]. However, the latter approach involves some issues on the estimation of the depth and area of the source. In the highly underdetermined inverse problems, such as F/BDOT reconstruction, it is often requested to reconstruct tomographic images from a small number of the measurement data. In this case, it is difficul for the conventional methods to recover the sources deeply located in the medium as shown by the simulation in this study. A new method for the reconstruction from a small number of the measurement data is desired for extensive practical implementation of F/BDOT. It is preferable that the solving method is fuss-free and not so much dependent on prior information, regularization parameters or initial guess. For effective diagnoses of breast cancer in early stages, other methods characterized by high accuracy and spatial resolution will be required.

In this paper, we propose a new method for reconstructing fluorescence/bioluminescenc source distributions by use of a new spatial filte and successive updating process of the forward model. A conventional spatial filte, which was originally employed for radar and sonar technologies, distinguishes and extracts desired signals from interferences by exploiting the differences at locations where they are generated [29]. Recently, the existing spatial filte is applied to the inverse problem in functional brain imaging, such as electroencephalography (EEG) and magnetoencephalography (MEG) [30–33]. By using the forward model, the new spatial filte proposed in this study transforms the measurement data to an estimation of a source, fluorescent bioluminescence light intensity at a location of interest. The new spatial filte is formulated to indicate the correct depth of the emitting light source by use of a small number of the measurement data even in the underdetermined inverse problem.

The estimated source is used to modify the forward model. From the result of the source estimation by the spatial filte, the proposed algorithm determines which sources are dispensable in the forward model and which are not. The dispensable sources to produce the measurement data are removed gradually from the forward model. The iterative estimation process which consists of spatial filterin and updating the forward model improves the localization of the estimated source distribution [34]. It has been reported that a similar iterative refinemen of the model achieves high-quality reconstruction in 3D [35].

Additionally, the reduction of the artifacts induced by noises accompanying the spatial filte is introduced. The effects of the random noises added to the measurement data and the mismatch in the optical properties between the forward model and the real object are reduced by the method based on the singular value decomposition (SVD). The estimation by the spatial filte is used to determine the noise space. The noise space is removed from the measurement data, and as the result, the reconstruction is improved and becomes robust.

In the following sections, we introduce the new reconstruction scheme mentioned above. Some numerical experiments are performed to show how the algorithm works and how the proposed method is different from the conventional methods. Some numerical case studies intend to provide the readers with the quantitative evaluation of the method and quantitative guidelines for practical uses.

2. Methods

2.1. Forward model of light propagation

Light propagation in biological media is essentially the phenomenon of the transport of radiant energy. The conservation of the radiant energy is rigorously described by the radiative transfer equation (RTE). This complex integro-partial differential equation can be numerically solved by the discrete ordinate method deterministically or by a Monte Carlo simulation stochastically. The solutions will be accurate but computationally costful.

PDE is a useful and widely used approximation of RTE. By taking the firs order PN approximation with the spherical harmonic expansion of the quantities in RTE, PDE is obtained.

Let us assume that the emitting fluorescence/bioluminescenc sources are reconstructed from continuous wave (CW) measurement data which are the fluorescenc light intensities measured at the surface of the medium. Then the propagation of the fluorescenc light in bulk tissues can be described by the steady-state photon diffusion equation (1) [3],

where r is the position vector, D = 1/(3µ′_s) is the diffusion coefficien with the reduced scattering coefficien µ′_s, µ_a is the absorption coefficient Φ is the fluenc rate which is the integral of the photon radiance over all directions, and q ₀ is the fluorescenc light source. Equation (1) can be solved under the boundary condition expressed by Eq. (2),

where A = (1 + R_f)/(1 − R_f), R_f is the internal reflectio ratio, and n is the outward normal to the surface ∂Ω of the medium. The light intensities measured at the surface, Γ(r), are the flu es of the fluenc rates at the surface given by the left hand side of Eq. (2).

The source term of the fluorescenc light, q ₀(r), is described as q ₀ = ηµ_afΦ_x, where η is the quantum yield of the fluorophore µ_af is the absorption coefficien of the fluorophor and Φ_x is the fluenc rate of the excitation light. The mechanisms of the emissions are different between bioluminescence and fluorescence However, propagation of light emitted by either bioluminescence or fluorescenc probes is described by Eq. (1) commonly. It is possible to reconstruct the fluorescenc properties such as η and µ_af after the reconstruction of q ₀ when Φ_x is given. Therefore, this study is focused on the reconstruction of q ₀ which can be converted to either fluorescenc or bioluminescence light sources.

In order to solve Eq. (1) by FEM, Eq. (1) is usually discretized as Eq. (3),

where Φ and Q are the vectors of the discretized fluenc rates and the source strengths at N nodes of FEM, and the N × N matrices K and C are calculated with the basis functions and the optical properties [36]. The matrix B is for the boundary condition. The entries of the matrices are, K_ij = ∫_Ω D∇ψ_i·∇ψ_jdr, C_ij = ∫_Ω µ_aψ_iψ_jdr and B_ij = 1/(2A) ∫_∂Ω ψ_iψ_jdγ, where i and j denote the indexes assigned to the nodes, B_ij is the integral over the boundary surface, and ψ_i is a basis function.

If the source term Q is given, the vector-matrix equation (3) can be solved for Φ,

where G = {K(D) + C(µ_a) + B}⁻¹.

The measurement data are the flu es of the fluenc rates at the surface, −D∇_nΦ, and they can be readily calculated as 1/(2A)Φ due to the boundary condition, Eq. (2). Since the flu es can be measured only at a limited number of positions, M, at the surface of the medium, the forward model in F/BDOT is formulated as follows,

where m is a vector composed of the measurable flu es at M positions, and the M × N matrix L consists of the M row vectors corresponding to the measurement positions among the rows of the matrix G.

2.2. Inverse problem to reconstruct the light source distribution

2.2.1. Spatial filte

The inverse problem for the image reconstruction in this study is to determine the source term Q. The components of Q are linearly mixed into m by L. To reconstruct the source term, we design a spatial filte which extracts a single component in the source vector from the measurement data. The spatial filte exploits the differences in the contributions of spatially different sources to the measurement data. The spatial filte is constructed for every node with the source strength to reconstruct the whole distribution of the source strengths. The spatial filte w_k we employ in this study transforms the set of the measurements m to the estimation of a source strength q̂_k at a specifi location. The estimated source distribution Q̂ = [q̂ ₁, q̂ ₂, ⋯, q̂_N]^T is obtained with the M × N matrix W which consists of the spatial filter w_k,

where q̂_k is the k-th source estimated by w_k, which is the k-th column vector of W, represents the spatial filte targeting on the k-th source strength q_k, and l_i is the i-th column vector of L. l_i represents the contribution of the single source with a unit strength in the i-th position to the measurement data. The spatial filter are constructed for all nodes of FEM model where the sources can exist.

The estimation of the source strength in the k-th position agrees with the true one when the spatial filte and the column vectors l_j hold w^T_kl_j, where δ_kj is the Kronecker’s delta. However, in the ill-posed inverse problem such as the source reconstruction for F/BDOT in this paper, it is not possible to obtain such a spatial filte for every l_j since the matrix L is singular. For successful estimation, w^T_k l_j (k ≠ j) appearing in Eq. (7) should be close to zero, and w^T_k l_k should be as close to unity as possible. In other words, the spatial filte is the solution of the following maximization problem,

We can obtain the spatial filte as the solution of Eq. (8) with the Lagrange method of multipliers,

When (LL^T) is singular, its inverse matrix in Eq. (9) can be substituted by a pseudo-inverse matrix. The coefficien in Eq. (9), {l^T (LL^T)⁻¹ l_i}^−1/2, differentiates the spatial filter from the Moore-Penrose pseudo inversion.

When we solve the minimization problem, ${min}_{w_k}{\parallel w_k^{T}L\parallel }^2$ subject to (w^T_kl_k) = 1, w_k = (LL^T)⁻¹ l_k is obtained as a spatial filte. In this case, the reconstructed distribution of the sources coincides with the minimum norm solution by the Moore-Penrose matrix inverse, which is difficul to reconstruct the emitting sources deeply embedded inside the medium in the underdetermined inverse problem. On the other hand, the proposed spatial filte has an advantage that the maximum of the estimated source strengths is located at the position of the source expected to make the largest contribution to the measurement data, regardless of the distance from the detectors. This comes from Eq. (8), which means that the spatial filte is the vector which maximizes q̂_k estimated from the measurement data generated by the single source with the unit strength at the targeted k-th position.

It is obvious that the proposed spatial filter satisfy (w^T_kl_k)² ≥ (w^T_jl_k)², since the spatial filte w_k maximizes the squared product of (w^Tl_k)², where w is the arbitrary vector satisfying ∥w^TL∥² = 1. If the source q_k contributes to the measurement data to a great extent, i.e. m≃l_kq_k, the estimated source strength q̂_k = w^T_k m at the k-th position takes a larger value than those at the other positions. Thus we can expect that the spatial filter identify the correct locations of strong sources, even if the spatial filte can not obtain w^T_kl_k sufficientl close to unity and as the result the estimated source strength cannot agree with the true one.

In practice, the measurement data generated by multiple sources which are closely located each other are so similar that it is not possible to completely identify their true positions of the emitting sources. However, the spatial filter surely indicate the positions with the dominantly contributing sources.

Conventional methods inverting Eq. (5) pay attention to the characteristics of Q such as the smoothness, the sparseness, etc. The regularization technique is employed to obtain the solution with the characteristics expected for the true solution. On the other hand, the method proposed in this paper pays attention to the relation among the components of Q. The spatial filte investigates which source is large and which is not.

The most important advantage of the proposed method is that the spatial filte allows us the estimation of the depth of the emitting fluorescen source. Generally, even though the true emitting source is located at a deep position in the object, the largest value of the source strength is often reconstructed close to the detector, which is the minimum-norm solution in underdetermined problems. On the other hand, the proposed spatial filte does not obey the minimum-norm criterion but tends to output the larger value of the source strength at the position where the true source is expected, regardless of the distances from the detectors.

The method proposed in this paper does not depend on whether the reconstruction is 2D or 3D. The system matrix L depends on the shape of the medium, the optical properties and the positions of the detectors. The spatial filte depends on L.

2.2.2. Updating of the forward model

We propose a method updating the forward model which utilizes the results of spatial filterin mentioned above. In the animal experiments for instance, we expect that the true light sources emitted by probes labeled with fluorophore or luciferase enzymes, are usually concentrated within small regions, although the reconstructed source distributions have generally wider regions than the true regions. By reducing the number of the possible source positions in the object, the ill-posed nature of the F/BDOT inverse problem will be alleviated. As the result, the reconstructed image of the source distribution will be improved.

In this updating process, the result of the source estimation via spatial filterin is fed back to the forward model. The spatial filter output large values at the source positions where the source strengths make great contributions to the measurement data. Therefore the positions where the reconstructed source strengths are small can be judged dispensable. In this context, once we get the information of the source strengths by the spatial filte, the “weight” can be define to represent the importance of each source position as follows,

where ∥ · ∥_∞ denotes the maximum norm. Then the forward model is updated by use of the weight, p_i, as the following,

This process transforms our inverse problem formulated by Eq. (5) to that reformulated as m = LPQ, where P is the diagonal matrix with p_i. The weight, p_i, takes the value in the range of [0, 1]. p_i = 1 means that the i-th source position is necessary for the forward model, while p_i = 0 means that the source position can be ignored. The spatial filte is recalculated with the updated matrix L ← LP.

By ignoring the sources at the dispensable source positions and enhancing the sources at the indispensable positions in the forward model with using the above manner for the next iteration, the localization of the source distribution is improved. The process is iterated until the source distribution is refine satisfactorily.

Finally, the source distribution is adjusted to fi the measurement data in the least-squares sense. While the proposed spatial filte promotes to reconstruct the true source strength, it is not guaranteed that the set of the spatial filter provide the solution Q̂ of m = LQ̂. However, the solution can relatively compares the source strengths among the source positions. Now we introduce a coefficien for the adjustment, α, by solving the optimization problem as follows,

Finally, the source distribution, Q̄, is estimated by the following equation,

2.2.3. Reduction of artifacts based on SVD

We apply SVD technique to reduce artifacts which is fictitiou light sources due to the measurement noises and the calculation errors. The measurement noises mentioned in this paper refer to the statistical random additive noises caused by several factors such as the electronics of the measurement system. The calculation errors include the errors induced by the finit numbers of elements in the FEM calculation and the errors caused by the mismatch between the forward model and the actually measured object, particularly the incorrect prior information about the optical properties of the medium. In the cases where the measurement noises and the calculation errors are not negligible, the measurement data are described as follows,

where ε represents the sum of the noises and the errors.

The spatial filter project the measurement data onto the source strengths without distinguishing signals and noises. Here, “signal” means the light emitted by the probes in the medium such as the fluorophores Even if the noises are small, updating the forward model does not eliminate the artifacts but localizes them. Fortunately, the feature of the spatial filte mentioned before can be of assistance to the reduction of the effects of the noises.

By using the source strength estimated by the spatial filte, we calculate the estimation of the contribution of the i-th source as l_iq̂_i, which is used in the updating process. The basis vectors which compose the estimated contributions are calculated by SVD of the updated L (or the eigen value decomposition of LL^T). The basis vectors corresponding to the large singular values (or the eigen values) are the dominant principle components of the estimated contributions. It means that the basis vectors are expected to contribute to the measurement data significantl. By assuming that the signal-to-noise ratio (SNR) of the data is high to some extent, the basis vectors with large singular values compose the signal space. And the uncorrelated basis vectors with small singular values can be regarded as the noise space.

In the process to reduces the artifacts, we calculate the singular value decomposition of the updated matrix L which consists of the products, l_iq̄_i, first L is decomposed as L = UDV^T where U and V are the M × M and N × N matrices respectively, whose columns form the orthonormal basis, and D is the M × N matrix which has the singular values s ₁, s ₂, ⋯, s_j, ⋯, s_M and (N − M) components with zero in the diagonal. The M × 1 basis vectors, u_j (j = 1, 2, ⋯ M) form the matrix U = [u ₁,u ₂,⋯,u_M]. The basis vector and the singular value with the identical index j correspond with each other. Note that the measurement data are represented by a linear combination of the basis vectors u_j, i.e. m = U_y, where y is a M × 1 vector determined uniquely by U. Therefore, by removing the basis vectors corresponding to the small singular values from the measurement data, i.e. by setting the corresponding components of y as zeros as the following procedure, ε in Eq. (14) can be reduced,

where H is the M × M diagonal matrix with the j-th diagonal components H_jj having 1 or 0 with depending on whether the j-th basis vector is the signal or the noise space. The number of the basis vectors selected as the signals depends on the situations. The number of the basis vectors for the signals increase as the SNR is high.

3. Results and discussions of numerical experiments

3.1. Conditions of numerical experiments

The performances of the reconstruction with spatial filterin and updating the forward model are demonstrated in this section. A 2D circular medium with a radius of 40 mm is the experimental object. Basically, the reduced-scattering and absorption coefficient of the background are given as 0.8 and 0.007 mm⁻¹, respectively, by assuming the application to breast cancer screening using fluorescen dye working in the near-infrared wavelength range [9]. CW excitation light illuminates the object at arbitrary locations, and the detectors measure the fluores cence emission light by cutting the excitation light with an optical filte. The parameter A in Eq. (2) is set as unity assuming no internal reflection Some regions with unit source strength are included in the medium.

16 detectors are placed circularly with an equal spacing at the surface of the medium, and the source distribution is reconstructed from 16 measurement data. The nodes at the surface of the medium are excluded from the source positions, and the source strengths of the nodes at the surface are set as zero. For generation of simulated measurement data, the numbers of the nodes and the triangular elements of FEM are 6561 and 12800, respectively.

Other conditions are changed from one case to another for the particular purpose of the numerical experiments in the following subsections.

3.2. Improvement in localization by updating the forward model

The effect of updating the forward model is demonstrated in this subsection. In this simulation, the FEM mesh for inverse calculation is identical with that for generating the measurement data. The generated measurement data contain no noise and the forward model has no mismatch in the optical properties in this case.

Figures 1 (a) to (g) show the results of the firs to seventh iteration, and Fig. 1 (h) shows the true image for comparison. The true emitting source region has the radius of 2.5 mm at the position of (x,y) =(20 mm, 0 mm) with the source strength of unity. From Fig. 1, it is understood that updating the forward model improves the localization of the reconstructed emitting source regions.

 

Fig. 1. The reconstructed source distributions: (a)–(g) are in order of the number of the iteration. (h) is the true source distribution; i.e., true source at (x, y) =(20 mm, 0 mm) with the radius of 2.5 mm and the source strength of unity.

Download Full Size | PDF

 

Fig. 2. The residual error ratio as the function of the updating iterations.

Download Full Size | PDF

Currently we do not have any criterion to terminate the iteration, but the sources can be quickly localized; only a few iterations are needed. Figure 2 shows that the estimation error ratio, define as ∥m−LQ̄∥²/∥m∥² ×100 (%), decreases with the increase of the number of iterations in successful cases. This feature can help users choose a reconstructed image which is recognized as an appropriate one from a physiological point of view.

3.3. Comparison with the conventional methods

The result obtained by the proposed method is compared with those by the typical conventional methods for the conditions with Fig. 1(h). Neither measurement noises nor calculation errors are considered.

After six iterations of the proposed reconstruction, the estimation error ratio is sufficientl small, and the reconstructed source region is well localized at the correct position with the almost correct source strengths as shown in Fig. 3(a).

 

Fig. 3. The reconstructed source distributions by (a) the proposed method, (b) the minimum-norm solution, (c) and (d) LSQR, and (e) and (f) Algebraic Reconstruction Technique. The true distribution is shown in Fig. 1 (h).

Download Full Size | PDF

As a typical conventional reconstruction method, we firs choose the least squares method which gives the least squares minimum norm solution [37]. The solution is obtained as Q = L^T (LL^T)⁻¹ m. To minimize the norm of the solution, the locations of the reconstructed emitting sources are shifted to the surface of the medium as shown in Fig. 3(b). It is obvious that the solution with weakly emitting sources provide small norm.

The results obtained by the solver such as LSQR which is based on the conjugate gradient method minimizing ∥m−LQ∥² [38] are presented in Figs. 3 (c) and (d). The solution depends on the initial guess. When we choose the initial guess having zero source strengths uniformly, the result shown in Fig. 3(c) becomes identical with the minimum norm solution. With the initial guess having the strengths of 0.005 uniformly, the reconstructed source strengths close to the detectors appear as spots to minimize the residual errors.

Another popular method is ART [39, 40]. The source distribution reconstructed by ART is shown in Fig. 3 (e) when the initial guess is given as zero uniformly. This solution is identical to the minimum-norm solution shown in Fig. 3 (b). Figure 3 (f) shows the image reconstructed by ART with the initial guess given as 0.01 uniformly. The region with non-zero source strength is broader than the true one, and the reconstructed source strengths are significantl smaller than the true value. It is difficul to obtain a well-localized source distribution by use of ART since the reconstruction is affected by the initial guess and the relaxation parameter [20–22].

Many reconstruction schemes based on the conventional methods have been proposed and have succeeded in reconstruction. However, these schemes involve various know-hows, such as the number of the measurement data. While the conventional methods cannot reconstruct the emitting sources deeply inside the medium, the proposed method depends less on the know-hows and localizes the sources in the highly underdetermined inverse problem.

3.4. Reduction of the artifacts induced by various noises and errors

This subsection demonstrates that the proposed method with the use of SVD technique can reduce the artifacts induced by the measurement noises and calculation errors, particularly (1) the random measurement noises, (2) the background emissions, and (3) the mismatch in the optical properties between the model and the experiment. Considering the inverse crime problem [41], the FEM meshing for solving the inverse problem is not identical to that for generating the simulated measurement data in the following simulations. The numbers of the node and the triangular elements are 1681 and 3200 respectively for the inverse problem, while those for generating the measurement data are 6561 and 12800, respectively.

3.4.1. Reduction of the effects of the random noises in the measurement data

The measurement data fluctuat due to many factors such as the electronics of the experimental instruments and the detector placement. In such cases, the measurement noises are usually modeled as the Gaussian noise. The conditions are the same as those is shown in Fig. 1(h) except that the Gaussian noises with the average of zero are added to the simulated measurement data. The number of the updating iterations is six, and the reduction of artifacts is carried out. The averages and the standard deviation (SD) of the reconstructed source strengths are estimated by a hundred trials of the reconstructions.

Figures 4 (a-1) and (a-2) compare the results with and without removing the noise space. Gaussian noises have SD of 1 percent of the maximum of the measurement data (SD_noise = 0.01). In Fig. 4 (a-1), the averages of the reconstructed source strengths clearly show the single emitting source region and most of SDs are less than 0.02, which are sufficientl small compared with the average values. In this case, after SVD of the updated matrix L, we select the larger singular values that yield more than 98 percent (λ_SVD = 0.98) of the total, and the corresponding basis vectors u_i are chosen as the signal space. The noise space is removed by Eq.(15) in every iteration. By using the proposed method of removing the noise space the emitting source region can be well identifie in most of the reconstructions with the noise of SD_noise=0.01. The artifacts tend to appear near the detectors far from the emitting sources. The light intensities measured by the detectors far from the emitting sources are weak with low SNR. Without removing the noise space, it is difficul to identify the emitting sources in most cases. The strong artifacts are reconstructed widely and randomly.

 

Fig. 4. The averages (top) and the standard deviations (bottom) of the reconstructed source strengths: (a-1) SD_noise = 0.01 and λ_SVD = 0.98, (a-2) SD_noise = 0.01 and λ_SVD = 1.0, (b) SD_noise = 0.001 and λ_SVD = 0.995, and (c) SD_noise = 0.1 and λ_SVD = 0.96.

Download Full Size | PDF

Figure 4 (b) shows that when the noise is sufficientl small, SD_noise=0.001, the proposed reconstruction scheme is quite robust, and the emitting sources can always be found at the correct positions. When the noise is significant SD_noise = 0.1, as shown in Fig. 4 (c), the reconstructed images can not provide the information of the fluorophore although removing the noise space is carried out.

3.4.2. Reduction of the artifacts induced by the background emissions

We evaluate the effectiveness of the proposed method for reduction of the artifacts induced by the background emissions due to thinly and widely distributed fluorophore and autoflu orescence. We assume that the source strengths of the background emissions are uniformly distributed in the model geometry of Fig. 1(h) with the same conditions unless otherwise specified The strength of the background emission, q_BG, at each node is assigned as (a) 0.0001, (b) 0.001 or (c) 0.01. Figure 5 plots the measured intensities from the target and background emissions. The number of the updating iterations is six.

 

Fig. 5. The measured intensities of (solid line) the target signal light and of (dashed lines) the background emissions.

Download Full Size | PDF

 

Fig. 6. The reconstructed source distributions using the proposed noise reduction with (a-1) q_BG = 0.0001 and λ_SVD = 0.995, (a-2) q_BG = 0.0001 and λ_SVD = 0.97, (b-1) q_BG = 0.001 and λ_SVD = 0.97, (b-2) q_BG = 0.001 and λ_SVD = 0.80, and (c-1) q_BG = 0.01 and λ_SVD = 0.97, (c-2) q_BG = 0.01 and λ_SVD = 0.8.

Download Full Size | PDF

In the reconstructed image shown in Fig. 6(a-1) for the case (a) with λ_SVD = 0.995, the noise reduction is found to be insufficient The uniformly distributed background emissions are measured equally by all the detectors and appear as the artifact emitting source localized at the center of the circular medium. When the background emission is not uniform, the reconstructed artifact sources will be distributed corresponding to the nonuniformity. The image shown in Fig. 6 (a-2) for the case with λ_SVD = 0.97 demonstrates that the influenc of the background is completely removed.

The reconstructed image for the case (b) with λ_SVD = 0.97 is shown in Fig. 6 (b-1) and is accompanied by the artifact. When λ_SVD = 0.8, the single emitting source region is well localized as shown in Fig. 6 (b-2). It is difficul to reconstruct the single target source region at the correct position for the case (c) as shown in Figs. 6 (c-1) and (c-2). Particularly, in Fig. 6(c-2) with λ_SVD = 0.8 the target source region disappears due to the low SNR.

3.4.3. Influenc of the mismatch in the optical properties

The mismatch in the optical properties between the forward model and the real object induces the artifacts. We evaluate the effect of the proposed method on the mismatches in µ′_s and µ_a.

The graph of Fig. 7 plots the simulated light intensities measured by the 16 detectors at different angles in the homogeneous and inhomogeneous media. The homogeneous medium has µ′_s = 0.8 mm⁻¹. And the inhomogeneous one has a particular sector regions in the range of (3/4)π ≤ θ ≤ (5/4)π with different µ′_s, (a) 0.72, (b) 0.4 or (c) 0.08 mm⁻¹. The measured light intensities in the particular region increase with the decrease in µ′_s.

The measurement data are generated by the inhomogeneous distribution, and the reconstructions are attempted with homogeneous µ′_s=0.8 mm⁻¹. The number of the updating iterations is six. The true source distribution is shown in Fig. 1 (h).

Figures 8(a) and (b) show that the reconstructions with the appropriate λ_SVD are successful in the cases (a) and (b), respectively. The artifact reduction alleviates the mismatches. On the other hand, it is not possible to reconstruct the images without the artifacts as shown in Fig. 8 (c-1) and (c-2) in the case (c). µ′_s of the particular region is so small that the measured light intensities increase to around 10% of the maximum of the measurement data as shown in Fig. 7. The artifacts appear to compensate the mismatch in the optical properties. The artifact is readily removed when the correct distribution of µ′_s is used for the reconstruction as shown in Fig. 8 (c-3). In the case of µ′_s = 8.0 mm⁻¹ in the particular region which is 10 times larger than µ′_s in the rest of the medium, the emitting source region is reconstructed well (Fig. 8(d)). The light intensities measured by the detectors far from the emitting source region are essentially so small that the strong attenuation by the particular region hardly affects the measurement data.

 

Fig. 7. The measured light intensities in the inhomogeneous medium (dashed curves) and in the homogeneous medium with µ′_s = 0.8 mm⁻¹ (solid curve). µ′_s in the region of (3/4)π ≤ θ ≤ (5/4)π of the inhomogeneous medium is 0.72 mm⁻¹, 0.4 mm⁻¹ or 0.08 mm⁻¹.

Download Full Size | PDF

 

Fig. 8. The source distributions reconstructed by use of incorrect homogeneous distribution of µ′_s. λ_SDV used in the noise reduction and the µ′_s in the particular region in the inhomogeneous medium are (a) λ_SVD = 0.995, µ′_s =0.72 mm⁻¹, (b) λ_SVD = 0.98, µ′_s = 0.4 mm⁻¹, (c-1) λ_SVD = 0.95, µ′_s =0.08 mm⁻¹, and (c-2) λ_SVD = 0.99, µ′_s =0.08 mm⁻¹. For (c-3), the correct distribution of µ′_s and λ_SVD = 0.99 are used. µ′_s =0.08 mm⁻¹ in the particular region. And for (d), homogeneous µ′_s and λ_SVD = 0.995 are used. The particular region has µ′_s =8.0 mm⁻¹.

Download Full Size | PDF

Then the mismatch in µ_a is studied by the same manner as that in µ′_s. Figure 9 shows the measured intensities for the homogeneous medium with µ_a = 0.007 mm⁻¹ and the inhomogeneous media with µ_a of (a) 0.0035 or (b) 0.0007 mm⁻¹ in the particular region. µ′_s is given as 0.8 mm⁻¹ uniformly. The smaller µ_a is, the larger intensities measured in the particular region are.

Figures 10 (a) and (b-1) show the source distributions reconstructed with incorrect homogeneous µ_a = 0.007 mm⁻¹ in the cases (a) and (b), respectively. The appropriate λ_SVD of 0.99 and 0.97 recovers the emitting source regions. Insufficien removing of the noise space with λ_SVD = 0.99, leads to the artifacts as shown in Fig. 10 (b-2). Figure 10 (b-3) shows that the reconstruction with the true distribution of µ_a does not produce any artifacts. The large µ_a in the particular region hardly affects the image as shown in Fig. 10(c).

 

Fig. 9. The detected light intensities (dashed line) in the inhomogeneous medium and (solid) in the homogeneous medium with µ_a = 0.007 mm⁻¹. The absorption coefficien in (3/4)π ≤ θ ≤ (5/4)π of the inhomogeneous medium varied.

Download Full Size | PDF

 

Fig. 10. The source distributions reconstructed by use of incorrect homogeneous distribution of µ_a. λ_SDV used in the noise reduction and the µ′_s in the particular region in the inhomogeneous medium are (a) λ_SVD = 0.99, µ_a =0.0035 mm⁻¹, (b-1) λ_SVD = 0.97, µ_a =0.0007 mm⁻¹, (b-2) λ_SVD = 0.99, µ_a =0.0007 mm⁻¹. For (b-3), the correct distribution of µ_a and λ_SVD = 0.99 were used. µ_a =0.0007 mm⁻¹ in the particular region. For (c), λ_SVD = 0.995 and the homogeneous µ_a are used. µ_a of the particular region is 0.07 mm⁻¹.

Download Full Size | PDF

3.5. Spatial resolution of the reconstructed images

The spatial resolution of the proposed method is discussed in this section. We assume two circular emitting source regions with the radius of 2.5 mm where the unit source strengths are placed at the nodes. Gaussian random noises with SD_noise = 0.01 are added to the simulated measurement data.

When the centers of the true source regions are at (x,y) = (−10, 0) and (10, 0), the reconstruction with λ_SVD = 0.98 and seven updating iterations separates two emitting sources clearly as shown in Fig. 11 (a). The number of the updating iterations to obtain well localized source regions increases when the true sources are embedded deeply. When the centers of the two sources are at (x,y) = (−7.5, 0) and (7.5, 0), the reconstruction with λ_SVD of 0.97 and seven updating iterations localizes two emitting sources separately (Fig. 11(b)).

 

Fig. 11. The reconstructed (top) and true (bottom) distributions of the source strengths by the proposed method. The true centers of the two emitting source regions are at (a) (x, y) = (−10, 0) and (10, 0), (b) (x, y) = (−7.5, 0) and (7.5,0), (c) (x, y) = (−5, 0) and (5, 0), and (d) (x, y) = (−7.5, 20) and (7.5, 20).

Download Full Size | PDF

When the distance between the two sources decreases further, with their centers at (x,y) = (−5.0, 0) and (5.0, 0), the two emitting sources are not clearly separated although two peaks are recognizable as shown in Fig. 11 (c). The number of the iteration is seven, and the separation is not improved by the updating iteration anymore. Figure 11 (d) shows the reconstructed image with λ_SVD = 0.99 and six iterations when the centers of the emitting regions are (x,y) = (−7.5, 20) and (7.5, 20). Two emitting source regions in Fig. 11(d) are separated as well as in Fig. 11 (b). The proposed method has a potential to localize two emitting source regions separately when the two source regions are separated more than about 10 mm.

4. Advanced simulations of FDOT

In the previous sections, we have solved the PDE only for the emitting light with neglecting the excitation light. This is applicable to BDOT, but unrealistic for FDOT where exciting fluo rophores is inevitable. This section describes the simulation of FDOT starting from solving the PDE for the excitation light.

4.1. 2D reconstruction

4.1.1. Generation of the measurement data

A circular medium with the radius of 40 mm is divided into three regions with different optical properties. One is the sector region of the circular medium in the range of (3/4)π ≤ θ ≤ (5/4)π with µ_a=0.0070 mm⁻¹ and µ′_s = 0.72 mm⁻¹ at the excitation wavelength. In the second sector region in −(1/4)π ≤ θ ≤ (1/4)π, µ_a=0.0056 mm⁻¹ and µ′_s =0.80 mm⁻¹, and in the third section µ_a=0.0070 mm⁻¹ and µ′_s=0.80 mm⁻¹. Figures 12(a) and (b) show the true distributions of µ_a and D = 1/(3µ′_s), respectively.

The fluorophor assumed in this simulation is Indocyanine-green (ICG), of which the extinction coefficient ε, is 0.013 µM⁻¹mm⁻¹ and the quantum yield, η, is 0.016. We assume two regions with high concentrations of ICG. One is a circular region having a radius of 2.5 mm with its center at (x,y) = (0, 20) and with the concentration of ICG, N, of 5µM. Therefore, µ_a of this region increases by µ_af (= εN) = 0.065 mm⁻¹ to 0.072 mm⁻¹ due to ICG. The other is a circular region having a radius of 5 mm with its center at (x,y) = (0,−20) and with N of 1µM leading to µ_a in this circular region of 0.020 mm⁻¹. Besides, ICG is distributed uniformly in the medium with the concentration of 0.001µM which causes the background emission. The true distribution of µ_a + µ_af is shown in Fig. 12(c). µ_a and µ′_s at the emission wavelength are assumed to be equal to those at the excitation wavelength.

 

Fig. 12. (a) The distribution of µ_a of the circular medium, (b) the distribution of D = 1/(3µ′_s), (c) the distribution of (µ_a + µ_af) where µ_af is the absorption coefficien of ICG, (d) the calculated fluenc rate of the excitation light in the logarithmic scale, log₁₀ Φ_x, and (e) the fluorescenc light sources in the logarithmic scale, log₁₀ q ₀ = log₁₀ ηµ_afΦ_x

Download Full Size | PDF

The excitation CW light sources with unit strengths illuminate the medium at 16 positions simultaneously. The fluenc rate of the excitation light, Φ_x, is shown in Fig. 12(d). The fluenc rate of the fluorescenc light is calculated with the true fluorescenc light sources, q ₀ = ηµ_afΦ_x which is shown in Fig. 12(e). There exist strong emission from the two regions containing the concentrated ICG and the background emissions observed clearly near the excitation light sources. The detectors for the fluorescenc light are placed at 16 positions in the middle of the neighboring two excitation light sources. The excitation light is filtere out and does not contaminate the measurement data. Gaussian noises with SD_noise= 0.01 are added to the measurement data.

4.1.2. Reconstruction of the fluorescenc light sources

In the reconstruction, the uniform distributions of µ_a and µ′_s are assumed in the medium with the values of 0.007 mm⁻¹ and 0.8 mm⁻¹, respectively. Figures 13(a) to (e) show the reconstructed images when λ_SVD are (a) 0.995, (b) 0.9, and (c), (d) and (e) 0.8, respectively. The numbers of the updating iterations are six for cases (a) to (c), f ve for case (d), and seven for case (e).

Due to the uniform distribution of ICG, a strong artifact caused by the background emission appears at the center of the medium in Fig. 13(a) with the largest λ_SVD. Less strong artifacts appear in the left half of the circular medium due to the mismatch in µ′_s. The influenc of the mismatch in µ′_s in the left half is more significan than that in µ_a in the right half. Although the two emitting source regions are reconstructed, the effect of removing the noise space is not sufficient

 

Fig. 13. The reconstructed images of the fluorescen source strengths with (a) λ_SVD =0.995 and six iterations, (b) λ_SVD = 0.9 and six iterations, (c) λ_SVD = 0.8 and six iterations, (d) λ_SVD = 0.80 and f ve iterations, and (e) λ_SVD = 0.8 and seven iterations.

Download Full Size | PDF

Images shown in Figs. 13(b) and (c) are reconstructed with λ_SVD smaller than that in Fig. 13(a). Comparing with Fig. 13 (a), the artifact at the center of the medium is eliminated, and the localization of two emitting sources is improved. From these results, λ_SVD = 0.80 seems more preferable.

The effect of the number of the updating iterations is demonstrated by comparing Figs. 13(c), 13(d) and 13(e). The reconstructed image in Fig. 13(e) has no artifact and shows well-localized emitting regions. However, the sizes of the emitting regions in the reconstructed image with λ_SVD = 0.8 and six iterations shown in Fig. 13(c) are the closest to the true ones among the three cases (c), (d) and (e).

From these results it can be said that the proposed method can robustly reconstruct the image of the source strengths for FDOT, although λ_SVD and the number of iterations are somewhat influential

4.2. 3D reconstruction

Finally, three dimensional reconstruction is demonstrated. The experimental object is a cylinder with a radius of 15 mm and a height of 30 mm consisting of annular regions which simulate various tissues in a mouse as shown in Fig. 14. Each region has µ′_s and µ_a listed in Table 1 [42,43]. The fluoropho (ICG) is concentrated in the two spherical regions having the centers at (x,y) = (−6, 0, 12) and (6, 0, 18) with the radius of 1.5 mm. The ICG concentrations in the regions are given as 1 µM. The illuminating and the detecting positions are placed in the x−y planes of z = 9 mm and 21 mm. Each plane has eight illuminating and eight detecting positions with an equal spacing. The detectors are placed in the middle of the two neighboring illuminating positions. The measurement data are generated by the same manner as in the 2D reconstruction in the previous section. The reconstructions are conducted from the measurement data of the 16 measured fluorescenc light intensities. The FEM calculation employs 5631 nodes and 24000 elements for the generation of the measurement data.

 

Fig. 14. The cross section in the plane of y = 0 mm of the cylindrical object consisting of annular regions which simulate various tissues of a mouse, i.e., (i) heart, (ii) lung, (iii) liver and (iv) fat. The illuminating and detecting positions are placed in the planes of z = 9, and 21 mm (dashed lines).

Download Full Size | PDF

Table 1. The optical properties of the tissues in the cylindrical object [mm⁻¹].

View Table

Figure 15(a) shows the reconstructed image when the measurement data contain neither background emission nor Gaussian noise and when there exist no mismatch in the optical properties and no artifact reduction is carried out. The FEM mesh identical to that for generating the measurement data is used for image reconstruction. The emitting regions are clearly reconstructed at the correct positions, and no artifact is observed. The recovered source strengths are compared well with the maximum of the true strength of 6.40 × 10⁻⁶.

In the following cases (b), (c) and (d), 2969 nodes and 12288 elements are used for image reconstruction. The averages of the true optical properties, µ_a = 0.0018 mm⁻¹ and µ′_s =0.996 mm⁻¹, are used in the reconstruction. In addition to the two emitting regions, ICG is uniformly distributed with the concentration of 0.0001µM in the whole medium. The measurement data contain the Gaussian noise with SD_noise = 0.01. Figure 15(b) shows the images of the source strength reconstructed from the noisy measurement data with seven iterations and λ_SVD = 0.995. Strong artifacts are making the emitting regions unclear.

When λ_SVD =0.9 and 0.8, the emitting sources are successfully reconstructed from the noisy data as shown in Figs. 15 (c) and 15(d), respectively, although some artifacts are observed. Also the positions of the emitting sources are slightly shifted from the true positions partly due to the coarser mesh than that for Fig. 15(a).

From these results, the proposed reconstruction scheme is validated for 3D geometries as well as 2D geometries. We should note that the efforts to reduce the noises and the back ground emissions in data acquisition and to use fin meshes in calculation are important for successful reconstructions of real objects.

 

Fig. 15. The cross sections of the reconstructed 3D source distributions (a) from the data without background emission, noise and mismatch in the optical properties, and no artifact reduction is carried out, (b) to (d) from the data with background emissions, noises and mismatches in the optical properties, λ_SVD = (b) 0.995, (c) 0.9, and (d) 0.8. The number of the updating iteration is seven.

Download Full Size | PDF

5. Conclusions

We have proposed a new method for improved reconstruction of the source distributions by use of a spatial filte, successive updating of the forward model and reduction of artifacts for fluorescence/bioluminescenc diffuse optical tomography.

The spatial filte transforms a set of the measurement data to a single source strength at a position of interest by using the differences in the contribution of the sources at different positions to the measurement data. The estimated source strengths indicate the importance of the source positions to produce the measurement data. Utilizing these characteristics of the spatial filte, the forward model is updated. This process ignores the dispensable source positions from the forward model according to the reconstructed source distribution, and as the result, the illposedness of the inverse problem is alleviated. The spatial resolution of the reconstructed image is also improved. The proposed reduction of artifacts by removing the noise space based on the singular value decomposition technique also exploits the features of spatial filtering Assuming sufficien signal-to-noise ratio, the signal space can be distinguished from the noise space by using the estimated source strengths and the forward model. The effects of several influencin factors such as the background emissions, measurement noises and mismatch in the optical properties between the forward model and the real object are reduced effectively.

Numerical experiments show the differences in the results between the proposed method and some conventional ones. The advantages of the proposed method are the capability of localizing the deeply embedded sources and the robustness to the measurement noises and the calculation errors. The quantitative criteria for truncating the singular values and selecting the number of the updating iterations are exemplifie by the simulation. The limits of the proposed reconstruction scheme are also revealed by the simulation. It is concluded that the proposed method employing the spatial filte, updating of the forward model and removing of the noise space based on the singular value decomposition can effectively improve the images of fluores cence/bioluminescence diffuse optical tomography.