DeepIGeoS: A Deep Interactive Geodesic Framework for Medical Image Segmentation

Accurate medical image segmentation is essential for diagnosis, surgical planning and many other applications. Convolutional Neural Networks (CNNs) have become the state-of-the-art automatic segmentation methods. However, fully automatic results may still need to be refined to become accurate and robust enough for clinical use. We propose a deep learning-based interactive segmentation method to improve the results obtained by an automatic CNN and to reduce user interactions during refinement for higher accuracy. We use one CNN to obtain an initial automatic segmentation, on which user interactions are added to indicate mis-segmentations. Another CNN takes as input the user interactions with the initial segmentation and gives a refined result. We propose to combine user interactions with CNNs through geodesic distance transforms, and propose a resolution-preserving network that gives a better dense prediction. In addition, we integrate user interactions as hard constraints into a back-propagatable Conditional Random Field. We validated the proposed framework in the context of 2D placenta segmentation from fetal MRI and 3D brain tumor segmentation from FLAIR images. Experimental results show our method achieves a large improvement from automatic CNNs, and obtains comparable and even higher accuracy with fewer user interventions and less time compared with traditional interactive methods.


INTRODUCTION
S EGMENTATION of anatomical structures is an essential task for a range of medical image processing applications such as image-based diagnosis, anatomical structure modeling, surgical planning and guidance. Although automatic segmentation methods [1] have been investigated for many years, they can rarely achieve sufficiently accurate and robust results to be useful for many medical imaging applications. This is mainly due to poor image quality (with noise, artifacts and low contrast), large variations among patients, inhomogeneous appearances brought by pathology, and variability of protocols among clinicians leading to different definitions of a given structure boundary. Interactive segmentation methods, which take advantage of users' knowledge of anatomy and the clinical question to overcome the challenges faced by automatic methods, are widely used for higher accuracy and robustness [2].
Although leveraging user interactions helps to obtain more precise segmentation [3], [4], [5], [6], the resulting requirement for many user interactions increases the burden on the user. A good interactive segmentation method should require as few user interactions as possible, leading to interaction efficiency. Machine learning methods are commonly used to reduce user interactions. For example, GrabCut [7] uses Gaussian Mixture Models to represent color distributions. It requires the user to provide a bounding box around the object of interest for segmentation and allows additional scribbles for refinement. SlicSeg [8] employs Online Random Forests to segment a Magnetic Resonance Imaging (MRI) volume by learning from user-provided scribbles in only one slice. Active learning is used in [9] to actively select candidate regions for querying the user.
Recently, deep learning techniques with convolutional neural networks (CNNs) have achieved increasing success in image segmentation [10], [11], [12]. They can find the most suitable features through automatic learning instead of manual design. By learning from large amounts of training data, CNNs have achieved state-of-the-art performance for automatic segmentation [12], [13], [14]. One of the most widely used CNNs is the Fully Convolutional Network (FCN) [11]. It outputs the segmentation directly by computing forward propagation only once at the testing time.
Recent advances of CNNs for image segmentation mainly focus on two aspects. The first is to overcome the problem of reduction of resolution caused by repeated combination of max-pooling and downsampling. Though some upsampling layers can be used to recover the resolution, this easily leads to blob-like segmentation results and low accuracy for tiny structures [11]. In [12], [15], dilated convolution is proposed to replace some downsampling layers and it allows exponential expansion of the receptive field without the loss of resolution. However, the CNNs in [12], [15] keep three layers of pooling and downsampling therefore their output resolution is still reduced eight times compared with the input. The second aspect is to enforce inter-pixel dependency to get a spatially regularized result. This helps to recover edge details and reduce noise in pixel classification. DeepLab [16] and DeepMedic [14] used fully connected Conditional Random Fields (CRFs) as a post-processing step. However, the parameters of these CRFs rely on manual tuning which is time consuming and may not ensure optimal values. It was shown in [17] that the CRF can be formulated as a Recurrent Neural Network (RNN) so that it can be trained end-to-end utilizing the back-propagation algorithm. However, this CRF constrains the pairwise potentials as Gaussian functions, which may be too restrictive for some complex cases, and this method does not apply automatic learning to all its parameters. Thus, using more freeform learnable pairwise potential functions and allowing automatic learning of all the parameters can potentially achieve better results.
This paper aims to integrate user interactions into CNN frameworks to obtain accurate and robust segmentation of 2D and 3D medical images, and at the same time, we aim to make the interactive framework more efficient with a minimal number of user interactions by using CNNs. With the good performance of CNNs shown for automatic image segmentation tasks [10], [11], [13], [14], [16], we hypothesize that they can reduce the number of user interactions for interactive image segmentation. However, only a few works have been reported to apply CNNs to interactive segmentation tasks [18], [19], [20], [21].
The contributions of this work are four-fold. 1). We propose a deep CNN-based interactive framework for 2D and 3D medical image segmentation. We use one CNN to get an initial automatic segmentation, which is refined by another CNN that takes as input the initial segmentation and user interactions; 2). We present a new way to combine user interactions with CNNs based on geodesic distance maps that are used as extra channels of the input for CNNs. We show that using geodesic distance can lead to improved segmentation accuracy compared with using Euclidean distance; 3). We propose a resolution-preserving CNN structure which leads to a more detailed segmentation result compared with traditional CNNs with resolution loss, and 4). We extend the current RNN-based CRFs [17] for segmentation so that the back-propagatable CRFs can use user interactions as hard constraints and all the parameters of potential functions can be trained in an end-to-end way. We apply the proposed method to 2D placenta segmentation from fetal MRI and 3D brain tumor segmentation from fluid attenuation inversion recovery (FLAIR) images.

