如何将 3D xyz 矢量几何投影到与三角表面的表面法线正交的平面上?右

myf*_*son 5 r mesh vector-graphics triangulation

我正在使用 OLS 求解方程,该方程返回 x/y/z 表面网格上的计时点。我可以按照下面概述的方法求解速度和方向性。问题是,我的矢量与表面网格不平行,并且我不确定使其朝正确方向所需的正确变换。

文献指出:

“为了在 3D 几何中实现矢量的可视化......如果在矢量计算中仅使用表面数据,则传导速度矢量会以几何方式投影到与表面法线正交的平面中。此约束修正了指向内部或内部的矢量。表面之外”

我不知道该怎么做。

代表数据:

library(plotly)
library(dplyr)
Run Code Online (Sandbox Code Playgroud)
tri_list <- list(surf_tri = structure(c(1L, 1L, 4L, 8L, 11L, 12L, 17L, 6L, 
19L, 13L, 3L, 5L, 22L, 1L, 26L, 18L, 18L, 8L, 21L, 14L, 27L, 
4L, 14L, 7L, 6L, 29L, 22L, 16L, 23L, 26L, 3L, 28L, 2L, 13L, 12L, 
31L, 33L, 35L, 32L, 28L, 33L, 31L, 20L, 9L, 36L, 12L, 29L, 20L, 
31L, 27L, 17L, 26L, 8L, 15L, 9L, 29L, 9L, 35L, 15L, 11L, 16L, 
16L, 38L, 31L, 37L, 32L, 24L, 35L, 23L, 23L, 23L, 37L, 10L, 35L, 
2L, 4L, 5L, 7L, 10L, 13L, 16L, 18L, 20L, 2L, 2L, 7L, 3L, 3L, 
25L, 17L, 27L, 17L, 7L, 2L, 18L, 27L, 28L, 5L, 7L, 1L, 30L, 29L, 
30L, 31L, 13L, 5L, 5L, 12L, 14L, 33L, 10L, 25L, 25L, 20L, 12L, 
12L, 28L, 36L, 15L, 19L, 4L, 7L, 37L, 4L, 27L, 13L, 36L, 36L, 
7L, 16L, 10L, 32L, 16L, 15L, 23L, 30L, 16L, 34L, 34L, 37L, 35L, 
24L, 13L, 24L, 26L, 11L, 11L, 38L, 3L, 2L, 2L, 6L, 9L, 14L, 15L, 
8L, 9L, 14L, 13L, 21L, 23L, 22L, 24L, 8L, 17L, 15L, 20L, 28L, 
6L, 6L, 19L, 4L, 4L, 22L, 16L, 22L, 22L, 32L, 23L, 21L, 28L, 
31L, 19L, 34L, 34L, 32L, 26L, 19L, 9L, 33L, 21L, 11L, 11L, 9L, 
1L, 9L, 32L, 29L, 29L, 31L, 9L, 8L, 8L, 17L, 33L, 38L, 38L, 38L, 
39L, 23L, 39L, 37L, 10L, 38L, 39L, 25L, 26L, 39L, 24L, 38L, 37L, 
39L), dim = c(74L, 3L), dimnames = list(NULL, c("V1", "V2", "V3"
))), coordinates = structure(c(1.9194540977478, 11.1435489654541, 
5.42491292953491, -1.66259050369263, 7.44240760803223, -6.69524192810059, 
-0.0333123430609703, -4.23227882385254, 21.1540908813477, 28.65553855896, 
26.0386810302734, 24.6995525360107, 38.495361328125, 22.120174407959, 
-2.17373514175415, 0.292593091726303, -11.4042282104492, -10.3343534469604, 
17.5129470825195, 11.7786016464233, 9.42895698547363, -1.22890889644623, 
23.1779632568359, 31.624174118042, 38.5983734130859, 45.1263847351074, 
-11.1577320098877, 12.0059661865234, -2.91624927520752, 0.309472501277924, 
30.9889793395996, 42.6134490966797, 23.7120399475098, 26.9033088684082, 
30.3514823913574, 7.96694564819336, 31.1637763977051, 31.1527690887451, 
25.4455451965332, -61.4243927001953, -56.5476036071777, -45.9739837646484, 
-66.2752914428711, -66.6171951293945, -84.5177154541016, -78.7427749633789, 
-92.5044403076172, -86.8224487304688, -90.4450073242188, -100.261276245117, 
-64.494514465332, -54.8481636047363, -57.9494781494141, -109.206504821777, 
-76.7904586791992, -95.8042831420898, -88.341796875, -64.2513122558594, 
-69.0683898925781, -68.2359237670898, -69.7075653076172, -76.6936798095703, 
-76.4758529663086, -78.4939727783203, -64.2876815795898, -86.8634643554688, 
-61.3102798461914, -68.8540496826172, -74.4825897216797, -79.063835144043, 
-81.4954528808594, -86.3758010864258, -87.7453536987305, -83.0354080200195, 
-109.458961486816, -90.4418487548828, -90.170280456543, -80.8295745849609, 
189.259994506836, 207.376663208008, 187.388595581055, 200.771896362305, 
215.514175415039, 214.221862792969, 217.531005859375, 215.624313354492, 
211.211242675781, 205.350524902344, 203.870376586914, 210.061508178711, 
195.598449707031, 210.308349609375, 195.291015625, 177.159469604492, 
194.172241210938, 210.24040222168, 213.330307006836, 216.293930053711, 
216.293441772461, 180.884536743164, 163.47248840332, 163.422592163086, 
168.801071166992, 186.711502075195, 207.275924682617, 212.505798339844, 
190.655349731445, 176.77978515625, 204.541122436523, 191.065948486328, 
209.147216796875, 206.617050170898, 172.961532592773, 211.9609375, 
203.223968505859, 179.055145263672, 170.081588745117), dim = c(39L, 
3L)), activation_timing_ms = c(-25.9116878216266, -35, 17.8156572141554, 
-32.3085435482249, -30.425790486056, -22.7921350518616, -19.7547757364543, 
-17.1422476398632, -5.21201605679312, 6.529324221754, 6.37415402817896, 
-18.6397688908753, 6.16765189118496, -23.9073664920415, -6.74190321661786, 
-17, -15.3633306040704, -20.0486637091315, -23.098251364298, 
-23.5797990311573, -23.6085020691451, -27.8322123439857, 11.4184099588406, 
12.5422915301849, 9.07414737038971, 5.42139612786991, -21.8214591939045, 
-29.3013020135356, -33.236757970705, -19.3552220222891, -5.4698807256143, 
-2.27373675443232e-13, -3.44126922278861, -0.401211436169433, 
17.9550207172763, 1, 5.84672565064739, -2.04852001824679, 26.5347286690835
))
Run Code Online (Sandbox Code Playgroud)
# Polynomial Surface Fitting 

