- © 2014 American Society of Plant Biologists. All Rights Reserved.

## Abstract

Root system traits are important in view of current challenges such as sustainable crop production with reduced fertilizer input or in resource-limited environments. We present a novel approach for recovering root architectural parameters based on image-analysis techniques. It is based on a graph representation of the segmented and skeletonized image of the root system, where individual roots are tracked in a fully automated way. Using a dynamic root architecture model for deciding whether a specific path in the graph is likely to represent a root helps to distinguish root overlaps from branches and favors the analysis of root development over a sequence of images. After the root tracking step, global traits such as topological characteristics as well as root architectural parameters are computed. Analysis of neutron radiographic root system images of lupine (*Lupinus albus*) grown in mesocosms filled with sandy soil results in a set of root architectural parameters. They are used to simulate the dynamic development of the root system and to compute the corresponding root length densities in the mesocosm. The graph representation of the root system provides global information about connectivity inside the graph. The underlying root growth model helps to determine which path inside the graph is most likely for a given root. This facilitates the systematic investigation of root architectural traits, in particular with respect to the parameterization of dynamic root architecture models.

Crucial factors for plant development are light quantity and quality as well as water and nutrient availability in soils. Regarding water and nutrient uptake, root architecture is the main aspect of plant productivity (Lynch, 2007; Smith and De Smet, 2012) and needs to be accurately considered when describing root processes. Currently, understanding the impact of roots and rhizosphere traits on plant resource efficiency is of highest relevance (Hinsinger et al., 2011). Development in this area will increase food security by enabling more sustainable production with reduced fertilizer input by improving cropping systems and cultivars for resource-limited environments (de Dorlodot et al., 2007).

Root architectural development includes architectural, morphological, anatomical, as well as physiological traits. For the systematic investigation of such complex biological systems, mathematical modeling is inevitable (Roose and Schnepf, 2008). Ideally, experiments and theoretical models are developed to mutually support each other. In this way, models are created that include state-of-the-art knowledge and have significant parameters. There are various root architectural models incorporating a multitude of processes (Dunbabin et al., 2013) that are originally based on Pagès et al. (1989) and Diggle (1988). Generally, the parameterization of such models is difficult and demands elaborate experimental effort. In this work, we present a novel approach for recovering root system parameters based on image-analysis techniques. In this way, we simplify the systematic investigation of root architectural traits, in particular with respect to the parameterization of root system models.

Imaging techniques for the visualization of soil-grown root systems in two and three dimensions include x-ray computed tomography (Heeraman et al., 1997; Tracy et al., 2010; Mooney et al., 2012), neutron radiography (NR; Oswald et al., 2008), and magnetic resonance imaging (Pohlmeier et al., 2008). NR is one of the most suitable techniques to investigate roots grown in soil, because it allows a high throughput, provides a strong contrast between roots and soil, and therefore requires little effort for image processing. A major advantage of NR as well as magnetic resonance imaging is the possibility to monitor water distribution and roots simultaneously (Menon et al., 2007; Oswald et al., 2008; Moradi et al., 2009; Carminati et al., 2010; Stingaciu et al., 2013). This is especially useful as water is a crucial factor ruling root allocation in soil (Hodge, 2010).

Images of root architecture contain a huge amount of information, and image analysis helps to recover parameters describing certain root architectural and morphological traits. The majority of imaging systems for root systems are designed for two-dimensional images, such as RootReader2D (Clark et al., 2013), GiA Roots (Galkovskyi et al., 2012), SmartRoot (Lobet et al., 2011), EZ-Rhizo (Armengaud et al., 2009), and Growscreen (Nagel et al., 2012). See also Le Bot et al. (2010) for a review of available software. A starting point for image analysis is commonly a grayscale image of a root system. The first step is to create a binary image by segmentation. Further steps include skeletonization, root tracking, and data analysis. The most common segmentation method is some form of thresholding (e.g. RootReader2D, GiA Roots, SmartRoot, EZ-Rhizo; Stingaciu et al. [2013]). Other methods include the livewire algorithm (Basu and Pal, 2012) or the levelset method (RootTrak; Mairhofer et al., 2012) that determine the borders of each root. The creation of a root system skeleton is either done manually (e.g. DART; Le Bot et al., 2010) or based on morphological operators such as thinning and closing (GiA Roots, RootReader2D, EZ-Rhizo), sometimes with options for the user to correct skeleton points (EZ-Rhizo). The root tracking step can be performed manually (DART) or based on creating a graph representation of the root system combined with Dijkstra’s algorithm, a search algorithm that finds the shortest path between two nodes inside a graph (RootReader2D; Stingaciu et al. [2013]). Furthermore, algorithms can operate on the skeleton (EZ-Rhizo) or directly on the image source (SmartRoot). In SmartRoot, the user selects a root in the original image with a mouse click and then a skeletonization algorithm determines the skeleton of the selected root. The output of all root tracking algorithms is a data structure of a set of roots that stores information such as connectivity between roots and their position in space.

