Re: [R] mutlidimensional in.convex.hull(wasmultidimensionalpoint.in.polygon??)
Syntax suggestions implemented. Inhull's original author consulted. Function submitted for potential inclusion in geometry package. Seasons greetings to all. Keith Jewell - "baptiste auguie" wrote in message news:de4e29f50912180238m45c4b7c3g792113d5b6c4c...@mail.gmail.com... Hi, Excellent, thanks for doing this! I had tried the 2D case myself but I was put off by the fact that Octave's convhulln had a different ordering of the points to R's geometry package (an improvement to Octave's convhulln was made after it was ported to R). I'm not sure how you got around this but it's good to see that the result was worth the effort! Below are a few minor syntax suggestions, # p <- dim(calpts)[2] # columns in calpts # and other similar lines could be replaced with ncol(calpts) nrow(testpts) nrow(hull) # length(degenflag[degenflag]) # can probably be written sum(degenflag) # center = apply(calpts, 2, mean) # more efficient colMeans(calpts) Would you consider submitting this function to the maintainer of the geometry package, after checking it's OK with inhull's original author? Best regards, baptiste __ 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.
Re: [R] mutlidimensional in.convex.hull (wasmultidimensionalpoint.in.polygon??)
Hi, Excellent, thanks for doing this! I had tried the 2D case myself but I was put off by the fact that Octave's convhulln had a different ordering of the points to R's geometry package (an improvement to Octave's convhulln was made after it was ported to R). I'm not sure how you got around this but it's good to see that the result was worth the effort! Below are a few minor syntax suggestions, # p <- dim(calpts)[2] # columns in calpts # and other similar lines could be replaced with ncol(calpts) nrow(testpts) nrow(hull) # length(degenflag[degenflag]) # can probably be written sum(degenflag) # center = apply(calpts, 2, mean) # more efficient colMeans(calpts) Would you consider submitting this function to the maintainer of the geometry package, after checking it's OK with inhull's original author? Best regards, baptiste 2009/12/18 Keith Jewell : > 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(len
Re: [R] mutlidimensional in.convex.hull (wasmultidimensionalpoint.in.polygon??)
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.340.071.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.090.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@r-project.org mailing list ht