source: project/release/4/tensor-bspline/trunk/tensor-bspline.scm @ 30772

Last change on this file since 30772 was 30772, checked in by Ivan Raikov, 6 years ago

tensor-bspline: added 2d interpolation routines

File size: 6.7 KB
Line 
1;;
2;; Chicken Scheme bindings for the DTENSBS package for interpolation
3;; of two and three dimensional gridded data using tensor products of
4;; B-spline basis functions.
5;;
6;; Chicken Scheme bindings created by Ivan Raikov.
7;;
8;; The DTENSBS package was created by Ronald Boisvert.
9;;
10;; This program is free software: you can redistribute it and/or
11;; modify it under the terms of the GNU General Public License as
12;; published by the Free Software Foundation, either version 3 of the
13;; License, or (at your option) any later version.
14;;
15;; This program is distributed in the hope that it will be useful, but
16;; WITHOUT ANY WARRANTY; without even the implied warranty of
17;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
18;; General Public License for more details.
19;;
20;; A full copy of the GPL license can be found at
21;; <http://www.gnu.org/licenses/>.
22;;
23
24
25(module tensor-bspline
26
27        (compute2d evaluate2d compute3d evaluate3d)
28
29        (import scheme chicken)
30        (import srfi-4 foreign)
31        (require-extension lolevel)
32
33
34;; DB2INK(X,NX,Y,NY,FCN,LDF,KX,KY,TX,TY,BCOEF,WORK,IFLAG)
35
36        (foreign-declare #<<EOF
37extern int db2ink_ (double*,int*,
38                    double*,int*,
39                    double*,
40                    int*, 
41                    int*, int*,
42                    double*, double*, 
43                    double*,
44                    double*, int*);
45EOF
46)
47
48;; DB3INK(X,NX,Y,NY,Z,NZ,FCN,LDF1,LDF2,KX,KY,KZ,TX,TY,TZ,BCOEF,WORK,IFLAG)
49
50        (foreign-declare #<<EOF
51extern int db3ink_ (double*,int*,
52                    double*,int*,
53                    double*,int*,
54                    double*,
55                    int*, int*,
56                    int*,int*,int*,
57                    double*, double*, double*,
58                    double*,
59                    double*, int*);
60EOF
61)
62
63;; DB2VAL(XVAL,YVAL,IDX,IDY,TX,TY,NX,NY,KX,KY,BCOEF,WORK)
64
65        (foreign-declare #<<EOF
66extern double db2val_ (double*,double*,
67                       int*,int*,
68                       double*,double*,
69                       int*,int*,
70                       int*,int*,
71                       double*,double*);
72EOF
73)
74             
75
76;; DB3VAL(XVAL,YVAL,ZVAL,IDX,IDY,IDZ,TX,TY,TZ,NX,NY,NZ,KX,KY,KZ,BCOEF,WORK)
77
78        (foreign-declare #<<EOF
79extern double db3val_ (double*,double*,double*,
80                       int*,int*,int*,
81                       double*,double*,double*,
82                       int*,int*,int*,
83                       int*,int*,int*,
84                       double*,double*);
85EOF
86)
87             
88
89
90(define f:db2ink
91  (foreign-lambda* int 
92                   ((f64vector x) (int nx) 
93                    (f64vector y) (int ny) 
94                    (f64vector fcn) 
95                    (int ldf) 
96                    (int kx) (int ky) 
97                    (f64vector tx) 
98                    (f64vector ty) 
99                    (f64vector bcoef)
100                    (f64vector work)
101                    )
102#<<EOF
103   int iflag;
104   iflag = 0;
105   db2ink_ (x,&nx,y,&ny,fcn,&ldf,
106            &kx,&ky,tx,ty,bcoef,work,&iflag);
107   C_return(iflag);
108EOF
109))
110
111(define f:db3ink
112  (foreign-lambda* int 
113                   ((f64vector x) (int nx) 
114                    (f64vector y) (int ny) 
115                    (f64vector z) (int nz) 
116                    (f64vector fcn) 
117                    (int ldf1) (int ldf2)
118                    (int kx) (int ky) (int kz)
119                    (f64vector tx) 
120                    (f64vector ty) 
121                    (f64vector tz) 
122                    (f64vector bcoef)
123                    (f64vector work)
124                    )
125#<<EOF
126   int iflag;
127   iflag = 0;
128   db3ink_ (x,&nx,y,&ny,z,&nz,fcn,&ldf1,&ldf2,
129            &kx,&ky,&kz,tx,ty,tz,bcoef,work,&iflag);
130   C_return(iflag);
131EOF
132))
133
134
135
136(define f:db2val
137  (foreign-lambda* double
138                   ((double x) (double y) 
139                    (int idx) (int idy) 
140                    (f64vector tx) (f64vector ty) 
141                    (int nx) (int ny) 
142                    (int kx) (int ky) 
143                    (f64vector bcoef) (f64vector work)
144                    )
145#<<EOF
146   double n;
147
148   n = db2val_ (&x,&y,&idx,&idy,tx,ty,
149                &nx,&ny,&kx,&ky, 
150                bcoef,work);
151
152   C_return (n);
153EOF
154))
155
156(define f:db3val
157  (foreign-lambda* double
158                   ((double x) (double y) (double z)
159                    (int idx) (int idy) (int idz)
160                    (f64vector tx) (f64vector ty) (f64vector tz)
161                    (int nx) (int ny) (int nz)
162                    (int kx) (int ky) (int kz)
163                    (f64vector bcoef) (f64vector work)
164                    )
165#<<EOF
166   double n;
167
168   n = db3val_ (&x,&y,&z,&idx,&idy,&idz,tx,ty,tz,
169                &nx,&ny,&nz,&kx,&ky,&kz, 
170                bcoef,work);
171
172   C_return (n);
173EOF
174))
175
176
177(define (compute2d kx ky x y fval ) 
178
179  (let (
180        (nx (f64vector-length x))
181        (ny (f64vector-length y))
182        )
183
184    (assert (and (>= nx 2) (>= ny 2) ))
185    (assert (= (f64vector-length fval) (* nx ny)))
186
187    (assert (and (and (>= kx 2) (< kx nx))
188                 (and (>= ky 2) (< ky ny))))
189
190    (let* (
191           (bcoef (make-f64vector (* nx ny)))
192           (tx    (make-f64vector (+ nx kx)))
193           (ty    (make-f64vector (+ ny ky)))
194           (work  (make-f64vector (+ (* nx ny)
195                                     (max (* 2 kx (+ 1 nx))
196                                          (* 2 ky (+ 1 ny))))))
197           )
198
199    (let ((status (f:db2ink x nx y ny fval nx kx ky
200                            tx ty bcoef work)))
201
202      (values tx ty bcoef status))
203
204  )))
205
206
207(define (compute3d kx ky kz x y z fval ) 
208
209  (let (
210        (nx (f64vector-length x))
211        (ny (f64vector-length y))
212        (nz (f64vector-length z))
213        )
214
215    (assert (and (>= nx 2) (>= ny 2) (>= nz 2)))
216    (assert (= (f64vector-length fval) (* nx ny nz)))
217
218    (assert (and (and (>= kx 2) (< kx nx))
219                 (and (>= ky 2) (< ky ny))
220                 (and (>= kz 2) (< kz nz))))
221
222    (let* (
223           (bcoef (make-f64vector (* nx ny nz)))
224           (tx    (make-f64vector (+ nx kx)))
225           (ty    (make-f64vector (+ ny ky)))
226           (tz    (make-f64vector (+ nz kz)))
227           (work  (make-f64vector (+ (* nx ny nz)
228                                     (max (* 2 kx (+ 1 nx))
229                                          (* 2 ky (+ 1 ny))
230                                          (* 2 kz (+ 1 nz))))))
231           )
232
233    (let ((status (f:db3ink x nx y ny z nz fval nx ny kx ky kz
234                            tx ty tz bcoef work)))
235
236      (values tx ty tz bcoef status))
237
238  )))
239
240
241(define (evaluate2d nx ny kx ky dx dy x y tx ty bcoef)
242   
243    (let ((work (make-f64vector (+ (* 3 (max kx ky)) ky))))
244
245      (f:db2val x y dx dy tx ty nx ny kx ky bcoef work)
246   
247      ))
248
249
250(define (evaluate3d nx ny nz kx ky kz dx dy dz x y z tx ty tz bcoef)
251   
252    (let ((work (make-f64vector (+ (* ky kz) (+ (max kx ky) kz)))))
253
254      (f:db3val x y z dx dy dz tx ty tz nx ny nz kx ky kz bcoef work)
255   
256      ))
257
258 
259)
Note: See TracBrowser for help on using the repository browser.