Dense Transformer Networks for Brain Electron Microscopy Image Segmentation

1 Introduction

In recent years, deep convolution neural networks (CNNs) have achieved promising performance on many artificial intelligence tasks, including image recognition [LeCun et al., 1998], object detection [Sermanet et al., 2014], and segmentation [Chen et al., 2015]. Among these tasks, dense prediction tasks take images as inputs and generate output maps with similar or the same size as the inputs. For example, in image semantic segmentation, we need to predict a label for each pixel on the input images [Long et al., 2015; Noh et al., 2015]. Other examples include depth estimation [Laina et al., 2016], image super-resolution [Dong et al., 2016], and surface normal prediction [Eigen and Fergus, 2015]. These tasks can be generally considered as image-toimage translation problems in which inputs are images, and outputs are label maps [Isola et al., 2017].

Given the success of deep learning methods on imagerelated applications, numerous recent attempts have been made to solve dense prediction problems using CNNs. A central idea of these methods is to extract a square patch centered on each pixel and apply CNNs on each of them to compute the label of the center pixel. The efficiency of these approaches can be improved by using fully convolutional or encoder-decoder networks. Specifically, fully convolutional networks [Long et al., 2015] replace fully connected layers with convolutional layers, thereby allowing inputs of arbitrary size during both training and test. In contrast, deconvolution networks [Noh et al., 2015] employ an encoderdecoder architecture. The encoder path extracts high-level representations using convolutional and pooling layers. The decoder path uses deconvolutional and up-pooling layers to recovering the original spatial resolution. In order to transmit information directly from encoder to decoder, the UNet [Ronneberger et al., 2015] adds skip connections [He et al., 2016] between the corresponding encoder and decoder layers. A common property of all these methods is that the label of any pixel is determined by a regular (usually square) patch centered on that pixel. Although these methods have achieved considerable practical success, there are limitations inherent in them. For example, once the network architecture is determined, the patches used to predict the label of each pixel is completely determined, and they are commonly of the same size for all pixels. In addition, the patches are usually of a regular shape, e.g., squares.

∗Contact Author

In this work, we propose the dense transformer networks to address these limitations. Our method follows the encoderdecoder architecture in which the encoder converts input images into high-level representations, and the decoder tries to make pixel-wise predictions by recovering the original spatial resolution. Under this framework, the label of each pixel is also determined by a local patch on the input. Our method allows the size and shape of every patch to be adaptive and datadependent. In order to achieve this goal, we propose to insert a spatial transformer layer [Jaderberg et al., 2015] in the encoder part of our network. We propose to use nonlinear transformations, such as these based on thin-plate splines [Bookstein, 1989]. The nonlinear spatial transformer layer transforms the feature maps into a different space. Therefore, performing regular convolution and pooling operations in this space corresponds to performing these operations on irregular patches of different sizes in the original space. Since the nonlinear spatial transformations are learned automatically from data, this corresponds to learning the size and shape of each patch to be used as inputs for convolution and pooling operations.

There has been prior work on allowing spatial transformations or deformations in deep networks [Jaderberg et al., 2015; Dai et al., 2017], but they do not address the spatial correspondence problem, which is critical in dense prediction tasks. The difficulty in applying spatial transformations to dense prediction tasks lies in that the spatial correspondence between input images and output label maps needs to be preserved. A key innovation of this work is that we provide a new technical solution that not only allows data-dependent learning of patches but also enables the preservation of spatial correspondence. Specifically, although the patches used to predict pixel labels could be of different sizes and shapes, we expect the patches to be in the spatial vicinity of pixels whose labels are to be predicted. By applying the nonlinear spatial transformer layers in the encoder path as described above, the spatial locations of units on the intermediate feature maps after the spatial transformation layer may not be preserved. Thus a reverse transformation is required to restore the spatial correspondence.

In order to restore the spatial correspondence between inputs and outputs, we propose to add a corresponding decoder layer. A technical challenge in developing the decoder layer is that we need to map values of units arranged on input regular grid to another set of units arranged on output grid, while the nonlinear transformation could map input units to arbitrary locations on the output map. We develop a interpolation method to address this challenge. Altogether, our work results in the dense transformer networks, which allow the prediction of each pixel to adaptively choose the input patch in a data-dependent manner. The dense transformer networks can be trained end-to-end, and gradients can be back-propagated through both the encoder and decoder layers. Experimental results on biological images demonstrate the effectiveness of the proposed dense transformer networks.

