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.

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