Nov 27, 2018 - volumetric image processing. Additionally .... by warping, or registering an image to a labeled atlas image. ... Using Tensorflow, we c...

0 downloads 4 Views 504KB Size

arXiv:1811.11226v1 [cs.NE] 27 Nov 2018

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible.

2

CT organ segmentation using GPU data augmentation, unsupervised labels and IOU loss Blaine Rister, Darvin Yi, Kaushik Shivakumar, Tomomi Nobashi and Daniel L. Rubin

Abstract—Fully-convolutional neural networks have achieved superior performance in a variety of image segmentation tasks. However, their training requires laborious manual annotation of large datasets, as well as acceleration by parallel processors with high-bandwidth memory, such as GPUs. We show that simple models can achieve competitive accuracy for organ segmentation on CT images when trained with extensive data augmentation, which leverages existing graphics hardware to quickly apply geometric and photometric transformations to 3D image data. On 3 mm3 CT volumes, our GPU implementation is 2.6-8×faster than a widely-used CPU version, including communication overhead. We also show how to automatically generate training labels using rudimentary morphological operations, which are efficiently computed by 3D Fourier transforms. We combined fully-automatic labels for the lungs and bone with semi-automatic ones for the liver, kidneys and bladder, to create a dataset of 130 labeled CT scans. To achieve the best results from data augmentation, our model uses the intersection-over-union (IOU) loss function, a close relative of the Dice loss. We discuss its mathematical properties and explain why it outperforms the usual weighted cross-entropy loss for unbalanced segmentation tasks. We conclude that there is no unique IOU loss function, as the naive one belongs to a broad family of functions with the same essential properties. When combining data augmentation with the IOU loss, our model achieves a Dice score of 78-92% for each organ. The trained model, code and dataset will be made publicly available, to further medical imaging research. Index Terms—computer vision, image segmentation, deep learning, data augmentation, computed tomography (CT)

I. I NTRODUCTION Fully-convolutional 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 application-specific 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, memory-efficient model trained with a new loss Manuscript received 11/27/2018. The authors are with the Departments of Electrical Engineering, Radiology and Biomedical Data Science at Stanford University, Stanford, CA USA 94305 (email: [email protected], [email protected], [email protected], [email protected], [email protected]).

