Interpolate points along a 3d line

Question

Although I outline the problem in R, answers in R and Python are welcome.

Suppose we have a set of points in x,y,z that define a discrete path (or set of connected line segments) over a vector of time t. These paths are not monotonic in z.

path <- data.frame(x = c(245, 233, 270, 400, 380),
                   y = c(245, 270, 138, 225, 300),
                   z = c(0, 1.2, 5, 3, 9),
                   t = 1:5)


plot3D::scatter3D(path$x, path$y, path$z, type = "b", bty = "g", 
                  phi = 0, col = "red", pch = 20, ticktype = "detailed")

How can I interpolate along the path at some arbitrary resolution in z?

For example, say I wanted to retain a resolution in z of 1 unit. The points along z are 0, 1.2, 5, 3, 9. Therefore, one possible solution that satisfies this constraint is to interpolate at 1, 2, 3, 4, 4, 4, 5, 6, 7, 8 in the z direction, yielding the blue points in the figure below (labels indicate z location):

Ultimately, I would like the coordinates of the blue dots. We can brute force the solution by sequentially solving for the equation of the line in 3d for each pair of points z, and then interpolating along each line segment. However, I want to make sure I'm not missing some existing implementation or clever hack.

Answer 1

I used purrr::map just to tie it all together into a single tibble, but you could easily do them independently. The only problem is if your x or y is NOT monotonic over z. Then you'll need to do a fourth variable indicating row order, and it gets even more complicated from there.

library(dplyr)
library(purrr)
library(plotly)

path3d <- data.frame(
  x = c(245, 233, 270, 400, 380),
  y = c(245, 270, 138, 225, 300),
  z = c(0, 1.2, 3, 5, 9)
  )

path3d_interp <- list(
  x = approxfun(path3d$z, path3d$x),
  y = approxfun(path3d$z, path3d$y)
) %>% 
  map(~.x(1:8)) %>% as_tibble() %>% 
  mutate(z = 1:8)

      x     y     z
  <dbl> <dbl> <int>
1  235   266.     1
2  249.  211.     2
3  270   138      3
4  335   182.     4
5  400   225      5
6  395   244.     6
7  390   262.     7
8  385   281.     8

plot_ly(path3d) %>% 
  add_paths(x = ~x, y = ~y, z = ~z) %>% 
  add_markers(
    x = ~x, y = ~y, z = ~z,
    data = path3d_interp
    )

---

Updated for non-monotonic in z:

I can't think of an especially elegant way to predict x and y exactly if there isn't exactly 1 solution for each desired z. The general outline I'd suggest is creating a funX, funY, and funZ of each variable, as predicted by your time t. Then using a very high-resolution vector of new t values and subsetting it for funZ(new_t_values). You'll never get exactly to the values you're looking for, but you can approximate them to the desired arbitrary precision:

path3d <- data.frame(
  x = c(245, 233, 270, 400, 380),
  y = c(245, 270, 138, 225, 300),
  z = c(0, 1.2, 5, 3, 9),
  t = 1:5
  )

Just to get a good visual of what's going on here with respect to t:

library(ggplot2)
ggplot(path3d) + 
  geom_path(aes(t, x), color = "blue") +
  geom_path(aes(t, y), color = "red") +
  geom_path(aes(t, z*50), color = "orange") +
  labs(y = "x, y, z*50")

This is looping over each column of path3d (x, y, z, and t) and creating a linear interpolation function for each variable against t as a predictor.

path3d_interp_funs <- 
  map(path3d, ~approxfun(path3d$t, .x))

---

Inspect the full range of t:

Now we can make a high resolution vector across the entire range of t, here 1 million elements. You can increase this as much as your precision requires and your memory allows.

new_t_values <- seq(min(path3d$t), max(path3d$t), length.out = 1e6)

 1.000000 1.000004 1.000008 1.000012 1.000016 1.000020 ...

Generate the full path of z:

Now we see what the z values are for each possible t in the range.

z_candidates <- path3d_interp_funs$z(new_t_values)

0.000000e+00 4.800005e-06 9.600010e-06 1.440001e-05 1.920002e-05 2.400002e-05 ...

Test where approximate-z is closest to your desired z:

So now we take each desired value of z (1:8), and ask which element of the z_candidates vector has the smallest absolute deviation from it. This returns indices we can use to subset new_t_values.

t_indices <- map_dbl(1:8, ~which.min(abs(z_candidates-.x)))

208334 302632 750000 434211 500001 875000 916667 958333

A sanity check: do those chosen t values result in the zs you want?

