2019-04-06

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