# Core LiDAR processing package
if (!requireNamespace("lidR", quietly = TRUE)) install.packages("lidR")

# Optional ground classification algorithm support
if (!requireNamespace("RCSF", quietly = TRUE)) install.packages("RCSF")

# Spatial handling and raster processing
if (!requireNamespace("sf", quietly = TRUE)) install.packages("sf")
if (!requireNamespace("terra", quietly = TRUE)) install.packages("terra")
if (!requireNamespace("stars", quietly = TRUE)) install.packages("stars")
if (!requireNamespace("alphashape3d", quietly = TRUE)) install.packages("alphashape3d")
if (!requireNamespace("knitr", quietly = TRUE)) install.packages("knitr")

# Load libraries
library(lidR)
library(sf)
library(terra)
library(stars)
library(alphashape3d)
library(knitr)

# Load the pre-classified LAS/LAZ file
las_path <- "./classified_3dsegformer.laz"
classified <- readLAS(las_path)

# Verify CRS and structure
print(classified)

# Check data integrity and attributes
las_check(classified) # We have duplicates since the tiles of laz files overlap to a certain degree (flight paths of the plane with the lidar sensor onboard overlap). In these overlapping zones, you can have dubble the amount of points. You can remove them before analysis (for example, remove the overlapping zones), but most functions inside the package will detect them automatically and remove them (but for for statistical analysis/metrics, it is indeed better to remove duplicates). The degenerated ground points can be caused by small errors in the sensor, or by "moving" objects on the ground during lidar scanning (For example, if there is wind, tree crowns move slightly, causing different z-values for the same locations).

# Visualization colored by elevation
plot(classified, color = "Z")

# Color by intensity
#plot(classified, color = "Intensity")

# Color by AI classification values
plot(classified, color = "Classification")

# Color by RGB (if available in point cloud attributes)
plot(classified, color = "RGB")

# Height normalization using a Triangular Irregular Network (TIN) ground model
norm <- normalize_height(classified, tin())

# Define transect endpoints (Coordinates in local projection, e.g., Belge 1972 / Lambert 72)
X1 <- 149800; Y1 <- 180514
X2 <- 150000; Y2 <- 180514

# Slice 2-meter wide transects for raw and normalized point clouds
slice <- clip_transect(classified, p1 = c(X1, Y1), p2 = c(X2, Y2), width = 2)
slice_norm <- clip_transect(norm, p1 = c(X1, Y1), p2 = c(X2, Y2), width = 2) 

# Compare raw vs. normalized profile plots
par(mfrow = c(2, 1))
plot(slice$X, slice$Z, col = slice$Z, pch = 19, main = "Raw Elevation Profile", xlab = "Distance", ylab = "Elevation (m)")
plot(slice_norm$X, slice_norm$Z, col = slice_norm$Z, pch = 19, main = "Normalized Height Profile", xlab = "Distance", ylab = "Height above ground (m)")

par(mfrow = c(1, 1))

# Filter normalized point cloud to high vegetation (ASPRS class 5)
high_veg <- filter_poi(norm, Classification == 5)

# Generate a 0.5m Canopy Height Model using Pit-free or Point-to-Raster algorithms
chm <- rasterize_canopy(high_veg, res = 0.5, algorithm = p2r(subcircle = 0.15))

# Locate individual tree tops using a Local Maximum Filter (LMF)
# Variable window size (ws) can also be used: e.g., function(z) { z * 0.1 + 1 }
ttops <- locate_trees(high_veg, lmf(ws = 5)) # This algorithm work both on normalized point clouds, as dsm (but results can be slightly different).

# Plot CHM with detected tree tops
plot(chm, main = "Canopy Height Model with Detected Treetops")
plot(sf::st_geometry(ttops), add = TRUE, col = "red", pch = 3, cex = 0.8)

# Generate a 0.5m Canopy Height Model using Pit-free or Point-to-Raster algorithms
chm2 <- rasterize_canopy(norm, res = 0.5, algorithm = p2r(subcircle = 0.15))

