New Step-by-Step Tutorial for Velodyne Drone LiDAR from Snoopy by LidarUSA

The folks from Harris Aerial gave us LiDAR data from a test-flight of one of their drones, the Carrier H4 Hybrid HE (with a 5kg maximum payload and a retail price of US$ 28,000), carrying a Snoopy A series LiDAR system from LidarUSA in the countryside near Huntsville, Alabama. The laser scanner used by the Snoopy A series is a Velodyne HDL 32E that has 32 different laser/detector pairs that fire in succession to scan up to 700,000 points per second within a range of 1 to 70 meters. You can download the raw LiDAR file from the 80 second test flight here. As always, the first thing we do is to visualize the file with lasview and to generate a textual report of its contents with lasinfo.

lasview -i Velodyne001.laz -set_min_max 680 750

It becomes obvious that the drone must have scanned parts of itself (probably the landing gear) during the flight and we exploit this fact in the later processing. The information which of the 32 lasers was collecting which point is stored into the ‘point source ID’ field which is usually used for the flightline information. This results in a psychedelic look in lasview as those 32 different numbers get mapped to the 8 different colors that lasview uses for distinguishing flightlines.

The lasinfo report we generate computes the average point density with ‘-cd’ and includes histograms for a number of point attributes, namely for ‘user data’, ‘intensity’, ‘point source ID’, ‘GPS time’, and ‘scan angle rank’.

lasinfo -i Velodyne001.laz ^
        -cd ^
        -histo user_data 1 ^
        -histo point_source 1 ^
        -histo intensity 16 ^
        -histo gps_time 1 ^
        -histo scan_angle_rank 5 ^
        -odir quality -odix _info -otxt

You can download the resulting report here and it will tell you that the information which of the 32 lasers was collecting which point was stored both into the ‘user data’ field and into the ‘point source ID’ field. The warnings you see below have to do with the fact that the double-precision bounding box stored in the LAS header was populated with numbers that have many more decimal digits than the coordinates in the file, which only have millimeter (or millifeet) resolution as all three scale factors are 0.001 (meaning three decimal digits).

WARNING: stored resolution of min_x not compatible with x_offset and x_scale_factor: 2171988.6475160527
WARNING: stored resolution of min_y not compatible with y_offset and y_scale_factor: 1622812.606925504
WARNING: stored resolution of min_z not compatible with z_offset and z_scale_factor: 666.63504345017589
WARNING: stored resolution of max_x not compatible with x_offset and x_scale_factor: 2172938.973065129
WARNING: stored resolution of max_y not compatible with y_offset and y_scale_factor: 1623607.5209975131
WARNING: stored resolution of max_z not compatible with z_offset and z_scale_factor: 1053.092674726669

Both the “return number” and the “number of returns” attribute of every points in the file is 2. This makes it appear as if the file would only contain the last returns of those laser shots that produced two returns. However, as the Velodyne HDL 32E only produces one return per shot we can safely conclude that those numbers should all be 1 instead of 2 and that this is just a small bug in the export software. We can easily fix this with las2las.

reporting minimum and maximum for all LAS point record entries ...
[...]
 return_number 2 2
 number_of_returns 2 2
[...]

The lasinfo report lacks information about the coordinate reference system as there is no VLR that stores projection information. Harris Aerial could not help us other than telling us that the scan was near Huntsville, Alamaba. Measuring certain distances in the scene like the height of the house or the tree suggests that both horizontal and vertical units are in feet, or rather in US survey feet. After some experimenting we find that using EPSG 26930 for NAD83 Alabama West but forcing the default horizontal units from meters to US survey feet gives a result that aligns well with high-resolution Google Earth imagery as you can see below:

lasgrid -i flightline1.laz ^
        -i flightline2.laz ^
        -merged ^
        -epsg 26930 -survey_feet ^
        -step 1 -highest ^
        -false -set_min_max 680 750 ^
        -o testing26930usft.png

Using EPSG code 26930 but with US survey feet instead of meters results in nice alignment with GE imagery.

