Understanding the Impact of Evaluation Metrics in Kinetic Models for Consensus-based Segmentation

We now consider the bidimensional case, $d=2$, and we specify the interaction function based on the so-called bounded confidence model. In more detail, we consider $N\geq 2$ agents and define their opinion variable through a vector $\mathbf{x}=(x_{i}(t),y_{i}(t))\in\mathbf{R}^{2}$, characterised by initial states $\{\mathbf{x}_{1}(0),\dots,\mathbf{x}_{N}(0)\}$. Agents will modify their opinion as a result of the interaction with other agents only if $|\mathbf{x}_{i}-\mathbf{x}_{j}|\leq\Delta$ , where $\Delta\geq 0$ is a given confidence level. Hence, we can write (1) as follows

where $P_{\Delta}(\mathbf{x}_{i},\mathbf{x}_{j})=\chi(|\mathbf{x}_{i}-\mathbf{x}_{j}|\leq\Delta):\mathbb{R}^{2}\to\{0,1\}$, being $\chi(A)$ the characteristic function of the set $A\subseteq\mathbb{R}^{2}$. We can easily observe that the mean position of the ensemble of agents is conserved in time, indeed

thanks to the symmetry of the considered bounded-confidence interaction function. The bounded confidence model converges to a steady configuration, meaning that the systems reaches consensus in finite time. The structure of the steady state depends on the value of $\Delta$, see [43].

Furthermore, to account for random fluctuations given by external factors in the opinion of agents we may consider a diffusion component as follows

where $\{\mathbf{W}_{i}\}_{i=1}^{N}$ is a set of independent Wiener processes. The impact of the diffusion is weighted by the variable $\sigma^{2}>0$. To visualise the interplay between consensus forces and diffusion we depict in Figure 1 the steady configuration of the model (4) for different combinations of the model parameters. For $\sigma^{2}=0$, the system forms a finite number of clusters depending on the value of $\Delta>0$, as illustrated in Figure 1(a). For values of the diffusion coefficient $\sigma^{2}>0$, the number of clusters of the system varies as depicted in Figure 1(b). The right panel of Figure 1(b) shows the scenario in which the diffusion effect becomes comparable to the tendency of agents to cluster. Finally, in Figure 1(c), for $\sigma=0.05$, the diffusion effect dominates the grouping tendency, resulting in a homogeneous steady-state distribution.

In the limit $N\to+\infty$ it can be shown that the empirical density

of the system of particles (4) converges to a continuous density $f(\textbf{x},t):\mathbb{R}^{2}\times\mathbb{R}_{+}\to\mathbb{R}_{+}$ solution to the following mean-field equation

where $\Xi[f](\mathbf{x},t)$ is defined as follows

see e.g. [11].

We can derive (5) using a kinetic approach by writing $\textbf{x}:=\textbf{x}_{i}(t)$ and $\textbf{x}_{*}:=\textbf{x}_{j}(t)$ for a generic pair $(i,j)$ of interacting agents/particles and we approximate the time derivative in (4) in a time step $\epsilon=\Delta t>0$, through an Euler-Maryuama approach, in the same spirit as [10, 45]. Hence, we recover the binary interaction rule

where $\mathbf{x}^{\prime}=\mathbf{x}_{i}(t+\epsilon)$, $\mathbf{x}_{*}^{\prime}=\mathbf{x}_{j}(t+\epsilon)$ and $\eta,\eta_{*}$ are two independent 2D centered Gaussian distribution random variable such that

where $\langle\cdot\rangle$ denotes the integration with respect to the distribution $\eta$. Furthermore, in (7) we shall consider $P(\mathbf{x},\mathbf{x}_{*})=\chi(|\mathbf{x}-\mathbf{x_{*}}|<\Delta)$. We can remark that, if $\sigma=0$, since $P_{\Delta}\in[0,1]$ and $\epsilon\in(0,1)$ we get

since the interaction function $P_{\Delta}$ is symmetric, consistently with (3). This shows that the mean position is conserved at every interaction. Finally, we have

