Majority voting and mask averaging are simple and easy-to-understand mechanisms to choose a consensus. Yet they suffer from the fact that this decision is purely local without any influence from the neighboring pixels. This can lead to situations where the hard consensus includes some isolated voxels or has very irregular boundaries. This is especially true for mask averaging, which does not have any mechanisms to enforce inter-rater consistency and that relies on the implicit assumption that the neighboring voxels of a segmented voxel are likely to be segmented, which is not the case on the boundaries. Another limitation of majority voting is the case where the number of raters $K$ is even and therefore many decisions are ambiguous with as many foreground than background voxels. Finally, those simple models assume that all raters’ contributions to the consensus are equal which may not be the case. In particular, an underperforming rater will bias the soft consensus with mask averaging.

The consensus prior probability is here supposed to factorize as the product of voxel priors $w_{n}$ values $p(T)=\prod_{n=1}^{N}P(T_{n})=\prod_{n=1}^{N}w_{n}$. The original STAPLE paper (Warfield et al., 2004) also introduced an Ising Markov random field model as a prior consensus probability to enforce that a voxel prior value depends on that of its neighbors. However, this approach leads to solving iteratively graph cuts problems and is not available in most widely used STAPLE implementations. Instead, the original paper assumes simple independent priors that lead to closed-form updates. Choosing $w_{n}=w=\frac{1}{2}$ is a non-informative prior but another common choice is to have a spatially uniform value $w_{n}=w=\frac{1}{NK}\sum_{n,k}S_{n}^{k}$ which is the average relative size of the foreground object in the observed segmentation masks. We further consider more general priors of the form $w=\frac{A}{N^{\alpha}}$, with $A$ a constant independent of the image size, and $\alpha\in\mathbb{N}$ an exponent. The non-informative case $w_{n}=0.5$ corresponds to $\alpha=0$ while the average object size to $\alpha=1$.

The likelihood of the observed data simply writes as $\mathcal{L}(T,\theta)=\prod_{k=1}^{K}p_{k}^{\mathrm{TP}_{k}}(1-p_{k})^{\mathrm{FN}_{k}}q_{k}^{\mathrm{TN}_{k}}(1-q_{k})^{\mathrm{FP}_{k}}$ and does not involve the prior on the consensus. There is no closed-form expression for the estimation of the rater parameters ($p_{k},q_{k}$) and the hard consensus ($T$) maximizing the likelihood. But an iterative maximization of the likelihood is possible by setting its derivatives to zero which leads to the update equation :

The marginal likelihood or evidence writes as $p(\mathcal{S}|\theta)=\prod_{n=1}^{N}(w_{n}\prod_{k}p_{k}^{S_{n}^{k}}(1-p_{k})^{1-S_{n}^{k}}+(1-w_{n})\prod_{k}q_{k}^{1-S_{n}^{k}}(1-q_{k})^{S_{n}^{k}})$ and is only a function of the rater parameters $\theta_{k}$. Its maximization is not tractable in closed form but the expectation-maximization algorithm provides a way to estimate some local maxima. The E-step consists in evaluating the posterior probability from Bayes law with the current estimated sensitivities and specificities :

.

The M-step updates the parameters $p_{k}$ and $q_{k}$ as follows :

where $\mathrm{sTP}_{k}$, $\mathrm{sTN}_{k}$, $\mathrm{sFP}_{k}$, $\mathrm{sFN}_{k}$ are the ”soft extension” of the number of true positive, true negative, false positive, and false negative voxels from rater $k$.

We formulate the estimation of a hard consensus $T$ as the minimization of the sum of the square distance between the consensus $T$ and each observed binary mask $S^{k}$ :

where $d(T,S^{k})$ is a distance as defined in Deza and Deza (2016) between the two masks $S^{k}$ and $T$. This is equivalent to estimating the consensus as a maximum likelihood where the likelihood can be written as $p(S^{k}|T)\propto\exp(-\lambda d(T,S^{k})^{2})$. We can note that the squared sum $\sum_{k=1}^{K}d(M,S^{k})^{2}$ also corresponds to the definition of a Fréchet variance. Based on this interpretation, $T$ appears as the Fréchet mean of $\mathcal{S}$ i.e. its centroid in the metric space defined by $d$.

