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.
To create an arrow R, we can use plot function and arrows function. We just need to understand all the coordinate values that should be passed inside the arrows function. For example, if we have two vectors that contains values from 1 to 10 then the arrow can be created by using arrows function as arrows(1,1,10,10).
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))
}
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With