Star-convex Polyhedra for 3D Object Detection and Segmentation in Microscopy

TL;DR

This paper introduces StarDist-3D, a novel method using star-convex polyhedra for accurate 3D cell nuclei detection and segmentation in microscopy images, overcoming challenges of shape complexity and anisotropic voxel sizes.

Contribution

It extends 2D star-convex shape prediction to 3D, addressing shape representation, anisotropic data adaptation, and efficient intersection computation for improved nuclei segmentation.

Findings

Outperforms classical and deep learning methods on challenging datasets

Accurately models complex nuclei shapes with star-convex polyhedra

Handles anisotropic voxel sizes effectively

Abstract

Accurate detection and segmentation of cell nuclei in volumetric (3D) fluorescence microscopy datasets is an important step in many biomedical research projects. Although many automated methods for these tasks exist, they often struggle for images with low signal-to-noise ratios and/or dense packing of nuclei. It was recently shown for 2D microscopy images that these issues can be alleviated by training a neural network to directly predict a suitable shape representation (star-convex polygon) for cell nuclei. In this paper, we adopt and extend this approach to 3D volumes by using star-convex polyhedra to represent cell nuclei and similar shapes. To that end, we overcome the challenges of 1) finding parameter-efficient star-convex polyhedra representations that can faithfully describe cell nuclei shapes, 2) adapting to anisotropic voxel sizes often found in fluorescence microscopy datasets, and 3) efficiently computing intersections between pairs of star-convex polyhedra (required for non-maximum suppression). Although our approach is quite general, since star-convex polyhedra include common shapes like bounding boxes and spheres as special cases, our focus is on accurate detection and segmentation of cell nuclei. Finally, we demonstrate on two challenging datasets that our approach (StarDist-3D) leads to superior results when compared to classical and deep learning based methods.

Tables

Table 1. Table 1 : Accuracy (average over 5 5 5 independent trials for trained models) for several IoU thresholds τ 𝜏 \tau for datasets Worm and Parhyale .

Threshold $\tau$	0.10	0.20	0.30	0.40	0.50	0.60	0.70	0.80	0.90
Worm
IFT-Watershed	$0.794$	$0.771$	$0.708$	$0.601$	$0.472$	$0.364$	$0.222$	$0.074$	$0.005$
U-Net	$0.873$	$0.861$	$0.814$	$0.706$	$0.570$	$0.418$	$0.255$	$0.116$	$0.027$
U-Net+	$0.920$	$0.905$	$0.872$	$0.807$	$0.700$	$0.593$	$0.406$	$0.144$	$0.005$
StarDist-3D	$0.936$	$0.926$	$0.905$	$0.855$	$0.765$	$0.647$	$0.460$	$0.154$	$0.004$
Parhyale
IFT-Watershed	$0.467$	$0.426$	$0.356$	$0.245$	$0.161$	$0.096$	$0.036$	$0.000$	$0.000$
U-Net	$0.547$	$0.501$	$0.423$	$0.330$	$0.247$	$0.171$	$0.091$	$0.021$	$0.000$
U-Net+	$0.592$	$0.552$	$0.481$	$0.372$	$0.280$	$0.198$	$0.097$	$0.010$	$0.000$
StarDist-3D	$0.766$	$0.757$	$0.741$	$0.698$	$0.593$	$0.443$	$0.224$	$0.038$	$0.000$

Equations

$$z_k$$

$$y_k$$

$$x_k$$

$$L(p,\hat{p},d_k,{\hat{d}}_k)=L_{obj}(p,\hat{p})+{\lambda }_d{L}_{dist}(p,\hat{p},d_k,{\hat{d}}_k).$$

$$L_{obj}(p,\hat{p})=-p\log \hat{p}-(1-p)\log (1-\hat{p}).$$

$$\begin{gathered} L_{dist}(p,\hat{p},d_k,{\hat{d}}_k)=p\cdot \text{mathbbm}{1}_{p>0}\cdot \frac{1}{n}{\sum }_k|d_k-{\hat{d}}_k|+ \\ {\lambda }_{reg}\cdot \text{mathbbm}{1}_{p=0}\cdot \frac{1}{n}{\sum }_k|{\hat{d}}_k|. \end{gathered}$$

Full text

11footnotetext: Equal contribution

1 Introduction

