Mathematical morphology

Mathematical morphology is a methodology for transforming geometrical struc-tures to other geometrical strucstruc-tures with certain central operators. Very general theory can be formed; see e.g. [42, pp. 13–35]. For our purposes, we can concentrate on morphology of subsets ofR²; see e.g. [12, pp. 6–8]. Usually the morphological transformations are combined with measurements of the transformed set. For ex-ample, a set is transformed several times, and every time, the area, boundary length, or the number of connected components of the transformed set is calculated. This

dilation erosion

Figure 4.2: Examples of the morphological transformations erosion and dilation performed on a planar set with the black disc.

procedure could be compared to extending a linear regression model with trans-formed covariates to gain additional information on the process being modelled.

As an example, Arns et al. [43] used measures of transformed sets to characterize different types of disordered materials.

The morphological operators are built upon structuring elements K ⊂ ^R2. Dif-ferent shapes of K lead to different types of transformations and hence different interpretations of measurements made from the transformed sets. To simplify the situation we only consider the origin centred disc B(o,r) = {x : ||x|| ≤ r} ^{as a} structuring element. Theerosionof setA⊂^R2with the discB(o,r)is the set

A⊖B(o,r) ={x∈ A:B(x,r)⊂A}^{. (4.12)} In other words, the erosion of a set A with B(o,r) is the set of those points in A where the centre point of a disc of radius rcan be located so that the disc is fully contained inA. Thedilationof setAwith the discB(o,r)is the set

A⊕B(o,r) ={x∈^R2:B(x,r)∩A6=^∅}^{. (4.13)} In other words, the dilation ofAwithB(o,r)contains all the points where the centre point of a disc of radius r can be located in such a way that the disc comes into contact withA. Erosion and dilation of a planar set with a disc have been presented in Figure 4.2.

5 HORVITZ–THOMPSON-LIKE ESTIMATORS

Horvitz and Thompson [9] presented an estimator for population totals when sam-pling without replacement from a finite population. The estimator is a sum of the sampled attributes of interest, over which the estimation of a total is wanted, weighted with the probabilities of sampling the subjects that were sampled. For the Horvitz–Thompson estimator, the sampling probabilities are known beforehand.

However, when ITD is used for tree detection, and in a sense for sampling the under-lying tree population, the probability of detecting a certain tree, the detectability, is not known. The probabilities must be estimated somehow. To make the difference clear, estimators that use estimated detectabilities are called Horvitz–Thompson-like e.g. in distance sampling of animal population sizes [11]. These estimators were also calledHorvitz–Thompson typeinI. There are different approaches to calculating the detection probabilities, but our approach is based on interpreting forest as a re-alization of a marked point process, the assumption that the detection of a tree is affected by other trees, and certain assumptions relating to mathematical morphol-ogy. Other approaches have been studied e.g. by [44] and [45].

LetX = {(xi,mi,si)}_i=1^N be a realization of a marked point process onW ⊂ ^R2 with attributes of interestmi and secondary markssi that affect the detectability of pointi. We assume that the point locations and marks do not contain measurement errors. We are interested in estimating a population totalτ=∑_i=1^N mi. The Horvitz–

Thompson-like estimator of this total can be written as

b τ=

∑

miDi

p(s_i)^{, (5.1)}

where Di is an indicator variable of detection – 1 if pointi has been detected and zero otherwise – and p is the detectability that depends on s_i. During estimation, the sum is calculated over the detected points, and hence the marks or detectabilities of undetected points are not needed.

To calculate the detectabilities, it is assumed that the detection of points has a sequentialnature. The sequence emerges from the simple idea that in the application cases considered in this thesis, a tree that is closer to the remote sensing device has an ability to affect the detection of a tree that is further away – for terrestrial remote sensing the closeness relates to the actual distance of the tree from the device, and for airborne remote sensing it relates to the height of the tree. More generally, this means that point i affects the detection of point j only if i < j. Nothing affects the detection of the first point, and hence it always has detectability 1. For point j>1 the pointsi =1, . . . ,j−1 produce areasAi⊂^R2with properties that depend on the location xi and mark si that affect the detectability of point j. To be more precise, a point with a mark sj would not be detected if xj ∈ G(^S_i=1^j−1Ai), where Gis a morphological transformation that depends on sj. If we assume that xj is a single random point that is uniformly distributed inW the probability of not being

detected can be written as

and hence the detectability is

p(sj) =1−P

It should be noted that the detectabilitiesp(s_j)that we have formulated are con-ditional probabilities, conditioned on the points i < j, in one sense the observed

”history” of the process. Assuming that the observationsi< jare fixed and x_j is a single uniformly distributed random point, the detectabilityp(s_j)is not an estimate, but a known fixed probability.

