2.1. Multitask learning for membrane and nuclei detection

The membrane cue is invariably necessary to segment cells from immuno-fluorescent images⁽Reference Cao, Wong, Zhao and Yan^{4, Reference Mosaliganti, Noche, Xiong, Swinburne and Megason7–Reference Azuma and Onami10)} but as discussed in Ref. (Reference Azuma and Onami10) combining the nuclei and membrane is more conducive for cell segmentation from tissue samples. This is specifically more relevant for cardiomyocyte segmentation which demonstrates significant variation in cell morphology and noncellular structures. In the convolutional neural network community, the efficacy of multitask learning has been demonstrated in the literature⁽Reference Ramesh and Tasdizen^12, Reference Sarkar, Mukherjee, Labruyère and Olivo-Marin¹³⁾, for various image analysis applications. We exploit the multitask neural network approach to combine the membrane and nuclei information at an initial detection stage for more robust detection.

The neural network is trained on 2D slices of nuclei and membrane images. The network is designed with a partially shared encoder between the nuclei and membrane channels to mutually benefit from one another. Further, due to low resolution, spatially delineating two nuclei is a challenging task. In this scenario, the membrane cue is advantageous to spatially resolve nuclei. An additional constraint for a robust delineation of the nuclei is employed via detection of intermembranous (inverted membrane signal) regions, which would ideally denote the cytoplasmic regions of the myocardial cells. An initial convolution block encodes the features from nuclei and membrane, which are then concatenated and input to the common encoder. Two separate decoder-blocks input the original features and output from the encoder to provide pixel-wise prediction of the nuclei and membrane region. The cytoplasmic region detection is performed via another branch which uses input from the membrane features and nuclei-membrane encoder output. The detailed network architecture is shown in Figure 4. The neural network is trained on the combined loss function from the three channels,

(1) $$ {L}_m\left({I}_m,{\hat{I}}_m\right)+{L}_d\left(\left({I}_{im}\right),\left({\hat{I}}_{im}\right)\right)+{L}_d\left({I}_n,{\hat{I}}_n\right). $$

Here $ {I}_m $ , $ {I}_n $ are the ground truth for membrane and nuclei, whereas $ {\hat{I}}_m $ , $ {\hat{I}}_n $ are the corresponding predictions. $ {I}_{im} $ is the intermembranous regions denoted by $ \left(1-{I}_m\right) $ , which is precomputed and augmented (nonmyocardial regions are removed) during the ground truth generation process.

Dice loss $ {L}_d $ is used for nuclei and inverted membrane detection whereas mean squared loss, $ {L}_m $ is used for membrane detection. The final refined 2D nuclei prediction is obtained as an overlap of nuclei and intermembranous region prediction. This aids in discarding the nonmyocardial cell nuclei and noisy detection. Sample prediction results of the nuclei and membrane are shown in Figure 3. As can be seen, the enhanced membrane signal and the desired nuclei are obtained with significantly good precision.

Figure 3. Deep learning results of 2D detection of membrane (green) and nuclei (blue) (from three different ventricular regions) are shown in (b,d,f). The corresponding original images are shown in (a,c,e).

Figure 4. Proposed neural network architecture is shown in the figure. The size of convolution kernel and the channel (height × width × channel size) for each layer are shown in the above figure. The channel size in each layer is also marked above each block in the diagram.

2.2. 3D nuclei detection via association

The 3D nuclei structure can be obtained from the 2D detections. This can be achieved by associating the 2D detections axially to construct the 3D nuclei. A naïve approach is to employ local heuristics (such as position, size, etc.)⁽Reference Stegmaier, Amat and Lemon^9, Reference Wang, Sarkar, Aziz, Vaccari, Gahlmann and Acton¹⁴⁾ to solve the association problem, but often lead to irregular and erroneous 3D structures. In Refs. (Reference Mukherjee, Basu, Condron and Acton15–Reference Xu, Su, Zhu, Guan and Zhang17), the authors present more sophisticated formulations for 3D segmentation by association. Xu et al. ⁽Reference Xu, Su, Zhu, Guan and Zhang¹⁷⁾ employ a graph-based global association to associate 2D detections for 3D neuron segmentation. In Ref. (Reference Mukherjee, Basu, Condron and Acton15), Mukherjee et al. used a graph-based method which solves the minimum spanning tree problem to associate disjoint objects and achieve 3D neuron segmentation. An attraction force-based association of disjoint objects in 3D is designed by Ref. (Reference Mukherjee, Condron and Acton16). These aforementioned methods⁽Reference Mukherjee, Basu, Condron and Acton^15, Reference Mukherjee, Condron and Acton^16, Reference Pop, Dufour and Le Garrec¹⁸⁾ have been designed for segmenting neurons, which demonstrate a tubular tree-like geometry unlike cell nuclei that are more convex in structure. A robust data association for 3D nuclei segmentation forms the core of our method since this serves as the initialization of an active mesh model for 3D cell segmentation. Here, we design a novel graph-based association approach to detect 3D nuclei while maintaining the convex geometry and certain biological constraints. In this section, we formally define the mathematical formulation of the association problem.

