I Introduction
Fullyconvolutional neural networks (FCNs) have recently become the de facto standard for medical image segmentation. These models simultaneously detect and segment objects of interest from raw image data, eliminating the need to design applicationspecific features. Current trends suggest that FCN models are capable of performing a wide variety of organ segmentation tasks, but they are limited by the scarcity of training data and the computational issues inherent in volumetric image processing. Additionally, medical images typically exhibit severe class imbalance, for which the usual segmentation loss functions are poorly suited. First, this paper addresses the data issue by creating a new dataset, using manual annotation for some organs and unsupervised methods for others, then performing 3D augmentations on this dataset using the GPU. Next, the computational issues are addressed with a small, memoryefficient model trained with a new loss function which excels at unbalanced segmentation. The end result is a fast and robust model for multiorgan segmentation.
Typical FCNs contain millions of trained parameters across dozens of layers. Training is almost always accelerated by graphics processing units (GPUs), which provide massive parallelism and high memory bandwidth. However, complex FCNs consume a prohibitive amount of highbandwidth memory when extended to process 3D images. This has led to a variety of clever approaches for applying 2D models to 3D data. In contrast, we argue that simpler 3D models can achieve competitive accuracy when combined with extensive data augmentation. By randomly transforming images at each training iteration, the network optimizes over a practically infinite distribution of training images. We leverage GPU texture sampling hardware to perform 3D data augmentation at negligible computational cost. By combining geometric and photometric transformations, our implementation minimizes memory usage and accelerates training. As part of our data augmentation, we show how GPUs can quickly generate random imaging noise. Our code is available as a small, easytouse library^{1}^{1}1https://github.com/bbrister/cudaImageWarp written in Python, C++ and CUDA, which is 2.68faster than a SciPy implementation, including communication overhead [1].
Besides the input data, we found that the choice of loss function is important for achieving high accuracy in unbalanced segmentation tasks. In our experiments, data augmentation always decreases the validation loss, but this only corresponds to an improvement in binary segmentation accuracy for certain loss functions. Weighted crossentropy, the most common loss function for segmentation, performs poorly when the classes are significantly unbalanced, as is often the case in medical images. This is because it assigns a lesser penalty for false positives than for false negatives. Instead, we propose the intersectionoverunion (IOU) loss function, which is closely related to the Dice loss, but has the advantage of being a metric [2, 3]. We analyze the mathematical properties of the IOU and Dice losses, showing that they have approximately balanced penalties for each type of error, and that each belongs to an infinite family of functions having essentially the same properties. In our experiments, the IOU loss significantly outperforms weighted crossentropy, and results are further improved by data augmentation.
Finally, we created a new dataset and applied these technologies to CT organ segmentation, which is the task of detecting and delineating organs in a CT scan. Manual annotation of large datasets is impractical for complex organs, such as the skeleton. To address this challenge, we show how to automatically generate training labels using simple image morphology, which is accelerated by 3D Fourier transforms. We automatically generate labels for lungs and bone, which are then combined with manual annotations for the liver, kidneys and bladder. Our experiments suggest that the FCN is able to learn the “average appearance” of the organ, and is not distracted by intermittent failures in the crude morphological algorithm. In addition to fullyautomated organ labels, we manually delineated the liver, kidneys and bladder in a dataset of 130 CT scans from the Liver Tumor Segmentation Challenge (LiTS), to create a free organ segmentation dataset with all five organ classes [4]. We then used this dataset to train a multiclass FCN. Our model segments the lung, liver, bone, kidneys and bladder with respective Dice scores of 91.8%, 92.2%, 79.3%, 78.3% and 83.7%, consistent with results from the literature which depend on more complex models [5, 6]. We also evaluated the liver segmentation on the LiTS challenge test set, with similar results. The trained model, code and dataset will be made publicly available, to further medical imaging research.
Our main contributions, along with the structure of the paper, are summarized as follows:

Unsupervised algorithms for labeling the bones and lungs in CT scans (section III)

An efficient GPU implementation of 3D data augmentation (section IV)

A new loss function which improves results for unbalanced segmentation (section VC)

