Multi-dimensional topological loss for cortical plate segmentation in fetal brain MRI

During in utero development, the human fetal cortical plate (CP) that is the future cortex undergoes drastic changes [Tierney and Nelson, 2009]. Indeed, the brain moves from a smooth surface at 10 weeks of gestation to a convoluted one at 35 weeks thanks to the appearance of the cortical gyrifications (see Figure 1). Jointly, the surface area and the volume of the future cortex are respectively increased 50 and 40 times during the 2nd and 3rd trimesters of pregnancy [Vasung et al., 2016, Vasung et al., 2020].

Cortical folding consists in the chronologic appearance of primary, secondary and lastly tertiary sulci [Garcia et al., 2018]. Nearly all gyri are in place at birth, even though the gyrification process carries on afterward [Lenroot and Giedd, 2006]. Cortical gyrification is considered to be a relevant marker of fetal brain maturation, as the chronological sequence of appearance of sulci is well known during the fetal period [Garel et al., 2003]. Conversely, abnormal fetal sulcation and gyration indicate disruption of one of the three main fetal stages of normal cortical development (i.e. cell proliferation, migration and cortical organization) [Barkovich et al., 2005, Leibovitz et al., 2022].

Altered cerebral cortex formation, induced either by genetic mutations, vascular injuries, metabolic diseases, fetal infection or teratogenic causes, may lead to malformations of cortical development. An updated classification of this group of heterogenic disorders has been recently published in a consensus statement [Severino et al., 2020]. Those rare disorders usually manifest with developmental delay, seizures, and motor and sensory deficits [Leventer et al., 2008]. Given the consequences of abnormal brain gyrification, early diagnosis is crucial, for which the analysis of cortical maturation is an asset.

Advanced neuroimaging techniques can greatly benefit the characterization of neurotypical and pathological fetal brain development. Ultrasound (US) is the first prenatal imaging modality for fetal screening. Unfortunately, its sensitivity to the surrounding maternal tissues easily alters the image quality. Conversely, magnetic resonance imaging (MRI) does not share these limitations and hence comes as a key complementary imaging tool to look for additional information in equivocal situation or to confirm or rule out eventual US pathological findings [Salomon et al., 2006, Griffiths et al., 2010, Prayer et al., 2017, Griffiths et al., 2019]. MRI is a non-invasive and reliable imaging method for the monitoring and follow-up of fetal brain development that relies on the tissue properties to generate image contrasts [Gholipour et al., 2014]. MRI is a valuable technique highly suitable for evaluating fetal brain morphometry and connectivity in neurotypical and pathological cases [Garel, 2004, Prayer et al., 2017]. In fetal brain MRI, T2-weighted (T2w) sequences evidence soft-tissue contrast. In this respect, they are essentially used to assess morphological and structural development [Gholipour et al., 2014].

In clinical practice, fast MR sequences are run in order to “freeze” the unpredictable fetal motion in the plane of acquisition [Prayer et al., 2017]. While the series presents an excellent in-plane resolution (1 to 2 $mm^{2}$), slice thickness has a sub-optimal resolution (3 to 5 mm) [Gholipour et al., 2014]. In a clinical examination, several low-resolution (LR) series are acquired in three orthogonal planes to offer complementary structural information from the three dimensions. However, the strong anisotropy of the volumes and the remaining inter-slice motion hampers a correct realignment of the images in the anatomical planes and corrupts any three-dimensional based measurements. Over the last decades, novel advanced image processing algorithms, based on super-resolution (SR) methods, were developed for the reconstruction of 3D volumes [Rousseau et al., 2006, Gholipour et al., 2010, Kuklisova-Murgasova et al., 2012, Tourbier et al., 2015, Ebner et al., 2020]. The underlying reconstruction principle is to combine a set of clinical LR acquisitions into a single high-resolution isotropic motion-free volume [Uus et al., 2022]. In addition to improved visualization, the availability of SR reconstructed volumes opens up to more accurate 3D-based quantitative analysis of the fetal brain anatomy such fetal brain biometry [Pier et al., 2016, Kyriakopoulou et al., 2017, Khawam et al., 2021].

