Nondestructive automated workflow for analyzing diverse leaf morphologies using computed tomography

Automated plant analysis methods can inform basic and applied plant research. Computed tomography (CT) has been used to image plants in canopies and individual plants; however this method is underutilized for dynamic analysis of plant growth. In this work we present a workflow and associated algorithm to nondestructively extract leaf area from 120 individual CT scans for plants with flat [soybean, Glycine max (L.) Merr.], textured (tomato, Solanum lycopersicum L.), and grassy (wheat, Triticum aestivum L.) leaves. Under low water conditions, we see significant changes in leaf area depending on leaf type. This work enables future automated phenotyping using CT scanning of whole plants.


INTRODUCTION
Automated extraction of quantitative plant metrics represents a bottleneck to understanding genetic and environmental influences on plant growth and yield (Furbank & Tester, 2011). Experimental and analytical methods that provide high-resolution plant structural data can enhance our understanding of plant biology. X-ray computed tomography (CT) provides unique benefits in that it does not suffer from occlusion caused by optically opaque plant tissues, such as is observed from using light-based approaches such as laser scanning or three-dimensional (3D) reconstructions using multiple camera images. Automatic methods to utilize CT reconstructions are required to harness the full potential of this technology.
There are multiple technologies and diverse imaging approaches for phenotyping. Field-based phenotyping methods include stand-off imaging approaches, which use hyperspectral sensors and algorithms to compare reflection of specific wavelengths and link to physiological responses (Meacham-Hensold et al., 2019). Threedimensional reconstructions have previously been used to study tree canopies (Ni, Burks, & Lee, 2016). Similar technologies include laser scanning or light detection and ranging (LiDAR) technologies (Gage et al., 2019), and point cloud analysis methods to extract leaf area and leaf metrics (Ziamtsov & Navlakha, 2019) or for fruit phenotyping (Feldmann, Tabb, & Knapp, 2019) have recently been developed. Finally, hand-held leaf scanners can be used for non-destructive leaf area estimation, although these require manual access to individual leaves (Hammami et al., 2017).
Emerging technologies for high-throughput phenotyping include chlorophyll fluorescence and imaging integrated into plant growth facilities. Magnetic resonance imaging (MRI) has been used to acquire high-precision plant morphology measures, such as imaging in the root zone (Tripodi et al., 2018). Positron emission tomography has also been applied to plants using radioactive tracers (Garbout et al., 2012). Interestingly, non-imaging approaches such as volatile sampling are being developed as a method to characterize plant phenotypes, and could be adapted for phenotyping greenhouses (Niederbacher, Winkler, & Schnitzler, 2015). In imaging approaches, plant segmentation can be a problem and algorithms to identify leaves and plants structures are actively being developed (Gehan et al., 2017). Automated imaging approaches and analytical techniques that can reliably segment plant tissues and provide quantitative metrics will enable continued improvements in high-throughput plant phenotyping.
Computed tomography scanning operates by subjecting samples to X-ray radiation and reconstructing images from the collected X-rays. One benefit of this approach is the ability to make repeated scans on the same sample, allowing for dynamic measurements (Lafond, Han, & Dutilleul, 2015). This approach can be implemented using medical-grade scanners with resolutions of ∼1 mm (Dutilleul, Lontoc-Roy, & Prasher, 2005), micro-CT with resolution down to 2 μm, and emerging nano-CT systems that approach 100 nm resolution. Diverse plants have been studied using available X-ray CT systems. Sorghum [Sorghum bicolor (L.) Moench] stem properties (length, diameter, pith) were obtained using medical-grade scanners as a method for high-throughput quantification of phenotypes (Gomez, Carvalho, Shi, Muliana, & Rooney, 2018), and starch reserves in grapevine (Vitis spp.) were inferred from micro-CT scans (Earles et al., 2018). Micro-CT was also used to study wheat (Triticum aestivum L.) grain traits (Hughes et al., 2017) and mesophyll structure in maize (Zea mays L.) (Dorca-Fornell et al., 2013;Du et al., 2017). Interestingly, root properties can be inferred by observing soil characteristics in growing plant systems (Keyes et al., 2016), micro-CT scanning has been used to investigate fluid dynamics in the rhizosphere (Yang, Varga, Liu, & Scheibe, 2017), and soil compaction in response to root growth has been quantified using X-ray computed tomography methods (Tracy et al., 2012). Furthermore, micro-CT scanning has been used to study the effects of multi-cropping on root systems (Mairhofer, Sturrock, Bennett, Mooney, & Pridmore, 2015).
In this study, we implemented a workflow to CT scan whole plants in order to quantify dynamic leaf area and demonstrate the workflow under high and low water treatments for plants with diverse leaf morphologies. We developed an associated automated image processing algorithm that requires little user input to extract leaf areas. We measure and report leaf areas throughout a water stress for soybean [Glycine max (L.) Merr.], tomato (Solanum lycopersicum L.), and wheat plants.

Core Ideas
• We performed repeated, non-destructive CT scanning of whole plants with a low water treatment. • We developed an algorithm to automate calculation of leaf area. • This method could be used in high-throughput phenotyping systems for high-resolution plant structural analyses.

Plant scanning and leaf area determination
Soybean, tomato, and wheat plants were selected to provide examples of flat, textured, and grassy leaves. A total of 10 plants per species across two watering regimes were scanned weekly on a clinical whole body CT scanner, specifically the Lightspeed VCT by GE. Using MATLAB (Mathworks), we developed a process to analyze CT scans to extract leaf area, volume, and projected area. The algorithm followed the general workflow: (a) preprocessing, (b) segmentation, and (c) calculation. In preprocessing, data were converted from digital imaging and communications in medicine (DICOM) to Neuroimaging Informatics Technology Initiative (NIfTI) format for 3D processing, manual separation of individual plants that had been scanned in a single run, size standardization of voxels between images, and orientation adjustment of plants. Segmentation involves two substeps: (a) thresholding and (b) feature engineering. The goal of the first substep of thresholding is to identify by image intensity connected regions that may be identified as leaves. To accomplish this, a segmentation was created using an empirically determined threshold of 80 (i.e., -944 Hounsfield units, HU). In addition, a morphologically dilated segmentation created using an empirically determined threshold of 400 (-624 HU) was also subtracted to remove bright objects like the table, pot, and stems. The end result is a leaf segmentation with some noise (Figure 1b). The goal of the second step of feature engineering is to determine features that would be useful in distinguishing between leaf regions and noise regions. We used two features to create our final segmentation. The first feature was volume as we expected leaves to be larger than a certain size (125 voxels or 125 mm 3 ). The second feature was termed dilated volume ratio. To calculate this feature, each region was dilated under a thresholded mask and divided the F I G U R E 1 Algorithm workflow. (a) pre-processing workflow and example preprocessed image. (b) Segmentation workflow and example segmented image. (c) Calculation workflow and example rendered image resulting dilated volume over the original volume. Leaves were expected to have small ratios (<1.1) since they are connected to stems of small sizes. The end result is a segmentation with all leaf voxels of the plant labeled. The final step of the pipeline was calculation, where three metrics were calculated: (a) volume, (b) surface area, and (c) projected area. Volume and surface area were calculated using MATLAB's regionprops3 function, and projected area was calculated by counting the number of pixels included in a top-down orthographic two-dimensional (2D) projection of the 3D segmentation. This process enabled us to identify regions of scans that belong to leaves while omitting stem, pot, and other confounding factors (e.g., image artefacts). An example illustrating specific algorithm steps is shown for a soybean plant in Figure 1.

Algorithm performance on different leaf types
Our algorithm was next implemented on the final scan of the study for the different plant types (Figure 2), which allows us to compare the extracted leaf area to the manual method of removing leaves, flattening, and scanning. A tunable parameter in the algorithm is the threshold used in the thresholding substep of segmentation, which affects the size of the leaf regions segmented as well as the number of false positive regions. We observe that a threshold setting of -944 HU minimized the difference between the CT scanning method and manual method for soy plants, but this parameter was increased to a value of -934 HU to capture area of wheat leaves (Supplemental Figures S1; Figure 2c, 2f, 2i). Optimization of the threshold would be ideal for additional plant species.
Final leaf area was decreased for both soy and wheat plants in response to low water condition, whereas tomato plants showed increased leaf area at the final time point. This response likely indicates the high water condition represented the stress treatment for tomato.