Detection and segmentation of cell nuclei in volumetric (3D) fluorescence microscopy images is a ubiquitous problem in developmental biology and often constitutes the first step when studying cellular expression patterns, or when tracing cell lineages in developing organisms [17, 22]. The task of nuclei detection is to roughly locate all individual nuclei inside a 3D volume, e.g. by enumerating their center points or bounding boxes. On the other hand, semantic segmentation aims to label each pixel with a semantic class (e.g., nucleus or background), but is not concerned with discerning individual nuclei. Finally, instance segmentation is more ambitious since it combines these tasks by seeking a separate label mask for each individual nucleus. As modern microscopes produce increasingly large 3D datasets, many automated instance segmentation methods have been proposed over the years [17]. These include classical thresholding approaches with pixel-grouping via connected component, morphological methods based on the watershed transform [4, 16, 7], and optimization via graph-cuts [5]. More recently, methods based on deep learning have been shown to vastly improve results for natural and biomedical images alike [11, 24, 25].

In general, deep learning based instance segmentation can be roughly categorized into 1) methods that first perform semantic segmentation followed by grouping of pixels into distinct objects (e.g. U-Net [9, 6]), and 2) methods that first predict axis-aligned bounding boxes of individual objects with a subsequent semantic segmentation step for each found object (e.g. [11, 27, 26]). Despite the advances made by these methods, they often still underperform due to the low signal-to-noise ratios and dense packing of nuclei in typical fluorescence microscopy datasets. In particular, methods of category 1) are prone to erroneously fuse touching nuclei, and those of category 2) may fail to discern objects that are poorly approximated with bounding boxes.

These problems have recently been highlighted by Schmidt et al. [21] for the case of 2D fluorescence microscopy images. To alleviate these issues, [21] proposed a method called StarDist, which uses a neural network that directly predicts an appropriate shape representation (star-convex polygons) for cell nuclei and demonstrated improved results. Concretely, for every pixel inside an object (nucleus) they predict the distance to the object boundary along several radial directions, thus defining a star-convex polygon. Furthermore, they also predict an object probability to determine which pixels are part of cell nuclei and thus are allowed to vote for an object shape. Since every pixel is predicting a polygon to represent the entire shape of the object it belongs to, they perform non-maximum suppression to prune redundant shapes that likely represent the same object. Note that [21] sits somewhere in between object detection and instance segmentation because the predicted shapes are of relatively high fidelity, but are not pixel-accurate.

In this paper, we adopt and extend the approach of [21] to the case of 3D volumes and use star-convex polyhedra as shape representations for cell nuclei and similar shapes. We directly predict the polyhedra parameters densely for each pixel and then use non-maximum suppression (NMS) to prune the highly redundant set of obtained polyhedron shapes to ideally retain only one predicted shape for each true object in the image. Please see Fig. 1 for an overview of our approach. Note that we keep the benefits of [21], first and foremost to accurately disambiguate densely packed objects in images with low signal-to-noise ratios. Furthermore, star-convex polygons/polyhedra are a superset of convex shapes in 2D/3D and thus include common shapes like bounding boxes and circles/spheres as special cases.

Contributions

The extension of [21] from 2D to 3D is challenging and our main contribution in this paper. First, computing the intersection of two star-convex polyhedra (as required for NMS) efficiently is non-trivial (see Section 2.3 and Fig. 1c), but highly necessary to make this approach practical for large 3D volumes. Second, while [21] used $32$ radial directions to represent 2D nuclei shapes, a naive extension to 3D would require $32^{2}=1024$ directions which is not feasible due to the excessive amount of computation and memory required for large 3D volumes. We show that a more judicious selection of radial directions (Section 2 and Fig. 1a) enables faithful shape representations with as little as $64$ values. Third, microscopy volumes are commonly acquired with anisotropic voxel sizes that result in squeezed nuclei shapes along the axial (Z) direction. We find that it is critical to adapt the star-convex representation to account for this anisotropy of the data to achieve good results (Sections 2 and 3). Finally, we demonstrate on two challenging datasets that our proposed method (StarDist-3D) leads to superior results when compared to a classical watershed method and U-Net baselines.

2 Method

2.1 Star-convex polyhedra

We describe the 3D shape of a single object (cell nucleus) with a star-convex polyhedron. Concretely, for each pixel inside an object we compute the distances $d_{k}$ to the object boundary along a fixed set of $n$ unit rays $\vec{r}_{k}$. To obtain a faithful shape model, we use rays that are approximately evenly distributed on an ellipsoid representative of the objects in a dataset. To that end, we first compute the points $(x_{k},y_{k},z_{k})_{k=0\ldots n-1}$ of a spherical Fibonacci lattice [10]

[TABLE]

