Deep Learning–Based Detection of Endothelial Tip Cells in the Oxygen-Induced Retinopathy Model

Annotated Data Set of Retina Flat Mount Images

Data sets

We prepared a data set of 220 retina flat mount images with labeled tip cells for training of the tip cell detector. All images were from either different mice or different eyes of the same mouse (140 mice with flat mount images from only one eye and 40 mice with flat mount images from both eyes). A skilled person (N.Z.) manually labeled 17,994 tip cells in these images. N.Z. is a biomedical scientist with years-long experience in the field of angiogenesis. Additionally, 2 smaller validation and test data sets were built, consisting of 16 and 25 annotated flat mount images, correspondingly. The validation data set was used to tune parameters related to the training process and the system architecture. The separated test data set was used for the estimation of the final system performance. All the images are gray-scale images of 5120 × 5120 pixel size (see an example in Figure 1). There was no overlap between mice used for test and train (or validation) data sets. Since our detector is designed to work with images of 512 × 512 size, the images, before being processed, are tiled to patches of this size. Correspondingly, the coordinates of labeled tip cells were updated to match the created patches. For training, we generated patches with a small overlap of 12.5% to avoid loss of labeled tip cells located at patch borders.

The labeling work was done with the use of the Zeiss Efficient Navigation (ZEN) software (Carl Zeiss). Zeiss Efficient Navigation is an image acquisition, processing, and analysis software for digital microscopy developed by Zeiss. The software allows image manipulation and annotation of image regions with various graphical tools. In our case, tip cell locations (tip cell ends) were labeled with point markers (see Figure 2, top-left). We designed a software that collects the labels saved by ZEN in xml files and feeds them to our training pipeline. In most cases, tip cells are directed into the direction of avascular area, have pointed ends due to filopodial extensions, and are located at tips of vessel sprouts. The typical range of the tip cell sizes is between 25 and 60 µm, though usually we do not see the whole 3-dimensional cell in the 2D images.

OPEN IN VIEWER

Figure 2.

An example of a retina flat mount image with tip cells labeled by the reference expert (A, B) and automatically detected by the system (C, D). Magnified regions of ground truth labels and detections are shown on the right. Scale bar: 500 µm.

Flat mount masks

We use binary masks of flat mount areas in order to exclude background regions (black in Figure 3) from being used in a few stages of the detection pipeline (sections Redesigned SSD: DNN Model, Hard Negative Mining, Post-processing). To generate flat mount masks, images were first smoothed by an edge-preserving smoothing algorithm ²⁴ and scaled down to 25% of the original size. The region growing algorithm was applied after image gray values were converted to the [0-255] range. The region growing started from 8 points taken from the vicinity of the image corners. At any point in the process of the region growing, the mean gray value of the current region is calculated. New pixels were added to the region if they differ from the current region mean by less than 3 gray levels. The complement of the generated image was further smoothed with morphological closing, which also eliminated small holes. The resulting image mask was scaled back to the original image size.

OPEN IN VIEWER

Figure 3.

Retina flat mount mask of the image in Figure 1.

Redesigned Single-Shot-Detector (SSD)

Application of DNNs to the object detection task has shown remarkable advances in the last years. ^25,26 The first DNN detectors were based on patch image classifiers. Later, they evolved to more computationally effective detectors with dedicated architectures. Such DL-based detectors outperform all the other types of detectors with a large margin when evaluated on large data sets of natural images.

There are 2 types, one and two stage, DNN object detectors. The two-stage object detectors include a fast object proposal unit as a first stage, where an object of any type is first detected. It is then categorized at the second, more computationally demanding object classification stage. Since in our case we detect only a single type of objects, tip cells versus a vascular background, we use an one-stage object detector, where the object of interest is directly detected.

Redesigned SSD: DNN Model

We have chosen one of the successful (state-of-the-art 2016) one-stage DNN object detectors, the SSD. ²⁷ The SSD has a relatively simple architecture that generalizes fairly well on different types of images, which over a few recent years allowed its extensive use in a broad range of applications. It is one of the fast detectors and has a few open-source high-quality implementations available online, upon one of which we built our system. ²⁸ Based on the principles of transfer learning, the SSD uses a truncated VGG-16 network ²⁹ pretrained on the ImageNet data set. ³⁰ In the SSD, the fully connected VGG layers were converted to convolutional layers and extra convolutional layers were added to generate several multiscale feature maps.

