OSA's Digital Library

Biomedical Optics Express

Biomedical Optics Express

  • Editor: Joseph A. Izatt
  • Vol. 1, Iss. 5 — Dec. 1, 2010
  • pp: 1259–1277
« Show journal navigation

Bioluminescence tomography with Gaussian prior

Hao Gao, Hongkai Zhao, Wenxiang Cong, and Ge Wang  »View Author Affiliations


Biomedical Optics Express, Vol. 1, Issue 5, pp. 1259-1277 (2010)
http://dx.doi.org/10.1364/BOE.1.001259


View Full Text Article

Acrobat PDF (1406 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

Parameterizing the bioluminescent source globally in Gaussians provides several advantages over voxel representation in bioluminescence tomography. It is mathematically unique to recover Gaussians [Med. Phys. 31(8), 2289 (2004)] and practically sufficient to approximate various shapes by Gaussians in diffusive medium. The computational burden is significantly reduced since much fewer unknowns are required. Besides, there are physiological evidences that the source can be modeled by Gaussians. The simulations show that the proposed model and algorithm significantly improves accuracy and stability in the presence of Gaussian or non- Gaussian sources, noisy data or the optical background mismatch. It is also validated through in vivo experimental data.

© 2010 OSA

1. Introduction

Bioluminescence tomography (BLT) [1

1. C. H. Contag and B. D. Ross, “It’s not just about anatomy: in vivo bioluminescence imaging as an eyepiece into biology,” J. Magn. Reson. Imaging 16(4), 378–387 (2002). [PubMed]

3

3. G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. 31(8), 2289–2299 (2004). [PubMed]

] is an emerging molecular imaging modality, which can be used to monitor physiological and pathological activities at molecular levels, specially visualize primary tumor growth and tumor cell metastasis. In BLT [4

4. 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). [PubMed]

13

13. Y. Lv, J. Tian, W. Cong, and G. Wang, “Experimental study on bioluminescence tomography with multimodality fusion,” Int. J. Biomed. Imaging 2007, 86741 (2007). [PubMed]

], one attempts to recover both location and intensity of bioluminescent sources from boundary measurements.

Conventionally the bioluminescent source is represented in voxels. For example, in a triangulation {τi, i≤N} of the domain Ω, the source in piecewise constants is
q=i=1Nqi1τi,
(1)
where 1i is the basis function supported on the subdomain τi that is 1 on τi, and 0 otherwise, and qi is the average intensity value on τi to be recovered.

The outcome of the Gaussian model is formulated as a minimization problem with respect to the number of Gaussians, intensity, the center, the radius and the aligned direction of each Gaussian.

Just to clarify, although the number of Gaussians may be fixed in a specific reconstruction process, the exact number of actual Gaussian sources is not assumed. If the number of Gaussians in reconstruction is larger than the actual number, the source can still be recovered; otherwise, the actual number can be approached from a smaller initial guess in a combinatorial fashion through minimizing the data discrepancy. This can be illustrated in the following example. Suppose there are three Gaussians. In the reconstruction, we can set the number of Gaussians to four. With the proposed reconstruction algorithm, three recovered Gaussians would correspond to actual ones while the other one has the zero or nearly zero intensity values. On the other hand, if the assumed number is smaller than the actual number, e.g., only one Gaussian, the data discrepancy cannot be effectively reduced, and in this case the problem can be formulated as combinatorial optimization problem in which the number of Gaussians is increased iteratively until three so that data discrepancy can be reduced to the desired tolerance. On the other hand, a rough estimate of the actual number of Gaussians does improve both accuracy and efficiency of the reconstruction. In summary, the actual number of Gaussians is not assumed and it is practically sufficient to fix the number of Gaussians in reconstructions since the sources are localized and limited in quantity.

Last, although we model the source in Gaussians in the reconstruction, the actual source distribution in simulations is not necessarily Gaussian. However the reconstruction is better if the actual source is in Gaussians as well. Please refer to the result section for details.

Here is the outline of the paper. In section 2, we introduce the new model of bioluminescent source by Gaussians, various intrinsic geometric constraints, and then the minimization algorithm; in section 3, the proposed model and algorithm is validated through simulations and experimental data.

2. Methods

2.1. Forward Modeling

In the following DA with Robin boundary condition, we use f for the light intensity, q for the bioluminescent source, μa for the absorption coefficient, μs' for the reduced scattering coefficient, D for the diffusion coefficient, i.e., D = 1/ [3(μa + μs')], n^ for the outgoing normal at the medium boundary, A for the boundary constant [24

24. A. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).

].

(Dφ)+μaφ=qφ+2ADφn^=0
(2)

Here we choose piecewise linear finite element method (FEM) to solve Eq. (2). After discretization, Eq. (2) is reduced to a linear system with the discretized light intensity Φ, the discretized source Q, and the system matrix F dependent on μa and D, i.e.,

FΦ=Q
(3)

We refer the readers to [24

24. A. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).

] for the details of FEM solution of Eq. (2). After solving Eq. (3), we have the following boundary fluxes

f=Dφn^=φ/2A.
(4)

2.2. Source Representation by Gaussians

As a quantitative imaging method, BLT is expected to provide the source intensity or power accurately. However, there are some intrinsic aforementioned difficulties in achieving high-resolution reconstruction with the voxel-based BLT. Therefore, instead of local voxel representation of the bioluminescent source, we consider the global representation by characteristic shapes of the tumor, e.g., the Gaussians in this study. A motivation is that because the light of interest is highly diffusive, a localized bioluminescent tumor does look like a Gaussian if it is sufficiently distant from the boundary. It turns out that this new representation with significantly reduced degrees of freedom is more robust than the voxel-based representation.

Mathematically, it was proven in [3

3. G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. 31(8), 2289–2299 (2004). [PubMed]

]: the solution of BLT is non-unique without the incorporation of effective a priori knowledge on the source distribution; in particular the uniqueness can be established when source is assumed to be composed of solid balls with the known intensities. However, using solid balls is not effective in representing general source distribution, especially features like anisotropy and fast decay of source intensity which are quite common in practice. In this study, we use summation of Gaussians to model general source distribution. This structural representation is flexible and improves the reconstruction stability for BLT. In particular the anisotropic Gaussian representation is an effective mathematical approximation in the sense that it parameterizes the major quantitative information of interest, such as the intensity, the center and the size of the source.