It is also possible to condition the detectability, if it is known that xj must be-long to a certain subset of W, by changingW to that subset in the aforementioned formulas. The subset can be of lower dimensionality thanW, for example a line or a circle – a situation encountered inIII. In this case|^.|measures length. It is useful to think Xas a point process in a larger window Sthat containsW, so that all the points xi that lie outsideW but have such an area Ai that comes into contact with W can be taken into account during detectability calculations, as these points can affect the detectability of points inW. Depending on the application and available data the detectability can depend on all points i = 1, . . . ,j−1, if it is possible to derive Ai also for undetected points from the data, or only those points that have been detected. In the latter case Ai = ∅ for the undetected trees in Equations 5.2 and 5.3.

Horvitz and Thompson [9] showed that the Horvitz–Thompson estimator is un-biased, i.e. E[bτ] =τ. Similarly, Horvitz–Thompson-like estimator is unbiased if the assumptions related to the detection process and underlying marked point process hold; an example proof is presented inIIIfor data generated by the Poisson process.

For the Horvitz–Thompson estimator, formula for the variance is known [10, p. 70]:

var(_bτ) =

This theoretical formula depends on all of the members of the population N and cannot actually be used for variance calculation. However, an unbiased estimator of the variance, calculated with the detected members of population, is also known [10, p. 70]:

This estimator of variance can be used with the Horvitz–Thompson-like estimator and should hold when the aforementioned model assumptions hold.

An important special case of the Horvitz–Thompson-like estimator is the esti-mator of the population size, i.e. the number of pointsN, which is given by settingb mi=1 for alli:

Nb =

∑

Di

p(s_i) =

∑

idetected

1

p(s_i)^{. (5.6)}

For the detected pointsithe estimatorNb can be divided into componentsNbi= _p(s¹_i), so that Nb = ∑Nbi. Nbi can be interpreted as an estimator of the number of points that are similar to point i. This leads to an estimate of mark distribution, and an estimate of the cumulative distribution function can be given for s– assuming the mark is one-dimensional – as

Fb(s) = ¹ Nb

∑

s_i≤s

Nb_i, (5.7)

where the indexing is again over the detected points. This formulation can also be extended to give distribution estimates over other observed marks of interest by replacings_iwith the mark in question in Equation 5.7. TheNb_ialso make it possible to estimate the mean of a mark as a weighted mean

1 Nb

∑

idetected

Nbimi. (5.8)

5.1 DERIVING DETECTABILITIES FROM AIRBORNE REMOTE SENS-ING DATA

InI, the problem of estimating stem density – that is, the total number of trees in a forest area – from ALS based ITD was considered. For stem density estimation mi =1. Although the algorithm of Lähivaara et al. [23] which was used for ITD in the paper produces 3D tree objects, the situation was simplified by considering the problem only in two dimensions. This led to si being the maximum crown radius of treei, andAi= B(xi,si). It was assumed that only the trees with larger value of s_icould affect the detection of treejwiths_j <s_i. This led to the sequence of points being ordered according tos_i, so thats1>s2>. . .>sN. It was assumed that a tree would not be detected if a certain proportion of its crown was covered by the larger trees. This proportion was controlled by introducing a tuning parameter α ≥ ^{0 to} the morphological transformation, and the detectability became

p(sj) =1− h_Sj−1

i=1A_i

⊖B(o,αs_j)^iTW

|W| ^{. (5.9)}

The interpretation of parameter α is as follows. If α = 1, it is assumed that a tree will not be detected if the crown of the tree is fully covered by the larger trees. If α = 0 – meaning no erosion is performed – it is assumed that a tree will not be detected if the centre point of the crown is covered by the larger trees. A value of αbetween 0 and 1 means that at least α per cent of the crown radius has to be covered for detection to fail. The limiting case here is the centre point being covered. Naturally, this construction could be extended to allow also trees with crown centre points outside of the set of larger tree crowns to be left undetected by

Figure 5.1: An example of the morphological aspects of detectability calculation.

On the left, an example forest area with detected trees. Detectability is calculated for the tree marked as a black disc. On the center, the data are limited to contain only trees that are larger than the black one. On the right, an appropriate morphological transformation is performed, in this case erosion withα=1, producing a set where the tree could have been located in such a way that it would have been completely covered by the larger tree crowns.

replacing the erosion with a dilation. However, this was not considered inIorII. In this case the effect of undetected trees to the calculation of the detectability was not considered, andAi =∅for undetected trees. The derivation of the detectability has been illustrated in Figure 5.1.

5.2 DERIVING DETECTABILITIES FROM TERRESTRIAL REMOTE