Our task is more restricted than generic object detection. We can therefore refine the original SSD architecture, thereby reducing the number of parameters to be learned. The refined SSD-512 architecture is shown in Figure 4. We experimented with the number of multiscale feature maps and concluded that for our task layers on the top of the truncated VGG base network do not improve the performance of our system. This is explained by the fact that the structures of interest are of approximately similar size, which is relatively small. Detection of objects of such size requires the use of the feature map at highest resolution. We therefore use only a single-scale feature map of 64 × 64 spatial size and exclude both SSDs extra convolutional layers and VGGs fully connected layers, which were converted to convolutional layers in the original SSD architecture. We also use only a single square default bounding box (64 pixels size), which matches the size of bounding boxes around labeled tip cells. Thus, for each location in a 64 × 64 feature map, the detector computes a pair of confidence score values c^obj, c^back and a pair of positional offset l^cx, l^cy values (see Equation 3). After a softmax operation, the resulting object’s confidence value c ^ o b j is thresholded at a suitable level to determine whether an object is present at the corresponding offset location in the input image.

OPEN IN VIEWER

Figure 4.

Refined Single Shot Detector (SSD) architecture consists of a modified VGG-16 backbone that generates a feature map at a single resolution and a predictor of confidence and position offset for each of 64 × 64 spatial locations. VGGs first layer filters are shrunk to 1 × 3 × 3 size to allow processing gray-tone images. VGGs fully connected layers are cut out as well as the last 3 convolutional layers.

We use a SSD architecture designed for images with a spatial resolution of 512 × 512 size (SSD512). Since retina flat mount images were available as one channel images, we converted the filters in the first convolutional layer of the pretrained VGG network from 3 × 3 × 3 size to 1 × 3 × 3 size, where the first dimension of the filters corresponds to color channels of an RGB image. Such conversion can be carried out either by summing up filter coefficients along the first dimension or by using only the coefficients from a single color dimension, for example, dimension that works on a green channel of an RGB image. We used the former approach, though both generated comparable results in our experiments.

Every input gray-scale image is preprocessed by means of subtraction of a mean value computed over the whole training data set. Values outside of the flat mount mask were not taken into account.

Training Procedure

The original SSD training procedure was modified to be more suitable to the type of images we deal with. For optimization with stochastic gradient descent, we set the learning rate to 5 × 10⁻⁴, which was reduced after 65,000 iterations to 5 × 10⁻⁵. We stopped training after 85,000 iterations, when performance improvement was saturating. However, the system has reached 99% of the final performance already after 30,000 iterations. As in the original SSD, we used momentum and weight decay equal 0.9 and 5 × 10⁻⁴, respectively. We used a warmup strategy ³¹ that sets a smaller learning rate in the beginning of the training procedure, which makes learning more stable. In our case, we started with learning rate 2.5 × 10⁻⁴ for the first 300 iterations and then switched to 5 × 10⁻⁴ as mentioned above. The batch size was set to 16 (due to GPU memory considerations only). Aforementioned parameter adjustments were based on performance evaluations on the validation data set (section Data sets).

Training Objective

We used a slightly modified objective function as in the original SSD consisting of the confidence and the localization losses

L ( c , l ) = 1 N ( L conf ( c ) + L loc ( l ) ) ,

1

for N predictions of confidence score c and spatial offset l made by the detector. Since objects of our interest are of approximately the same size, we do not regress the height and width of the bounding boxes, which simplifies the localization loss compared to the original SSD as follows

L loc ( l ) = ∑ i , j ∑ m ϵ { c x , c y } x i , j smooth L 1 ( l i m − g ^ j m ) ,

2

g ^ j c x = g j c x − d i c x s , g ^ j c y = g j c y − d i c y s ,

where i goes over all indexes of default bounding boxes (detector’s default locations) and j goes over all ground truth bounding boxes (object locations), x_i,j is an indicator function matching default bounding box i to the ground truth bounding box j. smooth L 1 is a robust L₁ loss used in fast R-CNN detector, ³² g and d are the coordinates of the ground truth locations and default (bounding box) locations, s denotes the predefined size of all bounding boxes (64 pixels). Eliminating regression of bounding box sizes speeded up convergence of the training procedure.

The confidence loss is the cross entropy with the softmax function computed for 2 categories only, the background and the object (tip cell):

L conf ( c ) = − ∑ i ∈ Pos log ( c ^ i obj ) − ∑ i ∈ Neg log ( c ^ i obj ) ,

3

c ^ i obj = exp ( c i obj ) exp ( c i obj ) + exp ( c i back ) , c ^ i back = exp ( c i back ) exp ( c i obj ) + exp ( c i back ) ,

where Pos and Neg are the sets of indexes of default bounding boxes that matched ground truth object or remained unmatched, respectively. Only a part of unmatched default bounding boxes (negatives) are used in objective function as described in section Hard Negative Mining. c ^ i obj and c ^ i back are the confidences of an object (tip cell) and background, respectively, computed by the detector for default bounding box i.

Data Augmentation

