Efficient and robust cell detection: A structured regression approach.

Since microscopy images usually contain a large portion of cell clusters, it is critical to accumulate large context information for separating between cell clusters from background noise. Meanwhile, it is also important to capture fine, local information for accurately splitting touched cells. In the previous conference paper ( Xie et al., 2015c ), we use the convolutional neural network to conduct structured regression. However, we observed that it is difficult to train an ideal model that is capable of producing a large-size target proximity mask (e.g. the same size as the input image patch). Although it is reasonable to apply max-pooling operations to reduce the feature dimension and increase the receptive filed for image classification and recognition tasks, this operation actually results in massive loss of high resolution information contained in the input image, which is crucial for dense prediction problem. In this paper, we present a fully residual convolutional neural network that is capable of achieving both of those two objectives (large receptive filed and high resolution information). Instead of doing patch-wise classification, our method produces dense proximity mask that is of the same size to the input image. Our model encodes the topological structured information exhibited in the training data and explicitly forces the pixels near cell centers to get higher values than their neighbor pixels. 4.0.1.

Fully residual convolutional neural network

The detailed network architecture is illustrated in Fig. 2 . It consists of one contracting path (upper side) which encodes the input to high-level features and one expanding path (lower side) that decodes those features to the output mask. The contracting path consists of a repeated stack of 3×3 convolution followed by a residual block ( Fig. 3 ) and a 2×2 down-sampling layer. Before we down-sample the feature maps, we double the feature map channels using 1×1 convolution until it reaches 256. All the down-sampling operations used in our method are mean-pooling. The biggest difference between the expanding and contracting paths is that the down-sampling in the contracting path is replaced with up-sampling. In this paper, we use bilinear interpolation to up-sample the feature maps. To compensate the massive loss of high-resolution information during down-sampling layers and improve the localization accuracy, we concatenate the higher resolution feature maps from the contracting path to the corresponding up-sampled feature maps in the expanding path. Residual connection An identity mapping in residual blocks can be expressed as: (1) x l + 1 = x l + ℱ ( x l , 𝒲 l ) where, ℱ denotes the residual function, 𝒲 l the parameters of the l –th residual block. This type of residual identity mapping has a nice back-propagation property that the gradient does not vanish even when the weights are arbitrarily small. In He et al. (2015) , the ℱ that achieves highest classification accuracy consists of a sequence of layers: BN-ReLu-Conv-BN-ReLu-Conv, where BN, ReLu, Conv stands for Batch Normalization, rectified linear activation, and convolution respectively. During the model training, we find that adding batch normalization in the residual blocks makes the training slower and does not have positive effects on the performance, especially on breast cancer dataset. One of the possible reasons is that we already use extensively dropout in the model. However, as is observed in Shah et al. (2016) , without batch normalization, the gradient tends to explode for large network. To stabilize the training Szegedy et al. (2016) , we scale down the activation of the last convolution layer in ℱ before adding them to the input, we pick 0.3 as the scaling factor in all of our experiments. All the convolution kernels are 3×3. To prevent overfiting, we add dropout between two convolutional layers. We also replace the ReLu activation function with the more recent ELU Clevert et al. (2015) , and thus the resulting residual function can be summarized as: ELU-Conv-Dropout-ELU-Conv-Scaling . This type of residual block is illustrated in Fig 3b . The residual connection require that the input should have the identical dimension to the output. In this paper, we use 1×1 convolution to map the feature channels to the desired output’s dimension when the dimensions of input and the output mismatch with each other. Feature map concatenation During the expanding path, up-sampling layer is used to expand the feature map size, however, the up-sampling layer and down-sampling layer are not exactly ’symmetric’. For instance, the size of 25×25 feature maps after down-sampling operation becomes 12×12. However, when we up-sample 12×12 feature maps back, we can only get 24×24 feature maps. In order to preserve the size during this process, we pad the up-sampled feature maps to make match the size of the corresponding feature maps in the contracting path. Those two set of feature maps are then concatenated and passed to the following convolutional layers to inference the output. This process is illustrated using brown arrows in Fig 2 . 4.1. Structured regression for cell detection 4.1.1.

Data preprocessing Denote x ∈

