R/hull.R

In gMOIP: Tools for 2D and 3D Plots of Single and Multi-Objective Linear/Integer Programming Models

## Functions related to the convex hull

#' Return the dimension of the convex hull of a set of points.
#'
#' @param pts A matrix/data frame/vector that can be converted to a matrix with a row for each
#'   point.
#' @param dim The dimension of the points, i.e. assume that column 1-dim specify the points. If
#'   NULL assume that the dimension are the number of columns.
#'
#' @return The dimension of the object.
#' @export
#'
#' @examples
#' ## In 1D
#' pts <- matrix(c(3), ncol = 1, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix(c(1,3,4), ncol = 1, byrow = TRUE)
#' dimFace(pts)
#'
#' ## In 2D
#' pts <- matrix(c(3,3,6,3,3,6), ncol = 2, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix(c(1,1,2,2,3,3), ncol = 2, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix(c(0,0), ncol = 2, byrow = TRUE)
#' dimFace(pts)
#'
#' ## In 3D
#' pts <- c(3,3,3,6,3,3,3,6,3,6,6,3)
#' dimFace(pts, dim = 3)
#' pts <- matrix( c(1,1,1), ncol = 3, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix( c(1,1,1,2,2,2), ncol = 3, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix(c(2,2,2,3,2,2), ncol=3, byrow= TRUE)
#' dimFace(pts)
#' pts <- matrix(c(0,0,0,0,1,1,0,2,2,0,5,2,0,6,1), ncol = 3, byrow = TRUE)
#' dimFace(pts)
#' pts <- matrix(c(0,0,0,0,1,1,0,2,2,0,0,2,1,1,1), ncol = 3, byrow = TRUE)
#' dimFace(pts)
#'
#' ## In 4D
#' pts <- matrix(c(2,2,2,3,2,2,3,4,1,2,3,4), ncol=4, byrow= TRUE)
#' dimFace(pts,)
dimFace<-function(pts, dim = NULL) {
   if (is.null(dim)) dim <- ncol(pts)
   if (is.vector(pts)) pts <- matrix(pts, ncol = dim, byrow = TRUE)
   x <- as.matrix(unique(pts[,1:dim,drop=FALSE]))
   if (nrow(x)==1) return(0)
   m1 <- x[2:dim(x)[1],]
   m2 <- matrix(rep(x[1,], times = dim(x)[1]-1), ncol = dim(x)[2], byrow = TRUE)
   x <- m1 - m2
   return(Matrix::rankMatrix(x)[1])
}


