Three-Dimensional Reconstruction of Soybean Canopies Using Multisource Imaging for Phenotyping Analysis

Abstract

Geometric three-dimensional (3D) reconstruction has emerged as a powerful tool for plant phenotyping and plant breeding. Although laser scanning is one of the most intensely used sensing techniques for 3D reconstruction projects, it still has many limitations, such as the high investment cost. To overcome such limitations, in the present study, a low-cost, novel, and efficient imaging system consisting of a red-green-blue (RGB) camera and a photonic mixer detector (PMD) was developed, and its usability for plant phenotyping was demonstrated via a 3D reconstruction of a soybean plant that contains color information. To reconstruct soybean canopies, a density-based spatial clustering of applications with noise (DBSCAN) algorithm was used to extract canopy information from the raw 3D point cloud. Principal component analysis (PCA) and iterative closest point (ICP) algorithms were then used to register the multisource images for the 3D reconstruction of a soybean plant from both the side and top views. We then assessed phenotypic traits such as plant height and the greenness index based on the deviations of test samples. The results showed that compared with manual measurements, the side view-based assessments yielded a determination coefficient (R²) of 0.9890 for the estimation of soybean height and a R² of 0.6059 for the estimation of soybean canopy greenness index; the top view-based assessment yielded a R² of 0.9936 for the estimation of soybean height and a R² of 0.8864 for the estimation of soybean canopy greenness. Together, the results indicated that an assembled 3D imaging device applying the algorithms developed in this study could be used as a reliable and robust platform for plant phenotyping, and potentially for automated and high-throughput applications under both natural light and indoor conditions.

1. Introduction

Soybeans are one of the main cash crops worldwide. To meet the needs of the growing human population, plant scientists and breeders must increase the productivity and yield of soybean crops, which is a substantial challenge [1]. High-throughput phenotyping platforms are essential for tracking the growth of soybean plants in the field and the contributions of these plants to both the food supply and the generation of bioenergy from their biomass [2]. Plant phenotyping involves the comprehensive assessment of plant characteristics that result from genetic and environmental factors [3,4,5,6]. Phenotypic characteristics include external morphological parameters such as plant height, size, petiole length, and initiation angle, as well as internal properties such as chlorophyll and nutrient (nitrogen, phosphorus, and potassium) contents [7,8]. The measurement and evaluation of these phenotypic traits can provide guidelines for plant breeding. However, most of the plant phenotyping methods depend on observations and manual measurements with contact sensors; these methods are considered low throughput, costly, and labor-intensive [9]. Furthermore, these traditional techniques typically require the destruction of plant organs, which negatively affects the normal growth of the measured plants [10]. Alternatively, plant phenotyping, which involves noninvasive measurements, is an emerging technology that has recently attracted the attention of researchers in the plant science and agricultural fields.

A current bottleneck that restricts phenotyping analysis lies in the development of the noninvasive and automated calculation methods that quantify the phenotypic traits of plants. In this regard, computer vision technologies show great potential for noninvasive measurement methods [11,12]. Currently, by recording whole three-dimensional (3D) datasets and reconstructing entire plants, 3D laser scanners can be used to acquire highly detailed plant phenotypic information on geometric traits, such as plant height. Photogrammetric information can be obtained via several different 3D techniques and laser devices, such as handheld laser scanners [13]. 3D plant models can then be generated by multiview stereo algorithms or by structure from motion [9,14]. Having very high resolutions, these techniques also present disadvantages, such as sensitivity to occlusion, a lack of color information, and failure to accurately reflect important phenotypic traits. Moreover, for the automatic phenotypic analysis of plant structure, 3D point clouds generated by laser scanners must be properly extracted from a large amount of 3D data, and must be classified. In addition, the high cost and limited availability of laser scanning devices have prevented their widespread application.

In addition to geometric traits, spectral reflectance in visible spectrum is another important characteristic for phenotyping analysis, especially for diagnosing nutrient conditions. Significant correlations between “soil and plant analyzer development (SPAD)” readings, leaf N concentration (LNC), and image color indices in RGB color space were observed for rice in natural light conditions [15]. Additionally, the strong relationship between the deficiency of nitrogen, phosphorus, potassium, and image color indices in RGB color space was further proved for soybean plants in an outdoor environment [16]. Although light intensity is one of the major factors that lead to color distortion, this kind of distortion can be registered by reference-based approaches to a great degree [17].

The above-mentioned studies showed that the nutrient conditions of plants of different cultivars and at different growth stages can be estimated in an effective way, that is, by using R, G, and B channel values that can be obtained by low-cost RGB cameras, without a wet lab, field-based sensors, or expensive image acquisition equipment.

Currently, geometric and spectral characteristics are acquired separately by different devices. With respect to agricultural applications, the spectral reflectance of plants’ canopies can be obtained by using an RGB-based camera or a multispectral camera [18,19], whereas the geometric information therein, such as plant height and canopy width, can be calculated via 3D cloud points generated by using 3D imaging devices such as terrestrial laser scanners or handheld laser scanners. Although some devices such as Kinect sensors [20] can capture image information (i.e., RGB images) and 3D distance information simultaneously, these two types of information cannot be acquired within the same coordinate system. Consequently, to obtain 3D plants with color information, a coordinate transformation is further required to be performed on these two types of information, thereby fusing multisource images. However, few techniques integrated with image fusion have been developed for the phenotyping analysis of plants, especially for soybean plants. Thus, there is still an urgent need in 3D reconstruction for plant canopies with both geometric and color information.

