Creating a Better DTM from Photogrammetic Points by Avoiding Shadows

At INTERGEO 2016 in Hamburg, the guys from Aerowest GmbH shared with us a small photogrammetric point cloud from the city of Soest that had been generated with the SURE dense-matching software from nFrames. We want to test whether – using LAStools – we can generate a decent DTM from these points that are essentially a gridded DSM. If this interest you also see this, this, this, and this story.

soest_00_google_earth

Here you can download the four original tiles (tile1, tile2, tile3, tile4) that we are using in these experiments. We first re-tile them into smaller 100 meter by 100 meter tiles with a 20 meter buffer using lastile. We use option ‘-flag_as_withheld’ that flags all the points falling into the buffer using the withheld flag so they can easily be removed on-the-fly later with the ‘-drop_withheld’ filter (see the README for more on this). We also add the missing projection with ‘-epsg 32632’.

lastile -i original\*.laz ^
        -tile_size 100 -buffer 20 ^
        -flag_as_withheld ^
        -epsg 32632 ^
        -odir tiles_raw -o soest.laz

Below is a screenshot from one of the resulting 100 meter by 100 meter tiles (with 20 meter buffer) that we will be focusing on in the following experiments.

The tiles 'soest_437900_5713800.laz'

The tile ‘soest_437900_5713800.laz’ used in all experiments.

Next we use lassort to reorder the points ordered along a coherent space-filling curve as the existing scan-line order has the potential to cause our triangulation engine to slow down. We do this on 4 cores.

lassort -i tiles_raw\*.laz ^
        -odir tiles_sorted -olaz ^
        -cores 4

We then use lasthin to lower the number of points that we consider as ground points (see the README for more info on this tool). We do this because the original 5 cm spacing of the dense matching points is a bit excessive for generating a DTM with a resolution of, for example, 50 cm. Instead we only use the lowest point in each 20 cm by 20 cm cell as a candidate for becoming a ground point, which reduces the number of considered points by a factor of 16. We achieve this by classifying these lowest point to a unique classification code (here: 9) and later tell lasground to ignore all other classifications.

lasthin -i tiles_sorted\*.laz ^
        -step 0.2 -lowest -classify_as 9 ^
        -odir tiles_thinned -olaz ^
        -cores 4
Then we run lasground on 4 cores to classify the ground points with options ‘-step 20’, ‘-bulge 0.5’, ‘-spike 0.5’ and ‘-fine’ that we selected after two trials on a single tile. There are several other options in lasground to play with that may achieve better results on other data sets (see README file for more on this). The ‘-ignore_class 0’ option instructs lasground to ignore all points that are not classified so that only those points that lasthin was classifying as 9 in the previous step are considered.
lasground -i tiles_thinned\*.laz ^
          -step 20 -bulge 0.5 -spike 0.5 -fine ^
          -ignore_class 0 ^
          -odir tiles_ground -olaz ^
          -cores 4
However, when we scrutinize the resulting ground classification notice that there are bumps in the corresponding ground TIN that seem to correspond to areas where the original imagery has dark shadows that in turn seem to significantly affect the geometric accuracy of the points generated by the dense-matching algorithm.
Looking a the bump from below we identify the RGB colors of the points have that form the bump that seem to be of much lower accuracy than the rest. This is an effect that we have noticed in the past for data generated with other dense-matching software and maybe an approach similar to the one we take here could also help with this low noise. It seems that points that are generated from shadowed areas in the input images can have a lot lower accuracy than those from well-lit areas. We use this correlation between RGB color and geometric accuracy to simply exclude all points whose RGB colors indicate that they might be from shadow areas from becoming ground points.
The RGB colors of low-accuracy points are mostly from very dark shadowed areas.

The RGB colors of low-accuracy points are mostly from very dark shadowed areas.

We use las2las with the relatively new ‘-filtered_transform’ option to reclassify all points whose RGB color is close to zero to yet classification code 7 (see README file for more on this). We do this for all points whose red value is between 0 and 30, whose green value is between 0 and 35, and whose blue value is between 0 and 40. Because the RGB values were stored with 16 bits in these files we have to multiply these values with 256 when applying the filter.
las2las -i tiles_thinned\*.laz ^
        -keep_RGB_red 0 7680 ^
        -keep_RGB_green 0 8960 ^
        -keep_RGB_blue 0 10240 ^
        -filtered_transform ^
        -set_classification 7 ^
        -odir tiles_rgb_filtered -olaz ^
        -cores 4