lm_dat <- cbind(tri_list$coordinates, activation_timing_ms = tri_list$activation_timing_ms) |>
  as.data.frame() |>
  rename(x = V1, y = V2, z = V3)

lm <- lm(activation_timing_ms~polym(x,y,z, degree=2, raw = TRUE), lm_dat)

lm_coef <- coef(lm)
a <- lm_coef[["polym(x, y, z, degree = 2, raw = TRUE)2.0.0"]]
b <- lm_coef[["polym(x, y, z, degree = 2, raw = TRUE)0.2.0"]]
c <- lm_coef[["polym(x, y, z, degree = 2, raw = TRUE)0.0.2"]]
d <- lm_coef[["polym(x, y, z, degree = 2, raw = TRUE)1.1.0"]]
e <- lm_coef[["polym(x, y, z, degree = 2, raw = TRUE)1.0.1"]]
f <- lm_coef[["polym(x, y, z, degree = 2, raw = TRUE)0.1.1"]]
g <- lm_coef[["polym(x, y, z, degree = 2, raw = TRUE)1.0.0"]]
h <- lm_coef[["polym(x, y, z, degree = 2, raw = TRUE)0.1.0"]]
i <- lm_coef[["polym(x, y, z, degree = 2, raw = TRUE)0.0.1"]]
j <- lm_coef[["(Intercept)"]]

