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[phull2nrp])-nrp
+ done - FALSE
+ while(!done){
+ phull3 - convhulln(rbind(ps,xs[-(outside),]))
+ also.outside - (1:nrx)[-outside][unique(phull3[phull3nrp])-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
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 k.jew...@campden.co.uk:
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[phull2nrp])-nrp
+ done - FALSE
+ while(!done){
+ phull3 - convhulln(rbind(ps,xs[-(outside),]))
+ also.outside - (1:nrx)[-outside][unique(phull3[phull3nrp])-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]
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 baptiste.aug...@googlemail.com 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.