Physiologically, the Gaussian shape is indeed a natural choice. According to http://www.humpath.com/tumorigenesis, tumorigenesis is a collection of complex genetic diseases characterized by multiple defects in the homeostatic mechanisms that regulate cell growth, proliferation and differentiation. Cancer is caused by uncontrolled proliferation and the inappropriate survival of damaged cells, which results in tumor formation. Cells have developed several safeguards to ensure that cell division, differentiation and death occur correctly and in a coordinated fashion, both during development and in the adult body. Many regulatory factors switch on or off genes that direct cellular proliferation and differentiation. Damage to these genes, which are referred as tumor-suppressor genes and oncogenes, is selected for in cancer. Most tumor-suppressor genes and oncogenes are first transcribed from DNA into RNA, and are then translated into protein to exert their effects. Recent evidence indicates that small non-protein-coding RNA molecules, called microRNAs (miRNAs), might also function as tumour suppressors and oncogenes. Therefore, in the tumorigenesis stage, cancerous cells are localized, limited in quantity, irregular in shape, and with some decaying behavior away from the center. As a result, they can be modeled effectively as Gaussian sources after they are bioluminescently labeled. Other inflammation foci can be similarly modeled at earlier stages.

Therefore we represent the source q as the summation of Gaussians with the intensity ρ, the center (xc, yc, zc), anisotropic radius (rx, ry, rz) and Euler angles (j, f, q), i.e.,
q(x,y,z)=i=1nρiExp(XiTRiTΣiRiXi).
(5)
where Xi=[xxc,iyyc,izzc,i],Σi=[1/rx,i20001/ry,i20001/rz,i2]andRi=RiθRiϕRiφwithRiθ=[cosθisinθi0sinθicosθi0001],Riϕ=[1000cosϕisinϕi0sinϕicosϕi] and Riφ=[cosφi0sinφi010sinφi0cosφi]. For simplicity, we restrict the discussion on the rotation variables in 2D, which is simplified to q only, i.e., Ri=[cosθisinθisinθicosθi].

Computationally, the Gaussian source representation Eq. (5) reduces the burden significantly since we solve BLT with at most n´(1 + d)(2 + d)/2 unknowns in d dimensions, which is in general much smaller than that from voxel representation, e.g., N (total number of voxels) unknowns for Eq. (1).

2.3. Gaussian-based BLT

Modeling bioluminescent source in Gaussians, we formulate BLT as the minimization problem with respect to variables XG = {ρj, xc,j, yc,j, zc,j, rx,j, ry,j, rz,j, jj, fj, qj, j≤n},

XG=argminXG12i=1M[PiTφ(XG)fi]2+R(XG).
(6)

Here we use the set {fi, i≤M} for the measured data, the set {Pi, i≤M} of column vectors for the discretized measuring operator, and geometric constraints R on XG, which can be regarded as the regularization term. Please note that f is now nonlinearly dependent on XG in this new formulation, although it linearly depends on q in Eq. (1) in the voxel-based representation.

In Eq. (6), we do not treat n as a variable, mainly because there is generally no practical need for that since the tumors are usually localized and limited in quantity. We can just assign a practical estimate to n that is larger enough than or comparable with the actual number.

However, in the case that such an estimate is unavailable, we formulate Eq. (6) as a combinatorial optimization problem with respect to the variable n. That is we start from an initial guess n0 of n, which can be smaller than the actual number, and solve Eq. (6) iteratively until the data discrepancy d reaches certain tolerancee. Here the update of n is simply nk + 1 = nk + 1. For the notation convenience, we will use X for XG whenever it is appropriate. Then we have the following algorithm.

  • Algorithm 1 : Combinatorial optimization
  • Given: initial guess n0 and e = 0.05.
  • Repeat: 1. Computing Xk through Eq. (6) with the fixed nk
  • 2. d = ||fk ( Xk ) f||/||f||;
  • 3. Stopping criterion: Quit if d <e.
  • 4. nk + 1 = nk + 1

In Algorithm 1, e = 0.05 is an empirical number that can identify those n’s that are sufficient for the reconstruction; fk denotes the measurements with the source Xk; || × || can be the simple summation of the absolute values.

Next, we turn to the geometric constraints R, which can be imposed naturally from the known geometry of the medium. Please note that we do not assume any a priori knowledge on the anatomical structure of the medium, e.g., from other imaging modalities, although they can be easily incorporated into the scheme to improve the reconstruction. All the geometric constraints are either from the shape of the medium or from some commonsense assumptions on the shape of the source to avoid some non-uniqueness in geometric representation.

The first type of constraints is the min-max constraint, the minimum and maximum for each variable of XG. For example, we can use the size of the medium as the min-max constraint for the center, and realistic values for anisotropic radius or intensity such as zeros as the lower bounds. That is

XminG<X<XmaxG.
(7.1)

Notice that although we constraint the rotation angles, e.g., £q<180°, the representation of ellipsoid is not unique, e.g., by simultaneously exchanging rx and ry and considering the supplementary angle of q in 2D. However, this does not affect the reconstruction result.

Secondly, we restrict the shape variation along different directions so that the source is not too “narrow”, which is physiologically reasonable, i.e.,

rx<cxyry,ry<cyzrz,rz<czxrx.
(7.2)

Here the parameter cxy, cyz and czx control the size difference between any two directions for each source.

In the case with multiple Gaussians, it is important to impose the following type of constraints, minimal-separation constaints, on centers and radiuses between any two Gaussians so that they will not overlap or get too close to each other. That is for any two Gaussian i and j,

cx(rx,i+rx,j)2+cy(ry,i+ry,j)2+cz(rz,i+rz,j)2<(xixj)2+(yiyj)2+(zizj)2.
(7.3)