# Locate individual tree tops using a Local Maximum Filter (LMF)
# Variable window size (ws) can also be used: e.g., function(z) { z * 0.1 + 1 }
ttops2 <- locate_trees(norm, lmf(ws = 5))

# Plot CHM with detected tree tops
plot(chm2, main = "Canopy Height Model with Detected Treetops")
plot(sf::st_geometry(ttops2), add = TRUE, col = "red", pch = 3, cex = 0.8)

nrow(ttops)
nrow(ttops2)

# Point-to-point tree segmentation using Li et al. (2012)
segmented_las <- segment_trees(high_veg, li2012()) # +/- 3 min. computation time

# Plot the segmented point cloud where colors represent individual tree IDs
plot(segmented_las, color = "treeID")

ttops_seg <- tree_metrics(segmented_las, func = .stdmetrics) %>% st_as_sf()
kable(summary(ttops_seg[,c( "zmax", "zmean","zsd","zq95","area")]%>% st_drop_geometry()))


hist(ttops_seg$zmax,plot = T, main = "Height Distribrution (m)",xlab =   "Height (m)",ylab = "Frequency")

# Clip transect of segmented trees
slice_tree <- clip_transect(segmented_las, p1 = c(X1, Y1), p2 = c(X2, Y2), width = 2)

# Plot profile colored by segmented tree IDs
plot(slice_tree$X, slice_tree$Z, col = slice_tree$treeID, pch = 19, 
     xlab = "Distance", ylab = "Height", main = "Individual Tree Segmentation Cross-Section")

# Extract a single tree
tree_1 <- filter_poi(segmented_las, treeID == 1)
plot(tree_1)

# Extract maximum tree height directly
max_height <- max(tree_1$Z)
cat("Tree 1 Max Height:", max_height, "m\n")

# Standard crown metrics
crown_m <- crown_metrics(tree_1, .stdtreemetrics)
summary(crown_m)

# Standard height metrics
cloud_m <- cloud_metrics(tree_1, .stdmetrics_z)
print(cloud_m)

get_crown_attributes <- function(X, Y, Z, alpha = c(Inf, 1)) {
  if (length(X) <= 3 || length(Y) <= 3 || length(Z) <= 3) {
    message("Cannot compute a 3D hull from 3 or fewer points")
    return(NULL)
  }
  
  # Prepare centered matrix for stable alpha shape computation
  mean_x <- mean(X)
  mean_y <- mean(Y)
  
  a3d_centered <- cbind(X - mean_x, Y - mean_y, Z)
  
  top_x <- X[which.max(Z)]
  top_y <- Y[which.max(Z)]
  
  # Compute alpha shape
  ashape <- alphashape3d::ashape3d(x = a3d_centered, alpha = alpha, pert = TRUE)
  
  # Calculate crown metrics
  df <- data.frame(
    Zmax = max(Z),
    Zq999 = as.numeric(quantile(Z, 0.999)),
    Zq99 = as.numeric(quantile(Z, 0.990)),
    Z_mean = mean(Z),
    n_points = length(Z),
    vol_convex = alphashape3d::volume_ashape3d(ashape, indexAlpha = 1),
    vol_concave = alphashape3d::volume_ashape3d(ashape, indexAlpha = 2),
    CV_Z = sd(Z) / mean(Z),
    CRR = (mean(Z) - min(Z)) / (max(Z) - min(Z)),
    X = top_x,
    Y = top_y
  )
  
  # Return both metrics and the fitted hull object along with shift offsets
  return(list(
    metrics = df, 
    ashape = ashape, 
    centered_points = a3d_centered
  ))
}

# Run function
result <- get_crown_attributes(
  tree_1@data[["X"]],
  tree_1@data[["Y"]],
  tree_1@data[["Z"]],
  alpha = c(Inf, 1)
)

# Extract components
metrics <- result$metrics
ashape_obj <- result$ashape
pts_centered <- result$centered_points

# --- Plot 1: Convex Hull (alpha = Inf) ---
plot(ashape_obj, indexAlpha = 1, transparency = 0.3, axes = TRUE)
rgl::points3d(pts_centered, color = "black", size = 4)

