DETERMINATION OF MATERIAL PROPERTIES: OPTICAL TOMOGRAPHY APPLICATIONS
1. INTRODUCTION
During the last decade, increasing developments of new optical techniques of clinical controls and medical diagnostics have been achieved. These developments were made possible thanks to advances in the theoretical understanding of the interaction between light and optical properties of semitransparent media such as human tissues. On the other hand, the improvement of light-detection systems as well as light sources that irradiate the tissues was necessary in order to put these techniques into practice. One of the most spectacular of these advances is the possibility to acquire a picture of the spatial distribution of the optical properties of different parts of the body from measurements at the surface. This technology of imaging is known as optical tomography (Arridge, 1999). It consists of identifying the parameters of a numerical model based on light transport through tissues from boundary measurements by using a minimization of the mean square discrepancy between the measurements and the predictions of the numerical model at detectors locations. In the following, an algorithm is presented together with some results of retrieval of optical properties of human-like media.
2. FORWARD MODEL
Light transport in biological tissues is described by a Boltzmann-type integro-differential equation called the radiative transfer equation. This equation is difficult to solve, and analytical solutions are available only for simple cases. Different approximations are used, among which the diffusion approximation is well suited for biological tissues due to their high scattering (Modest, 1993). While this approximation leads to an equation that is easy to solve, it fails to describe light transport near boundaries, sources, and in void-like regions where the mean free path is very large (Hielscher et al., 1998). Nowadays, major improvements are carried out in order to obtain rapid and accurate algorithms. Among them, the use of the complete radiative transfer equation (RTE) as a forward model to overcome the limitations of the diffuse approximation in optical tomography has recently received great interest (Klose et al., 2002; Charette et al., 2007). The RTE is used in different forms (time, steady, and frequency domains):
(1) |
where n is the refractive index of the medium, c_{0} is the light speed, κ and σ are, respectively, the absorption and scattering coefficients, I is the light intensity at the spatial position, , in the direction , t is the time and ω is the modulation frequency. S(, , t) represents the emission and all additional source terms. Φ(, ) is the scattering phase function. Generally, the scattering in tissues is given by the Henyey-Greenstein phase function (Modest, 1993)
(2) |
where Θ is the scattering angle between directions and , and g is the anisotropic factor. Note that in Eq. (1) the wavelength dependency is not considered as it is generally assumed that light scattering does not affect the wavelength (scattering coherence). For optical tomography applications, the measurable quantities used are
(3) |
where ∂D is the boundary of the domain D being considered, and is the outward normal unit vector. In the following, the development will be made with the frequency domain formulation where a single modulation frequency is used. However, the ideas can be applied to steady and time domains as well with some modifications.
3. INVERSION
In the inversion procedure, the aim is to recover the optical properties of the medium through the minimization of an objective function. We present here a gradient type optimization where the objective function gradient is computed through an adjoint formulation
3.1. Objective Function
The objective function J to be minimized in the inversion procedure is the mean square discrepancy between the measurements and the predictions of the forward model such that
(4) |
where θ is the vector of parameters [θ = (κ,σ), i.e., the absorption and scattering coefficients], P(r_{d}, ω) and M_{d} are, respectively, the dth detector prediction and measurement, and r_{d} is the position of the dth detector. N_{d} is the number of detectors, ||.|| is an l_{2} norm on the solution space. We suppose that κ and σ are constant per mesh element, i.e., θ = (κ, σ) ∈ R^{2Ne} where N_{e} is the number of mesh elements. Then, the minimization problem writes as
(5) |
where R(θ,I) is the residual form of the forward problem, which is written in the frequency domain as
(6) |
3.2 Gradient of the Objective Function
Gradient computation for large-scale problems is highly time consuming. The adjoint formulation gives an easy and fast way to compute the gradient by solving an additional equation for the adjoint variable and evaluating the gradient through a simple scalar product. Here, we present a Lagrangian formulation. The Lagrangian formulation of the minimization problem given by Eq. (5) takes the following form:
(7) |
where λ(, , ω) is a vector that represents the adjoint variable and [.|.] is an inner product associated to the norm used in Eq. (4). As the residual R(θ,I) is always null when I is the solution of the forward problem, then we have L(θ,I,λ) = J(θ) for any value of λ. The differentiation of the Lagrangian L(θ,I,λ) with respect to θ in the direction δθ writes as
(8) |
Taking into account that the residual equation is null [see Eq. (6)], we rewrite Eq. (8) as
(9) |
In order to compute the gradient of the objective function J(θ) without having to compute the sensitivity of I with respect to θ in the direction δθ, the adjoint variable λ is chosen such that
(10) |
This leads to the following equation for the adjoint variables (Favennec et al., 2006; Balima et al., 2007; Balima, 2011):
(11) |
where the subscript * is relative to the transposition of the residual operator. Using the expression of the residual equations, we obtain the following adjoint equation to be solved for λ(, , ω):
(12) |
Using Eqs. (8) and (10), the differentiation of the objective function J(θ) with respect to θ in the direction δθ is rewritten as
(13) |
Since the objective function does not depend explicitly on θ, we have , and therefore Eq. (13) reduces to
(14) |
which is the expression for evaluating the objective function gradient.
4. EXAMPLES IN THE FREQUENCY DOMAIN
4.1 Example Using a Finite Volume Forward Model
The test medium is a 2D domain of 2 cm × 2 cm that contains three objects in the background medium with scattering and absorption coefficients of 50 cm^{-1} and 0.5 cm^{-1}, respectively, see Fig. 1. Two of them (white) have optical properties of σ = 60 cm^{-1} and κ = 0.6 cm^{-1}, while the third (black) is assigned optical properties of σ = 40 cm^{-1} and κ = 0.4 cm^{-1}. Four zero-phased sources with a frequency of 600 MHz are placed at mid-center of the medium boundaries, and 80 detectors are used around the medium, which yields a total of 320 source-detector pairs to be used in the inversion. The synthetic data are generated on a 81 × 81 grid with 80 discrete ordinates (S_{8} quadrature) by using the forward model based on the finite volume method, and the generated complex-valued intensities are used as the input data for the inversion. The medium is assumed to be forward scattering with a Henyey-Greenstein anisotropic factor g = 0.9. Note that the reduced scattering coefficient σ' = σ(1 - g) is used in the figure. The inversion calculations are performed on a 41 × 41 grid, coarser than the grid used to obtain the input data and the Limited memory BFGS method (Liu and Nocedal, 1989) is employed. The quality is measured with the error norm:
Figure 1. (a) and (b) depict the original phantom in absorption and scattering, respectively; the plain circles indicate the location where the beam is applied; (c) and (d) represent the numerical retrieval of absorption and scattering coefficients, respectively, with the L-BFGS method; error between reconstructed and original images is 0.01523 for scattering and 0.01858 for absorption; the color shades (in cm^{-1}) are implicitly defined by figures (a) and (b) (from Charette et al., 2008).
(15) |
where superscripts r and 0 refer to the reconstruction and original images, respectively, and subscripts i and j are the parameter mesh nodes indices for the corresponding optical property θ, which can be either the absorption or the scattering coefficient.
4.2 Example Using a Finite Elements Forward Model
We have done some work to handle complex geometries for tomography application with a least square finite element formulation (Balima et al., 2010a, 2011). Some preliminary results are given in Fig. 2. All details on the accuracy of the forward model are given in Balima et al. (2010b). The test problem used here is the same as above, but with two circular inclusions, and the sources are collimated. The synthetic data are generated on a regular triangle mesh of 958 elements, with 2 × 520 spatial degrees of freedom (P1 elements) and 48 discrete ordinates (S_{6} quadrature) by using the forward model (Balima et al., 2010b). The generated complex-valued intensities are either noise free or noised, with a uniform Gaussian of -20 dB level and used as the input data for the inversion. The number of detectors is 74 for each source. The inversion calculations are performed on the same mesh with the L-BFGS method.
Figure 2. (a) and (b) are the original absorption and scattering coefficients, respectively; (c) and (d) are, respectively, the identified absorption and scattering coefficients with L-BFGS on noise-free data; error between noise-free reconstructed and original images is 0.0174867 for scattering and 0.0228672 for absorption. (e) and (f) are, respectively, the identified absorption and scattering coefficients with L-BFGS on noised data (-20 dB); error between reconstructed with noised data and original images is 0.0329933 for scattering and 0.0487816 for absorption.
REFERENCES
Arridge, S. R., Optical Tomography in Medical Imaging, Inverse Prob., vol. 15, pp. R41-R93, 1999.
Balima, O., Favennec, Y., and Petit, D., Model Reduction on Heat Conduction with Radiative Boundary Conditions Using the Modal Identification Method, Numer. Heat Transfer Part B, vol. 52, pp. 107-130, 2007.
Balima, O., Charette, A., and Marceau, D., Comparison of Light Transport Models in View of Optical Tomography Applications, J. Comput. Appl. Math., vol. 234, no. 7, pp. 2259–2271, 2010a.
Balima, O., Pierre, T., Charette, A., and Marceau, D., A Least Square Finite Element Formulation of the Collimated Irradiation in Frequency Domain for Optical Tomography Applications, J. Quant. Spectrosc. Radiat. Transfer, vol. 111, no. 2, pp. 280-286, 2010b.
Balima, O., Boulanger, J., Charette, A., and Marceau, D., New developments in frequency domain Optical tomography. Part I, Part II, J. Quant. Spectrosc. Radiat. Transfer, vol. 112, no. 7, pp. 1229–1240, 2011.
Charette, A., Boulanger, J., and Kim, H., Optical Tomography as an Inverse Radiation Problem, Fifth International Symposium on Radiation Transfer, Bodrum, Turkey, 2007.
Charette, A., Boulanger, J., and Kim, H. K., An Overview on Recent Radiation Transport Algorithm Development for Optical Tomography Imaging, J. Quant. Spectrosc. Radiat. Transfer, vol. 109, no. 17-18, pp. 2743-2766, 2008.
Favennec, Y., Girault, M., and Petit, D., The Adjoint Method Coupled with the Modal Identification Method for Nonlinear Model Reduction, Inverse Prob. Sci. Eng., vol. 14, pp. 153-170, 2006.
Hielscher, A., Alcouffe, R., and Barbour, R., Comparison of Finite Difference Transport and Diffusion Calculations for Photon Migration in Homogeneous and Heterogeneous Tissues, Phys. Med. Biol., vol. 43, pp. 1285-1302, 1998.
Klose, A., Netz, U., Beuthan, J., and Hielscher, A., Optical Tomography Using the Time-Independent Equation of Radiative Transfer--Part 1: Forward Model, J. Quant. Spectrosc. Radiat. Transfer, vol. 72, pp. 691-713, 2002.
Liu, D. and Nocedal, J., On the Limited Memory BFGS Method for Large Scale Optimization, Math. Prog., vol. 45, no. 3, pp. 503-528, 1989.
Modest, M. F., Radiative Heat Transfer, Mechanical Engineering, McGraw Hill, New York, 1993.
References
- Arridge, S. R., Optical Tomography in Medical Imaging, Inverse Prob., vol. 15, pp. R41-R93, 1999.
- Balima, O., Favennec, Y., and Petit, D., Model Reduction on Heat Conduction with Radiative Boundary Conditions Using the Modal Identification Method, Numer. Heat Transfer Part B, vol. 52, pp. 107-130, 2007.
- Balima, O., Charette, A., and Marceau, D., Comparison of Light Transport Models in View of Optical Tomography Applications, J. Comput. Appl. Math., vol. 234, no. 7, pp. 2259â€“2271, 2010a.
- Balima, O., Pierre, T., Charette, A., and Marceau, D., A Least Square Finite Element Formulation of the Collimated Irradiation in Frequency Domain for Optical Tomography Applications, J. Quant. Spectrosc. Radiat. Transfer, vol. 111, no. 2, pp. 280-286, 2010b.
- Balima, O., Boulanger, J., Charette, A., and Marceau, D., New developments in frequency domain Optical tomography. Part I, Part II, J. Quant. Spectrosc. Radiat. Transfer, vol. 112, no. 7, pp. 1229â€“1240, 2011.
- Charette, A., Boulanger, J., and Kim, H., Optical Tomography as an Inverse Radiation Problem, Fifth International Symposium on Radiation Transfer, Bodrum, Turkey, 2007.
- Charette, A., Boulanger, J., and Kim, H. K., An Overview on Recent Radiation Transport Algorithm Development for Optical Tomography Imaging, J. Quant. Spectrosc. Radiat. Transfer, vol. 109, no. 17-18, pp. 2743-2766, 2008.
- Favennec, Y., Girault, M., and Petit, D., The Adjoint Method Coupled with the Modal Identification Method for Nonlinear Model Reduction, Inverse Prob. Sci. Eng., vol. 14, pp. 153-170, 2006.
- Hielscher, A., Alcouffe, R., and Barbour, R., Comparison of Finite Difference Transport and Diffusion Calculations for Photon Migration in Homogeneous and Heterogeneous Tissues, Phys. Med. Biol., vol. 43, pp. 1285-1302, 1998.
- Klose, A., Netz, U., Beuthan, J., and Hielscher, A., Optical Tomography Using the Time-Independent Equation of Radiative Transfer--Part 1: Forward Model, J. Quant. Spectrosc. Radiat. Transfer, vol. 72, pp. 691-713, 2002.
- Liu, D. and Nocedal, J., On the Limited Memory BFGS Method for Large Scale Optimization, Math. Prog., vol. 45, no. 3, pp. 503-528, 1989.
- Modest, M. F., Radiative Heat Transfer, Mechanical Engineering, McGraw Hill, New York, 1993.