Here cx, cy and cz are the control parameters, which can for example be set to ln2, corresponding to full width at half maximum (FWHM). Without using constraints from Eq. (7.2), multiple inclusions may not be separated. An example is given in the result section.

All the above constraints Eq. (7) can be represented by the following abstract inequalities

gk(XG)<0.
(8)

R(XG)=kln[gk(XG)].
(9)

Please note that when anatomical priors are available from other imaging modalities, such as X-ray CT or MRI, they can be incorporated naturally into the scheme Eq. (6). For example, the min-max constraints Eq. (7.1) can be tightened for centers or radiuses with the measured values from anatomical images; more accurate constraints can be imposed by assigning the appropriate constants in Eq. (7.2) and Eq. (7.3). In contrast, it is difficult to impose those global geometric features locally in the voxel-based BLT. On the other hand, the straightforward co-registration of the images from X-ray CT or MRI may not improve BLT satisfactorily since there is usually no simple one-to-one anatomical correspondence between them due to the fundamentally distinct imaging mechanisms.

2.4. Minimization by Barrier Method

As the standard approach for minimizing nonlinear least squares, we linearize the first term (the data fidelity) in Eq. (6), and solve for the incremental change dX iteratively via the following outer loop of iterations

Xn+1=Xn+argminδX12||JδX(ffn)||22+R(Xn+δX).
(10)

During each iteration in the outer loop for Eq. (10), we minimize the following with b = f-fn and x = dX

x=argminx12||Jxb||22+R(Xn+x).
(11)

minxL(x,t)=minx12||Jxb||221tkln[gk(Xn+x)].
(12)

Minimizing Eq. (12) strictly enforces the geometric constraint R since the value of logarithmic penalty functions would otherwise become infinitely large. During each step the solution xn is no more than K/t-suboptimal with K as the total number of constraints. This implies xn converges to the optimal point of Eq. (11) such that Jx = b with strictly enforced geometry constraints as t→∞. The Eq. (12) can be minimized with Newton’s method via the following inner loop

x2L(xn,tn)δx=xL(xn,tn)xn+1=xn+sδxtn+1=μtn.
(13)

2.5. Algorithm Implementation

The flowchart of the algorithm for solving Eq. (6) is as follow: the outer loop comes from the linearization Eq. (10) of the data fidelity term while the inner loop solves each linearization step using barrier method via Eq. (13).

  • Algorithm 2 : minimizing BLT Eq. (6) with Gaussian representation
  • Outer loop: Linearization of the data fidelity via Eq. (10)
  • Given: initial guess X0 and eo = 0.01.
  • Repeat: 1. Compute Jacobian J from Xn.
  • 2. Inner loop: solve x = dX using barrier method via Eq. (13)
  • Given: t0 = −2R(Xn)/||b||2, μ = 2, initial guess x0 = 0.
  • Repeat: 2.1. Compute the descent direction dx;
  • 2.2. Compute the moving step s via backtracking line search;
  • Given: a = 0.01, s = 1, b = 0.5.
  • While: L(xn + sdx)>L(xn) + α(ÑLx)T × (sdx)
  • s = bs.
  • 2.3. Update xn + 1 = xn + sdx and tn + 1 = μtn.
  • 2.4. Stopping criterion: Quit if K/tn<ei.
  • 3. Stopping criterion: Quit if |En + 1-En|/ En <eo.

In this flowchart, we use the ratio difference of the total source power E = òqdΩ as the stopping criterion for the outer loop, i.e., assuming that the total power E is stable when the algorithm converges. As the stopping criterion for the inner loop, we use the sub-optimal ratio K/t, which naturally measures the difference between the iterative solution of Eq. (12) and the original solution with strictly enforced constraints. We set ei empirically to 0.0001t0. Please note that the formula for t0 is for balancing the linearized data fidelity and the penalty function.

Next, we give the details of computing Jacobian and the descent direction in Algorithm 2 for the completeness.

Let us start with the computation of Jacobian J = {Jij, i£M, j£n}, which can be derived from the Jacobian J0 = {J0,ij, i£M, j£N} with respect to the piecewise-constant voxel representation via Eq. (1). The connection between J and J0 is through the following coordinate transformation
Jij=PiTφxj=k=1NqkxjPiTφqk=k=1NqkxjJ0,ik,
(14)
where qk represents the kth voxel value in Eq. (1), xj is the jth parameter in global representation by Gaussians in Eq. (5).

An efficient way to compute J0 is through the adjoint method [24

24. A. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).

]. After computing J0, J follows immediately from (14). However, noticing that there are only a few parameters in the new source representation by Gaussians, i.e., n†M†N, we directly compute J without computing J0 via the direct method as follow.

From Eq. (3), we have F(∂ϕ/∂qk) = ∂Q/∂qk, then

Jij=PiT(k=1Nqkxjφqk)=PiTk=1Nqkxj[F1(Qqk)]=PiT[F1(k=1NqkxjQqk)].
(15)

Through Eq. (15), we only need to compute the forward solver n times, which would be at least M times if we compute J0 first via Eq. (14). Please note that ∂Q/∂qk is a sparse vector, which is nonzero only at the nodes in the kth voxel; on the other hand ∂qk/∂xj can be analytically derived from Eq. (5). Therefore, we can now compute the Jacobian J merely by computing n forward solvers, which is much more efficient than that in the voxel-based BLT. The trade-off is that we need to compute the Jacobian more than once due to the nonlinear dependence of J on XG, which needs to be computed only once in the voxel-based BLT due to the linear dependence of J0 on X. However, since n†M, we still achieve a considerable gain in speed, since the new algorithm usually converges in less than 20 iterations while the ratio M/n is usually at least about a hundred. That is we usually achieve at least 5 gains in speed on the computation of Jacobians.

In the rare cases when n>M, the adjoint method is preferred, i.e.,
Jij=[(FT)1Pi]T(k=1NqkxjQqk).
(16)
That is we compute M forward solvers instead.

2Rk=(Rk)T(Rk).
(17)