Although soybean plants are globally important crops for providing oils, most studies on reconstruction and phenotyping have focused on a single plant in a laboratory circumstance, such as a plant grown in a greenhouse [21]. These laboratory-based systems can be useful for obtaining certain phenotypic traits such as plant height, but they also have limitations. The main shortcoming is that certain traits calculated from an individual plant can not reflect the true cases in the field due to the artificial laboratory environment, which can significantly affect plant growth. Thus, phenotyping methods under natural light conditions are highly desirable for plant breeding. However, in contrast with progress made in laboratory-based research [22], the lack of a well-constructed platform has hindered the development of plant phenotyping under natural light conditions in large plantations. Therefore, there is a great opportunity for innovation with respect to the development of 3D-based technology for plant phenotyping under natural growing conditions.

To address these issues, it is demonstrated in this paper that it is feasible to reconstruct the 3D geometric characteristics containing the color information of soybean plants under natural light conditions by using a relatively low-cost multisource imaging system. This reconstruction was carried out through the efficient extraction and registration of multisource images by using density-based spatial clustering of applications with noise (DBSCAN), principal component analysis (PCA), and iterative closest point (ICP) algorithms. Our approach resulted in the successful characterization of soybean plant phenotypic traits (plant height and greenness index) acquired from 3D reconstructed images. This technique provides an alternative for the large-scale characterization of plant phenotypes under natural conditions.

2. Materials and Methods

2.1. Overall Process Flow for 3D Reconstruction

The overall processes and methods for performing a 3D reconstruction of the soybean canopy and the calculation of phenotypic traits are shown in Figure 1. First, raw 3D point cloud data including the soybean canopy information were captured from the side view (Figure 1a) and top view (Figure 1b), respectively, by using a multisource imaging system consisting of a photonic mixer detector (PMD; model: Camcube 3.0, PMDTech Company, Siegen, Germany) and an RGB camera (model: C270, Logitech Company, Lausanne, Switzerland). Second, the soybean canopies (Figure 1d,f) were extracted from the raw 3D cloud point data by executing the DBSCAN algorithm. Third, the RGB image and associated 3D point cloud information were fused together to reconstruct the geometric morphology of the soybean canopies (Figure 1g,h) by executing the PCA algorithm for a rough registration and the ICP algorithm for optimal registration between the input point cloud and reference point cloud. Last, phenotypic traits such as plant height and the greenness index were calculated. The phenotypic data of the soybean canopies were compared with the data collected by manual measurements, such that the accuracy of the developed system can be evaluated.

2.2. Experimental Treatments and Measurement of Phenotypic Traits

The plantation experiment was conducted in 2016 at the Heilongjiang Bayi Agricultural University. A total of 70 soybean plants, including four varieties (namely, Kennong23, Kennong29, Kennong30, and Kennong33), were cultivated in pots. The experiment was conducted in accordance with a randomized complete block design. Each variety was replicated 17 times and planted in a polyvinyl chloride (PVC) pot (25-cm diameter and 40-cm height) prior to disinfection and germination treatments of soybean seeds. Twenty kilograms of soil and sand (2:1 w/w) were mixed together and added to each pot. The nitrogen, phosphorus, and potassium nutrient compositions were 50 mg/kg, 30 mg/kg, and 30 mg/kg, respectively. One, two, or three soybean plants were grown in an individual pot. Plant height was measured with a ruler.

2.3. Multisource Imaging System

The multisource imaging system (Figure 2) was composed of a PMD camera and an RGB camera. The PMD is a developer of complementary metal oxide semiconductor (CMOS) 3D time-of-flight (TOF) components, and a provider of engineering support in the field of digital 3D imaging. With its active imaging system, this depth camera can irradiate an active light source onto objects. The light reflected from the objects is then used to generate the depth image by calculating the TOF between the emission and reception [11]. In addition to generating depth images, the PMD camera was used to construct three additional multisource images: an intensity image, an amplitude image, and a flag image. The intensity image was recorded to illustrate the average intensity of incident light, natural light, and near-infrared light mixed together, and the amplitude image was used to indicate the reflecting ability of the objects. The flag image reflected the quality of the image pixels. Although this imaging system has the powerful capability of acquiring multisource images, there is no color information therein. This deficiency can be overcome by using the spectral reflectance acquired from an RGB camera at a spatial resolution of 640 × 480 pixels. Moreover, the PMD camera can be operated at a high frame rate (40 frames/s), but a resolution thereof (200 × 200 pixels) is relatively low.

Figure 3 shows examples of the multisource images, including an RGB image, a distance image, an amplitude image, a flag image, and an intensity image. The multisource images played an important role in the 3D reconstruction of soybean canopies, and the RGB images provided rich color information of the 3D canopies. The distance images were used to generate accurate 3D distance information, including the information with respect to x, y, and z-axes. Additionally, the amplitude images aided in the removal of the background of the RGB images. The flag images could be used to verify the quality of the multisource images when the valid images were selected. Last, the intensity images and RGB images were used together to calibrate the binocular cameras via the manual selection of feature points. Therefore, the combination of these two cameras is advantageous for the 3D reconstruction of soybean canopies with color information therein. Figure 2 shows an integration of the RGB camera and the PMD camera. In this way, these two cameras (integrated as a multisource imaging system) could collect images simultaneously via the software preinstalled thereon.