# --- Plot 2: Concave Hull (alpha = 1) ---
plot(ashape_obj, indexAlpha = 2, transparency = 0.4, axes = TRUE)
rgl::points3d(pts_centered, color = "black", size = 4)

# Using example dataset for multi-tree crown geometries
LASfile <- system.file("extdata", "MixedConifer.laz", package = "lidR")
trees <- readLAS(LASfile, filter = "-drop_z_below 0")

# Convex Hull
metrics_convex <- crown_metrics(trees, .stdtreemetrics, geom = "convex")

# Concave Hull
metrics_concave <- crown_metrics(trees, .stdtreemetrics, geom = "concave")

# Bounding Box
metrics_bbox <- crown_metrics(trees, .stdtreemetrics, geom = "bbox")

# Visualize concave hulls
par(mfrow=c(1,3))
plot(sf::st_geometry(metrics_concave))
plot(sf::st_geometry(metrics_convex))
plot(sf::st_geometry(metrics_bbox))

par(mfrow=c(1,1))     

# Define custom height & intensity metric function
myMetrics <- function(z, i) {
  metrics <- list(
    zwimean = sum(z * i) / sum(i), # Intensity-weighted mean elevation
    zimean  = mean(z * i),         # Mean product of height and intensity
    zsqmean = sqrt(mean(z^2))      # Quadratic mean height (RMS)
  )
  return(metrics)
}

# Apply custom metrics over a stars spatial grid template
template_stars <- stars::st_as_stars(st_bbox(segmented_las), dx = 10, dy = 10)
m_custom <- template_metrics(segmented_las, myMetrics(Z, Intensity), template_stars)
plot(m_custom)

plot(chm)

# 1. Basic approach
first_returns <- filter_poi(segmented_las, ReturnNumber == 1)
m_pixel_basic <- pixel_metrics(first_returns, ~mean(Z), res = 20)

# 2. Optimized approach (Inline filtering - memory efficient)
m_pixel_opt <- pixel_metrics(segmented_las, ~mean(Z), res = 20, filter = ~ReturnNumber == 1)

# 3. Echo-type breakdown
echo_types <- c("all", "first", "multiple")
func_echo <- ~list(avgI = mean(Intensity))
m_echo <- pixel_metrics(segmented_las, func_echo, res = 20, by_echo = echo_types)
plot(m_echo)

# Load sample forest inventory plot data
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile, filter = "-keep_random_fraction 0.5")
plot(las)
shpfile <- system.file("extdata", "efi_plot.shp", package = "lidR")
inventory <- sf::st_read(shpfile, quiet = TRUE)
inventory # contains an ID and a Value Of Interest (VOI) per plot
plot(inventory[,1]) # Point shapefile dataset.

# Calculate metrics per circular plot (radius = 11.28m, ~400 m² plot)
fun_q85 <- ~list(q85 = quantile(Z, probs = 0.85))
m_plot <- plot_metrics(las, fun_q85, inventory, radius = 11.28)
plot(header(las))
plot(m_plot["q85"], pch = 19, cex = 3, add = TRUE)

# Buffer inventory points into explicit polygon objects
inventory_poly <- sf::st_buffer(inventory, 11.28)
m_poly <- plot_metrics(las, .stdmetrics_z, inventory_poly)

# Works with polygons as well
inventory <- sf::st_buffer(inventory, 11.28)
plot(header(las))
plot(sf::st_geometry(inventory), add = TRUE)
m <- plot_metrics(las, .stdmetrics_z, inventory)
plot(m["zq85"], pch = 19, cex = 3, add = TRUE)

# Generate 8m voxel counts and mean intensity
v_count <- voxel_metrics(las, length(Z), res = 8)
v_intensity <- voxel_metrics(las, mean(Intensity), res = 8)

# 3D Voxel rendering
plot(v_intensity, color = "V1", colorPalette = heat.colors(50), trim = 60, voxel = TRUE)

# 3D Voxel rendering
plot(v_count, color = "V1", colorPalette = heat.colors(50), trim = 60, voxel = TRUE)