#' Efficient test for points inside a convex hull in p dimensions.
#'
#' @param pts A \eqn{nxp} array to test, \eqn{n} data points, in dimension \eqn{p}. If you have
#'   many points to test, it is most efficient to call this function once with the entire set.
#' @param vertices A \eqn{mxp} array of vertices of the convex hull. May contain redundant
#'   (non-vertex) points.
#' @param hull Tessellation (or triangulation) generated by `convhulln` (only works if the dimension
#'   of the hull is \eqn{p}). If hull is \code{NULL}, then it will be generated.
#' @param tol Tolerance on the tests for inclusion in the convex hull. You can think of `tol` as the
#'   difference a point value may be different from the values of the hull, and still be perceived
#'   as on the surface of the hull. Because of numerical slop nothing can ever be done exactly here.
#'   In higher dimensions, the numerical issues of floating point arithmetic will probably suggest a
#'   larger value of `tol`. `tol` is not used if the dimension of the hull is larger than one and
#'   not equal \eqn{p}.
#'
#' @note Some of the code are inspired by the [Matlab
#'   code](https://www.mathworks.com/matlabcentral/fileexchange/10226-inhull) by
#'   John D'Errico and
#'   [how to find a point inside a hull](https://stat.ethz.ch/pipermail/r-help/2009-December/415377.html).
#'   If the dimension of the hull is below \eqn{p} then PCA may be used to check (a
#'   warning will be given).
#'
#' @return An integer vector of length \eqn{n} with values 1 (inside hull), -1 (outside hull) or 0
#'   (on hull to precision indicated by `tol`).
#'
#' @author Lars Relund \email{lars@@relund.dk}
#' @export
#'
#' @examples
#' ## In 1D
#' vertices <- matrix(4, ncol = 1)
#' pt <- matrix(c(2,4), ncol = 1, byrow = TRUE)
#' inHull(pt, vertices)
#' vertices <- matrix(c(1,4), ncol = 1)
#' pt <- matrix(c(1,3,4,5), ncol = 1, byrow = TRUE)
#' inHull(pt, vertices)
#'
#' ## In 2D
#' vertices <- matrix(c(2,4), ncol = 2)
#' pt <- matrix(c(2,4, 1,1), ncol = 2, byrow = TRUE)
#' inHull(pt, vertices)
#' vertices <- matrix(c(0,0, 3,3), ncol = 2, byrow = TRUE)
#' pt <- matrix(c(0,0, 1,1, 2,2, 3,3, 4,4), ncol = 2, byrow = TRUE)
#' inHull(pt, vertices)
#' vertices <- matrix(c(0,0, 0,3, 3,0), ncol = 2, byrow = TRUE)
#' pt <- matrix(c(0,0, 1,1, 4,4), ncol = 2, byrow = TRUE)
#' inHull(pt, vertices)
#'
#' ## in 3D
#' vertices <- matrix(c(2,2,2), ncol = 3, byrow = TRUE)
#' pt <- matrix(c(1,1,1, 3,3,3, 2,2,2, 3,3,2), ncol = 3, byrow = TRUE)
#' inHull(pt, vertices)
#'
#' vertices <- matrix(c(2,2,2, 4,4,4), ncol = 3, byrow = TRUE)
#' ini3D()
#' plotHull3D(vertices)
#' pt <- matrix(c(1,1,1, 2,2,2, 3,3,3, 4,4,4, 3,3,2), ncol = 3, byrow = TRUE)
#' plotPoints3D(pt, addText = TRUE)
#' finalize3D()
#' inHull(pt, vertices)
#'
#' vertices <- matrix(c(1,0,0, 1,1,0, 1,0,1), ncol = 3, byrow = TRUE)
#' ini3D()
#' plotHull3D(vertices)
#' pt <- matrix(c(1,0.1,0.2, 3,3,2), ncol = 3, byrow = TRUE)
#' plotPoints3D(pt, addText = TRUE)
#' finalize3D()
#' inHull(pt, vertices)
#'
#' vertices <- matrix(c(2,2,2, 2,4,4, 2,2,4, 4,4,2, 4,2,2, 2,4,2, 4,2,4, 4,4,4), ncol = 3,
#'             byrow = TRUE)
#' ini3D()
#' plotHull3D(vertices)
#' pt <- matrix(c(1,1,1, 3,3,3, 2,2,2, 3,3,2), ncol = 3, byrow = TRUE)
#' plotPoints3D(pt, addText = TRUE)
#' finalize3D()
#' inHull(pt, vertices)
#'
#' ## In 5D
#' vertices <- matrix(c(4,0,0,0,0, 0,4,0,0,0, 0,0,4,0,0, 0,0,0,4,0, 0,0,0,0,4, 0,0,0,0,0),
#'             ncol = 5, byrow = TRUE)
#' pt <- matrix(c(0.1,0.1,0.1,0.1,0.1, 3,3,3,3,3, 2,0,0,0,0), ncol = 5, byrow = TRUE)
#' inHull(pt, vertices)
inHull <- function(pts, vertices, hull=NULL,
                   tol=mean(mean(abs(as.matrix(vertices))))*sqrt(.Machine$double.eps)) {
   # ensure arguments are matrices (not data frames) and get sizes
   pts <- .checkPts(pts)
   vertices <- .checkPts(vertices)

   p <- ncol(vertices)   # dimension of points
   cx <- nrow(pts)   # points to test
   d <- dimFace(vertices, dim = p)
   if (d == 0) return(-as.integer(apply(pts, 1, FUN = function(x) any(abs(x - vertices[1, ]) > tol))))
   if (p == 1) { #1D space - check if between max and min
      m <- min(vertices)
      M <- max(vertices)
      val <- rep(0, cx)  # first set all on hull
      val[pts > M + tol | pts < m - tol] <- -1
      val[pts < M - tol & pts > m + tol] <- 1
      return(val)
   }
   if (d == 1) { # a line (i.e. no interior points)
      denom <- vertices[2,]-vertices[1,]  # the direction d of the line p = d * t + p0
      # case: if d_i = 0 => p_i = p0_i otherwise not on line
      zeros <- denom == 0
      if (any(zeros)) {
         val <- apply(pts[, zeros, drop=F], 1, FUN = function(x) {any(x != vertices[1, zeros])})  # are some points not on the line
         # reduce dimension where d_i = 0
         vertices <- vertices[, !zeros, drop = F]
         pts <- pts[, !zeros, drop=F]
         res <- inHull(pts, vertices, tol = tol)
         res[res == 1] <- 0   # set interior pts to true value
         res[res == -1] <- 1
         res <- as.logical(res)
         return(-as.integer(res | val))
      } else {
         l <- apply(pts, 1, FUN = function(x) {(x-vertices[1,])/denom})  # calc t values t = d_i/(p_i-p0_i)
         tRng <- range(apply(vertices, 1, FUN = function(x) {(x[1] - vertices[1,1])/denom[1]})) # t range for vertices
         eq <- apply(l, 2, FUN = function(x) abs(max(x) - min(x)) < tol) # check if values in each col are equal (i.e. on line)
         l <- apply(l, 2, FUN = function(x) max(x) < tRng[2] + tol & min(x) > tRng[1] - tol) # check if t values in range
         l <- eq & l
         return(-as.integer(!l))
      }
   }
   if (p == d) {
      if (is.null(hull)) hull <- geometry::convhulln(vertices)
      nt <- nrow(hull)    # number of simplexes in hull
      # find normal vectors to each simplex
      nrmls <- matrix(NA, nt, p)         # predefine each nrml as NA, degenerate
      degenflag <- matrix(TRUE, nt, 1)
      for (i in  1:nt) {
         nullsp <-
            t(MASS::Null(t(vertices[hull[i,-1],] - matrix(vertices[hull[i, 1],], p - 1, p, byrow = TRUE))))
         if (dim(nullsp)[1] == 1) {
            nrmls[i, ] <- nullsp
            degenflag[i] <- FALSE
         }
      }
      # Warn of degenerate faces, and remove corresponding normals
      if (sum(degenflag) > 0)
         warning(sum(degenflag), " degenerate faces in convex hull")
      nrmls <- nrmls[!degenflag, ]
      nt <- nrow(nrmls)
      # find center point in hull, and any (1st) point in the plane of each simplex
      center = colMeans(vertices)
      a <- vertices[hull[!degenflag, 1], ]
      # scale normal vectors to unit length and ensure pointing inwards
      nrmls <- nrmls / matrix(apply(nrmls, 1, function(x) sqrt(sum(x ^ 2))), nt, p)
      dp <- sign(apply((matrix(center, nt, p, byrow = TRUE) - a) * nrmls, 1, sum))
      nrmls <- nrmls * matrix(dp, nt, p)
      # if  min across all faces of dot((x - a),nrml) is
      #      +ve then x is inside hull
      #      0   then x is on hull
      #      -ve then x is outside hull
      # Instead of dot((x - a),nrml)  use dot(x,nrml) - dot(a, nrml)
      aN <- diag(a %*% t(nrmls))
      val <- apply(pts %*% t(nrmls) - matrix(aN, cx, nt, byrow = TRUE), 1, min)
      # code  values inside 'tol' to zero, return sign as integer
      val[abs(val) < tol] <- 0
      return(as.integer(sign(val)))
   }
   if (p != d) {  # note p>2
      outside <- rep(2,cx)
      row.names(pts) <- 1:nrow(pts)
      nrp <- nrow(vertices)
      done <- FALSE
      while(!done){
         pt <- pts[which(outside != -1),,drop=F]
         set <- rbind(vertices,pt)
         dS <- dimFace(set, dim = p)
         if (p == dS) {
            val <- geometry::convhulln(set)
            idx <- as.integer(row.names(pt[unique(val[val>nrp])-nrp,,drop=F]))
            outside[idx] <- -1
            done <- all(outside != 2)
         }
         if (dS != p) {
            idx <- as.integer(row.names(pt))
            warning(
               "Using PCA to transform. Results for pts index ",
               paste(idx, collapse = ","),
               " might not be accurate!"
            )
            set <- stats::princomp(set)
            idx <- which(set$sdev > tol)
            set <- stats::predict(set)[,idx]
            # set <- prcomp(set, scale. = T, retx = T, rank = dS)$x
            val <- inHull(set[(nrp+1):nrow(set),], set[1:nrp,], tol = 5*tol)
            val[val == 1] <- 0
            outside[as.integer(row.names(pt))] <- val
            return(outside)
         }
      }
      return(outside)
   }
}


