Unsupervised segmentation of greenhouse plant images based on modified Latent Dirichlet Allocation

Introduction

Plant phenotype analysis based on image processing has been a popular application field of agricultural computer vision in recent years. An automatic, high-throughput, accurate and rapid imaging technique and processing method for plant phenotypic analysis will not only monitor the growth of plants, but also lay a visual foundation for the optimization of an intelligent greenhouse environment control system. It is a benefit to gene breeding, environment regulation, cultivation management, and the estimation of yield and quality. In the literature, the convolutional neural network (CNN) has been applied to agriculture, however, a large amount of reliable training data is needed which becomes an urgent problem. For the analysis of greenhouse plants, it is an important process to get enough high-quality labelled images. In this regard, a well-segmented result of the plant and leaves can help in labelling the image quickly and accurately.

There are various imaging methods in the literature, among which a large amount of segmentation methods are developed from different aspects. Currently, there are several common methods to collect information of plants: (1) typical 2D color camera, (2) RGB-D camera, including stereo camera, ToF camera, structured light Kinect and laser scanner, (3) spectrometer, and (4) airborne laser radar. It is known that 2D imaging technologies are simple, low-cost and non-destructive, which can be suitable for plant phenotype analysis such as studying the growth of plant leaves, tracing their position and orientation, classifying new leaves from old ones, examining growth regulation and evaluating environmental stress if an accurate result of leaf segmentation is obtained (Scharr et al., 2016). The 2D imaging technologies lay a visual foundation for the analysis and warning of leaf disease as well.

For plant segmentation, the use of color in combination with depth images or multi-view images for supervised or unsupervised segmentation is a popular practice (Alenya, Dellen & Torras, 2011; Song et al., 2007; Teng, Kuo & Chen, 2011). Additionally, researchers have proposed some methods of segmenting plants from their background, in order to handle leaves with lesions (Zou et al., 2015; Valliammal & Geethalakshmi, 2012) and those without lesions (Zhang, Zhang & Guo-Hong, 2017; Wang, Wang & Cui, 2011). The background of the plants is artificial or relatively brief. One method for the detection of leaf disease spots is to directly extract the lesion part from the leaf of interest (Ma et al., 2017), while the other method is to separate a whole leaf from the complex background and extract the lesion part subsequently (Fang, Lu & Lisi, 2014). The latter can provide a better view to observe the site of lesions and analyze the severity of the disease. For single leaf segmentation, methods can be roughly categorized into four classes: shape model constraints (Manh et al., 2001; Mezzo et al., 2007; Cerutti et al., 2011; Sogaard, 2005; Persson & Astrand, 2008), boundary information detection (Valliammal & Geethalakshmi, 2011; Noble & Brown, 2008; Tang et al., 2009; Lee & Slaughter, 2004), depth information integration (Guo & Xu, 2017; Alenya et al., 2013; Xia et al., 2015; Sanz-Cortiella et al., 2011) and machine learning methods (Zheng, Zhang & Wang, 2009; Meunkaewjinda et al., 2008; Hernández-Rabadán, Guerrero & Ramos-Quintana, 2012; Pape & Klukas, 2015). Recently, the end-to-end recurrent neural network (RCNN) architecture (Ren & Zemel, 2016; Romera-Paredes & Torr, 2016), using an attention mechanism to model a human-like counting process, has provided a new way to handle occlusion by proceeding sequentially. Some of the methods above are mainly utilized in the classification of plants by extracting one complete leaf from a plant, while other methods segment all leaves in the image. Compared with the former, the latter methods have advantages in applications like genome-environment-phenotype synergistic analysis which cannot only rely on the information of one single leaf.

In recent years, more and more statistical methods have been applied in the image processing field. The Latent Dirichlet Allocation (LDA) model (Blei, Ng & Jordan, 2003) introduces a parameter θ which obeys Dirichlet distribution on the basis of Probabilistic Latent Semantic Analysis (pLSA) to establish the probability distribution of latent topic variable z. Li & Perona (2005) first applied the LDA model to image segmentation and many improved models came into view later. For example, a Spatial Latent Dirichlet Allocation (SLDA) topic model (Wang & Grimson, 2008) encodes the spatial structure of visual words as a random hidden variable in a better way with LDA’s generative procedure. A spatially coherent latent topic model (Spatial-LTM) (Cao & Li, 2011) provides a unified representation for a spatially coherent bag of word topic models, which can simultaneously segment and classify objects in the case of occlusion and multiple instances. Russell et al. (2006) partition a set of segmented objects into visual object classes using LDA and the visual object classes are further used to assess the accuracy of a segmentation. Reddy, Singhal & Krishna (2014) propose an algorithm that jointly infers the semantic class and motion labels of an object, integrating the semantic, geometric and optical flow based constraints into a Dense-CRF model.

In this paper, we first analyze the difficulties of applying a natural language processing model to image segmentation, and improve the spatial structure encoding of LDA through the word-document assignment strategy. Considering the problems of leaf disease spots and mosses, a spatial LDA image segmentation algorithm based on Mean-shift document bandwidths searching (MSBS-LDA) is proposed and applied to plant segmentation. In view of the problems with tomato fruits and the complicated illumination environment, an improved leaf segmentation method based on watershed is proposed to further segment the leaves. Experimental results show that the proposed methods can achieve a high accuracy of image segmentation.

Latent Dirichlet Allocation (LDA)

The Latent Dirichlet Allocation (LDA) was first proposed by Blei, Ng & Jordan (2003). It is a topic generation model which contains a three-level structure of word, topic and document (Griffiths & Steyvers, 2004; David, 2010). Both the topics of a document and the words of a topic obey polynomial distribution. It uses the Bag of Words (BoW) model, where each document is regarded as a word frequency vector, such that textual information is transformed into easy-to-model digital information. Therefore, LDA can be used to identify latent topic information in a large document collection or corpus. In recent years, LDA has been widely used in the field of machine vision as an unsupervised machine learning technology, such as target discovery, scene classification, behavior detection and visual surveillance (Wang & Grimson, 2008). However, spatial structures among visual words, which are vital in computer vision issues, are ignored in the language model since BoW forms a huge gap between low-level visual words and high-level semantics.

LDA gives the topic of each document in a document collection or corpus as a probability distribution. It is an unsupervised learning algorithm, which does not need manual annotation of the training set. Only a corpus and the number of topics are demanded. Note that a document can contain more than one topic, and each word in the document is generated by one of the topics. As depicted in the LDA Bayesian network structure (Fig. 1), the generative procedure of a document in LDA can be described as follows:

(1) For a document i, its topic distribution θ_i is sampled from Dirichlet priori $Dir\left(\alpha \right)$;

(2) For a word j in document i, its topic label z_i,j is sampled from topic multinomial distribution θ_i;

(3) For a topic z_i,j, its word distribution ϕ_zi,j is sampled from Dirichlet priori $Dir\left(\beta \right)$;

(4) The value w_i,j of word j in document i is sampled from word multinomial distribution ϕ_zi,j;

Thus, the joint distribution of all visible variables and latent variables in the model is: (1)${\displaystyle p\left(w_i,z_i,{\theta }_i,\Phi \left|\alpha ,\beta \right.\right)=\prod _{j=1}^{N}p\left({\theta }_i\left|\alpha \right.\right)p\left(z_{i,j}\left|{\theta }_i\right.\right)p\left(\Phi \left|\beta \right.\right)p\left({\omega }_{i,j}\left|{\phi }_{z_{i,j}}\right.\right).}$

The maximum likelihood estimation for a document’s word distribution can be obtained by integrating θ_i, Φ and summing z_i, as shown by: (2)${\displaystyle p\left(w_i\left|\alpha ,\beta \right.\right)={\int }_{{\theta }_i}{\int }_{\Phi }\sum _{z_i}p\left(w_i,z_i,{\theta }_i,\Phi \left|\alpha ,\beta \right.\right).}$

According to the maximum likelihood estimation of $p\left(w_i\left|\alpha ,\beta \right.\right)$, the parameters in the model can finally be estimated by Gibbs sampling.

Spatial LDA

When applied to image segmentation, the modified LDA needs to meet the following two requirements: (1) it should have a good description of local features and (2) it should be equipped with a good spatial summarization ability. The former concentrates on the construction of visual words while the latter focuses on the design of visual documents, making up for the defect of the BoW model. The two processes of local feature description and global structure summarization need to complement each other, leading to the improvement of spatial LDA algorithm.

MSBS-LDA

For the SLDA algorithm (Wang & Grimson, 2008), each document is represented by a point in the image, assuming that its region covers the whole image. If an image patch is close to a document, it has a high probability to be assigned to that document. To describe the word-document mechanism, several parameters are introduced as follows:

(1) d_j: a hidden variable indicates which document word j is assigned to;

(2) $c_i^d$: a hyper-parameter for document i known as a priori;

(3) $g_i^d$: the index of the image where document i is placed;

(4) $\left(x_i^d,y_i^d\right)$: the location of document j;

(5) $c_j=\left(g_j,x_j,y_j\right)$: storing location $\left(x_j,y_j\right)$ and image index g_j of word j;

In the generation procedure of SLDA, for a word j, a random variable d_j is sampled from prior $\mathrm{p}\left(d_j|\eta \right)$ indicating which document word j is assigned to and a uniform prior is used. The image index and location of word j is sampled from distribution $p\left(c_j|c_{d_j}^d,\sigma \right)$, and a Gaussian kernel is chosen: (5)${\displaystyle p\left(\left(g_j,x_j,y_j\right)\left|\left(g_{d_j}^d,x_{d_j}^d,y_{d_j}^d\right)\right.,\sigma \right)\propto {\delta }_{g_{d_j}^d}\left(g_j\right)exp\left\{-\frac{{\left(x_{d_j}^d-x_j\right)}^2+{\left(y_{d_j}^d-y_j\right)}^2}{{\sigma }^2}\right\}}$ where $p\left(c_j|c_{d_j}^d,\sigma \right)=0$ if the word and the document are not in the same image. In the procedure of parameter estimation, z_j and d_j are sampled through a Gibbs sampling procedure integrating out ϕ_k and θ_i. The conditional distribution of z_j given d_j is the same as in LDA, which is given by: (6)${\displaystyle p\left(z_j=k|d_j=i,d_{-j},z_{-j},w,\alpha ,\beta \right)\propto \frac{n_{-j,w_j}^{\left(k\right)}+{\beta }_{w_j}}{\sum _{w=1}^W\left(n_{-j,w}^{\left(k\right)}+{\beta }_w\right)}\cdot \frac{n_{-j,k}^{\left(i\right)}+{\alpha }_k}{\sum _{k^{\prime }=1}^K\left(n_{-j,k^{\prime }}^{\left(i\right)}+{\alpha }_{k^{\prime }}\right)}}$

The method attempts to borrow the language model to image segmentation, for which a uniform prior is applied to determine which document a word is assigned to and a Gaussian kernel is adopted to describe the location of the word. Inspired by this idea, we consider the local maximum points of probability density function (P.D.F) of an image as documents, and the density estimation should be estimated with a nonparametric method due to the fact that the distribution of image data has no fixed pattern.

As a feature space analysis method, Mean-shift (Comaniciu & Meer, 2002; Comaniciu, Ramesh & Meer, 2000) attempts to find the local maximum points of P.D.F in a joint space, and applies nonparametric kernel density estimation (KDE) with smooth effect to density estimation, which provides a new way for word-document assignment of LDA. Thus we propose a modified LDA segmentation algorithm, namely Mean-shift Bandwidths Searching Latent Dirichlet Allocation (MSBS-LDA), for which Mean-shift is adopted to determine which document a word is assigned to. Moreover, the documents are represented by the modes of P.D.F and instead of applying Gibbs sampling to parameter estimation, it applies kernel density estimation to encode image information.

The pixels of each image include two types of information in coordinate space (x^s, $\left[px,py\right]$) and color space (x^r, $\left[l,u,v\right]$) such that the five-dimensional joint feature space $\left[px,py,l,u,v\right]$ is constituted. For a point x, Mean-shift iteratively searches its mode y in the joint space and assigns the color value of the mode to itself, that is x^r = y^r. If we set the footprint of x_i to $\left\{y_{i,0},y_{i,1},y_{i,2},\dots ,y_{i,k},\dots ,y_{i,c}\right\}$ then y_i,0 = x_i at the beginning and it converges to the mode y_i,c. The procedure of mode detection is as follow:

(1) Screen points close to $y_{i,k}^s$ in the coordinate space to the next step. Note that h_s is the kernel function bandwidth in coordinate space.

(2) Use the points survived to calculate the center of gravity and move towards it according to Eq. (7). (7)${\displaystyle y_{i,k+1}^s=\frac{\sum _{n=1}^N{x}_n^{s}g\left({∥\frac{x_n^r-y_{i,k}^r}{h_r}∥}^2\right)}{\sum _{n=1}^{N}g\left({∥\frac{x_n^r-y_{i,k}^r}{h_r}∥}^2\right)}}$ where h_r is the kernel function bandwidth in color space.

