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]))
}

Answer 2

The rgl:::showNormals() function does this. It's an internal function, used for debugging, and it plots the normals that are included in the mesh object, at each vertex.

There are a couple of issues for using this function with your data. First, being an internal function, it's subject to change without notice. So I'd make a copy of it and work with that.

Second, it plots the normals component of the mesh, and your mesh doesn't have one. You can use the rgl::addNormals() function to add normals, and then it will work. For example,

qmn <- addNormals(qm)
rgl:::showNormals(qmn)

If you want the normals to point up, you can say qmn$normals <- -qmn$normals before plotting.

The third issue is that it uses the ugly "lines" type of arrow. But if you've got a copy of the function, you could modify that.

The fourth issue is that the arrows are at the vertices, not the centers. If you really want them at the centers, you'll need to calculate normals there. Normals to a pair of vectors are calculated using the cross product (e.g. rgl:::xprod()), but you'll have to decide which pair of vectors to use. rgl doesn't force quads to be planar; I don't know if quadmesh does. This means that you might have more than one choice for the normal to the quad.