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 (nucleuslevel) 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, membranecentered 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.5min interval for 240 time points, using a Leica SP8 confocal microscope with a 70slice 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 zaxis 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 zaxis.
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 I^{en} is the stack image after enhancement.
Region filter
Preliminary experiment shows membrane based EDT (in “Membranecentered segmentation” section) is highly dependent on the quality of binary membrane image. The region filter is designed to screen noise regions in I^{en}. After suppressing noise and enhancing membrane signal, we choose a threshold to convert I^{en} into binary image I^{bn}. 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 I^{bn}. 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 I^{fm} 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 I^{fm}. We defined a distance matrix as the same size as one slice. For the upper half surface of the segmented embryonic surface S^{eu}, the distance matrix delineated the vertical distance between S^{eu} and membrane signal I^{fm}. 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 R^{cavity}. Positive masks in R^{cavity} indicated the location where the membrane signal should be modified with S^{eu}. 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 I^{fm}. Partial surfaces with positive mask were added into I^{fm}, shown as gray points in Fig. 13c.
Membranecentered 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 oversegmentation where one region is split; whereas, absent seeds lead to undersegmentation with two regions combined together. The terrain map plays a dominant role in generating region boundaries. In 3DMMS, a welldefined 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 mergeorsplit mistakes. Generally, the terrain map is the linear combination of membrane intensity in nucleuscentered 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 membranecentered watershed. The nucleus stack was processed by AceTree to generate the nucleus matrix. The nucleus matrix I^{n} was constructed as
$$ I^{\mathrm{n}}=l_{i} $$
(6)
where (x_{i},y_{i},z_{i}) and l_{i} were the nucleus location and label in the lineage, respectively. We noted D^{m} as the membranecentered EDT on I^{fm}. Then D^{m} was reversed and normalized by
$$ D^{\mathrm{m}}=\frac{\max(D^{\mathrm{m}})D^{\mathrm{m}}}{\max(D^{\mathrm{m}})} $$
(7)
The nucleus matrix I^{n}, plus a background minimum, were used as seeds for the watershed segmentation on new terrain map D^{m}. 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 (membranecentered EDT). Channelconnected 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 membranecentered 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.