In section 2.2.2, we have seen that when the sensitivity and specificity are equal, the maximization of the STAPLE model leads to the majority voting algorithm. In this case, we can write the likelihood $p(S^{k}|T)=\gamma_{k}^{\mathrm{TP}_{k}+\mathrm{TN}_{k}}(1-\gamma_{k})^{\mathrm{FP}_{k}+\mathrm{FN}_{k}}$ (where $\gamma_{k}$ is the accuracy parameter) which is a product of $N$ independent Bernoulli distributions. Since the Bernoulli distribution is a member of the exponential family (Dai et al., 2013), it can be also written as $p(S^{k}|T)\propto\exp(-\lambda_{k}(\mathrm{FP}_{k}+\mathrm{FN}_{k}))$ where $\lambda_{k}=\operatorname*{logit}(\gamma_{k})$. The number of false positives or false negatives $\mathrm{FP}_{k}+\mathrm{FN}_{k}$ is the number of elements of symmetric difference between the two sets $S^{k}$ and $T$ : $\mathrm{FP}_{k}+\mathrm{FN}_{k}=|T\Delta S^{k}|=|(T\cup S^{k})\setminus(T\cap S^{k})|$ and is also called the Hamming distance in information theory. Thus, by choosing $d(T,S^{k})=\sqrt{|T\Delta S^{k}|}$, the maximum likelihood leads to majority voting consensus (as detailed in Appendix B).

On the baseline models, soft consensuses were obtained as posterior probabilities of having a consensus from the observed binary masks. However, from the likelihoods $p(S^{k}|\tilde{T})\propto\exp(-\lambda d(\tilde{T},S^{k})^{2})$, the computation of the posterior $p(\tilde{T}|\mathcal{S})$ may not be tractable due to the difficulty of computing the normalization constant. Instead, we propose to approximate $p(\tilde{T}_{n}|\mathcal{S})$ by the quantity $\tilde{U}_{n}\in[0,1]$ such that $\tilde{U}\in[0,1]^{N}$ minimizes the quantity :

where $d^{s}(\tilde{X},S^{k})$ is a distance between the probabilistic array $\tilde{X}$ and the binary mask $S^{k}$. More precisely, the distances $d^{s}(\tilde{X},S^{k})$ considered are soft surrogate of the distance between binary sets $d(\tilde{X},S^{k})$ such that $d^{s}(\tilde{X},S^{k})^{2}=d(\tilde{X},S^{k})^{2}$ when $\tilde{X}\in\{0,1\}^{N}$. For instance, the distance $d(\tilde{X},S^{k})=\|\tilde{X}-S^{k}\|$ is a soft surrogate of the Hamming distance since $|\tilde{X}\Delta S^{k}|=\|\tilde{X}-S^{k}\|^{2}$. Besides it is clear that the mask averaging (MA) method is a soft consensus minimizing the following squared sum $\sum_{k=1}^{K}\|\tilde{U}-S^{k}\|^{2}$.

The estimation of the soft and hard consensus is independent of the background size if the distance $d(T,S^{k})$ is invariant to the number of true negatives. Besides, unlike the MV and MA algorithms, the optimization cannot be performed at the voxel level when the distance cannot be split voxelwise. Instead of optimizing the whole foreground object, we chose to consider each connected component separately from each other to obtain more coherent results. Finally, we further split the optimization into subcrowns with various heuristics to speed up the computation.

The Jaccard coefficient (aka IoU) between binary masks $A$ and $B\in\{0,1\}^{N}$ is defined as : $\text{Jac}(A,B)=\frac{|A\cap B|}{|A\cup B|}$. In Kosub (2019), it is shown that its complementary to 1 $\mathrm{dist}_{J}(A,B)=1-\text{Jac}(A,B)=\frac{|A\Delta B|}{|A\cup B|}$ is a metric between binary sets following the triangular inequality. Several formulations of soft surrogates exist that extend the Jaccard distance. We focused specifically on two of them : the Soergel metric (Späth, 1981; Deza and Deza, 2016) ${d}_{\text{Sg}}(x,y)=\frac{\sum_{i}\max(x_{i},y_{i})-\min(x_{i},y_{i})}{\sum_{i}\max(x_{i},y_{i})}$ which follows the triangular inequality but is not differentiable, and the widely-used Tanimoto distance (Willett et al., 1998; Deza and Deza, 2016; Leach and Gillet, 2007) ${d}_{\text{Tan}}(x,y)=1-\frac{\sum_{i}x_{i}y_{i}}{\sum_{i}x_{i}^{2}+y_{i}^{2}-x_{i}y_{i}}=\frac{||x-y||^{2}}{||x-y||^{2}+<x,y>}$.

It is defined as $\text{DSC}(A,B)=\frac{2|A\cap B|}{|A|+|B|}$ and is widely used in image segmentation as a performance index. Indeed, the Dice index is equal to the F1-score and corresponds to the harmonic mean of the sensitivity and positive predictive value. It is closely related to the Jaccard coefficient as $\text{DSC}(A,B)=\frac{2\text{Jac}(A,B)}{1+\text{Jac}(A,B)}$. The Dice distance $\mathrm{dist}_{D}(A,B)=1-\text{DSC}(A,B)$ is a near-metric i.e. it respects a relaxed form of the triangular inequality (Gragera and Suppakitpaisarn, 2018). Soft surrogates of the Dice distance have been developed especially as a loss function in deep learning. We consider in the remainder two main extensions of the Dice distance (Ma et al., 2021) on non-binary sets defined as ${d}_{\text{pSD}}(x,y)=1-\frac{2\sum_{i}x_{i}y_{i}}{\sum_{i}x_{i}^{p}+\sum_{i}y_{i}^{p}}$ where $p\in\{1,2\}$.

