Spatial and Texture Analysis of Root System Architecture with Earth Mover's Distance (STARSEED)

Purpose: Root system architectures are complex, multidimensional, and challenging to characterize eﬀectively for agronomic and ecological dis-covery. Methods: We propose a new method, Spatial and Texture Analysis of Root System architEcture with Earth mover’s Distance (STARSEED), for comparing root architectures that incorporate spatial information through a novel application of the Earth Mover’s Distance (EMD). Results: We illustrate that the approach captures the response of sesame root systems for diﬀerent genotypes and soil moisture levels. STARSEED provides quantitative and visual insights into changes that occur in root architectures across experimental treatments. Conclusion: STARSEED can be easily generalized to other plants and provides insight into root system architecture development and response to varying growth conditions not captured by existing root architecture metrics and models. The code and data for our experiments are publicly available: https://github.com/GatorSense/STARSEED.


Introduction
Studying plant roots, a key to achieving the second Green Revolution (Lynch, 2007), requires effective characterization and comparisons of root growth and architecture. This includes understanding how genetic and environmental factors impact early root development, not only in terms of biomass and root length but also in terms of architecture and spatial exploration. Current methods, such as WinRhizo, use 2D images of roots to measure individual parameters related to morphology and topology. Several open-source software packages have also been developed to further make accessible the analysis of individual root traits (Pierret et al, 2013;Armengaud et al, 2009). These tools capture information relevant to root system architecture (RSA) characterization. However, they typically provide different parameters each only capturing one specific aspect of RSA such as total root length, surface area, branching angle, link magnitude, among others. (Regent Instruments Inc., Quebec, Canada). These packages are lacking methods for holistic architectural analysis of root systems.
There is a need to develop new techniques that are able to provide indepth spatially explicit RSA characterization from 2D images using a limited number of parameters. In the literature, deep learning has been investigated as methods to study root architectures (Chung et al, 2020;Pound et al, 2017;Yasrab et al, 2019;Yu et al, 2020;Xu et al, 2020). Deep learning approaches have provided mechanisms for automating root detection and segmentation in imagery (Chung et al, 2020;Yasrab et al, 2019;Yu et al, 2020;Xu et al, 2020) as well as localizing and identifying unique features for root phenotyping (Pound et al, 2017). Despite the success of deep learning models, there are several issues including computational costs, the need for very large amounts of labeled data, and a lack of explainability (Gunning, 2017;Adadi and Berrada, 2018).
As an alternative to help in addressing these shortcomings, we propose the Spatial and Texture Analysis of Root System architEcture with Earth mover's Distance (STARSEED) approach, to better characterize RSA (e.g., root architecture and soil exploration). Our approach uses the Earth Mover's Distance (EMD) (Rubner et al, 1998) to quantify the differences among root structures. EMD is intuitive and explainable by providing spatially explicit information on changes in root architecture. The proposed STARSEED approach allows for us to etablish meaningful biological connections between the treatments such as cultivar and moisture level, and the observed RSA as well as provide detailed qualitative and quantitative analysis. STARSEED jointly 1) considers the spatial arrangements of roots in local (i.e., smaller regions of the image) and global (i.e., whole image) contexts, 2) extracts informative features from the root architectures, and 3) gives precise insight for biological and agronomic interpretation.

Related Work
Texture Features for RSA Texture (i.e., spatial arrangement of the pixel values in a raster grid) is a powerful cue that can be used to identify patterns tied to RSA characterization. A straightforward approach to compute texture from an image are counts (i.e., histograms) of local pixel intensities (Cula and Dana, 2001). In addition to this baseline approach, one can compute more complex texture features such as fractal dimension and lacunarity. Fractal dimension is used to quantify the "roughness" in an image (Sarkar and Chaudhuri, 1992) and has been used for root analysis to capture the way roots develop to fill space (Li et al, 2020;Wang et al, 2009;TATSUMI et al, 1989). Fractal dimension is inspired by the concept of self-similarity (i.e., idea that objects are comprised of smaller copies of the same structure) for which roots may develop in a consistent pattern invariant to the environment (Mandelbrot, 1967).
A related feature to fractal dimension is lacunarity, which captures the "gaps" or spatial arrangement in images (Mandelbrot and Mandelbrot, 1982;Plotnick et al, 1993). Since some visually distinct images can have the same fractal dimension, lacunarity provides a features that can aid in discriminating between these distinct textures (Keller et al, 1989;Voss, 1991;Mandelbrot and Van Ness, 1968). Lacunarity will have small values when root architectures are dense and larger values with gaps and coarse arrangement of roots (Keller et al, 1989;Mandelbrot and Mandelbrot, 1982).

