Kim Milferstedt
2007-Apr-03 15:15 UTC
[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)*sin(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)*cos(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
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. >