# Automatic Segmentation for Plant Leaves via Multiview Stereo Reconstruction

Authors: Jingwei Guo; Lihong Xu
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

This paper presented a new method for automatic plant point cloud acquisition and leaves segmentation. Quasi-dense point cloud of the plant is obtained from multiview stereo reconstruction based on surface expansion. In order to overcome the negative effects from complex natural light changes and to obtain a more accurate plant point cloud, the Adaptive Normalized Cross-Correlation algorithm is used in calculating the matching cost between two images, which is robust to radiometric factors and can reduce the fattening effect around boundaries. In the stage of segmentation for each single leaf, an improved region growing method based on fully connected conditional random field (CRF) is proposed to separate the connected leaves with similar color. The method has three steps: boundary erosion, initial segmentation, and segmentation refinement. First, the edge of each leaf point cloud is eroded to remove the connectivity between leaves. Then leaves will be initially segmented by region growing algorithm based on local surface normal and curvature. Finally an efficient CRF inference method based on mean field approximation is applied to remove small isolated regions. Experimental results show that our multiview stereo reconstruction method is robust to illumination changes and can obtain accurate color point clouds. And the improved region growing method based on CRF can effectively separate the connected leaves in obtained point cloud.


---

## Body

## 1. Introduction

In greenhouse crop cultivation, the three-dimensional structures of plants convey the real-time continuous information of the whole process of crop growth in a most direct way. Also, it is the most intuitive feature for plant phenotype analysis [1]. Therefore acquiring the 3D structure of plants by an automated, high-throughput, accurate, and rapid method not only helps to monitor crop growth, but also provides visual basis for optimizing the intelligent greenhouse environment control system; and finally, it has a significant meaning for plant phenotypic and genome-environment-phenotype synergistic analysis [2–6]. Owing to noninvasive and high-throughput characteristics, automated 3D-digitization technology has greatly promoted the accurate modeling of 3D structures of plants.To date, there are several approaches for 3D structure reconstruction. For example, Microsoft Kinect provides depth images by utilizing coded infrared light patterns. However, Kinect has limited range and can only be used indoors. Laser Scanning can obtain highly accurate and complete 3D structure but the devices are quite expensive for agricultural usage. Therefore the passive stereo vision technology, which puts low requirements on equipment and environment [7–10], has received wide attention. Biskup et al. [8] used a binocular stereo vision system to measure the structural parameters of plant canopies.Binocular stereo vision system needs to manually calibrate the camera and its imaging angle is relatively fixed. Hence the method is sometimes inefficient and it is difficult for it to reveal the complex 3D structure of plants. Compared to binocular stereo, the multiview stereo reconstruction based on an image sequence is used more widely because it can perform automatic calibration of the camera and generate more intact 3D structure information. In the literature, multiview stereo algorithms are roughly categorized into three classes. The first class operates by first computing a cost function on a 3D volume and then extracting a surface from this volume [11–14]. The second class is an image-space methodology that operates by fusing multiple depth maps [15–18]. The third class is surface expansion method which builds a quasi-dense point cloud from a set of matches [19–22]. The accuracy of volume-based method depends on the density of the grid, and therefore the computational and memory costs of volumetric grids are very high. The following two methods are both able to produce good results. In this paper, we propose a multiview stereo reconstruction method based on surface expansion which iteratively expands a sparse set of initial matches into a quasi-dense point cloud that represents surfaces of the scene. In order to overcome the negative effects from complex natural light changes and obtain more accurate plant point clouds, the Adaptive Normalized Cross-Correlation algorithm (ANCC) [23] is used in calculating the matching cost between points of interest in two images, which is robust to radiometric factors such as illumination direction, illuminant color, and imaging device changes. The ANCC method also can reduce the fattening effect that object boundaries are not reconstructed correctly; and this effect clearly deteriorates results of the Zero-mean Normalized Cross-Correlation (ZNCC)—the method used in comparison.We then segment the plant leaves after getting the plant point cloud from the above multiview stereo reconstruction method. Generally, the segmentation of 3D point cloud can be roughly categorized into five classes [24]: edge-based methods, region-based methods, attributes-based methods, model-based growing methods, and graph-based methods. References [25–27] segmented a plant into two sets, one for the stems and the other for the leaves according to surface feature histograms, followed by a classification or a clustering method. Yin et al. [28] proposed an approach to accurately segment and reconstruct plant. However their approach requires manually cutting the plant leaves and scanning each leaf individually. Hétroy-Wheeler et al. [29] used spectral clustering to realize the segmentation of leaf point cloud, but it needed to manually specify the number of leaves before segmentation.Based on the characteristics of the subject plants, an improved region growing method based on conditional random field (CRF) is proposed to segment the leaf point cloud. In a canopy point cloud, there may exist some connections between leaves, which easily disturb the segmentation. Therefore we perform some erosion operations around the leaf edge before segmentation. In the initial stage of segmentation, considering that different leaves usually have different surface normals and the points on the same leaf surface usually have similar normals, we use region growing algorithm based on local surface normal and curvature to segment leaves [30]. Usually the segmentation result contains some small isolated regions; therefore a fully connected CRF model defined with the complete set of points is built to further classify these regions to where they ought to belong. A highly efficient approximate inference algorithm based on mean field approximation is applied to the CRF model in which the pairwise edge potentials are defined by Gaussian kernel [31].