#' Find segments (lines) of a face.
#'
#' @param vertices A `m x p` array of vertices of the convex hull, as used by
#'   [geometry::convhulln()].
#' @param hull Tessellation (or triangulation) generated by [geometry::convhulln()] If hull is
#'   left empty or not supplied, then it will be generated.
#' @param tol Tolerance on the tests for inclusion in the convex hull. You can
#'   think of `tol` as the distance a point may possibly lie outside the hull, and
#'   still be perceived as on the surface of the hull. Because of numerical slop
#'   nothing can ever be done exactly here. I might guess a semi-intelligent
#'   value of `tol` to be
#'
#'   `tol = 1.e-13*mean(abs(vertices(:)))`
#'
#'   In higher dimensions, the numerical issues of floating point arithmetic
#'   will probably suggest a larger value of `tol`.
#'
#' @return A matrix with segments.
#' @author Lars Relund \email{lars@@relund.dk}
hullSegment <- function(vertices, hull=geometry::convhulln(vertices),
                        tol=mean(mean(abs(vertices)))*sqrt(.Machine$double.eps)) {
   vertices <- as.matrix(vertices)
   p <- ncol(vertices)   # columns in vertices
   nt <- nrow(hull)    # number of simplexes in hull
   # find normal vectors to each simplex
   nrmls <- matrix(NA, nt, p)         # predefine each nrml as NA, degenerate
   degenflag <- matrix(TRUE, nt, 1)
   for (i in  1:nt) {
      nullsp <-
         t(MASS::Null(t(
            vertices[hull[i,-1],] - matrix(vertices[hull[i, 1],], p - 1, p, byrow =
                                              TRUE)
         )))
      if (dim(nullsp)[1] == 1) {
         nrmls[i, ] <- nullsp
         degenflag[i] <- FALSE
      }
   }
   nt <- nrow(nrmls)
   # find center point in hull, and any (1st) point in the plane of each simplex
   center = colMeans(vertices)
   a <- vertices[hull[, 1], ]
   # scale normal vectors to unit length and ensure pointing inwards
   nrmls <- nrmls / matrix(apply(nrmls, 1, function(x) sqrt(sum(x ^ 2))), nt, p)
   dp <- sign(apply((matrix(center, nt, p, byrow = TRUE) - a) * nrmls, 1, sum))
   nrmls <- nrmls * matrix(dp, nt, p)
   idx <- !duplicated(round(nrmls,10))
   nrmls <- nrmls[idx,]
   a <- a[idx,]
   aN <- diag(a %*% t(nrmls))
   res <- NULL
   m <- nrow(vertices)
   for( i in 1:m )
      for( j in 1:m )
         if( i < j ) {
            # Middle of the segment
            p <- .5 * vertices[i,] + .5 * vertices[j,]
            # Check if it is at the intersection of two planes
            tmp <- p %*% t(nrmls) - aN #matrix(aN, 1, nrow(nrmls), byrow = TRUE)
            if(sum( abs( tmp ) < 1e-6 ) >= 2) res <- rbind(res, c(i,j))
         }
   return(as.matrix(res))
}


