OSA's Digital Library

Virtual Journal for Biomedical Optics

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

  • Editor: Gregory W. Faris
  • Vol. 5, Iss. 10 — Jul. 19, 2010

Optics InfoBase > Virtual Journal for Biomedical Optics > Volume 5 > Issue 10 > Page 13151

« Show journal navigation

Reconstruction of fluorescence/bioluminescence sources in biological medium with spatial filter

Shinpei Okawa and Yukio Yamada  »View Author Affiliations


Optics Express, Vol. 18, Issue 12, pp. 13151-13172 (2010)
http://dx.doi.org/10.1364/OE.18.013151


View Full Text Article

Acrobat PDF (9789 KB)


Advanced Search

Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

We propose a new method for reconstruction of emitting source distributions by use of a spatial filte and a successive updating process of the forward model 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, and the forward model is updated by use of the estimated source strengths. This updating process ignores the dispensable source positions from reconstruction according to the reconstructed source distribution, and the spatial resolution of the reconstructed image is improved. The estimated sources are also used for the reduction of artifacts induced by noises based on the singular value decomposition. Some numerical experiments show the advantages of the proposed method by comparing the present results with those obtained by the conventional methods of the least squares method and Algebraic Reconstruction Technique. Finally the criteria for practical use of the method are quantitatively presented by the simulations for 2D and 3D geometries.

© 2010 Optical Society of America

1. Introduction

Fluorescence/ bioluminescence diffuse optical tomography (F/BDOT) has been studied for some medical applications [1

V. Ntziachristos, C-H. Yung, C. Bremerand, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo ,” Nat. Med. 8 (7), 757–760 (2002). [CrossRef] [PubMed]

,2

V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005). [CrossRef] [PubMed]

]. 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

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Prob. 15, R41–R93 (1999). [CrossRef]

]. 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

V. Ntziachristos, C-H. Yung, C. Bremerand, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo ,” Nat. Med. 8 (7), 757–760 (2002). [CrossRef] [PubMed]

]. 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

D. Grosenick, H. Wabnitz, H. H. Rinneberg, T. Moesta, and P. M. Schlag, “Development of a time-domain optical mammography and fir t in vivo applications,” Appl. Opt. 38(13), 2927–2943 (1999). [CrossRef]

]. 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

R. Weissleder, “Molecular Imaging in Cancer,” SCIENCE 321, 1168–1171 (2006). [CrossRef]

]. 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

K. Vishwanath, B. Pogue, and M-A. Mycek, “Quantitative fluore cence lifetime spectroscopy in turbid media: comparison of theoretical, experimental and computational method,” Phys. Med. Biol. 47, 3387–3405 (2002). [CrossRef] [PubMed]

].

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

J. Wu, Y. Wang, L. Perelman, I. Itzkan, R. R. Dasari, and M. S. Feld, “Time-resolved imaging of fluore cent objects embedded in turbid media,” Opt. Lett. 20(5), 489–491 (1995). [CrossRef] [PubMed]

]. 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

L. Zhang, F. Gao, H. He, and H Zhao, “Three-dimensional scheme for time-domain fluore cence molecular tomography based on Laplace transforms with noise-robust factors,” Opt. Express 16(10), 7214–7223 (2008). [CrossRef] [PubMed]

]. 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

F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, “A self-normalized, full time-resolved method for fluore cence diffuse optical tomography,” Opt. Express 16(17), 13104–13121 (2008). [CrossRef] [PubMed]

, 25

M. A. O’Leary, D. A. Boas, X. D. Li, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21(2), 158–160 (1996). [CrossRef] [PubMed]

]. 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

B. D. van Veen and K. M. Buckly, “Beamforming: A versatile approach to spatial filtering” IEEE ASSP Mag. 15, 4–23 (1988). [CrossRef]

]. Recently, the existing spatial filte is applied to the inverse problem in functional brain imaging, such as electroencephalography (EEG) and magnetoencephalography (MEG) [30–33

S. Baillet, J. C. Mosher, and R. M. Leahy, “Electromagnetic brain mapping,” IEEE Signal Process Mag. 18, 14–30 (2001). [CrossRef]

]. 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

S. Okawa and Y. Yamada, “Source estimation with spatial filte for fluore cence diffuse optical tomography,” in Biomedical Optics Topical Meeting Technical Digest (OSA, Miami, Florida, 2008) BSuE41.

]. It has been reported that a similar iterative refinemen of the model achieves high-quality reconstruction in 3D [35

C. Kuo, O. Coquoz, T. L. Troy, H. Xu, and B. W. Rice, “Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging,” J. Biomed. Opt. 12, 024007 (2007). [CrossRef] [PubMed]

].

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

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Prob. 15, R41–R93 (1999). [CrossRef]

],