Let $ \mathcal{I} $ be the 3D volumetric image of nuclei and $ {\mathcal{I}}_z $ be a 2D image ( $ z th $ slice of the volumetric image $ \mathcal{I} $ ) in the domain $ \left[0,1\right] $ . Let $ \Phi ={\left\{{o}_z^j\right\}}_{z=1\dots {N}_s,\hskip0.24em j=1\dots {n}_z} $ a set consisting of binary object on 2D images and $ {o}_z^j\cap {o}_z^k=\mathrm{\varnothing} $ _, where $ j, k\in \left\{1\dots {n}_z\right\} $ . Here, N_s denotes the number of 2D slices and n_z denotes the number of 2D objects in zth slice. $ \Psi ={\left\{{\psi}_p\right\}}_{p=1\dots {N}_3} $ , where N ₃ is the total number of 3D objects in the volumetric image $ \mathcal{I} $ and $ {\psi}_p\cap {\psi}_q=\mathrm{\varnothing} $ . For 3D nuclei detection, we seek to find a function $ \mathcal{F}:\Phi \mapsto \Psi $ . The mapping, $ \mathcal{F} $ , is accomplished by solving a constrained optimization problem, which simultaneously selects a 2D object and a link connecting two 2D objects from consecutive z-slices. The problem is defined as follows:

(2) $$ {\displaystyle \begin{array}{l}\underset{Y_{i, j}}{\min}\sum \limits_{i, j}\left\{{C}_{i, j}{Y}_{i, j}+{C}_{s, i}{Y}_{s, i}+{C}_{t, i}{Y}_{t, i}+{C}_i{Y}_i+\sum \limits_k{Q}_{i, j, k}{Y}_{i, j}{Y}_{j, k}\right\}\\ {} s. t.\sum \limits_{i, j\in \varepsilon}{Y}_{i, j}=\sum \limits_{j, k\in \varepsilon}{Y}_{j, k};{Y}_{i, j}\in \left\{0,1\right\}\hskip0.24em {Y}_i\in \left\{0,1\right\}.\end{array}} $$

Here a graph structure is imposed on the ensemble of 2D detections. Each 2D nuclei detection $ {o}_z^i $ is considered as a node $ {x}_i $ in the graph (Figure 5a). The connectivity between nodes is defined such that $ {x}_i $ in stack $ z $ is connected to a node $ {x}_j $ in the subsequent stack $ z+1 $ . Each node in the graph is connected to a source, s and a sink, t. The edge between two nodes x_i and x_j is denoted by an indicator function $ {Y}_{i, j} $ , that is, Y_i,j = 1 if an edge is present between two 2D objects and 0 otherwise. $ {Y}_{i, j} $ determines the connectivity between 2D detections to reconstruct the 3D nuclei. The nucleus connectivity is penalized by a weight $ {C}_{i, j}, $ which is computed based on how likely $ \mathcal{F}\left({o}_z^i\right)=\mathcal{F}\left({o}_{z+1}^j\right) $ . The link between a source and a node is denoted by $ {Y}_{\hskip-1pt s, i} $ and $ {Y}_{i, t} $ with connectivity weights $ {C}_{s, i} $ and $ {C}_{i, t} $ respectively. $ {C}_i $ is defined as the likelihood of a 2D object belonging to any flow path. The second order cost $ {Q}_{i, j, k} $ employed in the optimization regulates if nuclei from three subsequent stacks belong to the same 3D object. The optimal solution is obtained by sending the maximum flow from source to sink while minimizing overall cost. The 2D objects in one flow path from source to sink can be inferred as one 3D object. A detailed description of the optimization problem solution is given in Section 3.