We use the fact that the drone has scanned itself to extract an (approximate) trajectory by isolating those LiDAR returns that have hit the drone. Via a visual check with lasview (by hovering with the cursor over the lowest drone hits and pressing hotkey ‘i’) we determine that the lowest drone hits are all above 719 feet. We use two calls to las2las to split the point cloud vertically. In the same call we also change the resolution from three to two decimal digits, fix the return number issue, and add the missing geo-referencing information:

las2las -i Velodyne001.laz ^
        -rescale 0.01 0.01 0.01 ^
        -epsg 26930 -survey_feet -elevation_survey_feet ^
        -set_return_number 1 ^
        -set_number_of_returns 1 ^
        -keep_z_above 719 ^
        -odix _above719 -olaz

las2las -i Velodyne001.laz ^
        -rescale 0.01 0.01 0.01 ^
        -epsg 26930 -survey_feet -elevation_survey_feet ^
        -set_return_number 1 ^
        -set_number_of_returns 1 ^
        -keep_z_below 719 ^
        -odix _below719 -olaz

We then use the manual editing capabilities of lasview to change the classifications of the LiDAR points that correspond to drone hits from 1 to 12, which is illustrated by the series of screen shots below.

lasview -i Velodyne001_above719.laz

The workflow illustrated above results in a tiny LAY file that is part of the LASlayers functionality of LAStools. It only encodes the few changes in classifications that we’ve made to the LAZ file without re-writing those parts that have not changed. Those interested may use laslayers to inspect the structure of the LAY file:

laslayers -i Velodyne001_above719.laz

We can apply the LAY file on-the-fly with the ‘-ilay’ option, for example, when running lasview:

lasview -i Velodyne001_above719.laz -ilay

To separate the drone-hit trajectory from the remaining points we run lassplit with the ‘-ilay’ option and request to split by classification with this command line:

lassplit -i Velodyne001_above719.laz -ilay ^
         -by_classification -digits 3 ^
         -olaz

This gives us two new files with the three-digit appendices ‘_001’ and ‘_012’. The latter one contains those points we marked as being part of the trajectory. We now want to use lasview to – visually – find a good moment in time where to split the trajectory into multiple flightlines. The lasinfo report tells us that the GPS time stamps are in the range from 418,519 to 418,602. In order to use the same trick as we did in our recent article on processing LiDAR data from the Hovermap Drone, where we mapped the GPS time to the intensity for querying it via lasview, we first need to subtract a large number from the GPS time stamps to bring them into a suitable range that fits the intensity field as done here.

lasview -i Velodyne001_above719_012.laz ^
        -translate_gps_time -418000 ^
        -bin_gps_time_into_intensity 1
        -steps 5000

The ‘-steps 5000’ argument makes for a slower playback (press ‘p’ to repeat) to better follow the trajectory.

Hovering with the mouse over a point that – visually – seems to be one of the turning points of the drone and pressing ‘i’ on the keyboard shows an intensity value of 548 which corresponds to the GPS time stamp 418548, which we then use to split the LiDAR point cloud (without the trajectory) into two flightlines:

las2las -i Velodyne001_below719.laz ^
        -i Velodyne001_above719_001.laz ^
        -merged ^
        -keep_gps_time_below 418548 ^
        -o flightline1.laz

las2las -i Velodyne001_below719.laz ^
        -i Velodyne001_above719_001.laz ^
        -merged ^
        -keep_gps_time_above 418548 ^
        -o flightline2.laz

Next we use lasoverlap to check how well the LiDAR points from the flight out and the flight back align vertically. This tool computes the difference of the lowest points for each square foot covered by both flightlines. Differences of less than a quarter of a foot are mapped to white, differences of more than half a foot are mapped to saturated red or blue depending on whether the difference is positive or negative:

lasoverlap -i flightline1.laz ^
           -i flightline2.laz ^
           -faf ^
           -min_diff 0.25 -max_diff 0.50 -step 1 ^
           -odir quality -o overlap_025_050.png

We then use a new feature of the LAStools GUI (as of version 180429) to closer inspect larger red or blue areas. We want to use lasmerge and clip out any region that looks suspect for closer examination with lasview. We start the tool in the GUI mode with the ‘-gui’ command and the two flightlines pre-loaded. Using the new PNG overlay roll-out on the left we add the ‘overlap_025_050_diff.png’ image from the quality folder created in the last step and clip out three areas.

lasmerge -i flightline1.laz -i flightline2.laz -gui

You can also clip out these three areas using the command lines below:

lasmerge -i flightline1.laz -i flightline2.laz ^
         -faf ^
         -inside 2172500 1623160 2172600 1623165 ^
         -o clip2500_3160_100x005.laz

lasmerge -i flightline1.laz -i flightline2.laz ^
         -faf ^
         -inside 2172450 1623450 2172550 1623455 ^
         -o clip2450_3450_100x005.laz

lasmerge -i flightline1.laz -i flightline2.laz ^
         -faf ^
         -inside 2172430 1623290 2172530 1623310 ^
         -o clip2430_3290_100x020.laz

A closer inspection of the three cut out slices explains the red and blue areas in the difference image created by lasoverlap. We find a small systematic error in two of the slices. In slice ‘clip2500_3160_100x005.laz‘ the green points from flightline 1 are on average slightly higher than the red points from flightline 2. Vice-versa in slice ‘clip2450_3450_100x005.laz‘ the green points from flightline 1 are on average slightly lower than the red points from flightline 2. However, the error is less than half a foot and only happens near the edges of the flightlines. Given that our surfaces are expected to be “fluffy” anyways (as is typical for Velodyne LiDAR systems), we accept these differences and continue processing.

The strong red and blue colors in the center of the difference image created by lasoverlap is easily explained by looking at slice ‘clip2430_3290_100x020.laz‘. The scanner was “looking” under a gazebo-like open roof structure from two different directions and therefore always seeing parts of the floor in one flightline that were obscured by the roof in the other.

While working with this data we’ve also implemented a new feature for lastrack that computes the 3D distance between LiDAR points and the trajectory and allows storing the result as an additional per point attribute with extra bytes. Those can then be visualized with lasgrid. Here an example:

lastrack -i flightline1.laz ^
         -i flightline2.laz ^
         -track Velodyne001_above719_012.laz ^
         -store_xyz_range_as_extra_bytes ^
         -odix _xyz_range -olaz ^
         =cores 2

lasgrid -i flightline*_xyz_range.laz -merged ^
        -drop_attribute_below 0 1 ^
        -attribute0 -lowest ^
        -false -set_min_max 20 200 ^
        -o quality/closest_xyz_range_020ft_200ft.png

lasgrid -i flightline*_xyz_range.laz -merged ^
        -drop_attribute_below 0 1 ^
        -attribute0 -highest ^
        -false -set_min_max 30 300 ^
        -o quality/farthest_xyz_range_030ft_300ft.png

Below the complete processing pipeline for creating a median ground model from the “fluffy” Velodyne LiDAR data that results in the hillshaded DTM shown here. The workflow is similar to those we have developed in earlier blog posts for Velodyne Puck based systems like the Hovermap and the Yellowscan.

Hillshaded DTM with a resolution of 1 foot generated via a median ground computation by the LAStools processing pipeline detailed below.

lastile -i flightline1.laz ^
        -i flightline2.laz ^
        -faf ^
        -tile_size 250 -buffer 25 -flag_as_withheld ^
        -odir tiles_raw -o somer.laz

lasnoise -i tiles_raw\*.laz ^
         -step_xy 2 -step 1 -isolated 9 ^
         -odir tiles_denoised -olaz ^
          -cores 4

lasthin -i tiles_denoised\*.laz ^
        -ignore_class 7 ^
        -step 1 -percentile 0.05 10 -classify_as 8 ^
        -odir tiles_thinned_1_foot -olaz ^
        -cores 4

