[R] rgl how to plot a cylinder like arrow3d?

[Duncan Murdoch]

On 11-08-09 5:22 AM, René Mayer wrote:
> Dear List,
> I'm trying to draw vector in XYZ with rgl under use of a cylinder3d.
> Therefore I scale and rotate a basis-cylinder).
> However, somehow the rotation is wrong as
> verified by overplotting arrow3d().
> Where is my mistake?

I would guess it is in assuming that theta can be computed on the 
original vector, not on the rotated one:  rotating about the Y axis 
changes end[1].

I usually use a Gram-Schmidt type calculation to do this sort of thing. 
  The undocumented rgl:::GramSchmidt function does this for 3x3 
matrices, orthogonalizing by row.  So you could replace all 3 of your 
rotations (including the first one) by this:

rotation <- 
rgl:::GramSchmidt(end-start+c(1,0,0),end-start+c(0,1,0),end-start, 
order=c(3,1,2))
c <- rotate3d(c, matrix=rotation)

(I notice you also forgot to translate the cylinders to start; this 
doesn't matter in your example, but would if you had a non-zero origin.)

One problem with this kind of display is that it assumes the coordinates 
are all of equal range.  If you try plotting one of these vectors when 
the range of x, y and z are drastically different, the cylinders will 
look really strange.

Duncan Murdoch

>
> library(heplots)
> library(rgl)
>
> # ... 2 vectors
> data=data.frame(row.names=c('X','Y','Z'), x1=c(2,1,5),y=c(4,3,2))
>
> vector3D=function(start=c(0,0,0), end, mycol='green'){
> # ... cylinder as basis-vector
> c=cylinder3d(rbind(      c(start), # start
> 			 c(0,0,1)), # end
> 		         radius = 0.2,
>                            e1=cbind(0, 0, 1),
>                            e2=cbind(1, 0, 0),
> 		         sides=10
>                      )
>
> # ... rotate cylinder horizontally and scale it
> len=sqrt(sum(abs(end)*abs(end)))
> c=scale3d(c,0.1,0.1,len)
> c=rotate3d(c,-(pi/2),0,1,0)
>
> # ... rotation angles
> theta = atan2(end[2],end[1]) # angle in yx-plane for Z-achses rotation
> phi = atan2(end[3],sqrt(end[1]^2+end[1]^2)) # angle in zx for Y-achses
> rotation
>
> # ... rotation
> c=rotate3d(c,phi,0,1,0) # Y-achses rotation
> c=rotate3d(c, -theta,0,0,1) # Z-achses rotation
> shade3d(addNormals(c), col=mycol)
> }
>
> open3d()
> vector3D(start=c(0,0,0),end=data[,1], 'red');
> vector3D(start=c(0,0,0),end=data[,2],'green')
> arrow3d(c(0,0,0),data[,1],color='red');arrow3d(c(0,0,0),data[,2],
> color='green')
> axes3d(c('x','y','z'));title3d('main','sub','X','Y','Z');box3d()
>
> thanks,
> René
>
> ______________________________________________
> R-help at r-project.org 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.