## 2. Multiview Stereo Reconstruction

The proposed multiview stereo approach consists of two steps: structure from motion (SFM: from images to sparse points) and dense matching (from sparse points to quasi-dense points). The SFM step is to obtain 3D sparse feature points location of the scene and camera poses (location and orientation). In the step of dense matching, we use the ANCC algorithm as a robust and accurate correspondence measure to match points of interest, which is robust to radiometric factors and can reduce the fattening effect.


### 2.1. Structure from Motion

An overview of the SFM process is presented in Figure1. The first step is to find feature points in each image, where the SIFT key point detector [32] is used because of its invariance to image transformations. Other feature detectors could also be used instead. For instance, several detectors were evaluated in [33]. In addition to the key point locations themselves, SIFT provides a local descriptor for each key point. For each pair of images, key point descriptors are matched by using the approximate nearest neighbors package [34]. Then we robustly estimate a fundamental matrix for the pair using RANSAC. During each iteration of RANSAC, a candidate fundamental matrix is computed using the eight-point algorithm [35], followed by nonlinear refinement based on Levenberg–Marquardt algorithm. Afterwards, we convert the fundamental matrix to the essential matrix and then obtain camera poses [35]. At last, the 3D positions of features points can be acquired through the triangulation calculation. (The scale factor is not considered during triangulation calculation because in the following stage leaf segmentation can be performed without the knowledge of a scale factor.) [35].Figure 1
A simplified overview of the SFM process.After generation of the point cloud of features points between each pair of images, we should merge these point clouds because they are not in the same coordinate system, as shown in Figure2. Finally, we get 3D positions of all the feature points in a single coordinate system with respect to all point clouds.Figure 2
Merging point clouds in one coordinate system.


### 2.2. Dense Matching