lasthin -i tiles_thinned_1_foot\*.laz ^
        -ignore_class 7 ^
        -step 2 -percentile 0.05 10 -classify_as 8 ^
        -odir tiles_thinned_2_foot -olaz ^
        -cores 4

lasthin -i tiles_thinned_2_foot\*.laz ^
        -ignore_class 7 ^
        -step 4 -percentile 0.05 10 -classify_as 8 ^
        -odir tiles_thinned_4_foot -olaz ^
        -cores 4

lasthin -i tiles_thinned_4_foot\*.laz ^
        -ignore_class 7 ^
        -step 8 -percentile 0.05 10 -classify_as 8 ^
        -odir tiles_thinned_8_foot -olaz ^
        -cores 4

lasground -i tiles_thinned_8_foot\*.laz ^
          -ignore_class 1 7 ^
          -town -extra_fine ^
          -odir tiles_ground_lowest -olaz ^
          -cores 4

lasheight -i tiles_ground_lowest\*.laz ^
          -classify_between -0.05 0.5 6 ^
          -odir tiles_ground_thick -olaz ^
          -cores 4

lasthin -i tiles_ground_thick\*.laz ^
        -ignore_class 1 7 ^
        -step 1 -percentile 0.5 -classify_as 2 ^
        -odir tiles_ground_median -olaz ^
        -cores 4

las2dem -i tiles_ground_median\*.laz ^
        -keep_class 2 ^
        -step 1 -use_tile_bb ^
        -odir tiles_dtm -obil ^
        -cores 4

blast2dem -i tiles_dtm\*.bil -merged ^
          -step 1 -hillshade ^
          -o dtm_hillshaded.png

We thank Harris Aerial for sharing this LiDAR data set with us flown by their Carrier H4 Hybrid HE drone carrying a Snoopy A series LiDAR system from LidarUSA.

First Look with LAStools at LiDAR from Hovermap Drone by CSIRO

Last December we had a chance to visit the team of Dr. Stefan Hrabar at CSIRO in Pullenvale near Brisbane who work on a drone LiDAR system called Hovermap. This SLAM-based system is mainly developed for the purpose of autonomous flight and exploration of GPS-denied environments such as buildings, mines and tunnels. But as the SLAM algorithm continuously self-registers the scan lines it produces a LiDAR point cloud that in itself is a nice product. We started our visit with a short test flight around the on-site tower. You can download the LiDAR data and the drone trajectory of this little survey here:

The Hovermap system is based on the Velodyne Puck Lite (VLP-16) that is much cheaper and more light-weight than many other LiDAR systems. One interesting tidbit in the Hovermap setup is that the scanner is installed such that the entire Puck is constantly rotating as you can see in this video. But  the Velodyne Puck is also known to produce somewhat “fluffy” surfaces with a thickness of a few centimeters. In a previous blog post with data from the YellowScan Surveyor system (that is also based on the Puck) we used a “median ground” surface to deal with the “fluff”. In the following we will have a look at the LiDAR data produced by Hovermap and how to process it with LAStools.

LiDAR data of CSIRO tower acquired during test flight of Hovermap system.

As always we start with a lasinfo report that computes the average density ‘-cd’ and histograms for the intensity and the GPS time:

lasinfo -i CSIRO_Tower\results.laz ^
        -cd ^
        -histo intensity 16 -histo gps_time 2 ^
        -odir CSIRO_Tower\quality -odix _info -otxt

A few excerpts of the resulting lasinfo report that you can download here are below:

lasinfo (180409) report for 'CSIRO_Tower\results.laz'
[...]
 number of point records: 16668904
 number of points by return: 0 0 0 0 0
 scale factor x y z: 0.0001 0.0001 0.0001
 offset x y z: -5.919576153930379 22.785394470724583 9.535698734939086
 min x y z: -138.6437 -125.2552 -34.1510
 max x y z: 126.8046 170.8260 53.2224
