3D Morphological Processing for Wheat Spike Phenotypes Using Computed Tomography Images

Abstract

Wheat is the main food crop today world-wide. In order to improve its yields, researchers are committed to understand the relationships between wheat genotypes and phenotypes. Compared to progressive technology of wheat gene section identification, wheat trait measurement is mostly done manually in a destructive, labor-intensive and time-consuming way. Therefore, this study will be greatly accelerated and promoted if we can automatically discover wheat phenotype in a nondestructive and fast manner. In this paper, we propose a novel pipeline based on 3D morphological processing to detect wheat spike grains and stem nodes from 3D X-ray micro computed tomography (CT) images. We also introduce a set of newly defined 3D phenotypes, including grain aspect ratio, porosity, Grain-to-Grain distance, and grain angle, which are very difficult to be manually measured. The analysis of the associations among these traits would be very helpful for wheat breeding. Experimental results show that our method is able to count grains more accurately than normal human performance. By analyzing the relationships between traits and environment conditions, we find that the Grain-to-Grain distance, aspect ratio and porosity are more likely affected by the genome than environment (only tested temperature and water conditions). We also find that close grains will inhibit grain volume growth and that the aspect ratio 3.5 may be the best for higher yield in wheat breeding.

1. Introduction

The global agriculture demand is expanding rapidly, not only because of the global population growth, but also because of climate change and resource depletion. Agriculture is always facing challenges to provide adequate food. Wheat is one of the most popular crops [1], which supports 20% of the daily protein and calories for 4.5 billion people [2]. And it is also one of the most important crops in the world.

In order to provide sufficient production, new technologies are required to speed up breeding, the most important one being genetic breeding [3]. By analyzing the inner relationships between the genotype and phenotype of a crop, scientists are able to modify the gene sections and thereby improve yields. Crop genetic breeding is a currently a significant research field of interest, and great results have thus far been achieved [4,5,6,7,8]. Compared to genotypes, the phenotype methods have developed much slower [9]. Obtaining wheat phenotypes information and studying wheat growth, performance and composition, are considered as a bottleneck in wheat breeding process. Presently, the collection of phenotypes depends on labor measurement, including counting the number of grains of a spike, measuring the length and weight of grains, and so on. However, manual measurement is time consuming, subjective, and lacking in accuracy [10]; therefore a fast and accurate phenotypic extraction method is imminent.

Some spike phenotypes are very important to wheat domestication [11], such as the quantity and weight of grains are highly correlated with the yields [12], and the variation in grain shape and size, which may also affect grain quality and yields [13]. There have been few studies focusing on the relationship between the spike and grain phenotypes. The main reason may be because of lacking an effective method to measure grain size, shape and other traits. To link genetic variation with plant phenotypes, we need an automated phenotypes extraction tool to accurately measure plant traits.

In this paper we introduce a 3D morphological image processing pipeline to automatically detect 3D phenotypes of wheat spikes from layer by layer scanned computed tomography (CT) images. It automatically and accurately detects grains and stem nodes from CT images and reconstructs 3D structure of wheat spikes. Illustrated in Figure 1, this overall automatic pipeline has three steps: (1) CT image pre-processing, (2) grain and stem detection and (3) 3D phenotypes extraction. Compared to the state-of-the-art methods, the proposed methods offer a more comprehensive 3D structure of wheat spike and a more accurate detection result. We do not process CT images layer by layer and then combine their results. We firstly combine all layers into 3D voxel images and then directly process the 3D images. This innovation enables us to identify grains and stems in a lossless way, by fully applying all CT images at the same time.

1.1. Related Works

Computer Vision has been applied to study plant traits for many years [15]. A variety of image techniques are used to capture trait information from different plant tissues. Visible light and laser scanning have good performances on measuring growth status, color and other external spike traits [16]. Visible light images are used to extract features and classify leaves [17]. Lidar is applied for surface trait extraction and the classification of 3D objects [18,19,20,21]. However, these techniques do not perform well on occluded or inner parts [18]. Therefore, we are considering X-ray imaging technology to collect plant surface and inner structure information. Compared with Positron Emission Computed Tomography (PET) and Magnetic Resonance Imaging (MRI) [22,23,24,25,26], Computed Tomography (CT) combines structural tomography and functional imaging, thereby more precisely scanning the physiological activity of plant tissue [27]. Presently, CT is widely used to assess tissue density [28], measure tiller numbers [29], grain quality [30], etc.

Since its introduction in the 1970s, CT has become an important tool for medical imaging to supplement X-rays and medical ultrasonography. As a non-invasive advanced imaging technology for high-resolution measurement of living tissue, CT provides inherent high-contrast resolution. It is able to distinguish tissues with physical density differences of less than 1% [31]. Modern optical 3D tomography and imaging techniques have been applied to reconstruct plant tissue in 3D for visualization and analysis, include plant vasculature [32], stem properties of sorghum [33], maize developing seeds [34] and plant root [35].