Since our annotated data set is relatively small (see section Annotated Dataset of Retina Flat-Mount Images for details), it was crucial to perform extensive augmentation of the data. In addition to a set of photometric and spatial augmentations carried out by the authors of the SSD, we augment the data with rotation at a random angle in the range [0-360). Since the size of tip cells does not vary significantly, we reduced the range of the random crop sizes suggested for the SSD to [0.66-1] and used zoom-out strategy that reduces objects sizes by the factor in the [1-1.3] range.

Hard Negative Mining

In the SSD approach, negative examples are default bounding boxes that do not match object labels. Since usually there are many more negative examples than labeled positives, the authors of the original SSD used a procedure called “hard negative mining” that reduces such an imbalance. According to this approach, only a part of the negatives is selected, namely the most informative cases for the training (ie, the “hard” negative cases that may resemble a tip cell). The selection of hard cases is based on the loss function. The details on the exact procedure of the hard negative mining can be found in the original SSD paper. ²⁷ We adopt the hard negatives mining procedure using the same ratio of 3:1 between negative and positive examples. We, however, also exclude large black parts outside of the cloverleaf-like flat mount area (see Figure 3) from the pool of potential locations for negative examples. This way, we put an emphasis on differentiation between various vascular structures in the retina and tip cells, rather than between the black background and tip cells.

For this purpose, we used the binary flat mount masks generated as described in section Annotated Dataset of Retina Flat-Mount Images. Additionally, border regions of retina flat mount areas were excluded from the pool of potential locations for negatives. Due to the preparation of retina flat mounts, these regions contain cut vessels that can resemble tip cells. We therefore applied morphological erosion to flat mount masks with a square structuring element of 129 pixels per side.

Postprocessing

The trained tip cell detector may also identify tip cells at the ends of leaves of the cloverleaf-like flat mount region. These tip cells are of no interest for us, because they reflect developmental retinal angiogenesis instead of revascularizarion of the damaged retina. We, therefore, add a postprocessing filter that suppresses detections that are too far away from the retina center.

For this purpose, we calculate centroid and radius of the circle that encloses retina in the flat mount masks. To speed up computations on the large masks of 5120 × 5120 pixel size, we make calculations from its edges, which are obtained by means of internal morphological gradient. ³³ The radius is estimated from the histogram of the edge distances from the centroid. Namely, we take a radius that equals the distance giving us a maximum value of the histogram. Finally, we suppress all detected tip cells that are at the larger distances than the estimated radius reduced by a constant, which was set to 300 pixels (375 μm).

System Performance

The system performance was evaluated using a separate test set consisting of 25 flat mount images of 5120 × 5120 pixel size and reference labels, which were generated by a human expert (N.Z). N.Z. is the same expert who also provided labels for the training data set (see section Data Sets for details). We will refer to this expert as the reference expert.

Figure 2 illustrates how our system performs qualitatively in comparison to the reference expert. Figure 2A shows an example of retina flat mount image with tip cells labeled by the reference expert. Figure 2C shows the same image with automatically detected tip cells. It is evident that the system identifies tip cell within the same regions as the reference expert does. The zoomed-in regions of the retina are shown on the right and visualize details of reference expert labels (Figure 2B) and automatic detections (Figure 2D). In the majority of cases, detections match the reference labels. Additionally, the examples show a few cases of tip cells that were not identified by the system, and conversely, a few structures that were not labeled by the reference expert but were detected by the system. An analysis conveyed by the reference expert has pointed out that a majority of the latter detections (not being a part of labeled tip cells) are in a gray-zone, that is, cannot not be decisively assessed. In addition, some of the novel detections were actually overlooked during the labeling process.

The detector outputs locations of tip cells and corresponding detection confidences in the [0-1] range. We quantitatively evaluated its performance using a precision-recall (PR) curve, which is commonly used for performance estimation of such detectors. The curve shows the trade-off between precision and recall values with varying detection threshold. Precision indicates the number of true detections relative to the number of all detections, while recall (or sensitivity) indicates the number of correctly detected objects relative to the overall number of objects of interest.

Figure 5 shows the PR curve for our case of the tip cell detector, when we defined a localization tolerance, equal to 64 pixels (80 µm), that sets an allowed maximal distance between predicted and ground-truth tip cell locations. A detector can operate at any point on the PR curve, which is defined by a chosen threshold on detector’s confidence. We empirically set the confidence threshold to 0.35, which corresponds to the operating point shown in red in Figure 5. An average precision is frequently used to summarize the PR curve with a single value in order to compare systems in a simplified way or to track improvements in performance at every step of a training process. Average precision equals to the mean of precision values over the range of recall values, which in our case reaches 0.90.

OPEN IN VIEWER

Figure 5.

