Regularization Techniques for Inverse Problem in DOT Applications

Diffuse optical tomography (DOT) is an emerging diagnostic technique which uses near–infra–red light to investigate the optical coefficients distribution in biological tissues. The surface of the tissue is illuminated by light sources, then the outgoing light is measured by detectors placed at various locations on the surface itself. In order to reconstruct the optical coefficients, a mathematical model of light propagation is employed: such model leads to the minimization of the discrepancy between the detected data and the corresponding theoretical field. Due to severe ill–conditioning, regularization techniques are required: common procedures consider mainly ℓ 1–norm (LASSO) and ℓ 2–norm (Tikhonov) regularization. In the present work we investigate two original approaches in this context: the elastic–net regularization, previously used in machine learning problems, and the Bregman procedure. Numerical experiments are performed on synthetic 2D geometries and data, to evaluate the performance of these approaches. The results show that these techniques are indeed suitable choices for practical applications, where DOT is used as a cheap, first–level and almost real–time screening technique for breast cancer detection.


Introduction
Diffuse optimal tomography (DOT) is a technique that received a growing interest in recent years. It consists in a non-invasive and non-ionizing imaging procedure that enables to investigate absorption and scattering optical properties of biological tissue. It uses near-infra-red light (600-900 nm): in this particular spectral window the absorption properties of water, oxygenated hemoglobin (HbO) and de-oxygenated hemoglobin (Hb) allow the light to penetrate deeper to examine the tissue up to several centimetres in depth [10,24]. The concentration of HbO and Hb produce a significant difference in optical coefficients between healthy and diseased tissue: such difference provides an imaging contrast tool for tissue diagnostic.
The technical arrangement of DOT consists in several detectors recording the light intensity from differently located sources. The light emerging from the sample is collected by an array of photodetectors or a CCD camera, depending on the specific industrial technology (see [10] for more details). The focus of this work is on the inverse problem in DOT, assuming that the image pre-processing and filtering has already been performed by the recording hardware and/or specialized software. The aim is reconstructing the distribution of the tissue optical coefficient. This optical distribution reveals to be the solution of a minimization problem, where the objective functional measures the discrepancy between the observed data and the corresponding optical field given by a light propagation model. The ill-conditioning of the problem, due to scattering properties of the tissue, requires the use of regularization techniques [1,2] in order to obtain a reliable distribution.
Natural choices for the regularization functional are 2 norm and 1 norm. The former, known also as Tikhonov regularization or ridge regression, has been widely used in DOT [25,12,18]. The 1 -norm [27,15,9] promotes sparse solutions: it forces some components on the solution to strongly emerge and at the same time pushes others to very small values or even to zero. In DOT application the 1 norm has been used to improve sparsity and sharp edges detection [26]. On one hand, the 1 -norm is indeed a suitable choice, since it allows to emphasize the contrast between the homogeneous and healthy tissue and pathological formations, usually small and localized. On the other hand, this norm selects only a small number of variables from a group of correlated ones. Moreover, the presence of noise on the recorded data amplifies these aspects: thus, the use of the sole 1 -norm could lead to imprecisions on localizing the inclusions and on reconstructing their shape.
In this work two novel regularization approaches are proposed. The first one is based on the elastic-net regularization [29], which consist in a linear combination of 1 and 2 norms and it was firstly proposed in the context of machine learning problems. This functional selects groups of correlated variables ( 2 term), yielding thus a good localization, and produces a realistic sparse pattern in the solution ( 1 term). Whilst this technique has been largely used in machine learning [13], at the best of our knowledge, use in DOT applications is original [11]: beside the already mentioned 2 [28] and 1 [23] norms, the interested reader can find in the literature other several approaches, such as algebraic (ART, SIRT) and subspaces (TSVD) techniques, TCG methods [17], but the combination of Tikhonov and 1 is a novelty. The second approach in this context is the usage of the inexact Bregman procedure [4,6,30]: such procedure consists in substituting the regularization functional with its Bregman distance, triggering an iterative procedure which usually yields to contrast enhancement [7] and sometimes a minor computational time [4]. Both aspects are crucial: the former satisfy the objective of improving the visual inspection of the different absorption coefficients, while the latter contributes in creating an almost real-time screening technique. Some authors [21] had applied the Linearized Bregman technique, but actually this approach amounts to solve a different problem than the original [30].
The paper is organized as follows. In Section 2 DOT forward and backward problems are presented. Section 3 is dedicated to the regularization approaches. Section 4 is devoted to the numerical experimentation, where the performance of the proposed techniques is evaluated. Section 5 draws the conclusions of the present work.
2. Mathematical models for light propagation 2.1. Forward problem An exhaustive but computationally intensive mathematical model of light propagation in turbid media is represented by the Radiative Transfer Equation (RTE), an integro-differential equation which ensures the conservation of energy of the light radiance ([W/cm 2 sr]). Performing an expansion in spherical harmonics of the radiance provides a series of approximating equations: its truncation to the first N terms yields the so-called P N approximation.
Let r ∈ Ω ⊂ R 3 be the position in the domain Ω (linear unit in [cm]) and S(r) [W/cm 3 ] the isotropic volume source. Considering a P 1 expansion of the RTE and under the additional hypotheses of isotropic sources and slow variation of the photonic flux, the following diffusive equation (DE) is obtained for the fluence rate at steady-state (see [14] for a complete derivation).
where the photon fluence Φ [W/cm 2 ] represents the integral of the radiance over the solid angle where µ s (r) = µ s (r)(1 − g) [1/cm] is the reduced scattering coefficient, µ s being the scattering coefficient, g the anisotropy scattering factor and µ a = µ a (r) [1/cm] is the absorption coefficient, with µ s µ a . From the microscopic viewpoint, in the DE modeling approach, the photons are supposed to move throughout the sample along random paths. Each photon travels along straight segments with a sudden discontinuity when the photon changes direction or is absorbed. The average length of rectilinear tracts l tr is the mean free path, l tr ≈ 1 µs .