2 Spatial Transformer Networks Based on Thin-Plate Spline

Spatial transformer networks [Jaderberg et al., 2015] are deep models containing spatial transformer layers. These layers explicitly compute a spatial transformation of the input feature maps. They can be inserted into convolutional neural networks to perform explicit spatial transformations. The spatial transformer layers consist of three components; namely, the localization network, grid generator and sampler.

The localization network takes a set of feature maps as input and generates parameters to control the transformation. If there are multiple feature maps, the same transformation is applied to all of them. The grid generator constructs transformation mapping between input and output grids based on parameters computed from the localization network. The sampler computes output feature maps based on input feature maps and the output of grid generator. The spatial transformer layers are generic and different types of transformations, e.g., affine transformation, projective transformation, and thin-plate spline (TPS), can be used. Our proposed work is based on the TPS transformation, and it is not described in detail in the original paper [Jaderberg et al., 2015]. Thus we provide more details below.

2.1 Localization Network

When there are multiple feature maps, the same transformation is applied to all of them. Thus, we assume there is only one input feature map below. The TPS transformation is determined by 2K fiducial points among which K points lie on the input feature map and the other K points lie on the output feature map. On the output feature map, the K fiducial points, whose coordinates are denoted as ˜ F = [ ˜ f1, ˜ f2, · · · , ˜ fK] ∈ R2×K , are evenly distributed on a fixed regular grid, where ˜ fi = [˜xi, ˜ yi]T denotes the coordinates of the ith point. The localization network is used to learn the K fiducial points F = [f1, f2, · · · , fK] ∈ R2×K on the input feature map. Specifically, the localization network, denoted as floc(·), takes the input feature maps U ∈ RH×W×C as input, where H, W and C are the height, width and number of channels of input feature maps, and generates the normalized coordinates F as the output as F = floc(U).

A cascade of convolutional, pooling and fully-connected layers is used to implement floc(·). The output of the final fully-connected layer is the coordinates F on the input feature map. Therefore, the number of output units of the localization network is 2K. In order to ensure that the outputs are normalized between -1 and 1, the activation function tanh(·) is used in the fully-connected layer. Since the localization network is differentiable, the K fiducial points can be learned from data using error back-propagation.

2.2 Grid Generator

For each unit lying on a regular grid on the output feature map, the grid generator computes the coordinate of the corresponding unit on the input feature map. This correspondence is determined by the coordinates of the fiducial points F and ˜ F. Given the evenly distributed K points ˜ F = [ ˜ f1, ˜ f2, · · · , ˜ fK] on the output feature map and the K fiducial points F = [f1, f2, · · · , fK] generated by the localization network, the transformation matrix T in TPS can be expressed as follows:

where ∆ ˜ F ∈ R(K+3)×(K+3) is a matrix determined only by ˜ F as

where R ∈ RK×K , and its elements are defined as ri,j = d2i,j lnd2i,j , and di,j denotes the Euclidean distance between ˜ fi and ˜ fj .

Through the mapping, each unit (˜xi, ˜ yi) on the output feature map corresponds to unit (xi, yi) on the input feature map. To achieve this mapping, we represent the units on the regular output grid by {˜ pi} ˜H× ˜W i=1 , where ˜ pi = [˜xi, ˜ yi]T is the (x, y)- coordinates of the ith unit on output grid, and ˜H and ˜W are the height and width of output feature maps. Note that the fiducial points { ˜ fi}Ki=1 are a subset of the points {˜ pi} ˜H× ˜W i=1 , which are the set of all points on the regular output grid.

To apply the transformation, each point ˜ pi is first extended from R2 space to RK+3 space as ˜ qi = [1, ˜xi, ˜ yi, si,1, si,2, · · · , si,K]T ∈ RK+3 , where si,j = e2i,j lne2i,j , and ei,j is the Euclidean distance between ˜ pi and ˜ fj . Then the transformation can be expressed as

