[R] boundary check

[Keith Jewell]

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.auguie at googlemail.com> wrote in message 
news:AANLkTikf3+higWyHWyXEU2bRWqS0x9mtYWz9xzyQ_OLW at 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