Even though there is a reasonable amount of research on analyzing medical CT images, the research on analyzing plant CT images is becoming increasingly popular in recent years. Hughes [14] introduces a CT image analysis pipeline based on morphology processing to extract spike and grain traits. This method is able to effectively and non-destructively analyze wheat traits, thereby enabling scientists to understand the environment effects and the gene functions.

Hughes’ method, however, cannot accurately identify wheat grains from a noisy background and tightly connected tissues, because it is based on 2D image processing and does not fully use the 3D information. We believe that the association information between CT image layers is very important. It is difficult to distinguish a grain and the background when only considering a single layer, especially when this layer only has a few pixels for the grain. This grain could be identified with another image if it has enough grain pixels. One may combine results from multi-layers after processing all images layer by layer. This way inevitably misses grain pixels, and it is indirect and complex. Inspired by Hughes’ method, we propose an effective but simple method to use all 3D voxels directly. All morphological processing is done in the 3D space, and pixels/voxels from all image layers are used simultaneously. This novel method significantly increases detection accuracy of wheat grain and stem nodes. The 3D morphology in image processing is a natural extension of the 2D version. Pardàs segments image sequencing by 3D morphological segmentation [36]. Maire uses the 3D morphological method to process X-ray tomography images and detect cellular ceramics [37]. We creatively introduce a CT image processing pipeline based on 3D morphology to measure plant traits automatically and nondestructively.

1.2. Contributions

(1)

We propose a novel CT image processing pipeline based on 3D morphology analysis. The proposed novel method is a fully automatic, highly accurate method that allows for rapid and nondestructive phenotypes extraction. Compared with the state-of-the-art algorithm, the proposed method fully uses all 3D information of multi-layers CT images in a straightforward manner.

(2)

We define a set of new 3D phenotypes of wheat spikes. The new 3D phenotypes include Grain-to-Grain distance, aspect ratio, porosity, angles between grains and stem. These 3D phenotypes enable breeding scientists to find the relationship between phenotypes and genotypes precisely.

(3)

We analyze the correlation among wheat grain traits and distinguish the traits that are more likely to be controlled by genome than by environments. The aspect ratio, Grain-to-Grain distance and porosity are slightly affected by the environments, we speculate that they may be mainly genetically controlled. We also find close grains will inhibit grain volume growth and the aspect ratio 3.5 may be the best for higher yields in wheat breeding.

2. Materials

2.1. Physical Equipment and Plant Materials

This experimental data comes from the data set published by Hughes’ work [14]. The CT images were collected by a μCT scanner developed by Scanco Medical, Switzerland, with model number μCT100. This scanner equipped with a cone beam X-ray source has power ranging from 20 to 100 kVp and a detector consisting of 3072 × 400 elements and a maximum resolution of 1.25 μm. The selected spring wheat (Triticum aestivum cv Paragon) was grown as single plants, and it was bred under various temperatures (25 °C, 30 °C, 35 °C and 40 °C) and water conditions (80% field capacity as high water and 40% field capacity as low water) until the six leaf stage. The number of wheat spikes is 90 and the number of wheat grains is around 3364.

For each scanning treatment, twelve representative spikes were selected and placed in a plastic holder. The maximum height of the holder is 70 mm. A few spikes were cut into two parts to fit the holder and each part is scanned separately. Sample preparation and loading into scanner cost about 30 min per 12 samples, and each scanned around 40 min. Figure 2 shows a slice image sample, from which we are able to clearly identify the wheat structures including grains, pedicel, floret and stems from the background.

2.2. Manual Labeling of Wheat Grains

In order to evaluate the proposed method, we manually labeled the data set as ground trues. There are 188 sets of unprocessed wheat spike CT images published, containing approximately 3000 grains. We selected the first 58 sets (only 43 sets contain grains, other 15 sets have not grain) of CT images and manually labeled all grains layer by layer. The manually labeled treatments were grown in the conditions of 40 °C high water, 40 °C low water, 35 °C low water. From visual checking, the intensity of the grain region is significantly higher than the background and has a sharp edge. An operator used Abode Photoshop software to do manual label work. Each layer image was manually given a binary image with the same size of the CT image. The grain region is highlighted in the binary image. The labeled 43 sets contain about 500 grains. Figure 3 shows three CT images and the corresponding labeled binary images. The labeling work took around 3 months to complete by an operator. Because the manual labeling is time consuming, we only processed about a quarter of the data set published by Hughes [14].

3. Method