A new dataset consisting of five organs labeled in 130 CT scans (section VIA)
Ii Related work
This section describes the related work supporting each of our main contributions–first organ segmentation in general, then using deep learning, then data augmentation, then weak supervision, and finally loss functions.
Organ segmentation has long been an active area of medical imaging research. Approaches can be divided into two main categories: semiautomatic and fullyautomatic. The former group requires a usergenerated starting shape, which grows or shrinks into the organ of interest. These approaches typically take into account intensity distributions, sharp edges and even the shape of the target organ, and their success depends greatly on the initialization. See Freiman et al. for an example of semiautomatic liver segmentation [7].
The other category, fullyautomatic methods, require no user input, as they directly detect the object of interest in addition to delineating its boundaries. These techniques can be divided into two main areas, pattern recognition and atlasbased. Pattern recognition methods are currently dominated by artificial neural networks such as FCNs, but other other deep learning methods have been tried
[8]. For example, Qin et al. used a CNN to classify superpixels, while Dong et al. combined an FCN with a generative adversarial network (GAN)
[9, 10]. In contrast to pattern recognition, atlasbased methods work by warping, or registering an image to a labeled atlas image. While atlasbased methods can achieve high accuracy, interpatient registration is computationally expensive and extremely sensitive to changes in imaging conditions, for example the absence of an organ. These methods are best suited to stationary objects of consistent size and shape, such as the brain [11]. For wholebody imaging, pattern recognition methods, and in particular deep learning, have rapidly grown in popularity. For the sake of completeness, it is worth noting that the distinction between these categories can be vague. For example, Okada et al. proposed a method which automatically generates seeds for semiautomatic segmentation methods, based on the size and location of a userprovided liver segmentation [12].As in our work, many researchers are now applying FCNs to the task of CT organ segmentation[5, 6, 3, 13]. However, most of these models are limited to a specific organ or body region, around which they assume the images will be cropped. In addition, some multiorgan projects are based on privately held data, and thus can neither be replicated nor extended beyond the obvious limitations of finetuning [5]. In contrast, our simple model naturally applies to a wide variety of organs, and requires no manual intervention in preparing the images. Leveraging a memoryefficient model and resampling the images to a coarse mm
resolution allows us to process large sections of the body at once. Using Tensorflow, we created a computationallyefficient model which can be deployed on a wide variety of devices, for use by practitioners of other disciplines
[14].Data augmentation has become standard practice in applying deep learning models to images. As in our work, several others have applied random affine transformations to 3D volumes [3]
. However, to our knowledge, we are the first to develop a computationally efficient system for augmenting 3D images using GPU texture sampling, and the first to leverage GPUs for random noise generation for the purpose of data augmentation. Some of the current deep learning frameworks, such as Keras, provide builtin operations for image manipulation
[15]. However, at the time of writing these do not support 3D geometric operations, and they are probably not implemented in the most efficient way.
Automatic generation of training labels is an active area of machine learning research, sometimes called “weak supervision.” Ghafoorian et al. used a regiongrowing algorithm to train a deep neural network for brain ventricle segmentation
[16]. We apply an even simpler method to the lungs and bone, which is accelerated via 3D Fourier transforms, and unlike the previous work, requires no user input. This is essential for bones, which are too numerous to annotate manually.The final contribution of our work is the IOU loss, both its proposal and mathematical analysis. This is closely related to the Dice loss which was proposed by Milltari et. al [3]. Sudre et al. noted that the Dice loss, as well as a variant based on fuzzy set theory, outperform weighted crossentropy for unbalanced classes [17]. However, to our knowledge, no one has yet expounded on the mathematical properties of these loss functions, their relationship to the binary scores for which they are named, nor the reason for their superior performance to the usual pervoxel loss functions [17].
Iii Training label generation
This section explains how we used image morphology to automatically detect and segment organs, which generates training labels for our neural network. First we describe the basics of dimensional image morphology, and how we accelerated these operations using Fourier transforms. Then we describe the specific algorithms used to segment the lungs and bones.
Iiia Morphology basics, acceleration by Fourier transforms
Let denote a binary image, and the structuring element. Then dilation is
That is, we first convolve with , treating the two as realvalued functions on . Then, we convert back to a binary image by setting zerovalued pixels to black, all others to white. Erosion is computed similarly: if is the binary complement of , then erosion is just . Similarly, the opening and closing operations are compositions of erosion and dilation.
Written this way, all of the basic operations in dimensional binary morphology reduce to a mixture of complements and convolutions, the latter of which can be computed by fast Fourier transforms, due to the identity , where denotes the Fourier transform of . This allows us to quickly generate training labels by morphological operations, especially when the structuring elements are large. In what follows, we describe how morphology is used to generate labels for two organs of interest, the skeleton and lungs. These organs were chosen because their size and intensity contrast enable detection by simple thresholding.
IiiB Morphological detection and segmentation of CT lungs
The lungs were detected and segmented based on the simple observation that they are the two largest air pockets in the body. First, we extract the air pockets from the CT scan by removing all voxels greater than 150 Hounsfield units (HU). The resulting mask is called the thresholded image. Then, we remove small air pockets by morphologically eroding the image using a spherical structuring element with a diameter of 1 cm. Next, we remove any air pockets which are connected to the boundary of any axial slice in the image^{2}^{2}2That is, the maximal and minimal extent of the image in the and axes.. This removes air outside of the body, while preserving the lungs. From the remaining air pockets, we take the two largest connected components, which are almost certainly the lungs. Finally, we undo the effect of erosion by taking the components of the thresholded image which are connected to the two detected lungs. An example output in shown in figure 1.
A quick web search reveals that similar algorithms have previously been used to segment the lungs [19]. The novelty in our approach lies in using this algorithm to train a deep neural network, which can be more robust than the morphological algorithm from which it was trained, as discussed in section VIC.
IiiC Morphological detection and segmentation of CT skeleton
Bone segmentation proceeds similarly to lung segmentation, by a combination of thresholding, morphology and selection of the largest connected components. This time we define two intensity thresholds, and HU. These were selected so that almost all bone tissue is greater than , while the hard exterior of each bone is usually greater than . First, we select the exteriors of all the bones by thresholding the image by . This step inevitably also includes some unwanted tissues, such as the aorta, kidneys and intestines, especially in contrastenhanced CTs. To remove these unwanted tissues, we select only the largest connected component, which is the skeleton. Next, we fill gaps in the exteriors of the bones by morphological closing, using a spherical structuring element with a diameter of 2.5 cm. This step has the unwanted side effect of filling gaps between bones as well, so we apply the threshold to remove most of this unwanted tissue.
At this stage, there could be holes in the center of large bones, such as the pelvis and femurs. To fill these, we note that, when the patient is reclined on the exam table, large bones almost always lie parallel to the axis of the image. Accordingly, we process each plane (axial slice) in the image, filling in any holes which are not connected to the boundaries. This fills in the centers of large bones in most slices, which suffices for our purposes of training a deep neural network. See figure 2 for a 3D visualization of the resulting skeleton segmentation. The accuracy of this scheme is evaluated in section VIC.
Iv GPUaccelerated data augmentation
Deep neural networks can be viewed as tabulae rasae on which a wide variety of classifiers can be inscribed, depending on the training data. For example, we typically expect that image classifiers should be invariant to rotation and scaling, but the basic units of convolutional neural networks (CNNs) need not be. Given a basic dataset, data augmentation applies random transformations to generate new training data. Theoretically, this can be justified as modifying the data distribution, which alters the expected loss for each choice of parameters. The expanded dataset also helps to combat overfitting.
Our organ segmenter was trained using a variety of data transformations including affine warping, intensity windowing, additive noise and occlusion. These common operations are expensive to compute over large 3D volumes. In particular, affine warping requires random access to a large buffer of image data, with little reuse, which is inefficient for the cacheheavy memory hierarchy of CPUs. A typical CT scan consists of hundreds of slices of 12bit data. When arranged into a 3D volume, a CT scan is hundreds of times larger than a typical lowresolution photograph used in conventional computer vision applications.
Although 3D data augmentation is slow on conventional CPUs, the operations involved are actually very efficient when implemented on more specialized hardware. All of the aforementioned operations are common in computer graphics, so GPUs have been designed to handle them efficiently. In particular, GPU texture memory is optimized for parallel access to 2D or 3D images, and includes special hardware for handing outofbounds coordinates and interpolation between neighboring pixels. This is especially wellsuited to accelerating affine warping. Photometric operations such as noise generation, windowing and cropping also benefit from the massive parallelism offered by GPUs.
Since these operations involve little reuse of data, each output pixel is drawn by its own CUDA thread. Each thread computes the affine coordinate map, samples the input image at that coordinate, applies the photometric transformations, then writes the final intensity value to the output volume. In this way, each output requires only a single access to texture memory.. To mitigate the cost of transferring data between memory systems, we designed our data augmentation library with a firstin firstout (FIFO) queue to pipeline jobs. The concept is illustrated in figure 3. While one image is being processed, the next has already begun transferring from main memory to graphics memory, effectively hiding its transfer latency. The FIFO programming model naturally matches the intended use case of augmenting an entire training batch at once. Our experiments in section VIB show this scheme accelerates processing by .
In what follows, we describe our image manipulation pipeline, with the details of each operation.
First, we warp the image using a random affine transformation. Sampling coordinates are computed as , where and . The matrix is generated by composing a variety of geometric transformations drawn uniformly from userspecified ranges, including rotation, scaling, shearing, reflection and generic affine warping^{3}^{3}3For reflection, the user specifies the probability that the image is reflected about each of the three axes.. Finally, a random displacement
is drawn from a uniform distribution according to userspecified ranges. Taking
to be the coordinates of the center of the volume, we then compute according to the formula , which guarantees . That is, the center of the image is displaced by units. The output image is defined by , where denotes the input image volume from the training set. The discrete image data is sampled from texture memory using trilinear interpolation, whereas the labels are sampled according to the nearest neighboring voxel.Second, we apply random occlusion, also known as cropping. To encourage robustness to missing anatomy, we randomly set part of the input volume to zero, while ensuring that every voxel has an equal chance of being occluded. The occluded region is a rectangular prism, with height uniformly distributed in the interval, and starting coordinate , where is the maximum possible coordinate. Then, we compute
Since we are already applying an affine transformation to the image, removing an axisaligned prism from the output effectively removes a randomlyoriented prism from the input. For efficiency, we evaluate occlusion prior to sampling the image texture. If the value is negative, all future operations are skipped, including the texture fetch.
Third, we introduce additive Gaussian noise as a simplistic model of the artifacts introduced in image acquisition. The operation is simply , where
is drawn from an independent, identically distributed (IID) Gaussian process with zero mean and standard deviation
. The sole parameter is drawn from a uniform distribution for each training example. In this way, some images are severely corrupted by noise, while others are hardly changed. We used the cuRand library to quickly generate noise on the GPU [20]. We initialize a separate random number generator (RNG) for each GPU thread, with one thread per output voxel. To reduce the initialization overhead, each thread uses a copy of the same RNG, starting at a different seed. This sacrifices the guarantee of independence between RNGs, but the effect is not perceivable in practice.Finally, we apply a random window/level transformation to the image intensities. To increase image contrast, radiologists always view CT scans within a certain range of Hounsfield units (HU). For example, bones might be viewed with a window of 10001500 HU, while abdominal organs would be viewed with a narrower window of 150230 HU. To simulate this, we randomly draw limits according to userspecified ranges. Then we compute
In other words, the intensity values are clamped to the range , and then affinely mapped to . This is straightforward to implement on GPUs.
V Neural network architecture
This section describes the design of our predictive model; both the inputs and outputs, pre and postprocessing, the neural network itself, and the training loss function.
Va Pre and postprocessing
Our neural network takes as input a image volume, and outputs a
probability map, where each voxel is assigned a class probability distribution. This becomes a
prediction map by taking the probability for each voxel. To reduce memory requirements, we resample all image volumes to mm before they are fed into the model. Resampling consists of Gaussian smoothing, which serves as a lowpass filter to avoid aliasing artifacts, followed by interpolation at the new resolution. Since each CT scan has its own millimeter resolution for each dimension , we adjust the Gaussian smoothing kernel according to the formula where the smoothing factors are computed from the desired resolution according toThis heuristic formula is based on the fact from digital signal processing that, in order to avoid aliasing, the cutoff frequency should be placed at
, the ratio of sampling rates, on a frequency scale.After the neural network, we resample the prediction map to the original image resolution by nearest neighbor interpolation. One difficulty with this scheme is that CT scans vary in resolution and number of slices, and at mm we are unlikely to fit the whole scan in our network. For training, we address this by selecting a subregion from the scan uniformly at random. For inference, we cover the scan by partiallyoverlapping subregions, averaging predictions where overlap occurs. Others have addressed the issue of limited viewing resolution by training multiple FCNs, each working at different scales [5]. However, we found that a single 3 mm network achieves competitive performance.
VB Model
We chose a relatively simple model which balances speed and memory consumption with accuracy. Our model is based on GoogLeNet, but its convolution and pooling operators work in 3D instead of 2D [21]. Like other fullyconnected networks (FCNs), our model essentially consists of two parts: one for decimation and one for interpolation [8]
. The decimation network is similar to the usual convolutional neural network (CNN) used for image classification, with maxpooling and strided convolution layers which reduce the image size by a factor of 8 in each dimension. The interpolation layer restores the feature maps to their image dimensions, essentially through convolutional interpolation. Unlike most of the literature, we chose not to use any “skip connections” which forward feature maps in the decimation part to later layers in the interpolation part
[8, 3, 6]. We did this two reasons: first, to save memory which is precious when dealing with 3D models. Second, to disambiguate the advantages of our loss function and data augmentation from any improvements in model architecture.Purpose  Type  Filter size  Outputs  Stride 

