6

I am using the rgl package to create 3D plots of my data. For some reasons (3D PCA biplots) I need vectors -- a line segment with an arrow. And I'm stuck, because I want to have 3D cones as arrow heads.

Somehow, I cannot wrap my senile mind around the geometry of the problem. Say, I would draw the vector with

segments3d( rbind( c( 0, 0, 0 ), c( 3, 3, 3 ) ) )

that is, a vector from the origin of the user coordinate system to [3,3,3].

I would like to create a cone with the tip at [3,3,3]. The base of the cone can be formed with a circle. Drawing a circle on the xz plane (perpendicular to the y plane) with radius r is easy:

n <- 10
sin.t <- sin( seq( 0, 2 * pi, len= n ) )
cos.t <- cos( seq( 0, 2 * pi, len= n ) )
r <- 0.1 
xv <- x + r * sin.t
yv <- rep( y, n )
zv <- z + r * cos.t

but how do I now transform these points such that the circle is now perpendicular to the vector? And its center 0.2 from the tip along the vectors direction? Once I have this transformation, I will draw triangles with the triangles3d function, each triangle having one corner at the tip and two vertices within the points of the circle.

This is basic maths, and I know the 18 year old me would not have a problem (or even a 28 year old me). Any hook (as opposed to fish) would be appreciated.

0

1 Answer 1

8

In the demos for rgl there is a cone3d function. It takes the base and the tip seperaately. In any event you could do something like this:

vec=rbind( c( 0, 0, 0 ), c( 3, 3, 3 ) )
segments3d( vec )


cone3d(base=vec[2,]-(vec[1,]+vec[2,]/6), 
     #this makes the head go 1/6th the length of the arrow
       rad=0.5,
       tip=vec[2,],
       col="blue",
       front="lines",
       back="lines")

Here is the cone3d function:

   cone3d <- function(base=c(0,0,0),tip=c(0,0,1),rad=1,n=30,draw.base=TRUE,qmesh=FALSE,
                    trans = par3d("userMatrix"), ...) {
   ax <- tip-base
   if (missing(trans) && !rgl.cur()) trans <- diag(4)
   ### is there a better way?
   if (ax[1]!=0) {
     p1 <- c(-ax[2]/ax[1],1,0)
     p1 <- p1/sqrt(sum(p1^2))
     if (p1[1]!=0) {
       p2 <- c(-p1[2]/p1[1],1,0)
       p2[3] <- -sum(p2*ax)
       p2 <- p2/sqrt(sum(p2^2))
     } else {
       p2 <- c(0,0,1)
     }
   } else if (ax[2]!=0) {
     p1 <- c(0,-ax[3]/ax[2],1)
     p1 <- p1/sqrt(sum(p1^2))
     if (p1[1]!=0) {
       p2 <- c(0,-p1[3]/p1[2],1)
       p2[3] <- -sum(p2*ax)
       p2 <- p2/sqrt(sum(p2^2))
     } else {
       p2 <- c(1,0,0)
     }
   } else {
     p1 <- c(0,1,0); p2 <- c(1,0,0)
   }
   degvec <- seq(0,2*pi,length=n+1)[-1]
   ecoord2 <- function(theta) {
     base+rad*(cos(theta)*p1+sin(theta)*p2)
   }
   i <- rbind(1:n,c(2:n,1),rep(n+1,n))
   v <- cbind(sapply(degvec,ecoord2),tip)
   if (qmesh) 
     ## minor kluge for quads -- draw tip twice
     i <- rbind(i,rep(n+1,n))
   if (draw.base) {
     v <- cbind(v,base)
     i.x <- rbind(c(2:n,1),1:n,rep(n+2,n))
     if (qmesh)  ## add base twice
       i.x <-  rbind(i.x,rep(n+2,n))
     i <- cbind(i,i.x)
   }
   if (qmesh) v <- rbind(v,rep(1,ncol(v))) ## homogeneous
   if (!qmesh)
     triangles3d(v[1,i],v[2,i],v[3,i],...)
   else
     return(rotate3d(qmesh3d(v,i,material=...), matrix=trans))
 }     

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.