Statistical Analysis and Modeling of the Geometry and Topology of Plant Roots

TL;DR

This paper introduces a statistical framework for analyzing and modeling the geometry and topology of plant roots, aiding in understanding growth patterns and environmental effects.

Contribution

The authors develop a novel tree-shape space representation and metric for root structures, enabling statistical analysis, geodesic computation, and applications like classification.

Findings

Framework effectively captures root variability

Application to wheat roots demonstrates environmental influence

Framework supports classification and regression tasks

Abstract

The root is an important organ of a plant since it is responsible for water and nutrient uptake. Analyzing and modelling variabilities in the geometry and topology of roots can help in assessing the plant's health, understanding its growth patterns, and modeling relations between plant species and between plants and their environment. In this article, we develop a framework for the statistical analysis and modeling of the geometry and topology of plant roots. We represent root structures as points in a tree-shape space equipped with a metric that quantifies geometric and topological differences between pairs of roots. We then use these building blocks to compute geodesics, i.e., optimal deformations under the metric between root structures, and to perform statistical analysis on root populations. We demonstrate the utility of the proposed framework through an application to a dataset of wheat roots grown in different environmental conditions. We also show that the framework can be used in various applications including classification and regression.

Figures

Tables

Table 1. TABLE I: Comparison of the computation time (in seconds) between the proposed approach and the approaches of Feragen et al.   [ 22 ] and Wang et al.   [ 24 ] . ”Size” refers to the number of lateral roots.

	Size	Approach	Alignment	Correspondence	Geodesic	Overall
Fig. 4 - Example 1	$23$	This paper	$147.0$	N/A	$0.04$	$147.04$
	$23$	Feragen et al. [22]	$147.0$	N/A	$10.98$	$157.98$
	$23$	Wang et al. [24]	$147.0$	$0.0050$	$11.40$	$158.40$
Fig. 4 - Example 2	$14$	This paper	$137.0$	N/A	$0.04$	$137.04$
	$14$	Feragen et al. [22]	$137.0$	N/A	$10.60$	$147.60$
	$14$	Wang et al. [24]	$137.0$	$0.0028$	$11.31$	$148.31$

Equations

