Segmentation of 3D images of plant tissues at multiple scales using the level set method

Background Developmental biology has made great strides in recent years towards the quantification of cellular properties during development. This requires tissues to be imaged and segmented to generate computerised versions that can be easily analysed. In this context, one of the principal technical challenges remains the faithful detection of cellular contours, principally due to variations in image intensity throughout the tissue. Watershed segmentation methods are especially vulnerable to these variations, generating multiple errors due notably to the incorrect detection of the outer surface of the tissue. Results We use the level set method (LSM) to improve the accuracy of the watershed segmentation in different ways. First, we detect the outer surface of the tissue, reducing the impact of low and variable contrast at the surface during imaging. Second, we demonstrate a new edge function for a level set, based on second order derivatives of the image, to segment individual cells. Finally, we also show that the LSM can be used to segment nuclei within the tissue. Conclusion The watershed segmentation of the outer cell layer is demonstrably improved when coupled with the LSM-based surface detection step. The tool can also be used to improve watershed segmentation at cell-scale, as well as to segment nuclei within a tissue. The improved segmentation increases the quality of analysis, and the surface detected by our algorithm may be used to calculate local curvature or adapted for other uses, such as mathematical simulations. Electronic supplementary material The online version of this article (10.1186/s13007-017-0264-5) contains supplementary material, which is available to authorized users.


Segmentation of 3D images of plant tissues at multiple scales using the level set method
Kiss et al.

Mathematical details of the level set method
We have adopted the distance regularised level set evolution (DRLSE) [1]. In this formulation of the LSM a regularising term is included in the energy function in order to avoid the necessity of periodically reinitialising the LSF function. Simultaniously, despite the use of simple finite difference implementation and of relatively large time steps sufficient numerical accuracy is ensured by this formulation. More precisely, let φ : Ω → R be an LSF defined on a domain Ω having its zero level the evolving contour. Assuming that the embedding LSF function takes positive values inside the contour and negative values outside, the inward normal vector can be expressed as n = ∇φ/|∇φ|, while the scalar curvature is κ = ∇n.
The evolution of the LSF minimizes the energy functional a linear combination of the following terms.
The image term L g in this energy functional guides the evolving contour to the desired structure on the image, where δ is the Dirac delta function, so that the energy L g (φ) is the integral of the edge function g along the zero level contour of φ.
The accelerating term A g is based on a weighted area (volume in 3D) of the domain enclosed by the contour, with H the Heaviside function. Lastly, we also added a simple regularising term with the potential p(|∇φ|) = 1 2 (|∇φ| − 1) 2 , in order to prevent the numerical degradation of the LSF during evolution by maintaining its gradient close to unity without any reinitialisation following [1] To summarize, the global energy functional in Equation (8) is the linear combination of all these terms, with the Dirac delta function δ and the Heaviside function H appearing in L g and A g respectively approximated by the functions and in which we set the parameter ε to 0.5, so that the width of the approximation of δ is one pixel. So given an initial LSF φ(x, 0) = φ 0 (x), the energy functional (1) can be minimized by solving the gradient flow ∂φ ∂t = − δEε δφ , which is of the form [1]

Parameters and initialisation
In Equation (8), we fix the parameter of the image term (λ = 10), and all the other weight parameters are determined relative to this one. In particular, the regularisation parameter µ is chosen fixed, µ = 0.001 worked suitably for all our 3D tests. The most important parameters are thus the remaining weight parameters in the energy functional, α for the accelerating term and β for the smoothing term. Different choices of these parameters lead to (slightly) different solutions, their choice is made by the averted user, depending on the specific question and situation. In Section 2.1 we already illustrated the effect of these parameters on the solution of the level set method used to detect the tissue surface. Here we study the sensitivity with respect to the initialisation process.
First we created an artificial test image stack by taking mirror symmetries in a real tissue image (Supplementary Figure 1A), so that the object in the image doesn't touches the boundary of the image. We considered five different initialisations for the LSF, all computed as thresholds of the input image stack at different values, so that the initial contour is in the exterior of the tissue. Then we followed the evolution of the contours for the three representative parameter sets (α,β) already presented in Section 2.1 (Supplementary Figure 1B-D).
Our results are briefly summarised in Supplementary Figure 2E. Given a parameter set (α,β), notice the convergence to the same tissue volume regardless the initialisation. Furthermore, notice the acceleration property of α = 1: when the accelerating term is turned on, the evolution is faster and the final object volume is lower. Notice also the effect of the smoothing term β = 1: the contour is smoothed out, the evolution is decelerated at the beginning, but the final tissue volume is still lower then the case when this term is turned off. Author details