---
title: "Linear transformations and matrix inverse in 3D"
author: "Michael Friendly"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Linear transformations and matrix inverse in 3D}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r nomessages, echo = FALSE}
knitr::opts_chunk$set(
  warning = FALSE,
  message = FALSE,
  echo = TRUE,
  collapse = TRUE,
  fig.height = 5,
  fig.width = 5
)
options(digits=4)
par(mar=c(4,4,1,1)+.1)
```
```{r setuprgl, echo=FALSE}
library(rgl)
library(knitr)
knit_hooks$set(webgl = hook_webgl)
```

```{r}
library(matlib)   # use the package
library(rgl)      # also using rgl
```

This vignette illustrates how linear transformations can be represented in 3D using
the [`rgl` package](https://dmurdoch.github.io/rgl/).
It is explained more fully in this StackExchange discussion
[3D geometric interpretations of matrix inverse](https://math.stackexchange.com/questions/1948184/3d-geometric-interpretations-of-matrix-inverse).
It also illustrates the determinant, $\det(\mathbf{A})$, as the volume of the transformed
cube, and the relationship between $\mathbf{A}$ and $\mathbf{A}^{-1}$.


## The unit cube
Start with a unit cube, which represents the identity matrix $\mathbf{I}_{3 \times 3}$ = `diag(3)` in 3 dimensions.
In `rgl` this is given by `cube3d()` which returns a `"mesh3d"` object containing the vertices from -1 to 1
on each axis. I translate and scale this to make each vertex go from 0 to 1.
These are represented in homogeneous coordinates to handle perspective and transformations.
See `help("identityMatrix")` for details.

```{r unit-cube}
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)
str(c3d)
```

The vertices are the 8 columns of the `vb` component of the object.
```{r vertices}
c3d$vb
```

I want to show the transformation of $\mathbf{I}$ by a matrix $\mathbf{A}$ as the corresponding transformation of the cube.

```{r matrix-A}
# matrix A: 
A <- matrix(c( 1, 0, 1, 
               0, 2, 0,  
               1, 0, 2), 3, 3) |> print()
det(A)

# A can be produced by elementary row operations on the identity matrix
I <- diag( 3 )

AA <- I |>
  rowadd(3, 1, 1) |>   # add 1 x row 3 to row 1
  rowadd(1, 3, 1) |>   # add 1 x row 1 to row 3
  rowmult(2, 2)        # multiply row 2 by 2

all(AA==A)
```

## Define some useful functions

I want to be able to display 3D mesh objects, with colored points, lines and faces.
This is a little tricky, because in rgl colors are applied to vertices, not faces. The faces are colored by interpolating the colors at the vertices. See the StackOverflow discussion
[drawing a cube with colored faces, vertex points and lines](https://stackoverflow.com/questions/39730889/rgl-drawing-a-cube-with-colored-faces-vertex-points-and-lines).
The function `draw3d()` handles this for `"mesh3d"` objects.
Another function, `vlabels()` allows for labeling vertices.

```{r draw3d}
# draw a mesh3d object with vertex points and lines
draw3d <- function(object, col=rep(rainbow(6), each=4), 
                   alpha=0.6,  
                   vertices=TRUE, 
                   lines=TRUE, 
                   reverse = FALSE,
                   ...) {
  if(reverse) col <- rev(col)
	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 $\mathbf{I}$ and $\mathbf{A}$ all together in one figure

We can show the cube representing the identity matrix $\mathbf{I}$ using `draw3d()`.
Transforming this by $\mathbf{A}$ is accomplished using `transform3d()`.

```{r 3d-demo1, webgl=TRUE}
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

```{r 3d-demo2, webgl=TRUE}
# 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}$


The inverse of `A` is found using `solve()
```{r AI}
AI <- solve(A) |> print()
```

The determinant of the matrix `A` here is 2. The determinant of $\mathbf{A}^{-1}$ is the reciprocal,
$1/ \det(\mathbf{A})$. Geometrically, this means that the larger $\mathbf{A}^{-1}$ is the smaller is
$\mathbf{A}^{-1}$.

```{r det-AI}
det(A)
det(AI)
```

You can see the relation between them by graphing both together in the same plot.
It is clear that $\mathbf{A}^{-1}$ is small in the directions $\mathbf{A}$ is large,
and vice versa.

```{r 3d-demo3, webgl=TRUE}
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

The `rgl` graphic can be animated by spinning it around an axis using `spin3d()`.
This can be played on screen using `play3d()` or
saved to a `.gif` file using `movie3d()`.

```{r 3d-demo4, webgl=TRUE}
play3d(spin3d(rpm=15),  duration=4)

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