Quantitative analysis of imaging biomarkers has defined cortical development for the typically developing brain [Rajagopalan et al., 2011, Clouchoux et al., 2012, Wright et al., 2014, Xia et al., 2019], while other works evidenced discrepancies in cortical volumes and sulcal patterns in the pathological brain [Clouchoux et al., 2013, Egaña-Ugrinovic et al., 2013, Im et al., 2013, Tarui et al., 2018]. Nevertheless, these analyses require prior additional segmentation processing steps. While expert manual image annotation is considered the gold standard, it is time-consuming, tedious, and subject to inter- and intra- expert variability. For this reason, manual segmentation is not an enduringly reliable method.

Thus robust automated fetal brain MRI segmentation methods are key for further analysis. Nonetheless, in contrast to adult brain segmentation, fetal brain segmentation remains challenging as to provide an age-invariant method [Makropoulos et al., 2018].

Segmentation of the cortical plate is particularly challenging as it undergoes significant changes throughout gestation due to brain growth and maturation, respectively modifying the morphology and the image contrast (see Figure 1). Furthermore, the cortex being a thin structure (from 1 to 2 mm for fetuses between 21 and 40 weeks of gestational age [Vasung et al., 2016]) strongly affected by partial volume effects, anatomical topology is prone to be incorrectly captured by automatic segmentation methods.

In Caldairou et al. [Caldairou et al., 2011], they introduced the first topological-based segmentation of the fetal cortex based on geometrical constraints along with anatomical and topological priors. However, the sample size in that study was small (i.e. six fetuses) and topological correctness was not evaluated for the provided segmentations. More recently, deep learning (DL) methods have also focused on fetal brain MRI cortical gray matter segmentation. Based on a neonatal segmentation framework, a recent study introduced a hybrid segmentation process that minimizes the need for human interaction in the segmentation of the developing cortex [Fetit et al., 2020]. Also with DL models, a multiple-predictions approach with a test-time augmentation to improve the robustness of the method was suggested [Hong et al., 2020]. Finally, a two-stage segmentation framework with an attention refinement was proposed [Dou et al., 2021]. Nevertheless, while the segmentation accuracy of these recent DL-based methods is promising, none of these CP segmentation works includes topological constraints nor assesses the topology. Overall, these studies of automatic segmentation frameworks report high similarity, as for overlap (Dice similarity coefficient of 0.87 in [Dou et al., 2021] and 0.90 in [Hong et al., 2020]), and low difference, as for the boundary distance-based metrics (average symmetric surface distance of 0.28 $mm$ in [Dou et al., 2021] and mean surface distance of 0.18 $mm$ in [Hong et al., 2020]), compared to the ground truth segmentation, but illustrated results show a lack of topological consistency with notably discontinuous/broken cortical ribbons (see Fig. 5 in  [Dou et al., 2021], Figures 5 and 6 in [Fetit et al., 2020]).

Here, we propose to integrate a topological constraint in a deep image segmentation framework to overcome the limitation of disjoint CP segmentation in fetal MRI. To our knowledge, only two works previously explored the topological fidelity of the semantic medical image segmentation with DL. First, Hu et al. [Hu et al., 2019] proposed using a topological loss for neuronal membrane segmentation. Second, a study presented topological constraints for MR cardiac image segmentation [Clough et al., 2021] based on prior topological knowledge, such as the number of connected components or handles present in the structure of interest. Although theoretical CP topological features are known, such prior information could only be valid in a whole 3D volume analysis. Therefore, we inspired from [Hu et al., 2019] that does not share this limitation to build a prior-free framework.

In this work, we incorporate topological constraints and assess the topology of the CP segmentation in fetal brain MRI. To this end, our first contribution is the generalization of a topological loss function based on persistent homology into a multi-dimensional (in topological spaces) formulation. Using a publicly available dataset of SR reconstructed fetal brain MRI [Payette et al., 2021] along with a spatiotemporal atlas [Gholipour et al., 2017], we compare our topology-integrative optimization method (TopoCP) to other widely used loss functions (Baseline and Hybrid) [Payette et al., 2022] with a state-of-the-art U-Net segmentation network. We further assess our proposed method compared to semi-automatic manual annotations. For the first time, the topological correctness of the segmentation is assessed in the evaluation of automatic CP segmentation. In that respect, our second contribution is a new topology-based metric for the quantification of the CP segmentation defects not only in number but also in size. We quantitatively compare automatic and semi-automatic segmentation with complementary metrics on two independent pure testing sets. As a complement, three fetal brain MRI experts further visually compare automatic segmentation on an additional out-of-distribution clinical dataset. The results evidence an overall significant improvement in the segmentation using our proposed topological loss function.

