How to plot vectors orthogonal to quadmesh?

Question

I am trying to use the rgl package in R to plot unit vectors orthogonal to each cell in a quadmesh. I've used this page for a bit of guidance, and currently have plotted a quadmesh surface and some placeholder vectors using arrow3d, but I'd like to get these vectors to be orthogonal to the cell they're contained within and to all be the same length (e.g. length of 1). Does anyone know how to get the proper endpoint to put into arrow3d to satisfy this condition, or have a different approach?

library(raster)
library(rgl)
library(quadmesh)

m <- matrix(c(seq(0, 0.5, length = 5), 
              seq(0.375, 0, length = 4)), 3)
x <- seq(1, nrow(m)) - 0.5
y <- seq(1, ncol(m)) - 0.5
r <- raster(list(x = x, y = y, z = m))

qm <- quadmesh(r)
image(r)
op <- par(xpd = NA)
text(t(qm$vb), lab = 1:ncol(qm$vb)) #Plot index numbers for vertices

vrts<- list(c(9,10,14,13),
            c(10,11,15,14),
            c(11,12,16,15),
            c(5,6,10,9),
            c(6,7,11,10),
            c(7,8,12,11),
            c(1,2,6,5),
            c(2,3,7,6),
            c(3,4,8,7)) #Index for vertices of each raster cell starting from bottom left and moving to the right

shade3d(qm, col = "firebrick")
axes3d()
title3d(xlab="X", ylab="Y", zlab="Z")
for (i in 1:9) {
  row_number<- floor((i-1)/3)+1
  col_number<- ((i-1)%%3)+1
  p<- apply(qm$vb[1:3,  vrts[[i]]], mean, MARGIN = 1) #get current xyz position of raster cell
  rgl.spheres(x = p[1], y = p[2], z = p[3], r=0.1) #Plot points
  p0<- c(x[col_number], y[row_number],p[3]) #arrow start point
  p1<- c(x[col_number],y[row_number],p[3]+1) #arrow end point
  arrow3d(p0, p1) #Plot arrow
}

Answer 1

This is tricky. The fundamental insights are that:

1. Each square in the mesh can be thought of as a plane
2. Each plane can be defined by a formula z = ax + by + c
3. The coefficients a, b and c can be found by running a linear regression on the 4 vertices than lie at the corners of each square.
4. The vector [-a, -b, 1] is normal to the plane
5. The vector [-a, -b, 1] / sqrt(a² + b² + 1²) has length 1
6. The arrow bases are at the midpoint of each square
7. The arrow tips are at the bases of each square plus the vector [-a, -b, 1] / sqrt(a² + b² + 1²)

To implement this, we first define a couple of helper functions:

midpoints <- function(x) diff(x)/2 + x[-length(x)]

centers_from_vertices <- function(v) {
  apply(apply(v, 1, midpoints), 1, midpoints)
}

Now we can create lists of x, y, z points for the vertices and the center points like this:

vertices <- lapply(asplit(qm$vb, 1), `dim<-`, value = dim(m) + 1)
centres  <- lapply(vertices, centers_from_vertices)

Rather than having to count all the vertices and create a list belonging to each square, we can automate the process and have the indices of the vertices for each square in a list.

index <- matrix(seq(prod(dim(m) + 1)), nrow = nrow(m) + 1, byrow = TRUE)

indices <- unlist(lapply(seq(dim(m)[1]), function(i) {
  lapply(seq(dim(m)[2]), function(j) {
    index[0:1 + i, 0:1 + j]
    })
  }), recursive = FALSE)

Now we can get the end-points of each arrow using the insights above:

ends <- lapply(asplit(sapply(indices, function(i) {
  co <- coef(lm(z~x+y,  as.data.frame(lapply(vertices, c))[c(i),]))
  co <- -c(co[2], co[3], z = -1)
  co/sqrt(sum(co^2))
  }), 1), c)

ends <- Map(`+`, centres[1:3], ends)

Finally, we can draw the result:

shade3d(qm, col = "firebrick")
axes3d()
title3d(xlab="X", ylab="Y", zlab="Z")
rgl.spheres(x = centres$x, y = centres$y, z = centres$z, r = 0.1)

for(i in 1:9) {
  arrow3d(c(centres$x[i], centres$y[i], centres$z[i]),
          c(ends$x[i], ends$y[i], ends$z[i]))
}