1 files changed, 90 insertions, 0 deletions
diff --git a/linterp.scm b/linterp.scm
new file mode 100644
index 0000000..5be2b36
--- /dev/null
+++ b/linterp.scm
@@ -0,0 +1,90 @@
+;;; "linterp.scm" Interpolate array access.
+;Copyright 2005 Aubrey Jaffer
+;
+;Permission to copy this software, to modify it, to redistribute it,
+;to distribute modified versions, and to use it for any purpose is
+;granted, subject to the following restrictions and understandings.
+;
+;1.  Any copy made of this software must include this copyright notice
+;in full.
+;
+;2.  I have made no warranty or representation that the operation of
+;this software will be error-free, and I am under no obligation to
+;provide any services, by way of maintenance, update, or otherwise.
+;
+;3.  In conjunction with products arising from the use of this
+;material, there shall be no use of my name in any advertising,
+;promotional, or sales literature without prior written consent in
+;each case.
+
+;;@code{(require 'array-interpolate)}
+
+(require 'array)
+(require 'subarray)
+(require 'array-for-each)
+
+;;@args ra x1 ... xj
+;;
+;;@1 must be an array of rank j containing numbers.  @0 returns a
+;;value interpolated from the nearest j-dimensional cube of elements
+;;of @1.
+;;
+;;@example
+;;(interpolate-array-ref '#2A:fixZ32b((1 2 3) (4 5 6)) 1 0.1)
+;;                              ==> 4.1
+;;(interpolate-array-ref '#2A:fixZ32b((1 2 3) (4 5 6)) 0.5 0.25)
+;;                              ==> 2.75
+;;@end example
+(define (interpolate-array-ref ra . xs)
+  (define (mix rat x1 x2) (+ (* (- 1 rat) x1) (* rat x2)))
+  (define (iar ra xs dims)
+    (define x1 (car xs))
+    (define b1 (car dims))
+    (define idx (inexact->exact (floor (car xs))))
+    (define dim-1 (+ -1 (car dims)))
+    (set! xs (cdr xs))
+    (set! dims (cdr dims))
+    (cond ((<= x1 0) (if (null? xs)
+			 (array-ref ra 0)
+			 (iar (subarray ra idx) xs dims)))
+	  ((>= x1 dim-1) (if (null? xs)
+			     (array-ref ra dim-1)
+			     (iar (subarray ra dim-1) xs dims)))
+	  ((integer? x1) (if (null? xs)
+			     (array-ref ra idx)
+			     (iar (subarray ra idx) xs dims)))
+	  ((null? xs) (mix (- x1 idx)
+			   (array-ref ra idx)
+			   (array-ref ra (+ 1 idx))))
+	  (else (mix (- x1 idx)
+		     (iar (subarray ra idx) xs dims)
+		     (iar (subarray ra (+ 1 idx)) xs dims)))))
+  (if (null? xs)
+      (array-ref ra)
+      (iar ra xs (array-dimensions ra))))
+
+;;@args ra1 ra2
+;;
+;;@1 and @2 must be numeric arrays of equal rank.  @0 sets @1 to
+;;values interpolated from @2 such that the values of elements at the
+;;corners of @1 and @2 are equal.
+;;
+;;@example
+;;(define ra (make-array (A:fixZ32b) 2 2))
+;;(resample-array! ra '#2A:fixZ32b((1 2 3) (4 5 6)))
+;;ra              ==>  #2A:fixZ32b((1 3) (4 6))
+;;(define ra (make-array (A:floR64b) 3 2))
+;;(resample-array! ra '#2A:fixZ32b((1 2 3) (4 5 6)))
+;;ra              ==>  #2A:floR64b((1.0 3.0) (2.5 4.5) (4.0 6.0))
+;;@end example
+(define (resample-array! ra1 ra2)
+  (define scales (map (lambda (rd1 rd2)
+			(if (<= rd1 1)
+			    0
+			    (/ (+ -1 rd2) (+ -1 rd1))))
+		      (array-dimensions ra1)
+		      (array-dimensions ra2)))
+  (array-index-map! ra1
+		    (lambda idxs
+		      (apply interpolate-array-ref ra2
+			     (map * scales idxs)))))