In this stage, we expand a sparse set of initial matches into a quasi-dense point cloud representing the surfaces of the scene based on the best-first expansion strategy by Lhuillier and Quan [36], which is robust to initial seed match outliers and is efficient in terms of spatial-temporal complexity. In order to overcome the negative effects from complex natural light environments and to obtain a more accurate plant point cloud, the ANCC algorithm is utilized in the process as the matching cost to match points of interest instead of ZNCC. It is noteworthy that the ANCC method is robust to radiometric factors such as illumination direction, illuminant color, and imaging device changes. Moreover, it can reduce the fattening effect from which the Zero-mean Normalized Cross-Correlation (ZNCC) suffers. The best-first expansion strategy is presented in the table of Algorithm 1.Algorithm 1:The multiview best-first match expansion algorithm.
Input: Seed Output: Map (1)Map←⏀ (2)while Seed!=⏀do (3) pull the ANCC-best matchx,x′ from Seed (4)Local←⏀ (5)for eachu,u′inNx,x′do (6)ifsu>t and su′>t and ANCCu,u′>z (7)then store u,u′ in Local (8)end-if (9)end-for (10)while Local!=⏀ (11)  pull the ZNCC-best matchu,u′ from Local (12)ifu,∗ and ∗,u′ are not in Map (13)then store u,u′ in Map and Seed (14)end-if (15)end-while (16)end-whileIn Algorithm1, su is a difference-based confidence measure to avoid propagating too uniform areas. The definition is (1)sx=max⁡Ix+Δ-Ix,Δ∈1,0,-1,0,0,1,0,-1.Ix means the pixel grayscale value at x. Nx,x′ means the neighbors of the seed. Let N(x)={u,u-x∈-N,N2}, Nx′={u′,u′-x′∈-N,N2} denote all 2N+1×2N+1 neighboring pixels of pixels x and x′, respectively. The possible matches limited by the discrete 2D disparity gradient are given as(2)Nx,x′=u,u′,u∈Nx,u′∈Nx′,u′-u-x′-x≤ε.Figure 3 illustrates a case for N=2, ε=1.Figure 3
Possible matchesu,u′ and v,v′ around a seed match x,x′ come from its 5×5 neighbor N(x) and Nx′ as the smallest size for discrete 2D disparity gradient limit. The match candidates for u (resp., v′) are within the 3×3 (black framed) centered at u′ (resp., v).In Algorithm1, ANCC(u,u′) is similarity measure to match points of interest in two images. Generally, the color image formation model at pixel p can be represented as follows:(3)RpGpBp⟶R~pG~pB~p=ρpaRpγρpbGpγρpcBpγ,where the brightness factor ρ(p) depends on the angle between the direction of light and the direction of surface normal at that point. The global scale factors a,b,c depend on illumination color and γ depends on camera. The well-known similarity measure ZNCC that uses raw stereo images does not generate satisfactory result for two reasons. One is that various radiometric changes caused by ρ,a,b,c, and γ are not taken into consideration. The other is that the object boundaries as the supporting windows from the left and right images do not appear exactly because of the view changes.To keep the color values away from the impact of radiometric factors, we transform the nonlinear relationship between the corresponding pixels to the affine-transformed relationship that can be solved effectively by the NCC framework. First, we take logarithm to transform this nonlinear relationship into a linear one. Then, the color value can be represented by (use the red channel to describe)(4)RL′p=log⁡ρLp+log⁡aL+γLRLp,RR′p′=log⁡ρRp′+log⁡aR+γRRRp′.To delete the item log⁡ρL(p), we subtract the average of the transformed color values in R, G, B channels. Then each color value becomes(5)RL′′p=log⁡aLaLbLcL3+γLlog⁡RLpRLpGLpBLp3≜αL+γLKLp,RR′′p′≜αR+γRKRp′.We can see that KL(p) and KR(p′) are not dependent on ρ,a,b,c, and γ.In order to reduce the fattening effect caused by outliers in the matching windows, each pixelt in m×m window Wp around the center pixel p is allocated with different weights. The weight is computed by using the bilateral filter as follows:(6)ωt=exp⁡-p-t22σd2-Ip-It22σs2.The weighted sum using the bilateral filter for the center pixel p is defined by the following equation, in which Zp is the normalizing constant:(7)Sp=∑t∈WpωtItZp.Then we subtract S(p) for each channel, thereby resulting in the removal of α.(8)RL′′′p=RL′′p-SLp=γLKLp-∑t∈WpωLtKLtZp.The similarity betweenp and p′ is defined by the following equation that is close to NCC’s framework:(9)ANCClog⁡Chrom_Rp,p′=∑i=1MωLtiωRtiRL′′′tiRR′′′ti∑i=1MωLtiRL′′′ti2×∑i=1MωRtiRR′′′ti2.Because the intensity information of each channel is normalized in log-chromaticity, the discriminability is lower than the original color. We therefore combine log-chromaticity and original color for keeping the discriminability. The final similarity equation is defined as follows:(10)ANCCp,p′=θ∑ξANCCξp,p′3+1-θ∑kANCCkp,p′3,where ξ∈{log⁡Chrom_R,log⁡Chrom_G,log⁡Chrom_B},k∈{R,G,B}. θ is a balancing factor between the log-chromaticity color and the original RGB color.To obtain the disparity map of the two images in different lighting conditions, we define stereo matching as a minimization problem of the following energy in the MAP-MRF framework:(11)Ef=Edataf+Esmoothf,Edataf=∑pDpfp,Esmoothf=∑p∑q∈NpVpqfp,fq,where N(p) is the neighborhood of pixel p and f is the disparity value we are trying to find for each pixel p. Dp(fp) is the data cost function which is defined by Dp(fp)=1-ANCC(p,p+fp). Vpq(fp,fq) is the smoothness cost function based on truncated quadratic model which is defined by Vpq(fp,fq)=λ∗min⁡fp-fq2,Vmax. The total energy E(f) was optimized using the Graph Cuts (GC) algorithm [37] in order to find the disparity map f. Under different lighting conditions, it can be seen from Figure 4 that the algorithm using ANCC as the matching cost function is very robust to illumination changes compared with results of two classic algorithms SGBM and Graph Cuts by OpenCV library.Figure 4
Results of test stereo algorithms on Aloe image pair in different lighting conditions. (a) The left image in lighting condition 1. (b) The right image in lighting condition 2. (c) The ground truth disparity map. (d) The disparity map of SGBM by OpenCV. (e) The disparity map of GC by OpenCV. (f) The disparity map of ANCC+GC.
(a) (b) (c) (d) (e) (f)The point clouds from multiview stereo reconstruction in different matching method are illustrated in Figures5 and 6. From the results we can see that many points are not reconstructed correctly (marked in red circles) based on the ZNCC matching method because of the fattening effect. And more accurate points are acquired with our reconstruction method based on ANCC in Figure 6.Figure 5
Multiview stereo reconstruction based on ZNCC: (a) bust; (b)Scindapsus; (c)Pachira macrocarpa.
(a) (b) (c)Figure 6
Multiview stereo reconstruction based on ANCC: (a) bust; (b)Scindapsus. (c)Pachira macrocarpa.
(a) (b) (c)


## 2.1. Structure from Motion

An overview of the SFM process is presented in Figure1. The first step is to find feature points in each image, where the SIFT key point detector [32] is used because of its invariance to image transformations. Other feature detectors could also be used instead. For instance, several detectors were evaluated in [33]. In addition to the key point locations themselves, SIFT provides a local descriptor for each key point. For each pair of images, key point descriptors are matched by using the approximate nearest neighbors package [34]. Then we robustly estimate a fundamental matrix for the pair using RANSAC. During each iteration of RANSAC, a candidate fundamental matrix is computed using the eight-point algorithm [35], followed by nonlinear refinement based on Levenberg–Marquardt algorithm. Afterwards, we convert the fundamental matrix to the essential matrix and then obtain camera poses [35]. At last, the 3D positions of features points can be acquired through the triangulation calculation. (The scale factor is not considered during triangulation calculation because in the following stage leaf segmentation can be performed without the knowledge of a scale factor.) [35].Figure 1
A simplified overview of the SFM process.After generation of the point cloud of features points between each pair of images, we should merge these point clouds because they are not in the same coordinate system, as shown in Figure2. Finally, we get 3D positions of all the feature points in a single coordinate system with respect to all point clouds.Figure 2
Merging point clouds in one coordinate system.


