Classification of LiDAR data

Professional software available in the domain of classification of LiDAR data as well as the literature available for feature extraction cum classification reveal that a certain hierarchy has to be followed in the process. Before proceeding to extract the features present on the terrain through the process of classification, we should recall that in a typical terrain there are ground, buildings, trees, roads and vehicles. The errors present in capturing the data often results in outliers also being present in the dataset.

One typically begins with the elimination of outliers. Outliers are those points which lay isolated within a sphere of radius $\varepsilon$. If the sphere contains less than $n$ points, the points are labelled as outliers. Both $\varepsilon$ and $n$ are thresholds that are provided by the user of the software or the classification algorithm. To find the outliers, each of the points in the dataset is considered as the center of a sphere of radius $\varepsilon$, and then the number of points within the sphere are counted. If the number is less than $n$, the said point is labelled as an outlier.

The popular software named TerraSolid (which works on the CAD platform Microstation), also has a concept of finding low lying points in addition to the first step of finding outliers. This routine is present in the utility named TerraScan. To find a low lying point it assumes a cylindrical neighbourhood of a point, where the axis of the cylinder is the z-direction. If for a given threshold $h$, which forms the height of the cylinder, there is no other point contained in the cylinder, then the point is labelled as low-lying.

The next step after the classification of the outliers and low-lying points, is to classify the ground points. Since 1999, several algorithms have been researched on for labelling the ground points. Sithole and Vosselman (2004) have provided an interesting review of the algorithms for ground classification. In their paper, the authors review and test the performance of several algorithms on the ISPRS test dataset. It has been reported later in the literature that these algorithms were not suitable for all the terrains, and therefore some interesting additional algorithms have been developed. The issues with ground point classification have been reviewed and addressed in a paper by Meng, Currit and Zhao (2010). The process of classification of the ground points has been referred in the literature as “filtering”. In TerraSolid, in addition to the slope of the ground, information regarding the longest edge of a building is also sought. This information is required in order to avoid classifying a pretty long building as a ground cluster.

After the ground points have been “filtered out”, there are trees and buildings to be detected. Some of the urban areas do not even contain trees, but some do! Sometimes, the trees are too close to the buildings. If the intensity information is not used, the tree points and the building points appear to be in the same cluster.

There could be multiple strategies for building extraction from the unclassified datasets. TerraSolid people first classify the low-vegetation and high-vegetation points just by their height from the ground. The remaining points are then classified as buildings using their own set of algorithms. Although this sounds pretty crude, it does help. The buildings can be then reconstred into a CAD model as TerraSolid sits on a Microstation enviroment. We shall deal with building extraction from LiDAR data in a separate post.

Road points could be classified from the ground points themselves. Researchers have reported the use of intensity values from LiDAR data to separate the road and other points. However, the problem becomes different when bridges and flyovers have to detected. Apart from the intensity values, the height values also need to be used.

Trees could be detected using a template matching procedure. A botanist usually knows the shape of a tree. The property that LiDAR data can capture multiple storeys from the trees, comes in handy here. Tree templates are available as RPC (Rational Polynomial Coefficient) models for purchase. However, this database is pretty limited. There is a research opportunity to create these RPC models by scanning different forests. The Indian biodiversity is pretty high, and an initiative to create tree models for the different species of trees available in India (at least) will be an excellent direction for research and development.

1. Sithole and Vosselman (2004), Experimental comparion of filter algorithms for bare-earth extraction from airborne laser scanning point clouds. doi:10.1016/j.isprsjprs.2004.05.004
2. Meng, Currit and Zhao (2010), Ground Filtering Algorithms for Airborne LiDAR Data: A Review of Critical Issues , doi:10.3390/rs2030833 

Delaunay Triangulations - Algorithms and APIs

The triangulation of a surface is a topological process. The concept of triangulation is well explained in Algebraic topology and followed up with enormous amount of research in Computational Geometry.

Delaunay Triangulation is a different ball game because in the process of triangulation, a certain rule has to be followed - Any point must not lie in the circumcircle formed by three other points. Put in pretty simple words, wow, what a rule! I will attempt to simplify the explanation of the process in the following paragraphs.

