如何将向量绘制到QuadMesh的正交？

Question

我正在尝试使用R中的rgl软件包来绘制QuadMesh中每个单元格的单位向量。我已经使用此页面对于一些指导，目前使用arrow3d绘制了Quadmesh Surface和某些占位符向量，但是我想让这些向量与它们内部和所有的单元格是正交的长度相同（例如长度为1）。有谁知道如何将适当的端点放入arrow3d以满足此情况或具有不同的方法中？

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

这很棘手。基本的见解是：

1. 网格中的每个正方形可以被视为平面
2. 每个平面，每个平面都可以由公式z = a a x + b b y + c
3. 系数 a ， b 和 c 可以通过在4个顶点上进行线性回归而不是躺在每个正方形的角落。
4. 向量[ -a ， -b ， 1 ]与矢量
5. [ -a ， -b ， 1 ]/sqrt（a² + b² + + 1²）长度1
6. 箭头底座位于每个正方形的中点，
7. 箭头尖端位于每个正方形 plus 的底部[ -a -a ， -b ， 1 ]/sqrt（a² + b² + 1²）

要实现此目标，我们首先定义A几个辅助功能：

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

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

现在我们可以创建顶点和类似中心点的x，y，z点的列表：

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

而不是必须计数所有顶点并创建一个属于每个正方形的列表，我们可以自动化该过程和过程具有列表中每个正方形的顶点索引。

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)

现在，我们可以使用上面的见解来获取每个箭头的终点：

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)

最后，我们可以绘制结果：

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

rgl ::: showmals（）函数可以做到这一点。它是一个内部功能，用于调试，并绘制每个顶点的网格对象中包含的正态。

在您的数据中使用此功能有一些问题。首先，作为内部功能，它需要更改，恕不另行通知。因此，我会复制它并与之合作。

其次，它绘制了网格的formals组件，而您的网格没有一个。您可以使用rgl :: addNormals（）函数添加正常功能，然后它将起作用。例如，

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

如果您希望正态指向，则可以在绘图之前说qmn $ normals＆lt; - -qmn $ normals。

第三个问题是它使用丑陋的“ lines”箭头类型。但是，如果您有该功能的副本，则可以对其进行修改。

第四个问题是箭头在顶点，而不是中心。如果您真的希望它们在中心，则需要计算那里的正态。使用交叉产品（例如rgl ::: xprod（））计算出一对矢量的正态，但是您必须决定要使用哪对矢量。 rgl不迫使四方面为平面；我不知道QuadMesh是否确实如此。这意味着您可能有多个选择的正常选择。