SINGLE TREE DETECTION FROM AIRBORNE LASER SCANNING DATA USING A MARKED POINT PROCESS BASED METHOD

Tree detection and reconstruction is of great interest in large-scale city modelling. In this paper, we present a marked point process model to detect single trees from airborne laser scanning (ALS) data. We consider single trees in ALS recovered canopy height model (CHM) as a realization of point process of circles. Unlike traditional marked point process, we sample the model in a constraint configuration space by making use of image process techniques. A Gibbs energy is defined on the model, containing a data term which judge the fitness of the model with respect to the data, and prior term which incorporate the prior knowledge of object layouts. We search the optimal configuration through a steepest gradient descent algorithm. The presented hybrid framework was test on three forest plots and experiments show the effectiveness of the proposed method.


INTRODUCTION
There has been a great deal of interest in the 3D modelling of urban and suburban environments. Vegetation detection and reconstruction plays an important role in the task, as they present in forms of row trees, tree groups and in largely, parks and urban forests. Serving as basis of tree modelling and many other applications, single tree detection is a topic which has received significant interest from the remote sensing community over the past decades. The predominant studies on single tree detection are mainly in forest inventory, however, the methods developed in forest environments can be migrated to detect trees in urban forests or forests near power line corridorswhich are of great significance in 3D urban and infrastructure modelling, planning and management.

Low-level Image Processing Techniques
Numerous methods have been proposed to detect single trees from ALS data. Most of them focus on the reconstruction of the canopy height model (CHM), which provides a representation of outer geometry of tree canopy. The peaks and valleys on CHM are assumed to be better estimations of treetop positions and crown edges than that can be obtained from aerial photos or satellite imageries. For this reason, many studies have extended methods developed for passive optical imageries to detect single trees from ALS data. Those methods include but not limited to: local maxima filtering (Popescu et al., 2002), region growing (Solberg et al., 2006), valley following (Leckie et al., 2003), template matching (Korpela et al., 2007), watershed segmentation and its variance marker-controlled watershed segmentation (Chen et al., 2006;Pyysalo and Hyyppa, 2002), and multi-scale segmentation (Brandtberg et al., 2003).
Among the proposed methods, local maxima filtering and marker-controlled watershed segmentation are most commonly used and ready for operational application, for their fast implementation, while being able to produce relatively plausible results (Kaartinen et al., 2012). Hyyppa et al. (2001) detected treetops from CHMs using local maxima filtering with fixed window size. Popescu et al. (2002) then tested a variable window local maxima filtering on CHMs attempting to overcome errors of omission and commission associated with fixed window local maxima filtering. This approach requires a priori knowledge of the relationship between tree height and crown size and the detection accuracy can greatly influenced by the specification of the relationship. In many cases, this relationship is either hard to obtain or different from study to study. Moreover, Falkowski et al. (2006) pointed out the relationship between tree height and crown size can be weak in some forest conditions. Further complicating the issue is that, in those applications, CHMs were usually smoothed with Gaussian filter before local maxima detection. Inappropriate specification of the size of the low pass filter and the degree of smoothing can also decrease the detection accuracy even before local maxima filtering is carried out in the subsequent step.
Once the treetops are detected, marker-controlled watershed segmentation is well suited to delineate the tree crown segments from the CHM. Marker-controlled watershed segmentation, introduced by Meyer and Beucher (1990), imposes the advantages of other segmentation methods of region growing and edge detection, thus being able to overcome the oversegmentation problem of ordinary watershed segmentation. In marker-controlled watershed segmentation, user specified markers will replace the local maxima that usually spread all over the image to perform the segmentation. In the resultant segmentation, one segment will correspond to each marker: in the case of single tree detection, one tree crown will be captured by one treetop. That also means the detection accuracy of marker-controlled watershed segmentation subjects to the accuracy of pre-determined local maxima as true treetops in the previous stage. Most local maxima based approaches tend to produce high commission errors in coniferous forests, as spurious treetops are detected within tree crowns from large branches. Otherwise, in cases local maxima filtering produces low commission error, omission error often goes up as small tree crowns are more likely to be undetected (Gebreslasie et al., 2011;Kaartinen et al., 2012).