Image Segmentation based on CNNs
Typical CNNs such as AlexNet [22], GoogleNet [23], VGG [24] and ResNet [25] were originally designed for image classification tasks. Some early works adapted such networks for pixel labeling with patch or region-based methods [10], [13]. Such methods achieved higher accuracy than traditional methods that relied on hand-crafted features. However, they suffered from inefficiency for testing. FCNs [11] take an entire image as input and give a dense segmentation. In order to overcome the problem of loss of spatial resolution due to multi-stage max-pooling and downsampling, it uses a stack of deconvolution (a.k.a. upsampling) layers and activation functions to upsample the feature maps. Inspired by the convolution and deconvolution framework of FCNs, a U-shape network (U-Net) [26] and its 3D version [18] were proposed for biomedical image segmentation. A similar network (V-Net) [27] was proposed to segment the prostate from 3D MRI volumes.
To overcome the drawbacks of successive max-pooling and downsampling that lead to a loss of feature map resolution, dilated convolution [12], [15] was proposed to preserve the resolution of feature maps and enlarge the receptive field to incorporate larger contextual information. In [28], a stack of dilated convolutions was used for object tracking and semantic segmentation. Dilated convolution has also been used for instance-sensitive segmentation [29] and action detection from video frames [30].
Multi-scale features extracted from CNNs have been shown to be effective for improving segmentation accuracy [11], [12], [15]. One way of obtaining multi-scale features is to pass several scaled versions of the input image through the same network. The features from all the scales can be fused for pixel classification [31]. In [13], [14], the features of each pixel were extracted from two concentric patches with different sizes. In [32], multi-scale images at different stages were fed into a recurrent convolutional neural network. Another widely used way to obtain multi-scale features is exploiting the feature maps from different levels of a CNN. For example, in [33], features from intermediate layers are concatenated for segmentation and localization. In [11], [12], predictions from the final layer are combined with those from previous layers.

Interactive Image Segmentation
Interactive image segmentation has been widely used in various applications [20], [34], [35]. There are many kinds of user interactions, such as click-based [36], contour-based [4] and bounding box-based methods [7]. Drawing scribbles is user-friendly and particularly popular, e.g., in Graph Cuts [3], GeoS [6], [37], and Random Walks [5]. However, most of these methods rely on low-level features and require a relatively large amount of user interactions to deal with images with low contrast and ambiguous boundaries. Machine learning methods [8], [38], [39] have been proposed to learn from user interactions. They can achieve higher segmentation accuracy with fewer user interactions. However, they are limited by hand-crafted features that depend on the user's experience.
Recently, using deep CNNs to improve interactive segmentation has attracted increasing attention due to CNNs' automatic feature learning and high performance. For instance, 3D U-Net [18] learns from sparsely annotated images and can be used for semi-automatic segmentation. Scribble-Sup [19] also trains CNNs for semantic segmentation supervised by scribbles. DeepCut [20] employs user-provided bounding boxes as annotations to train CNNs for the segmentation of fetal MRI. However, these methods are not fully interactive for testing since they do not accept further interactions for refinement. In [21], a deep interactive object selection method was proposed where user-provided clicks are transformed into Euclidean distance maps and then concatenated with the input of FCNs. However, the Euclidean distance does not take advantage of image context information. In contrast, the geodesic distance transform [6], [37], [40] encodes spatial regularization and contrast-sensitivity but it has not been used for CNNs.