2.4. Calibration of Multisource Imaging System

The multisource imaging system was a binocular vision system composed of a PMD camera and an RGB camera (Figure 4). The camera calibration toolbox of MATLAB (version: 8.2) was used to calibrate both cameras. The position of any point P in the space could be approximated as a pinhole imaging model. There were four coordinate systems (Figure 4): the camera coordinate system (OcL-X_LY_LZ_L and OcR-X_RY_RZ_R), the imaging plane coordinate system (OL-X_LY_L and O_R-X_RY_R), the image coordinate system (O1-uLvL and O2-uRvR), and the world coordinate system.

The coordinates of P were (X₁, Y₁, Z₁), (x₁, y₁), and (u₁, v₁) in the PMD camera’s coordinate system, imaging plane coordinate system, and image coordinate system, respectively. The coordinates of P were (X₂, Y₂, Z₂), (x₂, y₂), and (u₂, v₂) in the RGB camera’s coordinate system, imaging plane coordinate system, and image coordinate system, respectively.

The relationship between the PMD and RGB coordinate systems is expressed by Equation (1):

$$\left[\begin{array}{l}{X}_2\\ Y_2\\ Z_2\end{array}\right]=R\cdot \left[\begin{array}{l}{X}_1\\ Y_1\\ Z_1\end{array}\right]+T$$

(1)

where R is a 3 × 3 orthogonal matrix and T is a translation vector, which was the external parameter determined by the position of the cameras.

R and T are solved by using the calibration method provided by Zhou [23], and the final results of R and T are as follows:

$$R=\left[\begin{array}{ccc}1.0000& 0.0064& 0.0023\\ -0.0061& 0.9930& -0.1178\\ -0.0031& 0.1178& 0.9930\end{array}\right]T=\left[\begin{array}{c}-22.89743\\ 290.58127\\ -22.86262\end{array}\right]$$

The imaging system (the PMD and RGB cameras) achieved the optimal calibration effect via the optical calibration method mentioned above. However, due to the complexity of soybean canopies, we further calibrated the PMD and RGB cameras for fusion by using the registration method of image feature points based on the optical calibration method, so as to accurately reconstruct 3D soybean plants (Section 2.6.2).

2.5. Data Collection and RGB Image Preprocessing

The images were acquired by using the multisource imaging system between 10 am and 12 pm between 4–8 July 2016. The images were captured under natural environmental conditions at 28 °C and a light intensity of 1000 lux, and then stored in the laptop computer in a native raw format. Two shooting methods—namely, shooting from side view and shooting from top view—were compared based on an assessment of their ability to reflect plant height and greenness index. In addition, the effectiveness of 3D reconstruction by both methods was evaluated. Regardless of which method was used, when the images were collected, the whole plant or canopy had to be fitted within the field of view (FOV) of the cameras.

Since the imaging effect of RGB images is easy to be affected by natural light, thus, the values of the R, G, and B channels acquired by a RGB camera need to be corrected. A calibration chart was used for correction so as to obtain RGB images that have uniform spectral reflectance [17,24].

2.6. 3D Reconstruction

A three-step framework was proposed to reconstruct an individual soybean plant in the side view and top view. We compared the two data acquisition methods in terms of their ability to calculate two phenotypic traits. The comparison results were used to further evaluate the accuracy of the 3D reconstruction.

2.6.1. DBSCAN Algorithm for Point Cloud Filtering

We started with a raw point cloud without considering geometric and color information. Thus, the first important step in 3D reconstruction was to classify the useful points from the raw point clouds. We proposed extraction tasks with the DBSCAN [25] algorithm for side-view and top-view image acquisition methods (Figure 5). The information embedded in the raw data required for 3D reconstruction was not only reconstituted from the soybean plant, but also related to the background (building, sky, or people) within seven meters.

In this study, P was one set of 3D points in the raw data, and could be written as follows:

$$\mathrm{P}=\left\{P_1,P_2,P_3,\cdots ,P_i|1\le i\le n\right\}$$

(2)

where $P_i=\left\{\left(x,y,z\right)|x\in \mathbf{r},y\in \mathbf{r},z\in \mathbf{r}\right\}$, $\left(x,y,z\right)$ are the coordinates of $P_i$, and r represents the real number set. NH represents the neighborhood of $P_i$, which belongs to the following aggregate:

$$\mathrm{NH}\left(P_i\right)=\left\{P_j|\sqrt{\left(P_i-P_j\right){\left(P_i-P_j\right)}^T}\le NH\right\}$$

(3)

