# Introducing shape constraints into object-based traveltime tomography

Gil Gaullier, Pierre Charbonnier, Fabrice Heitz, Philippe Côte

► To cite this version:

Gil Gaullier, Pierre Charbonnier, Fabrice Heitz, Philippe Côte. Introducing shape constraints into object-based traveltime tomography. Inverse Problems, 2016, 32 (9), pp.095002. 10.1088/0266-5611/32/9/095002. hal-02177235

HAL Id: hal-02177235

https://hal.science/hal-02177235v1

Submitted on 5 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Introducing shape constraints into object-based
traveltime tomography

G. Gaullier¹, P. Charbonnier¹, F. Heitz² and P. Côte³

¹ CEREMA/DTerEst/Laboratoire de Strasbourg, Strasbourg, France

² Lab. ICube, UMR 7357 CNRS, Université de Strasbourg, France

³ Dept. GERS, IFSTTAR, Bouguenais, France

E-mail: pierre.charbonnier@cerema.fr

**Abstract.** Traveltime tomography is a difficult, ill-posed reconstruction problem due to the nonlinearity of the forward model and the limited number of measurements usually available. In such an adverse situation, pixel-based regularization methods are generally unable to provide satisfactory reconstructions. In this paper we propose a novel object-based reconstruction method that introduces prior information about the shape of the structures to be reconstructed, which yields high quality geoacoustic inversion. The proposed method approaches the forward model by a series of linear problems, leading to a sequence of minimizations during which the shape prior is introduced. The method is demonstrated on synthetic and real data, collected on a specific bench dedicated to non-destructive testing of civil engineering structures.

*Keywords:* Traveltime tomography, Level sets method, Geoacoustic inversion, Shape-constrained reconstruction

Submitted to: *Inverse Problems*
# 1. Introduction