CRFs for Spatial Regularization
Graphical models such as CRFs [12], [41], [42] have been widely used to enhance segmentation accuracy by introducing spatial consistency. In [41], spatial regularization was obtained by minimizing the Potts energy with a min-cut/maxflow algorithm. In [42], the discrete max-flow problem was mapped to its continuous optimization formulation. Such methods encourage segmentation consistency between adjacent pixel pairs with high similarity. In order to better model long-range connections within the image, a fully connected CRF was used in [43] to establish pairwise potentials on all pairs of pixels in the image. To make the inference of this CRF efficient, the pairwise edge potentials were defined by a linear combination of Gaussian kernels in [44]. The parameters of CRFs in these works were manually tuned or inefficiently learned by grid search. In [45], a maximum margin learning method was proposed to learn CRFs using Graph Cuts. Other methods including structured output Support Vector Machines [46], approximate marginal inference [47] and gradient-based optimization [48] were also proposed to learn parameters in CRFs. They treat the learning of CRFs as an independent step after the training of classifiers.
The CRF-RNN network [17] formulated dense CRFs as RNNs so that the CNNs and CRFs can be jointly trained in an end-to-end system for segmentation. However, the pairwise potentials in [17] are limited to weighted Gaussians and not all the parameters are trainable due to the Permutohedral lattice implementation [49]. In [50], a Gaussian Mean Field (GMF) network was proposed and combined with CNNs where all the parameters are trainable. More freeform pairwise potentials for a pair of super-pixels or image patches were proposed in [31], [51], but such CRFs have a low resolution. In [52], a generic CNN-CRF model was proposed to handle arbitrary potentials for labeling body parts in depth images. However, it has not yet been validated with other segmentation applications.

METHOD
The proposed deep interactive segmentation method based on CNNs and geodesic distance transforms (DeepIGeoS) is depicted in Fig. 1. To minimize the number of user interactions, we propose to use two CNNs: an initial segmentation proposal network (P-Net) and a refinement network (R-Net). P-Net takes as input a raw image with C I channels and gives an initial automatic segmentation. Then the user checks the segmentation and provides some interactions (clicks or scribbles) to indicate mis-segmented regions. R-Net takes as input the original image, the initial segmentation and the user interactions to provide a refined segmentation. P-Net and R-Net use a resolution-preserving structure that captures high-level features from a large receptive field without loss of resolution. They share the same structure except the difference in the input dimensions. Based on the initial automatic segmentation obtained by P-Net, the user might give clicks/scribbles to refine the result more than one time through R-Net. Differently from previous works [53] that re-train the learning model each time when new user interactions are given, the proposed R-Net is only trained with user interactions once since it takes a considerable time to re-train a CNN model with a large training set.  To make the segmentation result more spatially consistent and to use scribbles as hard constraints, both P-Net and R-Net are connected with a CRF, which is modeled as an RNN (CRF-Net) so that it can be trained jointly with P-Net/R-Net by back-propagation. We use freeform pairwise potentials in the CRF-Net. The way user interactions are used is presented in 3.1. The structures of 2D/3D P-Net and R-Net are detailed in 3.2. In 3.3, we describe the implementation of our CRF-Net. Training details are described in 3.4.

User Interaction-based Geodesic Distance Maps
In our method, scribbles are provided by the user to refine an initial automatic segmentation obtained by P-Net. A scribble labels a set of pixels as the foreground or back-  Fig. 3. The CNN structure of 2D/3D P-Net with CRF-Net(f). The parameters of convolution layers are (kernel size, output channels, dilation) in dark blue rectangles. Block 1 to block 6 are resolution-preserving. 2D/3D R-Net uses the same structure as 2D/3D P-Net except its input has three additional channels shown in Fig. 2 and the CRF-Net(f) is replaced by the CRF-Net(fu) (Section 3.3).
ground. Interactions with the same label are converted into a distance map. In [21], the Euclidean distance was used due to its simplicity. However, the Euclidean distance treats each direction equally and does not take the image context into account. In contrast, the geodesic distance helps to better differentiate neighboring pixels with different appearances, and improve label consistency in homogeneous regions [6]. GeoF [40] uses the geodesic distance to encode variable dependencies in the feature space and it is combined with Random Forests for semantic segmentation. However, it is not designed to deal with user interactions. We propose to encode user interactions via geodesic distance transforms for CNN-based segmentation.
Suppose S f and S b represent the set of pixels belonging to foreground scribbles and background scribbles, respectively. Let i be a pixel in an image I, then the unsigned geodesic distance from i to the scribble set S(S ∈ {S f , S b }) is: where P i,j is the set of all paths between pixel i and j. p is one feasible path and it is parameterized by s ∈ [0,1]. u(s) = p (s)/ p (s) is a unit vector that is tangent to the direction of the path. If no scribbles are drawn for either the foreground or background, the corresponding geodesic distance map is filled with random numbers. Fig. 2 shows an example of geodesic distance transforms of user interactions. The geodesic distance maps of user interactions and the initial automatic segmentation have the same size as I. They are concatenated with the raw channels of I so that a concatenated image with C I +3 channels is obtained, which is used as the input of the refinement network R-Net.