R d × d × c as one local image patch extracted from image I at location ( u, v ), in which d and c represent the patch size and image channel, respectively. For simplicity, we only use square local image patches, and x can be identified by one quintuple { u, v, d, c, I }. Please note that image patches have the same image channel with the original image For each image I with human annotations, we compute the corresponding proximity map M using the following function: (2) ℳ ( u , v ) = { e α ( 1 − D ( u , v ) d ) − 1 e α − 1 if D ( u , v ) ≤ d , 0 otherwise , in which D ( u, v ) is the Euclidean distance from pixel ( u, v ) to the closest human annotation. d is a distance threshold and α is the decay ration, and both of them are used to control the shape of this exponential function. This function is a normalized version of the one used in ( Kainz et al., 2015 ; Sironi et al., 2014 ), in practice, we can choose a scaling factor to scale up this proximity value. After obtaining the proximity map M for image I , we can define the proximity patch s ∈ R d × d for x . s can be viewed as the structured label of x , and can be identified as a quintuple { u, v, d , ℳ}. This data generation process is illustrated in Fig. 3a . 4.1.2. Inference in structured regression We define { f l } l = 1 L as the transformation of each of the L layers parameterized by { θ l } l = 1 L , respectively. Our goal is to learn a mapping function ψ = f L ◦ f L −1 ◦ ⋯ f 1 , which maps the input local image patch to a proximity patch. Please note that, f i is a general notation for the transformation of i -th layer, the corresponding θ i has distinct forms for different types of f i . For example, θ i is given as [ W i , b i ] if f i denotes a conventional fully connected layer. Given one input x i , the network computes the output o i as ψ ( x i ; θ 1 , …, θ L ). In order to evaluate the model’s parameters, we formulate the structured regression as the following optimization problem: (3) arg min θ 1 , … , θ L 1 𝒩 ∑ i = 1 𝒩 ℒ ( ψ ( x i ; θ 1 , … , θ L ) , y i ) , in which ℒ is the loss function. A widely adopted loss function is mean square error, but in our case, a dominant portion of the value in the proximity patch is zeros and only a small portion of pixels has positive response. This might bias our model to produce trivial outputs containing all zeros. To solve this problem, we adopt a weighting strategy to allow the model to assign different weights to the loss coming from different regions of the proximity patches. More specifically, it is defined as: (4) ℒ ( ψ ( x i ; θ 1 , … , θ L ) , y i ) = 1 2 ∑ j = 1 p ( β y i j + λ y ¯ i ) ( y i j − o i j ) 2 , in which o i j denotes the j -th element of o i , ȳ i represents the mean value of y i , β, b and λ are predefined constants and used to tune the weights of the losses coming from different parts of the model’s output. This loss function does not use a fixed weight for every training sample; instead, it allows the model to determine based on the mean value of the training proximity patch. Denote a i as the inputs to the last layer for training sample x i . We can obtain a i = ψ( x i ; θ 1 , …, θ L −1 ). We denote the j -th element of y i , and a i as y i j and a i j , respectively. In order to back propagate the gradients from the last layer (structured regression layer) to the lower layers, we need to calculate the partial derivative of (4) with respect to the input to the last layer. More specifically, if the activation function is chosen to be sigmoid in the last layer, it can be given by (5) ∂ ℒ ( o i , y i ) ∂ a i j = ∂ ℒ ( o i , y i ) ∂ o i j ∂ o i j ∂ a i j = ( β y i j + λ y ¯ i ) ( o i j − y i j ) a i j ( 1 − a i j ) . After computing the value of (5) , we can calculate the gradients of (3) with respect to { θ l } l = 1 L recursively using the chain-rule based back-propagation algorithm. Our algorithm uses a more complicated output layer since the regression is performed on proximity patches y i ∈ 𝒴 that encode topological information. The output of the proposed model is explicitly computed as quantitative predictions related to the locations of cell centers and thus produce much more precise and robust local maxima for locating cell centers than conventional sliding-window based classification methods. After obtaining the proximity prediction map (denoted as 𝒫), a small threshold ξ ∈ [0, 1] is applied to remove the values smaller than ξ · max(𝒫). The final procedure for nucleus localization is to find all the local maximum locations in 𝒫.