Global traits of the root system are obtained directly from the segmented image or the skeleton (GiA Roots, RooTrak). Global traits include convex hull, network depth, network length distribution, maximum number of roots, maximum width of the root system, network length, and specific root length. Furthermore, the data structure from a root tracking procedure is used to obtain individual, local root parameters (DART, RootReader2D, EZ-Rhizo, SmartRoot). The latter are able to obtain root architectural parameters that can be used for model parameterization. An additional aspect is the dealing with dynamic data (i.e. images of the same root system taken at several times). Analysis of such sequences may lead to better insight into the development of the root system (e.g. DART, SmartRoot) or even reveal growth zones and their local growth velocities (Basu and Pal, 2012).

Analysis software for two-dimensional images of soil-grown root systems currently work in a semiautomated way with respect to tracing individual roots. This requires considerable user input for larger root systems. We present a new, fully automated approach for recovering root architectural parameters from two-dimensional images of root systems. The software Root System Analyzer is, to our knowledge, the first algorithm for two-dimensional analysis of soil-grown root systems that features fully automated root tracking. Only primary roots have to be initiated manually by the user. The user is also free to initiate any laterals, but this is not mandatory. Further growth of primary roots and laterals is then tracked in a fully automated way. In addition, there is a user interface that allows for manual correction of individual roots if required. In this work, we do not go into the details about the segmentation step but focus on the root tracking step and the parameterization of a root system model (Leitner et al., 2010). The described algorithm starts with a sequence of segmented two-dimensional images showing the dynamic development of a root system. For each image, morphological operators are used for skeletonization. Based on this, a graph representation of the root system is created. A dynamic root architecture model helps to determine which edges of the graph belong to an individual root. The algorithm elongates each root at the root tip and simulates growth confined within the already existing graph representation. The increment of root elongation is calculated assuming constant growth. For each root, the algorithm finds all possible paths and elongates the root in the direction of the optimal path. In this way, each edge of the graph is assigned to one or more coherent roots. The algorithm considers the fact that new branches can only emerge after the apical zone has developed, which helps in the decision of whether the root is branching or two roots are crossing or overlapping. Image sequences of root systems are handled in such a way that the previous image is used as a starting point for the current image. This is helpful in the analysis of complex root systems as well as for retrieving dynamic parameters such as elongation rates. The algorithm is implemented in a set of Matlab m-files, which makes the code flexible so that it can easily be adjusted to specific experimental setups or mathematical models.

We exemplify the approach with two-dimensional neutron radiography images of lupine (*Lupinus albus*) root systems grown in mesocosms filled with a sandy soil. Furthermore, we compare our approach with the approaches of SmartRoot and RootReader2D and demonstrate how our approach can be used to analyze large root systems.

## RESULTS

### Root Tracking in Neutron Radiography Image Sequences

