Calculate distances between a line and all points on an intersecting plane in r

Question

I have a regular, rectangular grid of x, y, z values (an image).
I also have an x, y, z line that intersects the grid.

For each plane (z-level) in the grid, I want to compute the Euclidean distance between each point in the grid and the intersection point of the line with that grid.

Example:

# create regular 3D array of values
vec <- array(1:21,c(21,21,21))
dimnames(vec) = list(seq(-10,10), seq(-10,10), seq(-10,10))

# convert to data.frame (with names x, y, z, value)
grid <- melt(vec, varnames=c("x","y","z"))

# and a set of points along a line
line <- data.frame(
  x = seq(-10, 10),
  y = seq(-10, 10),
  z = seq(-10, 10)
)

I tried a few things and settled on using a for loop on the z-values.

Solution:

# loop through each z-level to compute the euclidean distance between 
# all points on the plane at that level, and 
# the point on the line at that level.
tmp = data.frame()
for(i in line$z) { 
  point <- subset(line, z == i)
  plane <- subset(grid, z == i)

  plane$euclidean = (plane$x - point$x)^2 + (plane$y - point$y)^2
  if(nrow(tmp) == 0) {
    tmp = plane
  } else {
    tmp = rbind(tmp, plane)
  }
}

I wasn't particularly happy with this solution even though it was relatively quick. I didn't like the idea that I had to split and recombine the dataset resulting in a new ordering/sorting, and feel that I'm doing something wrong when I resort to a for loop in r.

I have a strong feeling that this way is somewhat inefficient and that there may be other (better?) ways to do this using one or more of the following:

1. without the for loop, subset, and rbind method.
2. using linear algebra
3. using one of the apply functions
4. using spatial data types and the function sp::spDistsN1()

Another solution was to merge the grid and line on the z-value and then perform direct calculation, but the merge step was extremely slow and the x, y, z columns in the resulting data.frame get renamed due to the duplicate column names.

---

Updates for clarification:

1. Although the value at each point (pixel) of the image is not important for the calculation of distance, it is needed later and should be brought along.
2. As in my own solution, I need the distance values assigned to each x,y,z point. So the output should include (x, y, z, value, euclidean) but doesn't necessarily need to be a data.frame.

Answer 1

Here's my tidyverse attempt, though not obvious if it's actually clearer than your approach. Conceptually I think it's basically the same: for each z calculate the distance between the line point for that z with all the other points.

1. Coerce the array to a dataframe (this way is basically the same as melt)
2. Make a function that uses dist to calculate the distances between a point and a matrix. dist actually calculates between all rows of a matrix, so we only want to keep the bottom row of the resulting triangle of distances. However, it's probably still faster than doing Euclidean distance by hand.
3. nest the data so we have one row per z, with the x and y as dataframes in a list column
4. left_join the line points on and then use pmap to apply our new function
5. unnest back out so we have columns for x, y, z points in the grid, px and py which are the line intersection point, and distance which is the distance to the point. one row per point in the grid.

vec <- array(1:21,c(21,21,21))
dimnames(vec) = list(x = seq(-10,10), y = seq(-10,10), z = seq(-10,10))

library(tidyverse)
grid <- vec %>% # same thing as melt basically
  as.tbl_cube(met_name = "value") %>%
  as_tibble() 

line <- data.frame(
  px = seq(-10, 10),
  py = seq(-10, 10),
  pz = seq(-10, 10)
)

my_dist <- function(point_x, point_y, mat){
  point_mat <- rbind(c(point_x, point_y), mat)
  dist_mat <- as.matrix(dist(point_mat))
  dist_vec <- dist_mat[nrow(dist_mat), 1:(ncol(dist_mat) - 1)]
  attributes(dist_vec) <- NULL
  return(dist_vec)
}

grid %>%
  select(-value) %>% 
  nest(x, y) %>% # One row per z
  left_join(line, by = c("z" = "pz")) %>%
  mutate(distance = pmap(list(px, py, data), my_dist)) %>%
  unnest() # Expand back out to one row per point
#> # A tibble: 9,261 x 6
#>        z    px    py distance     x     y
#>    <int> <int> <int>    <dbl> <int> <int>
#>  1   -10   -10   -10     28.3   -10   -10
#>  2   -10   -10   -10     28.3    -9   -10
#>  3   -10   -10   -10     27.6    -8   -10
#>  4   -10   -10   -10     26.9    -7   -10
#>  5   -10   -10   -10     26.2    -6   -10
#>  6   -10   -10   -10     25.6    -5   -10
#>  7   -10   -10   -10     25      -4   -10
#>  8   -10   -10   -10     24.4    -3   -10
#>  9   -10   -10   -10     23.9    -2   -10
#> 10   -10   -10   -10     23.3    -1   -10
#> # ... with 9,251 more rows

Created on 2018-09-18 by the reprex package (v0.2.0).