source: project/release/4/parametric-curve/trunk/parametric-curve.scm @ 30416

Last change on this file since 30416 was 30416, checked in by Ivan Raikov, 7 years ago

parametric-curve: improved numeric stability in tests

File size: 7.2 KB
Line 
1
2;; An implementation of parametric curves.
3;;
4;; This code is inspired by the Haskell rsagl library.
5;;
6;; Copyright 2012-2014 Ivan Raikov and the Okinawa Institute of
7;; Science and Technology.
8;;
9;; This program is free software: you can redistribute it and/or
10;; modify it under the terms of the GNU General Public License as
11;; published by the Free Software Foundation, either version 3 of the
12;; License, or (at your option) any later version.
13;;
14;; This program is distributed in the hope that it will be useful, but
15;; WITHOUT ANY WARRANTY; without even the implied warranty of
16;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17;; General Public License for more details.
18;;
19;; A full copy of the GPL license can be found at
20;; <http://www.gnu.org/licenses/>.
21;;
22
23(module parametric-curve
24
25        (
26         parametric-curve?
27         simple-curve linear-curve line-segment
28         map-curve translate-curve scale-curve
29         sample-curve sample-curve*
30         iterate-curve bbox-curve
31         arc-length
32         )
33       
34        (import scheme chicken data-structures)
35 
36        (require-library srfi-1 srfi-4 extras bvsp-spline)
37
38        (import (only srfi-1 fold every list-tabulate zip concatenate)
39                (only srfi-4 f64vector make-f64vector list->f64vector f64vector->list f64vector-length f64vector-ref f64vector-set!)
40                (only extras fprintf pp)
41                (prefix bvsp-spline bvsp-spline:)
42                )
43
44;; A boundary-value preserving spline
45(define-record-type spline (make-spline n k x y d d2)
46  spline? 
47  (n spline-n)
48  (k spline-k)
49  (x spline-x)
50  (y spline-y)
51  (d spline-d)
52  (d2 spline-d2)
53  )
54
55
56;; A parametric equation with an associated interval, interpolating
57;; spline.
58(define-record-type peq (make-peq fn xmin xmax ymin ymax spl)
59  peq? 
60  (fn      peq-fn )
61  (xmin    peq-xmin)
62  (xmax    peq-xmax)
63  (ymin    peq-ymin)
64  (ymax    peq-ymax)
65  (spl     peq-spl)
66  )
67
68
69(define-record-printer (peq x out)
70  (fprintf out "#(peq fn=~A xmin=~A xmax=~A ymin=~A ymax=~A spl=~A)"
71           (peq-fn x)
72           (peq-xmin x)
73           (peq-xmax x)
74           (peq-ymin x)
75           (peq-ymax x)
76           (peq-spl x)
77           ))
78
79
80
81;; A N-dimensional parametric curve is a set of N parametric equations.
82(define (parametric-curve? x) (every peq? x))
83
84
85;; Defines a simple parametric equation.
86(define (simple-peq n k fn xmin xmax)
87  (if (> xmin xmax)
88      (simple-peq n k fn xmax xmin)
89      (let ((dx (/ (- xmax xmin) n)))
90        (let* ((x (list-tabulate (+ 1 n) (lambda (i) (+ xmin (* i dx)))))
91               (y (map fn x))
92               (ymin (fold min +inf.0 y))
93               (ymax (fold max -inf.0 y))
94               (xv (list->f64vector x))
95               (yv (list->f64vector y))
96               )
97
98          (let ((spl
99                 (let-values (((d d2 constr errc diagn)
100                               (bvsp-spline:compute n k xv yv)))
101                   (if (zero? errc)
102                       (make-spline n k xv yv d d2) 
103                       (error 'simple-peq "unable to compute interpolating spline"))
104                   ))
105                )
106
107
108
109              (make-peq fn xmin xmax ymin ymax spl)
110              ))
111          ))
112      )
113
114
115;; Samples a parametric equation at the given point of interest.
116(define (sample-peq c)
117  (let ((s (peq-spl c))
118        (min (peq-xmin c))
119        (max (peq-xmax c)))
120    (lambda (xp)
121      (and (>= xp min) (<= xp max)
122           (let ((v (bvsp-spline:evaluate (spline-n s) (spline-k s) 
123                                          (spline-x s) (spline-y s) 
124                                          (spline-d s) (spline-d2 s) 
125                                          (f64vector xp) 0)))
126             (and v (f64vector-ref v 0)))
127           ))
128    ))
129
130
131;; Samples a parametric equation at the given points of interest.
132(define (sample-peq* c)
133  (let ((s (peq-spl c)))
134    (lambda (xps)
135      (let ((xpsv (if (list? xps) (list->f64vector xps) xps)))
136        (bvsp-spline:evaluate (spline-n s) (spline-k s) 
137                              (spline-x s) (spline-y s) 
138                              (spline-d s) (spline-d2 s) 
139                              xpsv 0))
140      )))
141
142
143
144;; Samples a parametric equation at regular intervals in the range xmin..xmax inclusive.
145(define (iterate-peq c n)
146  (let* (
147         (f     (sample-peq* c))
148         (xmin  (peq-xmin c))
149         (xmax  (peq-xmax c))
150         (delta (- xmax xmin))
151         (dx    (if (zero? delta) 0
152                    (if (< n 2)
153                        (error 'iterate-peq "number of iterations must be >= 2")
154                        (/ (- xmax xmin) (- n 1)))))
155         )
156    (f (list-tabulate n (lambda (i) (+ xmin (* dx i)))))
157    ))
158
159
160;; Transforms an equation using the given function.
161(define (map-peq f x)
162  (let ((fx (sample-peq x)) (splx (peq-spl x)))
163    (let ((f1 (lambda (u) (f (fx u)))))
164      (simple-peq (spline-n splx) (spline-k splx) f1 (peq-xmin x) (peq-xmax x))
165      )))
166
167
168;; Combines two equations using the given function.
169(define (compose-peq f x y)
170  (let ((fx   (sample-peq x)) 
171        (splx (peq-spl x)) 
172        (fy   (sample-peq y)))
173    (let ((f1 (lambda (u) (f (fx u) (fy u)))))
174      (simple-peq (spline-n splx) (spline-k splx) f1 (peq-xmin x) (peq-xmax x))
175      )))
176
177
178;; Translates a parametric equation.
179(define (translate-peq x p) (map-peq (lambda (v) (+ x v)) p))
180
181
182;; Scales a parametric equation. 
183(define (scale-peq s p) (map-peq (lambda (v) (* s v)) p))
184
185
186;; Defines a simple parametric curve.
187(define (simple-curve n k fs tmin tmax)
188  (if (null? fs) '()
189      (let ((c (simple-peq n k (car fs) tmin tmax)))
190        (cons c (simple-curve n k (cdr fs) tmin tmax))
191        )))
192 
193 
194;; Samples a curve at the given point
195(define (sample-curve s)
196  (let ((scs (map sample-peq s)))
197    (lambda (t) (map (lambda (sc) (sc t)) scs))
198    ))
199
200;; Samples a curve at the given points
201(define (sample-curve* s)
202  (let ((scs (map sample-peq* s)))
203    (lambda (ts) (map (lambda (sc) (sc ts)) scs))
204    ))
205
206
207;; Linear curve of the form c1 * x + c2
208;; Argument coeffs supplies c1 and c2 for the different dimensions
209(define (linear-curve n coeffs tmin tmax )
210  (simple-curve n 1 (map (lambda (s) 
211                           (let ((c1 (car s)) (c2 (cadr s)))
212                             (cond ((and (zero? c2) (zero? c1)) (lambda (x) 0.0))
213                                   ((zero? c1) (lambda (x) c2))
214                                   (else (lambda (x) (+ (* c1 x) c2))))))
215                           coeffs)
216                tmin tmax))
217
218
219;; Line segment curve of the form (x1,xn) defined on the parameter range 0.0 .. 1.0
220(define (line-segment n coeffs )
221  (simple-curve n 1 (map (lambda (s) (lambda (x) (* x s))) coeffs) 0.0 1.0))
222
223
224;; Maps the given functions to the parametric curve.
225(define (map-curve fs c)
226  (map (lambda (f p) (map-peq f p)) fs c))
227
228
229;; Translates a parametric curve.
230(define (translate-curve xs c)
231  (map (lambda (x p) (translate-peq x p)) xs c))
232
233                         
234;; Scales a parametric curve.
235(define (scale-curve xs c)
236  (map (lambda (x p) (scale-peq x p)) xs c))
237
238                         
239;; Obtain the bounding box of a curve
240(define (bbox-curve c) (append (map peq-ymin c) (map peq-ymax c)))
241
242
243;; Samples a parametric curve at regular intervals in the range xmin..xmax inclusive.
244(define (iterate-curve c n) (map (lambda (p) (iterate-peq p n)) c))
245
246;; Computes the arc length of the parametric curve given step dx
247(define (arc-length c dx)
248  (let* ((n (inexact->exact (round (/ 1.0 dx))))
249         (v (iterate-curve c n)))
250    (let recur ((i 1) (l 0.0) 
251                (s (map (lambda (x) (f64vector-ref x 0)) v)))
252      (if (< i n)
253          (let ((s1 (map (lambda (x) (f64vector-ref x i)) v)))
254            (recur (+ 1 i)
255                   (+ l (sqrt (fold (lambda (x x1 l) (+ l (expt (- x1 x) 2))) 0.0 s s1)))
256                   s1))
257          l))
258    ))
259           
260   
261 
262
263)
Note: See TracBrowser for help on using the repository browser.