Resolution-Preserving CNNs using Dilated Convolution
CNNs in our method are designed to capture high-level features from a large receptive field without the loss of resolution of the feature maps. They are adapted from VGG-16 [24] and made resolution-preserving. Fig. 3 shows the structure of 2D and 3D P-Net. In 2D P-Net, the first 13 convolution layers are grouped into five blocks. The first and second blocks have two convolution layers respectively, and each of the remaining blocks has three convolution layers. The size of the convolution kernel is fixed as 3×3 in all these convolution layers. 2D R-Net uses the same structure as 2D P-Net except that its number of input channels is C I +3 and it employs user interactions in the CRF-Net. To obtain an exponential increase of the receptive field, VGG-16 uses a max-pooling and downsampling layer after each block. However, this implementation would decrease the resolution of feature maps exponentially. Therefore, to preserve resolution through the network, we remove the max-pooling and downsampling layers and use dilated convolution in each block.
Let I be a 2D image of size W × H, and let K rq be a square dilated convolution kernel with a size of (2r+1)×(2r+1) and a dilation parameter q, where r ∈ Z and q ∈ Z. The dilated convolution of I with K rq is defined as: For 2D P-Net/R-Net, we set r to 1 for block 1 to block 5, so the size of a convolution kernel becomes 3×3. The dilation parameter in block i is set to: where d ∈ Z is a system parameter controlling the base dilation parameter of the network. We set d=1 in experiments.
The receptive field of a dilated convolution kernel K rq is (2rq+1)×(2rq+1). Let R i × R i denote the receptive field of block i. R i can be computed as: where τ j is the number of convolution layers in block j, with a value of 2, 2, 3, 3, 3 for the five blocks respectively. When r=1, the receptive field size of each block is R 1 =4d+1, Thus, these blocks capture features at different scales. The stride of each convolution layer is set to 1. The number of output channels of convolution in each block is set to a fixed number C. In order to use multi-scale features, we concatenate the features from different blocks to get a composed feature of length 5C. This feature is fed into a classifier that is implemented by two additional layers as shown in block 6 in Fig. 3(a). These two layers use convolution kernels with size of 1×1 and dilation parameter of 0. Block 6 gives each pixel an initial score of belonging to the foreground or background class. In order to get a more spatially consistent segmentation and add hard constraints when scribbles are given, we apply a CRF on the basis of the output from block 6. The CRF is implemented by a recurrent neural network (CRF-Net, detailed in 3.3), which can be jointly trained with P-Net or R-Net. The CRF-Net gives a regularized prediction for each pixel, which is fed into a cross entropy loss function layer.
Similar network structures are used by 3D P-Net/R-Net for 3D segmentation, as shown in Fig. 3(b). To reduce the memory consumption for 3D images, we use one downsampling layer before the resolution-preserving layers and compress the output features of block 1 to 5 by a factor four via 1×1×1 convolutions before the concatenation layer.

Back-propagatable CRF-Net with Freeform Pairwise Potentials and User Constraints
In [17], a CRF based on RNN was proposed and it can be trained by back-propagation. Rather than using Gaussian functions, we extend this CRF so that the pairwise potentials can be freeform functions and we refer to it as CRF-Net(f). In addition, we integrate user interactions in our CRF-Net(f) in the interactive refinement context, which is referred to as CRF-Net(fu). The CRF-Net(f) is connected to P-Net and the CRF-Net(fu) is connected to R-Net.
Let X be the label map assigned to an image I with a label set L = {0, 1, ..., L -1}. The Gibbs distribution P (X = x|I) = 1 Z(I) exp(−E(x|I)) models the probability of X given I in a CRF, where Z(I) is the normalization factor known as the partition function, and E(x) is the Gibbs energy: where the unary potential ψ u (x i ) measures the cost of assigning label x i to pixel i, and the pairwise potential is the cost of assigning labels x i , x j to a pixel pair i, j. N is the set of all pixel pairs. In our method, the unary potential is obtained from P-Net or R-Net that gives classification scores at each pixel. The pairwise potential is: where d ij is the Euclidean distance between pixels i and j. µ(x i , x j ) is the compatibility between the label of i and that of j represented by a matrix of size L × L.
where f i and f j represent the feature vectors of i and j, respectively. The feature vectors can either be learned by a network or be derived from image features such as spatial location with intensity values. For experiments we used the latter one, as in [3], [17], [44] for simplicity and efficiency. Fig. 4. The Pairwise-Net for pairwise potential function f (f ij , d ij ).f ij is the difference of features between a pixel pair i and j. d ij is the Euclidean distance between them.
f (·) is a function in terms off ij and d ij . Instead of defining f (·) as a single Gaussian function [3] or a combination of several Gaussian functions [17], [44], we set it as a freeform function represented by a fully connected neural network (Pairwise-Net) which can be learned during training. The structure of Pairwise-Net is shown in Fig. 4. The input is a vector composed off ij and d ij . There are two hidden layers and one output layer.
Graph Cuts [3], [45] can be used to minimize Eq. (6) when ψ p (·) is submodular [54] such as when the segmentation is binary with µ(·) being the delta function and f (·) being positive. However, this is not the case for our method since we learn µ(·) and f (·) where µ(·) may not be the delta function and f (·) could be negative. Continuous max-flow [42] can also be used for the minimization, but its parameters are manually designed. Alternatively, meanfield approximation [17], [44], [50] is often used to efficiently minimize the energy while allowing learning parameters by back-propagation. Instead of computing P (X|I) directly, an approximate distribution Q(X|I) = i Q i (x i |I) is computed so that the KL-divergence D(Q||P ) is minimized. This yields an iterative update of Q i (x i |I) [17], [44], [50].
where L is the label set. i and j are a pixel pair. For the proposed CRF-Net(fu), with the set of user-provided scribbles S f b = S f ∪ S b , we force the probability of pixels in the scribble set to be 1 or 0. The following equation is used as the update rule for each iteration: where s i denotes the user-provided label of a pixel i that is in the scribble set S f b . We follow the implementation in [17] to update Q through a multi-stage mean-field method in an RNN. Each mean-field layer splits Eq. (8) into four steps including message passing, compatibility transform, adding unary potentials and normalizing [17].