Below you see all points that will no longer be considered because their classification was set to 7 by the command above.
Points whose RGB values indicate they might lie in the shadows.

Points whose RGB values indicate they might lie in the shadows.

We then re-run lasground with the same options ‘-step 20’, ‘-bulge 0.5’, ‘-spike 0.5’ and ‘-fine’ as before but now we ignore all points that are still have classification 0 because they were not classified as 9 by lasthin earlier and we also ignore all points that have been assigned classification 7 by las2las in the previous step.
lasground -i tiles_thinned\*.laz ^
          -step 20 -bulge 0.5 -spike 0.5 -fine ^
          -ignore_class 0 7 ^
          -odir tiles_ground -olaz ^
          -cores 4
The situation has significantly improved for the bumb we saw before as you can see in the images below.

Finally we create a DTM with blast2dem (see README) and a DSM with lasgrid (see README) by merging all points into one file but dropping the buffer points that were marked as withheld by the initial run of lastile (see README).

blast2dem -i tiles_ground\*.laz -merged ^
          -drop_withheld -keep_class 2 ^
          -step 0.5 ^
          -o dtm.bil

lasgrid -i tiles_ground\*.laz -merged ^
        -drop_withheld ^
        -step 0.5 -average ^
        -o dsm.bil
 As our venerable lasview (see README) can directly read BIL rasters as points (just like all the other LAStools), so we can compare the DTM and the DTM by loading them as two flight lines (‘-faf’) and then switch back and forth between the two by pressing ‘0’ and ‘1’ in the viewer.
lasview -i dtm.bil dsm.bil -faf

Above you see the final DTM and the original DSM. So yes, LAStools can definitely create a DTM from point clouds that are the result of dense-matching photogrammetry. We used the correlation between shadowed areas in the source image and geometric errors to remove those points from consideration for ground points that are likely to be too low and result in bumps. For dense-matching algorithms that also output an uncertainty value for each point there is the potential for even better results as our range of eliminated RGB colors may not cover all geometrically uncertain points while at the same time may be too conservative and also remove correct ground points.

LASmoons: Jakob Iglhaut

Jakob Iglhaut (recipient of three LASmoons)
Program for Geospatial Information Management
Carinthia University of Applied Sciences, Villach, AUSTRIA

Background:
As part of the EU LIFE programme two river stretches in Carinthia, Austria have recently been subject to restoration measures. The LIFE-project aims at protecting valuable riverine flora and fauna while improving flood protection. By remodelling the river beds, the construction of groynes and still water bodies the river environment was directed to more natural morphology and state. The joint R&D project “Remotely Piloted Aircraft Multi Sensor System (RPAMSS)” aims at capturing multi-dimensional environmental data in order to monitor the development of these rivers stretches in a holistic way. Flights with an RTK capable fixed wing UAV are carried out at a particular section of the rivers Gail and Drau respectively. The project site at the Upper-Drau is located in the area of Obergottesfeld, Austria (560m ASL), with an area currently remotely monitored by the RPAmSS of approximately 3.5km². The second study area is located close to Feistritz at the river Gail (550m ASL) with an area of approx. 0.9km². Apart from being addressed by the LIFE project both study areas are also defined as NATURA 2000 nature protection sites. In both areas frequent UAV flights are carried out collecting high-resolution multi-spectral imagery. Structure from Motion photogrammetry enables the creation of high-density multi-spectral point clouds.

lasmoons_jakob_iglhaut_0

Goal:
The aim of the project is to assess the morphology and related temporal changes of the described riverine environment based on SfM point clouds. A full processing chain will be developed to take full advantage of the high-density data. Particular interest lies in the extraction of ground points underneath vegetation in leaf-on/leaf-off. Ground points will be gridded to generate DTMs. The qualitative performance of the data will be held against an ALS acquired DTM. Furthermore forest metrics will be extracted for the riparian zone in order to quantify their current state and changes.

Data:
+
High-density multi-spectral (R,G,B,NIR) SfM derived point clouds (UAS imagery)
+ Variable point densities, GSD ~3cm.

