Deleting Trees data points

Introduction

In this project we are going to detect the section in the point cloud data set that defines and specifies the trees. We start with the point cloud pruned with the algorithm in the DeleteGround project. Then we apply two different algorithms to detect and delete trees. The first one is based on the return number of each point and rectangular metric. The second one is based on density clustering and intensity. We also use the rectangular metric and density clustering to clean the data set.

Loading packages, Data set and basic preparations

First we load the R packages that we need for this project.

Now we load the data set. Unfortunately, since the size of the data set is too big, we can not render the 3D plot of the data set. Instead, we show a aerial snapshot of the data set.

main=readLAS("main.las")
m=as.data.frame(main)
qmain=kmeans(m[,1:2],30,iter.max = 30)
m$color=qmain$cluster
l=list()
for(i in 1:30){
      l[[i]]=m[m$color==i,]
}
dim(m)

## [1] 1364499       6

aerial snapshot

As you can see in the dimension of \(m\), we have \(1364499\) points in our cloud. Now we are going to dlete the ground points set by using the algorithms we presented in the previous project.

a=list()
for(i in 1:30){
      ww=l[[i]]
      a[[i]]=ww[ww$Z>=delGh(l[[i]],threshold = 2),]
}
for(i in 1:30){
      a[[i]]=delg(a[[i]],ndeletes = 3)
}
res=a[[1]]
for(i in 2:30){
      res=rbind(res,a[[i]])
}
lmain=l
l=a
dim(res)

## [1] 375050      6

As you can see, we have \(375050\) points after deleting ground data points. Let’s take a look at the aerial snapshot of the result.

Result

Now we need to cluster the result by using k-mean clustring on XY-coordinate.

Deleting Trees

Deleting by using return number-rectangular metric

In this method, first we use the return number of each point to identify the location of the trees. Usually the maximum returning numbers are corresponded to trees. After detecting trees locations, we build a rectangular cylinder by using rectangular metric around the each location and delete all the points in the cylinder. Let’s focus on the first two clusters. First let’s have a 3D plot of each cluster. You can click on the plot, rotate it or zoom in.

q1=kmeans(l[[1]][,1:2],10,iter.max = 30)
q2=kmeans(l[[2]][,1:2],10,iter.max = 30)
open3d()

## wgl 
##   1