Figure 5. The schematic for the directed graph connecting the nuclei is shown in (a). Subfigure (b) shows a schematic for the possible conflicts which need to be handled in the mathematical formulation.

This optimization problem satisfies three main properties necessary to reconstruct 3D nuclei from 2D detection.

Here $ \mathcal{N} z $ denotes the neighboring 2D slice. The constraint $ {\sum}_{i, j\in \varepsilon}{Y}_{i, j}={\sum}_{j, k\in \varepsilon}{Y}_{j, k} $ in equation (2) is essential in satisfying this property. The constraint indicates that the number of incoming edges to a node equals the number of outgoing edges. When the total number of incoming edges equals one it satisfies the property that a 2D detection can be connected to only one 2D detection in the subsequent stack.

The second-order cost, $ {Q}_{i, j, k} $ determines when Y_i,j = 1 and Y_j,k = 1 simultaneously and thus $ {x}_i $ , $ {x}_j $ , and $ {x}_k $ of stacks z, z + 1, and z + 2 belong to the same 3D object. The optimization details to solve for $ \mathcal{F} $ using equation (2) are detailed in Section 3.

2.3. Active mesh cell segmentation

The 3D nuclei detection provides an approximate location of nuclei. In order to obtain final cell segmentation, the membrane information is essential. Edge-based deformable models⁽Reference Pop, Dufour and Le Garrec^11, Reference Dufour, Thibeaux, Labruyere, Guillen and Olivo-Marin¹⁹⁾ are more appropriate choice for segmentation in this scenario. This is specifically advantageous in our case, since this formulation allows the mesh to coincide with image edge features and simultaneously forbid overlapping of meshes⁽Reference Dufour, Thibeaux, Labruyere, Guillen and Olivo-Marin¹⁹⁾. The active mesh model implemented in the bioimage analysis software Icy⁽Reference De Chaumont, Dallongeville and Chenouard²⁰⁾ is used for this purpose. The active contour evolution parameters are estimated in an automated manner as described in Refs. (Reference Pop, Dufour and Le Garrec11) and (Reference Dufour, Thibeaux, Labruyere, Guillen and Olivo-Marin19). We initialize the active mesh at the detected nuclei and evolve using the 3D membrane image $ {\mathcal{I}}_M $ . Each 2D slice in $ {\mathcal{I}}_M $ is multitask neural network membrane prediction $ {\hat{I}}_m $ .

3.1. Update Z

The solution for $ \mathbf{Z} $ can be obtained by first taking the derivative of the augmented Lagrangian with respect to $ \mathbf{Z} $ and setting it to zero. The derivative of cost function in equation (5) with respect to Z is given as

(6) $$ \underset{\mathbf{Z}}{\min }{\mathbf{Z}}^T\mathbf{QY}+{\unicode{x03BB}}_1\mid \mathbf{Z}-\mathbf{Y}\mid +{\unicode{x03BB}}_2\mid \mathbf{Z}-\mathbf{1}\mid +{\unicode{x03BB}}_3\mid \mathbf{Z}-\nu \mid +\frac{\mu}{2}\left(\parallel \mathbf{Z}-\mathbf{1}{\parallel}^2+\parallel \mathbf{Z}-\mathbf{Y}{\parallel}^2+\parallel \mathbf{Z}-\nu {\parallel}^2\right). $$

By setting the derivative of equation (6) to 0 and setting $ \beta =\frac{1}{3\mu} $ , we obtain a closed-form solution of Z as

(7) $$ \mathbf{Z}=\beta \left(\left(\mu -\mathbf{Q}\right)\mathbf{Y}+\mu \left(1+\nu \right)-{\unicode{x03BB}}_1-{\unicode{x03BB}}_2-{\unicode{x03BB}}_3\right). $$

3.2. Update Y

The solution for Y is obtained by solving the augmented Lagrangian with the constraints on Y as in equation (3):