{ − ∇ ·D (r) ∇+ μa (r)}Φ (r)= q0 (r),
(1)

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 0 is the fluorescenc light source. Equation (1) can be solved under the boundary condition expressed by Eq. (2),

−n ·D ∇Φ= 1 2AΦ on ∂Ω,
(2)

where A = (1 + Rf )/(1 − Rf ), Rf 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 0(r), is described as q 0 = ηµ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 0 when Φ x is given. Therefore, this study is focused on the reconstruction of q 0 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),

{K (D)+C ( μa)+B}Φ=Q,
(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

M. Schweiger, S. R. Arridge, and D. T. Delpy, “Application of the Finite-Element Method for the Forward and Inverse Model in Optical Tomography,” J. Math. Imaging Vis. 3, 263–283 (1993). [CrossRef]

]. The matrix B is for the boundary condition. The entries of the matrices are, Kij = ∫Ω D∇ψi ·∇ψjdr, Cij = ∫Ω µaψiψjdr and Bij = 1/(2A) ∫∂Ω ψiψjdγ, where i and j denote the indexes assigned to the nodes, Bij 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 Φ,

Φ=GQ,
(4)

where G = {K(D) + C(µa ) + B}−1.

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,

m=LQ,
(5)

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 wk 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̂ 1, q̂ 2, ⋯, q̂N ] T is obtained with the M × N matrix W which consists of the spatial filter wk ,

Q̂= WTm,
(6)
q̂k= wkTm= Σ i=1N wkT li qi, (k=1,2, ⋯,N),
(7)

where q̂k is the k-th source estimated by wk , which is the k-th column vector of W, represents the spatial filte targeting on the k-th source strength qk , and li is the i-th column vector of L. li 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 lj hold wT klj , 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 lj since the matrix L is singular. For successful estimation, wT k lj (k ≠ j) appearing in Eq. (7) should be close to zero, and wT k lk should be as close to unity as possible. In other words, the spatial filte is the solution of the following maximization problem,

max wk ( wkT lk)2,subject to ∥ wkTL∥2 (= Σ i=1N ( wk li)2)=1.
(8)

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

wk= { lkT ( LLT) −1 lk} − 12 ( LLT) −1 lk.
(9)

When (LLT ) is singular, its inverse matrix in Eq. (9) can be substituted by a pseudo-inverse matrix. The coefficien in Eq. (9), {lT (LLT )−1 li }−1/2, differentiates the spatial filter from the Moore-Penrose pseudo inversion.

When we solve the minimization problem, min wk ∥ wkTL∥2 subject to (wT klk ) = 1, wk = (LLT )−1 lk 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 (wT klk )2 ≥ (wT jlk )2, since the spatial filte wk maximizes the squared product of (wTlk )2, where w is the arbitrary vector satisfying ∥wTL∥2 = 1. If the source qk contributes to the measurement data to a great extent, i.e. m≃lkqk , the estimated source strength q̂k = wT 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 wT klk 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,

pi= q̂i ∥ Q̂∥∞,
(10)

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

li← li pi.
(11)

This process transforms our inverse problem formulated by Eq. (5) to that reformulated as m = LPQ, where P is the diagonal matrix with pi . The weight, pi , takes the value in the range of [0, 1]. pi = 1 means that the i-th source position is necessary for the forward model, while pi = 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,

α=arg minα ∥ m−αL Q̂∥2= mTL Q̂ Q̂T LTL Q̂,
(12)

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

Q¯=α · Q̂.
(13)

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,

m=LQ+ε,
(14)

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 liq̂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 LLT ). 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, liq̄i , first L is decomposed as L = UDVT 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 1, s 2, ⋯, sj , ⋯, sM and (N − M) components with zero in the diagonal. The M × 1 basis vectors, uj (j = 1, 2, ⋯ M) form the matrix U = [u 1,u 2,⋯,uM ]. 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 uj , i.e. m = Uy , 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,

m← UHUTm (=UHy),
(15)