Dynamic leaf area
Leaf area and projected area were extracted using the automated workflow and followed the same trends as height and leaf number (Figure 3, Supplemental Figure S2). We observed extreme dehydration in two of five low water soybean plants during the 4th week of treatment. Tomato plants under low water conditions had increased leaf area, height, and leaf number as compared to controls, further supporting the conclusion that tomato plants grew more readily under the low water treatment. In both soy and tomato, leaf area showed similar trends between low water and control groups as other nondestructive measures of height and leaf number. In wheat plants, height showed little variation over time, whereas leaf area and leaf number showed similar separations.

DISCUSSION
We developed and implemented a CT scanning workflow with an associated algorithm that nondiscriminately extracts leaf area from plants with very different leaf morphologies. We present an image analysis algorithm to quantify leaf area, and demonstrate the use of the algorithm to quantify leaf area under high and low water conditions for flat, textured, and grassy leaved plants.

F I G U R E 2 Comparison of CT with traditional leaf area for the final time point. Pictures of grown plants soy (a), tomato (b), and wheat (c).
Algorithm segmentation of leaves for soy (d), tomato (e), and wheat (f). Leaf area by CT scanning (gray bars) or traditional destructive leaf area measure (white bars) for soy (g), tomato (h), and wheat (i) This approach is an example of a technology that could be readily integrated into high-throughput phenotyping greenhouses. Single-pass scans coupled to automatic image processing algorithms results in detailed data that can be analyzed offline. However, we note that plants not tested in this study may provide unique challenges that require further customization of such algorithms. Leaf area is a common metric that is sensitive to plant stress (Banerjee, Krishnan, & Mridha, 2017;Gómez-del-Campo, Ruiz, & Lissarrague, 2002). Methods to evaluate leaf area response to stress include canopy imaging, point cloud reconstruction from RGB (Zhou et al., 2019) or LiDAR imagery (Gage et al., 2019), handheld leaf scanning technologies (Hammami et al., 2017) and destructive methods where leaves are harvested from plants and analyzed by scanning (Timm et al., 2018). The method we present here provides high-resolution data for individual plants that can be used to advance our understanding of plant stress response. Extensions of this method include higher resolution CT scanning that allows analysis of internal plant structures (Wang, Verboven, & Nicolai, 2017), analytical methods to analyze plant growth and branching patterns (Bessonov, Morozova, & Volpert, 2008), and optimization to image root structures (Subramanian, Han, Dutilleul, & Smith, 2015), all factors that are important for understanding plant stress responses.
One practical complication we encountered was the slow reconstruction time for scanned data. Although active scan time of any individual plant was less than a minute, reconstructions could take up to 10 min, potentially causing a backlog in processing. Data management is thus a consideration for systems in which scans can be rapidly collected, and efficient data management would mean scan-to-answer time on the order of minutes for any single plant. A major advantage of the approach is the detailed data that can be collected and stored for ongoing analysis of phenotypes of interest.
Further practical considerations for our study are the orientation and multiplexing of plants. Initially, plants were scanned individually and vertically passed upright through the scan area. To increase throughput, we found that larger scans with multiple plants reduced overall scan time. However, as plants grew through our treatments, they eventually outgrew the vertical limits of the scan area. Our approach to resolve this issue was to lay plants horizontally on their side. For implementation in phenotyping applications plants could be dropped vertically through a scanning ring. Plant type and width would thus dictate the exact implementation of this approach. Further improvements of this method include development into additional leaf metrics, such as leaf count, shape, size, and plant metrics such as projected area, and branching structure. In this work, leaf thickness approached the resolution limit of the scanner. Recent methods have demonstrated utility of 3D reconstructions by combining multiple images from different perspectives relative to the plant (Gibbs et al., 2018), and the translatability of algorithm presented in this work to 3D point clouds is of interest.
In this work, we presented a scanning approach and associated algorithm to extract leaf metrics dynamically for soybean, tomato, and wheat plants. Our goal was to develop a single algorithm that performs on multiple leaf types and to show how leaf metrics change with stress, in this case a water stress. The methods presented here represent how CT scanning could be implemented for highthroughput phenotyping. This approach can be implemented into current or new phenotyping greenhouses to generate detailed data on plant growth metrics.

