Open AccessArticle

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

Biao Xiong

Bo Wang

Shengwu Xiong

^*,

Chengde Lin

and

Xiaohui Yuan

School of Computer Science and Technology, Wuhan University of Technology, Wuhan 430070, China

Author to whom correspondence should be addressed.

Remote Sens. 2019, 11(9), 1110; https://doi.org/10.3390/rs11091110

Submission received: 6 April 2019 / Revised: 28 April 2019 / Accepted: 6 May 2019 / Published: 9 May 2019

(This article belongs to the Special Issue Advanced Imaging for Plant Phenotyping)

Download

Browse Figures

Graphical abstract
"> Figure 1
The pipeline of the proposed 3D morphological processing. (a) An optical image of a wheat spike [<a href="#B14-remotesensing-11-01110" class="html-bibr">14</a>]; (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. The whole pipeline is done automatically while some parameters are given manually. "> Figure 2
A slice image sample and representations of wheat spike tissue. "> Figure 3
Labeled grains images. The top row images are sets of raw slices after holder clearance processing, and the bottom three are the corresponding three labeled binary slices. "> Figure 4
The pipeline of pro-processing of CT image. The parts are all automated. "> Figure 5
Stack 2D image sequence into 3D image. (a) The original group of 2D images; (b) Three example images; (c) The 2D matrices of each image are stacked along the vertical direction. Their horizontal coordinates remain unchanged. The coordinate values in the vertical direction are their numbers in the image group; (d) The result of stacked 3D images. "> Figure 6
Holder clearance processing. (a) The raw image; (b) Image after removing holder using Region of Interest (ROI). "> Figure 7
Image intensity equalization. Top three are images of different wheat spikes with different intensities, while the bottom three are equalized images correspondingly. "> Figure 8
CT image histogram of a wheat spike. (a) The complete histogram; (b) the histogram between 50 and 255 intensity. "> Figure 9
Foreground detection by GMM. (a) The raw CT image; (b) Image after threshold segmentation. "> Figure 10
The 3D morphological erosion for separating connected objects. Left one is before erosion processing, we can see the adhesion between the right two grains; Middle and right images are the grains after one and two steps of erosion. The top row are three 2D images and the bottom row are the 3D views, correspondingly. In the red circles, we can see the adhesion has been removed, and each connected component is one grain. "> Figure 11
Detected grains and stem nodes. Stem nodes have similar intensity with grains, while stem nodes and grains are spatially close. "> Figure 12
The 3D region growing method to recover lost voxels. (a) The raw CT image; (b) Detected grain region by erosion processing; (c) The yellow voxels are candidate voxels by the first iteration; (d) The green voxels are newly grown grains; (e) The newly grown voxels are set to seeds; (f) Candidate voxels (in yellow) and newly grown grain voxels (in green) by the second growing iteration; (g) Six-neighboring voxels in 3D view; (h) An example of an voxel on grain surface and its six-neighboring voxels. "> Figure 13
Computation of length and radius of grains by PCA. "> Figure 14
The nine unique surface voxel classes. "> Figure 15
Stem detection and grain angle computation. (a) 3D images of a wheat spike; (b) Detected stem by 3D morphology; (c) Combined stem nodes and grains; (d) Computed grain angle. "> Figure 16
The detection results of the top part of wheat spike with treatment No. 1178. (a) Grains detected by Hughes’ method; (b) Grains detected by our method; (c) Manually labeled grains. The difference between our and Hughes results is marked with a red box. "> Figure 17
Spike treatment No. 1172 (upper one) and No. 1188 (lower one) contain one grain but some background too. The green grain is detected by our method, the red one is from Hughes. We can see that their grain is much bigger than actual value, which cause the precision loss. "> Figure 18
The details of the grain detection results for wheat spike (the bottom part of No. 1178). (a) Detected grains by Hughes’ method, visualized in red; (b) Detected grains by our method, visualized in green; (c) the two results are stacked together for evaluation; (d) Zoomed-in picture of (c). We can see that tips of grain (lower one) and some other voxels (middle one) are not detected using Hughes’ method, some other tissues (upper one) are erroneously detected as grains. "> Figure 19
The raw images of the three groups of wheat spikes and the counted grains. Some treatments contain 2 set of images, we all set for one wheat spike in the figure. The treatment No. 0116C contains 25 grains but Hughes counted 24; 01161 contains 23 grains but Hughes counted 22. The treatment No. 01179 has no grain both counted by us and Hughes, but one grain is counted using their method. "> Figure 20
The proposed method failed to find one grain in treatment No. 1180 spike because the grain is very close to another grain. "> Figure 21
Detected stem nodes and grains. (a) The whole wheat spike; (b) the detected stem nodes; (c) the detected grain and stem node. Upper row is spike No. 1170, middle row is bottom part of spike No. 1165 and bottom row is top part of spike No. 1178. "> Figure 22
The detected grains for wheat spike No. 1170. The grain number is sorted with their gravity height in descending order. We can see that No. 5 grain has bigger angle, No. 11-14 grains have bigger volume and No. 13 grain has biggest surface area. "> Figure 23
The volume distribution of grains in 43 groups of spikes; the black short line means average value of the spike. The volume distribution of all wheat spikes please refer to Additional file II: All Spike Volume.pdf. "> Figure 24
Spike yields according to different temperature and water conditions. Temperature and water have a great impact on production. The yields at 25 °C is nearly twice that of 35 °C. The yields in low water conditions will be slightly higher than the one in high water. "> Figure 25
(a–d) Single grain morphology under different temperature and water conditions. "> Figure 26
Correlation analysis of individual volume of grain between the key traits, including Grain-to-Grain distance, grain angle and grain aspect ratio. One sub-graph only considers one environment condition. (a,c,e) consider the effect of the Grain-to-Grain distance, angle and aspect ratio on volume at different temperatures. And (b,d,f) consider the effect of the Grain-to-Grain distance, angle and aspect ratio on volume under different water conditions. "> Figure 27
Correlation matrix of wheat phenotypic traits. The upper triangle is the graphical representation, and the lower triangle is the corresponding correlation coefficient. ">

Versions Notes

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.

Keywords:

3D morphological processing; 3D phenotypes; X-ray micro computed tomography; image analysis; grain and stem detection

Graphical Abstract

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_{o u t (x, y, z)} = \frac{I_{i n (x, y, z)} - \min (I_{i n})}{\max (I_{i n}) - \min (I_{i n})}

(1)

where

I_{i n (x, y, z)}

means the intensity of image

I

at the coordinate

x

y

z

\min (I_{i n})

and

\max (I_{i n})

are the minimum and maximum intensity of all CT images of each wheat spike

I_{i n}

I_{o u t}

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{matrix} L_{C_{i}} = G r a i n i f s_{C_{i}} > 10000 \\ L_{C_{i}} = S t e m N o d e i f s_{C_{i}} > 50 & & s_{C_{i}} < 2000 \\ L_{C_{i}} = O t h e r i f s_{C_{i}} < 50 \end{matrix}

(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_{m i n} > λ_{1} \cdot i_{m e a n} o r i_{m i d} > λ_{2} \cdot i_{m e a n}

(3)

where, for all the candidate voxel and its six neighboring voxels,

i_{m i n}

is the minimum intensity,

i_{m e a n}

is mean intensity, and

i_{m i d}

is the median intensity.

λ_{1}

and

λ_{2}

are adjustment coefficients, whose values may differ with crop types. In our experiments,

λ_{1}

and

λ_{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

μ_{1}

μ_{2}

, and

μ_{3}

(

μ_{1} \geq μ_{2} \geq μ_{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}

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

{N_{i}}

\hat{S} = \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 – 3}

associated with voxels in classes

S_{1 – 3}

. These weights are

W_{1} \approx 0.894

W_{2} \approx 1.3409

W_{3} \approx 1.5879

. The weight

W_{4 – 9}

associated with voxel classes

S_{4 – 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 – 3}

usually appear in some planes, types

S_{4 – 6}

often appear in curved border regions. Voxel types

S_{7 – 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 – 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:

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

(5)

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

(6)

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

(7)

Table 1 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 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.

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.

Author Contributions

Data set labelling, B.W. and C.L.; Funding acquisition, B.X., S.X. and X.Y.; Methodology, B.X., B.W. and C.L.; Supervision, S.X.; Validation, S.X.; Visualization, B.W.; Writing—original draft, B.X. and B.W.; Writing—review and editing, B.X. and B.W.

Funding

This work was in part supported by the National Key Research and Development Program of China (Grant No. 2016YFD0101900), the Defense Industrial Technology Development Program (Grant No. JCKY2018110C165), Hubei Provincial Natural Science Foundation of China (Grant No. 2017CFA012), and the Fundamental Research Funds for the Central Universities (WUT: 2019IVA043).

Acknowledgments

We would like to acknowledge the work and dataset of Nathan Hughes, etc. Without their inspiring work and generous dataset, our work is not possible to be done. We would also like to thank all the reviewers, whose valuable comments and suggestions greatly improve this article.

Conflicts of Interest

The authors declare no conflict of interest.

References

Glémin, S.; Bataillon, T. A comparative view of the evolution of grasses under domestication. New Phytol. 2009, 183, 273–290. [Google Scholar] [CrossRef]
Reynolds, M.; Foulkes, M.J.; Slafer, G.A.; Berry, P.; Parry, M.A.J.; Snape, J.W.; Angus, W.J. Raising yield potential in wheat. J. Exp. Bot. 2009, 60, 1899–1918. [Google Scholar] [CrossRef]
Tardieu, F.; Tuberosa, R. Dissection and modelling of abiotic stress tolerance in plants. Curr. Opin. Plant Biol. 2010, 13, 206–212. [Google Scholar] [CrossRef]
Fang, C.; Ma, Y.; Wu, S.; Liu, Z.; Wang, Z.; Yang, R.; Hu, G.; Zhou, Z.; Yu, H.; Zhang, M.; et al. Genome-wide association studies dissect the genetic networks underlying agronomical traits in soybean. Genome Biol. 2017, 18, 161. [Google Scholar] [CrossRef]
Gao, J.; Yang, S.; Cheng, W.; Fu, Y.; Leng, J.; Yuan, X.; Jiang, N.; Ma, J.; Feng, X. GmILPA1, Encoding an anaphase-promoting complex-like Protein, affects Leaf Petiole Angle. Plant Physiol. 2017, 174, 1167–1176. [Google Scholar] [CrossRef] [PubMed]
Lu, S.; Zhao, X.; Hu, Y.; Liu, S.; Nan, H.; Li, X.; Fang, C.; Cao, D.; Shi, X.; Kong, L.; et al. Natural variation at the soybean J locus improves adaptation to the tropics and enhances yield. Nat. Genet. 2017, 49, 773–779. [Google Scholar] [CrossRef] [PubMed]
Li, Z.; Hu, G.; Liu, X.; Zhou, Y.; Li, Y.; Zhang, X.; Yuan, X.; Zhang, Q.; Yang, D.; Wang, T.; et al. Transcriptome Sequencing Identified Genes and Gene Ontologies Associated with Early Freezing Tolerance in Maize. Front. Plant Sci. 2016, 7, 1–14. [Google Scholar] [CrossRef] [PubMed]
Zhao, X.; Cao, D.; Huang, Z.; Wang, J.; Lu, S.; Xu, Y.; Liu, B.; Kong, F.; Yuan, X. Dual functions of GmTOE4a in the regulation of photoperiod-mediated flowering and plant morphology in soybean. Plant Mol. Biol. 2015, 88, 343–355. [Google Scholar] [CrossRef]
Fiehn, O. Metabolomics–the link between genotypes and phenotypes. Plant Mol. Biol. 2002, 48, 155–171. [Google Scholar] [CrossRef] [PubMed]
Kirigwi, F.M.; Van Ginkel, M.; Brown-Guedira, G.; Gill, B.S.; Paulsen, G.M.; Fritz, A.K. Markers associated with a QTL for grain yield in wheat under drought. Mol. Breed. 2007, 20, 401–413. [Google Scholar] [CrossRef]
Granier, C.; Vile, D. Phenotyping and beyond: Modelling the relationships between traits. Curr. Opin. Plant Biol. 2014, 18, 96–102. [Google Scholar] [CrossRef]
Leilah, A.A.; Al-Khateeb, S.A. Statistical analysis of wheat yield under drought conditions. J. Arid Environ. 2005, 61, 483–496. [Google Scholar] [CrossRef]
Gegas, V.C.; Nazari, A.; Griffiths, S.; Simmonds, J.; Fish, L.; Orford, S.; Sayers, L.; Doonan, J.H.; Snape, J.W. A Genetic Framework for Grain Size and Shape Variation in Wheat. Plant Cell 2010, 22, 1046–1056. [Google Scholar] [CrossRef]
Hughes, N.; Askew, K.; Scotson, C.P.; Williams, K.; Sauze, C.; Corke, F.; Doonan, J.H.; Nibau, C. Non-Destructive, high-Content analysis of wheat grain traits using X-ray micro computed tomography. Plant Methods 2017, 13, 1–16. [Google Scholar] [CrossRef]
Furbank, R.T.; Tester, M. Phenomics—Technologies to relieve the phenotyping bottleneck. Trends Plant Sci. 2011, 16, 635–644. [Google Scholar] [CrossRef]
Fahlgren, N.; Gehan, M.A.; Baxter, I. Lights, camera, action: High-throughput plant phenotyping is ready for a close-up. Curr. Opin. Plant Biol. 2015, 24, 93–99. [Google Scholar] [CrossRef]
Wang, B.; Gao, Y. Fast and effective retrieval of plant leaf shapes. Lect. Notes Comput. Sci. 2013, 7725 LNCS, 475–486. [Google Scholar]
Xiong, B.; Jancosek, M.; Oude Elberink, S.; Vosselman, G. Flexible building primitives for 3D building modeling. ISPRS J. Photogramm. Remote Sens. 2015, 101, 275–290. [Google Scholar] [CrossRef]
Xiong, B.; Oude Elberink, S.; Vosselman, G. A graph edit dictionary for correcting errors in roof topology graphs reconstructed from point clouds. ISPRS J. Photogramm. Remote Sens. 2014, 93, 227–242. [Google Scholar] [CrossRef]
Xiong, B.; Liu, Q.; Xiong, J.; Li, S.; Wang, S.; Liang, D. Field-of-Experts Filters Guided Tensor Completion. IEEE Trans. Multimed. 2018, 20, 2316–2329. [Google Scholar] [CrossRef]
Wang, Y.; Wen, W.; Wu, S.; Wang, C.; Yu, Z.; Guo, X.; Zhao, C. Maize Plant Phenotyping: Comparing 3D Laser Scanning, Multi-View Stereo Reconstruction, and 3D Digitizing Estimates. Remote Sens. 2019, 11, 63. [Google Scholar] [CrossRef]
Borisjuk, L.; Rolletschek, H.; Neuberger, T. Surveying the plant’s world by magnetic resonance imaging. Plant J. 2012, 70, 129–146. [Google Scholar] [CrossRef]
Hubeau, M.; Steppe, K. Plant-PET Scans: In Vivo Mapping of Xylem and Phloem Functioning. Trends Plant Sci. 2015, 20, 676–685. [Google Scholar] [CrossRef]
Jahnke, S.; Menzel, M.I.; Van Dusschoten, D.; Roeb, G.W.; Bühler, J.; Minwuyelet, S.; Blümler, P.; Temperton, V.M.; Hombach, T.; Streun, M.; et al. Combined MRI-PET dissects dynamic changes in plant structures and functions. Plant J. 2009, 59, 634–644. [Google Scholar] [CrossRef] [PubMed]
Xu, X.; Liu, Y.; Liu, Q.; Lu, H.; Zhang, M. Gradient-Based Low Rank Method for Highly Undersampled Magnetic Resonance Imaging Reconstruction. J. Shanghai Jiaotong Univ. 2018, 23, 384–391. [Google Scholar] [CrossRef]
Wang, S.; Tan, S.; Gao, Y.; Liu, Q.; Ying, L.; Xiao, T.; Liu, Y.; Liu, X.; Zheng, H.; Liang, D. Learning Joint-Sparse Codes for Calibration-Free Parallel MR Imaging. IEEE Trans. Med. Imaging 2018, 37, 251–261. [Google Scholar] [CrossRef]
Rahaman, M.M.; Chen, D.; Gillani, Z.; Klukas, C.; Chen, M. Advanced phenotyping and phenotype data analysis for the study of plant growth and development. Front. Plant Sci. 2015, 6, 1–15. [Google Scholar] [CrossRef] [PubMed]
Aerts, H.J.W.L.; Velazquez, E.R.; Leijenaar, R.T.H.; Parmar, C.; Grossmann, P.; Cavalho, S.; Bussink, J.; Monshouwer, R.; Haibe-Kains, B.; Rietveld, D.; et al. Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach. Nat. Commun. 2014, 5, 4006. [Google Scholar] [CrossRef] [PubMed]
Yang, W.; Guo, Z.; Huang, C.; Duan, L.; Chen, G.; Jiang, N.; Fang, W.; Feng, H.; Xie, W.; Lian, X.; et al. Combining high-throughput phenotyping and genome-wide association studies to reveal natural genetic variation in rice. Nat. Commun. 2014, 5, 1–9. [Google Scholar] [CrossRef] [PubMed]
Neethirajan, S.; Karunakaran, C.; Jayas, D.S.; White, N.D.G. X-ray Computed Tomography Image Analysis to explain the Airflow Resistance Differences in Grain Bulks. Biosyst. Eng. 2006, 94, 545–555. [Google Scholar] [CrossRef]
Repesa, M.; Sofic, A.; Jakupovic, S.; Tosum, S.; Kazazic, L.; Dervisevic, A. Comparison of Results of Measurement of Dimensions of the Placed Dental Implants on Cone Beam Computed Tomography with Dimensions of the Producers of the Implants. Acta Inform. Medica 2017, 25, 116. [Google Scholar] [CrossRef]
McElrone, A.J.; Choat, B.; Parkinson, D.Y.; MacDowell, A.A.; Brodersen, C.R. Using High Resolution Computed Tomography to Visualize the Three Dimensional Structure and Function of Plant Vasculature. J. Vis. Exp. 2013, e50162. [Google Scholar] [CrossRef]
Gomez, F.E.; Carvalho, G.; Shi, F.; Muliana, A.H.; Rooney, W.L. High throughput phenotyping of morpho-anatomical stem properties using X-ray computed tomography in sorghum. Plant Methods 2018, 14, 59. [Google Scholar] [CrossRef]
Rousseau, D.; Widiez, T.; Di Tommaso, S.; Rositi, H.; Adrien, J.; Maire, E.; Langer, M.; Olivier, C.; Peyrin, F.; Rogowsky, P. Fast virtual histology using X-ray in-line phase tomography: Application to the 3D anatomy of maize developing seeds. Plant Methods 2015, 11, 55. [Google Scholar] [CrossRef]
Metzner, R.; Eggert, A.; van Dusschoten, D.; Pflugfelder, D.; Gerth, S.; Schurr, U.; Uhlmann, N.; Jahnke, S. Direct comparison of MRI and X-ray CT technologies for 3D imaging of root systems in soil: Potential and challenges for root trait quantification. Plant Methods 2015, 11, 17. [Google Scholar] [CrossRef]
Pardas, M.; Salembier, P.; Torres, L. 3D morphological segmentation for image sequence processing. In Proceedings of the IEEE Winter Workshop on Nonlinear Digital Signal Processing, Tampere, Finland, 17–20 January 1993. [Google Scholar] [CrossRef]
Maire, E.; Colombo, P.; Adrien, J.; Babout, L.; Biasetto, L. Characterization of the morphology of cellular ceramics by 3D image processing of X-ray tomography. J. Eur. Ceram. Soc. 2007, 27, 1973–1981. [Google Scholar] [CrossRef]
Xiong, B.; Zhang, X.; Jiang, W. Semi-Supervised Classification based on Gaussian Mixture Model for remote imagery. Sci. China Technol. Sci. 2010, 53, 85–90. [Google Scholar] [CrossRef]
Mangan, A.P.; Whitaker, R.T. Partitioning 3D surface meshes using watershed segmentation. IEEE Trans. Vis. Comput. Graph. 1999, 5, 308–321. [Google Scholar] [CrossRef]
Wold, S.; Esbensen, K.; Geladi, P. Principal component analysis. Chemom. Intell. Lab. Syst. 1987, 2, 37–52. [Google Scholar] [CrossRef]
Windreich, G.; Kiryati, N.; Lohmann, G. Voxel-based surface area estimation: From theory to practice. Pattern Recognit. 2003, 36, 2531–2541. [Google Scholar] [CrossRef]
Mullikin, J.C.; Verbeek, P.W. Surface area estimation of digitized planes. Bioimaging 1993, 1, 6–16. [Google Scholar] [CrossRef]
Marshalla, D.R.; Maresa, D.J.; Mossb, H.J.; Ellison, F.W. Effects of grain shape and size on milling yields in wheat. II. Experimental studies Effects of Grain Shape and Size on Milling Yields in Wheat. II *. Experimental Studies. Aust. J. Agric. Res. 1986, 37, 331–342. [Google Scholar] [CrossRef]
Karimi, M.; Kheiralipour, K.; Tabatabaeefar, A.; Khoubakht, G.M.; Naderi, M.; Heidarbeigi, K. The Effect of Moisture Content on Physical Properties of Wheat. Pakistan J. Nutr. 2009, 8, 90–95. [Google Scholar]

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. The whole pipeline is done automatically while some parameters are given manually.

Figure 2. A slice image sample and representations of wheat spike tissue.

Figure 3. Labeled grains images. The top row images are sets of raw slices after holder clearance processing, and the bottom three are the corresponding three labeled binary slices.

Figure 4. The pipeline of pro-processing of CT image. The parts are all automated.

Figure 5. Stack 2D image sequence into 3D image. (a) The original group of 2D images; (b) Three example images; (c) The 2D matrices of each image are stacked along the vertical direction. Their horizontal coordinates remain unchanged. The coordinate values in the vertical direction are their numbers in the image group; (d) The result of stacked 3D images.

Figure 6. Holder clearance processing. (a) The raw image; (b) Image after removing holder using Region of Interest (ROI).

Figure 7. Image intensity equalization. Top three are images of different wheat spikes with different intensities, while the bottom three are equalized images correspondingly.

Figure 8. CT image histogram of a wheat spike. (a) The complete histogram; (b) the histogram between 50 and 255 intensity.

Figure 9. Foreground detection by GMM. (a) The raw CT image; (b) Image after threshold segmentation.

Figure 10. The 3D morphological erosion for separating connected objects. Left one is before erosion processing, we can see the adhesion between the right two grains; Middle and right images are the grains after one and two steps of erosion. The top row are three 2D images and the bottom row are the 3D views, correspondingly. In the red circles, we can see the adhesion has been removed, and each connected component is one grain.

Figure 11. Detected grains and stem nodes. Stem nodes have similar intensity with grains, while stem nodes and grains are spatially close.

Figure 12. The 3D region growing method to recover lost voxels. (a) The raw CT image; (b) Detected grain region by erosion processing; (c) The yellow voxels are candidate voxels by the first iteration; (d) The green voxels are newly grown grains; (e) The newly grown voxels are set to seeds; (f) Candidate voxels (in yellow) and newly grown grain voxels (in green) by the second growing iteration; (g) Six-neighboring voxels in 3D view; (h) An example of an voxel on grain surface and its six-neighboring voxels.

Figure 13. Computation of length and radius of grains by PCA.

Figure 14. The nine unique surface voxel classes.

Figure 15. Stem detection and grain angle computation. (a) 3D images of a wheat spike; (b) Detected stem by 3D morphology; (c) Combined stem nodes and grains; (d) Computed grain angle.

Figure 16. The detection results of the top part of wheat spike with treatment No. 1178. (a) Grains detected by Hughes’ method; (b) Grains detected by our method; (c) Manually labeled grains. The difference between our and Hughes results is marked with a red box.

Figure 17. Spike treatment No. 1172 (upper one) and No. 1188 (lower one) contain one grain but some background too. The green grain is detected by our method, the red one is from Hughes. We can see that their grain is much bigger than actual value, which cause the precision loss.

Figure 18. The details of the grain detection results for wheat spike (the bottom part of No. 1178). (a) Detected grains by Hughes’ method, visualized in red; (b) Detected grains by our method, visualized in green; (c) the two results are stacked together for evaluation; (d) Zoomed-in picture of (c). We can see that tips of grain (lower one) and some other voxels (middle one) are not detected using Hughes’ method, some other tissues (upper one) are erroneously detected as grains.

Figure 19. The raw images of the three groups of wheat spikes and the counted grains. Some treatments contain 2 set of images, we all set for one wheat spike in the figure. The treatment No. 0116C contains 25 grains but Hughes counted 24; 01161 contains 23 grains but Hughes counted 22. The treatment No. 01179 has no grain both counted by us and Hughes, but one grain is counted using their method.

Figure 20. The proposed method failed to find one grain in treatment No. 1180 spike because the grain is very close to another grain.

Figure 21. Detected stem nodes and grains. (a) The whole wheat spike; (b) the detected stem nodes; (c) the detected grain and stem node. Upper row is spike No. 1170, middle row is bottom part of spike No. 1165 and bottom row is top part of spike No. 1178.

Figure 22. The detected grains for wheat spike No. 1170. The grain number is sorted with their gravity height in descending order. We can see that No. 5 grain has bigger angle, No. 11-14 grains have bigger volume and No. 13 grain has biggest surface area.

Figure 23. The volume distribution of grains in 43 groups of spikes; the black short line means average value of the spike. The volume distribution of all wheat spikes please refer to Additional file II: All Spike Volume.pdf.

Figure 24. Spike yields according to different temperature and water conditions. Temperature and water have a great impact on production. The yields at 25 °C is nearly twice that of 35 °C. The yields in low water conditions will be slightly higher than the one in high water.

Figure 25. (a–d) Single grain morphology under different temperature and water conditions.

Figure 26. Correlation analysis of individual volume of grain between the key traits, including Grain-to-Grain distance, grain angle and grain aspect ratio. One sub-graph only considers one environment condition. (a,c,e) consider the effect of the Grain-to-Grain distance, angle and aspect ratio on volume at different temperatures. And (b,d,f) consider the effect of the Grain-to-Grain distance, angle and aspect ratio on volume under different water conditions.

Figure 27. Correlation matrix of wheat phenotypic traits. The upper triangle is the graphical representation, and the lower triangle is the corresponding correlation coefficient.

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

Treatment No.	Grain Number	Our Method			Hughes’ Method
Treatment No.	Grain Number	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.

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

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)	(mm³)	(mm³)	(mm²)	(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

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Xiong, B.; Wang, B.; Xiong, S.; Lin, C.; Yuan, X. 3D Morphological Processing for Wheat Spike Phenotypes Using Computed Tomography Images. Remote Sens. 2019, 11, 1110. https://doi.org/10.3390/rs11091110

AMA Style

Xiong B, Wang B, Xiong S, Lin C, Yuan X. 3D Morphological Processing for Wheat Spike Phenotypes Using Computed Tomography Images. Remote Sensing. 2019; 11(9):1110. https://doi.org/10.3390/rs11091110

Chicago/Turabian Style

Xiong, Biao, Bo Wang, Shengwu Xiong, Chengde Lin, and Xiaohui Yuan. 2019. "3D Morphological Processing for Wheat Spike Phenotypes Using Computed Tomography Images" Remote Sensing 11, no. 9: 1110. https://doi.org/10.3390/rs11091110

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

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

Abstract

1. Introduction

1.1. Related Works

1.2. Contributions

2. Materials

2.1. Physical Equipment and Plant Materials

2.2. Manual Labeling of Wheat Grains

3. Method

3.1. Holder Clearance

3.2. Intensity Equalization

3.3. Foreground Detection by GMM

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

3.4.1. Detection of Grain and Stem Nodes

3.4.2. 3D Growth to Recover Grains

3.5. 3D Phenotypes Calculation

3.5.1. Grain Length and Radius

3.5.2. Grain Volume and Surface Area

3.5.3. Aspect Ratio and Porosity

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

4. Results and Discussion

4.1. Results of Grain Detection

4.2. Grain Number Counting

4.3. Stem Node Detection

4.4. Grain Size and Angle

4.5. Analysis of 3D Phenotypes

4.5.1. Volume Distribution

4.5.2. Factors Affecting Spike Yields

4.5.3. Factors Affecting Grain Phenotype

4.5.4. Correlation Between Volume and Shapes

5. Conclusions

Author Contributions

Funding

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI