[R] mutlidimensional in.convex.hull (was multidimensional point.in.polygon??)

[baptiste auguie]

2009/12/10 Charles C. Berry <cberry at tajo.ucsd.edu>:
[snipped]
> Many?
>
>
>> set.seed(1234)
>> ps <- matrix(rnorm(4000),ncol=4)
>> phull <-  convhulln(ps)
>> xs <- matrix(rnorm(1200),ncol=4)
>> 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]
> +     print(length(also.outside))
> +     outside <- c(outside,also.outside)
> +     done <- length(also.outside)==0
> + }
> [1] 3
> [1] 0
>>
>
> phull2 was evaluated once, phull3 twice.
>
> Any point that is in the convex hull of rbind(ps,xs) is either in or outside
> the convex hull of ps. Right? So, just recursively eliminate points that are
> in the convex hull of the larger set.
>

If I'm not mistaken this method is efficient only because the two
point distributions are very similar (drawn from rnorm, so they look
like two concentric balls). If one of the convex hulls is very
distorted along one axis, say, I believe the method will involve many
more iterations and in the limit will require computing a convex hull
for each test point as Duncan suggested.

Such a pathological of test points example might be,

    xs <- matrix(0,ncol=4,nrow=100)
    xs[,1] <- seq(1,100)

Or did I completely miss something? (quite possible)


----

Regarding the "inhull" Matlab code, I came to the opposite conclusion:
it should be easily ported to R. 1) it is a very short piece of code
(even more so if one disregards the various checks and handling of
special cases), with no Matlab-specific objects (only integers,
booleans, matrices and vectors). 2) The core of the program relies on
the qhull library, and the same applies to R I think. 3) Matlab and R
use very similar indexing for matrices and similar linear algebra in
general.

That said, I'm a bit short on time to give it a go myself. I think the
open-source Octave could run this code too, so it might help in
checking the code step-by-step.


All the best,

baptiste