(8) $$ {\displaystyle \begin{array}{l}\underset{\mathbf{Y}}{\min}\hskip0.24em {\mathbf{C}}^T\mathbf{Y}+{\mathbf{Z}}^T\mathbf{QY}+{\unicode{x03BB}}_1\mid \mathbf{Z}-\mathbf{Y}\mid +\frac{\mu}{2}\parallel \mathbf{Z}-\mathbf{Y}{\parallel}^2\\ {}\hskip0.24em \mathrm{s}.\mathrm{t}.\sum \limits_{i, j\in \varepsilon}{Y}_{i, j}=\sum \limits_{j, k\in \varepsilon}{Y}_{j, k};{Y}_{i, j}\in \left\{0,1\right\};{Y}_i\in \left\{0,1\right\}.\end{array}} $$

Since Y_i,j is a binary variable, one can assume $ {Y}_{i, j}^2={Y}_{i, j} $ . We can then rewrite $ \parallel {Z}_{i, j}-{Y}_{i, j}{\parallel}_2^2={Y}_{i, j}-2{Y}_{i, j}{Z}_{i, j}+{Z}_{i, j}^2 $ . By imposing this approximation⁽Reference Chari, Lacoste-Julien, Laptev and Sivic²⁵⁾, $ \parallel \mathbf{Z}-\mathbf{Y}{\parallel}^2\simeq {\left(1-2\mathbf{Z}\right)}^T\mathbf{Y}+\parallel \mathbf{Z}{\parallel}^2 $ . The optimization problem for Y can now be written as

(9) $$ {\displaystyle \begin{array}{l}\underset{\mathbf{Y}}{\min}\hskip0.24em \left({\mathbf{C}}^T+{\mathbf{Z}}^T\mathbf{Q}-{\unicode{x03BB}}_1+\mu {\left(1-2\mathbf{Z}\right)}^T\right)\mathbf{Y}\\ {}\hskip0.24em s. t.\sum \limits_{i, j\in \varepsilon}{Y}_{i, j}=\sum \limits_{j, k\in \varepsilon}{Y}_{j, k};{Y}_{i, j}\in \left\{0,1\right\};{Y}_i\in \left\{0,1\right\}.\end{array}} $$

The optimal flow Y, is obtained by solving equation (9) as a min-cost network flow optimization using the push-relabel algorithm⁽Reference Goldberg and Tarjan^22, Reference Zhang, Li and Nevatia²³⁾. The pseudo-algorithm of the optimization problem is given in Algorithm 1.

3.3. Graph weight selection

In addition to formulating the optimization method, a robust design of the link costs is an important aspect of the data association problem. For our problem, the link cost connecting the 2D nuclei detections $ {o}_z^i $ and $ {o}_{z+1}^j $ is designed taking into account three different criteria: spatial distance between 2D objects, the change in nuclei size, and a biologically inspired membrane-overlap condition. The spatial distance is simply obtained as the Euclidean distance between two nuclei given as $ {c}_{i j}^c=\frac{\parallel c\left({o}_z^i\right)- c\left({o}_{z+1}^j\right){\parallel}_2^2}{2{\sigma}_c^2} $ , where $ c\left({o}_z^i\right) $ is the centroid of nuclei $ {o}_z^i $ and $ {\sigma}_c $ constrains the spatial distance (axially) between two nuclei centers. The size constraint is imposed by the change in area of two nuclei given by $ {c}_{i j}^a=\frac{\parallel a\left({o}_z^i\right)- a\left({o}_{z+1}^j\right){\parallel}_2^2}{2{\sigma}_a^2} $ , $ a\left({o}_z^i\right) $ is the area of nuclei detection $ {o}_z^i $ and $ {\sigma}_a^2 $ controls the allowed change in the area. This criterion imposes that a significant increase (or decrease) in 2D nuclei detection size denote they are part of different 3D objects, that is, $ \mathcal{F}\left({o}_z^i\right)\ne \mathcal{F}\left({o}_{z+1}^j\right) $ (or false positives in 2D detection).

The biologically inspired membrane-overlap criterion states that even if $ {o}_z^i $ and $ {o}_{z+1}^j $ are spatially close, $ \mathcal{F}\left({o}_z^i\right)\ne \mathcal{F}\left({o}_{z+1}^j\right) $ , if cell membrane is present between them. We ensure this criterion by designing a membrane cross-over cost defined by,

