Deep nested level sets: Fully automated segmentation of cardiac MR images in patients with pulmonary hypertension
Jinming Duan Jo Schlemper Wenjia Bai Timothy J W Dawes Ghalib Bello Georgia Doumou Antonio De Marvao Declan P O'Regan Daniel Rueckert



Deep nested level set method



J Duan, et al.
^1Biomedical Image Analysis Group, Imperial College London, London, UK

^2MRC London Institute of Medical Sciences, Imperial College London, London, UK

j.duan@imperial.ac.uk



Deep nested level sets: Fully automated segmentation of cardiac MR images in patients with pulmonary hypertension
    Jinming Duan^1,2, Jo Schlemper^1, Wenjia Bai^1, Timothy J W Dawes^2, Ghalib Bello^2, Georgia Doumou^2, Antonio De Marvao^2, Declan P O'Regan^2, Daniel Rueckert^1
    December 30, 2023
=====================================================================================================================================================================









In this paper we introduce a novel and accurate optimisation method for segmentation of cardiac MR (CMR) images in patients with pulmonary hypertension (PH). The proposed method explicitly takes into account the image features learned from a deep neural network. To this end, we estimate simultaneous probability maps over region and edge locations in CMR images using a fully convolutional network. Due to the distinct morphology of the heart in patients with PH, these probability maps can then be incorporated in a single nested level set optimisation framework to achieve multi-region segmentation with high efficiency. The proposed method uses an automatic way for level set initialisation and thus the whole optimisation is fully automated. We demonstrate that the proposed deep nested level set (DNLS) method outperforms existing state-of-the-art methods for CMR segmentation in PH patients.







§ INTRODUCTION


Pulmonary hypertension (PH) is a cardiorespiratory syndrome characterised by increased blood pressure in pulmonary arteries. It typically follows a rapidly progressive course. As such, early identification of PH patients with elevated risk of a deteriorating course is of paramount importance. For this, accurate segmentation of different functional regions of the heart in CMR images is critical.  

Numerous methods for automatic and semi-automatic CMR image segmentation have been proposed, including deformable models <cit.>, atlas-based image registration models <cit.> as well as statistical shape and appearance models <cit.>. More recently, deep learning-based methods have achieved state-of-the-art performance in the CMR domain <cit.>. However, the above approaches for CMR image segmentation have multiple drawbacks. First, they tend to focus on left ventricle (LV) <cit.>. However, the prognostic importance of the right ventricle (RV) is a broad range of cardiovascular disease and using the coupled biventricular motion of the heart enables more accurate cardiac assessment. Second, existing approaches rely on manual initialisation of the image segmentation or definition of key anatomical landmarks <cit.>. This becomes less feasible in population-level applications involving hundreds or thousands of CMR images. Third, existing techniques have been mainly developed and validated using normal (healthy) hearts <cit.>. Few studies have focused on abnormal hearts in PH patients.

To address the aforementioned limitations of current approaches, in this paper we propose a deep nested level set (DNLS) method for automated biventricular segmentation of CMR images. More specifically, we make three distinct contributions to the area of CMR segmentation, particularly for PH patients: First, we introduce a deep fully convolutional network that effectively combines two loss functions, i.e. softmax cross-entropy and class-balanced sigmoid cross-entropy. As such, the neural network is able to simultaneously extract robust region and edge features from CMR images. Second, we introduce a novel implicit representation of PH hearts that utilises multiple nested level lines of a continuous level set function. This nested level set representation can be effectively deployed with the learned deep features from the proposed network. Furthermore, an initialisation of the level set function can be readily derived from the learned feature. Therefore, DNLS does not need user intervention (manual initialisation or landmark placement) and is fully automated. Finally, we apply the proposed DNLS method to clinical data acquired from 430 PH patients (approx. 12000 images), and compare its performance with state-of-the-art approaches.











§ MODELLING BIVENTRICULAR ANATOMY IN PATIENTS WITH PH




ro0.4



    < g r a p h i c s >



Short-axis images of a healthy subject (left) and a PH subject (right), including the anatomical explanation of both LV and RV. The desired epicardial contours (red) and endocardial contours (yellow) from both ventricles are plotted.




To illustrate cardiac morphology in patients with PH, Fig <ref> shows the difference in CMR images from a representative healthy subject and a PH subject. In health, the RV is crescentic in short-axis views and triangular in long-axis views, wrapping around the thicker-walled LV. In PH, the initial hypertrophic response of the RV increases contractility but is followed invariably by progressive dilatation and failure heralding clinical deterioration and ultimately death. During this deterioration, the dilated RV pushes onto the LV to deform and lose its roundness. Moreover, in PH the myocardium around RV become much thicker than a healthy one, allowing PH cardiac morphology to be modelled by a nested level set. Next, we incorporate the biventricular anatomy of PH hearts into our model for automated segmentation of LV and RV cavities and myocardium.