To apply DBSCAN clustering for the extraction of plant information, the raw points were divided into core points, (density) reachable points, and outliers, and are defined as follows: a core point $P_i$ must be within distance k (k was the maximum radius of the neighborhood from $P_i$); a point $Q$ is reachable from $P_i$ if there is a path T₁, …, T_n, with T₁ = $P_i$ and T_n = Q, whereas each T_i + 1 is directly reachable from T_i; the outliers are the points that are not reachable from any other point. If $P_i$ was a core point, then it formed a cluster together with all of the points (core or non-core) that were reachable from it. Each cluster contained at least one core point; non-core points could be part of a cluster, but they formed its “edge” because they could not be used to reach additional points.

In Figure 6, minPts = 4. The point $P_i$ and the other blue points were core points, as the area surrounding these points in an NH radius contained at least four points (including the point itself); since they were all reachable from one another, they formed a single cluster. Points N and Q were not core points, but were reachable from P_i (via other core points); thus, points N and Q belonged to the cluster as well. Point M is a noise point that is neither a core point nor directly reachable.

For running DBSCAN, in addition to parameter k, a minimum number of points was needed to build a dense region, starting with a random point that had not been accessed before. A k-neighborhood of this starting point was then retrieved. In this research, if it contained at least 30 points within 0.06 m, a new cluster was generated. Otherwise, the point was marked as a noise point. If a point was a dense part of a cluster, then its k-neighborhood was also classified as that cluster. Hence, all of the points in the k-neighborhood were added to the cluster. This process continued until the density-connected cluster was thoroughly retrieved. All of the unvisited points were then traversed and processed for the purpose of detecting an additional cluster or noise point. After the DBSCAN algorithm was run, at least 32,690 points were filtered from the raw point cloud consisting of 200 × 200 points. Thus, 7310 effective points were extracted to form the 3D shape of a soybean plant.

2.6.2. Fusion of Multisource Images

The fusion of multisource images is a recent development and remains an active research field [26]. A major fusion problem concerning these multisource sensors resulted from the merging of a low-spatial resolution distance image and a high-spatial resolution RGB image. During this process, the PMD amplitude image and its corresponding RGB image were fused together to construct a 3D soybean plant containing color information. The whole procedure of this fusion involved coordinate selection, affine transformation, and the assignment of color information for 3D points.

First, coordinate transformation was implemented based on the PMD camera coordinate system to generate two different coordinate systems that had uniform structures. A control point selection tool called cpselect (Figure 7) in the MATLAB (version: 8.2) environment was used to select three pairs of control points between the PMD amplitude image and its related RGB image, after which the coordinates of the control points were returned to a cpstruct structure.

After obtaining the cpstruct structure, the coordinate system of the RGB camera was transformed to that of the PMD camera via affine transformation, which is a linear mapping method that strongly preserves shapes for points, straight lines, and planes in an RGB image. Furthermore, to assign color information to its corresponding 3D point, a color index matrix was defined, which stored the corresponding relations between the coordinates of 3D points and the values of R, G, and B channels. After the images were fused, a new distance image with color information was constructed in a 3D coordinate system. In the resultant RGB image, the size was changed to 200 × 200 pixels, which was a one-to-one match with its original distance image according to the color index matrix.

2.6.3. Registration of 3D Point Clouds between Front and Back Sides

The whole canopy image could be obtained simultaneously by cameras in the top view. As such, the multisource images acquired in the top view were not subjected to image registration between two sides (front and back), while this process was needed for those images collected in the side view. In this study, to show the 3D reconstruction of a soybean plant in the side view, the images acquired from both sides (front and back) of the plant were processed by registering the 3D point clouds.

In this research, $C^k$ was the aggregate of the point clouds, $C^k=\left\{x_i^k|i=1,2,\dots ,m\right\}$; $\overline{x}$ represents the center of a point cloud, $\overline{x}=\frac{1}{m}{\displaystyle \sum _{i=1}^m(x_i)}$; the variance was ${\sigma }_x^2=\frac{1}{m}{\displaystyle \sum _{i=1}^m{(x_i-\overline{x})}^2}$; and the covariance matrix was $Co{v}_x^{}=\frac{1}{m}{\displaystyle \sum _{i=1}^m(x_i-\overline{x}){(x_i-\overline{x})}^T}$.

For the input point cloud $x_i$ and the reference point cloud $y_i$, the purpose of the registration was to determine the optimum similarity transformation $T=\left[sR|t\right]$ and apply this transformation on $x_i$, which transformed the coordinate of $x_i$ to the counterpart of $y_i$.

$$y_i=sR{x}_i+t$$

(4)

where R, t, and s are the rotation matrix, translation vector, and scale factor, respectively.

The ICP, which is an algorithm [27] that is capable of maximally reducing the difference between the two point clouds, was used to perform the two abovementioned processes of point cloud registration. Many variants of ICP, which has been widely used for the geometric alignment of 3D surfaces generated by various scan methods, have been developed. The two most popular variants are point-to-point and point-to-projection algorithms. The former performs better at registration tasks primarily because the PCA algorithm [28] is applied in advance to obtain a rough registration of the two components of the point clouds. This work was achieved by the combination of the PCA and ICP algorithms. R and t of Equation (4) could ultimately be solved by the singular value decomposition (SVD) of the ICP algorithm.

(1) PCA for rough registration

The PCA algorithm was applied to determine the principal axis direction between $x_i$ and $y_i$, and the rough registration results were used as inputs for the ICP algorithm for exact registration by the ICP algorithm. The processes were as follows:

