[R] which points within an ellipsoid? Sorting data in 3d

[Greg Snow]

I don't understand how the x,y,z vectors define the ellipsoid, but one
approach would be:

Find a 3x3 matrix that represents the ellipsoid as a variance matrix
(some of the code in the ellipse package may be informative for this).

Subtract x.coord.point, y.coord.point, and z.coord.point from your data.

Post multiply the data (in an nx3 matrix) with the inverse of the
cholesky decomposition of the above variance matrix (this should convert
your data on the scale of the ellipsoid to be now based on the unit
sphere).

Find x^2+y^2+z^2 (where x,y,z are the transformed coordinates of your
points) which will be the distance of each point from 0, any values < 1
correspond to points originally within the ellipsoid, >1 outside the
ellipsoid.

Hope this helps,

-- 
Gregory (Greg) L. Snow Ph.D.
Statistical Data Center
Intermountain Healthcare
greg.snow at intermountainmail.org
(801) 408-8111
 
 

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch 
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Kim Milferstedt
> Sent: Tuesday, April 03, 2007 9:16 AM
> To: r-help at stat.math.ethz.ch
> Subject: [R] which points within an ellipsoid? Sorting data in 3d
> 
> Hello,
> 
> in a three dimensional coordinate system, I'd like to find 
> all my experimental data points that fall within an ellipsoid 
> around a fixed coordinate. The fixed point is defined by 
> (x.coord.point, y.coord.point, z.coord.point). The 
> coordinates of the ellipsoid are given by the three vectors x,y,z.
> 
> In a previous version of my code, I simply used a box instead 
> of an ellipsoid to sort my data, which was really easy as I 
> only had to compare each data point to three coordinates. I 
> used something like this...
> 
> XYZ.within <- which( x.data  < x.coord.point + multipl &
>                               x.data  > x.coord.point - multipl &
>                               y.data  < y.coord.point + multipl &
>                               y.data  > y.coord.point - multipl &
>                         z.data  < z.coord.point + multipl &
>                               z.data  > z.coord.point - multipl
>                  )
> 
> But now I have many more coordinate sets to compare my data 
> to. How can I find out in an efficient way which of the data 
> points lie within the ellipsoid?
> 
> Thanks!
> 
> Kim
> 
> #mock-data for the display with rgl
> data.1 <- xyz.coords(5,-2,17)
> data.2 <- xyz.coords(15, -10,18)
> data.3 <- xyz.coords(-19, 13,9)
> 
> #code I use to construct the ellipsoid
> x.coord.point <- 4
> y.coord.point <- -7
> z.coord.point <- -3
> 
> radius <- 8
> x.radius.multiplier <- 1
> y.radius.multiplier <- 2
> z.radius.multiplier <- 3
> 
> endpoint <- 350
> interval <- 2
> alpha <- rep(seq(0,endpoint, by
> =interval),rep((endpoint+interval)/interval,(endpoint+interval
> )/interval))
> beta <- c(rep(seq(0,endpoint, by 
> =interval),(endpoint+interval)/interval))
> 
> x <-
> x.coord.point+(x.radius.multiplier*radius)*cos(alpha*pi/180)*s
> in(beta*pi/180)
> y <- y.coord.point+(y.radius.multiplier*radius)*sin(alpha*pi/180)
> z <-
> z.coord.point+(z.radius.multiplier*radius)*cos(alpha*pi/180)*c
> os(beta*pi/180)
> 
> # using rlg to visualize the mock-data and the ellipsoid 
> plot.lim <- c(-20,20) plot3d(x, y, z,
>          xlim = plot.lim,
>          ylim = plot.lim,
>          zlim = plot.lim
>          )
> points3d(data.1, col ="green", size = 6) points3d(data.2, col 
> ="red", size = 6) points3d(data.3, col ="blue", size = 6)
> 
> __________________________________________
> 
> Kim Milferstedt
> University of Illinois at Urbana-Champaign Department of 
> Civil and Environmental Engineering
> 4125 Newmark Civil Engineering Laboratory
> 205 North Mathews Avenue MC-250
> Urbana, IL 61801
> USA
> phone: (001) 217 333-9663
> fax: (001) 217 333-6968
> email: milferst at uiuc.edu
> http://cee.uiuc.edu/research/morgenroth
> 
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>