In Section 2 we introduce our multi-dimensional topological loss function, the overall segmentation framework and our original topology hole ratio assessment metrics; in Section 3, we describe the experimental design, along with materials and description of the experiments performed; in Section 4, we describe and discuss the results; and finally, we conclude on this work in Section 5.

In semantic image segmentation, conventional optimization loss functions (e.g. the cross-entropy) often proceed to a pixel-wise comparison of the class-prediction that is summarized in a likelihood function $f$ to the one-hot encoded target vector. In this work, we aim to integrate the analysis of the global shape correctness of the prediction through the study of topology during the model optimization.

Our topological loss is an architecture-agnostic optimization function. In other words, it is independent of the deep-learning network architecture. In our segmentation framework, we use a state-of-the-art architecture, U-Net [Ronneberger et al., 2015], to compare different optimization methods. Two reference loss functions (Baseline and Hybrid) are implemented to evaluate the added value of our topological loss function (TopoCP) (see configurations in Section 2.2.3).

Recent works [Yeghiazaryan and Voiculescu, 2018, Maier-Hein et al., 2022] evidenced the importance of considering complementary metrics for the assessment of semantic segmentation methods. Specifically, in this work, segmentation should be assessed for its closeness to the target topology.

We presented in Section 2.1.1 the BN$k$ that define global features of a 3D binary image. To quantitatively compare topology of binary images, we introduce the $k$-dimensional Betti number error (BNE$k$) as the absolute difference of the ground truth expected value and the prediction measures. As the CP segmentation is filtered for its biggest connected component (see Section 2.2.4), BNE$0$ is incidentally irrelevant to consider. Additionally, this work specifically focuses on the presence of 1-dimensional holes. In this way, we only consider BNE$1$. Besides, considering that the human cortex is a closed structure with no obstruction, its ground truth expected BN$1$ is 0.

While BNE$1$ focuses on the count of occlusions on the CP surface, this score is not providing any information on the holes themselves. Indeed, we often observe that the number of holes that is a discrete value is not necessarily correlated to the size of the broken connections (i.e. the obstruction) nor to the size of the structure of interest (i.e. the fetal CP).

In that respect, we introduce a new metric that aims to quantify the size of the hole. The hole ratio (HR) is the ratio of false negative voxels that are connected to a hole ($FN_{holes}$) over the true voxels of the region of interest, which are represented by the true positives ($TP$) and false negatives ($FN$).

Our in-house implementation (available online, see Section 2.4) illustrated in Figure 5 successively identifies the location of one voxel per 1-dimensional hole, propagates these voxels into the mask of holes candidates, i.e. the FN, and finally, computes the volume ratio presented in Equation 5. Let us note that this measure strongly relies on the topological correctness of the ground truth.

Figure 6(a) shows (B) and (C), two 3D rendering of the same image patch CP segmentation. (A) shows the 3D rendering of the ground truth, enlighting the region of interest. One can easily observe that the (B) segmentation has a main hole compared to the (C) segmentation that presents multiple medium-size holes. Additionally, quantitative results (see Table 6(b)) confirm the discrepancy between the quantity of holes as a number and the quantity of holes as a ratio of the region of interest.

Albeit BNE1 must be used with caution, it is still a relevant score to assess the CP topology in the absence of topologically accurate ground truth segmentation. Nonetheless, we promote the use of an additional quantitative metric relative to the size of the occlusions to undertake a robust quantitative analysis of the broken connections in the CP segmentation.

Both topology-based measures (BNE1 and HR) rely on the cubical complex implementation of the GUDHI library (v3.5.0) [noa,].

