Re: [R] mutlidimensional in.convex.hull (wasmultidimensionalpoint.in.polygon??)

From: Keith Jewell <k.jewell_at_campden.co.uk>
Date: Fri, 18 Dec 2009 10:15:51 +0000

Hi All,

I couldn't resist doing this the "right" way! A colleague explained the vector algebra to me (thanks Martin!) and I've followed the structure of the Matlab code referenced below, but all errors are mine!

I don't pretend to great R expertise, so this code may not be optimal (in time or the memory issue addressed by the Matlab code), but it is a lot faster than the other algorithm discussed in this thread. These timings are on the real example data described in my previous message in this thread.

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

   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.

list of date sections of archive