Earth Mover's Distance (EMD)
EMD, also known as the Wasserstein-1 distance (Gulrajani et al, 2017), is used to determine quantitative differences between two distributions (Rubner et al, 1998). In STARSEED, we use EMD to measure distances between the spatial distribution of root textures, thus, quantifying root architecture. EMD has several advantages such as allowing for partial matching (i.e., comparisons can be made between representations of different length such as comparing smaller and larger root architectures) and matches our human perception when the chosen ground distances (i.e., distance between feature vectors) is meaningful (Rubner et al, 1998). The computation for EMD is based on the solution to the engineering problem for transportation corridors (Hitchcock, 1941). Essentially, we want to find the minimal amount of "work" to transform one root architecture to another.
EMD can be used to not only compare 1-dimensional distributions, but EMD can also be extended to multi-dimensional distributions such as images.
Images are comprised of pixels and these pixels can be "clustered" or assigned to meaningful groups based on shared characteristics such as spatial location. The set of these clusters are used to form a "signature" or a more compact representation of the image to increase computational efficiency (Rubner et al, 1998). Given an image with C clusters, the signature representation, P , is given by Equation 1 where p i is the cluster representative and w pi is the weight of the cluster: Typically, the cluster representative is a feature vector and the weight is the percentage of pixels in a cluster expressing that feature (Rubner et al, 1998). The selection of the cluster representative is application dependent, but by defining the cluster representative as the information/descriptors (i.e., features) and importance (i.e., weight) of a part of an image provides a clear interpretation for EMD. Once the signature representations are constructed, one can compute EMD. After the transportation problem is solved to find the minimal flow, F, between two clusters (please refer to (Rubner et al, 1998) for more details), EMD can be computed. Given two image signatures, P and Q, with S and T representatives respectively, EMD is formulated in Equation 2 where d ij is the ground distance between the centroid of regions p i and q j and f ij is the optimal flow between p i and q j : EMD captures dissimilarity between distributions (i.e., larger values indicate more dissimilarity or 'work' to move the defined 'earth' feature to the cluster representative). EMD allows one to measure the global change between two distributions (i.e., distance measure) as well as local changes between the two sources of information through the flow matrix, F .