#' Add all points on the bounding box hit by the rays.
#'
#' @param pts A data frame with all points
#' @param m Minimum values of the bounding box.
#' @param M Maximum values of the bounding box.
#' @param direction Ray direction. If i'th entry is positive, consider the i'th column of the `pts`
#'   plus a value greater than on equal zero. If negative, consider the i'th column of the `pts`
#'   minus a value greater than on equal zero.
#'
#' @note Assume that `pts` has been checked using [.checkPts()].
#' @return The points merged with the points on the bounding box. The column `pt` equals 1 if
#'   points from `pts` and zero otherwise.
#' @export
#'
#' @examples
#' pts <- genNDSet(3,10)[,1:3]
#' addRays(pts)
#' addRays(pts, dir = c(1,-1,1))
#' addRays(pts, dir = c(-1,-1,1), m = c(0,0,0), M = c(100,100,100))
#' pts <- genSample(5,20)[,1:5]
#' addRays(pts)
addRays <- function(pts, m = apply(pts,2,min)-5, M = apply(pts,2,max)+5, direction = 1) {
   pts <- as.data.frame(pts)
   p <- ncol(pts)
   if (length(direction) != p) direction = rep(direction[1],p)
   v <- purrr::map_dbl(1:p, function(i) if (sign(direction[i]) > 0) M[i] else m[i])
   names(v) <- colnames(pts)
   set <- cbind(pts, pt = 1) # point in (1=original set, 0=artificial)

   for (i in 1:nrow(pts)) {
      x <- dplyr::bind_rows(pts[i,,drop=FALSE],v)
      x <- as.list(x)
      x <- expand.grid(x)
      x$pt <- 0
      set <- rbind(set,x)
   }
   # colnames(set)[1:p] <- paste0("z",1:p)
   # dplyr::distinct(set, V1, .keep_all = TRUE)
   # rownames(set) <- 1:dim(set)[1]
   set <- set[rownames(unique(set[,1:p,drop=FALSE])),, drop=FALSE]
   # rownames(set)[set$pt==0][6] <- (nrow(pts)+1):nrow(set)  # may give error if number in first rows
   return(set)
}


