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 的范围,因此我将其省略。
| 归档时间: |
|
| 查看次数: |
241 次 |
| 最近记录: |