(15) $$ {c}_{ij}^m=\frac{1}{N}\sum \limits_{k= p\left( i, j\right)-\frac{N}{2}}^{p\left( i, j\right)+\frac{N}{2}}{e}^{-\frac{{\left( k- p\left( i, j\right)\right)}^2}{2{\sigma}_m^2}} f\left({M}_{k, z},{M}_{k, z+1}\right)\Big). $$

Here $ f\left({M}_{., z},{M}_{., z+1}\right) $ is the interpolation of membrane stack between z and z + 1. N is the number of pixels within the window centered at $ p(.) $ . $ p\left( i, j\right) $ is the mid-point of the line joining centroids of nuclei $ {o}_z^i $ and $ {o}_{z+1}^j $ . The cost of linking two nuclei, C_i,j is given as the sum of $ {c}_{ij}^c $ , $ {c}_{ij}^a $ , and $ {c}_{ij}^m $ . The pairwise cost, Q_i,j,k is formulated as the sum of C_i,j and C_j,k to ensure the co-occurrence of three nodes in a single flow. This implies that when the link costs concur with the above criteria for three subsequent nuclei $ {o}_z^i $ , $ {o}_{z+1}^j $ , and $ {o}_{z+1}^k $ , then $ \mathcal{F}\left({o}_z^i\right)=\mathcal{F}\left({o}_{z+1}^j\right)=\mathcal{F}\left({o}_{z+2}^k\right) $ . This also ensures a local smoothness in the linked path.

Sample data association results for 3D nuclei reconstruction are shown in Figure 6. It can be observed from the figure, the nuclei association (overlaid on neural network detection results, in magenta) can delineate nuclei axially even in scenarios where visually identifying the separation is challenging.

Figure 6. Nuclei association results. First and third column shows the original images. Second and fourth column shows nuclei association results (magenta) overlaid on neural network detection of the membrane (green) and nuclei (blue).

4.1. Data description

The mouse hearts were dissected from embryos at 10.5 days after fertilization (E10.5). They were fixed in 4% paraformaldehyde and labeled as a whole mount by immuno-staining with anti-Scrib and anti-Cadh2 antibodies to detect cardiomyocyte membranes. Nuclei were counter-stained with Hoechst. Whole-mount hearts were cleared with CUBIC approach⁽Reference Le Garrec, Domnguez and Desgrange²⁶⁾ and mounted on 0.4% agarose in Reagent-2 in a capillary. Multichannel 8-bit images were acquired on a Z.1 (Zeiss) light-sheet microscope with a 20×/1.0 objective with voxel resolution is 0.63 × 0.63 × 1.26 μm.

4.2. Neural network training

We obtain 56—2D cropped images (dimension 255 × 255) from the 3D scan of an embryo ventricular region. They are manually annotated to create nuclei and membrane ground truth. The 2D ground truth was generated in a semi-automated manner using Icy, a bio-image analysis software⁽Reference De Chaumont, Dallongeville and Chenouard²⁰⁾. Ground truth for nuclei was generated by drawing manual contours. The membrane images are first enhanced using a hessian enhancement scheme and manual thresholding to obtain binary detection. An in-painting tool⁽Reference De Chaumont, Dallongeville and Chenouard²⁰⁾ is used to connect (disconnect) membrane regions to preserve continuity (remove noisy structures) of the structures. Finally, we compute the inverse of membrane ground truth and apply area-threshold to discard nonmyocardial regions. This forms the ground truth for the intermembrane prediction branch. We perform a qualitative and quantitative comparison with four state-of-the-art methods to demonstrate the efficacy of our method. A dataset consisting of 500 images is created by employing various augmentations, such as horizontal and vertical flips, rotation, variation in brightness.

The neural network is trained for 100 epochs with a learning rate of 10⁻³ and a learning rate decay of 0.9 at every 50 epochs. Adam Optimizer⁽Reference Kingma and Ba²⁷⁾ is used to train the network. Details of the network architecture are shown in Figure 4. The filter dimensions (height and width) and the number of channels (Channel Size) at each convolution layer are shown in the aforementioned figure (above each convolution step). The channel size at each step is also marked above each block in the schematic. The dotted blocks denote the concatenation of channels. Each convolution layer consists of ReLU activation followed by Maxpool (size 2 × 2) for the encoder section and an upsampling (scale factor 2) in the decoder section. Sigmoid activation layer is used in the final convolutional layer to obtain a smooth and scaled output.

