• uniform, centripedal and chordal variants
• thanks to mika rantanen on whose information the code is based on

• note that map with multiple lists as arguments works similar to vector notation. it maps all first values of each list, then all second values and so on
• p0d and similar are the individual values for each dimension of points - x, y, z, and so on. since any number of dimensions is allowed, they are not specifically named
• point vectors are given as lists for simplicity as they are so easy to handle in scheme

; helpers

(define (euclidean-distance p0 p1)
  "(number ...) (number ...) -> number
   unlimited dimensions"
  (sqrt (apply + (map (lambda (p0d p1d) (expt (- p0d p1d) 2)) p0 p1))))

(define (map-segments size f a)
  "integer procedure:{any ... -> any} list -> list
   map over each overlapping segment of length len.
   each segment is one step apart.
   example: for (1 2 3 4) size 2 maps (1 2) (2 3) (3 4)"
  (let loop ((rest a) (buffer (list)) (result (list)) (count size))
    (if (null? rest) (reverse (if (null? buffer) result (cons (apply f buffer) result)))
      (if (< count 1)
        (loop (cdr rest) (append (cdr buffer) (list (car rest))) (cons (apply f buffer) result) 0)
        (loop (cdr rest) (append buffer (list (car rest))) result (- count 1))))))

(define (map-integers count f)
  "integer procedure:{integer -> any} -> list
   map over integers from 0 to count - 1"
  (let loop ((n 0)) (if (= n count) (list) (cons (f n) (loop (+ 1 n))))))

; main algorithm

(define* (catmull-rom-interpolate-f p0 p1 p2 p3 #:optional (alpha 0.5) (tension 0))
  "(number ...) (number ...) (number ...) (number ...) [real:0..1 real:0..1] -> procedure:{real:0..1 -> (number ...)}
   return a function that gives points on a catmull-rom spline segment between p1 and p2 at fractional offsets.
   the returned function is called as (f time) where time is a real number between 0 and 1.
   to draw paths this function can be called with overlapping segments of a points series
   * alpha: 0 uniform, 0.5 centripetal, 1 chordal
   * tension: 0: smooth, 1: linear
   * no limit on the number of point dimensions"
  ; some coefficients constant for the segment are pre-calculated
  (let*
    ( (t01 (expt (euclidean-distance p0 p1) alpha))
      (t12 (expt (euclidean-distance p1 p2) alpha))
      (t23 (expt (euclidean-distance p2 p3) alpha))
      (m1
        (map
          (lambda (p0d p1d p2d)
            (*
              (- 1 tension)
              (+
                (- p2d p1d)
                (* t12
                  (-
                    (/ (- p1d p0d) t01)
                    (/ (- p2d p0d) (+ t01 t12)))))))
          p0 p1 p2))
      (m2
        (map
          (lambda (p1d p2d p3d)
            (*
              (- 1 tension)
              (+
                (- p2d p1d)
                (* t12
                  (-
                    (/ (- p3d p2d) t23)
                    (/ (- p3d p1d) (+ t12 t23)))))))
          p1 p2 p3))
      (a (map (lambda (p1d p2d m1d m2d) (+ (* 2 (- p1d p2d)) m1d m2d)) p1 p2 m1 m2))
      (b (map (lambda (p1d p2d m1d m2d) (- (* -3 (- p1d p2d)) m1d m1d m2d)) p1 p2 m1 m2))
      (c m1)
      (d p1))
    (lambda (t) (map (lambda (a b c d) (+ (* a t t t) (* b t t) (* c t) d)) a b c d))))

; example path drawing

(define (catmull-rom-path alpha tension resolution points)
  "real real integer ((number ...):point ...) -> ((number ...):point ...)
   create a smooth interpolated path from intermediate points.
   example: (catmull-rom-path 0.5 0 100 (quote ((-0.72 -0.3) (0 0) (1 0.8) (1.1 0.5) (2.7 1.2) (3.4 0.27))))"
  (let
    ( (points
        ; add one before and one after the series to interpolate between all given points,
        ; because the interpolation is only always between two of four points
        (let
          ( (first-point
              (map (lambda (p0d p1d) (- (* 2 p0d) p1d))
                (car points)
                (car (cdr points))))
            (last-point
              (map (lambda (p0d p1d) (- (* 2 p0d) p1d))
                (car (reverse points))
                (car (cdr (reverse points))))))
          (append (list first-point) points (list last-point)))))
    (apply append
      (map-segments 4
        (lambda (p0 p1 p2 p3) "map points to point lists"
          (let ((interpolate-f (catmull-rom-interpolate-f p0 p1 p2 p3 alpha tension)))
            (map-integers resolution
              (lambda (t) (interpolate-f (/ t resolution))))))
        points))))

(define (display-example-path)
  "write points to standard output in a gnuplot compatible format.
   can be plotted for example with:
   guile example.scm | gnuplot --persist -e \"set key off; plot '-' lc rgb 'blue'\""
  (let
    ( (result-points
        (catmull-rom-path 0.5 0
          100 (quote ((-0.72 -0.3) (0 0) (1 0.8) (1.1 0.5) (2.7 1.2) (3.4 0.27))))))
    (for-each
      (lambda (a)
        (display (car a))
        (display " ")
        (display (car (cdr a)))
        (display "\n"))
      result-points)))

(display-example-path)