In this paper, we introduce an image processing pipeline based on 3D morphology to automatically detect 3D phenotypes of wheat spikes from CT images. Illustrated in Figure 1, this method has three folders: (1) CT image pre-processing, (2) Block-like tissue detection, include grains and stem nodes, by 3D morphological processing and (3) 3D phenotypes extraction. We combine all layers into 3D voxel images and detect grain and stem nodes by the 3D morphological processing. This method simultaneously uses all pixels/voxels from all slices and is able to identify grains and stems very precisely. Before detecting block-like tissue, we do holder clearance and intensity equalization to increase the efficiency of recognition, as shown in Figure 4.

3.1. Holder Clearance

The CT images are scanned slice by slice along the Z-direction. By stacking all slices along the Z-direction, the 3D CT images are easily obtained. We do all 3D morphological processing and object detection directly on the 3D images, the stack processing is as shown in Figure 5.

The wheat spike is scanned with a supporting holder. Thereby, each CT image has a circle around the wheat spike. The holder should be removed for further processing of the images. Because the holder is a fixed size plastic container, a simple Region of Interest (ROI) method is good enough to remove the holder. In our experiments, we only kept the images inside the holder area. The holder clearance processing is shown in Figure 6.

3.2. Intensity Equalization

The CT image intensities diverse across sample groups because the scanning parameters differ from different time, which leads to a diversity of image intensities in different sets. We use a linear stretch mapping to equalize images of all sets. The linear stretch mapping follows the equation:

$$I_{out\left(x,y,z\right)}=\frac{I_{in\left(x,y,z\right)}-\min \left(I_{in}\right)}{\max \left(I_{in}\right)-\min \left(I_{in}\right)}$$

(1)

where $I_{in\left(x,y,z\right)}$ means the intensity of image $I$ at the coordinate $x$, $y$, $z$. $\min \left(I_{in}\right)$ and $\max \left(I_{in}\right)$ are the minimum and maximum intensity of all CT images of each wheat spike $I_{in}$. $I_{out}$ is the equalized image. The image equalization process enables us to process all image sets with consistent parameters. The process of intensity equalization is shown in Figure 7.

3.3. Foreground Detection by GMM

After the holder clearance, there are only spike tissues and noise left in the CT images. The background is composed of noise, artifacts and near black pixels. The wheat structure has a stronger signal than the background, while the background has many more pixels than the foreground, as illustrated in Figure 8a. In the CT image of plants, the image intensity of one class follows Gaussian distributions, and the whole CT image follows GMM. Therefore, we use Gaussian Mixture Model (GMM) to identify background and foreground [38].

The GMM is a probabilistic model used to represent K sub-distributions in the overall distribution, where each sub-distribution follows Gaussian distribution. After image equalization, the background pixels are near black (not zero intensity for most cases). The intensity of noise is usually around 50, between the near-black pixels and plant tissues. The background has two dominant Gaussian models. We use GMM with three components to fit each image and detect the largest two Gaussian models, which are identified as the background. By filtering out the background, the histogram of foreground has two to three Gaussian models (shown in Figure 8b); however, this is not enough to identify grains, leaves and stem, because their intensity overlaps. We need to use their 3D spatial connections to identify grains and stem nodes.

By finding the two biggest Gaussian models, we can get the threshold to filter out the background. In our experiments, the found threshold is around 40 to 60, for different image set. We use threshold-based segmentation to remove the portion of lower intensity values. The remaining are of high intensity tissues, including grains and stem nodes, as shown in Figure 9.

3.4. Block-Like Tissue Detection by 3D Morphological Processing

After foreground detection, there are only patches with high intensity values from the image, including grains, stem nodes and thick leaves. Different tissues vary in size, thereby in an ideal situation, we can detect the grain by size. However, during the foreground detection, grains are not accurately separated from other tissues, e.g., leaves and stems. The nearby grains may stick together, which influence the detection of individual grains.

Hughes [14] detects grains in each CT slice image and identifies individual grains in the 3D space by 3D connected component analysis. This method cannot accurately identify wheat grains from a noisy background and tightly connected tissues. E.g., several morphological erosion steps are applied to split neighboring grains. A grain will be missed on the image which only has a few pixels of the grain. The lost grain cannot be recovered because the morphological processing is done independently on each slice image.

We believe that the association information between CT image layers is very important. The CT images are scanned slice by slice along the Z-direction. By stacking all slices along the Z-direction, the 3D CT images are easily obtained, we do all 3D morphological processing and object detection directly on the 3D images, as shown in Figure 1. The pipeline of the proposed 3D morphological processing. (a) An optical image of a wheat spike [14]; (b) Layer by layer scanned CT images of the wheat spike; (c) The 3D image from stacked CT images; (d) Detected grains and stems by the proposed 3D morphology analysis; (e) Extracted 3D phenotypes of the wheat spike. This innovation enables us to precisely identify grains and stem nodes, by fully applying all CT images at the same time. Here, we only consider block-like tissue including grains and stem nodes. Paper-thin tissue such as leaves and florets cannot be identified by this method.