Therefore, we have
xL(x,t)=JTb+1tkRkx2L(x,t)=JTJ+1tk(Rk)T(Rk).
(18)
Please note each constraint Rk is evaluated at Xn + x rather than x.

The overall algorithm for BLT with Gaussian representation is more efficient than voxel-based BLT, mainly because the minimization is now with respect to only a few parameters of Gaussians instead of at least thousands of voxels in the conventional voxel-based BLT. For the current BLT algorithm with new model by Gaussians, minimization through barrier method in Step 2 is extremely fast although the algorithm seems complicated, which is understandable due to the significantly reduced number of variables; the major computation cost comes from the Jacobian J in Step 1, although it takes less computational time than that in voxel-based BLT, which is again due to the reduced number of unknowns.

3. Results

In this section, we first compare the performance of Gaussian-based BLT, “G-BLT”, with that of the conventional voxel-based BLT, “V-BLT”, through simulated data in single-inclusion and multiple-inclusion cases, and then validate the proposed method in mouse study.

3.1. Reconstruction with Single Inclusion

We first evaluate the proposed algorithm with single anisotropic Gaussian inclusion, and then perform the single-inclusion reconstructions with various shapes, different noise levels, and different approximation errors in optical background to compare the performance of G-BLT and V-BLT. For the presentation purpose, we only show 2D comparisons. The results and the conclusions are similar for the 3D case. In particular, 3D results will be demonstrated with in vivo data.

All the 2D reconstructions are performed on a circular phantom with the diameter of 20 millimeters. 120 data are collected all around the phantom at the boundary with equal angular distance in between. The reconstruction mesh has 2162 triangular elements and 1130 nodes. To avoid the inverse crime, the data are generated with different meshes that are about four times as large as the reconstruction mesh. For example, there are 8664 triangular elements and 4425 nodes for the case with circular inclusion as shown in Fig. 1(a)
Fig. 1 Reconstructions with single Gaussian inclusion. (a)-(d) are the phantoms with 0°, 45°, 90° and 135° rotations respectively; (e)-(h) are from the voxel-based BLT; (i)-(l) are from Gaussian-based BLT. Please see Table 1 for the true and recovered parameters. Please note that the maximum intensity recovered via the voxel-based BLT is only up to 20% of the true value. The results show that Gaussian-based BLT not only localizes the inclusion better than, but also quantitatively more accurate than the voxel-based BLT.
.

We assume the homogeneous optical background with μa = 0.01mm−1 and μs' = 1mm−1. Without further mentioning, we use mm as the unit for the length, mm−1 as the unit for absorption or scattering coefficients, nano Watts (nW) as the unit for the total source power, and nW/mm2, nW/mm3 as the unit for the source intensity in 2D and 3D respectively. The unit of angle is in degree.

For the V-BLT, we solve the problem via the standard Levenberg-Marquardt method with L2 regularization [24

24. A. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).

]. That is to minimize
X=argminX12||J0Xf||22+λ||X||22.
(19)
where the Jacobian J0 is with respect to the piecewise-constant voxel representation by Eq. (1), and l is the regularization parameter.

Please notice that one may optimize the choice of regularization term in Eq. (19). However we found out that the results did not improve drastically when using other regularization strategies [14

14. H. Gao and H. K. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization,” Opt. Express 18(3), 1854–1871 (2010). [PubMed]

,15

15. H. Gao and H. K. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation part 2: total variation and l1 data fidelity,” Opt. Express 18(3), 2894–2912 (2010). [PubMed]

], which may be due to the severely illposed nature of the inverse problem in the optical regime modeled by diffusion approximation Eq. (2). In the following, we restrict our discussion based on L2 regularization as in Eq. (19).

3.1.1. Inclusions with Gaussian Shapes

Since it is physiologically reasonable to approximate the bioluminescence sources by Gaussians, we first consider reconstructions where we assume the source is a single Gaussian. All the tests have the same Gaussian up to a rotation, with parameters specified in Table 1

Table 1. Tests with single Gaussian source. The Gaussian phantom has ρ = 1, xc = 5, yc = 0, rx = 2, ry = 1. The rotation angles are 0°, 45°, 90°, 135°respectively for Test 1-4. In all tests, the total energy of the phantom E is approximately 6.285. Please notice that the result from Case 2 can be understood according to the non-unique representation of the Gaussian

table-icon
View This Table
| View All Tables
. The reconstructed results with both G-BLT and V-BLT are plotted in Fig. 1 and the recovered parameter values are in Table 1 as well.

From Fig. 1, in all tests G-BLT can recover the Gaussian phantom when using the representation of source in Gaussians while V-BLT gives relatively much poorer results. Moreover, G-BLT is quantitatively better than V-BLT in the sense it gives more accurate total power E and maximal intensity as documented in Table 1.

3.1.2. Inclusions with Various Shapes

We evaluate the algorithm through the simulated data with single inclusions of various shapes. The parameters for various single inclusions are specified in Table 2

Table 2. The parameters for single-object simulations with inclusions of various shapes. E is the total source energy; ρ is the maximal source intensity; xc and yc are x- and y-coordinate of the center of the inclusion; rx and ry are maximum distances from the center to the boundary along x and y direction respectively

table-icon
View This Table
| View All Tables
. The reconstructed results with both G-BLT and V-BLT are plotted in Fig. 2
Fig. 2 Single-inclusion reconstructions with various shapes. (a)-(d) are the phantoms with circular, square, triangular and rectangular inclusion respectively; (e)-(h) are from the voxel-based BLT; (i)-(l) are from Gaussian-based BLT. Please see Table 2 and 3 for the true and recovered parameters. Please note that the maximum intensity recovered via the voxel-based BLT is only up to 20% of the true value. The results show that Gaussian-based BLT not only localizes the inclusion better than, but also quantitatively more accurate than the voxel-based BLT.
and the parameter values are in Table 3

Table 3. The reconstructed parameters for single-object simulations with inclusions of various shapes. Please see Table 2 for the corresponding true values

table-icon
View This Table
| View All Tables
.