and the mean energy is dissipated at each interaction, since $P_{\Delta}\geq 0$. Hence, we consider the distribution function $f=f(\mathbf{x},t):R^{2}\times R_{+}\rightarrow R_{+}$, such that $f(\mathbf{x},t)d\mathbf{x}$ represents the fraction of agents/particles in $[{x}_{1},{x}_{1}+d{x}_{1})\times[{x}_{2},{x}_{2}+d{x}_{2}]$ at time $t\geq 0$. The evolution of $f$ as a result of binary-interaction scheme (7) is obtained by a Boltzmann-type equation which reads in weak form

being $\varphi(\cdot)$ a test function. As observed in [54], when $\Delta t=\epsilon\to 0^{+}$ we can observe that the binary scheme (7) becomes quasi-invariant and we can introduce the following expansion

being $R_{\epsilon}(\mathbf{x},\mathbf{x}_{*})$ a reminder term and $H[\varphi]$ the Hessian matrix. Hence, scaling $\tau=\epsilon t$ and the distribution $f_{\epsilon}(\mathbf{x},\tau)=f(\mathbf{x},\tau/\epsilon)$, we may plug (12) in (11) to get

Following [9], see also [42], we can prove that

as $\epsilon\to 0^{+}$. Hence, integrating back by parts the first two terms we obtain (5). In more detail, we can prove that $f_{\epsilon}$ converges, up to extraction of a subsequence to a probability density $f(\mathbf{x},\tau)$ that is weak solution to the nonlocal Fokker-Planck equation (5).

An application of the Hegselman-Krause model for data-clustering problems has been proposed in [26]. The idea is to extend the 2D model by characterising each particle with an internal feature $c_{i}\in[0,1]$ that represents the gray color of the $i$th pixel. Therefore, we interpret each pixel in the image as a particle characterized by a position vector and the static feature $c$ as shown in Figure 2.

To address the segmentation task, we can define a dynamic feature for the system of pixels through an interaction function that accounts for alignment processes among pixels with sufficiently similar features. In particular, let us consider the following

Therefore, the time-continuous evolution for the system of pixels is given by

In this case we introduced two confidence bounds $\Delta_{1}\geq 0$, $\Delta_{2}\geq 0$ taking into account the position and the gray level of the pixels, respectively. In this way, the interactions between the pixels will generate a large time distribution which is characterized by several clusters, depending on the value of $\Delta_{1}$ and $\Delta_{2}$. Hence, coherently with k-means methods, see e.g. [37], a pixel belongs to a cluster $\mathcal{C}_{\mu}=\{\mathbf{x}_{i}:\|\mathbf{x}_{i}-\mu\|\leq\alpha\}$, being $\alpha>0$ the pixel size, if it is sufficiently close to the local quantity $\mu\in\mathbb{R}^{2}$. We highlight how we are only interested in clustering with respect to the space variable. This dynamics is represented in Figure 3.

Biomedical images are often subject to ambiguities arising from various sources of uncertainty related to clinical factors and potential bottlenecks in data acquisition processes [5, 32]. These uncertainties can be broadly categorized into aleatoric uncertainty, stemming from inherent stochastic variations in the data collection process, and epistemic uncertainty, relates to uncertainties in model parameters and can lead to deviations in the results. Aleatoric uncertainties poses significant challenges in image segmentation, as image processing models must contend with limitations in the raw acquisition data. Addressing these uncertainties is critical, and the study of uncertainty quantification in image segmentation is an expanding field aimed at developing robust segmentation algorithms capable of mitigating erroneous outcomes. To this end, in [9], it has been proposed an extension of (14) to consider segmentation of biomedical images. In particular, the particle model (15) has been integrated a nonconstant stochastic part to take into account aleatoric uncertainties arising from the data acquisition process. These uncertainties may include factors such as motion artifacts or field inhomogeneities in magnetic resonance imaging (MRI). They modified equation 14 as follows:

where $\{\mathbf{W}_{i}\}_{i=1}^{N}$ is set of independent Wiener processes, $P_{\Delta_{1},\Delta_{2}}(\cdot,\cdot,\cdot,\cdot)\in[0,1]$ is the interaction function defined in (13), and $D(c)\geq 0$ quantifies the impact of diffusion related to the value of the feature $c\in[0,1]$. Since the aleatoric uncertainties are expected to appear far away from the static feature’s boundaries, only diffusion functions that are maximal at the center and satisfy $D(0)=D(1)=0$ are considered. Similarly to (7), we may introduce the following binary interaction scheme by writing $(\mathbf{x},\mathbf{x}_{*}):=(\mathbf{x}_{i}(t),\mathbf{x}_{j}(t))$ a random couple of pixels having features $(c,c_{*}):=(c_{i}(t),c_{j}(t))$. We get

where $(\mathbf{x}^{\prime},\mathbf{x}_{*}^{\prime}):=(\mathbf{x}_{i}(t+\Delta t),\mathbf{x}_{j}(t+\Delta t))$ and $(c^{\prime},c_{*}^{\prime}):=(c_{i}(t+\Delta t),c_{j}(t+\Delta t))$. At the statistical level, as in [9], we may follows the approach described in Section 2.2. Hence, we introduce the distribution function $f=f(\mathbf{x},c,t):R^{2}\times[0,1]\times R_{+}\rightarrow R_{+}$, such that $f(\mathbf{x},t)d\mathbf{x}$ represents the fraction of agents/particles in $[x_{1},x_{1}+dx_{1})\times[x_{2},x_{2}+dx_{2}]$ characterized by a feature $c\in[0,1]$ at time $t\geq 0$. The evolution of $f$ whose interaction follow the binary scheme (16) is given by the following Boltzmann-type equation

Hence, since the feature is not evolving in time, we can proceed as in Section 2.2 to derive in the quasi-invariant limit for $\epsilon\to 0^{+}$ the corresponding Fokker-Planck-type PDE

where

The numerical approximation of Boltzmann-type equations has been deeply investigated in the recent decades, see e.g. [19, 41]. The approximation of this class of equations is particularly challenging due to the curse of dimensionality brought up by the multidimensional integral of the collision operator, and the presence of multiple scales. Furthermore, the preservation of relevant physical quantities are essential for a correct description of the underlying physical problem [44].

In view of its computational efficiency, in the following we will adopt a DSMC approach. Indeed the computational cost of this method is $O(N)$ where $N$ is the number of particles. Next, we describe the DSMC method based on a Nanbu-Bavosky scheme [41]. We begin by randomly selecting $N/2$ pairs of particles and making them evolve following the binary scheme presented in (7). We consider a time interval $[0,T]$ which we divide in $N_{t}$ intervals of size $\Delta t>0$. The DSMC approach for the introduced kinetic equation is based on a first order forward time discretization. In the following we will always consider the case $\Delta t=\epsilon>0$ such that all the particles are going to interact, see [41] for more details. We introduce the stochastic rounding of a positive real number x as:

where $\lfloor x\rfloor$ is the integer part of $x$. The random variable $\eta$ is sampled from a 2D Gaussian Distribution centered at zero and a diagonal covariance matrix.

In this section we present the procedure to estimate the Segmentation Mask of Brain Tumor Images. The procedure described in this section closely follows the methodology presented in [9]. For a given image, we define the feature’s values in relation to the gray level of each pixel. In more detail, for a given pixel $i\in\{1,\dots,N\}$ we define