3.4.1. Detection of Grain and Stem Nodes

Erosion is one of the fundamental operations in morphological image processing, which was firstly used on binary images, later being extended to grayscale images and it is currently used on 3D CT images. We use this method mainly for noise elimination and segmentation of independent grains. It is observed that the adhesion between the grains is relatively thin and can be separated by some conventional segmentation methods, such as Watershed [39]. We chose 3D morphological erosion, a simpler and better-performing method, to separate individual grains. Morphological erosion could be understood as the minimum value filter that takes the minimum value in a certain field of voxel instead of the original value. We use a structural element with six neighbors in the 3D space for the 3D erosion. The adhesion is usually at the edge of the grain, so the minimum intensity of its adjacency matrix is generally 0, the area of adhesion will be smaller after completing an erosion processing every time. After the erosion processing, we detect connected components, classify the components to be grains, stem nodes and other according to their size, as follows:

$$\{\begin{array}{c}{L}_{C_i}=Grainif{s}_{C_i}>10000\\ L_{C_i}=StemNodeif{s}_{C_i}>50\&\&s_{C_i}<2000\\ L_{C_i}=Otherif{s}_{C_i}<50\end{array}$$

(2)

where $C_i$ is connected component i, $L_{C_i}$ is the label of $C_i$, and $s_{C_i}$ is the size of $C_i$. The size unit is voxel.

Then we can get the accurate numbers of grains in a spike based on patches. The erosion processing is shown in Figure 10. Stem nodes have similar intensity as shown in Figure 11, so we can detect it in same way.

3.4.2. 3D Growth to Recover Grains

After the 3D morphological erosion, we get the exact position of each grain. However, the erosion process damages the grain, resulting in the detected grains smaller than actual ones. The 3D morphological dilation process is not able to recover the lost grain voxels as the ‘hit’ operation is too simple to handle complex situations. We use a 3D region growing method to restore grains, illustrated in Figure 12.

The detected grains by 3D morphological erosion are taken as seeds for the 3D region growing. The region growing starts from any seed voxel, and checks its six neighbors in the 3D space (see Figure 12g). If one neighbor voxel meets growing criteria, it will be assigned to the seed region, and used as seed for the next searching iteration. The 3D region growing follows breath-first strategy and stops when no new seed is found. The growing criteria is defined based on the following observation. Intensity of grains and stem nodes is usually higher than background, leaves and floret. Intensity of grain border and stem nodes is higher than one of grain inner part. Therefore, a candidate voxel is being assigned as a newly found seed only if all of its neighbor voxels meet the following criteria:

$$i_{min}>{\lambda }_1·i_{mean}or{i}_{mid}>{\lambda }_2·i_{mean}$$

(3)

where, for all the candidate voxel and its six neighboring voxels, $i_{min}$ is the minimum intensity, $i_{mean}$ is mean intensity, and $i_{mid}$ is the median intensity. ${\lambda }_1$ and ${\lambda }_2$ are adjustment coefficients, whose values may differ with crop types. In our experiments, ${\lambda }_1$ and ${\lambda }_2$ are set to 0.87 and 0.99. The minimum and mean intensities are to verify whether there are low-intensity voxels in the neighbor. If low-intensity voxel exists, the candidate voxels should be considered as outside of grains.

3.5. 3D Phenotypes Calculation

After detection of grains and stem nodes, we accurately obtain the voxels set which belong to grains and stem nodes. The volume of the grain can be easily calculated by counting the number of voxels in the point set. We can also find the pores inside the grain and then calculate the volume of the internal pores. Some phenotypic information such as length, radius and surface area seem more difficult to be calculated. We also compute the angle of grains between stems. Those wheat traits are important to evaluate wheat grain features and crop yields.

3.5.1. Grain Length and Radius

We use Principle Component Analysis (PCA) [40] to compute the length and radius of grains. For a 3D object block with an irregular shape, it is difficult to calculate its length directly. Because this grain may be tilted, we cannot use its distance on the z-axis to indicate the length. The grain’s shape is irregular, we cannot change its tilt angle and z-axis projection distance to length. We use PCA to detect the three principle components of a grain, and the grain length and radius are computed in the PCA space.

We use the 3D coordinates of all voxel center of each grain as input for PCA. The voxels are from stacking layered CT images in the 3D space. PCA applies orthogonal transformation to convert the block voxels into a set of linearly uncorrected variables, each is a principle component. By sorting the eigenvalues in descending order, we get three eigenvalues ${\mu }_1$, ${\mu }_2$, and ${\mu }_3$ (${\mu }_1\ge {\mu }_2\ge {\mu }_3$), and three corresponding eigenvectors $v_1$, $v_2$, and $v_3$. By projecting grain voxels onto each eigenvector, we will get length $l_1$, $l_2$, and $l_3$ correspondingly on each direction. Therefore, the grain length $L_g$ is $l_1$ and grain radius $R_g$ is the mean of $l_2$, and $l_3$. Figure 13 shows the eigenvalues as well as grain length and radius.

3.5.2. Grain Volume and Surface Area

Grain volume is calculated by simply counting voxel number for each grain. After detecting a single grain, we obtained a voxel point set. All the points in the point set form a complete grain, and we can count the number of points to get the volume of the grain. There are pores in some of the spikes. We fill the internal pores by morphology closing and calculate the volume again. The pore volume is the difference value between the two volumes. Grain surface area, however, cannot be simply calculated in this way. The edge voxels are jagged and not smooth, therefore we cannot only count the outermost voxels as exposed surface. So, we estimate the surface area by assigning different weights to the number of surfaces so that the surface voxels are exposed to the background.

Windreich [41] extended the 2D perimeter estimation theory to 3D surface area estimation. Their algorithm begins by detecting all surface voxels and their six-connected voxels. The surface voxels are divided into nine categories. The surface area is estimated as a linear combination of class member values $\left\{N_i\right\}$:

$$\widehat{S}={\displaystyle \sum }_{i=1}^9{W}_i{N}_i$$

(4)

$W_i$ means weights for 9 surface voxel types, shown in Figure 14. By calculating the number of these nine types of voxels, we can then accurately obtain the surface area of the grain. They defined the voxel classification scheme, and Mullikin [42] determined the first three weights $W_{1\text{–}3}$ associated with voxels in classes $S_{1\text{–}3}$. These weights are $W_1\approx 0.894$, $W_2\approx 1.3409$, $W_3\approx 1.5879$. The weight $W_{4\text{–}9}$ associated with voxel classes $S_{4\text{–}9}$ are determined by Windreich [41]. These weights are $W_4=2$, $W_5=8/3$, $W_6=10/3$, $W_7=1.79$, $W_8=2.68$ and $W_9=4.08$. Voxels of types $S_{1\text{–}3}$ usually appear in some planes, types $S_{4\text{–}6}$ often appear in curved border regions. Voxel types $S_{7\text{–}9}$ only exist in extreme situations, such as plane, line or point alone. In our experiments, we only count the number of voxel types of $S_{1\text{–}6}$ because of no extreme situations.

3.5.3. Aspect Ratio and Porosity

After computing grain length, radius and volume, we use these parameters to describe grain shape. We firstly define grain aspect ratio to be the ratio of length to radius. The larger aspect ratio indicates that the grain is closer to the slender shape, and the smaller value indicates a thicker and shorter grain. The aspect ratio of grain is related to the milling yields of wheat [43]. The porosity is defined as the ratio of the pore volume to the total volume of the grain which may be related to the moisture content of the grain [44]. The aspect ratio and porosity of the grain describe the shape of individual grains.

3.5.4. Grain Angle and Grain-to-Grain Distance (GGD)

The above grain phenotypes are for individual grains. Another group of phenotypes that affect yields are the distribution of grains along the stem. We defined two features, grain angle and distance, to describe the grain distributions shown in Figure 15. Therefore, grain angle is defined as the angle between the stem node and the grain direction. The stem nodes are extracted by 3D morphology process described in Section 3.4. The connection line between stem nodes describes the position of the straw, the size of the node describes the thickness of the straw. We define the Grain-to-Grain distance (GGD) of one grain to be the average distances between the grain to its two adjacent grains. The GGD presents the grain number in a unit distance along the stem.

4. Results and Discussion

4.1. Results of Grain Detection

After morphological erosion and growth processing, we can separate all the grains from a single wheat spike one by one. We compared the grains that we detected, the ones we labeled and the ones detected by the Hughes’ method. To evaluate the accuracy critically, we used the experiment results that they published. The three datasets are shown in Figure 16, where we can clearly see that the results by the proposed methods fit the ground truth data in a better manner. Hughes’ method erroneously detects many artifacts, and fails to detect the grain pit.

We use some quantitative indicators to evaluate our method as follows:

• TP (True Positive): the number of voxels which belong to ground truth grains among the detection results.
• FP (False Positive): the number of voxels which do not belong to ground truth grains among the detection results.
• FN (False Negative): the number of voxels which belong to ground truth grains but have not been detected.

Based on the above definitions of TP, FP and FN, we calculate three metrics, the precision, recall and F-Score to evaluate our method and Hughes’ method, the three metrics are explained below:

$$\mathrm{precision}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}\times 100\%$$