path3d_interp_funs$z(new_t_values[t_indices])

0.9999994 1.9999958 3.0000020 3.9999986 4.9999960 5.9999970 7.0000060 7.9999910

---

Generate your desired data:

So let's loop over each of those approximating functions, evaluating each one at our new chosen values of t:

path3d_interp <- 
  path3d_interp_funs %>%
  map(~.x(new_t_values[t_indices])) %>% 
  as_tibble()

# A tibble: 8 x 4
      x     y     z     t
  <dbl> <dbl> <dbl> <dbl>
1  235.  266. 1.000  1.83
2  241.  242. 2.00   2.21
3  400.  225. 3.00   4.00
4  260.  173. 4.00   2.74
5  270.  138. 5.00   3.00
6  390.  262. 6.00   4.50
7  387.  275. 7.00   4.67
8  383.  287. 8.00   4.83

---

Visualize your results

You can check to make sure the points do indeed fall on the correct paths:

ggplot(path3d) + 
  geom_path(aes(t, x), color = "blue") +
  geom_path(aes(t, y), color = "red") +
  geom_path(aes(t, z*50), color = "orange") +
  geom_point(data = path3d_interp, aes(t, x), color = "blue") +
  geom_point(data = path3d_interp, aes(t, y), color = "red") +
  geom_point(data = path3d_interp, aes(t, z*50), color = "orange") +
  geom_text(data = path3d_interp, aes(t, z*50, label = round(z))) +
  labs(y = "x, y, z*50")

And see them in 3D:

plot_ly(path3d) %>% 
  add_paths(x = ~x, y = ~y, z = ~z) %>% 
  add_markers(
    x = ~x, y = ~y, z = ~z,
    data = path3d_interp
    )  %>% 
  add_text(
    x = ~x, y = ~y, z = ~z, text = ~round(z),
    data = path3d_interp
  )

---

Updated again:

At the risk of being way too long-winded, I was reminded of the uniroot function, which does a pretty good job of this backsolving:

t_solutions <- map(1:8,
  ~uniroot(
    function(x) path3d_interp_funs$z(x) - .x, 
    interval = range(path3d$t)
  )
) %>% map_dbl("root") 

  1.833333 2.210526 2.473684 2.736842 4.333333 4.500000 4.666667 4.833333

However, you may notice that these solutions aren't the same as the ones from the previous method!

uniroot found solutions closer to the extremes of the interval, and not in the local region where the function changes direction. But that brings up the issue that there may be multiple t values for each of your desired z values. So a more robust solution:

root_finder <- function(f, zero, range, cuts) {

  endpts <- seq(range[1], range[2], length.out = cuts+1)

  range_list <- map2(endpts[-(cuts+1)], endpts[-1], c)

  safe_root <- possibly(uniroot, otherwise = NULL)
  f0 <- function(x) f(x) - zero

  map(range_list, ~safe_root(f0, interval = .x, maxiter = 100)$root) %>% 
    compact() %>% 
    unlist() %>% 
    unique()

}

This function takes a function, a new desired "zero" for uniroot to solve for, a range of values to test the function over, and how many buckets to divide that range into. Then it tests for a solution in each bucket, and returns NULL if none is found. It then throws out the NULLs and removes any duplicates (if a solution is exactly at the border of a bucket, for example).

root_finder(path3d_interp_funs$z, zero = 4, range = range(path3d$t), cuts = 10)

2.736842 3.500000 4.166667

Then you can loop over all your desired z values to find the value(s) of t that satisfy it.

t_solutions <- map(
    1:8, 
    ~root_finder(path3d_interp_funs$z, zero = .x, range = range(path3d$t), cuts = 100)
  ) %>% unlist()

1.833333 2.210526 2.473684 4.000000 2.736842 3.500000 4.166667 3.000000 4.333333 4.500000 4.666667 4.833333

Then again passing these t values into each of the functions you created earlier, you can make a dataframe of all of them.

map(path3d_interp_funs, ~.x(t_solutions)) %>% 
  as_tibble()

       x     y     z     t
   <dbl> <dbl> <dbl> <dbl>
 1  235   266.  1.    1.83
 2  241.  242.  2.    2.21
 3  251.  207.  3.    2.47
 4  400   225   3     4   
 5  260.  173.  4     2.74
 6  335   182.  4     3.5 
 7  397.  238.  4.    4.17
 8  270   138   5     3   
 9  393.  250.  5.00  4.33
10  390   262.  6     4.5 
11  387.  275   7.    4.67
12  383.  288.  8.00  4.83