# Differentiable Soft Morphological Filters for Medical Image Segmentation

Lisa Guzzi, Maria A. Zuluaga, Fabien Lareyre, Gilles Di Lorenzo, Sébastien Goffart, Andrea Chierici, Juliette Raffort, Hervé Delingette

► To cite this version:

Lisa Guzzi, Maria A. Zuluaga, Fabien Lareyre, Gilles Di Lorenzo, Sébastien Goffart, et al. *Differentiable Soft Morphological Filters for Medical Image Segmentation*. MICCAI 2024 - Medical Image Computing and Computer Assisted Intervention, Oct 2024, Marrakesh, Morocco. hal-04592007

HAL Id: hal-04592007

https://hal.science/hal-04592007

Submitted on 29 May 2024

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.
# Differentiable Soft Morphological Filters for Medical Image Segmentation

Lisa Guzzi, Maria A. Zuluaga, Fabien Lareyre, Gilles Di Lorenzo, Sébastien Goffart, Andrea Chierici, Juliette Raffort, Hervé Delingette

► To cite this version:

Lisa Guzzi, Maria A. Zuluaga, Fabien Lareyre, Gilles Di Lorenzo, Sébastien Goffart, et al. *Differentiable Soft Morphological Filters for Medical Image Segmentation*. MICCAI 2024 - Medical Image Computing and Computer Assisted Intervention, Oct 2024, Marrakesh, Morocco. hal-04592007v2

HAL Id: hal-04592007

https://hal.science/hal-04592007v2

Submitted on 10 Jul 2024

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
# Differentiable Soft Morphological Filters for Medical Image Segmentation

Lisa Guzzi<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, Maria A. Zuluaga<fnref n="2" />, Fabien Lareyre<fnref n="4" />,<fnref n="5" />, Gilles Di Lorenzo<fnref n="4" />, Sébastien Goffart<fnref n="1" />,<fnref n="3" />, Andrea Chierici<fnref n="3" />, Juliette Raffort<fnref n="3" />,<fnref n="5" />, and Hervé Delingette<fnref n="1" />

**Abstract.** Morphological operations such as erosion, dilation, and skeletonization offer valuable tools for processing and analyzing segmentation masks. Several studies have investigated the integration of differentiable morphological operations within deep segmentation neural networks, particularly for the computation of loss functions. However, those methods have shown limitations in terms of reliability, versatility or applicability to different types of operations and image dimensions. In this paper, we present a novel framework that provides differentiable morphological filters on probabilistic maps. Given any morphological filter defined on 2D or 3D binary images, our approach generates a soft version of this filter by translating Boolean expressions into multilinear polynomials. Moreover, using proxy polynomials, these soft filters have the same computational complexity as the original binary filter. We demonstrate on diverse biomedical datasets that our method can be easily integrated into neural networks either as a loss function or as the final morphological layer in a segmentation network. In particular, we show that the proposed filters for mask erosion, dilation or skeletonization lead to competitive solutions compared to the state-of-the-art.

**Keywords:** Image Segmentation · Morphological Operations · Deep Learning

## 1 Introduction

Morphological operations are often required for analyzing images based on object shape and structure [19, 21]. These operations include morphological erosion and dilation for edge detection [6, 23, 27], morphological opening for noise reduction [2, 4], and skeletonization for analyzing vessel trajectory and topology [9].

Recent works have explored the integration of morphological operations either as layers of a neural network or as loss functions. In all cases, morphological

<footnote n="1">¹ Université Côte d'Azur, Inria, Epione Team, Sophia Antipolis, France</footnote>

<footnote n="2">² Data Science Department, EURECOM, Sophia Antipolis, France</footnote>

<footnote n="3">³ University Hospital of Nice, Nice, France</footnote>

<footnote n="4">⁴ Department of Vascular Surgery, Hospital of Antibes Juan-les-Pins, Antibes, France</footnote>

<footnote n="5">⁵ Université Côte d'Azur, CNRS, UMR7370, LP2M, Nice, France
lisa.guzzi@inria.fr</footnote>
# Differentiable Soft Morphological Filters for Medical Image Segmentation

Lisa Guzzi<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, Maria A. Zuluaga<fnref n="2" />, Fabien Lareyre<fnref n="4" />,<fnref n="5" />, Gilles Di Lorenzo<fnref n="4" />, Sébastien Goffart<fnref n="1" />,<fnref n="3" />, Andrea Chierici<fnref n="3" />, Juliette Raffort<fnref n="3" />,<fnref n="5" />, and Hervé Delingette<fnref n="1" />

lisa.guzzi@inria.fr

**Abstract.** Morphological operations such as erosion, dilation, and skeletonization offer valuable tools for processing and analyzing segmentation masks. Several studies have investigated the integration of differentiable morphological operations within deep segmentation neural networks, particularly for the computation of loss functions. However, those methods have shown limitations in terms of reliability, versatility or applicability to different types of operations and image dimensions. In this paper, we present a novel framework that provides differentiable morphological filters on probabilistic maps. Given any morphological filter defined on 2D or 3D binary images, our approach generates a soft version of this filter by translating Boolean expressions into multilinear polynomials. Moreover, using proxy polynomials, these soft filters have the same computational complexity as the original binary filter. We demonstrate on diverse biomedical datasets that our method can be easily integrated into neural networks either as a loss function or as the final morphological layer in a segmentation network. In particular, we show that the proposed filters for mask erosion, dilation or skeletonization lead to competitive solutions compared to the state-of-the-art.

**Keywords:** Image Segmentation · Morphological Operations · Deep Learning

# 1 Introduction

Morphological operations are often required for analyzing images based on object shape and structure [19, 21]. These operations include morphological erosion and dilation for edge detection [6, 23, 27], morphological opening for noise reduction [2, 4], and skeletonization for analyzing vessel trajectory and topology [9].

Recent works have explored the integration of morphological operations either as layers of a neural network or as loss functions. In all cases, morphological

<footnote n="1">¹ Université Côte d'Azur, Inria, Epione Team, Sophia Antipolis, France</footnote>

<footnote n="2">² Data Science Department, EURECOM, Sophia Antipolis, France</footnote>

<footnote n="3">³ University Hospital of Nice, Nice, France</footnote>

<footnote n="4">⁴ Department of Vascular Surgery, Hospital of Antibes Juan-les-Pins, Antibes, France</footnote>

<footnote n="5">⁵ Université Côte d'Azur, CNRS, UMR7370, LP2M, Nice, France</footnote>
operations must be differentiable for gradient-based optimization during training [8]. Some works have investigated the use of supervised Convolutional Neural Networks (CNN) to emulate specific morphological operations. While differentiable, these methods [1, 12, 17] require substantial training time, datasets and fine-tuning, and they only apply to a given operator. As morphological operations were originally designed for binary images and subsequently extended to grayscale images with min-max filters [21], other works have considered soft min and maxpooling layers (based on the counter-harmonic mean) to replace erosion and dilation in a differentiable manner [5]. Some studies aim at replacing some convolutional layers with non-linear feature-extracting morphological layers leading to the so-called *Deep Morphological Neural Networks* [3, 11, 13, 20]. Other authors combine traditional segmentation CNNs with loss functions requiring morphological operators. For instance, Jurdi et al. [5] compute the perimeter of a segmented structure using min and max pooling. The clDice loss [22] introduces a similarity measure based on the interpolation of the segmentation masks and their associated skeletons to preserve the topology of tubular structures, while computing the soft skeleton of a segmentation by iterating between min and maxpooling layers. However, this method assumes a fixed structuring element and may not preserve an object's topology by creating discontinuities in the skeleton, leading to topological errors. Menten et al. [10] propose a skeletonization operation compatible with gradient-based optimization by using a binary morphological operator on binary images sampled from probability maps using the reparameterization trick. However, this approach requires specifying several sensitive hyperparameters and is restricted to the extraction of skeletons.

In this paper, we propose a novel extension of binary morphological operations on probability maps that can be seamlessly integrated into neural networks either as a loss function or as a final morphological layer. A first methodological contribution consists of the translation of any morphological operation based on Boolean expressions into a single multilinear or proxy polynomial. The resulting soft morphological filters are differentiable, have no hyperparameters to be tuned, and can be generated from any binary morphological filter.

