如何绘制与四边形网格正交的向量？

Question

我正在尝试使用 R 中的 rgl 包来绘制与四边形网格中每个单元格正交的单位向量。我已经使用 this page 进行了一些指导，目前已经使用 arrow3d 绘制了一个四边形网格表面和一些占位符向量，但我想让这些向量与它们的单元格正交re contained within 和 to 都具有相同的长度（例如长度为 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

这很棘手。基本见解是：

网格中的每个方块都可以看作是一个平面
每个平面都可以用公式 z = ax + by + c
系数a、b和c可以通过运行对 4 个顶点的线性回归而不是位于每个正方形的角上。
矢量[-a,-b,1]正常飞机
向量[-a, -b, 1] / sqrt(a² + b² + 1²) 的长度为 1
箭基在每个正方形的中点
箭头尖端位于每个正方形的底部加上向量 [-a, -b, 1] / sqrt(a² + b² + 1² )

为了实现这一点，我们首先定义了几个辅助函数：

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)

现在我们可以使用上面的见解得到每个箭头的 end-points：

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:::showNormals() 函数执行此操作。这是一个内部函数，用于调试，它在每个顶点绘制网格对象中包含的法线。

将此函数用于您的数据有几个问题。首先，作为一个内部函数，它如有更改，恕不另行通知。所以我会复制它并使用它。

其次，它绘制了网格的 normals 组件，而您的网格没有。可以使用rgl::addNormals()函数添加法线，然后就可以了。例如，

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

如果你想让法线向上，你可以在绘图前说qmn$normals <- -qmn$normals。

第三个问题是它使用了丑陋的 "lines" 类型的箭头。但是如果你有函数的副本，你可以修改它。

第四个问题是箭头在顶点，而不是中心。如果你真的想要它们在中心，你需要在那里计算法线。一对向量的法线是使用叉积计算的（例如 rgl:::xprod()），但您必须决定使用哪对向量。 rgl 不强制四边形是平面的；我不知道 quadmesh 是否如此。这意味着对于四边形的法线，您可能有不止一种选择。

如何绘制与四边形网格正交的向量？

How to plot vectors orthogonal to quadmesh?

r

rgl