(5)

$$\mathrm{recall}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}\times 100\%$$

(6)

$$\mathrm{F}1-\mathrm{Score}=\frac{\mathrm{precision}\times \mathrm{recall}}{\mathrm{precision}+\mathrm{recall}}\times 2\times 100\%$$

(7)

Table 1. Detection results of the proposed method and Hughes’ method in 11 wheat spikes.

Treatment No.	Grain Number	Our Method			Hughes’ Method
		Precision	Recall	F1 Score	Precision	Recall	F1 Score
B1165	22	96.96%	98.64%	97.79%	95.78%	98.60%	97.17%
1189	9	95.36%	99.26%	97.27%	84.11%	95.50%	89.44%
T1180	8	95.40%	99.29%	97.31%	78.81%	99.72%	88.04%
B1178	36	96.37%	98.83%	97.59%	94.46%	98.27%	96.33%
1175	3	96.26%	99.11%	97.67%	83.23%	99.74%	90.74%
T1178	31	94.75%	98.93%	96.80%	88.10%	94.57%	91.22%
1170	15	98.72%	96.78%	97.74%	90.29%	98.52%	94.23%
1188	1	92.61%	98.01%	95.23%	68.08%	99.16%	80.73%
T1183	2	96.56%	92.90%	94.70%	83.03%	93.20%	87.82%
1172	1	89.22%	99.85%	94.24%	72.63%	99.76%	84.06%
Total	128	96.36%	98.49%	97.42%	90.49%	97.65%	93.93%

