3D U-Net Segmentation Improves Root System Reconstruction from 3D MRI Images in Automated and Manual Virtual Reality Work Flows

Magnetic resonance imaging (MRI) is used to image root systems grown in opaque soil. However, reconstruction of root system architecture (RSA) from 3-dimensional (3D) MRI images is challenging. Low resolution and poor contrast-to-noise ratios (CNRs) hinder automated reconstruction. Hence, manual reconstruction is still widely used. Here, we evaluate a novel 2-step work flow for automated RSA reconstruction. In the first step, a 3D U-Net segments MRI images into root and soil in super-resolution. In the second step, an automated tracing algorithm reconstructs the root systems from the segmented images. We evaluated the merits of both steps for an MRI dataset of 8 lupine root systems, by comparing the automated reconstructions to manual reconstructions of unaltered and segmented MRI images derived with a novel virtual reality system. We found that the U-Net segmentation offers profound benefits in manual reconstruction: reconstruction speed was doubled (+97%) for images with low CNR and increased by 27% for images with high CNR. Reconstructed root lengths were increased by 20% and 3%, respectively. Therefore, we propose to use U-Net segmentation as a principal image preprocessing step in manual work flows. The root length derived by the tracing algorithm was lower than in both manual reconstruction methods, but segmentation allowed automated processing of otherwise not readily usable MRI images. Nonetheless, model-based functional root traits revealed similar hydraulic behavior of automated and manual reconstructions. Future studies will aim to establish a hybrid work flow that utilizes automated reconstructions as scaffolds that can be manually corrected.

(sequence, resolution). The undistorted center slice of the scanned MRI phantom stack was used as reference slice. Using affine transformation, each horizontal slice up and down the center slice 731 was registered to the reference slice and dewarped. The generated transformation matrices were 732 subsequently used to de-warp the MRI lupine scans slice-by-slice. As mentioned above, each column 733 was scanned in three sections with an overlap between 5 and 10 mm. These scans were subsequently 734 stitched together using the following approach: first, we normalized the voxel intensities of the 735 different sections using the intensities of the marker tube as a reference. Next, we determined 736 the two overlapping slices of two neighboring sections using the criteria of maximal correlation. 737 The overlapping slice of the lower section was then mapped to the overlapping slice of the upper 738 section using an affine transformation. All remaining slices of the lower section were then equally 739 registered with the slice above using the same transformation matrix. Finally, the three sections were 740 concatenated. Subsequently, the connected structure of the marker tube was masked and excluded 741 from the MRI scans. The images pre-processed in this way were the starting point for the three root To allow a more complete extraction of roots from imperfect data (i.e. data with gaps), Dijkstra's 746 shortest path algorithm [Dijkstra1959] is modified. The path-cost threshold is momentarily ignored 747 to allow the exploration of high-cost voxels within the maximum gap length. Allowing such traversals 748 requires a new cost-map, as using the initial cost-map (Fig. 3b) would result in the emergence of 749 multiple high-cost paths from disconnected segments of the same root. Hence, the initial cost-map is 750 adapted when using the gap closing option: voxels below the given minimum intensity threshold are 751 not excluded from the cost map, but their cost is increased by a factor of ten (Fig. 3c). Subsequently, 752 all voxels above half of the maximum voxel cost are considered as potential gaps. This preserves low-753 intensity information while simultaneously enhancing the contrast between gap and no-gap voxels.

754
Penalizing these low-intensity voxels also ensures that the shortest path algorithm explores no-gap

766
A radius estimate map is extracted from the largest connected component (Fig. 3d). This is used as basis for another cost map used for the Dijkstra's shortest path algorithm [Dijkstra1959]. Larger is calculated based on the container volume V c : For the MRI experiment, V c was 455 cm 3 . The half-mean-distance between roots, HMD (cm), is also 786 based on the root length per container volume and approximated following the classical approach 787 proposed by Newman [Newman]: We calculate the mean radius of a root system as where r i is the radius of a root segment, l i is the length of a root segment, and n is the total number 790 of root segments of a given root system.

791
The total number of root tips per tracing, as well as the number of lateral root tips is given as 792 a topology measure. Here, we specify lateral tips up to the highest order observed in the manual 793 reconstructions (3 rd order lateral roots). Although lateral roots above 3 rd order are recorded in the 794 automatic tracings, we refrain from stating them in the root measures results section for improved 795 readability. Laterals above 3 rd order are qualitatively accessible in the visual comparisons. 796 We compute the equivalent conductance of the root system, K rs (cm 2 d −1 ), according to Meunier where T act (cm 3 day −1 ) is the actual transpiration rate, Ψ sr (cm) is the mean soil water potential 799 at the soil-root interfaces, and Ψ collar (cm) is the water potential at the root collar. We chose a 800 scenario where Ψ collar is set to -15000 cm, applied as Dirichlet boundary condition, and Ψ sr is set 801 to -500 cm at the soil surface, while assuming a hydrostatic equilibrium in the soil domain. K rs 802 reflects the ability of a root system to take up a certain water volume under a given water potential

807
The standard uptake fraction of a root segment, SUF i (-), indicates its relative contribution to 808 T act and is calculated via where J r is the radial water flux into a root segment (cm 3 day −1 ). To obtain an aggregated (scalar) 810 metric, we determine the mean depth of standard root water uptake, zSUF (cm), by where z i (cm) is the respective depth of a root segment.

812
Finally, we investigate in how far the root hydraulic properties affect the impact of differences in 813 tracings of the three reconstruction methods on root system function metrics. In the first scenario, 814 called constant scenario, we calculate the metrics by applying the same fixed axial and radial con-815 ductivities to all roots. Based on Zarebanadkouki et al.
[Zarebanadkouki], we apply a fixed axial 816 conductivity of k z = 4.32e −2 (cm 3 d −1 ) and a fixed radial conductivity of k r = 1.73e −4 d −1 to all 817 roots. For the second scenario, called variable scenario, we apply order and age-dependent root hy-818 draulic properties. As we did not reconstruct time-series, we linearly interpolate root segment ages 819 as a function of root length while assuming a daily growth rate of 1 cm. A resulting age distribution 820 for an exemplary 14-day old root system is shown in Fig. S1 Figure S1: Age-distribution used to assign age-dependent root hydraulic conductivities to manual tracings M (left), manual tracings based on segmented images M+ (middle), and automated tracings A (right). Shown is an exemplary 14-day old Lupinus albus root system of the MRI dataset. Figure S2: Age-dependent root hydraulic conductivities applied in the variable simulation scenario.
For the M and M+ reconstructions, we also report the respective reconstruction speed, v r (cm 829 root min −1 ), which is calculated as where t r (min) is the total time required to trace a root system. Division by the respective RL is 831 done to mitigate differences in total root length on t a between M and M+ reconstructions.

832
To give an estimate of the quality of the MRI images, we calculate an exemplary contrast-to-