dev/3d-demo.R

In friendly/matlib: Matrix Functions for Teaching and Learning Linear Algebra and Multivariate Statistics

#' ---
#' title: "3D demonstrations of linear transformations and matrix inverse"
#' author: "Michael Friendly"
#' date: "30 Sep 2016"
#' ---

#' Start with a unit cube, representing the identity matrix. Show its transformation
#' by a matrix $A$ as the corresponding transformation of the cube.
#'
#' This also illustrates the determinant, det(A), as the volume of the transformed
#' cube, and the relationship between $A$ and $A^{-1}$.


library(rgl)
library(matlib)
# cube, with each face colored differently
colors <- rep(2:7, each=4)
c3d <- cube3d()
# make it a unit cube at the origin
c3d <- scale3d(translate3d(c3d, 1, 1, 1),
               .5, .5, .5)

# matrix A: 
A <- matrix(c( 1, 0, 1, 0, 2, 0,  1, 0, 2), 3, 3)
det(A)

# same as the elementary row operations: 2*y, 2*z, x+y
I <- diag( 3 )
AA <- rowmult(rowadd(rowadd(I, 3, 1, 1), 1, 3, 1), 2, 2)
all(AA==A)

#' ## Define some useful functions

# draw a mesh3d object with vertex points and lines
# see: http://stackoverflow.com/questions/39730889/rgl-drawing-a-cube-with-colored-faces-vertex-points-and-lines
draw3d <- function(object, col=rep(rainbow(6), each=4), alpha=0.6,  vertices=TRUE, lines=TRUE, ...) {
	shade3d(object, col=col, alpha=alpha, ...)
	vert <- t(object$vb)
	indices <- object$ib
	if (vertices) points3d(vert, size=5)
	if (lines) {
	for (i in 1:ncol(indices))
		lines3d(vert[indices[,i],])
		}
}

# label vertex points
vlabels <- function(object, vertices, labels=vertices, ...) {
	text3d( t(object$vb[1:3, vertices] * 1.05), texts=labels, ...)
}


#' ## show I and A all together in one figure

open3d()
draw3d(c3d)
vlabels(c3d, c(1,2,3,5))

axes <- rbind( diag(3), -diag(3) )
rownames(axes) <- c("x", "y", "z", rep(" ", 3))
vectors3d(axes, frac.lab=1.2, headlength = 0.2, radius=1/20, lwd=3)

c3t<- transform3d(c3d, A) 
draw3d(c3t)
vlabels(c3t, c(1,2,3,5))

#' ## Same, but using separate figures, shown side by side

# NB: this scales each one separately, so can't see relative size
open3d()
mfrow3d(1,2, sharedMouse=TRUE)
draw3d(c3d)
vectors3d(axes, frac.lab=1.2, headlength = 0.2, radius=1/20, lwd=3)

next3d()	
draw3d(c3t)
vectors3d(axes, frac.lab=1.2, headlength = 0.2, radius=1/20, lwd=3)


#' ## A and A^{-1}

open3d()
draw3d(c3t)
c3Inv <- transform3d(c3d, solve(A))
draw3d(c3Inv)
vectors3d(axes, frac.lab=1.2, headlength = 0.2, radius=1/20, lwd=3)
vlabels(c3t, 8, "A", cex=1.5)
vlabels(c3t, c(2,3,5))

vlabels(c3Inv, 4, "Inv", cex=1.5)
vlabels(c3Inv, c(2,3,5))

#' Animate
play3d(spin3d(rpm=15),  duration=4)

movie3d(spin3d(rpm=15), duration=4, movie="inv-demo", dir=".")