Individual Tree Segmentation segmenting each tree into many pieces

I am a first time user trying to help my community inventory trees over a threshold height using public Lidar data. My goal is to extract hulls and/or approximate GPS coordinates of each tree exceeding the target height.

I followed along with the individual tree segmentation examples (1,2) but in my results trees are generally segmented into 4-8 parts based on the colorization in the image (included).

I first thought this had to do with the watershed parameters, but I didn’t hit on the right combo in the workflow. Later I tried the li2012 algorithm and my results were no better.

Can someone suggest either the key parameters I should be changing, or a different workflow I should be following?

print(Sys.time())
library(lidR)

las = readLAS("/Users/brianvanvoorst/Desktop/USGS_LPC_MI_GrandTraverseCO_2015_380522_LAS_2017.las")
las = classify_ground(las, csf())
print ("Classify done")
las = normalize_height(las, tin())
print ("Normalize done")
algo = pitfree(thresholds = c(0,10,20,30,40,50), subcircle = 0.2)
print ("Pitfree done")
chm  = grid_canopy(las, 0.5, algo)

plot(chm, col = height.colors(50))
# smoothing post-process (e.g. two pass, 3x3 median convolution)
ker = matrix(1,3,3)
chm = focal(chm, w = ker, fun = median)
chm = focal(chm, w = ker, fun = median)

plot(chm, col = height.colors(50)) # check the image
algo = watershed(chm, th = 4)
las  = segment_trees(las, algo)

# remove points that are not assigned to a tree
trees = filter_poi(las, !is.na(treeID))

plot(trees, color = "treeID", colorPalette = pastel.colors(100))
print(Sys.time())

---

Okay here is my latest source and result

require(lidR)
require(rlas)
require(rgdal)
require(tictoc)

las <- readLAS("/Users/brianvanvoorst/Desktop/USGS_LPC_MI_GrandTraverseCO_2015_380522_LAS_2017.las", filter="-keep_class 1L")

#dtm <- grid_terrain(las, algorithm = knnidw(k = 8, p = 2))
# Error: No ground points found. Impossible to compute a DTM.
#las_normalized <- normalize_height(las, dtm)

# Create a filter to remove points above 95th percentile of height
filter_noise = function(las, sensitivity)
{
  p95 <- grid_metrics(las, ~quantile(Z, probs = 0.95), 10)
  las <- merge_spatial(las, p95, "p95")
  las <- filter_poi(las, Z < p95*sensitivity)
  las$p95 <- NULL
  return(las)
}

las_denoised <- filter_noise(las, sensitivity = 1.2)

chm <- grid_canopy(las_denoised, 0.5, pitfree(c(0,2,5,10,15), c(3,1.5), subcircle = 0.2))

plot_dtm3d(chm)

ker <- matrix(1,5,5)
chm_s <- focal(chm, w = ker, fun = median)

algo <- watershed(chm_s, th = 4)
las_watershed  <- segment_trees(las_denoised, algo)

# remove points that are not assigned to a tree
trees <- filter_poi(las_watershed, !is.na(treeID))

# View the results
plot(trees, color = "treeID", colorPalette = pastel.colors(100))

las = readLAS("/USGS_LPC_MI_GrandTraverseCO_2015_380522_LAS_2017.laz", filter = "-drop_y_above 523580 -drop_x_below 19381500")
las = normalize_height(las, tin())
thresholds = round(c(0,5,10,15,20,25)/0.3048,0) # The highest point is ~80 feet ~= 25 m
algo = pitfree(thresholds = thresholds, max_edge = c(0, 2), subcircle = 0.2/0.3048)
chm = grid_canopy(las, 2, algo)
plot(chm, col = height.colors(50))