[R] orthographic projection of ellipsoids

[baptiste auguie]

Dear list,


I've generated a list of 3D coordinates representing ellipsoids in  
arbitrary orientations. I'm now trying to obtain a 2D projection of  
the scene, that is to draw the silhouette of each object on a plane  
(x,y). The only way I could think of is to compute the convex hull of  
the (x,y) coordinates of each object and use this as the outline of  
the object. This is clearly not very efficient or satisfying.

I think I'm on the wrong track from the start. Is there an obvious  
analytical parametrisation of such projections? Any comments are  
welcome.

Many thanks,

baptiste

>
> rotM3d <- function(theta=0, phi=0, psi=0){ # 3D rotation matrix
> 	a11 <- cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi)
> 	a12 <- cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi)
> 	a13 <- sin(psi)*sin(theta)
> 	a21 <- -sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi)
> 	a22 <- -sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi)
> 	a23 <- cos(psi)*sin(theta)
> 	a31 <- sin(theta)*sin(phi)
> 	a32 <- -sin(theta)*cos(phi)
> 	a33 <- cos(theta)
> 	matrix(c(a11, a12, a13, a21, a22, a23, a31, a32, a33), ncol=3)
> }
> rotM3d() # I
>
> ellipsoid <-  # idea borrowed from a post in the R-mailing list  
> (John Fox i think)
> function(x=0, y=0, z=0, radius=1, shape=diag(c(10, 2, 2)),theta=0,  
> phi=0, psi=0,  segments=11) {
> 			 angles <- (0:segments)*2*pi/segments
> 			ecoord2 <- function(p) {
> 			    c(cos(p[1])*sin(p[2]), sin(p[1])*sin(p[2]), cos(p[2]))
> 			}
> unit.sphere <- t(apply(expand.grid(angles, angles), 1, ecoord2))
> xyz <- t(c(x, y, z) + radius * rotM3d(theta, phi, psi)%* 
> %t(unit.sphere %*% chol(shape)))
> chull(x=xyz[, 1], y=xyz[, 2])->points
> mdf <- data.frame(x=xyz[points, 1], y=xyz[points, 2])
> polygon(mdf, col=hcl(h = 0, c = 35, l = 85, 0.5))
> invisible(xyz)
> }
>
>
> xx <- seq(-5, 5, len=10)
> xy <- expand.grid(xx, xx)
>
> xy.jit <- apply(xy, 2, jitter, amount=0.4)
>
> par(mar=c(0, 0, 0, 0))
> plot(xy.jit, t="n", axes=F, xlab="", ylab="")
>
> x <- xy.jit[, 1]
> y <- xy.jit[, 2]
>
> twist <- pi*y/max(abs(y)) * rep(1, length(y))
> tilt <- pi*x/max(abs(x)) * rep(1, length(x))
> b.quiet <- mapply(ellipsoid,
> 	theta=twist, psi=tilt,x=x, y=y,  SIMPLIFY=F, radius=0.15)

_____________________________

Baptiste Auguié

School of Physics
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK

Phone: +44 1392 264187

http://newton.ex.ac.uk/research/emag