In this paper, we consider traveltime tomography, which aims at reconstructing wave speed distributions from first-arrival time measurements (for a survey, see [1]). Such a problem is highly nonlinear because the ray paths that define the forward model depend on the unknown speed distribution. Standard pixel-based techniques linearize the energy criterion using a first-order approximation and apply algorithms such as Simultaneous Iterative Reconstruction Technique (SIRT to perform the optimization [2]. They yield poor reconstructions at objects boundaries, especially when a limited number of measurements is available, which is a common situation in practice.

To overcome the limitations of standard regularization, more sophisticated priors may be included into penalty terms based e.g. on sparsity [3] or bounded variation [4] criteria. Half-quadratic edge-preserving methods (see, e.g., [5]) introduce an auxiliary variable, called line process, to break the smoothing constraint at the location of discontinuities. This yields better reconstructions at object boundaries, but local discontinuity models are generally not sufficient at high noise levels or for a small number of projections. In contrast to such pixel-based methods, shape-based or object-based reconstruction explicitly models objects boundaries using evolving closed curves, namely active contours [6, 7, 8].

<figure><img src="image_2.png" /><figcaption>**Figure 1.** The object is modeled by a shape in $\mathbb{R}^2$.</figcaption></figure>

Shape-based approaches are in the spirit of the classic model of Mumford and Shah [9]. In this framework, the unknown field (here, we consider the reciprocal of the velocity, namely the slowness, $s$) is modeled as a piecewise continuous function related to a region $\Omega$ in $\mathbb{R}^2$ or $\mathbb{R}^3$ (Fig. 1). The inverse problem can be set up as follows: given a collection of noisy measurements, estimate a slowness map, $s(\cdot, \Omega)$, defined from a shape $\Omega$, as

$$s(x, \Omega) = s_{\Omega}(x)1_{\Omega}(x) + s_{\Omega^c}(x)1_{\Omega^c}(x), \quad (1)$$

where $1_{\Omega}$ denotes the indicator or characteristic function of $\Omega$. In (1), the functions $s_{\Omega}(x)$ and $s_{\Omega^c}(x)$ denote the slowness of (respectively) the object and the medium in which it is immersed. They depend on the shape $\Omega$, which is unknown and must therefore be optimized when solving the inverse problem.

Let us notice, as illustrated Fig. 1, that this model does not constrain the object to be connected and is thus able to deal with arbitrary topologies. However, its non-differentiable nature prevents from applying conventional linearization techniques, so
shape optimization schemes must be sought. Contrary to pixel-based methods, which aim at reconstructing a function defined on the whole investigated region $D$, shape-based approaches explicitly model the geometry and the functional form of the object of interest. This enables the introduction of global prior information about the slowness distribution of foreground and background regions, as well as on the objects shape. This has already proven to be beneficial in tomography for reconstruction problems associated to linear forward problems [7, 8, 10].

The object-based approach to inverse problems has first been used in [6] by Santosa in the framework of the level-set method to solve the problems of deconvolution and diffraction screen reconstruction. Since then, many efforts have been made to extend the ideas developed in [6, 11, 12] for linear reconstruction applications to a larger class of problems governed by partial differential equations (PDE), see for example [8, 13, 14]. We address here the inverse problem associated to a nonlinear direct problem governed by curve integrals, namely traveltime tomography. Object-based approaches for traveltime tomography have been introduced independently in two recent papers. In [15], a reconstruction method, derived from the adjoint-state theory and based on the level-set algorithm [16], is proposed. It assumes that the slowness of the background region $\Omega^c$ is known. While the method is based on piecewise constant or piecewise continuous slowness models, it makes no assumption about the geometry of the solution. In contrast, a stochastic method that optimizes a population of ellipses to reconstruct high velocity structures is proposed in [17]. In order to reduce the complexity of the direct problem, which is intensively used, due to the stochastic nature of the approach, all objects are supposed to have the same velocity and a (quite restrictive) high-contrast assumption is made.

In this paper, we introduce an object-based reconstruction method with shape constraints for nonlinear traveltime tomography. As in the above-mentioned approaches, a piecewise homogeneous model (1) is considered. Contrary to [17], our approach is deterministic and does not limit the shape, topology or physical properties of the object to reconstruct. To overcome the difficulties related to the nonlinear and non-differentiable nature of the problem, we consider a decoupling approach that leads to a series of linear functional minimizations, implemented using a level-set algorithm. Although more sophisticated solvers might be considered, we compute the direct problem with a parametric ray tracing (bending) algorithm, extracted from an existing SIRT software [2]. Using this versatile and efficient code leads to a rather fast implementation and allows a fair and direct comparison with a proven, well-known pixel-based reconstruction technique. Another contribution of this paper, that clearly differentiates our work from [15], is the introduction of high-level shape constraints in traveltime tomographic reconstruction. The proposed method is based on a hierarchical representation of shapes, so it may be adapted for invariance to geometric transformations. We assess the reconstruction performance of the method using simulated and real data, acquired in controlled conditions with an ultrasound test bench. Reconstruction results demonstrate the efficiency of our approach, even in the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Sources ($\alpha$) and sensors ($\beta$) are located on the boundary of the investigated domain, $D$.</figcaption></figure>

case of sparse data.

The rest of the paper is organized as follows. In Section 2, we introduce the traveltime tomography reconstruction problem. The proposed constrained shape evolution method is introduced in Section 3. Reconstruction results obtained on synthetic and real data are presented in Section 4 and compared to pixel-based SIRT reconstructions for different source/receiver configurations.

# 2. Traveltime tomographic reconstruction

## 2.1. *Principle and experimental setup*

Traveltime tomographic reconstruction is based on noisy measurements given by the first-arrival time of acoustic waves, propagated from a number of sources to a set of sensors, located on the boundary of the investigated domain, $D$ (Figure 2). Traveltime tomography aims at reconstructing at each spatial location, $x$, the wave propagation speed or its inverse, called *slowness* and denoted here by $s$. This cartography may be related to the physical properties of the sample in order, for example, to detect objects of interest, such as cavities or bedrocks.

In this paper, a source location is denoted by $\alpha$, whereas a sensor, or receiver location is denoted by $\beta$. The traveltime of the wavefront from $\alpha$ to $\beta$, associated to the scalar slowness $s$ will be denoted $T(\alpha, \beta, s)$. The set of measured traveltimes is stacked in a vector of noisy observation, $T_{obs}$.

In order to obtain real data in controlled conditions for comparing different reconstruction methods, we have used a specific test bench [18], dedicated to small-scale geophysical measurements (Figure 3). We consider a 3D object composed of an aluminium half-cylinder embedded in a $300 \text{ mm} \times 600 \text{ mm} \times 350 \text{ mm}$ block of resin. Traveltime measurements are made on the central slice of the object, in the following manner. A piezoelectric, ultrasonic source is successively placed at 35 locations distributed on three sides of the structure. The wavefronts traverse the aluminium at about $6,000 \text{ m.s}^{-1}$ and the resin, at approximately $3,000 \text{ m.s}^{-1}$. For each source location, an interferometric laser sensor is placed at 26 successive positions on the fourth side of the block. The sensor records the displacements of the surface induced by the
<figure><img src="image_4.png" /><figcaption>Figure 3. Experimental setup (a), the real 3D object on the test bench (b).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Example of the signal emitted by the source (left) and recorded by the sensor (right).</figcaption></figure>

wave, with angstrom accuracy. The first-arrival times of the wavefront are then pointed manually by an operator using dedicated software. An example of record of both the emitted and received signals is shown on Fig. 4, at an arbitrary scale.

## 2.2. Forward problem

The forward problem consists in computing the traveltimes $T(\alpha, \beta, s)$ between known end points $\alpha$ and $\beta$ given a slowness field, $s$. Considering the high-frequency approximation, we may use ray theory to describe the propagation of wave fronts [19]. Fermat's principle states that the ray path followed by the propagating wave minimizes the traveltime. In this paper, a shape-based approach, in which a discontinuous slowness field $s(\cdot, \Omega)$ defined by (1) is considered. Hence, the forward operator $T$, that depends on $\alpha, \beta$ and $s(\cdot, \Omega)$, reads:

$$T(\alpha, \beta, s(\cdot, \Omega)) = \inf_{\psi \in \Gamma_{\alpha\beta}} \int_{\psi} s(x, \Omega) dl \quad (2)$$

where $\Gamma_{\alpha\beta}$ denotes the set of all paths joining $\alpha$ to $\beta$. Equivalently, (2) may be written as:

$$T(\alpha, \beta, s(\cdot, \Omega)) = \int_{\gamma_{\alpha\beta}} s(x, \Omega) dl \quad (3)$$
where the integration path $\gamma_{\alpha\beta}$ satisfies:

$$ \gamma_{\alpha\beta} = \arg \inf_{\psi \in \Gamma_{\alpha\beta}} \int_{\psi} s(x, \Omega) dl. \qquad (4) $$

As a consequence, contrary to computerized tomography [20] where the integration paths are fixed, straight lines, here they depend on the slowness field $s(\cdot, \Omega)$, which makes the forward problem nonlinear with respect to $s$.

## 2.3. Inverse problem

Traveltime inversion aims at reconstructing the propagation velocity (or the slowness $s$) at each spatial location, using a collection of noisy first-arrival time measurements, $T_{obs}$. In the shape-based approach, which consists in determining the discontinuous slowness field defined in (1), the unknowns of the problem are the functions $s_{\Omega}(x), s_{\Omega^c}(x)$ and the boundaries of the shape $\Omega$.

The reconstruction problem is cast as the minimization of a nonlinear objective function. In the shape-based framework, we use the Least-Squares (LS) criterion

$$ \mathcal{J}(s_{\Omega}, s_{\Omega^c}, \Omega) = \frac{1}{2} \|T(\cdot, \cdot, s(\cdot, \Omega)) - T_{obs}\|^2. \qquad (5) $$

where $T(\cdot, \cdot, s(\cdot, \Omega))$ denotes the vector that collects the traveltimes (predicted by the model) between all sources $\alpha$ and sensors $\beta$.

A particularity of our application is that the direct problem is governed by curve integrals (3)-(4), in which the unknown slowness field, $s$, is involved in a nonlinear manner, through the optimal integration path. Hence, contrary to the case of linear tomography, it is not possible to derive a closed-form expression for the shape derivative of the criterion, which makes the optimization problem much harder. In the next section, we introduce a decoupling technique which alleviates this difficulty.

# 3. A shape-constrained optimization method

In this section, we first introduce a method to solve the nonlinear, unconstrained optimization problem (5), first with respect to the shape, $\Omega$, then w.r.t. the parameters of the slowness model (1). Then, we show how to introduce high level constraints on the geometry and topology of the sought object in the criterion. Finally, we describe the proposed algorithm.

## 3.1. Optimization with respect to shape

To circumvent the nonlinearity problem in the context of the piecewise homogeneous slowness distribution model (1), we resort to a decoupling technique, which leads to a series of linear optimization problems, that can be solved in the shape evolution framework.

The nonlinearity of the forward operator stems from the fact that the unknown shape, $\Omega$, appears both in the integrand of (2) and in the (minimal) path of integration.
We propose to separate the forward operator in two sub-steps: estimation of the path and traveltime computation. To this end, we define a sequence of shape estimates, $(\Omega_k)_{k \ge 0}$ which aims at approximating the optimal shape $\Omega^*$ for $k$ becoming large. Let us define the function $T^{\gamma^k}$ depending on $\alpha, \beta$, on the current shape estimate $\Omega_k$ and on the shape $\Omega$ to be estimated at step $k + 1$ as follows

$$T^{\gamma^k}(\alpha, \beta, \Omega) = \int_{\gamma_{\alpha\beta}^k} s(\cdot, \Omega)dl \qquad (6)$$

where the path of integration $\gamma_{\alpha\beta}^k$ satisfies Fermat's principle related to the slowness field (1) defined by the current shape $\Omega_k$, i.e.

$$\gamma_{\alpha\beta}^k = \arg\inf_{\psi \in \Gamma_{\alpha\beta}} \int_{\psi} s(\cdot, \Omega_k)dl. \qquad (7)$$

The shape $\Omega_{k+1}$, i.e. the shape at the next iteration of the procedure, is defined as the solution of the following shape optimization problem based on $T^{\gamma^k}$:

$$\Omega_{k+1} = \arg\inf_{\Omega} \frac{1}{2} \|T^{\gamma^k}(\cdot, \cdot, s(\cdot, \Omega)) - T_{\text{obs}}\|^2 \qquad (8)$$

where, as before, $T^{\gamma^k}(\cdot, \cdot, s(\cdot, \Omega))$ and $T_{\text{obs}}$ are vectors collecting respectively the predicted and observed travel times for all pairs of source-receiver locations. We will denote by $n_s$ the number of source locations and $n_r$ the number of receiver locations, which leads to vectors of size $n_s \times n_r$.

The main interest of (8) is that it involves a linear operator. Indeed, once the integration path $\gamma_{\alpha\beta}^k$ has been estimated using the current estimate $\Omega_k$ in (7), the curve integral in (6) is performed along a fixed curve. Since the latter is independent from the variations of the slowness field defined by the evolving shape $\Omega$, the operator defined by (6) is linear with respect to $s(\cdot, \Omega)$. Therefore, for a collection of $n_s$ source locations and $n_r$ receiver locations, the vector $T^{\gamma^k}(\cdot, \cdot, \Omega)$ of traveltimes can be written as a matrix-vector product:

$$T^{\gamma^k}(\cdot, \cdot, s(\cdot, \Omega)) = \mathcal{R}^k s(\cdot, \Omega). \qquad (9)$$

In this expression, $s(\cdot, \Omega)$ is the vector of slowness values at each pixel location and $\mathcal{R}^k$ is the so-called matrix of ray lengths [1, p 113], whose elements correspond to the distance covered by each ray in each pixel of the investigated domain $\mathcal{D}$. Finally, equation (8) can be written in its simplest form:

$$\Omega_{k+1} = \arg\inf_{\Omega} \frac{1}{2} \|\mathcal{R}^k s(\cdot, \Omega) - T_{\text{obs}}\|^2. \qquad (10)$$

The interest of this formulation, which is akin to usual tomography, is that it is possible to differentiate the criterion in closed form (see Appendix B) and hence, to compute $\Omega_{k+1}$ by shape evolution.

The decoupling technique we propose is a heuristic that is reminiscent of alternate optimization techniques, since we alternately fix the shape to compute the optimal raypaths, and fix the raypaths to optimize the shape. The difference is that alternate optimization considers a single criterion while here, there are two different
nested criteria (Fermat's principle and Least-Squares). A consequence is that the application of Fermat's principle may cause temporary increases of the LS energy. Our experimentation (see Sec. 4.1) shows that we indeed observe such temporary jumps in energy when recomputing the optimal raypaths. We also have observed in all our experiments that the algorithm yields an overall decrease of the LS energy and thus converges (experimentally). However, contrary to alternate optimization, establishing the convergence properties of the algorithm is a difficult mathematical question, that remains open for the moment.

## 3.2. Optimal slowness distributions

The least-squares criterion $||R^k s(\cdot, \Omega) - T_{obs}||^2$ must, simultaneously, be optimized w.r.t. the parameters of the slowness distributions inside and outside the shape. We consider the case where $s_\Omega$ and $s_{\Omega^c}$ are linear combinations of spatial basis functions, as proposed in [7] in the case of linear tomography. In other words, we write (1) as

$$s(x, \Omega) = \sum_{i=1}^{K} a_i s_i^\Omega(x) 1_\Omega(x) + \sum_{i=1}^{L} b_i s_i^{\Omega^c}(x) 1_\Omega^c(x) \quad (11)$$

where the basis functions, $s_i^\Omega(x)$ and $s_i^{\Omega^c}(x)$, may be polynomials or trigonometric functions, for example. Let us denote by $A_i^k$ (resp. $B_i^k$) the result of the forward operator with fixed rays (6) applied to function $s_i^\Omega 1_\Omega$ (resp. to function $s_i^{\Omega^c} 1_\Omega^c$):

$$A_i(\alpha, \beta, \Omega) = \int_{\gamma_{\alpha\beta}^k} s_i^\Omega 1_\Omega, \quad B_i(\alpha, \beta, \Omega) = \int_{\gamma_{\alpha\beta}^k} s_i^{\Omega^c} 1_\Omega^c. \quad (12)$$

We stack the $A_i^k$'s (resp. the $B_i^k$'s) into a vector $A_i^k$ (resp. $B_i^k$) of length $(n_s \times n_r)$ and collect the $a_i$'s and $b_i$'s into a parameter vector of length $(K + L)$:

$$\theta = (a_1, \dots, a_K, b_1, \dots, b_L)^T. \quad (13)$$

Under these notations and exploiting linearity, we express the vector of arrival times predicted by the forward model, $T^{\gamma^k}$ as a function of $\theta$

$$T^{\gamma^k} = R^k s = M^k \theta \quad (14)$$

where $M^k$ is a matrix whose K first columns are the $A_i^k$'s and the last L columns are the $B_i^k$'s. Hence, the optimal values of the slowness model parameters are given by

$$\hat{\theta} = \arg\inf_{\theta} \frac{1}{2} \|M^k \theta - T_{obs}\|^2. \quad (15)$$

Since $M^{k^T} M^k$ is semipositive definite and rarely singular [7], $\hat{\theta}$, is solution of the linear system

$$M^{k^T} M^k \hat{\theta} = M^{k^T} T_{obs} \quad (16)$$

which may be obtained through inversion (or pseudo-inversion). A standard case occurs when the background and object of interest may both be characterized by a constant slowness. Then, the parameter vector reduces to the pair $(s_\Omega, s_{\Omega^c})$ and $\hat{\theta}$ is solution of a 2 × 2 system.
### 3.3. Incorporating shape prior

The shape prior used to constrain the evolution of $\Omega$ was first introduced by the authors in [21], in the context of segmentation and used in [10] for linear tomography. It is based on the Legendre moments of the indicator function of a shape, truncated at the order $N$: $\lambda = \{\lambda_{p,q}, p + q \le N\}$. We recall that the geometric or regular moments of the indicator function of shape $\Omega$ are given for $u \ge 0, v \ge 0$ by:

$$M_{u,v}(\Omega) = \int_{\mathcal{D}} \mathbb{1}_{\Omega}(x) x_1^u x_2^v dx = \int_{\Omega} x_1^u x_2^v dx \quad (17)$$

where $x_1^u$ and $x_2^v$ are powers of the first and second coordinate of a point $x = (x_1, x_2)$, $x \in \mathcal{D} \subset \mathbb{R}^2$. From (17), Legendre moments can be considered by substituting the usual monomials by the Legendre polynomials, which have the property to be orthogonal in $[-1, 1]$:

$$P_p(x) = \frac{1}{2^p p!} \times \frac{d^p}{dx^p} (x^2 - 1)^p, \quad (18)$$

where $p$ is the degree of the polynomial. Moreover, Legendre moments can be made intrinsically invariant to geometric transformations by exploiting the linear relationship between Legendre and regular moments. We consider here scale and translation invariant moments‡ as follows [22]:

$$\lambda_{p,q}(\Omega) = C_{pq} \sum_{0 \le u \le p} \sum_{0 \le v \le q} a_{pu} a_{qv} \eta_{u,v}(\Omega) \quad (19)$$

where $a_{pu}, a_{qv}$ are the coefficient of $x_1^u$ and $x_2^v$ in the series expansion of the Legendre polynomials of order $p$ and $q$, respectively, i.e.

$$P_p(x_1) = \sum_{0 \le u \le p} a_{pu} x_1^u \quad \text{and} \quad P_q(x_2) = \sum_{0 \le v \le q} a_{qv} x_2^v, \quad (20)$$

$C_{pq} = (2p + 1)(2q + 1)/4$ and $\eta_{u,v}$'s are normalized central moments, defined as

$$\eta_{u,v}(\Omega) = \int_{\Omega} H_{u,v}(x, \Omega) dx \quad (21)$$

with

$$H_{u,v}(x, \Omega) = \frac{(x_1 - \bar{x}_1)^u (x_2 - \bar{x}_2)^v}{(S_F |\Omega|)^{\frac{u+v+2}{2}}}. \quad (22)$$

In (22), $(\bar{x}_1, \bar{x}_2)$ is the centroid of $\Omega$, $|\Omega|$ denotes its area and $S_F$ is a scale factor which is fixed so that the modeled shape fits in the $[-1, 1]^2$ interval (see [21] for more explanations). Note that these quantities are readily available once the first three regular moments have been computed, since $|\Omega| = M_{0,0}$, $\bar{x}_1 = M_{1,0}/M_{0,0}$ and $\bar{x}_2 = M_{0,1}/M_{0,0}$.

Shapes may be reconstructed from their invariant Legendre moments, as illustrated on Figure 5. In this example, a half-disk object is reconstructed from its moments up to different orders. As can be seen, the higher the order $N$, the sharper the reconstruction.

<footnote>‡ Invariance to rotations may also be incorporated at the expense of a harder optimization problem (see [21])</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 5. Reconstruction of a reference shape (half disk) from its invariant Legendre moments on [-1, 1]², for different maximal orders of the Legendre moments, N (bottom), with $S_F = 1$. The corresponding orders are marked by blue circles on the reconstruction error curve (top). Numerical instabilities arise for N = 50.</figcaption></figure>

All the reconstructions are centered, and their normalized scale is given by $1/S_F$. The upper part of the figure shows the corresponding error curve. Such curves are used in practice to choose the value of the maximal order of representation N in order to avoid numerical instabilities [22].

The shape representation by moments is hierarchical, compact and does not depend on shape topology, since it is based on the characteristic function. Since it is invariant with respect to scaling and translation, comparing the evolving shape with a set (or alphabet) of a priori known reference shapes $\Omega_j^{ref}$ (Figure 6) is made possible.

Denoting by $\lambda_j^{ref}$ (for $j = 1, \dots, N_{ref}$) the invariant Legendre descriptors of these reference shapes, the shape prior is defined as a mixture of equiprobable isotropic Gaussians, with variance $\sigma$, centered on each reference shape descriptor. This leads to the following prior energy:

$$J_{prior}(\Omega) = -\log \left( \sum_{j=1}^{N_{ref}} e^{-\frac{\|\lambda(\Omega) - \lambda_{(j)}^{ref}\|^2}{2\sigma^2}} \right). \quad (23)$$

In our implementation, $\sigma$ is chosen so that the classification error probability, between the two closest reference shapes (in terms of descriptors) is less than 3% [21]. Note that, in the spirit of [23], this model may be seen as a Parzen-window estimator, with Gaussian kernels of bandwidth $\sigma$, that approximates the probability density function of the reference shape descriptors when $N_{ref}$ becomes large.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Sample set (dictionary) of reference shapes used for reconstruction.</figcaption></figure>

In the proposed optimization process, the data term is defined at each step (see (8) and (10)) as:

$$J_k^{\text{data}}(\Omega) = \frac{1}{2} \|T^{\gamma^k}(\cdot, \cdot, s(\cdot, \Omega)) - T_{\text{obs}}\|^2 = \frac{1}{2} \|\mathcal{R}_k s(\cdot, \Omega) - T_{\text{obs}}\|^2. \quad (24)$$

We incorporate the shape prior by a combination with the data fidelity term at each step k:

$$J_k(\Omega) = (1 - \alpha_p) J_k^{\text{data}}(\Omega) + \alpha_p J_{\text{prior}}(\Omega). \quad (25)$$

where $\alpha_p$ is a trade-off parameter which tunes the balance between the data term and the prior.

### 3.4. A shape constraint reconstruction algorithm

Starting from an initial slowness field, $s(\cdot, \Omega_0)$, the proposed algorithm iterates the following two steps until convergence:

* Step 1. Compute the set of ray paths satisfying Fermat's principle related to the slowness field $s(\cdot, \Omega_k)$. This defines $\mathcal{R}^k$ (eq. 9).

* Step 2. Minimize (25) by gradient descent to obtain $\Omega_{k+1}$, while periodically updating the parameters $\theta$ (defined in eq. 13) of the slowness model, using (16).

In the following paragraphs, we present each step in more details, in the particular case of constant slowness values: $s_\Omega$ for the object and $s_\Omega^c$ for the background. This corresponds to the experimental setup that we have used.

#### 3.4.1. Initialization
We initialize $\Omega_0$, for instance with a simple disk, and set the initial values of $s_\Omega$ and $s_{\Omega^c}$ using the back-projection of the data. More specifically, we superimpose the mask given by the disk on the back-projection image and compute the mean values in terms of slowness inside and outside the mask. As the example of Figure 7 shows, we obtain a rather crude approximation of the shape, and the corresponding slowness fields, $s(\cdot, \Omega_0)$ is in general poorly contrasted. However, we will see in Sec. 4 that this initialization yields a satisfactory solution at convergence. Other initialization strategies (e.g. seed snakes) may be considered (see for instance [21]).

#### 3.4.2. Step 1: minimal path computation
At Step 1, different methods may be used to implement the computation of the traveltimes [1]. Ray tracing methods proceed either by shooting rays from the source with varying orientations until the ray trajectory
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Initialization from real data. Initial mask (a), back-projection of the data (units: m.s⁻¹) (b), initial slowness field, s(·, Ω₀) (c).</figcaption></figure>

accurately targets the receiver position or by bending a curve between the source and the receiver until it satisfies Fermat's principle. Graph-based methods implement Fermat's principle by finding minimal paths on a graph. Wavefront tracking methods simulate the propagation of waves using grid-based eikonal equation solvers such as *Fast Marching* [24] or *Fast Sweeping* [25].

In this paper, we consider a fast variant of the bending technique, which models rays as circular arcs [2]. More specifically, s(·, Ωₖ) being fixed, a parametric, hence fast) exhaustive search is performed over the values of the bending parameter of the circular arc that joins α and β, to select the one that leads to the minimal traveltime (i.e. the one that satisfies at best Fermat's principle, eq.(2)). Repeating this procedure for each source-sensor pair, approximations of ℝᵏ and its adjoint, ℝᵏ*, related to the slowness field s(·, Ωₖ), are made readily available.

Figure 8 shows the ray paths for two different configurations of sources and receivers, and for three different ray path families. In our experiments, the object and the background are homogeneous mediums. As a consequence, the slowness is piecewise constant, and the raypaths should be polygonal lines. As expected, straight lines (Fig. 8c) are a poor approximation of the actual raypaths (Fig. 8a). We will see in Sec. 4 that the circular arc approximation (Fig. 8b) provides traveltimes which are close to the optimal and to the times measured on real data.

### 3.4.3. Step 2: shape evolution
To obtain the shape Ωₖ₊₁ (Step 2), we optimize (25) with respect to the shape by using a gradient descent algorithm [26]

$$ \left\{ \begin{array}{l} \frac{\partial \Gamma}{\partial \tau}(x, \tau) = -\nabla_{\Omega} J_k \triangleq V \cdot \mathcal{N} \\ \Gamma(x, 0) = \Gamma_k(x) \end{array} \right. \qquad (26) $$

where Γₖ = ∂Ωₖ is the boundary of the shape at the current iteration, V is called the speed function and N is the inward normal to ∂Ω. Applying (26) defines, at convergence, a new shape Ωₖ₊₁ which locally minimizes the shape functional Jₖ. In (26), ∇ΩJₖ is the
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Computed raypaths for two different configurations of sources and receivers, superimposed on the 2D digital model of the object, computed using (a) *Fast Marching* and back-propagation (optimal), (b) circular arc approximation (good), (c) straight line approximation (poor). Top: 5 sources and 10 receivers; Bottom: 6 sources and 16 receivers.</figcaption></figure>

shape gradient of $J_k$ which requires the calculation of the shape derivative of the data fidelity term and of the shape prior. This may be done in the Eulerian framework [26, 27] depicted in Appendix A.

The shape derivative of $J_k^{data}$, which is linear, is calculated as in [10] (see Appendix B), leading to:

$$ \nabla_{\Omega} J_k^{data} = \left( \sum_i a_i s_i^{\Omega} - \sum_i b_i s_i^{\Omega^c} \right) \mathcal{R}^{k*} [T_{obs} - \mathcal{R}^k s(\cdot, \Omega)]\mathcal{N} \quad (27) $$

In our current implementation we have considered a piecewise constant slowness model. In this case, (11) becomes:

$$ s(x, \Omega) = s_{\Omega} \mathbb{1}_{\Omega}(x) + s_{\Omega^c} \mathbb{1}_{\Omega^c}(x), \quad (28) $$

leading to the following shape derivative:

$$ \nabla_{\Omega} J_k^{data} = (s_{\Omega} - s_{\Omega^c}) \mathcal{R}^{k*} [T_{obs} - \mathcal{R}^k s(\cdot, \Omega)]\mathcal{N} \quad (29) $$

As can be seen, the resulting speed function is the product of two factors. The first one reflects the competition (in the sense of the term *Region Competition* coined by Ronfard [28]) between the background and object slowness models. The second one is akin to a back-projection of the reconstruction error. According to the sign of the product, the evolution tends to be in the direction of $\mathcal{N}$ or opposite to it. Moreover, the evolution speed vanishes when the reconstruction error becomes small enough.

Differentiating $J_{prior}$ leads to the following shape derivative of the shape prior (see Appendix C):

$$ \nabla_{\Omega} J_{prior} = - \sum_{0 \le u+v \le N} A_{uv} (H_{uv} + G_{uv}) \mathcal{N}, \quad (30) $$
where $H_{uv}$ is defined in (9), and $G_{uv}$ is given in (C.4). The $A_{uv}$ factor is defined as the weighted average of individual factors $A_{(j)uv}$:

$$A_{uv} = \left( \sum_{j=1}^{N_{ref}} A_{(j)uv} w_{(j)} \right) / (2\sigma^2 \sum_{j=1}^{N_{ref}} w_{(j)}) \quad (31)$$

The expression of $A_{(j)uv}$ (given in Appendix C) involves the distance, in terms of Legendre moments, between the evolving shape $\Omega$ and the reference shapes, $\Omega_j$. The weights, $w_{(j)}$ are exponentially decreasing functions of this distance

$$w_{(j)} = e^{\frac{-\|\lambda(\Omega) - \lambda_{(j)}^{\text{ref}}\|^2}{2\sigma^2}} \quad (32)$$

Hence, the evolution speed induced by the minimization of $J_{prior}$ is the weighted average of the forces that would attract the evolving shape towards each individual reference shape. Let us notice that monitoring the values of the weights allows the recognition of the object of interest among the reference shapes, simultaneously to its reconstruction.

In practice, the evolution of $\Omega$ according to (26) is implemented using the Level-Set algorithm, with variable step-size. A narrow-band technique is used to speed up computations [16]. The narrow band is periodically reinitialized using Sussman's equation [29]. The parameters of the slowness model are updated by solving (16) at each reinitialization of the narrow-band.

In order to stop the algorithm, we monitor the sequence of successive energies, $J_k(\Omega_k)$. The minimization is stopped when $k$ satisfies:

$$|J_{k+1}(\Omega_{k+1}) - J_k(\Omega_k)| < 0.05 J_0(\Omega_0). \quad (33)$$

## 4. Shape reconstruction results

In this section, we assess the performance of the method in different conditions. A synthetic model of the object will be used for the first tests, in § 4.1, but we will mainly use real data (§ 4.2) acquired in controlled conditions using the test bench, as described in Sec. 2.1. A schematic view of the central slice of the physical object of interest is shown on Figure 9. The materials used to make the object are specified to be homogeneous, so the piecewise constant slowness model, i.e. $\theta = (s_{\Omega}, s_{\Omega}^c)$, is physically plausible and will be considered in what follows. In § 4.3, we will present a second synthetic experiment, with a more complex object, to illustrate the capacity of the proposed approach to deal with more complex shape topologies.

### 4.1. Reconstruction from synthetic data

For the purpose of testing the behaviour of the algorithm, a digital model of the object is considered at first. We discretize the 600 mm × 350 mm central slice of the 3D object with a 5mm step size to obtain a 71 × 121 binary image which depicts a vertical half-disk over a constant background. We set the velocity of the aluminium inclusion to 6,000 m.s⁻¹ and the velocity of the resin to 3,000 m.s⁻¹, respectively. These values
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Configuration of sources and receivers (left). Traveltimes (in µs): (a) real data, (b) synthetic data computed by Fast Marching, (c) circular arc approximation, (d) straight line approximation (right).</figcaption></figure>

are chosen close to the values measured on the bench. We consider 35 locations for the source and 26 for the receiver, as shown on the left of Figure 9.

On Figure 9 (right), we show a matrix view of the synthetic traveltime data obtained for different travel path models, compared to the real traveltimes data acquired on the bench. We can notice that the traveltimes computed under the circular arc approximation [2] are close to those provided by the *Fast Marching* method, which is obviously not the case for the straight line approximation. In order to add noise to the synthetic data obtained under the circular arc approximation, we use a mixture of a Gaussian model for the noise and a uniform distribution for the outliers. The parameters of the mixture are estimated from a subset of real data, that corresponds to straight rays traversing only the resin, using an EM algorithm described in [30]. Figure 10 shows a sample of the noisy data generated using the circular arc approximation and the mixture model, which are very close to the corresponding real data. The root mean square (RMS) difference between the two series is around 4.4 µs, while the arrival times vary from 17 µs to 172 µs, with a mean of 114 µs.

The reconstruction results are shown on Figure 11. We use a *continuation* strategy that exploits the hierarchical nature of Legendre descriptors. More specifically, we first apply the proposed algorithm (as described in § 3.4), without shape prior, i.e. with $\alpha_p = 0$ (second column of Figure 11). One may notice that the reconstruction result is already close to the ideal object. Then, we start the algorithm again, initializing it with the first object-based reconstruction, introducing the shape constraint gradually by considering invariant moment representation of increasing orders, and adjusting the weight parameter $\alpha_p$ accordingly. The set of reference shapes used for this experiment contains duplicates of the object with different orientation and two other shapes, as shown on Figure 6. We successively fix the order of the invariant moment representation
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Comparison of real measurements and synthetic data (computed using (2) and corrupted by mixture noise). For clarity, only 1 traveltime out of 5 is displayed (the total number of traveltimes in this experiment is 35×26 = 910). Times (ordinates) are given in µs.</figcaption></figure>

to $N = 13$ (third column) and to $N = 27$ (fourth column), while the weight parameter is set to $\alpha_p = 0.3$ and $\alpha_p = 0.9$, respectively. The reconstruction results (fifth column) are visually satisfactory. The estimated velocities are very close to the specified values (see Table 1): the relative error is about 5 % for the aluminium and 2 % for the resin.

On Figure 12, we present the evolution of the reconstruction energy $J_k(\Omega)$ (eq. 25) along the iterations of the algorithm. The blue curves correspond to the shape evolutions by gradient descent (Step 2 of the algorithm). We use red circles to identify the values of $J_k(\Omega_k)$ obtained after updating $\mathcal{R}^k$ (Step 1). Even if this operation induces energy increases, these jumps tend to become smaller at each iteration $k$ and one can check that the series $(J_k(\Omega_k))_{k \ge 0}$ is decreasing. Let us notice that we obtained similar curves in all our experiments.

<table><caption>**Table 1.** Details of the constrained object-based reconstruction from noisy numerical data for 35 source locations and 26 receiver locations. Recall that the specified velocities are 6,000 m.s$^{-1}$ for the aluminium inclusion and 3,000 m.s$^{-1}$ for the resin.</caption><thead><tr><th>Number of steps</th><th>N</th><th>α<sub>p</sub></th><th>Velocity of the aluminium (m.s<sup>-1</sup>)</th><th>Velocity of the resin (m.s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Initial values</td><td>-</td><td>-</td><td>3,582</td><td>3,365</td></tr><tr><td>4</td><td>0</td><td>0</td><td>4,987</td><td>3,056</td></tr><tr><td>4</td><td>13</td><td>0.3</td><td>5,600</td><td>3,053</td></tr><tr><td>2</td><td>27</td><td>0.9</td><td>5,685</td><td>3,052</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Constrained object-based reconstruction from noisy numerical data for 35 source locations and 26 receiver locations: reconstructed velocity field (top); evolving contour (red) superimposed on the ideal object shown in yellow (bottom).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Energy curves (eq. 25) for the constrained object-based reconstruction from noisy numerical data for 35 source locations and 26 receiver locations. Blue curves show the energy $J_k(\Omega)$ minimized by gradient descent (step 2) and red circles show the values of $J_k(\Omega_k)$ after each update of $\mathcal{R}_k$ (step 1).</figcaption></figure>

## 4.2. Reconstruction from real data

We now consider shape reconstruction from real data. For comparison purposes, we first present the results of a pixel-based SIRT algorithm [2], which shares the same ray path model (circular arc approximation) as the proposed algorithm. Figure 13 (a) shows the reconstruction of the velocity field. As can be seen, the pixel-based method yields poor quality reconstructions.

On Figure 14, we show the results of our constrained shape evolution method for real data. As previously (Fig. 11), we first perform an optimization for $\alpha_p = 0$. The
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Left: Pixel-based reconstructions of the velocity field (in $m.s^{-1}$) from real data, using a SIRT algorithm. Right: the boundary of the ideal shape (in yellow) is presented, superimposed on the reconstruction (gray scale). Standard setup (a). Reduced data setup (b).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Constrained object-based reconstruction from real data for 35 source locations and 26 receiver locations: reconstructed velocity field (top); evolving contour (red) superimposed on the ideal object shown in yellow (bottom).</figcaption></figure>

<table><caption>**Table 2.** Details of the constrained object-based reconstruction from real data for 35 source locations and 26 receiver locations.</caption><thead><tr><th>Number of steps</th><th>N</th><th>α<sub>p</sub></th><th>Velocity of the aluminium (m.s<sup>-1</sup>)</th><th>Velocity of the resin (m.s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Initial values</td><td>-</td><td>-</td><td>3,522</td><td>3,330</td></tr><tr><td>3</td><td>0</td><td>0</td><td>5,081</td><td>3,039</td></tr><tr><td>2</td><td>13</td><td>0.7</td><td>5,390</td><td>3,025</td></tr><tr><td>2</td><td>27</td><td>0.8</td><td>5,500</td><td>3,024</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Constrained object-based reconstruction for the reduced data setup (real data, 18 source location and 13 receiver locations): reconstructed velocity field (top); evolving contour (red) superimposed on the ideal object shown in yellow (bottom).</figcaption></figure>

result is shown on the second column of Figure 14. Then, we introduce the shape prior. We first adjust the weight parameter $\alpha_p$ in (25) for a moment-based representation up to the order $N = 13$. The result is shown on the third column of Figure 14. The final result, obtained for $N = 27$ is shown on the fourth column. As can be observed, the proposed method gives good results, both in terms of shape reconstruction and estimated velocity values (Table 2). The energy curves (which are not shown here) are similar to those obtained in the previous experiment.

### 4.2.1. Reduced data setup
To assess the method when a significant portion of data is missing, we retain only one source/sensor out of two, in a second experimental setup. The result for the pixelwise SIRT approach appears Figure 13 (b). As expected, it is of low quality. One can see that the constrained shape reconstruction method (using the same values for $N$ and $\alpha_p$ as for the full data set) is able to recover properly the shape (Figure 15) and the velocity values (Table 3) even in such adverse conditions.

### 4.2.2. Handling rotations
In the above experiments, the ideal object belongs to the set of reference shapes up to a scaling and translation. In the following experiment, performed with the standard setup (35 sources, 25 receivers), we consider rotations (from -45° to +45° every 2°, see figure 16), such that the actual object does not belong to the set of reference shapes. As can be seen on figure 17, the reconstruction only slightly degrades.

These results show that it is possible to take into account rotations through the
<table><caption>**Table 3.** Details of the constrained object-based reconstruction for the reduced real data setup (real data, 18 source location and 13 receiver locations).</caption><thead><tr><th>Number of steps</th><th>N</th><th>α<sub>p</sub></th><th>Velocity of the aluminium (m.s<sup>-1</sup>)</th><th>Velocity of the resin (m.s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Initial values</td><td>-</td><td>-</td><td>3,526</td><td>3,311</td></tr><tr><td>3</td><td>0</td><td>0</td><td>4,818</td><td>3,028</td></tr><tr><td>3</td><td>13</td><td>0.7</td><td>5,237</td><td>2,992</td></tr><tr><td>2</td><td>27</td><td>0.8</td><td>5,238</td><td>2,993</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** 12 reference shapes corresponding to rotations of the ideal object up to scaling and translation (from -45° to +45° every 2°).</figcaption></figure>

alphabet of reference. Indeed, while Legendre moments can be made invariant to rotation, we showed in [21] that this makes the optimization more computationally intensive and that it complexifies the energetic landscape (the number of local minima). Accounting for rotations via the set of reference shapes is an interesting alternative because this can be done at a constant complexity, and because potential ambiguities between rotated reference shapes can be avoided.

Moreover, this experiment illustrates the behaviour of the algorithm when the ideal object does not belong to the set of reference shapes. On Figure 17, one can see that, for high values of $\alpha_p$, the solution is biased towards the closest shape which is 1° away from the ideal solution. By decreasing the value of the regularization parameter, the data term is given more importance and the solution is closest to the observed shape, at the price of a lower quality.

### 4.2.3. Impact of the set of reference shapes
To further illustrate the impact of the reference alphabet, we show on Figure 18 a set of experimental results obtained from real data (35 sources and 26 receivers) with different sets of reference shapes. More specifically, we add or substitute ambiguous shapes to the alphabet of Figure 6: a crescent and a crude polygonal approximation of the half-disk.

On the three first lines of the table in Figure 18, it can be seen that when the ideal shape is present in the reference alphabet, the reconstruction is constrained towards it,
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Boundary of the recovered shape (red) superimposed on the ideal object (yellow) for the constrained object-based reconstruction for the set of reference shapes of Fig. 16 (successive rotations of angle 2°).</figcaption></figure>

as expected. However, if we substitute the crescent or the polygon to the ideal shape (fourth and fifth lines), the solution evolves towards the lure, which corresponds to the closest approximation to the ideal object. Note that the same set of parameters was used for all these experiments.

### 4.3. Dealing with topological changes

We finally present an illustration of the capacity of the method to deal with more complex object topologies. The synthetic object considered here is shown on Figure 19. It is made up of two components: a square with an ellipsoid hole and a disk located in the top-right part of the image. As for the previous series of experiments, the speed value for the object and for the background are respectively set to 6,000 m.s⁻¹ and 3,000 m.s⁻¹. We consider 35 source locations and 40 receiver locations, as shown on Figure 19, left. The sets of raypaths shown on Figure 20, give an idea of the coverage of the domain. A matrix view of the simulated traveltimes, for the complete set of ray paths, is given in Figure 19 (right), in false colors. In this case, we use additive Gaussian noise, with RMS = 2.8 µs. We first perform pixel-based reconstruction using the SIRT algorithm described in [2]. The best result we obtained with this method remains of a poor quality, as shown on Figure 21.

We next apply the object-based method proposed in the present paper, using the
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Reconstruction results from real data (35 sources and 26 receivers), for different sets of reference shapes, with lures. When the ideal shape belongs to the reference set, the algorithm converges towards it; else, it converges to the closest shape (polygon or crescent).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Experimental setup (left) and noisy numerical traveltimes, in µs (right).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Circular ray paths for two sample sets of source/receiver pairs: left, 18 sources and 13 receivers; right, 18 sources and 7 receivers.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 21.** Pixel-based reconstruction with the SIRT algorithm (a), boundaries of the object superimposed on the reconstructions (b).</figcaption></figure>

same continuation strategy as previously: without shape constraint, then with shape prior for increasing orders of Legendre moments, $N$. Figure 22 shows the first 16 iterations of the reconstruction process without shape prior (i.e. $\alpha_p = 0$), for $k = 1$. As can be seen, the initial circular tends to split, in order to capture the two components of the object. However, the separation is incomplete at the end of this reconstruction stage
<figure><img src="image_23.png" /><figcaption>**Figure 22.** Evolution (first 15 iterations) of the shape reconstructed from numerical data with a Gaussian additive noise of 20 dB, without shape prior ($\alpha_p$ = 0) and for $k = 1$ (red curve), superimposed on the ideal object (yellow).</figcaption></figure>

(Figure 24, second column). Then, we introduce the shape constraint by considering the set of reference shapes presented on Figure 23, for $N = 21$ and then, for $N = 45$. The reconstruction results are shown on Figure 24. As can be seen, the two components of the object are already recovered at the end of the second stage ($N = 21$). The last stage only brings minor improvements. Table 4 shows that the speed constants are nicely estimated.

## 5. Conclusion

In this paper, we have investigated the case of nonlinear shape reconstruction from traveltime data. A shape evolution method based on the resolution of a series of linear optimization problems, and incorporating high-level shape constraints, has been proposed. The method simultaneously performs shape reconstruction and
<figure><img src="image_24.png" /><figcaption>**Figure 23.** Set of reference shapes.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 24.** Constrained object-based reconstruction results from numerical data with a Gaussian additive noise of 20 dB: reconstructed velocity field (top), evolving contour, in red, superimposed on the ideal object, in yellow (bottom).</figcaption></figure>

<table><caption>**Table 4.** Details of the constrained object-based reconstruction of Figure 24 from numerical data with a Gaussian additive noise of 20 dB for 35 source locations and 40 receiver locations. Recall that the modeled velocities are 6,000 m.s⁻¹ for the aluminium inclusion and 3,000 m.s⁻¹ for the resin.</caption><thead><tr><th>Number of steps</th><th>N</th><th>α<sub>p</sub></th><th>Velocity of the aluminium (m.s<sup>-1</sup>)</th><th>Velocity of the resin (m.s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Initial values</td><td>-</td><td>-</td><td>3,584</td><td>3,466</td></tr><tr><td>5</td><td>0</td><td>0</td><td>5,536</td><td>2,995</td></tr><tr><td>3</td><td>21</td><td>0.5</td><td>5,897</td><td>2,996</td></tr><tr><td>3</td><td>45</td><td>0.9</td><td>5,992</td><td>2,988</td></tr></tbody></table>

recognition. It leads to promising reconstruction results on both numeric and real data, even for reduced sets of data. The method handles non-trivial topologies by exploiting the topological change capacities of the level-sets implementation. These
results demonstrate the feasibility and suitability of introducing shape constraints into traveltime tomographic reconstruction.

Our implementation of the direct problem is certainly perfectible. In particular, using the Fast Marching [24] technique would be a straightforward extension of this work. Another interesting direction would be the introduction of an adjoint-state solver, similar to [15]. More sophisticated physical models, including multiple arrival times, or even the full waveform [31], could be considered as well, at the price of an increased complexity. We note that the proposed model is limited to the reconstruction of a single object, belonging to a set of known reference shapes. Extending it to the case of multiple independent objects, as well as to a non-uniform background, is another direction for our future work. Besides, even if we observed a decrease of the energy in all our experiments, with real and synthetic data, the theoretical study of the convergence of the algorithm and the characterization of its solutions remains a non-trivial and open issue. The application of this kind of methods to complex data acquired in real site conditions is another future prospect.

## Appendix A. Shape gradients of domain integrals (Eulerian framework)

The framework of shape derivative [26, 27] is based on shape variations $\Omega_{\tau}$ from a shape $\Omega$. We model shape variations by the speed method introduced in [26], namely $\Omega_{\tau}$ is the set of $x(\tau) \in \mathcal{D}$ satisfying

$$ \left\{ \begin{array}{l} \frac{dx}{d\tau} = \mathcal{V}(x(\tau)) \quad \text{for } \tau > 0 \\ x(0) \in \Omega \end{array} \right. \qquad (A.1) $$

For the conditions on the function $\mathcal{V}$ to ensure that the shape $\Omega_{\tau}$ remains contained into the whole domain $\mathcal{D}$, we refer to [26]. Then, the shape derivative of a shape functional is defined as follows:

**Definition 1** Let $J$ a functional depending on the shape $\Omega$. Then, if it exists, the shape derivative of $J$ in the direction $\mathcal{V}$ is given by

$$ J'(\Omega) \cdot \mathcal{V} = \lim_{\tau \searrow 0} \frac{J(\Omega_{\tau}) - J(\Omega)}{\tau} \qquad (A.2) $$

For a domain integral, i.e. a shape functional $J$ of the form

$$ J(\Omega) = \int_{\Omega} f(x, \Omega)dx, \qquad (A.3) $$

one has the following theorem [26]:

**Theorem 1** The shape derivative of a domain integral $J = \int_{\Omega} f(x, \Omega)$ is given by

$$ J'(\Omega) \cdot \mathcal{V} = \int_{\Omega} f'(x, \Omega) \cdot \mathcal{V}dx - \int_{\Gamma} f(x, \Omega)(\mathcal{V}(x), \mathcal{N}(x))d\Gamma(x) \qquad (A.4) $$

where $(\cdot, \cdot)$ is the usual dot product in $\mathbb{R}^2$, $\Gamma = \partial\Omega$ is the boundary of $\Omega$ and $\mathcal{N}$ is the inward normal to $\Gamma$.
When the integration holds on $\Omega^c$, a similar formula is obtained, except that the integration on $\Gamma$ has a minus sign.

Furthermore, related to the shape derivative, one can define the shape gradient as in [26]. The shape gradient is supported on $\Gamma = \partial\Omega$ and depends on the considered inner product. We provide the shape gradient of a rather general class of shape derivative in the next definition.

**Definition 2** Let the shape derivative of $J$ be of the form

$$J'(\Omega) \cdot \mathcal{V} = \int_{\Gamma} g(x) (\mathcal{V}(x), \mathcal{N}(x)). \qquad (A.5)$$

Then, considering the $L^2$ inner product, the shape gradient $\nabla_{\Omega} J$ of $J$ is defined as

$$\nabla_{\Omega} J = -g\mathcal{N} \qquad (A.6)$$

In [32], other inner products are considered, which leads to other forms of the shape gradient.

## Appendix B. Shape derivative of the data fidelity term

In this Appendix, we calculate the shape derivative of

$$J_k^{\text{data}}(\Omega) = \frac{1}{2} \| T^{\gamma^k}(\cdot, \cdot, s(\cdot, \Omega)) - T_{\text{obs}} \|^2 = \frac{1}{2} \| \mathcal{R}^k s(\cdot, \Omega) - T_{\text{obs}} \|^2 = \frac{1}{2} \| \mathcal{M}^k \theta - T_{\text{obs}} \|^2.$$

One may notice that, in this expression, $\theta$ is, itself, a function of $\Omega$. Therefore, a term of the form $(\partial J_k^{\text{data}}(\hat{\theta})/\partial\theta).(\theta'(\Omega) \cdot \mathcal{V})$ should appear in the shape derivative of $J_k^{\text{data}}$. However, as noticed in [33], this term vanishes due to the necessary optimality conditions satisfied by $\hat{\theta}$. So, what remains to be calculated is the shape derivative of $J_k^{\text{data}}$ for $\theta$ fixed. Omitting parameters for clarity, we write $J_k^{\text{data}}$ as:

$$J_k^{\text{data}} = \frac{1}{2} \sum_{j=1}^{j=n_s \times n_r} (T_{\text{obs}} - \mathcal{R}^k s)_j^2, \qquad (B.1)$$

$\mathcal{R}^k s = s_{\Omega} \mathcal{R}^k 1_{\Omega} + s_{\Omega^c} \mathcal{R}^k 1_{\Omega^c}$, the shape derivative of $J_k^{\text{data}}$ is given by

$$\nabla_{\Omega} J_k^{\text{data}} = - \sum_{j=1}^{j=n_s \times n_r} (T_{\text{obs}} - \mathcal{R}^k s)_j \left( \sum_{i=1}^{i=K} a_i \nabla_{\Omega} \mathcal{A}_{i,j}^k + \sum_{i=1}^{i=L} b_i \nabla_{\Omega} \mathcal{B}_{i,j}^k \right) \qquad (B.2)$$

where we recall (see Sec. 3.2) that $\mathcal{A}_{ij}$ is the $j$-th component of $\mathcal{R}^k s_i^{\Omega^c}$:

$$\mathcal{A}_{ij} = \int_{\gamma_j^k} s_i^{\Omega} 1_{\Omega}. \qquad (B.3)$$

and $\mathcal{B}_{ij}$ is the $j$-th component of $\mathcal{R}^k s_i^{\Omega^c}$. If we denote by $l_j^k(x)$ the projection kernel [7] that defines ray path $\gamma_j^k$, then we may write $\mathcal{A}_{ij}$ as a domain integral:

$$\mathcal{A}_{ij} = \int_{\Omega^k} s_i^{\Omega}(x) l_j^k(x) dx, \qquad (B.4)$$
Since $l_j^k(x)$ does not depend on $\Omega$, applying Theorem 1 and Definition 2, we have:

$$ \nabla_{\Omega} A_{ij} = -s_i^{\Omega} l_j^k \mathcal{N}. \qquad (B.5) $$

Similarly, we have that

$$ \nabla_{\Omega} B_{ij} = s_i^{\Omega^c} l_j^k \mathcal{N}, \qquad (B.6) $$

where the minus sign indicates that the normal to $\Omega^c$ is opposite to $\mathcal{N}$. Collecting these results, we obtain:

$$ \nabla_{\Omega} J_{k}^{\text{data}} = \left( \sum_{i=1}^{i=K} a_i s_i^{\Omega} - \sum_{i=1}^{i=L} b_i s_i^{\Omega^c} \right) \sum_{j=1}^{j=n_s \times n_r} (T_{\text{obs}} - \mathcal{R}^k s)_j l_j^k \mathcal{N} \qquad (B.7) $$

where the last sum may be interpreted as the back-projection of $(t - \mathcal{R}^k s)$ to the image domain, i.e.:

$$ \nabla_{\Omega} J_{k}^{\text{data}} = \left( \sum_{i=1}^{i=K} a_i s_i^{\Omega} - \sum_{i=1}^{i=L} b_i s_i^{\Omega^c} \right) \mathcal{R}^{k*} (T_{\text{obs}} - \mathcal{R}^k s) \mathcal{N} \qquad (B.8) $$

## Appendix C. Shape derivative of the shape prior

In this appendix, we sketch the calculus of the shape derivative of the prior, $J_{\text{prior}}$, defined in (23) (more details may be found in [34]). Applying the chain rule, it is sufficient to determine the shape derivative of the invariant Legendre moments $\lambda_{p,q}(\Omega)$. Moreover, since expansion (19) is linear, the shape derivative of $\lambda_{p,q}(\Omega)$ is straightforwardly obtained by the calculation of the shape derivative of the invariant regular moments $\eta_{u,v}$, for all $u$ and $v$ such that $0 \le u \le p, 0 \le v \le q$. We recall that

$$ \eta_{u,v}(\Omega) = \int_{\Omega} H_{u,v}(x, \Omega) dx \qquad (C.1) $$

where $H_{u,v}(x, \Omega)$ is given by (22). Then, using Theorem 1, one has:

$$ \eta'_{u,v}(\Omega) \cdot \mathcal{V} = \int_{\Omega} H'_{u,v}(x, \Omega) \cdot \mathcal{V} dx - \int_{\Gamma} H_{u,v}(x, \Omega)(\mathcal{V}(x), \mathcal{N}(x)) d\Gamma(x) \qquad (C.2) $$

where $H'_{u,v}(x, \Omega)$ is the derivative of $H_{u,v}(x, \Omega)$ with respect to the shape $\Omega$. To perform the calculation of $H'_{u,v}(x, \Omega)$, one has to note that $H_{u,v}(x, \Omega)$ depends on $\Omega$ through the components $\bar{x}_1$ and $\bar{x}_2$ of the centroid and through $|\Omega|$, the area of the shape $\Omega$. Using Theorem 1 and the chain rule, we get

$$ \nabla_{\Omega} J_{\text{prior}} = \sum_{0 \le u+v \le N} A_{u,v} (H_{u,v}(x, \Omega) + G_{u,v}(x)) \mathcal{N}(x) \qquad (C.3) $$

with $\mathcal{N}$ the outward normal to $\Gamma = \partial\Omega$ and where $G_{u,v}$ can be calculated by the shape differentiation of $\bar{x}_1, \bar{x}_2$ and $|\Omega|$ with respect to $\Omega$. We obtain:

$$ G_{u,v}(x) = \frac{u(x_1 - \bar{x}_1)}{S_F^{1/2} |\Omega|^{3/2}} \eta_{u-1,v} + \frac{v(x_2 - \bar{x}_2)}{S_F^{1/2} |\Omega|^{3/2}} \eta_{u,v-1} + \frac{u+v+2}{2|\Omega|} \eta_{u,v} \qquad (C.4) $$
The factor $A_{u,v}$ is the weighted average

$$A_{u,v} = \left( \sum_{j=1}^{N_{ref}} A_{(j)uv} w_{(j)} \right) / \left( 2\sigma^2 \sum_{j=1}^{N_{ref}} w_{(j)} \right), \qquad (C.5)$$

of the expansion

$$A_{(j)uv} = 2 \sum_{0 \le p+q \le N} \left( \lambda_{p,q} - \lambda_{(j)p,q}^{\text{ref}} \right) C_{pq} a_{pu} a_{qv} \qquad (C.6)$$

with weights

$$w_{(j)} = \exp \left( - \frac{\| \boldsymbol{\lambda}(\Omega) - \boldsymbol{\lambda}_{(j)}^{\text{ref}} \|^2}{2\sigma^2} \right). \qquad (C.7)$$

The weights are exponentially decaying functions of the distance (in terms of invariant Legendre moments) between the evolving shape $\Omega$ and the $j$-th reference shape.

## References

[1] N. Rawlinson and M. Sambridge. Seismic traveltime tomography of the crust and lithosphere. *Advances in Geophysics*, 46:81–198, 2003.

[2] O. Abraham, K. Ben Slimane, and Ph. Côte. Factoring anisotropy into iterative geometric reconstruction algorithms for seismic tomography. *International Journal of Rock Mechanics and Mining Sciences*, 35(1):31–41, 1998.

[3] J. Charléty, S. Voronin, G. Nolet, I. Loris, F. Simons, K. Sigloch, and I. Daubechies. Global seismic tomography with sparsity constraints: Comparison with smoothing and damping regularization. *Journal of Geophysical Research: Solid Earth*, 118(9):4887–4899, 2013.

[4] G. Aubert and P. Kornprobst. *Mathematical Problems in Image Processing: Partial Differential Equation and Calculus of Variations*, volume 147 of *Applied Mathematical Sciences*. Springer-Verlag, 2nd edition, 2006.

[5] P. Charbonnier, L. Blanc-Féraud, G. Aubert, and M. Barlaud. Deterministic edge-preserving regularization in computed imaging. *IEEE Transactions on Image Processing*, 6(2):298–311, February 1997.

[6] F. Santosa. A level-set approach for inverse problems involving obstacles. *The European Series in Applied and Industrial Mathematics: Control, Optimization and Calculus of Variations*, 1:17–33, 1996.

[7] H. Feng, W. Karl, and D. Castañon. A curve evolution approach to object-based tomographic reconstruction. *IEEE Transactions on Image Processing*, 12(1):44–57, January 2003.

[8] O. Dorn and D. Lesselier. Level set methods for inverse scattering. *Inverse Problems*, 22(4): R67–R131, 2006.

[9] D. Mumford and J. Shah. Optimal approximations by piecewise smooth functions and associated variational problems. *Com. pure & Appli. Math.*, 42(5):577–685, 1989.

[10] G. Gaullier, P. Charbonnier, and F. Heitz. Introducing shape priors in object-based tomographic reconstruction. In *ICIP 2009, 16th International Conference on Image Processing*, pages 1077–1080. IEEE, 2009.

[11] J.-P. Bruandet, F. Peyrin, J.-M. Dinten, and M. Barlaud. 3D tomographic reconstruction of binary images from cone beam projections: a fast level set approach. In *IEEE International Symposium on Biomedical Imaging*, pages 677–680, 2002.

[12] C. Alvino and A. Yezzi. Tomographic reconstruction of piecewise smooth images. In *IEEE Computer Society Conference on Computer Vision and Pattern Recognition*, pages 576–581, 2004.
[13] A. Litman, D. Lesselier, and F. Santosa. Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set. *Inverse Problems*, 14(3): 685–706, 1998.

[14] M. Burger. A level set method for inverse problems. *Inverse Problems*, (17):1327–1355, 2001.

[15] W. Li and S. Leung. A fast local level set adjoint state method for first arrival transmission traveltime tomography with discontinuous slowness. *Geophysical Journal International*, 195(1):582–596, October 2013.

[16] S. Osher and J.A. Sethian. Fronts propagating with curvature-dependant speed: algorithms based on Hamilton-Jacobi formulations. *Journal of Computational Physics*, 79(1):12–49, November 1988.

[17] Y. Lin and A. Ortega. Object-based high contrast travel time tomography. *Submitted to SEG Geophysics*, arXiv:1303.3052 [physics.data-an], 2013.

[18] F. Bretaigneau, D. Leparoux, O. Durand, and O. Abraham. Small-scale modeling of onshore seismic experiment: A tool to validate numerical modeling and seismic imaging methods. *Geophysics*, 76(5):T101–T112, January 2011.

[19] V. Červený. *Seismic Ray Theory*. Cambridge University Press, Cambridge, 2001.

[20] F. Natterer. *The Mathematics of Computerized Tomography*. Teubner, Stuttgart, 1986.

[21] A. Foulonneau, P. Charbonnier, and F. Heitz. Multi-reference shape priors for active contours. *International Journal on Computer Vision*, 81(1):68–81, January 2009.

[22] A. Foulonneau, P. Charbonnier, and F. Heitz. Affine-invariant geometric shape priors for region-based active contours. *IEEE Transactions On Pattern Analysis and Machine Intelligence*, 28(8):1352–1357, August 2006.

[23] D. Cremers, S. Osher, and S. Soatto. Kernel density estimation and intrinsic alignment for shape priors in level set segmentation. *International Journal of Computer Vision*, 69(3):335–351, 2006.

[24] J.A. Sethian. A fast marching level set method for monotonically advancing fronts. *Proc. Nat. Acad. Sci.*, 93(4):1591–1595, 1996.

[25] R. Tsai Y.-H., L.T. Cheng, S. Osher, and H.-K. Zhao. Fast sweeping algorithms for a class of hamilton-jacob equations. *SIAM Journal on Numerical Analysis*, 2(41):673–694, 2003.

[26] J. Sokolowski and J.P. Zolesio. *Introduction to shape optimization: shape sensitivity analysis*, volume 16 of Springer Series in Computational Mathematics. Springer Verlag, Berlin, Heidelberg, New York, July 1992.

[27] G. Aubert, M. Barlaud, O. Faugeras, and S. Jehan-Besson. Image segmentation using active contours: calculus of variations or shape gradients? *SIAM, Journal on Applied Mathematics*, 63(6):2128–2154, September 2003.

[28] R. Ronfard. Region based strategies for active contour models. *International Journal of Computer Vision*, 13(2):229–251, October 1994.

[29] M. Sussman, P. Smereka, and S. Osher. A level set approach for computing solutions to incompressible two-phase flow. *Journal of Computational Physics*, 114:146–159, 1994.

[30] P. Torr and A. Zisserman. Mlesac: a new robust estimator with application to estimating image geometry. *Computer Vision and Image Understanding*, 78(1):138–156, April 2000.

[31] F. Bretaigneau, R. Brossier, D. Leparoux, O. Abraham, and J. Virieux. 2d elastic full-waveform imaging of the near-surface: application to synthetic and physical modelling data sets. *Near Surface Geophysics*, 11(3):307–316, 2013.

[32] M. Burger. A framework for the construction of level set methods for shape optimization and reconstruction. *Interfaces and Free Boundaries*, 5:301–329, 2003.

[33] R. Ramlau and W. Ring. A Mumford-Shah level-set approach for the inversion and segmentation of X-ray tomography data. *Journal of Computational Physics*, 221(2):539–557, 2007.

[34] A. Foulonneau, P. Charbonnier, and F. Heitz. Multi-reference affine-invariant geometric shape priors for region-based active contours. Technical Report RR-AF01-08, LRS ERA 27 LCPC / LSIIT UMR 7005 CNRS, May 2008. An abridged version (without the appendices) of this preliminary report was published in [21].