## 2.2. Dense Matching

In this stage, we expand a sparse set of initial matches into a quasi-dense point cloud representing the surfaces of the scene based on the best-first expansion strategy by Lhuillier and Quan [36], which is robust to initial seed match outliers and is efficient in terms of spatial-temporal complexity. In order to overcome the negative effects from complex natural light environments and to obtain a more accurate plant point cloud, the ANCC algorithm is utilized in the process as the matching cost to match points of interest instead of ZNCC. It is noteworthy that the ANCC method is robust to radiometric factors such as illumination direction, illuminant color, and imaging device changes. Moreover, it can reduce the fattening effect from which the Zero-mean Normalized Cross-Correlation (ZNCC) suffers. The best-first expansion strategy is presented in the table of Algorithm 1.Algorithm 1:The multiview best-first match expansion algorithm.
Input: Seed Output: Map (1)Map←⏀ (2)while Seed!=⏀do (3) pull the ANCC-best matchx,x′ from Seed (4)Local←⏀ (5)for eachu,u′inNx,x′do (6)ifsu>t and su′>t and ANCCu,u′>z (7)then store u,u′ in Local (8)end-if (9)end-for (10)while Local!=⏀ (11)  pull the ZNCC-best matchu,u′ from Local (12)ifu,∗ and ∗,u′ are not in Map (13)then store u,u′ in Map and Seed (14)end-if (15)end-while (16)end-whileIn Algorithm1, su is a difference-based confidence measure to avoid propagating too uniform areas. The definition is (1)sx=max⁡Ix+Δ-Ix,Δ∈1,0,-1,0,0,1,0,-1.Ix means the pixel grayscale value at x. Nx,x′ means the neighbors of the seed. Let N(x)={u,u-x∈-N,N2}, Nx′={u′,u′-x′∈-N,N2} denote all 2N+1×2N+1 neighboring pixels of pixels x and x′, respectively. The possible matches limited by the discrete 2D disparity gradient are given as(2)Nx,x′=u,u′,u∈Nx,u′∈Nx′,u′-u-x′-x≤ε.Figure 3 illustrates a case for N=2, ε=1.Figure 3
Possible matchesu,u′ and v,v′ around a seed match x,x′ come from its 5×5 neighbor N(x) and Nx′ as the smallest size for discrete 2D disparity gradient limit. The match candidates for u (resp., v′) are within the 3×3 (black framed) centered at u′ (resp., v).In Algorithm1, ANCC(u,u′) is similarity measure to match points of interest in two images. Generally, the color image formation model at pixel p can be represented as follows:(3)RpGpBp⟶R~pG~pB~p=ρpaRpγρpbGpγρpcBpγ,where the brightness factor ρ(p) depends on the angle between the direction of light and the direction of surface normal at that point. The global scale factors a,b,c depend on illumination color and γ depends on camera. The well-known similarity measure ZNCC that uses raw stereo images does not generate satisfactory result for two reasons. One is that various radiometric changes caused by ρ,a,b,c, and γ are not taken into consideration. The other is that the object boundaries as the supporting windows from the left and right images do not appear exactly because of the view changes.To keep the color values away from the impact of radiometric factors, we transform the nonlinear relationship between the corresponding pixels to the affine-transformed relationship that can be solved effectively by the NCC framework. First, we take logarithm to transform this nonlinear relationship into a linear one. Then, the color value can be represented by (use the red channel to describe)(4)RL′p=log⁡ρLp+log⁡aL+γLRLp,RR′p′=log⁡ρRp′+log⁡aR+γRRRp′.To delete the item log⁡ρL(p), we subtract the average of the transformed color values in R, G, B channels. Then each color value becomes(5)RL′′p=log⁡aLaLbLcL3+γLlog⁡RLpRLpGLpBLp3≜αL+γLKLp,RR′′p′≜αR+γRKRp′.We can see that KL(p) and KR(p′) are not dependent on ρ,a,b,c, and γ.In order to reduce the fattening effect caused by outliers in the matching windows, each pixelt in m×m window Wp around the center pixel p is allocated with different weights. The weight is computed by using the bilateral filter as follows:(6)ωt=exp⁡-p-t22σd2-Ip-It22σs2.The weighted sum using the bilateral filter for the center pixel p is defined by the following equation, in which Zp is the normalizing constant:(7)Sp=∑t∈WpωtItZp.Then we subtract S(p) for each channel, thereby resulting in the removal of α.(8)RL′′′p=RL′′p-SLp=γLKLp-∑t∈WpωLtKLtZp.The similarity betweenp and p′ is defined by the following equation that is close to NCC’s framework:(9)ANCClog⁡Chrom_Rp,p′=∑i=1MωLtiωRtiRL′′′tiRR′′′ti∑i=1MωLtiRL′′′ti2×∑i=1MωRtiRR′′′ti2.Because the intensity information of each channel is normalized in log-chromaticity, the discriminability is lower than the original color. We therefore combine log-chromaticity and original color for keeping the discriminability. The final similarity equation is defined as follows:(10)ANCCp,p′=θ∑ξANCCξp,p′3+1-θ∑kANCCkp,p′3,where ξ∈{log⁡Chrom_R,log⁡Chrom_G,log⁡Chrom_B},k∈{R,G,B}. θ is a balancing factor between the log-chromaticity color and the original RGB color.To obtain the disparity map of the two images in different lighting conditions, we define stereo matching as a minimization problem of the following energy in the MAP-MRF framework:(11)Ef=Edataf+Esmoothf,Edataf=∑pDpfp,Esmoothf=∑p∑q∈NpVpqfp,fq,where N(p) is the neighborhood of pixel p and f is the disparity value we are trying to find for each pixel p. Dp(fp) is the data cost function which is defined by Dp(fp)=1-ANCC(p,p+fp). Vpq(fp,fq) is the smoothness cost function based on truncated quadratic model which is defined by Vpq(fp,fq)=λ∗min⁡fp-fq2,Vmax. The total energy E(f) was optimized using the Graph Cuts (GC) algorithm [37] in order to find the disparity map f. Under different lighting conditions, it can be seen from Figure 4 that the algorithm using ANCC as the matching cost function is very robust to illumination changes compared with results of two classic algorithms SGBM and Graph Cuts by OpenCV library.Figure 4
Results of test stereo algorithms on Aloe image pair in different lighting conditions. (a) The left image in lighting condition 1. (b) The right image in lighting condition 2. (c) The ground truth disparity map. (d) The disparity map of SGBM by OpenCV. (e) The disparity map of GC by OpenCV. (f) The disparity map of ANCC+GC.
(a) (b) (c) (d) (e) (f)The point clouds from multiview stereo reconstruction in different matching method are illustrated in Figures5 and 6. From the results we can see that many points are not reconstructed correctly (marked in red circles) based on the ZNCC matching method because of the fattening effect. And more accurate points are acquired with our reconstruction method based on ANCC in Figure 6.Figure 5
Multiview stereo reconstruction based on ZNCC: (a) bust; (b)Scindapsus; (c)Pachira macrocarpa.
(a) (b) (c)Figure 6
Multiview stereo reconstruction based on ANCC: (a) bust; (b)Scindapsus. (c)Pachira macrocarpa.
(a) (b) (c)


