Optical appearance of cell membrane is variable due to different size, number, and position of fluorescent signals on the focal plane. In our method, a membrane image is preprocessed with multiple steps. A fluorescent microscope produces membrane stack (red) and nucleus stack (blue) simultaneously. While nucleus channel is used to generate (nucleus-level) seeds matrix by existing methods, we obtain the cellular shapes by leveraging the membrane channel. The framework of 3DMMS can be divided into three parts, membrane image preprocessing, membrane-centered segmentation and division correction, as illustrated in Fig. 9.
Data
C. elegans was first stained with dual labelling in cell nucleus and membrane. All the animals were maintained on NGM plates seeded with OP50 at room temperature unless stated otherwise. Membrane marker and lineaging marker were rendered homozygous for automated lineaging. To improve the overall resolution, 4D imaging stacks were sequentially collected on both green and red fluorescent protein (mCherry) channels at a 1.5-min interval for 240 time points, using a Leica SP8 confocal microscope with a 70-slice resonance scanner. All images were acquired with resolutions of 512×712×70 stack (with voxel size 0.09×0.09×0.43μm). All the images were deconvoluted and resized into 205×285×70 before analysis.
Membrane image preprocessing
Statistical intensity normalization
Fluorescent images are often corrupted by noise, such as Poisson distributed incoming photos. Besides, signal intensity decreases along the z-axis because of the attenuation of laser energy. To achieve parameter generalization through the whole stack, Gaussian smoothed membrane image was adjusted by statistical intensity normalization, which balanced the intensity distribution of symmetrical slices in each stack. First, pixel intensity histogram of each slice was embedded into an intensity distribution matrix as a row. Background pixels were ignored for computational stability. An example of Gaussian smoothed intensity distribution matrix is shown in Fig. 10a. A threshold of the pixel number was applied, thus a threshold line (red in Fig. 10a) was formed across all slices. Slices at the deeper half of the stack were multiplied by the ratio of this slice’s intensity on the red line to that of its symmetrical slice. The stack intensity distribution after the adjustment is shown in Fig. 10b.
Additionally, the membrane stack was resampled to 205×285×134 with linear interpolation on the z-axis.
Hessian matrix enhancement
Cell surfaces are composed of plane components. Membrane signals can be enhanced by selecting all pixels that belong to a plane structure. We took the associate quadratic form to exploit intensity changes surrounding a pixel, and further determined its structure components. By diagonalizing the quadratic form, the Hessian descriptor is defined as
$$ {{} \begin{aligned} H\,=\,\left[ \begin{array}{ccc} \frac{\partial^{2}{I^{\mathrm{m}}}}{\partial{x^{2}}} &\frac{\partial^{2}{I^{\mathrm{m}}}}{\partial{xy}} &\frac{\partial^{2}{I^{\mathrm{m}}}}{\partial{xz}}\\ \frac{\partial^{2}{I^{\mathrm{m}}}}{\partial{yx}} &\frac{\partial^{2}{I^{\mathrm{m}}}}{\partial{y^{2}}} &\frac{\partial^{2}{I^{\mathrm{m}}}}{\partial{yz}}\\ \frac{\partial^{2}{I^{\mathrm{m}}}}{\partial{zx}} &\frac{\partial^{2}{I^{\mathrm{m}}}}{\partial{zy}} &\frac{\partial^{2}{I^{\mathrm{m}}}}{\partial{z^{2}}} \end{array} \right]\,=\,\left[\begin{array}{ccc} \vec{e_{1}} &\vec{e_{2}} &\vec{e_{3}} \end{array} \right]\left[\begin{array}{ccc} \lambda_{1} &0 &0\\ 0 &\lambda_{2} &0\\ 0 &0 &\lambda_{3} \end{array} \right]\!\left[\begin{array}{c} \vec{e_{1}}\\\vec{e_{2}}\\\vec{e_{3}} \end{array} \right] \end{aligned}} $$
(2)
where λ1,λ2,λ3 are eigenvalues with |λ1|<|λ2|<|λ3|, and \(\vec {e_{1}}, \vec {e_{2}}, \vec {e_{3}}\) are the corresponding eigenvectors. Pixels could be allocated to three structures regarding the eigenvalues: (1) when |λ1|,|λ2|<1 and |λ3|≥1, the pixel locates on a plane; (2) when |λ1|<1 and |λ2|,|λ3|≥1, the point locates on a stick; and (3) when |λ1|,|λ2|,|λ3|≥1, the point locates in a ball. So membrane surface signal can be enhanced with
$$ I^{\text{en}}(x,y,z)=\frac{|\lambda_{3}(x,y,z)|}{\max\left(|\lambda_{3}(x,y,z)|x,y,z\in{\text{stack voxels}}\right)} $$
(3)
where Ien is the stack image after enhancement.
Region filter
Preliminary experiment shows membrane based EDT (in “Membrane-centered segmentation” section) is highly dependent on the quality of binary membrane image. The region filter is designed to screen noise regions in Ien. After suppressing noise and enhancing membrane signal, we choose a threshold to convert Ien into binary image Ibn. It is composed of disconnected regions, denoted as Φ={ϕi}, some of which are noise spots. The largest connected region ϕi belongs to valid cell surface signal χ, but other regions need to be screened. Keeping noise spots would introduce erroneous cell boundaries, whereas missing valid signal results in segmentation leakages.
Herein, principal component analysis (PCA) was employed to analyze the location relationship between ϕmax and small regions in {Φ∖ϕmax}. Noise and valid regions had different influence on the Euclidean distance transformation (EDT) of the membrane surface ϕmax. The flow chart of the region filter is shown in Fig. 11. Cell surface signal was initialized as χ={ϕmax}. Following steps were repeatedly used to update χ:
-
1.
Construct zero matrix L with the same size as Ibn. Points already in ϕmax are set as 1 in L. DL denotes the EDT results on L. Similarly, after another region ϕi (green or yellow region in Figs. 11b and d) in {ϕ∖χ} is combined into L, EDT is also used to generate DL′.
-
2.
We use
$$ R=\left\{(x,y,z)|DL(x,y,z)\neq DL'(x,y,z)\right\} $$
(4)
to obtain the influenced EDT region R when we add ϕi into L.
-
3.
Use PCA to analyze the polarization features of R. Variance percentage on three directions are γ1,γ2,γ3 and γ1<γ2<γ3. The coefficient for adding ϕi into χ is measured by \(\frac {\gamma _{1}}{\gamma _{1}+\gamma _{2}+\gamma _{3}}\). Our experiments shows that if this coefficient is larger than 0.1, ϕi can be regarded as membrane signal and should be grouped into χ. Otherwise, ϕi will be ignored.
An example result is shown in Fig. 12. Filtered membrane stack Ifm is a binary image whose points in χ is positive.
Surface regression
The embryonic surface cannot be imaged completely because of a balance between the phototoxicity and signal intensity. Moreover, the stain concentration is much lower at the boundary where only one layer of the membrane exists. Incomplete surface degrades the performance of 3DMMS because of the leakage between different targets, as shown in Fig. 13b. We use surface regression to recover the boundary surface signal around the missing embryonic surface area, noted as surface cavity. In surface regression, we only modify surfaces in the cavities and this is different from the embryonic region segmentation in BCOMS.
We apply the active surface first to obtain the initial surface of the entire embryo. The smooth factor is tuned to be a large value to prevent segmented surface dropping into the cavity. From Fig. 14, we know that cavity surface can be found according to the vertical distance between the segmented embryo surface and the membrane signal Ifm. We defined a distance matrix as the same size as one slice. For the upper half surface of the segmented embryonic surface Seu, the distance matrix delineated the vertical distance between Seu and membrane signal Ifm. The distance was set to zero when there were no corresponding signals. Distance matrix was smoothed, and further thresholded using Ostu’s method [30], to construct a binary mask Rcavity. Positive masks in Rcavity indicated the location where the membrane signal should be modified with Seu. We used
$$ I^{\text{fm}}\left(x,y,S^{eu}(x,y)\right)= \begin{cases} 1, &\text{if }R^{\text{cavity}}(x,y)=1\\ 0, &\text{if }R^{\text{cavity}}(x,y)\neq 1 \end{cases} $$
(5)
to repair Ifm. Partial surfaces with positive mask were added into Ifm, shown as gray points in Fig. 13c.
Membrane-centered segmentation
Watershed segmentation is a fast algorithm to group points with different labels according to specific terrain map based on image intensity. Along the steepest descent, all pixels are classified into different catchment basin regions by tracing points down to the corresponding local minima [31], which are also termed as seeds. After watershed transformation, each region consists of points whose geodesic descent paths terminate at the same seed. The number of seeds controls the number of regions. Redundant seeds result in over-segmentation where one region is split; whereas, absent seeds lead to under-segmentation with two regions combined together. The terrain map plays a dominant role in generating region boundaries. In 3DMMS, a well-defined terrain map, combined with nucleus channel, accommodates the difficulty of lost information and membrane perception.
The nucleus image is simultaneously acquired with the membrane image, which can be used as seeds to eliminate merge-or-split mistakes. Generally, the terrain map is the linear combination of membrane intensity in nucleus-centered watershed segmentation [21,32–34]. However, it is difficult to make a tradeoff between two sources of influence on the final region boundary, as shown in Fig. 15 (combination of EDT and membrane). To overcome this problem, we combined nucleus and membrane stacks in a different way, noted as membrane-centered watershed. The nucleus stack was processed by AceTree to generate the nucleus matrix. The nucleus matrix In was constructed as
$$ I^{\mathrm{n}}=l_{i} $$
(6)
where (xi,yi,zi) and li were the nucleus location and label in the lineage, respectively. We noted Dm as the membrane-centered EDT on Ifm. Then Dm was reversed and normalized by
$$ D^{\mathrm{m}}=\frac{\max(D^{\mathrm{m}})-D^{\mathrm{m}}}{\max(D^{\mathrm{m}})} $$
(7)
The nucleus matrix In, plus a background minimum, were used as seeds for the watershed segmentation on new terrain map Dm. This map can, to a certain extent, relieve the segmentation leakage by building a ridge at the holes of the binary membrane signal, as demonstrated in Fig. 15 (membrane-centered EDT). Channel-connected cells were well separated with each other. It produces reasonable boundaries in both the blurry area and surface cavities.
Cell division revision
Two nuclei in a dividing cell would lead to a split, indicated with red lines in Fig. 16b. We resolved this problem by considering the membrane signal distribution of the interface between two cells. First, we analyzed nucleus lineage information and found out the daughter cells (or nuclei). Details on the rules of finding daughter cells can be found in [“Additional file 1”]. For each pair of daughter cells, the intensity of their interface is examined to determine whether the division has finished. The membrane-centered segmentation yields cell boundaries with the membrane signal or ridges in EDT. We calculated the average intensity of two cells’ interface to determine whether this interface located at ridges with a hole. If the interface includes a hole, the division is in process and two cells should be merged. The average intensity threshold is experimentally determined to be 40. Segmentation results after cell division correction is shown in Fig. 16c.