LAStools processing:
1) 
fix SfM owing incoherence [lassort]
2) create 100m tiles (10m buffer) for parallel processing [lastile]
3) noise removal introduced by the SfM algorithm [lasnoise]
4).extract ground points [lasground_new]
5) generate normalized above heights [lasheight]
6) classify based on height-above-ground (low veg, high veg) [lasheight]
7) create DSM and DTM [blast2dem]
8) 
generate a Canopy Height Model (CHM) using the pit-free method of Khosravipour et al. (2014) with the workflow described here [lasthin, las2dem, lasgrid]
9) 
sub-sample the point clouds for other (spectral) analyses [lassplit, lasthin, lasmerge]

Reference:
Westoby, M. J., et al. “Structure-from-Motion photogrammetry: A low-cost, effective tool for geoscience applications.” Geomorphology 179 (2012): 300-314.
Fonstad, Mark A., et al. “Topographic structure from motion: a new development in photogrammetric measurement.” Earth Surface Processes and Landforms 38.4 (2013): 421-430.
Khosravipour, A., Skidmore, A.K., Isenburg, M., Wang, T.J., Hussin, Y.A., 2014. Generating pit-free Canopy Height Models from Airborne LiDAR. PE&RS = Photogrammetric Engineering and Remote Sensing 80, 863-872.
Javernick, L., J. Brasington, and B. Caruso. “Modeling the topography of shallow braided rivers using Structure-from-Motion photogrammetry.” Geomorphology 213 (2014): 166-182.

Removing low noise from Semi-Global Matching (SGM) Points

At PhoWo and INTERGEO 2015 rapidlasso was spending quality time with VisionMap who make the A3 Edge camera that provides fine resolution images from high altitudes and can quickly cover large areas. Under the hood of their LightSpeed software is the SURE dense-matching algorithm from nframes that turns images into photogrammetric point clouds. We were asked whether LAStools is able to create bare-earth DTM rasters from such points. If you have read our first, second, or third blog post on the topic you know that our asnwer was a resounding “YES!”. But we ran into an issue due to the large amount of low noise. Maybe the narrow angle between images at a high flying altitude affects the semi-global matching (SGM) algorithm. Either way, in the following we show how we use lascanopy and lasheight to mark low points as noise in a preprocessing step.

We obtained a USB stick containing a 2.42 GB file called “valparaiso_DSM_SURE_100.las” containing about 100 million points spaced 10 cm apart generated by SURE and stored with an (unnecessary high) resolution of millimeters (aka “resolution fluff”) as the third digit of all coordinates was always zero:

las2txt -i F:\valparaiso_DSM_SURE_100.las -stdout | more
255991.440 6339659.230 89.270
255991.540 6339659.240 89.270
255991.640 6339659.240 88.660
255991.740 6339659.230 88.730
[...]

We first compressed the bulky 2.42 GB LAS file into a compact 0.23 GB LAZ to our local hard drive – a file that is 10 times smaller and that will be 10 times faster to copy:

laszip -i F:\valparaiso_DSM_SURE_100.las ^
       -rescale 0.01 0.01 0.01 ^
       -o valparaiso_DSM_SURE_100.laz ^

Then we tiled the 100 million points into 250 meter by 250 meter tiles with 25 meter buffer using lastile. We use the new option ‘-flag_as_withheld’ to mark all buffer points with the withheld flag so they can be easily removed on-the-fly via the ‘-drop_withheld’ command-line filter (also see the README file).

lastile -i valparaiso_DSM_SURE_100.laz ^
        -tile_size 250 -buffer 25 ^
        -flag_as_withheld ^
        -odir valparaiso_tiles_raw -o val.laz
250 meter by 250 meter tiling with 25 meter buffer

250 meter by 250 meter tiles with 25 meter buffer

Before processing millions to billions of points we experiment with different options to find what works best on a smaller area, namely the tile “val_256750_6338500.laz” pointed to above. Using the workflow from this blog posts did not give perfect results due to the large amount of low noise. Although many low points were marked as noise (violett) by lasnoise, too many ended up classified as ground (brown) by lasground as seen here:
excessive low noise affects ground classification

excessive low noise affects ground classification

We use lascanopy – a tool very popular with forestry folks – to compute four BIL rasters where each 5m by 5m grid cell contains the 5th, 10th, 15th, and 20th percentile of the elevation values from all points falling into a cell (also see the README file):
lascanopy -i val_256750_6338500.laz ^
          -height_cutoff -1000 -step 5 ^
          -p 5 10 15 20 ^
          -obil