## 3. Leaves Segmentation from Point Clouds

To segment each individual leaf from point clouds, an improved region growing method based on fully connected CRF is applied in this section. The method can be divided into three steps: boundary erosion, initial segmentation, and refinement. At the beginning, the edge of each leaf point cloud is eroded to remove the connectivity between leaves. Then leaves will be initially segmented by region growing algorithm based on local surface normal and curvature. Finally an efficient CRF inference method based on mean field approximation is used to remove small isolated regions. The details are further explained as follows.


### 3.1. Boundary Erosion and Initial Segmentation

Because of the similarity of the colors of leaves, the color information is not reliable for segmentation. Considering the fact that the points on the same leaf surface usually have similar normals, and different leaves usually have different mean normals, the region growing algorithm is used. This method can find smoothly connected areas in point cloud based on point normals. First, the normal for each point is estimated by fitting a plane to some neighboring points. The neighborhood can be specified in two different methods, namely,K nearest neighbors (KNN) and fixed distance neighbors (FDN). The next is the region growing algorithm based on normal and curvature, which is presented in Algorithm 2.Algorithm 2:The region growing algorithm for point cloud segmentation.
Input: Point cloud = P, point normals N, residuals r, neighbor finding function Ω(⋅), residual threshold γth, angle threshold θth Output: Leaf point cloud R (1)Initialize Region List R←⏀, Available points list A←{1⋯Pcount} (2)whileA is not emptydo (3) Current regionRc←⏀, Current seeds Sc←⏀ (4) Point with minimum residual inA→Pmin (5)Pmin→insert{Sc}&{Rc} (6)Pmin→removeA (7)fori=0 to size ({Sc}) do (8)  Find nearest neighbors of current seed pointBc←Ω(Sc{i}) (9)forj=0 to size ({Bc}) do (10)    Current neighbor pointPj←Bcj (11)ifA contains Pj and cos-1⁡(NSci,NPj)<θth (12)thenPj→insertRc Pj→removeA (13)ifrPj<γth (14)thenPj→insert{Sc} (15)end-if (16)end-if (17)end-for (18)end-for (19) Add current region to global segment listRc→insertR (20)end- while (21) Sort{R} according to the size of the region. (22)ReturnRThis algorithm is based on two constraints—local connectivity and surface smoothness. For local connectivity, the points in a segment should be locally connected. This can be realized by using only the neighboring points (through KNN or FDN) during region growing. For surface smoothness, the points in a segment should create a local smooth surface whose normals do not vary much. This constraint can be realized with a thresholdθth on the angles between the current seed point and the points added to the region. Additionally, a threshold on curvature values γth can ensure that smooth areas are broken on the edges. The result of region growing segmentation without boundary erosion for a canopy point cloud is shown in Figure 7.Figure 7
Region growing segmentation without boundary erosion for the point cloud fromPachira macrocarpa. θth=3.25, γth=1.0. The red points are small isolated regions because of too little points, which we will remove in the future refinement stage.As it can be seen from Figure7 some leaves are segmented into one part, which indicates relatively poor performance because some leaves are connected (marked by red circles) and these leaves have similar normals. To solve the problem, some erosion operations around the leaf edge are performed before segmentation. First, the boundary of the point cloud is extracted by detecting the boundary of the corresponding depth map. The boundary of the point cloud is illustrated in Figure 8. Then we constructk-d tree based on all points including the boundary ones. For every boundary point, we find the k nearest points nearby and then delete them. After erosion, we carry out the region growing algorithm for the point cloud and the new result is shown in Figure 9. It can be seen that the result is much better than previous result without boundary erosion and also leaves that were originally connected are now separated.Figure 8
The boundary of the point cloud fromPachira macrocarpa leaves.Figure 9
Region growing segmentation after boundary corrosion.θth=3.25, γth=1.0.