(3) Decide if the stopping condition is met or if the number of iterations exceed the maximum limit. If so, stop the search and turn to next step, otherwise, return to step (1) and start from y_k+1. The stopping condition for search is determined by: (8)${\displaystyle \left\{\begin{array}{c}{y}_{i,k+1}^s=y_{i,k}^s\hfill \\ ∥y_{i,k+1}^r-y_{i,k}^r∥\le thr\hfill \end{array}\right.}$

(4) Assign the color $y_{i,c}^r$ of the mode y_i,c to the starting point x_i∕y_i,0, namely $x_i^r\leftarrow y_{i,c}^r$.

For density estimation at point x, the non-parametric density estimation is computed as: (9)${\displaystyle \hat{f}\left(x\right)=\frac{1}{N{h}^d}\sum _{n=1}^{N}K\left(\frac{x_n-x}{h}\right)=\frac{1}{N{h}^d}\sum _{n=1}^N{c}_{k,d}k\left({∥\frac{x_n-x}{h}∥}^2\right)}$ where $\mathrm{K}\left(x\right)$ is the kernel function, for which the radially symmetric function is used here such that $\mathrm{K}\left(x\right)=c_{k,d}k\left({∥x∥}^2\right)$. Note that, c_k,d is a normalization constant, making ${\int }_{R^d}K\left(x\right)dx=1$, and $\mathrm{k}\left(x\right)$ is the profile for $\mathrm{K}\left(x\right)$. The Epanechnikov Kernel is adopted and its profile is given by: (10)${\displaystyle \begin{array}{c}k\left(x\right)=\left\{\begin{array}{cc}1-x\hfill & 0\le x\le 1\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right.\hfill \\ K_E\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2}\left(d+2\right)\left(1-{∥x∥}^2\right)\hfill & ∥x∥\le 1\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right.\hfill \end{array}}$

In addition, the density gradient estimation for point x is computed as: (11)${\displaystyle \hat{\nabla }{f}_{h,K}\left(x\right)=\frac{2c_{k,d}}{N{h}^{\left(d+2\right)}}\sum _{n=1}^N\left(x_n-x\right)\left[-k^{\prime }\left({∥\frac{\left(x_n-x\right)}{h}∥}^2\right)\right]}$

We know where to move and then the step length is determined by: (12)${\displaystyle \begin{array}{ccc}\hat{\nabla }{f}_{h,K}\left(x\right)\hfill & =\hfill & \frac{2c_{k,d}}{N{h}^{\left(d+2\right)}}\sum _{n=1}^N\left(x_n-x\right)\left[-k^{\prime }\left({∥\frac{\left(x_n-x\right)}{h}∥}^2\right)\right]\hfill \\ \hfill & =\hfill & \frac{2c_{k,d}}{N{h}^{\left(d+2\right)}}\sum _{n=1}^N\left(x_n-x\right)g\left({∥\frac{\left(x_n-x\right)}{h}∥}^2\right)\hfill \\ \hfill & =\hfill & \frac{2c_{k,d}}{N{h}^{\left(d+2\right)}}\left[x_n\sum _{n=1}^{N}g\left({∥\frac{\left(x_n-x\right)}{h}∥}^2\right)-x\sum _{n=1}^{N}g\left({∥\frac{\left(x_n-x\right)}{h}∥}^2\right)\right]\hfill \\ \hfill & =\hfill & \underset{C}{\frac{2c_{k,d}}{h^2{c}_{g,d}}}\underset{{\hat{f}}_{h,G}\left(x\right)}{\left[\frac{c_{g,d}}{N{h}^d}\sum _{n=1}^{N}g\left({∥\frac{\left(x_n-x\right)}{h}∥}^2\right)\right]}\underset{m_{h,G}\left(x\right)}{\left[\frac{\sum _{n=1}^N{x}_{n}g\left({∥\frac{\left(x_n-x\right)}{h}∥}^2\right)}{\sum _{n=1}^{N}g\left({∥\frac{\left(x_n-x\right)}{h}∥}^2\right)}-x\right]}\hfill \end{array}}$ where ${\mathrm{m}}_{h,G}\left(x\right)$ is the mean shift vector. It is a directional vector, which corresponds to the gradient direction. For a point x moving to x’ in the gradient direction, the new coordinate is: (13)${\displaystyle x^{\prime }=x+m_{h,G}=\frac{\sum _{n=1}^N{x}_{n}g\left({∥\frac{x_n-x}{h}∥}^2\right)}{\sum _{n=1}^{N}g\left({∥\frac{x_n-x}{h}∥}^2\right)}.}$

As we can see from the above steps, the document assignment of word i is hugely affected by the bandwidths h_s and h_r of multivariate kernel density estimation. When we fix h_s and increase h_r monotonically within a certain range, the number of clusters decreases. A large h_r leads to inadequate documents (regions) while a small h_r results in excessive documents (regions). Similarly, the number of clusters decreases when h_r is fixed and h_s increases. In this regard, we propose to use the word frequency statistics and structural features of the image to find the required parameters $\left(h_s,h_r\right)$. According to the meaning of the coordinate scale parameter h_s, we count the word frequency w_fren of the quantified visual words in accordance with the LDA dictionary. Then h_s is calculated on the basis of the proportion of the first major word by Eq. (14) where C is a constant in [1.41, 1.67]. For the color scale parameter h_r, we iteratively call the LDA segmentation procedure based on Mean-shift documents to obtain an optimal result. The mean value of LUV distance between the K classes is calculated by Eq. (15) and the last three $\Delta \bar{E}$ serve as a cut-off condition for the iteration of LDA. (14)${\displaystyle h_s=10^C\times \frac{w_{fren}^{\left(1\right)}}{N}}$ (15)${\displaystyle \Delta \bar{E}=\frac{2\sum _{i=1}^{K-1}\sum _{j=i+1}^K\sqrt{{\left({\bar{L}}_i-{\bar{L}}_j\right)}^2+{\left({\bar{U}}_i-{\bar{U}}_j\right)}^2+{\left({\bar{V}}_i-{\bar{V}}_j\right)}^2}}{K\left(K-1\right)}.}$

Under the framework of LDA, the number of topics K is determined manually. After setting h_s, we initialize h_r and decrease it with a step of one. In the beginning, Mean-shift obtains few clusters where the structure information is not obvious and changes a lot. Hence, LDA based on these documents leads to mis-segmentation, and the resulting $\Delta \bar{E}$ between the K classes is unstable. During the implementation process, h_r gradually decreases and the number of clusters increases, and under-segments change to over-segments as a result. Thereafter, the mis-segmentation situation of LDA is significantly improved, and $\Delta \bar{E}$ tends to stabilize after a big jump. If h_r is quite small, Mean-shift turns up excessive segments and the structural information is destroyed. Therefore, we aim to find a set of over-segmented Mean-shift documents, which should maintain certain structural description. The qualitative estimate of the value of h_r should fall in a middle range.

To obtain the final segmentation result, we need to consider the stopping criterion for the iteration of LDA. Note that, when under-segments change to over-segments, the Mean-shift documents have a great impact on LDA segmentation result. When h_r falls in a suitable middle range, the LDA segmentation result tends to be stable, and if h_r continues to decrease, the segmentation result is adversely effected. Therefore, we need to find the segmentation result immediately when the stable state is reached. We take it as the condition of convergence that the difference of three consecutive $\Delta \bar{E}$ is less than one. But sometimes there may occur special cases where two consecutive $\Delta \bar{E}$ are close and then spread out or the optimal range of the LDA segmentation result is limited. At this time, we take the segmentation corresponding to the first large jump of $\Delta \bar{E}$ as the optimal result to output. In addition, when h_r decreases within 10, the number of Mean-shift clusters increases rapidly, which leads to an “overleaping” on the number of spatial documents, thus affecting the judgement of convergence. Therefore, when h_r is less than 10, we change the step length to 0.5 to reduce the adverse effect of sudden changes in the document number on the segmentation results.

Based on the calculation above, the topic probability distribution for each pixel can be obtained by LDA given the word-document assignment. The generation process is as follows:

(1) For a document i, multinomial parameter θ_i for K topics is sampled from Dirichlet priori, ${\theta }_i\sim Dir\left(\alpha \right)$;

(2) For a topic k, multinomial parameter ϕ_k is sampled from Dirichlet priori, ${\varphi }_k\sim Dir\left(\beta \right)$;

(3) For a word j, its document d_j is determined by the search of kernel function bandwidth $h=\left(h_s,h_r\right)$ of Mean-shift;

(4) For a word j, its topic label z_j is sampled from discrete distribution of document d_j, $z_j\sim Discrete\left({\theta }_{d_j}\right)$;

(5) The value w_j of word j is sampled from discrete distribution, $w_j\sim Discrete\left({\varphi }_{z_j}\right)$.

Based on the analysis above, MSBS-LDA makes an improvement to Assumption 2 and proposes Assumption 3: If two image patches are from the same object class, they often appear in one image, and they belong to the pre-classified document in image space. That is to say, word i is likely to be assigned to document j if satisfying the following conditions: (1) they are in the same image, (2) word i belongs to the same document as its mode, and (3) the mode is found by Mean-shift combined with word frequency statistics and LDA iterative search. The assumption it adds is that a word tends to have the same topic label as other words in the pre-classified document of image space.

To overcome the limitations of two word-document assignments proposed in Sect. ‘Visual documents’, MSBS-LDA makes some breakthroughs in the two major difficulties of document design: image data modelling and document tolerance.

(1) Spatial information of the image is modelled by nonparametric estimation and the filtering effect of Mean-shift enhances the tolerance of the documents to feature differences, leading to a good global summary ability of the algorithm.

(2) Consider that the bandwidth parameters need to be manually adjusted for spatial documents, which is a non-trivial issue in case of real-world problems without domain knowledge. We propose to use word frequency statistics to determine the coordinate space bandwidth parameter h_s and iteratively search the optimal color space bandwidth parameter h_r in combination with the LDA segmentation algorithm of Mean-shift documents. The optimization of word-document assignment can reduce the time and resource consumption and greatly improve the accuracy of LDA segmentation result.

With the help of the Gaussian filter bank, we increase the diversity of words to ensure the local description ability of the LDA algorithm along with a rapid convergence of the unsupervised learning dictionary. The word-document assignment based on Mean-shift improves the image data modelling of documents and the tolerance to different features of the same object, thereby improving the global summary ability of the spatial LDA algorithm. Based on LDA’s word frequency statistics and structural features of segmentation result, Mean-shift is guided to automatically find the optimal bandwidth parameters. The combination of the three processes improves the ability of LDA to solve image segmentation problems. Fig. 2 shows the flowchart of the MSBS-LDA algorithm.

Plant Segmentation and Leaf Segmentation

The segmentation of greenhouse plant images is divided into two processes: (1) plant extraction from background and (2) single leaf segmentation. The superiority of MSBS-LDA to describe local details in a complex scene and its tolerance to differences in features of the same object make it suitable for the problem of plant segmentation under the greenhouse scene, especially in the case of complex lighting, leaf lesions and mosses. Therefore, we first adopt MSBS-LDA to handle plant segmentation, extracting plants from the complex scene.

The illumination in greenhouses is ever-changing and the images captured during the fruits bearing period sometimes include tomato fruits of various growth stages, which have a negative impact on leaf segmentation. In these regards, we first compute color difference (R-B) on plant foreground with fruit and then perform threshold segmentation to remove the fruit part. Thereafter, the lighting environment is investigated based on the gray distribution of the leaf foreground histogram, as shown in Fig. 3. The grayscale distribution rate α′ is calculated as: (16)${\displaystyle {\alpha }^{\prime }=\frac{C_{hist}\left(M_s-R_s,M_s+R_s\right)}{C_{hist}\left(1,255\right)},R_s=255\times {\beta }^{\prime }∕2}$ where C_hist is a grayscale statistical function, M_s is the maximum gray value of the accumulated grayscale, R_s is a grayscale radius, and β′ represents a grayscale distribution range. We set β′ and the threshold of α′ as 0.2 and 0.8, respectively. Equation (16) denotes that if more than 80% of the pixels are distributed within 20% of the grayscale range, the lighting environment of the foreground leaf image is considered to be ideal, otherwise, illumination correction of homomorphic filtering is required to weaken its effect.

Then, leaf segmentation is carried out on the extracted leaf part after illumination response in three steps: strong edge detection, leaf centroid location, and marker-controlled watershed segmentation. Removal of strong edges from foreground leaf mask is a preprocessing step for leaf center detection, which can significantly improve the efficiency and accuracy of local maximum filtering. In the literature, some scholars have provided new ideas for closed contour effects (Valliammal & Geethalakshmi, 2011; Sampath & Shan, 2007; Dollár & Zitnick, 2015; Ming, Li & He, 2012). To decrease algorithm complexity and meet the demand of rapidity, the Structured Edge (SE) detector (Dollár & Zitnick, 2015) is applied to obtain leaf-leaf boundaries. Then we subtract the detected edges from the foreground leaf mask and compute a Euclidean distance map on it. Afterwards, the local maximum of the distance map is searched to locate each leaf’s centroid based on the dilation operation. Further segmentation by watershed algorithm is conducted as the post-process step. The schematic diagram of the segmentation process is shown in Fig. 4 and the overall flowchart of the proposed plant segmentation and leaf segmentation method is shown in Fig. 5.

Experiments and Analysis

In our research, we first test MSBS-LDA on the images from the Microsoft Research Cambridge Object Recognition Image Database v2 (MSRC-v2) (Microsoft, 2000) to evaluate its performance on the general dataset. For the segmentation of plant and leaf, all the tomato plant images were taken under real greenhouse conditions from three Venlo greenhouses of the Chongming Base of National Facility Agricultural Engineering Technology Research Center (abbreviated as CM), the Sunqiao Modern Agricultural Development Zone in Shanghai (abbreviated as SQ), and the Jiading Experimental Greenhouse in Tongji University (abbreviated as JD). The main characteristic of the CM images is that they have uniform illumination and the plant is in vegetal stage with smooth leaves; while for the SQ images, the plant is in the stage of blossoming and bearing fruits, and the complicated light environment causes reflection as well as shadows on the leaves. The leaves with light or serious disease spots of the JD images increase the diversity of the plant images. It deserves pointing out that all the ground-truths and training set for comparison experiments were labelled manually by the first author. In addition to the tomato plant images, we took 13 images (among which, nine images contain obvious mosses) of Arabidopsis from the Computer Vision Problems in Plant Phenotyping (CVPPP) dataset A1 subset (Scharr et al., 2014; Minervini, Abdelsamea & Tsaftaris, 2014) to test the ability of MSBS-LDA to handle mosses. All experiments were conducted on a PC HP-g4-1059TX machine (Shanghai, China) with 2.10 GHz CPU and 6GB RAM.

Evaluation metrics

In order to quantitatively evaluate the accuracy of the spatial LDA segmentation algorithms including NR-LDA, OR-LDA, SLIC-LDA and MSBS-LDA, we define the following three metrics:

(1) SA (Segmentation Accuracy) measures the area of overlap between ground-truth and algorithm result: (17)${\displaystyle SA=\frac{1}{K}\sum _{i=1}^K\frac{\left|P_i^{gt}\cap P_i^{ar}\right|}{\left|P_i^{gt}\right|}}$

(2) OR (Over-segmentation Rate) measures the area of algorithm result not in ground-truth: (18)${\displaystyle OR=\frac{1}{K}\sum _{i=1}^K\frac{\left|P_i^{ar}\right|-\left|P_i^{gt}\cap P_i^{ar}\right|}{\left|P_i^{gt}\right|}\epsilon \left(\left|P_i^{ar}\right|-\left|P_i^{gt}\cap P_i^{ar}\right|\right)}$

(3) UR (Under-segmentation Rate) measures the area of ground-truth not in algorithm result: (19)${\displaystyle UR=\frac{1}{K}\sum _{i=1}^K\frac{\left|P_i^{gt}\right|-\left|P_i^{gt}\cap P_i^{ar}\right|}{\left|P_i^{gt}\right|}\epsilon \left(\left|P_i^{gt}\right|-\left|P_i^{gt}\cap P_i^{ar}\right|\right)}$ where $P_i^{gt}$ and $P_i^{ar}$ denote the ground-truth and the algorithm result of the ith class, and ε(⋅) is the step function. All of these metrics are counted in pixels, with larger values of SA, smaller values of OR and UR representing higher agreement between ground-truth and algorithmic results.

The metrics above are suitable for the evaluation of general multi-class classifiers, but for plant and leaf segmentation, more targeted quantitative evaluation metrics are required. In this regard, some evaluation criteria proposed in Scharr et al. (2016) are adopted here, which are, respectively, defined as follows:

(1) FBD (Foreground-Background Dice) is Dice score of foreground plant: (20)${\displaystyle Dice\left(\text{\%}\right)=\frac{2\left|P^{gt}\cap P^{ar}\right|}{\left|P^{gt}\right|+\left|P^{ar}\right|}}$ which measures the degree of overlap between ground-truth P^gt and segmentation result P^ar.

(2) SBD (Symmetric Best Dice) is symmetrical mean Dice of all leaves, as computed by: (21)${\displaystyle SBD\left(L^{ar},L^{gt}\right)=min\left\{BD\left(L^{ar},L^{gt}\right),BD\left(L^{gt},L^{ar}\right)\right\}}$ where BD (Best Dice) is defined as: (22)${\displaystyle BD\left(L^a,L^b\right)=\frac{1}{M}\sum _{i=1}^M\underset{1\le j\le N}{max}\frac{2\left|L_i^a\cap L_j^b\right|}{\left|L_i^a\right|+\left|L_j^b\right|}}$

(3) Dic (Difference in Count) is difference in leaf number between ground-truth and algorithm result: (23)${\displaystyle Dic=\#L^{ar}-\#L^{gt}}$

(4) —Dic— is absolute value of Dic,

(5) MHD (modified Hausdorff distance) measures accuracy of shapes and boundaries by comparing two sets of points, A and B, on the edge of a leaf: (24)${\displaystyle MHD\left(A,B\right)=max\left\{D\left(A,B\right),D\left(B,A\right)\right\}}$ with (25)${\displaystyle D\left(A,B\right)=\frac{1}{\left|A\right|}\sum _{p\in A}\underset{q\in B}{min}∥p-q∥.}$

Spatial LDA evaluation

The experiment is carried out on the MSRC-v2 dataset to test the spatial LDA segmentation algorithms including NR-LDA, OR-LDA, SLIC-LDA and MSBS-LDA.

First, we investigate the influence of different word-document assignment strategies on the performance of spatial LDA. The segmentation results are shown in Figs. 6 and 7. We can figure out that the traditional LDA takes local descriptors as words and the whole image as one document, then it clusters the visual words that often appear in the same image into one object class, which exists two problems: (1) the segmentation result is noisy because spatial information is not considered and (2) since the whole image is treated as one document, if an object is dominant in the image, then other non-dominant objects could be labelled as the dominant one. NR-LDA improves the effect of noise, but documents of non-overlapping regions make rectangle edges noticeable. This kind of borderline situation is further improved by OR-LDA. Compared with spatial structure, SLIC-LDA more meticulously depicts the homogeneity of the atomic regions and the differences between them, so that weakens the tolerance of the documents to the different features in the same object. As a result, some unnecessary edges and noise in the segmentation results are generated. With the description of color and location information of the documents, MSBS-LDA tends to get better segmentation results in both details and global structure. The smooth effect of kernel density estimation and filtering effect of Mean-shift greatly reduce the noise of LDA caused by the BoW model. The results are no longer limited by the size of the rectangular documents, also, the shapes and boundaries are more natural. The evaluation results are shown in Table 1, where MSBS-LDA achieves the highest segmentation accuracy (0.9899), the lowest under-segmentation rate (0.0787) and the lowest over-segmentation rate (0.0101).

Table 1:

Segmentation evaluation of spatial LDA algorithms for images in Fig. 7.

	SA	OR	UR
LDA	0.7483	0.7072	0.2517
NR-LDA	0.8582	0.7885	0.1419
OR-LDA	0.9509	0.7998	0.0491
SLIC-LDA	0.8976	0.7342	0.1024
MSBS-LDA	0.9899	0.0787	0.0101

DOI: 10.7717/peerj.5036/table-1

Afterwards, we compare MSBS-LDA with some unsupervised learning segmentation methods. The segmentation results are shown in Fig. 8. It can be seen that MSBS-LDA shows superior performance in general. Compared with 2D-OTSU and ICM, MSBS-LDA can significantly improve the noise situation, and it also has better tolerance to feature differences of the same object. Meanwhile, compared with Normalized-cut and Co-segmentation, the description of the details of MSBS-LDA has not been weakened by the tolerance.

Plant segmentation evaluation

For three kinds of tomato plant images with distinct characteristics, we adopt MSBS-LDA to fulfill the plant segmentation task. Meanwhile, we test some other segmentation algorithms and features in comparison: (1) BP with 3D color and texture features proposed in Minervini, Abdelsamea & Tsaftaris (2014), (2) OTSU with a* channel color of SLIC super-pixels, (3) co-segmentation, (4) dense-CRF based on texton-boost (Krähenbühl & Koltun, 2011; Yang & Li, 2014). The results of the five foreground segmentation methods on three tomato plant image sets are illustrated in Fig. 9, and the quantitative evaluation is shown in Tables 2 and 3. For Arabidopsis images from CVPPP A1 subset, the results of MSBS-LDA are illustrated in Fig. 10, and the quantitative evaluation is shown in Table 4.

According to the experimental results on tomato plant images, for the CM images, MSBS-LDA obtains the segmentation result comparable to BP and SLIC-OTSU, the stems and other details can be well separated from the background. The performance of MSBS-LDA is not affected by the complex background of the SQ image, and there is almost no noise in the result. Moreover, MSBS-LDA shows superior segmentation performance on the JD image which is most difficult for plant extraction. From Fig. 9 it is clear that not only the plant and the background can be effectively separated, but also the leaf disease spots can be accurately included in the foreground plant. When compared with the other four segmentation methods, MSBS-LDA achieves better segmentation accuracy, leaf shape and edge precision in general. In fact, through experiments we also find that gmentation performance of some methods, such as BP and co-segmentation which rely on the training set or the same type of objects in a group, could be greatly affected if the characteristic between images are quite different or the design of the training set is not appropriate. In addition, from Table 5 we can assume that the time consumption of MSBS-LDA is competitive and practicable for most greenhouse applications. From the experimental results on CVPPP A1 subset, we can see that MSBS-LDA is able to separate mosses and plants effectively by adjusting the topic number. The appearance of mosses does not affect the performance of the algorithm much.

Table 2:

FBD of five methods on three tomato plant image sets.

FBD (for plant, %)	BP	SLIC-OTSU	Co-segmentation	Dense-CRF	MSBS-LDA
CM1	97.90	98.50	93.00	92.76	98.08
CM2	97.83	97.65	89.44	80.65	97.54
CM3	97.31	97.12	82.04	90.47	98.02
SQ1	96.32	96.03	92.44	85.78	96.42
SQ2	97.07	96.93	88.99	92.63	97.44
SQ3	98.80	99.13	97.15	96.18	99.18
JD1	89.17	86.65	89.07	88.85	90.78
JD2	94.00	92.02	92.89	92.62	94.31
JD3	84.93	90.82	84.31	72.19	96.54

DOI: 10.7717/peerj.5036/table-2

Table 3:

MHD of five methods on three tomato plant image sets.

MHD (for plant)	BP	SLIC-OTSU	Co-segmentation	Dense-CRF	MSBS-LDA
CM1	0.71	0.48	5.92	7.04	1.48
CM2	0.46	0.72	1.89	3.14	0.44
CM3	1.51	2.51	6.90	6.41	2.03
SQ1	4.88	6.88	9.90	8.22	3.51
SQ2	6.43	4.17	14.79	5.05	2.10
SQ3	3.83	2.23	4.66	2.13	0.45
JD1	8.80	9.93	10.36	8.12	8.75
JD2	3.92	5.17	3.71	4.36	3.02
JD3	7.67	4.69	6.86	32.55	2.35

DOI: 10.7717/peerj.5036/table-3

Table 4:

Scores of MSBS-LDA on CVPPP A1 subset.

Samples	FBD (%)	MHD
13	95.99	0.51
59	95.69	0.62
91	95.70	0.98
96	95.55	0.70
115	95.24	0.70
153	94.54	0.88

DOI: 10.7717/peerj.5036/table-4

Table 5:

Time complexity of five methods for each image.

Methods	BP	SLIC-OTSU	Co-segmentation	Dense-CRF	MSBS-LDA
Time consumption (s)	22.51	5.51	213.14	166.90	141.64

DOI: 10.7717/peerj.5036/table-5

Leaf segmentation evaluation

For the plant segmentation results containing fruit part, the separation results of leaves and fruits are shown in Fig. 11. The leaf segmentation results on three tomato plant image sets and the CVPPP A1 subset are shown in Figs. 12 and 13, respectively. The evaluation scores for all the testing images are listed in Table 6. For more complicated cases, like the JD images with lesions and the SQ images with both complex lighting environment and seriously overlapped leaves, our method reaches segmentation accuracy ranging from 50% to 60%. For the CM images with more uniform illumination and smooth leaves, the accuracy becomes more than 65%. The leaf segmentation scores of CVPPP A1 images are listed in Table 7. For Arabidopsis images with mosses and weak boundaries between overlapping leaves, the segmentation accuracy is up to 70%.

In this study, we find that if centroids of leaves are not detected correctly owing to serious overlapping, complicated lighting environment or side-effect of lesions, the performance of watershed algorithm will be limited. Besides, low contrast, weak boundaries and complex posture can lead to low segmentation accuracy. The method can also be fine-tuned to achieve a higher accuracy on one particular image, such as utilizing a more accurate edge detector, a more appropriate illumination normalization module or a better local maximum filter.

Table 6:

Leaf segmentation evaluation on three tomato plant image sets.

Image sets	FBD (for leaves, %)	SBD (%)	|DiC |	DiC
CM	96.77 (1.83)	65.37 (6.25)	1.80 (1.32)	−1.00 (2.05)
SQ	96.75 (1.71)	55.86 (4.00)	2.60 (2.17)	−2.60 (2.17)
JD	93.73 (2.79)	51.88 (7.16)	1.25 (0.89)	1.25 (0.89)
All	95.89 (2.46)	58.12 (8.03)	1.93 (1.63)	−0.93 (2.37)

DOI: 10.7717/peerj.5036/table-6

Notes:

Average values are shown for metrics described in ‘Evaluation metrics’ and in parenthesis standard deviation.

Table 7:

Leaf segmentation evaluation on CVPPP A1 subset.

	FBD (%)	SBD (%)	—DiC—	DiC
Mean	94.42 (1.48)	73.81 (5.32)	2.62 (1.33)	−2.62 (1.33)
Median	95.24	74.71	3.00	−3.00
Max	95.99	82.40	5.00	−1.00
Min	92.03	66.12	1.00	−5.00

DOI: 10.7717/peerj.5036/table-7

Discussion and Conclusion

In this study, we propose a modified statistical model of LDA, namely MSBS-LDA, to segment greenhouse tomato plants in an unsupervised way, and leaf segmentation is carried out subsequently. Through our experiments in different cases, some conclusions are drawn as follows:

(1) After analyzing the difficulties of using natural language processing model in image segmentation, it is proposed to improve the LDA algorithm in the spatial structure encoding of images through the word-document assignment. The diversity of visual vocabulary is guaranteed by constructing visual words, which can improve the ability of LDA to describe image details. Moreover, the single document mapping is changed by designing visual documents, ameliorating the spatial summarization ability of LDA. The comparison experiments show that the spatial LDA algorithms with word-document assignments outperform the traditional one in image segmentation.

(2) In regards to leaf lesion spots and complicated backgrounds of greenhouses, a spatial LDA segmentation algorithm based on bandwidths searching of Mean-shift documents (MSBS-LDA) is proposed and applied to plant segmentation. The non-parametric estimation is adopted to space modelling, and the modes of Mean-shift are mapped as documents. The comparison experiments show that the proposed MSBS-LDA algorithm can simultaneously give good expression to image details and the global structure information, so that it can accurately include the lesion part in the plant and separate the plant from mosses of similar color through the adjustment of the topic number.

(3) In regards to tomato fruits and ever-changing illumination in greenhouses, an improved watershed segmentation method is proposed and applied to leaf segmentation. The fruit separation, illumination normalization, strong edge detection, leaf centroid location, and marker-controlled watershed segmentation are carried out sequentially to complete leaf segmentation. The experiments show that the method can guarantee a certain segmentation accuracy for three greenhouse tomato plant image sets of different characteristics.

Supplemental Information