The four resulting rasters can be visually inspected and compared with lasview:
lasview -i val_256750_6338500_*.bil -files_are_flightlines
comparing 5th and 10th elevation percentiles

comparing the 5th and the 10th elevation percentiles

By pressing the hot keys <0>, <1>, <2> and <3> to switch between the different percentiles and <t> to triangulate them into a surface, we can see that for this example the 10th percentile works well while the 5th percentile is still affected by the low noise. Hence we use the 10th percentile elevation surface and classify all points below it as noise with lasheight (also see the README file).
lasheight -i val_256750_6338500.laz ^
          -ground_points val_256750_6338500_p10.bil ^
          -classify_below -0.5 7 ^
          -odix _denoised -olaz
We visually confirm that all low points where classified as noise (violett). Note the obvious “edge artifact” along the front boundary of the tile. This is why we always recommend to use a buffer in tile-based processing.
points below 10th percentile surface marked as noise

points below 10th percentile surface marked as noise

At the end of the blog post we give the entire command sequence that first computes a 10th percentile raster with 5m resolution for the entire file with lascanopy and then marks all points of each tile below the10th percentile surface as noise with lasheight. When we classify all points into ground and non-ground points with lasground we ignore all points classified as noise. Here are the results:

DTM extracted from SGM points despite low noise

DTM extracted from dense-matching points despite low noise

corresponding DSM with all buildings and vegetaion included

corresponding DSM with all buildings and vegetaion included

Above you see the generated DTM and the corresponding DSM. So yes, LAStools can create DTMs from points that are result of dense-matching photogrammetry … even when there is a lot of low noise. There are many other ways to mix and match the modules of LAStools for more refined workflows. Sometimes declaring all points below the 10th percentile surface as noise may be too agressive. In a future blog post we will look how to combine lascanopy and lasnoise for a more adaptive approach.

:: compute 10th percentile for entire area
lascanopy -i valparaiso_DSM_SURE_100.laz ^
          -height_cutoff -1000 -step 5 ^
          -p 10 ^
          -obil

:: tile input into 250 meter tiles with buffer
lastile -i valparaiso_DSM_SURE_100.laz ^
        -tile_size 250 -buffer 25 ^
        -flag_as_withheld ^
        -odir valparaiso_tiles_raw -o val.laz