Greenhouse Setup and Data Collection
The data and images collected from a root development experiment in a controlled setting was used to develop the STARSEED approach. Example images collected from the experiment are shown in Figure 1. For the experiment, custom-made rhizoboxes were constructed with dimensions of 35.6 × 20.3 × 5.1cm, and each had one side made of clear plastic to allow for observation of root development over time. Each of 64 rhizoboxes were filled with 1550g of inert calcine clay (Turface Athletics, Buffalo Grove, IL, USA), henceforth referred to as soil. A single seed from one of four non-dehiscent sesame cultivars, all provided by Sesaco Inc. (Sesaco32, Sesaco35, Sesaco38 and Sesaco40), was planted per box. In addition, four different soil water content treatments were implemented: 60%, 80%, 100% and 120% of the soil water holding capacity, corresponding to 688, 918, 1147 and 1376mL of water per rhizobox, respectively. The soil and water were mixed thoroughly together before filling the rhizoboxes to promote homogeneity of the soil water content throughout the rhizobox. The top of each rhizobox was covered with Press'N Seal Cling Wrap to reduce water evaporation. A small hole was pierced in the plastic wrap upon seedling emergence to allow for leaf and stem growth. These four cultivars and water treatments were arranged in a complete randomized design with 4 repetitions and were set up in a temperature-controlled greenhouse between 25-35 degrees Celsius. Four runs of the 64 rhizobozes were completed to complete the experiment. Run 1 and 2 were prepared with soil and water only, while run 3 and 4 were fertilized with 1.51mg of 15-5-15 + Ca + Mg Peters Excel mix fertilizer (ICL Specialty Fertilizers, Summerville SC) and 0.29g of ammonium sulfate were dissolved in the water of each rhizobox prior to mixing with the soil. Baking soda was added as needed to neutralize the fertilizer solution to a pH of 6. The rhizoboxes were set up on benches with a 30-degree incline from vertical to promote root growth along the clear plastic side of the box. Plants were grown for a duration of 21 days after planting (DAP) for the runs without fertilizer and 16 DAP for the runs with fertilizer. This duration corresponded to the time it took the first 5 plants in Run 1 without fertilizer and Run 3 with fertilizer to reach the bottom of the rhizobox.
The last day of the experiment, all rhizoboxes were scanned clear side down with a Plustek OpticSlim 1180 A3 flatbed scanner (Plustek Inc., Santa Fe, CA, USA). For each run, images from seeds that did not germinate or seedlings that died during the experiment were taken out of the images set. Eleven images were thus removed from Run 1, 10 from Run 2, 11 from Run 3 and 5 from Run 4. In total, there were 107 images across run 1 and 2, and 113 images for Run 3 and 4. The roots from all remaining images were hand-traced with WinRHIZO Tron (Regent Instruments Inc., Quebec, Canada) and total root length (TRL) was obtained. Fig. 2: Overview of STARSEED. The binary images are reconstructed according to the hand-traced bounding boxes. The images are then downsampled by eight and a spatial grid is overlaid on the image to create equally sized spatial bins. Within each local region, a feature is computed to generate a texture signature of each image. We then compute the Earth Mover's Distance (EMD) between each pair of images to generate a distance matrix from the image chosen as the reference. To visualize the distance matrix and access the selection of scale (i.e., number of spatial regions) for our analysis, we project the roots using multi-dimensional scaling (Kruskal, 1964) for qualitative and quantitative assessment through the Calinski-Harabasz (Caliński and Harabasz, 1974) index to validate our method. Each step in the process corresponds to the same step in Algorithm 1.

Spatial and Texture Analysis of Root System architEcture with Earth mover's Distance (STARSEED)
The overall pipeline of our proposed STARSEED method is illustrated in Figure 2. We provide the details and rational behind our approach in the following subsections. We also provide pseudocode for the overall process in Algorithm 1.

Image Preprocessing
The initial step of the proposed approach is to perform preprocessing of the data to isolate the root pixels in the image. We manually traced roots using WinRHIZO Tron software. Using the WinRHIZO Tron tracings, binary masks were extracted in which pixels that correspond to root were assigned a value of 1 in the mask and the non-root pixels were assigned a value of 0. In addition to root tracing, each image was downsampled by a factor of eight through average pooling and smoothed by a Gaussian filter (σ = 1) to mitigate noise before extracting texture features.

Algorithm 1 STARSEED Overall Process
Require: N input images X ∈ R H×W , N corresponding treatment labels y ∈ R 1 , B number of spatial bins, g θ (·) preprocessing function, f θ (·) texture feature extraction, d embedding dimension 1: Preprocess images, X processed ← g θ (X)

2:
Partition each image into B regions 3: Compute feature and centroid location, X f eatures ← f θ (X processed ) for each region to get a signature, P ∈ R B×3 , for each image 4: Compute EMD score and flow matrix between each image signature using Equation 2 to generate pairwise distance matrix, D ∈ R N ×N

5:
Project D using MDS (Kruskal, 1964) into d-dimensional feature space to compute embedding, E ∈ R N ×d 6: Compute Calinski-Harabasz (CH) score using E and y 7: return EMD Scores, CH Index, Flow Matrices

