Tree_segmentation function integrated to tutorial, modifications based on sf package

R/tree.detection.Rmd

@@ -51,7 +51,7 @@ Otherwise you can load your own data provided positions and heights are measured
 ```{r prepareTreeInventory, eval=FALSE}
 # import field inventory
 fichier <- "chablais3_listeR.csv"
-tree.inventory <- read.csv(file = fichier, sep = ";", header = F, stringsAsFactors = TRUE)
+tree_inventory_chablais3 <- read.csv(file = fichier, sep = ";", header = F, stringsAsFactors = TRUE)
 names(tree_inventory_chablais3) <- c("x", "y", "d", "h", "n", "s", "e", "t")
 # save as rda for later access
 # save(tree_inventory_chablais3,file="tree_inventory_chablais3.rda")
@@ -95,13 +95,17 @@ ggplot(tree_inventory_chablais3, aes(x = x, y = y, group = s)) +
   labs(color = "Species") # titre de la légende
 ```
 
-We define the region of interest (ROI) to crop ALS data to corresponding extent before further processing. ROI is set on the extent of tree inventory, plus a 10 meter buffer.
+We define the region of interest (ROI) to crop ALS data to corresponding extent before further processing. ROI is set on the extent of tree inventory, plus a 10 meter buffer. The tree table is converted to a `sf` object to make spatial processing easier in the following steps.
 
 ```{r roi, include = TRUE}
+# duplicate coordinates to ensure they remain in the data.frame
+tree_inventory_chablais3[, c("X", "Y")] <- tree_inventory_chablais3[, c("x", "y")]
+# convert to spatial sf object
+tree_inventory_chablais3 <- sf::st_as_sf(tree_inventory_chablais3, coords = c("X", "Y"), crs = 2154)
 # buffer to apply around ROI (meters)
 ROI_buffer <- 10
-# ROI limits
-ROI_range <- data.frame(round(apply(tree_inventory_chablais3[, c("x", "y")], 2, range)))
+# ROI limits: bounding box of trees
+ROI_range <- round(sf::st_bbox(tree_inventory_chablais3))
 ```
 
 ### Airborne Laser Scanning data
@@ -126,13 +130,10 @@ lazdir <- "../data/tree.detection"
 cata <- lidR::readALSLAScatalog(lazdir)
 # set coordinate system
 lidR::projection(cata) <- 2154