:: mark points below as noise
lasheight -i valparaiso_tiles_raw/*.laz ^
          -ground_points valparaiso_DSM_SURE_100_p10.bil ^
          -classify_below -0.5 7 ^
          -odir valparaiso_tiles_denoised -olaz ^
          -cores 4

:: ground classify while ignoring noise points
 lasground -i valparaiso_tiles_denoised\*.laz ^
          -ignore_class 7 ^
          -town -bulge 0.5 ^
          -odir valparaiso_tiles_ground -olaz ^
          -cores 4 

:: create 50 cm DTM rasters in BIL format
las2dem -i valparaiso_tiles_ground\*.laz ^
        -keep_class 2 ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir valparaiso_tiles_dtm -obil ^
        -cores 4 

:: average 50 cm DTM values into single 1m DTM 
lasgrid -i valparaiso_tiles_dtm\*.bil -merged ^
        -step 1.0 -average ^
        -o valparaiso_dtm.bil

:: create hillshade adding in UTM 19 southern
blast2dem -i valparaiso_dtm.bil ^
          -hillshade -utm 19M ^
          -o valparaiso_dtm_hill.png

:: create DSM hillshade with same three steps
las2dem -i valparaiso_tiles_raw\*.laz ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir valparaiso_tiles_dsm -obil ^
        -cores 4
lasgrid -i valparaiso_tiles_dsm\*.bil -merged ^
        -step 1.0 -average ^
        -o valparaiso_dsm.bil
blast2dem -i valparaiso_dsm.bil ^
          -hillshade -utm 19M ^
          -o valparaiso_dsm_hill.png

Creating DTMs from dense-matched points of UAV imagery from SenseFly’s eBee

Tim Sutton and his team at Kartoza work on flood modelling and risk assessment using Inasafe. They have been trying to generate a DTM from point cloud data derived via dense-matching from UAV imagery collected by an eBee of SenseFly in the “unplanned developments” or “slums” North West of Dar es Salaam, the capital city of Tanzania. Tim’s team was stuck after “other software” produced this result:

results for ground points classification with other software

poor ground classification of “Tandale” with “other software”

Tim reached out to us at rapidlasso asking whether LAStools could handle this better. After all, we had published two blog articles – namely this one and that one – showing how to generate DTMs from the point clouds generated by the dense-matching photogrammetry software of Pix4D. Below the workflow we devised and the results we produced for Tim and his team.

We obtained 3 different data sets of areas called “Tandale”, “Borahatward”, and “Bugurunni”. We added one new option to our lasground software called ‘-bulge 1.0’ (see README) to improve the removal of smaller buildings and got this result for “Tandale”.

ground classification with LAStools

DTM of “Tandale” from ground points classified with LAStools

Before you point out the “facetted” look of this DTM keep in mind that “Tandale” is a densely populated poor area. A first hand account of the rough life in this area can be found here. Most dense-matching points are on corrugated roofs that become voids that need to be interpolated across in the DTM. Take a look at the corresponding DSM where all objects are still present.

original data

DSM of “Tandale” from all dense-matching points

Below we give a detailed description at the example of the “Bugurunni” data set of the workflow that was used to generate DTMs for the three data sets. At the end of this article you will see some more results.

We first use lassort to quantize, sort, and compress on 4 cores the seven spatially incoherent LAS files of the “Bugurunni” data set (totalling 4.5 GB with excessive resolution of millimeters) into LASzip-compressed files with a more reasonable resolution of centimeters and points ordered along a space-filling curve. We also add the missing projection information with ‘-utm 37M’. The resulting 7 LAZ files occupy only 0.7 GB meaning we get a compression of 9 : 1. The option ‘-odir’ specifies the output directory.

lassort -i bugurunni_densified_*.las ^
        -rescale 0.01 0.01 0.01 ^
        -utm 37M ^
        -odir bugurunni_strips -olaz ^
        -cores 4

Next we tile the sorted strips into 500 meter by 500 meter tiles with 50 meter buffer using lastile. We use the new option ‘-flag_as_withheld’ to mark all buffer points with the withheld flag so they can easily be removed on-the-fly with the ‘-drop_withheld’ command-line filter (see the README file for more on this).

lastile -i bugurunni_strips\*.laz ^
        -files_are_flightlines ^
        -tile_size 500 -buffer 50 ^
        -flag_as_withheld ^
        -o bugurunni_raw\bugu.laz
Using lasnoise on 4 cores we classify isolated points that might hinder ground-classification as noise (class 7). The parameters ‘-isolated 15’ means that all points surrounded by less than 15 other points in their 3 by 3 by 3 = 27 cells neighborhood in a 3D grid are considered isolated. The size of each grid cell is specified with ‘-step_xy 2 -step_z 1’  as 2 meter by 2 meter by 1 meter. These parameters were found experimentally (see the README file for more on this).
lasnoise -i bugurunni_raw\*.laz ^
         -step_xy 2 -step_z 1 ^
         -isolated 15 ^
         -odir bugurunni_noise -olaz ^
         -cores 4
Then we run lasground on 4 cores to classify the ground points with options ‘-metro’ and ‘-bulge 1.0’. The option ‘-metro’ is a convenient short-hand for ‘-step 50’ that will remove all objects on the terrain (e.g. large buildings) that have an extend of 50 meters or less. The option ‘-bulge 1.0’ instructs lasground to be conservative and only add points that are 1 meter or less above a smoothed version of the initial ground estimate (see the README file for more on this)..
lasground -i bugurunni_noise\*.laz ^
          -ignore_class 7 ^
          -metro -bulge 1.0 ^
          -odir bugurunni_ground -olaz ^
          -cores 4
Now we use las2dem to raster a DTM from only those points that were classified as ground. The option ‘-step 0.5’ sets the output grid resolution to 0.5 meters, ‘-kill 200’ interpolates across voids of up to 200 meters, and ‘-use_tile_bb’ rasters only the original 500 meter by 500 meter tile interior but not the 50 meter buffer. This assures that the resulting raster tiling aligns without artifacts across tile boundaries. The option ‘-obil’ chooses BIL as the output raster format.
las2dem -i bugurunni_ground\*.laz ^
        -keep_class 2 ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir bugurunni_dtm -obil ^
        -cores 4
As a simply form of anti-aliasing we average each four pixels of 0.5 meter resolution into one pixel of 1.0 meter resolution with lasgrid as all LAStools can read BIL files via on-the-fly conversion to points.
lasgrid -i bugurunni_dtm\*.bil -merged ^
        -step 1.0 -average ^
        -o bugurunni_dtm.bil

Finally we create a hillshade of the DTM adding back the projection that was “lost” in the BIL file generation so that blast2dem – the extremely scalable BLAST version of las2dem – can automatically produce a KML file for display in Google Earth.

blast2dem -i bugurunni_dtm.bil ^
          -hillshade -utm 37M ^
          -o bugurunni_dtm_hill.png

For comparison we also create a DSM with the same three steps but using all points.

las2dem -i bugurunni_raw\*.laz ^
        -step 0.5 -kill 200 -use_tile_bb ^
        -odir bugurunni_dsm -obil ^
        -cores 4
lasgrid -i bugurunni_dsm\*.bil -merged ^
        -step 1.0 -average ^
        -o bugurunni_dsm.bil
blast2dem -i bugurunni_dsm.bil ^
          -hillshade -utm 37M ^
          -o bugurunni_dsm_hill.png
DTM of "Bugurunni" from ground points classified with LAStools

DTM of “Bugurunni” from ground points classified with LAStools

Above you see the generated DTM and below the corresponding DSM. So yes, LAStools can create DTMs from points that are result of dense-matching photogrammetry … under one assumption: there is not too much vegetation.

DSM of "Bugurunni" from all dense-matching points

DSM of “Bugurunni” from all dense-matching points

Below also the results for the “Borahatward” data. In a future blog post we will see how to deal with the excessive low noise sometimes present in dense-matching points.

DTM of "Bo"

DTM of “Borahatward” from ground points classified with LAStools

DSM of "Borahatward" from all dense-matching points

DSM of “Borahatward” from all dense-matching points

LAStools Toolbox for ERDAS IMAGINE 2014 Released

PRESS RELEASE (for immediate release)
April 8, 2015
rapidlasso GmbH, Gilching, Germany

The latest release of LAStools from rapidlasso GmbH contains a new toolbox for ERDAS IMAGINE®2014, allowing the users of Hexagon Geospatial’s powerful remote sensing software to utilize the popular rapidlasso LiDAR processing modules. The tools have been fully integrated into the software so that they are also available as operators within the IMAGINE Spatial Modeler framework.

This new ERDAS IMAGINE 2014 toolbox is included in the latest version of the LAStools distribution. This release instantly augments the existing image analysis tools of ERDAS IMAGINE 2014 with the widely-popular point processing capabilities of LAStools. These new tools empower users to validate, quality-check, clean, classify, thin, raster, contour, and compress LiDAR pointclouds in LAS or LAZ formats as well as directly consume and further analyze the resulting raster or vector products with the rich functionality of ERDAS IMAGINE.

Sample LAStools workflow created with the IMAGINE Spatial Modeler.

Sample LAStools workflow created with the Spatial Modeler of ERDAS IMAGINE.

About rapidlasso GmbH:
Technology powerhouse rapidlasso GmbH specializes in efficient LiDAR processing tools that are widely known for their high productivity. They combine robust algorithms with efficient I/O and clever memory management to achieve high throughput for data sets containing billions of points. The company’s flagship product – the LAStools software suite – has deep market penetration and is heavily used in industry, government agencies, research labs, and educational institutions. Visit http://rapidlasso.com for more information.

About Hexagon Geospatial:
Hexagon Geospatial helps you make sense of the dynamically changing world. Hexagon Geospatial provides the software products and platforms to a large variety of customers through direct sales, channel partners and other Hexagon businesses. For more information, visit http://hexagongeospatial.com.

Hexagon Geospatial is part of Hexagon, a leading global provider of information technologies that drive quality and productivity improvements across geospatial and industrial enterprise applications. Hexagon’s solutions integrate sensors, software, domain knowledge and customer workflows into intelligent information ecosystems that deliver actionable information, automate business processes and improve productivity. They are used in a broad range of vital industries. Hexagon (Nasdaq Stockholm: HEXA B) has more than 15,000 employees in 46 countries and net sales of approximately 3.1bn USD. Learn more at http://hexagon.com.

© 2015 Hexagon AB and/or its subsidiaries and affiliates. All rights reserved. Hexagon and the Hexagon logo are registered trademarks of Hexagon AB or its subsidiaries. All other trademarks or servicemarks used herein are property of their respective owners. Hexagon Geospatial believes the information in this publication is accurate as of its publication date. Such information is subject to change without notice.

Processing Large Dense-Matching Point Clouds

It is easy to generate Digital Surface Models (DSMs) from the point clouds generated by dense-matching photogrammetry. The blast2dem tool that is part of the BLAST extension of LAStools can seamlessly triangulate and raster up to 2 billion points. In an earlier blog article we showed how to also create Digital Terrain Models (DTMs) from photogrammetric points by combining lassort, lasground, and las2dem. We were since asked how to do this for larger inputs. Here the answer.

pix4d_large1_raw.png

In the following we describe a LAStools pipeline that turns the point cloud shown above into a DTM, a DSM, and 1 meter elevation contours using tile-based multi-core batch processing. Below you see the result.

macmahanfarm_dsm_dtm_hill_small

The ‘macmahanfarm.laz’ file has 140,616,508 RGB-colored points and was provided by Pix4D as a LAZ file of 1.2 GB that decompresses to 4.5 GB of LAS. The unusually low LASzip compression ratio of “only” 1:4 comes from the incoherent order of the points in the file and their excessive coordinate resolution.

First we use lastile to tile the points of ‘macmahanfarm.laz’ into 500 meter by 500 meter tiles with 25 buffer. This buffer will help to avoid edge artifacts along tile boundaries. A lasinfo report tells us that point coordinates are stored with millimeter resolution in the original file. Because the data is not that precise we rescale it on-the-fly to centimeters, which improves compression and eases I/O bottlenecks.

lastile -i macmahanfarm.laz ^
        -rescale 0.01 0.01 0.01 ^
        -tile_size 500 -buffer 25 ^
        -o poop\mmfarm.laz -olaz

One oddity in the data obtained from Pix4D is the high spatial incoherence in the ordering of points in the original file. Here an illustration of the point order in ‘mmfarm_266500_6175500.laz’, one of the 65 tiles generated in the last step.

Such incoherence in point order slows down subsequent processing with lasground and las2dem (which benefit from spatial locality). We use lassort to reorder the tiles on 4 cores in parallel and add appendix ‘_s’ to each file name:

lassort -i poop\mmfarm*.laz ^
        -odix _s -olaz ^
        -cores 4

This sort the points along a standard Hilbert or Morton space filling curve. Coherence in point order also improves compression: the sorted tile files are about 30 percent.smaller than the unsorted ones.

Next we use lasground to classify the 65 sorted tiles on 4 cores in parallel. A step size of 20 was chosen after neither 10 nor 15 were able to remove a few dense patches of vegetation. Alternatively, use a smaller step size and reclassify points manually with lasview.

lasground -i poop\mmfarm*_s.laz ^
          -step 20 -extra_fine ^
          -odix g -olaz ^
          -cores 4

Below a closer look at the result for tile ‘mmfarm_266500_6175000_sg.laz’. Note that our data is mostly open terrain. Your milage will vary (a lot!!!) depending on the amount of vegetation in the scene. There will generally not be any elevation data in vegetated areas. For larger forested patches you will need LiDAR to get a DTM. The same is true if you are particularly interested in bare-earth terrain features beneath trees and shrubs.

Then we use las2dem to raster the ground points of each tile as a DTM with 0.5 meter resolution in the BIL format, again on 4 cores in parallel. The ‘-use_tile_bb’ limits rasterization to the original 500 meter by 500 meters of the tile, clipping the 25 meter buffer along the tile boundaries from the output.

las2dem -i poop\mmfarm*_sg.laz ^
        -keep_class 2 ^
        -step 0.5 -use_tile_bb ^
        -ocut 3 -odix _dtm -obil ^
        -cores 4

And then merge them into a single DTM using lasgrid, which (like all other LAStools) can read raster data in BIL, ASC, and DTM format and convert it on-the-fly to points:

lasgrid -i poop\mmfarm*_dtm.bil -merged ^
        -step 0.5 ^
        -o macmahanfarm_dtm.tif

The DSM is created the same way but without filtering out points. Here a trick how to quickly generate contours that are not too complex (but also not of catographic quality). We first use lasgrid to merge and average the DTMs into a smoother and coarser DTM from which we extract 1 meter elevation contours with blast2iso.

lasgrid -i poop\mmfarm*_dtm.bil -merged ^
        -step 1.0 -average ^
        -o poop\dtm_coarse.bil
blast2iso -i poop\dtm_coarse.bil ^
          -iso_every 1.0 ^
          -smooth 2 -simplify_length 0.5 -simplify_area 0.5 -clean 15 ^
          -o macmahanfarm_iso.shp

pix4d_large1_qgis_dtm

You can find an example batch script that summarizes all the above steps into one workflow in the latest LAStools distribution.

pix4d_large1_qgis_dtm_contours

Finally, a word on compression: lowering the coordinate resolution to centimeters, sorting the points of every tile into a coherent spatial order, and merging the tiles back into one, results in a ‘macmahanfarm_cm.laz’ file that is only 0.56 GB: less than half the size of the ‘macmahanfarm.laz’ file we obtained from Pix4D. Now the LASzip compression ratio is 1:8.

Generating a DTM from Dense-Matching Photogrammetry

Point clouds from dense-matching photogrammetry are popular. Stored in the LAS format – or its compressed LAZ twin – these “near-LiDAR” points are easily ingested into LiDAR processing software to, for example, generate a Digital Surface Model (DSM). “Can LAStools also create a Digital Terrain Model (DTM) from such data?” is a question we often get asked by private email or via the LAStools user forum. At INTERGEO Euroasia 2014 we gave a demo of this at the UltraCam booth. We repeated this demo with different data at the Pix4D booth during the Geospatial World Forum 2014. Now the folks from Visionmap have also become curious …

pix4d_raw

Above the “topography_cadastre_switzerland_densified.las” point cloud we got at the Pix4D booth. You may have already seen this colorized point cloud at some other Pix4D event. The original points were in a highly incoherent spatial order that negativly affected all subsequent processing. Therefore we first use lassort to reorder the points into a space-filling Hilbert curve. At the same time we rescale the coordinate resolution to centimeters and compress the file into the LAZ format.

lassort -i topography_cadastre_switzerland_densified.las ^
        -rescale 0.01 0.01 0.01 -olaz
Next we use lasground to classify the ground points with options ‘-city’ and ‘-ultra_fine’. The option ‘-city’ is a convenient short-hand for ‘-step 25’ that will remove all objects on the terrain that have an extend of 25 meters or less. The option ‘-ultra_fine’ instructs lasground to spend more time on finding a good initial ground estimate. The option ‘-odix _g’ add the appendix ‘_g’ to the output file name.
lasground -i topography_cadastre_switzerland_densified.laz ^
          -city -ultra_fine ^
          -odix _g -olaz
pix4d_ground
Above you see only the points that were classified as ground (classification 2) and below you see how triangulating these points interpolates the voids where buildings and vegetation used to be.
pix4d_triangulated
Now we use las2dem to raster a DTM from only those points that were classified as ground. The option ‘-step 0.5’ sets the output grid resolution to 0.5 meters, ‘-ocut 2’ removes the ‘_g’ from the file name, ‘-odix _dtm’ adds the appendix ‘_dtm’ to the file name, and ‘-obil’ chooses BIL as the output raster format.
las2dem -i topography_cadastre_switzerland_densified_g.laz ^
        -keep_class 2 -thin_with_grid 0.125 -extra_pass ^
        -step 0.5 ^
        -ocut 2 -odix _dtm -obil
For comparison we also use las2dem to create the corresponding DSM. The option ‘-thin_with_grid 0.125’ makes the density of the points that are triangulated into a TIN more “compatible” with the 0.5 by 0.5 grid that the TIN is rastered with by keeping only the first point that falls within each 0.125 by 0.125 area or a maximum of 16 points per 0.5 by 0.5 raster cell. The option ‘-extra_pass’ counts all surviving points in a first pass to minimize the memory footprint needed for Delaunay TIN construction.
las2dem -i topography_cadastre_switzerland_densified_g.laz ^
        -thin_with_grid 0.125 -extra_pass ^
        -step 0.5 ^
        -ocut 2 -odix _dsm -obil

pix4d_dtm

Did you know that lasview can visualize BIL files via on-the-fly conversion from grids to points? Above you see the generated DTM and below the corresponding DSM. So yes, LAStools can create DTMs from points that are result of dense-matching photogrammetry … under one assumption: there is not too much vegetation.

pix4d_dsm