where $\varphi=\tfrac{1+\sqrt{5}}{2}$ denotes the golden ratio. To account for anisotropy of the data we generate intermediate, anisotropically scaled vectors \vec{u}_{k}=\bigl{(}\tfrac{x_{k}}{s_{x}},\tfrac{y_{k}}{s_{y}},\tfrac{z_{k}}{s_{z}}\bigr{)}. The respective anisotropy factor $\vec{s}=(s_{x},s_{y},s_{z})$ is calculated as the median bounding box size of all objects in the training images. The final unit rays $\vec{r}_{k}$ are then computed via normalization $\vec{r}_{k}=\frac{\vec{u_{k}}}{|\vec{u_{k}}|}$. The surface of a star-convex polyhedron represented by the distances $d_{k}$ is then given by its vertices $d_{k}\cdot\vec{r}_{k}$ and triangulation induced by the convex hull facets of the unit rays $\vec{r}_{k}$ (which is a convex set by definition).

We generally find that a sufficiently accurate reconstruction of the labeled 3D cell nuclei in our ground-truth (GT) images can be obtained with as few as $64$ rays. Fig. 2 shows the reconstruction fidelity for a dataset with highly anisotropic images (Parhyale, cf. Section 3) and highlights the importance of using an appropriate anisotropy factor $\vec{s}$. Note that $\vec{s}$ is automatically computed from the GT images and does not have to be chosen manually. Furthermore, Fig. 2 shows that our ray definition (Fibonacci) is more accurate than using equidistant (polar/azimuthal) distributed rays.

2.2 Model

Following [21], we use a convolutional neural network (CNN) to densely predict the star-convex polyhedron representation and a value that indicates how likely a pixel is part of an object. Concretely, for each pixel $(x,y,z)$, the CNN is trained to predict the $n$ radial distances $\{d_{k}(x,y,z)\}_{k=0\ldots n-1}$ to the object boundary as defined above and additionally an object probability $p(x,y,z)$ defined as the (normalized) Euclidean distance to the nearest background pixel (Fig. 1a). To save computation and memory we predict at a grid of lower spatial resolution than the input image, since a dense (i.e., per input pixel) output is often not necessary (this is similar to the concept of bounding box anchors in object detection approaches [19, 18]).

We use a slightly modified 3D variant of ResNet [12] as a neural network backbone111We find that using a U-Net [9] backbone leads to very similar results. to predict both the radial distances and object probabilities (Fig. 1b). In particular, we use residual blocks with $3$ convolution layers of kernel size $3\mkern-2.0mu\times\mkern-2.0mu3\mkern-2.0mu\times\mkern-2.0mu3$. Similar to [12], we start with two convolution layers of kernel sizes $7\mkern-2.0mu\times\mkern-2.0mu7\mkern-2.0mu\times\mkern-2.0mu7$ and $3\mkern-2.0mu\times\mkern-2.0mu3\mkern-2.0mu\times\mkern-2.0mu3$, but without strides to avoid downsampling. This is followed by $m$ residual blocks, where each block only performs downsampling if the spatial resolution is still higher than the prediction grid (see above); we double the number of convolution filters after each downsampling. After the last residual block, we use a single-channel convolution layer with sigmoid activation to output the per-pixel222To improve readability we will drop from now on the explicit pixel coordinate $(x,y,z)$ for both $p(x,y,z)$ and $d_{k}(x,y,z)$. object probabilities $p$. The last residual block is additionally connected to an $n$-channel convolution layer to output the radial distances $d_{k}$. Our code based on Keras/TensorFlow [8, 1] and documentation is available at https://github.com/mpicbg-csbd/stardist.

Training

Given the pixel-wise object probabilities and distances of the prediction $(\hat{p},\hat{d}_{k})$ and ground-truth $(p,d_{k})$, we minimize the following loss function (averaged over all pixels) during training:

[TABLE]

For the object loss $L_{\mathit{obj}}$ we use standard binary cross-entropy

[TABLE]

For the distance loss $L_{\mathit{dist}}$ we use the mean absolute error weighted by the object probability (active only on object pixels, i.e. $p>0$) and add a regularization term (active only on background pixels, i.e. $p=0$):

[TABLE]

This specific form was chosen to promote increased accuracy for points closer to the object centers (which eventually will be used as polyhedra center candidates).

Prediction

