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

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

initial import of tensor-bspline, a wrapper for the DTENSBS library

File size: 4.1 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        (compute3d evaluate3d)
28
29        (import scheme chicken)
30        (import srfi-4 foreign)
31        (require-extension lolevel)
32
33
34;; DB3INK(X,NX,Y,NY,Z,NZ,FCN,LDF1,LDF2,KX,KY,KZ,TX,TY,TZ,BCOEF,WORK,IFLAG)
35
36        (foreign-declare #<<EOF
37extern int db3ink_ (double*,int*,
38                    double*,int*,
39                    double*,int*,
40                    double*,
41                    int*, int*,
42                    int*,int*,int*,
43                    double*, double*, double*,
44                    double*,
45                    double*, int*);
46EOF
47)
48
49;; DB3VAL(XVAL,YVAL,ZVAL,IDX,IDY,IDZ,TX,TY,TZ,NX,NY,NZ,KX,KY,KZ,BCOEF,WORK)
50
51        (foreign-declare #<<EOF
52extern double db3val_ (double*,double*,double*,
53                       int*,int*,int*,
54                       double*,double*,double*,
55                       int*,int*,int*,
56                       int*,int*,int*,
57                       double*,double*);
58EOF
59)
60             
61
62
63(define f:db3ink
64  (foreign-lambda* int 
65                   ((f64vector x) (int nx) 
66                    (f64vector y) (int ny) 
67                    (f64vector z) (int nz) 
68                    (f64vector fcn) 
69                    (int ldf1) (int ldf2)
70                    (int kx) (int ky) (int kz)
71                    (f64vector tx) 
72                    (f64vector ty) 
73                    (f64vector tz) 
74                    (f64vector bcoef)
75                    (f64vector work)
76                    )
77#<<EOF
78   int iflag;
79   iflag = 0;
80   db3ink_ (x,&nx,y,&ny,z,&nz,fcn,&ldf1,&ldf2,
81            &kx,&ky,&kz,tx,ty,tz,bcoef,work,&iflag);
82   C_return(iflag);
83EOF
84))
85
86
87
88(define f:db3val
89  (foreign-lambda* double
90                   ((double x) (double y) (double z)
91                    (int idx) (int idy) (int idz)
92                    (f64vector tx) (f64vector ty) (f64vector tz)
93                    (int nx) (int ny) (int nz)
94                    (int kx) (int ky) (int kz)
95                    (f64vector bcoef) (f64vector work)
96                    )
97#<<EOF
98   double n;
99
100   n = db3val_ (&x,&y,&z,&idx,&idy,&idz,tx,ty,tz,
101                &nx,&ny,&nz,&kx,&ky,&kz, 
102                bcoef,work);
103
104   C_return (n);
105EOF
106))
107
108
109(define (compute3d kx ky kz x y z fval ) 
110
111  (let (
112        (nx (f64vector-length x))
113        (ny (f64vector-length y))
114        (nz (f64vector-length z))
115        )
116
117    (assert (and (>= nx 2) (>= ny 2) (>= nz 2)))
118    (assert (= (f64vector-length fval) (* nx ny nz)))
119
120    (assert (and (and (>= kx 2) (< kx nx))
121                 (and (>= ky 2) (< ky ny))
122                 (and (>= kz 2) (< kz nz))))
123
124    (let* (
125           (bcoef (make-f64vector (* nx ny nz)))
126           (tx    (make-f64vector (+ nx kx)))
127           (ty    (make-f64vector (+ ny ky)))
128           (tz    (make-f64vector (+ nz kz)))
129           (work  (make-f64vector (+ (* nx ny nz)
130                                     (max (* 2 kx (+ 1 nx))
131                                          (* 2 ky (+ 1 ny))
132                                          (* 2 kz (+ 1 nz))))))
133           )
134
135    (let ((status (f:db3ink x nx y ny z nz fval nx ny kx ky kz
136                            tx ty tz bcoef work)))
137
138      (values tx ty tz bcoef status))
139
140  )))
141
142
143(define (evaluate3d nx ny nz kx ky kz dx dy dz x y z tx ty tz bcoef)
144   
145    (let ((work (make-f64vector (+ (* ky kz) (+ (max kx ky) kz)))))
146
147      (f:db3val x y z dx dy dz tx ty tz nx ny nz kx ky kz bcoef work)
148   
149      ))
150
151 
152)
Note: See TracBrowser for help on using the repository browser.