Before doing the reconstruction, we simply compare the measurements from Gaussian and non-Gaussian sources as shown in Fig. 3
Fig. 3 Comparison of boundary measurements from Gaussian sources and non-Gaussian sources. (a)-(c) are from non-Gaussian sources with side length ratio of 1:1, 1:2 and 2:1. (d)-(f) are from Gaussian sources with the radius ratio of 1:1, 1:2 and 2:1. (g)-(i) plot the boundary measurements originating counterclockwise from x-axis, in which the red curve are from Gaussian sources and the blue curves are from non-Gaussian sources. Please note that the data are normalized so that the absolute sum is 1.
. Here we embed inclusions deeply at the center of the medium. The plotted measurements are normalized with the absolute sum 1. The figure shows that the measurement profiles from Gaussian sources are almost the same as those from the non-Gaussian ones. It suggests that various shapes can be characterized by anisotropic Gaussians when the sources are deeply embedded into the diffusive medium.

From Fig. 2, visually we can tell that G-BLT is consistently better in localizing the inclusions than V-BLT. Moreover, G-BLT is quantitatively better than V-BLT in the sense it gives more accurate total power E and maximal intensity as documented in Table 3.

One minor issue regarding G-BLT is that the reconstructed radiuses along x- and y- axis seem not to match those in true phantoms, whish is possibly due to the Gaussian approximation error of non-Gaussian phantom. For example, with rectangular inclusion, the recovered rx = 1.752 and ry = 1.797, while the true rx = 1 and ry = 2. Although rx<ry in the reconstruction, this scale difference is much smaller than the true difference.

Notice that in this section and the rest of reconstructions for non-Gaussian sources, we did not consider q as reconstruction variable. The reason is that G-BLT does not recover the shapes for non-Gaussian sources as exactly as for Gaussian sources even with rotation.

3.1.3. Different Noise Level

We evaluate the algorithm with respect to different noise level. 5%, 10%, 20%, 30% percentage of Gaussian noise is added to the measurements f respectively, e.g., f(1 + 5% × Randn), where Randn represents the random Gaussian distribution with zero mean and unit variance. The reconstructed results with both G-BLT and V-BLT are plotted in Fig. 4
Fig. 4 Single-inclusion reconstructions with different Gaussian noise levels. The phantom is as the same as (a) in Fig. 2 with single circular inclusion; (a)-(d) are from the voxel-based BLT with 5%, 10%, 20% and 30% noise level respectively; (e)-(h) are correspondingly from Gaussian-based BLT. Please see Table 2 and 4 for the true and recovered parameters. The results show that Gaussian-based BLT is more robust to the noise than the voxel-based BLT.
and the parameter values are in Table 4

Table 4. The reconstructed parameters for single-circular-object simulations with different Gaussian noise levels. Please see Table 2 for the corresponding true values

table-icon
View This Table
| View All Tables
. Here the simulations are performed on the same circular phantom with single circular inclusion as in Fig. 2(a).

Again, from Fig. 4, visually G-BLT localizes much better than V-BLT. In particular, the result from G-BLT changes a little for the noise up to 20%. On the other hand, it is difficult to tell even the center of object in V-BLT when the noise is increased beyond 10%. Moreover, from Table 4, G-BLT gives more accurate quantitative information than V-BLT.

3.1.4. Mismatch of Optical Background

We evaluate the algorithm with different mismatch in optical background. The background is perturbed with ± 10%, ± 30%, ± 50%, ± 70% percentage of error respectively, e.g., μa(1 + 10%) and μs'(1 + 10%). Notice that we either increase or decrease both absorption and scattering coefficients since empirically we find that the mismatch may cancel if one is increased while the other is decreased. Therefore, we change both in the same direction to simulate the worst-case errors. Please also note that the reconstruction results with V-BLT are not presented since it does not localize the inclusion satisfactorily even in the case without optical background error.

The results from G-BLT are presented in Fig. 5
Fig. 5 Single-inclusion reconstructions with different mismatch in optical background. The phantom is as the same as (a) in Fig. 2 with single circular inclusion; (a)-(d) are from Gaussian-based BLT with −10%, −30%, −50%, −70% error in optical background; (e)-(h) are from Gaussian-based BLT with + 10%, + 30%, + 50%, + 70% error in optical background. Please see Table 2 and 5 for the true and recovered parameters. Please also note that the reconstruction results with the voxel-based BLT are not presented since it does not localize the inclusion satisfactorily even in the case without optical background error. The results show that Gaussian-based BLT is robust to the mismatch error of optical background in the sense that the inclusion can still be localized well despite of the shifting of centers due to the systematic discrepancy between the true values and the values used in the reconstruction.
and Table 5

Table 5. Gaussian-based BLT for single-circular-object simulations with different mismatch in optical background. Please see Table 2 for the corresponding true values

table-icon
View This Table
| View All Tables
, which show that G-BLT is robust to the mismatch error of optical background in the sense that the inclusion can still be localized well despite of the shifting of centers due to the systematic discrepancy between the true values and the values used in the reconstruction.

3.2. Reconstruction with Multiple Inclusions

Now we evaluate the algorithm in the presence of multiple inclusions by considering first Gaussian inclusions and then non-Gaussian shapes.

Despite of the fact that there are two inclusions in the test on Gaussian sources (Fig. 6
Fig. 6 Multiple-inclusion reconstructions with 5% Gaussian noise level. (a) is the phantom with two elliptical inclusions.; (b) is from the voxel-based BLT; (c) is from Gaussian-based BLT. Please see Table 6 for the true and recovered parameters. The results show again that Gaussian-based BLT not only localizes the inclusion better than, but also quantitatively more accurate than the voxel-based BLT.
) and three inclusions in the test on non-Gaussian sources (Fig. 6), we assume four inclusions as the initial guess in G-BLT, since the exact number of inclusions is usually not known in practice unless the correct anatomical priors are available from other imaging modalities. However, since the bioluminescent source is localized and limited in quantity as discussed previously, we only need to set the number of sources n to a reasonable larger value than the actual number, in order to faithfully reconstruct the source distribution without considering n as an independent variable. Particularly, as the initial guess, we assume there are four Gaussian inclusions, i.e., n = 4, with ρ1 = ρ2 = ρ3 = ρ4 = 0.1, x1,c = 5, y1,c = 5, x2,c = 5, y2,c = −5, x3,c = −5, y3,c = −5, x4,c = −5, y4,c = 5, r1,x = r2,x = r3,x = r4,x = 1, r1,y = r2,y = r3,y = r4,y = 1 and q1 = q2 = q3 = q4 = 90°. Notice that for the same reason mentioned above, we do not consider q as reconstruction variable for non-Gaussian sources for simplicity.