§ METHODOLOGY




Nested level set approach: We view image segmentation in PH as a multi-region image segmentation problem. Let I : Ω→ℝ^d denote an input image defined on the domain Ω⊂ℝ^2. We segment the image into a set of n pairwise disjoint region Ω_i, with Ω  =  ∪ _i = 1^nΩ _i, Ω _i∩Ω _j = ∅ ∀ i  j. The segmentation task can be solved by computing a labelling function l(x):Ω→{1,...,n} that indicates which of the n regions each pixel belongs to: Ω _i = {x| l( x ) = i.}. The problem is then formulated as an energy minimisation problem consisting of a data term and a regularisation term

    min_Ω _1,...,Ω _n{∑_i = 1^n ∫_Ω _if_i( x )dx  + λ∑_i = 1^n Per_g( Ω _i,Ω)}.


ro0.4



    < g r a p h i c s >



An example of partitioning the domain Ω into 4 disjoint regions (right), using 3 nested level lines {x|ϕ(x)= c_i , i=1,2,3} of the same function ϕ (left). The intersactions between the 3D smooth surface ϕ and the 2D plans correspond to the three nested curves on the right.





The data term, f_i:Ω→ℝ is associated with region that takes on smaller values if the respective pixel position has stronger response to region. In a Bayesian MAP inference framework, f_i( x ) =  - log P_i ( I( x )|Ω _i) corresponds to the negative logarithm of the conditional probability for a specific pixel color at the given location x within region Ω_i. Here we refer to f_i as region feature. The second term, Per_g( Ω _i,Ω) is the perimeter of the segmentation region Ω_i, weighted by the non-negative function g. This energy term alone is known as geodesic distance, the minimisation over which can be interpreted as finding a geodesic curve in a Riemannian space. The choice of g can be an edge detection function which favours boundaries that have strong gradients of the input image I. Here we refer to g as edge feature. 

