Tree stem detection

Published: December 7, 2018 / Last updated: February 24, 2020

Tested on Matlab r2019b, GNU Octave 5.2.0

Introduction

In this tutorial, you will learn how to detect tree stems in a 3D point cloud. Stem detection is an alternative approach to tree top detection and is likely to perform better on irregularly shaped canopies (in particular, if they were surveyed in leaf-off conditions).

We will use Airborne Laser Scanning (ALS) data acquired over the state of Geneva (Switzerland) in February 2017 (leaf-off conditions) with the Riegl LMS-Q1560 sensor. The data is openly available from Geneva's Official Geodata Portal.

A small forest located about 10 km to the North of Geneva city (see map below) will serve as our study case. Within this area, the main tree species are pedunculate oak (Quercus robur) and Scots pine (Pinus sylvestris).

The stem detection method we will apply is based on a simple observation: the space in a tree gets more cluttered as you approach the stem. Consequently, considering a homogeneous point sampling of the tree geometry, the point density (i.e. the number of points within a vertical column) tends to be higher near the stem (see figure 1 left). To reduce the effect of inhomogenous sampling and to simplify the detection of density peaks, the point cloud can be rasterized (see figure 1 right). The locations of local density peaks may subsequently serve as a basis for tree crown delineation, as proposed in Rahman and Gorte (2008) or more recently in Ayrey et al. (2017).

Figure 1 - The point density is higher near the stem, for a tree growing vertically. Left: Point counts. Right: Raster cell counts.

Setup

  1. Download and uncompress the Digital Forestry Toolbox (DFT) Zip or Tar archive
  2. Download the ge_2017_a.laz file from the DFT Zenodo repository and uncompress it with LASzip
  3. Start Matlab/Octave
  4. Delete any previous versions of the toolbox
  5. Add the DFT folders to the Matlab/Octave search path using addpath(genpath('path to DFT main folder'))
  6. Open dft_tutorial_3.m (located in the tutorials folder)

Step 1 - Reading the LAS file

We start by reading the LAS file using LASread():

pc = LASread('ge_2017_a.las');

Note that you may have to adjust the file path in the code above, depending on where you are storing the LAS file. Also note that the point classification uses the following scheme:

Class Description
0 Created, never lassified
1 Unclassified
2 Terrain
3 Low vegetation (< 50 cm)
5 High vegetation (≥ 50 cm)
6 Buildings
7 Outliers and incorrect measurements (Low noise)
9 Water
13 Bridges
15 Terrain (additional points)
16 Outliers and incorrect measurements (High noise)
19 Points measured outside of the area of interest (unclassified)

Step 2 - Normalize the point cloud elevation

The stem detection algorithm uses the planimetric coordinates and height of the points above ground as an input. To compute this height, we start by building a 0.5 m resolution raster terrain model from the classified 3D point cloud using elevationModels():

cellSize = 0.5;
[models, refmat] = elevationModels([pc.record.x, pc.record.y, pc.record.z], ...
    pc.record.classification, ...
    'classTerrain', [2, 15], ...
    'classSurface', [4,5], ...
    'cellSize', cellSize, ...
    'closing', 5, ...
    'interpolation', 'idw', ...
    'searchRadius', inf, ...
    'weightFunction', @(d) d^-3, ...
    'smoothingFilter', fspecial('gaussian', [2, 2], 0.8), ...
    'outputModels', {'terrain'}, ...
    'fig', false, ...
    'verbose', true);

We then subtract the terrain elevation from the point cloud elevation, to obtain the height of the points above the ground:

P = round([pc.record.x - refmat(3,1), pc.record.y - refmat(3,2)] / refmat(1:2,:));
ind = sub2ind([nrows, ncols], P(:,1), P(:,2));
xyh = [pc.record.x, pc.record.y, pc.record.z - models.terrain.values(ind)];

Step 3 - Filter the points by classification and return number

When compared to small branches, stems and primary branches are more likely to fully intercept the laser beam and generate a last return (see figure 2). Since we are only interested in detecting stems, we can apply a filter which keeps only last returns (and also removes terrain points or other unwanted classes). To do this, we create a logical (boolean) vector: 

idxl_last = pc.record.return_number == pc.record.number_of_returns;
idxl_veg = ismember(pc.record.classification, [4, 5]);
idxl_filter = idxl_veg & idxl_last;

Note that while using only the last returns may increase the detection reliability for large stems, it may also reduce the detection rate for smaller stems. You can modify this filtering criteria or add new ones (e.g. acquisition date, intensity, colour) to suit your specific needs.

Figure 2 - Cross section of the normalized point cloud. Top: all points. Bottom: only last returns and vegetation classes.

Step 4 - Detect the stems

Using only the filtered points, we can now detect the stems with the treeStems() function:

[label, xyh_stem] = treeStems(xyh, ...
    idxl_filter, ...
    'cellSize', 0.4, ...
    'bandWidth', 1.5, ...
    'verticalStep', 0.25, ...
    'searchRadius', 2, ...
    'minLength', 5, ...
    'verbose', true, ...
    'fig', true);

The parameters have the following meaning:

By modifying the parameters values, you can adjust the tradeoff between the sensitivity (recall) and the reliability (precision) of the detection. The interactive map below illustrates the detected stem positions.

Step 5 - Export the stem attributes to a CSV file

You can write the stem attributes to a Comma Separated Values text file (.csv) with:

fid = fopen('ge_2017_a_stems.csv', 'w+'); % open file
fprintf(fid, 'X, Y, H\n'); % write header line
fprintf(fid, '%.2f, %.2f, %.2f\n', xyh_stem'); % write records
fclose(fid); % close file

After exporting the CSV file, you can try opening it in any text editor (e.g. Notepad++).

Step 6 - Export the stem attributes to an ESRI shapefile


Matlab users, the shapewrite() function included here is currently not compatible with Matlab. Please use the shapewrite() function from the offical Matlab mapping toolbox instead.

To export the stem attributes to an ESRI shapefile, you first have to format them into a non-scalar structure and add a "Geometry" field:

S = struct('Geometry', repmat({'Point'}, size(xyh_stem,1),1), ...
      'X', num2cell(xyh_stem(:,1)), ...
      'Y', num2cell(xyh_stem(:,2)), ...
      'BoundingBox', [], ...
      'H', num2cell(xyh_stem(:,3)));

Specify the output filepath (by default, the file is created in the current working directory) and use shapewrite() to write the file:

shapewrite(S, 'ge_2017_a_stems.shp');

After exporting the shapefile, you can try opening it in any GIS software (e.g. Quantum GIS).

References