From: Keith Jewell <k.jewell_at_campden.co.uk>

Date: Fri, 18 Dec 2009 10:15:51 +0000

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)

*# DEFAULT: tol = 1e-6
*

*#
*

*# 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)

#--------------------------------------------------------

R-help_at_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. Received on Fri 18 Dec 2009 - 10:19:48 GMT

Date: Fri, 18 Dec 2009 10:15:51 +0000

> system.time(inhull(xs,ps)) # the "right" way

user system elapsed

1.34 0.07 1.41

> system.time({phull <- convhulln(ps) # the other algorithm

+ phull2 <- convhulln(rbind(ps,xs)) + nrp <- nrow(ps) + nrx <- nrow(xs) + outside <- unique(phull2[phull2>nrp])-nrp + done <- FALSE + while(!done){ + phull3 <- convhulln(rbind(ps,xs[-(outside),])) + also.outside <- (1:nrx)[-outside][unique(phull3[phull3>nrp])-nrp] + outside <- c(outside,also.outside) + done <- length(also.outside)==0 + }})

user system elapsed

15.09 0.09 15.22

Although I really must move on now, if anyone has comments, criticisms, suggestions for improvement I would be interested.

Enjoy!

Keith Jewell

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

# In this R implementation default

tol=mean(mean(abs(calpts)))*sqrt(.Machine$double.eps)

# 0 = on hull (to precision indicated by tol)

#--------------------------------------------------------

require(geometry, quietly=TRUE) # for convhulln
require(MASS, quietly=TRUE) # for Null

# ensure arguments are matrices (not data frames) and get sizes

calpts <- as.matrix(calpts)

testpts <- as.matrix(testpts)

p <- dim(calpts)[2] # columns in calpts
cx <- dim(testpts)[1] # rows in testpts
nt <- dim(hull)[1] # 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(length(degenflag[degenflag]) > 0) warning(length(degenflag[degenflag])," degenerate faces in convex hull")

nrmls <- nrmls[!degenflag,]

nt <- dim(nrmls)[1]

# find center point in hull, and any (1st) point in the plane of each

simplex

center = apply(calpts, 2, mean)

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))

}

R-help_at_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. Received on Fri 18 Dec 2009 - 10:19:48 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Fri 18 Dec 2009 - 11:30:30 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*