We apply the variational level set method <cit.> to (<ref>) in this study. Because a PH heart can be implicitly represented by two nested level lines of a continuous level set function ({x|ϕ(x)=c_i, i=1,2} in Fig <ref>). Note that the nested level set idea present here is inspired from previous work <cit.>. Our approach uses features learned from many images while previous work only consider single image. With the idea, we are able to approximate the multi-region segmentation energy (<ref>) by using only one continuous function. The computational cost is thus small. Now assume that the contours in the image I can be represented by level lines of the same Lipschitz continuous level set function ϕ:Ω→ℝ. With n-1 distinct levels {c_1 < c_2 < ··· <c_n-1}, the implicit function ϕ partitions the domain Ω into n disjoint regions, together with their boundaries (see Fig <ref> right). We can then define the characteristic function χ_iϕ for each region Ω_i as


    χ _iϕ(x)  = {[                    H( c_i - ϕ(x) )                                i=1; H( ϕ(x)  - c_i - 1)H( c_i - ϕ(x) )                        2 ≤ i ≤ n-1;                H( ϕ(x)  - c_i - 1)                                i=n ].,

where H is the one-dimensional Heaviside function that takes on either 0 or 1 over the whole domain Ω. Due to the non-differentiate nature of H it is usually approximated by its smooth version H_ϵ for numerical calculation <cit.>. Note that in (<ref>) ∑_i = 1^n χ _iϕ = 1 is automatically satisfied, meaning that the resulting segmentation will not produce a vacuum or an overlap effect. That is, by using (<ref>) Ω  =  ∪ _i = 1^nΩ _i and Ω _i∩Ω _j = ∅ hold all the time. With the definition of χ _iϕ, we can readily reformulate (<ref>) in the following new energy minimisation problem 

    min_ϕ(x){∑_i = 1^n ∫_Ωf_i( x )χ _iϕ( x )dx  + λ∑_i = 1^n - 1∫_Ωg( x )| ∇ H( ϕ( x ) - c_i )|dx}.

Note that (<ref>) differs from (<ref>) in multiple ways due to the use of the smooth function ϕ and characteristic function (<ref>). First, the variable to be minimised is the n regions Ω _1,...,Ω _n in (<ref>) while the smooth function ϕ in (<ref>). Second, the minimisation domain is changing from over Ω_i in (<ref>) to over Ω in (<ref>). Third, (<ref>) uses an abstract Per_g( Ω _i,Ω) for the weighted length of the boundary between two adjacent regions, while (<ref>) represents the weighted length with the co-area formula, i.e. ∫_Ωg| ∇ H( ϕ - c_i)|dx. Finally, the upper limit of summation in the regularisation term of (<ref>) is n while n-1 in that of (<ref>). So far, the region features f_i and the edge feature g have not been defined. Next, we will tackle this problem.






Learning deep features using fully convolutional network: We propose a deep neural network that can effectively learn region and edge features from many labelled PH CMR images. Learned features are then incorporated to (<ref>). Let us formulate the learning problem as follows: we denote the input training data set by S={(U_p, R_p, E_p), p=1,...,N}, where sample U_p={u^p_j,j=1,...,|U_p|} is the raw input image, R_p={r^p_j,j=1,...,|R_p|}, r^p_j ∈{1,...,n} is the ground truth region labels (n regions) for image U_p, and E_p={e^p_j,j=1,...,|E_p|}, e^p_j ∈{0,1} is the ground truth binary edge map for U_p. We denote all network layer parameters as W and  propose to minimise the following objective function via the (back-propagation) stochastic gradient descent

    W^* = argmin (L_R(W) + α L_E(W) ),

where L_R(W) is the region associated cross-entropy loss that enables the network to learn region features, while L_E(W) is the edge associated cross-entropy loss for learning edge features. The weight α balances the two losses. By minimising (<ref>), the network is able to output joint region and edge probability maps simultaneously. In our image-to-image training, the loss function is computed over all pixels in a training image U={u_j,j=1,...,|U|}, a region map R={r_j,j=1,...,|R|}, r_j ∈{1,...,n} and an edge map E={e_j,j=1,...,|E|}, e_j ∈{0,1}. The definitions of L_R(W) and L_E(W) are given as follows.

    L_R(W) = -∑_jlogP_so(r_j|U,W),
 
where j denotes the pixel index, and P_so(r_j|U,W) is the channel-wise softmax probability provided by the network at pixel j for image U. The edge loss is 

    L_E(W) = -β∑_j ∈ Y_+log P_si(e_j=1|U,W) - (1-β)∑_j ∈ Y_-logP_si(e_j=0|U,W).

For a typical CMR image, the distribution of edge and non-edge pixels is heavily biased. Therefore, we use the strategy  <cit.> to automatically balance edge and non-edge classes. Specifically, we use a class-balancing weight β. Here, β  =  |Y_ -|/|Y| and 1-β=|Y_ +|/|Y|, where |Y_-| and |Y_+| respectively denote edge and non-edge ground truth label pixels. P_si(e_j=1|U,W) is the pixel-wise sigmoid probability provided by the network at non-edge pixel j for image U. 




In Fig <ref>, we show the network architecture for automatic feature extraction, which is a fully convolutional network (FCN) and adapted from the U-net architecture <cit.>. Batch-normalisation (BN) is used after each convolutional layer, and before a rectified linear unit (ReLU) activation. The last layer is however followed by the softmax and sigmoid functions. In the FCN, input images have pixel dimensions of 160 × 160. Every layer whose label is prefixed with `C' performs the operation: convolution → BN → ReLU, except C17. The (filter size/stride) is (3×3/1) for layers from C1 to C16, excluding layers C3, C5, C8 and C11 which are (3×3/2). The arrows represent (3×3/1) convolutional layers (C14a-e) followed by a transpose convolutional (up) layer with a factor necessary to achieve feature map volumes with size 160 × 160 × 32, all of which are concatenated into the red feature map volume. Finally, C17 applies a (1×1/1) convolution with a softmax activation and a sigmoid activation, producing the blue feature map volume with a depth n+1, corresponding to n (3) region features and an edge feature of an image.

After the network is trained, we deploy it on the given image I in the validation set and obtain the joint region and edge probability maps from the last convolutional layer

    (P_R,P_E)=CNN(I,W^*),

where CNN(·) denotes the trained network. P_R is a vector region probability map including n (number of regions) channels, while P_E is a scalar edge probability map. These probability maps are then fed to the energy (<ref>), in which f_i = -logP_Ri, i={1,...,n} and g=P_E. With all necessary elements at hand, we are ready to minimise (<ref>) next. 


Optimisation: The minimisation process of (<ref>) entails the calculus of variations, by which we obtain the resulting Euler-Lagrange (EL) equation with respect to the variable ϕ. A solution (ϕ^*) to the EL equation is then iteratively sought by the following gradient descent method

    ∂ϕ/∂ t =  - ∑_i = 1^n f_i∂χ _iϕ/∂ϕ  + λκ _g∑_i = 1^n - 1δ_ϵ( ϕ  - c_i),

where κ _g = div( g∇ϕ/.-| ∇ϕ|) is the weighted curvature that can be numerically implemented by the finite difference method on a half-point grid <cit.>. δ_ϵ is the derivative of H_ϵ, which is defined in <cit.>. 

At steady state of (<ref>), a local or global minimiser of (<ref>) can be found. Note that the energy (<ref>) is nonconvex so it may have more than one global minimiser. To obtain a desirable segmentation result, we need a close initialisation of the level set function (ϕ^0) such that the algorithm converges to the solution we want. We tackle this problem by thresholding the region probability map P_R3 and then computing the signed distance function (SDF) from the binary image using the fast sweeping algorithm. The resulting SDF is then used as ϕ^0 for (<ref>). In this way, the whole optimisation process is fully automated. 




§ EXPERIMENTAL RESULTS


Data: Experiments were performed using short-axis CMR images from 430 PH patients. For each patient 10 to 16 short-axis slices were acquired roughly covering the whole heart. Each short-axis image has resolution of 1.5 × 1.5 × 8.0 mm^3. Due to the large slice thickness of the short-axis slices and the inter-slice shift caused by respiratory motion, we train the FCN in a 2D fashion and apply the DNLS method to segment each slice separately. The ground truth region labels were generated using a semi-automatic process which included a manual correction step by an experienced clinical expert. Region labels for each subject contain the left and right ventricular blood pools and myocardial walls for all 430 subjects at end-diastolic (ED) and end-systolic (ES) frames. The ground truth edge labels are derived from the region label maps by identifying pixels with label transitions. The dataset was randomly split into training datasets (400 subjects) and validation datasets (30 subjects). For image pre-processing, all training images were reshaped to the same size of 160 × 160 with zero-padding, and image intensity was normalised to the range of [0, 1] before training. 

 Parameters: The following parameters were used for the experiments in this work: First, there are six parameters associated with finding a desirable solution to (<ref>). They are the weighting parameter λ (1), regularisation parameter ϵ (1.5), two levels c_1 (0) and c_2 (8), time step t (0.1), and iteration number (200). Second, for training the network, we use Adam SGD with learning rate (0.001) and batch size (two subjects) for each of 50000 iterations. The weight α in (<ref>) is set to 1. We perform data augmentation on-the-fly, which includes random translation, rotation, scaling and intensity rescaling of the input images and labels at each iteration. In this way, the network is robust against new images as it has seen millions of different inputs by the end of training. Note that data augmentation is crucial to obtain better results. Training took approx. 10 hours (50000 iterations) on a Nvidia Titan XP GPU, while testing took 5s in order to segment all the images for one subject at ED and ES. 




 Comparsion: The segmentation performance was evaluated by computing the Dice overlap metric between the automated and ground truth segmentations for LV & RV cavities and myocardium. We compared our method with the vanilla CNN proposed in <cit.>, the code of which is publicly available. DNLS was also compared with the vanilla CNN with a conditional random field (CRF) <cit.> refinement (CRF-CNN). In Fig <ref>, visual comparison suggests that DNLS provides significant segmentation improvements over CNN and CRF-CNN. For example, at the base of the right ventricle both CNN and CRF-CNN fail to retain the correct anatomical relationship between endocardium and epicardium portraying the endocardial border outside the epicardium. CRF-CNN by contrast retains the endocardial border within the epicardium, as described in the ground truth. In Table <ref>, we report their Dice metric of ED and ES time frames in the validation dataset and show that our DNLS method outperforms the other two methods for all the anatomical structures, especially for the myocardium. CNN is the fastest method as it was deployed with GPU, and DNLS is the most computationally expensive method due to its complex optimisation processes. 









§ CONCLUSION


In this paper, we proposed the deep nested level set (DNLS) approach for segmentation of CMR images in patients with pulmonary hypertension. The main contribution is that we combined the classical level set method with the prevalent fully convolutional network to address the problem of pathological image segmentation, which is a major challenge in medical image segmentation. The DNLS inherits advantages of both level set method and neural network, the former being able to model complex geometries of cardiac morphology and the latter providing robust features. We have shown the derivation of DNLS in detail and demonstrated that DNLS outperforms two state-of-the-art methods. 


  §.§.§ Acknowledgements.
 
The research was supported by the British Heart Foundation (NH/17/1/32725); National Institute for Health Research (NIHR) Biomedical Research Centre based at Imperial College Healthcare NHS Trust and Imperial College London; and the Medical Research Council, UK. We would like to thank Dr Simon Gibbs, Dr Luke Howard and Prof Martin Wilkins for providing the CMR image data. The Titan Xp for this research was donated by NVIDIA. 

splncs