Implementation Details
The raster-scan algorithm [6] was used to compute geodesic distance transforms by applying a forward pass scanning and a backward pass scanning with a 3×3 kernel for 2D and a 3×3×3 kernel for 3D. It is fast due to accessing the image memory in contiguous blocks. For the proposed CRF-Net with freeform pairwise potentials, two observations motivate us to use pixel connections based on local patches instead of full connections within the entire image. First, the permutohedral lattice implementation [17], [44] allows efficient computation of fully connected CRFs only when pairwise potentials are Gaussian functions. However, a method that relaxes the requirement of pairwise potentials as freeform functions represented by a network (Fig. 4) cannot use that implementation and therefore would be inefficient for fully connected CRFs. Suppose an image with size M × N , a fully connected CRF has M N (M N -1) pixel pairs. For a small image with M =N =100, the number of pixel pairs would be almost 10 8 , which requires not only a huge amount of memory but also long computational time. Second, though long-distance dependency helps to improve segmentation in most RGB images [12], [17], [44], this would be very challenging for medical images since the contrast between the target and background is often low [55]. In such cases, long-distance dependency may lead the label of a target pixel to be corrupted by the large number of background pixels with similar appearances. Therefore, to maintain a good efficiency and avoid longdistance corruptions, we define the pairwise connections for one pixel within a local patch centered on that. In our experiments, the patch size is set to 7×7 for 2D images and 5×5×3 for 3D images. We initialize µ(·) as µ( is the Iverson Bracket [17]. A fully connected neural network (Pairwise-Net) with two hidden layers is used to learn the freeform pairwise potential function (Fig. 4). The first and second hidden layers have 32 and 16 neurons, respectively. In practice, this network is implemented by an equivalent fully convolutional neural network with 1×1×1 kernels. We use a pre-training step to initialize the Pairwise-Net with an approximation of a contrast sensitive function [3]: where F is the dimension of the feature vectors f i and f j , and ω and σ are two parameters controlling the magnitude and shape of the initial pairwise function respectively. In this initialization step, we set σ to 0.08 and ω to 0.5 based on experience. Similar to [16], [17], [44], we set f i and f j as values in input channels (i.e, image intensity in our case) of P-Net for simplicity of implementation and for obtaining contrast-sensitive pairwise potentials. use a Stochastic Gradient Descent (SGD) algorithm with a quadratic loss function to pre-train the Pairwise-Net. For pre-processing, all the images are normalized by the mean value and standard variation of the training set. We apply data augmentation by vertical or horizontal flipping, random rotation with angle range [-π/8, π/8] and random zoom with scaling factor range [0.8, 1.25]. We use the cross entropy loss function and SGD algorithm for optimization with minibatch size 1, momentum 0.99 and weight decay 5×10 −4 . The learning rate is halved every 5k iterations. Since a proper initialization of P-Net and CRF-Net(f) is helpful for a faster convergence of the joint training, we train the P-Net with CRF-Net(f) in three steps. First, the P-Net is pre-trained with initial learning rate 10 −3 and maximal number of iterations 100k. Second, the Pairwise-Net in the CRF-Net(f) is pre-trained as described above. Third, the P-Net and CRF-Net(f) are jointly trained with initial learning rate 10 −6 and maximal number of iterations 50k.
After the training of P-Net with CRF-Net(f), we automatically simulate user interactions to train R-Net with CRF-Net(fu). First, P-Net with CRF-Net(f) is used to obtain an automatic segmentation for each training image. It is compared with the ground truth to find mis-segmented regions. Then the user interactions on each mis-segmented region are simulated by randomly sampling n pixels in that region. Suppose the size of one connected under-segmented or over-segmented region is N m , we set n for that region to 0 if N m < 30 and N m /100 otherwise based on experience. Examples of simulated user interactions on training images are shown in Fig. 5. With these simulated user interactions on the initial segmentation of training data, the training of R-Net with CRF-Net(fu) is implemented through SGD, which is similar to the training of P-Net with CRF-Net(f).
We implemented our 2D networks by Caffe 1 [56] and 3D networks by Tensorflow 2 [57] using NiftyNet 3 [58]. Our training process was done via two 8-core E5-2623v3 Intel Haswells and two K80 NVIDIA GPUs and 128GB memory. The testing process with user interactions was performed on a MacBook Pro (OS X 10.9.5) with 16GB RAM and an Intel Core i7 CPU running at 2.5GHz and an NVIDIA GeForce GT 750M GPU. A Matlab and PyQt GUI were developed for 2D and 3D interactive segmentation tasks, respectively. (See supplementary videos) 1

Comparison Methods and Evaluation Metrics
We compared our P-Net with FCN [11] and DeepLab [16] for 2D segmentation and DeepMedic [14] and High-Res3DNet [58] for 3D segmentation. Pre-trained models of FCN 4 and DeepLab 5 based on ImageNet were fine-tuned for 2D placenta segmentation. Since the input of FCN and DeepLab should have three channels, we duplicated each of the gray-level images twice and concatenated them into a three-channel image as the input. DeepMedic and High-Res3DNet were originally designed for multi-modality or multi-class 3D segmentation. We adapted them for single modality binary segmentation. We also compared 2D/3D P-Net with 2D/3D P-Net(b5) that only uses the features from block 5 (Fig. 3) instead of the concatenated multi-scale features.
The proposed CRF-Net(f) with freeform pairwise potentials was compared with: 1). Dense CRF as an independent post-processing step for the output of P-Net. We followed the implementation in [14], [16], [44]. The parameters of this CRF were manually tuned based on a coarse-to-fine search scheme as suggested by [16], and 2). CRF-Net(g) which refers to the CRF that can be trained jointly with CNNs by using Gaussian pairwise potentials [17].
We compared three methods to deal with user interactions. 1). Min-cut user-editing [7], where the initial probability map (output of P-Net in our case) is combined with user interactions to solve an energy minimization problem with min-cut [3]; 2). Using the Euclidean distance of user interactions in R-Net, which is referred to as R-Net(Euc), and 3). The proposed R-Net with the geodesic distance of user interactions.
We also compared DeepIGeoS with several other interactive segmentation methods. For 2D slices, DeepIGeoS was compared with: 1). Geodesic Framework [37] that computes a probability based on the geodesic distance from userprovided scribbles for pixel classification; 2). Graph Cuts [3] that models segmentation as a min-cut problem based on user interactions; 3). Random Walks [5] that assigns a pixel with a label based on the probability that a random walker reaches a foreground or background seed first, and 4). SlicSeg [8] that uses Online Random Forests to learn from the scribbles and predict the remaining pixels. For 3D images, DeepIGeoS was compared with GeoS [6] and ITK-SNAP [59]. Two users (an Obstetrician and a Radiologist) respectively used these interactive methods to segment every test image until the result was visually acceptable.
For quantitative evaluation, we measured the Dice score and the average symmetric surface distance (ASSD).
where R a and R b represent the region segmented by the algorithm and the ground truth, respectively.  where S a and S b represent the set of surface points of the target segmented by the algorithm and the ground truth, respectively. d(i, S b ) is the shortest Euclidean distance between i and S b . We used the Student's t-test to compute the p-value in order to see whether the results of two algorithms significantly differ from each other.