4.3. 3D segmentation and validation

3D segmentation using Ca3d pipeline is performed on six crops of size 255 × 255 × 70, extracted from the ventricular region of the volumetric image of two different embryos. From each crop, the 2D slices of nuclei and membrane channel are first used to obtain 2D nuclei, membrane and intermembranous segmentation using the multitask neural network. As discussed earlier, the final refined 2D nuclei prediction is obtained as an overlap of nuclei and intermembranous region predictions. The 3D nuclei detections are obtained via the proposed association approach. Before applying active contour, the 3D nuclei of relatively small size are identified and discarded. This ensures errors due to nuclei over-segmentation are reduced in the final step. Due to the size and shape variations (resulting from intensity inhomogeneity), the terminal regions of 3D nuclei (axially) may lead to multiple associated fragments causing nuclei over-segmentation. As we seek to employ a single nuclei detection during active contour initialization for individual cell segmentation, we employ, a simple morphology-based threshold approach to discard relatively smaller 3D nuclei. Only the associated 3D nuclei which consist of 2D detections comprising of seven or more Z-slices are used for active contour segmentation. This is determined from expert annotation of nuclei. 2D nuclei regions (from neural network) of the cardiomyocytes are merged along the axial direction manually by a developmental biologist to obtain 3D nuclei structures (this is also used to compute semi-manual 3D cell segmentation for validation as discussed later in Section 4.4). From these annotations, minimum size (elongation) of the nuclei was determined as 12 μm. With Z-resolution 1.26 μm and keeping into account possible errors ensuing from neural network and the 3D association, 9 μm(1.26 × 7 μm) is selected as threshold of nuclei size. 3D nuclei objects below this size are discarded before employing active contour segmentation. This confirms smaller fragments, specifically resulting from over-segmentation of nuclei, are discarded before active contour initialization.

Due to the lack of ground truth in 3D, we first perform manual validation on a part of the segmented data. The segmentation result on the rest of the data is statistically compared with the validation for different cell morphological parameters. We also compare our method with state-of-the-art methods both qualitatively and quantitatively. The details of the validation and comparison methods are given in the following subsections.

4.4. Comparison with competing methods

4.4.2. 3D segmentation validation

As we have mentioned earlier, obtaining ground truth for individual cells from the tissue is a time consuming and laborious task even for experts in the field. Consequently, we compare the cell morphological values with two types of validation methods—manual and semi-manual cell segmentation.

1. i. Manual evaluation: The manual evaluation was performed by looking at individually segmented cells and original images by experts to determine the segmentation quality. The segmentation accuracy was measured based on two different criteria: (a) cell-overlap: axially and spatially nonoverlapping cells and (b) approximate cell size. The nonoverlapping criterion considers if the segmentation contour overlaps multiple cells either axially or spatially. The cell size criterion looks at the number of Z-slices comprising a cell and cell elongation. Each of the segmented cells is given a high-confidence or low confidence score based on the above two criteria. Manual evaluation is performed on randomly selected 300 cells for three different volumetric crops. The frequency distribution of the evaluated cells based on high-confidence value using cell size and cell overlap criterion are individually provided in Figure 7. The distribution of the cell morphological parameters, when both criteria have high confidence value simultaneously, is also shown in Figure 7 (in green). Although the distributions are almost similar to individual criteria, high confidence in both provides a more precise range for cell parameters. We employ this combined high-confidence range to evaluate the efficiency of Ca3D and compare it with state-of-the-art methods.

2. ii. Semi-manual 3D cell segmentation: A semi-manual method is used to generate cell segmentation for comparison and evaluation. In this process, the 2D nuclei detection (from neural network) of the cardiomyocytes are merged along the axial direction manually by a developmental biologist. The shape and size change of detected nuclei, spatial displacement, presence of cell membrane are evaluated manually during the merging procedure. The merged nuclei act as the initialization of active meshes which is evolved using the detected membrane image are to reconstruct individual cells in 3D. This also provides a confidence measure for the 3D nuclei reconstruction problem⁽Reference Sarkar, Darby, Foucambert, Meilhac and Olivo-Marin³⁰⁾.

Figure 7. The distribution of cells within a valid range for each cell morphological parameter using different validation criteria.

