# Space-Dependent Sobolev Gradients as a Regularization for Inverse Radiative Transfer Problems

Authors: Y. Favennec; F. Dubot; D. Le Hardy; B. Rousseau; D. R. Rousse
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Diffuse optical tomography problems rely on the solution of an optimization problem for which the dimension of the parameter space is usually large. Thus, gradient-type optimizers are likely to be used, such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, along with the adjoint-state method to compute the cost function gradient. Usually, theL
2-inner product is chosen within the extraction procedure (i.e., in the definition of the relationship between the cost function gradient and the directional derivative of the cost function) while alternative inner products that act as regularization can be used. This paper presents some results based on space-dependent Sobolev inner products and shows that this method acts as an efficient low-pass filter on the cost function gradient. Numerical results indicate that the use of Sobolev gradients can be particularly attractive in the context of inverse problems, particularly because of the simplicity of this regularization, since a single additional diffusion equation is to be solved, and also because the quality of the solution is smoothly varying with respect to the regularization parameter.


---

## Body

## 1. Introduction

Diffuse optical tomography (DOT) consists in reconstructing the spatial radiative property distributions of an optically thick participating (absorbing and scattering) medium,D, from optical measurements collected at the border of D [1–3]. This inverse problem has been a surge of interest since the middle of the nineties, due to better theoretical understanding in radiative transfer; improvements in light detection and control of light sources; and the increase of the computing resource performances. This optical diagnostic technique is mainly driven by potential applications in biological tissues imaging although it can be used for nondestructive testing in industrial processes [4, 5] or the characterization of complex semitransparent media such as cellular foams with spatially varying porosities or fiber reinforced thermoplastic composites with inhomogeneity in the fibers distribution. In this domain, the characterization has often been limited to the reconstruction of constant absorption and scattering coefficients and anisotropy factor of the phase function, considering the unidimensional radiative transfer equation (RTE) as forward model [6].Biomedical applications of the DOT include imaging of finger joints [7], breast imaging [8, 9], neuroimaging of humans [10–12], or small-animal imaging [13] to name but a few. In this context, near-infrared light is used in order to take advantage of the so-called therapeutic window (600–900 nm) in which tissues have low absorption. According to the application, RTE has to be considered as forward model [14–19], diffuse approximation (DA) may be applied [8, 12, 20–23], or hybrid model is employed [24–28]. This paper deals with two very different forward models, namely, the bidimensional steady-state RTE and the three-dimensional frequency DA, in order to gauge the potential ability of the introduced regularization with the different models.The evaluation of space-distributed physical properties from the knowledge of measurements is an inverse problem as opposed to forward problems where the state is computed from the knowledge of properties, boundary conditions, and so forth. Such inverse problems are usually difficult to solve due to their ill-posed nature [29].The solution of the inverse problem is obtained through the minimization of a cost function that measures the misfit between the predictions (solution of the forward problem) and some experimental measurements. Though zero-order optimization methods may be used for relatively similar problems [31, 32], gradient-type methods are dealt with in this paper because of the large dimension of the parameter space. Due to the ill-posed nature of the inverse problem, unstable solutions may be obtained if no regularization is performed.In order to cope with such instabilities, the ordinary Tikhonov approach consists in assuming that the properties are subject to a minimal norm in a given vectorial space which leads to enhancing thea priori regularity of the properties [33, 34]. This approach relies on the penalization of the cost function adding prior information on the properties to be recovered. However, the weight associated with the penalization function must be chosen carefully. Indeed, a too small value for the penalization function leads to a small cost function but with high property variations, while too high values lead to small property variations about priors but with a high resulting cost function value at the optimum. The compromise minimizing both the cost function and the property variations may be performed according to the construction of the so-called L-curve as suggested by Hansen and O’Leary [35]. Even though, for appropriate penalization, the optimal Tikhonov parameter has been found to be quasi-independent of the parameter space dimension [36], the construction of such curve, searching its convex corner location, is likely to necessitate a large number of inverse problem solutions for plenty of penalization function weights. Furthermore, it is now well known that the Tikhonov regularization is to be used only together with matrix-inversion-based optimizers, but not with BFGS-like optimizers because these do not rely on matrix inversion [19, 23].Appropriate parameterization for the space-dependent properties may also affect the solution regularity. As a matter of fact, the amount of measurements data is likely to be insufficient to retrieve the large amount of parameters when considering, at the end, the discrete optimization problem on a fine grid. Thus, regularization by reduction of the dimension of the parameter space could be employed, which consists in searching the physical properties projected on a coarse finite element grid, as suggested for instance in [37], and applied for optical tomography applications in [23].This paper focuses on a different strategy to improve the reconstruction of space-distributed parameters when the data is corrupted by noise. The proposed strategy filters the effect of the noise, not on the data as it is done for instance when using the mollification method [29, 38], but on the cost function gradient itself. This means that the optimization is not modified at all since the original cost function remains the same. It also means that the data is not modified: the whole information contained within the experimental data is integrated within the cost function to be minimized. As a matter of fact, several inner products are to be defined when considering the inverse optimization problem. According to Protas [39], each of them (cost function definition, adjoint identity, and gradient extraction after defining the space used in the definition relationship between the cost function derivative and its gradient) may somehow yield to a particular regularization opportunity. The noise present in the measurements data propagates through the adjoint problem and then to the cost function gradient from which the directions of descent are computed to supply the optimizer. The propagation of the noise induces high variations of the cost function gradient with respect to what the gradient would be without noise. The idea is thus to somehow apply a low-pass filter to the cost function gradient, exactly where the noise is the highest. This is done choosing the space-dependent weighted Sobolev inner product definition when extracting the cost function gradient rather than the ordinary Hilbert inner product. Note that space-independent Sobolev inner products have been used in other situations, in image treatment (e.g., [40]), in partial differential equations (PDEs), in energy minimization (e.g., [41]), and in boundary conditions reconstruction with accurate final time [42, 43].The underlying application of this study is optical tomography, that is, the reconstruction of radiative properties from light intensity measurements on some parts of the boundary. The paper is therefore organized as follows. Section2 gives the theoretical bases for extracting the cost function gradients when the cost function is expressed in terms of a state which is solution of some PDEs. This section is written in a concise way, following [44], such that the results can be easily transposed to other practical applications. Also, only the continuous equations are dealt with, for the states, the adjoints, and the gradients. This “differentiate-then-discretize” approach has the advantage of versatility when considering afterwards different discretization choices (one for the states, one for the physical properties, etc.). The cost function gradient is then extracted firstly in a general framework using the ordinary L
2-inner product and then using space-independent and finally space-dependent Sobolev inner products. Section 3 presents application of space-dependent Sobolev gradients for the inverse problem of optical tomography in which the forward model is the RTE. The continuous cost function gradient is fully derived and numerical tests are performed for several levels of noise. As the Sobolev filter is chosen to follow an exponential function, only one parameter fully describes the regularization. Next, in Section 4, a very different model, namely, the frequency DA in the three-dimensional case, is considered as forward model of light propagation. For the two models, the quality of the reconstructions is measured quantitatively with the aid of both deviation and correlation factors. Without contest, reconstructions are greatly improved with the space-dependent Sobolev gradients. Finally, Section 5 naturally summarizes the main conclusions and perspectives of the current work.


## 2. Space-Dependent Sobolev Cost Function Gradient