where T is defined in Eq. (1). By this transformation, each coordinate (˜xi, ˜ yi) on the output feature map corresponds to a coordinate (xi, yi) on the input feature map. Note that the transformation T is defined so that the points ˜ F map to points F .

2.3 Sampler

The sampler generates output feature maps based on input feature maps and the outputs of grid generator. Each unit ˜ pi on the output feature map corresponds to a unit pi on the input feature map as computed by Eq. (3). However, the coordinates pi = (xi, yi)T computed by Eq. (3) may not lie exactly on the input regular grid. In these cases, the output values need to be interpolated from input values lying on regular grid. For example, a bilinear sampling method can be used to achieve this. Specifically, given an input feature map U ∈ RH×W , the output feature map V ∈ R ˜H× ˜W can be obtained as

for i = 1, 2, · · · , ˜H × ˜W, where Vi is the value of pixel i, Unm is the value at (n,m) on the input feature map, pi = (xi, yi)T , and pi is computed from Eq. (3). By using the transformations, the spatial transformer networks have been shown to be invariant to some transformations on the inputs. Other recent studies have also attempted to make CNNs to be invariant to various transformations [Jia et al., 2016; Henriques and Vedaldi, 2016; Cohen and Welling, 2016; Dieleman et al., 2016].

3 Dense Transformer Networks

The central idea of CNN-based method for dense prediction is to extract a regular patch centered on each pixel and apply CNNs to compute the label of that pixel. A common property of these methods is that the label of each pixel is determined by a regular (typically square) patch centered on that pixel. Although these methods have been shown to be effective on dense prediction problems, they lack the ability to learn the sizes and shapes of patches in a data-dependent manner. For a given network, the size of patches used to predict the labels of each center pixel is determined by the network architecture. Although multi-scale networks have been proposed to allow patches of different sizes to be combined [Farabet et al., 2013], the patch sizes are again determined by network architectures. In addition, the shapes of patches used in CNNs are invariably regular, such as squares. Ideally, the shapes of patches may depend on local image statistics around that pixel and thus should be learned from data. In this work, we propose the dense transformer networks to enable the learning of patch size and shape for each pixel.

As illustrated in Figure 2, features extracted by the original convolutional operation corresponds to a regular path on input feature maps. We employ a nonlinear TPS transformation in our model. Therefore, a regular patch ˜ f1 ˜ f2 ˜ f3 ˜ f4 corresponds to an area f1f2f3f4 with different shape and size. The parameters of nonlinear transformation are learned by a network based on inputs. Therefore, our model can learn an appropriate transformation to achieve better performance. To tackle the spatial correspondence problem, a reverse transformation is used to restore the spatial correspondence. The reverse TPS transformation share parameters with the TPS transformation in the encoder network.

3.1 An Encoder-Decoder Architecture

In order to address the above limitations, we propose to develop a dense transformer network model. Our model employs an encoder-decoder architecture in which the encoder path extracts high-level representations using convolutional and pooling layers and the decoder path uses deconvolution and un-pooling to recover the original spatial resolution [Noh et al., 2015; Ronneberger et al., 2015; Badrinarayanan et al., 2015]. To enable the learning of size and shape of each patch automatically from data, we propose to insert a spatial transformer module in the encoder path in our network. As has been discussed above, the spatial transformer module transforms the feature maps into a different space using nonlinear transformations. Applying convolution and pooling operations on regular patches in the transformed space is equivalent to operating on irregular patches of different sizes in the original space. Since the spatial transformer module is differentiable, its parameters can be learned with error backpropagation algorithms. This is equivalent to learning the size and shape of each patch from data.

Although the patches used to predict pixel labels could be of different sizes and shapes, we expect the patches to include the pixel in question at least. That is, the patches should be in the spatial vicinity of pixels whose labels are to be predicted. By using the nonlinear spatial transformer layer in encoder path, the spatial locations of units on the intermediate feature maps could have been changed. That is, due to this nonlinear spatial transformation, the spatial correspondence between input images and output label maps is not retained in the feature maps after the spatial transformer layer. In order to restore this spatial correspondence, we propose to add a corresponding decoder layer, known as the dense transformer decoder layer. This decoder layer transforms the intermediate feature maps back to the original input space, thereby reestablishing the input-output spatial correspondence.