Novel Earth Mover's Distance Application
In order to apply EMD to characterize and compare RSAs, we construct a "signature" to represent each root image. The signatures used in STARSEED are a composition of textures features computed spatially across each image. Since texture is undefined at a single point (Tuceryan and Jain, 1993), a local neighborhood needs to be identified to compute the texture features comprising the signature. In our proposed approach, we simply divide each image into grid of equally sized regions with each region serving as the local neighborhood. Within each grid square, the collection of texture features being used to locally characterize the RSA is computed. The size of the grid trades off between computational efficiency and localization. A larger grid square corresponds to less computation needed but each texture feature is computed over a larger area, resulting in a loss of localization.
Each region of the image is represented through a cluster representative. The cluster representative is comprised of spatial coordinates (i.e., horizontal and vertical position of region) and a weight to create an R 3 vector to represent a local spatial region of the image. The weight given to each region is computed by the texture feature extracted from the pixels in the region. The final representation of the image will be a "signature" which is the set of clusters from the image. By constructing our signature in this manner, we are incorporating spatial and texture information to represent the root architectures.
We used EMD to compare the signature representations of each root architecture as described in step 4 of Algorithm 1. Only regions containing roots were used for our calculations to mitigate the impact of background on our results. Therefore, we are comparing root signatures of different sizes. In order to not favor smaller signatures, the normalization factor is added to the EMD calculation as shown in Equation 2. Any distance can be used to define the ground distance, d ij (Rubner et al, 1998). In order to consider the spatial arrangement of the root architectures, we compared the center location of each spatial bin by selecting the Euclidean distance as our ground distance to compute d ij . We computed the EMD between each pair of images within Run 1-2 and Run 3-4 separately to generate a pairwise distance matrix, D, as shown in Algorithm 1.