In a second contribution, we demonstrate their impact on two segmentation applications, while comparing them with state-of-the-art morphological filters. First, we define a soft 2D skeletonization and integrate it into the clDice loss [22] for retinal vessel segmentation. Second, we integrate erosion/dilation filters as the final layer of a 2D and 3D U-Net to integrate a post-processing step within the network for improving the detection of calcification plaques in computed tomography angiography (CTA) of the lower limbs. To foster reproducibility, our code is available at https://github.com/lisaGUZZI/Soft-morph

# 2 Method

## 2.1 Probabilistic Definition of Soft Morphological filters

We consider a generic binary segmentation problem of an image $X$ consisting of $N$ voxels $\{X_n\}, n = 1...N$ for any image $X$. We want to formalize the definition
operations must be differentiable for gradient-based optimization during training [8]. Some works have investigated the use of supervised Convolutional Neural Networks (CNN) to emulate specific morphological operations. While differentiable, these methods [1, 12, 17] require substantial training time, datasets and fine-tuning, and they only apply to a given operator. As morphological operations were originally designed for binary images and subsequently extended to grayscale images with min-max filters [21], other works have considered soft min and maxpooling layers (based on the counter-harmonic mean) to replace erosion and dilation in a differentiable manner [5]. Some studies aim at replacing some convolutional layers with non-linear feature-extracting morphological layers leading to the so-called *Deep Morphological Neural Networks* [3, 11, 13, 20]. Other authors combine traditional segmentation CNNs with loss functions requiring morphological operators. For instance, Jurdi et al. [5] compute the perimeter of a segmented structure using min and max pooling. The clDice loss [22] introduces a similarity measure based on the interpolation of the segmentation masks and their associated skeletons to preserve the topology of tubular structures, while computing the soft skeleton of a segmentation by iterating between min and maxpooling layers. However, this method assumes a fixed structuring element and may not preserve an object's topology by creating discontinuities in the skeleton, leading to topological errors. Menten et al. [10] propose a skeletonization operation compatible with gradient-based optimization by using a binary morphological operator on binary images sampled from probability maps using the reparameterization trick. However, this approach requires specifying several sensitive hyperparameters and is restricted to the extraction of skeletons.

In this paper, we propose a novel extension of binary morphological operations on probability maps that can be seamlessly integrated into neural networks either as a loss function or as a final morphological layer. A first methodological contribution consists of the translation of any morphological operation based on Boolean expressions into a single multilinear or proxy polynomial. The resulting soft morphological filters are differentiable, have no hyperparameters to be tuned, and can be generated from any binary morphological filter.

In a second contribution, we demonstrate their impact on two segmentation applications, while comparing them with state-of-the-art morphological filters. First, we define a soft 2D skeletonization and integrate it into the clDice loss [22] for retinal vessel segmentation. Second, we integrate erosion/dilation filters as the final layer of a 2D and 3D U-Net to integrate a post-processing step within the network for improving the detection of calcification plaques in computed tomography angiography (CTA) of the lower limbs. To foster reproducibility, our code is available at https://github.com/lisaGUZZI/Soft-morph

# 2 Method

## 2.1 Probabilistic Definition of Soft Morphological filters

We consider a generic binary segmentation problem of an image $X$ consisting of $N$ voxels $\{X_n\}, n = 1...N$ for any image $X$. We want to formalize the definition
<figure><img src="image_1.png" /><figcaption>Fig. 1: (a) Definition of the 4 + 1 neighborhood variables;(b) Probabilistic map created from a binary retinal blood vessel;(c) Match between the output of 2D thinning operators based on multilinear vs proxy polynomials.</figcaption></figure>

of a soft morphological filter $\mathcal{F}^*(\cdot)$ applied on a probabilistic segmentation image $\mathcal{Y} = \{y_n\} \in [0, 1]^N$ and generate a new probabilistic image $\mathcal{Z} = \{z_n\} \in [0, 1]^N = \mathcal{F}^*(\mathcal{Y})$. This filter should generalize a given binary filter $Z = F(Y) \in \{0, 1\}^N$, $Y \in \{0, 1\}^N$ such that both give the same result when the input probabilistic image is binary $Z = F(Y) = \mathcal{F}^*(Y)$.

More precisely, the input probabilities $y_n$ correspond to the posterior probability $y_n = p(Y_n = 1|X) \in [0, 1]$ of the binary variables $Y_n \in \{0, 1\}$, resulting from a segmentation algorithm, typically the output of a neural network. The binary morphological operator $F(\cdot)$ is applied on the (unknown) binary image $Y$ resulting in a new binary segmentation $Z$. This operator $F(\cdot)$ takes as input $k$ binary variables $Y_{\mathcal{N}(1,n)}, \dots, Y_{\mathcal{N}(k,n)}$ and outputs a binary variable $Z_n = F(Y_{\mathcal{N}(1,n)}, \dots, Y_{\mathcal{N}(k,n)}) \in \{0, 1\}$. The neighborhood function $\mathcal{N}(i, n)$ provides the index of the $i$th neighbor of voxel $n$ and in mathematical morphology, it provides the domain where the structuring element is defined. For instance, considering the typical neighborhoods of 2D or 3D images leads to the choice of $k = 4 + 1$ or $k = 8 + 1$ in 2D and $k = 6 + 1, 18 + 1, 26 + 1$ in 3D.

The binary operator $F(\cdot)$ is a Boolean function that has $2^k$ different input values and outputs a binary variable. The exhaustive list of those values $F(a) \in \{0, 1\}$, for $a \in \{0, 1\}^k$ is called the *truth table* of $F(\cdot)$ and can be established for small values of $k$. Besides, it can be shown that any Boolean function can be written as a propositional formula involving the $k$ binary variables with the logical operator AND ($\land$), OR ($\lor$) and NOT ($\neg$). For instance, based on the notation defined in Fig. 1a, the dilation operator acting on a $4+1$-neighborhood of a 2D image can be written as $f_{\text{Dil}} = Y0 \lor Y1 \lor Y2 \lor Y3 \lor Y4$ whereas the erosion operator is $f_{\text{Ero}} = Y0 \land Y1 \land Y2 \land Y3 \land Y4$.

The morphological operator $F(\cdot)$ is deterministic but its input is a binary image $Y$ only known through its posterior probability $\mathcal{Y} = p(Y|X)$. Therefore we aim at estimating the posterior $z_n = p(Z_n = 1|X)$ of the final segmentation $Z$ knowing that it results from the application of the morphological operation $Z = F(Y)$. Formally, this posterior can be obtained through the law of total
<figure><img src="image_2.png" /><figcaption>Fig. 1: (a) Definition of the 4 + 1 neighborhood variables;(b) Probabilistic map created from a binary retinal blood vessel;(c) Match between the output of 2D thinning operators based on multilinear vs proxy polynomials.</figcaption></figure>

of a soft morphological filter $\mathcal{F}^*(\cdot)$ applied on a probabilistic segmentation image $\mathcal{Y} = \{y_n\} \in [0, 1]^N$ and generate a new probabilistic image $\mathcal{Z} = \{z_n\} \in [0, 1]^N = \mathcal{F}^*(\mathcal{Y})$. This filter should generalize a given binary filter $Z = F(Y) \in \{0, 1\}^N$, $Y \in \{0, 1\}^N$ such that both give the same result when the input probabilistic image is binary $Z = F(Y) = \mathcal{F}^*(Y)$.

More precisely, the input probabilities $y_n$ correspond to the posterior probability $y_n = p(Y_n = 1|X) \in [0, 1]$ of the binary variables $Y_n \in \{0, 1\}$, resulting from a segmentation algorithm, typically the output of a neural network. The binary morphological operator $F(\cdot)$ is applied on the (unknown) binary image $Y$ resulting in a new binary segmentation $Z$. This operator $F(\cdot)$ takes as input $k$ binary variables $Y_{\mathcal{N}(1,n)}, \dots, Y_{\mathcal{N}(k,n)}$ and outputs a binary variable $Z_n = F(Y_{\mathcal{N}(1,n)}, \dots, Y_{\mathcal{N}(k,n)}) \in \{0, 1\}$. The neighborhood function $\mathcal{N}(i, n)$ provides the index of the $i$th neighbor of voxel $n$ and in mathematical morphology, it provides the domain where the structuring element is defined. For instance, considering the typical neighborhoods of 2D or 3D images leads to the choice of $k = 4 + 1$ or $k = 8 + 1$ in 2D and $k = 6 + 1, 18 + 1, 26 + 1$ in 3D.