Clinical Background and Experiments Setting
Fetal MRI is an emerging diagnostic tool complementary to ultrasound due to its large field of view and good soft tissue contrast. Segmenting the placenta from fetal MRI is important for fetal surgical planning such as in the case of twin-to-twin transfusion syndrome [60]. Clinical fetal MRI data are often acquired with a large slice thickness for good contrast-to-noise ratio. Movement of the fetus can lead to inhomogeneous appearances between slices. In addition, the location and orientation of the placenta vary largely between individuals. These factors make automatic and 3D segmentation of the placenta a challenging task [61]. Interactive 2D slice-based segmentation is expected to achieve more robust results [8], [53]. The 2D segmentation results can also be used for motion correction and high-resolution volume reconstruction [62]. We collected clinical MRI scans for 25 pregnancies in the second trimester. The data were acquired in axial view with pixel size between 0.7422 mm×0.7422 mm and 1.582 mm×1.582 mm and slice thickness 3 -4 mm. Each slice was resampled with a uniform pixel size of 1 mm×1 mm and cropped by a box of size 172×128 containing the placenta. We used 17 volumes with 624 slices for training, three volumes with 122 slices for validation and five volumes with 179 slices for testing. The ground truth was manually delineated by a experienced Radiologist. Fig. 6 shows the automatic segmentation results obtained by different networks. It shows that FCN is able to capture the main region of the placenta. However, the segmentation results are blob-like with smooth boundaries. DeepLab is better than FCN, but its blob-like results are similar to those of FCN. This is mainly due to the downsampling and upsampling procedure employed by these methods. In contrast, 2D P-Net(b5) and 2D P-Net obtain more detailed results. It can be observed that 2D P-Net achieves better results than the other three networks. However, there are still some obvious mis-segmented regions by 2D P-Net. Table 1 presents a quantitative comparison of these networks based on all the testing data. 2D P-Net achieves a Dice score of 84.78±11.74% and an ASSD of 2.09±1.53 pixels, and it performs better than the other three networks.