$$q(t)=\{\begin{array}{ll}\frac{\dot{\beta }(t)}{\sqrt{\parallel \dot{\beta }(t)\parallel }}& \text{ if }\dot{\beta }(t)\text{ exists and is nonzero,}\\ 0& \text{otherwise.}\end{array}$$

$$d_{\mathcal{C}}({\mathbf{q}}^1,{\mathbf{q}}^2)={\lambda }_m\parallel q_0^1-q_0^2{\parallel }^2+{\lambda }_s\sum _{k=1}^N\parallel q_k^1-q_k^2{\parallel }^2+{\lambda }_p\sum _{k=1}^N(s_k^1-s_k^2{)}^2,$$

$$d({\mathbf{\beta }}_1,{\mathbf{\beta }}_2)=\underset{O\in SO(2),\gamma \in \Gamma }{\min }{d}_{\mathcal{C}}({\mathbf{q}}^1,O({\mathbf{q}}^2,\gamma )),$$

$$\mathbf{\alpha }(r)=(1-r){\mathbf{q}}^1+r{\mathbf{q}}^2,0\le r\le 1,$$

$$\mathbf{\mu }={\mathrm{argmin}}_{\mathbf{q}\in \mathcal{S},O\in SO(2),\gamma \in \Gamma }{d}_{\mathcal{C}}(\mathbf{q},O({\mathbf{q}}^i,\gamma )).$$

$$\mathbf{q}={\text{Exp}}_{\mathbf{\mu }}\left(\sum _{i=1}^m{b}_i\sqrt{{\lambda }_i}{\Lambda }_i\right).$$

$$\text{M}[p_1,p_2,\dots ,p_l,1{]}^{\top }=\text{b}.$$

$$\text{M}=\text{B}{\text{P}}^{+},$$

Full text

**Keywords: ** geodesics, tree-shape space, root shape alignment, shape analysis, principal component analysis, shape modeling, root shape synthesis

1 Introduction

Roots are the primary site of nutrient and water uptake for plants and play a crucial role in plant growth. Understanding structural (*i.e., *geometric and topological) similarities and dissimilarities between roots, and capturing and modeling the structural variability in root populations can help in assessing the plant’s health, understanding its interaction with the surrounding environment, and modelling growth processes.

Existing techniques for modelling structural variability are limited to objects that only bend and stretch [1, 2, 3, 4, 5]. Early works, *e.g., * [1], use morphable models, which represent the shape of an object using a dense set of landmarks and thus it can be seen as a point in a high-dimensional Euclidean space. Statistical analysis can then be performed using standard techniques such as principal component analysis [6]. However, these techniques are limited to objects which only slightly bend or stretch. Recent works such as [2, 3, 7, 4] proposed new formulations that are suitable for large elastic deformations. These techniques, however, do not handle topological variations such as those encountered in natural objects that have a tree structure111By a tree structure or a tree-shape, we refer to the mathematical concept of tree, *i.e., *a graph of nodes and edges with no cycles. Each edge can have attribute that define its shape. As such, although a botanical tree and a root are two different things, mathematically they have a tree-shape or a tree structure. , *e.g., *plant’s shoots or roots, airway trees, and neuronal structures in the human brain. In these types of objects, growth processes, disease progression, and environmental effects affect not only their geometry but also their branching structures.

In this article, we propose a framework for the statistical analysis of the geometry and topology of objects that have a tree structure. We focus on plant roots characterised by a main root and first order lateral roots. Higher order lateral roots (sub laterals and off laterals) are not catered for in the present form, although the model can be extended to treat these as well. In the meantime, we foresee that actual and more complex root systems are dissected to attain the form applicable to the analysis presented here. The framework we propose builds on the representation suggested by Duncan et al. [8] for the analysis of simple neuronal structures. It allows one to:

• •

Compute correspondences and geodesic paths between plant roots even when the latter undergo large bending, stretching, and topological deformations.

• •

Compute statistical atlases, *i.e., *the means and principal modes of variation, of plant root collections.

• •

Characterize the geometric and structural variability within a collection of roots using probability distributions.

• •

Develop a mechanism a program in which plant roots are synthesized either randomly, using random sampling, or in a controlled manner using regression from biologically motivated parameters.

The proposed framework has a wide range of applications. We show its utility in root classification and root synthesis. It also has multiple applications in plant biology including (1) quantifying differences in root morphology, and (2) comparing root systems of either different genotypes grown under the same environmental conditions, or root systems of the same genotype plant grown under different environmental conditions, including different nutrient concentrations, or soils at different levels of moisture content or different toxicity.

The remaining parts of this article are organised as follows; Section 2 overviews the related work. Section 3 lays down the mathematical formulation of the concepts tree-shape space and its associated metric and geodesics. These concepts are then used to perform statistical analysis of collections of plant roots (Section 4), and to synthesize and simulate 3D plant roots (Section 5). Section 6 presents the results, and compares the performance of the proposed approach and the quality of its results with the state-of-the-art. Section 7 summarizes the main findings of this article.

2 Related works

The framework developed in this article can be seen as the generalization to objects with tree-like structures (*e.g., *plant roots and shoots) of the statistical analysis techniques that have been proposed for manifold 3D shapes. These types of objects vary not only in geometry but also in topology. Thus, we focus our survey on the techniques which have been proposed for capturing root morphology, and on techniques for statistical shape analysis.

2.1 Root morphology analysis

Quantitative characterization of root shapes can reveal fine differences between various phenotypes, signals, and exogenous regulatory substances [9]. This can facilitate molecular-genetic and physiological studies of root development. Some approaches extract shape descriptors to capture the important geometric and morphological properties of root shapes. Examples of such properties include the angular deviation of the root tip from vertical axis [10, 11, 12], the vertical growth index [9, 13], which is measured as the ratio between the root tip ordinate and the root length, and the horizontal growth index defined as the ratio between the root tip abscissa and the root length. Schultz et al. [14] employed straightness, wave density, and horizontal growth index to describe the root shape.

These shape descriptors have been used to represent and compare root shapes in the descriptor space. They, however, only represent spatial information about roots. Also, variability in the descriptor space often does not correspond to variability in root shapes. For example, the average of two root descriptors is not guaranteed to correspond to the average root, not even to a valid root shape.

2.2 Statistical Shape analysis and modeling

Instead of using descriptors to characterize the shape of objects, several recent papers treat the shape of an object as a point in a shape space. Equipping the shape space with a proper metric allows for a comparison of objects based on their shapes, computing geodesics, and performing statistical analysis including regressions and shape synthesis. These concepts have been extensively investigated in the case of 3D manifold objects such as human faces [1], human bodies [15], arbitrary natural objects [2, 3], and man-made objects [16]. These methods, however, are limited to 3D models with fixed topology. They cannot capture and model structural variabilities such as those present in plant roots.

More closely related to our approach are techniques based on tree statistics. The seminal work of Billera et al. [17] proposed the notion of continuous tree-space and its associated tools for computing summary statistics. Some variants of this idea have been introduced for the statistical analysis of tree-structured data, *e.g., * [18, 19]. These works, however, only consider the topological structure of trees and ignore the geometric attributes of edges, which limits their usage. To overcome this limitation, several extended methods have been proposed for defining a more general tree-space. This includes Feragen et al.’s framework [20, 21, 22, 23], which proposed a tree-shape space for computing statistics of airway trees, and its extension to complex botanical trees [24, 25].

Despite their efficiency and accuracy in certain situations, these techniques exhibit three main fundamental limitations. First, they use the Quotient Euclidean Distance (QED), which is not suitable for capturing large elastic deformation, *i.e., *bending and stretching, of the branches. Second, they represent tree shapes as father-child branching structure, which leads to a significant shrinkage along the geodesics between trees that exhibit large topological differences. Third, branch-wise correspondences need to be manually specified, especially when dealing with complex tree-like structures, which restricts its utility in practical applications.

**Contributions: ** We propose in this article a statistical framework that is more suitable for analyzing plant root shapes. It builds upon and extends the recent work of Duncan et al. [8]. First, we show that the proposed main-side branching representation is more efficient for capturing topological changes than the father-child branching employed in the frameworks of Feragen et al. [20, 21, 22, 23] and Wang et al. [24, 25]. Second, we use the new representation to develop tools for computing, jointly, one-to-one correspondences and geodesic deformations between plant roots that significantly differ in geometry and topology. Finally, we demonstrate the utility of the framework in high-level applications such computing statistical summaries and atlases, synthesising root shapes, either randomly or in a controlled manner, and classifying plant roots based on their geometry and topology. To the best of our knowledge, this is the first approach that deals with statistical modeling of plants root based on their shapes.

3 Mathematical representation and shape space metric

The input to our framework is a collection of 2D plant roots. We first skeletonize each root and convert it into a set of two-layer planar curves: ${\bm{\beta}}=(\beta_{0},\{\beta_{i}\}_{i=1}^{n})$. Here, $\beta_{0}$ represents the main root branch and $\{\beta_{i}\}_{i=1}^{n}$ represent the finite (possibly empty) set of lateral branches, each one is attached to the main root at some locations, see Figure 1. Each branch $\beta_{i}$ is a continuous curve in $\mathbb{R}^{2}$, i.e., $\beta_{i}\ :[0,1]\rightarrow\mathbb{R}^{2}$. Let $t_{i}\in[0,1]$ such that $\beta_{0}(t_{i})=\beta_{i}(0)$ is the location on $\beta_{0}$ to which the $i-$th lateral branch $\beta_{i}$ is attached.

In order to perform statistical analysis of the geometry and topology of a collection of such roots, we need to define a tree-shape space222A tree is a graph with no cycles. It is defined as a set of nodes and edges that connect them. Each edge can also be augmented with shape attributes. , a metric on this space, and a mechanism for computing correspondences and geodesics in the tree-shape space. Due to their arbitrary structure, two roots will rarely have the same number of lateral roots. To simplify the search for correspondences, and following [8], we augment each root $\bm{\beta}$ by adding virtual lateral branches, *i.e., *branches of length zero, according to the $t$ value of each non-virtual lateral branch. For example, for two roots ${\bm{\beta}_{1}}=(\beta_{0}^{1},\{\beta_{i}^{1}\}_{i=1}^{n_{1}})$ and ${\bm{\beta}_{2}}=(\beta_{0}^{2},\{\beta_{i}^{2}\}_{i=1}^{n_{2}})$, we add $n_{2}$ null branches to $\bm{\beta}_{1}$, whose locations on $\beta_{0}^{1}$ are computed as $t^{1}_{n_{1}+i}=t^{2}_{i},i=1,2,\dots,n_{2}$. Then, we add $n_{1}$ null branches to $\bm{\beta}_{2}$ in the same manner. In what follows, and for simplicity of notation, we will also use the symbol $\bm{\beta}$ to indicate an augmented root.

We represent each root $\bm{\beta}$ using its Square Root Velocity Function Tree (SRVFT) [26], $\bm{q}$, which consists of a collection of Square Root Velocity Functions, one for each branch $\beta$ in $\bm{\beta}$. That is, given one branch $\beta$, its SRVF representation is defined as its derivative scaled by the square-root of the $\mathcal{L}_{2}$ norm of the derivative:

[TABLE]

As such, the SRVF is translation invariant. Thus, when converting an entire root into its SRVFT representation, we lose track of the location of the bifurcation points. In order to preserve them, we represent each side branch using an ordered pair $(q_{i},s_{i})$ where $s_{i}\in[0,1]$ is the starting location on $\beta_{0}$. Then, the SRVFT representation of the entire root $\bm{\beta}$ becomes $\bm{q}=(q_{0},\{(q_{i},s_{i})\}_{i=1}^{N})$.

Let $\mathcal{C}$ denote the space of all root shapes, which are represented by their SRVFTs. $\mathcal{C}$ is called pre-shape space [26]. We now need to define the metric and the mechanism for computing geodesics between a pair of roots.

3.1 Metric for tree-like shapes

Our goal is to define a metric, which quantifies the geometric, *i.e., *bending and stretching, and topological deformations. We follow the approach of Duncan et al. [8]. Let $\bm{q}^{i}=(q_{0}^{i},\{(q_{k}^{i},s_{k}^{i})\}_{k=1}^{N}),i\in\{1,2\}$ be the SRVFT representations of two roots $\bm{\beta}_{1}$ and $\bm{\beta}_{2}$, respectively. We define the distance, or dissimilarity, between such roots in $\mathcal{C}$ as:

[TABLE]

which is a weighted sum of three terms. The first term is the $\mathcal{L}_{2}$ distance between the main roots. Since the $\mathcal{L}_{2}$ distance in the SRVF space is equivalent to the full elastic metric in the space of curves [26], then the first term measures the amount of bending and stretching needed to align one curve onto the other. The second term is the $\mathcal{L}_{2}$ distance between the SRVFs of the lateral roots. The last term measures the distance between the bifurcation points and is used to capture topological changes. The parameters $\bm{\lambda}=(\lambda_{m},\lambda_{s},\lambda_{p})$ control the relative contributions of the three terms.

A good metric for statistical shape analysis should be invariant to shape-preserving transformations, *i.e., *translation, scaling, rotation, and re-parameterization. The SRVF representation is translation-invariant by construction since it uses derivatives. Invariance to scale can be efficiently handled by scaling each root by the length of its main root. Note that the latter might not be needed depending on the application at hand. For instance, growth analysis requires preserving the scales of the roots.

The remaining variabilities, *i.e., *rotation and reparameterization, are handled algebraically following [4, 8]. The action of all possible rotations $O\in SO(2)$ and reparameterizations $\gamma\in\Gamma$ on a root shape $\bm{\beta}$ forms a set of roots that have the same shape ($\Gamma$ here is the space of all orientation-preserving diffeomorphisms of $[0,1]$ to itself). Similar to [4], we redefine the dissimilarity between two roots $\bm{\beta}_{1}$ and $\bm{\beta}_{2}$ as the length of the shortest geodesic between $\bm{q}^{1}$ and $O(\bm{q}^{2},\gamma)$ (here, $(\bm{q},\gamma)$ is the SRVFT of $\bm{\beta}\circ\gamma$):

[TABLE]

where $d_{\mathcal{C}}(\bm{q}^{1},O(\bm{q}^{2},\gamma))$ is the length of the geodesic between $\bm{q}^{1}$ and the rotated and re-parameterized version of $\bm{q}^{2}$. The optimization over $O$ and $\gamma$ is the registration process, which consists of searching for optimal alignment and correspondence between $\bm{\beta}_{1}$ and $\bm{\beta}_{2}$. It is solved as a linear assignment problem following [8].

3.2 Geodesics

A geodesic is the optimal deformation (bending, stretching, and topological changes), under the metric, from one root to another. Its length quantifies the minimum amount of deformation that one needs to apply to one root in order to align it onto the other. For simplicity, let us also denote by $\bm{q}^{2}$ the version of $\bm{q}^{2}$ after applying the optimal rotation and reparameterization found by optimizing Equation (3). Since the distance between $\bm{q}^{1}$ and $\bm{q}^{2}$ is a weighted sum of Euclidean distances, see Equation (2), then the geodesic $\bm{\alpha}$ connecting them is the linear path that connects the two points, *i.e., *:

[TABLE]

which is defined in the SVFT space. The geodesic between $\bm{\beta}_{1}$ and $\bm{\beta}_{2}$ is then given by mapping the path $\bm{\alpha}$ back to the space of trees.

Figure 2 shows an example of the geodesic between two simple tree-like shapes that differ in topology. For this example, we also show the results obtained by the approach of Feragen et al. [22] and the approach of Wang et al. [24]. Each row in the figure is one geodesic between the leftmost and rightmost roots. Branchwise correspondences are indicated using colors. In this example, we can see that the intermediate shapes along the geodesic obtained using [22] are not natural; they exhibit unnatural shrinkage and are unable to find the correct correspondences. The approach of Wang et al. [24] is less prune to global shrinkage. However, as seen in the last row of Figure 2, it can match branches located in different sides. As such, the lateral branches unnaturally shrink and expand along the geodesic.

4 Statistical atlas

For a given set of root samples, denoted as $\{\bm{\beta}_{i},i=1,2,...,m\}$, the mean root is defined as the root that is as close as possible to all the roots in the set. Mathematically, it is the root that minimizes the sum of the geodesic distances to all the roots. Similar to geodesic computation, we first add to every root $\bm{\beta}_{i}$ virtual lateral branches to equalize the number of lateral branches across the roots in the dataset. Let $\bm{\mu}$ be the SRVFT representation of the mean root. It can be computed as:

[TABLE]

Here, $\mathcal{S}$ is the tree-shape space. The solution to Equation (5) is known as the Karcher mean. We employ the same gradient descent approach described in [4] to solve this optimization problem.

Figure 3(a) shows the mean root of three roots generated by this approach. Figure 3(b) compares between the mean roots generated by our approach, the approach of Feragen et al. [22], and the approach of Wang et al. [24]. It can be clearly seen that the mean root generated by our approach is more plausible than those generated by the other two methods.

In addition to the mean root, which can be regarded as a template that characterizes the main morphological properties of the shape of roots in the dataset, one would also like to (1) quantify how the samples in the dataset deviate from the mean root, and (2) analyze the distribution of roots around the mean root. This can be done by computing the covariance matrix and the modes of variation. The tree-shape space $\mathcal{S}$ is a non-linear manifold, thus we employ its tangent space $T_{\bm{\mu}}(\mathcal{S})$ instead of $\mathcal{S}$ to perform statistical analysis. Similar to [26, 4], we first map each $\bm{q}^{i}$ onto the tangent space using the inverse exponential map: $\bm{q}^{i}\rightarrow\bm{v^{i}}=\log_{\bm{\mu}}(\bm{q}^{i})$. Then, we calculate the covariance matrix: $C=\frac{1}{m-1}\sum_{i=1}^{m}\bm{v^{i}}\bm{(v^{i})^{t}}$. After calculating the eigenvalues $\{\lambda_{i}\}$ and eigenvectors $\{\Lambda_{i}\}$ of $C$, we can obtain the principal modes of shape variation for the root collection. We also can generate a sample $v_{k}\in T_{\bm{\mu}}(\mathcal{S})$ along each eigenvector $\Lambda_{i}$ as: $v_{k}=\alpha\sqrt{\lambda_{i}}\Lambda_{i},\alpha\in\mathbb{R}$ on the tangent space, where $\lambda_{i}$ is the eigenvalue associated with the eigenvector $\Lambda_{i}$.

5 Plant roots synthesis

5.1 Root synthesis by random sampling

The eigenvectors $\{\Lambda_{i}\}$ form an orthonormal basis of a Euclidean vector space. One can characterize the distribution of the input roots by fitting probability models either from a non-parametric family by computing their probability densities directly from the data, or from a parametric family, such as multivariate Gaussians. In this article, for the sake of simplicity, we fit to the population a multivariate Gaussian with mean $\bm{\mu}$ and a diagonal covariance matrix $C$ whose diagonal elements are the square roots of the eigenvalues $\lambda_{i}$.

Once we have the probability model, a new botanical tree can be generated by randomly sampling from the distribution. In the case of a Gaussian model, a random root can be synthesized by first randomly generating $m$ real values $b_{1},b_{2},...,b_{n}\in\mathbb{R}$ and then setting

[TABLE]

Here, $\text{Exp}_{\bm{\mu}}$ is the exponential map, which maps elements in the tangent space to $\mathcal{S}$ at $\bm{\mu}$ to the SRVFT space. A root $x$ can be obtained by mapping $\bm{q}$ back to the tree-shape space. In this article, we only consider the $m$-leading eigenvalues such that $\frac{\sum_{i=1}^{n}\lambda_{i}}{\sum_{i}\lambda}>0.99$. In order to generate plausible roots, one can restrict $b_{i}$’s to be within a certain range, e.g., $[-1,1]$.

5.2 User-controlled root synthesis

While random sampling from the multivariate Gaussian distribution provides an easy way of synthesizing new roots, it lacks control. In fact, sometimes, users would like to generate roots by just adjusting a few parameters, *e.g., *biologically motivated parameters. In this article, we formulate this problem as a regression in the tree-shape space [25]. In particular, let $\textbf{p}\in\mathbb{R}^{l}$ be the vector of $l$ parameters. Let $x$ be a point in the tree-shape space and $\bm{q}$ its SRVFT representation. From Equation (6), $\bm{q}$ can be represented as a real valued vector $\textbf{b}=(b_{1},b_{2},\dots,b_{m})^{\top}$. If we assume that the relation between the biological parameters p and the root to be linear, then the mapping can be represented using $m\times(l+1)$ matrix M such that:

[TABLE]

After assembling all the parameter vectors into an $(l+1)\times m$ matrix P and all the vectors b into an $n\times m$ matrix B, the mapping matrix M can be calculated as follows;

[TABLE]

where $\textbf{P}^{+}$ is the pseudoinverse of P.

There are several biologically-motivated parameters that can be used. In this article, we consider (1) the main root length, (2) the mean length of the side-roots, and (3) the standard deviation of the side root lengths. We first extract these three parameters from a training dataset and use them to learn the regression model. At runtime, the user can specify these parameters and the system automatically synthesizes new root models.

6 Results and discussion

In this section, we consider a collection of wheat roots that exhibit different degrees of geometric and structural variations, and demonstrate the results of the proposed approach in computing (1) geodesics and (2) statistical summaries such as means and modes of variations. We also report the timing and compare our results to the state-of-the-art. Additional results are included in the supplementary material. For this, we have collected $53$ wheat plant roots with rich structural variations. First, wheat plants have been taken from soil, washed, and scanned using a flatbed scanner. Their images have then been binarized and automatically skeletonized and converted into a two-layer structure representation as described in Section 3.

6.1 Geodesics

In this experiment, we consider pairs of plant roots, each pair is composed of a source and a target root, and generate new roots by computing the geodesic (optimal deformation) that connects the source to the target. Figure 4 shows two examples of such geodesics. For comparison, we also show geodesics between the same pair of roots but computed using the approaches of Feragen et al. [22] and of Wang et al. [24]. In each example, the colors indicate the branch-wise correspondences.

As one can see, the geodesics generated with the approach of Feragen et al. [22] exhibit unnatural shrinkage: the intermediate roots along the geodesic shrink and then expand. Also, the approach of Wang et al. [24] fails to find correct branch-wise correspondences, which significantly affects the quality of the geodesics obtained with this approach. In comparison, the approach proposed in this article is able to produce smooth geodesics with plausible-looking intermediate roots. More results are included in the supplementary material.

In addition, we check the influence of the weight $\lambda_{m},\lambda_{s}$, and $\lambda_{p}$ of Equation (2) on the quality of the geodesics. For the purpose of this experiment, we vary the values of these parameters and observe the results. This is illustrated in Figure 5 on a toy example. In the first three results where the triplet $(\lambda_{m},\lambda_{s},\lambda_{p})$ is set to $(1.0,0.01,0.01),(0.01,1.0,0.01),(0.01,0.01,1.0)$, respectively, the three geodesics look very similar. However, in the fourth result where $\lambda_{m}=0.01$, $\lambda_{s}=0.00001$, and $\lambda_{p}=1.0$, the approach favours the creation of new lateral branches rather than sliding the existing one. This is predictable since the last term, which measures distance between bifurcation points, is highly penalized.

The supplementary material also includes additional complex examples that show the effect of varying the parameters of Equation (2) on the quality of the geodesics.

6.2 Summary statistics

Figure 6 shows the mean root computed from the entire data set. The 53 input roots are shown on the left while the generated mean root is shown on the right. Note that the correspondences between each pair in the dataset, as well as between each sample in the dataset and the computed mean root, are automatically computed and do not require any user input or interaction. Also, the average of a pair of plant roots is a by-product of the geodesic computation process. In fact, the root that is equidistant to the source and the target roots is exactly their mean. Thus, the root that lies exactly at the middle of each geodesic in Figure 4 is actually the average of the leftmost and rightmost roots.

Figure 7 shows the first two principal modes of variation for the root dataset of Figure 6. From these results, we clearly see that the first two leading modes capture the main geometric and topological variations in the root dataset.

6.3 Random root synthesis

Figure 8 shows plant roots that have been automatically synthesized by random sampling from the statistical model fitted to the root dataset of Figure 6. These roots have been synthesized without any biological knowledge or interaction from the users. For a better visualization, we group these synthesized plant roots into three clusters depending on their distances to the mean root. Note that those that are close to the mean share some similarities with the input roots but are not the same. The further the synthesized roots deviate from the mean the less plausible they become.

6.4 User-controlled root synthesis

We show a few examples of plant roots generated with intuitive controls such as main root length, mean length of side-roots, and standard deviation of side roots length. We take the plant roots of Figure 6, compute their mean root and their modes of variation, and then learn a regression model that maps the three control parameters into points in the tree-shape space. This provides a direct way to explore the range of root shapes. Figure 9, 10, 11 show a few representative results produced with this approach.

6.5 Timing

The propsoed framework has been implemented in Matlab and runs on a desktop PC with Intel(R) Core(TM) i5-4460 CPU@3.20GHz and 8GB of RAM. Table I reports the computation time for each geodesic shown in Figure 4. It also reports the side root number of the roots after adding virtual side roots.

To compute each of the geodesics of Figure 4, we first align the source and target roots to one another (alignment step), and then compute their geodesic using the method proposed in this paper, the method of Feragen et al. [22], and the method of Wang et al. [24]. Table I shows that most of the time is consumed in the alignment process. Once the roots are aligned, the method proposed in this article for computing geodesics is significantly faster than the approaches of Feragen et al. and Wang et al. ($40$ms for ours vs. $11$s for the other two methods).

The computation of the mean root takes $9$ hours $4$ minutes. Once the mean root has been computed, computing the principal modes of variation takes $1.1$ seconds. It only takes $0.22$ seconds to generate one random root. For regression, it only tales $0.01$ seconds to compute the mapping matrix $M$ and about $0.001$ seconds (on average) to generate a new root. This suggests that the approach proposed in this article is well suited for the generation and synthesis of root structures

6.6 Application to plant root classification

Finally, the length of the geodesic between two roots is a measure of dissimilarity between these roots. To demonstrate the utility of this measure for plant root classification, we computed the pairwise geodesic distance matrix for the entire dataset and performed a hierarchical binary clustering using MATLAB’s linkage function. Figure 13 shows the root clustering result according to the hierarchical binary clustering result as is shown in Figure 12. As one can see, the metric clusters together roots that have similar geometry and topology.

7 Conclusion and Future Work

In this article, we have developed a full framework for the statistical analysis and modeling of the geometry and topology of plant roots. Central to the framework is an elastic curve-based tree-shape space and an elastic metric that enables us to find optimal correspondences between roots, even in the presence of large topological deformations. As shown in the experimental results, the metric is suitable for comparing roots based on their shape, and for computing geodesics between them. We showed that the representation and metric can be used to compute statistical summaries of root collections and to synthesize new roots either randomly or interactively by varying biologically-motivated parameters. This can lead to many practical applications. For instance, the proposed framework can be used alongside with root imaging techniques, which are recently receiving a growing attention, see for example [27] for a repository of root datasets. The proposed geodesic computation algorithm and geodesic distance can be used to compare Root System Architectures (RSA) and to understand and model similarities and differences of RSAs within and across plant species. The tools for computing statistics (means and modes of variations) can be used to model the variability in geometry and topology of RSAs within species. It can also be used for probabilistic classification. Finally, regression tools are fundamental to the simulation and modelling of RSAs from biological or environmental parameters. To the best of our knowledge, this is the first paper that provides all the building blocks for these high-level tasks.

Although effective as evidenced by the results shown in this article, there are still some limitations that can be addressed in future work. First, we only considered roots that have two layers of branches. However, roots have complex tree structures often composed of more than two layers. In the future, we plan to extend the current framework to processing multi-layer roots. Second, the current formulation does not take into account the thickness of the roots. As such, the framework does not differentiate between coarse and fine roots. While root thickness can be easily taken into account by adding a fourth term to Equation (2), it will result in an additional parameter to tune and an increase in computation time. As such, we plan in the future to explore more efficient techniques. Third, the type of geodesics that our framework generates between a pair of roots depends on the weights of the three terms of the metric defined in Equation (2). While in this article, we manually set these weights, in practice they depend on the application and on the plant species being analysed. As such, it would be interesting to learn these parameters from data. Third, we have demonstrated that this framework can be used to regress plant root shapes from a few parameters. Our framework is general and can be easily extended to incorporate biological knowledge and environmental effects such as soil humidity and nutrient content. In the future, it will be interesting to incorporate these factors into the regression process and explore more regression tools beyond the linear ones.

Finally, we plan in the future to use the proposed framework to quantify differences in 3D root morphology, and comparing 3D root systems of either different genotypes grown under the same environmental conditions, or 3D root systems of the same genotype plant grown under different environmental conditions, including different nutrient concentrations, or soils at different levels of moisture content, or soils of different toxicity. This could be achieved either simply by calculating the geodesic distance between mean roots, or by using advanced statistical measures that capture variability within each group of plants.

**Acknowledgement. ** The authors would like to thank Adam Duncan for sharing the code for analyzing simple neuronal trees, and Aasa Feragen for making the code for tree-like shape analysis publicly available. Guan Wang would like to thank the China Scholarship Council for funding his visit and stay at Murdoch University, Australia.