Figure 1, A to D, shows the segmented images of four lupine root systems that were grown in mesocosms in a sandy soil under the same homogenous soil moisture conditions. The root systems were imaged by neutron radiography; the different colors indicate three different measurement times. The images of each sequence were registered using the plugin stackreg of ImageJ (http://rsb.info.nih.gov/ij/disclaimer.html). The segmentation algorithm was based on a matched filter response (Hoover et al., 2000). The new root tracking algorithm of Root System Analyzer was applied to each of these images. Figure 1, E to H, illustrates the sequential root tracking for the three measurements of the first image (Fig. 1A). In order to detect coherent roots, the algorithm is based on two assumptions. First, roots are elongated by growing root tips. Second, new root tips emerge in the branching zone and form new lateral roots. In this way, the root tracking algorithm uses a dynamic root architectural model to decide whether a detected root is valid from a root developmental point of view.

The first step in the detection procedure is the skeletonization of the segmented image by morphological operations. This process reduces each root to its center line. From the skeleton, a graph is created where the nodes *N* are the branching points, crossing points, and root tips at each measurement time. On the skeleton, neighboring nodes are connected by an edge stored in an adjacency matrix *A*. For each edge, the corresponding coordinates are stored in a list *E*. Figure 1E shows the skeleton of the above image together with the nodes of the graph. The same graph representation of the root system is also used by RootReader2D. The method of Stingaciu et al. (2013) is also based on a graph representation; however, in that case, the nodes represent voxels of the three-dimensional image.

In the second step, we apply the root tracking algorithm to the first measurement (indicated by the red color in Fig. 1, A–D). The algorithm can initiate a tap root system automatically by selecting the node with the largest *z*-coordinate. Alternatively, the user can initiate one or more roots manually, then all other roots are traced in a fully automated way. The edges of the graph are assigned to roots in the following way. For each existing root tip, further root growth is calculated in small time steps according to the underlying root architectural model. In each time step, all possible growth paths in the graph of a certain length are evaluated. The optimal growth path is chosen dependent on straightness and average diameter. Furthermore, it is penalized if the edge is already assigned to a root. The increment of growth for each root tip occurs along the optimal path, and corresponding edges are assigned to belong to the root. New root tips are created in all other possible directions. These tips start to grow after a time delay and form new lateral roots. The method enables the distinction between crossings and branches. Analyzing all possible paths increases the scale on which decisions are made and, therefore, makes it more likely to find the correct solution. Additionally, we developed a graphical user interface that enables a visual check and manual correction based on Dijkstra’s algorithm if needed. Similar to that, RootReader2D and Stingaciu et al. (2013) apply Dijkstra’s algorithm for the root tracking procedure.

The second step is repeated for all measurements using the edge assignment of the previous measurement as initialization. Figure 1, F to H, represents the assignment of the edges to roots for each measurement time. The color red denotes the tap root, blue denotes first order laterals, green denotes second order laterals, and magenta denotes higher order laterals. In this way, our algorithm can handle images with temporal information, like DART or SmartRoot. This benefits, on the one hand, the root tracking, and on the other hand, it enables the extraction of dynamic root traits like elongation rates of different root types.

The way the roots are tracked is new in that the underlying dynamic root architecture model is used for deciding which paths in the graph potentially belong to an individual root. The algorithm starts at the tips of user-provided or previously detected roots and simulates further root growth according to the root architectural model with assumed parameters but confined within the already existing graph representation of the root system. The increment of root elongation is calculated for a small time step. Then the algorithm finds all possible paths of this length that the root could take in the graph. “Possible” considers the fact that a new branch can only emerge once the length of the parent root is at least equal to the sum of the basal and apical zones. This helps in the decision of whether there is a branching or a crossing. The dynamic root assignment is illustrated in Supplemental Movie S1.

As in Le Bot et al. (2010), the output information for each root is an identification number, the branching order, the time of emergence, the parent identification number, the distance between the branching point to the parent root base, and the root length at the observation time. In addition, we store the area of the root in the image as well as the nodes of the graph that belong to each individual root. Together with the adjacency matrix *A* and the list of coordinates *E*, this gives us all information about the position of the root in the source images. Additional examples, as well as software and documentation, can be found at http://www.csc.univie.ac.at/rootbox/rsa.html or in Supplemental Data S1.

Certain global parameters can be determined from this result without further analysis. Table I shows means and sd of number of roots and total length over the four root systems for the three measurement times. Further postprocessing is necessary to retrieve parameters for root architecture models.

### Parameterizing the Root Architectural Model

For each of the four lupine root systems, we derived a data structure as described above (following Le Bot et al. [2010]). Since the growth conditions were the same, we view the four images as replicates. Thus, we merge the data structures and use these data to parameterize the dynamic root architectural model of Leitner et al. (2010). This model needs 12 input parameters for each root type, as shown in Table II. Figure 2 outlines the work flow for parameterizing the root architectural model.

We assume that the root system is not strictly organized in terms of root orders but in terms of different root types that do not necessarily coincide with root orders (following Pagès et al. [2004]). For lupine, we assume that there is one tap root and different types of lateral roots that grow from any predecessor with a certain probability. To distinguish these groups of laterals, we perform a cluster analysis. We use the Matlab-implemented algorithm k-means. The algorithm requires previous knowledge of the number of clusters. From visual comparison, we decided that we need two clusters of laterals, long and short. The algorithm assigns each root to one of the clusters using the observed root length. Outliers and random starting points could result in wrong or unintuitive clustering of roots. We validated the clustering by examining the clusters in the histogram of root lengths.

In a further postprocessing step, root architectural parameters are retrieved from the data structure and averaged over each root type. The estimated parameters are shown in Table II.

The values for *l*_{a}, *l*_{b}, and *l*_{n} can directly be retrieved from the data structure. The great variation of the parameters *l*_{a}, *l*_{b}, and *l*_{n} is reflected in large sd values. This is not a measurement error but is the true variability that can be observed in the original images. The root radius *a* is computed from the stored values of area and length and averaged over each root type. The branching angle Θ is measured from the corresponding coordinates of a root and its predecessor.

Root elongation is described in the model by the root growth function λ, which is dependent on the maximal root length *k* and the initial elongation rate *r*. These two parameters cannot be observed directly but are obtained by fitting the root growth function λ to the data of root age versus root length for each type (Fig. 3). Root age is only known exactly for the zero order roots. For higher root orders, root age is estimated based on the growth rate and length of the apical zone of the predecessor. Note that errors in the calculation amplify in higher orders. However, the sd values of *r* and *k*, which are given by the diagonal elements of the covariance matrix, are small. The maximal number of branches per root *nob* is computed from the values of *k*, *l*_{n}, *l*_{a}, and *l*_{b}:Variation from growth direction is described by two parameters, σ and *N*. The deflection parameter σ describes the expected angular deviation from growth direction per cm of root length. It is computed along the individual roots and averaged per root type. The parameter *N* describes the strength that keeps the root in its initial growth direction. For the tap root, this describes gravitropism, and for lateral roots, this describes exotropism. This is the only parameter that we obtained by visual comparison. Different tropisms are hard to disentangle without specific experiments (e.g. hydrotropism and gravitropism; Tsutsumi et al., 2004).

### Modeling Case Study

Models of root architecture and function increase our mechanistic understanding of soil-plant interactions. For example, Schnepf et al. (2012) investigated phosphorus uptake by a root system as affected by root growth and exudation. The parameterization of the root architectural model is of utmost importance for a realistic description of root surface development in soil. In this work, we use the parameter set of young lupine plants to create simulated root architectures based on the model of Leitner et al. (2010) and to investigate the dynamic development of global root system properties.

We performed a topological characterization of 100 simulated root systems according to Fitter (1987). The analysis reveals that, after 25 d, the mean altitude is 73.60 ± 3.28, the mean magnitude is 240.09 ± 36.95, and the mean external path length is 8,376.65 ± 1,667.71. Table III shows the dynamic development of theses characteristic values. According to Fitter et al. (1991), the topological index is a number between zero and one and correlates positively with the exploitation efficiency of the root system. Thus, these numbers provide an estimate of the exploitation efficiency and enable quick comparisons (e.g. between cultivars).

Simulations were performed in three dimensions and constrained within the mesocosm as in the experimental setup. This is illustrated in Figure 4A. Figure 4B shows further analysis in terms of root length densities. Such information can be further used in plant-soil interaction models by calculating, for example, sink terms for plant water or nutrient uptake. The Matlab implementation of the root architecture model facilitates both postprocessing and coupling to other models. The root growth inside the mesocosm is visualized in Supplemental Movie S2.

The dynamic development of total root length and number of roots is illustrated in Figure 5. The solid and dashed lines represent the total root length and number of roots of a simulated lupine root system based on the parameter set of Table II. The asterisks with error bars show the corresponding measured root lengths and number of roots according to Table I. We can observe that, initially, the average root length is approximately 1 cm; this increases over time to an average length of approximately 2 cm.

### Comparison with SmartRoot and RootReader2D

Several image-analysis tools for two-dimensional root systems are available. Our algorithm combines features that are partly present in other softwares. It shares the skeletonization and graph representation of the root system with RootReader2D, while it is similar to SmartRoot with respect to the use of image sequences and also with respect to postprocessing and parameter retrieval. We chose the first lupine image shown in Figure 1A and analyzed it with both softwares for comparison.

In Table IV, we compare the results of Root System Analyzer with the results obtained by analyzing the same image with SmartRoot. The resulting tracings according to both tools are shown in Figure 6.

The main parameters, such as overall root length, internodal distance on zero order roots, and branching angles, agree well. Thus, we suggest that both produce satisfactory results. However, there are two main differences in the results. First, our automated algorithm detects more very small roots than SmartRoot (i.e. roots are detected at an earlier stage). This means that the overall root length is still similar; however, the number of roots (in this case, second order) is larger, and thus the internodal distance (in this case, on the first order roots) is smaller. This has implications for the parameterization of root architectural models. Second, the image was analyzed independently with both imaging tools, and this resulted in the fact that the zero order root took a different path in the two cases. The underlying algorithms are very different between the two tools. Our automated procedure is based on skeletonization and graph representation of the root system coupled with a dynamic root architectural model for automatic root tracking. It does feature manual correction should the user find it necessary. SmartRoot works on the (grayscale) image itself; the user can follow individual roots by mouse clicks. SmartRoot does feature an automated algorithm for individual roots that is triggered by a mouse click inside the root, and this algorithm is based on finding the midline of the root and progressing forward and backward until the base and tip of the root are reached. The user input required by our automated tool is considerably less; thus, it is more suitable for the complete analysis of large root systems. Although the underlying algorithms for root detection are different, the handling in the graphical user interface is similar, in particular with respect to handling sequences of images of the same root system. Both use the tracings from the previous time step as starting points for current root tracking. In SmartRoot, the user is required to determine anchor points in each image of a sequence in order to be able to compensate for little shifts between images. The postprocessing in SmartRoot is user friendly, with links to SQL databases or output produced as txt files and images. When using Root System Analyzer, one needs to be familiar with Matlab. However, this offers a lot of postprocessing options within Matlab, and several Matlab functions for parameterization of root architectural models are already provided for.

Like Root System Analyzer, RootReader2D is based on an algorithm that first creates the center line of each root and then a graph representation where the root branching points, crossing points, and root tips are the nodes of the graph and the edges are the links between the nodes. Root System Analyzer and RootReader2D produce the same skeleton and graph with only tiny variations of the position of the center line due to differences in the skeletonization procedure. Root tracking in RootReader2D is based on Dijkstra’s algorithm. However, the handling for the user is similar to that of SmartRoot in that each root is tracked individually by one or more mouse clicks. Neither RootReader2D nor SmartRoot offers a fully automated procedure for tracking all the roots of the root system, like Root System Analyzer. Thus, this tool is more suitable for the extensive analysis of larger root systems. Results of analyzing a large root system are shown in the next section.

### Application to a Large Maize Root System

Root System Analyzer is so far the only algorithm based on imaging techniques that was applied to root systems as complex as that of a 78-d-old maize (*Zea mays*) root system. There is the potential to use this algorithm to extract more information from the large number of existing two-dimensional images of excavated root systems such as those from the “Wurzelatlas” series.

To test the algorithm on a large fibrous root system, we used the software to track roots in a two-dimensional drawing of a large maize root system (Kutschera et al., 2009). Assuming the sowing time to be May 1, the root system is 78 d old. It goes down to a depth of approximately 60 cm and has a maximum width of approximately 120 cm, thus being much larger than the laboratory root systems typically used for imaging. The original hand drawing was scanned at high resolution (1,200 dots per inch), yielding an image size of 116 megapixels. The image was turned into a black-and-white image by thresholding and then analyzed by Root System Analyzer.

In the tracking procedure, we encountered two problems. First, it is hard to distinguish individual roots in the highly dense part of the root system just below the stem. Thus, the node points assigned in this area are questionable. Second, the algorithm does not yet detect multiple primary roots automatically. Therefore, we assigned several primary roots by hand and only then started the automated procedure. In this way, we picked 21 initial roots and the algorithm detected 3,457 roots automatically.

The tracking results agreed well with data from the literature. Roots of cereal plants, such as maize, typically have three or four orders (Roose et al., 2001), and this has been confirmed by our analysis. The resulting root system has a maximum root order of four; however, most of the roots belong to orders zero to three. The original image and the tracked root system are shown in Figure 7. The number and length of roots in each order that were detected in this two-dimensional view of the root system are provided in Table V. Selected parameters of orders zero to three are presented in Table VI; dynamic parameters *r* and *k* could not be estimated from only one image. Average root architectural parameters for the first three orders compare well with those presented by Pagès et al. (1989) with the exception of the length of the apical zone of zero order roots.

## CONCLUSION

We present a novel approach for root tracking from two-dimensional images of root systems. The combination of the graph representation of the root system together with an underlying root architectural model results in a reliable method for root assignment. In this way, the decisions are made based on global information about the connectivity within the graph, and it is assured that root assignment is valid from a root system development point of view. Particularly, this is beneficial for dealing with overlaps and intersections in two-dimensional images. Comparison with SmartRoot shows similar capabilities for postprocessing, parameter retrieval, and handling of image sequences but differences in the underlying algorithm for root tracing. Our algorithm uses the global information contained in the graph representation of the root system and traces all roots of the root system in a fully automated way. Comparison with RootReader2D shows similarities with respect to the automated skeletonization and graph creation procedure. However, also in RootReader2D, individual roots have to be traced manually by at least one mouse click. Thus, our new algorithm is more suitable for larger root systems, as demonstrated by tracing the roots of a large maize root system. This is promising for retrieving more information from the many two-dimensional hand drawings of excavated root systems (e.g. in the Wurzelatlas book series).

If an image sequence is available that shows root system development, the root assignment becomes more robust, since the detected roots of the previous image act as initial roots for the current image. Furthermore, the dynamic development of the root system can be analyzed, like elongation rates for individual soil-grown roots.

Our approach works on the segmented images of root systems. In contrast, Lobet et al. (2011) and Naeem et al. (2011) work directly on the original image data. These data contain more information and, therefore, can better seize small-scale image features and thus work with high precision on small root systems. However, they lack the ability to use global information about connectivity, which is the strength of our graph-based technique. A combination of these approaches could further increase the quality of root tracking algorithms.

While the literature for automatic root assignment is limited, there is a vast literature for biomedical applications, like vessel tracking (Lesage et al., 2009). Ideas from this area could further benefit image analysis for root architecture.

Plant root experiments are extremely costly regarding labor and time. Therefore, we suggest that the proposed image-analysis method can help to extract the most information from root system images. In particular, we propose this approach for image-based parameterization of root architectural models.

## MATERIALS AND METHODS

The starting point for the automated root tracking algorithm is a segmented two-dimensional image. In order to annotate coherent roots from this image, we make use of morphological and graph theoretic methods to perform a fully automated or semiautomated root tracking. The methods are implemented in Matlab. The corresponding m-files are freely available at http://www.csc.univie.ac.at/rootbox/rsa.html and as well in Supplemental Data S1.

### Binary Two-Dimensional Neutron Radiography Images of Root Systems of Lupine

Neutron radiography is one of the few noninvasive and nondestructive techniques available for studying plant root systems in situ (Moradi et al., 2009). Although this technique is too time consuming for high-throughput screening of plants, it is capable of imaging plant root systems in opaque soil. Providing a soil environment to the plant roots is necessary for studying plastic traits of root systems such as root length. Furthermore, this technique can be used to recover complete root systems.

We demonstrated the applicability of Root System Analyzer to detect root systems in neutron radiography images. Segmented neutron radiography images were provided by B.F., who performed an experiment on the effects of phosphorus and water inhomogeneity on root system growth of lupine (*Lupinus albus*; B. Felderer and P. Vontobel, unpublished data). Briefly, the experimental setup, filtering, and segmentation are outlined below.

Single white lupine plantlets were grown in Teflon-coated aluminum containers of 27 × 27 × 1.4 cm over a period of 26 d. The aluminum containers were filled with sandy soil extracted from the forefield of an open-cast mine near Cottbus (Welzow Süd, Germany). The soil contained almost no organic matter. Plants were grown in a climate chamber at a humidity of 60% with a 16-h/8-h day/night cycle with 21°C/16°C temperature. The root system was visualized with neutron radiography at 12, 19, and 26 d after germination (at the Paul Scherrer Institute in Villigen, Switzerland), which is highly suitable to visualize roots and water in soil (Menon et al., 2007; Moradi et al., 2009). The field of view was 18 × 18 cm with a nominal resolution of 0.176 mm; thus, four scans per plant container were needed to record the whole sample. Image sequences were registered using the ImageJ plugin stackreg.

The neutron radiography images were first corrected for their flat field, for sample scattering, and for variation of beam intensity between the scans using the Quantitative Neutron Imaging Software (Hassanein et al., 2005). After stitching, a matched filter response was applied for root segmentation using a software developed by Anders Kästner (Paul Scherrer Institut), which amplifies thin structures and creates one binary image per measurement.

For this study, the images of four containers with initially homogenous water content were selected (Fig. 1, A–D) in order to demonstrate the new root tracking approach.

### Automatic Root Tracking

#### From the Binary Image to a Graph

In a first step, we derive an approximated center line from each binary image *im*_{B}. Therefore, we create a skeleton *im*_{S} from *im*_{B} by iteratively applying two morphological operations: thinning, which removes boundary pixels without disconnecting any domains (Lam et al., 1992); and closing (i.e. dilation followed by erosion).

In the next step, *im*_{S} is used to create a graph representation of the root system. For each pixel on the skeleton in *im*_{S}, we count the number of neighboring skeleton pixels in an eight-pixel neighborhood. This yields a matrix *im*_{C}. If a pixel of *im*_{C} has exactly one neighbor, it represents a leaf in the corresponding graph. If it has more than two neighbors, it represents a branching or crossing point. Therefore, all pixels that correspond to nodes in the graph are given byThe connected components in *im*_{N} represent the nodes in the graph. They are labeled with the Matlab function bwlabel, and the coordinates of the nodes *N* are exactly the centers of each connected component. Two nodes are connected by an edge if they are connected by the skeleton exclusive of the nodes (i.e. *im*_{S}\*im*_{N}). For each edge, the corresponding pixel coordinates from the skeleton *im*_{S} are stored in a list *E*. Node connectivity is represented by an adjacency matrix *A* of the undirected graph, where the entries represent the edge indices of *E*. In a final step, results are improved by applying a one-dimensional Gauss filter to the edge coordinates and, optionally, by removing small terminal edges from the graph. The above approach is implemented in the file image2graph.m. For the root tracking algorithm, we use the following edge weights. First, the length of each edge is calculated from the coordinates *E* using the Euclidean distance. Second, the approximate radius of the root at each coordinate is obtained by calculating the distance function of the binary image *u*_{b}. In the resulting distance matrix *D*, each pixel value holds the distance to the closest boundary of the root. Therefore, the values of *D* at the center line are exactly the root radii. The distance matrix *D* is calculated by the Matlab function bwdist (described by Maurer et al., 2003).

#### Root Tracking in the Graph

In the root tracking algorithm, we assign each edge of the graph to one or more individual roots and retrieve information about the connectivity of these individual roots. We use basic knowledge about root system development to detect coherent roots. Roots emerge from growing root tips. Each root consists of a basal zone succeeded by a branching zone followed by an apical zone. Lateral roots emerge only in the branching zone, after the apical zone has developed. The following algorithm is implemented in the Matlab file trackRoots.m.

We store a list *T* of growing root tips. Initially, these are set manually, or if there is only one tap root, the algorithm picks the node with the largest *z*-coordinate.

The time that has passed at each measurement time is subdivided into smaller time steps typically used for root growth modeling. While there are any tips in the list *T*, the growth increment for each root during each time step is calculated according to the underlying root growth model. For each root, we add edges until the growth increment is reached. This is done in the following way. (1) Find all possible growth paths in the graph from the node that represents the growing root tip (getPath.m). This is performed in a recursive way. All paths from the node with a length smaller than a maximal path length *R*_{S} are retrieved. Each path is a list of edges with corresponding nodes and edge weights such as root length and radius. (2) Choose the optimal path by evaluating the quality *q*_{j} of each path with index *j* as described in the following section. The optimal path has index *j* = *argmax*_{i∈I}*q*_{i} with *q*_{j} = max_{i∈I}*q*_{i}. The heuristic is based on the path length, radius, straightness, and information on previous edge assignments. (3) If the optimal path *j* fulfills the quality criterion *q*_{j} > 1, the first edge of the optimal path is added to the growing root, and the root tip moves to the next node. The assigned edge is denoted as visited in order to penalize multiple visits. Otherwise, the root stops growing and is removed from the list *T*. (4) New laterals emerge into all other possible directions. These new laterals are added to the list *T* and start to grow after a time delay in order to let the apical zone of the base root develop. The procedure (steps 1–4) is repeated for all time steps until no growing root tips are left.

The algorithm dynamically assigns each edge in the graph to one or more roots. Following Le Bot et al. (2010), the output information for each root is an identification number *number*, the branching order *order*, the time of emergence *ct*, the parent identification number *predecessor*, the distance between branching point to the parent root base *prelength*, and the root length at the observation time *length*. In addition, we store the area of the root in the image *area*, as well as the nodes that belong to each individual root *path*.

#### Path Evaluation

A path in the graph is evaluated by the following heuristic quality criterion (implemented in the Matlab file quality.m). The quality criterion is based on the average root radius and the straightness of the root.

The algorithm starts at a growing root, which is typically thicker than its lateral branches. Therefore, we want to follow the path with the largest average root radius *r* to find a coherent root. The computation of the average path radius *r* is described above. Additionally, we take into account if the edge is already assigned to another root. In this case, we subtract its average root radius from the average radius of the edge.

We assume that roots grow preferably straight. Thus, we want to favor the straightest path in the graph. The straightness *s* between two coordinates *x*_{1} and *x*_{2} is given by the ratio of the length of the straight line between the points and the length along the edges. The coordinate *x*_{1} is located at a distance *R*_{a} along the root in front of the current root tip location, and *x*_{2} is located at *R*_{a} along the new path after the current root tip. Therefore, the straightness is given byWe define the quality *q* of each path as the product of the radius *r* and straightness *s*:The exponent *e* describes the strength of influence of the straightness (according to experience, *e* = 2 performs well).

#### Image Sequences

Image sequences enable the calculation of elongation rates of individual roots. Furthermore, our algorithm uses this dynamic information to improve the root tracking. Young root systems are still small and easier to track automatically. After the tracking of the initial measurement, the root assignment of the previous measurement always acts as an initial root system for the current image. If the time difference between the measurements is small, the root development is easily manageable automatically by the tracking algorithm. This procedure facilitates correct root tracking even in large and complex root systems.

Manual correction is hardly necessary. However, we provide a graphical user interface (RSA_GUI.m) to enable manual root corrections if required. This includes the removal of wrongly assigned roots as well as manual assignment by picking the starting and end points of the root. Figure 8 shows a screen shot of the graphical user interface after automatic detection of the roots of the lupine root system at the first measurement time. If a path is found to be incorrect, the user can remove the root by selecting “remove” and clicking at any point on this root. A new root can be added by selecting “add root” and clicking on start and end points of the root. The path of minimal length between those nodes is found by Dijkstra’s algorithm. Dijkstra’s algorithm is a graph search algorithm that solves the single-source shortest path problem for a graph with positive edge path costs. In our application, these costs are defined as edge lengths or edge radii. Thus, for any wrongly determined root, one mouse click for its removal and two mouse click for creation of a new path are required.

#### Processing Efficiency

On a standard personal computer, the analysis of young root systems takes only a few minutes. The most time-consuming step is after the automated analysis if the user wishes to manually correct some of the assigned roots. The automatic detection of roots in a large root system such as a large maize root system (11,658 nodes, 14,783 edges, and 3,478 detected roots) may take a few hours on a standard personal computer. In mature root systems, there are limitations with respect to the complexity of the root system. The part just below the stem is so densely rooted that individual roots can no longer be recognized. As a result, this area appears as a white or black area in the segmented image, and detection fails. This can be overcome if sequences of the same root system at earlier times are available.

Note that Root System Analyzer can work with all two-dimensional imaging techniques that offer a sufficient resolution such that a segmented binary image can be created, where roots are represented by white pixels and the background by black pixels. Furthermore, all pixels that belong to the roots must be connected. Currently, there is no automatic segmentation provided in the software. Generally, this kind of segmentation is a nontrivial task for most imaging techniques due to different types of artifacts. Depending on the imaging technique that is used, there might be specialized software available for the segmentation step.

### Model Parameterization

We use the resulting data structure to parameterize the dynamic root architectural model of Leitner et al. (2010); however, any other root architectural model could be parameterized as well. This model needs 12 input parameters, as shown in Table II. They are computed as averages over root order or root type (Pagès et al., 2004). In order to find root types, we perform a cluster analysis regarding root length using the Matlab function kmeans. Alternatively, root types could be distinguished according to the developmental stage of a plant (e.g. a maize plant). This is not implemented in the current version of Root System Analyzer, but any user is free to replace the kmeans function by any other function that groups the roots into different types. The computation of the root architectural parameters for each root type is implemented in the Matlab file analyseRS.m.

The values for *l*_{b}, *l*_{a}, and *l*_{n} can be retrieved from the data structure as described below

For each root *j* with laterals *I*, the length of the basal zone is given by *l*_{b,j} = min_{i∈I}(*prelength*_{i}). The apical zone is given by *l*_{a,j} = *length*_{j} − max_{i∈I}(*prelength*_{i}). Apical and basal zones can only be determined if the root has at least one branch. The internodal distances are given by *l*_{n,j,i} = *prelength*_{i+1} − *prelength*_{i}, where the distances between branching points and the parent root base are in ascending order. Internodal distances can only be determined if the root has at least two branches.

The root radius *a* is computed from the stored values of area and length and averaged over each root type. The branching angle Θ is calculated as the angle between two straight lines approximating the root and its lateral in the branching point.

Root elongation is described in the model by the growth function λ, which is dependent on the maximal root length *k* and the initial elongation rate *r*:These two parameters are obtained by fitting the root growth function λ to the data of root age versus root length for each type. The sd values of *r* and *k* are the diagonal elements of the covariance matrix.

Variation from growth direction is described by two parameters, σ and *N*. The deflection parameter σ describes the expected angular deviation from the growth direction per cm of root length. This value is obtained by calculating the angular change along the roots. All parameters except for *N* are computed along the individual roots and then averaged per root type. The parameter *N* describes the strength with which the root will keep its initial growth direction (i.e. gravitropism for tap roots, exotropism for later roots). This is the only parameter that we obtained by visual comparison.

### Supplemental Data

The following materials are available in the online version of this article.

**Supplemental Movie S1.**Illustration of the dynamic assignment of roots by the tracking algorithm.**Supplemental Movie S2.**Visualization of the simulated root growth.**Supplemental Data S1.**Software, documentation, and additional examples.

## Footnotes

The author responsible for distribution of materials integral to the findings presented in this article in accordance with the policy described in the Instructions for Authors (www.plantphysiol.org) is: Daniel Leitner (daniel.leitner{at}univie.ac.at).

↵1 This work was supported by the Austrian Science Fund (FWF) (grant no. V220–N13) and by an APART fellowship of the Austrian Academy of Sciences at the Computational Science Center, University of Vienna (to D.L.).

↵[C] Some figures in this article are displayed in color online but in black and white in the print edition.

↵[OPEN] Articles can be viewed online without a subscription.

↵[W] The online version of this article contains Web-only data.

### Glossary

- NR
- neutron radiography

- Received September 3, 2013.
- Accepted November 6, 2013.
- Published November 11, 2013.