Note: B and T ahead of Name means the bottom and top of the spike, because some spikes are too long to scan so be cut into two parts. The spike treatment No. 1172 only contains one grain but also contains many other tissues, so its precision is much lower than others. The image of this spike is shown in Figure 17.

lists the evaluation results of the proposed method and Hughes’ method for 11 wheat spikes. The manual label data is used as ground truth data, and both methods are compared with the ground truth data. From Table 1, we can see that our method is better than Hughes’ method no matter the precision, recall and F1-Score. Our mean precision, recall and F1-Score for the 11 wheat spikes are 96.36%, 98.49% and 97.42%, while Hughes’ are 90.49%, 97.65%, and 93.93%. The quantitative evaluation matches well with the visualization evaluation. Hughes’ method only reaches precision 78.81%, 68.08%, 72.63% for spike treatment No. T1180, 1188, and 1172. While our method reaches precision 95.40%, 92.61% and 89.22%, respectively. In the most extreme case, 1188, our method increases the precision 24.53%, from 68.08% to 92.61%.

Compared with Hughes’ method, we detect grains directly on the 3D image rather than on 2D images. Because there is a strong correlation between adjacent 2D layers, the detection on the 2D layer will lose this part of the association, so the result will be biased. For example, it is hard to judge whether some discrete color patches with small area are noise or grains on a single layer. The 2D erosion will directly delete these pixels, which cause the grain tips to be deleted (Figure 18). While 3D erosion more precisely detects grains on 3D image as it directly uses all 3D voxels. Additionally, the region growing method is able to recover the eroded voxels. Therefore, the shape of the grains detected by our method is more regular; there is no obvious jaggedness, matching the real situation well.

4.2. Grain Number Counting

Grain number for each wheat spike is crucial to evaluate spike yield. Hughes also evaluates grain number and provides ground truth of grain number for each grain spike. However, after manually counting grain by ourselves, we find that our manually counted number is different to Hughes’ manual counting. Therefore, we list four grain numbers for each grain spike in

Table 2. The number of grains counted by the proposed method. The numbers in red show the difference between the ground truth data between ours and Hughes’, the numbers in green mean the detected grain is different from our ground truth data.

Treatment No.	Our Method	Our Ground Truth	Hughes’ Method	Hughes’ Ground Truth	Treatment No.	Our Method	Our Ground Truth	Hughes’ Method	Hughes’ Ground Truth
1170	15	15	15	15	118A	0	0	2	0
1171	0	0	1	0	118B	2	2	2	2
1172	1	1	1	1	118C	1	1	2	1
1173	0	0	1	0	118D	36	36	36	35
1175	3	3	3	3	118E	4	4	4	4
1176	1	1	1	1	1189	9	9	8	8
1177	1	1	2	2	1160	14	14	14	14
1178	67	67	64	62	1161	23	23	24	22
1179	0	0	2	0	1162	22	22	22	22
117A	5	5	5	5	1163	24	24	25	24
117B	0	0	1	0	1164	22	22	23	22
117C	0	0	1	0	1165	26	26	27	26
1180	43	44	44	44	1166	30	30	32	30
1182	31	31	31	31	1167	18	18	19	18
1183	23	23	24	23	1168	20	20	20	20
1184	13	13	14	13	1169	14	14	15	14
1185	0	0	1	0	0116A	14	14	15	14
1188	1	1	2	2	0116C	25	25	27	24
Total	508	509	530	502