mfrow3d(1,2)
plot3d(l[[1]]$X,l[[1]]$Y,l[[1]]$Z,aspect = F,col = q1$cluster,xlab = "X",ylab = "Y",zlab = "Z")
next3d()
plot3d(l[[2]]$X,l[[2]]$Y,l[[2]]$Z,aspect = F,col = q2$cluster,xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

Now let’s look at those points at which we have \(returnnumber=4\), in both clusters.

open3d()

## wgl 
##   3

mfrow3d(1,2)
plot3d(l[[1]]$X,l[[1]]$Y,l[[1]]$Z,aspect = F,col = q1$cluster,xlab = "X",ylab = "Y",zlab = "Z",size = 1)
points3d(l[[1]][l[[1]]$ReturnNumber==4,]$X,l[[1]][l[[1]]$ReturnNumber==4,]$Y,l[[1]][l[[1]]$ReturnNumber==4,]$Z,aspect=F,add=TRUE,col="brown",size=6)
next3d()
plot3d(l[[2]]$X,l[[2]]$Y,l[[2]]$Z,aspect = F,col = q1$cluster,xlab = "X",ylab = "Y",zlab = "Z",size = 1)
points3d(l[[2]][l[[2]]$ReturnNumber==4,]$X,l[[2]][l[[2]]$ReturnNumber==4,]$Y,l[[2]][l[[2]]$ReturnNumber==4,]$Z,aspect=F,add=TRUE,col="brown",size=6)

You must enable Javascript to view this page properly.

In these plots, the points with \(returnnumber=4\) are shown by brown rectangulars. As you can see, all the brown points are pointing at trees or bushes. Now we need to build rectangular cylinder by using rectangular metric around the each brown point and delete all points in those cylinders. We specify the raduis of cylinders to be 10 (the raduis is deefined by rectangular metric).

delettrees=function(w,raduis=1,returnnumber=4){
      ww=w
      l=c()
      a=which(w$ReturnNumber==returnnumber)
      for(i in 1:length(a)){
            x=w$X[a[i]]
            y=w$Y[a[i]]
            dx=abs(w$X-x)
            dy=abs(w$Y-y)
            q=which(dx<=raduis & dy<=raduis)
            if(length(q)!=0){
                  l=c(l,q)
            }
      }
      ww=ww[-l,]
      ww
}
w1=delettrees(l[[1]],raduis = 10,returnnumber = 4)
w2=delettrees(l[[2]],raduis = 10,returnnumber = 4)
dim(l[[1]])[1]-dim(w1)[1]

## [1] 720

dim(l[[2]])[1]-dim(w2)[1]

## [1] 2120

As you can see, we are deleting \(720\) points from the first cluster and \(2120\) points from the second cluster by using this algorithm. No wlet’s visualize the results:

q1=kmeans(w1[,3],10,iter.max = 30)
q2=kmeans(w2[,3],10,iter.max = 30)
open3d()

## wgl 
##   5

mfrow3d(1,2)
plot3d(w1$X,w1$Y,w1$Z,aspect = F,col = q1$cluster,xlab = "X",ylab = "Y",zlab = "Z")
next3d()
plot3d(w2$X,w2$Y,w2$Z,aspect = F,col = q2$cluster,xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

Now let’s focus on the points at wich the return number is \(3\).

open3d()

## wgl 
##   7

mfrow3d(1,2)
plot3d(w1$X,w1$Y,w1$Z,aspect = F,col = q1$cluster,xlab = "X",ylab = "Y",zlab = "Z",size = 1)
points3d(w1[w1$ReturnNumber==3,]$X,w1[w1$ReturnNumber==3,]$Y,w1[w1$ReturnNumber==3,]$Z,aspect=F,add=TRUE,col="brown",size=6)
next3d()
plot3d(w2$X,w2$Y,w2$Z,aspect = F,col = q1$cluster,xlab = "X",ylab = "Y",zlab = "Z",size = 1)
points3d(w2[w2$ReturnNumber==3,]$X,w2[w2$ReturnNumber==3,]$Y,w2[w2$ReturnNumber==3,]$Z,aspect=F,add=TRUE,col="brown",size=6)

You must enable Javascript to view this page properly.

In these plots, the points with \(returnnumber=3\) are shown by brown rectangulars. Now we need to build rectangular cylinders delete all points in those cylinders. We specify the raduis of cylinders to be \(8\) (the raduis is deefined by rectangular metric).

w11=delettrees(w1,raduis = 8,returnnumber = 3)
w22=delettrees(w2,raduis = 8,returnnumber = 3)
dim(w1)[1]-dim(w11)[1]

## [1] 2253

dim(w2)[1]-dim(w22)[1]

## [1] 3125

As you can see, we are deleting \(2253\) points from the first cluster and \(3125\) points from the second cluster by using this algorithm. No wlet’s visualize the results:

q1=kmeans(w11[,3],10,iter.max = 30)
q2=kmeans(w22[,3],10,iter.max = 30)
open3d()

## wgl 
##   9

mfrow3d(1,2)
plot3d(w11$X,w11$Y,w11$Z,aspect = F,col = q1$cluster,xlab = "X",ylab = "Y",zlab = "Z")
next3d()
plot3d(w22$X,w22$Y,w22$Z,aspect = F,col = q2$cluster,xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

Cleaning Point Cloud by using Rectangular Metric

As you can see in the previous plots, previous algorithm leaves us with a point cloud containing a lot of dirts. We need to clean that. The first method that we are using is based on rectangular metric. In this method, we pick a point and build a rectangular cylinder around that point. Then we focus on the points contained in that cylinder. If the distance (both the distance and the raduis of rectangular will be calculated by rectangular metric) between the points and the cylinder shell is more than a threshold, we can be sure that those points construct an isolated mass. If we choose the raduis of the rectangular to be small enough, we can be sure that all those isolated masses are not buildings and we can delete them. Let’s focus on the first cluster.

clean1=function(w,maxraduis=1,minraduis=0.5,n=10,maxpoint=0){
      l=dim(w)[1]
      ww=w[,1:3]
      ind=c()
      fi=c()
      while(l>0){
            x=ww[1,1]
            y=ww[1,2]
            dx=abs(w$X-x)
            dy=abs(w$Y-y)
            q=which(dx<=maxraduis & dy<=maxraduis) 
            qq=which(dx<=maxraduis-minraduis & dy<=maxraduis-minraduis) 
            if(length(q)==length(qq)){
                  fi=c(fi,q)
            } 
            ind=unique(c(ind,q))
            ww=w[-ind,]
            l=dim(ww)[1]
      }
      if(length(fi)!=0){
            w=w[-fi,]
      }
      w
}

a11=clean1(w11,maxraduis =6,minrad = .2)
w11$color=rownames(w11)%in%rownames(a11)

open3d()

## wgl 
##  11

mfrow3d(1,2)
plot3d(w11$X,w11$Y,w11$Z,aspect = F,col=as.integer(w11$color)+1,xlab = "X",ylab = "Y",zlab = "Z")
plot3d(a11$X,a11$Y,a11$Z,aspect = F,col="blue",xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

The left side plot is the original data set and the left side plot corespond to the deleted daa set. In the left side plot, the black points are those which are targeted by this algorithm and deleted in the right side plot.Let’s try another raduis.

a111=clean1(a11,maxraduis =15,minrad = .5)
a11$color=rownames(a11)%in%rownames(a111)

open3d()

## wgl 
##  13

mfrow3d(1,2)
plot3d(a11$X,a11$Y,a11$Z,aspect = F,col=as.integer(a11$color)+1,xlab = "X",ylab = "Y",zlab = "Z")
plot3d(a111$X,a111$Y,a111$Z,aspect = F,col="blue",xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

We can continue this process by using different raduis. Now let’s implement the algorithm on the second cluster.

a22=clean1(w22,maxraduis =6,minrad = .2)
w22$color=rownames(w22)%in%rownames(a22)

open3d()

## wgl 
##  15

mfrow3d(1,2)
plot3d(w22$X,w22$Y,w22$Z,aspect = F,col=as.integer(w22$color)+1,xlab = "X",ylab = "Y",zlab = "Z")
plot3d(a22$X,a22$Y,a22$Z,aspect = F,col="blue",xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

Now let’s use different raduis.

a222=clean1(a22,maxraduis =5,minrad = .2)
a22$color=rownames(a22)%in%rownames(a222)

open3d()

## wgl 
##  17

mfrow3d(1,2)
plot3d(a22$X,a22$Y,a22$Z,aspect = F,col=as.integer(a22$color)+1,xlab = "X",ylab = "Y",zlab = "Z")
plot3d(a222$X,a222$Y,a222$Z,aspect = F,col="blue",xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

We can continue this process by using different raduis. Now let’s implement the algorithm on the second cluster.

Cleaning by Density Dlustering

Now that we have a point cloud data set in which the masses are far enough from each other that would enable us to use density clustering. Density clustering tries to put each separate mass into a cluster. The masses must be far enough so that the algorithm works properly. To apply density clustering, we need to specify the minimum distance between masses. To do that, we use k-tree neighborhod algorithm and plot the denisty function of freqiencies.

library(dbscan)
par(mfrow=c(1,2))
kNNdistplot(a111[,1:2])
abline(h=1.1)
kNNdistplot(a222[,1:2])
abline(h=.9)

As you can see, for the first cluster, the best raduis would be \(1.5\) and for the second one the best raduis would be \(1.4\). Now we apply density clustering by using these raduis.

qa111=dbscan(a111[,1:2],eps = 1.5)
a111$color=qa111$cluster
qa222=dbscan(a222[,1:2],eps = 1.4)
a222$color=qa222$cluster

open3d()

## wgl 
##  19

mfrow3d(1,2)
plot3d(a111$X,a111$Y,a111$Z,aspect = F,col=a111$color+1,xlab = "X",ylab = "Y",zlab = "Z")
plot3d(a222$X,a222$Y,a222$Z,aspect = F,col=a222$color+1,xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

Now let’s take a look at the cluster table for the first cluster (a111):

d111=as.data.frame(table(qa111$cluster))
d111=d111[order(d111$Freq,decreasing = T),]
d111

##    Var1 Freq
## 3     2 2471
## 6     5  975
## 2     1   77
## 5     4   75
## 4     3   68
## 8     7   64
## 10    9   39
## 1     0   29
## 7     6   11
## 9     8    5

And here is the cluster table for the first cluster (a222):

d222=as.data.frame(table(qa222$cluster))
d222=d222[order(d222$Freq,decreasing = T),]
d222

##    Var1 Freq
## 8     7  495
## 11   10  478
## 20   19  466
## 6     5  414
## 5     4  400
## 38   37  362
## 7     6  330
## 17   16  322
## 37   36  322
## 33   32  294
## 21   20  157
## 4     3  153
## 29   28  137
## 25   24  119
## 13   12   98
## 22   21   97
## 39   38   84
## 16   15   80
## 23   22   80
## 26   25   80
## 18   17   70
## 40   39   65
## 14   13   61
## 3     2   49
## 1     0   47
## 2     1   36
## 27   26   36
## 19   18   30
## 9     8   22
## 30   29   18
## 15   14   13
## 24   23   11
## 35   34   11
## 10    9   10
## 34   33    9
## 31   30    8
## 32   31    8
## 12   11    6
## 36   35    6
## 28   27    5

As you can see in both cases, there some clusters containing less than \(200\) points. We are sure that any cluster containing buildings must have more than \(200\) points. So we can delete those clusters with less than \(200\) points.

c1=d111[d111$Freq>=200,1]
b1=a111[a111$color%in%c1,]
c2=d222[d222$Freq>=200,1]
b2=a222[a222$color%in%c2,]

open3d()

## wgl 
##  21

mfrow3d(1,2)
plot3d(b1$X,b1$Y,b1$Z,aspect = F,col="red",xlab = "X",ylab = "Y",zlab = "Z")
next3d()
plot3d(b2$X,b2$Y,b2$Z,aspect = F,col="blue",xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

Deleteing Trees by Density Clustering and Z-Distribution

Now that the data sets are clean enough, we can apply density clustring to have a clean set of clusters.

library(dbscan)
qb2=dbscan(b2[,1:2],eps = 1.5)
b2$color=qb2$cluster
open3d()

## wgl 
##  23

plot3d(b2$X,b2$Y,b2$Z,aspect = F,col=b2$color,xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

As you can see, we have \(10\) different clusters. Some of them consists of buildings, some of consists of trees and some of consists of both. To detect the clusters consisting of trees, we need to focus on the distribution of Z-coordinate of clusters. Theoretically, if a cluster consists of trees, the distribution must be left skewed since we already delete the ground and so the frequency of the points is more around the tree crown. Let’s look at denity plot of Z-coordinate for each cluster.

library(e1071)
par(mfrow=c(2,5))
for(i in 1:10){
      bb=b2[b2$color==i,]
      plot(density(bb$Z),main=skewness(bb$Z))
}

The number on top of each plot is the skewness of Z-coordinate correspond to that cluster. The more positive is the skewness, the more the plot is left skewed. Since the skewness of \(7\)th and \(9\)th plot are \(0.81\) and \(0.59\), according to our theory, the \(7\)th and \(9\)th clusters must be trees. Let’s check that

open3d()

## wgl 
##  25

mfrow3d(1,2)
plot3d(b2[b2$color==7,]$X,b2[b2$color==7,]$Y,b2[b2$color==7,]$Z,aspect = F,col="red",xlab = "X",ylab = "Y",zlab = "Z")
next3d()
plot3d(b2[b2$color==9,]$X,b2[b2$color==9,]$Y,b2[b2$color==9,]$Z,aspect = F,col="red",xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

The left side plot is the \(7\)th cluster and the right side plot is \(9\)th cluster. As you can see, our theory is true. So we can detect those clusters consisting of trees and then delete them. Let’s go through the same process for the first cluster.

qb2=dbscan(b1[,1:2],eps = 1.5)
b1$color=qb2$cluster
open3d()

## wgl 
##  27

plot3d(b1$X,b1$Y,b1$Z,aspect = F,col=b1$color,xlab = "X",ylab = "Y",zlab = "Z")

You must enable Javascript to view this page properly.

As you can see we only have two clusters. Let’s look at the skewness of each cluster.

par(mfrow=c(1,2))
for(i in 1:2){
      bb=b1[b1$color==i,]
      plot(density(bb$Z),main=skewness(bb$Z))
}

As you can see, the skewness of both clusters is so small. That means none of the clusters are left skewed.