Changeset 30772 in project


Ignore:
Timestamp:
04/28/14 04:04:20 (5 years ago)
Author:
Ivan Raikov
Message:

tensor-bspline: added 2d interpolation routines

Location:
release/4/tensor-bspline/trunk
Files:
2 added
4 edited

Legend:

Unmodified
Added
Removed
  • release/4/tensor-bspline/trunk/tensor-bspline.meta

    r30747 r30772  
    2424 (author "Ivan Raikov")
    2525
    26  (synopsis "An interface to the DTENSBS package (Interpolation of two and three dimensional gridded data)."))
     26 (synopsis "An interface to the DTENSBS package (Interpolation of 2D and 3D gridded data)."))
  • release/4/tensor-bspline/trunk/tensor-bspline.scm

    r30747 r30772  
    2525(module tensor-bspline
    2626
    27         (compute3d evaluate3d)
     27        (compute2d evaluate2d compute3d evaluate3d)
    2828
    2929        (import scheme chicken)
     
    3232
    3333
     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
    3448;; DB3INK(X,NX,Y,NY,Z,NZ,FCN,LDF1,LDF2,KX,KY,KZ,TX,TY,TZ,BCOEF,WORK,IFLAG)
    3549
     
    4660EOF
    4761)
     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             
    4875
    4976;; DB3VAL(XVAL,YVAL,ZVAL,IDX,IDY,IDZ,TX,TY,TZ,NX,NY,NZ,KX,KY,KZ,BCOEF,WORK)
     
    6087             
    6188
     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))
    62110
    63111(define f:db3ink
     
    86134
    87135
     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
    88156(define f:db3val
    89157  (foreign-lambda* double
     
    105173EOF
    106174))
     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  )))
    107205
    108206
     
    141239
    142240
     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
    143250(define (evaluate3d nx ny nz kx ky kz dx dy dz x y z tx ty tz bcoef)
    144251   
  • release/4/tensor-bspline/trunk/tensor-bspline.setup

    r30747 r30772  
    77         dbsplin/dbintk.c  dbsplin/dbnslv.c  dbsplin/dbvalu.c
    88         dbsplin/dbnfac.c  dbsplin/dbspvn.c  dbsplin/dintrv.c 
    9          db3ink.c  db3val.c  dbknot.c   dbtpcf.c 
     9         db3ink.c  db3val.c  db2ink.c  db2val.c  dbknot.c   dbtpcf.c 
    1010         -j tensor-bspline -lf2c -lm)
    1111(compile -O3 -s tensor-bspline.import.scm)
  • release/4/tensor-bspline/trunk/tests/run.scm

    r30747 r30772  
    1010(print "machine-epsilon = " machine-epsilon)
    1111
    12 ;;  Input data
    1312
    1413(let* (
     
    2322       (yv  (list->f64vector ys))
    2423       (zv  (list->f64vector zs))
    25        (g   (f64vector-meshgrid xv yv zv))
    26        (xx  (car g))
    27        (yy  (cadr g))
    28        (zz  (caddr g))
    29        (fv  (make-f64vector n))
    30        (f   (lambda (n)
    31               (let ((x (f64vector-ref xx n))
    32                     (y (f64vector-ref yy n))
    33                     (z (f64vector-ref zz n)))
     24       (g2   (f64vector-meshgrid xv yv))
     25       (xx2  (car g2))
     26       (yy2  (cadr g2))
     27       (g3   (f64vector-meshgrid xv yv zv))
     28       (xx3  (car g3))
     29       (yy3  (cadr g3))
     30       (zz3  (caddr g3))
     31       (fv3 (make-f64vector n))
     32       (f3  (lambda (n)
     33              (let ((x (f64vector-ref xx3 n))
     34                    (y (f64vector-ref yy3 n))
     35                    (z (f64vector-ref zz3 n)))
    3436                (- (* x x) y z))))
     37       (fv2 (make-f64vector (* nx ny)))
     38       (f2  (lambda (n)
     39              (let ((x (f64vector-ref xx2 n))
     40                    (y (f64vector-ref yy2 n)))
     41                (- (* x x) y))))
    3542      )
    3643
     
    3845    (if (> n 0)
    3946        (begin
    40           (f64vector-set! fv (- n 1) (f (- n 1)))
     47          (f64vector-set! fv3 (- n 1) (f3 (- n 1)))
    4148          (recur (- n 1)))
    4249        ))
    4350
    44   (print "xx = " xx)
    45   (print "yy = " yy)
    46   (print "zz = " zz)
    47   (print "fv = " fv)
     51  (let recur ((n (* nx ny)))
     52    (if (> n 0)
     53        (begin
     54          (f64vector-set! fv2 (- n 1) (f2 (- n 1)))
     55          (recur (- n 1)))
     56        ))
     57
     58  (print "xx3 = " xx3)
     59  (print "yy3 = " yy3)
     60  (print "zz3 = " zz3)
     61  (print "fv3 = " fv3)
     62
     63  (print "xx2 = " xx2)
     64  (print "yy2 = " yy2)
     65  (print "fv2 = " fv2)
    4866
    4967  ;;  Set the degree and the class of continuity of the spline.
    5068  (let ((kx 2) (ky 2) (kz 2))
    5169   
    52     (let-values (((tx ty tz bcoef errc)
    53                   (compute3d kx ky kz xv yv zv fv)))
     70    (let-values (((tx3 ty3 tz3 bcoef3 errc3)
     71                  (compute3d kx ky kz xv yv zv fv3))
     72                 ((tx2 ty2 bcoef2 errc2)
     73                  (compute2d kx ky xv yv fv2)))
    5474
    55       (print "errc = " errc)
     75      (print "errc3 = " errc3)
     76      (print "errc2 = " errc2)
    5677
    57       (assert (= errc 1))
     78      (assert (= errc3 1))
     79      (assert (= errc2 1))
    5880
    59       (let* ((x -0.5)
     81      (let* (
     82             (x -0.5)
    6083             (y -1.5)
    6184             (z -2.5)
    62              (ans (evaluate3d nx ny nz kx ky kz 0 0 0 x y z tx ty tz bcoef))
     85             (ans3 (evaluate3d nx ny nz kx ky kz 0 0 0 x y z tx3 ty3 tz3 bcoef3))
     86             (ans2 (evaluate2d nx ny kx ky 0 0 x y tx2 ty2 bcoef2))
    6387             )
    6488       
    65         (print "ans = " ans))
     89        (print "ans3 = " ans3)
     90        (print "ans2 = " ans2)
    6691
    6792        ))
    6893
    69     )
     94    ))
    7095
    7196   
Note: See TracChangeset for help on using the changeset viewer.