### 3.2. Segmentation Refinement

After initial segmentation, there still exist some small isolated regions (small red points region). Therefore we need to classify these regions and determine the actual regions they belong to. A common solution is to transform this problem as maximum a posteriori (MAP) inference in a conditional random field (CRF) defined over pixels or image patches [38–41]. The CRF model contains pairwise potentials, which can simultaneously remove the noise and resolve ambiguous classification results.In this subsection, a fully connected CRF is built, establishing pairwise potentials on all pairs of points from the canopy point cloud. Considering a random fieldX defined over a set of variables X1,…,XN, the domain of each variable is a set of labels L=l1,l2,…,lk. Also a random field I is defined over variables I1,…,IN. In our setting, I ranges over possible input point cloud of size N and X ranges over possible point-level labelings. Ij is the position vector of point j and Xj is the label assigned to point j. Hence, the corresponding Gibbs energy of the fully connected pairwise CRF model is formulated as(12)Ex=∑iψuxi+∑i<jψpxi,xj,where i and j both range from 1 to N. A conditional random field I,X is characterized by a Gibbs distribution P(x∣I)=exp-E(x)/Z(I). The maximum a posteriori (MAP) labeling of the random field is x∗=argmaxx∈LNP(x∣I). The unary potential is computed independently for each point with a classifier that produces a distribution over the label assignment xi given point cloud features. The following is the function of the pairwise potentials in our model:(13)ψpxi,xj=∑mμmxi,xjkmfi,fj,where k(m) is a kernel function and μ(m) is a label compatibility function. A simple label compatibility function μ is given by the Potts model μ(xi,xj)=[xi≠xj]. We limit the choice of the kernel function to Gaussian kernel:(14)kmfi,fj=ωmexp⁡-12fi-fjTΛmfi-fj,where ω(m) is a linear combination weight and Λ(m) is a symmetric, positive-definite precision matrix, defining the shape of the kernel. The vectors fi and fj are feature vectors for points i and j in an arbitrary feature space, respectively. For our points labeling problem, the smoothness kernel k(fi,fj)=exp⁡[-pi-pj2/2θγ2] is used, defined in terms of point positions pi and pj to remove isolated regions.Usually the high complexity of inference in fully connected models will limit its application. Hence, a highly efficient inference algorithm based on a mean field approximation to the CRF distribution [31] is adopted. This approximation yields an iterative message passing algorithm for approximate inference. As message passing in the presented model can be performed using Gaussian filtering in feature space, highly efficient approximations for high-dimensional filtering are utilized to reduce the complexity of message passing from quadratic to linear.


## 3.1. Boundary Erosion and Initial Segmentation

