The phenotyping of three-dimensional root system architecture traits is necessary for effective genetic studies that will identify connections between root traits and their genotypic control, which in turn may unlock the potential of roots for crop improvement. However, the measurement of these traits proves difficult; unlike in above-ground work on shoots where manual measurements can be made, any manual measurement of root system architecture (RSA) traits results in destruction of the plant because it is displaced from the soil. That being said, it remains that for many contexts excavation for manual measurement, the ‘shovelomics’ method, still remains the most reliable method for RSA [1, 5, 23].
Other approaches include growing plants in gels [4, 28] and hydroponic systems ; the root shape is extracted using a turntable system with a shape from silhouette approach. These systems allow the plant to continue growing, and its root system to be examined over time [21, 22]. One concern with acquiring RSA traits from plants grown in gels or liquids in that these materials influence the RSA in ways that soil does not , and there are constraints to using gel-based systems, such as the size and maturity of the plant. In addition, the frictional and nutritional properties of the gel growing medium does not exactly simulate the environment plant roots experience in soil-based media.
For these reasons, the use of X-ray computed tomography (CT) has emerged as a way to visualize the structure of actively-growing plant root systems in soil or soil-like media [9, 13, 14, 16], with the possibility of capturing root growth over time because the technique is nondestructive. In order to quantify the root’s shape, the root regions of the X-ray CT regions must be segmented from the non-root regions, which is extremely arduous to do manually for complex root systems. Once the segmentation is complete, features of the root systems can be measured using metrics within the root community and this data is then used for various physiological and genomic studies.
In this paper, we provide a three-dimensional level set method adapted to the root segmentation problem in X-ray CT, specifically when the volumes are large (, or voxels). Our method requires minimal, non-exact initialization of a selection of image slices, of various orientations, from the X-ray CT volume. The level set formulation is derived from Rousson and Deriche [18, 19], but instead of modeling the whole CT volume, only regions close to the root class voxels are represented by the root and non-root distributions. We introduce some modifications specific for dealing with root systems as well as the size of the datasets we encounter. Using a narrow-band approach to level sets, we denote active regions of the contour using an occupancy grid, where the dimension of each grid cube edge is at least the size of the band. This grid is used for the updating of the distributions as well as the distance map computations. Each voxel is permitted to be a member of the active region of the contour until a threshold is exceeded, at which point the voxel becomes a permanent member of the static contour. This process restricts computational resources to the active region, and since the roots are thin, prevents previously discovered root regions from being smoothed away. The method is demonstrated on five datasets, with three plant species: soybean (Glycine max), maize (Zea mays), and cassava (Manihot esculenta) of various different maturity levels and root thicknesses, and grown in three different media: expanded clay, turface, and potting mix. Specifically, our contributions are:
General algorithm to segment root systems from non-root regions from a range of media with a varied texture and size characteristics in X-ray CT volumes.
A level set method that focuses computational resources on the active region of the contour through use of an occupancy grid and contour age parameter.
Fast exact narrow band distance map computations within the positive occupancy grid regions by an alteration of Meijster et al.  distance transform algorithm.
Level sets have been in use for some time in the medical community to segment items of interest in the human body for quantification or further review from CT image volumes, such as skulls, lung nodules, and portions of the heart ([10, 8, 20, 27] respectively). In the medical community, segmentations are often desired in a clinical setting, so one concern is the computational time involved with level set implementations. Lefohn et al.  implemented level sets on the GPU to mitigate this concern. Wang et al. [24, 25] proposed a front evolution strategy called coherent propagation, in which the oscillation of the front is reduced resulting in reduced computational time. The coherent propagation strategy was implemented in a bone segmentation setting in .
The prior work on segmenting X-ray CT images of roots fall into two groups: those based on threshold and/or morphological techniques within established software packages [9, 16] and those that use a two-dimensional level set formulation to track root regions throughout the scan volume [13, 14]. All of the prior work operations on X-ray micro-computed tomography. The method described in Metzner et al. 2015  specifies thresholds for small, medium, and large roots, and then uses connected components and morphological operations to locate root regions and a filter to remove air or non-root regions. Flavel et al. 2017  created a plugin for ImageJ, Root1, that normalizes the greylevels between two scans, and then performs adaptive thresholding, filtering and morphological operates to locate the root regions. Mairhofer et al. 2012 with RooTrak  use a level-set variation to track roots from the top of the pot to the bottom. This is done by the user initializing the root regions on the top layer of the image volume. Then, the segmentation at the current layer is determined using a level set method. Then, the subsequent layer down is initialized with the segmentation from the current layer, and the segmentation into root and non-root regions is again found, and so on until the bottom layer. In , the authors supplemented the system with the ability to recover upward-growing roots as well. Finally, the key ideas of RooTrak were extended in  to track, and differentiate, two or more plants whose roots are growing in the same pot.
This work is most similar to Mairhofer et al.’s RooTrak [13, 14] in that both employ level set formulations. However, our method incorporates the three-dimensional growth of root systems in a different way by directly implementing a three-dimensional level set segmentation algorithm. In addition, we allow for different kinds of orientations and initializations.
2 Background and related work on level sets
We formulate the segmentation problem using a level set approach, following the general form of Rousson and Deriche 2002 [18, 19] in their formulation for scalar images (a review of level sets for segmentation is , and much of the basic notation here is consistent with that work). Let be the image domain, and in this paper consists of all the pixels forming the images of the X-ray CT volume. The segmentation problem is to find the division of into two sets and . The boundary dividing the two sets of pixels is the contour . Using the implicit representation of a contour, the set of pixels is the zero level set of the embedding function , so that . then is the signed distance function from the contour: at the contour, positive for pixels in and negative for pixels in .
In the Rousson and Deriche formulation, it is assumed that the two classes of pixels are drawn from Gaussian distributions,, . Given an initial division of pixels into and , with Gaussian parameters and , the energy function to minimize is
where is the greyvalue of a pixel and is a regularization parameter that favors smaller contour surfaces .
The gradient descent equation to evolve is
The process in many level-set methods, per iteration, given these preliminaries can be represented roughly in four steps as shown in Algorithm 1. Because of the presence of in Equation 2, we use the narrow band hypothesis, and the band distance is a parameter in our method. To compute the curvature , we use the difference of normals method [26, 11], and as in  approximate because is the signed distance function.
3 Level sets for segmenting root systems
We modified many aspects of the approach presented in [18, 19] such that segmentations of large volumes with small, thin roots would be successful. There are two major ideas, and both are implemented through the use of occupancy grids. The first idea is that the computational resources should be focused on the active portions of the contour (Section 3.2). The second idea is that the image domain is not divided into two sets and , but rather three: , , and , and by we indicate those voxels that have not been explored yet and are unlabeled. By doing so, we are able to take into account some of the physical characteristics of the root CT volumes (Section 3.3). Before we get into the details of these ideas and how they affect the steps presented in Section 2 of the Rousson and Deriche level set implementation, we describe the initialization procedure in the context of root segmentation in Section 3.1.
As mentioned in Section 2, and are initialized. In the context of root segmentation, this must be done manually. In [13, 14], the user initializes the root regions from the top layer of the image volume, corresponding to the top layer of the material in the pot. In our system, the user initializes a selection of images by marking the root regions, which will make up , with red color. Examples are shown in Figure 8. The choice of red is arbitrary, and can be changed to better match user preferences. Since these regions will serve as the initial values of the root class segmentation, it is important that only roots are marked, but it is not necessary that all roots are marked in an image. Our practice is to indicate root regions with squiggles or dots on the larger root regions, and leave regions where we are unsure about the presence or absence of a root unmarked. In particular, in the turface medium, it can be difficult for the user to determine root versus non-root regions (Figure (b)b).
The orientation of initialization images is chosen according to the ease of locating the root region for that particular plant and medium. For example, in maize the plant has a great deal of downward growth, so it may be easier to localize roots in image slices from the CT volume that are parallel to the ground plane than in other orientations, particularly for younger plants, as in Figure (b)b. On the other hand, more mature plants, such as that shown in Figure (c)c have large root regions that may be initialized from slices perpendicular to the ground plane, as with plants that grow in a variety of directions, such as cassava (Figure (a)a). Multiple images are marked per dataset, possibly from multiple orientations as appropriate, with the number of images being roughly proportional to the complexity of the root system or dataset size. Since we use the narrow band assumption, the distance from the initialization regions and the root regions to be discovered will influence the number of iterations of steps 1-4 in Section 2. For this reason, we typically space the initialization images evenly throughout the volume.
3.2 Active regions of the contour
The modifications we make to the general level set algorithm in Algorithm 1 has at its base the splitting of the contour into two parts: those regions representing the active region of the contour and those representing the static regions of the contour . This splitting operation is done for two reasons: 1) to reduce runtime, 2) prevent root regions discovered in early iterations from being smoothed away in later iterations. Others have pursued strategies that aim to set portions of the contour with the aim of reducing runtime of level set methods, such as the coherent propagation strategy [24, 25] mentioned in Section 1. In this section, we will explain our methodology.
A voxel is classified as part of one or the other contour sets in the following way. Given that the algorithm is at iteration, a function records the number of times a voxel has been a member of (the root class) from iteration to . A parameter represents the maximum value of such that voxel is part of . Consequently, we define the two sets according to this function, , and for any voxel : and . Figure 9 shows an example of the regions versus the regions.
With this definition of the active region of the contour, we will now the introduce our use of occupancy grids into the level set method. Let be the occupancy grid indicating the presence or absence of active regions of the contour . Each cube of covers voxels, where and is the band size for the narrow-band distance computation of (Line 1 from Algorithm 1). To mark per iteration, first all entries are . Given a voxel in , the grid cube location is found and marked , and its 26-connected neighbors within the grid coordinate system are found, and also marked . This is repeated for every voxel in . The result is that the grid cubes marked are guaranteed to contain the tube of voxels within voxels of , provided .
’s guarantee of containment of the -voxel wide tube of voxels around is important, because we use an exact linear time algorithm for computing the distance map only within this region. Use of such an algorithm provides an efficient narrow-band distance transform implementation. The algorithm we use for distance transform (DT) computation is that of Meijster et al. , a grid sweep algorithm whose complexity is linear in the number of elements in the grid. Another advantage of the algorithm in  is that it is easily parallelizable. Suppose the smallest Euclidean distance from voxel to the contour is . Then, the truncated Euclidean distance is , so that all voxels within the narrow band will have . We alter  to calculate the truncated Euclidean distance, and this altered version is applied within the regions of the image volume where . We notate this set of voxels as , or . As a result, the complexity of our distance transform step is linear in .
The modification of the process is summarized in Algorithm 2 for one iteration. First, all entries of are set to zero. Line 2 is as before from Algorithm 1. The active and static portions of the contour are determined in lines 3-4. Then, grid cubes of that include active regions of the contour are updated in lines 5-8. The modified distance transform of Meijster et al. is applied to those voxels in in 10. We note that is a signed distance transform; so voxels have . For those voxels within the narrow band, the gradient and new value of is computed (lines 12-14). We compute distributions by use of histograms of the image values; if the class labels change, we need to update the histogram with the change (line 16). Finally, updated distributions can be quickly computed by way of examining the histograms for each class (line 17).
3.3 Occupancy grid for gradual exploration of image volumes
Large pots are characterized not only by correspondingly large image volumes, but also variations within the pot as a result of compaction, water density variation, and root density variation. The initialization procedure outlined in Section 1, which produces a set , and following the classical procedure of dividing into two classes, may not produce a set with distribution parameters that are representative of the regions surrounding . To deal with this issue, we also use the occupancy grid introduced in Section 3.2 to incrementally add new regions to . To do so, we require a second occupancy grid , where the grid cubes where represents regions that have been explored, or the history of the regions that the method has visited in the evolution of the front.
The algorithm for this procedure is shown in Algorithm 3. Any grid entries where are checked against ; if has already been explored, or , then nothing is updated. Otherwise, is set to explored, and its contents are added to the class histogram for , non-root regions. Grid cubes where represent , or the unlabeled region.
Note that we set the same grid edge length for and , but that it is possible to alter this scheme by making one grid cube edge length a multiple of the other, for instance, and the implementation is still straightforward. One may want to pursue such a strategy to reflect some differences between the band size for the narrow band implementation of level sets and the assumptions about the distance around root regions () needed to properly model it.
3.4 Termination criterion and extraneous region removal
The termination of the method occurs when , where is a constant. In other words, the number of voxels on the active portion of the contour has become small enough that we terminate.
Following termination of the iteration level set procedure, the connected components are determined from the set of the voxels that are in . Those components that are connected to the initialization set are retained, the voxels in the remaining components are marked as , the non-root class.
X-ray CT images were acquired by a North Star Imaging (NSI) X5000 X-ray CT instrument. Five datasets were acquired with a range of characteristics, and they are listed in Table 1. In some cases, only a selection of the image slices were selected from the entire volume because the root system had not yet grown through the entire pot. In this case, the truncated volume dimensions are listed in the table.
|Dataset||Plant species||Medium||Voxel size||Dimensions||No. init. images||No. iter.s||Time|
|1||Maize||expanded clay||66 m||10||1||6||883||4.79 hr.s|
|2||Cassava||Berger 45 potting medium||55 m||10||1||11||302||25.26 min.s|
|3||Soybean||expanded clay||62 m||10||1.2||1||171||8.03 min.s|
|4||Maize||expanded clay||100 m||20||1.5||17||514||4.69 hr.s|
|5||Maize||turface||1 mm||10||1.05||19||598||3.02 hr.s|
4.1 Implementation details
The parameters described in the method are also shown in Table 1. The maximum contour membership threshold was set to for all datasets, and the termination constant for all datasets as well. The absence of material, such as outside of a pot, manifests in the images as a dark color. We specify a minimum greylevel value for members of and to prevent voxels from being added to either set if their greylevel is below the minimum. In addition, we specify a wide band of values that , the root class, may take. A voxel whose greylevel is outside of this band is not permitted to change its class label from to . Since the root and non-root colors have such a high degree of overlap, this simple implementation detail prevents erroneous front evolutions which can be easily predicted a priori. We implemented our method in C/C++, and all results in this paper were produced on a machine with a 12 core Intel Xeon(R) 2.7 GHz processor and 256 GB RAM
4.2 Results and discussion
The datasets represent a range of characteristics in terms of the texture and greylevels of the medium as well as root thickness, as defined in terms of pixels. The reconstruction results of using the proposed method are shown in Figures 13-21. The maize plant from Dataset 1 and shown in Figure 13 represents the most mature root system of the set. Qualitatively, the method performs well here, recovering both the large and small features of the root system. False positives are present mainly on the lower side of the plant, where some artifacts in imaging resulted in blurring near the edges.
The cassava root system from Dataset 2 (Figure 14) is also characterized by large and small roots, but with a different hierarchy than maize. We found that in this dataset, better results were obtained when the iterative growing adaption in Section 3.3 was not used; as a result, we set for all at the start of the method. It may be that larger regions than are needed to model , since the root regions are so large in Dataset 2 compared to other datasets. The method was able to recover very small roots which were not obvious to the human eye when viewing the dataset slice by slice, and reconstructs the large roots well.
Dataset 3 is the only example of soybean in our experiments, and that species’ root system is characterized by the traditional thin, elongated structures as well as rounded regions called nodules. The parameter was increased to in this dataset to prevent leaking from the nodules, in the class, to the interior of the expanded clay particles, which had similar greylevels. Very small roots were not recovered, and in the images the appearance of the very small roots was very faint. However, larger roots and all of the nodules were recovered. Unusually, this dataset required only one initialization image and had a reasonable run time, at eight minutes.
Dataset 4 was an example of maize planted in expanded clay. Here, the small size of the roots (defined by the number of pixels) in the middle to lower sections of the pot combined with the expanded clay medium resulted in some challenges for the proposed method. Figure 20 shows the recovered root system. As is evident, sections are recovered but the root system is much more extensive than was discovered by the method. In addition, in this dataset alteration of the parameter did not prevent the leaking phenomenon into expanded clay particles, shown in Figure (b)b near the top of the root system.
The last dataset again shows an example of maize, and this time the medium was turface, and is shown in Figure 21. In this dataset, the differentiation of medium from root can be quite difficult visually, and the root regions are very noisy. The method was able to recover parts of the root system without false positives, however the root reconstruction was truncated and smaller roots not reconstructed.
The method was able to recover larger roots on the majority of datasets, with few false positives. Small root recovery was problematic, though that may be a function of the parameters of the X-ray CT data acquisition. Figure 24 shows, as an example, the scale of the small roots in Dataset 4 and the corresponding detection by the proposed method.
The use of active regions of the contour resulted in a particular pattern concerning the number of active regions () versus the number of regions explored so far (). Figure 25 shows a graph of these two values over all iterations for Dataset 1. Consequently, the computational resources are focused on smaller and smaller regions near the end of the algorithm, as the total number of regions to be explored nears convergence. An illustration of the active regions on a reconstruction is shown in Figure 9.
4.2.1 Comparison to other approaches
We did asses the performance of another approach, RooTrak [13, 14], on two of the datasets here, Datasets 2 and 3. In RooTrak, there is a similarity measure enforced between neighboring layers. We found in both of these cases that RooTrak terminated early and hypothesize that the reason is that the size of the root regions between layers differed by too great a value in the types of situations that we consider; for Dataset 2, the maximum number of iterations was . For Dataset 3, the result was better with the maximum number of iterations at , but termination was still very early. The result is shown in Figure 26. The implementations of another method, the Root1 plugin for ImageJ  is not available at this time.
We presented a methodology for segmenting plant roots from X-ray CT images of them growing in soil. The method is intended to be used for a variety of media and plant root maturity levels and sizes.
While this was a general method, some incorporation of features of particular media might be useful to extracting better segmentations. For instance, in the expanded clay medium sometimes the front propagation produces a leak from the root to a section of an expanded clay granule because of the similar greylevels between the two; we mitigate this situation by increasing the value of the smoothness parameter . However, if the imaging of plant roots in a particular medium is necessary on a large scale for an experiment, modifying this method by considering specialized solutions to removing known non-root shapes given the plant species, X-ray CT imaging parameters, root size, and medium may be more efficient.
In addition, while the region exploration adaption in Section 3.2 reduces some issues associated with regional variation of texture and greylevels, it may be that distributions may need to be defined regionally to take into account physical realities.
-  S. Abiven, A. Hund, V. Martinsen, and G. Cornelissen. Biochar amendment increases maize root surface areas and branching: a shovelomics study in zambia. Plant and Soil, 395(1):45–55, Oct 2015.
-  M. Chowdhury, D. Jörgens, C. Wang, Ö. Smedby, and R. Moreno. Segmentation of cortical bone using fast level sets. In Proc. SPIE, volume 10133, pages 1013327–1 – 1013327–7, 2017.
-  P. Cignoni, M. Callieri, M. Corsini, M. Dellepiane, F. Ganovelli, and G. Ranzuglia. MeshLab: an Open-Source Mesh Processing Tool. In V. Scarano, R. D. Chiara, and U. Erra, editors, Eurographics Italian Chapter Conference. The Eurographics Association, 2008.
-  R. T. Clark, R. B. MacCurdy, J. K. Jung, J. E. Shaff, S. R. McCouch, D. J. Aneshansley, and L. V. Kochian. Three-dimensional root phenotyping with a novel imaging and software platform. Plant Physiology, 156(2):455–465, 2011.
-  T. Colombi, N. Kirchgessner, C. A. Le Marié, L. M. York, J. P. Lynch, and A. Hund. Next generation shovelomics: set up a tent and rest. Plant and Soil, 388(1):1–20, Mar 2015.
-  D. Cremers, M. Rousson, and R. Deriche. A review of statistical approaches to level set segmentation: Integrating color, texture, motion and shape. International Journal of Computer Vision, 72(2):195–215, Apr 2007.
-  M. Doube, M. M. Kłosowski, I. Arganda-Carreras, F. P. Cordelières, R. P. Dougherty, J. S. Jackson, B. Schmid, J. R. Hutchinson, and S. J. Shefelbine. Bonej: Free and extensible bone image analysis in imagej. Bone, 47(6):1076 – 1079, 2010.
-  A. A. Farag, H. E. A. E. Munim, J. H. Graham, and A. A. Farag. A novel approach for lung nodules segmentation in chest ct using level sets. IEEE Transactions on Image Processing, 22(12):5202–5213, Dec 2013.
-  R. J. Flavel, C. N. Guppy, S. M. R. Rabbi, and I. M. Young. An image processing and analysis tool for identifying and analysing complex plant root systems in 3d soil using non-destructive analysis: Root1. PLOS ONE, 12(5):1–18, 05 2017.
-  S. Ghadimi, H. A. Moghaddam, R. Grebe, and F. Wallois. Skull segmentation and reconstruction from newborn ct images using coupled level sets. IEEE Journal of Biomedical and Health Informatics, 20(2):563–573, March 2016.
-  A. E. Lefohn, J. M. Kniss, C. D. Hansen, and R. T. Whitaker. A streaming narrow-band algorithm: interactive computation and visualization of level sets. IEEE Transactions on Visualization and Computer Graphics, 10(4):422–433, July 2004.
-  S. Mairhofer, J. Johnson, C. J. Sturrock, M. J. Bennett, S. J. Mooney, and T. P. Pridmore. Visual tracking for the recovery of multiple interacting plant root systems from x-ray ct images. Machine Vision and Applications, 27(5):721–734, Jul 2016.
-  S. Mairhofer, S. Zappala, S. Tracy, C. Sturrock, M. J. Bennett, S. J. Mooney, and T. P. Pridmore. Recovering complete plant root system architectures from soil via x-ray -computed tomography. Plant Methods, 9(1):8, Mar 2013.
-  S. Mairhofer, S. Zappala, S. R. Tracy, C. Sturrock, M. Bennett, S. J. Mooney, and T. Pridmore. Rootrak: Automated recovery of three-dimensional plant root architecture in soil from x-ray microcomputed tomography images using visual tracking. Plant Physiology, 158(2):561–569, 2012.
-  A. Meijster, J. B. Roerdink, and W. H. Hesselink. A general algorithm for computing distance transforms in linear time. In Mathematical Morphology and its applications to image and signal processing, pages 331–340. Springer, 2002.
-  R. Metzner, A. Eggert, D. van Dusschoten, D. Pflugfelder, S. Gerth, U. Schurr, N. Uhlmann, and S. Jahnke. Direct comparison of mri and x-ray ct technologies for 3d imaging of root systems in soil: potential and challenges for root trait quantification. Plant Methods, 11(1):17, Mar 2015.
-  M. A. Piñeros, B. G. Larson, J. E. Shaff, D. J. Schneider, A. X. Falcão, L. Yuan, R. T. Clark, E. J. Craft, T. W. Davis, P.-L. Pradier, N. M. Shaw, I. Assaranurak, S. R. McCouch, C. Sturrock, M. Bennett, and L. V. Kochian. Evolving technologies for growing, imaging and analyzing 3d root system architecture of crop plants. Journal of Integrative Plant Biology, 58(3):230–241, 2016.
M. Rousson and R. Deriche.
A variational framework for active and adaptative segmentation of vector valued images.Technical Report 4515, INRIA Sophia Antipolis, July 2002.
-  M. Rousson and R. Deriche. A variational framework for active and adaptative segmentation of vector valued images. In Workshop on Motion and Video Computing, 2002. Proceedings., pages 56–61, Dec 2002.
-  R. Shahzad, S. Gao, Q. Tao, O. Dzyubachyk, and R. van der Geest. Automated Cardiovascular Segmentation in Patients with Congenital Heart Disease from 3D CMR Scans: Combining Multi-atlases and Level-Sets, pages 147–155. Springer International Publishing, Cham, 2017.
-  O. Symonova, C. N. Topp, and H. Edelsbrunner. Dynamicroots: A software platform for the reconstruction and analysis of growing plant roots. PLOS ONE, 10(6):1–15, 06 2015.
-  C. N. Topp, A. S. Iyer-Pascuzzi, J. T. Anderson, C.-R. Lee, P. R. Zurek, O. Symonova, Y. Zheng, A. Bucksch, Y. Mileyko, T. Galkovskyi, B. T. Moore, J. Harer, H. Edelsbrunner, T. Mitchell-Olds, J. S. Weitz, and P. N. Benfey. 3d phenotyping and quantitative trait locus mapping identify core regions of the rice genome controlling root architecture. Proceedings of the National Academy of Sciences, 110(18):E1695–E1704, 2013.
-  S. Trachsel, S. M. Kaeppler, K. M. Brown, and J. P. Lynch. Shovelomics: high throughput phenotyping of maize (zea mays l.) root architecture in the field. Plant and Soil, 341(1):75–87, Apr 2011.
-  C. Wang, H. Frimmel, and Ö. Smedby. Level-set based vessel segmentation accelerated with periodic monotonic speed function. In In Proc. SPIE, volume 7962, pages 7962 – 7962 – 7, 2011.
-  C. Wang, H. Frimmel, and Ö. Smedby. Fast level-set based image segmentation using coherent propagation. Medical Physics, 41(7):073501–n/a, 2014. 073501.
-  R. T. Whitaker and X. Xue. Variable-conductance, level-set curvature for image denoising. In Proceedings 2001 International Conference on Image Processing (Cat. No.01CH37205), volume 3, pages 142–145 vol.3, 2001.
-  C. Yang, W. Wu, Y. Su, and S. Zhang. Left ventricle segmentation via two-layer level sets with circular shape constraint. Magnetic Resonance Imaging, 38(Supplement C):202 – 213, 2017.
-  Y. Zheng, S. Gu, H. Edelsbrunner, C. Tomasi, and P. Benfey. Detailed reconstruction of 3d plant root shape. In 2011 International Conference on Computer Vision, pages 2026–2033, Nov 2011.