Automated segmentation of objects or equivalently boundary surfaces plays a very import role in quantitative image analysis. In several years, deep learning based method for semantic segmentation has become very popular in computer vision and medical imaging. The fully convolutional networks (FCN) , and then U-net 
for medical image segmentation have been proposed. All these methods model the segmentation problem as a pixel-wise or voxel-wise classification problem. Popular loss functions for training these networks include cross entropy (CE), weighted cross entropy (WCE) for imbalanced classes, and dice similarity coefficient (DSC). As there is no explicit constraint among the labeling of different pixels, usually post processing, such as morphological operation, is required to get reasonable prediction.
It is well known that prior knowledge, like shape priors, can help significantly improve segmentation performance, especially in medical scenario the target objects often lack strong boundaries, and/or present similar intensity profiles with multiple crowed nearby tissues. Some studies have investigated along this line. For example, Milletari et al. 
proposed to integrate into CNN the statistical shape priors, which were obtained from principal components analysis (PCA). Although the robustness of segmentation was improved, the mean DSC did not increase. Ravishankaret al.  utilized a shape regularization network to project the preliminary prediction into shape space, and a shape loss was added to make the prediction more fit into the desired shape space. Chen et al.  proposed an end-to-end approach with CNN and Conditional Random Fields (CRFs) for segmentation. However, the inference of CRFs was approximate and not globally optimal.
Inspired by the graph search method , Shah et al.  first modeled the terrain-like surfaces segmentation as direct surface identification using regression. The network consists of an FCN followed by fully-connected layers. The network is very light weighted and no post processing is required. Surprisingly the results are very promising . It is well known that U-net outperforms FCN because U-net has an additional expansive path such that features of high resolution can be learned and then better prediction accuracy can be achieved. To improve segmentation accuracy further, it is natural to think to replace the FCN with a U-net. However, it is not reasonable to concatenate a U-net with fully-connected layers, as the invariance of feature maps in the original resolution is supposed to be much less than that in the low resolution, such that there would be much more chance that the fully-connected layers heavily overfit to the training data. Another drawback of Shah et al’s work  is that the surface smoothness is implicitly learned within the whole network (mainly within fully-connected layers) as a black box. It is hard to decode after training.
To resolve problems mentioned above, we propose to explicitly model the surface segmentation problem as a quadratic programming, which can be solved with guaranteed global optimality. The whole network can be trained end-to-end. Our contributions are in three folds: 1) The first time to parameterize the output of the U-net as Gaussians (the mean represents surface position prediction from the U-net and the standard deviation encodes the prediction confidence), which converts the description from the discrete to the continuous space; 2) The parameter that controls smoothness can be learned and has a clear meaning; 3) The method works in the continuous space and enables sub-pixel segmentation.
We first define the surface segmentation problem. A -D image can be viewed as a
-D tensor. A terrain-like surface in is oriented and shown in Fig. 1. Let , and denote the image sizes in three dimensions, respectively. Let all column index be a set . The surface is defined by , . Thus any feasible surface in intersects with each column exactly once. Generally in surface segmentation , the problem is formulated as minimizing the energy function
where the unary term is the energy when considering each column independently, and the pairwise energy term penalizes discontinuity of surface position among adjacent columns. The design of and will be detailed in Section 2.1.
2.1 Proposed Inference Pipeline
One should note that the proposed method can be applied in both -D and -D. For the purpose of proof-of-concept and clear explanation, the detailed description of the proposed method and all experiments are done in -D. The inference pipeline of the proposed method is demonstrated in Fig. 1(a). The trained U-net takes in the original image , and outputs the discrete probability map . Ideally, for each image column, the probability for the target surface position is high, and it is gradually reduced on the positions away from the one on the surface, as demonstrated in Fig. 1(b). We thus propose a block to convert the discrete probability map to a Gaussian parameterization , where specifies the mean surface position on each column and is the corresponding standard deviation. The Gaussian parameterization is then fed into the trained smoothing block (SB) , which incorporate the learned surface smoothness priors to infer the optimal target surface. Next, we detail the novel blocks in our deep optimal surface segmentation neural network.
D2C block The D2C block is basically designed to convert the discrete probability map of each column , which is output from U-net, to a continuous representation (Gaussian is utilized in our design). This enables optimizing directly on the surface position and sub-pixel accuracy prediction. The proposed conversion is realized by fitting a continuous Gaussian function to the discrete probability map
, which can be thought as discrete samples of a continuous Gaussian probability density function. Recall that one dimensional Gaussian function has the formulaand then where , , and . In our setting, for each column, we have samples of . We can define an error function , namely Then minimizing the weighted mean square error (MSE)
, one can get the estimates ofby solving a set of three linear equations and then . The problem is very similar to least square problem 
. In our implementation, a linear transform is utilized to normalize the probability map for each column to the range [0,1], then we can ignore the magnitude. As the computation of each column is independent, it is straightforward to be extended to -D.
Smoothing Block (SB) To integrate the surface segmentation model (Eq. 1) with smoothness priors, we define the energy function as
where and are defined as and is the set of neighbor columns.
For simplicity, the nearest neighbor pairs, i.e. , are considered as the set of neighbor columns. The whole energy in Eq. 2 can be reformulated as the standard quadratic form It can be proved that the Hessian matrix is positive definite by using Gershgorin circle theorem and then the energy function is convex. The gradient is Let the gradient to be zero, we have the global optimal solution
Another advantage of the proposed energy formulation is the optimal solution can be calculated in one step, i.e. we do not need to make use of a recurrent neural network (RNN) to implement SB. It is also straightforward to implement the smoothing block in 3D.
2.2 Training Strategy
2.2.1 U-net Pre-training
A common U-net architectures is utilized to generate the discrete probability map for the input image (). The detailed architecture is demonstrated in Fig. 10
. In the proposed method, the softmax layer works on eachcolumn, not on each pixel. The rational is that we assume the target surface intersects with each column by exactly once, and so the probabilities are normalized within each column. Also we assume the U-net should output a Gaussian shaped probability map for each column, which mimics the Bayesian learning for each column and shares merits with knowledge distillation  and distillation defense 
. To encourage the U-net outputs reasonable probability maps, the Kullback–Leibler divergence (KLD) loss is utilized for the U-net pre-training. KLD is a measure of how one probability distribution is different from a reference probability distribution. It is equivalent to the Cross Entropy when the reference is a Kronecker delta function. We propose to relax the delta function to a Gaussian distribution such that the proposed D2C block can work properly. We pickto be around 0.1 times of the column length and our method is insensitive to . Then we will have one Gaussian distribution ground truth for each column, denoted as . One illustration is demonstrated in Fig. 1(b). We denote the output from the U-net as and the loss for the pre-training is formulated as
2.2.2 Fine Tuning
During the fine tunning phase, the mean square surface position error (MSE) is utilized as the loss function, which is formulated as where denotes the ground truth surface positions. The fine tuning of the whole network proceeds in an alternation fashion (Fig. 4).
The validation data is utilized to train the SB, and the training data is for U-net training. As SB only has one parameter- to train, the overfitting chance of it is very low. Also the U-net is not trained on validation data, the learned should be more representative in the wild. Otherwise if fine tuning U-net+SB simultaneously on the training data, the learned is generally smaller than necessary, as the pre-trained U-net generally has fit the training data well and SB does not play an important role.
3.1 SD-OCT Retinal Layer Segmentation
The proposed method was applied to retinal layer segmentation in SD-OCT images, which were obtained from the public dataset . Since the manual tracings were only available for a region centered at the fovea, subvolumes of size were extracted around the fovea. The dataset was randomly divided into 3 sets: 1) Training set - 263 volumes (79 normal, 184 with age-related macular degeneration (AMD)); 2)Validation set - 59 volumes (18 normal, 41 AMD); 3) Testing set - 58 volumes (18 normal, 40 AMD). The surfaces considered are S2-Inner Aspect of Retinal Pigment Epithelium Drusen Complex (IRPE) and S3-Outer Aspect of Bruch Membrane (OBM) as shown in Fig. 5, which are very challenging to segment.
Pre-processing and Augmentation The intensity of each slice was normalized to the range [-1, 1]. No additional pre-processing methods were utilized. For the purpose of pre-training the U-Net, the standard deviation of the Gaussian model of the surface position on each column was set . We augmented the training data by applying mirroring along the horizontal dimension, and 2 random translations along the axial dimension on the original image and the mirrored image.
Hyperparameters All training utilized Adam optimizer. For U-net pre-training, the learning rate was , and the training ran for 2000 epochs. For fine-tuning, the learning rate of U-net () was , the learning rate of SB ( was , the total epochs for the two phases were: , and . The initial smoothness parameter was set to .
|Surface||W/O, normal||W/, normal||W/O, AMD||W/, AMD|| normal|| AMD|
Results Unsigned mean surface positioning error (UMSP)  is utilized for evaluation of segmentation accuracy. The quantitative results are summarized in Table 1. We compare to another deep learning based method , which is the state-of-the-art on this dataset. As for the IRPE surface, the proposed method without SB outperforms . With SB plugged in, the performance can be improved further: for the normal cases, the UMSP can be improved by 39%; for the AMD cases, we can achieve a similar improvement (40%). The segmentation of the OBM surface is more challenging than that of the IRPE surface. One can notice that the learned SB can improve the segmentation accuracy consistently. As for the OBM surface, compared to , the proposed method achieves 41% improvement on the normal cases and performs comparably on the AMD cases. An interesting observation is that the learned are around 0.05 and 1.0 for IRPE and OBM, respectively, which coincides with the human observation that generally less context information exists for OBM. Sample segmentation results are illustrated in Fig. 5, Fig. 6 and Fig. 7.
3.2 IVUS Vessel Wall Segmentation
The data used for this experiment was obtained from the standardized evaluation of IVUS image segmentation database . In this experiment, the dataset B was used. This dataset consists of 435 images with a size of , as well as the respective expert manual tracings of lumen and media surfaces. It comprises two groups - a training set (109 slices) and a testing set (326 slices).
The experiment with the proposed method was conducted in conformance with the directives provided for the IVUS challenge. In our experiment, we randomly split the 109 training slices into 100 slices for training and 9 slices for validation. Pre-processing and Augmentation Each slice was transformed to be represented in the polar coordinate system with a size of , as illustrated in Fig. 8. The intensity of each slice was normalized to have a zero mean and a unit standard deviation. The Gaussian truth was generated for each column using , as it has a shorter column (128 vs 512 in the SD-OCT data). As the number of training data was limited, we augmented the data on the fly by random combinations of various operations including mirroring, circulation shifting along the polar dimension, adding Gaussian noises (mean=0, std=0.1), adding Salt and Pepper noises (5%), and cropping (90% of the original size) and then resizing (). All training parameters were the same as those for the SD-OCT segmentations, except that the number of epochs for the pre-training was 4000.
is the state-of-the-art method for this IVUS dataset. It is an expectation maximization based method and issemi automated. VGG-U-net denotes a deep learning based method . The state-of-the-art fully automation method is a graph search based method working in irregularly sampled space 
, the unary energy of which is learned with random forests. From Table2, one can find that the proposed method outperforms the sub-voxel graph based method  for the lumen surface and the performance on the media surface is comparable. Compared to the state-of-the-art semi-automate method P3, the performance of the proposed is comparable for the lumen surface and is marginally inferior on the media surface. The proposed method outperforms the VGG-U-net by a big margin. Sample segmentation results are illustrated in Fig. 8 and Fig. 9.
As to the inference computation time, the proposed method needs more overhead than the VGG-U-net method (0.21 vs 0.09 sec/slice). The overhead is mainly from the Hessian matrix computation, as well as two separate runs of the program for two surfaces. While compared to P3 [3, 4] (8.6 sec/slice) and the graph based method  (187.4 sec/slice), the proposed method is highly efficient.
In this paper, a novel segmentation model based on a convolutional neural network (CNN) and a learnable surface smoothing block is proposed to tackle the surface segmentation problem with end-to-end training. To the best of our knowledge, this is the first study for surface segmentation which can achieve guaranteed globally optimal solutions using deep learning. Experiments on SD-OCT retinal layer segmentation and IVUS vessel wall segmentation demonstrated very promising results. The proposed method is applicable to -D and -D.
-  (2018) Automatic detection of lumen and media in the ivus images using u-net with vgg16 encoder. arXiv preprint arXiv:1806.07554. Cited by: §3.2, Table 2.
-  (2014) Standardized evaluation methodology and reference database for evaluating ivus image segmentation. Comput Med Imaging Graph 38 (2), pp. 70–90. Cited by: §3.2.
-  (2006) Intravascular ultrasound image segmentation: a three-dimensional fast-marching method based on gray level distributions. IEEE Trans. Med. Imag 25 (5), pp. 590–601. Cited by: §3.2, §3.2, Table 2.
-  (2010) Fast-marching segmentation of three-dimensional intravascular ultrasound images: a pre-and post-intervention study. Med. Phys. 37 (7Part1), pp. 3633–3647. Cited by: §3.2, §3.2, Table 2.
-  (2019) An end-to-end approach to semantic segmentation with 3d cnn and posterior-crf in medical images. In Medical Imaging meets NIPS Workshop at NIPS 2018, Cited by: §1.
-  (2014) Quantitative classification of eyes with and without intermediate age-related macular degeneration using optical coherence tomography. Ophthalmology 121 (1), pp. 162–172. Cited by: §3.1.
-  (2009) Automated 3-d intraretinal layer segmentation of macular spectral-domain optical coherence tomography images. IEEE Trans. Med. Imag 28 (9), pp. 1436–1447. Cited by: §3.1.
-  (2011) A simple algorithm for fitting a gaussian function [dsp tips and tricks]. IEEE Signal Process Mag 28 (5), pp. 134–137. Cited by: §2.1.
-  (2015) Distilling the knowledge in a neural network. In Deep Learning and Representation Learning Workshop at NIPS 2014, Cited by: §2.2.1.
-  (2006) Optimal surface segmentation in volumetric images-a graph-theoretic approach. IEEE Trans Pattern Anal Mach Intell 28 (1), pp. 119–134. Cited by: §1, §2.
-  (2015) Fully convolutional networks for semantic segmentation. In CVPR 2015, pp. 3431–3440. Cited by: §1.
-  (2017) Integrating statistical prior knowledge into convolutional neural networks. In MICCAI 2017, pp. 161–168. Cited by: §1.
-  (2016) Distillation as a defense to adversarial perturbations against deep neural networks. In 2016 IEEE Symposium on Security and Privacy (SP), pp. 582–597. Cited by: §2.2.1.
-  (2017) Learning and incorporating shape models for semantic segmentation. In MICCAI 2017, pp. 203–211. Cited by: §1.
-  (2015) U-net: convolutional networks for biomedical image segmentation. In MICCAI 2015, pp. 234–241. Cited by: §1.
-  (2019) Optimal surface segmentation with convex priors in irregularly sampled space. Med. Image Anal.. Cited by: §3.2, §3.2, Table 2.
-  (2018) Multiple surface segmentation using convolution neural nets: application to retinal layer segmentation in oct images. Biomed. Opt. Express 9 (9), pp. 4509–4526. Cited by: §1, §3.1, Table 1.