Inverse problem
We consider the reduced scattering coefficient µ s to be the same inside the target, so that D ≈ D 0 = const. In this case, the goal of DOT is to reconstruct µ a (r) from a set of boundary measurements of Φ(r), assuming D is known. In the following, the Rytov perturbative approach is adopted. This choice is acceptable in the hope of disposing of fast algorithms to provide "almost real time" diagnostic results. Namely where a series solution for ψ(r) can be derived, ψ 0 (and Φ 0 ) being the so-called background conditions (in absence of perturbations) and ψ 1 the perturbed field induced by the heterogeneities. After some algebra (see [20], Vol.II, Ch.17 for a detailed derivation), the following equation for where µ a (r) = µ a,0 + δµ a (r), δµ a being a perturbation from background conditions. A linearized solution ψ 1 ∼ ψ 1,0 for problem (3) is obtained from G(r − r ) being the Green's function of (3) evaluated in r for a Dirac-δ source located in r . Notice that ψ 1 = log(Φ/Φ 0 ) represents the logarithmic amplitude fluctuation of the light intensity and it is the datum used as a rhs in the inverse problem solution.
The integral relation (4) provides an explicit model solution and is usually dealt with by discretizing the domain into a large number N of volume elements (voxels) and evaluating the rhs of the equation at the detector positions. Using for example a midpoint quadrature rule (voxel centroid r j , voxel volume ∆V j , j = 1, . . . , N ) to discretize Eq. (4), one obtains where the dependence on the source and detector position r s and r d , respectively, has been made explicit. Considering now to dispose of M = N s × N d measurements (that is, M couples source-detector, N s is the number of the sources, N d the number of detectors) the following system of equations naturally arises where ∈ R (M ×1) is obtained from measurements or from numerical simulations of the forward model. The inversion of system (7) is a nontrivial task, since the system can be overdetermined (M > N ) or underdetermined (M < N ) in the case of strong data deficiency and the matrix J is ill conditioned. In the forthcoming sections, we discuss the issue of regularization procedures applied to the solution of (7). To do this, it is important to observe that problem (7)

Elastic Net Penalization
The functional in Eq. (8) is regularized via the elastic-net [29] min The parameter λ is a regularization parameter and α controls the influence of the penalties given by the 1 and 2 norms on the computed solution. More precisely, λ balances the tradeoff between the Least Square functional and the regularization, while α controls the convex combination of the two norms, allowing to decide which one has the greater influence (α = 0 sets the Thikonov regularization, α = 1 leads to 1 regularization). By varying α ∈ (0, 1) sparsity is traded for the inclusion of more solution components contributions. The choice of λ is quite delicate. For pure sparsity-promotion via 1 , there is still not any golden rule for an automatic selection for λ: some specific strategies refer to particular instances [3,22]. For the pure 2 penalization there exist several more consolidated methods for automatic choice of λ [25]. A general possibility, implemented in advanced software packages as the glmnet [16] package, is to use k−fold cross validation techniques to establish optimal values of λ.

