Return to Revisions

I tried naively to use nearest neighbors to separate the ground and the trees. I iteratively set a point to be a tree, ground, or other type of point. A point is a tree or ground point if a percentage of the k=10 nearest neighbors are tree or ground points. The percentage decreases with each iteration. The initial guess of trees are points with intensity < 600 (I also tried points where the convex hull of the nearest neighbors has volume > 1e6, but this produces worse results).

import numpy as np
import pylas
from sklearn.neighbors import NearestNeighbors
import sys


def summary(i, istree):
    print("{0}: trees: {1}, ground: {2}, other: {3}"
          .format(i,
                  (istree == 1).sum() / len(istree),
                  (istree == -1).sum() / len(istree),
                  (istree == 0).sum() / len(istree)))


def main(fname1, fname2):
    with pylas.open(fname1) as f:
        outf = f.read()
    pts = outf.points
    data = np.asarray([pts[i] for i in ("X",
                                        "Y",
                                        "Z")]).transpose().copy()
    coords = data[:, :3]

    Nnbs = 10
    nbrs = NearestNeighbors(n_neighbors=Nnbs,
                            algorithm='ball_tree').fit(coords)
    distances, indices = nbrs.kneighbors(coords)
        
    tree_pts = np.logical_or(pts["intensity"] < 600, False)
    ground_pts = True ^ tree_pts
    
    istree = tree_pts * 1 + -1 * ground_pts
    
    for i in range(3):
        istree = solve(istree, indices,
                       bound=0.2)
        summary(i, istree)

    niter = 10
    for i in range(niter):
        istree = solve(istree, indices,
                       bound=0.2 * (niter - 1 - i) / (niter - 1))
        summary(i, istree)
            
    classification = np.zeros(len(istree), dtype=np.uint8)
    classification += np.uint8(3) * (istree == 1)
    classification += np.uint8(2) * (istree == -1)
    classification += np.uint8(1) * (istree == 0)
    
    outf.points["raw_classification"] = classification
    outf.write(fname2)

    
def solve(istree, indices, bound):
    n_istree = np.zeros(shape=(istree.shape[0],), dtype=np.int8)
    
    for pt, inds in enumerate(indices):
        s = istree[inds].sum() / len(inds)
        if s > bound:
            n_istree[pt] = 1
        elif s < -bound:
            n_istree[pt] = -1
        else:
            n_istree[pt] = 0
            
    return n_istree


if __name__ == "__main__":
    main(sys.argv[1], sys.argv[2])

Unfortunatley this code produces worse results than pdal.

The .las file for testing is available here.