function which excels at unbalanced segmentation. The end result is a fast and robust model for multi-organ 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 high-bandwidth 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, easy-to-use library1 written in Python, C++ and CUDA, which is 2.6-8×faster 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 cross-entropy, 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 intersection-over-union (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 cross-entropy, 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 1 https://github.com/bbrister/cudaImageWarp

3

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 fully-automated 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 multi-class 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 V-C) • A new dataset consisting of five organs labeled in 130 CT scans (section VI-A) II. R ELATED 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: semi-automatic and fully-automatic. The former group requires a user-generated 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, fully-automatic 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 atlas-based. 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 super-pixels, while Dong et al. combined an FCN with a generative adversarial network (GAN) [9], [10]. In contrast to pattern recognition, atlas-based methods work by warping, or registering an image to a labeled atlas image. While atlas-based 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 whole-body 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 semi-automatic segmentation methods, based on the size and location of a user-provided liver segmentation [12]. As in our work, many researchers are now applying FCNs to the task of CT organ segmentation [3], [5], [6], [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 multi-organ projects are based on privately held data, and thus can neither be replicated nor extended beyond the obvious limitations of fine-tuning [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 memory-efficient model and resampling the images to a coarse 3 mm3 resolution allows us to process large sections of the body at once. Using Tensorflow, we created a computationally-efficient 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 built-in 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 region-growing 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 cross-entropy 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 per-voxel loss functions [17].

4

III. T RAINING 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 n-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. A. Morphology basics, acceleration by Fourier transforms Let f : Zn → F2 denote a binary image, and k : Zn → F2 the structuring element. Then dilation is ( 1, (f ∗ k)(x) > 0 D(f, k)(x) = (1) 0, otherwise. That is, we first convolve f with k, treating the two as real-valued functions on Zn . Then, we convert back to a binary image by setting zero-valued pixels to black, all others to white. Erosion is computed similarly: if f¯ is the binary complement of f , then erosion is just E(f, k)(x) = D(f¯, k). Similarly, the opening and closing operations are compositions of erosion and dilation. Written this way, all of the basic operations in ndimensional 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 f[ ∗k = fb· b k, where fb denotes the Fourier transform of f . 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. B. 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 image2 . 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 VI-C. 2 That

is, the maximal and minimal extent of the image in the x and y axes.

Figure 1. Example of CT lung detection and segmentation using image morphology. Lung mask overlaid in blue. Rendered by 3D Slicer [18].

C. 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, τ1 = 0 and τ2 = 200 HU. These were selected so that almost all bone tissue is greater than τ1 , while the hard exterior of each bone is usually greater than τ2 . First, we select the exteriors of all the bones by thresholding the image by τ2 . This step inevitably also includes some unwanted tissues, such as the aorta, kidneys and intestines, especially in contrast-enhanced 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 τ1 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 z-axis of the image. Accordingly, we process each xy-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 VI-C. IV. GPU- ACCELERATED

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

5

Figure 3. Pipelining jobs for efficient batch processing. With this scheme, a batch of training images is augmented with minimal latency from data transfers.

Figure 2. Example of CT skeleton detection and segmentation using image morphology. Skeleton mask overlaid in red. Rendered by 3D Slicer [18].

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 512 × 512 slices of 12-bit data. When arranged into a 3D volume, a CT scan is hundreds of times larger than a typical low-resolution 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 out-of-bounds coordinates and interpolation between neighboring pixels. This is especially well-suited 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 first-in first-out (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 VI-B show this scheme accelerates processing by 1.69×. 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 x′ = Ax + b, where x, b ∈ R3 and A ∈ R3×3 . The matrix A is generated by composing a variety of geometric transformations drawn uniformly from user-specified ranges, including rotation, scaling, shearing, reflection and generic affine warping3 . Finally, a random displacement d ∈ R3 is drawn from a uniform distribution according to user-specified ranges. Taking c ∈ R3 to be the coordinates of the center of the volume, we then compute b according to the formula b = c + d − Ac, which guarantees Ac + b = c + d. That is, the center of the image is displaced by d units. The output image is defined by Iaffine (x) = Iin (Ax + b), where Iin : R3 → R 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δ ∈ [0, δmax ], and starting coordinate z ∈ [−δmax , zmax ], where zmax is the maximum possible zcoordinate. Then, we compute ( 0, z ≥ x3 ≥ z + δ Iocc (x) = (2) Iaffine (x), otherwise. Since we are already applying an affine transformation to the image, removing an axis-aligned prism from the output effectively removes a randomly-oriented 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 Inoise (x) = Iocc (x) + n(x), where n(x) 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 3 For reflection, the user specifies the probability that the image is reflected about each of the three axes.

6

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 -1000-1500 HU, while abdominal organs would be viewed with a narrower window of -150-230 HU. To simulate this, we randomly draw limits −∞ < a < b < ∞ according to user-specified ranges. Then we compute I (x) − a , 0 , 1 . (3) Iwindow (x) = min max noise b−a In other words, the intensity values are clamped to the range [a, b], and then affinely mapped to [0, 1]. This is straightforward to implement on GPUs. V. N EURAL

NETWORK ARCHITECTURE

This section describes the design of our predictive model; both the inputs and outputs, pre- and post-processing, the neural network itself, and the training loss function. A. Pre- and post-processing Our neural network takes as input a 120 × 120 × 160 image volume, and outputs a 120 × 120 × 160 × 6 probability map, where each voxel is assigned a class probability distribution. This becomes a 120 × 120 × 160 prediction map by taking the arg max probability for each voxel. To reduce memory requirements, we resample all image volumes to 3 mm3 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 u = (u1 , u2 , u3 ), we adjust the Gaussiansmoothing kernel according to the formula P3 2 2 g(x) ∝ exp − k=1 xk /σk where the smoothing factors are computed from the desired resolution r = 3 according to σk = 13 max(r/uk − 1, 0). This heuristic formula is based on the fact from digital signal processing that, in order to avoid aliasing, the cutoff frequency should be placed at r/uk , the ratio of sampling rates, on a [0, 1] frequency scale. After the neural network, we resample the 120 × 120 × 160 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 3 mm3 we are unlikely to fit the whole scan in our network. For training, we address this by selecting a 120 × 120 × 160 subregion from the scan uniformly at random. For inference, we cover the scan by partially-overlapping sub-regions, averaging predictions where overlap occurs. Others have addressed the

Table I L IST OF LAYERS IN OUR FULLY- CONVOLUTIONAL NEURAL NETWORK . Purpose

Decimation

Interpolation Output

Type Convolution Convolution Convolution Pooling Inception Inception Inception Inception Inception Inception Pooling Inception Inception Interpolation Softmax

Filter size 7 1 4 2 1-5 1-5 1-5 1-5 1-5 1-5 2 1-5 1-5 8 n/a

Outputs 64 64 192 192 256 480 512 512 528 832 832 832 1024 6 6

Stride 2 1 1 2 1 1 1 1 1 1 2 1 1 1/8 1

issue of limited viewing resolution by training multiple FCNs, each working at different scales [5]. However, we found that a single 3 mm3 network achieves competitive performance. B. 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 fully-connected 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 max-pooling 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 [3], [6], [8]. 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. 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 7 × 7 × 7 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 memory-efficient 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.

7

class. This approach seems simpler and more naturally motivated than the multi-class 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 LDice = 1 −

Figure 4. Example of segmentation errors with unbalanced classes. The blue color represents a prediction of class A, the pink a prediction of class B. The ground truth is the filled in pentagon, while the two circles are prediction errors. With weighted cross-entropy, the pink circle receives more weight than the blue circle, since the weight of each pixel is determined by the size of the ground truth class.

C. IOU Loss The most common loss function used for image segmentation is weighted cross-entropy. 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 p ∈ [0, 1]n denote the vector of n output probabilities from the model, where n is the number of voxels, let y ∈ {0, 1}nPdenote the binary ground truth labels, and let w = (1/n) nk=1 yk denote the weight of the class yk = 1. Let Lk (pk ) denote the cross-entropy loss function for voxel k. Then weighted cross-entropy loss is n

LCE (p) =

1X ((1 − yk )(1 − w) + yk w) Lk (pk ). n

(4)

k=1

The problem with this loss function is that it weights errors unequally: a false positive receives the weight 1 − w, whereas a false negative receives w. 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 cross-entropy, 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 0 or 1. Since {0, 1}n ⊂ [0, 1]n , we can consider y as a vector in [0, 1]n consisting only of probabilities 0 and 1. Then, the smooth IOU loss is kpyk1 (5) LIOU = 1 − kpk1 + kyk1 − kpyk1 P where kpk1 = k pk , since pk ≥ 0 for all k. To extend this to multi-class segmentation, simply average the LIOU for each

2kpyk1 . kpk1 + kyk1

(6)

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 p = (0, 1) and y = (1, 0). Then it can be verified that LDice (p, y) ≥ LDice (p, p ∪ y) + LDice (y, p ∪ y), 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 p and labels y, we can define a restriction by the taking the intersection with a binary vector s. According to the triangle inequality, L(p, g) ≤ L(p, p ∩ s) + L(g ∩ s, p ∩ s) + L(g, g ∩ s), (7) so the loss decrease is bounded by L(p, p ∩ s) + L(g, g ∩ s), which simplifies to kp ∩ sk1 /kpk1 + kg ∩ sk1 /kgk1. 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: 1) They are equal to the desired binary loss, either Dice or IOU, when p is binary. 2) They are strictly increasing in each pk when yk = 1, decreasing when yk = 0. 3) They are maximized only when p = y, minimized only when p = 1 − y. 4) They are smooth functions, if we define the loss to be 1 at p = y = 0, which is otherwise undefined. Properties 1-3 ensure that minimizing the continuous loss corresponds to maximizing the binary score of the trained model, while properties 2-4 encourage the training optimization problem to be well-behaved. 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 pk . However, this property does not extend to an average of per-image 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 p is nonbinary. In fact, these properties are not even unique among the rational functions. For example, consider the family of loss functions Pn m m k=1 pk yk P P LIOU = Pn (8) n n m m k=1 pk + k=1 yk − k=1 pk yk

8

satisfies properties 1-4. 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 lowest-order 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ǫ = kp−yk1 denote the misprediction error, and let N = kyk1 . A false negative corresponds to pk ≤ yk for all k, whence the loss is LF N = 1 −

N −ǫ (N − ǫ) + N − (N − ǫ)

(10)

N (N + ǫ) + N − N

(11)

ǫ . N Similarly, a false positive has pk ≥ yk , which yields =

LF P = 1 −

ǫ . N +ǫ So long as N ≫ ǫ, the two penalties are approximately equal. This shows that, unlike weighted cross-entropy, the IOU loss strikes a reasonable balance between each type of error. The assumption N ≫ ǫ essentially means that the model is performing well on the training data, which we should expect towards the end of training.

700

Time per volume (ms)

for any power m > 0. It is easily verified that these functions also satisfy all of the above properties, but differ from the m = 1 case for non-binary values of p. More generally, if f1 , . . . fn is a collection of smooth increasing functions on [0, 1] with fk (0) = 0 and fk (1) = 1 for all k ∈ {1, . . . , n}, then the function Pn k f k=1 fk (pk )yP P (9) LIOU = Pn n n k=1 fk (pk ) + k=1 yk − k=1 fk (pk )yk

600 500 400 300 200 100 0 0

20 CPU

40 GPU

60

80

100

120

140

Batch size

Figure 5. Comparison of our GPU data augmentation to a CPU implementation using SciPy.

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 ITK-SNAP. For the bones, we saved time by starting with the automatic labels, and then manually correcting the errors and omissions.

=

VI. E XPERIMENTAL R ESULTS 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. A. 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 full-body. 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 semi-automated 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

B. 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 3 mm3 . 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 i7-6900K 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 multi-level cache-based memory hierarchy. Execution times range from 125-74 ms per CT scan on the GPU to 323-592 ms on the CPU. Accordingly, the GPU offers a speedup of 2.6-8.1×, depending on the batch size. The largest GPU batch size is 1.69× faster per volume than the smallest, due to job pipelining. Comparing the fastest times for each processor results in a 4.4× speedup for the GPU.

9

Table II AVERAGE DICE SCORES PER CT VOLUME FOR EACH ORGAN SEGMENTER . Method Neural Nework Loss Cross entropy IOU Data augmentation No Yes No Yes Lung 87.5 86.8 90.8 91.8 Liver 88.8 85.5 90.9 92.2 Bone 73.6 65.4 78.9 79.3 Kidney 71.8 65.2 77.5 78.3 Bladder 70.4 58.4 80.1 83.7 Other 97.7 96.4 98.6 98.7 ∗ Includes all organs besides lung and bone

Morphology n/a n/a 97.8 n/a 93.2 n/a n/a 99.7∗

2.5

Loss

2 1.5

1

Figure 7. Results on a CT scan from the validation set, for each of the four training schemes. IOU loss yields tighter boundaries than cross-entropy, and IOU with augmentation is the only version yielding a reasonable bladder.

0.5 0 0

10000

20000

30000

40000

Iteration Training, cross entropy

Validation, cross enttopy

Training, IOU

Validation, IOU

Figure 6. Learning curves for cross-entropy and IOU loss. The jump in training loss corresponds to the introduction of data augmentation.

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. C. CT organ segmentation accuracy In order to compare our proposed methods to the default variant, we trained the same neural network with two different loss functions, cross-entropy 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 sub-volumes 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 cross-entropy. Surprisingly, data augmentation increases validation accuracy with the IOU loss, but decreases

the accuracy with cross entropy. This occurs even though the validation cross-entropy loss is decreasing, as shown in figure 6. This illustrates the issue described in section V-C, where weighted cross-entropy 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 cross-entropy without data augmentation, since the model can over-fit 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 post-process 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 air-filled 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, organ-specific solutions. D. 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 liver-only 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

10

Table III R ESULTS OF EACH LIVER SEGMENTER ON THE L I TS CHALLENGE . Loss Data augmentation Avg. Dice Global Dice VOE RVD ASSD MSSD RMSD

Cross entropy No Yes 0.866 0.840 0.872 0.846 0.233 0.273 -0.193 -0.248 8.506 10.407 0.767 173.722 18.299 21.908

IOU No Yes 0.896 0.905 0.904 0.911 0.185 0.168 -0.024 0.033 5.876 3.767 88.843 54.832 13.533 8.380

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 top-scoring methods, as it operates at a coarse 3mm3 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 post-processing, 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 top-performing models from incorporating our suggested improvements. VII. C ONCLUSION 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 cross-entropy 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. R EFERENCES [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,” Operations-Research-Spektrum, vol. 3, pp. 91–94, Jun 1981.

[3] F. Milletari, N. Navab, and S. Ahmadi, “V-net: 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. Ben-Cohen, 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 multi-organ segmentation on abdominal ct with dense v-networks,” 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, “Superpixel-based and boundary-sensitive 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 image-to-image network,” in Medical Image Computing and Computer-Assisted Intervention (MICCAI 2017) (M. Descoteaux, L. Maier-Hein, 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 multi-organ segmentation from ct images using conditional shape-location 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, “H-denseunet: 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: Large-scale 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 Subject-Specific 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/segment-lungs-from-3-d-chest-mri-data.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

11

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, “User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability,” Neuroimage, vol. 31, no. 3, pp. 1116–1128, 2006.