How to geometrically project 3D xyz vectors onto a plane orthogonal to the surface normals of a triangulated surface? R

Question

I'm solving an equation with OLS which returns a timing point on an x/y/z surface mesh. I can solve for speed and directionality per the method outlined below. The problem is, my vectors are not parallel to the surface mesh and I'm uncertain of the correct transforms needed to get it in the correct direction.

The literature states:

"For the visualization of the vectors in the 3D geometry ... If only surface data were used in the calculation of the vectors, the conduction velocity vectors were geometrically projected into the plane orthogonal to the surface normals. This constraint corrected vectors pointing inside or outside the surface"

I have no idea how to do that.

Reprex Data:

library(plotly)
library(dplyr)

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

# 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
  )

Relevant Sections of the literature:

Calculating the weights is beyond the scope of MRE and I left it out.

@Chris I hope you find my updated post more accessible/reproducible. — myfatson, Jan 22 '23 at 19:19
Yes, now reproducible at home for us cut and pasters, though I get spheres rather than desired 'cones` that would seem more indicative of direction (and I'll try to jigger it around to cones). The lit is useful. Between (2.9) and (2.10), I think it is best read as 'Due to the orthogonality of the wavefront surface, the wavefront surface and propagation direction held true", as likely you did. And now back to your equations, and I removed my above comment so we don't get kicked into chat. — Chris, Jan 22 '23 at 20:48
I don't fully understand what you're trying to do but if you want to project onto an orthonormal basis of a matrix A's columns or rows you can do that with a bunch of different algorithms. My favorite is SVD which uses eigenvalues of colspace and rowspace. The gram schmidt algorithm is quite simple - if you want to post my code for it I can but you'd probably want to check complexity and that assumptions are met for your problem before relying on it. — Casey Jayne, Feb 06 '23 at 23:24
@CaseyJayne Given a 3D concave surface, a fixed known origin, and timing measurements along the entire surface, we can say for coordinate xyz, it took x seconds for some signal to get there. Given enough data, you can use OLS regression to solve for the best-fit relationship between timing and position, then using some calculus we can solve for the instantaneous velocity at any point along the surface. This returns the uvw components of a vector. What I can't do is get the arrows to points along the surface plane. I hope that helps clarify! I'm a bio-scientist and certainly beyond my depth! — myfatson, Feb 07 '23 at 21:23

How to geometrically project 3D xyz vectors onto a plane orthogonal to the surface normals of a triangulated surface? R

0 Answers0