Project5: Working with LiDAR data

Lidar data


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 
This will read the lidar point cloud zurich.txt, classify it and render it.

[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

The standard format for Lidar data is the las format (.las) and the compressed las format (.laz). If you download lidar data from the web it it will come as a .las or .laz file. For complete specification of the las format check out ASPRS las format.

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
cd bin
Some of the modules in bin are licensed, some are free. Check out the doc page for las2txt.
bin/las2txt -h
I ran it as
bin/las2txt -i data/zurich.las -o zurich.txt -parse xyznrc
This tells las2txt to write for each point: it's x,y,z, number of returns, return number, and code, in this order. When you convert lidar data to txt it's important to use the same -parse xyznrc command, because lidarview expects to read the values in this order.

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

To avoid converting to .txt and work with files in .las/.laz format directly you could use LAStools/LASlib, which is a free and open source C++ API for reading and writing LIDAR data. If you downloaded LAStools you already have it (or, you can find it on GitHub here).
cd LAStools
cd LASlib
cd example
The example folder has examples of how to read and write a .las file and includes a Makefile.

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]


  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:
    //convert (x,y,z) to [-1,1]x[-1,1]x[-1,1]
    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:
    //convert (x,y,z) to [-1,1]x[-1,1]x[-1,1]
    //set color of this vertex 
    //convert (x,y,z) to [-1,1]x[-1,1]x[-1,1]
    //set color of this vertex 
    //convert (x,y,z) to [-1,1]x[-1,1]x[-1,1]
    //set color of this vertex 
    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:

    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: 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

Push the code and the paper, also bring a hard copy of the paper to class.


Total 25 points [somewhat tentative, I'll need feedback].


Last modified: Sun Nov 5 10:37:29 EDT 2017