First, the point cloud $x_i$ was decomposed by SVD, from which the left singular vector U, eigenvalue matrix D, and right singular vector V could be obtained at the same time. The SVD operation for $x_i$ was calculated by Equation (2):

$$UD{V}^T=svd(Co{v}_x)$$

(5)

Second, the PCA coordinate was established, and the direction of axis was determined by $U$ and $V$. Afterward, a similarity transformation $T_{PCA}$ was calculated by Equation (6), during which $x_i$ was transformed from the original coordinate to the PCA coordinate.

$$T_{PCA}=U\left[I|-\overline{x}\right]$$

(6)

There were four candidate transformations because of the uncertainty of the principal axis direction. Thus, the optimal transformation must be selected from four transformations, which could minimize the registration error and maximize the normal vector consistency for $x_i$ and $y_i$.

The PCA algorithm ultimately achieved basic registration between the input and reference point clouds via rotation, translation, and zoom operations. The basic registration result was subsequently used as an input for the ICP algorithm.

(2) ICP for optimal registration

We defined P and Q as having the same points, although they were from two point cloud aggregates ($x_i$ and $y_i$, respectively). The coordinates of P and Q were $(x_1,y_1,z_1)$ and $(x_2,y_2,z_2)$; any point in aggregate $y_i$ could be expressed by $T_k$. Thus, the Euclidean distance between P and $Q_i$ was then calculated by the following:

$$\mathrm{d}\left(P,Q_i\right)=||\mathrm{P}-{\mathrm{Q}}_i||=\sqrt{{\left(x_1-x_{2i}\right)}^2+{\left(y_1-y_{2i}\right)}^2+{\left(z_1-z_{2i}\right)}^2}$$

(7)

Essentially, the steps of the ICP algorithm for optimal registration are described as follows:

First, we obtained the source point cloud $x_i$ and the reference point cloud $y_i$.

Second, we searched the nearest point $B_k$ in input point cloud $x_i$ for any point $A_k$ in reference point cloud $y_i$ using Formula (7), which could best align each source point to its match after weighting and rejecting outlier points.

Next, we solved the rotation transform matrix $R$, the translation vector t and the deviation $s$ via SVD in accordance with Equations (9) to (15). In addition, Equation (8) served as the objective function with the minimum square error [29].

$$e^2(R,t,s)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\Vert y_i-(sR{x}_i+t)\Vert }^2}$$

(8)

$$Co{v}_{yx}=\frac{1}{n}{\displaystyle \sum _i(y-\overline{y}){(x-\overline{x})}^T}$$

(9)

$$UD{V}^T=svd(Co{v}_x)$$

(10)

