-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathspline.rkt
197 lines (178 loc) · 8.42 KB
/
spline.rkt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
#lang typed/racket/base
;; spline.rkt -- construct spline interpolation functions from data points
;;
;; This file is part of data-frame -- https://github.com/alex-hhh/data-frame
;; Copyright (c) 2018 Alex Harsányi <[email protected]>
;;
;; This program is free software: you can redistribute it and/or modify it
;; under the terms of the GNU Lesser General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or (at your
;; option) any later version.
;;
;; This program is distributed in the hope that it will be useful, but WITHOUT
;; ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
;; FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
;; License for more details.
;;
;; You should have received a copy of the GNU Lesser General Public License
;; along with this program. If not, see <http://www.gnu.org/licenses/>.
;; Construct a spline interpolation function based on a set of data points.
;; The implementation is based on
;; https://en.wikipedia.org/wiki/Spline_interpolation
(require math/matrix
racket/match)
(define-type Data-Point (Vector Real Real))
(define-type Delta-Point (Vector Real Real))
(define-type Data-Series (U (Listof Data-Point) (Vectorof Data-Point)))
(define-type Delta-Series (Listof Delta-Point))
;; Calculate the delta series for a set of data points. The delta is the
;; difference between adjacent points in the data series. It will have
;; one-less element than data series.
(: ->delta-series (-> Data-Series Delta-Series))
(define (->delta-series data-series)
(cond ((list? data-series)
(for/list ([d0 (in-list data-series)] [d1 (in-list (cdr data-series))])
(match-define (vector x0 y0) d0)
(match-define (vector x1 y1) d1)
(vector (- x1 x0) (- y1 y0))))
((vector? data-series)
(for/list ([d0 (in-vector data-series)] [d1 (in-vector data-series 1)])
(match-define (vector x0 y0) d0)
(match-define (vector x1 y1) d1)
(vector (- x1 x0) (- y1 y0))))
(#t
(raise-argument-error "->delta-series -- unknown data-series type" data-series))))
;; Construct a B term for solving the equation system.
(: b-term (-> Delta-Point Real))
(define (b-term delta)
(match-define (vector dx dy) delta)
(* 3 (/ dy (* dx dx))))
;; Construct the B terms for the equation system. From a delta series, we add
;; together the b-term of each delta, except for the first and last element
;; which stand by themselves.
(: ->b-terms (-> Delta-Series (Listof Real)))
(define (->b-terms delta-series)
(for/list ([d0 (in-sequences (in-value #f) delta-series)]
[d1 (in-sequences delta-series (in-value #f))])
(+ (if d0 (b-term d0) 0) (if d1 (b-term d1) 0))))
;; Construct the column matrix of the B terms
(: b-matrix (-> Delta-Series (Matrix Real)))
(define (b-matrix delta-series)
(->col-matrix (->b-terms delta-series)))
;; Construct an A term for solving the equation system
(: a-term (-> Delta-Point Real))
(define (a-term delta)
(match-define (vector dx dy) delta)
(/ 1 dx))
;; Construct the A matrix
(: a-matrix (-> Delta-Series (Matrix Real)))
(define (a-matrix delta-series)
(define n (+ (length delta-series) 1))
(build-matrix
n n
(lambda ((x : Integer) (y : Integer))
(cond ((eqv? x y) ; diagonal
(cond ((eqv? x 0) ; first element
(* 2 (a-term (list-ref delta-series 0))))
((eqv? x (- n 1)) ; last element
(* 2 (a-term (list-ref delta-series (- n 2)))))
(#t ; other diagonal element
(* 2 (+ (a-term (list-ref delta-series (- x 1)))
(a-term (list-ref delta-series x)))))))
((or (eqv? x (- y 1)) (eqv? y (- x 1)))
(a-term (list-ref delta-series (min x y))))
(#t
0)))))
;; Given the K-TERMS and delta-series, compute the a_i, b_i terms for each
;; interpolation polynomial
(: ab-terms (-> (Listof Real) Delta-Series (Listof (Vector Real Real))))
(define (ab-terms k-terms delta-series)
(for/list ([k0 k-terms] [k1 (cdr k-terms)] [delta delta-series])
(match-define (vector dx dy) delta)
(vector
(- (* k0 dx) dy)
(+ (* (- k1) dx) dy))))
;; Construct the spline terms necessary for doing the interpolation. Given
;; the data-series, k-terms and delta series, we calculate the a_i, b_i terms
;; for each interpolation polynomial and return a list of vectors containing
;; the X, Y coordinates of the two points between which the polynomial will
;; interpolate; the k0 and k1 tangents at the start and end of the
;; interpolated range, followed by the a and b term for the polynomial.
(: spline-terms (-> Data-Series (Listof Real) Delta-Series
(Listof (Vector Real Real Real Real Real Real Real Real))))
(define (spline-terms data-series k-terms delta-series)
(cond ((list? data-series)
(for/list ([p0 (in-list data-series)]
[p1 (in-list (cdr data-series))]
[k0 (in-list k-terms)]
[k1 (in-list (cdr k-terms))]
[delta (in-list delta-series)])
(match-define (vector dx dy) delta)
(match-define (vector x0 y0) p0)
(match-define (vector x1 y1) p1)
(vector x0 y0 x1 y1 k0 k1 (- (* k0 dx) dy) (+ (* (- k1) dx) dy))))
((vector? data-series)
(for/list ([p0 (in-vector data-series)]
[p1 (in-vector data-series 1)]
[k0 (in-list k-terms)]
[k1 (in-list (cdr k-terms))]
[delta (in-list delta-series)])
(match-define (vector dx dy) delta)
(match-define (vector x0 y0) p0)
(match-define (vector x1 y1) p1)
(vector x0 y0 x1 y1 k0 k1 (- (* k0 dx) dy) (+ (* (- k1) dx) dy))))
(#t
(raise-argument-error "spline-terms: unknown data-series type" data-series))))
;; Compute the a, b terms for each interpolation polynomial for a list of data
;; points in DATA-SERIES
(: poly-terms (-> Data-Series (Listof (Vector Real Real Real Real Real Real Real Real))))
(define (poly-terms data-series)
(define delta-series (->delta-series data-series))
(define A (a-matrix delta-series))
(define B (b-matrix delta-series))
(define k-terms (matrix->list (matrix-solve A B)))
(spline-terms data-series k-terms delta-series))
;; Build a function that interpolates between points in the DATA-SERIES.
;; DATA-SERIES is a list of 2 or more data points, each data point is a
;; (vector/c real? real?). If only 2 data points are given, the resulting
;; spline will be a line.
;;
;; The returned function accepts one argument, X and returns the interpolated
;; value, Y. A spline function is only defined between the first and last
;; point in data-series and produces a function that follows these points
;; smoothly. Our function will extend this over the entire domain, but
;; outside the data points the spline will turn into a line with the same
;; slope as the first and last data point.
(: spline (-> Data-Series (-> Real Real)))
(define (spline data-series)
(when (or (and (vector? data-series) (< (vector-length data-series) 3))
(and (list? data-series) (< (length data-series) 3)))
(raise-argument-error
'spline "data series needs 2 or more points" data-series))
(define spline-terms (poly-terms data-series))
(lambda ((x : Real))
(let loop ([sterms spline-terms])
(match-define (vector x0 y0 x1 y1 k0 k1 a b) (car sterms))
(cond
;; At each of the endpoints, the value of the spline is y0 or y1
;; respectively (the spline goes through the specified points.
((= x x0) y0)
((= x x1) y1)
;; Before the first point, spline turns into a straight line, with k0
;; as the slope
((< x x0)
(+ (* k0 x) (- y0 (* k0 x0))))
;; Between these points, we interpolate according to the polynomial.
((<= x0 x x1)
(let* ((t (/ (- x x0) (- x1 x0)))
(^t (- 1 t)))
(+ (* t y1) (* ^t y0)
(* t ^t (+ (* a ^t) (* b t))))))
;; After last point, spline turns into a straight line, with k1 as the
;; slope.
((and (> x x1) (null? (cdr sterms)))
(+ (* k1 x) (- y1 (* k1 x1))))
(#t
(loop (cdr sterms)))))))
;;............................................................. provides ....
(provide spline)