where H is the M × M diagonal matrix with the j-th diagonal components Hjj 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−1, respectively, by assuming the application to breast cancer screening using fluorescen dye working in the near-infrared wavelength range [9

J. C. Hebden, H. Veenstra, H. Dehghani, E. M. C. Hillman, M. Schweiger, S. R. Arridge, and D. T. Delpy, “Three-dimensional time-resolved optical tomography of a conical breast phantom,” Appl. Opt. 40(19), 3278–32887 (2001). [CrossRef]

]. 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.
Fig. 2. The residual error ratio as the function of the updating iterations.

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̄∥2/∥m∥2 ×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).

As a typical conventional reconstruction method, we firs choose the least squares method which gives the least squares minimum norm solution [37

C. R. Vogel,Computational Methods for Inverse Problems (Frontiers in Applied Mathematics) (SIAM, Philadelphia, 2002). [CrossRef]

]. The solution is obtained as Q = LT (LLT )−1 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∥2 [38

C. C. Paige and M. A. Saunders, “LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares,” ACM trans. on Math. Software 8(1), 43–71 (1982). [CrossRef]

] 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

C. W. Groetsch, Inverse Problems in the Mathematical Sciences 1. Auflage (Friedr. Vieweg & Sohn Verlagsge-sellschaft mbH, Braunschweig/Wiesbaden, 1993).

, 40

C. W. Groetsch, Inverse problem (the Mathematical Association of America, Washington, 1999).

]. 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

F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, “A self-normalized, full time-resolved method for fluore cence diffuse optical tomography,” Opt. Express 16(17), 13104–13121 (2008). [CrossRef] [PubMed]

].

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

S. Holder, Electrical Impedance Tomography: Methods, History and Applications (Institute of Physics Publishing, Bristol and Philadelphia, 2005).

], 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 (SDnoise = 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 ui 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 SDnoise =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) SDnoise = 0.01 and λSVD = 0.98, (a-2) SDnoise = 0.01 and λSVD = 1.0, (b) SDnoise = 0.001 and λSVD = 0.995, and (c) SDnoise = 0.1 and λSVD = 0.96.

Figure 4 (b) shows that when the noise is sufficientl small, SDnoise =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 SDnoise = 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, qBG , 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.
Fig. 6. The reconstructed source distributions using the proposed noise reduction with (a-1) qBG = 0.0001 and λSVD = 0.995, (a-2) qBG = 0.0001 and λSVD = 0.97, (b-1) qBG = 0.001 and λSVD = 0.97, (b-2) qBG = 0.001 and λSVD = 0.80, and (c-1) qBG = 0.01 and λSVD = 0.97, (c-2) qBG = 0.01 and λSVD = 0.8.

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−1. 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−1. 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−1. 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−1 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−1 (solid curve). µ′ s in the region of (3/4)π ≤ θ ≤ (5/4)π of the inhomogeneous medium is 0.72 mm−1, 0.4 mm−1 or 0.08 mm−1.
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−1, (b) λSVD = 0.98, µ′ s = 0.4 mm−1, (c-1) λSVD = 0.95, µ′ s =0.08 mm−1, and (c-2) λSVD = 0.99, µ′ s =0.08 mm−1. For (c-3), the correct distribution of µ′ s and λSVD = 0.99 are used. µ′ s =0.08 mm−1 in the particular region. And for (d), homogeneous µ′ s and λSVD = 0.995 are used. The particular region has µ′ s =8.0 mm−1.

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−1 and the inhomogeneous media with µa of (a) 0.0035 or (b) 0.0007 mm−1 in the particular region. µ′ s is given as 0.8 mm−1 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−1 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−1. The absorption coefficien in (3/4)π ≤ θ ≤ (5/4)π of the inhomogeneous medium varied.
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−1, (b-1) λSVD = 0.97, µa =0.0007 mm−1, (b-2) λSVD = 0.99, µa =0.0007 mm−1. For (b-3), the correct distribution of µa and λSVD = 0.99 were used. µa =0.0007 mm−1 in the particular region. For (c), λSVD = 0.995 and the homogeneous µa are used. µa of the particular region is 0.07 mm−1.

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 SDnoise = 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).

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−1 and µ′ s = 0.72 mm−1 at the excitation wavelength. In the second sector region in −(1/4)π ≤ θ ≤ (1/4)π, µa =0.0056 mm−1 and µ′ s =0.80 mm−1, and in the third section µa =0.0070 mm−1 and µ′ s =0.80 mm−1. 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−1mm−1 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−1 to 0.072 mm−1 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−1. 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, log10 Φ x , and (e) the fluorescenc light sources in the logarithmic scale, log10 q 0 = log10 ηµaf Φ x

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 0 = ηµ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 SDnoise = 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−1 and 0.8 mm−1, 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.

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

