summaryrefslogtreecommitdiffstats
path: root/linterp.scm
diff options
context:
space:
mode:
Diffstat (limited to 'linterp.scm')
-rw-r--r--linterp.scm90
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)))))