diff --git a/NAMESPACE b/NAMESPACE
index 49cc9c5fc6c45b0526080ddea2f7b361114cb73d..116c12b39690c01744085c1a39497579692e1aa2 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -13,6 +13,7 @@ export(cimg2Raster)
 export(circle2Raster)
 export(cleanPredictionRaster)
 export(cloudMetrics)
+export(cloudTreeMetrics)
 export(coregistration)
 export(createDisk)
 export(demFiltering)
diff --git a/R/common.R b/R/common.R
index 8d9b1d6fccd2c6279134bb1640a619150987ed08..1283c5d24faf965c24176521dd56b5ca3e3a6106 100644
--- a/R/common.R
+++ b/R/common.R
@@ -372,8 +372,8 @@ plotTreeInventory <-function(xy, height=NULL, diam=NULL, species=NULL, add=FALSE
 #' raster::res(r) <- 1
 #' 
 #' # xy positions
-#' xy <- data.frame(c(10,20,31.25,15),
-#'   c(10,20,31.25,25))
+#' xy <- data.frame(x=c(10,20,31.25,15),
+#'   y=c(10,20,31.25,25))
 #' # compute mask
 #' mask1 <- rasterXYMask(xy, c(5, 8, 5, 5), r)
 #' mask2 <- rasterXYMask(xy, c(5, 8, 5, 5), r, binary=FALSE)
@@ -390,7 +390,7 @@ plotTreeInventory <-function(xy, height=NULL, diam=NULL, species=NULL, add=FALSE
 rasterXYMask <- function(xy, buff, r, binary=TRUE)
 {
   # convert vector to data.frame in case only one pair of coodinates is provided
-  if (length(xy) == 2) {xy <- data.frame(matrix(coord, 1, 2))}
+  if (is.null(dim(xy)) & length(xy)==2) {xy <- data.frame(matrix(xy, 1, 2))}
   # compute squared buffers
   buff2 <- buff^2
   # compute XY coordinates of cell centers
diff --git a/R/metrics.R b/R/metrics.R
index 3e9ecaacca5e1203e73614585c8b9341fe3033f2..006a5e894ceceb7bbb84996d443260acc648c684 100644
--- a/R/metrics.R
+++ b/R/metrics.R
@@ -155,7 +155,7 @@ ABAmodelMetrics <- function(z, i, rn, c, hmin = 2, breaksH=NULL)
 #'
 stdTreeMetrics <- function(x, area.ha=NA)
 {
-  data.frame(Tree.meanH=mean(x$h), Tree.sdH=stats::sd(x$h), Tree.giniH=reldist::gini(x$h), Tree.density=length(x$h)/area.ha, TreeInf10.density=sum(x$h<=10)/area.ha, TreeSup10.density=sum(x$h>10)/area.ha, TreeSup20.density=sum(x$h>20)/area.ha, TreeSup30.density=sum(x$h>30)/area.ha, Tree.meanCrownSurface=mean(x$s), Tree.CanopyCover=sum(x$sp)/area.ha/10000, Tree.meanCrownV=mean(x$v), Tree.meanCanopyH=sum(x$vp)/area.ha/10000, Tree.meanCrownH=sum(x$v)/sum(x$s))
+  data.frame(Tree.meanH=mean(x$h), Tree.sdH=stats::sd(x$h), Tree.giniH=ifelse(is.null(x), NA, reldist::gini(x$h)), Tree.density=length(x$h)/area.ha, TreeInf10.density=sum(x$h<=10)/area.ha, TreeSup10.density=sum(x$h>10)/area.ha, TreeSup20.density=sum(x$h>20)/area.ha, TreeSup30.density=sum(x$h>30)/area.ha, Tree.meanCrownSurface=mean(x$s), Tree.CanopyCover=sum(x$sp)/area.ha/10000, Tree.meanCrownV=mean(x$v), Tree.meanCanopyH=sum(x$vp)/area.ha/10000, Tree.meanCrownH=sum(x$v)/sum(x$s))
   }
 #
 ###############################
@@ -219,24 +219,50 @@ points2terrainStats <- function(p, centre=NULL, r=NULL)
     round(data.frame(altitude=altitude, azimut.gr=azimut, slope.gr=slope, adjR2.plane=summary(modlin)$adj.r.squared * 100),1)
   } else { NULL}
 }
-
 ################################
-#' compute tree metrics for list of LAS objects (should be normalized clouds, with buffer around area of interest). Call the treedetection function on each object and then arranges the results in a data.frame
+#' Computation of tree metrics
+#' 
+#' This function performs tree detection in each of LAS objects contained in a list and then computes summary statistics based on compute tree metrics for list of LAS objects (should be normalized clouds, with buffer around area of interest).
+#' - calls the tree detection and extraction functions on each object
+#' - remove the detected trees located outside of the regions of interest defined by their centers and radius
+#' - computes summary statistics of remaining detected trees.
 #'
 #' @param llasn list of LAS objects (from lidR package)
 #' @param XY a dataframe or matrix with XY coordinates of plot centers
-#' @param plot.radius a numeric plot radius in m
-#' @param res CHM resolution for tree detection (m)
-#' @param hmin min height for a tree (m)
-#' @param h.points height for CHM thresholding (m)
+#' @param plot.radius numeric. plot radius in same unit as coordinates of LAS objects
+#' @param res numeric. resolution of canopy height model computed with \code{\link{points2DSM}} before tree segmentation
 #' @param func a function to be applied to the attributes of extracted trees (return from internal call to treeExtraction function to compute plot level metrics
+#' @param ... other parameters to be passed to \code{\link{treeSegmentation}}
 #' @return a dataframe with tree metrics in columns corresponding to LAS objects of the list (lines)
+#' @seealso \code{\link{treeSegmentation}}, \code{\link{treeExtraction}}, \code{\link{stdTreeMetrics}}
+#' @examples
+#' data(laschablais3)
+#' 
+#' # extract three point clouds from LAS object
+#' llas <- list()
+#' llas[["A"]] <- lidR::lasclipCircle(laschablais3, 974350, 6581680, 10)
+#' llas[["B"]] <- lidR::lasclipCircle(laschablais3, 974390, 6581680, 10)
+#' llas[["C"]] <- lidR::lasclipCircle(laschablais3, 974350, 6581640, 10)
+#' # normalize point clouds
+#' llas <- lapply(llas, function(x) {lidR::lasnormalize(x, lidR::tin())})
+#' 
+#' # compute tree metrics
+#' cloudTreeMetrics(llas, 
+#'                  cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)),
+#'                  10, res=0.5)
+#' 
+#' # compute metrics with user-defined function
+#' # number of detected trees between 20 and 30 meters and their mean height
+#' user.func <- function(x, area.ha=NA)
+#' {
+#'   dummy <- x$h[which(x$h>10 & x$h<20)]
+#'   data.frame(Tree.between.20.30=length(dummy)/area.ha, Tree.meanH=mean(dummy))
+#' }
+#'  cloudTreeMetrics(llas, cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)), 10, res=0.5, func=user.func)
 #' @export
 #'
-cloudTreeMetrics <- function(llasn, XY, plot.radius, res=0.5, hmin=5, h.points=60, func=stdTreeMetrics)
+cloudTreeMetrics <- function(llasn, XY, plot.radius, res=0.5, func=stdTreeMetrics, ...)
 {
-  # ADD example
-  # ALLOW ... parameters to be passed to treeSegmentation
   plot.area.ha <- pi*plot.radius^2/10000
   xy.list <- split(XY, seq(nrow(XY)))
   # LOOP on list <- replace by lapply
@@ -248,16 +274,18 @@ cloudTreeMetrics <- function(llasn, XY, plot.radius, res=0.5, hmin=5, h.points=6
     # compute dsm
     dummy <- points2DSM(x,res=res)
     # replace NA, low and high values
-    dummy[is.na(dummy) | dummy<0 | dummy>h.points] <- 0
+    dummy[is.na(dummy) | dummy<0] <- 0
     # tree detection
-    dummy <-  treeSegmentation(dummy, hmin=hmin)
+    dummy <-  treeSegmentation(dummy, ...)
     # tree extraction
-    lsegms[[i]] <- treeExtraction(dummy$filled.dem, dummy$local.maxima, dummy$segments.id, rasterXYMask(coord, plot.radius, dummy$local.maxima, binary=T))
+    dummy <- treeExtraction(dummy$filled.dem, dummy$local.maxima, dummy$segments.id, rasterXYMask(coord, plot.radius, dummy$local.maxima, binary=T))
+    lsegms[[names(llasn)[i]]] <- dummy
   }
-  # GERER les cas ou pas d'arbre detecte ??
-  # GERER les NA 
-  # compute metrics
-  tree.metrics  <- lapply(lsegms, FUN=function(x){func(x, area.ha=plot.area.ha)})
-  # bind data.frames
-  tree.metrics <- do.call(rbind,tree.metrics)
-}
+  #
+  if (length(lsegms)>0)
+  {
+    tree.metrics  <- lapply(lsegms, FUN=function(x){func(x, area.ha=plot.area.ha)})
+    tree.metrics <- do.call(rbind,tree.metrics)
+  } else {tree.metrics <- NULL}
+  return(tree.metrics)
+}
\ No newline at end of file
diff --git a/R/tree.detection.R b/R/tree.detection.R
index 62767ba4f2380ea9b2ec04f10d7d3d2bda7ebb28..6a35f33417b8ed79e35f425bfd5112fb3ba194c4 100644
--- a/R/tree.detection.R
+++ b/R/tree.detection.R
@@ -541,6 +541,11 @@ treeExtraction <- function(r.dem.nl, r.maxi, r.dem.w, r.mask=NULL)
   #
   # extract segment id, height and dom.radius
   cells <- which(raster::values(r.maxi)>0)