The binary operator $F(\cdot)$ is a Boolean function that has $2^k$ different input values and outputs a binary variable. The exhaustive list of those values $F(a) \in \{0, 1\}$, for $a \in \{0, 1\}^k$ is called the *truth table* of $F(\cdot)$ and can be established for small values of $k$. Besides, it can be shown that any Boolean function can be written as a propositional formula involving the $k$ binary variables with the logical operator AND ($\land$), OR ($\lor$) and NOT ($\neg$). For instance, based on the notation defined in Fig. 1a, the dilation operator acting on a $4+1$-neighborhood of a 2D image can be written as $f_{\text{Dil}} = Y0 \lor Y1 \lor Y2 \lor Y3 \lor Y4$ whereas the erosion operator is $f_{\text{Ero}} = Y0 \land Y1 \land Y2 \land Y3 \land Y4$.

The morphological operator $F(\cdot)$ is deterministic but its input is a binary image $Y$ only known through its posterior probability $\mathcal{Y} = p(Y|X)$. Therefore we aim at estimating the posterior $z_n = p(Z_n = 1|X)$ of the final segmentation $Z$ knowing that it results from the application of the morphological operation $Z = F(Y)$. Formally, this posterior can be obtained through the law of total
probability as the expectation of filtered binary segmentation $F(a)$:

$$ 
\begin{aligned}
z_n = p(Z_n = 1|X) &= \sum_{Y_1=0}^{1} \dots \sum_{Y_N=0}^{1} p(Z_n = 1|Y) p(Y|X) \\
&= \sum_{a \in \{0,1\}^N} F(a) p(a|X) = \mathbb{E}_{a \sim p(Y|X)} F(a) \quad (1)
\end{aligned}
$$

Eq.(1) defines implicitly the relationship $z_n = \mathcal{F}^*(\mathcal{Y})$ of the soft morphological filter and it is furthermore required to make the filter differentiable, i.e. to estimate the derivatives $\frac{\partial \mathcal{F}^*(\mathcal{Y})}{\partial y_m}$

## 2.2 Soft Morphological Filters as multilinear Polynomials

To provide a closed-form expression of a soft filter defined from Eq.(1), we propose to adopt a polynomial representation of the Boolean function $F(a)$, $a \in \{0, 1\}^k$. Indeed we can associate with any Boolean function $F(a)$ a multilinear polynomial $\mathcal{F}^*(x)$, $x = (x_1, \dots, x_k)^T \in \mathbb{R}^k$ defined as:

$$ \mathcal{F}^*(x) = \sum_{a \in \{0,1\}^k} F(a) \prod_{i=1}^{k} x_i^{1-a_i} (1-x_i)^{a_i} \quad (2) $$

It is easy to see that $F(a) = \mathcal{F}^*(a)$, i.e. that the polynomial $\mathcal{F}^*(x)$ coincides by construction with the Boolean function on the hypercube $\{0, 1\}^k$. Each monomial $F(a) \prod_{i:a_i=0} x_i \prod_{i:a_i=1} (1-x_i)$ is of degree $k$ and linear with each variable $x_j$, making the polynomial multilinear. This property leads to the following result:

**Theorem** The expectation of a Boolean function $F(a) \in \{0, 1\}$, $a \in \{0, 1\}^k$ over a set of $k$ independent variables with $a \sim \text{Bernoulli}(p)$, $p \in [0, 1]^k$ is $\mathcal{F}^*(p)$

*Proof.* It is easy to show that $\mathbb{E}_{a_i \sim \text{Bernoulli}(p_i)}(\alpha + \beta a_i) = \alpha + \beta p_i$ using the linearity of expectation. Thus, we have

$$ \mathbb{E}_{a_1} (F(a) \prod_{j=1}^{k} a_j^{1-a_j} (1-a_j)^{a_j}) = F(a) p_1^{1-a_1} (1-p_1)^{a_1} \prod_{i=2}^{k} a_i^{1-a_i} (1-a_i)^{a_i}. $$

By taking the expectation over each variable $a_i$, $i > 1$, we get $\mathbb{E}_a(F(a)) = \sum_{a \in \{0,1\}^k} F(a) \prod_{i=1}^{k} p_i^{1-a_i} (1-p_i)^{a_i} = \mathcal{F}^*(p)$. $\square$

Therefore, assuming that the marginal posteriors $y_n = p(Y_n|X)$ are independently distributed (which is the case when dealing with the output of segmentation neural networks or mean field approximations), we define the soft morphological filter associated with the binary filter $F()$ as the polynomial value $\mathcal{F}^*(y_{\mathcal{N}(1,n)}, \dots, y_{\mathcal{N}(k,n)}) \in [0, 1]$. Besides, this result can be simply interpreted as the multilinear interpolation of the Boolean function over the hypercube $\{0, 1\}^k$.

Thus, the soft erosion operator is simply $\mathcal{Z}_{\text{Ero}} = y_0 y_1 y_2 y_3 y_4$ whereas the dilation operator is $\mathcal{Z}_{\text{Dil}} = \sum_{a \in \{0,1\}^5, a \neq 0} \prod_{i=1}^{5} y_{i-1}^{1-a_i} (1-y_{i-1})^{a_i}$
probability as the expectation of filtered binary segmentation $F(a)$ :

$$ 
\begin{aligned}
z_n = p(Z_n = 1|X) &= \sum_{Y_1=0}^{1} \dots \sum_{Y_N=0}^{1} p(Z_n = 1|Y) p(Y|X) \\
&= \sum_{a \in \{0,1\}^N} F(a) p(a|X) = \mathbb{E}_{a \sim p(Y|X)} F(a)
\end{aligned} 
\quad (1) $$

Eq.(1) defines implicitly the relationship $z_n = \mathcal{F}^*(\mathcal{Y})$ of the soft morphological filter and it is furthermore required to make the filter differentiable, i.e. to estimate the derivatives $\frac{\partial \mathcal{F}^*(\mathcal{Y})}{\partial y_m}$

## 2.2 Soft Morphological Filters as multilinear Polynomials

To provide a closed-form expression of a soft filter defined from Eq.(1), we propose to adopt a polynomial representation of the Boolean function $F(a)$, $a \in \{0, 1\}^k$. Indeed we can associate with any Boolean function $F(a)$ a multilinear polynomial $\mathcal{F}^*(x)$, $x = (x_1, \dots, x_k)^T \in \mathbb{R}^k$ defined as:

$$ \mathcal{F}^*(x) = \sum_{a \in \{0,1\}^k} F(a) \prod_{i=1}^{k} x_i^{a_i} (1-x_i)^{1-a_i} \quad (2) $$

It is easy to see that $F(a) = \mathcal{F}^*(a)$, i.e. that the polynomial $\mathcal{F}^*(x)$ coincides by construction with the Boolean function on the hypercube $\{0, 1\}^k$. Each monomial $F(a) \prod_{i:a_i=1} x_i \prod_{i:a_i=0} (1-x_i)$ is of degree $k$ and linear with each variable $x_j$, making the polynomial multilinear. This property leads to the following result:

**Theorem** The expectation of a Boolean function $F(a) \in \{0, 1\}$, $a \in \{0, 1\}^k$ over a set of $k$ independent variables with $a \sim \text{Bernoulli}(p)$, $p \in [0, 1]^k$ is $\mathcal{F}^*(p)$

*Proof.* It is easy to show that $\mathbb{E}_{a_i \sim \text{Bernoulli}(p_i)}(\alpha + \beta a_i) = \alpha + \beta p_i$ using the linearity of expectation. Thus, we have

$$ \mathbb{E}_{a_1} (F(a) \prod_{j=1}^{k} a_j^{a_j} (1-a_j)^{1-a_j}) = F(a) p_1^{a_1} (1-p_1)^{1-a_1} \prod_{i=2}^{k} a_i^{a_i} (1-a_i)^{1-a_i}. $$

By the taking the expectation over each variable $a_i$, $i > 1$, we get $\mathbb{E}_a(F(a)) = \sum_{a \in \{0,1\}^k} F(a) \prod_{i=1}^{k} p_i^{a_i} (1-p_i)^{1-a_i} = \mathcal{F}^*(p)$. $\square$