Decimation  Convolution  7  64  2 
Convolution  1  64  1  
Convolution  4  192  1  
Pooling  2  192  2  
Inception  15  256  1  
Inception  15  480  1  
Inception  15  512  1  
Inception  15  512  1  
Inception  15  528  1  
Inception  15  832  1  
Pooling  2  832  2  
Inception  15  832  1  
Inception  15  1024  1  
Interpolation  Interpolation  8  6  1/8 
Output  Softmax  n/a  6  1 
Table I lists the layers of our neural network, in order from the input image data to the final probability maps. The exact details of each layer type are beyond the scope of this paper, but should be familiar to practitioners. Filter sizes and strides apply to all three dimensions, for example a filter size of 7 implies a
isotropic filter. All convolutions are followed by constant “bias” addition, batch normalization and rectification
[22]. Pooling always refers to taking neighborhood maxima. An inception module consists of a multitude of convolution layers of sizes 1, 3 and 5, along with a pooling layer, which are concatenated to form four heterogeneous output paths
[21]. The inception module seems to be a memoryefficient way to construct very deep neural networks, since it employs relatively inexpensive operations of heterogeneous sizes. For simplicity, we report the total number of outputs of the inception module, rather than the number of filters of each type. The final softmax layer outputs class probabilities for each voxel.
VC IOU Loss
The most common loss function used for image segmentation is weighted crossentropy. This method assigns a separate loss function to each voxel, minimizing the weighted average of all losses. To compensate for extreme class imbalances encountered in medical imaging, a separate weight is assigned to each voxel so that all objects weigh the same, regardless of size. We found this apporach to give poor results, and instead formulated a different loss function, the IOU loss.
Imagine a binary classification scenario, as shown in figure 4. Let
denote the vector of
output probabilities from the model, where is the number of voxels, let denote the binary ground truth labels, and let denote the weight of the class . Let denote the crossentropy loss function for voxel . Then weighted crossentropy loss isThe problem with this loss function is that it weights errors unequally: a false positive receives the weight , whereas a false negative receives . As a result, the model learns to exaggerate the boundaries of small objects, since false positives almost always receive the weight of the large “background” class. This is a limitation not just of crossentropy, but of any loss function which is a weighted average over each voxel.
In order to address this limitation, we adopted a loss function which minimizes the intersection over union (IOU) of the output segmentation, also called the Jaccard index. Since the IOU depends on binary values, it cannot be optimized by gradient descent, and thus it is not directly suitable as a loss function for deep learning. Instead, our strategy is to define a smooth function which is equal to the IOU when the output probabilities are all either
or .Since , we can consider as a vector in consisting only of probabilities and . Then, the smooth IOU loss is
where , since for all . To extend this to multiclass segmentation, simply average the for each class. This approach seems simpler and more naturally motivated than the multiclass Dice scheme of Sudre et al. [17].
The IOU loss is closely related to the Dice loss, which was unceremoniously proposed by Milletari et al. [3]. The Dice loss is
Since the IOU loss is similar to the Dice loss, one might wonder why we bothered to propose a new loss function. One advantage is that, when restricted to binary scores, the IOU loss obeys the triangle inequality, and thus qualifies as a metric [2]. This is not true of the Dice loss. For example, let and . Then it can be verified that , which violates the triangle inequality.
One useful application of the triangle inequality is to bound the amount by which we can decrease the IOU loss by restricting the experiment to a subset of the domain. For example, for any set of binary predictions and labels , we can define a restriction by the taking the intersection with a binary vector . According to the triangle inequality,
(1) 
so the loss decrease is bounded by , which simplifies to .
To our knowledge, no one has yet described the mathematical properties of the Dice loss in detail. Since it has essentially the same properties as the IOU loss, both being continuous interpolations of binary or set functions, we treat both losses simultaneously. They have the following properties, which are easily verified:

They are equal to the desired binary loss, either Dice or IOU, when is binary.

They are strictly increasing in each when , decreasing when .

They are maximized only when , minimized only when .

They are smooth functions, if we define the loss to be at , which is otherwise undefined.
Properties 13 ensure that minimizing the continuous loss corresponds to maximizing the binary score of the trained model, while properties 24 encourage the training optimization problem to be wellbehaved. For example, property 2 implies that the function has no strict local minima, since these would also be local minima when restricted to a single input probability . However, this property does not extend to an average of perimage loss functions over a dataset, since this could be a sum of both increasing and decreasing functions, which need not be monotone.
Importantly, these properties do not suffice for uniqueness of the loss function, since there are other possible behaviors when is nonbinary. In fact, these properties are not even unique among the rational functions. For example, consider the family of loss functions
for any power . It is easily verified that these functions also satisfy all of the above properties, but differ from the case for nonbinary values of . More generally, if is a collection of smooth increasing functions on with and for all , then the function
satisfies properties 14. This shows that there is nothing like a unique choice of “IOU loss” or “Dice loss,” unless deeper properties are discovered. However, we can at least say that the proposed loss functions are the lowestorder rational functions satisfying our properties, and thus the most computationally efficient.
We now analyze the penalty of each type of error, false positives and false negatives, for the IOU loss. Let denote the misprediction error, and let . A false negative corresponds to for all , whence the loss is
(2)  
Similarly, a false positive has , which yields
(3)  
So long as , the two penalties are approximately equal. This shows that, unlike weighted crossentropy, the IOU loss strikes a reasonable balance between each type of error. The assumption essentially means that the model is performing well on the training data, which we should expect towards the end of training.
Vi Experimental Results
This section describes our experiments validating the methods on real data. First describe how we created our CT organ dataset. Then, we measure the speed of our GPU data augmentation against a CPU baseline. Finally, we train various organ segmentation networks, evaluating their performance alongside the unsupervised label generators, first on our own dataset, and then on a public challenge.
Via Dataset
In order to train and evaluate our organ segmentation system, we annotated a set of 130 anonymized CT volumes exhibiting a wide variety of imaging conditions, both with and without contrast, abdominal, thoracic and fullbody. The original images came from the Liver Tumor Segmentation (LiTS) Challenge organized by Christ et al. [4]. The LiTS challenge provides liver masks for almost all of the data, which were extracted using a semiautomated segmentation method. To complete the dataset, we added a few liver masks ourselves, and cropped some of the image volumes to eliminate undesired artifacts. To this we added our own masks for the kidneys and bladders. All of these annotations were created using ITKSNAP, a free tool for volumetric image segmentation, using a mixture of active contours and manual correction [23].
Since training labels were automatically generated for the bones and lungs, we must evaluate the model performance on a different dataset. However, manually labeling these organs is extremely tedious, which was the original motivation for developing the automatic labelers. As a compromise, we manually labeled the lungs and bones in 10 out of the 130 CT scans, again using ITKSNAP. For the bones, we saved time by starting with the automatic labels, and then manually correcting the errors and omissions.
ViB Data augmentation speed
After developing our GPU data augmentation program, we implemented the same operations using SciPy, a popular library for scientific computing [1]. We then compared the speed of the two programs on our CT organ dataset. In order to maintain consistency with our organ segmentation experiment, all images and labels were resampled to a resolution of mm. To remove the effects of file I/O, we wrote the whole dataset into main memory at the start of each experiment. Execution times were averaged over 5 different image batches. Our data augmentation code was written in C++ and CUDA and controlled by a Python wrapper. All experiments were performed on a single machine with an Intel Core i76900K CPU and four NVIDIA GeForce Titan X Pascal GPUs.
The average time per image volume is reported in figure 5. The GPU implementation utilizes the memory transfer scheme from section IV, so it benefits from increased batch size which hides the communication overhead. Each batch is evenly distributed across the four GPUs in our server, which allows for a larger total batch size, albeit at marginal benefit over a single GPU with a batch size of 32. On the other hand, the CPU grows slower with increasing batch size, which is probably due to its multilevel cachebased memory hierarchy. Execution times range from 12574 ms per CT scan on the GPU to 323592 ms on the CPU. Accordingly, the GPU offers a speedup of 2.68.1, depending on the batch size. The largest GPU batch size is faster per volume than the smallest, due to job pipelining. Comparing the fastest times for each processor results in a speedup for the GPU. All GPU times include the overhead of sending images to and from graphics memory, which could theoretically be elided if deep learning frameworks supported direct access to CUDA objects created by other programs. This experiment shows that GPUs offer a considerable speedup over generic CPUs for 3D data augmentation.
ViC CT organ segmentation accuracy
Method  Neural Nework  Morphology  