The spatial transformer module can be inserted after any layer in the encoder path while the dense transform decoder module should be inserted into the corresponding location in decoder path. In our framework, the spatial transformer module is required to not only output the transformed feature maps, but also the transformation itself that captures the spatial correspondence between input and output feature maps. This information will be used to restore the spatial correspondence in the decoder module. Note that in the spatial transformer encoder module, the transformation is computed in the backward direction, i.e., from output to input feature maps (Figure 1). In contrast, the dense transformer decoder module uses a forward direction instead; that is, a mapping from input to output feature maps. This encoder-decoder pair can be implemented efficiently by sharing the transformation parameters in these two modules.

A technical challenge in developing the dense transformer decoder layer is that we need to map values of units arranged on input regular grid to another set of units arranged on regular output grid, while the decoder could map to units at arbitrary locations on the output map. That is, while we need to compute the values of units lying on regular output grid from values of units lying on regular input grid, the mapping itself could map an input unit to an arbitrary location on the output feature map, i.e., not necessarily to a unit lying exactly on the output grid. To address this challenge, we develop a sampler method for performing interpolation. We show that the proposed samplers are differentiable, thus gradients can be propagated through these modules. This makes the entire dense transformer networks fully trainable. Formally, assume that the encoder and decoder layers are inserted after the i-th and j-th layers, respectively, then we have the following relationships:

where U i is the feature map of the i-th layer, p is the coordinate of a point, T is the transformation defined in Eq. (1), which maps from the coordinates of the (i+1)-th layer to the i-th layer, Sampling(·) denotes the sampler function.

From a geometric perspective, a value associated with an estimated point in bilinear interpolation in Eq. (4) can be interpreted as a linear combination of values at four neighboring grid points. The weights for linear combination are areas of rectangles determined by the estimated points and four neighboring grid points. For example, in Figure 1, when a point is mapped to P on input grid, the contributions of points A, B, C, and D to the estimated point P is determined by the areas of the rectangles S1, S2, S3, and S4. However, the interpolation problem needs to be solved in the dense transformer decoder layer is different with the one in the spatial transformer encoder layer, as illustrated in Figure 1. Specifically, in the encoder layer, the points A, B, C, and D are associated with values computed from the previous layer, and the interpolation problem needs to compute a value for P to be propagated to the next layer. In contrast, in the decoder layer, the point P is associated with a value computed from the previous layer, and the interpolation problem needs to compute values for A, B, C, and D. Due to the different natures of the interpolation problems need to be solved in the encoder and decoder modules, we propose a new sampler that can efficiently interpolate over decimal points in the following section.

3.2 Decoder Sampler

In the decoder sampler, we need to estimate values of regular grid points based on those from arbitrary decimal points, i.e., those that do not lie on the regular grid. For example, in Figure 1, the value at point P is given from the previous layer. After the TPS transformation in Eq. (3), it may be mapped to an arbitrary point. Therefore, the values of grid points A, B, C, and D need to be computed based on values from a set of arbitrary points. If we compute the values from surrounding points as in the encoder layer, we might have to deal with a complex interpolation problem over irregular quadrilaterals. Those complex interpolation methods may yield more accurate results, but we prefer a simpler and more efficient method in this work. Specifically, we propose a new sampling method, which distributes the value of P to the points A, B, C, and D in an intuitive manner. Geometrically, the weights associated with points A, B, C, and D are the area of the rectangles S1, S2, S3, and S4, respectively (Figure 1). In particular, given an input feature map V ∈ R ˜H× ˜W , the output feature map U ∈ RH×W can be obtained as

where Vi is the value of pixel i, pi = (xi, yi)T is transformed by the shared transformation T in Eq. (1), Unm is the value at the (n,m)-th location on the output feature map, Snm is a normalization term that is used to eliminate the effect that different grid points may receive values from different numbers of arbitrary points, and n = 1, 2, · · · ,N, m = 1, 2, · · · ,M. More details are given in Figure 2.

In order to allow the backpropagation of errors, we define the gradient with respect to Unm as dUnm. Then the gradient with respect to Vnm and xi can be derived as follows:

A similar gradient can be derived for dyi. This provides us with a differentiable sampling mechanism, which enables the gradients flow back to both the input feature map and the sampling layers.

4 Experimental Evaluation

