From 237c6e380aebdcbc70bd1c9ecf7d3f6effca2752 Mon Sep 17 00:00:00 2001 From: Bryan Newbold Date: Mon, 20 Feb 2017 00:05:34 -0800 Subject: Import Upstream version 3a3 --- linterp.scm | 90 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) create mode 100644 linterp.scm (limited to 'linterp.scm') 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))))) -- cgit v1.2.3