High-throughput volumetric reconstruction for 3D wheat plant architecture studies

Wei Fang*, Hui Feng*, Wanneng Yang*†‡, Lingfeng Duan, Guoxing Chen, Lizhong Xiong and Qian Liu*¶ *Britton Chance Center for Biomedical Photonics Wuhan National Laboratory for Optoelectronics Huazhong University of Science and Technology 1037 Luoyu Rd. Wuhan 430074, P. R. China National Key Laboratory of Crop Genetic Improvement and National Center of Plant Gene Research Huazhong Agricultural University Wuhan 430070, P. R. China College of Engineering Huazhong Agricultural University Wuhan 430070, P. R. China MOA Key Laboratory of Crop Ecophysiology and Farming System in the Middle Reaches of the Yangtze River Huazhong Agricultural University Wuhan 430070, P. R. China ¶qianliu@mail.hust.edu.cn


Introduction
In recent decades, it has been estimated that we will need to double the quantity of food produced to meet the demands of a rapidly growing population. 1 Increasing crop yields has become a major challenge for modern agriculture. To meet this challenge, more productive crop varieties must be bred based on the identi¯cation of new plant characteristics. Plant architecture (PA) has been proven to significantly a®ect crop yield, where the yield is a®ected primarily by the PA components of plant height, 2 number of tillers, 3 tiller diameter and tiller angle. 4,5 The characteristics of the ideal plant architecture (IPA) include small numbers of tillers, few unproductive tillers and thick and sturdy stems. 6 To measure these traits in multiple tiller crops (e.g., rice, wheat and barley), the usual method consists of many laborious manual measurements. Clearly, considering the rapid development of functional genomics and gene technologies over the past decade, it is impractical to perform manual phenotypic analysis. Consequently, it is imperative to develop better measurement technologies to automate the tedious task of collecting large amounts of phenotypic data with a high-throughput for agricultural crop improvement.
Many approaches have been developed to facilitate plant phenotypic research. Based on two-dimensional (2D) image processing, visible light imaging technology has been widely used in automated plant phenotypic analysis, 7 such as leaf phenotypic analysis (to determine the leaf shape, size, width, length, area and perimeter), [8][9][10][11] yieldrelated trait analysis (to determine the number of tillers, the panicle length and the grain size, length, width and thickness), [12][13][14] root system analysis (to determine the average radius, area, horizontal width and length distribution) [15][16][17] and general plant analysis and measurement (to determine the plant height, plant width, center of gravity, projected area and biovolume). 18,19 However, the information provided by these 2D approaches is very limited; many features cannot be extracted using 2D techniques, such as the angle, thickness, bending and orientation. Therefore, three-dimensional (3D) approaches have been proposed to overcome the loss of spatial information (thickness, volume and angles) in 2D image processing. The 3D plant model could be represented as point cloud, surface or volume. And it could be obtained directly from special devices or reconstructed from plant images. Point cloud data of plant could be obtained directly from devices which could measure the depth of the scene, such as laser scanner, 20 light detection and ranging (LIDAR) laser, 21 and time of°ight (TOF) scanners. 22 The equipment used in the system, however, could be expensive and often places restrictions on the environments in which it works.
Based on low-cost camera and the techniques developed in computer vision¯eld, image based approaches for plant modeling attempt to extract information from images and reconstruct 3D model of the plant. This kind of method will be more e±cient with the continued increase of computing power, which means that the method has great potential in plant phenotyping. Existing 3D modeling approaches from images could be classied as surface-based and volume-based. In surface modeling approaches, the reconstructed plants are described as a set of planar sections. As the measurement of leaf area could be easily made in surface model, surface model is suitable for plant canopy photosynthesis e±ciency and leaf morphology research. [23][24][25] While, volume-based modeling approach has been always developed to measure those traits of solid or thick organ with complex structure, such as root system (RootReader 3D software) 26 or tillers. That's because, in volumetric model, plants are represented as many voxels, which makes the measurements of height and volume easy.
Image-based 3D reconstruction requires a set of images captured from di®erent viewpoints. And the calibration in the imaging system could be implemented in various ways. Clark et al. 26 developed an axis of rotation (AOR) calibration technique to determine the orientation of the AOR in relation to the camera. Bellasio et al. 24 used coded light illumination to obtain the 3D information. However, a complicated calibration technique would make the imaging system unaccessible. Therefore, there is a need for e±cient, accessible techniques for acquiring the 3D PA of wheat.
In our previous work, we developed a highthroughput rice phenotyping facility for rice plant parameter extraction using 2D image processing technology. 27 In this study, we propose a method based on automatic volumetric reconstruction to nondestructively acquire 3D PA and the corresponding phenotypic traits (including the plant fresh weight, plant height, number of tillers, stem diameter and tiller angle) with high e±ciency (20 s per plant) for wheat.

