Catmull-Rom Splines in R

# based on  https://github.com/vz-risk/VCDB/commits/master

tj <- function(ti, Pi, Pj, alpha) {
  #xi, yi = Pi
  #xj, yj = Pj
  xi <- Pi[1]
  yi <- Pi[2]
  xj <- Pj[1]
  yj <- Pj[2]
  return( ( ( (xj-xi)^2 + (yj-yi)^2 )^0.5 )^alpha + ti)
}
catmullrom <- function(data, x, y, alpha=0.5, nPoints=100) {
  x <- data[[rlang::ensym(x)]]
  y <- data[[rlang::ensym(y)]]
  if (length(x) != length(y)) {stop("x and y vectors must be the same length.")}
  # calculate start and end points
  sz <- length(x) # number of points
  dom <- max(x) - min(x)
  rng <- max(y) - min(y)
  prex <- x[1]+sign(x[1]-x[2])*dom*0.01 # 0.01 = pctdom
  prey <- (y[1]-y[2])/(x[1]-x[2])*(prex-x[1])+y[1]
  endx <- x[sz]+sign(x[sz]-x[sz-1])*dom*0.01 # 0.01 = pctdom
  endy <- (y[sz]-y[sz-1])/(x[sz]-x[sz-1])*(endx-x[sz])+y[sz]
  x <- c(prex, x, endx)
  y <- c(prey, y, endy)
  ## The curve C will contain an array of (x,y) points.
  #C = []
  #for i in range(sz-3):
  #  c = CatmullRomSpline(P[i], P[i+1], P[i+2], P[i+3],alpha)
  #C.extend(c)
  #
  #return C
  points <- matrix(c(x, y), ncol=2)
  ret <- do.call(rbind, lapply(1:(sz-1), function(i) {
    P0 <- matrix(points[i, ], nrow=1)
    P1 <- matrix(points[i+1, ], nrow=1)
    P2 <- matrix(points[i+2, ], nrow=1)
    P3 <- matrix(points[i+3, ], nrow=1)
    # Convert the points to numpy so that we can do array multiplication
    #P0, P1, P2, P3 = map(numpy.array, [P0, P1, P2, P3])
    t0 = 0
    t1 = tj(t0, P0, P1, alpha)
    t2 = tj(t1, P1, P2, alpha)
    t3 = tj(t2, P2, P3, alpha)
    # Only calculate points between P1 and P2
    #t = numpy.linspace(t1,t2,nPoints)
    t_seq <- seq(t1, t2, length.out=nPoints)
    # Reshape so that we can multiply by the points P0 to P3
    # and get a point for each value of t.
    #t = t.reshape(len(t),1)
    t_seq <- matrix(t_seq, ncol=1)
    #  Barry and Goldman's pyramid algorithm for cubic Catmull-Rom curves?
    A1 = ((t1-t_seq)/(t1-t0)) %*% P0 + ((t_seq-t0)/(t1-t0)) %*% P1
    A2 = ((t2-t_seq)/(t2-t1)) %*% P1 + ((t_seq-t1)/(t2-t1)) %*% P2
    A3 = ((t3-t_seq)/(t3-t2)) %*% P2 + ((t_seq-t2)/(t3-t2)) %*% P3
    B1 = as.numeric((t2-t_seq)/(t2-t0)) * A1 + as.numeric((t_seq-t0)/(t2-t0)) * A2
    B2 = as.numeric((t3-t_seq)/(t3-t1)) * A2 + as.numeric((t_seq-t1)/(t3-t1)) * A3
    C  = as.numeric((t2-t_seq)/(t2-t1)) * B1 + as.numeric((t_seq-t1)/(t2-t1)) * B2
    return(C)
  }))
  return(data.frame(x=ret[, 1], y=ret[, 2]))
}

stat_catmullrom <- function(mapping = NULL, data = NULL, geom = "path",
                      position = "identity", 
                      ...,
                      catmullrom_alpha = 0.5,
                      na.rm = FALSE, 
                      show.legend = NA, inherit.aes = TRUE) {
  layer(
    data = data,
    mapping = mapping,
    stat = StatCatmullRom,
    geom = geom,
    position = position,
    show.legend = show.legend,
    inherit.aes = inherit.aes,
    params = list(
      catmullrom_alpha = catmullrom_alpha,
      na.rm = na.rm,
      ...
    )
  )
}

#' @rdname ggplot2-ggproto
#' @format NULL
#' @usage NULL
#' @export
StatCatmullRom <- ggproto("StatCatmullRom", Stat,
                    required_aes = c("x", "y"),
                    
                    #default_aes = aes(size = 0.5),
                    
                    setup_params = function(data, params) {
                      params$flipped_aes <- has_flipped_aes(data, params, main_is_orthogonal = FALSE, main_is_continuous = TRUE)
                      
                      params
                    },
                    
                    extra_params = c("catmullrom_alpha", "na.rm"),
                    
                    compute_group = function(data, scales, catmullrom_alpha = 0.5, flipped_aes = FALSE) {
                      data <- flip_data(data, flipped_aes)
 
                      data_catmullrom <- catmullrom(data, x, y)
                      
                      data_catmullrom$flipped_aes <- flipped_aes
                      flip_data(data_catmullrom, flipped_aes)
                    }
)