S. Okawa and Y. Yamada, “3D Light Source Reconstruction with Spatial Filter for Fluorescence/ Bioluminescence Diffuse Optical Tomography,ffi Diffuse Optical Imaging II, R. Cubeddu and A. H. Hielscher, eds., Proc. SPIE 7369, 736916.

,43

Y. Lv, J. Yian, W. Cong, G. Wang, W. Yang, C. Qin, and M. Xu, “Spectrally resolved bioluminescence tomography with adaptive finit element analysis: methodology and simulation,” Phys. Med. Biol. 52, 4497–4512 (2007). [CrossRef] [PubMed]

]. 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).
Table 1.  The optical properties of the tissues in the cylindrical object [mm−1].
heartlungliverfat
µ′ s 0.822.040.61.14
µa 0.00140.00280.0030.0012

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−6.

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−1 and µ′ s =0.996 mm−1, 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 SDnoise = 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.

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.

References and links

1.

V. Ntziachristos, C-H. Yung, C. Bremerand, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo ,” Nat. Med. 8 (7), 757–760 (2002). [CrossRef] [PubMed]

2.

V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005). [CrossRef] [PubMed]

3.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Prob. 15, R41–R93 (1999). [CrossRef]

4.

A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). [CrossRef] [PubMed]

5.

M. S. Patterson and B. W. Pogue, “Mathematical model for time-resolved and frequency-domain fluore cence spectroscopy in biological tissues,” Appl. Opt. 33(10), 1963–1974 (1994). [CrossRef] [PubMed]

6.

D. Grosenick, H. Wabnitz, H. H. Rinneberg, T. Moesta, and P. M. Schlag, “Development of a time-domain optical mammography and fir t in vivo applications,” Appl. Opt. 38(13), 2927–2943 (1999). [CrossRef]

7.

D. Grosenick, K. T. Moesta, M. Möller, J. Mucke, H. Wabnitz, B. Gebauer, C. Stroszczynski, B. Wassermann, P. M. Schlag, and H. Rinneberg,“Time-domain scanning optical mammography: I. Recording and assessment of mammograms of 154 patients,” Phys. Med. Biol. 50, 2429–2449 (2005). [CrossRef] [PubMed]

8.

D. Grosenick, H. Wabnitz, K. T. Moesta, J. Mucke, P. M. Schlag, and H. Rinneberg, “Time-domain scanning optical mammography: II. Optical properties and tissue parameters of 87 carcinomas,” Phys. Med. Biol. 50, 2451–2468 (2005). [CrossRef] [PubMed]

9.

J. C. Hebden, H. Veenstra, H. Dehghani, E. M. C. Hillman, M. Schweiger, S. R. Arridge, and D. T. Delpy, “Three-dimensional time-resolved optical tomography of a conical breast phantom,” Appl. Opt. 40(19), 3278–32887 (2001). [CrossRef]

10.

T. Yates, C. Hebdan, A. Gibson, N. Everdell, S. R. Arridge, and M. Douek, “Optical tomography of the breast using a multi-channel time-resolved imager,” Phys. Med. Biol. 50, 2503–2517 (2005). [CrossRef] [PubMed]

11.

G. Boverman, E. L. Miller, A. Li, Q. Zhang, T. Chaves, D. H. Brooks, and D. A. Boas “Quantitative spectroscopic optical tomography of the breast guided by imperfect a priori structural information,” Phys. Med. Biol. 50, 3941–3956 (2005). [CrossRef] [PubMed]

12.

R. Weissleder, “Molecular Imaging in Cancer,” SCIENCE 321, 1168–1171 (2006). [CrossRef]

13.

D. J. Hawrysz and E. M. Sevick-Muraca, “Developments Toward Diagnostic Breast Cancer Imaging Using Near-Infrared Optical Measurements and Fluorescent Contrast Agents,” Neoplasia 2 (5), 388–417 (2000). [CrossRef]

14.

L. Bakker, M. van der Mark, M van Beek, M. van der Voort, G. Hooft, T. Nielsen, T. Koehler, R. Ziegler, K. Licha, and M. Pessel, “Optical Fluorescence Imaging of Breast Cancer,” in Proceedings of Biomedical Optics Topical Meeting (OSA, Miami, Florida, 2006) SH56.

15.

