Project5: Working with LiDAR data

• lidar-data
• On dover you can find LIDAR data at: /mnt/research/gis/DATA/LIDAR/

Overview

The goal of this project is to work with Lidar data and take a shot at classifying it, that is, identifying the ground, buildings and vegetation.

Since understanding Lidar would be pretty hard, if not impossible, without visualizing, I provided a skeleton that is able to read a lidar file and visualize it based on return number and code; you will find this skeleton in your GitHub repositories once you accept the project. Your code will build on this skeleton so you probably want to start by reading it and having an idea how it works. You can ignore all the OpenGL stuff, and focus on how the data is filtered. LiDAR data comes classified or partially classified from the data provider.

The interface: Your code will take as input the name of a lidar dataset, in txt form. For example,

[ltoma@dover:\~] ./lidarview  zurich.txt 

[Note: do we want to save the coded points? This may be necessary to compare outputs. It's easy to add it once it's necessary, so for now let's say its not needed]

Converting lidar .las/.laz files to .txt files with las2txt

For simplicity, we'll convert .las/.laz data to .txt using las2txt tool from LAStools. LAStools is a suite of tools for working with Lidar data by Martin Isenburg. You can download LAStools from its GitHub site; or, if you prefer, you can download the LAStools library from Martin's rapidlasso site.

Once you download LAStools:

cd LAStools
make 
cd bin
ls

bin/las2txt -h

bin/las2txt -i data/zurich.las -o zurich.txt -parse xyznrc

Note that LAStools comes with a folder with data (you can also find it here on Github), and you'll recognize zurich.las in the list.

Or, you could read .las/.laz files directly using LASlib

cd LAStools
cd LASlib
ls
cd example

In fact, glancing at the example LAStools/LASlib/example/lasexample_simple_classification.cpp , you will see that it has code for how to create a min_z and a max_z grid for the points, which is very similar to what you all did in class; the difference is that the .las data has a header, which stores minx and maxx, miny and maxy values. That's convenient! Perhaps using text lidar data wasn's such a simplification afterall. [Note: next time work with .las/.laz directly]

Outline

1. Create a grid corresponding to the first returns: The first task is to create a grid corresponding to the first return points (and render it, see below).

   If multiple points fall in one cell, the cell should store the maximum value.

   The resolution of the grid: Ideally it needs to match closely the resolution of the lidar data. You don't want too many points per cell, not too many empty cells. Since lidar data is usually approx. 1m resolution, a good starting point is step=1 or step=2. The step shoud be either specified by the user on the command line, or pre-defined at the top of the code. Your function should print some info on how many points fall in one cell, and how many empty cells, etc, so that the user can infer whether step should be increased or decreased.

   It seems that this grid should be created when the data is loaded, and it should exist as a global along with the array of lidar points.

2. Render a grid with the color based on height and using a triangle strip: The current code renders the lidar points by iterating through the array and for each point:

   glBegin(GL_POINT); 
   //convert (x,y,z) to [-1,1]x[-1,1]x[-1,1]
   glVertex3f(...)
   glEnd();
   

   You'll need to add code to render a grid. We'll assume the grid is diagonalized in the obvious way and you'll iterate through the points, and render each triangle:

   glBegin(GL_TRIANGLE);
   
   //convert (x,y,z) to [-1,1]x[-1,1]x[-1,1]
   //set color of this vertex 
   glVertex3f(...)
   
   //convert (x,y,z) to [-1,1]x[-1,1]x[-1,1]
   //set color of this vertex 
   glVertex3f(...)
   
   //convert (x,y,z) to [-1,1]x[-1,1]x[-1,1]
   //set color of this vertex 
   glVertex3f(...)
   
   glEnd();
   

   The color: Assign the color by height. You could divide the z-range in a couple of buckets (3 or so), [z1,z2,z3,z4] and assign a color to each bucket. When you want to render a point, you choose a color based on which bucket its height falls in.

3. Render a grid with hill shading: Same as above, but add in hill shading. Hill shading is particularly relevant for archaeology because its shows the relief. [Note: this is tentative until someone gets it to work]

4. Label the ground points: The three main types of points that need to be distinguished in Lidar data are ground, vegetation and building. Your task is to label the ground points. To this end, the structure lidarpoint contains a field mycode which is initialized to 0. After you identify the ground, all ground points shoudl have mycode=2 and the non-ground points should stay at mycode=0 (or be labeled as vegetation or building or noise). For e.g. if you decide that a point is noise, you could set it as mycode=7.

   The ASPRS standard codes are below:

   

   Note that structure lidarpoint also includes a field called code: this stores the code read from the las file, and it is the classification code given by the data provider. Sometimes lidar data comes unclassified (all codes are 0), sometimes it's partially or fully classified by the data provider; the qualiy of the classification can vary and the data may need additional rounds of filtering with refined methods. Anyways, if the codes are available, you can use them for comparison with your code.

   Note: Pressing 'c' cycles through the color maps, one of which is by mycode. This makes it easy to test.

   Idea: To find the ground points you will want to start from the last returns: these are the points with either a single return (return_number = 1 and nb_returns=1), or points part of a multiple return, and have their return_number equal to the nb_returns (return_number = nb_returns). Denote the last return points by P'.

   Take a look at a classified dataset (e.g. zurich) and see how the last returns look, and how they are classified. Not all last return points are ground: some are ground, some are buildings, some are vegetation and some are noise. In other words, you can;t just take all points in P' and label them as ground. You will need to implement some smarter filtering that filters out the non-ground points, i.e. the roof tops and other noise (note that there isn't much vegetation in the last return).

   Some possible ideas:

   • a point of min elevation must be ground, so you can start the search (e.g.bfs) there
   • use slope: a steep increase is unlikely to be natural earth, so it's non-ground
   • use height difference: if the height difference between two points is too high, its unlikely that both are ground
   • height threshold: (most) buildings are higher than x meters (note: lidar measuremenst are in meters)
   • erosion: traverse the set of points and set each point to the minimum z-value in its neighborhood. Here we think of the non-ground points as features, because their z-value is larger than the z-value of the ground and they "stand out". The effect of erosion is that features larger that the window will disappear and in general features will get eroded because their perimeter will be set to ground. Several round of erosion will get rid of buildings, but also of ground that may be high (like natural hills).
   • dilation: traverse the set of points and set each point to the maximum z-value in its neighborhood.
   • Erosion followed by dilation is commonly used in computer vision to eliminate noise and strengthen/emphasize the features.

   Create a grid of last returns: No matter which ideas you use, you will need to be able to find the neighbors of a point. I suggest that you store the set P' (last returns) in a grid. If more than one point falls in the same cell, you could set the value of the cell to the lowest point; or you could have each cell in the grid store, instead of a value, an array of the points inside it.

5. Create a grid corresponding to the ground points: Let G be the set of points labeled as ground above. Create a grid corresponding to G. The holes that correspond to point that are not-ground (like under a roof ) are assigned a ground height by interpolating from ground neighbors (e.g. nearest neighbor). Add a key to keypress so when the user presses, this ground grid is rendered.

6. Paper:Write a paper/report that summarizes what you did:
   • Your name(s) and github repo for this project
   • Creating a grid for the first returns: what spacing is default, and how it can be changed.
   • Render a grid:what key to press, etc.
   • Label ground points: Summarize your approach to label the ground, and how it works. Run it on zurich.txt and include pictures. How close can you get to the codes that come with zurich.txt? Running it on more datasets is a plus, and may actually be easier to detect ground when there are no buildings. The approach does not to be perfect (it can't be), so focus on what you learned.
   • Create a grid corresponding to the ground points: Sumarize your approach, and what one needs to do to see it (what keys to press, etc).
   Keep the paper concise, think of it as a report: its goal is to summarize what you have done.

And finally, this is an open ended project, so feel free to explore and take it to any direction you want to. If you have any suggestions, let me know!

What to turn in

Grading

• Grid of first returns
• Render grid
• Hill shading
• Label ground
• Ground grid
• Paper

Enjoy!