2023-02-27

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