Loss  Cross entropy  IOU  n/a  
Data augmentation  No  Yes  No  Yes  n/a 
Lung  87.5  86.8  90.8  91.8  97.8 
Liver  88.8  85.5  90.9  92.2  n/a 
Bone  73.6  65.4  78.9  79.3  93.2 
Kidney  71.8  65.2  77.5  78.3  n/a 
Bladder  70.4  58.4  80.1  83.7  n/a 
Other  97.7  96.4  98.6  98.7  99.7 
Includes all organs besides lung and bone
In order to compare our proposed methods to the default variant, we trained the same neural network with two different loss functions, crossentropy and IOU. In each case, we began without data augmentation, only introducing it after the training loss ceased to decline. This strategy saves time, since the network first learns the basic appearance of each organ before adapting to the various augmentations. The learning curves are shown in figure 6. Each training iteration consists of a batch of 4 image volumes. The training loss is averaged over all the random crops in 100 training iterations, while the validation loss is averaged over all subvolumes of each image in the validation set. In each case, training loss spikes upward when augmentation is introduced, and then slowly settles back down. Interestingly, data augmentation always increases the training loss, as the training data becomes more difficult, but decreases the loss on the unmodified validation data.
The final results are shown in table II. For all organs, the IOU loss with data augmentation outperforms all other variants. Even without augmentation, the IOU loss outperforms both variants of crossentropy. Surprisingly, data augmentation increases validation accuracy with the IOU loss, but decreases the accuracy with cross entropy. This occurs even though the validation crossentropy loss is decreasing, as shown in figure 6. This illustrates the issue described in section VC, where weighted crossentropy encourages false positives over false negatives for minority classes, so the model tends to overestimate organ boundaries. This is especially evident with the bladder, as shown in figure 7. The effect is less pronounced for crossentropy without data augmentation, since the model can overfit the small training set.
On our dataset, the neural network failed to match the accuracy of the unsupervised morphological algorithm for bones and lungs. For the largest and most distinctive organs, image morphology is simpler and probably more precise than a neural network. However, the two approaches are not mutually exclusive, as unsupervised methods can postprocess the output of a neural network. The main disadvantage of morphology is that it is prone to catastrophic failure in ways which are difficult to anticipate. For example, the morphological lung segmenter can be thrown off by other airfilled objects in the scan, such as an exam table. In contrast, a neural network captures the visual appearance of each organ, which is often more reliable. Finally, when segmenting a large variety of organs, a single network has obvious conceptual and computational advantages over a litany of brittle, organspecific solutions.
ViD Liver segmentation challenge
Seeking external validation of our previous results, we ran the four models on the LiTS challenge test dataset, without any additional training for the liveronly task [4]. Our results are shown in table III. Our scores are similar to the liver scores in table II, where IOU loss with data augmentation outperforms all the other models.
Our best model achieved a mean Dice per case of 90.5%, which at the time of writing would place us in rank 49 on that section of the challenge. The liver segmentation leaderboard is highly competitive, as the top score is 96.6% out of a possible 100%. Differences between such high scores are likely attributable to subtle discrepancies in organ boundaries, to which most applications are insensitive. Our model is considerably simpler and more general than the topscoring methods, as it operates at a coarse 3mm
resolution, uses no additional data, and identifies five different organs simultaneously. Our averge processing time was 5.98s per volume, consuming a modest 2429 MB of graphics memory, which is well within the capabilities of commodity graphics cards. In contrast, other submissions focused only on the liver, and leveraged significantly more complex and expensive methods, including higher processing resolution, cascades of 2D and 3D models, sophisticated postprocessing, and even transfer learning from other datasets
[13]. Our work shows that, with the right training, a relatively simple model suffices for many applications. Furthermore, nothing precludes the topperforming models from incorporating our suggested improvements.Loss  Cross entropy  IOU  

Data augmentation  No  Yes  No  Yes 
Avg. Dice  0.866  0.840  0.896  0.905 
Global Dice  0.872  0.846  0.904  0.911 
VOE  0.233  0.273  0.185  0.168 
RVD  0.193  0.248  0.024  0.033 
ASSD  8.506  10.407  5.876  3.767 
MSSD  0.767  173.722  88.843  54.832 
RMSD  18.299  21.908  13.533  8.380 
Vii Conclusion
We delineated a dataset of 130 abdominal CT scans using a mixture of manual annotation and automated morphological segmentation. We used this dataset to train a deep neural network for CT organ segmentation, using GPUs to accelerate data augmentation. Our model uses the IOU loss to improve segmentation accuracy. We explained mathematically why there is no unique IOU loss, and why the various IOU loss functions outperform weighted crossentropy for unbalanced segmentation tasks. The code, data and trained model will be made publicly available. We hope that the dataset will enable the development and evaluation of more accurate organ segmentation methods. We invite others to annotate new organs, which could easily be incorporated into the existing system.
Acknowledgments
This work was supported in part by grants from the National Cancer Institute, National Institutes of Health, 1U01CA190214 and 1U01CA187947.
References
 [1] E. Jones, T. Oliphant, P. Peterson, et al., “SciPy: Open source scientific tools for Python.” http://www.scipy.org/, 2001–. Accessed October 15, 2018.
 [2] H. Späth, “The minisum location problem for the jaccard metric,” OperationsResearchSpektrum, vol. 3, pp. 91–94, Jun 1981.
 [3] F. Milletari, N. Navab, and S. Ahmadi, “Vnet: Fully convolutional neural networks for volumetric medical image segmentation,” in 2016 Fourth International Conference on 3D Vision (3DV), pp. 565–571, Oct 2016.
 [4] P. Christ, J. Holch, C. Jacobs, G. Chartrand, A. BenCohen, J. Moreau, R. Vivanti, E. Carmichael, et al., “LiTS  liver tumor segmentation challenge.” https://competitions.codalab.org/competitions/17094, 2017. Accessed October 13, 2018.
 [5] H. R. Roth, H. Oda, X. Zhou, N. Shimizu, Y. Yang, Y. Hayashi, M. Oda, M. Fujiwara, K. Misawa, and K. Mori, “An application of cascaded 3d fully convolutional networks for medical image segmentation,” Computerized Medical Imaging and Graphics, vol. 66, pp. 90 – 99, 2018.
 [6] E. Gibson, F. Giganti, Y. Hu, E. Bonmati, S. Bandula, K. Gurusamy, B. Davidson, S. P. Pereira, M. J. Clarkson, and D. C. Barratt, “Automatic multiorgan segmentation on abdominal ct with dense vnetworks,” IEEE Transactions on Medical Imaging, vol. 37, pp. 1822–1834, Aug 2018.
 [7] M. Freiman, O. Eliassaf, Y. Taieb, L. Joskowicz, Y. Azraq, and J. Sosna, “An iterative bayesian approach for nearly automatic liver segmentation: algorithm and validation,” International Journal of Computer Assisted Radiology and Surgery, vol. 3, p. 439, Jul 2008.
 [8] E. Shelhamer, J. Long, and T. Darrell, “Fully convolutional networks for semantic segmentation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 39, pp. 640–651, April 2017.
 [9] W. Qin, J. Wu, F. Han, Y. Yuan, W. Zhao, B. Ibragimov, J. Gu, and L. Xing, “Superpixelbased and boundarysensitive convolutional neural network for automated liver segmentation,” Physics in Medicine & Biology, vol. 63, no. 9, p. 095017, 2018.
 [10] D. Yang, D. Xu, S. K. Zhou, B. Georgescu, M. Chen, S. Grbic, D. Metaxas, and D. Comaniciu, “Automatic liver segmentation using an adversarial imagetoimage network,” in Medical Image Computing and ComputerAssisted Intervention (MICCAI 2017) (M. Descoteaux, L. MaierHein, A. Franz, P. Jannin, D. L. Collins, and S. Duchesne, eds.), (Cham), Springer International Publishing, 2017.
 [11] P. Kalavathi and V. B. S. Prasath, “Methods on skull stripping of mri head scan images—a review,” Journal of Digital Imaging, vol. 29, pp. 365–379, Jun 2016.
 [12] T. Okada, M. G. Linguraru, M. Hori, R. M. Summers, N. Tomiyama, and Y. Sato, “Abdominal multiorgan segmentation from ct images using conditional shapelocation and unsupervised intensity priors,” Medical Image Analysis, vol. 26, no. 1, pp. 1 – 18, 2015.
 [13] X. Li, H. Chen, X. Qi, Q. Dou, C. Fu, and P. Heng, “Hdenseunet: Hybrid densely connected unet for liver and tumor segmentation from ct volumes,” IEEE Transactions on Medical Imaging, pp. 1–1, 2018.
 [14] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, “TensorFlow: Largescale machine learning on heterogeneous systems.” https://www.tensorflow.org/, 2015.
 [15] F. Chollet et al., “Keras.” https://keras.io, 2015.
 [16] M. Ghafoorian, J. Teuwen, R. Manniesing, F.E. de Leeuw, B. van Ginneken, N. Karssemeijer, and B. Platel, “Student beats the teacher: deep neural networks for lateral ventricles segmentation in brain mr,” 2018.
 [17] C. H. Sudre, W. Li, T. Vercauteren, S. Ourselin, and M. Jorge Cardoso, “Generalised dice overlap as a deep learning loss function for highly unbalanced segmentations,” in Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support, pp. 240–248, Springer International Publishing, 2017.
 [18] R. Kikinis, S. D. Pieper, and K. G. Vosburgh, 3D Slicer: A Platform for SubjectSpecific Image Analysis, Visualization, and Clinical Support, pp. 277–289. New York, NY: Springer New York, 2014.
 [19] Mathworks, “Segment lungs from 3d chest scan and calculate lung volume.” https://www.mathworks.com/help/images/segmentlungsfrom3dchestmridata.html, 2018. Accessed October 13, 2018.
 [20] NVIDIA, “cuRAND.” "https://developer.nvidia.com/curand", 2018. Accessed November 8, 2018.
 [21] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “Going deeper with convolutions,” in The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2015.
 [22] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” in Proceedings of the 32nd International Conference on Machine Learning (F. Bach and D. Blei, eds.), vol. 37 of Proceedings of Machine Learning Research, (Lille, France), pp. 448–456, PMLR, Jul 2015.
 [23] P. A. Yushkevich, J. Piven, H. Cody Hazlett, R. Gimpel Smith, S. Ho, J. C. Gee, and G. Gerig, “Userguided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability,” Neuroimage, vol. 31, no. 3, pp. 1116–1128, 2006.
Comments
There are no comments yet.