Baseline, Hybrid and TopoCP models implementation, including the optimization loss functions can be found in our Github repository ²²2https://github.com/Medical-Image-Analysis-Laboratory/FetalCP˙segmentation. The weights of the trained model are made available to perform inference. The implementation of the topology-based evaluation metric is also made available at this link.

A summary of clinical and atlas datasets is shown in Table 1.

Although our segmentation framework infers CP segmentation in a 2.5D multi-view strategy (2D image patches from the three orthogonal planes), we proceed to the quantitative evaluation in 3D that is of the whole cortical volume. Automatic medical image segmentation requires the use of complementary metrics for the assessment of different aspects of the segmentation [Yeghiazaryan and Voiculescu, 2018, Maier-Hein et al., 2022]. Most commonly, overlap-based (e.g. Dice similarity coefficient, Jaccard similarity coefficient, Intersection over union) and distance-based (e.g. X^th percentile Hausdorff distance, average symmetric surface distance) metrics are reported [Yeghiazaryan and Voiculescu, 2018]. However, in this work we aim to assess the segmentation not only in terms of overlap and distance accuracy to the ground truth but also in terms of shape correctness. Therefore, we also consider topology-based metrics. Table 2 summarizes the metrics used in the training (for learning monitoring) and testing (for final evaluation) phases.

The Dice Similarity Coefficient (DSC) [Dice, 1945] is an overlap-based similarity metric. Robust to outliers, it is a widely-used metric to assess medical image segmentation accuracy. Average Symmetric Surface Distance (ASSD) is the mean of the directed average surface distances [Yeghiazaryan and Voiculescu, 2018]. The latter is defined as the average of the distances of points from one surface to their closest points on the other object boundary. The ASSD is computed using the python MedPy ³³3https://loli.github.io/medpy/ implementation (v0.4.0).

In the case of absence of topologically accurate ground truth (i.e. see Section 3.3.1 for cross-validation details), we consider BNE1 to quantitatively assess the topology, while our proposed hole ratio HR is used in the pure testing phase (see details on topology metrics in Section 2.3).

Arrows in Table 2 indicate whether each metric is better maximized or minimized. Taking values between 0 and 1, DSC is a similarity metric that is better maximized ($\uparrow$). Difference metrics (ASSD, BNE1 and HR) must be minimized ($\downarrow$).

Table 3 shows the averaged validation scores of all three configurations, and specifically for each $\lambda_{topo}$ assessed in the TopoCP configuration. Our first observation is that, regardless of the value of the $\lambda_{topo}$ parameter, TopoCP is better performing than both reference methods, as we reach the state of the art performances in all the three complementary metrics. Overall, all $TopoCP_{n}$ give similar DSC (mean: 0.76) and ASSD (mean $\pm$ standard deviation: 0.27 $\pm$ 0.01) performances, although $\lambda_{topo}=0.01$ is of the highest rank for both overlap and boundary-distance based scores. An increased inter-$TopoCP_{n}$ variability appears for the topology-based metric (BNE$1$) with mean score from 20.6 to 22.5. Counting the number of bores in the CP segmentation, $\lambda_{topo}=0.005$ is giving the best performances. We observe large BNE$1$ standard deviation for all $TopoCP_{n}$. Nonetheless, the finest topology-relative $\lambda_{topo}$ is not only giving the minimum averaged BNE$1$, but is also noticeably presenting a smaller BNE$1$ standard deviation of 7.8 (BNE$1$ range: from 7.8 to 10). Therefore, $\lambda_{topo}=0.005$ is the most accurate and precise of the $\lambda_{topo}$ assessed as for the topology fidelity. The substantial fluctuation in the topological metric shows the importance of the choice of the $\lambda_{topo}$ hyper-parameter.

Finally, our global ranking that is derived from metric-wise rankings evidences the ideal value $\lambda_{topo}=0.005$. We observe that none of the extreme values considered (i.e. $0.0002$ and $0.2$) are in the Top 3 best performing $\lambda_{topo}$. Therefore, we can say that although $\lambda_{topo}=0.005$ might not be the exact optimal $\lambda_{topo}$, it certainly falls in a relevant range and in the right order of magnitude.