being $C_{i}$, $i=1,\dots,N$, the gray value of the original image. Therefore, the value $c_{i}=1$ represents a white pixel and $c_{i}=0$ represents black pixel. In particular, for this work we used the brain tumor dataset that consists of 3D in multi-parametric MRI of patients affected by glioblastoma or lower-grade glioma, publicly available in the context of the Brain Tumor Image Segmentation Challenge http://medicaldecathlon.com/. The acquisition sequences include $T_{1}$-weighted, post-Gadolinium contrast $T_{1}$-weighted, $T_{2}$-weighted and $T_{2}$ Fluid-Attenuated Inversion Recovery volumes. Each MRI scan is accompanied by corresponding ground truth segmentation mask, which is a binary image where anatomical regions of interest are highlighted as white pixels while all other areas are represented as black pixels. These ground truth segmentation masks were manually delineated by experienced radiologists and specifically identify three structures: ”tumor core”, ”enhancing tumor” and ”whole tumor”. We evaluate the performances of the DSMC algorithm for two different segmentation tasks: the ”tumor core” and the ”whole tumor” annotations. For the first task we use a single slice in the axial plane of the post-Gadolinium contrast $T_{1}$-weighted scans while for the second task we use a single slice in the axial plane of the $T_{2}$-weighted scans. The procedure to generate the segmentation masks is as follows:

1. 1.

   We begin by associating each pixel with a position vector $(x_{i},y_{i})$ and with static feature $c_{i}$. We scale the vector position to a domain $[-1,1]\times[-1,1]$ and the static feature to $[0,1]$.

2. 2.

   We apply a DSMC approach as described in Algorithm 1 to numerically approximate the large-time solution of the Boltzmann-type model defined in (11). This approach enables pixels to aggregate into clusters based on their Euclidean distance and gray color level.

3. 3.

   The segmentation masks are generated by assigning to the original position of each pixel the mean value of the clusters they belong to. In this way we generate a multi-level mask composed of a number of homogenous regions.

4. 4.

   Finally we obtain the binary mask by defining a threshold $\tilde{c}$ such that:

   $$c_{i}=\begin{cases}1\ if\ c\geq\tilde{c}\\ 0\ if\ c<\tilde{c}\end{cases}$$

   (22)

   For all the following experiments, $\tilde{c}$ is defined as the 10th percentile of pixels in the image that belong to the region of interest. This percentile was chosen as an optimal value for brain tumor images; however, it could also be considered as a parameter to be optimized within the process outlined in Section 3.2.1.

Following this procedure, we apply two morphological refinement steps to remove small regions that have been misclassified as foreground parts and to fill small regions that have been incorrectly categorized as background pixels. We begin by labeling all the connected pixels in the foreground and reassigning them to the background those whose number of pixels is less than a certain threshold. Then we repeat the same procedure but for the pixels in the background. To this end, we use the scikit-image python library that detects distinct objects of a binary image [55]. This allows us to obtain more precise segmentation masks by reducing small imperfections. This entire process is illustrated in Figure 4.

Next, we introduce the principal optimization metrics used for evaluating a binary segmentation mask. We define $\{S_{g}^{0},S_{g}^{1},S_{t}^{0},S_{t}^{1}\}$ where $S_{g}^{0}$ and $S_{g}^{1}$ represent the set of pixels that belong to the background and foreground of the ground truth segmentation mask respectively. Same applies for ${S_{t}^{0},S_{t}^{1}}$ but for the binary mask we want to evaluate. One could also wish to asses the validity of a segmentation mask with multiple labels, we refer to [53] for an introduction to the subject. Figure 5 presents a summary of the key terms used in the definitions of metrics.

In this section we study the impact of choosing different diffusion functions $D(c)$ in images consisting of a blurry background and a geometric shape in the center, as shown in Figure 7. The objective is to detect the shape of the Geometric Figure and to compare how the choice of different diffusion functions affect the value of the model parameters $\Delta_{1},\Delta_{2}$ and $\sigma^{2}$ where the optimization process is identical to the one introduce in Section 3.2.1. To this end, we chose the following diffusion functions:

We point the reader to Figure 6 for a summary of the various introduced diffusion functions in (34).