Therefore, assuming that the marginal posteriors $y_n = p(Y_n|X)$ are independently distributed (which is the case when dealing with the output of segmentation neural networks or mean field approximations), we define the soft morphological filter associated with the binary filter $F()$ as the polynomial value $\mathcal{F}^*(y_{\mathcal{N}(1,n)}, \dots, y_{\mathcal{N}(k,n)}) \in [0, 1]$. Besides, this result can be simply interpreted as the multilinear interpolation of the Boolean function over the hypercube $\{0, 1\}^k$.

Thus, the soft erosion operator is simply $\mathcal{Z}_{\text{Ero}} = y_0 y_1 y_2 y_3 y_4$ whereas the dilation operator is $\mathcal{Z}_{\text{Dil}} = \sum_{a \in \{0,1\}^5, a \neq 0} \prod_{i=1}^{5} y_{i-1}^{a_i} (1-y_{i-1})^{1-a_i}$
## 2.3 Use of Proxy Polynomials

The construction of the multilinear polynomial $F^*(\cdot)$ (Eq.(2)) requires the summation over non-zero elements $F(a)$ of the truth table of size $2^k$. For non-trivial truth tables, writing such polynomial requires the use of symbolic computation software such as SymPy or Maple. But the complexity of such polynomial grows exponentially with the number $k$ of variables. In practice, when $k > 10$ for non-trivial filters the number of monomial becomes often prohibitively large.

When the computation of $F^*$ becomes intractable, we propose to use proxy polynomials $F^\bullet$ that have the same complexity as the binary filter $F()$. To this end, we first write the binary filter $F(a)$, $a \in \{0,1\}^k$ as a proposition formula involving the AND, OR and NOT logical operators. For non-trivial filters, this is done by efficiently enumerating the binary decision diagram representing the filter. We then substitute the logical operators with their linear polynomial equivalents: AND($Y_0, Y_1, \dots, Y_l$) $\equiv Y_0Y_1\dots Y_l$, OR($Y_0, Y_1, \dots, Y_l$) $\equiv 1 - (1 - Y_0)(1 - Y_1)\dots(1 - Y_l)$ and NOT($Y_0$) $= 1 - Y_0$. This results in a polynomial $F^\bullet(a)$ that coincides with $F$ and $F^*(\cdot)$ on the hypercube $\{0,1\}^k$: $F(a) = F^*(a) = F^\bullet(a), \forall a \in \{0,1\}^k$. The proxy polynomial $F^\bullet(x)$ comes in a factorized form, can be of degree greater than $k$ and usually is not multilinear. The only difference between the multilinear and proxy polynomials is the application of the idempotence rule of Boolean variables $Y_0^i = Y_0, \forall i > 0$ on all variables. Proxy polynomials can be interpreted as a non-linear interpolation of the values over the hypercube as opposed to a (multi) linear one for the multilinear polynomials. It can also be interpreted as a fuzzy logic translation of the binary filter using the product logic [14]. The proxy soft filters $F^\bullet(y_{N(1,n)}, \dots, y_{N(k,n)})$ approximates the true soft filter $F^*(y_{N(1,n)}, \dots, y_{N(k,n)})$ but remains computationally tractable with a complexity similar to the one of the binary filter. For instance the proxy polynomial for the dilation operator on the 4-neighbourhood written as $Z_{\text{Dil}} = 1 - (1 - y_0)(1 - y_1)(1 - y_2)(1 - y_3)(1 - y_4)$ which is multilinear and much more efficient to compute than the expression given in section 2.2.

## 2.4 Thinning Filters

Skeletonization aims at extracting a thinned version of an object, while preserving its essential shape, topology and connectivity [7]. The Boolean expressions utilized for skeletonization filters were derived from sub-iterative thinning algorithms from [25], for 2D images, and [16], for 3D images. These thinning algorithms operate on base masks representing the 8+1 or 26+1 local neighbourhood configuration of foreground pixels that can be turned into background in a specific orientation. Each sub-iteration of these algorithms corresponds to a distinct direction along which thinning operations are applied. In 2D, these directions include North, South, East, West, while in 3D they are extended to additional Up and Down directions. The Boolean expressions are translated into multilinear filters for 2D images and into proxy polynomials for 2D and 3D images. They are applied iteratively to the image until no change is observed between 2 iterations.
## 2.3 Use of Proxy Polynomials

The construction of the multilinear polynomial $F^*(\cdot)$ (Eq.(2)) requires the summation over non-zero elements $F(a)$ of the truth table of size $2^k$. For non-trivial truth tables, writing such polynomial requires the use of symbolic computation software such as SymPy or Maple. But the complexity of such polynomial grows exponentially with the number $k$ of variables. In practice, when $k > 10$ for non-trivial filters the number of monomial becomes often prohibitively large.

When the computation of $F^*$ becomes intractable, we propose to use proxy polynomials $F^\bullet$ that have the same complexity as the binary filter $F()$. To this end, we first write the binary filter $F(a)$, $a \in \{0,1\}^k$ as a proposition formula involving the AND, OR and NOT logical operators. For non-trivial filters, this is done by efficiently enumerating the binary decision diagram representing the filter. We then substitute the logical operators with their linear polynomial equivalents: AND($Y_0, Y_1, \dots, Y_l$) $\equiv Y_0Y_1\dots Y_l$, OR($Y_0, Y_1, \dots, Y_l$) $\equiv 1 - (1 - Y_0)(1 - Y_1)\dots(1 - Y_l)$ and NOT($Y_0$) $= 1 - Y_0$. This results in a polynomial $F^\bullet(a)$ that coincides with $F$ and $F^*(\cdot)$ on the hypercube $\{0,1\}^k$: $F(a) = F^*(a) = F^\bullet(a), \forall a \in \{0,1\}^k$. The proxy polynomial $F^\bullet(x)$ comes in a factorized form, can be of degree greater than $k$ and usually is not multilinear. The only difference between the multilinear and proxy polynomials is the application of the idempotence rule of Boolean variables $Y_0^i = Y_0, \forall i > 0$ on all variables. Proxy polynomials can be interpreted as a non-linear interpolation of the values over the hypercube as opposed to a (multi) linear one for the multilinear polynomials. It can also be interpreted as a fuzzy logic translation of the binary filter using the product logic [14]. The proxy soft filters $F^\bullet(y_{N(1,n)}, \dots, y_{N(k,n)})$ approximates the true soft filter $F^*(y_{N(1,n)}, \dots, y_{N(k,n)})$ but remains computationally tractable with a complexity similar to the one of the binary filter. For instance the proxy polynomial for the dilation operator on the 4-neighbourhood written as $Z_{\text{Dil}} = 1 - (1 - y_0)(1 - y_1)(1 - y_2)(1 - y_3)(1 - y_4)$ which is multilinear and much more efficient to compute than the expression given in section 2.2.

## 2.4 Thinning Filters

Skeletonization aims at extracting a thinned version of an object, while preserving its essential shape, topology and connectivity [7]. The Boolean expressions utilized for skeletonization filters were derived from sub-iterative thinning algorithms from [25], for 2D images, and [16], for 3D images. These thinning algorithms operate on base masks representing the 8+1 or 26+1 local neighbourhood configuration of foreground pixels that can be turned into background in a specific orientation. Each sub-iteration of these algorithms corresponds to a distinct direction along which thinning operations are applied. In 2D, these directions include North, South, East, West, while in 3D they are extended to additional Up and Down directions. The Boolean expressions are translated into multilinear filters for 2D images and into proxy polynomials for 2D and 3D images. They are applied iteratively to the image until no change is observed between 2 iterations.
# 3 Results

## 3.1 Experimental Setup

**Datasets.** Our framework was evaluated on three datasets: the DRIVE dataset, comprising 2D retinal blood vessel images [24]; the VesSap dataset [15], with synthetic 3D brain vessels, to benchmark 3D skeletonization methods; and a private dataset of lower limb CTA scans from the Hospital of Nice, with manual annotations of calcified plaques, for the second use case.

**Evaluation Metrics.** We assessed performance with multiple metrics: Dice co-efficient for mask overlaps, the clDice [22] for topology preservation, the mean of absolute Betti number errors representing the number of connected components ($\beta0$) and holes ($\beta1$), and the absolute error of Euler's characteristics.

