Hi, Below Baptiste's message I attach the R code and the .Rd documentation I treated as 'final', it may be slightly different from that in the Dec 2009 post.
I did submit if for inclusion in the geometry package, but last time I checked it wasn't there. I have found (and others have reported via private e-mail) that high dimensional data sets cause crashes (R exits with no warnings or messages). I tentatively believe this is a 'bug' in convhulln, but tracking it takes me to a complicated package and compiled code. It isn't a problem for me, so I can't make time to follow it up. hth. Keith J ------------------------------------ "baptiste auguie" <baptiste.aug...@googlemail.com> wrote in message news:aanlktikf3+higwyhwyxeu2brwqs0x9mtywz9xzyq_...@mail.gmail.com... Hi, I remember a discussion we had on this list a few months ago for a better way to decide if a point is inside a convex hull. It eventually lead to a R function in this post, http://tolstoy.newcastle.edu.au/R/e8/help/09/12/8784.html I don't know if it was included in the geometry package in the end. HTH, baptiste --------------------------------------------- inhull <- function(testpts, calpts, hull=convhulln(calpts), tol=mean(mean(abs(calpts)))*sqrt(.Machine$double.eps)) { #++++++++++++++++++++ # R implementation of the Matlab code by John D'Errico 04 Mar 2006 (Updated 30 Oct 2006) # downloaded from http://www.mathworks.com/matlabcentral/fileexchange/10226-inhull # with some modifications and simplifications # # Efficient test for points inside a convex hull in n dimensions # #% arguments: (input) #% testpts - nxp array to test, n data points, in p dimensions #% If you have many points to test, it is most efficient to #% call this function once with the entire set. #% #% calpts - mxp array of vertices of the convex hull, as used by #% convhulln. #% #% hull - (OPTIONAL) tessellation (or triangulation) generated by convhulln #% If hull is left empty or not supplied, then it will be #% generated. #% #% tol - (OPTIONAL) 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(calpts(:))) #% #% In higher dimensions, the numerical issues of floating #% point arithmetic will probably suggest a larger value #% of tol. #% # In this R implementation default tol=mean(mean(abs(calpts)))*sqrt(.Machine$double.eps) # # VALUE: Matlab returns a vector of TRUE (inside/on) or FALSE (outside) # This R implementation returns an integer vector of length n # 1 = inside hull # -1 = inside hull # 0 = on hull (to precision indicated by tol) #-------------------------------------------------------- # ensure arguments are matrices (not data frames) and get sizes calpts <- as.matrix(calpts) testpts <- as.matrix(testpts) p <- ncol(calpts) # columns in calpts cx <- nrow(testpts) # rows in testpts 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(Null(t(calpts[hull[i,-1],] - matrix(calpts[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(calpts) a <- calpts[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(testpts %*% 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 as.integer(sign(val)) } ---------------------------------------- \name{inhull} \alias{inhull} %- Also NEED an '\alias' for EACH other topic documented here. \title{ Test if one set of N-D points is inside the convex hull defined by another set } \description{ R implementation of the Matlab code by John D'Errico 04 Mar 2006 (Updated 30 Oct 2006) downloaded from http://www.mathworks.com/matlabcentral/fileexchange/10226-inhull with some modifications and simplifications } \usage{ inhull(testpts, calpts, hull = convhulln(calpts), tol = mean(mean(abs(calpts))) * sqrt(.Machine$double.eps)) } %- maybe also 'usage' for other objects documented here. \arguments{ \item{testpts}{ numeric n x p array representing n points in p dimensions; converted to matrix by \code{as.matrix(testpts); each point is tested whether inside the convex hull (defined by calpts or hull). } } \item{calpts}{ numeric m x p array representing m points in p dimensions; converted to matrix by \code{as.matrix(testpts). If 'hull' is not given, \code{geometry::convhulln(calpts)} is used to define the convex hull } } \item{hull}{ (OPTIONAL) tessellation (or triangulation) generated by convhulln. If hull is left empty or not supplied, then it will be generated. } \item{tol}{ (OPTIONAL) 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. } } \details{ %% ~~ If necessary, more details than the description above ~~ } \value{ An integer vector of length n \tabular{rl}{ 1 \tab inside hull\cr -1 \tab inside hull\cr 0 \tab on hull (to precision indicated by tol) }} \references{ \url{http://www.mathworks.com/matlabcentral/fileexchange/10226-inhull} } \author{ Keith Jewell 2009 } \note{ submitted for inclusion in geometry package~ } %% ~Make other sections like Warning with \section{Warning }{....} ~ \seealso{ %% ~~objects to See Also as \code{\link{help}}, ~~~ } \examples{ set.seed(1234) ps <- matrix(rnorm(4000),ncol=4) xs <- matrix(rnorm(1200),ncol=4) inhull(xs, ps) } % Add one or more standard keywords, see file 'KEYWORDS' in the % R documentation directory. \keyword{ ~kwd1 } \keyword{ ~kwd2 }% __ONLY ONE__ keyword per line ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