$$S=\{\begin{array}{ll}I,& \left|U\right|\left|V\right|=1\\ diag(1,\cdots 1,-1),& \left|U\right|\left|V\right|=1\end{array}$$

(11)

$$R=US{V}^T$$

(12)

$$t=\overline{y}-sR\overline{x}$$

(13)

$$s=\frac{1}{{\sigma }_X{}^2}tr(DS)$$

(14)

$$T_{icp}=\left[sR|t\right]$$

(15)

Further, we transformed the source points $x_i$ using the obtained transformation in the previous step to obtain $A_{k+1}$.

We then analyzed whether the deviation $D_k$ was convergent or not. If it was, the final matrix of the coordinate transformation and $A_{k+1}$ were determined; otherwise, it was iterated from step 2 until $D_k$ was a convergence. The similarity transformation $T_{icp}$ could be obtained after the convergence of $S$.

The final transformation matrix T from $x_i$ to $y_i$ was as follows:

$$T=T_{pca}^{-1}(ref)T_{icp}(input,ref)T_{pca}(input)$$

(16)

When the algorithms above were integrated together, the 3D shape of the soybean plant was reconstructed by registering the two point clouds acquired from the front and back. Traditionally, missing data due to occlusions and misarranged 3D point positions are calibrated by hole-filling algorithms according to the curvature of the surrounding triangle mesh [30]. Although an entire 3D model could be built in accordance with this method, it was not applicable for reconstruction with color information primarily because the values of the R, G, and B channels of the point generated by the hole-filling algorithms could not represent the real spectral reflectance of soybean canopies. Therefore, the next step for phenotyping analysis was based on the 3D plant model built without hole-filing algorithms applied.

2.7. Methods of Calculating 3D Phenotypic Traits

Referring to the current 3D reconstruction approaches for soybean plants, we used a mathematical method to describe plant organs, aiming at simplifying data processing and providing a repeatable and objective parameterization of growth processes. We simulated relevant plant phenotypic traits such as plant height and canopy greenness index in the side view and top view separately. These traits are important parameters for evaluating plant quality, as they play an important role in the entire growth stage via photosynthesis. A total of 70 soybean plants (35 plants for the side view and 35 plants for the top view) were selected from the field for measuring these two traits.

2.7.1. Method of Calculating Plant Height

Plant height is defined as the shortest distance from the upper boundary of the main photosynthetic tissues (excluding inflorescences) on a plant to the ground level [31,32], expressed in meters. Among various phenotypes, plant height is an important geometric parameter; it can be used not only as an indicator of whole soybean growth, but also to quantify other advanced parameters such as yield and total biomass [33,34]. Soybean plant height is a particularly important factor used by breeders for screening and selecting improved varieties.

Plant height is traditionally measured manually with rulers or handheld devices such as laser rangefinders. However, these methods are time-consuming and labor-intensive, and thus are not applicable for large-scale phenotyping analysis. Advances in 3D imaging techniques allow the measurement of geometric traits such as plant height via the acquisition of accurate and efficient data concerning the 3D structures of plants.

For top-view depth images, the camera was installed on a tripod that had an adjustable height, and each plant was placed on the ground in the center of the camera’s FOV. According to the TOF principle, the pixel intensity reflected the distance between objects and the lens of the camera in millimeters. A pixel would be invalid if its intensity was zero because of no reflection from an object. According to Figure 5(right), soybean plant height was calculated via Equation (17) in the 3D space of the soybean reconstruction model.

$$H=H_c-h-min(D_p)p\in P\mathrm{and}{D}_p>0$$

(17)

where $H$ represents the height of the soybean plant in millimeters and $H_c$ is the distance from the lens of the camera to the ground level in millimeters; $H_c$ was 2000 mm when the images were captured via the top view-based method. To capture the whole canopy image, an adjusted value was used according to the actual plant height. In addition, in this equation, $h$, $P$, $p$, and $D$ indicate the pot height shown in the image, the pixel set in a distance image, one pixel in the pixel set $P$, and the distance value of pixel $p$, respectively. Consequently, a total of 35 plant height measurements were made for calculating the relative error between the 3D measurements and manual measurements.

Moreover, for the side-view depth images, plant height was derived as follows:

$$H=Y_{max}-h$$

(18)

where $H$ and $h$ have the same meaning as above, and $Y_{max}$ is the maximum value of the coordinate $Y$ of the depth image. Further, the height of 35 plants was subjected to correlation analysis with the manual measurements.

2.7.2. Method of Calculating Greenness Index

The greenness index is a primary phenotypic trait concerning the color characteristic of the leaf, which can be used for the classification of plant health status [35] and the inversion of both the nitrogen content and chlorophyll content [19]. Conventionally, visual inspection and chemical analysis have been the major methods used for the evaluation of plant health. Methods involving leaf color charts (LCCs) are widespread among farmers, but much subjectivity is involved in the results. The chemical analysis-based approach is destructive and not conducive to continuous plant growth. The image color indices in RGB color space can serve as a more objective yet nondestructive method for plant health assessment [36].

In this study, the values of the R, G, and B channels of each point cloud were extracted from the 3D plant canopy. The invalid points and object colors were omitted from the calculation according to the distance value of the objects if they were not derived from the canopy. All of the points of the 3D point cloud were traversed to extract their values of the R, G, and B channels for the purpose of calculating the greenness index.

An original input 3D image with color information was investigated in the RGB color space. The imaging effect therein was subjected to a white balance treatment [17] that was processed in its corresponding 2D image prior to the RGB decomposition in which the R, G, and B channels of the image were extracted for determining the greenness index. We applied the following scheme to obtain the greenness index, which was used to quantify the relative health of the leaves [16]:

$$GS=B_{xyz}/(R_{xyz}+G_{xyz}+B_{xyz})$$

(19)

where GS is the greenness index and $R_{xyz},G_{xyz}$ and $B_{xyz}$ represented the values of the $R,G$, and $B$ channels at coordinates $\left(x,y,z\right)$ respectively. Consequently, 35 values of canopy greenness from the side view and 35 values of canopy greenness from the top view were compared with those obtained manually. The manual measurement values of canopy greenness index were calculated via the following method: all of the 2D images of the soybean canopies acquired by the RGB camera of the multisource imaging system were extracted from their complex background, and the average values of the R, G, and B channels of the whole canopy image were calculated separately by Adobe Photoshop software (version: CS6). The manually measured greenness values of greenness were determined according to Equation (19).

3. Results

3.1. 3D Reconstruction

For illustration purposes in this paper, a single plant and multiple plants were selected as typical examples for showing the reconstruction results. Soybean plants studied via two different acquisition methods (side view and top view) can be distinguished by the color of their pots. During data collection, the multisource imaging system was placed 80 cm from the target potted soybean plants; the positions were not changed. Thus, each pot was rotated 180 degrees around the center of the pot to acquire the back images after obtaining the front images of a soybean plant.

For the side view, the target 3D points (Figure 8d,e) from the front and back sides were extracted from the complex background (Figure 8c) via the DBSCAN algorithm mentioned in Section 2.6.1, and were then fused with color information according to the corresponding control points between the RGB image and amplitude image. Figure 8f,g illustrates the rough registration and exact registration results via PCA and ICP algorithms, respectively. Final reconstruction results for a single plant are shown in Figure 8h, which demonstrates that the entire shape with color information was reconstructed well, although some scattered 3D points were not fully clustered. Moreover, the accuracy of correspondence between RGB and 3D coordinates needs to be improved further.

Although many more phenotypic traits such as stem diameter and leaf angle can be researched via the 3D reconstruction of a soybean plant in the side view compared with the top view, relatively complex and time-consuming algorithms such as DBSCAN and registration algorithms were used to reconstruct the 3D shape of soybean canopies. Thus, for the two phenotypic traits (plant height and greenness index) in this paper, the much simpler top-view acquisition method was used to extract information on plant height and greenness index; this information was not subjected to image registration between the front and back of the soybean canopies. In addition, for the purpose of avoiding the effects of wind, the plants were moved indoors for acquiring images from the top view. Similarly, the canopy images were segmented from the raw data (Figure 9b) based on the DBSCAN algorithm combined with distance thresholding (TH = 80 mm) to remove any miscellaneous points. Notably, the thresholding needed to be adjusted according to the height of the plants. The canopy depth and plant height indicated by color bar in Figure 9c were calculated by Equation (17), and the 3D canopy was ultimately reconstructed by the fusion of the distance image and the two-dimensional (2D) RGB image. As shown in Figure 9d, the pixels at the edge of the leaves were missing because of the weak reflection at the edges. In addition, some useless red pixels, which were part of the pot, were observed in the final reconstruction. This phenomenon occurred mainly because a few leaves covered some sections of the pot, which affected the threshold setting. However, these red pixels did not affect the greenness calculations, and did not negatively impact the plant height measurements.

The resolution of the 3D reconstruction of soybean plants is 200 × 200 pixels not only for the side view, but also for the top view. Each pixel pitch is 45 μm, which depended on the resolution of the PMD camera.

3.2. Accuracy of Plant Height Measurements in the Side and Top Views

The plant height data obtained from the side view-based and top view-based 3D reconstructions were compared with the manually measured data (Figure 10). The plant height was measured via the algorithm described in Section 2.7.1. The linear best-fit models comparing the side view-based and top view-based plant height measurements with the manual measurements yielded determination coefficient (R²) values of 0.9890 and 0.9936, respectively. However, there were some biases: the side view-based plant heights fluctuated to some degree (by −1.8 cm to 1.7 cm), and the top view-based errors in the plant height calculations ranged from −1.1 cm to 1.5 cm. The average error for the side view and top view was 0.6713 cm and 0.2600 cm, respectively.

The 3D measurements indicated that the multisource imaging system could measure soybean plant height with a high degree of accuracy under both natural light conditions and indoor conditions. These findings provided useful guidelines for soybean plant height measurements under field conditions. Plant height phenotypes could be applied to guide the screening of soybean genotypes. For example, plant height is an indicator of early maturing soybean cultivars, which could reduce yield losses due to pests and diseases. Additionally, plant height can also be used to predict leaf area index and yield.

3.3. Accuracy of Greenness in the Side and Top Views

A crucial focus of this study was to produce a practical method for determining the greenness index, which provides valuable data to crop experts for diagnosing the diseases and nutrient conditions of soybean plants. The greenness index was calculated using the procedure described in Section 2.6.2. Figure 11 shows a comparison of the results of the 3D measurements and manual measurements. The greenness index was highly correlated with that assessed manually (R² of 0.8864 and an average error of 0.0117) for top view-based data; the minimum and maximum deviations were −0.03 and 0.03, respectively. The side view-based measurements yielded a correlation of R² = 0.6059, with an average error of 0.0386, and the deviation of calculation fluctuated between −0.14 and 0.07.

Although the top view-based plant height and greenness measurements were more accurate than their side view-based counterparts, other phenotypic traits such as leaf angle, branch angle, and stem diameter could not be measured effectively from the top view. Consequently, we performed 3D soybean reconstructions in the side view. The above-mentioned phenotypic traits will be considered in our future research. Furthermore, side-view 3D reconstruction will be improved by considering weather factors during data collection.

4. Discussion

4.1. Analysis of Experimental Results

In this paper, algorithms of the 3D reconstruction of soybean canopies were proposed for use in phenotyping analysis based on two acquisition methods (side view and top view). Additionally, plant height and greenness index were calculated.

In terms of accuracy, there was no significant difference in plant height between the side view and top view; the R² values from the side view and top view were 0.9890 and 0.9936, respectively. However, the side view-based greenness was less accurate than its top-view counterpart because of the random environmental factors under the natural light conditions affecting the 3D reconstructions of soybean canopies. The R² values from the side view and top view were 0.6059 and 0.8864, respectively.

Although the present experimental results meet the requirements for calculating the plant height and greenness index of soybean canopies, some aspects should be considered for obtaining more accurate greenness index data in the side view under natural conditions. The first aspect involves environmental factors, especially weather conditions. The ideal weather conditions involve a sunny day with no wind. However, the actual environmental conditions, such as high winds or rainy weather, cannot be controlled, but can be avoided. The second aspect involves the algorithms; both the DBSCAN and registration algorithms should be improved to acquire the optimal 3D reconstruction results of soybean canopies. The use of a modified DBSCAN algorithm [37] to cluster a 3D point cloud for extracting soybean canopies from complex backgrounds will be studied. In addition, compared with classic algorithms, improved ICP algorithms [38] can be used for more accurate and efficient exact registration between two pieces of 3D point clouds.

4.2. Evaluation of Algorithm Robustness

Robustness is a key indicator in algorithm evaluation. The algorithms in this research met the requirements needed for the 3D reconstruction and phenotyping analysis of soybean plants. The robustness of the algorithms could be evaluated from the aspects described below.

First, it is well known that light intensity will impact imaging. Therefore, to prove the applicability of the algorithms for 3D reconstruction, the multisource images were acquired both indoors and under natural conditions. Image acquisition under natural light was performed between 10 am and 12 pm, and the maximum light intensity reached 1000 lux. Under these circumstances, the algorithms worked well for the 3D reconstruction of soybean canopies.

Moreover, wind is another factor that affects 3D reconstruction. The presence of wind during data collection might result in the vibration or color distortion of the canopy, leading to inaccurate fusion between the 3D points and the RGB image. Some results proved that there was obvious image distortion at the edges and thickness of leaf organs when the wind speed varied from 0.9 m/s to 2.4 m/s. Therefore, to improve the stability of the algorithms, weather conditions with no or low wind (less than 2.4 m/s) are the best environments in which to obtain optimal 3D reconstructions of soybean canopies for phenotyping analysis [39,40].

In addition, background complexity is the third factor that affects imaging. During image acquisition, the background consisted of buildings, plants, sky, and ground, which could be considered a complex environment (Figure 3a). Furthermore, we acquired images under the conditions indicated above (sunny day, gentle breeze, and complex background).

Our correlation results yielded an R² value of 0.9890 for plant height and an R² value of 0.6059 for the greenness index under natural light conditions, both of which could be used to evaluate 3D reconstruction results. In addition, these experimental results, which were based on existing soybean samples, were accurate. Thus, the algorithms proposed in this paper exhibited excellent robustness for plant phenotyping analysis.

4.3. Advantages of Multisource Imaging Systems

This study has shown the capability of using both the fusion of PMD and RGB camera images and the proposed algorithm to measure the plant height and greenness index of soybean canopies rapidly and accurately. From plant science and breeding perspectives, the plant height trait could be used for selecting soybean genotypes. For example, plant height is used as an indicator of early maturing plant cultivars [41,42] to avoid yield losses that result from diseases and insect–pest complexes. In addition, measuring the greenness index over time at every growth stage could be used to indicate the effectiveness of plant energy consumption, which can be potentially used for the diagnosis of nutrient deficiencies. From a technical perspective, plant height and greenness can also be accurately measured with other imaging sensors such as light detection and ranging (LiDAR) [43] and more accurate RGB cameras in hue, saturation, and intensity (HSI) or hue, saturation, and value (HSV) color spaces [44], respectively. The use of PMDs to reconstruct the 3D geometry of plants has been explored, and the results of the present study demonstrated that a PMD camera operating in tandem with an RGB camera can provide accurate measurements of plant height and greenness, which creates new opportunities for field-based soybean plant phenotyping.

4.4. Future Work

In the present framework of phenotyping, plants are usually monitored one after another with noninvasive measurement devices. As an initial proof of feasibility, we have studied soybean plants cultivated in pots. However, a promising strategy could be developed to increase the throughput at larger observation scales by capturing multiple soybean plants in a single image. Additionally, under field conditions, the spatial arrangement of soybean plants is less regular than that in well-controlled environments. Thus, it is challenging to quantify the effects of crop spatial arrangement.

In addition, both views (side view and top view) are important for phenotyping analysis. Plant height, the leaf area index, and the greenness index can be measured accurately from the top view, but this view is not suitable for measuring other phenotypic traits such as leaf angle, branch angle, and stem diameter, which can be calculated effectively from the side view. Thus, future work will focus on improving the algorithms for 3D reconstruction from not only the side view but also the top view to acquire additional elaborate phenotypic traits. Notably, environmental factors—especially weather conditions—should be considered when data are collected in natural environments. Sunny days with low wind are considered the best conditions for phenotyping analysis.

5. Conclusions

In this paper, we developed a streamlined method to measure the specific phenotypic traits of soybean plants based on a 3D reconstruction containing color information. The main achievements are summarized as follows:

(1)

An active imaging system consisting of a PMD camera and an RGB camera was used to collect multi-images of soybean plants. First, the DBSCAN algorithm was used to extract soybean plant information from the complex raw dataset. Next, the multisource images were fused together for the purpose of constructing 3D images that contain color information. Last, 3D points from the front and back sides were registered using the ICP algorithm. The proposed methodology can be used to reconstruct a 3D soybean plant for a phenotyping analysis that includes measurements of plant height and greenness.

(2)

By combining this multisource imaging system and the proposed algorithms, we can accurately measure soybean plant height. Correlation analysis between the estimated and manual measurements yielded R² values of 0.9890 and 0.9936 for the side view and top view, respectively, and their average errors were 0.6713 cm and 0.2600 cm, respectively. From a plant breeding perspective, this finding could be especially useful for rapidly predetecting a subset of soybean genotypes that are of suitable height for expected yields and machine harvesting.

(3)

Compared with the side view-based greenness, the top view-based greenness was much more accurate. The greenness index estimated from the top view-based data was highly correlated with the manually assessed greenness index: the R² value was 0.8864, and the average error was 0.0117. However, the R² value decreased to 0.6059 (average error of 0.0386) for the side view-based results. This result was primarily due to the impact of the natural environment, such as wind and sunlight, which led to some fusion and registration deviations between the 3D points and their corresponding RGB images. The algorithm itself needs to be improved.