Assessment of Relationships for RSA
Once the distances between each image signature representation is computed, we perform qualitative and quantitative analysis of the results using D and the associated labels of the root structures. Given the pairwise distance matrix, we can use different dimensionality reduction approaches (Oza and Tumer, 2001;Becht et al, 2019;Maaten and Hinton, 2008) to project the matrix into two dimension to visually identify trends in our data as differences among cultivars and moisture treatments. We used multi-dimensional scaling (MDS) (Kruskal, 1964) because this method has been shown to work well with EMD to identify patterns among images that share characteristics (Rubner et al, 1998(Rubner et al, , 1997. In addition to visualizing trends, we want to access the inter-and intrarelationships among the data. Ideally, we expect similar root structures (e.g., roots from the same cultivar) to have a lower EMD values and dissimilar root structures to have larger EMD scores. One measure that can be used to compute these relationships is the Calinksi-Harabasz (CH) index (Caliński and Harabasz, 1974). The CH index is quick to compute and considers both intra-and inter-relationships among the root structures. An analysis of variance (ANOVA) was performed on the CH scores within each treatment for each pair of runs across features. Significantly different means (p < .05) were separated using Tukey's honestly significant difference (HSD) test.

Results
We performed quantitative and qualitative analysis for local and global RSA supported by the STARSEED method. Our dataset was comprised of a total of 220 images across four different runs (2 unfertilized, 2 fertilized, four treatments for cultivar (Sesaco32, Sesaco35, Sesaco38 and Sesaco40) and four water levels (60%, 80%, 100%, and 120%). To measure the intra-cluster dissimilarity (i.e., RSA from the same cultivar or water level) and inter-cluster dissimilarity (i.e., RSA from different cultivar or water level) across treatment levels in each pair of runs, we calculated the CH index (Caliński and Harabasz, 1974) for each level within a fertilizer condition. Ideally, samples that belong to the same cluster should be "close" (i.e., smaller distances or intra-cluster variance) while samples from different clusters should be "well-separated" (i.e., larger distances or inter-cluster variance). When the clusters are dense and well-separated, the CH index is larger. For our global analysis, we generated representative images for each treatment level by averaging the feature representations of each individual image within the same level to carefully explore the root architecture spatial differences across the different cultivars and moisture contents.  For image processing and down-scaling, we investigated the impact of a coarser (i.e., small number of spatial bins) and finer (i.e., large number of spatial bins) scale. The number of spatial bins was varied from 100 to 2000 in steps of 100 (total of 20 values). We calculated three feature values in each local region: percentage of root pixels, fractal dimension (Mandelbrot and Mandelbrot, 1982), and lacunarity (Mandelbrot and Mandelbrot, 1982). To calculate the percentage of root pixels, we simply computed the total number of root pixels divided by the area of the spatial bin. The percentage of root pixels texture provides insight into the distribution of roots in the image. If a region has more dense roots, the percentage of root pixels will be larger.

Local Analysis: Relationships to Treatment Level Representatives
The average and standard deviation of the CH indices across all number of spatial bins for Runs 1 and 2 are shown in Table 1 while the scores for Runs 3 and 4 are shown in Table 2. The small standard deviations indicate that scores are generally stable across the different scales. The results show that our method is robust across the number of spatial bins (i.e., scale) used for our analysis. The ANOVA on the CH scores showed that % root pixels tended to be the feature with the most discriminating power (i.e., the highest CH score) out of the three features considered, except for fractal dimension on cultivar in Runs 1 and 2. For the rest of this study, we will thus display the EMD results for the % root pixels feature, except when presenting results for the cultivar differences in Runs 1 and 2 where we will use the fractal dimension feature. In Runs 1 and 2, the CH indices scores for water levels (1.50 overall average) were not as high (p < .05) as for cultivars (3.96 overall average), indicating more overall variability of root architecture without fertilizer between cultivars than between the moisture levels (Table 1). When looking specifically at each treatment level's score, both the S32 cultivar and the 80% water level have the largest CH indices within cultivar and water level treatments respectively, indicating these treatments as expressing the most distinctly similar root architectures. For Runs 3 and 4 with added fertilizer, the trends were inverse in that the average CH indices were significantly higher across all features and spatial bins for water level (14.70 overall average) than for cultivar (1.62 overall average). The highest average CH score was observed for 120% water level, indicating an extremely distinct RSA for this treatment compared to the roots grown in the three other water treatments ( Table 2).
The present study goes one step further in exploring the EMD results visually. To achieve this, we chose the number of spatial bins corresponding to the largest CH index for each feature as shown in the Table 3. To better capture the changes of RSA, we visualized the root images in a 2D space by projecting the EMD matrix through MDS (Kruskal, 1964). Since the root architectures have the most variance for water levels in Run 3 and 4, we showed a example visualization as shown in Figure 3. Here, we used the percentage of root pixels as the feature and 2000 spatial bins for the images. From the bottom left to the top right, the water level of the root architectures changed from 120% to 60%. Our proposed STARSEED can clearly capture the trend that the main roots (i.e., tap roots) are shorter and have less sub-roots as the moisture level changed.

Global Analysis: EMD Between Treatments
In order to identify global trends in RSA for each treatment, we aggregated the signature representations of each centered root image to create a representative for each treatment. For each root image, we computed the feature representation as shown in Figure 2 from the average feature values within (a) Root Images (b) Water Levels Fig. 3: Example of qualitative results produced by projecting EMD matrix through MDS. The result shown is from the percentage of root pixels as the feature and 2000 spatial bins for images from Runs 3 and 4. The root architectures are arranged in a meaningful way in the 2D space (i.e., images that are similar are near one another). By using MDS, we are able to qualitatively assess if EMD is capturing the relationships between the various root architectures. In Figure 3a, we can see that the MDS projection of the EMD distance matrix groups the root architectures well based off shared characteristics. In Figure 3b, we also observe that the root architectures are well clustered based on moisture level. STARSEED computes the EMD distances between the root architectures without considering the treatment value. Through the proposed approach, we observe that there is some relationship between root architecture differences and treatment levels.
each spatial bin for roots that belong to the same cultivar or water level.  Figure 4e shows the EMD score overall. Figure 4f, Figure 4g, and Figure 4h represent the top 20% value changes between S32 and the three other cultivars; red spatial bins are unique to S32, green spatial bins are unique to S35, S38 or S40, and yellow spatial bins are common to both cultivars.
After we obtained these representatives, we computed EMD between each cultivar or water level. By using the EMD approach, we not only get a score for the differences between each treatment level, but we can also highlight a) the magnitude and b) location of the RSA dissimilarities. The EMD calculations were performed first between each pair of cultivars, then between each pair of moisture levels, and separately between Run 1-2 and Run 3-4. This process generated a large number of different results, thus we chose to display and discuss only a sub-sample of these EMD figures. For each figure, only the top 20% value changes between the reference treatment level (Figures 4-8a) and the other levels (Figures 4-8b, c, d) are displayed as white arrows. We also show the EMD scores in Figures 4-8e between the reference treatment level and the other levels. Figure 4 shows the processed RSA images for each cultivar and EMD results with the reference to S32. The EMD score (Figure 4e) was smallest between S32 and S35, indicating that the two architectures were most similar. This is confirmed by comparing Figure 4f Figure  5e shows the EMD score overall. Figure 5f, Figure 5g, and Figure 5h represent the top 20% value changes between S40 and the three other cultivars; red spatial bins are unique to S40, green spatial bins are unique to S32, S35 or S38, and yellow spatial bins are common to both cultivars. between these two cultivars are relatively small, and for the most part are concentrated in the deeper section of the root system. Specifically, S35 seemed to have a more laterally spread out root distribution towards the deepest portion of the image, and a more concentrated root mass along the tap root right above that. Out of the three other cultivars, S32 was most different from S38, and Figure 4g shows that the most variation in architecture between these two cultivars was located in the middle part of the root system. Indeed, S38 tended to spread out more laterally in the mid-section of the image compared to S32.

Differences between cultivars
The differences between S32 and S40 were located in the upper part of the root system, and Figure 4h shows that S32 tended to have more a more spread out root distribution towards the top of the soil surface compared to S40. According to the Table 1, the EMD results for S40 without fertilizer ( Figure 5 showed that the architecture for S40 was most different from S32, and most similar to S35 (Figure 5e). STARSEED also highlighted the fact that S40 tended to have a pronounced shallow lateral root growth, while the other cultivars tended to be denser and less spread out in the very shallowest depth. When fertilizer was added, S40 showed a tendency to produce more shallow lateral roots than the three other cultivars (Figure 6). S40 was most dissimilar from S32, and the differences consisted of S40 having a stronger lateral growth than S32, which was more vertical. Similar observations can be  Figure 6e shows the EMD score overall. Figure 6f, Figure 6g, and Figure 6h represent the top 20% value changes between S40 and the three other cultivars; red spatial bins are unique to S40, green spatial bins are unique to S32, S35 or S38, and yellow spatial bins are common to both cultivars. made between S40 and S35, though S35 seemed to have more lateral growth than S32. With S38 the differences between the two architectures were not focused on one specific area and tended to be smaller overall, indicative of a more similar architecture also evidenced by the low EMD score between S40 and S38 (Figure 6e).

Differences between moisture levels
For Runs 1 and 2, when considering water levels, the root architectures were overall very similar as shown in Table 1. The architecture that was most distinguished from the others was that of the 80% moisture treatment ( Figure  7). It was overall most different from 100%, and was closest to 120% and the 60% treatment architecture (Figure 7e). Specifically, the 80% treatment had a visibly denser central region than the three other moisture levels (Figure 7a through Figure 7d). This is reflected more clearly in Figure 7f through h: the 80% treatment showed a more concentrated root distribution closer to the top of the root system compared to the other treatments, and so the arrows tend to point toward the upper center of the root system. The fertilizer added in Run 3 and 4 substantially changed the plant development curve from sigmoidal to exponential. The architecture for the 120% moisture level was appeared extremely different from the other moisture levels (Table 2, Figure 8). The  Figure 7e shows the EMD score overall. Figure  7f, Figure 7g, and Figure 7h represent the top 20% value changes between 80% and the three other mositure levels; red spatial bins are unique to 80%, green spatial bins are unique to 60%, 100% or 120%, and yellow spatial bins are common to both water levels.
roots barely grew beyond the flooded layer of soil, and tended to spread out laterally to a greater extent than in the other soil moisture conditions. The roots for 120% also tended to be more dense right above the flooded layer of soil, as evidenced by higher color scale value in Figure 8a.

Avoidance vs Tolerance in waterlogged soil conditions
When looking at the EMD results, one interesting observation is that the three features that were considered (i.e., % root pixels, fractal dimension and lacunarity) had very similar discriminating power between cultivars and between moisture levels (Table 1, Table 2) though the ANOVA/Tukey HSD test showed that % root pixels tended to be a better feature than fractal dimension and lacunarity. Both fractal dimension and lacunarity have been used successfully in RSA analysis (Li et al, 2020;Walk et al, 2004). However, the main differences between these previous studies and the current one is that in this study these features were not calculated across the entire image but for each spatial bin. 120% and the three other moisture levels; red spatial bins are unique to 120%, green spatial bins are unique to 60%, 80% or 100%, and yellow spatial bins are common to both water levels.
We also calculated the fractal dimension and lacunarity values for whole images, but we observed limited differences among treatments for the CH index (e.g., S32 = 0.06 and S38 = 0.03 for Runs 3 and 4 with fractal dimension). By using these texture features to describe the individual regions of an image, we improved the RSA characterization and comparisons for each treatment. Results from Tables 1 and 2 indicate that for future similar work, % root pixels, an easily calculated parameter, should be preferentially used rather than fractal dimension and lacunarity.
More specifically, the results for Runs 1-2 and Runs 3-4 (Table 1 and 2) indicate that without fertilizer, cultivar is the dominant driver of root architecture, and that water only plays a minimal role. This can be seen in the highest CH index score on average across all cultivars compared with the averages for water level (Table 1). S32 was the cultivar with the most distinct architecture, and through Figure 4, we are able to pinpoint what made this architecture different from the others such as having most changes in the deeper section of the root architecture through a less concentrated root mass along the right of the tap root (S35), the most variation in the lateral spread in the middle part of the root system (S38), and the most difference in the upper part of the root system by a more spread out root distribution towards the top of the soil surface (S40).
The situation is reversed when fertilizer is added, and water level becomes the factor with the biggest impact on root architecture, as evidenced by the much higher CH index value for moisture level. This high score is actually driven by a single treatment, 120% moisture level, which had a very different architecture compared to the other soil moisture treatment (Figure 8). The roots seemed to not be able to grow much into the flooded layer of soil, and tended to grow laterally more so than in the other treatments. This very distinct RSA in flooded conditions when fertility is adequate in Run 3 and 4 contrasts strongly with that of Run 1 and 2. Sesame is known to be highly sensitive to waterlogged soil (Sarkar et al, 2016), and can activate a variety of morpho-anatomical and physiological responses to cope with flooding stress (Wei et al, 2013). Here, we can highlight two different morpho-anatomical strategies employed by the crop depending on the soil fertility status. Without adequate nutrient availability, the crop seems to opt for a tolerance strategy to flooding, developing roots into the saturated soil (Figure 7). We can assume that this apparent tolerance strategy is accompanied but other changes not captured in this study such as the formation of aerenchyma and changes in the enzymatic activities (Wei et al, 2013). When there are enough nutrients in the soil however, the plants seemed to adopt an avoidance strategy, and did not grow roots down into the waterlogged soil, but tended to proliferate laterally ( Figure 8). These observations being constant across cultivars, we can thus conclude that fertility conditions the response of sesame early RSA to flooding stress.

Validity of the EMD Method
Although of interest in estimating early root vigor and biomass, it is clear here that only considering TRL does not provide any insight on RSA. Increasing our knowledge and understanding of root architectural development in response to environmental stresses is of capital importance for global agricultural production, and has been defined as the pillar of the Second Green Revolution (Lynch, 2007). Many methods are now available to observe RSA, yet some of the more advanced techniques such as CT imaging remain expensive and impractical, underlining the need for improved trans-disciplinary phenotyping approaches that can capture and quantify the complexity of RSA Wu et al, 2018;Zhu et al, 2011).
Our proposed STARSEED method could be the answer to some of these issues. It provides a fast, reproducible, quantitative, and spatially explicit characterization of RSA, also allowing for comparisons between root systems. A valuable aspect of the analysis is the introduction of the EMD class score, which summarizes all of the differences between two RSAs down to a single number which can be compared between these RSAs. The EMD scores also provide an interpretable measure that quantifies the visual observations between the root structures. EMD, combined with the CH index, can thus be used to very quickly know which treatment lead to the most distinct RSA, and the degree to which these RSAs relate to one another globally.
The actual spatial EMD result (i.e., Figures 4 through 8) can then be looked at to further elucidate specific spatial differences between RSAs. To the best of the authors' knowledge, this is the first time EMD has been used to characterize RSA. The method has been thoroughly validated and used in varied fields of science, including fluid mechanics (Benamou and Brenier, 2000), linguistics (Kusner et al, 2015), or, more related to the present study, image classification and comparison (Zhang et al, 2020). The CH index used in STARSEED serves as evaluation step for the method to not only quantify what we observe, but also to ensure the approach captures distinct features to effectively group root architectures based on a shared characteristic (e.g., cultivar, moisture treatment).
One of the most critical strengths of STARSEED and the way it is performed here is that it allows for the generation of an "average" root architecture visual representation as shown in Figures 4 through 8a through d. Indeed, by taking the overlap of all images within one treatment, partitioning the image into local regions, and calculating a feature value for each spatial bin, we can generate a meaningful and visually explicit "average" architecture that reflects differences between genotypes or environmental conditions. Recent studies have attempted to either develop new techniques to directly measure RSA (Armengaud et al, 2009), find a way to accurately represent average RSAs corresponding to specific growing conditions (Shahzad et al, 2018), and many more studies have developed or refined root architecture development models (Tron et al, 2015;Postma et al, 2017;Zhao et al, 2017;Pagès et al, 2020). One study in particular created similar 2D heat maps of "root frequency" observed on the four transparent surfaces of rhizoboxes but did not perform a spatially explicit analysis to the level of STARSEED (Jørgensen et al, 2014).

Going beyond the 2D images of rhizoboxes
Given the rhizobox set up of this study and the use of flatbed scanners for root images acquisition, the observations are confined in 2D when root systems always exist in 3D spaces. Imaging and characterizing a whole RSA in 3D in situ or even in vitro is usually low throughput and expensive if not nearly impossible. Most scientists thus resort to using models in lieu of direct observations, as previously mentioned. One of the main difficulty of using these 3D RSA models is correct parameterization, which is crucial to the accuracy and validity of a model's prediction . Recent work has shown that 2D measurements of root systems could be used to adequately inform the parameters for 3D models (in wheat ). It is thus a reasonable hypothesis to say that the current EMD results could be further refined and subsequently used to generate and interpret such 3D models. There is potential to go even further and directly use EMD on 3D data (Salti et al, 2014), although the issue of acquiring such data still remains.
Combining the average EMD visual maps with new root architecture quantification techniques (Wu and Guo, 2014) could really boost our understanding of RSA plasticity. This would allow for both global (i.e., EMD scores) and local (i.e., magnitude and direction of changes) comparisons between genotypes and environmental conditions. Additionally, combining our understanding of RSA plasticity derived from the EMD visual maps along with molecular and genetic work would most probably be of critical importance to breeders.

Conclusion
In this paper, we presented STARSEED, an approach to characterize root architecture from 2D images. Qualitative and quantitative analysis demonstrate the effectiveness of the proposed method. STARSEED successfully incorporated spatial and texture information to describe root architectures at both global and local contexts. The method is explainable and provides clear connections to biological aspects of each root image. STARSEED also allows for aggregating individual architectures to assess and average response for each environmental condition and genotype. Future work includes clustering based on root architecture (i.e., using the pairwise EMD matrix for relational clustering), applying our general framework over time to characterize how RSA develops in various conditions, along with scaling up to 3D characterization.

Declarations
Funding This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1842473.

Conflict of interest/Competing interests
The authors have no conflicts of interest to declare that are relevant to the content of this article.

Ethics approval Not applicable
Consent to participate Not applicable Consent for publication Not applicable experiments, and performed analysis. Joshua Peeples, Weihuang Xu, and Romain Gloagen wrote the manuscript. Diane Rowland, Alina Zare, and Zachary Brym reviewed and edited the manuscript. All authors read and approved the final manuscript.