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])

Unfortunately this code produces worse results than pdal.

The .las file for testing is available here.

Edit:
Combining the above code with the individual tree segmentation of lidR gives good results.

require(lidR)

las = readLAS("test.las")
las = lasground(las, csf())
algo = pitfree(thresholds = c(0,10,20,30,40,50), subcircle = 0.2)
chm  = grid_canopy(las, 0.5, algo)
algo = watershed(chm, th = 4)
las  = lastrees(las, algo)
writeLAS(las, "tmp2.las")

For each point classified to be part of a tree by lidR I vote if this is a tree or ground by taking a weighted average of the result of the nearest neighbor python code. This is done in the following python code (tmp.las is the output of the nearest neighbor code).

import numpy as np
import pylas
from collections import defaultdict

with pylas.open("tmp.las") as f:
    x = f.read()

with pylas.open("tmp2.las") as f:
    y = f.read()

xx = np.array([x[v] for v in ("X", "Y", "Z", "intensity")]).transpose().copy()
yy = np.array([y[v] for v in ("X", "Y", "Z", "intensity")]).transpose().copy()

valid_points = []
ii = 0
for i in range(len(x.points)):
    if np.count_nonzero(xx[i] - yy[ii]) == 0:
        valid_points.append(x.points[i])
        ii += 1
        if ii > len(yy): break

assert len(valid_points) == len(yy)
        
x.points = np.array(valid_points, dtype=x.points.dtype)
        

treeID = y["treeID"]
classification = x["classification"]
ground = y["classification"] == 2

istree = (-0.25 * (classification == 2) +
          1 * (classification == 3) +
          0 * (classification == 1))


d = defaultdict(int)
for i in range(len(x.points)):
    d[treeID[i]] += istree[i]

for i in range(len(y.points)):
    if d[treeID[i]] > 0 and not ground[i]:
        classification[i] = 3
    else:
        classification[i] = 2

y["classification"] = classification
y.write("tmp3.las")

Still there are some trees which are classified as ground. I suspect the reason is that the treeID of lidR is an 8 bit integer and there seems to be more than 256 trees.