# Differentiation
t_x_fun <- function(x,y,z){
  2*a*x + d*y + e*z + g
}

t_y_fun <- function(x,y,z){
  2*b*y + d*x + f*z + h
}

t_z_fun <- function(x,y,z){
  2*c*z + e*x + f*y + i
}

CV_x_fun <- function(t_x,t_y,t_z){
  t_x/(t_x^2 + t_y^2 + t_z^2)
}

CV_y_fun <- function(t_x,t_y,t_z){
  t_y/(t_x^2 + t_y^2 + t_z^2)
}

CV_z_fun <- function(t_x,t_y,t_z){
  t_z/(t_x^2 + t_y^2 + t_z^2)
}

fitted_activation_timing <- fitted(lm)

# Loop through the V1 coordinates

node_indices_tmp <- vector(length=nrow(tri_list$surf_tri))
v1_vertex_coord <- vector("list", 
                          length = nrow(tri_list$coordinates))

for(i in seq_along(1:nrow(tri_list$coordinates))) {
  
  v3_index_tmp <- tri_list$surf_tri[i,3]
  v2_index_tmp <- tri_list$surf_tri[i,2]
  v1_index_tmp <- tri_list$surf_tri[i,1]
  
  timing_sorted <- matrix(c(fitted_activation_timing[v1_index_tmp], v1_index_tmp,
                            fitted_activation_timing[v2_index_tmp], v2_index_tmp,
                            fitted_activation_timing[v3_index_tmp], v3_index_tmp),
                          nrow = 3,
                          ncol = 2,
                          byrow=T)
  timing_sorted <- timing_sorted[order(timing_sorted[, 1]), ]
  
  v3_index <- timing_sorted[3,2]
  v2_index <- timing_sorted[2,2]
  v1_index <- timing_sorted[1,2]
  
  v1_vertex_coord[[i]] <- tri_list$coordinates[v1_index,]
}

v1_vertex_coord <- matrix(unlist(v1_vertex_coord),
                          nrow = length(v1_vertex_coord),
                          ncol = 3, 
                          byrow=T)

CV_x <- vector(length=nrow(v1_vertex_coord))
CV_y <- vector(length=nrow(v1_vertex_coord))
CV_z <- vector(length=nrow(v1_vertex_coord))

for(i in seq_along(1:nrow(v1_vertex_coord))) {
  x <- v1_vertex_coord[i,1]
  y <- v1_vertex_coord[i,2]
  z <- v1_vertex_coord[i,3]
  
  t_x <- t_x_fun(x,y,z)
  t_y <- t_y_fun(x,y,z)
  t_z <- t_z_fun(x,y,z)

  CV_x[i] <- CV_x_fun(t_x,t_y,t_z)
  CV_y[i] <- CV_y_fun(t_x,t_y,t_z)
  CV_z[i] <- CV_z_fun(t_x,t_y,t_z)
}

CV <- cbind(lm_dat, CV_x, CV_y, CV_z, fitted_activation_timing) |> as.data.frame()

plot_ly() |>
    add_trace(
      name = "Faces",
      type = "mesh3d",
    x = tri_list$coordinates[,1], y = tri_list$coordinates[,2], z = tri_list$coordinates[,3],
    i = tri_list$surf_tri[,1]-1, j = tri_list$surf_tri[,2]-1, k = tri_list$surf_tri[,3]-1,
    opacity = 0.8,
    flatshading = TRUE # we don't want smoothing
    ) |>
  add_trace(
    type = "scatter3d",
    mode = "markers",
    x = tri_list$coordinates[,1],
    y = tri_list$coordinates[,2],
    z = tri_list$coordinates[,3]) |>
  add_trace(
    type="cone",
    name = "CV Vectors",
    x = CV[,"x"],
    y = CV[,"y"],
    z = CV[,"z"],
    u = CV[,"CV_x"],
    v = CV[,"CV_y"],
    w = CV[,"CV_z"],
    opacity=0.8,
    sizemode= 'absolute',
    sizeref= 1
  ) 
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述

文献的相关部分: 在此输入图像描述 在此输入图像描述 在此输入图像描述

计算权重超出了 MRE 的范围,因此我将其省略。