4.4.3. Quantitative comparison with cell morphology distribution profile

The methods are compared with the above manual and semi-manual validation techniques with respect to various cell morphological parameters: the cell elongation, cell-width, cell elongation ratio, surface area of cells. Any one of these cell parameters is not sufficient to assess the efficacy of a method, but if the performance is uniform over the different parameters, the efficacy is then definitely convincing. The mean and standard deviation of the parameters of segmented cells using Ca3D, RACE⁽Reference Stegmaier, Amat and Lemon⁹⁾, 3DMMS⁽Reference Cao, Wong, Zhao and Yan⁴⁾, ShapeMetrics⁽Reference Takko, Pajanoja, Kurtzeborn, Hsin, Kuure and Kerosuo⁸⁾, and Pop et al. ⁽Reference Pop, Dufour and Le Garrec¹¹⁾ are presented in Table 1. In addition to that, the percentage of over- and undersegmented cells are also presented in Table 1.

Table 1. Percentage of over and under segmentation of cells based on cell morphological parameters based on high-confidence validated cells.

It is observed that for different cell parameters, Ca3D performs evenly in comparison to the other methods. For cell elongation, although the mean value of 3DMMS and ShapeMetrics is closer to the manually validated value, a significantly high standard deviation is observed. The mean values for RACE are significantly higher and that of Pop et al. are lower compared to the validation data. This observation is also apparent from the distribution of cells in Figure 8a. The distribution of cell population pertaining to cell elongation parameter for Ca3D is comparable to the validation data in comparison to the other methods. A similar observation is made for cell width. While Ca3D and 3DMMS have similar mean values of cell width, 3DMMS demonstrates a higher standard deviation. Pop et al. and RACE on the other hand have significantly smaller population mean demonstrating a shift in the cell population density from that of the validation data (Figure 8b). The combined effect of the cell elongation and cell width is reflected in the cell elongation ratio (=cell elongation/cell width). The high cell elongation and low cell width lead to a high value of cell elongation ratio for RACE leading to a more flat distribution (Figure 8c). In comparison, Ca3D demonstrates a distribution, which is concurrent with the validation data and closely followed by Pop et al. But, Pop et al. demonstrate a lower cell population mean for surface area. Contrarily, 3DMMS has a similar population mean as the validation data but a high standard deviation leading to higher over-sized and under-sized segmented cells.

Figure 8. The distribution of the cells within a valid range for different cell morphological parameter-(a.) cell elongation, (b.) cell width, (c.) elongation ratio and (d.) surface area are shown in this figure.

To quantify the distribution presented in Figure 8b, we compute the KL divergence measure and Bhattacharya distance between the parameter distribution of each comparison method and that of the validation methods. A lower value of the distance measure implies that the distribution is similar. It is noted from Table 1, that Ca3D, for all the cell parameters concur with that of manual validation (7th and 8th column) as well as semi-automatic validation (9th and 10th column).

The percentage of oversegmented, undersegmented, and within-range cells are also presented in Table 1. The range for accurate cell size based on the aforementioned parameters is given as μ_v ± k_v × σ_v. Here, μ_v and σ_v are the mean and standard deviation of the validated cells. k_v is decided such that 90% of the validated cells are within this range. It is observed, with Ca3D we obtain highest percentages of cells, which lie within the valid range for three of the cell parameters. The number of cells within the correct range for surface area is highest for RACE⁽Reference Stegmaier, Amat and Lemon⁹⁾, but percentage of under-sized cells with respect to cell width is comparatively higher, which implies that detected cells are thinner and longer in structure. This is also evident from the elongation ratio values, where RACE⁽Reference Stegmaier, Amat and Lemon⁹⁾ demonstrates the highest percentage of over-sized cells. Lowest percentage of over-sized cells is obtained with Pop et al. ⁽Reference Pop, Dufour and Le Garrec¹¹⁾, but percentage of undersized cells is significantly high. This is pertaining to the fact that intensity-based clustering possibly detects fragments of nuclei leading to the detection of smaller-sized cells. 3DMMS⁽Reference Cao, Wong, Zhao and Yan⁴⁾ and ShapeMetrics⁽Reference Takko, Pajanoja, Kurtzeborn, Hsin, Kuure and Kerosuo⁸⁾ both demonstrate significantly high percentages of both over-sized and under-sized cells in comparison to Ca3D.