Automatic Segmentation by 2D P-Net with CRF-Net(f)
Based on 2D P-Net, we investigated the performance of different CRFs. A visual comparison between Dense CRF, CRF-Net(g) with Gaussian pairwise potentials and CRF-Net(f) with freeform pairwise potentials is shown in Fig. 7. In the first column, the placenta is under-segmented by 2D P-Net. Dense CRF leads to very small improvements on the result. CRF-Net(g) and CRF-Net(f) improve the result by preserving more placenta regions, and the later shows a better segmentation. In the second column, 2D P-Net obtains an over-segmentation of adjacent fetal brain and maternal tissues. Dense CRF does not improve the segmentation noticeably, but CRF-Net(g) and CRF-Net(f) remove more oversegmentated areas. CRF-Net(f) shows a better performance than the other two CRFs. The quantitative evaluation of these three CRFs is presented in Table 1, which shows Dense CRF leads to a result that is very close to that of 2D P-Net (pvalue > 0.05), while the last two CRFs significantly improve the segmentation (p-value < 0.05). In addition, CRF-Net(f) is better than CRF-Net(g). Fig. 7 and Table 1 indicate that large mis-segmentation exists in some images, therefore we use 2D R-Net with CRF-Net(fu) to refine the segmentation interactively in the following.    Fig. 8 shows examples of interactive refinement based on 2D R-Net with CRF-Net(fu) that uses freeform pairwise potentials and employs user interactions as hard constraints. The first column in Fig. 8 shows initial segmentation results obtained by 2D P-Net + CRF-Net(f). The user provides clicks/scribbles to indicate the foreground (red) or the background (cyan). The second to last column in Fig. 8 show the results for five variations of refinement. These refinement methods correct most of the mis-segmented areas but perform at different levels in dealing with local details, as indicated by white arrows. Fig. 8 shows 2D R-Net with geodesic distance performs better than min-cut user-editing and 2D R-Net(Euc) that uses Euclidean distance. CRF-Net(fu) can further improve the segmentation. For quantitative comparison, we measured the segmentation accuracy after the first iteration of user refinement (giving user interactions to mark all the main mis-segmented regions and applying refinement once), in which the same initial segmentation and the same set of user interactions were used by the five refinement methods. The results are presented in Table 2, which shows the combination of the proposed 2D R-Net using geodesic distance and CRF-Net(fu) leads to more accurate segmentations than the other refinement methods with the same set of user interactions. The Dice score and ASSD of 2D R-Net + CRF-Net(fu) are 89.31±5.33% and 1.22±0.55 pixels, respectively. Fig. 9 shows a visual comparison between DeepIGeoS and Geodesic Framework [37], Graph Cuts [3], Random Walks [5] and SlicSeg [8] for 2D placenta segmentation. The first row shows the initial scribbles and the resulting segmentation. Notice no initial scribbles are needed for DeepIGeoS. The second row shows refined results, where DeepIGeoS only needs two short strokes to get an accurate segmentation, while the other methods require far more scribbles to get similar results. Quantitative comparison of these methods based on the final segmentation given by the two users is presented in Fig. 10. It shows these methods