Bregman Procedure.
The Bregman procedure [4] finds its basis on the Bregman distance: given a convex function f with a non-empty subdifferential ∂f , its Bregman distance from x to y is D p f (x, y) = f (x) − f (y) − p, x − y , p ∈ ∂f (y), being ∂f (y) the subdifferential of f at y [8,9]. When f ∈ C 1 , ∂f (y) ≡ ∇f (y). In this work, the Least Square functional is coupled with 2 and 1 regularization, triggering the following scheme for the Tikhonov regularization while for LASSO regularization the algorithm reads as where soft γ (x) = sign(x) max(|x| − γ, 0). In [4,6] the convergence proof of the general method is stated and more details are given in [30] for the particular case of Eq. (11).

Numerical Tests
A 2D domain Ω is considered: it consists in a semicircle of radius 5 cm (Fig. 1). This geometry, albeit in a reduced dimensional setting and idealized version, represents a female breast supported on a solid plate on which embedded light sources are alternatively turned on for the DOT exam. The synthetic fluence observations in the detectors are generated by solving (1) with the MatLab R PDE Toolbox on a fine triangular mesh consisting of almost 37858 triangles and using an approximation of degree 2. The boundary conditions are ϕ = 0 on ∂Γ D , ϕ+2 A D ∂ϕ/∂ n = 0 on Γ R , n is the outward normal vector and the coefficient A keeps into account the different tissue and air refractive indices [14]. The first condition represents the fact that light does not escape from the solid plate. N s = 20 point-wise light sources are simulated, uniformly located on the horizontal boundary at a distance of 1 mm from Γ D , N d = 200 detectors are radially disposed 1 mm inside the domain along the boundary Γ R . We discretize the integral using N = 3882 voxels. The background values are set to µ s = 0.1 cm −1 , µ a,0 = 0.01 cm −1 for all the tests, and Gaussian noise at 1% level is added to the logarithmic fluctuation to simulate uncertainties in the data and the real behaviour of technical equipment. Three different cases are generated A) a single circular inclusion is placed at (x, y) = (0, 2) cm with radius .3 cm. The absorption coefficient of the inclusion is set to µ a,inc = 2µ a,0 B) two circular inclusions are placed at (x, y) = (±2, 2) cm both with radius 0.3 cm. The absorption coefficient of the inclusions is set to µ a,inc = 2µ a,0 C) as in case B) but the absorption coefficient of the left inclusion is set to µ a,inc = 5µ a,0 , while the right inclusion has µ a,inc = 2µ a,0 .
Tests are performed using the elastic-net regularization method implemented in the glmnet software package [16] and alternatively using the Bregman procedure. The results are depicted in Fig. 2, where the two approaches are compared with the state-of-the-art techniques.
As for the elastic-net approach, it allows to correctly locate the inclusion(s) with their absorption coefficient for tests A) and B). In test C) the localization is less precise, but the expected values of the absorption coefficient are well estimated. As for the Bregman procedure, choosing the ridge regression reveals again to be not a very suitable option: although this procedure increases the quality of the reconstruction by the contrast enhancement ( [4,5]), the inclusion is not clearly identified and the recovered halo can be confused with some artefacts. On the other hand, setting the LASSO regularization really benefits from the Bregman implementation: the results for test (A) show that the location of the inclusion is properly recognized and its absorption value µ a is correctly reconstructed. The reconstruction of test B)  Figure 2. Results for test cases A, B, C in column-wise order. Top to bottom: Tikhonov regularization, Bregman procedure coupled with 2 -norm, pure LASSO approach, Bregman coupled with 1 -norm and eventually elastic-net procedure with α = 0.5.
suffers from a small displacement, but the value of µ a is well recovered. Finally, in test C) both inclusions are not completely placed in the correct location, but the absorption value is close to the correct one. This misplacement is due to the different values of absorption coefficient, hence the regularization parameter promotes the higher one.

Conclusions
This work has investigated two original approaches to regularize the inverse problem arising in DOT reconstruction: the elastic-net regularization and the Bregman procedure. The experimental results attest that both techniques actually improve the quality of the estimation of the distribution coefficient, both in localizing the inclusion and in recovering its correct value. Future work will involve the combination of the two procedure and the employment of local and adaptive regularization techniques.