Material and Methods
This work included plant volumetric reconstruction from multiple images, plant model processing and phenotypic trait extraction and analysis. The approach was based on a computer vision system consisting of a digital image acquisition system and 3D reconstruction software. The process°ow is shown in Fig. 1. This section will introduce the digital image acquisition system.

Digital image acquisition system
We constructed volumetric models of plants for phenotypic trait measurement. The method used was based on taking a set of digital photographs of the plants from 20 angles. To acquire the photographs, an image acquisition system was developed.
The system consisted of a color camera (AVT Stingray F-504B/F-504C, Allied Vision Technologies Corporation, 2452 (H) Â 2056 (V) resolution) with a 12 mm¯xed focal lens (M1214-MP, Computar), a rotatable platform, a servomotor controller (MBDDT2210, Panasonic Corporation, China), several LED lamps and a computer workstation (HP xw6400, Hewlett-Packard Development Company, USA). As shown in Fig. 2(a), a pot-grown wheat plant with a calibration pattern was placed on a rotatable platform driven by the servomotor. The computer commanded the servomotor controller, precisely driving the platform to a given angle. One image was obtained for every 18 of plant rotation ( Fig. 2(c)). All images were saved in PNG format for subsequent analysis. The system was integrated into an automated high-throughput system for measuring rice tillers (the H-SMART system) developed in our previous work. 12

Digital image pre-processing
To reconstruct the 3D shape from a silhouette, each pixel in the images must be classi¯ed as plant or background. To enhance the image plant region extraction, the excess green (ExG) vegetation index 28 was used in the image processing. The main steps of image processing were as follows: (1) The red, green and blue (RGB) color planes were extracted. (2) The ExG value used to distinguish the pixels was calculated for each pixel using the equation below, and the image was transformed into a grayscale image.
where r, g and b indicate the pixel values of the red, green and blue channels, respectively.
(3) Several operations were performed on the grayscale image:¯rstly, the image segmentation was implemented with a¯xed threshold and the plant region was extracted. Then, the noise spots in the binary image were eliminated with a prede¯ned area threshold, which was 10 pixels in this work. Last, morphological closing operation was applied to the binary image to obtain an integrated silhouette (Fig. 3).