Clinical Background and Experiments Setting
Gliomas are the most common brain tumors in adults with little improvement in treatment effectiveness despite considerable research works [63]. With the development of medical imaging, brain tumors can be imaged by different MR protocols with different contrasts. For example, T1weighted images highlight enhancing part of the tumor and FLAIR acquisitions highlight the peritumoral edema. Segmentation of brain tumors can provide better volumetric measurements and therefore has enormous potential value for improved diagnosis, treatment planning, and followup of individual patients. However, automatic brain tumor segmentation remains technically challenging because 1) the size, shape, and localization of brain tumors have considerable variations among patients; 2) the boundaries between adjacent structures are often ambiguous.
In this experiment, we investigate interactive segmentation of the whole tumor from FLAIR images. We used the 2015 Brain Tumor Segmentation Challenge (BraTS) [63] training set with images of 274 cases. The ground truth were manually delineated by several experts. Differently from previous works using this dataset for multi-label and multimodality segmentation [14], [64], as a first demonstration of deep interactive segmentation in 3D, we only use FLAIR images in the dataset and only segment the whole tumor. We randomly selected 234 cases for training and used the remaining 40 cases for testing. All these images had been skull-stripped and resampled to size of 240×240×155 with isotropic resolution 1mm 3 . We cropped each image based on the bounding box of its non-zero region. The feature channel number of 3D P-Net and R-Net was C = 16. Fig. 11 shows examples of automatic segmentation of brain tumor by 3D P-Net, which is compared with DeepMedic [14], HighRes3DNet [58] and 3D P-Net(b5). In the first column, DeepMedic segments the tumor roughly, with some missed regions near the boundary. High-Res3DNet reduces the missed regions but leads to some over-segmentation. 3D P-Net(b5) obtains a similar result to that of HighRes3DNet. In contrast, 3D P-Net achieves a more accurate segmentation, which is closer to the ground truth. More examples in the second and third column in  Quantitative evaluation of these four networks is presented in Table 3. DeepMedic achieves an average dice score of 83.87%. HighRes3DNet and 3D P-Net(b5) achieve similar performance, and they are better than DeepMedic. 3D P-Net outperforms these three counterparts with 86.68±7.67% in terms of Dice and 2.14±2.17 pixels in terms of ASSD. Note that the proposed 3D P-Net has far fewer parameters compared with HighRes3DNet. It is more memory efficient and therefore can perform inference on a 3D volume in interactive time.

Automatic Segmentation by 3D P-Net with CRF-Net(f)
Since CRF-RNN [17] was only implemented for 2D, in the context of 3D segmentation we only compared 3D CRF-Net(f) with 3D Dense CRF [14] that uses manually tuned parameters. Visual comparison between these two types of CRFs working with 3D P-Net is shown in Fig. 12. It can be observed that CRF-Net(f) achieves more noticeable improvement compared with Dense CRF that is used as post-processing without end-to-end learning. Quantitative measurement of Dense CRF and CRF-Net(f) is listed in Table 3. It shows that only CRF-Net(f) obtains significantly better segmentation than 3D P-Net with p-value < 0.05.    Fig. 13 shows examples of interactive refinement of brain tumor segmentation using 3D R-Net with CRF-Net(fu). The initial segmentation is obtained by 3D P-Net + CRF-Net(f).
With the same set of user interactions, we compared the refined results of min-cut user-editing and four variations of 3D R-Net: using geodesic or Euclidean distance transforms with or without CRF-Net(fu). Fig. 13 shows that min-cut user-editing achieves a small improvement. It can be found that more accurate results are obtained by using geodesic distance than using Euclidean distance, and CRF-Net(fu) can further help to improve the segmentation. For quantitative comparison, we measured the segmentation accuracy after the first iteration of refinement, in which the same set of scribbles were used for different refinement methods. The quantitative evaluation is listed in Table 4, showing that the proposed 3D R-Net with geodesic distance and CRF-Net(fu) achieves higher accuracy than the other variations with a Dice score of 89.93±6.49% and ASSD of 1.43±1.16 pixels.

Comparison with Other Interactive Methods
Fig. 14 shows a visual comparison between GeoS [6], ITK-SNAP [59] and DeepIGeoS. In the first row, the tumor has a good contrast with the background. All the compared methods achieve very accurate segmentations. In the second row, a lower contrast makes it difficult for the user to identify the tumor boundary. Benefited from the initial tumor boundary that is automatically identified by 3D P-Net, DeepIGeoS outperforms GeoS and ITK-SNAP. Quantitative comparison is presented in Fig. 15. It shows DeepIGeoS achieves higher accuracy compared with GeoS and ITK-SNAP. In addition, the user time for DeepIGeoS is about one third of that for the other two methods. Supplementary video 2 shows more examples of DeepIGeoS for 3D brain tumor segmentation.

CONCLUSION
In this work, we presented a deep learning-based interactive framework for 2D and 3D medical image segmentation. We proposed a P-Net to obtain an initial automatic segmentation and an R-Net to refine the result based on user interactions that are transformed into geodesic distance maps and then integrated into the input of R-Net. We also proposed a resolution-preserving network structure with dilated convolution for dense prediction, and extended the existing RNN-based CRF so that it can learn freeform pairwise potentials and take advantage of user interactions as hard constraints. Segmentation results of the placenta from 2D fetal MRI and brain tumors from 3D FLAIR images show that our proposed method achieves better results than automatic CNNs. It requires far less user time compared with traditional interactive methods and achieves higher accuracy for 3D brain tumor segmentation. The framework can be extended to deal with multiple organs in the future.