Figure 10 illustrates the accuracy of the fetal CP segmentation for a pathological subject of 26.6 (Subject 1) and a neurotypical subject of 34.8 (Subject 2) weeks of GA. The topologically corrected ground truth and the three configurations segmentation with Baseline, Hybrid, and TopoCP, are compared. Qualitative 2D assessment (top rows) of the segmentation are presented as an overlay on the T2w image on an axial, respectively coronal, view for Subject 1, respectively Subject 2. Additionally, 3D rendering of the CP segmentation are presented in the bottom rows. Overall, we observe that all configurations generates a thinner ribbon than the corrected ground truth. Specifically, TopoCP presents fixed cortical connections that are broken in the Baseline and Hybrid segmentations (white arrows). The CP TopoCP segmentation 3D rendering seems to present less bores than the two reference configurations Baseline and Hybrid. In particular, the segmentation appears, equivalently for the young and the old fetuses, more challenging, for all methods, in the lower parts of the frontal and temporal lobes, although TopoCP seems to exhibit a more sensitive segmentation in these areas. TopoCP appears to be more sensitive to the complexity of the CP morphology in older fetuses. White circles evidence in Subject 2 an improved segmentation in the hippocampal area and the depth of a gyrification. White arrows show area where the topological correctness recovered with TopoCP, compared to the Baseline. In Subject 1, two connections are fixed in the frontal lobe, although one of them is already fixed in the Hybrid configuration.

Quantitative results for both test sets (total of 27 cases), FeTA and STA, are respectively presented in Table 4(a) and Table 4(b). Tables show the mean ± standard deviation of the CP segmentation for each testing metric (DSC, ASSD and HR) in each configuration. Overall, the performance of the segmentation framework is improved when trained with our optimal TopoCP configuration (i.e. $\lambda_{topo}=0.005$) for all metrics in both datasets. TopoCP is always performing significantly better than both Baseline and Hybrid methods.

Furthermore, while all analyzed aspects of the segmentation, namely the overlap, the boundary-distance and the topology are improved with TopoCP compared to the reference methods, we observe a drop in the performances between FeTA and STA evaluation. Indeed, DSC goes from 0.85 in FeTA to 0.79 in STA, ASSD from 0.19 to 0.41 and HR from 0.06 to 0.23. We believe this is due to the domain shift between FeTA and STA images (different reconstruction pipelines, different intensity-based processing, etc) generating an inter-dataset variation in the data distribution. Such domain gap between the two sets of images is not learned from the training data that are only composed of FeTA images. Still, DL can generalize to some extent. Further training with multi-dataset images or the use of domain adaptation strategies can partially fit the domain gap. Let us note that we do not address the data distribution generalization in this paper as it is beyond its scope. Furthermore, such performance drop occurs also in Baseline and Hybrid approach.

Taking advantage of the steady GA-distribution in the STA set, we perform an analysis of the metrics throughout gestation. Figure 11 shows the performance metrics as a function of the GA for 18 cases from 21 to 37 weeks of GA. Regardless of the configuration, we observe a trend in the performances based on the GA. Indeed, all metrics reach better performances for subjects younger than 30 weeks of GA. From week 23 to 31, all three methods (Basline, Hybrid and TopoCP) seem to give equivalent scores. Outside this range (i.e. GA$<23$ and GA$>30$ weeks), TopoCP is always performing better than the other two methods, except for one outlier subject of 38 weeks of GA.

We visually inspect the TopoCP segmentation mask of the STA 38-weeks-old subject to better understand the origin of this outlier. Specifically in the cerebellum, we observe the presence of false positives that are connected to the main cortical segmentation through brainstem false positives. The cerebellum is a ”little brain” composed of white matter encased in the cerebellar cortex. In terms of fetal brain T2w MR contrast and similarly to the CP, the cerebellar cortex expresses as a thin dark ribbon surrounding white matter. Therefore, it is a challenging area to accurately differentiate in the segmentation of the fetal CP at a patch-level. Similar mis-segmentations appear in younger fetuses, nevertheless our post-processing step to keep the biggest connected component filters out most of it. In this oldest subject, cerebellum errors are worsened with brainstem mis-segmentation. We hypothesize such false positive errors in the cerebellum are due to mis-leading contextual information, due to the reduced field of view of the patches. Therefore, we believe that increasing the patch size could help to overcome these mis-segmentation. Nevertheless, while the performances are particularly damaged in the distance metric, TopoCP still performs better in terms of DSC and HR compares to the other configurations.

