1 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 semiautomatic CMR image segmentation have been proposed, including deformable models [1], atlasbased image registration models [2] as well as statistical shape and appearance models [3]
. More recently, deep learningbased methods have achieved stateoftheart performance in the CMR domain
[4]. However, the above approaches for CMR image segmentation have multiple drawbacks. First, they tend to focus on left ventricle (LV) [1]. 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 [1, 2, 3]. This becomes less feasible in populationlevel applications involving hundreds or thousands of CMR images. Third, existing techniques have been mainly developed and validated using normal (healthy) hearts [1, 2, 4]. 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 crossentropy and classbalanced sigmoid crossentropy. 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 stateoftheart approaches.
2 Modelling biventricular anatomy in patients with PH
To illustrate cardiac morphology in patients with PH, Fig 1 shows the difference in CMR images from a representative healthy subject and a PH subject. In health, the RV is crescentic in shortaxis views and triangular in longaxis views, wrapping around the thickerwalled 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.
3 Methodology
Nested level set approach: We view image segmentation in PH as a multiregion image segmentation problem. Let denote an input image defined on the domain . We segment the image into a set of pairwise disjoint region , with , . The segmentation task can be solved by computing a labelling function that indicates which of the regions each pixel belongs to: . The problem is then formulated as an energy minimisation problem consisting of a data term and a regularisation term
(1) 
The data term, 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, corresponds to the negative logarithm of the conditional probability for a specific pixel color at the given location within region . Here we refer to as region feature. The second term, is the perimeter of the segmentation region , weighted by the nonnegative function . 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 can be an edge detection function which favours boundaries that have strong gradients of the input image . Here we refer to as edge feature.
We apply the variational level set method [5, 6] to (1) in this study. Because a PH heart can be implicitly represented by two nested level lines of a continuous level set function ( in Fig 2). Note that the nested level set idea present here is inspired from previous work [1, 7]. Our approach uses features learned from many images while previous work only consider single image. With the idea, we are able to approximate the multiregion segmentation energy (1) by using only one continuous function. The computational cost is thus small. Now assume that the contours in the image can be represented by level lines of the same Lipschitz continuous level set function . With distinct levels , the implicit function partitions the domain into disjoint regions, together with their boundaries (see Fig 2
right). We can then define the characteristic function
for each region as(2) 
where is the onedimensional Heaviside function that takes on either 0 or 1 over the whole domain . Due to the nondifferentiate nature of it is usually approximated by its smooth version for numerical calculation [7]. Note that in (2) is automatically satisfied, meaning that the resulting segmentation will not produce a vacuum or an overlap effect. That is, by using (2) and hold all the time. With the definition of , we can readily reformulate (1) in the following new energy minimisation problem
(3) 
Note that (3) differs from (1) in multiple ways due to the use of the smooth function and characteristic function (2). First, the variable to be minimised is the regions in (1) while the smooth function in (3). Second, the minimisation domain is changing from over in (1) to over in (3). Third, (1) uses an abstract for the weighted length of the boundary between two adjacent regions, while (3) represents the weighted length with the coarea formula, i.e. . Finally, the upper limit of summation in the regularisation term of (1) is while in that of (3). So far, the region features and the edge feature 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 (3). Let us formulate the learning problem as follows: we denote the input training data set by , where sample is the raw input image, , is the ground truth region labels ( regions) for image , and , is the ground truth binary edge map for . We denote all network layer parameters as W
and propose to minimise the following objective function via the (backpropagation) stochastic gradient descent
(4) 
where is the region associated crossentropy loss that enables the network to learn region features, while is the edge associated crossentropy loss for learning edge features. The weight balances the two losses. By minimising (4), the network is able to output joint region and edge probability maps simultaneously. In our imagetoimage training, the loss function is computed over all pixels in a training image , a region map , and an edge map , . The definitions of and are given as follows.
(5) 
where denotes the pixel index, and is the channelwise softmax probability provided by the network at pixel for image . The edge loss is
(6) 
For a typical CMR image, the distribution of edge and nonedge pixels is heavily biased. Therefore, we use the strategy [8] to automatically balance edge and nonedge classes. Specifically, we use a classbalancing weight . Here, and , where and respectively denote edge and nonedge ground truth label pixels. is the pixelwise sigmoid probability provided by the network at nonedge pixel for image .
In Fig 3
, we show the network architecture for automatic feature extraction, which is a fully convolutional network (FCN) and adapted from the Unet architecture
[9]. Batchnormalisation (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
. Every layer whose label is prefixed with ‘C’ performs the operation: convolution BNReLU, 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 (33/2). The arrows represent (33/1) convolutional layers (C14ae) 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 (11/1) convolution with a softmax activation and a sigmoid activation, producing the blue feature map volume with a depth , corresponding to (3) region features and an edge feature of an image.After the network is trained, we deploy it on the given image in the validation set and obtain the joint region and edge probability maps from the last convolutional layer
(7) 
where CNN() denotes the trained network.
is a vector region probability map including
(number of regions) channels, while is a scalar edge probability map. These probability maps are then fed to the energy (3), in which and . With all necessary elements at hand, we are ready to minimise (3) next.Optimisation: The minimisation process of (3) entails the calculus of variations, by which we obtain the resulting EulerLagrange (EL) equation with respect to the variable . A solution () to the EL equation is then iteratively sought by the following gradient descent method
(8) 
where is the weighted curvature that can be numerically implemented by the finite difference method on a halfpoint grid [10]. is the derivative of , which is defined in [7].
At steady state of (8), a local or global minimiser of (3) can be found. Note that the energy (3) 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 () such that the algorithm converges to the solution we want. We tackle this problem by thresholding the region probability map and then computing the signed distance function (SDF) from the binary image using the fast sweeping algorithm. The resulting SDF is then used as for (8). In this way, the whole optimisation process is fully automated.
4 Experimental results
Data: Experiments were performed using shortaxis CMR images from 430 PH patients. For each patient 10 to 16 shortaxis slices were acquired roughly covering the whole heart. Each shortaxis image has resolution of . Due to the large slice thickness of the shortaxis slices and the interslice 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 semiautomatic 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 enddiastolic (ED) and endsystolic (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 preprocessing, all training images were reshaped to the same size of
with zeropadding, and image intensity was normalised to the range of
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 (3). They are the weighting parameter (1), regularisation parameter (1.5), two levels (0) and (8), time step (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 (4) is set to 1. We perform data augmentation onthefly, 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.
Methods  LV&RV Cavities  Myocardium  Time 

Vanilla CNN [4]  0.9020.047  0.7030.091  0.06s 
CRFCNN [11]  0.9110.045  0.7120.082  2s 
Proposed DNLS  0.9250.032  0.7720.058  5s 
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 [4], the code of which is publicly available. DNLS was also compared with the vanilla CNN with a conditional random field (CRF) [11] refinement (CRFCNN). In Fig 4, visual comparison suggests that DNLS provides significant segmentation improvements over CNN and CRFCNN. For example, at the base of the right ventricle both CNN and CRFCNN fail to retain the correct anatomical relationship between endocardium and epicardium portraying the endocardial border outside the epicardium. CRFCNN by contrast retains the endocardial border within the epicardium, as described in the ground truth. In Table 1, 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.
5 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 stateoftheart 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.
References
 [1] Feng, C., Li, C., Zhao, D., Davatzikos, C., Litt, H.: Segmentation of the left ventricle using distance regularized twolayer level set approach. In: MICCAI. (2013) 477–484
 [2] Bai, W., Shi, W., O’Regan, D.P., Tong, T., Wang, H., JamilCopley, S., Peters, N.S., Rueckert, D.: A probabilistic patchbased label fusion model for multiatlas segmentation with registration refinement: application to cardiac mr images. IEEE transactions on medical imaging 32(7) (2013) 1302–1315
 [3] Albà, X., Pereañez, M., Hoogendoorn, C., et al.: An algorithm for the segmentation of highly abnormal hearts using a generic statistical shape model. IEEE transactions on medical imaging 35(3) (2016) 845–859
 [4] Bai, W., Sinclair, M., Tarroni, G., et al.: Humanlevel cmr image analysis with deep fully convolutional networks. Journal of Cardiovascular Magnetic (2018)
 [5] Duan, J., Pan, Z., Yin, X., Wei, W., Wang, G.: Some fast projection methods based on chanvese model for image segmentation. EURASIP Journal on Image and Video Processing (2014) 1–16
 [6] Tan, L., Pan, Z., Liu, W., Duan, J., Wei, W., Wang, G.: Image segmentation with depth information via simplified variational level set formulation. Journal of Mathematical Imaging and Vision 60(1) (2018) 1–17
 [7] Chung, G., Vese, L.A.: Image segmentation using a multilayer levelset approach. Computing and visualization in science 12(6) (2009) 267–285
 [8] Xie, S., Tu, Z.: Holisticallynested edge detection. In: ECCV. (2015) 1395–1403
 [9] Simonyan, K., Zisserman, A.: Very deep convolutional networks for largescale image recognition. CoRR abs/1409.1556 (2014)

[10]
Duan, J., Haines, B., Ward, W.O., Bai, L.:
Surface reconstruction from point clouds using a novel variational
model.
In: International Conference on Innovative Techniques and Applications of Artificial Intelligence. (2015) 135–146
 [11] Krähenbühl, P., Koltun, V.: Efficient inference in fully connected crfs with gaussian edge potentials. In: NIPS. (2011) 109–117