Recovering root system traits using image analysis Exemplified by 2-dimensional neutron radiography images of lupine

One-sentence summary: Image-based parametrisation of root architectural models is advanced by a new approach for the analysis of image sequences of plant root systems, exemplified by neutron radiographic images of root systems of soil-grown lupine plants. 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 skeletonised 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 favours 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 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 decide 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 parametrisation of dynamic root architecture models.


Introduction
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 a 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 modelling is inevitable (Roose and Schnepf, 2008).
Ideally, experiments and theoretical models are developed mutually supporting each other. In this way models are created which 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) which are originally based on  and Diggle (1988). Generally, the parametrisation 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 parametrisation of root system models.
Imaging techniques for the visualisation of soil-grown root systems in two and three dimensions include x-ray computed tomography (Mooney et al., 2012;Tracy et al., 2010;Heeraman et al., 1997), 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 (Oswald et al., 2008;Moradi et al., 2009;Carminati et al., 2010;Menon et al., 2007;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 comprehend 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 is designed for 2-dimensional images, e.g. RootReader2D (Clark et al., 2013), GiA Roots (Galkovskyi et al., 2012), SmartRoot (Lobet et al., 2011), EZ-Rhizo (Armengaud et al., 2009), Growscreen (Nagel et al., 2012). See also Le Bot et al. (2010) for a review of available software. Starting point for image analysis is commonly a grey scale image of a root system. The first step is to create a binary image by segmentation. Further steps include skeletonisation, root tracking, and data analysis. The most common segmentation method is some form of thresholding, e.g.
RootReader2D, GiA Roots, SmartRoot, EZ-Rhizo, or 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 boarders 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 (GiARoots, 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 and 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 skeletonisation algorithm determines the skeleton of the selected root. 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 root system, network length, or 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 parametrisation. 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 on the development of the root system e.g. DART, SmartRoot, or even reveal growth zones and their local growth velocities (Basu and Pal, 2012  Analysing 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 which enables a visual check and manual correction based on Dijkstra's algorithm if needed. Similar to that, RootReader2D and Stingaciu et al. (2013) apply the Dijkstra's algorithm for the root tracking procedure. This helps in the decision whether there is a branching or a crossing.
The dynamic root assignment is illustrated in the video provided in supplementary material S1. Certain global parameters can be determined from this result without further analysis. Table Table 1 shows mean and standard deviations 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.

Parametrising the root architectural model
For each of the four lupine root systems, we derived a data structure as described above ( We assume that the root system is not strictly organised 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 beforehand knowledge of the amount 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 Table 2.
The values for l b , l a , 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 standard deviations. 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 can not be observed directly but are obtained by fitting the root growth function λ to the data of root age versus root length for each type (see Fig. Error Reference source not found.). Root age is only known exactly for the 0 th roots. For higher root orders, root age is estimated based on the growth rate and length of apical zone of the predecessor. Note that errors in the calculation amplify in higher orders. However, the standard deviations 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,

Modelling case study
Models of root architecture and function increase our mechanistic understanding of soil-plant interactions. For example, Schnepf et al. We performed a topological characterisation of 100 simulated root systems according to Fitter (1987). The analysis reveals that, after 25 days, the mean altitude is 73.60 ± 3.28, the mean magnitude is 240.09 ± 36.95 and the mean external path length is 8376.65 ± 1667.71. Table   Table 3 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 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 3 dimensions and constrained within the mesocosm as in the experimental setup. This is illustrated in  Table Table   2. The asterisks with error bars show the corresponding measured root lengths and number of roots according to Table Table 1 observe that initially the average root length is approximately 1 cm; this increases over time to an average length of approximately 2 cm.

RootReader2D
Several image analysis tools for 2-dimensional root systems are available. Our algorithm combines features that are partly present in other softwares. It shares the skeletonisation 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 Fig. Error! Reference source not found.a and analysed it with both software for comparison.
In Table Table 4 The main parameters such as overall root length, internodal distance on 0 th -order root and branching angles agree well. Thus, we suggest that both produce satisfactory results. However, there are two main differences in the results: Firstly, 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 (in this case 2 nd -order) roots is larger and thus the internodal distance (in this case on the 1 st -order roots) is smaller. This has implications for parametrisation of root architectural models.  Thus, this tool is more suitable for the extensive analysis of larger root systems. Results of analysing a large root system are shown in the next section.

Application to a large maize root system
"Root System Analyser is so far the only algorithm based on imaging techniques that was applied to root systems as complex as that of a 78 days old maize root system. There is the potential to use this algorithm to extract more information from the large number of existing 2D 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 the 2-dimensional drawing of a large maize root system ( Kutschera et al., 2009, p. 164). Assuming the sowing time to be May 1 st , the root system is 78 days 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 (Lichtenegger, 2003) was scanned at high resolution (1200 dpi) yielding an image size of 116 megapixel. The image was turned into a black-and-white image by thresholding and then analysed by Root System Analyser.
In the tracking procedure, we encountered two problems: Firstly, 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. Secondly, 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 3457 roots automatically.
The tracking results agreed well with data from 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 4, however most of the roots belong to orders 0-3. The original image and the tracked root system are shown in Figure Error! Reference source not found.. Number and length of roots in each order that were detected in this 2-dimensional view of the root system are provided in Table Table 5. Selected parameters of orders 0-3 are presented in Table 6, dynamic parameters r and k could not be estimated from only one image. Average root architectural parameters for the first three orders compare well to those presented in  with the exception of the the length of the apical zone of 0-order roots.

Conclusions and outlook
We presented a novel approach for root tracking from 2-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 developmental point of view. Particularly, this is beneficial for dealing with overlaps and intersections in 2-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 skeletonisation 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 2-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 dynamical development of the root system can be analysed, like elongation rates for individual soil grown roots.
Our approach works on the segmented images of root systems. In Plant root experiments are extremely costly regarding labour 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 parametrisation of root architectural models.

Materials and methods
Starting point for the automated root tracking algorithm is a segmented 2-dimensional image. In order to annotate coherent roots from this image, we make use of morphological and graph theoretic methods to

Binary 2-dimensional neutron radiography images of root systems of Lupinus albus
Neutron radiography is one of the few non-invasive and non-destructive techniques available for studying plant root systems in situ (Moradi et al., 2009)

From the binary image to a graph
In a first step we derive an approximated centreline 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 neighbouring skeleton pixels in an 8-neighbourhood. This yields a matrix im C . If a pixel of im C has exactly one neighbour it represents a leaf in the corresponding graph. If it has more than two neighbours, it represents a branching or crossing point. Therefore all pixels that correspond to nodes in the graph are given by

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 largest z-coordinate.
The time that has passed at each measurement time is subdivided into smaller time steps typically used for root growth modelling. 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 which 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 fulfils 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.

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 in section 3.2.1. 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 favour 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, x 2 is located at R a along the new path after the current root tip. Therefore, the straightness is given by We 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 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 removal of wrongly assigned roots as well as manual assignment by picking the starting and end points of the root. Figure  that solves the single source shortest path problem for a graph with positive edge path costs. In our application this 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 PC, the analysis of young root systems takes only takes 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 (11658 nodes, 14783 edges, 3478 detected roots) may take a few hours on a standard PC. In mature root systems there are limitations with respect to the complexity of root system. The part just below the stem is so densely rooted that individual roots can no longer be recognised. As a results 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 Analyser 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 non-trivial task for most imaging techniques due to different types of artefacts. Depending on the imaging technique that is used there might be specialised software available for the segmentation step.

Model parametrisation
We use the resulting data structure to parametrise the dynamic root architectural model of Leitner et al. (2010); however, any other root architectural model could be parametrised as well. This model needs 12 input parameters as shown in Table Table 2. 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 Analyser, 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 in end of section 3.2.2. 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