After the CNN predicts the radial distances $\hat{d}_{k}$ and object probabilities $\hat{p}$, we collect a set of object candidates by only considering radial distances at pixels with object probabilities above a reasonably high threshold, i.e. we only retain shapes that very likely belong to an object. Since the set of object candidates is highly redundant, we use non-maximum suppression (NMS) to obtain only one shape for every actual object in the image, as is common in object detection (e.g., [11]). Thereby, the object candidate with the highest object probability suppresses all other remaining candidates if they overlap substantially. This process is repeated until there are no further candidates to be suppressed. All remaining (i.e. not suppressed) candidates yield the final set of predicted object shapes.

2.3 Efficient non-maximum suppression

The NMS step requires to assess the pairwise overlap of a potentially large set of polyhedron candidates ($>10^{4}$). Unfortunately, computing the exact intersection volume between two star-convex polyhedra efficiently is non-trivial (in contrast to convex polyhedra). To address this issue, we employ a filtering scheme that computes as needed successively tighter upper and lower bounds for the overlap of polyhedron pairs (cf. Fig. 1c). Concretely, we compute the intersection volume of the respective i) bounding spheres (upper bound), ii) inscribed spheres (lower bound), iii) convex hulls (upper bound), and iv) kernels333The (convex) set of all points that can serve as center of the star-convex polyhedron. (lower bound). Note that iii) and iv) involve only the intersection of convex polyhedra and can thus be computed efficiently [3]. If a computed bound is already sufficient to decide whether a candidate should be suppressed or not, no further computation is carried out. Otherwise, we eventually perform the expensive but exact intersection computation by rasterization of both polyhedra. We find that this NMS filtering scheme leads to a noticeable reduction in runtime that makes StarDist-3D practical (e.g. $9$ s for a stack of size $1141\mkern-2.0mu\times\mkern-2.0mu140\mkern-2.0mu\times\mkern-2.0mu140$ with $12000$ initial candidates).

3 Experiments

We consider two qualitatively different datasets (Fig. 3) to validate the efficacy of our approach:

Worm

A subset of $28$ images used in Long et al. [15], showing DAPI-stained nuclei of the first larval stage (L1) of C. elegans (Fig. 3 left). Stacks are of average size $1157\mkern-2.0mu\times\mkern-2.0mu140\mkern-2.0mu\times\mkern-2.0mu140$ pixels with semi-automatically annotated cell nucleus instances ($15148$ in total) that underwent subsequent manual curation. We randomly choose $18$/$3$/$7$ images for training/validation/testing. Note that the images have (near) isotropic resolution.

Parhyale

A subset of recording #04 of Alwes et al. [2], showing Parhyale hawaiensis expressing Histone-EGFP (Fig. 3 right). It contains $6$ images of $512\mkern-2.0mu\times\mkern-2.0mu512\mkern-2.0mu\times\mkern-2.0mu34$ pixels with manually annotated nucleus instances ($1738$ in total). We randomly choose $3$/$1$/$2$ images for training/validation/testing. In contrast to Worm, the images are highly anisotropic in the axial direction. This dataset is more challenging due its substantially lower signal-to-noise ratio (cf. Fig. 6). Furthermore, it contains much fewer labeled training images, more akin to what is typical in many biological datasets.

3.1 Methods and Evaluation

We compare our proposed StarDist-3D approach against two kinds of methods (IFT-Watershed [16] and 3D U-Net [9]) that are commonly used for segmentation of fluorescence microscopy images. First, a classical watershed-based method that does not use machine learning. Second, a variant of the popular U-Net with and without more sophisticated postprocessing.

To evaluate the performance of all methods, we use $\mathit{accuracy(\tau)}=\frac{\mathit{TP}}{\mathit{TP}+\mathit{FN}+\mathit{FP}}$ for several overlap thresholds $\tau$. $\mathit{TP}$ are true positives, which are pairs of predicted and ground-truth nuclei having an intersection over union (IoU) value $\geq\tau$. $\mathit{FP}$ are false positives (unmatched predicted instances) and $\mathit{FN}$ are false negatives (unmatched ground-truth instances). We use the Hungarian method [14] (implementation from SciPy [13]) to compute an optimal matching whereby a single predicted nucleus cannot be assigned to multiple ground-truth instances (and vice versa). Note that a suitable value of $\tau$ depends on the biological application. For example, one would likely use a smaller $\tau<0.5$ for the purpose of counting objects, whereas intensity measurements inside each object would require more accurate shapes and thus demand a higher value of $\tau$.

StarDist-3D