WARNING: full resolution of min_x not compatible with x_offset and x_scale_factor: -138.64370561381907
WARNING: full resolution of min_y not compatible with y_offset and y_scale_factor: -125.25518631070418
WARNING: full resolution of min_z not compatible with z_offset and z_scale_factor: -34.150966206894068
WARNING: full resolution of max_x not compatible with x_offset and x_scale_factor: 126.80455330595831
WARNING: full resolution of max_y not compatible with y_offset and y_scale_factor: 170.82597525215334
WARNING: full resolution of max_z not compatible with z_offset and z_scale_factor: -34.150966206894068
[...]
 gps_time 121.288045 302.983110
WARNING: 2 points outside of header bounding box
[...]
covered area in square units/kilounits: 51576/0.05
point density: all returns 323.19 last only 318.40 (per square units)
 spacing: all returns 0.06 last only 0.06 (in units)
WARNING: for return 1 real number of points by return is 16424496 but header entry was not set.
WARNING: for return 2 real number of points by return is 244408 but header entry was not set.
[...]
real max z larger than header max z by 0.000035
real min z smaller than header min z by 0.000035
[...]

Most of these warnings have to do with poorly chosen offset values in the LAS header that have many decimal digits instead of being nice round numbers. The points are stored with sub-millimeter resolution (scale factors of 0.0001) which is unnecessarily precise for a UAV flying a Velodyne Puck where the overall system error can be expected to be on the order of a few centimeters. Also the histogram of return numbers in the LAS header was not populated. We can fix these issues with one call to las2las:

las2las -i CSIRO_Tower\results.laz ^
        -rescale 0.01 0.01 0.01 ^
        -auto_reoffset ^
        -odix _fixed -olaz

If you create another lasinfo report on the fixed file you will see that all the warnings have gone. The file size is now also only 102 MB instead of 142 MB because centimeter coordinate compress much better than sub-millimeter coordinates.

The average density of 318 last return per square meter reported by lasinfo is not that useful for a UAV survey because it does account for the highly varying distribution of LiDAR returns in the area surveyed. With lasgrid we can get a much more clear picture of that.

lasgrid -i CSIRO_Tower\results_fixed.laz ^
        -last_only ^
        -step 0.5 -use_bb -density ^
        -false -set_min_max 0 1500 ^
        -o CSIRO_Tower\quality\density_0_1500.png

LiDAR density: blue is close to zero and red is 1500 or more last returns / sqr mtr

The red dot in the point density indicated an area with over 1500 last returns per square meter. No surprise that this is the take-off and touch-down location of the copter drone. Naturally this spot is completely over-scanned compared to the rest of the area. We can remove these points with the help of the timestamps by cutting off the start and the end of the recording.

The total recording time including take-off, flight around the tower, and touch-down was around 180 seconds or 3 minutes as the lasinfo report tells us. Note that the recorded time stamps are neither “GPS Week Time” nor “Adjusted Standard GPS Time” but an internal system time. By visualizing the trajectory of the UAV with lasview while binning the timestamps into the intensity field we can easily determine what interval of timestamps describes the actual survey flight. First we convert the drone trajectory from the textual ASCII format to the LAZ format with txt2las:

txt2las -i CSIRO_Tower\results_traj.txt ^
        -skip 1 ^
        -parse txyz ^
        -set_classification 12 ^
        -olaz

lasview -i CSIRO_Tower\results_traj.laz ^
        -bin_gps_time_into_intensity 1

Binning timestamps into intensity allows visually determining start and end of survey.

Using lasview and pressing <i> while hovering over those points of the trajectory that appear to be the survey start and end we determine visually that the timestamps between 164 to 264 correspond to the actual survey flight over the area of interest with the take-off and touch-down maneuvers excluded. We use las2las to cut out the relevant part and re-run lasgrid:

las2las -i CSIRO_Tower\results_fixed.laz ^
        -keep_gps_time 164 264 ^
        -o CSIRO_Tower\results_survey.laz

lasgrid -i CSIRO_Tower\results_survey.laz ^
        -last_only ^
        -step 0.5 -use_bb -density ^
        -false -set_min_max 0 1500 ^
        -o CSIRO_Tower\quality\density_0_1500_survey.png

LiDAR density after removing take-off and touch-down maneuvers.