## 3.2 Morphological filters Benchmarking

**Proxy polynomials analysis** We evaluate the correspondence between the output of multilinear and proxy polynomial-based filters for the soft skeletonization of a 2D probability map generated from one binary DRIVE dataset (Fig. 1b) as the sigmoid of a scaled distance map. Results reveal a strong correlation (Pearson correlation coefficient = 0.98, p-value < 0.001)(Fig. 1c) between both polynomials, with proxy polynomials being five times faster to compute than multilinear polynomials. Hence, proxy polynomials are efficient substitutes for multilinear functions in morphological filters with a small impact on the output.

**Benchmarking differentiable Morphological operations** Table 1 presents the results for erosion, dilation, and skeletonization operations on 2D and 3D images. We compare our dilation and erosion filters against the soft max and minpooling layer from Pytorch. For skeletonization, we compare our filters to a neural network trained for skeletonization [12], the skeletonization algorithm from [22] ("Soft Skeleton"), and the skeletonization method from Menten et al. [10]. The baseline methods are the non-differentiable filters implemented in scikit-image package [26]. Our method outperforms other methods across different metrics. In erosion and dilation operations, our approach achieves comparable Dice and topological scores with the baseline. For skeletonization, our method demonstrates improved performance in preserving vessel topology and connectivity compared to state-of-the-art methods, comparable with Menten et al. [10], but achieving a lower computation time in 2D and a better Dice in 3D. We acknowledge that the Dice may not reflect entirely the skeleton quality as different skeletonizations can produce different outputs. The soft-skeleton of [22] shows high topological errors due to discontinuities (Fig. 2).

## 3.3 Applications

**Use case 1: Skeleton extraction in clDice loss function for 2D segmentation.** We compared the integration of the skeleton extraction methods (based
# 3 Results

## 3.1 Experimental Setup

**Datasets.** Our framework was evaluated on three datasets: the DRIVE dataset, comprising 2D retinal blood vessel images [24]; the VesSap dataset [15], with synthetic 3D brain vessels, to benchmark 3D skeletonization methods; and a private dataset of lower limb CTA scans from the Hospital of Nice, with manual annotations of calcified plaques, for the second use case.

**Evaluation Metrics.** We assessed performance with multiple metrics: Dice coefficient for mask overlaps, the clDice [22] for topology preservation, the mean of absolute Betti number errors representing the number of connected components ($\beta0$) and holes ($\beta1$), and the absolute error of Euler's characteristics.

## 3.2 Morphological filters Benchmarking

**Proxy polynomials analysis** We evaluate the correspondence between the output of multilinear and proxy polynomial-based filters for the soft skeletonization of a 2D probability map generated from one binary DRIVE dataset (Fig. 1b) as the sigmoid of a scaled distance map. Results reveal a strong correlation (Pearson correlation coefficient = 0.98, p-value < 0.001)(Fig. 1c) between both polynomials, with proxy polynomials being five times faster to compute than multilinear polynomials. Hence, proxy polynomials are efficient substitutes for multilinear functions in morphological filters with a small impact on the output.

**Benchmarking differentiable Morphological operations** Table 1 presents the results for erosion, dilation, and skeletonization operations on 2D and 3D images. We compare our dilation and erosion filters against the soft max and minpooling layer from Pytorch. For skeletonization, we compare our filters to a neural network trained for skeletonization [12], the skeletonization algorithm from [22] ("Soft Skeleton"), and the skeletonization method from Menten et al. [10]. The baseline methods are the non-differentiable filters implemented in scikit-image package [26]. Our method outperforms other methods across different metrics. In erosion and dilation operations, our approach achieves comparable Dice and topological scores with the baseline. For skeletonization, our method demonstrates improved performance in preserving vessel topology and connectivity compared to state-of-the-art methods, comparable with Menten et al. [10], but achieving a lower computation time in 2D and a better Dice in 3D. We acknowledge that the Dice may not reflect entirely the skeleton quality as different skeletonizations can produce different outputs. The soft-skeleton of [22] shows high topological errors due to discontinuities (Fig. 2).

## 3.3 Applications

**Use case 1: Skeleton extraction in clDice loss function for 2D segmentation.** We compared the integration of the skeleton extraction methods (based
<table><caption>Table 1: Benchmark of morphological operation methods on 2D and 3D datasets</caption><thead><tr><th></th><th></th><th></th><th>Dice ↑</th><th>β0 ↓</th><th>β1 ↓</th><th>Euler ↓</th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="6">DRIVE<br>2D dataset</td><td rowspan="3">Dilation</td><td>Baseline</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.002</td></tr><tr><td>Ours</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.09</td></tr><tr><td>Soft maxpooling</td><td>0.95</td><td>0.1</td><td>6.65</td><td>6.75</td><td>0.003</td></tr><tr><td rowspan="3">Erosion</td><td>Baseline</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.002</td></tr><tr><td>Ours</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.08</td></tr><tr><td>Soft minpooling</td><td>0.85</td><td>214.1</td><td>7.6</td><td>211.8</td><td>0.003</td></tr><tr><td rowspan="5">Skeleton</td><td>Baseline</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.004</td></tr><tr><td>Ours</td><td>0.65</td><td>0</td><td>0</td><td>0</td><td>0.43</td></tr><tr><td>Neural Network</td><td>0.77</td><td>206.15</td><td>22.3</td><td>226.85</td><td>3.75</td></tr><tr><td>Soft-skeleton</td><td>0.65</td><td>1414.2</td><td>66.5</td><td>1480.7</td><td>0.005</td></tr><tr><td>Menten et al.</td><td>0.65</td><td>0</td><td>0</td><td>0</td><td>0.84</td></tr><tr><td rowspan="5">VesSap<br>3D dataset</td><td rowspan="5">Skeleton</td><td>Baseline</td><td>1</td><td>0</td><td>0</td><td>0</td><td>2.5</td></tr><tr><td>Ours</td><td>0.50</td><td>0</td><td>0</td><td>0</td><td>1800</td></tr><tr><td>Soft-skeleton</td><td>0.41</td><td>8362.2</td><td>8.8</td><td>8371</td><td>72</td></tr><tr><td>Menten et al.</td><td>0.48</td><td>0</td><td>0</td><td>0</td><td>880</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Fig. 2: Comparison of 2D skeletonization methods. Red arrows highlight significant differences.</figcaption></figure>

on proxy polynomials) into the clDice loss function [22] for 2D retinal vessel segmentation using the DRIVE dataset split in 80% training/ 20% testing sets. We used a 2D U-Net model [18] with a batch size of 16, 500 epochs and a learning rate of 1e-4. The clDice is used with the soft Dice loss, the weight α representing the weight of clDice in the loss function. Table 2 shows the performance metrics of different methods including the Binary Cross Entropy (BCE) and soft Dice losses alone. Our method yields the highest clDice scores, as well as the lowest number of β0 and Euler errors, indicating superior segmentation accuracy and topology preservation compared to baseline methods and existing techniques, demonstrating its effectiveness in preserving vessel continuity.

We also tested using morphological filters as the final layer of the 2D U-Net, after the sigmoid layer. We applied morphological dilation, erosion, closing and opening operations with a soft Dice loss (Table 2). The additional layers have an important impact on the topological fidelity of the segmentation, matching the effect of the clDice loss. Topological performance increases even more with the clDice loss added with our final morphological layer achieving the lowest number of β0 number errors and the best clDice score overall.
<table><caption>Table 1: Benchmark of morphological operation methods on 2D and 3D datasets</caption><thead><tr><th></th><th></th><th></th><th>Dice ↑</th><th>β0 ↓</th><th>β1 ↓</th><th>Euler ↓</th><th>Time (s)</th></tr></thead><tbody><tr><td rowspan="6">DRIVE<br>2D dataset</td><td rowspan="3">Dilation</td><td>Baseline</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.002</td></tr><tr><td>Ours</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.09</td></tr><tr><td>Soft maxpooling</td><td>0.95</td><td>0.1</td><td>6.65</td><td>6.75</td><td>0.003</td></tr><tr><td rowspan="3">Erosion</td><td>Baseline</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.002</td></tr><tr><td>Ours</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.08</td></tr><tr><td>Soft minpooling</td><td>0.85</td><td>214.1</td><td>7.6</td><td>211.8</td><td>0.003</td></tr><tr><td rowspan="5">Skeleton</td><td>Baseline</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0.004</td></tr><tr><td>Ours</td><td>0.65</td><td>0</td><td>0</td><td>0</td><td>0.43</td></tr><tr><td>Neural Network</td><td>0.77</td><td>206.15</td><td>22.3</td><td>226.85</td><td>3.75</td></tr><tr><td>Soft-skeleton</td><td>0.65</td><td>1414.2</td><td>66.5</td><td>1480.7</td><td>0.005</td></tr><tr><td>Menten et al.</td><td>0.65</td><td>0</td><td>0</td><td>0</td><td>0.84</td></tr><tr><td rowspan="5">VesSap<br>3D dataset</td><td rowspan="5">Skeleton</td><td>Baseline</td><td>1</td><td>0</td><td>0</td><td>0</td><td>2.5</td></tr><tr><td>Ours</td><td>0.50</td><td>0</td><td>0</td><td>0</td><td>1800</td></tr><tr><td>Soft-skeleton</td><td>0.41</td><td>8362.2</td><td>8.8</td><td>8371</td><td>72</td></tr><tr><td>Menten et al.</td><td>0.48</td><td>0</td><td>0</td><td>0</td><td>880</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Fig. 2: Comparison of 2D skeletonization methods. Red arrows highlight significant differences.</figcaption></figure>

on proxy polynomials) into the clDice loss function [22] for 2D retinal vessel segmentation using the DRIVE dataset split in 80% training/ 20% testing sets. We used a 2D U-Net model [18] with a batch size of 16, 500 epochs and a learning rate of 1e-4. The clDice is used with the soft Dice loss, the weight α representing the weight of clDice in the loss function. Table 2 shows the performance metrics of different methods including the Binary Cross Entropy (BCE) and soft Dice losses alone. Our method yields the highest clDice scores, as well as the lowest number of β0 and Euler errors, indicating superior segmentation accuracy and topology preservation compared to baseline methods and existing techniques, demonstrating its effectiveness in preserving vessel continuity.

We also tested using morphological filters as the final layer of the 2D U-Net, after the sigmoid layer. We applied morphological dilation, erosion, closing and opening operations with a soft Dice loss (Table 2). The additional layers have an important impact on the topological fidelity of the segmentation, matching the effect of the clDice loss. Topological performance increases even more with the clDice loss added with our final morphological layer achieving the lowest number of β0 number errors and the best clDice score overall.
<table><caption>Table 2: 2D U-Net segmentation performance on DRIVE dataset. (*) denotes statistically significant differences with soft-skeleton and Menten et al. methods having the same α value using the Wilcoxon rank test.</caption><thead><tr><th colspan="2">Loss \ Metrics</th><th>Dice ↑</th><th>clDice ↑</th><th>β0 ↓</th><th>β1 ↓</th><th>Euler ↓</th></tr></thead><tbody><tr><td colspan="2">BCE</td><td>0.806</td><td>0.853</td><td>48.55</td><td>40.75</td><td>89.3</td></tr><tr><td colspan="2">SoftDice (baseline)</td><td>0.82</td><td>0.863</td><td>45.35</td><td>35.15</td><td>80.5</td></tr><tr><td rowspan="3">clDice</td><td rowspan="3">α=0.7</td><td>Soft-Ske</td><td>0.789</td><td>0.885</td><td>17.9</td><td><b>27.2</b></td><td>44.6</td></tr><tr><td>Menten et al.</td><td>0.756</td><td>0.855</td><td>18.15</td><td>38.4</td><td>56.55</td></tr><tr><td>Ours</td><td>0.813</td><td><b>0.888</b></td><td><b>14.7</b></td><td>27.95</td><td><b>42.65</b></td></tr><tr><td rowspan="3"></td><td rowspan="3">α=0.5</td><td>Soft-Ske</td><td>0.823</td><td>0.879</td><td>24.7</td><td>32.05</td><td>56.75</td></tr><tr><td>Menten et al.</td><td>0.78</td><td>0.865</td><td>25.2</td><td>34.35</td><td>59.55</td></tr><tr><td>Ours</td><td><b>0.828</b></td><td>0.887*</td><td>16.5*</td><td>32.9</td><td>44.4</td></tr><tr><td rowspan="4">SoftDice + final layer</td><td rowspan="4"></td><td>Dilation</td><td><b>0.82</b></td><td><b>0.887</b></td><td><b>19.4</b></td><td>33.15</td><td><b>52.55</b></td></tr><tr><td>Erosion</td><td><b>0.82</b></td><td>0.866</td><td>26.9</td><td>35.55</td><td>62.45</td></tr><tr><td>Closing</td><td>0.8</td><td>0.874</td><td>30.55</td><td>29.55</td><td>60.1</td></tr><tr><td>Opening</td><td>0.78</td><td><b>0.886</b></td><td>28</td><td><b>25</b></td><td>52.9</td></tr><tr><td rowspan="3">clDice + final layer</td><td rowspan="3">Ours α=0.7</td><td>Dilation</td><td>0.8</td><td><b>0.899</b>*</td><td>11.3*</td><td>30.35</td><td>41.65</td></tr><tr><td>Erosion</td><td>0.821*</td><td>0.887</td><td><b>10.15</b>*</td><td>32.75</td><td>42.9</td></tr><tr><td>Ours α=0.5</td><td>0.807</td><td><b>0.897</b>*</td><td><b>11.1</b>*</td><td>30.55</td><td><b>41.55</b>*</td></tr><tr><td></td><td></td><td>Erosion</td><td><b>0.826</b></td><td>0.885</td><td><b>11.1</b>*</td><td>36.85</td><td>47.95</td></tr></tbody></table>

**Use case 2: Calcification plaque detection on CTAs of the lower limbs**

We tested the morphological filters as final layers of a 3D U-Net on the 3D dataset with a batch size of 2 and 250 epochs. The dataset contains 88 images split into 80% training/ 20% testing sets. In this dataset, we segment the calcification plaques in the arteries of the lower limbs. However, due to their small size, the Dice coefficient is not appropriate to evaluate the segmentation performances. Instead, we use detection metrics to compare the accuracy of detected calcification plaques using a threshold of 30% overlap between the prediction and ground truth. Without the morphological filters, detection performance achieve 71.6% average precision, 75.1% F1 score. Using as final layer the closing filter yields better results with 75.5% average precision, 77.1% F1 score, as well as the opening filter with 79.09% average precision and 78.58% F1 score. Optimal results are obtained with the erosion filter reaching an average precision of 91.3% and 80.8% F1 score.

## 4 Conclusion