The same constants as in the single-inclusion simulations are used. In order to separate the multiple objects, we find out the minimal-separation constraints are usually necessary. This is also intuitively correct since the source representation may otherwise be non-unique, given the severely illposed nature of BLT. Here we set cx = cy = ln2, corresponding to that the minimal distance between any two objects is at least the average of FWHM.

3.2.1. Inclusions with Gaussian Shapes

The detailed setting is given in Table 6

Table 6. Multiple-Gaussian-inclusion simulation with Inclusion 1 [the left ellipse in Fig. 6(a)] and Inclusion 2 [the right ellipse in Fig. 6(a)]. The total energy E of the phantom is 12.570; the recovered total energy from G-BLT is 12.596 while that from V-BLT is 11.582

table-icon
View This Table
| View All Tables
and Fig. 6(a). The results are plotted in Fig. 6, which clearly shows that G-BLT is able to accurately recover the values of intensity, centers, radiuses and rotation angle. On the other hand, neither shapes nor locations of inclusions can be not be captured from V-BLT.

3.2.2. Inclusions with Non-Gaussian Shapes

The detailed setting is given in Table 7

Table 7. Multiple-inclusion simulations. Three circular inclusions of 1mm diameter are located at (0, 0), (4, 0) and (8, 0) respectively with the unit source intensity. The total energy E is 9.364; the recovered total energy from G-BLT is 9.290 while that from V-BLT is 8.282. Please note that the presented values are from G-BLT with minimal-separation constraints as in Eq. (7.3). Please also notice that three inclusions are assumed in the initial guess for the reconstruction with G-BLT although there are actually three inclusions. After the reconstruction, 3 out of 4 inclusions have distinguished peak values

table-icon
View This Table
| View All Tables
and Fig. 7(a)
Fig. 7 Multiple-inclusion reconstructions with 5% Gaussian noise level. (a) is the phantom with three circular inclusions.; (b) is from the voxel-based BLT; (c) is from Gaussian-based BLT without minimal-separation constraint; (d) is from Gaussian-based BLT with the minimal-separation constraint. Please see Table 7 for the true and recovered parameters. The results show that Gaussian-based BLT with proper geometric constraints is in general able to separate multiple inclusions better than the voxel-based BLT.
. The results are plotted in Fig. 7, which clearly shows that G-BLT with proper constraints is able to separate the inclusions while V-BLT fails. The parameter values are documented in Table 7, which are quite consistent with the true values. One minor issue is that the recovered center is not as accurate as one from the single-inclusion cases, which is understandable due to mutual effect of inclusions and the illposedness of the problem.

3.2.3. Combinatorial optimization

In this section, we give one example of combinatorial optimization (Algorithm 1) with the number of Gaussians n as a variable.

In the phantom [Fig. 8(a)
Fig. 8 Combinatorial optimization with the varying n for G-BLT. (a) is the phantom with three Gaussian inclusions.; (b)-(f) are the reconstructed images from n = 1, 2, 3, 4, 5 respectively.
], there are three Gaussians with the parameters specified in Table 8

Table 8. Combinatorial optimization with the varying n for G-BLT. Three Gaussian inclusions are located at (0, 5), (4.33, −2.5) and (−4.33, −2.5) respectively with ρ = 1, rx = 1 and ry = 2. The number of Gaussians n is treated as a variable in this combinatorial problem. The value of n in the reconstruction is set to 1 as the initial guess and iteratively increased until the solution converges. The following table records the reconstructed values for n = 1, 2, 3, 4 and 5. Please notice that the reconstruction should have terminated at n = 3 according to the given criterion if it were not for illustration purpose. The total energy E is 18.843; the recovered total energy is 21.164, 20.861, 18.879, 18.908 and 18.912 respectively. The measurement difference d is 0.167, 0.132, 0.00315, 0.00322 and 0.00327 respectively

table-icon
View This Table
| View All Tables
. In the reconstruction, we start with the initial guess n0 = 1, and continue the combinatorial process until the data discrepancy d reaches the tolerance. Please notice that the process should have terminated at n = 3. For illustration purpose, we run until n = 5.

The detailed setting and the recovered parameters are given in Table 8. Figure 8 shows that this combinatorial process is able to reconstruct the source even the initial guess is smaller than the actual number, e.g., the results from n = 1 and n = 2 fail to give a satisfactory reconstruction while those from n = 3, n = 4 and n = 5 all successfully recover the desired distribution.

It is interesting to notice that although the reconstruction with more than or equal to the actual number of Gaussians is able to recover the actual number of Gaussians, the reconstruction with more Gaussians than the actual number may degrade the image quality due to the apparent non-unique representation of the geometry. For example, two recovered sources in Fig. 8(e) constitute the top Gaussian in Fig. 8(a) when using n = 4 and two recovered sources in Fig. 8(f) constitute the left Gaussian in Fig. 8(a) when using n = 5.

3.3. 3D in vivo validation

After evaluating the propose G-BLT in 2D with various settings, we have established the superiority of G-BLT over V-BLT. Now we validate G-BLT in 3D with in vivo experimental data. The details of experimental setups and mouse experiments are given in [4

4. 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). [PubMed]

]. We also refer the readers to the above reference for the reconstruction parameters, such as the used mesh and the values for optical parameters.

Now we use G-BLT, Gaussian-based BLT, to reconstruct instead of V-BLT. Please note that we do not select the permissible region.