K. Vishwanath, B. Pogue, and M-A. Mycek, “Quantitative fluore cence lifetime spectroscopy in turbid media: comparison of theoretical, experimental and computational method,” Phys. Med. Biol. 47, 3387–3405 (2002). [CrossRef] [PubMed]

16.

J. Wu, Y. Wang, L. Perelman, I. Itzkan, R. R. Dasari, and M. S. Feld, “Time-resolved imaging of fluore cent objects embedded in turbid media,” Opt. Lett. 20(5), 489–491 (1995). [CrossRef] [PubMed]

17.

R. Roy, A. Godavarty, and E. M. Sevick-Muraca, “Fluorescence-enhanced optical tomography using referenced measurements of heterogeneous media,” IEEE trans. Med. Imaging. 22(7), 824–836 (2003). [CrossRef] [PubMed]

18.

A. B. Milstein, J. J. Scott, S. Oh, D. A. Boas, R. P. Millane, C. A. Bouman, and K. J. Webb, “Fluorescence optical diffuse tomography using multiple-frequency data,” J. Opt. Soc. Am. A 21(6), 1035–1049 (2004). [CrossRef]

19.

D. Y. Paithankar, U. A. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluore cent yield and lifetime from multiply scattered light reemitted from random medium,” Appl. Opt. 36(10), 2260–2272 (1997). [CrossRef] [PubMed]

20.

F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, “A self-normalized, full time-resolved method for fluore cence diffuse optical tomography,” Opt. Express 16(17), 13104–13121 (2008). [CrossRef] [PubMed]

21.

A. Marjono, A. Yano, S. Okawa, F. Gao, and Y. Yamada, “Total light approach of time-domain fluore cence diffuse optical tomography,” Opt. Express 16(19), 15268–15285 (2008). [CrossRef] [PubMed]

22.

L. Zhang, F. Gao, H. He, and H Zhao, “Three-dimensional scheme for time-domain fluore cence molecular tomography based on Laplace transforms with noise-robust factors,” Opt. Express 16(10), 7214–7223 (2008). [CrossRef] [PubMed]

23.

F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, “A linear, featured-data scheme for image reconstruction in time-domain fluore cence molecular tomography,” Opt. Express , 14(16), 7109–7124 (2006). [CrossRef] [PubMed]

24.

M. J. Eppstein, D. E. Doughety, D. J. Hawrysz, and E. M. Sevick-Muraka, “Three-Dimensional Baysian Optical Image Reconstruction with Domain Decomposition,” IEEE trans. Med. Imaging 20(3), 147–163 (2001). [CrossRef] [PubMed]

25.

M. A. O’Leary, D. A. Boas, X. D. Li, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21(2), 158–160 (1996). [CrossRef] [PubMed]

26.

A. Soubret, J. Ripoll, D. Yessayan, and V. Ntziachristos, “Three-dimensional fluore cent tomography in presence of absorption: Study of the normalized Born approximation,” in 2004 Biomedical Optics Topical Meeting Technical Digest (OSA, Miami, Florida, 2004) WB6.

27.

A. X. Cong and G. Wang, “A finite-elemen reconstruction method for 3D fluore cence tomography,” Opt. Express 13(24), 9847–9857 (2006). [CrossRef]

28.

G. Wang, W. Cong, K. Durairaj, X. Qian, H. Shen, P. Sinn, E. Hoffman, G. McLennan, and M. Henry, “ In vivo mouse studies with bioluminescence tomography,” Opt. Express 14(17), 7801–7809 (2006). [CrossRef] [PubMed]

29.

B. D. van Veen and K. M. Buckly, “Beamforming: A versatile approach to spatial filtering” IEEE ASSP Mag. 15, 4–23 (1988). [CrossRef]

30.

S. Baillet, J. C. Mosher, and R. M. Leahy, “Electromagnetic brain mapping,” IEEE Signal Process Mag. 18, 14–30 (2001). [CrossRef]

31.

B. D. van Veen, W. van Drongelen, M. Yuchtman, and A. Suzuki, “Localization of brain electrical activity via linear constrained minimum variance spatial filte,” IEEE trans. Biomed. Eng. 44(9), 867–880 (1997). [CrossRef] [PubMed]

32.

S. Okawa and S. Honda, “MEG Analysis with Spatial Filtered Reconstruction,” IEICE Trans. on Fundam. Electron. Commun. Comput. Sci. 89-A(5), 1428–1436 (2006). [CrossRef]

33.

S. Okawa and S. Honda, “Dipole estimation with a combination of noise reduction and spatial filte,” International Congress Series 1300, 249–252 (2007). [CrossRef]