Let us define the domain of interestD
⊂
R
d and denote the cost function J
(
u
) expressed in terms of the state u
:
D
→
K. This cost function may be expressed as the half of the squared norm of a quadratic function v based on the discrepancy between some state predictions u and the related measurements u
d. The norm is expressed with a complex inner product (
·
,
·
)
X defined later on for specific applications. Cost function J
(
u
) is actually to be minimized with respect to some space-distributed parameters ψ
∈
P, with P
⊆
L
2
(
D
), for example (depending on the forward model considered). The generic cost function is given by(1)
j
ψ
≔
J
u
=
1
2
f
u
-
u
d
X
2
.The stateu
(
ψ
) is related to the parameters ψ through an operator (that is likely to be nonlinear) that combines the partial differential equations along with the boundary conditions, initial condition, and so forth. This operator is denoted as S for the state problem. To be concise, one writes down(2)
S
u
,
ψ
=
0
.In order to use efficient optimization algorithms for large-scale inverse problems encountered in distributed parameter estimation, the cost function derivative or rather the cost function gradient has to be computed. To do so, the definition of the directional derivative ofj
(
ψ
) towards the direction ϕ is used.Definition 1 (directional derivative).
Let a pointψ and a direction ϕ
∈
P. The directional derivative of j at point ψ in the direction ϕ is [45](3)
j
′
ψ
;
ϕ
≔
lim
t
→
0
+
⁡
j
ψ
+
t
ϕ
-
j
ψ
t
.The application of the directional derivative (3) on the cost function (1) leads to the equality(4)
j
′
ψ
;
ϕ
=
f
u
-
u
d
,
f
′
u
u
′
ψ
;
ϕ
X
∀
ϕ
∈
Pand, because of the linearity of j
′
(
ψ
;
ϕ
) with respect to ϕ,(5)
j
′
ψ
;
ϕ
=
∇
j
ψ
,
ϕ
Z
∀
ϕ
∈
P
.The inner product·
,
·
Z involved in the right-hand side of (5) will be also defined later on, when necessary. That is, in this inner product, regularization should be introduced.The differentiation of the state problem in the directionϕ yields the directional derivative u
′
(
ψ
;
ϕ
) that can further be injected into the inner product (4). The state equation being given by (2), the derivative of S
(
u
,
ψ
) towards the direction ϕ is written as(6)
S
u
′
u
,
ψ
u
′
ψ
;
ϕ
+
S
ψ
′
u
,
ψ
ϕ
=
0
∀
ϕ
∈
P
,where S
u
′ is the tangent operator, that is, the Jacobian of S with respect to u, and S
ψ
′ is the Jacobian of S with respect to ψ. Using this basic approach, one needs dim
⁡
P integrations of (6) to access the full gradient ∇
j
(
ψ
)
∈
R
dim
⁡
P for a fully discretized problem. When dim
⁡
P is large, which is likely to be the case when recovering space-distributed parameters, this approach becomes inefficient even though the computation of the linear tangent operator S
u
′
(
u
,
ψ
) can be reused for all directions ϕ of the canonical basis of ψ.Following Lions [46] and Céa [47], the adjoint-state method computes the whole cost function gradient through the integration of a single additional adjoint-state problem. A new variable (the adjoint-state variable u
∗
(
ψ
)) is then introduced so that the cost function derivative also satisfies the relationship(7)
j
′
ψ
;
ϕ
=
S
ψ
′
u
,
ψ
ϕ
,
u
∗
Y
∀
ϕ
∈
P
,where, again, the inner product ·
,
·
Y should be defined later on when necessary. Combining (6) and (7), the cost derivative is rewritten as(8)
j
′
ψ
;
ϕ
=
-
S
u
′
u
,
ψ
u
′
,
u
∗
Y
∀
ϕ
∈
P
.The adjoint-state problem is then identified through (4) and (8), such that it satisfies the equality(9)
-
S
u
′
u
,
ψ
u
′
,
u
∗
Y
=
f
u
-
u
d
,
f
′
u
u
′
ψ
;
ϕ
X
∀
ϕ
∈
P
.Next, if the adjoint problem (9) is satisfied (this means that the adjoint-state u
∗ is accessible), then the cost function gradient is given by the inner product (7). In practice, operators involved in (9) are transposed. The inner product property S
u
′
u
,
ψ
u
′
,
u
∗
Y
=
S
∗
(
u
,
ψ
)
u
∗
,
u
′
Y (with S
∗ being the transposed conjugate operator of S
u
′) gives the adjoint equation concisely written as(10)
S
∗
u
,
ψ
u
∗
+
f
′
u
f
u
-
u
d
=
0
∀
ϕ
∈
P
.As reported by Protas [39], the gradient is identified in a given space for which the metric had been selected. To do so, the two following definitions are needed.Definition 2 (L
2
(
D
) Hilbert space).
Consider(11)
L
2
D
=
u
:
D
⟶
R
∣
∫
D
u
x
u
x
d
x
<
∞
.This functional space is associated with the inner product:(12)
u
,
w
L
2
D
=
∫
D
u
x
w
x
d
x
.Definition 3 (H
1
(
D
) Sobolev space).
Consider(13)
H
1
D
≔
ψ
:
D
⟶
R
∣
ψ
∈
L
2
D
;
∇
ψ
∈
L
2
D
.This functional space is associated with the inner product:(14)
ψ
,
ϕ
H
1
D
≔
ψ
,
ϕ
L
2
D
+
∇
ψ
,
∇
ϕ
L
2
D
.Most often, theL
2
(
D
) real inner product for ·
,
·
Z is used to extract the gradient. Combining (5) and (7) with the inner product of Definition 2 gives, after transposing within the right-hand side inner product, the following result:(15)
∇
j
ψ
,
ϕ
L
2
D
=
S
ψ
′
u
,
ψ
u
∗
,
ϕ
Y
∀
ϕ
∈
L
2
D
.This way for writing down the cost function gradient is the ordinary one as first suggested by Lions [46] and Céa [47]. However, according to Protas [39], it is suggested to use other inner product definitions due to the poor scaling of the corresponding discrete optimization problem. As reported in [48], the incorporation of such derivatives into the inner product, instead of using the L
2 inner product, has the effect of scale-dependent filtering and allows one to extract smoother gradients, thereby preconditioning the optimization process. In some other studies, solving partial differential equations through optimization (see, e.g., the articles of Mahavier [49], Danaila and Kazemi [41], and Majid and Sial [50]), in image treatment (see, e.g., Renka [40]), in boundary conditions reconstruction (see, e.g., Nassiopoulos and Bourquin [42] and Bourquin and Nassiopoulos [43]), or in initial state reconstruction [39], the Sobolev space H
1
(
D
) is used for extracting the cost function gradient. This acts as a preconditioner modifying the direction of descent involved in the gradient-type optimization algorithm. It is shown, for instance, in the study of Danaila and Kazemi [41] that the cost function decrease when solving the PDEs problem can be much higher with the H
1 inner product than with the ordinary L
2 one.In the present study, the choice of another inner product comes from the fact that we are faced with dealing with noisy measurementsu
d. As a matter of fact, the measurements are noisy by nature. These measures are involved in the cost function definition (1), through a difference with the predictions. The noisy measurements are thus also involved in the forcing term f
′
(
u
)
f
(
u
-
u
d
). The noise then propagates through the adjoint equation (10), and then to the cost function gradient (15).In order to deal with the smoothing of the measurements noise that propagates in the adjoint system and then to the cost function gradient, the weighted Sobolev space is introduced.Definition 4 (H
1
l
(
D
) weighted Sobolev space inner product).
Consider(16)
ψ
,
ϕ
H
1
l
D
≔
ψ
,
ϕ
L
2
D
+
l
∇
ψ
,
l
∇
ϕ
L
2
D
.Remark 5.
l
∈
R
+ is the regularization parameter.This strategy has been applied in bidimensional optical tomography applications based on the radiative transfer equation [19] where it was shown how the diffusion of the noise could attenuate fluctuations and thus give better reconstructions with a moderate parameter l.The strategy that is developed here goes much further: it consists in filtering the cost function gradient where, and only where it is needed, the high frequency fluctuations exist, at the vicinity of the sensors.In applications considered here, that is, in optical tomography, sensors are located on the boundaries, so that the effect of the noise appears on the boundaries before being diffused through the adjoint-state equation. Thus, the idea is to choose a filtering function whose value is high on the boundary and that continuously decreases within the medium. To do so, one uses the distance function definition.Definition 6 (distance function).
Letx
∈
D and y
∈
∂
D; d
(
x
,
∂
D
) is the Euclidean distance function to the set ∂
D if(17)
d
x
,
∂
D
=
min
y
∈
∂
D
⁡
x
-
y
.Remark 7.
As an example, the filtering function can be written as(18)
l
v
2
x
=
exp
⁡
-
v
d
2
x
,
∂
D
.This makes the construction of the space-dependent weighted Sobolev space possible.Definition 8 (H
1
(
v
)
(
D
) space-dependent weighted Sobolev space inner product).
Consider(19)
ψ
,
ϕ
H
1
v
D
≔
ψ
,
ϕ
L
2
D
+
l
v
x
∇
ψ
,
l
v
x
∇
ϕ
L
2
D
.Remark 9.
With the choice of the filtering function (18), the regularization parameter becomes v.Using theH
1
(
v
)
(
D
) inner product for ·
,
·
Z and if ϕ
∈
H
1
(
D
), then there exists a unique ∇
H
1
(
v
)
j
(
ψ
)
∈
H
1
(
v
)
(
D
) such that(20)
j
′
ψ
;
ϕ
=
∇
H
1
v
j
ψ
,
ϕ
H
1
v
D
∀
ϕ
∈
H
1
D
.Taking into account that(21)
∇
∇
H
1
v
j
ψ
,
∇
ϕ
L
2
D
=
-
Δ
∇
H
1
v
j
ψ
,
ϕ
L
2
D
+
∇
∇
H
1
v
j
ψ
·
n
,
ϕ
L
2
∂
D
,one obtains by identifying terms(22)
1
-
l
v
2
Δ
∇
H
1
v
j
ψ
,
ϕ
L
2
D
=
S
ψ
′
u
,
ψ
u
∗
,
ϕ
Y
∀
ϕ
∈
H
1
Dwith a null flux boundary condition; that is, ∇
∇
H
1
(
l
)
j
(
ψ
)
·
n
=
0 on ∂
D. At this stage, the application being defined (measurements, state model, etc.), the Y-norm is given accordingly, and the cost function gradient can be extracted.Remark 10.
Due to the inclusionH
1
(
v
)
(
D
)
⊂
L
2
(
D
), the space-dependent weighted Sobolev gradient ∇
H
1
(
v
)
j
(
ψ
) is indeed an acceptable direction of descent for the optimization problem that consists in minimizing (1).