Precision-recall curve of the deep neural network (DNN) detector, chosen operating point of the detector, and comparative performance of human experts.

To the best of our knowledge, our system is the first tip cell detection system, and we therefore compare the performance of our tip cell detector to performance of the current state of the art: human experts. For this purpose, the reference expert (N.Z.), who has an extensive experience in the field, trained 2 biomedical experts (T.S. and S.T.) to manually label tip cells in flat mount images according to the same guidelines.

T.S. is a veterinarian with a proven track record in histological evaluation of organs including eye/retina. S.T. is a trained technician with experience in sprouting angiogenesis. The performance of biomedical experts in predicting the tip cell locations in the test set was measured in terms of precision and recall pair values, which are depicted as points in Figure 5 together with the PR curve for the trained system. The figure shows that our system is superior over the human experts in detection and localization of the tip cells.

Comparison of Tip Cell Counts

Our system was built to automate tip cell counting in retina flat mounts. We therefore analyzed the correlation between the number of automatically detected tip cells (orange markers in Figure 6) and counts calculated from reference expert labels (n_ref). Figure 6 shows that the counts are indeed highly correlated, which is quantitatively expressed with a Pearson correlation coefficient r = 0.94 and the average relative difference | n ref − n det | ( n ref + n det ) / 2 of only 11%. As we did for the performance estimation in the System Performance section, we also compared the correlation of our system with the reference expert labels versus the correlation of 2 further human experts (blue and green markers) with the reference expert labels. Again, the counts based on the automatic detections have a considerably higher correlation with n_ref then the human expert counts. The figure also clearly shows a high annotation variability between experts (interobserver variability). The counts of biomedical expert 1 and expert 2 are correlated with a Pearson coefficient r = 0.56, that is, 2 human experts show considerably higher variability than the DNN detector compared to the reference expert (r = 0.94). The average relative difference between 2 biomedical experts | n exp 1 − n exp 2 | ( n exp 1 + n exp 2 ) / 2 equals 33%, again much larger than between the DNN detector and the reference expert (11%).

OPEN IN VIEWER

Figure 6.

Correlation between the number of tip cells counted by a reference expert n_ref and detected automatically with deep neural network (DNN) detector (orange), as well as a few test experts (blue and green). Gray data points were obtained from the tip cells labeled by the reference expert in the repeated attempt (one month after the first attempt). The detector was trained on a separate (training) data set of retina flat mount images annotated by the reference expert. r denotes Pearson correlation calculated from 25 data points (the number of flat mount retina images in the test data set), solid lines are obtained with linear regression from the corresponding data points.

To assess the quality of the counts of the reference expert and the consistency of human labeling (intraobserver variability), we asked the reference expert to label the same test data set again, one month after the first attempt. Gray markers in Figure 6 show the resulting counts, which turned out to be highly correlated (r = 0.99) with the counts of the first labeling session. This indicates a high quality of the labels provided by the experienced reference expert. It should be noted that with the automatic system we have a counting consistency that would be quantified with the correlation r = 1.0, since it will score always identically (as long as the DNN model is kept constant). From Figure 6, we can observe that the regression lines for the counts from the repeated labeling effort of the reference expert and from the DNN detector largely coincide, which suggests that our tip cell detection system can substitute for human annotation.

We, additionally, did an experiment comparing detected tip cell counts in 32 flat mounts captured with a Zeiss LSM and an Opera Phenix HCS system. Flat mounts were prepared at different years in comparison to the training and test data sets. The flat mount preparation was performed according to the same protocol as previously. In contrast to the images acquired with the LSM, image acquisition with HCS system was performed after one year storage of the flat mounts, which might have reduced their fluorescence intensity. To cope with very different illumination conditions and exposure times, brightness of all the images was linearly adjusted to match average brightness of the images in the training data set. An analysis of the tip cell counts detected in images acquired by the 2 systems resulted in a Pearson correlation coefficient of r = 0.86 (n = 32), see Figure 7. This shows the robustness of our system in detecting tip cells using images acquired with different devices, or under different acquisition conditions. However, as shown in Figure 7, the system detects slightly more tip cells in the images from the HCS system compared to the LSM, presumably due to slight differences in focal depth or optical resolution.

OPEN IN VIEWER

Figure 7.

Correlation between the number of tip cells detected in the images acquired with the Zeiss laser scanning microscope (LSM) and Opera high-content screening (HCS) systems. The Pearson coefficient calculated from 32 data points was r = 0.86.

Therefore, in order to study changes of tip cell counts, for example, in an efficacy study, it is advisable to compare detections only on images of the same type (device, acquisition settings, and similar flat mount preparation method). Ideally, the detection system should be retrained on examples of labeled tip cells taken from the same type of device.