*Article* **Image Segmentation with a Priori Conditions: Applications to Medical and Geophysical Imaging**

**Guzel Khayretdinova 1,2,†, Christian Gout 2,\*,†, Théophile Chaumont-Frelet 3,† and Sergei Kuksenko 1,†**

<sup>2</sup> INSA Rouen, LMI, 76000 Rouen, France


**Abstract:** In this paper, we propose a method for semi-supervised image segmentation based on geometric active contours. The main novelty of the proposed method is the initialization of the segmentation process, which is performed with a polynomial approximation of a user defined initialization (for instance, a set of points or a curve to be interpolated). This work is related to many potential applications: the geometric conditions can be useful to improve the quality the segmentation process in medicine and geophysics when it is required (weak contrast of the image, missing parts in the image, non-continuous contour. . . ). We compare our method to other segmentation algorithms, and we give experimental results related to several medical and geophysical applications.

**Keywords:** image segmentation; a priori segmentation; level set method; geodesic active contour

#### **1. Introduction**

The problem of segmenting an image into its significant components has been studied for over 30 years in computers science, applied mathematics and more generally in computer vision.

Recently, from convolutional neural networks (see Fukushima [1], Waibel et al. [2], and LeCun et al. [3], among others. . . ) to recurrent neural networks (Hochreiter and Schmidhuber [4]), encoder–decoders (Badrinarayanan et al. [5]), or generative adversarial networks (Goodfellow et al. [6]), deep learning based techniques have achieved huge successes in the field of artificial intelligence and image segmentation (see [7] for a recent survey). These approaches leads to excellent results, including in medical applications (see for example Fantazzini et al. [8]). However, in general, the performance heavily depends on labeled data, and this point is a main difficulty on several applications when lacking labeled data such as in many medical and geophysical applications. The data augmentation is possible [9], but it is often complicated, and requires more CPU/GPU time than other segmentation techniques based on energy minimization or variational approaches (Kass et al. [10], Chand and Vese [11], Mumford and Shah [12], Vese and Le Guyader [13]. . . ) or geometrical ones (fast marching methods, see Sethian [14] or Forcadel et al. [15]).

In this work, we focus on a modeling based on a variational approach: it basically consists of evolving an initial contour subject to constraints towards the boundary of the object to be detected. This deformation is done by minimizing a functional depending on the curve and defined so that a local minimum is obtained at the boundary of the object. This energy-like functional minimization problem has led to many research works in the last 25 years. Caselles, Kimmel and Sapiro [16] have shown that by setting one of the regularization parameters to zero in the classical active contour model, one gets a problem equivalent to finding a geodesic curve in a Riemann space whose metric depends on the image contents. The issue was then no longer seen as an energy-like minimization problem

**Citation:** Khayretdinova, G.; Gout, C.; Chaumont-Frelet, T.; Kuksenko, S. Image Segmentation with a Priori Conditions: Applications to Medical and Geophysical Imaging. *Math. Comput. Appl.* **2022**, *27*, 26. https:// doi.org/10.3390/mca27020026

Academic Editor: Gianluigi Rozza

Received: 3 March 2021 Accepted: 7 March 2022 Published: 11 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

97

<sup>1</sup> TUSUR CSR, 634050 Tomsk, Russia; guzel.khayretdinova@insa-rouen.fr (G.K.); ksergp@tu.tusur.ru (S.K.)

but as a curve evolution one. Let us note that this approach is very suitable to our problems, and it could be possible to use deep learning (DL) based algorithm (as used by Le Guyader et al. [17]), but the a priori conditions are complicated to integrate in such DL models (and such DL models require a lot of hardware, much more that our proposed algorithm).

An edge can be viewed as the locus of connected points for which the image gradient varies abruptly. However, when data acquisition cannot be performed in an optimal manner (e.g., liver in medical imaging), this criterion can no longer be applied. In many applications, image data are missing or of poor quality, and occultation phenomena can often occur: two objects that are adjacent to one another can own intrinsic homogeneous textures so that it is hard to clearly identify the interface between them. It is then relevant to introduce geometric constraints in the modeling to help the segmentation process. In [18–21], the authors introduce the concept of a finite set of points to be interpolated during the segmentation process (like in Gout and Le Guyader [20]) or a variational approach to integrate higherlevel shape priors into level set based segmentation method (see Cremers et al. [18]). However, in some applications (in geophysics or in several medical applications), this is not sufficient and adding new a priori conditions to be satisfied is required both to improve the segmentation process and also to simplify the initial condition (in [20], the method requires a priori conditions plus an initial guess to start the segmentation process, which constitutes a drawback of their algorithm).