By construction, all those distances only depend on segmented pixels and are independent of the background size. Note that both distances are extended to get a null distance between two empty sets. Using those distances in the Fréchet variance computation, the inclusion of voxels segmented by a large number of raters (resp. a few raters) decreases (resp. increases) its value. The different formulations of the MACCHIatO framework are summarized in table 1.

Since the distances listed in the previous section are independent of the number of true negatives, their computations can be restricted to the union of all rater masks : $\mathcal{E}_{\mathcal{S}}=\{n|\sum_{k=1}^{K}S_{n}^{k}>0\}$. Furthermore, we consider that to decide whether a voxel belongs to the consensus, one should only take into account the regional context associated with the connected components surrounding that voxel, since far-away components may not be relevant. Therefore, we choose to minimize separately the Fréchet variances of Eqs. 6 and 7 for each connected component $St$ of the masks union $\mathcal{E}_{\mathcal{S}}$. Therefore, in practice, we minimize the Local Mean Squared Distance between $\mathcal{S}$ and the consensus : $\mathrm{LMSD}_{d}(\mathcal{S},M)=\sum_{{St}\subset\mathcal{E}_{\mathcal{S}}}\frac{1}{K}\sum_{k}d(S^{k}_{\|St},M_{\|St})^{2}$ where $S^{k}_{\|St}$ (resp. $M_{\|St}$) are the restriction of the binary masks $S^{k}$ (resp. $M$) to the connected component $St$. A benefit of this choice is that the determination of the Fréchet Mean behaves similarly to a structure-wise MV, as the Fréchet Mean of components segmented by less than half of the raters is the null set. However, contrary to MV, raters who do not segment a component kept by the majority of raters do not bias its consensus segmentation, as their contribution to the associated LMSD is $\frac{1}{K}\delta_{\emptyset}(M_{\|St})$ = 0 and does not impact the Fréchet mean. To lighten notations, we drop the $St$ index in the remainder. It is equivalent to considering that $\mathcal{E}_{\mathcal{S}}$ has only one single connected component.

The minimization of the Fréchet variance is a combinatorial problem with a complexity of $2^{|\mathcal{E}_{\mathcal{S}}|}$ for the naive approach. Furthermore, it may lead to several global minima when the number of raters $K$ is small. For those reasons, we propose instead to seek a local minimum of the Fréchet variance by introducing some heuristics in the optimization. With this approach, the local minimum has a lower complexity to compute and, by construction, is maximally connected to avoid isolated voxels.

More precisely, instead of a computationally expensive per voxel minimization of the Fréchet variance, we decompose the set $\mathcal{E}_{\mathcal{S}}$ into a set of subcrowns that take into account the global morphological relationships between each rater mask. The formal definition of subcrowns requires the specification of distance maps $Dm_{\mathcal{N}}(S^{k})$ to each binary mask $S^{k}$ on $\mathcal{E}_{\mathcal{S}}$ according to a chosen neighborhood $\mathcal{N}$. This one can be either the 4 or 8 (resp. 6 or 26) connectivity in 2D (resp. 3D). The distance $Dm_{\mathcal{N}}(S^{k})$ is set to 0 for all voxels inside the object $S^{k}$.

The global morphological distance map is the sum of those distance maps

for all raters on $\mathcal{E}_{\mathcal{S}}$. A crown $C_{td}^{\mathcal{N}}$ is defined as the set of voxels having the same global morphological distance $td$. Those crowns realize a partition of $\mathcal{E}_{\mathcal{S}}$ ($\mathcal{E}_{\mathcal{S}}=\coprod_{td}C_{td}^{\mathcal{N}}$), and the $0$-crown corresponds by construction to the intersection of all masks in $\mathcal{S}$.

We further split each crown as a set of subcrowns by grouping the voxels that have been produced by the same set of raters. In other words, a subcrown corresponds to a set of voxels located at the same morphological distance from the intersection of all rater masks and which have been segmented by the same group of raters, as seen in Fig. 2(a). Formally, a subcrown is noted ${(C_{td}^{\mathcal{N}})}^{g}$ where the superscript $g$ corresponds to a group of raters and subcrowns realize a partition of a crown :

where $\mathcal{P}(\llbracket 1,K\rrbracket)$ is the power set (i.e. the set of all subsets) of the first K integers.

The process for the construction of subcrowns is illustrated in Fig. 2(a)