Let us suppose that we are given a set of N two dimensional points. For the sake of simplicity let us assume that these points are unique (no two points have the same coordinates). We take out three points from this set of N points, and then check whether these points are collinear or not. If they are, we will have to choose another subset of three points. So let us suppose that they are not collinear. Create a triangle with these three points. Take out another point from the main set.  Check if this point falls within the circumcircle of the earlier triangle. If yes, "cancel" the earlier triangle, and create triangles such that the given rule is followed. If no, create another triangle with the nearest edge. This process goes on and on until the entire list of points is exhausted.

We have to understand that this is a computationally intensive process. Now as a person, working with LiDAR data and interested in processing it to extract features like buildings, trees etc. you might not be really interested in designing and coding the algorithm for creating the Delaunay Triangulation of a given set of points. Surely there are people in the domain of computer science who have already done the hard work for coding these algorithms. Two of the most popular resources are QHULL and CGAL.

QHULL is a divide and conquer algorithm (Barber, C.B., Dobkin, D.P., and Huhdanpaa, H.T., "The Quickhull algorithm for convex hulls," ACM Trans. on Mathematical Software, 22(4):469-483, Dec 1996), and is available as an open source project on the Internet (http://www.qhull.org). The binaries are available too! One just needs to spawn the binary executable with the right options to generate the delaunay triangulation. Till very recently, MATLAB was using the QHULL library but has switched over to CGAL (for better or worse). Apart from usual computations of Delaunay Triangles, the QHULL code is able to compute the convex hulls and the voronoi diagrams. A comprehensive help is given on their website.

The Computational Geometry Algorithms Library (CGAL) is a large collection of algorithms in the field and a small section of this is also the computation of convex hulls, alpha shapes, and Delaunay Triangulations. (We can talk about alpha shapes later on.) CGAL has an API with C++ to enable programmers to write codes for computation of triangulations and delaunay triangulations. It is important to note that in the jargon of CGAL, "triangulation" and "delaunay triangulation" are two different processes.

As on April 20, 2011, I received a mail from the CGAL project people [subscribe to the cgal-announce newsletter] that they have released CGAL 3.8, with new features. It supports 3D mesh generation, 2D and 3D triangulations (faster and efficient algorithms), and better computations of 3D alpha shapes. 

RANSAC Algorithm - Some Outputs

It is time to give you a glimpse of how the algorithms mentioned in the previous post perform. I would like to acknowledge Optech for providing this data, and I am taking liberties to show some output on a small "cutout" of the same.

From an aerial photograph, the object containing multiple planes is given below:

Screenshot of the object with multiple planes

Now, as we know, the real picture is never seen from an aerial photograph, so we also give you a screenshot of the data corresponding to this cutout.

LiDAR data for the corresponding object

The RANSAC algorithm basically requires an input regarding the threshold for finding the inliers to the planes. But are all detected planes or best_models acceptable? One who is working with LiDAR data would often say no. Therefore, we add another restriction to the RANSAC algorithm. We accept a result only if it contains more than a certain number of points.

We therefore, make some assumptions regarding the thresholds as follows:

1. The minimum number of points in a plane should be 40. 
2. The threshold for finding inliers would be say 0.20m. 

I am not mentioning here how we have selected these numbers. Maybe we could discuss them later. But the moot point of the discussion is that we should get similar results if we use the same thresholds.

Peter Kovesi's algorithm will always find out one best plane from the data set. So what I did is to follow an approach mentioned in Chapter 14 of the book titled Topographic Laser Scanning and Ranging which is authored by F. Bretar. To explain in brief, we find out one plane, remove the inliers from the dataset, find out a second plane, remove the inliers again and so on. We do this until we are not able to find a plane, or the number of remaining points is less than or equal to 3.

Output using Peter Kovesi's algorithm is shown in the following figure:

Detected planes (denoted by different colors) using the Peter Kovesi's algorithm

There were 9 planes detected, and  the program took about 10-15 seconds to execute. This code has been written in MATLAB and 10-15 seconds is a fairly good time for a code to execute. However, the results generated are not satisfactory, as it seems that the RANSAC algorithm preferably selects horizontal planes than the tilting ones. There are two mildly sloping planes in this dataset which do not touch but have parallel center lines. However, their slopes are not in the same direction. The algorithm detects the midriffs of both the planes as belonging to one!

Our next goal is to test the output using the Mobile Robot Programming Toolkit (MRPT), I have already provided you with the link of the demo code. MRPT is easily installable on the Ubuntu platforms using Synaptic. Just find mrpt and select all the installable software from the displayed list. I just customised the demo code to read the data (shown above) and execute the program.

The output given is as follows:

OpenGL screenshot (basic API of MRPT) to display the detected planes and point cloud

The program took 500 milliseconds to execute and it detected 12 planes (Remember! The thresholds are the same). However, the base display of planes is not proper and a study of the TPlane structure in the code is required.

Let us now go to Tim Zaman's code in MATLAB given on his website. What is appreciable is that Tim Zaman with all his humility does accept that his code is slow, but argues that it needs to be structured and readable. I agree with him. However, let us see how his code has performed with our data. Unfortunately Tim Zaman's code needs to be tuned further to be used for detecting multiple planes, and I am therefore providing you the elementary output (of detecting a single plane). Remember, the thresholds are always the same.

Output found with Tim Zaman's code (www.timzaman.nl)

Tim Zaman's code in MATLAB took more than 10 minutes to find out something like the above figure. I have been interacting with Tim to figure out the problems with the code.  

RANSAC Algorithm - Tools (3D plane extraction)

As I said in my previous post, that there were variants for the RANSAC algorithm and one might be interested in their codes. I would be catering to the MATLAB and C++ enthusiasts only as of now. I could post my own Java Code sometime at a later stage, but it basically is an extension of the MATLAB code proposed by Peter Kovesi.

1. Category MATLAB: Peter Kovesi (See here)
2. Category MATLAB: Tim Zaman (See here)
3. Category C++: Mobile Robot Programming Toolkit (See here)

In my next post I will talk about some of the problems encountered by me in extracting planes from the RANSAC algorithms. 

RANSAC Algorithm - Introduction

RANSAC stands for RANDOM SAMPLE CONSENSUS. This algorithm was established by Fischler and Bolles in 1981. Since then, there have been multiple variants proposed for using the RANSAC algorithm. Some of these variants were used by Nardinocci et al(2001) and Forlani et al (2003) for extraction of planar features or to put in simply to extract buildings using LiDAR data. 

A brief description of the RANSAC algorithm could be found at WikiPedia. You might sneer at me for giving a reference to Wikipedia, but believe me, it makes the understanding of the algorithm simpler. Can I make it more simpler? Yes, I could. Just keep on reading.

Let us assume that we are just talking about extracting planes. Further, let us assume that we are given a set of N points, where N > 3 . Let us also assume that our data contains at least one plane. Now, for a point to belong to a plane, it should be within a certain distance from it. We have to ask the user of the algorithm regarding that distance. Let this distance "threshold" be t.

It is a well known fact that to define a plane, we need just 3 points. So, what do we do? Pick 3 points randomly from the set. This will define the plane. But, is this the right plane? We will have to find the distances of the remaining points from this plane. The points that would "belong" to the plane (points that are within a distance of t from the plane) are called "inliers". Record the three points and its number of inliers. Let us call this record as best_model.

Let us do another iteration of this. Pick another set of 3 points, defining a plane. Find the inliers for the plane. If the number of inliers are greater than those in the best_model, delete the record maintained earlier and save this one.

Now, we shall keep on doing this for a certain number iterations k.

You might argue that if N is the given number of points, and every time we have to extract 3 points from the set, then the right number of iterations would be C(N,3)! Turns out to be a fairly huge number if one were to do this for all the points! Fear not, there is an algorithm coded by Peter Kovesi where these problems are solved! In fact, the Wikipedia page also gives a hint of the process.

Now that we have solved the problem of the number of iterations, we must now concentrate on how do determine the exact threshold t for fitting the plane. Fischler and Bolles (1981) suggest that this can be achieved either by experiment or by analytical methods. There is a book titled Topographic Laser Scanning and Ranging in which Chapter Number 14 deals with determination of the threshold for a set a set of points. This chapter has been written by Frederic Bretar. Bretar uses an icosahedron to "collect" the normals of the points to find out the right threshold. I am however, not convinced by the method but he has claimed in his paper Bretar and Roux (2005) that this works.