-# extract points in ROI plus additional buffer
-las_chablais3 <- lidR::clip_rectangle(
+# extract points in ROI plus additional 5m buffer
+las_chablais3 <- lidR::clip_roi(
   cata,
-  ROI_range$x[1] - ROI_buffer - 5,
-  ROI_range$y[1] - ROI_buffer - 5,
-  ROI_range$x[2] + ROI_buffer + 5,
-  ROI_range$y[2] + ROI_buffer + 5
+  ROI_range + (ROI_buffer + 5) * c(-1, -1, 1, 1)
 )
 # save as rda for easier access:
 # save(las_chablais3, file="las_chablais3.rda", compress = "bzip2")
@@ -141,15 +142,15 @@ las_chablais3 <- lidR::clip_rectangle(
 ## Data preparation
 ### Digital Elevation Models
 
-From the ALS point cloud, digital elevation models are computed [@Monnet11c, pp. 43-46]. The Digital Terrain Model (DTM) represents the ground surface, it is computed by bilinear interpolation of points classified as ground. The Digital Surface Model (DSM) represents the canopy surface, it is computed by retaining in each raster's pixel the value of the highest point included in that pixel. The Canopy Height Model (CHM) is the normalized height of the canopy. It is computed by subtracting the DTM to the DSM.
+From the ALS point cloud, digital elevation models are computed [@Monnet11c, pp. 43-46]. The Digital Terrain Model (DTM) represents the ground surface, it is computed by bilinear interpolation of points classified as ground. The Digital Surface Model (DSM) represents the canopy surface, it is computed by retaining in each raster's pixel the value of the highest point contained in that pixel. The Canopy Height Model (CHM) is the normalized height of the canopy. It is computed by subtracting the DTM to the DSM.
 
 ```{r computeDEMs, include=TRUE, warning = FALSE}
 # define extent and resolution of raster
 output_raster <- terra::rast(resolution = 0.5,
-                             xmin = ROI_range$x[1] - ROI_buffer,
-                             xmax = ROI_range$x[2] + ROI_buffer,
-                             ymin = ROI_range$y[1] - ROI_buffer,
-                             ymax = ROI_range$y[2] + ROI_buffer,
+                             xmin = ROI_range$xmin - ROI_buffer,
+                             xmax = ROI_range$xmax + ROI_buffer,
+                             ymin = ROI_range$ymin - ROI_buffer,
+                             ymax = ROI_range$ymax + ROI_buffer,
                              crs = sf::st_crs(las_chablais3)$wkt
 )
 # terrain model computed from points classified as ground
@@ -178,19 +179,15 @@ A plot mask is computed from the inventoried positions, using a height-dependent
 
 ```{r computeMask, include=TRUE, message=FALSE}
 # plot mask computation based on inventoried positions
-# convex hull of plot
-mask_chull <- lidaRtRee::raster_chull_mask(tree_inventory_chablais3[, c("x", "y")], dsm)
-# union of buffers around trees
-mask_tree <- lidaRtRee::raster_xy_mask(
-  tree_inventory_chablais3[, c("x", "y")],
-  2.1 + 0.14 * tree_inventory_chablais3$h, dsm
-)
-# union of convexHull and tree buffers
-mask_plot <- max(mask_chull, mask_tree)
-# set zeros values to NA
-mask_plot[mask_plot==0] <- NA
-# vectorize plot mask
-mask_plot_v <- terra::as.polygons(mask_plot)
+# convex hull of union of points geometry
+mask_chull <- sf::st_convex_hull(sf::st_union(sf::st_geometry(tree_inventory_chablais3)))
+# union of buffers around points geometry
+mask_tree <- sf::st_union(sf::st_buffer(sf::st_geometry(tree_inventory_chablais3),
+                          2.1 + 0.14 * tree_inventory_chablais3$h))
+# union of convex hull and tree buffers
+mask_plot_v <- sf::st_union(mask_chull, mask_tree)
+# rasterize mask
+mask_plot <- terra::rasterize(terra::vect(mask_plot_v), dsm)
 ```
 Displaying inventoried trees on the CHM shows a pretty good agreement of crowns visible in the CHM with trunk locations and sizes.
 
@@ -210,8 +207,16 @@ terra::plot(mask_plot_v, border = "red", add = TRUE)
 ```
 
 ## Tree delineation
+
+Tree delineation consists in two steps, which are exemplified in the following paragraph:
+
+* image processing for local maxima detection, selection and image segmentation (function `lidaRtRee::tree_segmentation`)
++ extraction of tree information (function `lidaRtRee::tree_extraction`).
+
+The function `lidaRtRee::tree_detection` is a shortcut that combines those two steps, and allows to process a catalog of files. It is presented in the batch processing paragraph.
+
 ### Segmentation
-Tree segmentation is performed on the Canopy Height Model by using a general function (`lidaRtRee::tree_segmentation`) which consists in the following steps:
+Tree segmentation is performed on the Canopy Height Model by using the function `lidaRtRee::tree_segmentation` which consists in the following steps:
 
 * image pre-processing (non-linear filtering and smoothing for noise removal),
 + local maxima filtering and selection for apex (local maximum) detection,
@@ -234,7 +239,7 @@ terra::plot(dummy %% 8, main = "Segments (random colors)", col = rainbow(8))
 ```
 
 ### Extraction of apices positions and attributes
-A `data.frame` of detected apices located in the plot mask is then extracted with the function `tree_extraction.` Segments can be converted from raster to polygons but the operation is quite slow. Attributes are :
+A `data.frame` of detected apices located in the plot mask is then extracted with the function `tree_extraction.` Crown polygons can be vectorized from the segments. Attributes are :
 
 * `id`: apex id
 + `x`: easting coordinate of apex
@@ -243,25 +248,24 @@ A `data.frame` of detected apices located in the plot mask is then extracted wit
 + `dom_radius`: distance of apex to nearest higher pixel of CHM
 + `s`: crown surface
 + `v`: crown volume
-+ `sp`: crown surface inside plot
-+ `vp`: crown volume inside plot
++ `sp` (if plot mask is provided): crown surface inside plot
++ `vp` (if plot mask is provided): crown volume inside plot
++ `crown` (optional): 2D crown polygon in WKT format
 
 
 ```{r plotSegmTrees, include=TRUE, fig.width = 5.5, fig.height = 4.5}
 # tree extraction only inside plot mask for subsequent comparison
-apices <- lidaRtRee::tree_extraction(
-  segms$filled_dem,
-  segms$local_maxima,
-  segms$segments_id, mask_plot
-)
+apices <- lidaRtRee::tree_extraction(segms, r_mask = mask_plot, crown = TRUE)
+# convert WKT field to polygons
+crowns <- sf::st_as_sf(sf::st_drop_geometry(apices), wkt = "crown")
+# remove WKT field from apices
+apices <- apices[, -which(names(apices)=="crown")]
 head(apices, n = 3L)
-# convert segments from raster to polygons
-segments_v <- terra::as.polygons(segms$segments_id)
 #
 # display initial image
 terra::plot(chm, col = gray(seq(0, 1, 1 / 255)), main = "CHM and detected positions")
 # display segments border
-terra::plot(segments_v, border = "white", add = T)
+terra::plot(sf::st_geometry(crowns), border = "white", add = T, col = NA)
 # display plot mask
 terra::plot(mask_plot_v, border = "red", add = T)
 # display detected apices
@@ -277,12 +281,12 @@ To assess detection accuracy, reference (field) trees should be linked to detect
 # match detected apices with field trees based on relative distance of apices
 matched <- lidaRtRee::tree_matching(
   tree_inventory_chablais3[, c("x", "y", "h")],
-  cbind(sf::st_coordinates(apices), apices$h)
+  apices[, c("x", "y", "h")]
 )
 # display matching results
 lidaRtRee::plot_matched(
   tree_inventory_chablais3[, c("x", "y", "h")],
-  cbind(sf::st_coordinates(apices), apices$h), matched, chm, mask_plot_v
+  apices[, c("x", "y", "h")], matched, chm, mask_plot_v
 )
 ```
 
@@ -292,36 +296,38 @@ lidaRtRee::plot_matched(
 # height histogram of detections
 detection_stats <- lidaRtRee::hist_detection(
   tree_inventory_chablais3[, c("x", "y", "h")],
-  cbind(sf::st_coordinates(apices), apices$h), matched
+  apices[, c("x", "y", "h")], matched
 )
 ```
 
-Detection accuracy is evaluated thanks to the number of correct detections (`r detection_stats$true_detections`), false detections (`r detection_stats$false_detections`) and omissions (`r detection_stats$omissions`). In heterogeneous stands, detection accuracy is height-dependent, it is informative to display those categories on a height histogram.
+Detection accuracy is evaluated thanks to the number of correct detections (`r detection_stats$true_detections`), false detections (`r detection_stats$false_detections`) and omissions (`r detection_stats$omissions`). In heterogeneous stands, detection accuracy is height-dependent, it is informative to display those categories on a height histogram with `lidaRtRee::hist_detection`.
 
 ```{r plotDetection, include=TRUE, out.width = '40%', fig.dim=c(4.5, 4.5)}
 # height histogram of detections
 detection_stats <- lidaRtRee::hist_detection(
   tree_inventory_chablais3[, c("x", "y", "h")],
-  cbind(sf::st_coordinates(apices), apices$h), matched
+  apices[, c("x", "y", "h")], matched
 )
 ```
 
 ### Height estimation accuracy
 
 ```{r heightRegression, include=FALSE}
-height_reg <- lidaRtRee::height_regression(tree_inventory_chablais3[, c("x", "y", "h")],
-  cbind(sf::st_coordinates(apices), apices$h),
+height_reg <- lidaRtRee::height_regression(
+  tree_inventory_chablais3[, c("x", "y", "h")],
+  apices[, c("x", "y", "h")],
   matched,
   species = tree_inventory_chablais3$s
 )
 ```
 
-For true detections, estimated heights can be compared to field-measured heights. Here, the root mean square error is `r round(height_reg$stats$rmse,1)`m, while the bias is `r round(height_reg$stats$bias,1)`m. The linear regression is displayed hereafter.
+For true detections, estimated heights can be compared to field-measured heights. Here, the root mean square error is `r round(height_reg$stats$rmse,1)`m, while the bias is `r round(height_reg$stats$bias,1)`m. The linear regression is displayed hereafter (`lidaRtRee::height_regression`).
 
 ```{r plotRegression, include=TRUE, out.width = '40%', fig.dim=c(4.5, 4.5)}
 # linear regression between reference height and estimated height
-height_reg <- lidaRtRee::height_regression(tree_inventory_chablais3[, c("x", "y", "h")],
-  cbind(sf::st_coordinates(apices), apices$h),
+height_reg <- lidaRtRee::height_regression(
+  tree_inventory_chablais3[, c("x", "y", "h")],
+  apices[, c("x", "y", "h")],
   matched,
   species = tree_inventory_chablais3$s
 )
@@ -345,11 +351,11 @@ lasn <- lidR::normalize_height(las_chablais3, lidR::tin())
 # add segment id in LAS object
 lasn <- lidR::merge_spatial(lasn, segms$segments_id, "seg_id")
 # put all seg_id values in ordered list
-liste_seg_id <- sort(unique(lasn$seg_id))
+list_seg_id <- sort(unique(lasn$seg_id))
 # set names of list equal to values
-names(liste_seg_id) <- liste_seg_id
+names(list_seg_id) <- list_seg_id
 # extract point cloud for each segment id in a list
-las_l <- lapply(liste_seg_id, function(x) {lidR::filter_poi(lasn, seg_id == x)})
+las_l <- lapply(list_seg_id, function(x) {lidR::filter_poi(lasn, seg_id == x)})
 ```
 
 ### Metrics computation
@@ -374,11 +380,12 @@ Computed metrics are merged with reference trees and detected apices, thanks to
 # associate each reference tree with the segment that contains its trunk.
 dummy <- terra::extract(
   segms$segments_id,
-  tree_inventory_chablais3[, c("x", "y")]
+  sf::st_coordinates(tree_inventory_chablais3)
 )
 tree_inventory_chablais3$seg_id <- dummy$segments_id
-# create new data.frame by merging metrics and inventoried trees based on segment id
-tree_metrics <- base::merge(tree_inventory_chablais3, metrics)
+# create new data.frame by merging metrics and inventoried trees (without geometry)
+# based on segment id
+tree_metrics <- base::merge(sf::st_drop_geometry(tree_inventory_chablais3), metrics)
 # remove non-tree segment
 tree_metrics <- tree_metrics[tree_metrics$seg_id != 0, ]
 # add metrics to extracted apices data.frame
@@ -521,7 +528,7 @@ tree_metrics_h$predicted_s <- predict(lda_MASS,
                                       tree_metrics_h[, c("maxZ", "meanZ", "sdZ", "meanI", "sdI")])$class
 ```
 
-To display the classification results, an image of the classified segments is created and the highest field trees in each segment are displayed on top of it. Detentions are correct when:
+To display the classification results, an image of the classified segments is created and the highest field trees in each segment are displayed on top of it. Results are correct when:
 
 * purple dots are on purple segments (correct ABAL classification),
 * blue dots are on blue segments (correct PIAB classification),
@@ -547,7 +554,7 @@ rat$col[is.na(rat$col)] <- "green"
 levels(species)[[1]] <- rat
 # display results
 terra::plot(species, col = rat$col)
-terra::plot(segments_v, add = TRUE, border = "white")
+terra::plot(sf::st_geometry(crowns), add = TRUE, border = "white")
 lidaRtRee::plot_tree_inventory(tree_metrics_h
   [, c("x", "y")],
   tree_metrics_h$h,
@@ -557,18 +564,18 @@ lidaRtRee::plot_tree_inventory(tree_metrics_h
 )
 ```
 ```{r saveGIS, include=FALSE, warning=FALSE, message=FALSE}
-# # extract all apices
-# apices_all <- lidaRtRee::tree_extraction(segms$filled_dem, segms$local_maxima, segms$segments_id)
-# # round height values to cm
-# apices_all$h <- round(apices_all$h, 2)
-# # save outputs
-# terra::writeRaster(chm, file = "./data/output/chm.tif", overwrite = TRUE)
-# terra::writeRaster(species, file = "./data/output/r_species.tif", overwrite = TRUE)
-# terra::writeRaster(segms$segments_id, file = "./data/output/r_segments.tif", overwrite = TRUE)
-# terra::writeVector(segments_v, "./data/output/v_segments.gpkg", overwrite = TRUE)
-# sf::st_write(apices_all, "./data/output/apices.gpkg", delete_layer = TRUE)
-# write.csv(tree_metrics_h, "./data/output/correct_detections.csv", row.names = FALSE)
-# write.csv(tree_inventory_chablais3, "./data/output/tree_inventory_chablais.csv", row.names = FALSE)
+# extract all apices
+apices_all <- lidaRtRee::tree_extraction(segms)
+# round height values to cm
+apices_all$h <- round(apices_all$h, 2)
+# save outputs
+terra::writeRaster(chm, file = "../data/output/chm.tif", overwrite = TRUE)
+terra::writeRaster(species, file = "../data/output/r_species.tif", overwrite = TRUE)
+terra::writeRaster(segms$segments_id, file = "../data/output/r_segments.tif", overwrite = TRUE)
+sf::st_write(crowns, "../data/output/crowns.gpkg", delete_dsn = TRUE)
+sf::st_write(apices_all, "../data/output/apices.gpkg", delete_dsn = TRUE)
+write.csv(tree_metrics_h, "../data/output/correct_detections.csv", row.names = FALSE)
+write.csv(sf::st_drop_geometry(tree_inventory_chablais3), "../data/output/tree_inventory_chablais.csv", row.names = FALSE)
 ```
 
 ## Display point cloud
@@ -692,194 +699,111 @@ for (i in 1:nrow(apices))
 
 ## Batch processing
 
-The following code exemplifies how use the catalog engine of package `lidR` to process numerous LAS files and extract apices for the whole area with parallel processing. The catalog engine provides the data as chunks to the segmentation algorithm, and results are aggregated afterwards. 
-In order to avoid border effects, a buffer zone is added to chunks. Chunk results are cropped to prevent the same tree from appearing twice in the final results. Chunk size and buffer size are important parameters :
+The function `lidaRtRee::tree_detection` encompasses both the segmentation and extraction steps, and can take as input:
+
+- a Canopy Height Model provided as SpatRaster object,
+- a LAS object,
+- a LAS-catalog object.
+
+If the LAS object or files are not normalized, the function has an argument (`normalization = TRUE`) to do it on-the-fly. 
+
+A region of interest (ROI) can be specified in order to output only trees contained in the ROI, as in the example below which reproduces the results of previous paragraphs.
+
+```{r tree_detection_LAS, include=TRUE, fig.width = 5.5, fig.height = 4.5}
+# perform tree detection
+apices <- lidaRtRee::tree_detection(chm, ROI = mask_plot_v)
+# display initial image
+terra::plot(chm, col = gray(seq(0, 1, 1 / 255)), main = "CHM and detected positions")
+# display plot mask
+terra::plot(mask_plot_v, border = "red", add = T)
+# display detected apices
+plot(apices["h"], col = "blue", cex = apices$h / 20, pch = 2, add = TRUE)
+```
+In case of a LAS-catalog object, the `lidR` catalog engine provides the data as chunks to the segmentation algorithm, and results are aggregated afterwards. In order to avoid border effects, a buffer zone is added to chunks. Chunk results are cropped to prevent the same tree from appearing twice in the final results. Chunk size and buffer size are important parameters.
 
-- `chunk_size` is a trade-off between the number of chunks to process and the amount of RAM required to process a single chunk ;
-- `chunk_buffer` size is a trade-off between redundant processing of the overlap area, and assuring that a tree which treetop is located at the border of a chunk has its full crown within the buffer size.
+- `chunk_size` is a trade-off between the number of chunks to process and the amount of RAM required to process a single chunk. It is advisable to set its value to 0 (processing by file) in case files are non-overlapping tiles and they can be loaded in RAM.
+- `chunk_buffer` size is a trade-off between redundant processing of the overlap area, and making sure that a tree which treetop is located at the border of a chunk has its full crown within the buffer area.
 
-Chunks can be processed with parallel computing, within limits of the cluster's RAM and number of cores.
+Chunks can be processed with parallel computing, within limits of the cluster's RAM and number of cores. If a ROI is provided, only files intersecting it are processed.
 
 The steps for processing a batch of las/laz files are :
 
-- build catalog of files and set chunk processing options
+- build catalog of files and set catalog processing options
 - provide segmentation parameters
 - provide output parameters
 - set cluster options for parallel processing
-- use the catalog engine to proceed with the following steps for each chunk:
-    - compute CHM
-    - segment and extract apices
-    - remove apices inside buffer zone
-- aggregate list of results
+- apply `lidaRtRee::tree_detection` function to catalog
+- the Canopy Height Model is not an output of the function, it can be computed separately by applying `lidR::rasterize_canopy` to the catalog (if LAS files are normalized; more information on this step in the [preprocessing tutorial](https://gitlab.irstea.fr/jean-matthieu.monnet/lidartree_tutorials/-/blob/master/R/ALS_data_preprocessing.Rmd)).
 
-Be aware that in case tree segments are vectorized, some obtained polygons might overlap. The segmentation algorithm might not be deterministic and borders are sometimes not consistent when adjacent polygons are pasted from different tiles.
+Be aware that in case tree segments are vectorized, some obtained crowns might overlap. The segmentation algorithm is not be deterministic and borders might not be consistent when adjacent polygons are pasted from different chunks.
 
-```{r batch, include=TRUE, warning = FALSE}
+```{r tree_detection_LAScatalog, include=TRUE, fig.width = 5.5, fig.height = 4.5}
 rm(list = ls())
 # BUILD CATALOG OF FILES AND SET CATALOG PROCESSING OPTIONS
 # folder containing the files
 lazdir <- "../data/forest.structure.metrics"
 # build catalog and set options
-# - progress: disable display of catalog processing
+# - progress: disable or not display of catalog processing
 # - select: read only xyzc attributes (coordinates, intensity,
 # echo order and classification) from height files
 # - chunk_size: tile size to split area into chunks to process,
 # trade-off between RAM capacity VS total number of tiles to process
 # here 70 for example purpose with small area
-# - buffer_size: around each tile to avoid border effects on segmentation results
+# - chunk_buffer: around each tile to avoid border effects on segmentation results
 # trade-off between making sure a whole tree crown is processed in case its top
 # is on the border VS duplicate processing
 # 5 m is minimum, 10 is probably better depending on tree size
 cata <- lidR::readALSLAScatalog(lazdir,
                                 progress = FALSE, 
-                                select = "xyz",
+                                select = "xyzc",
                                 chunk_size = 70,
                                 chunk_buffer = 10)
 # set coordinate system
 lidR::projection(cata) <- 2154
 #
-# TREE SEGMENTATION PARAMETERS
-# set raster resolution
-res <- 1
-#
-# OUTPUT PARAMETERS
-# option to vectorize crowns (set to FALSE if too slow)
-out_vectorize_apices <- TRUE
-# output canopy height models ? (set to FALSE if too much RAM used)
-out_chm <- TRUE
-# save chms on disk
-out_save_chm <- FALSE
-#
 # CLUSTER PARAMETERS
 # specify to use two parallel sessions
 future::plan("multisession", workers = 2L)
 # remove warning when using random numbers in parallel sessions
 options(future.rng.onMisuse = "ignore")
 #
-# FUNCTION TO APPLY TO CATALOG
-routine <- function(chunk, resolution = res)
-{
-  # empty list for output
-  output <- list()
-  # bounding box of region, without buffer
-  bbox <- sf::st_bbox(chunk)
-  # read chunk of data
-  las <- lidR::readLAS(chunk)
-  # return NULL if empty
-  if (lidR::is.empty(las)) return(NULL)
-  # normalization if required
-  # las <- lidR::normalize_height(las, lidR::tin())
-  # in this example LAS tiles are already normalized
-  # compute canopy height model
-  chm <- lidR::rasterize_canopy(las, resolution, algorithm = lidR::p2r(), pkg = "terra")
-  #
-  # check all chm is not NA
-  if (all(is.na(terra::values(chm)))) return(NULL)
-  #
-  # tree detection (default parameters)
-  segms <- lidaRtRee::tree_segmentation(chm)
-  # tree extraction
-  apices <- lidaRtRee::tree_extraction(
-    segms$filled_dem,
-    segms$local_maxima,
-    segms$segments_id
-  )
-  # remove apices in buffer area
-  apices <- sf::st_crop(apices, sf::st_bbox(chunk))
-  # add tile id
-  apices$tile <- paste0(bbox$xmin, "_", bbox$ymin)
-  # put apices into output slot
-  output$apices <- apices
-  #
-  # crop chm if requested output (global variable)
-  if (out_chm | out_save_chm)
-  {
-    chm <- terra::crop(chm, terra::ext(bbox$xmin, bbox$xmax, bbox$ymin, bbox$ymax))
-  }
-  # output chm if asked for 
-  if (out_chm)
-  {
-    output$chm <- terra::wrap(chm)
-  }
-  # if save on disk
-  if (out_save_chm) {
-    terra::writeRaster(chm,
-    file = paste0("chm_", bbox$xmin, "_", bbox$ymin, ".tif"),
-    overwrite = TRUE
-    )
-  }
-  # convert to vectors if option is TRUE
-  if (out_vectorize_apices) {
-    # vectorize
-    apices_v <- terra::as.polygons(segms$segments_id)
-    # remove polygons which treetop is in buffer
-    apices_v <- apices_v[is.element(apices_v$segments_id, apices$id), ]
-    # convert to sf
-    apices_v <- sf::st_as_sf(apices_v)
-    names(apices_v)[1] <- "id"
-    # add attributes
-    apices_v <- merge(apices_v, sf::st_drop_geometry(apices), all.x = TRUE)
-    # save in list
-    output$apices_v <- apices_v
-  }
-  output
-} # end of routine function
-#
-# APPLY FUNCTION TO CATALOG
-resultats <- lidR::catalog_apply(cata, routine)
-#
-# RESULTS AGGREGATION
-# apices
-apices <- lapply(resultats, function(x) x[["apices"]])
-# remove NULL elements
-apices[sapply(apices, is.null)] <- NULL
-# bind remaining elements
-apices <- do.call(rbind, apices)
+# DETECTION PARAMETERS
+resolution <- 1
 #
-# chm
-if (out_chm) {
-  # extract and unwrap chms
-  chm <- lapply(resultats, function(x) terra::rast(x[["chm"]]))
-  # merge chm
-  # no names in list otherwise do.call returns an error
-  chm_all <- do.call(terra::merge, chm)
-}
-# apices_v
-if (out_vectorize_apices) {
-  apices_v <- lapply(resultats, function(x) x[["apices_v"]])
-  # remove NULL elements
-  apices_v[sapply(apices_v, is.null)] <- NULL
-  apices_v <- do.call(rbind, apices_v)
-  # 1-pixel overlapping in apices_v might be present because image segmentation
-  # is not fully identical in overlap areas of adjacent tiles.
-}
+# PROCESSING
+# perform tree detection with crown polygon output
+apices <- lidaRtRee::tree_detection(cata, res = resolution, crown = TRUE)
+# create crown polygons object
+crowns <- sf::st_as_sf(sf::st_drop_geometry(apices), wkt = "crown", crs = sf::st_crs(apices))
+# compute canopy height model (apply lidR::rasterize_canopy to catalog, with same resolution
+chm <- lidR::rasterize_canopy(cata, res = resolution)
 ```
 
 The following image displays the results for the whole area.
 
 ```{r batch.plot, include=TRUE, fig.width = 8, fig.height = 4.2}
 # threshold outsiders in chm
-chm_all[chm_all > 40] <- 40
-chm_all[chm_all < 0] <- 0
+chm[chm > 40] <- 40
+chm[chm < 0] <- 0
 # display chm
-terra::plot(chm_all,
+terra::plot(chm,
   main = "Canopy Height Model and segments"
 )
 # display segments border
-plot(sf::st_geometry(apices_v), border = "white", add = T)
+plot(sf::st_geometry(crowns), border = "white", add = T)
 # add apices
 plot(sf::st_geometry(apices), cex = apices$h / 40, add = TRUE, pch = 2)
 ```
 
-The following lines save outputs to files.
+The following lines save outputs to files (and overwrite existing files).
 
 ```{r batch.export, eval = FALSE}
 # merged chm
-terra::writeRaster(chm_all, file = "chm.tif", overwrite = TRUE)
+terra::writeRaster(chm, file = "chm.tif", overwrite = TRUE)
 # apices
-sf::st_write(apices, "apices_points.gpkg")# , delete_layer = TRUE)
-# vectorized apices
-if (out_vectorize_apices) sf::st_write(apices_v, "v_apices_points.gpkg")#, delete_layer = TRUE)
+sf::st_write(apices, "apices.gpkg", delete_dsn = TRUE)
+# crowns
+sf::st_write(crowns, "crowns.gpkg", delete_dsn = TRUE)
 ```
 
 ## References
\ No newline at end of file