We introduced a novel framework for differentiable soft morphological filtering in medical image analysis. We have demonstrated that multilinear polynomials associated with Boolean operations correspond to the expectation of binary morphological filters. In addition, we have introduced proxy polynomials that lead to efficient soft filters acting on probability maps. Finally, we have shown on two use-case applications that the proposed framework could be used effectively
<table><caption>Table 2: 2D U-Net segmentation performance on DRIVE dataset. (*) denotes statistically significant differences with soft-skeleton and Menten et al. methods having the same α value using the Wilcoxon rank test.</caption><thead><tr><th colspan="2">Loss \ Metrics</th><th>Dice ↑</th><th>clDice ↑</th><th>β0 ↓</th><th>β1 ↓</th><th>Euler ↓</th></tr></thead><tbody><tr><td colspan="2">BCE</td><td>0.806</td><td>0.853</td><td>48.55</td><td>40.75</td><td>89.3</td></tr><tr><td colspan="2">SoftDice (baseline)</td><td>0.82</td><td>0.863</td><td>45.35</td><td>35.15</td><td>80.5</td></tr><tr><td rowspan="3">clDice</td><td rowspan="3">α=0.7</td><td>Soft-Ske</td><td>0.789</td><td>0.885</td><td>17.9</td><td><b>27.2</b></td><td>44.6</td></tr><tr><td>Menten et al.</td><td>0.756</td><td>0.855</td><td>18.15</td><td>38.4</td><td>56.55</td></tr><tr><td>Ours</td><td>0.813</td><td><b>0.888</b></td><td><b>14.7</b></td><td>27.95</td><td><b>42.65</b></td></tr><tr><td rowspan="3"></td><td rowspan="3">α=0.5</td><td>Soft-Ske</td><td>0.823</td><td>0.879</td><td>24.7</td><td>32.05</td><td>56.75</td></tr><tr><td>Menten et al.</td><td>0.78</td><td>0.865</td><td>25.2</td><td>34.35</td><td>59.55</td></tr><tr><td>Ours</td><td><b>0.828</b></td><td>0.887*</td><td>16.5*</td><td>32.9</td><td>44.4</td></tr><tr><td rowspan="4">SoftDice + final layer</td><td rowspan="4"></td><td>Dilation</td><td><b>0.82</b></td><td><b>0.887</b></td><td><b>19.4</b></td><td>33.15</td><td><b>52.55</b></td></tr><tr><td>Erosion</td><td><b>0.82</b></td><td>0.866</td><td>26.9</td><td>35.55</td><td>62.45</td></tr><tr><td>Closing</td><td>0.8</td><td>0.874</td><td>30.55</td><td>29.55</td><td>60.1</td></tr><tr><td>Opening</td><td>0.78</td><td><b>0.886</b></td><td>28</td><td><b>25</b></td><td>52.9</td></tr><tr><td rowspan="3">clDice + final layer</td><td rowspan="3">Ours α=0.7</td><td>Dilation</td><td>0.8</td><td><b>0.899</b>*</td><td>11.3*</td><td>30.35</td><td>41.65</td></tr><tr><td>Erosion</td><td>0.821*</td><td>0.887</td><td><b>10.15</b>*</td><td>32.75</td><td>42.9</td></tr><tr><td>Ours α=0.5</td><td>Dilation</td><td>0.807</td><td><b>0.897</b>*</td><td><b>11.1</b>*</td><td>30.55</td><td><b>41.55</b>*</td></tr><tr><td></td><td></td><td>Erosion</td><td><b>0.826</b></td><td>0.885</td><td><b>11.1</b>*</td><td>36.85</td><td>47.95</td></tr></tbody></table>

**Use case 2: Calcification plaque detection on CTAs of the lower limbs**
We tested the morphological filters as final layers of a 3D U-Net on the 3D dataset with a batch size of 2 and 250 epochs. The dataset contains 88 images split into 80% training/ 20% testing sets. In this dataset, we segment the calcification plaques in the arteries of the lower limbs. However, due to their small size, the Dice coefficient is not appropriate to evaluate the segmentation performances. Instead, we use detection metrics to compare the accuracy of detected calcification plaques using a threshold of 30% overlap between the prediction and ground truth. Without the morphological filters, detection performance achieve 71.6% average precision, 75.1% F1 score. Using as final layer the closing filter yields better results with 75.5% average precision, 77.1% F1 score, as well as the opening filter with 79.09% average precision and 78.58% F1 score. Optimal results are obtained with the erosion filter reaching an average precision of 91.3% and 80.8% F1 score.

## 4 Conclusion

We introduced a novel framework for differentiable soft morphological filtering in medical image analysis. We have demonstrated that multilinear polynomials associated with Boolean operations correspond to the expectation of binary morphological filters. In addition, we have introduced proxy polynomials that lead to efficient soft filters acting on probability maps. Finally, we have shown on two use-case applications that the proposed framework could be used effectively
either inside loss functions or as final segmentation layers in order to preserve
topology and connectivity in segmentation problems.

Possible improvements include the optimization of 2D and 3D soft skele-
tonization filters, and the definition of new loss functions and layers ultimately
paving the way for its widespread use in clinical segmentation applications.

**Acknowledgments.** This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002. The authors are grateful to the OPAL infrastructure from Université Côte d'Azur.

References

1. Aouad, T., Talbot, H.: Binary morphological neural network. In: 2022 IEEE Inter-national Conference on Image Processing (ICIP). pp. 3276–3280 (2022)

2. Asano, A., Kobayashi, Y., Chie Muraki, Muneyasu, M.: Optimization of gray scale morphological opening for noise removal in texture images. In: The 2004 47th Midwest Symposium on Circuits and Systems, 2004. MWSCAS '04. vol. 1, pp. 1.313–1.316. IEEE, Hiroshima, Japan (2004)

3. Franchi, G., Fehri, A., Yao, A.: Deep morphological networks. Pattern Recognition **102**, 107246 (Jun 2020)

4. González-Hidalgo, M., Massanet, S., Mir, A., Ruiz-Aguilera, D.: Improving salt and pepper noise removal using a fuzzy mathematical morphology-based filter. Applied Soft Computing **63**, 167–180 (Feb 2018)

5. Jurdi, R.E., Petitjean, C., Honeine, P., Cheplygina, V., Abdallah, F.: A Surprisingly Effective Perimeter-based Loss for Medical Image Segmentation. In: Proceedings of the Fourth Conference on Medical Imaging with Deep Learning. pp. 158–167. PMLR (Aug 2021)

6. Krishnapuram, R., Gupta, S.: Morphological methods for detection and classification of edges in range images. Journal of Mathematical Imaging and Vision **2**(4), 351–375 (Dec 1992)

7. Lam, L., Lee, S.W., Suen, C.: Thinning methodologies-a comprehensive survey. IEEE Transactions on Pattern Analysis and Machine Intelligence **14**(9), 869–885 (1992)

8. LeCun, Y., Bengio, Y., Hinton, G.: Deep learning. Nature **521** (7553), 436–444 (May 2015)

9. Lidayová, K., Frimmel, H., Wang, C., Bengtsson, E., Smedby, O.: Chapter 12 - Skeleton-based fast, fully automated generation of vessel tree structure for clinical evaluation of blood vessel systems. In: Skeletonization, pp. 345–382. Academic Press (Jan 2017)

10. Menten, M.J., Paetzold, J.C., Zimmer, V.A., Shit, S., Ezhov, I., Holland, R., Probst, M., Schnabel, J.A., Rueckert, D.: A skeletonization algorithm for gradient-based optimization. In: Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV). pp. 21394–21403 (October 2023)

11. Mondal, R., Dey, M.S., Chanda, B.: Image Restoration by Learning Morphological Opening-Closing Network. Mathematical Morphology - Theory and Applications **4**(1), 87–107 (Jan 2020)
either inside loss functions or as final segmentation layers in order to preserve
topology and connectivity in segmentation problems.

Possible improvements include the optimization of 2D and 3D soft skele-
tonization filters, and the definition of new loss functions and layers ultimately
paving the way for its widespread use in clinical segmentation applications.

**Acknowledgments.** This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002. The authors are grateful to the OPAL infrastructure from Université Côte d'Azur.

References

1. Aouad, T., Talbot, H.: Binary morphological neural network. In: 2022 IEEE Inter-national Conference on Image Processing (ICIP). pp. 3276–3280 (2022)

2. Asano, A., Kobayashi, Y., Chie Muraki, Muneyasu, M.: Optimization of gray scale morphological opening for noise removal in texture images. In: The 2004 47th Midwest Symposium on Circuits and Systems, 2004. MWSCAS '04. vol. 1, pp. 1.313–1.316. IEEE, Hiroshima, Japan (2004)

3. Franchi, G., Fehri, A., Yao, A.: Deep morphological networks. Pattern Recognition **102**, 107246 (Jun 2020)

4. González-Hidalgo, M., Massanet, S., Mir, A., Ruiz-Aguilera, D.: Improving salt and pepper noise removal using a fuzzy mathematical morphology-based filter. Applied Soft Computing **63**, 167–180 (Feb 2018)

5. Jurdi, R.E., Petitjean, C., Honeine, P., Cheplygina, V., Abdallah, F.: A Surprisingly Effective Perimeter-based Loss for Medical Image Segmentation. In: Proceedings of the Fourth Conference on Medical Imaging with Deep Learning. pp. 158–167. PMLR (Aug 2021)

6. Krishnapuram, R., Gupta, S.: Morphological methods for detection and classification of edges in range images. Journal of Mathematical Imaging and Vision **2**(4), 351–375 (Dec 1992)

7. Lam, L., Lee, S.W., Suen, C.: Thinning methodologies-a comprehensive survey. IEEE Transactions on Pattern Analysis and Machine Intelligence **14**(9), 869–885 (1992)