For both the square and circle images, the Surface Dice Coefficient was used to optimize the parameters with a tolerance equal to the length of 1 pixel. Both images have a shape of $(256,256)$ pixels. The final time was set to $T=200$ with $\Delta t=0.1$. The resulting binary mask was the same for all choices of diffusion functions, obtaining the same loss function value. The results are shown in Figure 7. In the case of the square Figure 7 (a) we can see from Table 1 that for $D_{1}(c)$ and $D_{3}(c)$ the values of $\Delta_{1}$ do not differ greatly for this two diffusion functions. In the case of $\Delta_{2}$ we obtain a slightly smaller value for $D_{1}(c)$ compared to the one obtained for $D_{3}(c)$ and a larger value of the parameter $\sigma^{2}>0$ for $D_{3}(c)$ compared to the one obtained for $D_{1}(c)$. If we look at Figure 6 we notice that $D_{1}(c)\geq D_{3}(c)$. Therefore, a larger value of the diffusion functions is balanced by a smaller value of $\sigma^{2}$ to obtain a similar diffusion effect. This holds also for $D_{1}(c)$ and $D_{3}(c)$ for the circle Figure 7 (b). Furthermore, comparing $D_{2}(c)$ and $D_{4}(c)$ for the square image we can see that the resulting parameters are smaller for $D_{2}(c)$ in contrast to the one obtained with $D_{4}(c)$. This is consistent because, again, we can see from Figure 6 that $D_{2}(c)\geq D_{4}(c)$. If we now compare $D_{2}(c)$ and $D_{4}(c)$ for the circle image we can see that the value of $\sigma^{2}$ is similar in this case. Nevertheless, in this case the difference is given by the values of $\Delta_{1}$ and $\Delta_{2}$, which are both smaller for $D_{2}(c)$. This indicates that, for different diffusion functions, the optimal parameters adjust to yield similar results. A very straight way is to obtain similar values of $\Delta_{1}$ and $\Delta_{2}$ and a lower value of $\sigma^{2}$ for the diffusion function that has a higher value as in the case of the square image. However, the example of the circle image shows that us that we can also obtain different combinations of parameters so as to counter the effect of a bigger diffusion function.

From Table 2 we can see the parameters obtained by minimizing three different optimization metrics using as a diffusion function $D_{1}(c)$ for the square image. For all the cases the resulting Surface Dice was equal to one indicating a perfect overlap between the computed and the ground truth segmentation masks. The resulting binary mask obtained were the same for the three examples and are equivalent to the ones shown in Figure 7. For the Volumetric and Surface Dice Coefficient we can see that the parameters obtained where identical. Nevertheless, for the Jaccard Index, the resulting parameters differed, being smaller in this case. The loss is null in both cases, coherently with the relationship (29).

In this section we specify the criterion that we implemented to determine the final time $T>0$. As defined in Section 3.1, we approximate the solution of (18) through a DSMC approach, even though we have no analytical insight on the form of the steady state. The objective is to find the values of the final time $T>0$ such that a numerical steady state can be defined. We stress that the time taken to reach the equilibrium state for different initial conditions is not the same so we need to determine the time parameters for all the images we want to analyze. To this end, if $f^{n}(\mathbf{x},c)$ is the approximation of the density at time $t^{n}=n\Delta t$, we define

which represents an index of variation between two successive time steps of the reconstructed kinetic density. As the solution evolves, this quantity decreases and tends to zero as the equilibrium state is reached, as illustrated in Figure 8 for the case of the Square Image with a blurry background. Hence, we may introduce a breaking criterion based on the condition $\mathcal{T}<\delta$, for some $\delta>0$. When this condition is satisfied, the reconstructed density is considered an approximation of the steady state.

The same procedure was done for all images presented in this work so as to fulfill with the condition presented in this section.