: the number counted by our method, by our manually counting, by Hughes’ method, and by Hughes’ manual counting. We can see in most cases, the ground trues data provided by Hughes is less than the ground truth provided by us, see spike No. 1178, 1188, 118D, 1161 and 0116C in Table 2 (shown in red). It is possible that the hidden grains are hard to find by operators. While in the spike No. 1177 and No. 1188. Hughes counts one more grain than us, we counted the grains several times very carefully and therefore we are sure that our manual counting is correct. Some of the differences in manual counting are shown in the Figure 19. The results of the comparison are shown in Table 2.

Table 2 shows that the number of grains we detected is almost the same as our ground true data. We can also see that our proposed method performs better than Hughes’ ground true data, which is the normal human performance. In the whole dataset, the proposed method only missed one grain, which was in wheat spike No. 1180 (Figure 20). This spike had a large number of grains, which were close and hard to be identified.

Compared with Hughes’ method, the accuracy of our method increases from 95% to 99.8%. In most wheat spikes, Hughes’ method erroneously counted 1 to 2 grains more than our ground truth data. In spike treatment No. 1178 and 1189, Hughes’ method missed several grains. However, our method successfully identified most grains, only missing 1 grain in 509 grains.

4.3. Stem Node Detection

The detected stem nodes are illustrated Figure 21. The stem node detection is for calculation of grain angle, no high detection accuracy is demanded. On the other hand, the manual labeling is time-consuming. Therefore, we do not manually label the stem node and do not evaluate the detection accuracy. From Figure 21 we can see most stem nodes are detected well, except for the small nodes on the top part of a grain.

4.4. Grain Size and Angle

We calculate the grain size and angle for all wheat spikes and all grains. The CT that we used has an ultra-high resolution (this data accuracy is 68.8 μm per voxel in our experiments), which is able to accurately record the internal structure of plant tissue without damaging the plant. We calculate volume, surface area, length and radius for each grain. Some grains contain pores; we also use CT images to calculate the volume of grains quite well. The volume of the grain describes the size of one individual grain, the volume of all grains on the spike indicates the yields of one wheat spike. We can find the connection point between the grain and the straw by the oblique angle. Figure 22 shows the grain serial number of spike No. 1170 and zoom-in visualization of grains. The calculated wheat grain traits are shown in

Table 3. The phenotypic information of grains in spike treatment No. 1170.

Grain	Angle	Length	Radius	Volume	Pore Volume	Surface Area	GGD
No.	(°)	(mm)	(mm)	(mm3)	(mm3)	(mm2)	(mm)
1	20.62	6.2	1.97	36.08	1.75	82.04	9.49
2	13.64	5.8	1.93	33.01	0.10	84.8	9.46
3	14.64	5.11	1.59	23.18	0.04	48.37	4.16
4	14.37	6.06	1.93	38.41	0.04	76.12	4.18
5	31.35	6.41	1.99	32.5	2.43	92.29	4.35
6	24.75	6.06	2.00	39.44	2.53	85.31	4.97
7	19.55	5.84	1.95	34.93	0.89	89.58	7.72
8	21.41	6.36	1.97	35.39	1.84	102.19	6.87
9	17.18	6.17	1.79	34.7	2.08	71.80	4.23
10	19.58	6.78	1.89	43.35	0.20	86.32	4.91
11	30.28	7.38	2.02	43.86	2.59	98.99	4.72
12	8.08	6.61	2.03	44.27	0.04	82.41	6.06
13	17	7.00	2.09	46.45	2.52	104.13	6.65
14	9.73	7.02	2.06	50.11	0.04	84.96	7.83
15	22.33	6.16	1.81	35.99	0.07	70.62	9.30

.

4.5. Analysis of 3D Phenotypes

We extracted phenotypic information for all data published by Hughes. We detected 3364 grains in 188 sets of images; they are raised under different temperature and water conditions. In this section we analyze the environment affection on the grain traits, thereby finding out which grain traits are significantly controlled by environment and which are controlled by gene. All phenotypic information will be included as an attachment at the end of the manuscript.

4.5.1. Volume Distribution

We first calculated grain volume (including internal pores), length, radius, tilt angle, surface area and so on for a total of 508 grains (manually labeled data by us). The volume distribution of all grains in the 43 spikes is plotted as a box diagram shown in the Figure 23. The volume value means the volume of a single grain. We can see that the grains of one spike have different volumes. However, the volume RMS for each spike is small comparing to the volume difference between spikes. We can see that the mean volume and sum volume are highly related to environments.

4.5.2. Factors Affecting Spike Yields