Because of the similarity of the colors of leaves, the color information is not reliable for segmentation. Considering the fact that the points on the same leaf surface usually have similar normals, and different leaves usually have different mean normals, the region growing algorithm is used. This method can find smoothly connected areas in point cloud based on point normals. First, the normal for each point is estimated by fitting a plane to some neighboring points. The neighborhood can be specified in two different methods, namely,K nearest neighbors (KNN) and fixed distance neighbors (FDN). The next is the region growing algorithm based on normal and curvature, which is presented in Algorithm 2.Algorithm 2:The region growing algorithm for point cloud segmentation.
Input: Point cloud = P, point normals N, residuals r, neighbor finding function Ω(⋅), residual threshold γth, angle threshold θth Output: Leaf point cloud R (1)Initialize Region List R←⏀, Available points list A←{1⋯Pcount} (2)whileA is not emptydo (3) Current regionRc←⏀, Current seeds Sc←⏀ (4) Point with minimum residual inA→Pmin (5)Pmin→insert{Sc}&{Rc} (6)Pmin→removeA (7)fori=0 to size ({Sc}) do (8)  Find nearest neighbors of current seed pointBc←Ω(Sc{i}) (9)forj=0 to size ({Bc}) do (10)    Current neighbor pointPj←Bcj (11)ifA contains Pj and cos-1⁡(NSci,NPj)<θth (12)thenPj→insertRc Pj→removeA (13)ifrPj<γth (14)thenPj→insert{Sc} (15)end-if (16)end-if (17)end-for (18)end-for (19) Add current region to global segment listRc→insertR (20)end- while (21) Sort{R} according to the size of the region. (22)ReturnRThis algorithm is based on two constraints—local connectivity and surface smoothness. For local connectivity, the points in a segment should be locally connected. This can be realized by using only the neighboring points (through KNN or FDN) during region growing. For surface smoothness, the points in a segment should create a local smooth surface whose normals do not vary much. This constraint can be realized with a thresholdθth on the angles between the current seed point and the points added to the region. Additionally, a threshold on curvature values γth can ensure that smooth areas are broken on the edges. The result of region growing segmentation without boundary erosion for a canopy point cloud is shown in Figure 7.Figure 7
Region growing segmentation without boundary erosion for the point cloud fromPachira macrocarpa. θth=3.25, γth=1.0. The red points are small isolated regions because of too little points, which we will remove in the future refinement stage.As it can be seen from Figure7 some leaves are segmented into one part, which indicates relatively poor performance because some leaves are connected (marked by red circles) and these leaves have similar normals. To solve the problem, some erosion operations around the leaf edge are performed before segmentation. First, the boundary of the point cloud is extracted by detecting the boundary of the corresponding depth map. The boundary of the point cloud is illustrated in Figure 8. Then we constructk-d tree based on all points including the boundary ones. For every boundary point, we find the k nearest points nearby and then delete them. After erosion, we carry out the region growing algorithm for the point cloud and the new result is shown in Figure 9. It can be seen that the result is much better than previous result without boundary erosion and also leaves that were originally connected are now separated.Figure 8
The boundary of the point cloud fromPachira macrocarpa leaves.Figure 9
Region growing segmentation after boundary corrosion.θth=3.25, γth=1.0.


## 3.2. Segmentation Refinement

After initial segmentation, there still exist some small isolated regions (small red points region). Therefore we need to classify these regions and determine the actual regions they belong to. A common solution is to transform this problem as maximum a posteriori (MAP) inference in a conditional random field (CRF) defined over pixels or image patches [38–41]. The CRF model contains pairwise potentials, which can simultaneously remove the noise and resolve ambiguous classification results.In this subsection, a fully connected CRF is built, establishing pairwise potentials on all pairs of points from the canopy point cloud. Considering a random fieldX defined over a set of variables X1,…,XN, the domain of each variable is a set of labels L=l1,l2,…,lk. Also a random field I is defined over variables I1,…,IN. In our setting, I ranges over possible input point cloud of size N and X ranges over possible point-level labelings. Ij is the position vector of point j and Xj is the label assigned to point j. Hence, the corresponding Gibbs energy of the fully connected pairwise CRF model is formulated as(12)Ex=∑iψuxi+∑i<jψpxi,xj,where i and j both range from 1 to N. A conditional random field I,X is characterized by a Gibbs distribution P(x∣I)=exp-E(x)/Z(I). The maximum a posteriori (MAP) labeling of the random field is x∗=argmaxx∈LNP(x∣I). The unary potential is computed independently for each point with a classifier that produces a distribution over the label assignment xi given point cloud features. The following is the function of the pairwise potentials in our model:(13)ψpxi,xj=∑mμmxi,xjkmfi,fj,where k(m) is a kernel function and μ(m) is a label compatibility function. A simple label compatibility function μ is given by the Potts model μ(xi,xj)=[xi≠xj]. We limit the choice of the kernel function to Gaussian kernel:(14)kmfi,fj=ωmexp⁡-12fi-fjTΛmfi-fj,where ω(m) is a linear combination weight and Λ(m) is a symmetric, positive-definite precision matrix, defining the shape of the kernel. The vectors fi and fj are feature vectors for points i and j in an arbitrary feature space, respectively. For our points labeling problem, the smoothness kernel k(fi,fj)=exp⁡[-pi-pj2/2θγ2] is used, defined in terms of point positions pi and pj to remove isolated regions.Usually the high complexity of inference in fully connected models will limit its application. Hence, a highly efficient inference algorithm based on a mean field approximation to the CRF distribution [31] is adopted. This approximation yields an iterative message passing algorithm for approximate inference. As message passing in the presented model can be performed using Gaussian filtering in feature space, highly efficient approximations for high-dimensional filtering are utilized to reduce the complexity of message passing from quadratic to linear.


## 4. Results

The proposed method was applied to four sets of point clouds acquired fromPachira macrocarpa,Scindapsus, strawberry, and tomato plant.


### 4.1. Initial Segmentation