We evaluate the proposed methods on two image segmentation tasks. The U-Net [Ronneberger et al., 2015] is adopted as our base model in both tasks, as it has achieved state-ofthe-art performance on image segmentation tasks. Specifically, U-Net adds residual connections between the encoder path and decoder path to incorporate both low-level and highlevel features. Other methods like SegNet [Badrinarayanan et al., 2015], deconvolutional networks [Zeiler et al., 2010] and FCN [Long et al., 2015] mainly differ from U-Net in the up-sampling method and do not use residual connections. Experiments in prior work [Ronneberger et al., 2015; Zeiler et al., 2010; Long et al., 2015] show that residual connections are important while different up-sampling methods lead to similar results. The network consists of 5 layers in the encoder path and another corresponding 5 layers in the decoder path. A stack of two 3×3 convolutional layers have the same receptive field as a 5×5 convolutional layer, but with less parameters [Simonyan and Zisserman, 2015]. Therefore, we use 3×3 kernels and one pixel padding to retain the size of feature maps at each level.

In order to efficiently implement the transformations, we insert the spatial encoder layer and dense transformer decoder layer into corresponding positions at the same level. Using spatial transformation layers at a deep position in the encoderdecoder network can increase the receptive field of the spatial transformation operation. Therefore, the layers are applied to the 4th layer, and their performance is compared to the basic U-Net model without spatial transformations. As for the transformation layers, we use 16 fiducial points that are evenly distributed on the output feature maps. In the dense transformer decoder layer, if there are pixels that are not selected on the output feature map, we apply an interpolation strategy over its neighboring pixels on previous feature maps to produce smooth results.

4.1 Brain ElectronMicroscopy Image Segmentation

We evaluate the proposed methods on brain electron microscopy (EM) image segmentation task [Lee et al., 2015; Ciresan et al., 2012], in which the ultimate goal is to reconstruct neurons at the micro-scale level. A critical step in neuron reconstruction is to segment the EM images. We use data set from the 3D Segmentation of Neurites in EM Images (SNEMI3D, http://brainiac2.mit.edu/SNEMI3D/). The SNEMI3D data set consists of 100 1024×1024 EM image slices. Since we perform 2D transformations in this work, each image slice is segmented separately in our experiments. The task is to predict each pixel as either a boundary (denoted as 1) or a non-boundary pixel (denoted as 0).

Our model can process images of arbitrary size. However, training on whole images may incur excessive memory requirement. In order to accelerate training, we randomly pick 224×224 patches from the original images and use it to train the networks. The experimental results in terms of ROC curves are provided in Figure 4. We can observe that the proposed DTN model achieves higher performance than the baseline U-Net model, improving AUC from 0.8676 to 0.8953. These results demonstrate that the proposed DTN model improves upon the baseline U-Net model, and the use of the dense transformer encoder and decoder modules in the U-Net architecture results in improved performance. Some example results along with the raw images and ground truth label maps are given in Figure 3.

4.2 Timing Comparison

Table 1 shows the comparison of training and prediction time between the U-Net model and the proposed DTN model on the two data sets. We can see that adding DTN layers leads to only slight increase in training and prediction time.

5 Conclusion

In this work, we propose the dense transformer networks to enable the automatic learning of patch sizes and shapes in dense prediction tasks. This is achieved by transforming the intermediate feature maps to a different space using nonlinear transformations. A unique challenge in dense prediction tasks is that, the spatial correspondence between inputs and outputs should be preserved in order to make pixel-wise predictions. To this end, we develop the dense transformer decoder layer to restore the spatial correspondence. The proposed dense transformer modules are differentiable. Thus the entire network can be trained from end to end. Experimental results show that adding the spatial transformer and decoder layers to existing models leads to improved performance. To the best of our knowledge, our work represents the first attempt to enable the learning of patch size and shape in dense prediction. The current study only adds one encoder layer and one decoder layer in the baseline models. We will explore the possibility of adding multiple encoder and decoder layers at different locations of the baseline model. In this work, we develop a simple and efficient decoder sampler for interpolation. A more complex method based on irregular quadrilaterals might be more accurate and will be explored in the future.

Acknowledgments

This work was supported in part by National Science Foundation grants IIS-1615035 and DBI1641223.