+  # check if positive values are present
+  if (length(cells)==0)
+  {
+    segms <- NULL
+  } else {
   coord <- raster::xyFromCell(r.maxi,cells)
   segms <- data.frame(coord,id=raster::values(r.dem.w)[cells],h=raster::values(r.dem.nl)[cells],dom.radius=raster::values(r.maxi)[cells])
   segms <- merge(segms, s, all.x = T)
@@ -552,6 +557,7 @@ treeExtraction <- function(r.dem.nl, r.maxi, r.dem.w, r.mask=NULL)
   }
   sp::coordinates(segms) <- ~ x+y
   segms@proj4string <- r.dem.nl@crs
+  }
   segms
 }
 #
diff --git a/man/cloudTreeMetrics.Rd b/man/cloudTreeMetrics.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..e408d24ddb0c7a5d837b5b3a54dae6de2bd897e0
--- /dev/null
+++ b/man/cloudTreeMetrics.Rd
@@ -0,0 +1,59 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/metrics.R
+\name{cloudTreeMetrics}
+\alias{cloudTreeMetrics}
+\title{Computation of tree metrics}
+\usage{
+cloudTreeMetrics(llasn, XY, plot.radius, res = 0.5,
+  func = stdTreeMetrics, ...)
+}
+\arguments{
+\item{llasn}{list of LAS objects (from lidR package)}
+
+\item{XY}{a dataframe or matrix with XY coordinates of plot centers}
+
+\item{plot.radius}{numeric. plot radius in same unit as coordinates of LAS objects}
+
+\item{res}{numeric. resolution of canopy height model computed with \code{\link{points2DSM}} before tree segmentation}
+
+\item{func}{a function to be applied to the attributes of extracted trees (return from internal call to treeExtraction function to compute plot level metrics}
+
+\item{...}{other parameters to be passed to \code{\link{treeSegmentation}}}
+}
+\value{
+a dataframe with tree metrics in columns corresponding to LAS objects of the list (lines)
+}
+\description{
+This function performs tree detection in each of LAS objects contained in a list and then computes summary statistics based on compute tree metrics for list of LAS objects (should be normalized clouds, with buffer around area of interest).
+- calls the tree detection and extraction functions on each object
+- remove the detected trees located outside of the regions of interest defined by their centers and radius
+- computes summary statistics of remaining detected trees.
+}
+\examples{
+data(laschablais3)
+
+# extract three point clouds from LAS object
+llas <- list()
+llas[["A"]] <- lidR::lasclipCircle(laschablais3, 974350, 6581680, 10)
+llas[["B"]] <- lidR::lasclipCircle(laschablais3, 974390, 6581680, 10)
+llas[["C"]] <- lidR::lasclipCircle(laschablais3, 974350, 6581640, 10)
+# normalize point clouds
+llas <- lapply(llas, function(x) {lidR::lasnormalize(x, lidR::tin())})
+
+# compute tree metrics
+cloudTreeMetrics(llas, 
+                 cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)),
+                 10, res=0.5)
+
+# compute metrics with user-defined function
+# number of detected trees between 20 and 30 meters and their mean height
+user.func <- function(x, area.ha=NA)
+{
+  dummy <- x$h[which(x$h>10 & x$h<20)]
+  data.frame(Tree.between.20.30=length(dummy)/area.ha, Tree.meanH=mean(dummy))
+}
+ cloudTreeMetrics(llas, cbind(c(974350, 974390, 974350), c(6581680, 6581680, 6581640)), 10, res=0.5, func=user.func)
+}
+\seealso{
+\code{\link{treeSegmentation}}, \code{\link{treeExtraction}}, \code{\link{stdTreeMetrics}}
+}
diff --git a/man/rasterXYMask.Rd b/man/rasterXYMask.Rd
index 19d6d693a2bb7c2adea3471b7e18e341161cdf8a..d30a33dcffb6d14314afcdb7f428214072016f3b 100644
--- a/man/rasterXYMask.Rd
+++ b/man/rasterXYMask.Rd
@@ -28,8 +28,8 @@ raster::extent(r) <- c(0,40,0,40)
 raster::res(r) <- 1
 
 # xy positions
-xy <- data.frame(c(10,20,31.25,15),
-  c(10,20,31.25,25))
+xy <- data.frame(x=c(10,20,31.25,15),
+  y=c(10,20,31.25,25))
 # compute mask
 mask1 <- rasterXYMask(xy, c(5, 8, 5, 5), r)
 mask2 <- rasterXYMask(xy, c(5, 8, 5, 5), r, binary=FALSE)