As the initial guess, we assume there are four Gaussian inclusions, i.e., n = 4, with ρ1 = ρ2 = ρ3 = ρ4 = 10, x1,c = x2,c = x3,c = x4,c = 0, y1,c = y2,c = y3,c = y4,c = 0, z1,c = 5, z2,c = 10, z3,c = 15, z4,c = 20, r1,x = r2,x = r3,x = r4,x = 1, r1,y = r2,y = r3,y = r4,y = 1 and r1,z = r2,z = r3,z = r4,z = 1.

The constants in the geometric constraints are: ρmax = 1000, ρmin = 0.1, xmax = ymax = 12, xmin = ymin = −12, zmax = 27, zmin = 0, rx,max = ry,max = 6, rx,min = ry,min = 0.5, and cxy = cyz = czx = 5, cx = cy = cz = ln2. Here the min-max values on the geometry correspond to the physical size of the domain; the maximum of the intensity is estimated from the data, which is roughly 50 times of the maximal value in the data.

As documented in Table 9

Table 9. In vivo study using Gaussian-based BLT. Please also notice that four inclusions are assumed in the reconstruction although there are actually two tumors from the histological verification. The reconstructed result also shows exactly two inclusions with the distinguished values, i.e. Inclusion 1 and 2 in the following table

table-icon
View This Table
| View All Tables
, 2 out of 4 recovered inclusions, with the locations just corresponding to the kidneys, have significantly larger values, i.e., the one with the peak intensity of 16.593nW/mm3 (right) and the other one with the peak intensity of 7.692nW/mm3 (left). In addition, the peak intensity of the other 2 inclusions (Inclusion 3 and 4) is less than 5% of the maximum (Inclusion 1), which can be from experimental or numerical noise. The reconstructed bioluminescent volumes (2rx × 2ry × 2ry) are 506 mm3 (right) and 901 mm3 (left) respectively. The discrepancy between the reconstructed volumes and the measured ones is not surprising since the anatomical volume may not correspond to the bioluminescent volume.

4. Conclusions and discussions

In this proof-of-concept study, we have introduced a stable way to recover bioluminescent source through novel source representation by Gaussians. It is mathematically unique and computationally cheap. Besides, there are strong physiological evidences that the actual bioluminescent sources can be modeled by Gaussians, which will be further studied in the future work. The resultant problem with the fixed number of Gaussians is formulated as a nonlinear least-square minimization with appropriate geometric constraints, which is solved in turn by barrier method; on the other hand, the problem with the varying number of Gaussians can be formulated as a combinatorial optimization problem. The superiority of the proposed method has been established through simulations and in vivo experimental data. In particular, when the source itself can be approximated by Gaussians, the proposed method is able to accurately recover the intensity, centers, radiuses and rotation angles of the Gaussians.

Due to the significantly reduced number of variables, the proposed method is robust to the measurement noise and the mismatch in optical background; the computational burden is also reduced to the minimization problem with respect to a few parameters instead of at least thousands of parameters in the conventional voxel-based BLT.

In the future work, we shall perform more mouse studies to further develop and optimize the current method into a practical in vivo bioluminescence imaging method, for which we may need to incorporate RTE for more accurate modeling, simplified representation for optical heterogeneity, such as piecewise constants, to be simultaneously reconstructed for correcting optical background, and multi-spectral data for better stability.

Although the source is represented as Gaussians in this study, other shape functions may be appropriate for other purposes. An interesting topic to pursuit is the multi-modality imaging since the structural priors can be easily incorporated into the proposed method for more accurate quantitative estimation.

Acknowledgments

The work is partially supported by National Science Foundation (NSF) grant DMS0811254, National Institutes of Health (NIH)/National Heart, Lung and Blood Institute (NHLBI) grant HL098912 and NIH/National Cancer Institute (NCI) grant CA135151.

References and links

1.

C. H. Contag and B. D. Ross, “It’s not just about anatomy: in vivo bioluminescence imaging as an eyepiece into biology,” J. Magn. Reson. Imaging 16(4), 378–387 (2002). [PubMed]

2.

G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. Meinel, “Development of the first bioluminescent CT scanner,” Radiology 299, 566 (2003).

3.

G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. 31(8), 2289–2299 (2004). [PubMed]

4.

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). [PubMed]

5.

W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. Wang, E. Hoffman, G. McLennan, P. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express 13(18), 6756–6771 (2005). [PubMed]

6.

G. Wang, X. Qian, W. Cong, H. Shen, Y. Li, W. Han, K. Durairaj, M. Jiang, T. Zhou, and J. Cheng, “Recent development in bioluminescence tomography,” Curr. Med. Imaging Rev. 2, 453–457 (2006).

7.

G. Wang, H. Shen, K. Durairaj, X. Qian, and W. Cong, “The first bioluminescence tomography system for simultaneous acquisition of multi-view and multi-spectral data,” Int. J. Biomed. Imaging 2006, 1–8 (2006).

8.

X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express 12(17), 3996–4000 (2004). [PubMed]

9.

A. J. Chaudhari, F. Darvas, J. R. Bading, R. A. Moats, P. S. Conti, D. J. Smith, S. R. Cherry, and R. M. Leahy, “Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging,” Phys. Med. Biol. 50(23), 5421–5441 (2005). [PubMed]

10.

H. Dehghani, S. C. Davis, S. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally resolved bioluminescence optical tomography,” Opt. Lett. 31(3), 365–367 (2006). [PubMed]

11.

C. Kuo, O. Coquoz, T. Troy, D. Zwarg, and B. Rice, “Bioluminescent tomography for in vivo localization and quantification of luminescent sources from a multiple-view imaging system,” Mol. Imaging 4, 370 (2005).

12.

G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(17), 4225–4241 (2005). [PubMed]

13.

Y. Lv, J. Tian, W. Cong, and G. Wang, “Experimental study on bioluminescence tomography with multimodality fusion,” Int. J. Biomed. Imaging 2007, 86741 (2007). [PubMed]

14.

H. Gao and H. K. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization,” Opt. Express 18(3), 1854–1871 (2010). [PubMed]

15.

H. Gao and H. K. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation part 2: total variation and l1 data fidelity,” Opt. Express 18(3), 2894–2912 (2010). [PubMed]

16.

J. Liu, A. Li, A. E. Cerussi, and B. J. Tromberg, “Parametric diffuse optical imaging in reflectance geometry,” IEEE Sel. Top. Quantum Electron. 16, 555–564 (2010).

17.

K. M. Case, and P. F. P. F. Zweifel, Linear Transport Theory (Addison-Wesley Educational Publishers Inc., 1967).

18.

A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic Press, 1978).

19.

E. E. Lewis, and W. F. Miller, Computational Methods of Neutron Transport (Wiley, 1984).

20.

H. Gao and H. K. Zhao, “A fast forward solver of radiative transfer equation,” Transp. Theory Stat. Phys. 38, 149–192 (2009).

21.

H. Shen and G. Wang, “A tetrahedron-based inhomogeneous Monte Carlo optical simulator,” Phys. Med. Biol. 55(4), 947–962 (2010). [PubMed]

22.

W. Cong, H. Shen, A. Cong, Y. Wang, and G. Wang, “Modeling photon propagation in biological tissues using a generalized Delta-Eddington phase function,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 76(5 Pt 1), 051913 (2007). [PubMed]

23.

A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys. 220, 441–470 (2006).

24.

A. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).

25.

S. Boyd, and L. Vandenberghe, Convex Optimization (Cambridge university press, 2004).

OCIS Codes
(100.3190) Image processing : Inverse problems
(110.6960) Imaging systems : Tomography
(170.3010) Medical optics and biotechnology : Image reconstruction techniques
(170.6280) Medical optics and biotechnology : Spectroscopy, fluorescence and luminescence

ToC Category:
Image Reconstruction and Inverse Problems

History
Original Manuscript: September 20, 2010
Revised Manuscript: October 16, 2010
Manuscript Accepted: October 27, 2010
Published: October 29, 2010

Citation
Hao Gao, Hongkai Zhao, Wenxiang Cong, and Ge Wang, "Bioluminescence tomography with Gaussian prior," Biomed. Opt. Express 1, 1259-1277 (2010)
http://www.opticsinfobase.org/boe/abstract.cfm?URI=boe-1-5-1259


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. C. H. Contag and B. D. Ross, “It’s not just about anatomy: in vivo bioluminescence imaging as an eyepiece into biology,” J. Magn. Reson. Imaging 16(4), 378–387 (2002). [PubMed]
  2. G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. Meinel, “Development of the first bioluminescent CT scanner,” Radiology 299, 566 (2003).
  3. G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. 31(8), 2289–2299 (2004). [PubMed]
  4. 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). [PubMed]
  5. W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. Wang, E. Hoffman, G. McLennan, P. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express 13(18), 6756–6771 (2005). [PubMed]
  6. G. Wang, X. Qian, W. Cong, H. Shen, Y. Li, W. Han, K. Durairaj, M. Jiang, T. Zhou, and J. Cheng, “Recent development in bioluminescence tomography,” Curr. Med. Imaging Rev. 2, 453–457 (2006).
  7. G. Wang, H. Shen, K. Durairaj, X. Qian, and W. Cong, “The first bioluminescence tomography system for simultaneous acquisition of multi-view and multi-spectral data,” Int. J. Biomed. Imaging 2006, 1–8 (2006).
  8. X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express 12(17), 3996–4000 (2004). [PubMed]
  9. A. J. Chaudhari, F. Darvas, J. R. Bading, R. A. Moats, P. S. Conti, D. J. Smith, S. R. Cherry, and R. M. Leahy, “Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging,” Phys. Med. Biol. 50(23), 5421–5441 (2005). [PubMed]
  10. H. Dehghani, S. C. Davis, S. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally resolved bioluminescence optical tomography,” Opt. Lett. 31(3), 365–367 (2006). [PubMed]
  11. C. Kuo, O. Coquoz, T. Troy, D. Zwarg, and B. Rice, “Bioluminescent tomography for in vivo localization and quantification of luminescent sources from a multiple-view imaging system,” Mol. Imaging 4, 370 (2005).
  12. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(17), 4225–4241 (2005). [PubMed]
  13. Y. Lv, J. Tian, W. Cong, and G. Wang, “Experimental study on bioluminescence tomography with multimodality fusion,” Int. J. Biomed. Imaging 2007, 86741 (2007). [PubMed]
  14. H. Gao and H. K. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization,” Opt. Express 18(3), 1854–1871 (2010). [PubMed]
  15. H. Gao and H. K. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation part 2: total variation and l1 data fidelity,” Opt. Express 18(3), 2894–2912 (2010). [PubMed]
  16. J. Liu, A. Li, A. E. Cerussi, and B. J. Tromberg, “Parametric diffuse optical imaging in reflectance geometry,” IEEE Sel. Top. Quantum Electron. 16, 555–564 (2010).
  17. K. M. Case, and P. F. P. F. Zweifel, Linear Transport Theory (Addison-Wesley Educational Publishers Inc., 1967).
  18. A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic Press, 1978).
  19. E. E. Lewis, and W. F. Miller, Computational Methods of Neutron Transport (Wiley, 1984).
  20. H. Gao and H. K. Zhao, “A fast forward solver of radiative transfer equation,” Transp. Theory Stat. Phys. 38, 149–192 (2009).
  21. H. Shen and G. Wang, “A tetrahedron-based inhomogeneous Monte Carlo optical simulator,” Phys. Med. Biol. 55(4), 947–962 (2010). [PubMed]
  22. W. Cong, H. Shen, A. Cong, Y. Wang, and G. Wang, “Modeling photon propagation in biological tissues using a generalized Delta-Eddington phase function,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 76(5 Pt 1), 051913 (2007). [PubMed]
  23. A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys. 220, 441–470 (2006).
  24. A. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999).
  25. S. Boyd, and L. Vandenberghe, Convex Optimization (Cambridge university press, 2004).

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