## 3. Application on 2D RTE-Based Optical Tomography

### 3.1. Inverse Problem Statement

The inverse problem consists mathematically in minimizing a cost function which quantitatively measures the discrepancy between some optical measurements and related predictions over a set of sources, sensors, and angles of directional integration. In optical tomography, taking into account the reflection at the interface and assuming that the medium is convex, the measurable quantity is the exitance, or portion of it in a given solid angleω
(
s
)
⊆
s
·
n
>
0 centered on direction s. The related prediction can be written as (23)
E
x
,
s
=
∫
ω
s
1
-
R
x
,
s
I
x
,
s
s
·
n
d
ω
s
,where I is the radiative intensity and the function R
(
x
,
s
)
∈
[
0,1
] is the Fresnel reflection factor. In what follows in this section, the border of the medium is assumed to be transparent; that is, R
=
0. Moreover, the solid angle ω
(
s
) is assumed to be sufficiently small so that the approximation E
(
x
,
s
r
)
≈
ω
r
I
(
x
,
s
r
)
s
r
·
n holds. Therefore, the radiative intensity measurements I
˘
(
x
,
s
r
) can be assessed and the cost function (24) mapping from P
=
L
∞
(
D
)
×
L
∞
(
D
) to R uses, according to the schematic description of the problem (Figure 1), the specific Euclidean norm:(24)
j
κ
,
σ
≔
J
I
=
1
2
s
·
n
I
˘
I
-
I
˘
X
2
with
u
X
≜
∑
k
=
1
K
∑
d
=
1
D
∑
r
=
1
R
u
k
x
d
,
s
r
2
1
/
2
,where the index k refers to the source (test) number, the index d refers to the detection node, and the index r refers to the direction. In (24), J
(
I
) is the cost function explicitly expressed in terms of the prediction I, while j
(
κ
,
σ
) is the cost function implicitly expressed in terms of the parameters. This latter function is often called the reduced cost function. According to (1), one has the straightforward correspondence f
(
u
-
u
d
)
↔
s
·
n
/
I
˘
(
I
-
I
˘
). Note that the scalar s
·
n in (24) is introduced in order to avoid infinite Dirichlet boundary condition in the adjoint-state problem. For a given source configuration, say the kth one, the prediction I
k
(
x
d
,
s
r
) at the sensor d and direction r is given through the solution of the forward steady-state radiative transfer problem formulated as(25)
Find
I
k
x
,
s
such  that
∀
s
∈
S
n
-
1
,
∀
k
=
1
,
…
,
K
:
s
·
∇
+
κ
+
σ
I
k
x
,
s
=
σ
∮
S
n
-
1
I
k
x
,
s
Φ
s
,
s
′
d
ω
s
′
,
∀
x
∈
D
,
I
k
x
,
s
=
I
k
0
x
×
1
x
∈
∂
D
k
0
×
1
s
=
s
k
0
,
∀
x
∈
∂
D
,where s is the considered direction of propagation, Φ
(
s
,
s
′
) is the phase function representing the probability that a photon arriving from the direction s
′ is scattered to the direction s, κ and σ are the absorption and scattering space-dependent functions, respectively, and 1
[
·
] denotes the indicator function.Figure 1
Schematic description of the experiment. Only one source is on at once, while others are off. For each source configuration numberk
=
1
,
…
,
K, the emerging intensity is measured on all sensors d
=
1
,
…
,
D in different directions r
=
1
,
…
,
R. This is a theoretical extension of the experimental directional-directional spectral transmittances device using a Fourier transform infrared spectrometer described in [30].


### 3.2. Cost Function Gradient Derivation

The computation of the cost function gradient follows the generic methodology detailed in the previous section. The cost function derivative towards the directionϕ is given differentiating (24):(26)
j
′
κ
,
σ
;
ϕ
=
∇
j
κ
,
σ
,
ϕ
Z
=
s
·
n
I
˘
2
I
-
I
˘
,
I
′
X
∀
ϕ
∈
P
.In order to derive the cost gradient for this specific application,Y-based inner products are defined.Definition 11.
Consider(27)
u
,
w
Y
≜
∑
k
∮
S
n
-
1
∫
D
u
w
d
x
d
ω
s
,
u
,
w
∂
Y
≜
∑
k
∮
S
n
-
1
∫
∂
D
u
w
d
x
d
ω
s
.TheY-inner product satisfies Properties 1 and 2, both being useful for the cost function gradient derivation.Property 1.
Consider(28)
s
·
∇
u
,
w
Y
=
-
s
·
∇
w
,
u
Y
+
s
·
n
w
,
u
∂
Y
.Property 2.
Consider(29)
σ
∮
S
n
-
1
u
x
,
s
′
ϕ
s
,
s
′
d
ω
s
′
,
w
x
,
s
Y
=
∑
k
∮
S
n
-
1
∫
D
σ
∮
S
n
-
1
u
x
,
s
′
ϕ
s
,
s
′
d
ω
s
′
w
x
,
s
d
x
d
ω
s
=
∑
k
∫
D
σ
∮
S
n
-
1
∮
S
n
-
1
w
x
,
s
u
x
,
s
′
ϕ
s
,
s
′
d
ω
s
′
d
ω
s
d
x
=
∑
k
∫
D
∮
S
n
-
1
σ
∮
S
n
-
1
w
x
,
s
ϕ
s
,
s
′
d
ω
s
u
x
,
s
′
d
x
d
ω
s
′
=
σ
∮
S
n
-
1
w
x
,
s
ϕ
s
,
s
′
d
ω
s
′
,
u
x
,
s
′
Y
.The generic equation for the cost gradient (15) applied on both space-distributed parameters κ and σ gives, for direction ϕ
=
ϕ
κ
;
0
∈
P,(30)
∇
j
κ
,
σ
,
ϕ
Z
=
I
x
,
s
I
∗
x
,
s
,
ϕ
Yand, for direction ϕ
=
0
;
ϕ
σ
∈
P,(31)
∇
j
κ
,
σ
,
ϕ
Z
=
I
x
,
s
I
∗
x
,
s
-
∮
S
n
-
1
I
x
,
s
ϕ
s
,
s
′
d
ω
s
′
I
∗
x
,
s
,
ϕ
Y
.In these two equations, the directionϕ
∈
P is split into two terms ϕ
κ and ϕ
σ in order to be able to compute the cost function gradient relatively to each of the functions κ and σ separately.On one hand, (30) is rewritten such that(32)
∇
j
κ
,
σ
,
ϕ
Z
=
ϕ
I
x
,
s
,
I
∗
x
,
s
Yand, on the other hand, Property 2 enables rewriting (31) such that(33)
∇
j
κ
,
σ
,
ϕ
Z
=
ϕ
I
x
,
s
-
ϕ
∮
S
n
-
1
I
x
,
s
ϕ
s
,
s
′
d
ω
s
′
,
I
∗
x
,
s
Y
.Besides, the differentiation of the forward model enables rewriting previous equations too, forϕ being either ϕ
κ
;
0 or 0
;
ϕ
σ:(34)
∇
j
κ
,
σ
,
ϕ
Z
=
-
s
·
∇
+
κ
+
σ
I
′
x
,
s
,
I
∗
Y
+
σ
∮
S
n
-
1
I
′
x
,
s
ϕ
s
,
s
′
d
ω
s
′
,
I
∗
Y
.Properties1 and 2 enable rewriting (34) such that(35)
∇
j
κ
,
σ
,
ϕ
Z
=
s
·
∇
-
κ
-
σ
I
∗
x
,
s
,
I
′
Y
-
s
·
n
I
∗
,
I
′
∂
Y
+
σ
∮
S
n
-
1
I
∗
x
,
s
ϕ
s
,
s
′
d
ω
s
′
,
I
′
Y
.The adjoint problem is eventually obtained identifying (26) with (35):(36)
Find
I
k
x
,
s
such  that
∀
s
∈
S
n
-
1
,
∀
k
=
1
,
…
,
k
:
-
s
·
∇
+
κ
+
σ
I
k
∗
x
,
s
=
σ
∮
S
n
-
1
I
k
∗
x
,
s
ϕ
s
,
s
′
d
ω
s
′
,
∀
x
∈
D
,
-
s
·
n
I
k
∗
x
,
s
=
s
·
n
I
˘
k
2
I
k
-
I
˘
k
x
,
s
×
1
x
∈
∂
D
d
x
,
∀
x
∈
∂
D
,where ∂
D
d denotes the set of the pointwise sensors.