The other set of point we are less interested in are those occasional hits far from the scanner that sample the area too sparsely to be useful for anything. We use lastrack to reclassify points as noise (7) that exceed a x/y distance of 50 meters from the trajectory and then use lasgrid to create another density image without the points classified as noise..

lastrack -i CSIRO_Tower\results_survey.laz ^
         -track CSIRO_Tower\results_traj.laz ^
         -classify_xy_range_between 50 1000 7 ^
         -o CSIRO_Tower\results_xy50.laz

lasgrid -i CSIRO_Tower\results_xy50.laz ^
        -last_only -keep_class 0 ^
        -step 0.5 -use_bb -density ^
        -false -set_min_max 0 1500 ^
        -o CSIRO_Tower\quality\density_0_1500_xy50.png

LiDAR density after removing returns farther than 50 m from trajectory.

We process the remaining points using a typical tile-based processing pipeline. First we run lastile to create tiling of 200 meter by 200 meter tiles with 20 buffers while dropping the noise points::

lastile -i CSIRO_Tower\results_xy50.laz ^
        -drop_class 7 ^
        -tile_size 200 -buffer 20 -flag_as_withheld ^
        -odir CSIRO_Tower\tiles_raw -o eta.laz

Because of the high sampling we expect there to be quite a few duplicate point where all three coordinate x, y, and z are identical. We remove them with a call to lasduplicate:

lasduplicate -i CSIRO_Tower\tiles_raw\*.laz ^
             -unique_xyz ^
             -odir CSIRO_Tower\tiles_unique -olaz ^
             -cores 4

This removes between 12 to 25 thousand point from each tile. Then we use lasnoise to classify isolated points as noise:

lasnoise -i CSIRO_Tower\tiles_unique\*.laz ^
         -step_xy 0.5 -step_z 0.1 -isolated 5 ^
         -odir CSIRO_Tower\tiles_denoised_temp -olaz ^
         -cores 4

Aggressive parameters assure most noise point below ground are found.

This classifies between 13 to 23 thousand point from each tile into the noise classification code 7. We use rather aggressive settings to make sure we get most of the noise points that are below the terrain. Getting a correct ground classification in the next few steps is the main concern now even if this means that many points above the terrain on wires, towers, or vegetation will also get miss-classified as noise (at least temporarily). Next we use lasthin to classify a subset of points with classification code 8 on which we will then run the ground classification. We classify each point that is closest to the 5th percentile in elevation per 25 cm by 25 cm grid cell given there are at least 20 non-noise points in a cell. We then repeat this while increasing the cell size to 50 cm by 50 cm and 100 cm by 100 cm.

lasthin -i CSIRO_Tower\tiles_denoised_temp\*.laz ^
        -ignore_class 7 ^
        -step 0.25 -percentile 5 20 -classify_as 8 ^
        -odir CSIRO_Tower\tiles_thinned_025 -olaz ^
        -cores 4

lasthin -i CSIRO_Tower\tiles_thinned_025\*.laz ^
        -ignore_class 7 ^
        -step 0.50 -percentile 5 20 -classify_as 8 ^
        -odir CSIRO_Tower\tiles_thinned_050 -olaz ^
        -cores 4

lasthin -i CSIRO_Tower\tiles_thinned_025\*.laz ^
        -ignore_class 7 ^
        -step 1.00 -percentile 5 20 -classify_as 8 ^
        -odir CSIRO_Tower\tiles_thinned_100 -olaz ^
        -cores 4

 

Then we ground classify the points that were classified into the temporary classification code 8 in the previous step using lasground.

lasground -i CSIRO_Tower\tiles_thinned_100\*.laz ^
          -ignore_class 7 0 ^
          -town -ultra_fine ^
          -odir CSIRO_Tower\tiles_ground -olaz ^
          -cores 4

The resulting ground points are a lower envelope of the “fluffy” sampled surfaces produced by the Velodyne Puck scanner. We use lasheight to thicken the ground by moving all points between 1 cm below and 6 cm above the TIN of these “low ground” points to a temporary classification code 6 representing a “thick ground”. We also undo the overly aggressive noise classifications above the ground by setting all higher points back to classification code 1 (unclassified).

