TITLE: Processing Terrestrial LiDAR with PDAL
DATE: 2020-12-15
AUTHOR: John L. Godlee
====================================================================


I used a Leica HDS6100 terrestrial laser scanner for my PhD work. 
The scans are processed initially in a program produced by Leica 
called Cyclone. I used Cyclone to georeference and stitch together 
scans from a single plot to create a single shadowless scene. From 
there they can be exported to .ptx files, each of which can contain 
multiple scans, each with their own local coordinate reference 
system. Below is an annotated .ptx file, showing the role of each 
line:

  [Cyclone]: 
https://leica-geosystems.com/en-gb/products/laser-scanners/software/
leica-cyclone

    20224  # Number of columns
    8615  # Number of rows
    482595.121831 8330769.987967 1254.138086  # Scanner registered 
position in real space (xyz)
    -0.990870 -0.134818 -0.000312  # Scanner registered axis 'X'
    0.134818 -0.990870 -0.000175  # Scanner registered axis 'Y'
    -0.000285 -0.000215 1.000000  # Scanner registered axis 'Z'
    -0.990870 -0.134818 -0.000312 0  # 4x4 tranformation matrix
    0.134818 -0.990870 -0.000175 0 
    -0.000285 -0.000215 1.000000 0
    482595.121831 8330769.987967 1254.138086 1
    0 0 0 0.500000  # Start of point coordinates
    0 0 0 0.500000  # Unreturned pulses
    0 0 0 0.500000  
    0 0 0 0.500000
    -0.000046 0.909775 -1.885635 0.010376  # First returned pulse
    -0.000046 0.903366 -1.870834 0.015015
    -0.000046 0.895859 -1.853836 0.019165
    -0.000046 0.894424 -1.849380 0.020874
    -0.000046 0.898849 -1.857010 0.024781

When a new scan starts the same header material will be repeated, 
but with different values depending on the scanner position.

My end goal was to have a voxelised point cloud with noise removed 
so that only foliage material remains.

First I needed to split the .ptx file into separate scans, based on 
the header material in each scan:

    #!/usr/bin/env sh

    if [ $# -ne 1 ]; then
        printf "Must supply one argument:\n  [1] input.ptx\n"
        exit 1
    fi

    # Get lines at which to split 
    lines=$(rg -n --no-encoding -M 10 "^[0-9]+\s+?$" $1 | 
        sed 's/:.*//g' | 
        awk 'NR%2!=0' | 
        tr '\n' ' ' | 
        sed 's/^[0-9]\s//g')

    # Get name of file without extension
    noext="${1%.ptx}"

    # Split file by scans using array dimension rows in header as 
line ref
    csplit -f $noext -b "_%d.ptx" $1 $lines

Then I converted each scan into a .laz file, with the coordinates 
transformed according to the transformation matrix in the header 
material:

    #!/usr/bin/env sh

    if [ $# -lt 2 ]; then
        printf "Must supply at least two arguments:\n  [1] 
input.ptx\n  [2] output.laz\n"
        exit 1
    fi

    # For each argument
    matrix=$(head -n 10 $1 | tail -4 | sed -r 's/0\s+?$/0.0/g' | 
awk -f transpose.awk)

    pdal pipeline pipelines/ptx_laz.json --readers.text.filename=$1 
\
        --filters.transformation.matrix="${matrix}" \
        --writers.las.filename=$2

Here is the PDAL pipeline:

    [
        {
            "type" : "readers.text",
            "filename" : "input.txt",
            "header" : "X Y Z I",
            "skip" : 10
        },
        {
            "type" : "filters.transformation",
            "matrix" : "0 -1  0  1  1  0  0  2  0  0  1  3  0  0  0 
 1"
        },
        {
            "type" : "writers.las",
            "compression" : "true",
            "minor_version" : "2",
            "dataformat_id" : "0",
            "forward" : "all",
            "filename" : "output.laz"
        }
    ]

Then I merged the .laz files back together using pdal merge 
file1.ptx file2.ptx ....

Then I voxelised the .laz file, again using a PDAL pipeline. This 
method creates voxels of 0.01 m^3 (1 cm^3), with a point at the 
center of each occupied voxel:

    [
        {
            "type" : "readers.las",
            "filename" : "input.laz"
        },
        {
            "type":"filters.voxeldownsize",
            "cell" : 0.01,
            "mode" : "center"
        },
        {
            "type" : "writers.las",
            "compression" : "true",
            "minor_version" : "2",
            "dataformat_id" : "0",
            "forward" : "all",
            "filename" : "output.laz"
        }
    ]

Finally I excluded noise with a PDAL pipeline. This method measures 
the mean distance of each point to its eight nearest neighbours, 
then excludes points with mean distances greater than the 95% 
confidence interval of the distribution:

    [
        {
            "type" : "readers.las",
            "filename" : "input.laz"
        },
        {
            "type" : "filters.outlier",
            "method" : "statistical",
            "mean_k" : 8,
            "multiplier" : 1.96
        },
        {
          "type" : "filters.range",
          "limits" : "Classification![7:7]"
        },
        {
            "type" : "writers.las",
            "compression" : "true",
            "minor_version" : "2",
            "dataformat_id" : "0",
            "forward" : "all",
            "filename" : "output.laz"
        }
    ]