### 3.3. Numerical Results

The optimization relies on the BFGS algorithm described by Liu and Nocedal [51] coupled with a quadratic interpolation line search [52]. Generally speaking, this latter has proved to be particularly effective in these recent years in the area of inversion (see [53–55] to name but a few) due to its effectiveness in minimizing nonlinear problems [56].The two-dimensional domain is a square of 2 cm length. Four collimated sources are located at the center of each side, while sensors are placed around sources towards twenty directions. The target properties to be reconstructed,κ
t and σ
t, are defined by the background for which κ
=
0.25 cm−1 and σ
=
20 cm−1 and two square inclusions, the former for which κ
=
0.15 cm−1 and σ
=
10 cm−1 and the latter for which κ
=
0.35 cm−1 and σ
=
30 cm−1. Both the forward and the adjoint models are solved using P
1 discontinuous Galerkin finite elements for which the mesh convergence has been validated. Three distinct meshes have been created, one in order to build the synthetic data before adding random noise following Klose [57] (3 721 nodes), one in order to project the forward and adjoint-state variables (1 681 nodes), and one associated with the property functions with the aim of implicitly regularizing the inverse problem (441 nodes). Seeking to gauge the accuracy of the reconstructions, the deviation factor e
1
∈
[
0
,
+
∞
[ and the correlation factor e
2
∈
[
-
1
,
+
1
] are used [58]. The former factor measures the deviation of the reconstructed image from the target image: a small value indicates a reconstructed image with high accuracy. The latter factor measures the linear correlation between the target and the reconstructed image: a large value shows a high detectability in the reconstructed image and indicates a reconstructed image with high accuracy. These errors are defined as in [15], for the parameter ϑ being either κ or σ:(37)
e
1
ϑ
=
100
×
∑
j
=
1
N
ϑ
j
-
ϑ
t
j
2
∑
j
=
1
N
ϑ
t
j
2
,
e
2
ϑ
=
∑
j
=
1
N
ϑ
j
-
ϑ
¯
ϑ
t
j
-
ϑ
t
¯
N
-
1
σ
ϑ
σ
ϑ
t
,where N is the dimension of the parameter space (the number of finite elements nodes), ϑ
j is the value of ϑ on the ith node, and ϑ
¯ and ϑ
t
¯ are the mean values of vectors ϑ and ϑ
t, and σ
ϑ, σ
ϑ
t are standard deviations of target and reconstructed images, respectively. Standard deviations are given by(38)
σ
ϑ
=
∑
j
=
1
N
ϑ
j
-
ϑ
¯
2
N
-
1
,
σ
ϑ
t
=
∑
j
=
1
N
ϑ
t
j
-
ϑ
t
¯
2
N
-
1
.Numerical tests were performed for different levels of noise. First, a small noise of 20 dB is added to the synthetic data (this corresponds to 1% of noise). Then, a moderate noise of 15 dB is added to the data (approximately 3% of noise). Finally, the highest level of noise of 10 dB is considered (10% of noise).Figure2 presents the evolution of dimensionless cost function values along with iterations. Note that the main used stopping criterion has been based on the cost function stabilization. As expected, it is seen from this figure that the higher the level of noise, the higher the cost function value at stabilization. Moreover, the convergence to the solution is slightly faster with Sobolev gradients than with ordinary ones, even though, far at the end, the value of the cost function is the lowest for the ordinary gradient.Figure 2
2D RTE-based optical tomography. Evolution of the cost function with respect to iterations, forv
=
+
∞, 200, 100, and 50 and a 20 dB, 15 dB, and 10 dB noise. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient. The algorithm stops at iteration k when j
k
-
j
k
-
1
/
j
k
-
1
<
1
0
-
2.Note that the BFGS optimizer is not based on matrix inversion, as it is the case for Gauss-Newton method and its derivative. It has been shown in other works that such a choice implicitly regularizes the inverse problem. Another regularization tool comes from the parameter mesh that is chosen sufficiently coarse to avoid the presence of too many local minima [37]. These tools make the reconstructions possible without the need of adding a penalization term to the cost function.Note also that the choice of the stopping criterion plays a crucial role in the effectiveness of the algorithm: a too weak stopping criterion may lead to useless solution, while a too severe stopping criterion may lead to a divergent solution. In fact, the stopping criterion based on cost function stabilization is another dummy regularization that is to be used in practical applications. The determination of an appropriate value for this criterion is a challenging task in theory because it highly depends on the application. For the needs of this paper, the value equal to1
0
-
2 has been determined empirically.Figures3, 4, and 5 present the evolutions, with respect to iterations, of deviation and correlation factors, respectively, to both parameters κ and σ, for the three levels of noise. These results show that when compared to ordinary gradients, the use of space-dependent Sobolev gradients decreases deviation factors and increases correlation factors. The optimum regularization parameter, v, is about 100 in this particular case: when v is too small, there is too much gradient diffusion, while when v is too high, this yields the ordinary L
2
(
D
) gradient.Figure 3
2D RTE-based optical tomography. 20 dB noise. Evolution of errorse
1 ((a), (b)) and e
2 ((c), (d)) for κ ((a), (c)) and σ ((b), (d)), with respect to iterations, for v
=
+
∞, 200, 100, and 50. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient. The algorithm stops at iteration k when j
k
-
j
k
-
1
/
j
k
-
1
<
1
0
-
2.
(a)
(b)
(c)
(d)Figure 4
2D RTE-based optical tomography. 15 dB noise. Evolution of errorse
1 ((a), (b)) and e
2 ((c), (d)) for κ ((a), (c)) and σ ((b), (d)), with respect to iterations, for v
=
+
∞, 200, 100, and 50. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient. The algorithm stops at iteration k when j
k
-
j
k
-
1
/
j
k
-
1
<
1
0
-
2.
(a)
(b)
(c)
(d)Figure 5
2D RTE-based optical tomography. 10 dB noise. Evolution of errorse
1 ((a), (b)) and e
2 ((c), (d)) for κ ((a), (c)) and σ ((b), (d)), with respect to iterations, for v
=
+
∞, 200, 100, and 50. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient. The algorithm stops at iteration k when j
k
-
j
k
-
1
/
j
k
-
1
<
1
0
-
2.
(a)
(b)
(c)
(d)Figures6 and 7 corroborate these conclusions. The Sobolev gradients filter the high frequency fluctuations and avoid the propagation of the noise to the cost function gradient and then to the reconstructions. It is observed that the regularization parameter v
=
100 gives the best reconstructions as regards smoothness and precision of reconstructions.Figure 6
2D RTE-based optical tomography. 20 dB noise. Reconstructions obtained forκ (a) and σ (b). 1st row: targets; then v
=
+
∞, 200, 100, and 50, respectively. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient.
(a)
(b)Figure 7
2D RTE-based optical tomography. 10 dB noise. Reconstructions obtained forκ (a) and σ (b). 1st row: targets; then v
=
+
∞, 200, 100, and 50, respectively. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient.
(a)
(b)


## 3.1. Inverse Problem Statement

The inverse problem consists mathematically in minimizing a cost function which quantitatively measures the discrepancy between some optical measurements and related predictions over a set of sources, sensors, and angles of directional integration. In optical tomography, taking into account the reflection at the interface and assuming that the medium is convex, the measurable quantity is the exitance, or portion of it in a given solid angleω
(
s
)
⊆
s
·
n
>
0 centered on direction s. The related prediction can be written as (23)
E
x
,
s
=
∫
ω
s
1
-
R
x
,
s
I
x
,
s
s
·
n
d
ω
s
,where I is the radiative intensity and the function R
(
x
,
s
)
∈
[
0,1
] is the Fresnel reflection factor. In what follows in this section, the border of the medium is assumed to be transparent; that is, R
=
0. Moreover, the solid angle ω
(
s
) is assumed to be sufficiently small so that the approximation E
(
x
,
s
r
)
≈
ω
r
I
(
x
,
s
r
)
s
r
·
n holds. Therefore, the radiative intensity measurements I
˘
(
x
,
s
r
) can be assessed and the cost function (24) mapping from P
=
L
∞
(
D
)
×
L
∞
(
D
) to R uses, according to the schematic description of the problem (Figure 1), the specific Euclidean norm:(24)
j
κ
,
σ
≔
J
I
=
1
2
s
·
n
I
˘
I
-
I
˘
X
2
with
u
X
≜
∑
k
=
1
K
∑
d
=
1
D
∑
r
=
1
R
u
k
x
d
,
s
r
2
1
/
2
,where the index k refers to the source (test) number, the index d refers to the detection node, and the index r refers to the direction. In (24), J
(
I
) is the cost function explicitly expressed in terms of the prediction I, while j
(
κ
,
σ
) is the cost function implicitly expressed in terms of the parameters. This latter function is often called the reduced cost function. According to (1), one has the straightforward correspondence f
(
u
-
u
d
)
↔
s
·
n
/
I
˘
(
I
-
I
˘
). Note that the scalar s
·
n in (24) is introduced in order to avoid infinite Dirichlet boundary condition in the adjoint-state problem. For a given source configuration, say the kth one, the prediction I
k
(
x
d
,
s
r
) at the sensor d and direction r is given through the solution of the forward steady-state radiative transfer problem formulated as(25)
Find
I
k
x
,
s
such  that
∀
s
∈
S
n
-
1
,
∀
k
=
1
,
…
,
K
:
s
·
∇
+
κ
+
σ
I
k
x
,
s
=
σ
∮
S
n
-
1
I
k
x
,
s
Φ
s
,
s
′
d
ω
s
′
,
∀
x
∈
D
,
I
k
x
,
s
=
I
k
0
x
×
1
x
∈
∂
D
k
0
×
1
s
=
s
k
0
,
∀
x
∈
∂
D
,where s is the considered direction of propagation, Φ
(
s
,
s
′
) is the phase function representing the probability that a photon arriving from the direction s
′ is scattered to the direction s, κ and σ are the absorption and scattering space-dependent functions, respectively, and 1
[
·
] denotes the indicator function.Figure 1
Schematic description of the experiment. Only one source is on at once, while others are off. For each source configuration numberk
=
1
,
…
,
K, the emerging intensity is measured on all sensors d
=
1
,
…
,
D in different directions r
=
1
,
…
,
R. This is a theoretical extension of the experimental directional-directional spectral transmittances device using a Fourier transform infrared spectrometer described in [30].


## 3.2. Cost Function Gradient Derivation

The computation of the cost function gradient follows the generic methodology detailed in the previous section. The cost function derivative towards the directionϕ is given differentiating (24):(26)
j
′
κ
,
σ
;
ϕ
=
∇
j
κ
,
σ
,
ϕ
Z
=
s
·
n
I
˘
2
I
-
I
˘
,
I
′
X
∀
ϕ
∈
P
.In order to derive the cost gradient for this specific application,Y-based inner products are defined.Definition 11.
Consider(27)
u
,
w
Y
≜
∑
k
∮
S
n
-
1
∫
D
u
w
d
x
d
ω
s
,
u
,
w
∂
Y
≜
∑
k
∮
S
n
-
1
∫
∂
D
u
w
d
x
d
ω
s
.TheY-inner product satisfies Properties 1 and 2, both being useful for the cost function gradient derivation.Property 1.
Consider(28)
s
·
∇
u
,
w
Y
=
-
s
·
∇
w
,
u
Y
+
s
·
n
w
,
u
∂
Y
.Property 2.
Consider(29)
σ
∮
S
n
-
1
u
x
,
s
′
ϕ
s
,
s
′
d
ω
s
′
,
w
x
,
s
Y
=
∑
k
∮
S
n
-
1
∫
D
σ
∮
S
n
-
1
u
x
,
s
′
ϕ
s
,
s
′
d
ω
s
′
w
x
,
s
d
x
d
ω
s
=
∑
k
∫
D
σ
∮
S
n
-
1
∮
S
n
-
1
w
x
,
s
u
x
,
s
′
ϕ
s
,
s
′
d
ω
s
′
d
ω
s
d
x
=
∑
k
∫
D
∮
S
n
-
1
σ
∮
S
n
-
1
w
x
,
s
ϕ
s
,
s
′
d
ω
s
u
x
,
s
′
d
x
d
ω
s
′
=
σ
∮
S
n
-
1
w
x
,
s
ϕ
s
,
s
′
d
ω
s
′
,
u
x
,
s
′
Y
.The generic equation for the cost gradient (15) applied on both space-distributed parameters κ and σ gives, for direction ϕ
=
ϕ
κ
;
0
∈
P,(30)
∇
j
κ
,
σ
,
ϕ
Z
=
I
x
,
s
I
∗
x
,
s
,
ϕ
Yand, for direction ϕ
=
0
;
ϕ
σ
∈
P,(31)
∇
j
κ
,
σ
,
ϕ
Z
=
I
x
,
s
I
∗
x
,
s
-
∮
S
n
-
1
I
x
,
s
ϕ
s
,
s
′
d
ω
s
′
I
∗
x
,
s
,
ϕ
Y
.In these two equations, the directionϕ
∈
P is split into two terms ϕ
κ and ϕ
σ in order to be able to compute the cost function gradient relatively to each of the functions κ and σ separately.On one hand, (30) is rewritten such that(32)
∇
j
κ
,
σ
,
ϕ
Z
=
ϕ
I
x
,
s
,
I
∗
x
,
s
Yand, on the other hand, Property 2 enables rewriting (31) such that(33)
∇
j
κ
,
σ
,
ϕ
Z
=
ϕ
I
x
,
s
-
ϕ
∮
S
n
-
1
I
x
,
s
ϕ
s
,
s
′
d
ω
s
′
,
I
∗
x
,
s
Y
.Besides, the differentiation of the forward model enables rewriting previous equations too, forϕ being either ϕ
κ
;
0 or 0
;
ϕ
σ:(34)
∇
j
κ
,
σ
,
ϕ
Z
=
-
s
·
∇
+
κ
+
σ
I
′
x
,
s
,
I
∗
Y
+
σ
∮
S
n
-
1
I
′
x
,
s
ϕ
s
,
s
′
d
ω
s
′
,
I
∗
Y
.Properties1 and 2 enable rewriting (34) such that(35)
∇
j
κ
,
σ
,
ϕ
Z
=
s
·
∇
-
κ
-
σ
I
∗
x
,
s
,
I
′
Y
-
s
·
n
I
∗
,
I
′
∂
Y
+
σ
∮
S
n
-
1
I
∗
x
,
s
ϕ
s
,
s
′
d
ω
s
′
,
I
′
Y
.The adjoint problem is eventually obtained identifying (26) with (35):(36)
Find
I
k
x
,
s
such  that
∀
s
∈
S
n
-
1
,
∀
k
=
1
,
…
,
k
:
-
s
·
∇
+
κ
+
σ
I
k
∗
x
,
s
=
σ
∮
S
n
-
1
I
k
∗
x
,
s
ϕ
s
,
s
′
d
ω
s
′
,
∀
x
∈
D
,
-
s
·
n
I
k
∗
x
,
s
=
s
·
n
I
˘
k
2
I
k
-
I
˘
k
x
,
s
×
1
x
∈
∂
D
d
x
,
∀
x
∈
∂
D
,where ∂
D
d denotes the set of the pointwise sensors.


## 3.3. Numerical Results

The optimization relies on the BFGS algorithm described by Liu and Nocedal [51] coupled with a quadratic interpolation line search [52]. Generally speaking, this latter has proved to be particularly effective in these recent years in the area of inversion (see [53–55] to name but a few) due to its effectiveness in minimizing nonlinear problems [56].The two-dimensional domain is a square of 2 cm length. Four collimated sources are located at the center of each side, while sensors are placed around sources towards twenty directions. The target properties to be reconstructed,κ
t and σ
t, are defined by the background for which κ
=
0.25 cm−1 and σ
=
20 cm−1 and two square inclusions, the former for which κ
=
0.15 cm−1 and σ
=
10 cm−1 and the latter for which κ
=
0.35 cm−1 and σ
=
30 cm−1. Both the forward and the adjoint models are solved using P
1 discontinuous Galerkin finite elements for which the mesh convergence has been validated. Three distinct meshes have been created, one in order to build the synthetic data before adding random noise following Klose [57] (3 721 nodes), one in order to project the forward and adjoint-state variables (1 681 nodes), and one associated with the property functions with the aim of implicitly regularizing the inverse problem (441 nodes). Seeking to gauge the accuracy of the reconstructions, the deviation factor e
1
∈
[
0
,
+
∞
[ and the correlation factor e
2
∈
[
-
1
,
+
1
] are used [58]. The former factor measures the deviation of the reconstructed image from the target image: a small value indicates a reconstructed image with high accuracy. The latter factor measures the linear correlation between the target and the reconstructed image: a large value shows a high detectability in the reconstructed image and indicates a reconstructed image with high accuracy. These errors are defined as in [15], for the parameter ϑ being either κ or σ:(37)
e
1
ϑ
=
100
×
∑
j
=
1
N
ϑ
j
-
ϑ
t
j
2
∑
j
=
1
N
ϑ
t
j
2
,
e
2
ϑ
=
∑
j
=
1
N
ϑ
j
-
ϑ
¯
ϑ
t
j
-
ϑ
t
¯
N
-
1
σ
ϑ
σ
ϑ
t
,where N is the dimension of the parameter space (the number of finite elements nodes), ϑ
j is the value of ϑ on the ith node, and ϑ
¯ and ϑ
t
¯ are the mean values of vectors ϑ and ϑ
t, and σ
ϑ, σ
ϑ
t are standard deviations of target and reconstructed images, respectively. Standard deviations are given by(38)
σ
ϑ
=
∑
j
=
1
N
ϑ
j
-
ϑ
¯
2
N
-
1
,
σ
ϑ
t
=
∑
j
=
1
N
ϑ
t
j
-
ϑ
t
¯
2
N
-
1
.Numerical tests were performed for different levels of noise. First, a small noise of 20 dB is added to the synthetic data (this corresponds to 1% of noise). Then, a moderate noise of 15 dB is added to the data (approximately 3% of noise). Finally, the highest level of noise of 10 dB is considered (10% of noise).Figure2 presents the evolution of dimensionless cost function values along with iterations. Note that the main used stopping criterion has been based on the cost function stabilization. As expected, it is seen from this figure that the higher the level of noise, the higher the cost function value at stabilization. Moreover, the convergence to the solution is slightly faster with Sobolev gradients than with ordinary ones, even though, far at the end, the value of the cost function is the lowest for the ordinary gradient.Figure 2
2D RTE-based optical tomography. Evolution of the cost function with respect to iterations, forv
=
+
∞, 200, 100, and 50 and a 20 dB, 15 dB, and 10 dB noise. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient. The algorithm stops at iteration k when j
k
-
j
k
-
1
/
j
k
-
1
<
1
0
-
2.Note that the BFGS optimizer is not based on matrix inversion, as it is the case for Gauss-Newton method and its derivative. It has been shown in other works that such a choice implicitly regularizes the inverse problem. Another regularization tool comes from the parameter mesh that is chosen sufficiently coarse to avoid the presence of too many local minima [37]. These tools make the reconstructions possible without the need of adding a penalization term to the cost function.Note also that the choice of the stopping criterion plays a crucial role in the effectiveness of the algorithm: a too weak stopping criterion may lead to useless solution, while a too severe stopping criterion may lead to a divergent solution. In fact, the stopping criterion based on cost function stabilization is another dummy regularization that is to be used in practical applications. The determination of an appropriate value for this criterion is a challenging task in theory because it highly depends on the application. For the needs of this paper, the value equal to1
0
-
2 has been determined empirically.Figures3, 4, and 5 present the evolutions, with respect to iterations, of deviation and correlation factors, respectively, to both parameters κ and σ, for the three levels of noise. These results show that when compared to ordinary gradients, the use of space-dependent Sobolev gradients decreases deviation factors and increases correlation factors. The optimum regularization parameter, v, is about 100 in this particular case: when v is too small, there is too much gradient diffusion, while when v is too high, this yields the ordinary L
2
(
D
) gradient.Figure 3
2D RTE-based optical tomography. 20 dB noise. Evolution of errorse
1 ((a), (b)) and e
2 ((c), (d)) for κ ((a), (c)) and σ ((b), (d)), with respect to iterations, for v
=
+
∞, 200, 100, and 50. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient. The algorithm stops at iteration k when j
k
-
j
k
-
1
/
j
k
-
1
<
1
0
-
2.
(a)
(b)
(c)
(d)Figure 4
2D RTE-based optical tomography. 15 dB noise. Evolution of errorse
1 ((a), (b)) and e
2 ((c), (d)) for κ ((a), (c)) and σ ((b), (d)), with respect to iterations, for v
=
+
∞, 200, 100, and 50. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient. The algorithm stops at iteration k when j
k
-
j
k
-
1
/
j
k
-
1
<
1
0
-
2.
(a)
(b)
(c)
(d)Figure 5
2D RTE-based optical tomography. 10 dB noise. Evolution of errorse
1 ((a), (b)) and e
2 ((c), (d)) for κ ((a), (c)) and σ ((b), (d)), with respect to iterations, for v
=
+
∞, 200, 100, and 50. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient. The algorithm stops at iteration k when j
k
-
j
k
-
1
/
j
k
-
1
<
1
0
-
2.
(a)
(b)
(c)
(d)Figures6 and 7 corroborate these conclusions. The Sobolev gradients filter the high frequency fluctuations and avoid the propagation of the noise to the cost function gradient and then to the reconstructions. It is observed that the regularization parameter v
=
100 gives the best reconstructions as regards smoothness and precision of reconstructions.Figure 6
2D RTE-based optical tomography. 20 dB noise. Reconstructions obtained forκ (a) and σ (b). 1st row: targets; then v
=
+
∞, 200, 100, and 50, respectively. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient.
(a)
(b)Figure 7
2D RTE-based optical tomography. 10 dB noise. Reconstructions obtained forκ (a) and σ (b). 1st row: targets; then v
=
+
∞, 200, 100, and 50, respectively. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient.
(a)
(b)


## 4. Application on 3D DA-Based Optical Tomography

### 4.1. Inverse Problem Statement

The radiative transfer equation (RTE) provides an equation of light propagation valid in most of participating (absorbing and scattering) media as long as the independent scattering regime is fulfilled [59]. For specific applications, approximations of the RTE could be applied. Solutions of the RTE being computationally expensive, especially if the third dimension has to be considered, approximated models of the RTE should be employed when valid. Among the different models approaching the RTE, the diffuse approximation (DA) is the privileged one. The latter, which describes the photon density field φ
(
x
)
=
∮
S
n
-
1
‍
I
(
x
,
s
)
d
ω
(
s
) inside the medium, is given in the frequency domain by [1, 60](39)
-
∇
·
D
x
∇
φ
^
k
x
,
f
0
+
κ
x
+
2
π
i
f
0
c
φ
^
k
x
,
f
0
=
X
^
k
x
,
f
0
,
∀
x
∈
D
,
2
γ
φ
^
k
x
,
f
0
+
A
D
x
∇
φ
^
k
x
,
f
0
·
n
=
2
x
^
k
x
,
f
0
,
∀
x
∈
∂
D
.Here, X
^
k
(
x
,
f
0
) is the Fourier transform at the frequency f
0 of the internal source of radiation, X
k [W cm−3]; x
^
k
(
x
,
f
0
) is the Fourier transform at the frequency f
0 of the temporal signal associated with the kth diffuse source, x
k [W cm−2]; φ
^
k
:
D
×
C
↦
R is the Fourier transform of the kth photon density field, φ
k
:
D
×
R
+
↦
R
[
W
c
m
-
2
]; D
(
x
)
=
n
D
κ
(
x
)
+
σ
′
(
x
)
-
1 is the macroscopic scattering coefficient [cm]; n
D
(
=
2
,
3
) is the dimension of D; κ and σ
′
=
(
1
-
g
)
σ are the absorption and reduced scattering coefficients (σ is the scattering coefficient as previously defined) [cm−1]; g is the asymmetry factor which equals the average cosine of the scattering angle; c is the speed of light in D (constant) [cm s−1]; n is the unit normal vector to the boundary of D; γ is a constant parameter depending on n
D (1
/
π for n
D
=
2, 1
/
4 for n
D
=
3); and A is a parameter which characterizes the reflection on the border ∂
D that can be derived from the Fresnel laws if specular reflection is considered [61, 62].Similarly to Section3.3, the inverse problem consists in solving the optimization problem:(40)
Find
κ
,
σ
′
∈
P
=
L
∞
D
×
L
∞
D
such  that
κ
,
σ
′
=
arg
min
κ
,
σ
′
∈
P
⁡
j
κ
,
σ
′
,where(41)
j
κ
,
σ
′
≔
J
φ
^
=
1
2
φ
^
-
φ
˘
^
φ
˘
^
X
2
with
u
X
≜
∑
k
=
1
K
∑
d
=
1
D
u
k
x
d
2
1
/
2and where the index k refers to the source number and the index d refers to the detection node. It should be noted that, under the assumption of the DA, the measurable quantity (23) is related to the photon density through the equality E
(
x
)
=
E
(
x
,
s
·
n
>
0
)
=
2
γ
φ
(
x
) [60]. In the remainder of the section, the following values of the source terms are considered: X
^
k
=
0 W cm−2 and x
^
k
=
0.1 W cm−2 for all k
=
1
,
…
,
K (which corresponds to a temporal signal of the form cos
⁡
(
2
π
f
0
t
)).


### 4.2. Cost Function Gradients Derivation

The cost function gradients with respect toκ and σ
′ have been computed in [23] with the help of a continuous Lagrangian formulation. In terms of (5), it can be shown that the directional derivatives of the cost function are equal to(42)
∇
j
ψ
,
ϕ
Z
=
∑
k
=
1
K
R
φ
^
k
φ
^
¯
k
∗
-
n
D
D
2
∇
φ
^
k
·
∇
φ
^
¯
k
∗
,
ϕ
Z
∀
ϕ
=
ϕ
κ
;
0
∈
P
,
∇
j
ψ
,
ϕ
Z
=
∑
k
=
1
K
R
-
n
D
D
2
∇
φ
^
k
·
∇
φ
^
¯
k
∗
,
ϕ
Z
∀
ϕ
=
0
;
ϕ
σ
′
∈
P
,where R is the real part operator and the adjoint variable φ
^
k
∗ involved in (15) is the solution to the adjoint problem:(43)
-
∇
·
D
∇
φ
^
k
∗
+
κ
-
2
π
i
f
0
c
φ
^
k
∗
=
0
∀
x
∈
D
,
2
γ
A
φ
^
k
∗
+
D
∇
φ
^
k
∗
·
n
=
-
1
φ
˘
^
k
2
φ
^
k
-
φ
˘
^
k
1
x
∈
∂
D
d
x
∀
x
∈
∂
D
.Therefore, if the inner productZ
=
L
2
(
D
) is used to extract the cost function gradient in (5), the expressions for the continuous cost function gradients are such that(44)
∇
κ
j
=
∑
k
=
1
K
R
φ
^
k
φ
^
¯
k
∗
-
n
D
D
2
∇
φ
^
k
·
∇
φ
^
¯
k
∗
,
∇
σ
′
j
=
∑
k
=
1
K
R
-
n
D
D
2
∇
φ
^
k
·
∇
φ
^
¯
k
∗
.Note that an adimensionalisation of the parameters bya priori radiative properties is used in the inversion process [23].


### 4.3. Numerical Results

The same algorithm as employed in Section3 is used to reconstruct the target properties from initial estimates equal to values of those of the background. The domain considered, D, is a square of 5 cm length (D
=
[
0
;
5
]
2 cm2). The target properties to be reconstructed, κ
t and σ
t
′, are defined by(45)
κ
t
x
=
0.1
+
0.05
×
1
x
∈
D
1
x
+
0.1
×
1
x
∈
D
2
x
cm
-
1
,
σ
t
′
x
=
10
+
5
×
1
x
∈
D
1
x
+
10
×
1
x
∈
D
2
x
cm
-
1
,where D
1
=
[
0.75
;
1.75
]
2 cm2 and D
2
=
[
3.25
;
4.25
]
2 cm2.The stopping criteria are based on the cost function stabilization, with a critical value chosen equal to1
0
-
3 and a total number of iterations equal to 100. As in Section 3.3, these two criteria are used as dummy regularizations in order to get a low enough cost function value still without getting a divergent solution. Again, these parameters have been tuned empirically.Physical properties involved in the forward model (39) are fixed to f
0
=
100 MHz, c
=
(
3
×
1
0
10
/
1.4
) cm s−1, and γ
=
0.25. Synthetic data are considered for these numerical tests. These data, φ
˘
^, which represent pseudoexperimental measurements, are built on a finer mesh (132 651 nodes) than the one that generates the predictions φ
^ (68 921 nodes) in order to avoid the inverse crime [60, 63]. Then, a 10 dB multiplicative white Gaussian noise is applied to φ
˘
^ at the nodes of the sensors to simulate the noise inherent to experimental devices; that is, φ
˘
^
n
o
i
s
y
=
φ
˘
^
1
+
0.1
×
m, where m
~
N
(
0,1
). Synthetic data, state, and parameter meshes can be seen in Figure 8. It can be seen from Figure 8 that 6 sources and 24 sensors, 0.5 cm2 each (containing 16 discretization points), are used for this particular test. A reduction of the parameter space dimension is employed to improve the quality of the reconstructions obtained with the BFGS algorithm [23]: the mesh associated with the radiative properties κ and σ
′ contains only 9 261 nodes. Note that the total number of complex measurements is equal to 2 304.Figure 8
3D DA-based optical tomography. Synthetic data mesh: 132 651 nodes (a); mesh of the state and adjoint variable: 68 921 nodes (b); mesh associated with radiative propertiesκ and σ
′: 9 261 nodes (c).
(a)
(b)
(c)Figure9 presents the evolutions, with respect to iterations, of deviation and correlation factors, respectively, to both parameters. It is seen from this figure that when compared to ordinary L
2 gradients, the use of space-dependent Sobolev gradients decreases deviation factors and increases correlation factors. The optimum regularization parameter is about 10 in this particular case. Figure 10 corroborates these conclusions: the best reconstructions, in terms of smoothness and accuracy, are obtained for a value of the regularization parameter, v, equal to 10. Finally, remark that Figure 9 allows observing the establishment of a cross talk effect [64, 65]: for v
=
+
∞ and v
=
70, it is seen that the deviation factor e
1
κ decreases with the iterations, while the deviation factor e
1
σ increases. In these two cases, values of the cost function decrease with the iterations, while an improvement of the absorption coefficient reconstruction and a deterioration of the reduced scattering coefficient reconstruction are observed. This effect is prevented by the use of space-dependent Sobolev gradients with an appropriate regularization parameter weight.Figure 9
3D DA-based optical tomography. 10 dB noise. Evolution of errorse
1 ((a), (b)) and e
2 ((c), (d)) for κ ((a), (c)) and σ
′ ((b), (d)), with respect to iterations, for v
=
+
∞, 70, 40, and 10. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient. The algorithms stop at iteration k when j
k
-
j
k
-
1
/
j
k
-
1
<
1
0
-
3.
(a)
(b)
(c)
(d)Figure 10
3D DA-based optical tomography. 10 dB noise. Reconstructions obtained in the plane passing through the center of the cube of which the normal vectorn equals 1
,
-
1,0
⊤, for κ (a) and σ
′ (b). 1st row: targets; then v
=
+
∞, 70, 40, and 10, respectively. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient.
(a)
(b)


## 4.1. Inverse Problem Statement

The radiative transfer equation (RTE) provides an equation of light propagation valid in most of participating (absorbing and scattering) media as long as the independent scattering regime is fulfilled [59]. For specific applications, approximations of the RTE could be applied. Solutions of the RTE being computationally expensive, especially if the third dimension has to be considered, approximated models of the RTE should be employed when valid. Among the different models approaching the RTE, the diffuse approximation (DA) is the privileged one. The latter, which describes the photon density field φ
(
x
)
=
∮
S
n
-
1
‍
I
(
x
,
s
)
d
ω
(
s
) inside the medium, is given in the frequency domain by [1, 60](39)
-
∇
·
D
x
∇
φ
^
k
x
,
f
0
+
κ
x
+
2
π
i
f
0
c
φ
^
k
x
,
f
0
=
X
^
k
x
,
f
0
,
∀
x
∈
D
,
2
γ
φ
^
k
x
,
f
0
+
A
D
x
∇
φ
^
k
x
,
f
0
·
n
=
2
x
^
k
x
,
f
0
,
∀
x
∈
∂
D
.Here, X
^
k
(
x
,
f
0
) is the Fourier transform at the frequency f
0 of the internal source of radiation, X
k [W cm−3]; x
^
k
(
x
,
f
0
) is the Fourier transform at the frequency f
0 of the temporal signal associated with the kth diffuse source, x
k [W cm−2]; φ
^
k
:
D
×
C
↦
R is the Fourier transform of the kth photon density field, φ
k
:
D
×
R
+
↦
R
[
W
c
m
-
2
]; D
(
x
)
=
n
D
κ
(
x
)
+
σ
′
(
x
)
-
1 is the macroscopic scattering coefficient [cm]; n
D
(
=
2
,
3
) is the dimension of D; κ and σ
′
=
(
1
-
g
)
σ are the absorption and reduced scattering coefficients (σ is the scattering coefficient as previously defined) [cm−1]; g is the asymmetry factor which equals the average cosine of the scattering angle; c is the speed of light in D (constant) [cm s−1]; n is the unit normal vector to the boundary of D; γ is a constant parameter depending on n
D (1
/
π for n
D
=
2, 1
/
4 for n
D
=
3); and A is a parameter which characterizes the reflection on the border ∂
D that can be derived from the Fresnel laws if specular reflection is considered [61, 62].Similarly to Section3.3, the inverse problem consists in solving the optimization problem:(40)
Find
κ
,
σ
′
∈
P
=
L
∞
D
×
L
∞
D
such  that
κ
,
σ
′
=
arg
min
κ
,
σ
′
∈
P
⁡
j
κ
,
σ
′
,where(41)
j
κ
,
σ
′
≔
J
φ
^
=
1
2
φ
^
-
φ
˘
^
φ
˘
^
X
2
with
u
X
≜
∑
k
=
1
K
∑
d
=
1
D
u
k
x
d
2
1
/
2and where the index k refers to the source number and the index d refers to the detection node. It should be noted that, under the assumption of the DA, the measurable quantity (23) is related to the photon density through the equality E
(
x
)
=
E
(
x
,
s
·
n
>
0
)
=
2
γ
φ
(
x
) [60]. In the remainder of the section, the following values of the source terms are considered: X
^
k
=
0 W cm−2 and x
^
k
=
0.1 W cm−2 for all k
=
1
,
…
,
K (which corresponds to a temporal signal of the form cos
⁡
(
2
π
f
0
t
)).


## 4.2. Cost Function Gradients Derivation

The cost function gradients with respect toκ and σ
′ have been computed in [23] with the help of a continuous Lagrangian formulation. In terms of (5), it can be shown that the directional derivatives of the cost function are equal to(42)
∇
j
ψ
,
ϕ
Z
=
∑
k
=
1
K
R
φ
^
k
φ
^
¯
k
∗
-
n
D
D
2
∇
φ
^
k
·
∇
φ
^
¯
k
∗
,
ϕ
Z
∀
ϕ
=
ϕ
κ
;
0
∈
P
,
∇
j
ψ
,
ϕ
Z
=
∑
k
=
1
K
R
-
n
D
D
2
∇
φ
^
k
·
∇
φ
^
¯
k
∗
,
ϕ
Z
∀
ϕ
=
0
;
ϕ
σ
′
∈
P
,where R is the real part operator and the adjoint variable φ
^
k
∗ involved in (15) is the solution to the adjoint problem:(43)
-
∇
·
D
∇
φ
^
k
∗
+
κ
-
2
π
i
f
0
c
φ
^
k
∗
=
0
∀
x
∈
D
,
2
γ
A
φ
^
k
∗
+
D
∇
φ
^
k
∗
·
n
=
-
1
φ
˘
^
k
2
φ
^
k
-
φ
˘
^
k
1
x
∈
∂
D
d
x
∀
x
∈
∂
D
.Therefore, if the inner productZ
=
L
2
(
D
) is used to extract the cost function gradient in (5), the expressions for the continuous cost function gradients are such that(44)
∇
κ
j
=
∑
k
=
1
K
R
φ
^
k
φ
^
¯
k
∗
-
n
D
D
2
∇
φ
^
k
·
∇
φ
^
¯
k
∗
,
∇
σ
′
j
=
∑
k
=
1
K
R
-
n
D
D
2
∇
φ
^
k
·
∇
φ
^
¯
k
∗
.Note that an adimensionalisation of the parameters bya priori radiative properties is used in the inversion process [23].


## 4.3. Numerical Results

The same algorithm as employed in Section3 is used to reconstruct the target properties from initial estimates equal to values of those of the background. The domain considered, D, is a square of 5 cm length (D
=
[
0
;
5
]
2 cm2). The target properties to be reconstructed, κ
t and σ
t
′, are defined by(45)
κ
t
x
=
0.1
+
0.05
×
1
x
∈
D
1
x
+
0.1
×
1
x
∈
D
2
x
cm
-
1
,
σ
t
′
x
=
10
+
5
×
1
x
∈
D
1
x
+
10
×
1
x
∈
D
2
x
cm
-
1
,where D
1
=
[
0.75
;
1.75
]
2 cm2 and D
2
=
[
3.25
;
4.25
]
2 cm2.The stopping criteria are based on the cost function stabilization, with a critical value chosen equal to1
0
-
3 and a total number of iterations equal to 100. As in Section 3.3, these two criteria are used as dummy regularizations in order to get a low enough cost function value still without getting a divergent solution. Again, these parameters have been tuned empirically.Physical properties involved in the forward model (39) are fixed to f
0
=
100 MHz, c
=
(
3
×
1
0
10
/
1.4
) cm s−1, and γ
=
0.25. Synthetic data are considered for these numerical tests. These data, φ
˘
^, which represent pseudoexperimental measurements, are built on a finer mesh (132 651 nodes) than the one that generates the predictions φ
^ (68 921 nodes) in order to avoid the inverse crime [60, 63]. Then, a 10 dB multiplicative white Gaussian noise is applied to φ
˘
^ at the nodes of the sensors to simulate the noise inherent to experimental devices; that is, φ
˘
^
n
o
i
s
y
=
φ
˘
^
1
+
0.1
×
m, where m
~
N
(
0,1
). Synthetic data, state, and parameter meshes can be seen in Figure 8. It can be seen from Figure 8 that 6 sources and 24 sensors, 0.5 cm2 each (containing 16 discretization points), are used for this particular test. A reduction of the parameter space dimension is employed to improve the quality of the reconstructions obtained with the BFGS algorithm [23]: the mesh associated with the radiative properties κ and σ
′ contains only 9 261 nodes. Note that the total number of complex measurements is equal to 2 304.Figure 8
3D DA-based optical tomography. Synthetic data mesh: 132 651 nodes (a); mesh of the state and adjoint variable: 68 921 nodes (b); mesh associated with radiative propertiesκ and σ
′: 9 261 nodes (c).
(a)
(b)
(c)Figure9 presents the evolutions, with respect to iterations, of deviation and correlation factors, respectively, to both parameters. It is seen from this figure that when compared to ordinary L
2 gradients, the use of space-dependent Sobolev gradients decreases deviation factors and increases correlation factors. The optimum regularization parameter is about 10 in this particular case. Figure 10 corroborates these conclusions: the best reconstructions, in terms of smoothness and accuracy, are obtained for a value of the regularization parameter, v, equal to 10. Finally, remark that Figure 9 allows observing the establishment of a cross talk effect [64, 65]: for v
=
+
∞ and v
=
70, it is seen that the deviation factor e
1
κ decreases with the iterations, while the deviation factor e
1
σ increases. In these two cases, values of the cost function decrease with the iterations, while an improvement of the absorption coefficient reconstruction and a deterioration of the reduced scattering coefficient reconstruction are observed. This effect is prevented by the use of space-dependent Sobolev gradients with an appropriate regularization parameter weight.Figure 9
3D DA-based optical tomography. 10 dB noise. Evolution of errorse
1 ((a), (b)) and e
2 ((c), (d)) for κ ((a), (c)) and σ
′ ((b), (d)), with respect to iterations, for v
=
+
∞, 70, 40, and 10. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient. The algorithms stop at iteration k when j
k
-
j
k
-
1
/
j
k
-
1
<
1
0
-
3.
(a)
(b)
(c)
(d)Figure 10
3D DA-based optical tomography. 10 dB noise. Reconstructions obtained in the plane passing through the center of the cube of which the normal vectorn equals 1
,
-
1,0
⊤, for κ (a) and σ
′ (b). 1st row: targets; then v
=
+
∞, 70, 40, and 10, respectively. Note that the case v
=
+
∞ corresponds to the ordinary L
2 gradient.
(a)
(b)


## 5. Conclusion

In this paper, inverse models based on the BFGS algorithm have been developed to solve optical tomography problems based on two different forward models: the bidimensional steady-state radiative transfer equation and the three-dimensional frequency diffuse approximation. A Sobolev filter function was defined and space-dependent Sobolev inner products were used when extracting the cost function gradients instead of theL
2-inner product. Two test cases consisting in reconstructing discontinuous radiative properties in a square and a cube have been considered. Numerical results obtained have shown that, for all levels of noise and all values of the regularization parameter tested, the use of Sobolev gradients drastically enhances the quality of the reconstructions. Also, it has been observed that the space-dependent Sobolev gradients can allow avoiding the establishment of a crosstalk effect in the case of the frequency three-dimensional approximation diffuse based inversion.Generally speaking, due to simplicity of implementation and effectiveness in filtering noise measurements locally that provides better reconstructions, space-dependent Sobolev gradients appear to be particularly attractive in the context of inverse problems. In the near future, the anisotropy factor will be integrated in the inverse problem based on the radiative transfer equation and the use of a logarithmic cost function will be considered as our more recent results suggest. Later, inversion numerical algorithms of the three-dimensional radiative transfer equation and diffuse approximation model will be integrated in a global numerical code. The latter will be based on the BFGS algorithm with the use of adimensionalisation of the radiative properties, reduction of the parameter space dimension, space-dependent Sobolev gradients, and a wavelet multiscale approach, as developed in another work.



---

*Source: 3543571-2016-03-29.xml*