In this work, we propose to consider a curve (or a set of curves) belonging to the searched contour of interest, integrated as an initial condition to be satisfied by the final segmentation contour obtained at the end of the segmentation process. This is an added value considering previous models (like [19]), where it was necessary to both give an initial condition and geometric conditions. In order to define this curve, it is possible to use *Dm*-spline functions (see Gout et al. [22]) obtained by minimizing an energy functional on a suitable Hilbert space (and so a set of points given by the user is sufficient to create such geometric conditions).

In this paper, we generalize and improve the approach of Gout and Le Guyader [20]: we combine the curve evolution approach developed by Caselles, Kimmel and Sapiro [16], with a geometrical approach consisting in a curve given by the user (a set of points can also be considered like in Gout and Le Guyader [20]), and the initial condition is then automatically generated (from the geometric conditions).

In Section 2, we give the modeling of our problem, based on a level set method (LSM) approach (Osher and Sethian [23], see also [24]). Using Euler-Lagrange variational principle, we obtain the PDE satisfied by the 3D LSM function Φ and the corresponding evolution problem is then given. The existence and uniqueness of the viscosity solution of the associated parabolic problem is also established: we give the main steps of the proof (based on Caselles et al. [16] and Gout and Le Guyader [20]). The discretization of the problem is tackled in Section 3, and is solved using an Additive Operator Splitting scheme (Weickert and Kühne [25]). Comparisons and experimental results with the proposed approach are given for several medical and geophysical imaging applications.

#### **2. Mathematical Modelling**

The well-known level set approach (see Osher and Sethian [23]) consists in considering the evolving active contour Γ = Γ(*t*) as the zero level set of a function Φ, which is a Lipschitz continuous function defined by:

$$\begin{cases} \Phi: & \Omega \times [0, +\infty[ \longrightarrow \longrightarrow \mathbb{R}, \\ & (x, t) \qquad \longmapsto \quad \Phi(x, t), \end{cases} \tag{1}$$

such that:

$$\Gamma(t) = \{ \mathbf{x} \in \Omega \: \mid \: \Phi(\mathbf{x}, t) = 0 \}, \tag{2}$$

and Φ(·, *t*) takes opposite signs on each side of Γ(*t*). This level set approach enables us to re-write the energy in terms of Φ:

$$E(\Phi) = \int\_{\Omega} g(|\nabla I(\mathbf{x})|) |\nabla H(\Phi(\mathbf{x}))| d\mathbf{x},\tag{3}$$

where *<sup>H</sup>* is the one-dimensional Heaviside function and where *<sup>g</sup>*(*s*) = <sup>1</sup> <sup>1</sup> <sup>+</sup> *<sup>s</sup>*<sup>2</sup> . We formulate the shape optimization problem as follows: "Find a shape such that the energy *E* is minimized":

$$\begin{cases} \text{Search for } \Phi \text{ such that:}\\ E(\Phi) = \min\_{\vec{\xi}} E(\vec{\xi}). \end{cases} \tag{4}$$

The functional *E* being not Gâteaux-differentiable, we regularize the problem by considering slightly regularized versions (*C*<sup>1</sup> or *C*<sup>2</sup> regularization) of the functions *H* and *δ* (the one dimensional Dirac measure) denoted, respectively, by *H* and *δ*, and considering that *δ* = *H* , the associated regularized functional *E* becomes:

$$E\_{\mathbf{c}}(\Phi) = \int\_{\Omega} g(|\nabla I(\mathbf{x})|) \delta\_{\mathbf{c}}(\Phi(\mathbf{x})) |\nabla \Phi(\mathbf{x})| dx,\tag{5}$$

where:

$$\int\_{\Omega} \delta\_{\mathfrak{c}}(\Phi(\mathfrak{x})) |\nabla \Phi(\mathfrak{x})| dx \tag{6}$$

is an approximation of the length of the zero level set of Φ.

#### *2.1. A Priori Conditions*

We now propose to include a priori condition in the process. We consider the a priori condition as a smooth arc curve C given by the user. The curve C is parametrised by *<sup>c</sup>* <sup>∈</sup> *<sup>C</sup>*1([0, 1], <sup>Ω</sup>). By construction, we impose that the curve <sup>C</sup> belongs to the level set {Φ<sup>0</sup> = 0}. The first step consists in constructing the initial condition from the geometric conditions given by the user. We have chosen this option, because it permits to avoid the stage of giving an initial condition plus the a priori conditions.

The first step consists in closing the arc to get a closed curve. We introduce a point *<sup>a</sup>* <sup>∈</sup> <sup>Ω</sup> and a vector *<sup>v</sup>* <sup>=</sup> (*v*1, *<sup>v</sup>*2) <sup>∈</sup> IR2. It is easy to construct two (unique) polynomials curves C<sup>1</sup> and C<sup>2</sup> ⊂ Ω that are parameterized using *c*<sup>1</sup> and *c*<sup>2</sup> ∈ *P*3([0, 1], Ω) such that:

$$\begin{array}{llll}c c(1) = c\_1(0), & c\_1(1) = c\_2(0) = a, & c\_2(1) = c(0), \\ -c'(1) = c'\_1(0), & c'\_1(1) = -c'\_2(0) = v, & c'\_2(1) = -c'(0). \end{array} \tag{7}$$

We set Γ<sup>0</sup> = C∪C<sup>1</sup> ∪ C2, and we consider Φ<sup>0</sup> as the signed distance to Γ<sup>0</sup> :

$$\Phi\_0(\mathbf{x}) = \begin{cases} \begin{array}{l} d(\mathbf{x}, \Gamma\_0) \\ -d(\mathbf{x}, \Gamma\_0) \end{array} & \text{if } \mathbf{x} \in Int\left(\Gamma\_0\right), \\\ -d(\mathbf{x}, \Gamma\_0) & \text{otherwise.} \end{cases} \tag{8}$$

To define the energy *EC*, we also introduce the mask *η*: Ω →IR:

$$\eta(\mathbf{x}) = \exp\left(-\frac{d(\mathbf{x}, \mathbf{C})^2}{\sigma}\right) : \forall \mathbf{x} \in \Omega,\tag{9}$$

where *<sup>σ</sup>* <sup>∈</sup>IR+, is a positive parameter which controls the width of the mask. Then, we can define the following energy functional:

$$E\_{\mathbb{C}}(\Phi) = \int\_{\Omega} \eta(\mathbf{x}) (\Phi(\mathbf{x}) - \Phi\_0(\mathbf{x}))^2 d\mathbf{x}.\tag{10}$$

The goal is to minimize the *L*<sup>2</sup> norm between Φ and Φ0: we impose the contour to be close to Γ0. The mask permits to impose such influence in (and only in) a local neighborhood of C.

#### *2.2. Minimization Problem and Evolution Equation*

The modeling of our problem consists in minimizing the energy:

$$E(\Phi) = E\_{\mathfrak{a}}(\Phi) + \frac{\mathfrak{a}}{2} E\_{\mathbb{C}}(\Phi) \tag{11}$$

with *α* > 0.

We now give the evolution problem. We minimize *E* with respect to Φ and deduce the associated Euler–Lagrange equation for Φ:

$$E'(\Phi) = \delta\_\varepsilon(\Phi) \text{div}\left(g(|\nabla I|) \frac{\nabla \Phi}{|\nabla \Phi|}\right) + a\eta (\Phi - \Phi\_0) \tag{12}$$

with the boundary conditions (homogeneous Neumann):

$$\frac{\delta\_{\mathfrak{c}}(\Phi)}{|\nabla\Phi|} \frac{\partial\Phi}{\partial\nu}\_{|\partial\Omega} = 0,\tag{13}$$

*ν* denoting the exterior normal to the boundary of Ω.

When a local minimum is reached, then the quantity *<sup>∂</sup>*<sup>Φ</sup> *<sup>∂</sup><sup>t</sup>* tends to 0, which means that the steady state is reached. A rescaling can be made so that the motion is applied to all level sets by replacing *δ* by |∇Φ| (It makes the flow independent of the scaling of Φ).

Let us also note that the speed of convergence of the model can be increased by adding the component *βg*(|∇*I*|)|∇Φ| in the evolution equation, *β* being a constant. This component can be seen as an area constraint. An analogy with the Balloon model developed by Cohen [26] can be made: this constant motion term prevents the curve from stopping on a non-significative local minimum and is also of importance when initializing the process with a curve inside the object to be detected.

Thus, the proposed parabolic problem with the associated boundary conditions *<sup>∂</sup>*<sup>Φ</sup> *∂ν* <sup>=</sup> <sup>0</sup> can be written:

$$\begin{cases} \Phi(x,0) &= \Phi\_0(x),\\ \frac{\partial \Phi}{\partial t} &= -|\nabla \Phi| \left[ \operatorname{div}(\mathcal{g}(|\nabla I|) \frac{\nabla \Phi}{|\nabla \Phi|}) \right] + a\eta (\Phi - \Phi\_0) \\ &+ \beta \mathcal{g}(|\nabla I|) |\nabla \Phi|, \end{cases} \tag{14}$$

#### *2.3. Existence, Uniqueness of the Solution*

This section is a theoretical part showing the existence and uniqueness of the solution of our segmentation model. We can show the existence and uniqueness of a viscosity solution to the evolution problem (14). We rewrite our problem as:

$$\begin{cases} \begin{array}{rcl} u\_t + \mathbb{P}(t, \mathbf{x}, u, Du, D^2 u) &=& 0 \quad \text{in} \quad ]0, T[ \times \Omega \\ B(\mathbf{x}, Du) &=& 0 \quad \text{in} \quad ]0, T[ \times \partial \Omega, \end{array} \end{cases} \tag{15}$$

where, for (*t*, *<sup>x</sup>*, *<sup>u</sup>*, *<sup>p</sup>*, *<sup>X</sup>*) <sup>∈</sup> [0, *<sup>T</sup>*] <sup>×</sup> <sup>Ω</sup>¯ <sup>×</sup> IR2 × S2:

$$\bar{F}(t, \mathbf{x}, u, p, X) = F(t, \mathbf{x}, p, X) + \mathbf{g}(\mathbf{x}, u), \quad B(\mathbf{x}, p) = p \cdot \vec{n}, \tag{16}$$

and:

$$\begin{array}{rcl} F(t, \mathbf{x}, p, \mathbf{X}) &=& -\text{trace}\left(\mathbf{g}(|\nabla I(\mathbf{x})|) \left(I - \frac{p \odot p}{|p|^2}\right) \mathbf{X}\right) \\ &- \langle \nabla (\mathbf{g}(|\nabla I(\mathbf{x})|)), p \rangle \\\\ g(\mathbf{x}, \boldsymbol{\mu}) &=& a \boldsymbol{\eta}(\mathbf{x}) (\boldsymbol{\mu} - \boldsymbol{\Phi}\_0(\mathbf{x})) .\end{array} \tag{17}$$

We assume that Ω is a bounded domain in IR<sup>n</sup> with a *C*<sup>1</sup> boundary. Under the following conditions, the Equation (15) has a unique viscosity solution.


$$(t, x, p, X) \in [0, T] \times \vec{\Omega} \times (\mathbb{R}^n - \{0\}) \times \mathbb{S}^n,$$

the function:

$$u \mapsto \mathbb{P}(t, \mathbf{x}, u, p, X) - \gamma u \tag{18}$$

is non-decreasing on IR.

3. For each *R* > 0, there exists a continuous function *wR* : [0, ∞[−→ [0, ∞[ satisfying *wR*(0) = 0 such that if *<sup>X</sup>*,*<sup>Y</sup>* <sup>∈</sup> *<sup>S</sup><sup>n</sup>* and *<sup>μ</sup>*1, *<sup>μ</sup>*<sup>2</sup> <sup>∈</sup> [0, <sup>∞</sup>[ satisfy:

$$
\mu \begin{pmatrix} X & 0 \\ 0 & Y \end{pmatrix} \le \mu\_1 \begin{pmatrix} I & -I \\ -I & I \end{pmatrix} + \mu\_2 \begin{pmatrix} I & 0 \\ 0 & I \end{pmatrix} \tag{19}
$$

then:

$$\begin{array}{l} \mathcal{F}(t, \mathbf{x}, \boldsymbol{\mu}, \boldsymbol{p}, \mathbf{X}) - \mathcal{F}(t, \boldsymbol{y}, \boldsymbol{\mu}, \boldsymbol{q}, -\boldsymbol{\chi}) \\ \geq -w \mathbb{R} \left( \mu\_1 (|\mathbf{x} - \boldsymbol{y}|^2 + \rho(\boldsymbol{p}, \boldsymbol{q})^2) \right) \\ + \mu\_2 + |\boldsymbol{p} - \boldsymbol{q}| + |\mathbf{x} - \boldsymbol{y}| (1 + \max(|\boldsymbol{p}|, |\boldsymbol{q}|)) ), \end{array} \tag{20}$$


All those conditions are clearly satisfied for *F*˜ = *F*, and *B* (see for instance [20]). We here have to extend the proof to *F*˜ = *F* + *g*.


The proof for conditions 4–6 are basic, and can be found for instance in [20]. Therefore, all the conditions are satisfied and the evolution problem (14) has a unique viscosity solution.

#### **3. Experimental Results**

For the discretization scheme, we use the additive operator splitting scheme (AOS) proposed by Weickert and Kuhne ([25], see also [20]). This corresponding algorithm requires a computational cost linear in *M* × *N* (size of the unknown vector Φ) at each step. This scheme is well-suited to our problem that differs only in the introduction of the function *η*(*x*).

We use the method proposed in this paper knowing that:


To illustrate the interest of our approach, we start with a simple test case, using a threeband image. Then, we run the method on blood vessel image. We obtain the following results. Without a priori condition, the segmentation process will of course give the entire white band as a result. This is the interest of adding such a priori condition to enforce the contour to arbitrary stop following the choice of the user.

We can see that the final contour correctly takes into account the curve given by the user imposing a constraint independent of the pixel values of the images, and correctly segment the white zone.

Of course, our approach allows us (to help) the segmentation process on "complicated" images, when ground truth informations are required like in Figure 1 or in many other cases, especially in geophysical or medical imaging. In Section 3.1, we illustrate the added value brought by our approach considering the initial guess choice improvement. Two examples are given to show that, compared to Gout et al. method [19], the approach given in this paper is much faster to get the same result since the initial guess is (strongly) optimized.

In Section 3.2, we propose to compare our approach with well-known algorithms (Chan-Vese [11], U-Net [27]).

In Sections 3.3 and 3.4, we give experimental results in medicine and geophysics.

**Figure 1.** Three-band image (**top left**), user-defined constraint (**top right**), initial contour automatically defined by our approach (**bottom left**) and final contour (**bottom right**).

#### *3.1. Impact of the Initial Guess on the Segmentation Process*

In this part, we focus on the added value brought by the automatic initial guess we propose, obtained from the geometrical condition given by the user. We compare our approach to the one introduced by Gout, Le Guyader and Vese [19], we take the same geometrical conditions, and then compare the number iteration to get the same final contour. The comparison is relevant because the segmentation algorithm is analogous in the two considered methods. In [19], note that we have to add a closed contour to the geometric conditions to start the segmentation process, while our proposed approach leads to an automatic initial guess (obtained from the geometric conditions given by the user). In the following two numerical examples, we take exactly the same interpolation conditions (3 points).

In Figures 2–4, we can see that our approach requires 5 times less iterations.

We repeat the same process in the following example (Figures 5–7). The initial contour for Gout et al. [19] algorithm is here chosen outside the region of interest. Of course, we take exactly the same interpolation conditions (3 points) in both cases.

In the example given on Figures 5–7, we can see that our approach requires 30 times less iterations.

In this subsection, we illustrate the gain of the automatic generated initial guess we propose from the geometric conditions (see Table 1).

**Figure 2.** We use Gout et al. [19] algorithm. The initial guess is given inside the region of interest, it is defined by <sup>Ψ</sup> = (*<sup>X</sup>* <sup>−</sup> <sup>39</sup>)<sup>2</sup> + (*<sup>Y</sup>* <sup>−</sup> <sup>47</sup>)<sup>2</sup> <sup>−</sup> <sup>18</sup>2, the time step is equal to 0.6, space step is equal to 0.1, distance is computed using the fast marching method (Sethian [14]) and the regularization term is equal to 0.8. **Top left**: Interpolation condition and initial closed contour. **Top right**: iteration 120. **Bottom left**: Iteration 240. **Bottom right**: Iteration 420.

**Figure 3.** Final result (iteration 420), with the interpolation conditions (3 points) using Gout et al. algorithm [19].

**Figure 4.** Using our proposed approach, the process does not require an initial guess, the interpolation conditions automatically initiate the process. **Left**: interpolation conditions. **Right**: final result after 80 iterations.

**Table 1.** Gain of the automatic generated initial guess.


**Figure 5.** We use Gout et al. [19] algorithm. The initial guess is given outside the region of interest, it is defined by <sup>Ψ</sup> <sup>=</sup> (*<sup>X</sup>* <sup>−</sup> <sup>30</sup>)<sup>2</sup> <sup>15</sup><sup>2</sup> <sup>+</sup> (*<sup>Y</sup>* <sup>−</sup> 30.5)<sup>2</sup> <sup>232</sup> <sup>−</sup> 1, the time step is equal to 0.3, space step is equal to 0.1, distance is computed using the fast marching method (Sethian [14]) and the regularization term is equal to 0.8. **Top left**: Interpolation condition and initial guess. **Top right**: iteration 80. **Bottom left**: Iteration 100. **Bottom right**: Iteration 160.

**Figure 6.** Final result (iteration 260), with the interpolation conditions (3 points) using Gout et al. algorithm [19].

**Figure 7. Left**: interpolation conditions. **Right**: final result after 8 iterations using our proposed approach.

To get the same final contour, our method clearly requires much less iterations. This result is analogous on all the tests we have done, especially when a closed contour is chosen outside the region of interest (in [19]), in this case, the total number of iterations is reduced by a factor of 5 to 60. If the initial contour in [19] is chosen inside the region of interest, the total number of iterations is reduced by a factor of 4 to 20. The only cases where the number of iterations is equivalent between the two different methods are when we choose as initial guess in [19] a closed curve very close to the final result and interpolating at least 2 given interpolation conditions, this is of course a major constraint.

#### *3.2. Quantitative Performance*

We use the BraTS Dataset [28], which was collected and shared by the MICCAI Brain Tumor Segmentation Challenge (see [28] for more details) several years ago. We use 270 MR scans, each with four MIR sequences: T1-weighted (native image, sagittal or axial 2D acquisitions with 1–6 mm slice thickness), T1-weighted with contrast-enhanced (Gadolinium) image, T2-weighted image (axial 2D acquisition, with 2–6 mm slice thickness), and FLAIR

(T2-weighted FLAIR image, axial, coronal, or sagittal 2D acquisitions, 2–6 mm slice thickness) (see [28] for more details). The training data have the size 240 × 240 × 155 pixels, and have the manual segmentation labels for many types of brain tumors. We trained the network to generate segmentation masks corresponding to the center slice of the input. We used the well-known Adam optimization method [29]. We trained the deep network using 238 training data, we set the initial learning rate as 10−4, and multiplied by 0.5 after every 20 epochs. Using a single GPU, we trained the models with batch size 4 for 40 epochs.

In order to evaluate the performance of our method, we give comparisons between several algorithms. We here recall several metrics:

– The Jaccard index [30], or Intersection over Union (IoU), is a commonly used metric in segmentation. It is defined as the area of intersection between the predicted segmentation map and the ground truth, divided by the area of union between the predicted segmentation map and the ground truth:

$$\text{IoU} = f(G, P) = \frac{|G \cap P|}{|G \cup P|} \, \tag{21}$$

where *G* is the ground truth and *P* denotes the predicted segmentation maps. It ranges between 0 and 1. Mean-IoU (mIoU) is defined as the average IoU over all classes. It is widely used in reporting the performance of segmentation algorithms.

– The Dice coefficient (Dice) is a popular metric for image segmentation, especially in medical imaging. This coefficient can be defined as twice the overlap area of predicted and ground-truth maps, divided by the total number of pixels in both images:

$$\text{Dice} = \frac{2|G \cap P|}{|G| + |P|} \tag{22}$$

When applied to binary segmentation maps, and referring to the foreground as a positive class, the Dice coefficient is essentially identical to the F1 score:

$$\text{Dice} = \frac{2TP}{2TP + FP + FN} = \text{F1}\_{\prime} \tag{2.3}$$

(The F1 score, which is defined as the harmonic mean of precision (Prec) and recall (Rec): F1<sup>=</sup> 2Prec Rec Prec + Rec , where Prec <sup>=</sup> *TP TP* <sup>+</sup> *FP* and Rec <sup>=</sup> *TP TP* <sup>+</sup> *FN* .), where TP refers to the true positive fraction, FP refers to the false positive fraction, and FN refers to the false negative fraction.

– The Hausdorff distance (Hd) evaluates the quality of the segmentation boundaries by computing the maximum distance between the prediction and its ground truth:

$$\text{Hd} = \max \left\{ \sup\_{p \in \partial P} \inf\_{z \in \partial G} ||p - z||\_{2^\prime} \sup\_{z \in \partial G} \inf\_{p \in \partial P} ||p - z||\_2 \right\} \tag{24}$$

where *∂P* and *∂G* are the surface point sets of *P* (segmentation prediction) and *G* (ground truth).

Of course, a large Dice coefficient or a small HD means an accurate segmentation result. Figure 8 shows several numerical example (here, we ramdomly sampled 33% of BRATS labeled dataset). In Figure 8, we give the dice score obtained for U-Net and our proposed algorithm. We can see that our approach clearly gives better results.

We also compare with the U-net algorithm ([27], see also [31]), tables show the evaluation results, please note that we have suppressed the best and worst 10% of the results, and that for our proposed approach, we have considered only 2 to 4 points as geometric conditions (the initial guess -curve- is automatically constructed from these points), all that leads to make the global results quite homogeneous. We also compare with the Chan-Vese algorithm [11]: in this case, the initial condition is a closed contour located inside the region of interest.

If we use 100% of the BraTS dataset, we obtain the Table 2.

**Figure 8.** Brain Tumor Segmentation (using BRATS dataset [28]). Several results from our method using labeled data on the BRATS dataset [28]. The tumor is underlined in yellow. First row: center slices of input. Second row: initial conditions for our algorithm. Third row: ground-truth labels. Fourth row: results from the supervised U-Nets learning method introduced by Ronnenberger et al. [27]. Fifth row: results from our proposed method. The scores in the bottom of each results denote Dice score.


**Table 2.** Accuracies of segmentations models on a third of labeled data from BraTS Dataset [28]. We use a Nvidia GeForce RTX 2080 Ti 11G Turbo Edition (Boost frequence: 1545 MHz, GPU memory: 11 GB).

We see that the GPU time is clearly better with our algorithm (and the Chan-Vese one too). The quality is globally equivalent: a little better with our algorithm when not considering the entire database, and equivalent if not. This is an important point because it validates our algorithm when we compare with MR scans segmented with deep learning algorithms with few labeled data.

#### *3.3. Applications to Medical Imaging*

Here, we apply our method to 2D medical images, the goal being to outline the cross-sectional area of a great thoracic vessel, namely the main pulmonary artery, in order to non-invasively assess pulmonary arterial hypertension. Figure 9 refers to a slice set perpendicularly to the MPA axis (arrow) of a 78 years old patient, who was suspected of pulmonary arterial hypertension, suffering from breathlessness. The presented magnitude image was acquired from a velocity encoded sequence (Venc = 1.5 m/s), during a 5 min time duration, with a 30 ms temporal resolution (1T Magnetom Expert imager; Siemens; Germany). The image quality is moderate, and in particular the vessel wall is irregular around the arrow, that is, not anatomically possible, and therefore an interpolation in the outlining by the physician is needed.

**Figure 9.** Example of an image sequence (courtesy of CHU Bordeaux, France). The arrow shows the pulmonary artery to segment.

Since the lack of annotated data is one of the most common limitations encountered in many deep learning approaches like the ones we will study in this subsection, deep leaning approaches are not suitable (at least for now!) for our considered applications.

An example of the segmentation of a pulmonary artery with a curve to interpolate is given in Figure 10. A part of this artery is difficult to segment since its bottom left is always located in a blurred zone.

**Figure 10.** Blood vessel image, user-defined constraint (**top right**), initial contour (**bottom left**) and final (**bottom right**) contours. Let us note that we have successfully segmented the complete image sequence.

Let us note that the MPA CSA variations versus time, throughout a cardiac cycle, obtained from the present automatic segmentation method are in agreement with previous results obtained from a manual outlining (Hôpital of Bordeaux, France). Moreover, the present method enables us to get a reliable assessment of the systolic-diastolic difference in the vessel CSA, which should lead to a further improvement in the accuracy of a non-invasive blood pressure estimation (using well-known methods, see [32,33]). In several medical applications, let us note that a registration process must be added to the segmentation process (see [34–37] for more details).

#### *3.4. Applications to Geophysical Imaging*

We now apply our algorithm to the segmentation of a non-continuous layer on a geophysical dataset (Figure 11 (right)) extracted from a 3D dataset (Figure 11 (left)). This is a complicated image to segment since there is a vertical fault (see Figure 12). As a priori condition, we here consider a set of point to interpolate (instead of a curve) in order to help the segmentation process (Figure 12). The algorithm starts with a polygonal line interpolating the set of points (to interpolate) as an initial condition, and then segment the image until the edge(s) of the image. Here, the final contour is a segment and not a closed curve, this is of great interest when one wants to track faults or cracks in segmentation processes.

The obtained result (Figure 13) has been validated by geophysicians, but the next step requires an algorithm for 3D applications in order to segment 3D geophysical datasets. We are also in the process of creating a free labeled geophysical dataset, in order to test deep learning techniques on these images.

**Figure 11. Left**: 3D seismic dataset (courtesy of TotalEnergies, Pau, France). Two strong and continuous reflectors (Horizon A and Horizon B) appear. **Right**: 2D image (extracted from the 3D seismic dataset) with vertical faults and layers. Considering the complexity of such dataset, it is impossible to segment it without a priori given conditions. The segmentation process consists in finding layers and/or faults.

**Figure 12. Left**: The red zone underlines the vertical fault (with dotted segments). **Right**: As a geometric conditions, we consider a set of points given by the user.

**Figure 13.** Final result: we can see that the searched layer is correctly segmented. On this specific dataset, the final contour has to reach the edge of the image: the contour goes beyond the last point to interpolate (right side) to join the edge of the image.

#### **4. Conclusions**

We have proposed a method to segment medical or geophysical images with a priori conditions. The efficiency of the method has been shown on 2D images, knowing that such images are difficult to segment with usual approaches since the searched contour is not continuous (geophysical dataset) or when the zone to segment is blurred (medical images).

Compared to Gout et al. [19], the added value brought by our proposed method (automatic initial guess) is materialized by a total number of iterations reduced by a factor of 4 to 60.

The GPU time is also a key point to underline about our algorithm, much lower than deep learning based approaches. About the sensitivity of the model, especially for geophysical data, it is not possible to segment correctly without a precise choice of the interpolation conditions: the choice of such condition strongly impacts the final result, this point is less obvious in medical applications where the choice of the initial condition improves the result only on very noisy areas and/or with a blurred edge, otherwise, the result (whatever the initial condition) is very close to classic algorithms (Chan Vese algorithm for example). The implementation of a 3D code is a work in progress, the main problem for a 3D segmentation process consists in choosing the initial (geometric) condition (a surface, a curve and/or points), which is a complicated task on complex 3D images (mainly because there are 3D data visualization difficulties "to see inside" the 3D dataset, especially in geophysics to give the initial condition—see Figure 14)—the mathematical results remain valid in 3D.

**Figure 14. Top**: Another 3D view of the 3D geological dataset of Figure 5. **Bottom**: From the 3D dataset, we use a set of points (well data given by geologists), and each surface/layer is then constructed with a *Dm*-spline operator (see [22] for more details on *Dm*-spline surface approximation). We here automatically construct 7 surfaces that will be considered as an initial condition. Next step consists in the 3D segmentation process from these initial conditions—this last part is a work in progress.

**Author Contributions:** Conceptualization, G.K., C.G., T.C.-F. and S.K.; methodology, G.K., C.G. and T.C.-F.; software, G.K. and T.C.-F.; validation, G.K., C.G. and S.K.; formal analysis, G.K.; investigation, G.K., C.G., T.C.-F. and S.K.; writing—original draft preparation, G.K.; writing—review and editing, G.K., C.G., T.C.-F. and S.K.; funding acquisition, C.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research was funded by the European Union with the European regional development fund (ERDF, 18P03390/18E01750/18P02733/21E05300) and by the Région Normandie Council via the M2SINUM and DEFHY3GEO projects, Agence Nationale de la Recherche via MEDISEG project (ANR-21-CE23-0013), and Maison Normande des Sciences du Numérique (MNSN, under DSG2-2021 project).

**Acknowledgments:** The authors thank the CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie, France) for providing computational resources.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