Stochastic Methods
Stochastic methods represent another branch of powerful tools in image analysis. They have proven to hold great promise in solving inverse problems including image segmentation, image restoration, and feature extraction (Descombes and Zerubia, 2002) .
Markov random field (MRF) was introduced into computer vision community by the seminal works of Besag (1986) and Geman and Geman (1984). MRF models an image as a realization of a collection of random variables associated with each pixel. The appeal of MRF is that it provides a probabilistic framework to encode contextual constrains into the prior probability, making it robust with respect to noise in the image (Li, 2009). However, the local definition of pixelwise constrains in MRF is difficult to incorporate more global and strong geometric constrains (Descombes and Zerubia, 2002). Marked point process models can be seen as an extension of MRF (Mallet et al., 2010), such that random variables are associated not with each pixels in the image but with random configurations of geometric objects or shapes describing the image. This means marked point process could model higher level geometrical primitives more naturally, while inheriting the merit of inclusion of priori knowledge on spatial pattern of features.
Several stochastic models have been proposed to detect tree crowns from remote sensing data. Descombes and Pechersky (2006) presented a three state MRF model to detect tree crown from aerial image. This approach addressed the problem as an image segmentation problem and worked on pixel level. Each pixel has three states as vegetation, background and centre of trees. Although the MRF was defined on pixel level, the label update was performed on object level using elliptical templates of crowns. Zhang and Sohn (2010) proposed a high-level MRF model for tree detection from ALS recovered CHM, by modelling trees as objects using low-level image features. Furthermore, Perrin et al. (Perrin et al., 2005(Perrin et al., , 2006 employed marked point processes to detect tree crowns in plantations from near infrared (NIR) aerial imageries. Tree crowns in the remote sensing image are modelled as a configuration of discs or ellipses. In all works, tree crowns were detected by maximizing a Bayesian criterion such as maximum a posteriori (MAP), which ended up as an energy minimization problem and solved in a simulated annealing framework.

Motivation
Marked point processes offer a powerful tool for the inclusion of constraints in the prior at object level. But the model is linked with images via a data term measured using hypothesis testing schema. This inherited property of stochastic model requires explorations of huge configuration spaces in order to find the optimal configuration, especially for non-data-driven model which does not make use of any low-level information that can be extracted from images. The optimization process of such model is usually computationally expensive and takes a long time to converge.
Methods have been proposed to overcome this drawback in traditional marked point process (Descombes et al., 2009). This paper proposes to combine the virtue of marked point process and traditional image process techniques, and presents a hybrid framework to detect single trees from ALS data.

MARKED POINT PROCESS
Marked point processes have been introduced in image processing by Baddeley and Van Lieshout (1993) as powerful stochastic tools to detect an unknown number of objects in an image.

Definition
Let be a marked point process in space , a bounded set of . Each object of is defined by its position in and associated with a mark from . is a measurable mapping from a probability space ( ) to configurations of points of . In another word, is a random variable whose realizations are random configurations of geometric objects. The configurations, noted as in the following, belongs to a space ( ( ) ( )), where contains all configuration of a finite number of object of .

A Model for Tree Crown
In remote sensing images, distribution of tree crowns in forests can be represented by a marked point process of disks. The associated space can be written as: where and are respectively the width and height of the image , and ( , ) the minimum and maximum radius of the disks in the configuration. Note ( ) as a tree object, where is its position and its radius. The configuration space of marked point process of tree crowns in a remote sensed image can be written as:

Energy of the Model
We consider the probability distribution ( ) of a homogeneous Poisson Process on , which give us a probability measure ( ) on . Then, for a marked point process , we can define a density ( ) with respect to the reference measure as: The density of a configuration can be written in Gibbs form: where ∫ ( ) ( ) is a normalizing constant.
Feature extraction or object detection from images with the model then turns into finding the configuration ̂ maximizing the posterior probability. This issue is then reduced to an energy minimization problem, finding the configuration minimizing the Gibbs energy ( ), i.e., ̂ ( ), which is equivalent to finding the maximum density estimator ̂ ( ).

PROPOSED MODEL
In this work, we seek to find a constraint configuration space in which the optimal or near optimal configuration resides. We will then limit the search for optimal configuration in the constraint space , which could greatly reduce the heavy computation demand of random sampling in in optimization process. This proposed model is realized by combining lowlevel image processing techniques of local maxima filtering and marker-controlled watershed segmentation with marked point process.

Constrained Configuration Space
The CHM image, representing the height of tree crowns above ground, is firstly constructed from classified ALS data. Then, we apply local maxima filtering with variable window size on the CHM to extract local maxima as potential treetops, which serves as good representation of tree locations. Let represents the set of extracted local maxima: where is the total number of local maxima extracted.
For a given subset of local maxima , they are used as markers in marker-controlled watershed segmentation to obtain a partition ( ) { ( ) } of the CHM, where is the corresponding segment of local maxima . ( ) is a low level presentation of the CHM image, and the set of segments are assumed to be a reasonable approximation of tree crowns with respect to the set of local maxima , where ( ) is number of local maxima in .
A tree object ( ) is then defined using location and radius on segment , where tree location is corresponding to the local maximum and radius is calculated as the average radius of segment . A configuration ( ) { ( ) } can be then constructed from the set of local maxima . The whole procedure is illustrated in Figure 1. We note all configurations generated from a subset of local maxima set as ⋃ ( ) . Apparently, is a discrete subspace of the whole configuration space . In this way, we build a constrained configuration space for the marked point process to sample the optimal configuration from.

Energy Formulation
A configuration of tree objects is measured in terms of goodness or cost by the global energy ( ), which is composed of a data term ( ) and a prior term ( ) such that: where [ ] is used to tune the tradeoff between the data term and the prior term.

Data Term:
This term is used to indicate the likelihood of the tree objects in relation with low-level segments obtained from the CHM image. Two energy functions are proposed to measure how well those features support the object as a plausible tree.

(i) Symmetric Function ( )
A symmetric ratio is defined as a measure of how well a treetop is located in the central part of the crown and the degree of the tree crown being of a symmetric circular shape (see Figure 2). For a given tree object with corresponding segment , the average radius of the segment and its standard derivation are calculated and noted as ( ) and ( ). The asymmetric ratio ( ) [ ] of object is calculated as: A sigmoid function is then used to define the symmetric function to penalize asymmetric tree crowns given by Eq. (7): where and are parameters set to control the position and slope of the sigmoid function respectively. The larger the asymmetric ratio ( ) [ ], the higher the symmetric function score ( ) [ ], which indicates the treetop is more likely to be a false treetop.

(ii) Area Ratio Function ( )
Another area ratio term ( ) is also included to reinforce the assessment of the geometric features of objects in the configuration.

Likewise, an area ratio
[ ] is firstly calculated. The ratio is computed as the proportion of the intersection of object and the underlying segments to the whole area of the segments ( ) by Eq. (8). The larger the area ratio, the higher the degree the geometric feature of object is in accordance with the hypothesis (see Figure 3). Based on the area ratio of the object, the area ratio function is defined as: where and are used to control the position and slope of the sigmoid function respectively. Figure 3: Area ratio calculation of (a) symmetric and (b) asymmetric tree crowns.

(iii) Area/Size Constraint ( )
In conventional marked point processes, a bounded set is specified for marks associated with points, which in our case are the minimum and maximum radii of the disks ( ) as shown in Eq. (4). In our applications, we propose a hard area/size constraint ( ) function to penalize tree crowns falling out of a specific size of [ ] with respect to the forest plots, which is defined as: This constraint is put above the two previous energy functions, and the data term is finally written as:

Prior Term:
This term introduces a priori knowledge concerning the objects layout. In most mature coniferous forest stands, tree crowns won't be overlapped too severely. A repulsive term is then defined as a soft penalize function to penalize severe overlaps in the configuration.

(i) Overlap Function ( )
To define the overlap function, we first introduce a symmetric neighbourhood relation between objects. We say two objects ( ) and ( ) are overlapping if ( ) and we write . Then an overlapping ratio [ ] is calculated as the ratio of the overlapping area between the two objects normalized by the area of the smaller object (see Figure 4): The overlap score ( ) on is then given as: where and are set to control the position and slope of the sigmoid function respectively.
The overlap function of configuration can be expressed as: Figure 4: Overlap ratio calculation of overlapping tree crowns.
In our application, a tree object occurs just once at a location of the pre-extracted local maxima, so there will be no multiple detection problems. The prior term contains only the overlapping function and is written as:

Parameter Setting:
Parameter setting in the model can be distinguished into three categories: physical parameters, weights and thresholds.
The physical parameters and are the size constraints specifying the range of tree crown radius in the forest plots, and are set as 1.0 m and 6.0 m according to our knowledge about the test sites.
Weights are assigned to tune the relative importance we wish to grant to different energy functions or terms in the combination. There are two weights in our model: a factor (see Eq. (5)) used to weight data term and prior term in the calculation of global energy, and the other (see Eq. 11) in data term, to tune importance of symmetric function and area ratio function. Both and are set to 0.5 as we place equal importance on those functions in all our experiments, and they generally gives good results.
Thresholds are set in three energy functions (Eq. (7), (9), (13)). The first group of threshold is set as tolerances on different constraints of geometric characteristic or pairwise interactions. The second group control slopes of the sigmoid functions. In our model, we proposed a Monte Carlo method to estimate those parameters for the sigmoid function.

Model Optimization
Typically, a reversible jump Markov Chain Monte Carlo (RJMCMC) coupled with a simulated annealing is adapted to find the optimal configuration in a classical marked point process. Although the algorithm is proven to be effective to sample the huge configuration space with variable dimensions, the computation cost is usually very expensive, and the process takes a long time to converge. An important benefit we get by constraining the configuration space using features extracted from the CHM image through low-level image processing techniques, is that the resulting configuration space is discrete. Therefore, we limit the search for the optimal configuration within the discrete space of possible subsets (N as the number of local maxima) and adapt the steepest gradient descent algorithm (Baddeley and Van Lieshout, 1993) for effective model optimization in our application.

Experimental Data
The study area is a boreal mature coniferous forest located in the Great  Figure 5(a)-(c)). The CHM images for each plot were produced with a resolution of 0.5m.

Quantitative Results
The proposed model is applied on the three forest plots and detection results are presented and compared with traditional local maxima filtering (LM) in Figure 5 and Table 1.
The images (a)-(c) in the first row of Figure 5 show the results obtained by LM. In those images, red circles with blue crosses in the centre stand for correctly detected tree crowns, while green and cyan circles represent for commission and omission errors respectively. As can be seen, LM is prone to produce a lot of commission errors. The corresponding images (d)-(f) show the detection results of the proposed model. As can be easily interpreted, most green circles have been successfully removed, which means the proposed model could effectively reduce commission errors. It is also noticed in the results that omission errors of the proposed model stay in the similar level with that of LM method. Table 1 depicts detailed quantitative assessment of the detection results of LM and the proposed model. There is an obvious improvement on the results of the proposed model over the LM method it based on. The commission errors dropped dramatically across the three plots, with the largest extent in Plot 1 from 36.2% to 10.3%, while omission errors didn't change much. Averagely, the overall detection quality increased about 15% comparing our proposed model and traditional LM method.

CONCLUSION AND FUTURE WORKS
We have proposed a hybrid framework to detect single trees from ALS data by combining of marked point processes and low-level image processing techniques of local maxima filtering and marker-controlled watershed segmentation. In this model, tree crowns in ALS recovered CHM are considered as a realization of a point process of circles, which enables to take into account both the geometric characteristics and pair-wise interactions of objects in the configuration. Local maxima filtering and marker-controlled watershed segmentation are employed to produce low-level representation of the image, which provides a constraint configuration space for the marked point process to sample for the optimal configuration. Experimental results on real forest plots showed the proposed model has an improvement in detection quality over the traditional local maxima filtering. Future research involves comparing the proposed method with traditional Marked Point Process approach, and incorporating more stochastic scheme in the last stage to recover missed detections in the first step.