Results of the initial segmentation are shown in Figure10. Some erosion operations were performed around the leaf edges beforehand in case that two leaves are divided into one part.Figure 10
Initial segmentation results. (a)Pachira macrocarpa. θth=3.25, γth=1.0. (b)Scindapsus. θth=4.25, γth=1.0. (c) Strawberry. θth=3.0, γth=1.0. (d) Tomato. θth=4.2, γth=1.0.
(a) (b) (c) (d)When a largeθth of angle threshold is set, the algorithm tends to segment the point cloud into connected subsets. The smaller θth, the higher probability that a leaf is segmented into several clusters (see Figure 11).Figure 11
Initial segmentation results under different values ofθth. (a) θth=2.75. (b) θth=3.00. (c) θth=3.50. (d) θth=3.75.
(a) (b) (c) (d)


### 4.2. Final Segmentation Results and Comparison with the Original Region Growing Algorithm

The results after inference of fully connected CRF are shown in the second row of Figure12. It can be seen that all the small isolated regions are classified into its nearest leaves and all leaves are individually segmented.Figure 12
Segmentation results. The original region growing algorithm: (a)Pachira macrocarpa. θth=3.50, γth=1.0. (b)Scindapsus. θth=4.25, γth=1.0. (c) Strawberry. θth=3.7, γth=1.0. (d) Tomato. θth=4.2, γth=1.0. Our algorithm: (e)Pachira macrocarpa. θth=3.25, γth=1.0. (f)Scindapsus. θth=4.25, γth=1.0. (g) Strawberry. θth=3.0, γth=1.0. (h) Tomato. θth=4.2, γth=1.0.
(a) (b) (c) (d) (e) (f) (g) (h)For the original region growing algorithm, no matter how we adjust the values ofθth and γth, there always exists the case that two or more leaves are segmented into one subset, which is shown in the first row (highlighted by red circles) of Figure 12. Meanwhile segmentation results contain many noise sections that are not classified. The results of our method are shown in the second row of Figure 12, and the results are superior than the first row.


## 4.1. Initial Segmentation

Results of the initial segmentation are shown in Figure10. Some erosion operations were performed around the leaf edges beforehand in case that two leaves are divided into one part.Figure 10
Initial segmentation results. (a)Pachira macrocarpa. θth=3.25, γth=1.0. (b)Scindapsus. θth=4.25, γth=1.0. (c) Strawberry. θth=3.0, γth=1.0. (d) Tomato. θth=4.2, γth=1.0.
(a) (b) (c) (d)When a largeθth of angle threshold is set, the algorithm tends to segment the point cloud into connected subsets. The smaller θth, the higher probability that a leaf is segmented into several clusters (see Figure 11).Figure 11
Initial segmentation results under different values ofθth. (a) θth=2.75. (b) θth=3.00. (c) θth=3.50. (d) θth=3.75.
(a) (b) (c) (d)


## 4.2. Final Segmentation Results and Comparison with the Original Region Growing Algorithm

The results after inference of fully connected CRF are shown in the second row of Figure12. It can be seen that all the small isolated regions are classified into its nearest leaves and all leaves are individually segmented.Figure 12
Segmentation results. The original region growing algorithm: (a)Pachira macrocarpa. θth=3.50, γth=1.0. (b)Scindapsus. θth=4.25, γth=1.0. (c) Strawberry. θth=3.7, γth=1.0. (d) Tomato. θth=4.2, γth=1.0. Our algorithm: (e)Pachira macrocarpa. θth=3.25, γth=1.0. (f)Scindapsus. θth=4.25, γth=1.0. (g) Strawberry. θth=3.0, γth=1.0. (h) Tomato. θth=4.2, γth=1.0.
(a) (b) (c) (d) (e) (f) (g) (h)For the original region growing algorithm, no matter how we adjust the values ofθth and γth, there always exists the case that two or more leaves are segmented into one subset, which is shown in the first row (highlighted by red circles) of Figure 12. Meanwhile segmentation results contain many noise sections that are not classified. The results of our method are shown in the second row of Figure 12, and the results are superior than the first row.


## 5. Discussion and Conclusion

We proposed an automatic method to acquire and segment the point cloud of plant leaves. A multiview stereo reconstruction method based on surface expansion is first applied to obtain the 3D structure of plants. We use the Adaptive Normalized Cross-Correlation to match points of interest in two images, and it is robust to radiometric factors and can reduce the fattening effect. In the stage of leaf segmentation on point clouds, an improved region growing method based on fully connected conditional random field (CRF) is proposed to separate connected leaves with similar color. Our method requires no prior botanical knowledge; therefore it can be applied in a wide variety of cases. Qualitative results on several kinds of plants show the effectiveness of our method. To prevent two leaves from being segmented into only one subset, some erosion operations around the leaf edge are performed before carrying out the region growing segmentation based on normal and curvature. A fully connected CRF model is built in which the pairwise edge potentials are defined by Gaussian kernel to remove small isolated regions after segmentation.In the future, we will focus on developing advanced multiview stereo reconstruction algorithms to obtain more accurate point clouds for plant phenotype analysis. In addition, we will also study advanced point cloud segmentation algorithms to obtain other organs such as stems and fruits from plants, for realizing the final goal of accurate measurement on plant organs.



---

*Source: 9845815-2017-02-20.xml*