We summarized the grain volumes of all single spikes under different temperature and water conditions, and found that the single spike had the highest yields at 25 °C in low water condition. Figure 24 shows that the average wheat yields decrease from 1900 mm³ to 200 mm³ from 25 °C to 40 °C in high water condition. The average wheat yields decrease from 2200 mm3 to 400 mm3 from 25 °C to 40 °C in low water condition. The experiments indicate that low water condition produces a better yield than high water condition in any template situation, while a low temperature is better than a high temperature.

4.5.3. Factors Affecting Grain Phenotype

In addition to measuring the yields of each wheat spike, we analyzed the relationships between grain shapes with temperature and water condition, shown in Figure 25. The grain shapes include grain angle, Grain-to-Grain distance, aspect ratio and porosity. The grain angle and aspect ratio are minimally affected by the external environments. We estimate that these traits are primarily controlled by genes. The Grain-to-Grain distance and porosity are greatly affected by the temperature; the influence of water conditions is minimal. We noted that although the yields of wheat grain is not high at 35 °C, the grain is the most full and the shape is the most stable. If we need more full grains, 35 °C may be a better choice.

4.5.4. Correlation Between Volume and Shapes

In this experiment, we mainly explored factors affecting yields of individual grains. In addition to external environment changes, we also explored the relationship between the volume and the extracted phenotypes, include Grain-to-Grain distance, grain angle as shown in Figure 27. Grain-to-Grain distance and grain angle indicate the distribution of the grain on the spike, the aspect ratio describes the shape of the grain. Illustrated in Figure 26, the relationship of phenotype to volume are grouped in box plots for each trait. One sub-graph only considers one environment condition. For example, sub-graph (a) considers temperature, so the grains grown in both high water and low water conditions are counted together.

As discussed in Section 4.5.2, we know temperature and water condition significantly affect the wheat yield. This fact is also confirmed by the Figure 26, where we can see the temperature and water condition have a similar effect on the relationship between volume and Grain-to-Grain distance, angle and aspect ratio. Therefore we only need to analyze the relationship of phenotype to volume in one condition, for example the distance to volume in temperature 30°.

Figure 26a,b show the relation of Grain-to-Grain distance to volume in 25°, 30°, 35° temperature conditions and high water and low water conditions. The distance for one grain is the average distances between the grain to its two nears grains. Usually one grain will not be 10 mm further away from other grains, only if the grains between them are not developed. If one only considers the distance below 10 mm, we can see the grains with smaller volume have shorter distances, meaning that grains that are too close will inhibit their growth (shown in Figure 26a,b).

It is interesting to see that the grain with a large angle has a relatively large volume in all temperature and water conditions (shown in Figure 26c,d). We speculate that larger grains are more inclined to be affected by gravity during growth. The heavier weight pulls the grain closer to the ground thereby resulting a bigger angle with the stem.

Figure 26e and f show the relation of aspect ratio to volume in the three temperature conditions and two water conditions. From the Figure 26e,f, we see the aspect ratio to volume curve follows the parabolic distribution, the volume reaches highest volume around aspect ratio 3.5. From the Figure 25c we know that the environment condition slightly affects the aspect ratio. It means that the aspect ratio is more likely affected by genome than environments, suggesting that aspect ratio 3.5 may be the best shape in wheat breeding. Figure 27 shows volume, volume, length and radius are highly related.

5. Conclusions

In this paper we introduce a novel CT image process pipeline based on 3D morphology analysis to detect wheat grains and stem nodes from CT images. This method is able to accurately extract the phenotypic information and traits of plant tissues. The experimental results show that our method is able to count grains more accurately than normal human performance. Our method increases the average accuracy from 90.5% to 96.4% and the average recall rate reach 98%. Our subsequent experiments on the wheat stem node also indicate that this method not only works on seed identification, but can also be applied to other similar plant tissues.

We also define a set of new 3D phenotypes of wheat spikes, including aspect ratio, porosity, Grain-to-Grain distance, and grain angle. By analyzing these traits with the environment conditions, we find that the distance, aspect ratio and porosity are fairly affected by the environment, the grain angle is linear corrected with grain volume. The grains with smaller volume have shorter distances, meaning that too close grains will inhibit their growth. The aspect ratio is more likely affected by genome, the aspect ratio 3.5 maybe the best for higher yields in wheat breeding.

The 3D morphology analysis detects block-like tissues include grains and stem nodes quite well, however, cannot extend to paper-thin tissues such as leaves and floret. Those tissues are also very important to plant trait study. We believe that CT is a good imaging technique for analyzing plant phenotypes. In future work, we will study other methods to detect paper thin tissues and more plant phenotypic traits based on CT Images. The phenotype–genotype analysis uncovers that Grain-to-Grain distances and aspect ratio are more likely to be controlled by gene. However, the number of wheat spike is only 90, which maybe not enough for concrete conclusion. We will extend the experiments to be a much larger number of spikes.