#' Find the convex hull of a set of points.
#'
#' @param pts A matrix with a point in each row.
#' @param addRays Add the ray defined by `direction`.
#' @param useRGLBBox Use the RGL bounding box when add rays.
#' @param direction Ray direction. If i'th entry is positive, consider the i'th column of `pts`
#'   plus a value greater than on equal zero (minimize objective $i$). If negative, consider the
#'   i'th column of `pts` minus a value greater than on equal zero (maximize objective $i$).
#' @param tol Tolerance on standard deviation if using PCA.
#' @param m Minimum values of the bounding box.
#' @param M Maximum values of the bounding box.
#'
#' @return A list with \code{hull} equal a matrix with row indices of the vertices defining each
#'   facet in the hull and \code{pts} equal the input points (and dummy points) and columns:
#'   \code{pt}, true if a point in the original input; false if a dummy point (a point on a ray).
#'   \code{vtx}, TRUE if a vertex in the hull.
#'
#' @examples
#' ## 1D
#' pts<-matrix(c(1,2,3), ncol = 1, byrow = TRUE)
#' dimFace(pts) # a line
#' convexHull(pts)
#' convexHull(pts, addRays = TRUE)
#'
#' ## 2D
#' pts<-matrix(c(1,1, 2,2), ncol = 2, byrow = TRUE)
#' dimFace(pts) # a line
#' convexHull(pts)
#' plotHull2D(pts, drawPoints = TRUE)
#' convexHull(pts, addRays = TRUE)
#' plotHull2D(pts, addRays = TRUE, drawPoints = TRUE)
#' pts<-matrix(c(1,1, 2,2, 0,1), ncol = 2, byrow = TRUE)
#' dimFace(pts) # a polygon
#' convexHull(pts)
#' plotHull2D(pts, drawPoints = TRUE)
#' convexHull(pts, addRays = TRUE, direction = c(-1,1))
#' plotHull2D(pts, addRays = TRUE, direction = c(-1,1), addText = "coord")
#'
#' ## 3D
#' pts<-matrix(c(1,1,1), ncol = 3, byrow = TRUE)
#' dimFace(pts) # a point
#' convexHull(pts)
#' pts<-matrix(c(0,0,0,1,1,1,2,2,2,3,3,3), ncol = 3, byrow = TRUE)
#' dimFace(pts) # a line
#' convexHull(pts)
#' pts<-matrix(c(0,0,0,0,1,1,0,2,2,0,0,2), ncol = 3, byrow = TRUE)
#' dimFace(pts) # a polygon
#' convexHull(pts)
#' convexHull(pts, addRays = TRUE)
#' pts<-matrix(c(1,0,0,1,1,1,1,2,2,3,1,1), ncol = 3, byrow = TRUE)
#' dimFace(pts) # a polygon
#' convexHull(pts) # a polyhedron
#' pts<-matrix(c(1,1,1,2,2,1,2,1,1,1,1,2), ncol = 3, byrow = TRUE)
#' dimFace(pts) # a polytope (polyhedron)
#' convexHull(pts)
#'
#' ini3D(argsPlot3d = list(xlim = c(0,3), ylim = c(0,3), zlim = c(0,3)))
#' pts<-matrix(c(1,1,1,2,2,1,2,1,1,1,1,2), ncol = 3, byrow = TRUE)
#' plotPoints3D(pts)
#' plotHull3D(pts, argsPolygon3d = list(color = "red"))
#' convexHull(pts)
#' plotHull3D(pts, addRays = TRUE)
#' convexHull(pts, addRays = TRUE)
#' finalize3D()
#'
#' @importFrom rlang .data
#' @export
convexHull <- function(pts, addRays = FALSE, useRGLBBox = FALSE, direction = 1,
                       tol = mean(mean(abs(pts)))*sqrt(.Machine$double.eps)*2,
                       m = apply(pts,2,min)-5, M = apply(pts,2,max)+5) {
   pts <- .checkPts(pts)
   p <- ncol(pts)
   if (length(direction) != p) direction = rep(direction[1],p)
   # print(set)
   if (addRays) {
      if (rgl::cur3d() > 0 & useRGLBBox) {
         limits <- rgl::par3d()$bbox
         for (i in 1:dim(pts)[1]) {
            pt <- as.vector(pts[i,])
            if (!(limits[1] < pt[1] && pt[1] < limits[2] &&
                  limits[3] < pt[2] && pt[2] < limits[4] &&
                  limits[5] < pt[3] && pt[3] < limits[6])){
               stop("The point is not in the interior of the current bounding box. Resize your axes.")
            }
         }
         m <- c(limits[1], limits[3], limits[5])
         M <- c(limits[2], limits[4], limits[6])
         set <- addRays(pts, m, M, direction)
      } else set <- addRays(pts, m, M, direction = direction)
   } else {
      set <- cbind(pts, pt = 1) # point in (1=original set, 0=artificial)
   }
   d <- dimFace(set[,1:p, drop = FALSE])
   # cat("Object of dimension",d,"\n")

   ## note p can be > 3 !!
   hull <- NULL
   if (d==0) { # a point
      hull <- 1
   } else if (d==p) { # same dimension as space
      if (p == 1) hull <- grDevices::chull(cbind(set[,1:p],0))
      if (p == 2) hull <- grDevices::chull(set[,1:p])
      if (p > 2) hull <- geometry::convhulln(set[,1:p], return.non.triangulated.facets = TRUE)
   } else if (d==1) { # a line (d<p)
      l <- nrow(set)
      n <- p
      comb <- t(utils::combn(n,2))
      for (i in 1:dim(comb)[1]) {  # simple projection down on each axis
         pt <- unique(set[,comb[i,]])
         if (l == dim(pt)[1]) {
            hull <- grDevices::chull(pt)
            if (length(hull)==2) {
               hull <- as.numeric(rownames(pt[hull,]))
               break
            }
         }
      }
   } else if (d==2) { # a plane (d<p)
      l <- nrow(set)
      tmp <- cbind(set, id = 1:l)
      comb <- t(utils::combn(p,2))
      for (i in 1:dim(comb)[1]) {  # simple projection down on each axis (keep the projection with highest number of vertices)
         pt <- unique(set[,comb[i,]])
         iL <- 0
         if (l == dim(pt)[1]) {
            hull <- grDevices::chull(pt)
            if (length(hull)>iL) {
               iL <- length(hull)
               hB <- rownames(pt[hull,])
            }
         }
      }
      hull <- tmp[hB,'id'] # return idx of set
   } else { # d > 2 & d < p
         warning(
            "Using PCA to transform. Results might not be accurate!"
         )
         set1 <- stats::prcomp(set[,1:p])
         # set1 <- stats::princomp(set[,1:p])
         idx <- which(set1$sdev > tol)
         set1 <- stats::predict(set1)[,idx, drop = FALSE]
         if (ncol(set1) == 1) hull <- grDevices::chull(cbind(set1,0))
         if (ncol(set1) == 2) hull <- grDevices::chull(set1)
         if (ncol(set1) > 2) hull <- geometry::convhulln(set1, return.non.triangulated.facets = TRUE)
   }
   if (is.vector(hull)) hull <- matrix(hull, nrow = 1)
   ## classify
   idx <- unique(as.vector(hull))
   idx <- idx[!is.na(idx)]
   idx <- sort(idx)
   set <- as.data.frame(set)
   set$vtx <- FALSE
   set$vtx[idx] <- TRUE
   # remove unwarnted dummy points
   set$oldId <- 1:nrow(set)
   tmp <- dplyr::filter(set, .data$vtx | .data$pt, .preserve = TRUE)
   if (nrow(set) != nrow(tmp)) {
      tmp$newId <- 1:nrow(tmp)
      mch <- dplyr::left_join(set, tmp, by = 'oldId') %>% dplyr::select("newId")
      hull <- apply(hull, c(1,2), function(id) {if (is.na(id)) return(NA) else return(mch[id,1])} )
      set <- dplyr::select(tmp, -"oldId", -"newId")
   } else set <- dplyr::select(tmp, -"oldId")
   return(list(hull = hull, pts = set))
   stop("Cannot find the vertices!")
}