In this section we study the impact that the different optimization metrics have on the resulting binary mask for the Core and Whole Tumor. We also analyze the parameters obtained for the different optimization metrics. Both brain tumor images consist of $N=(240,240)$ pixels. For the optimization procedure we determine $T=100$ and $\Delta t=0.01$. For each segmentation mask generated we evaluated $300$ different combinations of parameters. Figure 9 show the Segmentation Masks obtained for both the Whole and Core Tumor by optimizing the Jaccard Index and the Volumetric Dice Coefficient. In Table 3 the resulting parameters and the loss obtained for both optimization metrics are presented, in this case the loss is equal to 1 for a perfect overlap and 0 if the images are totally disjoint. First, we can observe that the loss obtained with both metrics satisfy (29) as expected. It can be noticed that for both segmentation masks the loss obtained is greater for the Volumetric Dice Coefficient. Furthermore, the parameter $\Delta_{1}$ obtained with both optimization metrics is similar for both the Core and Whole Tumor. Nevertheless, we can see that for the Whole Tumor the $\Delta_{2}$ parameter obtained with the Jaccard Index is bigger than the one obtained with the Volumetric Dice Coefficient. For the case of the Core Tumor instead the $\Delta_{2}$ parameter is bigger for the Volumetric Dice Coefficient. If we compare this to the values obtained for $\sigma^{2}$ in both cases for both metrics we can see that a bigger diffusion value is countered by a smaller value of $\Delta_{2}$ so as to obtain similar Segmentation Masks as seen from Figure 9.

For the Surface Dice Coefficient, the tolerance $\tau$ was set to the length of 1 pixel, both when used as the optimization loss and when used as the evaluation metric. In Figure 10 shows the resulting binary mask obtained with the Surface Dice Coefficient and the Volumetric Dice Coefficient for the Core and Whole Tumor. In the case of the Whole Tumor the loss obtained with Surface Dice Coefficient is smaller than the one obtained with the Jaccard Index and the Volumetric Dice Coefficient. For the Core Tumor the loss obtained with the Surface Dice Coefficient is similar to the one reported by the Jaccard Index and both are smaller than the obtained with the Volumetric Dice Coefficient. For the Whole Tumor we can see that the resulting parameters are similar for all the optimization metrics. Nevertheless, for the Core Tumor we can notice that the parameters obtained with the Surface Dice Coefficient differ compared to the ones obtained with the Jaccard Index and the Volumetric Dice Coefficient. In particular, we obtained a smaller value for $\sigma^{2}$ and slightly bigger value for $\Delta_{1}$. This indicates that a smaller value for the diffusion of the particles is compensated by allowing the particles to aggregate with others that are slightly more separated than in the case of the Volumetric Dice and Jaccard Index. Given that both the Volumetric Dice Coefficient and the Jaccard Index are a measure of the superposition between two volumes (in this case two surfaces) they do not represent the proximity between two surfaces making the Surface Dice Coefficient more suitable to use as a loss metric when comparing two different surfaces.

For the $F_{\beta}-$measure we can see in Figure 11 the binary masks obtained for different values of $\beta$ for the Core and Whole Tumor. For the case of the Core Tumor we can observe that for $\beta=0.25$ we obtain areas of misclassified pixels in the tumor region. This can also be seen from Table 4 where the number of False Negatives is bigger and the number of False Positives is smaller compared to the results obtained for bigger values of $\beta$. If we recall (31) we can see that for low values of $\beta$ the False Negatives are multiplied by a factor of $\beta^{2}$ thus having a smaller weight compared to the False Positives. As we increase the value of $\beta$ we can notice from both Table 4 and Figure 11 that modifying the value of $\beta$ has no impact on the resulting binary mask. This also holds true for the Whole Tumor as no difference can be noticed in the results obtained for different values of $\beta$. Finally in Figure 12 we see the loss reported for different values of $\beta$ where the loss equal to $1$ represents a perfect overlap. First, it can be noticed that we obtain the higher value of the loss for $\beta=0.25$, meaning that this should be the most accurate result, which is anyway balanced by the fact that we obtain the larger number of False Negatives. Again, we saw that this can be obtained from (31), where low values of $\beta$ reduce the impact of a large number of False Negatives on the resulting loss. Secondly, we observe that the loss decreases for larger values of $\beta$. This behavior arises because the loss is inversely proportional to $\beta$, while the resulting segmentation masks remain unchanged, as shown in Table 4. This shows that the $F_{\beta}$-measure may not be a reliable metric for these types of segmentation masks and this segmentation method, and that modifying the value of $\beta$ provides no advantage.