34.

S. Okawa and Y. Yamada, “Source estimation with spatial filte for fluore cence diffuse optical tomography,” in Biomedical Optics Topical Meeting Technical Digest (OSA, Miami, Florida, 2008) BSuE41.

35.

C. Kuo, O. Coquoz, T. L. Troy, H. Xu, and B. W. Rice, “Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging,” J. Biomed. Opt. 12, 024007 (2007). [CrossRef] [PubMed]

36.

M. Schweiger, S. R. Arridge, and D. T. Delpy, “Application of the Finite-Element Method for the Forward and Inverse Model in Optical Tomography,” J. Math. Imaging Vis. 3, 263–283 (1993). [CrossRef]

37.

C. R. Vogel,Computational Methods for Inverse Problems (Frontiers in Applied Mathematics) (SIAM, Philadelphia, 2002). [CrossRef]

38.

C. C. Paige and M. A. Saunders, “LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares,” ACM trans. on Math. Software 8(1), 43–71 (1982). [CrossRef]

39.

C. W. Groetsch, Inverse Problems in the Mathematical Sciences 1. Auflage (Friedr. Vieweg & Sohn Verlagsge-sellschaft mbH, Braunschweig/Wiesbaden, 1993).

40.

C. W. Groetsch, Inverse problem (the Mathematical Association of America, Washington, 1999).

41.

S. Holder, Electrical Impedance Tomography: Methods, History and Applications (Institute of Physics Publishing, Bristol and Philadelphia, 2005).

42.

S. Okawa and Y. Yamada, “3D Light Source Reconstruction with Spatial Filter for Fluorescence/ Bioluminescence Diffuse Optical Tomography,ffi Diffuse Optical Imaging II, R. Cubeddu and A. H. Hielscher, eds., Proc. SPIE 7369, 736916.

43.

Y. Lv, J. Yian, W. Cong, G. Wang, W. Yang, C. Qin, and M. Xu, “Spectrally resolved bioluminescence tomography with adaptive finit element analysis: methodology and simulation,” Phys. Med. Biol. 52, 4497–4512 (2007). [CrossRef] [PubMed]

OCIS Codes
(100.3190) Image processing : Inverse problems
(170.3880) Medical optics and biotechnology : Medical and biological imaging

ToC Category:
Image Processing

History
Original Manuscript: March 25, 2010
Revised Manuscript: May 27, 2010
Manuscript Accepted: May 28, 2010
Published: June 3, 2010

Virtual Issues
Vol. 5, Iss. 10 Virtual Journal for Biomedical Optics