Camera calibration
In a pinhole camera model, each 3D point in the world coordinate system can be projected onto the image plane using the perspective transformation in Eq. (2). 29 where (X; Y ; ZÞ are the coordinates of a 3D point in the world coordinate system; (u; vÞ are the coordinates of the projection point in the image plane; (c x ; c y Þ are the coordinates of the principal point of the camera in the image plane; (f x ; f y Þ are the focal lengths expressed in pixel units; (R; T Þ are the rotation and translation matrices, also called the extrinsic parameter matrices. They are used to describe the camera position in the world coordinate system. (c x ; c y Þ and (f x ; f y Þ are called the camera intrinsic parameters, and (R; T Þ are called the extrinsic parameters. Given the intrinsic and extrinsic parameters, for each voxel, we can¯nd the corresponding pixel in the image by projecting the voxel onto the image plane. Real lensed usually have some distortion, mostly radial distortion and slight tangential distortion. So the reality coordinates of the projection point in the image plane was di®erent from the theoretical coordinates, and they can be expressed as below: where (x; yÞ are the theoretical coordinates, and (x 0 ; y 0 Þ are the actual coordinates, (k 1 ; k 2 ; k 3 . . .) are radial distortion coe±cients, (p 1 ; p 2 Þ are tangential distortion coe±cients, r is the distance between the point and the principal point, and r is calculated as Eq. (4).
usually, higher order coe±cients are not considered. So, in this work the lens distortion coe±cients means (k 1 ; k 2 ; p 1 ; p 2 ). Before the algorithm used to determine the 3D shape from the silhouettes it could be applied to the image sequences, the intrinsic and extrinsic camera parameters had to be determined. To obtain the intrinsic matrix, the camera was calibrated using a calibration pattern designed with easily detectable features. 30 The pattern was printed and pasted on a plastic plate with a 10 cm diameter hole in the center. This calibration object was placed on the potted plant so that the plant passed through the hole. Before the plant image acquisition process was performed, we obtained approximately 15-20 images from various angles without the plant and used these images to determine the intrinsic matrix of the camera and the distortion coe±cients of the lens. When we took photographs of the plant with the calibration pattern, we could then calculate the extrinsic matrix (i.e., the camera position and orientation) for each image from the extracted pattern in the image.

Volume reconstruction
Reconstruction from multiple views was the method used in this work. The object space, including the plant to be reconstructed, was divided into cells, which are called voxels in computer graphics. Initially, all the voxels were assumed to be empty. Reconstruction was accomplished by projecting the voxels onto the images to¯nd the voxels belonging to the plant.
The center of the calibration pattern was the origin of the 3D coordinate system. Based on the coordinate system, a rectangular bounding box was constructed around the plant (with dimensions of x: 40 cm, y: 40 cm, z: 50 cm) to ensure that the entire plant was included. The bounding box was represented by an array of voxels. The coordinates of each voxel were de¯ned by the center of the voxel. The coordinates of the voxel at the origin of the system were (0, 0, 0); the other voxels were created sequentially, and the coordinates were computed based on the voxel size. The voxel size was user de¯ned, and the size determined the spatial resolution (SR). The SR depends on the image resolution. And the SR could be calculated as below: where, f is the focal length of the lens, D is the distance between the object and the camera, is the size of the image unit. The camera used in this work is AVT Stingray F-504C, is 3.45 m, f is 12 mm, D is between 1000$1400 mm, so the SR is between 0.29$0.40 mm, which means the voxel size could not be smaller than 0.4 mm. Each image was processed as described in the previous section. As shown in Fig. 3(c), the pixels belonging to the plant were set to white (1), and the background was set to black (0). Each voxel in the bounding box was then projected onto the image plane, and the corresponding pixel i in image i was obtained using a perspective transformation. In this way, for each voxel, a set of pixels (pixel 1 ; . . . ; pixel k ) were obtained. The photo-consistency of each voxel was determined from the pixel set using Eq. (6). Nonphoto-consistent voxels (i.e., voxels for which consistent ¼ 0) were removed, and the plant volume consisted of the remaining voxels.
where k is the number of images; pixel i indicates the value of the pixel projected onto image i , and T indicates the threshold, which was 0.8 in this work.

Phenotypic parameter extraction
For the phenotypic study, several phenotypic parameters of interest were extracted from the volume model, including the plant height, plant volume, number of tillers, tiller angle and stem diameter. This section introduces the process used to extract these parameters.

Volume
The plant model obtained from the volumetric reconstruction consisted of many voxels. Therefore, the volume (V c ) of the 3D model of the plant could be obtained from the number of voxels (N v ) and the voxel size (scale) using Eq. (7).

Number of tillers and tiller angle
To determine the tiller traits, a 3D thinning algorithm was used to obtain the skeleton of the plant from the volumetric model. The 3D skeletonization provided shape features that were then extracted from the model. In discrete spaces, the thinning process iteratively deleted all the simple points whose deletion did not alter the topology, until the skeleton was¯nally obtained (Fig. 4). The tiller traits that describe the PA include the number of tillers, the tiller angle and the stem diameter. To analyze these tiller traits, the plant skeleton model must be partitioned into several individual tillers using a 3D region-growing algorithm.
At the bottom of the plant skeleton model, a nonlabeled voxel was chosen as the starting point and labeled as a stem. Beginning at the starting point, the algorithm then searched upward through all the connected voxels and labeled them as a stem or a leaf according to the following steps (de¯nitions of the terms used are given in Fig. 5): (1) The voxel at the top of a stem was used as the current voxel. Nonlabeled voxels were then visited beginning at the current voxel until a branch point was visited; all the visited voxels were then labeled as a stem.
(2) A line was¯t to the stem voxels, and the direction vector (V stem Þ of the line was calculated. (3) For i ¼ 1; 2; 3 . . ., visited N voxels on branch i (where N was a parameter used to eliminate false branches; in this paper, if the length of the branch was shorter than 2 cm, which means N ¼ 20 for the voxel size of 1 mm, then the voxels were left unlabeled) and calculated the direction vector (V bi Þ as in step 2. (4) For i ¼ 1; 2; 3 . . ., if the angle between V stem and V bi was less than 3 , labeled the N voxels on branch i as stem, otherwise labeled them as leaf. (5) Steps 1 to 4 were repeated until all the connected voxels were visited.
After the steps above were performed, the stem and some of the leaf voxels were labeled; if the stem was long enough and had more than one leaf on it, the stem was considered to be a tiller. The individual tillers were extracted, and the number of tillers was obtained. The axis of each tiller was determined from the line¯tted to the stem above. The inclination of the tiller was expressed as the angle between its axis and the z-axis (i.e., upward) of the world coordinate system, which could be calculated from the direction vector of the line.

Stem diameter
To obtain the size of a tiller, its axis was mapped to the volumetric model. On the axis, we labeled M points as samples at intervals of length d from the origin of the axis (M and d were de¯ned by the user according to di®erent requirements). For each sample point, using the following steps, a planar radius R could be obtained to¯t the cross-section of the tiller to a circular area in the volumetric model: (1) We started with a radius R 0 ; R 0 was initialized to the minimum of the distances from the axis to each border point. (2) The radius was increased by one unit, and the gap fraction (i.e., the proportion of empty pixels) was computed for the increased annular zone.  5. De¯nitions of terms used in the tiller extraction and tiller angle calculation algorithms. A branch point is a voxel on a stem that has more than one connected voxel. N is a parameter used to eliminate false branches (e.g., branch 1 on tiller 2) in the model resulting from reconstruction. Additionally, N was used to avoid connected search errors caused by leaf contact between tillers. Usually, N was set to 20 for a voxel size of 1 mm. (3) Step 2 was repeated, and R was obtained when the gap fraction was greater than a given threshold (the threshold was 0.5 in this paper) (Fig. 6).
All the programs used in this study were developed in the Cþþ language using Visual Studio 2010 (Microsoft) development tools. The 3D model was reconstructed from the image sequence using the method described above.

Performance validation
To validate the automated method for trait measurement, the di®erences between the automated measurements and the manual measurements were calculated. The estimation errors were expressed as the mean absolute percentage error (MAPE) and the root mean square error (RMSE) for the plant height and plant fresh weight. The MAPE and RMSE were calculated according to Eqs. (9) and (10): where a i and m i denote the automated measurements and manual measurements, respectively, and n denotes the number of plants. The correlation coe±cient (R) was also calculated to test the correlation between the automated and manual measurements. Figure 7 shows the rendered results of three plant models with a voxel size of 1 mm. All 3D volume models were rendered using Amira version 5.3.2 (Visage Imaging).

Measurement results
All the parameter results obtained from the volumetric models produced by applying our method to the wheat plant images are presented in this section. , and (c) shows a tiller section with a connected leaf. The black pixels represent the axis and the gray pixels represent the plant.
As described in the previous section, regression models were built between the parameters extracted from the 3D models (i.e., the model volume V c and model height H c Þ and the manual measurements to estimate plant fresh weight and calibrate the plant height. The model results and the statistical parameters are shown in Table 1. The regression analysis was performed to test the relationships between the estimated values and the manual measurements for the plant height and plant fresh weight. Scatter plots of the estimated values and the manual measurements are displayed in Fig. 8. The MAPE of the plant height was 2.71% (1.08 cm with an average plant height of 40.07 cm), and the MAPE of the plant fresh weight was 10.06% (1.41 g with an average plant fresh weight of 14.06 g). The RMSE was 1.37 cm and 1.79 g for the plant height and plant fresh weight, respectively. The scatter plots and linear regressions displayed in Fig. 8(a) and 8(b) indicate a strong correlation between the estimated parameters and the manually measured parameters, with correlation coe±cients of 0.95 and 0.96 for the plant height and plant fresh weight, respectively. The number of tillers was also extracted from the model. The results of the tiller   number estimation are displayed in Fig. 8(c). As shown in the scatter plot, the absolute error in the number of tillers was within AE1 for most plants, with an RMSE of 1.01.

Computational optimization
The data acquisition and processing°ow were designed for a high-throughput phenotypic platform. Consequently, it was necessary to signi¯cantly reduce the processing time for each plant. In the reconstruction process, each voxel was projected to each image with Eq. (2), and then there was no data dependency among di®erent projections. Therefore, it was able to parallelize the calculation.
The pseudo-code of the GPU based reconstruction algorithm was shown as below: Graphics processing unit (GPU) parallel computing technology was used to accelerate the reconstruction module. The program was coded and executed on a central processing unit (CPU) and on a GPU device (on an computer equipped with an Intel Core2Quad Q8200 (2.33 GHz) processor with 3.00 GB main memory and an NVIDIA GeForce 9800 GT graphics card that supported parallel computing). The performance and results are shown in Table 2.
As shown in Table 2, using parallel computing on the GPU, a substantial increase was obtained in the computing performance for the plant volumetric model reconstruction. The speedup ratio ranged from 9.21 to 23.61 at di®erent voxel sizes (4-0.5 mm), and the di®erence between the volumetric models reconstructed using the CPU and the GPU was negligible. Additionally, it can be seen that for a voxel size of 1 mm, the time cost is almost 4 min (234594 ms) with the CPU but only 10 s (10281 ms) with the GPU, demonstrating the potential of real-time reconstruction using a GPU.

Discussion
As described above, this study presents a PA parameter acquisition method for wheat plants.
Through the development of a volumetric model processing and analysis methodology for plant phenotyping, this study demonstrates the automatic measurement of PA parameters, such as the plant height, plant fresh weight and number of tillers, and the potential to measure the tiller angle and stem diameter. Many phenotypic parameters of interest are not easily extracted from 2D images.  By providing more information on the plant structure, there is no doubt that 3D reconstruction could facilitate the measurement of most parameters with 3D information. However, despite its advantages in estimating parameters that are not easily extracted from 2D images, 3D reconstruction also has several limitations that should be further investigated. 7 Additionally, the poor e±ciency of the 3D reconstruction method severely limits its application in a high-throughput phenotyping platform. Several sensing techniques, such as laser scanning and the Microsoft Kinect system, have been used in the 3D imaging of plants for phenotyping. However, due to their ine±ciency in data acquisition, these systems are not suitable for a high-throughput platform. 20 In comparison, reconstruction methods based on camera imaging are simpler and more°exible for data acquisition. Coupled with the use of parallel computing technology for volumetric reconstruction, camera-based methods could have a great advantage in speed. An average computational time for model processing and parameter extraction of 20 s per plant was achieved in the present study, which means that thousands of plants could be processed in one day.
As described above, the reconstruction process was accelerated using parallel computing technology on a GPU. Although the calculation precision of the CPU and GPU was di®erent, the di®erence between the results obtained using the CPU and GPU was negligible, as seen in Table 2 and Fig. 9. The voxel size is another factor a®ecting the reconstruction result. A smaller voxel size provides higher spatial resolution but also requires more computation time. As shown in Fig. 9, the reconstruction model becomes coarser as the voxel size increases. A clear sawtooth e®ect can be observed at voxel sizes of 2 mm, 3 mm and 4 mm, and fracture of the model can be seen at a voxel size of 4 mm. Although we could obtain better results at voxel sizes of 0.5 mm and 1 mm, considering computational e±ciency, a voxel size of 1 mm provided su±ciently accurate reconstruction and was chosen in this study.

Conclusions
In this study, we present a methodology based on volumetric model processing, including volumetric model reconstruction from 20 plant images, volumetric model processing and data regression analysis, which was applied to the extraction of Triticum aestivum PA parameters. Experiments were performed to test the method with 183 wheat plants. The MAPE for the plant height and the plant fresh weight was 2.71% (1.08 cm with an average plant height of 40.07 cm) and 10.06% (1.41 g with an average plant fresh weight of 14.06 g), respectively. The correlation coe±cients for the plant height and plant fresh weight were 0.95 and 0.96, respectively. Using a GPU, the image analysis was accelerated by a factor of 23, reducing the average computation time for each plant to only 20 s. In conclusion, the 3D volumetric reconstruction method for PA presented in this work has the potential to be applied in a high-throughput plant phenotyping platform.