8. LeCun, Y., Bengio, Y., Hinton, G.: Deep learning. Nature **521**(7553), 436–444 (May 2015)

9. Lidayová, K., Frimmel, H., Wang, C., Bengtsson, E., Smedby, O.: Chapter 12 - Skeleton-based fast, fully automated generation of vessel tree structure for clinical evaluation of blood vessel systems. In: Skeletonization, pp. 345–382. Academic Press (Jan 2017)

10. Menten, M.J., Paetzold, J.C., Zimmer, V.A., Shit, S., Ezhov, I., Holland, R., Probst, M., Schnabel, J.A., Rueckert, D.: A skeletonization algorithm for gradient-based optimization. In: Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV). pp. 21394–21403 (October 2023)

11. Mondal, R., Dey, M.S., Chanda, B.: Image Restoration by Learning Morphological Opening-Closing Network. Mathematical Morphology - Theory and Applications **4**(1), 87–107 (Jan 2020)
12. Nguyen, N.H.: U-Net based skeletonization and bag of tricks. In: 2021 IEEE/CVF International Conference on Computer Vision Workshops (ICCVW). pp. 2105–2109. IEEE, Montreal, BC, Canada (Oct 2021)

13. Nogueira, K., Chanussot, J., Mura, M.D., Santos, J.A.D.: An introduction to deep morphological networks. IEEE Access 9, 114308–114324 (2021)

14. Novák, V., Perfilieva, I., Mockor, J.: Mathematical Principles of Fuzzy Logic. The Springer International Series in Engineering and Computer Science, Springer US (1999)

15. Paetzold, J., Schoppe, O., Al-Maskari, R., Tetteh, G., Efremov, V., Todorov, M.I., Cai, R., Mai, H., Rong, Z., Ertuerk, A., Menze, B.H.: Transfer learning from synthetic data reduces need for labels to segment brain vasculature and neural pathways in 3D (Apr 2019)

16. Palàgyi, K., Kuba, A.: A 3D 6-subiteration thinning algorithm for extracting medial lines. Pattern Recognition Letters 19(7), 613–627 (May 1998)

17. Panichev, O., Voloshyna, A.: U-Net Based Convolutional Neural Network for Skeleton Extraction. In: 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW). pp. 1186–1189. IEEE, Long Beach, CA, USA (Jun 2019)

18. Ronneberger, O., Fischer, P., Brox, T.: U-Net: Convolutional Networks for Biomedical Image Segmentation. In: Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015. pp. 234–241. Lecture Notes in Computer Science, Springer International Publishing, Cham (2015)

19. Serra, J., Soille, P., Viergever, M.A.: Mathematical Morphology and Its Applications to Image Processing, Computational Imaging and Vision, vol. 2. Springer Netherlands, Dordrecht (1994)

20. Shen, Y., Shih, F.Y., Zhong, X., Chang, I.C.: Deep Morphological Neural Networks. International Journal of Pattern Recognition and Artificial Intelligence 36(12), 2252023 (Sep 2022)

21. Shih, F.Y.: Image Processing and Mathematical Morphology: Fundamentals and Applications. CRC Press, Boca Raton (Jan 2017)

22. Shit, S., Paetzold, J.C., Sekuboyina, A., Ezhov, I., Unger, A., Zhylka, A., Pluim, J.P.W., Bauer, U., Menze, B.H.: clDice - a Novel Topology-Preserving Loss Function for Tubular Structure Segmentation. pp. 16555–16564. IEEE Computer Society (Jun 2021)

23. Singh, G., Singh, S.: Medical Images Edge Detection using Mathematical Based Morphology. Apeejay Journal Of Computer Science And Applications (2013)

24. Staal, J., Abramoff, M., Niemeijer, M., Viergever, M., Van Ginneken, B.: Ridge-Based Vessel Segmentation in Color Images of the Retina. IEEE Transactions on Medical Imaging 23(4), 501–509 (Apr 2004)

25. Wagner, M.G.: Real-Time Thinning Algorithms for 2D and 3D Images using GPU processors. Journal of real-time image processing 17(5), 1255–1266 (Oct 2020)

26. Van der Walt, S., Schönberger, J.L., Nunez-Iglesias, J., Boulogne, F., Warner, J.D., Yager, N., Gouillart, E., Yu, T.: scikit-image: image processing in python. PeerJ 2, e453 (2014)

27. Zhao Yu-qian, Gui Wei-hua, Chen Zhen-cheng, Tang Jing-tian, Li Ling-yun: Medical Images Edge Detection Based on Mathematical Morphology. In: 2005 IEEE Engineering in Medicine and Biology 27th Annual Conference. pp. 6492–6495. IEEE, Shanghai, China (2005)
12. Nguyen, N.H.: U-Net based skeletonization and bag of tricks. In: 2021 IEEE/CVF International Conference on Computer Vision Workshops (ICCVW). pp. 2105–2109. IEEE, Montreal, BC, Canada (Oct 2021)

13. Nogueira, K., Chanussot, J., Mura, M.D., Santos, J.A.D.: An introduction to deep morphological networks. IEEE Access 9, 114308–114324 (2021)

14. Novák, V., Perfilieva, I., Mockor, J.: Mathematical Principles of Fuzzy Logic. The Springer International Series in Engineering and Computer Science, Springer US (1999)

15. Paetzold, J., Schoppe, O., Al-Maskari, R., Tetteh, G., Efremov, V., Todorov, M.I., Cai, R., Mai, H., Rong, Z., Ertuerk, A., Menze, B.H.: Transfer learning from synthetic data reduces need for labels to segment brain vasculature and neural pathways in 3D (Apr 2019)

16. Palàgyi, K., Kuba, A.: A 3D 6-subiteration thinning algorithm for extracting medial lines. Pattern Recognition Letters 19(7), 613–627 (May 1998)

17. Panichev, O., Voloshyna, A.: U-Net Based Convolutional Neural Network for Skeleton Extraction. In: 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW). pp. 1186–1189. IEEE, Long Beach, CA, USA (Jun 2019)

18. Ronneberger, O., Fischer, P., Brox, T.: U-Net: Convolutional Networks for Biomedical Image Segmentation. In: Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015. pp. 234–241. Lecture Notes in Computer Science, Springer International Publishing, Cham (2015)

19. Serra, J., Soille, P., Viergever, M.A.: Mathematical Morphology and Its Applications to Image Processing, Computational Imaging and Vision, vol. 2. Springer Netherlands, Dordrecht (1994)

20. Shen, Y., Shih, F.Y., Zhong, X., Chang, I.C.: Deep Morphological Neural Networks. International Journal of Pattern Recognition and Artificial Intelligence 36(12), 2252023 (Sep 2022)

21. Shih, F.Y.: Image Processing and Mathematical Morphology: Fundamentals and Applications. CRC Press, Boca Raton (Jan 2017)

22. Shit, S., Paetzold, J.C., Sekuboyina, A., Ezhov, I., Unger, A., Zhylka, A., Pluim, J.P.W., Bauer, U., Menze, B.H.: clDice - a Novel Topology-Preserving Loss Function for Tubular Structure Segmentation. pp. 16555–16564. IEEE Computer Society (Jun 2021)

23. Singh, G., Singh, S.: Medical Images Edge Detection using Mathematical Based Morphology. Apeejay Journal Of Computer Science And Applications (2013)

24. Staal, J., Abramoff, M., Niemeijer, M., Viergever, M., Van Ginneken, B.: Ridge-Based Vessel Segmentation in Color Images of the Retina. IEEE Transactions on Medical Imaging 23(4), 501–509 (Apr 2004)

25. Wagner, M.G.: Real-Time Thinning Algorithms for 2D and 3D Images using GPU processors. Journal of real-time image processing 17(5), 1255–1266 (Oct 2020)

26. Van der Walt, S., Schönberger, J.L., Nunez-Iglesias, J., Boulogne, F., Warner, J.D., Yager, N., Gouillart, E., Yu, T.: scikit-image: image processing in python. PeerJ 2, e453 (2014)

27. Zhao Yu-qian, Gui Wei-hua, Chen Zhen-cheng, Tang Jing-tian, Li Ling-yun: Medical Images Edge Detection Based on Mathematical Morphology. In: 2005 IEEE Engineering in Medicine and Biology 27th Annual Conference. pp. 6492–6495. IEEE, Shanghai, China (2005)