[R] rgl: draw multiple ellipsoids

[baptiste auguie]

Dear list,

I'm trying to visualise some ellipsoidal shapes in 3D. Their position,
axes, and angular orientation can be arbitrary. I saw an ellipse3d
function in rgl; however it is heavily oriented towards the
statistical concept of ellipse of confidence, whilst I am just
concerned with the geometrical object.

Below is my current implementation. It is quite slow when creating
many shapes, so I'm asking for suggestions (I'm new to rgl),


## helping function, probably already in rgl but nevermind
euler <- function(alpha=0, beta=0, gamma=0){

  Ra <- matrix(c(cos(alpha), sin(alpha), 0, -sin(alpha), cos(alpha),
0, 0, 0, 1), ncol=3, byrow=T)
  Rb <- matrix(c(1, 0, 0, 0, cos(beta), sin(beta), 0, -sin(beta),
cos(beta)), ncol=3, byrow=T)
  Rc <- matrix(c(cos(gamma), sin(gamma), 0, -sin(gamma), cos(gamma),
0, 0, 0, 1), ncol=3, byrow=T)

  return(Rc%*%Rb%*%Ra)

}

require(rgl)

## create one ellipsoid at a given location, with axes a,b,c and Euler
angles phi, theta, psi
rgl.ellipsoid <- function (x=0,y=0,z=0, a = 1,b=1,c=1, phi=0,theta=0,psi=0,
                       subdivide = 3, smooth = TRUE, ...)
{

    sphere <- subdivision3d(cube3d(...), subdivide)
    class(sphere) <- c("mesh3d","shape3d")

    norm <- sqrt(sphere$vb[1, ]^2 + sphere$vb[2, ]^2 + sphere$vb[3,
        ]^2)
    for (i in 1:3) sphere$vb[i, ] <- sphere$vb[i, ]/norm
    sphere$vb[4, ] <- 1
    sphere$normals <- sphere$vb
    result <- scale3d(sphere, a,b,c)
    rotM <- euler(phi,theta,psi)
    result <- rotate3d(result,matrix=rotM)
    result <- translate3d(result, x,y,z)
    invisible(result)
}


## loop over the specification of a cluster
rgl.ellipsoids <- function(positions, sizes, angles,...){

  N <- NROW(positions)
  ll <- lapply(seq(1,N), function(ii)
               rgl.ellipsoid(positions[ii,1],positions[ii,2],positions[ii,3],
                         sizes[ii,1],sizes[ii,2],sizes[ii,3],
                         angles[ii,1],angles[ii,2],angles[ii,3], ...))

  shapelist3d(ll,...)

}


N <- 20
set.seed(123)
positions <- matrix(rnorm(3*N), ncol=3)
sizes <- matrix(runif(3*N, 0.01, 0.05), ncol=3)
angles <- matrix(runif(3*N, 0, 2*pi), ncol=3)

rgl.ellipsoids(positions, sizes, angles, col= 1:N)


Best regards,

baptiste

sessionInfo()
R version 2.12.1 Patched (2010-12-30 r53895)
Platform: i386-apple-darwin9.8.0 (32-bit)

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/C/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rgl_0.92.798

loaded via a namespace (and not attached):
[1] tools_2.12.1