- Email: [email protected]
Gauss–Newton reconstruction method for optical tomography using the finite element solution of the radiative transfer equation
ARTICLE IN PRESS Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
Contents lists available at ScienceDirect
Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt
Gauss–Newton reconstruction method for optical tomography using the finite element solution of the radiative transfer equation T. Tarvainen a,,1, M. Vauhkonen a, S.R. Arridge b a b
Department of Physics, University of Kuopio, P.O. Box 1627, 70211 Kuopio, Finland Department of Computer Science, University College London, Gower Street, London WC1E 6BT, UK
a r t i c l e in fo
abstract
Article history: Received 10 July 2008 Received in revised form 29 August 2008 Accepted 29 August 2008
The radiative transfer equation can be utilized in optical tomography in situations in which the more commonly applied diffusion approximation is not valid. In this paper, an image reconstruction method based on a frequency domain radiative transfer equation is developed. The approach is based on a total variation output regularized least squares method which is solved with a Gauss–Newton algorithm. The radiative transfer equation is numerically solved with a finite element method in which both the spatial and angular discretizations are implemented in piecewise linear bases. Furthermore, the streamline diffusion modification is utilized to improve the numerical stability. The approach is tested with simulations. Reconstructions from different cases including domains with low-scattering regions are shown. The results show that the radiative transfer equation can be utilized in optical tomography and it can produce good quality images even in the presence of low-scattering regions. & 2008 Elsevier Ltd. All rights reserved.
Keywords: Optical tomography Radiative transfer equation Inverse problems Finite element method Gauss–Newton method
1. Introduction Optical tomography (OT) is a relatively new imaging modality in which images of the optical properties of the medium are estimated based on measurements of near-infrared light on the surface of the object. The image reconstruction problem in OT is a non-linear ill-posed inverse problem. Thus, even small errors in the measurements or modelling can cause large errors in the reconstructions. To overcome the difficulties due to the ill-posed nature of the problem, regularization techniques need to be used [1]. Furthermore, computationally feasible forward models that describe light propagation within the medium accurately are needed. A widely accepted model for light propagation in tissues is the radiative transfer equation (RTE). However, it is computationally expensive, and therefore the most typical approach in OT has been to use the diffusion approximation (DA) to the RTE as the forward model. The DA is basically a special case of the first order spherical harmonics approximation to the RTE, and thus it has some limitations. Firstly, the medium must be scattering dominated, and secondly, light propagation cannot be modelled accurately close to the collimated light sources and boundaries [1]. Due to these limitations, the DA fails to produce accurate estimates for light propagation close to the sources and boundaries, and in cases in which the turbid medium contains low-scattering or non-scattering regions [2–5]. Therefore, in these situations, more accurate light transport models such as the RTE should be applied.
Corresponding author. Tel.: +358 40 3552310; fax: +358 17 162585.
E-mail addresses: Tanja.Tarvainen@uku.fi (T. Tarvainen), Marko.Vauhkonen@uku.fi (M. Vauhkonen), [email protected] (S.R. Arridge). Also with the Department of Computer Science, University College London, UK.
1
0022-4073/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2008.08.006
ARTICLE IN PRESS 2768
T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
Despite its computationally expensive nature, the RTE has successfully been applied in OT in a few applications. The RTE inverse problem was first addressed in [6] where a transport–backtransport method for a time-dependent RTE was developed. In the approach, the RTE was solved using a finite difference method for spatial discretization and discrete ordinates for angular discretization. Further, the forward and inverse problem of the steady-state RTE were addressed in [4,7] using the finite difference method for the spatial discretization and discrete ordinates for the angular discretization, and an algorithm using a finite element method (FEM) and discrete ordinates was developed in [8]. A method for the forward solution to the steady-state RTE which used the FEM for the spatial discretization and a spherical harmonics approximation for the angular discretization was introduced in [9]. In the frequency domain, numerical solutions to the RTE have been obtained using a finite volume—discrete ordinates method [10], a FEM [11–13], a finite element (FE)—spherical harmonics method [14], and a finite difference—discrete ordinates method [15]. For a recent review of OT reconstruction algorithms that are based on the RTE see [16]. In addition to medical optical imaging, the RTE has been applied in many areas of science such as neutron transport, transfer of charged particles, underwater visibility, marine biology, and radiative transfer in planetary and stellar atmospheres [17,18]. In these cases, the numerical solutions to the RTE have been obtained for example using the FEM [19,20] and the FE—discrete ordinates method [21,22]. In this paper, the frequency domain RTE is solved with the FEM. The FEM is generally regarded as more flexible than other numerical approaches when issues of implementing different boundary conditions and handling complex geometries are considered. In this paper, both the spatial and angular discretizations are implemented in piecewise linear bases. The piecewise linear representation of the angular basis can be regarded as more accurate approach than the discrete ordinates method which basically corresponds to using delta functions as angular bases. Thus, the radiation is represented in all directions instead of just in the discretization angles which is the case in the discrete ordinates method. The ray effect may disturb the standard FE-techniques when solving the RTE [23,24]. The ray effect can produce oscillating results [21] or it can visually be seen as ‘‘photon rays’’ radiating from the source into the direction of the discretization angles. To overcome these problems, a streamline diffusion modification can be utilized. The streamline diffusion modification has been found to stabilize numerical solutions of the RTE in situations in which standard techniques produce oscillating results [21,22,12]. One example of such a situation is the FE-solution to the RTE in a low-scattering medium. In this paper, we utilize the streamline diffusion modification in the FE-solution to the RTE. Most of the approaches to the RTE inverse problem are based on a non-linear least squares approach. In this approach, one minimizes a regularized least squares error between the measured data and data obtained as a solution to the RTE in order to solve the optical properties of the medium. This minimisation problem has been solved using gradient based methods [10,15] and Newton type methods [14]. In [13], a ‘‘special case’’ of the RTE inverse problem using shape boundaries was published. In this paper, we treat the more general problem estimating absorption and scattering distributions within the object. However, some of the notations and methodology are the same. We use a total variation output regularized least squares method for the RTE inverse problem. The total variation prior has been found to suit well for OT reconstruction problem, for examples see [25–27]. The minimisation problem is solved with a Gauss–Newton method which is equipped with a line search algorithm. Gauss–Newton method has successfully been utilized in OT image reconstruction with the DA based approaches, see e.g. [26–28]. The line search algorithm ensures an optimal step length on each iteration. Furthermore, in this paper, scaling of the data and solution spaces is applied in order to improve the numerical stability of the problem. Similar scaling has previously been applied in the image reconstruction problem of OT utilising the DA as the light transport model [28]. One of the most significant problems in OT is an accurate modelling of light transport within low-scattering regions. A typical low-scattering region is the cerebrospinal fluid which surrounds the brain and fills the ventricles. The diffusion theory has been found to fail in situations in which its approximations are not valid such as within low-scattering regions [3–5,29,12,30]. It has also been shown that these limitations in diffusion theory can lead to large errors in the reconstructed images [31,13]. In this paper, we solve the OT reconstruction problem in the presence of low-scattering regions using the RTE as the light transport model. Although utilising the RTE in OT has been relatively intensively studied, these are to our knowledge the first reconstructions in which absorbing inclusions have been estimated in a case in which the object contains also low-scattering regions. The rest of the paper is organized as follows. The OT forward problem and the RTE are described in Section 2. The inverse problem is described in Section 3 along with some transformations which improve the method. The corresponding FE implementations are given in Section 4. The results of simulations are shown in Section 5, and finally conclusions are given in Section 6.
2. Forward problem The forward problem in OT is to solve the measurable data when the optical properties of the medium and the input light sources are given. Light propagation in biological material is usually described through transport theory. The transport theory can be modelled through stochastic methods, such as Monte Carlo, and deterministic methods which are based on describing particle transport with integro-differential equations.
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
2769
2.1. Radiative transfer equation Let O Rn , n ¼ 2 or 3 denote the physical domain which is considered isotropic in the sense that the probability of scattering between two directions depends only on the relative angle between those directions, not on an absolute direction. Furthermore, let qO denote the boundary of the domain and s^ 2 Sn1 denote a unit vector in the direction of interest. The RTE is written in the frequency domain as Z io fðr; s^ Þ þ s^ rfðr; s^ Þ þ ðms þ ma Þfðr; s^Þ ¼ ms Yðs^ s^ 0 Þfðr; s^ 0 Þ ds^ 0 þ qðr; s^ Þ; r 2 O (1) c Sn1 where c is the speed of light in medium, i is the imaginary unit, o is the angular modulation frequency of the input signal, and ms ¼ ms ðrÞ and ma ¼ ma ðrÞ are the scattering and absorption coefficients of the medium, respectively. Furthermore, fðr; s^ Þ is the radiance, Yðs^ s^ 0 Þ is the scattering phase function, and qðr; s^ Þ is the source inside O. In this paper, there are no internal light sources, and thus qðr; s^ Þ ¼ 0. The scattering phase function Yðs^ s^ 0 Þ describes the probability that a photon with an initial direction s^ 0 will have a direction s^ after a scattering event. In OT, the most usual phase function for isotropic material is the Henyey–Greenstein scattering function [32] which in a three-dimensional case is of the form
Yðs^ s^ 0 Þ ¼
1 1 g2 4p ð1 þ g 2 2g s^ s^ 0 Þ3=2
(2)
and in a two-dimensional case is of the form
Yðs^ s^ 0 Þ ¼
1 1 g2 2p ð1 þ g 2 2g s^ s^ 0 Þ
(3)
The scattering shape parameter g defines the shape of the probability density and it gets values between 1ogo1. With the value g ¼ 0, the scattering probability density is a uniform distribution. For forward dominated scattering g40 and for backward dominated scattering go0. Several boundary conditions can be applied to the RTE [33,34]. In OT, we use the boundary condition which assumes that no photons travel in an inward direction at the boundary qO, thus
fðr; s^ Þ ¼ 0;
r 2 qO;
^ s^ no0
(4)
where n^ is a outward unit normal [1]. This boundary condition, also known as the free surface boundary condition and the vacuum boundary condition, implies that once a photon escapes the domain O it does not re-enter it. The boundary condition (4) can be modified to include a boundary source f0 ðr; s^ Þ at the source position ej qO and it can be written in the form [20] ( ^ f0 ðr; s^ Þ; r 2 [j ej ; s^ no0 fðr; s^ Þ ¼ (5) ^ 0; r 2 qOn[j ej ; s^ no0 In OT, the measurable quantity is the exitance GðrÞ on the boundary of the domain. Utilising the boundary condition (4), it can be written as [35] Z Z ^ fðr; s^ Þ ds^ ¼ ^ fðr; s^ Þ ds^ ; r 2 qO GðrÞ ¼ ðs^ nÞ ðs^ nÞ (6) Sn1
^ s^n40
Another quantity of interest is the photon density which is of the form [35] Z FðrÞ ¼ fðr; s^ Þ ds^ Sn1
(7)
3. Inverse problem We consider the solution to the inverse problem in a discrete framework where the optical parameters are represented in a space of dimension K. In general, let us use the following notation: ! x:¼
ma ms
where ma ¼ ðma;1 ; . . . ; ma;K ÞT 2 RK and ms ¼ ðms;1 ; . . . ; ms;K ÞT 2 RK are vectors of discretized absorption and scattering parameters. The actual choice of basis function used in this paper, is defined later in Eqs. (34) and (35). Further, 2 let the measurement vector be Z ¼ ðZ1;1 ; . . . ; Z1;m ; Z2;1 ; . . . ; Zm;m ÞT 2 Cm where m is the number of light sources and detectors. In the inverse problem of OT, absorption and scattering functions within the object are reconstructed. Most of the approaches to the OT inverse problem are based on the non-linear least squares approach. In this approach,
ARTICLE IN PRESS 2770
T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
a single data acquisition is used to estimate absolute values for the optical properties of the medium. The regularized non-linear least squares problem is to estimate the absorption and scattering distributions x which minimize the functional
C ¼ kLðZ FðxÞÞk22 þ BðxÞ
(8)
when the measured data Z are given. In functional (8), F is the forward model for light transport which maps the absorption and scattering parameters to the measurable data. In this study, we use the radiative transfer Eq. (1) as the forward model. Further, matrix L is a weighting matrix which basically corresponds to the Cholesky factor of the inverse of the noise covariance matrix. The term BðxÞ40 is a regularizing penalty functional. Regularization is needed to overcome the difficulties due to the ill-posed nature of the problem. In this paper, we use a total variation prior for regularization. The total variation norm is a good prior for piecewise constant images which consists of a few constant levels with relatively short boundary lines [26]. In this paper, it is constructed similarly as in [26,36] for the DA based image reconstruction problem. 3.1. Scaling In OT, where the dynamic range of the measured light intensities is often very large, scaling of the data may be needed in order to ensure numerical stability of the optimisation algorithm. Furthermore, transformation of the solution space may be used to constitute a right preconditioning. In this paper, we scale data and solution spaces similarly as in [28] for the DA based image reconstruction problem. In the data space, we use the following transformation: ! ! ~A ln jZj Re Z ~ ¼ Z¼ ~j ¼ (9) ln Z argðZÞ Im Z This corresponds to the choice of logarithm of amplitude and phase shift as the measurables. This choice of measurables is widely used in OT. In addition, the data are scaled by inverse values of constants cre and cim which are the average values of the logarithm of amplitude and phase shift of measured data, thus cre ¼
cim ¼
m X m 1 X ~A Z 2 m s¼1 d¼1 s;d
(10)
m X m 1 X ~j Z m2 s¼1 d¼1 s;d
(11)
The purpose of this scaling is to normalize between the logarithmic amplitude and phase [28], and it corresponds to the choice of the weighting matrix of the form !!1 C re (12) L ¼ diag C im 2
2
where C re ¼ ðcre ; . . . ; cre ÞT 2 Rm and C im ¼ ðcim ; . . . ; cim ÞT 2 Rm . In the solution space, we use the following transformation: x~ ¼ ln
x x¯
(13)
where x¯ is the mean of the optical parameters. In this case, it is the mean of the absorption and scattering values of the initial guess for the reconstruction. This scaling makes the units of the solution space ‘‘dimensionless’’. Furthermore, the logarithmic transformation ensures a positive solution. With transformations (9) and (13), the minimized functional of the least squares problem (8) becomes ! Z ~ AF ~ A ð¯xex~ Þ 2 ~ ¼ C þ Bðx~ Þ (14) L ~ j ~ j ð¯xex~ Þ Z F 2
~ is the forward model for light transport which maps the absorption and scattering parameters to the scaled where F measurables. The minimisation problem is solved with the Gauss–Newton method. Thus, the estimated parameters are solved iteratively as T T ~ F ~ ð¯xex~ ÞÞ 1g ðx~ Þ Þ x~ ðiþ1Þ ¼ x~ ðiÞ þ sðiÞ ð~JðiÞ W ~J ðiÞ þ 12HB ðx~ ÞðiÞ Þ1 ð~J ðiÞ WðZ ðiÞ 2 B
(15)
where i is the iteration index, sðiÞ is a step length parameter, ~J is the Jacobian, and W ¼ LT L, where L is as in Eq. (12). The step length sðiÞ can be determined by a line search algorithm in order to achieve a fast convergence of the minimisation problem (14). Furthermore, HB ðx~ Þ and g B ðx~ Þ are the Hessian and gradient of the penalty functional which are constructed
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
similarly as in [36]. The scaled Jacobian ~J is of the form ! A Re ~J ¼ L J diagðxÞ ¼ L diagðFÞ1 J diagðxÞ Im Jj
2771
(16)
where J is the unscaled Jacobian, which is described in detail in Section 4.2. For more information about the choice of the transformations for both the data and solution spaces, see [28,37]. 4. FE implementations 4.1. FE solution to the RTE In this paper, the FEM is used for the solution to the RTE. In the FEM, a variational formulation is derived for the original problem. Then, a finite dimensional approximation for the variational formulation is constructed using suitably chosen basis and test functions in the solution space. In this paper, both the spatial and angular discretizations are implemented in piecewise linear bases. Further, in this paper, the streamline diffusion modification is applied. In the streamline diffusion modification, the test function is of the form ðv þ ds^ rvÞ instead of the standard method in which the test function is v. The parameter d is the ‘‘smoothing’’ parameter which is a spatially varying constant that depends on local absorption and scattering [21]. To obtain the variational formulation for the RTE, Eq. (1) is multiplied with a test function ðv þ ds^ rvÞ and integrated over spatial domain O and angular directions Sn1 . Utilising the boundary condition (5) and taking into account that there are no sources inside the domain, the following variational formulation for the RTE can be derived: Z Z Z Z io fðr; s^ Þðvðr; s^ Þ þ ds^ rvðr; s^ ÞÞ ds^ dr s^ rvðr; s^ Þfðr; s^ Þ ds^ dr n1 O Sn1 c Z O ZS Z Z ^ þ fðr; s^ Þvðr; s^ Þ ds^ dS dðs^ rfðr; s^ ÞÞðs^ rvðr; s^ ÞÞ ds^ dr þ ðs^ nÞ þ n1 qO Sn1 ZO ZS ðms þ ma Þfðr; s^ Þðvðr; s^ Þ þ ds^ rvðr; s^ ÞÞ ds^ dr þ n1 Z ZO ZS m Yðs^ s^ 0 Þfðr; s^0 Þ ds^ 0 ðvðr; s^ Þ þ ds^ rvðr; s^ ÞÞ ds^ dr s O Sn1 Sn1 Z Z ^ f0 ðr; s^ Þvðr; s^ Þ ds^ dS ðs^ nÞ (17) ¼ qO
Sn1
^ denote the positive and negative parts of ðs^ nÞ. ^ where f0 ðr; s^ Þ is the boundary source as in Eq. (5) and ðs^ nÞ In the FEM, the solution fðr; s^ Þ of the variational formulation (17) is approximated in a piecewise linear basis
fðr; s^ Þ fh ðr; s^ Þ ¼
Nn X Na X
fi‘ ci ðrÞc‘ ðs^ Þ
(18)
i¼1 ‘¼1
where ci ðrÞ and c‘ ðs^ Þ are the nodal basis functions of the spatial and angular FE meshes, respectively, fi‘ is the radiance in spatial nodal point i and ‘‘direction’’ ‘; Nn is the number of spatial nodes, and Na is the number of angular directions. The choice of the basis functions, Eq. (18), allows the separation of spatial and angular RTE discretizations, and thus the spatial and angular integrals can be calculated separately and connected using Kronecker products. Therefore, the approach is flexible, e.g. the amount of angular directions can be increased or decreased easily. Inserting approximation (18) into Eq. (17) and choosing the nodal basis functions cj ðrÞ and cm ðs^ Þ as the test functions, the FE-approximation of the RTE can be derived. The FE-approximation of the RTE with the streamline diffusion modification can be written in the form Af ¼ b
(19)
where A ¼ A0 þ A1 þ A2 þ A3 þ A4
(20)
b ¼ b1 q0
(21)
where vector q ¼ ðq01;1 ; ; q0Nn ;Na ÞT 2 RNn Na is the source strength vector 0 j qO such that qi;‘ is the source power at the node i in the ‘‘direction’’ N n Na 0
e
C
having non-zero elements for spatial nodes r i 2 ‘, and f ¼ ðf1;1 ; . . . ; f1;Na ; . . . ; fNn ;1 . . . ; fNn ;Na ÞT 2 is the radiance. The components of Eqs. (20) and (21) are of the form Z Z io ci ðrÞcj ðrÞ dr c‘ ðs^ Þcm ðs^ Þ ds^ A0 ðh; pÞ ¼ c O Sn1 Z Z þ d (22) s^ rcj ðrÞcm ðs^ Þc‘ ðs^ Þ ds^ ci ðrÞ dr O
Sn1
ARTICLE IN PRESS 2772
T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
A1 ðh; pÞ ¼
Z Z n1
ZO SZ
d
þ O
A2 ðh; pÞ ¼
Z qO
A3 ðh; pÞ ¼
s^ rcj ðrÞcm ðs^ Þc‘ ðs^ Þ ds^ ci ðrÞ dr ðs^ rcj ðrÞcm ðs^ ÞÞðs^ rci ðrÞc‘ ðs^ ÞÞ ds^ dr
Sn1
ci ðrÞcj ðrÞ dS
Z Sn1
Z ðms þ ma Þci ðrÞcj ðrÞ dr c‘ ðs^ Þcm ðs^ Þ ds^ O Sn1 Z Z s^ rcj ðrÞcm ðs^ Þc‘ ðs^ Þ ds^ ci ðrÞ dr dðms þ ma Þ þ
(24)
Z
Sn1
O
A4 ðh; pÞ ¼
Z Z
O
O
b1 ðh; pÞ ¼
^ þ c‘ ðs^ Þcm ðs^ Þ ds^ ðs^ nÞ
(23)
Z qO
ms ci ðrÞcj ðrÞ dr dms
Z n1
S
Z
Z Sn1
Sn1
s^ rcj ðrÞcm ðs^ Þ
ci ðrÞcj ðrÞ dS
Z Sn1
(25)
Yðs^ s^ 0 Þc‘ ðs^ 0 Þ ds^ 0 cm ðs^ Þ ds^ Z Sn1
Yðs^ s^ 0 Þc‘ ðs^ 0 Þ ds^ 0 ds^ ci ðrÞ dr
^ c‘ ðs^ Þcm ðs^ Þ ds^ ðs^ nÞ
(26)
(27)
where h ¼ N a ðj 1Þ þ m; p ¼ Na ði 1Þ þ ‘ (j; i ¼ 1; . . . ; Nn , m; ‘ ¼ 1; . . . ; Na , h; p ¼ 1; . . . ; Nn Na ). The FE-solution to the RTE is obtained using Eq. (19). Thus, as a solution to the forward problem, the radiance f in nodal points of the spatial and angular discretizations is obtained. The exitance, which is the measurable quantity, can be calculated using Eq. (6) as
Gh ¼ Mf ¼ MA1 b
(28)
where M 2 RmNn Na is a measurement matrix which contains the discretized measurement operator and m is the number of detectors. For detailed derivation of the FE-approximation of the RTE see e.g. [11,38]. 4.2. Jacobian In the Gauss–Newton iteration algorithm (15), the Jacobian matrix J needs to be solved on each iteration round. This can be efficiently done by using an adjoint method. We first consider the analytic Fre´chet derivative of the forward mapping for the RTE given by [1,6] ! Z C f F0 Z ¼ (29) R 0 0 C ð Sn1 Yðs^ s^ Þf ds^ fÞ ds^ Sn1 where
denotes conjugate transpose, and C is the solution to the adjoint RTE Z Yðs^ s^ 0 ÞCðr; s^ 0 Þ ds^ 0 ; r 2 O
io Cðr; s^ Þ s^ rCðr; s^ Þ þ ðms þ ma ÞCðr; s^ Þ ¼ ms c
Sn1
(30)
with boundary condition ^ Z; Cðr; s^ Þ ¼ ðs^ nÞ
^ r 2 qO; s^ no0
(31)
indicating the redistribution of the measurement isotropically into incoming directions. If we consider one measurement position ei and one source position ej , we can derive the sensitivity functions Z rðma Þ ðrÞ ¼ Ci ðr; s^ Þfj ðr; s^ Þ ds^ Sn1
rðms Þ ðrÞ ¼
Z Sn1
Ci ðr; s^ Þfj ðr; s^ Þ ds^
Z
Z Sn1
Sn1
Yðs^ s^ 0 ÞCi ðr; s^ 0 Þ ds^ 0 fj ðr; s^ Þ ds^
(32)
(33)
Let us assume that the domain O is divided into K disjoint sub-domains Ok , and that the absorption and scattering coefficients are expressed in basis
ma ðrÞ
K X
mak wkðma Þ ðrÞ
(34)
msk wkðms Þ ðrÞ
(35)
k¼1
ms ðrÞ
K X k¼1
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
2773
where wk ðrÞ is a characteristic function of the sub-domain Ok . The kth column of the Jacobian J for the absorption or scattering is a vectorized matrix Z ðm Þ jk a ¼ wkðma Þ ðrÞrðma Þ ðrÞ (36) O
ðms Þ
jk
Z ¼ O
wðkms Þ ðrÞrðms Þ ðrÞ
(37)
We should note that this expression for the Jacobian, based on projection of the sensitivity functions into the finite dimensional space spanned by the characteristic functions wk , is not the same as would be obtained by differentiation of the discrete matrix A, because the latter would also include terms involving the streamline diffusion integrals. 5. Results The performance of the RTE based reconstruction method was tested with 2D simulations. In simulations, a highly scattering circular domain O which contained different inclusions was investigated. The radius of the circle was 10 mm, and 16 equally spaced sources and detectors located on the boundary of the domain. The modulation frequency of the input signal was chosen as 100 MHz. The refractive indices of the background medium, inclusion, and the surrounding medium were chosen as nin ¼ 1. Four types of cases were investigated: a domain with a highly absorbing circular blob and a highly scattering circular blob (case 1), a domain with a highly absorbing circular blob and a low-scattering circular blob located both far apart and close together (cases 2 and 3, respectively), and a domain with a highly absorbing blob and a thin low-scattering gap (case 4). The optical properties of the simulation domains are summarized in Table 1. In all of the simulations, the background absorption and scattering coefficients were ma ¼ 0:025 mm1 and ms ¼ 2 mm1 , respectively. Further, in all of the simulations the scattering shape parameter was g ¼ 0:2 throughout the domain. The simulated absorption and scattering distributions are shown on left columns of Figs. 3–6. The data were generated with the RTE using the FEM as described in Section 4.1. The RTE was solved using Eq. (19) where the FE-integrals were calculated as in Eqs. (22)–(27). Further, the exitance was solved using Eq. (28). The spatial FE-discretization of case 1 (two blobs) contained 8589 nodal points and 16 792 triangular elements. The spatial FE-discretization of the case 2 (two blobs far apart) contained 8589 nodal points and 16 792 triangular elements and the FE-discretization of case 3 (two blobs close together) contained 8437 nodal points and 16 488 triangular elements. Further, the FE-discretization of the case 4 (a blob and a gap) contained 8545 nodal points and 16 704 triangular elements. In all of the simulations, the angular discretization contained 16 angular directions. Gaussian distributed noise with standard deviation of 1% of corresponding amplitude and phase was added to the simulated data. 5.1. Error surfaces First, we investigated the effect of the scaling on the minimisation problem. Therefore, we computed the norms (8) where L was a unit matrix, and (14) where L was as in Eq. (12) using the simulated data Z of case 1 (two blobs in highly scattering medium) and case 4 (a highly absorbing blob and a low-scattering gap) and forward solution F calculated with different optical properties. The forward data were generated for all combinations of 10 homogeneous absorption values on the range ma 2 ½0:005; 0:05 mm1 and 10 homogeneous scattering values on the range ms 2 ½0:5; 5 mm1 . The obtained error maps are shown in Figs. 1 and 2 where the objective function (8) (unscaled data) is on the left and the objective function (14) (scaled data) is on the right. As can be seen from Fig. 1, the error map from unscaled data shows a different shape profile than the error map from the scaled data in the case of a highly scattering medium. The error map from the scaled data shows a more clearly defined minimum than the unscaled data. This indicates that the scaling provide a better means to separate absorption and scattering. In the case of a highly absorbing blob and a low-scattering gap (Fig. 2), the error map from the scaled data shows again a more clearly defined minimum than the error map from the unscaled data. However, the minimum is not as clear as it was in the case of a highly scattering medium. This indicates, that the reconstruction from a medium with a low-scattering gap is challenging and that the absorption and scattering may not be that well distinguished. Table 1 The absorption coefficient ma and the scattering coefficient ms of the background medium and the inclusions Background
Case 1 Cases 2–3 Case 4
Blob 1
Blob 2/gap
ma ðmm1 Þ
ms ðmm1 Þ
ma ðmm1 Þ
ms ðmm1 Þ
ma ðmm1 Þ
ms ðmm1 Þ
0.025 0.025 0.025
2 2 2
0.125 0.125 0.125
2 2 2
0.025 0.025 0.01
5 0.02 0.02
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
0.05
0.05
0.04
0.04
0.03
0.03
a
a
2774
0.02
0.02
0.01
0.01 1
2
s
3
4
5
1
2
s
3
4
5
0.05
0.05
0.04
0.04
0.03
0.03
a
a
Fig. 1. Objective function for case 1 (a highly scattering medium) as a function of scattering and absorption. On the left, the minimized functional was (8) where L was a unit matrix, and on the right, the minimized functional was (14).
0.02
0.02
0.01
0.01 1
2
s
3
4
5
1
2
s
3
4
5
Fig. 2. Objective function for case 4 (a highly absorbing blob and a low-scattering gap) as a function of scattering and absorption. On the left, the minimized functional was (8) where L was a unit matrix, and on the right, the minimized functional was (14).
5.2. Reconstructions The inverse problem was solved with an image reconstruction method described in Section 3 using the RTE as light transport model. The reconstruction mesh contained 3644 triangular elements. Note that the reconstruction mesh was coarser than the data generation meshes in order to mitigate against the possibility of inverse crimes. First, an initial guess for absorption and scattering values was estimated by minimising the least squares error between simulated measurement and forward solution obtained with homogeneous absorption and scattering. Then, the absorption and scattering reconstructions were computed using the total variation regularized output least squares scheme by minimising functional (14). The previously solved homogeneous absorption and scattering values were used as an initial point of the reconstruction algorithm. In the image reconstruction procedure, the data were scaled using transformation (9), and the weighting matrix L was chosen as in Eq. (12). In the solution space, we used transformation (13). The minimisation problem was solved with Gauss–Newton algorithm (15) which was equipped with a line search algorithm for the determination of the step length. In all of the cases, convergence was achieved in 15 iterations.
5.2.1. Case 1: a highly scattering medium As a first case, we investigated a situation in which the medium contained a highly absorbing blob and a highly scattering blob. The radii of the blobs were 2 mm. The absorption and scattering of the background medium were, as in all of the simulations, ma ¼ 0:025 mm1 and ms ¼ 2 mm1 . The absorption and scattering within the highly absorbing blob were ma ¼ 0:125 mm1 and ms ¼ 2 mm1 , and the absorption and scattering within the highly scattering blob were
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
2775
Fig. 3. Absorption coefficients (top row) and scattering coefficients (bottom row) within a domain with a highly absorbing blob and a highly scattering blob (case 1). Simulated distributions are on the left column and the reconstructions are on the right column.
ma ¼ 0:025 mm1 and ms ¼ 5 mm1 . The simulated absorption and scattering distributions can be seen on the left column of Fig. 3. The reconstructions can be seen on the right column of Fig. 3. The absorption distributions are on the top row and the scattering distributions are on the bottom row. As it can be seen from Fig. 3, the reconstructions are clear and the inclusions are well distinguished. Some crosstalk between the scattering and absorption appears in the reconstructions. This, however, is typical for OT reconstructions, and although it has been shown that utilising the frequency domain RTE the amount of crosstalk is significantly smaller than with the steady-state RTE, there still exists some crosstalk between the optical parameters [10,15].
5.2.2. Cases 2 and 3: two blobs far apart and close together Then we investigated a situation in which the highly scattering medium contained a highly absorbing blob and a low-scattering blob. Two situations were considered: one in which the blobs were located far apart (case 2) and an another in which the blobs were located close together (case 3). The radii of the blobs were 2 mm as in the previous simulation. The absorption and scattering within the highly absorbing blob were ma ¼ 0:125 mm1 and ms ¼ 2 mm1 , and the absorption and scattering within the low-scattering blob were ma ¼ 0:025 mm1 and ms ¼ 0:02 mm1 . The absorption and scattering values of the background medium and the blobs are summarized in Table 1 and the simulated distributions can be seen on the left columns of Figs. 4 and 5. The reconstructions from case 2 in which the domain contained two blobs far apart can be seen in Fig. 4, and the reconstructions from case 3 in which the domain contained two blobs close together can be seen in Fig. 5. Figures show absorption coefficients on the top row and scattering coefficients on the bottom row. The simulated distributions are on the left column and the reconstructions are on the right column. As it can be seen, both the low-scattering blob and the highly absorbing blob can be well distinguished. Even in the case in which the low-scattering blob is located close to the absorbing blob (Fig. 5), the reconstruction shows a good quality absorption distribution with a clear blob.
5.2.3. Case 4: a highly absorbing blob and a low-scattering gap As a fourth case we investigated a situation in which the highly scattering medium contained a highly absorbing blob and a thin low-scattering and low-absorbing gap. The radius of the blob was 2 mm and the thickness of the gap was 0:5 mm. The absorption and scattering within the highly absorbing blob were ma ¼ 0:125 mm1 and ms ¼ 2 mm1 , and the absorption and scattering within the gap were ma ¼ 0:01 mm1 and ms ¼ 0:02 mm1 . The absorption and scattering values are summarized in Table 1 and the simulated distributions are shown on the left column of Fig. 6.
ARTICLE IN PRESS 2776
T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
Fig. 4. Absorption coefficients (top row) and scattering coefficients (bottom row) within a domain with a highly absorbing blob and a low-scattering blob (case 2). Simulated distributions are on the left column and the reconstructions are on the right column.
Fig. 5. Absorption coefficients (top row) and scattering coefficients (bottom row) within a domain with a highly absorbing blob and a low-scattering blob (case 3). Simulated distributions are on the left column and the reconstructions are on the right column.
The reconstructions are shown on the right column of Fig. 6. As it can be seen from Fig. 6, the RTE can detect the location of the gap but not the absolute absorption and scattering values. The highly absorbing blob can also be detected, although the absolute values of absorption are not as good as in the previous simulations.
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
2777
Fig. 6. Absorption coefficients (top row) and scattering coefficients (bottom row) within a domain with a highly absorbing blob and a low-absorbing and low-scattering gap (case 4). Simulated distributions are on the left column and the reconstructions are on the right column.
6. Conclusions In this paper, the forward and inverse problem of the frequency domain radiative transfer equation were considered. An image reconstruction method based on a total variation output regularized least squares method was developed. In the approach, the frequency domain RTE was used as the light transport model. The minimisation problem was solved with a Gauss–Newton algorithm which was equipped with a line search algorithm. Furthermore, scaling of the data and solution spaces were applied to improve the numerical stability of the problem. The radiative transfer equation was numerically solved with the FEM. In the FE-solution to the RTE both the spatial and angular discretizations were implemented in piecewise linear bases. Furthermore, the streamline diffusion modification was utilized to improve the numerical stability of the problem. The approach was tested with simulations. Reconstructions from different cases including domains with low-scattering regions were shown. The results show that the frequency domain RTE can be used in optical tomography to produce good quality reconstructions.
Acknowledgements This work was supported by the Academy of Finland (project 122499) and by EPSRC grant EP/E034950/1. References [1] Arridge SR. Optical tomography in medical imaging. Inverse Probl 1999;15:R41–93. [2] Wang L, Jacques SL. Hybrid model of Monte Carlo simulation and diffusion theory for light reflectance by turbid media. J Opt Soc Am A 1993;10(8): 1746–52. [3] Firbank M, Arridge SR, Schweiger M, Delpy DT. An investigation of light transport through scattering bodies with non-scattering regions. Phys Med Biol 1996;41:767–83. [4] Hielscher AH, Alcouffe RE, Barbour RL. Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues. Phys Med Biol 1998;43:1285–302. [5] Okada E, Firbank M, Schweiger M, Arridge SR, Cope M, Delpy DT. Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head. Appl Opt 1997;36(1):21–31. [6] Dorn O. A transport–backtransport method for optical tomography. Inverse Probl 1998;14(5):1107–30. [7] Klose AD, Hielscher AH. Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer. Med Phys 1999;26(8): 1698–707. [8] Abdoulaev G, Hielscher A. Three-dimensional optical tomography with the equation of radiative transfer. JEI 2003;12(4):594–601.
ARTICLE IN PRESS 2778
T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
[9] Aydin E, de Oliveira C, Goddard A. A comparison between transport and diffusion calculations using a finite element-spherical harmonics radiation transport method. Med Phys 2002;29(9):2013–23. [10] Ren K, Bal G, Hielscher A. Frequency domain optical tomography based on the equation of radiative transfer. SIAM J Sci Comput 2006;28(4):1463–89. [11] Tarvainen T, Vauhkonen M, Kolehmainen V, Kaipio J. Hybrid radiative–transfer–diffusion model for optical tomography. Appl Opt 2005;44(6): 876–86. [12] Tarvainen T, Vauhkonen M, Kolehmainen V, Arridge S, Kaipio J. Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and non-scattering regions. Phys Med Biol 2005;50:4913–30. [13] Arridge SR, Dorn O, Kaipio JP, Kolehmainen V, Schweiger M, Tarvainen T, et al. Reconstruction of subdomain boundaries of piecewise constant coefficients of the radiative transfer equation from optical tomography data. Inverse Probl 2006;22:2175–96. [14] Wright S, Schweiger M, Arridge S. Reconstruction in optical tomography using the P N approximations. Meas Sci Technol 2007;18:79–86. [15] Kim HK, Charette A. A sensitivity function-based conjugate gradient method for optical tomography with the frequency-domain equation of radiative transfer. JQSRT 2007;104:24–39. [16] Klose AD, Hielscher AH. Optical tomography with the equation of radiative transfer. Int J Numer Methods Heat Fluid Flow 2008;18(3/4):443–64. [17] Case KM, Zweifel PF. Linear transport theory. Reading: Addison-Wesley; 1967. [18] Ishimaru A. Wave propagation and scattering in random media, vol. 1. New York: Academic Press; 1978. [19] Boman E, Tervo J, Vauhkonen M. Modelling the transport of ionizing radiation using the finite element method. Phys Med Biol 2005;50:265–80. [20] Tervo J, Kolmonen P, Vauhkonen M, Heikkinen L, Kaipio J. A finite-element model of electron transport in radiation therapy and a related inverse problem. Inverse Probl 1999;15:1345–61. [21] Kanschat G. A robust finite element discretization for radiative transfer problems with scattering. East–West J Numer Math 1998;6(4):265–72. [22] Richling S, Meinko¨hn E, Kryzhevoi N, Kanschat G. Radiative transfer with finite elements I. Basic method and tests. Astron Astrophys 2001;380(2):776–88. [23] Lathrop K. Ray effects in discrete ordinates equations. Nucl Sci Eng 1968;32:357–69. [24] Lathrop K. Remedies for ray effects. Nucl Sci Eng 1971;45:255–68. [25] Paulsen KD, Jiang H. Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization. Appl Opt 1996;35(19):3447–58. [26] Kolehmainen V, Vauhkonen M, Kaipio J, Arridge SR. Recovery of piecewise constant coefficients in optical diffusion tomography. Opt Express 2000;7(13):468–80. [27] Tarvainen T, Kolehmainen V, Vauhkonen M, Vanne A, Gibson A, Schweiger M, et al. Computational calibration method for optical tomography. Appl Opt 2005;44(10):1879–88. [28] Schweiger M, Arridge SR, Nissila¨ I. Gauss–Newton method for image reconstruction in diffuse optical tomography. Phys Med Biol 2005;50:2365–86. [29] Arridge SR, Dehghani H, Schweiger M, Okada E. The finite element model for the propagation of light in scattering media: a direct method for domains with nonscattering regions. Med Phys 2000;27(1):252–64. [30] Aydin ED. Three-dimensional photon migration through voidlike regions and channels. Appl Opt 2007;46(34):8272–7. [31] Dehghani H, Arridge SR, Schweiger M, Delpy DT. Optical tomography in the presence of void regions. J Opt Soc Am A 2000;17(9):1659–70. [32] Henyey LG, Greenstein JL. Diffuse radiation in the galaxy. Astrophys. J. 1941;93:70–83. [33] Duderstadt JJ, Martin WR. Transport theory. New York: Wiley; 1979. [34] Kaipio J, Somersalo E. Statistical and computational inverse problems. New York: Springer; 2005. [35] Arridge SR, Hebden JC. Optical imaging in medicine: II. Modelling and reconstruction. Phys Med Biol 1997;42:841–53. [36] Kolehmainen V. Novel approaches to image reconstruction in diffusion tomography. PhD thesis, University of Kuopio, Kuopio, Finland, 2001. [37] Schweiger M, Arridge SR. Application of temporal filters to time resolved data in optical tomography. Phys Med Biol 1999;44:1699–717. [38] Tarvainen T, Vauhkonen M, Kolehmainen V, Kaipio J. Finite element model for the coupled radiative transfer equation and diffusion approximation. Int J Numer Methods Eng 2006;65(3):383–405.
Contents lists available at ScienceDirect
Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt
Gauss–Newton reconstruction method for optical tomography using the finite element solution of the radiative transfer equation T. Tarvainen a,,1, M. Vauhkonen a, S.R. Arridge b a b
Department of Physics, University of Kuopio, P.O. Box 1627, 70211 Kuopio, Finland Department of Computer Science, University College London, Gower Street, London WC1E 6BT, UK
a r t i c l e in fo
abstract
Article history: Received 10 July 2008 Received in revised form 29 August 2008 Accepted 29 August 2008
The radiative transfer equation can be utilized in optical tomography in situations in which the more commonly applied diffusion approximation is not valid. In this paper, an image reconstruction method based on a frequency domain radiative transfer equation is developed. The approach is based on a total variation output regularized least squares method which is solved with a Gauss–Newton algorithm. The radiative transfer equation is numerically solved with a finite element method in which both the spatial and angular discretizations are implemented in piecewise linear bases. Furthermore, the streamline diffusion modification is utilized to improve the numerical stability. The approach is tested with simulations. Reconstructions from different cases including domains with low-scattering regions are shown. The results show that the radiative transfer equation can be utilized in optical tomography and it can produce good quality images even in the presence of low-scattering regions. & 2008 Elsevier Ltd. All rights reserved.
Keywords: Optical tomography Radiative transfer equation Inverse problems Finite element method Gauss–Newton method
1. Introduction Optical tomography (OT) is a relatively new imaging modality in which images of the optical properties of the medium are estimated based on measurements of near-infrared light on the surface of the object. The image reconstruction problem in OT is a non-linear ill-posed inverse problem. Thus, even small errors in the measurements or modelling can cause large errors in the reconstructions. To overcome the difficulties due to the ill-posed nature of the problem, regularization techniques need to be used [1]. Furthermore, computationally feasible forward models that describe light propagation within the medium accurately are needed. A widely accepted model for light propagation in tissues is the radiative transfer equation (RTE). However, it is computationally expensive, and therefore the most typical approach in OT has been to use the diffusion approximation (DA) to the RTE as the forward model. The DA is basically a special case of the first order spherical harmonics approximation to the RTE, and thus it has some limitations. Firstly, the medium must be scattering dominated, and secondly, light propagation cannot be modelled accurately close to the collimated light sources and boundaries [1]. Due to these limitations, the DA fails to produce accurate estimates for light propagation close to the sources and boundaries, and in cases in which the turbid medium contains low-scattering or non-scattering regions [2–5]. Therefore, in these situations, more accurate light transport models such as the RTE should be applied.
Corresponding author. Tel.: +358 40 3552310; fax: +358 17 162585.
E-mail addresses: Tanja.Tarvainen@uku.fi (T. Tarvainen), Marko.Vauhkonen@uku.fi (M. Vauhkonen), [email protected] (S.R. Arridge). Also with the Department of Computer Science, University College London, UK.
1
0022-4073/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.jqsrt.2008.08.006
ARTICLE IN PRESS 2768
T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
Despite its computationally expensive nature, the RTE has successfully been applied in OT in a few applications. The RTE inverse problem was first addressed in [6] where a transport–backtransport method for a time-dependent RTE was developed. In the approach, the RTE was solved using a finite difference method for spatial discretization and discrete ordinates for angular discretization. Further, the forward and inverse problem of the steady-state RTE were addressed in [4,7] using the finite difference method for the spatial discretization and discrete ordinates for the angular discretization, and an algorithm using a finite element method (FEM) and discrete ordinates was developed in [8]. A method for the forward solution to the steady-state RTE which used the FEM for the spatial discretization and a spherical harmonics approximation for the angular discretization was introduced in [9]. In the frequency domain, numerical solutions to the RTE have been obtained using a finite volume—discrete ordinates method [10], a FEM [11–13], a finite element (FE)—spherical harmonics method [14], and a finite difference—discrete ordinates method [15]. For a recent review of OT reconstruction algorithms that are based on the RTE see [16]. In addition to medical optical imaging, the RTE has been applied in many areas of science such as neutron transport, transfer of charged particles, underwater visibility, marine biology, and radiative transfer in planetary and stellar atmospheres [17,18]. In these cases, the numerical solutions to the RTE have been obtained for example using the FEM [19,20] and the FE—discrete ordinates method [21,22]. In this paper, the frequency domain RTE is solved with the FEM. The FEM is generally regarded as more flexible than other numerical approaches when issues of implementing different boundary conditions and handling complex geometries are considered. In this paper, both the spatial and angular discretizations are implemented in piecewise linear bases. The piecewise linear representation of the angular basis can be regarded as more accurate approach than the discrete ordinates method which basically corresponds to using delta functions as angular bases. Thus, the radiation is represented in all directions instead of just in the discretization angles which is the case in the discrete ordinates method. The ray effect may disturb the standard FE-techniques when solving the RTE [23,24]. The ray effect can produce oscillating results [21] or it can visually be seen as ‘‘photon rays’’ radiating from the source into the direction of the discretization angles. To overcome these problems, a streamline diffusion modification can be utilized. The streamline diffusion modification has been found to stabilize numerical solutions of the RTE in situations in which standard techniques produce oscillating results [21,22,12]. One example of such a situation is the FE-solution to the RTE in a low-scattering medium. In this paper, we utilize the streamline diffusion modification in the FE-solution to the RTE. Most of the approaches to the RTE inverse problem are based on a non-linear least squares approach. In this approach, one minimizes a regularized least squares error between the measured data and data obtained as a solution to the RTE in order to solve the optical properties of the medium. This minimisation problem has been solved using gradient based methods [10,15] and Newton type methods [14]. In [13], a ‘‘special case’’ of the RTE inverse problem using shape boundaries was published. In this paper, we treat the more general problem estimating absorption and scattering distributions within the object. However, some of the notations and methodology are the same. We use a total variation output regularized least squares method for the RTE inverse problem. The total variation prior has been found to suit well for OT reconstruction problem, for examples see [25–27]. The minimisation problem is solved with a Gauss–Newton method which is equipped with a line search algorithm. Gauss–Newton method has successfully been utilized in OT image reconstruction with the DA based approaches, see e.g. [26–28]. The line search algorithm ensures an optimal step length on each iteration. Furthermore, in this paper, scaling of the data and solution spaces is applied in order to improve the numerical stability of the problem. Similar scaling has previously been applied in the image reconstruction problem of OT utilising the DA as the light transport model [28]. One of the most significant problems in OT is an accurate modelling of light transport within low-scattering regions. A typical low-scattering region is the cerebrospinal fluid which surrounds the brain and fills the ventricles. The diffusion theory has been found to fail in situations in which its approximations are not valid such as within low-scattering regions [3–5,29,12,30]. It has also been shown that these limitations in diffusion theory can lead to large errors in the reconstructed images [31,13]. In this paper, we solve the OT reconstruction problem in the presence of low-scattering regions using the RTE as the light transport model. Although utilising the RTE in OT has been relatively intensively studied, these are to our knowledge the first reconstructions in which absorbing inclusions have been estimated in a case in which the object contains also low-scattering regions. The rest of the paper is organized as follows. The OT forward problem and the RTE are described in Section 2. The inverse problem is described in Section 3 along with some transformations which improve the method. The corresponding FE implementations are given in Section 4. The results of simulations are shown in Section 5, and finally conclusions are given in Section 6.
2. Forward problem The forward problem in OT is to solve the measurable data when the optical properties of the medium and the input light sources are given. Light propagation in biological material is usually described through transport theory. The transport theory can be modelled through stochastic methods, such as Monte Carlo, and deterministic methods which are based on describing particle transport with integro-differential equations.
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
2769
2.1. Radiative transfer equation Let O Rn , n ¼ 2 or 3 denote the physical domain which is considered isotropic in the sense that the probability of scattering between two directions depends only on the relative angle between those directions, not on an absolute direction. Furthermore, let qO denote the boundary of the domain and s^ 2 Sn1 denote a unit vector in the direction of interest. The RTE is written in the frequency domain as Z io fðr; s^ Þ þ s^ rfðr; s^ Þ þ ðms þ ma Þfðr; s^Þ ¼ ms Yðs^ s^ 0 Þfðr; s^ 0 Þ ds^ 0 þ qðr; s^ Þ; r 2 O (1) c Sn1 where c is the speed of light in medium, i is the imaginary unit, o is the angular modulation frequency of the input signal, and ms ¼ ms ðrÞ and ma ¼ ma ðrÞ are the scattering and absorption coefficients of the medium, respectively. Furthermore, fðr; s^ Þ is the radiance, Yðs^ s^ 0 Þ is the scattering phase function, and qðr; s^ Þ is the source inside O. In this paper, there are no internal light sources, and thus qðr; s^ Þ ¼ 0. The scattering phase function Yðs^ s^ 0 Þ describes the probability that a photon with an initial direction s^ 0 will have a direction s^ after a scattering event. In OT, the most usual phase function for isotropic material is the Henyey–Greenstein scattering function [32] which in a three-dimensional case is of the form
Yðs^ s^ 0 Þ ¼
1 1 g2 4p ð1 þ g 2 2g s^ s^ 0 Þ3=2
(2)
and in a two-dimensional case is of the form
Yðs^ s^ 0 Þ ¼
1 1 g2 2p ð1 þ g 2 2g s^ s^ 0 Þ
(3)
The scattering shape parameter g defines the shape of the probability density and it gets values between 1ogo1. With the value g ¼ 0, the scattering probability density is a uniform distribution. For forward dominated scattering g40 and for backward dominated scattering go0. Several boundary conditions can be applied to the RTE [33,34]. In OT, we use the boundary condition which assumes that no photons travel in an inward direction at the boundary qO, thus
fðr; s^ Þ ¼ 0;
r 2 qO;
^ s^ no0
(4)
where n^ is a outward unit normal [1]. This boundary condition, also known as the free surface boundary condition and the vacuum boundary condition, implies that once a photon escapes the domain O it does not re-enter it. The boundary condition (4) can be modified to include a boundary source f0 ðr; s^ Þ at the source position ej qO and it can be written in the form [20] ( ^ f0 ðr; s^ Þ; r 2 [j ej ; s^ no0 fðr; s^ Þ ¼ (5) ^ 0; r 2 qOn[j ej ; s^ no0 In OT, the measurable quantity is the exitance GðrÞ on the boundary of the domain. Utilising the boundary condition (4), it can be written as [35] Z Z ^ fðr; s^ Þ ds^ ¼ ^ fðr; s^ Þ ds^ ; r 2 qO GðrÞ ¼ ðs^ nÞ ðs^ nÞ (6) Sn1
^ s^n40
Another quantity of interest is the photon density which is of the form [35] Z FðrÞ ¼ fðr; s^ Þ ds^ Sn1
(7)
3. Inverse problem We consider the solution to the inverse problem in a discrete framework where the optical parameters are represented in a space of dimension K. In general, let us use the following notation: ! x:¼
ma ms
where ma ¼ ðma;1 ; . . . ; ma;K ÞT 2 RK and ms ¼ ðms;1 ; . . . ; ms;K ÞT 2 RK are vectors of discretized absorption and scattering parameters. The actual choice of basis function used in this paper, is defined later in Eqs. (34) and (35). Further, 2 let the measurement vector be Z ¼ ðZ1;1 ; . . . ; Z1;m ; Z2;1 ; . . . ; Zm;m ÞT 2 Cm where m is the number of light sources and detectors. In the inverse problem of OT, absorption and scattering functions within the object are reconstructed. Most of the approaches to the OT inverse problem are based on the non-linear least squares approach. In this approach,
ARTICLE IN PRESS 2770
T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
a single data acquisition is used to estimate absolute values for the optical properties of the medium. The regularized non-linear least squares problem is to estimate the absorption and scattering distributions x which minimize the functional
C ¼ kLðZ FðxÞÞk22 þ BðxÞ
(8)
when the measured data Z are given. In functional (8), F is the forward model for light transport which maps the absorption and scattering parameters to the measurable data. In this study, we use the radiative transfer Eq. (1) as the forward model. Further, matrix L is a weighting matrix which basically corresponds to the Cholesky factor of the inverse of the noise covariance matrix. The term BðxÞ40 is a regularizing penalty functional. Regularization is needed to overcome the difficulties due to the ill-posed nature of the problem. In this paper, we use a total variation prior for regularization. The total variation norm is a good prior for piecewise constant images which consists of a few constant levels with relatively short boundary lines [26]. In this paper, it is constructed similarly as in [26,36] for the DA based image reconstruction problem. 3.1. Scaling In OT, where the dynamic range of the measured light intensities is often very large, scaling of the data may be needed in order to ensure numerical stability of the optimisation algorithm. Furthermore, transformation of the solution space may be used to constitute a right preconditioning. In this paper, we scale data and solution spaces similarly as in [28] for the DA based image reconstruction problem. In the data space, we use the following transformation: ! ! ~A ln jZj Re Z ~ ¼ Z¼ ~j ¼ (9) ln Z argðZÞ Im Z This corresponds to the choice of logarithm of amplitude and phase shift as the measurables. This choice of measurables is widely used in OT. In addition, the data are scaled by inverse values of constants cre and cim which are the average values of the logarithm of amplitude and phase shift of measured data, thus cre ¼
cim ¼
m X m 1 X ~A Z 2 m s¼1 d¼1 s;d
(10)
m X m 1 X ~j Z m2 s¼1 d¼1 s;d
(11)
The purpose of this scaling is to normalize between the logarithmic amplitude and phase [28], and it corresponds to the choice of the weighting matrix of the form !!1 C re (12) L ¼ diag C im 2
2
where C re ¼ ðcre ; . . . ; cre ÞT 2 Rm and C im ¼ ðcim ; . . . ; cim ÞT 2 Rm . In the solution space, we use the following transformation: x~ ¼ ln
x x¯
(13)
where x¯ is the mean of the optical parameters. In this case, it is the mean of the absorption and scattering values of the initial guess for the reconstruction. This scaling makes the units of the solution space ‘‘dimensionless’’. Furthermore, the logarithmic transformation ensures a positive solution. With transformations (9) and (13), the minimized functional of the least squares problem (8) becomes ! Z ~ AF ~ A ð¯xex~ Þ 2 ~ ¼ C þ Bðx~ Þ (14) L ~ j ~ j ð¯xex~ Þ Z F 2
~ is the forward model for light transport which maps the absorption and scattering parameters to the scaled where F measurables. The minimisation problem is solved with the Gauss–Newton method. Thus, the estimated parameters are solved iteratively as T T ~ F ~ ð¯xex~ ÞÞ 1g ðx~ Þ Þ x~ ðiþ1Þ ¼ x~ ðiÞ þ sðiÞ ð~JðiÞ W ~J ðiÞ þ 12HB ðx~ ÞðiÞ Þ1 ð~J ðiÞ WðZ ðiÞ 2 B
(15)
where i is the iteration index, sðiÞ is a step length parameter, ~J is the Jacobian, and W ¼ LT L, where L is as in Eq. (12). The step length sðiÞ can be determined by a line search algorithm in order to achieve a fast convergence of the minimisation problem (14). Furthermore, HB ðx~ Þ and g B ðx~ Þ are the Hessian and gradient of the penalty functional which are constructed
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
similarly as in [36]. The scaled Jacobian ~J is of the form ! A Re ~J ¼ L J diagðxÞ ¼ L diagðFÞ1 J diagðxÞ Im Jj
2771
(16)
where J is the unscaled Jacobian, which is described in detail in Section 4.2. For more information about the choice of the transformations for both the data and solution spaces, see [28,37]. 4. FE implementations 4.1. FE solution to the RTE In this paper, the FEM is used for the solution to the RTE. In the FEM, a variational formulation is derived for the original problem. Then, a finite dimensional approximation for the variational formulation is constructed using suitably chosen basis and test functions in the solution space. In this paper, both the spatial and angular discretizations are implemented in piecewise linear bases. Further, in this paper, the streamline diffusion modification is applied. In the streamline diffusion modification, the test function is of the form ðv þ ds^ rvÞ instead of the standard method in which the test function is v. The parameter d is the ‘‘smoothing’’ parameter which is a spatially varying constant that depends on local absorption and scattering [21]. To obtain the variational formulation for the RTE, Eq. (1) is multiplied with a test function ðv þ ds^ rvÞ and integrated over spatial domain O and angular directions Sn1 . Utilising the boundary condition (5) and taking into account that there are no sources inside the domain, the following variational formulation for the RTE can be derived: Z Z Z Z io fðr; s^ Þðvðr; s^ Þ þ ds^ rvðr; s^ ÞÞ ds^ dr s^ rvðr; s^ Þfðr; s^ Þ ds^ dr n1 O Sn1 c Z O ZS Z Z ^ þ fðr; s^ Þvðr; s^ Þ ds^ dS dðs^ rfðr; s^ ÞÞðs^ rvðr; s^ ÞÞ ds^ dr þ ðs^ nÞ þ n1 qO Sn1 ZO ZS ðms þ ma Þfðr; s^ Þðvðr; s^ Þ þ ds^ rvðr; s^ ÞÞ ds^ dr þ n1 Z ZO ZS m Yðs^ s^ 0 Þfðr; s^0 Þ ds^ 0 ðvðr; s^ Þ þ ds^ rvðr; s^ ÞÞ ds^ dr s O Sn1 Sn1 Z Z ^ f0 ðr; s^ Þvðr; s^ Þ ds^ dS ðs^ nÞ (17) ¼ qO
Sn1
^ denote the positive and negative parts of ðs^ nÞ. ^ where f0 ðr; s^ Þ is the boundary source as in Eq. (5) and ðs^ nÞ In the FEM, the solution fðr; s^ Þ of the variational formulation (17) is approximated in a piecewise linear basis
fðr; s^ Þ fh ðr; s^ Þ ¼
Nn X Na X
fi‘ ci ðrÞc‘ ðs^ Þ
(18)
i¼1 ‘¼1
where ci ðrÞ and c‘ ðs^ Þ are the nodal basis functions of the spatial and angular FE meshes, respectively, fi‘ is the radiance in spatial nodal point i and ‘‘direction’’ ‘; Nn is the number of spatial nodes, and Na is the number of angular directions. The choice of the basis functions, Eq. (18), allows the separation of spatial and angular RTE discretizations, and thus the spatial and angular integrals can be calculated separately and connected using Kronecker products. Therefore, the approach is flexible, e.g. the amount of angular directions can be increased or decreased easily. Inserting approximation (18) into Eq. (17) and choosing the nodal basis functions cj ðrÞ and cm ðs^ Þ as the test functions, the FE-approximation of the RTE can be derived. The FE-approximation of the RTE with the streamline diffusion modification can be written in the form Af ¼ b
(19)
where A ¼ A0 þ A1 þ A2 þ A3 þ A4
(20)
b ¼ b1 q0
(21)
where vector q ¼ ðq01;1 ; ; q0Nn ;Na ÞT 2 RNn Na is the source strength vector 0 j qO such that qi;‘ is the source power at the node i in the ‘‘direction’’ N n Na 0
e
C
having non-zero elements for spatial nodes r i 2 ‘, and f ¼ ðf1;1 ; . . . ; f1;Na ; . . . ; fNn ;1 . . . ; fNn ;Na ÞT 2 is the radiance. The components of Eqs. (20) and (21) are of the form Z Z io ci ðrÞcj ðrÞ dr c‘ ðs^ Þcm ðs^ Þ ds^ A0 ðh; pÞ ¼ c O Sn1 Z Z þ d (22) s^ rcj ðrÞcm ðs^ Þc‘ ðs^ Þ ds^ ci ðrÞ dr O
Sn1
ARTICLE IN PRESS 2772
T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
A1 ðh; pÞ ¼
Z Z n1
ZO SZ
d
þ O
A2 ðh; pÞ ¼
Z qO
A3 ðh; pÞ ¼
s^ rcj ðrÞcm ðs^ Þc‘ ðs^ Þ ds^ ci ðrÞ dr ðs^ rcj ðrÞcm ðs^ ÞÞðs^ rci ðrÞc‘ ðs^ ÞÞ ds^ dr
Sn1
ci ðrÞcj ðrÞ dS
Z Sn1
Z ðms þ ma Þci ðrÞcj ðrÞ dr c‘ ðs^ Þcm ðs^ Þ ds^ O Sn1 Z Z s^ rcj ðrÞcm ðs^ Þc‘ ðs^ Þ ds^ ci ðrÞ dr dðms þ ma Þ þ
(24)
Z
Sn1
O
A4 ðh; pÞ ¼
Z Z
O
O
b1 ðh; pÞ ¼
^ þ c‘ ðs^ Þcm ðs^ Þ ds^ ðs^ nÞ
(23)
Z qO
ms ci ðrÞcj ðrÞ dr dms
Z n1
S
Z
Z Sn1
Sn1
s^ rcj ðrÞcm ðs^ Þ
ci ðrÞcj ðrÞ dS
Z Sn1
(25)
Yðs^ s^ 0 Þc‘ ðs^ 0 Þ ds^ 0 cm ðs^ Þ ds^ Z Sn1
Yðs^ s^ 0 Þc‘ ðs^ 0 Þ ds^ 0 ds^ ci ðrÞ dr
^ c‘ ðs^ Þcm ðs^ Þ ds^ ðs^ nÞ
(26)
(27)
where h ¼ N a ðj 1Þ þ m; p ¼ Na ði 1Þ þ ‘ (j; i ¼ 1; . . . ; Nn , m; ‘ ¼ 1; . . . ; Na , h; p ¼ 1; . . . ; Nn Na ). The FE-solution to the RTE is obtained using Eq. (19). Thus, as a solution to the forward problem, the radiance f in nodal points of the spatial and angular discretizations is obtained. The exitance, which is the measurable quantity, can be calculated using Eq. (6) as
Gh ¼ Mf ¼ MA1 b
(28)
where M 2 RmNn Na is a measurement matrix which contains the discretized measurement operator and m is the number of detectors. For detailed derivation of the FE-approximation of the RTE see e.g. [11,38]. 4.2. Jacobian In the Gauss–Newton iteration algorithm (15), the Jacobian matrix J needs to be solved on each iteration round. This can be efficiently done by using an adjoint method. We first consider the analytic Fre´chet derivative of the forward mapping for the RTE given by [1,6] ! Z C f F0 Z ¼ (29) R 0 0 C ð Sn1 Yðs^ s^ Þf ds^ fÞ ds^ Sn1 where
denotes conjugate transpose, and C is the solution to the adjoint RTE Z Yðs^ s^ 0 ÞCðr; s^ 0 Þ ds^ 0 ; r 2 O
io Cðr; s^ Þ s^ rCðr; s^ Þ þ ðms þ ma ÞCðr; s^ Þ ¼ ms c
Sn1
(30)
with boundary condition ^ Z; Cðr; s^ Þ ¼ ðs^ nÞ
^ r 2 qO; s^ no0
(31)
indicating the redistribution of the measurement isotropically into incoming directions. If we consider one measurement position ei and one source position ej , we can derive the sensitivity functions Z rðma Þ ðrÞ ¼ Ci ðr; s^ Þfj ðr; s^ Þ ds^ Sn1
rðms Þ ðrÞ ¼
Z Sn1
Ci ðr; s^ Þfj ðr; s^ Þ ds^
Z
Z Sn1
Sn1
Yðs^ s^ 0 ÞCi ðr; s^ 0 Þ ds^ 0 fj ðr; s^ Þ ds^
(32)
(33)
Let us assume that the domain O is divided into K disjoint sub-domains Ok , and that the absorption and scattering coefficients are expressed in basis
ma ðrÞ
K X
mak wkðma Þ ðrÞ
(34)
msk wkðms Þ ðrÞ
(35)
k¼1
ms ðrÞ
K X k¼1
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
2773
where wk ðrÞ is a characteristic function of the sub-domain Ok . The kth column of the Jacobian J for the absorption or scattering is a vectorized matrix Z ðm Þ jk a ¼ wkðma Þ ðrÞrðma Þ ðrÞ (36) O
ðms Þ
jk
Z ¼ O
wðkms Þ ðrÞrðms Þ ðrÞ
(37)
We should note that this expression for the Jacobian, based on projection of the sensitivity functions into the finite dimensional space spanned by the characteristic functions wk , is not the same as would be obtained by differentiation of the discrete matrix A, because the latter would also include terms involving the streamline diffusion integrals. 5. Results The performance of the RTE based reconstruction method was tested with 2D simulations. In simulations, a highly scattering circular domain O which contained different inclusions was investigated. The radius of the circle was 10 mm, and 16 equally spaced sources and detectors located on the boundary of the domain. The modulation frequency of the input signal was chosen as 100 MHz. The refractive indices of the background medium, inclusion, and the surrounding medium were chosen as nin ¼ 1. Four types of cases were investigated: a domain with a highly absorbing circular blob and a highly scattering circular blob (case 1), a domain with a highly absorbing circular blob and a low-scattering circular blob located both far apart and close together (cases 2 and 3, respectively), and a domain with a highly absorbing blob and a thin low-scattering gap (case 4). The optical properties of the simulation domains are summarized in Table 1. In all of the simulations, the background absorption and scattering coefficients were ma ¼ 0:025 mm1 and ms ¼ 2 mm1 , respectively. Further, in all of the simulations the scattering shape parameter was g ¼ 0:2 throughout the domain. The simulated absorption and scattering distributions are shown on left columns of Figs. 3–6. The data were generated with the RTE using the FEM as described in Section 4.1. The RTE was solved using Eq. (19) where the FE-integrals were calculated as in Eqs. (22)–(27). Further, the exitance was solved using Eq. (28). The spatial FE-discretization of case 1 (two blobs) contained 8589 nodal points and 16 792 triangular elements. The spatial FE-discretization of the case 2 (two blobs far apart) contained 8589 nodal points and 16 792 triangular elements and the FE-discretization of case 3 (two blobs close together) contained 8437 nodal points and 16 488 triangular elements. Further, the FE-discretization of the case 4 (a blob and a gap) contained 8545 nodal points and 16 704 triangular elements. In all of the simulations, the angular discretization contained 16 angular directions. Gaussian distributed noise with standard deviation of 1% of corresponding amplitude and phase was added to the simulated data. 5.1. Error surfaces First, we investigated the effect of the scaling on the minimisation problem. Therefore, we computed the norms (8) where L was a unit matrix, and (14) where L was as in Eq. (12) using the simulated data Z of case 1 (two blobs in highly scattering medium) and case 4 (a highly absorbing blob and a low-scattering gap) and forward solution F calculated with different optical properties. The forward data were generated for all combinations of 10 homogeneous absorption values on the range ma 2 ½0:005; 0:05 mm1 and 10 homogeneous scattering values on the range ms 2 ½0:5; 5 mm1 . The obtained error maps are shown in Figs. 1 and 2 where the objective function (8) (unscaled data) is on the left and the objective function (14) (scaled data) is on the right. As can be seen from Fig. 1, the error map from unscaled data shows a different shape profile than the error map from the scaled data in the case of a highly scattering medium. The error map from the scaled data shows a more clearly defined minimum than the unscaled data. This indicates that the scaling provide a better means to separate absorption and scattering. In the case of a highly absorbing blob and a low-scattering gap (Fig. 2), the error map from the scaled data shows again a more clearly defined minimum than the error map from the unscaled data. However, the minimum is not as clear as it was in the case of a highly scattering medium. This indicates, that the reconstruction from a medium with a low-scattering gap is challenging and that the absorption and scattering may not be that well distinguished. Table 1 The absorption coefficient ma and the scattering coefficient ms of the background medium and the inclusions Background
Case 1 Cases 2–3 Case 4
Blob 1
Blob 2/gap
ma ðmm1 Þ
ms ðmm1 Þ
ma ðmm1 Þ
ms ðmm1 Þ
ma ðmm1 Þ
ms ðmm1 Þ
0.025 0.025 0.025
2 2 2
0.125 0.125 0.125
2 2 2
0.025 0.025 0.01
5 0.02 0.02
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
0.05
0.05
0.04
0.04
0.03
0.03
a
a
2774
0.02
0.02
0.01
0.01 1
2
s
3
4
5
1
2
s
3
4
5
0.05
0.05
0.04
0.04
0.03
0.03
a
a
Fig. 1. Objective function for case 1 (a highly scattering medium) as a function of scattering and absorption. On the left, the minimized functional was (8) where L was a unit matrix, and on the right, the minimized functional was (14).
0.02
0.02
0.01
0.01 1
2
s
3
4
5
1
2
s
3
4
5
Fig. 2. Objective function for case 4 (a highly absorbing blob and a low-scattering gap) as a function of scattering and absorption. On the left, the minimized functional was (8) where L was a unit matrix, and on the right, the minimized functional was (14).
5.2. Reconstructions The inverse problem was solved with an image reconstruction method described in Section 3 using the RTE as light transport model. The reconstruction mesh contained 3644 triangular elements. Note that the reconstruction mesh was coarser than the data generation meshes in order to mitigate against the possibility of inverse crimes. First, an initial guess for absorption and scattering values was estimated by minimising the least squares error between simulated measurement and forward solution obtained with homogeneous absorption and scattering. Then, the absorption and scattering reconstructions were computed using the total variation regularized output least squares scheme by minimising functional (14). The previously solved homogeneous absorption and scattering values were used as an initial point of the reconstruction algorithm. In the image reconstruction procedure, the data were scaled using transformation (9), and the weighting matrix L was chosen as in Eq. (12). In the solution space, we used transformation (13). The minimisation problem was solved with Gauss–Newton algorithm (15) which was equipped with a line search algorithm for the determination of the step length. In all of the cases, convergence was achieved in 15 iterations.
5.2.1. Case 1: a highly scattering medium As a first case, we investigated a situation in which the medium contained a highly absorbing blob and a highly scattering blob. The radii of the blobs were 2 mm. The absorption and scattering of the background medium were, as in all of the simulations, ma ¼ 0:025 mm1 and ms ¼ 2 mm1 . The absorption and scattering within the highly absorbing blob were ma ¼ 0:125 mm1 and ms ¼ 2 mm1 , and the absorption and scattering within the highly scattering blob were
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
2775
Fig. 3. Absorption coefficients (top row) and scattering coefficients (bottom row) within a domain with a highly absorbing blob and a highly scattering blob (case 1). Simulated distributions are on the left column and the reconstructions are on the right column.
ma ¼ 0:025 mm1 and ms ¼ 5 mm1 . The simulated absorption and scattering distributions can be seen on the left column of Fig. 3. The reconstructions can be seen on the right column of Fig. 3. The absorption distributions are on the top row and the scattering distributions are on the bottom row. As it can be seen from Fig. 3, the reconstructions are clear and the inclusions are well distinguished. Some crosstalk between the scattering and absorption appears in the reconstructions. This, however, is typical for OT reconstructions, and although it has been shown that utilising the frequency domain RTE the amount of crosstalk is significantly smaller than with the steady-state RTE, there still exists some crosstalk between the optical parameters [10,15].
5.2.2. Cases 2 and 3: two blobs far apart and close together Then we investigated a situation in which the highly scattering medium contained a highly absorbing blob and a low-scattering blob. Two situations were considered: one in which the blobs were located far apart (case 2) and an another in which the blobs were located close together (case 3). The radii of the blobs were 2 mm as in the previous simulation. The absorption and scattering within the highly absorbing blob were ma ¼ 0:125 mm1 and ms ¼ 2 mm1 , and the absorption and scattering within the low-scattering blob were ma ¼ 0:025 mm1 and ms ¼ 0:02 mm1 . The absorption and scattering values of the background medium and the blobs are summarized in Table 1 and the simulated distributions can be seen on the left columns of Figs. 4 and 5. The reconstructions from case 2 in which the domain contained two blobs far apart can be seen in Fig. 4, and the reconstructions from case 3 in which the domain contained two blobs close together can be seen in Fig. 5. Figures show absorption coefficients on the top row and scattering coefficients on the bottom row. The simulated distributions are on the left column and the reconstructions are on the right column. As it can be seen, both the low-scattering blob and the highly absorbing blob can be well distinguished. Even in the case in which the low-scattering blob is located close to the absorbing blob (Fig. 5), the reconstruction shows a good quality absorption distribution with a clear blob.
5.2.3. Case 4: a highly absorbing blob and a low-scattering gap As a fourth case we investigated a situation in which the highly scattering medium contained a highly absorbing blob and a thin low-scattering and low-absorbing gap. The radius of the blob was 2 mm and the thickness of the gap was 0:5 mm. The absorption and scattering within the highly absorbing blob were ma ¼ 0:125 mm1 and ms ¼ 2 mm1 , and the absorption and scattering within the gap were ma ¼ 0:01 mm1 and ms ¼ 0:02 mm1 . The absorption and scattering values are summarized in Table 1 and the simulated distributions are shown on the left column of Fig. 6.
ARTICLE IN PRESS 2776
T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
Fig. 4. Absorption coefficients (top row) and scattering coefficients (bottom row) within a domain with a highly absorbing blob and a low-scattering blob (case 2). Simulated distributions are on the left column and the reconstructions are on the right column.
Fig. 5. Absorption coefficients (top row) and scattering coefficients (bottom row) within a domain with a highly absorbing blob and a low-scattering blob (case 3). Simulated distributions are on the left column and the reconstructions are on the right column.
The reconstructions are shown on the right column of Fig. 6. As it can be seen from Fig. 6, the RTE can detect the location of the gap but not the absolute absorption and scattering values. The highly absorbing blob can also be detected, although the absolute values of absorption are not as good as in the previous simulations.
ARTICLE IN PRESS T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
2777
Fig. 6. Absorption coefficients (top row) and scattering coefficients (bottom row) within a domain with a highly absorbing blob and a low-absorbing and low-scattering gap (case 4). Simulated distributions are on the left column and the reconstructions are on the right column.
6. Conclusions In this paper, the forward and inverse problem of the frequency domain radiative transfer equation were considered. An image reconstruction method based on a total variation output regularized least squares method was developed. In the approach, the frequency domain RTE was used as the light transport model. The minimisation problem was solved with a Gauss–Newton algorithm which was equipped with a line search algorithm. Furthermore, scaling of the data and solution spaces were applied to improve the numerical stability of the problem. The radiative transfer equation was numerically solved with the FEM. In the FE-solution to the RTE both the spatial and angular discretizations were implemented in piecewise linear bases. Furthermore, the streamline diffusion modification was utilized to improve the numerical stability of the problem. The approach was tested with simulations. Reconstructions from different cases including domains with low-scattering regions were shown. The results show that the frequency domain RTE can be used in optical tomography to produce good quality reconstructions.
Acknowledgements This work was supported by the Academy of Finland (project 122499) and by EPSRC grant EP/E034950/1. References [1] Arridge SR. Optical tomography in medical imaging. Inverse Probl 1999;15:R41–93. [2] Wang L, Jacques SL. Hybrid model of Monte Carlo simulation and diffusion theory for light reflectance by turbid media. J Opt Soc Am A 1993;10(8): 1746–52. [3] Firbank M, Arridge SR, Schweiger M, Delpy DT. An investigation of light transport through scattering bodies with non-scattering regions. Phys Med Biol 1996;41:767–83. [4] Hielscher AH, Alcouffe RE, Barbour RL. Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues. Phys Med Biol 1998;43:1285–302. [5] Okada E, Firbank M, Schweiger M, Arridge SR, Cope M, Delpy DT. Theoretical and experimental investigation of near-infrared light propagation in a model of the adult head. Appl Opt 1997;36(1):21–31. [6] Dorn O. A transport–backtransport method for optical tomography. Inverse Probl 1998;14(5):1107–30. [7] Klose AD, Hielscher AH. Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer. Med Phys 1999;26(8): 1698–707. [8] Abdoulaev G, Hielscher A. Three-dimensional optical tomography with the equation of radiative transfer. JEI 2003;12(4):594–601.
ARTICLE IN PRESS 2778
T. Tarvainen et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 109 (2008) 2767–2778
[9] Aydin E, de Oliveira C, Goddard A. A comparison between transport and diffusion calculations using a finite element-spherical harmonics radiation transport method. Med Phys 2002;29(9):2013–23. [10] Ren K, Bal G, Hielscher A. Frequency domain optical tomography based on the equation of radiative transfer. SIAM J Sci Comput 2006;28(4):1463–89. [11] Tarvainen T, Vauhkonen M, Kolehmainen V, Kaipio J. Hybrid radiative–transfer–diffusion model for optical tomography. Appl Opt 2005;44(6): 876–86. [12] Tarvainen T, Vauhkonen M, Kolehmainen V, Arridge S, Kaipio J. Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and non-scattering regions. Phys Med Biol 2005;50:4913–30. [13] Arridge SR, Dorn O, Kaipio JP, Kolehmainen V, Schweiger M, Tarvainen T, et al. Reconstruction of subdomain boundaries of piecewise constant coefficients of the radiative transfer equation from optical tomography data. Inverse Probl 2006;22:2175–96. [14] Wright S, Schweiger M, Arridge S. Reconstruction in optical tomography using the P N approximations. Meas Sci Technol 2007;18:79–86. [15] Kim HK, Charette A. A sensitivity function-based conjugate gradient method for optical tomography with the frequency-domain equation of radiative transfer. JQSRT 2007;104:24–39. [16] Klose AD, Hielscher AH. Optical tomography with the equation of radiative transfer. Int J Numer Methods Heat Fluid Flow 2008;18(3/4):443–64. [17] Case KM, Zweifel PF. Linear transport theory. Reading: Addison-Wesley; 1967. [18] Ishimaru A. Wave propagation and scattering in random media, vol. 1. New York: Academic Press; 1978. [19] Boman E, Tervo J, Vauhkonen M. Modelling the transport of ionizing radiation using the finite element method. Phys Med Biol 2005;50:265–80. [20] Tervo J, Kolmonen P, Vauhkonen M, Heikkinen L, Kaipio J. A finite-element model of electron transport in radiation therapy and a related inverse problem. Inverse Probl 1999;15:1345–61. [21] Kanschat G. A robust finite element discretization for radiative transfer problems with scattering. East–West J Numer Math 1998;6(4):265–72. [22] Richling S, Meinko¨hn E, Kryzhevoi N, Kanschat G. Radiative transfer with finite elements I. Basic method and tests. Astron Astrophys 2001;380(2):776–88. [23] Lathrop K. Ray effects in discrete ordinates equations. Nucl Sci Eng 1968;32:357–69. [24] Lathrop K. Remedies for ray effects. Nucl Sci Eng 1971;45:255–68. [25] Paulsen KD, Jiang H. Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization. Appl Opt 1996;35(19):3447–58. [26] Kolehmainen V, Vauhkonen M, Kaipio J, Arridge SR. Recovery of piecewise constant coefficients in optical diffusion tomography. Opt Express 2000;7(13):468–80. [27] Tarvainen T, Kolehmainen V, Vauhkonen M, Vanne A, Gibson A, Schweiger M, et al. Computational calibration method for optical tomography. Appl Opt 2005;44(10):1879–88. [28] Schweiger M, Arridge SR, Nissila¨ I. Gauss–Newton method for image reconstruction in diffuse optical tomography. Phys Med Biol 2005;50:2365–86. [29] Arridge SR, Dehghani H, Schweiger M, Okada E. The finite element model for the propagation of light in scattering media: a direct method for domains with nonscattering regions. Med Phys 2000;27(1):252–64. [30] Aydin ED. Three-dimensional photon migration through voidlike regions and channels. Appl Opt 2007;46(34):8272–7. [31] Dehghani H, Arridge SR, Schweiger M, Delpy DT. Optical tomography in the presence of void regions. J Opt Soc Am A 2000;17(9):1659–70. [32] Henyey LG, Greenstein JL. Diffuse radiation in the galaxy. Astrophys. J. 1941;93:70–83. [33] Duderstadt JJ, Martin WR. Transport theory. New York: Wiley; 1979. [34] Kaipio J, Somersalo E. Statistical and computational inverse problems. New York: Springer; 2005. [35] Arridge SR, Hebden JC. Optical imaging in medicine: II. Modelling and reconstruction. Phys Med Biol 1997;42:841–53. [36] Kolehmainen V. Novel approaches to image reconstruction in diffusion tomography. PhD thesis, University of Kuopio, Kuopio, Finland, 2001. [37] Schweiger M, Arridge SR. Application of temporal filters to time resolved data in optical tomography. Phys Med Biol 1999;44:1699–717. [38] Tarvainen T, Vauhkonen M, Kolehmainen V, Kaipio J. Finite element model for the coupled radiative transfer equation and diffusion approximation. Int J Numer Methods Eng 2006;65(3):383–405.