lasheight -i CSIRO_Tower\tiles_ground\*.laz ^
          -classify_between -0.01 0.06 6 ^
          -classify_above 0.06 1 ^
          -odir CSIRO_Tower\tiles_ground_thick -olaz ^
          -cores 4

Profile view for 25 centimeter wide strip of open terrain. Top: Green points are low ground. Orange points are thickened ground with 5 cm drop lines. Middle: Brown points are median ground computed from thick ground. Bottom: Comparing low ground points (in green) with median ground points (in brown).

From the “thick ground” we then compute a “median ground” using lasthin in a similar fashion as we used it before. A profile view for a 25 centimeter wide strip of open terrain illustrates the workflow of going from “low ground” via “thick ground” to “median ground” and shows the slight difference in elevation between the two.

lasthin -i CSIRO_Tower\tiles_ground_thick\*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.25 -percentile 50 10 -classify_as 2 ^
        -odir CSIRO_Tower\tiles_ground_median_025 -olaz ^
        -cores 4

lasthin -i CSIRO_Tower\tiles_ground_median_025\*.laz ^
        -ignore_class 0 1 7 ^
        -step 0.50 -percentile 50 10 -classify_as 2 ^
        -odir CSIRO_Tower\tiles_ground_median_050 -olaz ^
        -cores 4

lasthin -i CSIRO_Tower\tiles_ground_median_050\*.laz ^
        -ignore_class 0 1 7 ^
        -step 1.00 -percentile 50 10 -classify_as 2 ^
        -odir CSIRO_Tower\tiles_ground_median_100 -olaz ^
        -cores 4

Then we use lasnoise once more with more conservative settings to remove the noise points that are sprinkled around the scene.

lasnoise -i CSIRO_Tower\tiles_ground_median_100\*.laz ^
         -step_xy 1.0 -step_z 1.0 -isolated 5 ^
         -odir CSIRO_Tower\tiles_denoised -olaz ^
         -cores 4

While we classify the scene into building roofs, vegetation, and everything else with lasclassify we also move all (unused) classifications to classification code 1 (unclassified). You may play with the parameters of lasclassify (see README) to achieve better a building classification. However, those buildings the laser can peek into (either via a window or because they are gazebo-like structures) will not be classified correctly. unless you remove the points that are under the roof somehow.

lasclassify -i CSIRO_Tower\tiles_denoised\*.laz ^
            -ignore_class 7 ^
            -change_classification_from_to 0 1 ^
            -change_classification_from_to 6 1 ^
            -step 1 ^
            -odir CSIRO_Tower\tiles_classified -olaz ^
            -cores 4

A glimpse at the final classification result is below. A hillshaded DTM and a strip of classified points. Of course the tower was miss-classified as vegetation given that it looks just like a tree to the logic used in lasclassify.

The hillshaded DTM with a strip of classified points.

Finally we remove the tile buffers (that were really important for tile-based processing) with lastile:

lastile -i CSIRO_Tower\tiles_classified\*.laz ^
        -remove_buffer ^
        -odir CSIRO_Tower\tiles_final -olaz ^
        -cores 4

And publish the LiDAR point cloud as version 1.6 of Potree using laspublish:

laspublish -i CSIRO_Tower\tiles_final\*.laz ^
           -i CSIRO_Tower\results_traj.laz ^
           -only_3D -elevation -overwrite -potree16 ^
           -title "CSIRO Tower" ^
           -description "HoverMap test flight, 18 Dec 2017" ^
           -odir CSIRO_Tower\tiles_portal -o portal.html -olaz

Note that we also added the trajectory of the drone because it looks nice and gives a nice illustration of how the UAV was scanning the scene.

Via Potree we can publish and explore the final point cloud using any modern Web browser.

We would like to thank the entire team around Dr. Stefan Hrabar for taking time out of their busy schedules just a few days before Christmas.