The qualitative evaluation by visual inspection is shown in Figures 9 and 10. The 3D segmented image viewed from three viewing planes is presented in Figure 9. As observed from the images, RACE⁽Reference Stegmaier, Amat and Lemon⁹⁾, Pop et al. ⁽Reference Pop, Dufour and Le Garrec¹¹⁾, and ShapeMetrics⁽Reference Takko, Pajanoja, Kurtzeborn, Hsin, Kuure and Kerosuo⁸⁾ demonstrate significant erroneous cell segmentation involving partial detection, significantly low number of detected cells, ambiguous cell boundary. Cell segmentation obtained using 3DMMS⁽Reference Cao, Wong, Zhao and Yan⁴⁾ identifies cell which aligns with the cell-membrane but generates over-segmented cells spatially as well as axially. 3D visualization of cell segmentation on a sample volumetric image and 3D images of some individual cells obtained for each of these methods are shown in Figure 10. Significantly, nonconvex and noisy cell-boundary is observed for ShapeMetrics⁽Reference Takko, Pajanoja, Kurtzeborn, Hsin, Kuure and Kerosuo⁸⁾, 3DMMS⁽Reference Cao, Wong, Zhao and Yan⁴⁾, and RACE⁽Reference Stegmaier, Amat and Lemon⁹⁾. Pop et al. ⁽Reference Pop, Dufour and Le Garrec¹¹⁾ exhibit a smoother cell boundary owing to the application of active mesh model. The superiority in segmentation quality of Ca3D is observed from both Figures 9 and 10 where segmented cells are compact-sized, with smoother boundary complying with the cell membrane while maintaining approximately convex shape.

Figure 9. 3D visualization of segmented cells in planar view for two sample crops extracted from ventricular and interventricular regions. The original raw images are shown in (a.) The segmentation method for Ca3D (in b.) and state-of-the-art methods (c.) PoP et al. ⁽Reference Pop, Dufour and Le Garrec¹⁸⁾, (d.) RACE⁽Reference Stegmaier, Amat and Lemon⁹⁾, (e.) 3DMMS⁽Reference Cao, Wong, Zhao and Yan⁴⁾, and (f.) ShapeMetrics⁽Reference Takko, Pajanoja, Kurtzeborn, Hsin, Kuure and Kerosuo⁸⁾ are shown for visual comparison.

Figure 10. 3D visualization of segmented cells sample crops extracted from ventricular. Sample segmented cells are shown on the left side of each crop. The segmentation method for Ca3D (in b.) and state-of-the-art methods (b.) PoP et al. ⁽Reference Pop, Dufour and Le Garrec¹⁸⁾, (c.) RACE⁽Reference Stegmaier, Amat and Lemon⁹⁾, (d.) 3DMMS⁽Reference Cao, Wong, Zhao and Yan⁴⁾, and (e.) ShapeMetrics⁽Reference Takko, Pajanoja, Kurtzeborn, Hsin, Kuure and Kerosuo⁸⁾ are shown for visual comparison.

4.5. Ablation study

We also present a comparison for ablation study for the convolutional neural network. The multitask neural network has three different task (nuclei, membrane, and intermembranous region) detection objectives. The ablation study is performed by using different task combinations to train the neural network. We employ three different tasks combinations: (a) nuclei, (b) nuclei and membrane, and (c) nuclei and inverted membrane. For the first and third task combination, we use 3D enhanced membrane (as in Pop et al.) for the active contour evolution. Table 2 provides the mean and standard deviation values along with percentage of over- and under-sized cells obtained for each of the neural network task combinations with respect to the validation data. For cell elongation and cell width, a lower percentage of over-sized cells is obtained using nuclei and membrane tasks but demonstrates a higher percentage of under-sized cells. The same setting also demonstrates a high percentage of over-sized cells pertaining to ell-elongation ratio. The error possibly occurs due to the consideration of nonmyocardial cells which are discarded by Ca3D using the combination of nuclei and inverted-membrane function. It is observed from the table, that under different cell parameters, Ca3D demonstrates a higher percentage of correctly sized cells in comparison to the three other task settings of the neural network.

Table 2. Ablation study: % of over and under segmentation of cells based on cell morphological parameters and 90% of high-confidence validated cells for different training of the neural network.