Citation
Shinpei Okawa and Yukio Yamada, "Reconstruction of fluorescence/bioluminescence sources in biological medium with spatial filter," Opt. Express 18, 13151-13172 (2010)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-18-12-13151


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. V. Ntziachristos, C-H. Yung, C. Bremerand, and R. Weissleder, "Fluorescence molecular tomography resolves protease activity in vivo," Nat. Med. 8 (7), 757-760 (2002). [CrossRef] [PubMed]
  2. V. Ntziachristos, J. Ripoll, L. V. Wang, R. Weissleder, "Looking and listening to light: the evolution of wholebody photonic imaging," Nat. Biotechnol. 23 (3), 313-320 (2005). [CrossRef] [PubMed]
  3. S. R. Arridge, "Optical tomography in medical imaging," Inverse Prob. 15, R41-R93 (1999). [CrossRef]
  4. A. P. Gibson, J. C. Hebden, and S. R. Arridge, "Recent advances in diffuse optical imaging," Phys. Med. Biol. 50, R1-R43 (2005). [CrossRef] [PubMed]
  5. M. S. Patterson and B. W. Pogue, "Mathematical model for time-resolved and frequency-domain fluorescence spectroscopy in biological tissues," Appl. Opt. 33(10), 1963-1974 (1994). [CrossRef] [PubMed]
  6. D. Grosenick, H. Wabnitz, H. H. Rinneberg, T. Moesta, and P. M. Schlag, "Development of a time-domain optical mammography and first in vivo applications," Appl. Opt. 38(13), 2927-2943 (1999). [CrossRef]
  7. D. Grosenick, K. T. Moesta, M. Moller, J. Mucke, H. Wabnitz, B. Gebauer, C. Stroszczynski, B. Wassermann, P. M. Schlag, and H. Rinneberg, "Time-domain scanning optical mammography: I. Recording and assessment of mammograms of 154 patients," Phys. Med. Biol. 50, 2429-2449 (2005). [CrossRef] [PubMed]
  8. D. Grosenick, H. Wabnitz, K. T. Moesta, J. Mucke, P. M. Schlag, and H. Rinneberg, "Time-domain scanning optical mammography: II. Optical properties and tissue parameters of 87 carcinomas," Phys. Med. Biol. 50, 2451-2468 (2005). [CrossRef] [PubMed]
  9. J. C. Hebden, H. Veenstra, H. Dehghani, E. M. C. Hillman, M. Schweiger, S. R. Arridge, and D. T. Delpy, "Three-dimensional time-resolved optical tomography of a conical breast phantom," Appl. Opt. 40(19), 3278-3288 (2001). [CrossRef]
  10. T. Yates, C. Hebdan, A. Gibson, N. Everdell, S. R. Arridge, and M. Douek, "Optical tomography of the breast using a multi-channel time-resolved imager," Phys. Med. Biol. 50, 2503-2517 (2005). [CrossRef] [PubMed]
  11. G. Boverman, E. L. Miller, A. Li, Q. Zhang, T. Chaves, D. H. Brooks, and D. A. Boas "Quantitative spectroscopic optical tomography of the breast guided by imperfect a priori structural information," Phys. Med. Biol. 50, 3941-3956 (2005). [CrossRef] [PubMed]
  12. R. Weissleder, "Molecular Imaging in Cancer," Science 321, 1168-1171 (2006). [CrossRef]
  13. D. J. Hawrysz, and E. M. Sevick-Muraca, "Developments toward Diagnostic Breast Cancer Imaging using Near-Infrared Optical Measurements and Fluorescent Contrast Agents," Neoplasia 2 (5), 388-417 (2000). [CrossRef]
  14. L. Bakker, M. van der Mark, M van Beek, M. van der Voort, G. Hooft, T. Nielsen, T. Koehler, R. Ziegler, K. Licha, and M. Pessel, "Optical Fluorescence Imaging of Breast Cancer," in Proceedings of Biomedical Optics Topical Meeting (OSA, Miami, Florida, 2006) SH56.
  15. K. Vishwanath, B. Pogue, and M-A. Mycek, "Quantitative fluorescence lifetime spectroscopy in turbid media: comparison of theoretical, experimental and computational method," Phys. Med. Biol. 47, 3387-3405 (2002). [CrossRef] [PubMed]
  16. J. Wu, Y. Wang, L. Perelman, I. Itzkan, R. R. Dasari, and M. S. Feld, "Time-resolved imaging of fluorescent objects embedded in turbid media," Opt. Lett. 20(5), 489-491 (1995). [CrossRef] [PubMed]
  17. R. Roy, A. Godavarty, and E. M. Sevick-Muraca, "Fluorescence-enhanced optical tomography using referenced measurements of heterogeneous media," IEEE Trans. Med. Imaging. 22(7), 824-836 (2003). [CrossRef] [PubMed]
  18. A. B. Milstein, J. J. Scott, S. Oh, D. A. Boas, R. P. Millane, C. A. Bouman, and K. J. Webb, "Fluorescence optical diffuse tomography using multiple-frequency data," J. Opt. Soc. Am. A 21(6), 1035-1049 (2004). [CrossRef]
  19. D. Y. Paithankar, U. A. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, "Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random medium," Appl. Opt. 36(10), 2260-2272 (1997). [CrossRef] [PubMed]
  20. F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, "A self-normalized, full time-resolved method for fluorescence diffuse optical tomography," Opt. Express 16(17), 13104-13121 (2008). [CrossRef] [PubMed]
  21. A. Marjono, A. Yano, S. Okawa, F. Gao, and Y. Yamada, "Total light approach of time-domain fluorescence diffuse optical tomography," Opt. Express 16(19), 15268-15285 (2008). [CrossRef] [PubMed]
  22. L. Zhang, F. Gao, H. He, and H Zhao, "Three-dimensional scheme for time-domain fluorescence molecular tomography based on Laplace transforms with noise-robust factors," Opt. Express 16(10), 7214-7223 (2008). [CrossRef] [PubMed]
  23. F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, "A linear, featured-data scheme for image reconstruction in timedomain fluorescence molecular tomography," Opt. Express,  14(16), 7109-7124 (2006). [CrossRef] [PubMed]
  24. M. J. Eppstein, D. E. Doughety, D. J. Hawrysz, and E. M. Sevick-Muraka, "Three-Dimensional Baysian Optical Image Reconstruction with Domain Decomposition," IEEE Trans. Med. Imaging 20(3), 147-163 (2001). [CrossRef] [PubMed]
  25. M. A. O’Leary, D. A. Boas, X. D. Li, B. Chance and A. G. Yodh, "Fluorescence lifetime imaging in turbid media," Opt. Lett. 21(2), 158-160 (1996). [CrossRef] [PubMed]
  26. A. Soubret, J. Ripoll, D. Yessayan, and V. Ntziachristos, "Three-dimensional fluorescent tomography in presence of absorption: Study of the normalized Born approximation," in 2004 Biomedical Optics Topical Meeting Technical Digest (OSA, Miami, Florida, 2004) WB6.
  27. A. X. Cong and G. Wang, "A finite-element reconstruction method for 3D fluorescence tomography," Opt. Express 13(24), 9847-9857 (2006). [CrossRef]
  28. G. Wang, W. Cong, K. Durairaj, X. Qian, H. Shen, P. Sinn, E. Hoffman, G. McLennan, and M. Henry, "In vivo mouse studies with bioluminescence tomography," Opt. Express 14(17), 7801-7809 (2006). [CrossRef] [PubMed]
  29. B. D. van Veen and K. M. Buckly, "Beamforming: A versatile approach to spatial filtering," IEEE ASSP Mag. 15, 4-23 (1988). [CrossRef]
  30. S. Baillet, J. C. Mosher and R. M. Leahy, "Electromagnetic brain mapping," IEEE Signal Process Mag. 18, 14-30 (2001). [CrossRef]
  31. B. D. van Veen, W. van Drongelen, M. Yuchtman, and A. Suzuki, "Localization of brain electrical activity via linear constrained minimum variance spatial filter," IEEE Trans. Biomed. Eng. 44(9), 867-880 (1997). [CrossRef] [PubMed]
  32. S. Okawa, and S. Honda, "MEG Analysis with Spatial Filtered Reconstruction," IEICE Trans. on Fundam.Electron. Commun. Comput. Sci. 89-A(5), 1428-1436 (2006). [CrossRef]
  33. S. Okawa, and S. Honda, "Dipole estimation with a combination of noise reduction and spatial filter," International Congress Series 1300, 249-252 (2007). [CrossRef]
  34. S. Okawa, and Y. Yamada, "Source estimation with spatial filter for fluorescence diffuse optical tomography," in Biomedical Optics Topical Meeting Technical Digest (OSA, Miami, Florida, 2008) BSuE41.
  35. C. Kuo, O. Coquoz, T. L. Troy, H. Xu and B. W. Rice, "Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging," J. Biomed. Opt. 12, 024007 (2007). [CrossRef] [PubMed]
  36. M. Schweiger, S. R. Arridge, and D. T. Delpy, "Application of the Finite-Element Method for the Forward and Inverse Model in Optical Tomography," J. Math. Imaging Vis. 3, 263-283 (1993). [CrossRef]
  37. C. R. Vogel,Computational Methods for Inverse Problems (Frontiers in Applied Mathematics) (SIAM, Philadelphia, 2002). [CrossRef]
  38. C. C. Paige and M. A. Saunders, "LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares," ACM trans. on Math. Software 8(1), 43-71 (1982). [CrossRef]
  39. C. W. Groetsch, Inverse Problems in the Mathematical Sciences 1. Auflage (Friedr. Vieweg & Sohn Verlagsgesellschaft mbH, Braunschweig/Wiesbaden, 1993).
  40. C. W. Groetsch, Inverse Problem (the Mathematical Association of America, Washington, 1999).
  41. S. Holder, Electrical Impedance Tomography: Methods, History and Applications (Institute of Physics Publishing, Bristol and Philadelphia, 2005).
  42. S. Okawa and Y. Yamada, "3D Light Source Reconstruction with Spatial Filter for Fluorescence/ Bioluminescence Diffuse Optical Tomography,ffin Diffuse Optical Imaging II, R. Cubeddu and A. H. Hielscher, eds., Proc. SPIE 7369, 736916.
  43. Y. Lv, J. Yian, W. Cong, G. Wang, W. Yang, C. Qin and M. Xu, "Spectrally resolved bioluminescence tomography with adaptive finite element analysis: methodology and simulation," Phys. Med. Biol. 52, 4497-4512 (2007). [CrossRef] [PubMed]

Cited By

 Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited