Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Proc IEEE Int Symp Biomed Imaging. Author manuscript; available in PMC 2010 April 27.
Published in final edited form as:
Proc IEEE Int Symp Biomed Imaging. 2009 June 28 : 446–449.
doi:  10.1109/ISBI.2009.5193080
PMCID: PMC2860303



Recently, we generalized the Delta-Eddington phase function and applied it to the radiative transfer equation for modeling the photon propagation in biological tissue. The resultant phase approximation model was shown to be highly accurate with a wide range of optical properties, including the strongly absorbing and weakly scattering media. In this paper, we propose phase-approximation-based method for estimating the optical parameters. Specifically, we design an iterative algorithm to take advantage of both the global search ability of the differential evolution algorithm and the efficiency of the conjugate gradient method. Then, we demonstrate the feasibility and merits of the proposed method in both numerical simulation and phantom experiments.

Index Terms: Phase Approximation (PA), Bioluminescence Tomography (BLT)


Estimation of optical properties of biological tissue is an essential task for optical imaging applications [1, 2]. Over past years, a number of optical property measurement methods were proposed based on different solutions to the radiative transport equation (RTE) [3, 4, 5, 6]. Recently we generalized the Delta-Eddington phase function and applied it to the radiative transfer equation for modeling the photon propagation in biological tissue, which is referred to as the phase approximation (PA) model for brevity [7, 8, 9]. It is well known that the radiative transfer equation (RTE) can accurately model the photon propagation in turbid media. Our PA model approximates the phase function in RTE using the generalized Delta-Eddington phase function instead of the commonly used Henyay-Greenstein function, and essentially keeps the accuracy of RTE while maintaining a high efficiency. In this manuscript, we propsed a method for estimating the optical parameters for the PA model by matching the calculated and measured photon fluence rates iteratively. The iterative procedures based on gradient-guided optimization are sensitive to initial guesses. By combining the global search ability of the stochastic algorithm and the efficiency of the gradient-guided optimization method, we are able to obtain an accurate, efficient and robust parameter reconstruction.


The radiative transfer equation (RTE) is considered as the most accurate for modeling photon propagation in turbid media and has been adopted for tissue optics [10]. RTE is given as:

s·L(r,s)+(μa+μs)L(r,s)=μsS2L(r,s)p(s,s)ds,   rΩ

where μa and μs are the absorption and scattering coefficients respectively, r is the position vector in the domain Ω, L the photon radiance, and p (s′, s) the scattering phase function from direction s′ to s. Since the RTE is difficult to solve efficiently, an approximation is often used to achieve adequate computational efficiency, such as the diffusion approximation (DA). Our approach is to simplify the RTE by approximating the scattering phase function as a generalized Delta-Eddington function. Actually, the exact form of the phase function p (s′, s) is unknown in practice, and the Henyey-Greenstein (HG) phase function is a good approximation and is almost exclusively used in tissue optics. However, up to now RTE with the HG phase function could not be further reduced. Considering the strongly forward component of the photon scattering in biological tissue, we adopted the simplified Delta-Eddington phase function:


where υ = s′·s, and f is the anisotropy weight factor, which is a unique optical parameter to the PA model. With the Delta-Eddington phase function, we have the following phase approximation equation [7, 9]:


where [mu]s = (1 − fs is the reduced scattering coefficient, and [mu]t = μa + [mu]s the reduced attenuation coefficient. G(r,r)=1|rr|2exp(0|rr|μ˜t(r)(rtβ)dt),andβ=rr|rr|, S (r′) is a collimated illumination source on the boundary, and rd is the internal reflection coefficient that is approximated by rd = −1.4399η−2+0.7099η−1+0.6681+0.0636η, η = ntissue/nair.


To solve the photon fluence rate from Eq.(3), we discretize the domain Ω into tetrahedral finite elements. Then, the photon fluecence can be approximated as:


where [var phi] (ri) are the basis functions. Using Eq.(4), the PA model can be converted to a system of linear equations:


where {Φ} is the discretized photon flucence on the nodes of the finite element mesh. The components of M, B and S are defined as [7]:






The integrals are evaluated using the Gauss-Legendre quadrature rule. Singular points in the quadrature need to be handled to ensure the numerical accuracy of the solution. The singularities occur at ri = r′, which are the vertices of each tetrahedral element. We denote the vertices of a tetrahedron as left angle bracket(x1, y1, z1), (x2, y2, z2), (x3, y3, z3), (x4, y4, z4)right angle bracket. Without losing generality, let us only show the case ri = (x4, y4, z4). We integrate over an arbitrary tetrahedron by performing an integral over the unit orthogonal tetrahedron in volume coordinate (ξ, η, ζ) after the transformation

x(ξ,η,ζ)  =  x4+ξ(x1x4)+η(x2x4)+ζ(x3x4)y(ξ,η,ζ)  =  y4+ξ(y1y4)+η(y2y4)+ζ(y3y4)z(ξ,η,ζ)  =  z4+ξ(z1z4)+η(z2z4)+ζ(z3z4)

where the singular point (x4, y4, z4) is mapped to the origin of the volume coordinate system. To avoid the singular point, the integral is expressed in a spherical coordinate system:

mi,j=μ˜s|det J|4πV¯ϕjλ exp (μ˜tρλ)sinϕdρdθdϕ=μ˜s|det J|4π0π20π201τϕjλ exp (μ˜tρλ)sinϕdρdθdϕ



and λ =



With correct optical parameters, the calculated photon fluence rate from the PA model (Eq.(3)) will be a good approximation to the measurement. The reconstruction process is to minimize the objective function Ψ (μa, [mu]s) by updating the optical parameters iteratively. The minimization problem is expressed as:


where {Φm} is the measured photon fluence rate. We use the adjoint differentiation to effectively calculate the gradient of the objective function. The gradient of Ψ can be expressed as


We adopt the conjugate gradient (CG) optimizer to minimize the objective function using its gradient starting from an initial guess of optical parameters. However, the CG minimization procedure is sensitive to the initial optical parameters, and a poor selection of the initial guess will lead to a local minimum instead of the true solution [11]. To improve the robustness of the optimization, the CG method can be preceded by a global search procedure. The global search is used to provide the CG optimizer with a good initial solution which is close to the true solution. Thus, the the CG method will much more likely converge to the true solution. Differential evolution (DE) is a population-based stochastic method for global optimization [12]. In DE, the collection of solution vectors is called the population which is denoted as


where NP is the population size, and g the number of generations. The solution vector in this case consists of the optical parameters: x = (μa, [mu]s). The population is mutated using the unique differential operator [12]:


where F is the difference scale factor, which is the most important control parameter of DE. The base vector xr0,g and the difference vectors xr1,g and xr2,g are randomly selected but r0, r1, and r2 should be mutually exclusive. The mutated vector vi,g then exchanges information with a target vector yi,g which is selected sequentially from the population to generate a trial vector ui,g : each component in the trial vector inherits from that of the mutated vector if a binomial trial is a success, i.e. rand[0, 1) ≤ Cr, where Cr is the crossover rate; otherwise, that component will inherit from the target vector. If the trial vector is a better solution, i.e., it has a smaller objective function value, then the trial vector replaces the target vector in the population. Finally, the generation counter is advanced when the entire population is updated, and the solution vector with the smallest function value is used as the solution when the maximum generation is reached. The solution generate by the DE optimization in a small number of iterations is used as the initial starting point for the CG optimizer. Using this hybridization scheme, we can enhance the likelihood of finding the true solution when the a priori knowledge on the optical properties is uncertain.


5.1. Numerical Study

In this section, we demonstrate the effectiveness of our method and characterize the accuracy of our algorithm using simulated measurements. The simulations are performed on a digital phantom of 10mm × 10mm × 10mm cube with four different optical property configurations. An illumination source was set on one side of the cube and emitted a pencil beam of 1μW into the phantom. Monte Carlo (MC) simulation was used to generate the measurements on the opposite side of the source. Because of the intrinsic Poisson noise in MC simulated, no additional noise is added. The simulated measurements are presented in Fig. 1.

Fig. 1
Simulated measurements of four optical property configurations.

To apply the reconstruction algorithm, we generated a finite element model for the slab. The finite element model was composed of 1708 nodes and 8665 tetrahedron with average edge length of 1.0mm. We used a population size of 24, differential scale factor of 0.8, and crossover rate of 0.9 for the differential evolution optimization and allow maximum 20 generations and selected the single best candidate solution from the population as the starting point of the conjugate gradient optimizer. The reconstruction results are listed in Table 1. The relative error for the reconstruction of μa is less then 10% and the error for [mu]s is less then 6%.

Table 1
The reconstructed optical parameters.

5.2. Phantom Study

Using the optical property reconstruction algorithm, we estimated the optical parameters of a mouse-shaped phantom (Xenogen Corporation) for different light spectra. The phantom was made by a homogeneous material and a illumination light source (Xenogen XPM-2) was embedded into it. The light source can be treated as a point source and has a depth of 4.8mm to the bottom surface of the mouse phantom. Using two filters of 515–575nm and 575–650nm, we obtain two sets of measurement of the photons escaping from the flat surface of the phantom using Xenogen’s IVIS 100 series imaging system, respectively.

A finite element model was established from the microCT image volume of the mouse-shaped phantom. The finite element model had 3864 nodes and 14708 tetrahedra. The measurements were mapped to the corresponding surface of the finite element model, as in Fig. 2.

Fig. 2
The measurements mapped to the finite element model of the phantom.

Similar to the previous phantom experiments, we used a population size of 24, differential scale of 0.8, crossover rate of 0.5, and maximum generation of 20 for the differential evolution optimization before the CG optimizer. The reconstructed optical properties are listed in Table 2.

Table 2
The reconstructed optical parameters for different illumination spectra.


The generalized Delta-Eddington phase approximation greatly simplified RTE, making it more efficient to solve while maintaining a accuracy comparable to that of the original RTE. The PA model has an unique reduced scattering coefficient [mu]s that must be determined in order to be practically used in optical imaging applications. Optical property estimation that is based on the PA model will inherently yield accurate optical parameters for a wide range of tissue properties. The accuracy of the optical parameters characterization, however, does not solely depends on the accuracy of the forward model; the iterative algorithm could converge prematurely to an inaccurate solution if the initial guess is not close enough to the true solution. Fortunately, differential evolution is an effective global search method to find a good starting point quickly. In this paper, we focus on the optical property characterization of a homogeneous biological tissue under the PA model. The non-homogeneous counterpart will be developed in the future to apply the PA model in optical tomography for small animal and clinical imaging.

In conclusion, we have presented a model-based, iterative method for determining the optical properties of biological tissue based on the phase approximation model. We developed an robust iterative procedure to eliminate the need of initial guess of the optical parameters. The proposed optical property characterization method, along with the PA model will improve the accuracy for optical molecular imaging modalities.


*The authors are grateful for the support from the following grants : NIH/NIBIB EB001685, NIH/NIBIB EB006036, NIH/NIBIB EB008476, NIH/NIBIB EB001458, NIH/NCI 2U24 CA092865, and DE-FC02-02ER63520.

Contributor Information

A Cong, VT-WFU School of Biomedical Eng & Sci, Virginia Tech.

W Cong, VT-WFU School of Biomedical Eng & Sci, Virginia Tech.

H Shen, VT-WFU School of Biomedical Eng & Sci, Virginia Tech.

G Wang, VT-WFU School of Biomedical Eng & Sci, Virginia Tech.

Y Lu, David Geffen School of Medicine at UCLA University of California.

A Chatziioannou, David Geffen School of Medicine at UCLA University of California.


1. Wang G, Qian X, Cong W, Shen H, Li Y, Han W, Durairaj K, Jiang M, Zhou T, Cheng J, et al. Recent Development in Bioluminescence Tomography. Current Medical Imaging Reviews. 2006;2(4):453–457.
2. Weissleder R, Ntziachristos V. SCALING DOWN IMAGING: MOLECULAR MAPPING OF CANCER IN MICE. Nature Medicine. 2003;9:123–128.
3. Prahl S, Van Gemert M, Welch A. Determining the optical properties of turbid media by using the adding-doubling method. APPLIED OPTICS. 1993;32(4/1) [PubMed]
4. Roggan A, Dorschel K, Minet O, Wolff D, Muller G. The optical properties of biological tissue in the near infrared wavelength rangereview and measurements. Laser-Induced Interstitial Thermotherapy. 1995:10–44.
5. Hayakawa C, Spanier J, Bevilacqua F, Dunn A, You J, Tromberg B, Venugopalan V. Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues. Optics Letters. 2001;26(17):1335–1337. [PubMed]
6. Hielscher A, Klose A, Hanson K. Gradient-based iterative image reconstruction scheme for time-resolved optical tomography. Medical Imaging, IEEE Transactions on. 1999;18(3):262–271. [PubMed]
7. Cong W, Shen H, Cong A, Wang Y, Wang G. Modeling photon propagation in biological tissues using a generalized Delta-Eddington phase function. Physical Review E. 2007;76(5):51,913. [PubMed]
8. Cong W, Cong A, Shen H, Liu Y, Wang G. Flux vector formulation for photon propagation in the biological tissue. Optics Letters. 2007;32(19):2837–2839. [PubMed]
9. Cong W, Shen H, Cong A, Wang G. Integral equations of the photon fluence rate and flux based on a generalized Delta-Eddington phase function. Journal of Biomedical Optics. 2008;13:024,016. [PMC free article] [PubMed]
10. Welch A, van Gemert M. Optical-thermal Response of Laser-irradiated Tissue. Plenum Pub Corp. 1995
11. Hielscher A, Klose A, Beuthan J. Evolution strategies for optical tomographic characterization of homogeneous media. Optics Express. 2000;7(13):507–518. [PubMed]
12. Price K, Storn R, Lampinen J. Differential Evolution: A Practical Approach To Global Optimization. Springer. 2005