Overall, it is in the second half of the third trimester GAs (i.e. GA$>$30 weeks) that we observe an increased benefit from TopoCP, compared to other methods. Our topological loss has a stronger positive effect on the topological errors for old subjects with more complex topology, although the whole range of gestational age consistently presents benefits from TopoCP.

We derive two hypothesis on the variation of the performances throughout gestation. First, we recall that the training data present subjects in the range 20.9 to 34.8 weeks of GA with mean 28.2 and standard deviation 3.6. Therefore, variation of young and old fetal brains are less represented in the learning process. Additionally, third trimester subjects present advanced sulcal patterns, resulting in a substantially more complex topology. Therefore, we postulate this accentuate the unstable evolution of segmentation accuracy over gestation.

Figure 12 presents a comparison of HR at a lobe-level between the configurations. This boxplot evidences the significant benefits ($p<0.05$) of TopoCP in most areas (frontal, occipital, temporal and insula lobes) compared to the Baseline configuration. In the parietal lobe, TopoCP is on average performing better than the Baseline although without statistical significance. Compared to the Hybrid configuration, TopoCP presents a significantly lower HR in all brain lobes. Regardless of the configuration, the parietal lobe is always the better segmented lobe in terms of HR as opposed to the insula lobe.

Table 5 summarizes the performance metrics (mean $\pm$ standard deviation) in both neurotypical (NT) and pathological (PT) groups. We observe the same trend in this group-based analysis as in the overall method comparison (see Tables 4 in Section 4.2). Regrettably, precise information on the pathology are not available in the FeTA dataset meta data. Therefore, we cannot draw any conclusion relative to eventual cortical pathologies.

Figure 13 (top) shows a comparative T2w axial view of the ground truth topologically corrected segmentation (A), the original manual annotation provided in FeTA (B) and the TopoCP predicted segmentation (C). Overall, we observe an improved accuracy with the automatic method. Specifically, white arrows indicate cortical location where TopoCP fixes topological inconsistencies, compared to the original manual annotations. The white circle focuses on the hippocampal area where the manual annotations are confused. Figure 13 (bottom) shows 3D rendering of the true cortical volume (green). In the manual and automatic TopoCP segmentation, the 1-dimensional holes are evidenced with the false negatives connected to 1-dimensional holes (light camel).

Table 6 shows a comparison of the performance metrics of the original FeTA manual annotations on 9 subjects and our TopoCP method. TopoCP is significantly better (*) than the original manual annotations in all metrics (DSC, ASSD and HR).

The automatic TopoCP segmentation method is able to learn segmentation features from noisy annotations. This improvement is conveyed in all three similarity, boundary-distance and topology -based metrics.

Table 7 summarizes the classification results of the segmentation masks according to each expert into Worst, Medium or Best. A consensus of the three experts assessment is presented in the bottom row.

The estimated Gwet’s AC between the three experts was $0.68$ (standard error (SE): $0.10$) for Worst, $0.68$ (SE: $0.10$) for Medium and $1$ for Best segmentation classifications. According to Altman’s benchmarking scale, the estimated coefficients for Worst and Medium segmentations are considered to be either Moderate, Good or Very Good with a probability of $0.99$. The Best segmentation classification has a Very good agreement between the experts with a probability of 1. With a percentage agreement of 100 %, the consensus of experts classifies TopoCP as the Best segmentation method in 100% of the cases. Inter-rater discrepancies are observed in the choice of Worst and Medium between Baseline and Hybrid segmentation.

Although TopoCP is ranked as Best segmentation in all cases, predictions still present many segmentation errors. We emphasize that the distribution of this clinical set differs from the FeTA training set as they were generated with different SR methods. Nevertheless, while all three configurations present altered segmentation due to the domain shift, still, TopoCP remains the better performing method.