Plant growth and treatments
Soybean (Williams 82 PI518671), tomato (Heinz 1706 PI 652934), and wheat (Crystal PI 351960) seeds were obtained from the U.S. National Plant Germplasm System (NPGS-GRIN). Two to three seeds were sown in 7.6-cm (3-inch) pots filled with Fafard 2 potting mix (BFG Supply Co.) with 12-h days 31 • C/22 • C day/night temperatures at 400 pars of light, and then plants were thinned to 1 plant per pot for a total of 10 plants per species. Plants were maintained by watering three times a week, once with 0.66 g L -1 (100 ppm) Jacks 15-16-17 fertilizer (BFG Supply Co.) to saturation. After 2 wk of growth, plants were subjected to low or high water treatments (n = 5 per condition). The high water condition was maintained by watering 3× weekly, whereas the low water condition was implemented by watering 1× per week. Plant height and leaf number were manually measured weekly until Week 5 post-germination when tissues were harvested for destructive leaf area. At harvest, leaves were clipped from plants and photographed for leaf area calculation. Leaf area was quantified using thresholding and area measurements using ImageJ (Schneider, Rasband, & Eliceiri, 2012). Dry weights were measured by separating the leaf, stem, and belowground biomass at harvest and incubating tissues at 70 • C for 14 d before weighting.

CT scanning
All plants were scanned weekly, from Week 2 to 5 after germination, beginning concurrently with water level treatments using a clinical whole body CT scanner, specifically the GE Light Speed VCT 64 (GE Healthcare). Each slice thickness is 0.625 mm and the X-ray detector is 40 cm wide. Raw data from scans is available at CyVerse (https://doi.org/10.25739/6atr-5773). In Week 2, plants were scanned vertically through the CT scanner as independent scans. In Week 3, plants were scanned in groups of two to four, then later plants were first manually segmented prior to algorithm workflow. In Weeks 4-5, plants were scanned horizontally in pairs to accommodate for plant height exceeding vertical dimensions of the scan area.

Image analysis for leaf area
Our image analysis pipeline was developed primarily in MATLAB (Mathworks) to process the CT scans of these plants and consists of several steps: (a) preprocessing, (b) segmentation, and (c) calculation. Preprocessing involves four substeps: (a) conversion, (b) separation, (c) resizing, and (d) reorientation. The first substep of conversion involved converting the CT scan from DICOM to NIfTI file format. The DICOM (.dcm) format is standard for acquiring medical images like CT images, but given that it saves a 3D scan as multiple 2D files, it is not as convenient for 3D image processing as the NIfTI (.nii) format that saves a 3D scan as a single 3D file. As such, a batch script was created in MATLAB calling Convert3D functions (https://sourceforge.net/projects/c3d/) to automatically perform conversion from DICOM to NIfTI. The second substep of separation involved a procedure to manually separate individual plants in scans with multiple plants.
We manually defined cut-planes using ImageJ (https://fiji. sc/) to separate plants in these scans into their individual files. The next substep of resizing involved standardizing the voxel resolution and image dimensions across all scans. Because the voxel resolutions and image dimensions of the original raw scans were inconsistent, a MAT-LAB algorithm was written to make every scan isotropic at 1 × 1 × 1 mm using linear interpolation and consistently sized at 512 × 512 × 512 pixels with zero (i.e., -1,024 HU) padding. The last preprocessing substep of reorienting involved changing the orientation of some of the scans. In a fully automated phenotyping system, we imagine that every plant would be scanned in the same orientation.
To address the inconsistency in orientation during image processing, we manually identified the scans that needed reorientation and applied it to the scans using Convert3D. Segmentation involves two substeps: (a) thresholding and (b) feature engineering (Figure 1). The goal of the first substep of thresholding is to identify by image intensity connected regions that may be identified as leaves. To accomplish this, a segmentation was created using a threshold of 80 (i.e., -944 HU), but a dilated segmentation created using a threshold of 400 (-624 HU) was also subtracted to remove bright objects like the table, pot, and stems. We used two features to create our final segmentation. The first feature was volume as we expected leaves to be larger than a certain size (125 voxels or 125 mm 3 ). The second feature was termed dilated volume ratio. To calculate this feature, each region was dilated under a thresholded mask and divided the resulting dilated volume over the original volume. Leaves were expected to have small ratios (<1.1) since they are connected to stems of small sizes. The end result is a segmentation with all leaf voxels of the plant labeled.
The last step of the pipeline was calculation, where three metrics were calculated: (a) volume, (b) surface area, and (c) projected area. Volume and surface area were calculated using MATLAB's regionprops3 function, and projected area was calculated by counting the number of pixels included in a top-down orthographic 2D projection of the 3D segmentation. The source code for this approach is available at GitHub (https://github.com/nathanaelkuo/ CT_leaf_analysis).