We use StarDist-3D as explained in Section 2.2 with $n=96$ radial directions and $m=3$ residual blocks that start with $32$ convolution filters. We predict at a grid half the spatial resolution of the input image, except for the anisotropic Z axis of Parhyale. We use automatically computed anisotropy factors (cf. Section 2.1) of $\vec{s}=(1,1,1)$ for Worm and $\vec{s}=(1,1,7.1)$ for Parhyale. We use weights $\lambda_{d}=0.1$ and $\lambda_{reg}=10^{-4}$ for the loss function in  Eq. 1.

IFT-Watershed

The IFT-Watershed [16] is an efficient combination of maxima detection and 3D watershed segmentation. It represents an advanced classical image segmentation method that we know is being used in practice. Concretely, we use the Interactive Watershed plugin444https://imagej.net/Interactive_Watershed in Fiji [20] and perform extensive parameter tuning (such as Gaussian filter size during preprocessing and maxima detection thresholds) using the training images of each dataset.

U-Net

We train a 3D U-Net [9] to classify each pixel into background, nucleus, and also nucleus boundary, as this helps substantially to separate touching nuclei [6]. We threshold the predicted nucleus probabilities and group pixels in each connected component to obtain individual nuclei.

U-Net+

We use the same trained 3D U-Net as above, but apply more sophisticated postprocessing. Concretely, we observe improved performance by thresholding the nucleus probabilities to obtain seed regions that we grow (via 3D watershed [23]) until they reach pixels with background probability above a second threshold.

We apply random data augmentations during training, including flips, axis-aligned rotations, elastic deformations, intensity rescaling, and noise. After training, thresholds for all methods (as described above) are tuned on validation images to optimize accuracy averaged over $\tau\in\{0.3,0.5,0.7\}$.

3.2 Results

The results in Tables 1 and 4 show that StarDist-3D consistently outperforms all other methods that we compared to (note that we report the average result over $5$ trials for all trained models). The performance gap between StarDist-3D and the other methods is especially striking for dataset Parhyale, which may be explained by StarDist-3D’s shape model being especially helpful to disambiguate between neighboring nuclei in these challenging low-SNR images. In Fig. 6 we show lateral (XY) and axial (XZ) views of segmentation results for both datasets. Here, IFT-Watershed often produces imperfect boundaries and erroneous splits, particularly for dataset Parhyale. This is expected, as the watershed operation uses the input intensities alone without leveraging extracted features. U-Net tends to under-segment the image, generally producing object instances that are too small, as the use of a single threshold for the nucleus class leads to a trade-off between object size and avoidance of falsely merged objects. In contrast, U-Net+ exhibits slight over-segmentation, since a larger first threshold produces more (but smaller) objects that are then grown to yield the final instances. Finally, StarDist-3D produces superior segmentations, although it can sometimes fail to detect some nuclei (especially for dataset Parhyale). As an additional visualization, we show a 3D rendering of StarDist-3D segmentation results for both datasets in Fig. 7.

Note that we find (not shown) that the accuracy of StarDist-3D would drop dramatically (for example, from $0.593$ to $0.291$ for $\tau=0.5$) if we did not adapt the radial directions to account for the anisotropy of the nuclei shapes (Section 2.1) for Parhyale. While StarDist-3D’s lead is less pronounced for dataset Worm, this may be due to the higher signal quality of the input images and also the general abundance of labeled cell nuclei available for training and validation ($11387$ in total). In Fig. 5, we investigate how StarDist-3D and the other trained models cope with less annotated data by randomly selecting only a partial 3D image slice from each training and validation stack. Interestingly, we find that with only $4.15\%$ of the training and validation data ($472$ instances in total), StarDist-3D can for $\tau=0.5$ reach the same performance (accuracy of $0.7$) as U-Net+ with access to $100\%$ of the training and validation data.

4 Discussion

We presented StarDist-3D, an extension of [21] to detect and segment cell nuclei in volumetric fluorescence microscopy images, even when they exhibit substantial anisotropy. Our method outperformed strong watershed and U-Net baselines, yet is easy to train and use, and due to our star-convex polyhedra parameterization and efficient intersection implementation fast enough to process typical large 3D volumes. Furthermore, StarDist-3D should be generally applicable to segment objects whose shapes are well-represented with star-convex polyhedra.

Acknowledgments

We thank Frederike Alwes and Michalis Averof (IGFL) for providing Parhyale data, Dagmar Kainmüller (MDC Berlin) for worm annotation and curation, and Stephan Saalfeld and Philipp Hanslovsky (HHMI Janelia) for software assistance. Uwe Schmidt and Robert Haase were supported by the BMBF grant SYSBIO II (031L0044), and Ko Sugawara by the ERC Advanced Grant reLIVE (ERC-2015-AdG-694918).