Nonempty Intervals and Generalized Arrays
Bradley J. Lucier
This SRFI is currently in withdrawn status. Here is an explanation of each status that a SRFI can hold. To provide input on this SRFI, please send email to srfi-122@nospamsrfi.schemers.org. To subscribe to the list, follow these instructions. You can access previous messages via the mailing list archive.
The author has submitted SRFI 179 and SRFI 231, which are updated versions of this SRFI. SRFI 122 has been withdrawn in favor of SRFI 231.
This SRFI specifies an array mechanism for Scheme. Arrays as defined here are quite general; at their most basic, an array is simply a mapping, or function, from multi-indices of exact integers $i_0,\ldots,i_{d-1}$ to Scheme values. The set of multi-indices $i_0,\ldots,i_{d-1}$ that are valid for a given array form the domain of the array. In this SRFI, each array's domain consists of a rectangular interval $[l_0,u_0)\times[l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$, a subset of $\mathbb Z^d$, $d$-tuples of integers. Thus, we introduce a data type called intervals, which encapsulate the cross product of nonempty intervals of exact integers. Specialized variants of arrays are specified to provide portable programs with efficient representations for common use cases.
In a 1993 post to the news group comp.lang.scheme, Alan Bawden gave a simple implementation of multi-dimensional arrays in R4RS scheme. The only constructor of new arrays required specifying an initial value, and he provided the three low-level primitives array-ref, array-set!, and array?. His arrays were defined on rectangular intervals in $\mathbb Z^d$ of the form $[0,u_0)\times\cdots\times [0,u_{d-1})$. I'll note that his function array-set! put the value to be entered into the array at the front of the variable-length list of indices that indicate where to place the new value. He offered an intriguing way to "share" arrays in the form of a routine make-shared-array that took a mapping from a new interval of indices into the domain of the array to be shared. His implementation incorporated what he called an indexer, which was a function from the interval $[0,u_0)\times\cdots\times [0,u_{d-1})$ to an interval $[0,N)$, where the body of the array consisted of a single Scheme vector of length $N$. Bawden called the mapping specified in make-shared-array linear, but I prefer the term affine, as I explain later.
Mathematically, Bawden's arrays can be described as follows. We'll use the vector notation $\vec i$ for a multi-index $i_0,\ldots,i_{d-1}$. (Multi-indices correspond to Scheme values.) Arrays will be denoted by capital letters $A,B,\ldots$, the domain of the array $A$ (in Bawden's case $[0,u_0)\times \cdots\times [0,u_{d-1})$) will be denoted by $D_A$, and the indexer of $A$, mapping $D_A$ to the interval $[0,N)$ will be denoted by $I_A$. Initially, Bawden constructs $I_A$ such that $I_A(\vec i)$ steps consecutively through the values $0,1,\ldots,N-1$ as $\vec i$ steps through the multi-indices $(0,\ldots,0,0)$, $(0,\ldots,0,1)$, $\ldots$, $(0,\ldots,1,0)$, etc., in lexicographical order, which means that if $\vec i$ and $\vec j$ are two multi-indices, then $\vec i<\vec j$ iff the first coordinate $k$ where $\vec i$ and $\vec j$ differ satisfies $\vec i_k<\vec j_k$. In fact, $I_A(\vec i)=\vec v\cdot\vec i$ for some specially-constructed vector $\vec v$ that depends only on $D_A$, the domain of $A$, where $\vec v\cdot\vec i$ is the dot product of $\vec v$ and $\vec i$.
In make-shared-array, Bawden allows you to specify a new $r$-dimensional interval $D_B$ as the domain of a new array $B$, and a mapping $T_{BA}:D_B\to D_A$ of the form $T_{BA}(\vec i)=M\vec i+\vec b$; here $M$ is a $d\times r$ matrix of integer values and $\vec b$ is a $d$-vector. So this mapping $T_{BA}$ is affine, in that $T_{BA}(\vec i)-T_{BA}(\vec j)=M(\vec i-\vec j)$ is linear (in a linear algebra sense) in $\vec i-\vec j$. The new indexer of $B$ satisfies $I_B(\vec i)=I_A(T_{BA}(\vec i))$.
A fact Bawden exploits in the code, but doesn't point out in the short post, is that $I_B$ is again an affine map, and indeed, the composition of any two affine maps is again affine.
We incorporate Bawden-style arrays into this SRFI, but extend them in two relatively minor ways that we find quite useful.
First, we allow the intervals of multi-indices that form the domains of arrays to have nonzero lower bounds as well as upper bounds, so domains are rectangular, $d$-dimensional intervals $[l_0,u_0)\times\cdots\times[l_{d-1},u_{d-1})$.
Second, we introduce the notion of a storage class, an object that contains functions that manipulate, store, check, etc., different types of values. A generic-storage-class can manipulate any Scheme value, whereas,e.g., a u1-storage-class can store only the values 0 and 1 in each element of a body.
We also require that our affine maps be one-to-one, so that if $\vec i\neq\vec j$ then $T(\vec i)\neq T(\vec j)$. Without this property, modifying the $\vec i$th component of $A$ would cause the $\vec j$th component to change.
Requiring the transformations $T_{BA}:D_B\to D_A$ to be affine may seem esoteric and restricting, but in fact many common and useful array transformations can be expressed in this way. We give several examples below:
array-extract to define this common operation; it's like looking at a rectangular sub-part of a spreadsheet. We use it to extract the common part of overlapping domains of three arrays in an image processing example below. array-translate to provide this operation.array-permute for this operation. The only nonidentity permutation of a two-dimensional spreadsheet turns rows into columns and vice versa.array-curry for this operation, which returns an array whose domain is $\text{Int}_1$ and whose elements are themselves arrays, each of which is defined on $\text{Int}_2$. Currying a two-dimensional array would be like organizing a spreadsheet into a one-dimensional array of rows of the spreadsheet.#f and $i_j\to u_j+l_j-1-i_j$ if $F_j$ is #t.In other words, we reverse the ordering of the $j$th coordinate of $\vec i$ if and only if $F_j$ is true. $T_{BA}$ is an affine mapping from $D_B\to D_A$, which defines a new array $B$, and we can provide array-reverse for this operation. Applying array-reverse to a two-dimensional spreadsheet might reverse the order of the rows or columns (or both).interval-scale and array-sample for these operations.We make several remarks. First, all these operations could have been computed by specifying the particular mapping $T_{BA}$ explicitly, so that these routines are simply "convenience" procedures. Second, because the composition of any number of affine mappings are again affine, accessing or changing the elements of a restricted, translated, curried, permuted array is no slower than accessing or changing the elements of the original array itself. Finally, we note that by combining array currying and permuting, say, one can come up with simple expressions of powerful algorithms, such as extending one-dimensional transforms to multi-dimensional separable transforms, or quickly generating two-dimensional slices of three-dimensional image data. Examples are given below.
Bawden-style arrays are clearly useful as a programming construct, but they do not fulfill all our needs in this area. An array, as commonly understood, provides a mapping from multi-indices $(i_0,\ldots,i_{d-1})$ of exact integers in a nonempty, rectangular, $d$-dimensional interval $[l_0,u_0)\times[l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$ (the domain of the array) to Scheme objects. Thus, two things are necessary to specify an array: an interval and a mapping that has that interval as its domain.
Since these two things are often sufficient for certain algorithms, we introduce in this SRFI a minimal set of interfaces for dealing with such arrays.
Specifically, an array specifies a nonempty, multi-dimensional interval, called its domain, and a mapping from this domain to Scheme objects. This mapping is called the getter of the array, accessed with the procedure array-getter; the domain of the array (more precisely, the domain of the array's getter) is accessed with the procedure array-domain.
If this mapping can be changed, the array is said to be mutable and the mutation is effected
by the array's setter, accessed by the procedure array-setter. We call an object of this type a mutable array. Note: If an array does not have a setter, then we call it immutable even though the array's getter might not be a "pure" function, i.e., the value it returns may not depend solely on the arguments passed to the getter.
In general, we leave the implementation of generalized arrays completely open. They may be defined simply by closures, or they may have hash tables or databases behind an implementation, one may read the values from a file, etc.
In this SRFI, Bawden-style arrays are called specialized. A specialized array is an example of a mutable array.
Even if an array $A$ is not a specialized array, then it could be "shared" by specifying a new interval $D_B$ as the domain of a new array $B$ and an affine map $T_{BA}:D_B\to D_A$. Each call to $B$ would then be computed as $B(\vec i)=A(T_{BA}(\vec i))$.
One could again "share" $B$, given a new interval $D_C$ as the domain of a new array $C$ and an affine transform $T_{CB}:D_C\to D_B$, and then each access $C(\vec i)=A(T_{BA}(T_{CB}(\vec i)))$. The composition $T_{BA}\circ T_{CB}:D_C\to D_A$, being itself affine, could be precomputed and stored as $T_{CA}:D_C\to D_A$, and $C(\vec i)=A(T_{CA}(\vec i))$ can be computed with the overhead of computing a single affine transformation.
So, if we wanted, we could share generalized arrays with constant overhead by adding a single layer of (multi-valued) affine transformations on top of evaluating generalized arrays. Even though this could be done transparently to the user, we do not do that here; it would be a compatible extension of this SRFI to do so. We provide only the routine specialized-array-share, not a more general array-share.
Certain ways of sharing generalized arrays, however, are relatively easy to code and not that expensive. If we denote (array-getter A) by A-getter, then if B is the result of array-extract applied to A, then (array-getter B) is simply A-getter. Similarly, if A is a two-dimensional array, and B is derived from A by applying the permutation $\pi((i,j))=(j,i)$, then (array-getter B) is (lambda (i j) (A-getter j i)). Translation and currying also lead to transformed arrays whose getters are relatively efficiently derived from A-getter, at least for arrays of small dimension.
Thus, while we do not provide for sharing of generalized arrays for general one-to-one affine maps $T$, we do allow it for the specific functions array-extract, array-translate, array-permute, array-curry, array-reverse, and array-sample, and we provide relatively efficient implementations of these functions for arrays of dimension no greater than four.
Daniel Friedman and David Wise wrote a famous paper CONS should not Evaluate its Arguments. In the spirit of that paper, our procedure array-map does not immediately produce a specialized array, but a simple immutable array, whose elements are recomputed from the arguments of array-map each time they are accessed. This immutable array can be passed on to further applications of array-map for further processing, without generating the storage bodies for intermediate arrays.
We provide the procedure array->specialized-array to transform a generalized array (like that returned by array-map) to a specialized, Bawden-style array, for which accessing each element again takes $O(1)$ operations.
specialized-array-share is, conceptually, a multi-valued array.array-curry gets its name from the
curry operator in programming---we are currying the getter of the array and keeping careful track of the domains.
interval-projections can be thought of as currying the
characteristic function of the interval, encapsulated here as interval-contains-multi-index?.(interval-cross-product interval1 interval2 ...)interval-projections),
which don't seem terribly natural for arrays.(make-array ...), array-map, and array->specialized-array to construct arrays, and while there are several other ways to construct arrays, there is no really low-level interface given for constructing specialized arrays (where one specifies a body, an indexer, etc.). It was felt that certain difficulties, some surmountable (such as checking that a given body is compatible with a given storage class) and some not (such as checking that an indexer is indeed affine), made a low-level interface less useful. At the same time, the simple (make-array ...) mechanism is so general, allowing one to specify getters and setters as general functions, as to cover nearly all needs.This document refers to translations and permutations. A translation is a vector of exact integers. A permutation of dimension $n$ is a vector whose entries are the exact integers $0,1,\ldots,n-1$, each occurring once, in any order.
Procedure: translation? object
Returns #t if object is a translation, and #f otherwise.
Procedure: permutation? object
Returns #t if object is a permutation, and #f otherwise.
An interval represents the set of all multi-indices of exact integers $i_0,\ldots,i_{d-1}$ satisfying $l_0\leq i_0<u_0,\ldots,l_{d-1}\leq i_{d-1}<u_{d-1}$, where the lower bounds $l_0,\ldots,l_{d-1}$ and the upper bounds $u_0,\ldots,u_{d-1}$ are specified multi-indices of exact integers. The positive integer $d$ is the dimension of the interval. It is required that $l_0<u_0,\ldots,l_{d-1}<u_{d-1}$.
Intervals are a data type distinct from other Scheme data types.
Procedure: make-interval lower-bounds upper-bounds
Create a new interval; lower-bounds and upper-bounds
are nonempty vectors (of the same length) of exact integers that satisfy
(< (vector-ref lower-bounds i) (vector-ref upper-bounds i)) for
$0\leq i<{}$(vector-length lower-bounds). It is an error if
lower-bounds and upper-bounds do not satisfy these conditions.
Procedure: interval? obj
Returns #t if obj is an interval, and #f otherwise.
Procedure: interval-dimension interval
If interval is an interval built with
(make-interval lower-bounds upper-bounds)then interval-dimension returns (vector-length lower-bounds). It is an error to call interval-dimension
if interval is not an interval.
Procedure: interval-lower-bound interval i
Procedure: interval-upper-bound interval i
If interval is an interval built with
(make-interval lower-bounds upper-bounds)
and i is an exact integer that satisfies
$0 \leq i<$ (vector-length lower-bounds),
then interval-lower-bound returns
(vector-ref lower-bounds i) and interval-upper-bound returns
(vector-ref upper-bounds i). It is an error to call interval-lower-bound or interval-upper-bound
if interval and i do not satisfy these conditions.
Procedure: interval-lower-bounds->list interval
Procedure: interval-upper-bounds->list interval
If interval is an interval built with
(make-interval lower-bounds upper-bounds) then interval-lower-bounds->list returns (vector->list lower-bounds) and interval-upper-bounds->list returns (vector->list upper-bounds). It is an error to call
interval-lower-bounds->list or interval-upper-bounds->list if interval does not satisfy these conditions.
Procedure: interval-lower-bounds->vector interval
Procedure: interval-upper-bounds->vector interval
If interval is an interval built with
(make-interval lower-bounds upper-bounds) then interval-lower-bounds->vector returns a copy of lower-bounds and interval-upper-bounds->vector returns a copy of upper-bounds. It is an error to call
interval-lower-bounds->vector or interval-upper-bounds->vector if interval does not satisfy these conditions.
Procedure: interval-volume interval
If interval is an interval built with
(make-interval lower-bounds upper-bounds)then, assuming the existence of vector-map, interval-volume returns
(apply * (vector->list (vector-map - upper-bounds lower-bounds)))It is an error to call interval-volume if interval does not satisfy this condition.
Procedure: interval= interval1 interval2
If interval1 and interval2 are intervals built with
(make-interval lower-bounds1 upper-bounds1)and
(make-interval lower-bounds2 upper-bounds2)respectively, then interval= returns
(and (equal? lower-bounds1 lower-bounds2) (equal? upper-bounds1 upper-bounds2))It is an error to call interval= if interval1 or interval2 do not satisfy this condition.
Procedure: interval-subset? interval1 interval2
If interval1 and interval2 are intervals of the same dimension $d$, then interval-subset? returns #t if
(<= (interval-lower-bound interval1 j) (interval-lower-bound interval2 j))and
(<= (interval-upper-bound interval1 j) (interval-upper-bound interval2 j))for all $0\leq j<d$, otherwise it returns #f. It is an error if the arguments do not satisfy these conditions.
Procedure: interval-contains-multi-index? interval index-0 index-1 ...
If interval is an interval with dimension $d$ and index-0, index-1, ..., is a multi-index of length $d$,
then interval-contains-multi-index? returns #t if
(interval-lower-bound interval j)$\leq$index-j$<$(interval-upper-bound interval j)
for $0\leq j < d$, and #f otherwise.
It is an error to call interval-contains-multi-index? if interval and index-0,..., do not satisfy this condition.
Procedure: interval-projections interval right-dimension
Conceptually, interval-projections takes a $d$-dimensional interval
$[l_0,u_0)\times [l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$
and splits it into two parts
$[l_0,u_0)\times\cdots\times[l_{d-\text{right-dimension}-1},u_{d-\text{right-dimension}-1})$
and
$[l_{d-\text{right-dimension}},u_{d-\text{right-dimension}})\times\cdots\times[l_{d-1},u_{d-1})$
This function, the inverse of Cartesian products or cross products of intervals, is used to keep track of the domains of curried arrays.
More precisely, if interval is an interval and right-dimension is an exact integer that satisfies 0 < right-dimension < d then interval-projections returns two intervals:
(values
(make-interval
(vector (interval-lower-bound interval 0)
...
(interval-lower-bound interval
(- d right-dimension 1)))
(vector (interval-upper-bound interval 0)
...
(interval-upper-bound interval
(- d right-dimension 1))))
(make-interval
(vector (interval-lower-bound interval
(- d right-dimension))
...
(interval-lower-bound interval
(- d 1)))
(vector (interval-upper-bound interval
(- d right-dimension))
...
(interval-upper-bound interval
(- d 1)))))It is an error to call interval-projections if its arguments do not satisfy these conditions.
Procedure: interval-for-each f interval
This routine assumes that interval is an interval and f is a routine whose domain includes elements of interval. It is an error to call
interval-for-each if interval and f do not satisfy these conditions.
interval-for-each calls f with each multi-index of interval as arguments, all in lexicographical order.
Procedure: interval-dilate interval lower-diffs upper-diffs
If interval is an interval with
lower bounds l0, ..., ld-1 and
upper bounds u0, ..., ud-1, and lower-diffs is a vector of exact integers L0, ..., Ld-1 and upper-diffs is a vector of exact integers U0, ..., Ud-1, then interval-dilate returns a new interval with
lower bounds l0+L0, ..., ld-1+Ld-1 and
upper bounds u0+U0, ..., ud-1+Ud-1, as long as this is a
nonempty interval. It is an error if the arguments do not satisfy these conditions.
Examples:
(interval=
(interval-dilate (make-interval '#(0 0) '#(100 100))
'#(1 1) '#(1 1))
(make-interval '#(1 1) '#(101 101))) => #t
(interval=
(interval-dilate (make-interval '#(0 0) '#(100 100))
'#(-1 -1) '#(1 1))
(make-interval '#(-1 -1) '#(101 101))) => #t
(interval=
(interval-dilate (make-interval '#(0 0) '#(100 100))
'#(0 0) '#(-50 -50))
(make-interval '#(0 0) '#(50 50))) => #t
(interval-dilate
(make-interval '#(0 0) '#(100 100))
'#(0 0) '#(-500 -50)) => error
Procedure: interval-intersect interval-1 interval-2 ...
If all the arguments are intervals of the same dimension and they have a nonempty intersection,
the interval-intersect returns that intersection; otherwise it returns #f.
It is an error if the arguments are not all intervals with the same dimension.
Procedure: interval-translate interval translation
If interval is an interval with
lower bounds l0, ..., ld-1 and
upper bounds u0, ..., ud-1, and translation is a translation with entries T0, ..., Td-1, then interval-translate returns a new interval with
lower bounds l0+T0, ..., ld-1+Td-1 and
upper bounds u0+T0, ..., ud-1+Td-1.
It is an error if the arguments do not satisfy these conditions.
One could define (interval-translate interval translation) by (interval-dilate interval translation translation).
Procedure: interval-permute interval permutation
The argument interval must be an interval, and the argument permutation must be a valid permutation with the same dimension as interval. It is an error if the arguments do not satisfy these conditions.
Heuristically, this function returns a new interval whose axes have been permuted in a way consistent with permutation.
But we have to say how the entries of permutation are associated with the new interval.
We have chosen the following convention: If the permutation is $(\pi_0,\ldots,\pi_{d-1})$, and the argument interval represents the cross product $[l_0,u_0)\times[l_1,u_1)\times\cdots\times[l_{d-1},u_{d-1})$, then the result represents the cross product $[l_{\pi_0},u_{\pi_0})\times[l_{\pi_1},u_{\pi_1})\times\cdots\times[l_{\pi_{d-1}},u_{\pi_{d-1}})$.
For example, if the argument interval represents $[0,4)\times[0,8)\times[0,21)\times [0,16)$ and the
permutation is #(3 0 1 2), then the result of (interval-permute interval permutation) will be
the representation of $[0,16)\times [0,4)\times[0,8)\times[0,21)$.
Procedure: interval-scale interval scales
If interval is a $d$-dimensional interval $[0,u_1)\times\cdots\times[0,u_{d-1})$ with all lower bounds zero, and scales is a length-$d$ vector of positive exact integers, which we'll denote by $\vec s$, then interval-scale returns the interval $[0,\operatorname{ceiling}(u_1/s_1))\times\cdots\times[0,\operatorname{ceiling}(u_{d-1}/s_{d-1})$.
It is an error if interval and scales do not satisfy this condition.
Conceptually, a storage-class is a set of functions to manage the backing store of a specialized-array. The functions allow one to make a backing store, to get values from the store and to set new values, to return the length of the store, and to specify a default value for initial elements of the backing store. Typically, a backing store is a (heterogeneous or homogeneous) vector. A storage-class has a type distinct from other Scheme types.
Procedure: make-storage-class getter setter checker maker length default
Here we assume the following relationships between the arguments of make-storage-class. Assume that the "elements" of
the backing store are of some "type", either heterogeneous (all Scheme types) or homogeneous (of some restricted type).
(maker n value) returns an object containing n elements of value value.v is an object created by (maker n value) and 0 <= i < n, then (getter v i) returns the current value of the i'th element of v, and (checker (getter v i)) => #t.v is an object created by (maker n value), 0 <= i < n, and (checker val) => #t, then (setter v i val) sets the value of the i'th element of v to val.v is an object created by (maker n value) then (length v) returns n.If the arguments do not satisfy these conditions, then it is an error to call make-storage-class.
Note that we assume that getter and setter generally take O(1) time to execute.
Procedure: storage-class? m
Returns #t if m is a storage class, and #f otherwise.
Procedure: storage-class-getter m
Procedure: storage-class-setter m
Procedure: storage-class-checker m
Procedure: storage-class-maker m
Procedure: storage-class-length m
Procedure: storage-class-default m
If m is an object created by
(make-storage-class getter setter checker maker length default)
then storage-class-getter returns getter, storage-class-setter returns setter, storage-class-checker returns checker, storage-class-maker returns maker, and storage-class-length returns length, and storage-class-default returns default. Otherwise, it is an error to call any of these routines.
Variable: generic-storage-class
Variable: s8-storage-class
Variable: s16-storage-class
Variable: s32-storage-class
Variable: s64-storage-class
Variable: u1-storage-class
Variable: u8-storage-class
Variable: u16-storage-class
Variable: u32-storage-class
Variable: u64-storage-class
Variable: f32-storage-class
Variable: f64-storage-class
Variable: c64-storage-class
Variable: c128-storage-class
generic-storage-class is defined as if by
(define generic-storage-class
(make-storage-class vector-ref
vector-set!
(lambda (arg) #t)
make-vector
vector-length
#f))sX-storage-class is defined for X=8, 16, 32, and 64 (which have default values 0 and
manipulate exact integer values between -2X-1 and
2X-1-1 inclusive),
uX-storage-class is defined for X=1, 8, 16, 32, and 64 (which have default values 0 and manipulate exact integer values between 0 and
2X-1 inclusive),
fX-storage-class is defined for X= 32 and 64 (which have default value 0.0 and manipulate 32- and 64-bit floating-point numbers), and
cX-storage-class is defined for X= 64 and 128 (which have default value 0.0+0.0i and manipulate complex numbers with, respectively, 32- and 64-bit floating-point numbers as real and imaginary parts). Each of these
could be defined simply as generic-storage-class, but it is assumed that implementations with homogeneous vectors will give definitions
that either save space, avoid boxing, etc., for the specialized arrays.
Arrays are a data type distinct from other Scheme data types.
Procedure: make-array interval getter [ setter ]
Assume first that the optional argument setter is not given.
If interval is an interval and getter is a function from
interval to Scheme objects, then make-array returns an array with domain interval
and getter getter.
It is an error to call make-array if interval and getter
do not satisfy these conditions.
If now setter is specified, assume that it is a procedure such that getter and setter satisfy: If
(i1,...,in)$\neq$(j1,...,jn)
are elements of interval and
(getter j1 ... jn) => x
then "after"
(setter v i1 ... in)
we have
(getter j1 ... jn) => x
and
(getter i1,...,in) => v
Then make-array builds a mutable array with domain interval, getter getter, and
setter setter. It is an error to call make-array if its arguments do not satisfy these conditions.
Example:
(define a (make-array (make-interval '#(1 1) '#(11 11))
(lambda (i j)
(if (= i j)
1
0))))defines an array for which (array-getter a) returns 1 when i=j and 0 otherwise.
Example:
(define sparse-array
(let ((domain
(make-interval '#(0 0)
'#(1000000 1000000)))
(sparse-rows
(make-vector 1000000 '())))
(make-array
domain
(lambda (i j)
(cond ((assv j (vector-ref sparse-rows i))
=> cdr)
(else
0.0)))
(lambda (v i j)
(cond ((assv j (vector-ref sparse-rows i))
=> (lambda (pair)
(set-cdr! pair v)))
(else
(vector-set!
sparse-rows
i
(cons (cons j v)
(vector-ref sparse-rows i)))))))))
((array-getter sparse-array) 12345 6789) => 0.
((array-getter sparse-array) 0 0) => 0.
((array-setter sparse-array) 1.0 0 0) => undefined
((array-getter sparse-array) 12345 6789) => 0.
((array-getter sparse-array) 0 0) => 1.Procedure: array? obj
Returns #t if obj is an array and #f otherwise.
Procedure: array-domain array
Procedure: array-getter array
If array is an array built by
(make-array interval getter [setter])(with or without the optional setter argument) then array-domain returns interval and array-getter returns getter.
It is an error to call array-domain or array-getter if array is not an array.
Example:
(define a (make-array (make-interval '#(1 1) '#(11 11))
(lambda (i j)
(if (= i j)
1
0))))
((array-getter a) 3 3) => 1
((array-getter a) 2 3) => 0
((array-getter a) 11 0) => is an errorProcedure: array-dimension array
Shorthand for (interval-dimension (array-domain array)). It is an error to call this function if array is not an array.
Procedure: mutable-array? obj
Returns #t if obj is a mutable array and #f otherwise.
Procedure: array-setter array
If array is an array built by
(make-array interval getter setter)then array-setter returns setter. It is an error to call array-setter
if array is not a mutable array.
Procedure: specialized-array-default-safe? [ bool ]
With no argument, returns #t if newly-constructed specialized arrays check the arguments of setters and getters by default, and #f otherwise.
If bool is #t then the next call to specialized-array-default-safe? will return #t;
if bool is #f then the next call to specialized-array-default-safe? will return #f;
otherwise it is an error.
Procedure: make-specialized-array interval [ storage-class generic-storage-class ] [ safe? (specialized-array-default-safe?) ]
Constructs a specialized-array from its arguments.
interval must be given as a nonempty interval. If given, storage-class must be a storage class; if it is not given it defaults to generic-storage-class. If given, safe? must be a boolean; if it is not given it defaults to the current value of (specialized-array-default-safe?).
The body of the result is constructed as
((storage-class-maker storage-class)
(interval-volume interval)
(storage-class-default storage-class))
The indexer of the resulting array is constructed as the lexicographical mapping of interval onto the interval [0,(interval-volume interval)).
If safe is #t, then the arguments of the getter and setter (including the value to be stored) of the resulting array are always checked for correctness.
After correctness checking (if needed), (array-getter array) is defined simply as
(lambda multi-index
((storage-class-getter storage-class)
(array-body array)
(apply (array-indexer array) multi-index)))
and (array-setter array) is defined as
(lambda (val . multi-index)
((storage-class-setter storage-class)
(array-body array)
(apply (array-indexer array) multi-index)
val))
It is an error if the arguments of make-specialized-array do not satisfy these conditions.
Examples. A simple array that can hold any type of element can be defined with (make-specialized-array (make-interval '#(0 0) '#(3 3))). If you find that you're using a lot of unsafe arrays of unsigned 16-bit integers, one could define
(define (make-u16-array interval)
(make-specialized-array interval u16-storage-class #f))
and then simply call, e.g., (make-u16-array (make-interval '#(0 0) '#(3 3))).
Procedure: specialized-array? obj
Returns #t if obj is a specialized-array, and #f otherwise. A specialized-array is an array.
Procedure: array-storage-class array
Procedure: array-indexer array
Procedure: array-body array
Procedure: array-safe? array
array-storage-class returns the storage-class of array. array-safe? is true if and only if the arguments of (array-getter array) and (array-setter array) (including the value to be stored in the array) are checked for correctness.
(array-indexer array) is assumed to be a one-to-one, but not necessarily onto, affine mapping from (array-domain array) into (array-body array).
It is an error to call any of these routines if array is not a specialized-array.
Procedure: specialized-array-share array new-domain new-domain->old-domain
Constructs a new specialized-array that shares the body of the specialized-array array.
Returns an object that is behaviorally equivalent to a specialized array with the following fields:
domain: new-domain
storage-class: (array-storage-class array)
body: (array-body array)
indexer: (lambda multi-index
(call-with-values
(lambda ()
(apply new-domain->old-domain
multi-index))
(array-indexer array)))new-domain->old-domain must be an affine one-to-one mapping from new-domain to
(array-domain array).
Note: It is assumed that affine structure of the composition of new-domain->old-domain and (array-indexer array will be used to simplify:
(lambda multi-index
(call-with-values
(lambda ()
(apply new-domain->old-domain multi-index))
(array-indexer array)))It is an error if array is not a specialized array, or if new-domain is not an interval, or if new-domain->old-domain is not a one-to-one affine mapping with the appropriate domain and range.
Example: One can apply a "shearing" operation to an array as follows:
(define a
(array->specialized-array
(make-array (make-interval '#(0 0) '#(5 10))
list)))
(define b
(specialized-array-share
a
(make-interval '#(0 0) '#(5 5))
(lambda (i j)
(values i (+ i j)))))
;; Print the "rows" of b
(array-for-each (lambda (row)
(pretty-print (array->list row)))
(array-curry b 1))
;; which prints
;; ((0 0) (0 1) (0 2) (0 3) (0 4))
;; ((1 1) (1 2) (1 3) (1 4) (1 5))
;; ((2 2) (2 3) (2 4) (2 5) (2 6))
;; ((3 3) (3 4) (3 5) (3 6) (3 7))
;; ((4 4) (4 5) (4 6) (4 7) (4 8))
This "shearing" operation cannot be achieved by combining the procedures array-extract, array-translate, array-permute, array-translate, array-curry, array-reverse, and array-sample.
Procedure: array->specialized-array array [ result-storage-class generic-storage-class ] [ safe? (specialized-array-default-safe?) ]
If array is an array whose elements can be manipulated by the storage-class
result-storage-class, then the specialized-array returned by array->specialized-array can be defined by:
(let* ((domain
(array-domain array))
(getter
(array-getter array))
(result
(make-specialized-array domain
result-storage-class
safe?))
(result-setter
(array-setter result)))
(interval-for-each (lambda multi-index
(apply result-setter
(apply getter
multi-index)
multi-index))
domain)
result)It is guaranteed that (array-getter array) is called precisely once for each multi-index in (array-domain array) in lexicographical order.
It is an error if result-storage-class does not satisfy these conditions, or if safe? is not a boolean.
Procedure: array-curry array inner-dimension
If array is an array whose domain is an interval $[l_0,u_0)\times\cdots\times[l_{d-1},u_{d-1})$, and inner-dimension is an exact integer strictly between $0$ and $d$, then array-curry returns an immutable array with domain $[l_0,u_0)\times\cdots\times[l_{d-\text{inner-dimension}-1},u_{d-\text{inner-dimension}-1})$, each of whose entries is in itself an array with domain $[l_{d-\text{inner-dimension}},u_{d-\text{inner-dimension}})\times\cdots\times[l_{d-1},u_{d-1})$.
For example, if A and B are defined by
(define interval (make-interval '#(0 0 0 0)
'#(10 10 10 10)))
(define A (make-array interval list))
(define B (array-curry A 1))
so
((array-getter A) i j k l) => (list i j k l)then B is an immutable array with domain (make-interval '#(0 0 0) '#(10 10 10)), each
of whose elements is itself an (immutable) array and
(equal?
((array-getter A) i j k l)
((array-getter ((array-getter B) i j k)) l)) => #t
for all multi-indices i j k l in interval.
The subarrays are immutable, mutable, or specialized according to whether the array argument is immutable, mutable, or specialized.
More precisely, if
0 < inner-dimension < (interval-dimension (array-domain array))then array-curry returns a result as follows.
If the input array is specialized, then array-curry returns
(call-with-values
(lambda () (interval-projections (array-domain array)
inner-dimension))
(lambda (outer-interval inner-interval)
(make-array
outer-interval
(lambda outer-multi-index
(specialized-array-share
array
inner-interval
(lambda inner-multi-index
(apply values
(append outer-multi-index
inner-multi-index))))))))Otherwise, if the input array is mutable, then array-curry returns
(call-with-values
(lambda () (interval-projections (array-domain array)
inner-dimension))
(lambda (outer-interval inner-interval)
(make-array
outer-interval
(lambda outer-multi-index
(make-array
inner-interval
(lambda inner-multi-index
(apply (array-getter array)
(append outer-multi-index
inner-multi-index)))
(lambda (v . inner-multi-index)
(apply (array-setter array)
v
(append outer-multi-index
inner-multi-index))))))))Otherwise, array-curry returns
(call-with-values
(lambda () (interval-projections (array-domain array)
inner-dimension))
(lambda (outer-interval inner-interval)
(make-array
outer-interval
(lambda outer-multi-index
(make-array
inner-interval
(lambda inner-multi-index
(apply (array-getter array)
(append outer-multi-index
inner-multi-index))))))))It is an error to call array-curry if its arguments do not satisfy these conditions.
Example:
(define a (make-array (make-interval '#(0 0)
'#(10 10))
list))
((array-getter a) 3 4) => (3 4)
(define curried-a (array-curry a 1))
((array-getter ((array-getter curried-a) 3)) 4)
=> (3 4)Procedure: array-extract array new-domain
Returns a new array with the same getter (and setter, if appropriate) of the first argument, defined on the second argument.
Assumes that array is an array and new-domain is an interval that is a sub-interval of (array-domain array). If array is a specialized array, then returns
(specialized-array-share array
new-domain
values)
Otherwise, if array is a mutable array, then array-extract returns
(make-array new-domain
(array-getter array)
(array-setter array))
Finally, if array is an immutable array, then array-extract returns
(make-array new-domain
(array-getter array))
It is an error if the arguments of array-extract do not satisfy these conditions.
Procedure: array-translate array translation
Assumes that array is a valid array, translation is a valid translation, and that the dimensions of the array and the translation are the same. The resulting array will have domain (interval-translate (array-domain array) translation).
If array is a specialized array, returns a new specialized array
(specialized-array-share
array
(interval-translate (array-domain array)
translation)
(lambda multi-index
(apply values
(map -
multi-index
(vector->list translation)))))
that shares the body of array.
If array is not a specialized array but is a mutable array, returns a new mutable array
(make-array
(interval-translate (array-domain array)
translation)
(lambda multi-index
(apply (array-getter array)
(map -
multi-index
(vector->list translation))))
(lambda (val . multi-index)
(apply (array-setter array)
val
(map -
multi-index
(vector->list translation)))))
that employs the same getter and setter as the original array argument.
If array is not a mutable array, returns a new array
(make-array
(interval-translate (array-domain array)
translation)
(lambda multi-index
(apply (array-getter array)
(map - multi-index (vector->list translation)))))
that employs the same getter as the original array.
It is an error if the arguments do not satisfy these conditions.
Procedure: array-permute array permutation
Assumes that array is a valid array, permutation is a valid permutation, and that the dimensions of the array and the permutation are the same. The resulting array will have domain (interval-permute (array-domain array) permutation).
We begin with an example. Assume that the domain of array is represented by the interval $[0,4)\times[0,8)\times[0,21)\times [0,16)$, as in the example for interval-permute, and the permutation is #(3 0 1 2). Then the domain of the new array is the interval $[0,16)\times [0,4)\times[0,8)\times[0,21)$.
So the multi-index argument of the getter of the result of array-permute must lie in the new domain of the array, the interval $[0,16)\times [0,4)\times[0,8)\times[0,21)$. So if we define old-getter as (array-getter array), the definition of the new array must be in fact
(make-array (interval-permute (array-domain array)
'#(3 0 1 2))
(lambda (l i j k)
(old-getter i j k l)))
So you see that if the first argument if the new getter is in $[0,16)$, then indeed the fourth argument of old-getter is also in $[0,16)$, as it should be. This is a subtlety that I don't see how to overcome. It is the listing of the arguments of the new getter, the lambda, that must be permuted.
Mathematically, we can define $\pi^{-1}$, the inverse of a permutation $\pi$, such that $\pi^{-1}$ composed with $\pi$ gives the identity permutation. Then the getter of the new array is, in pseudo-code, (lambda multi-index (apply old-getter ($\pi^{-1}$ multi-index))). We have assumed that $\pi^{-1}$ takes a list as an argument and returns a list as a result.
Employing this same pseudo-code, if array is a specialized-array and we denote the permutation by $\pi$, then array-permute returns the new specialized array
(specialized-array-share array
(interval-permute (array-domain array) π)
(lambda multi-index
(apply values (π-1 multi-index))))The result array shares (array-body array) with the argument.
Again employing this same pseudo-code, if array is not a specialized array, but is
a mutable-array, then array-permute returns the new mutable
(make-array (interval-permute (array-domain array) π)
(lambda multi-index
(apply (array-getter array)
(π-1 multi-index)))
(lambda (val . multi-index)
(apply (array-setter array)
val
(π-1 multi-index))))which employs the setter and the getter of the argument to array-permute.
Finally, if array is not a mutable array, then array-permute returns
(make-array (interval-permute (array-domain array) π)
(lambda multi-index
(apply (array-getter array)
(π-1 multi-index))))It is an error to call array-permute if its arguments do not satisfy these conditions.
Procedure: array-reverse array flip?
We assume that array is an array and flip? is a vector of booleans whose length is the same as the dimension of array.
array-reverse returns a new array that is specialized, mutable, or immutable according to whether array is specialized, mutable, or immutable, respectively. Informally, if (vector-ref flip? k) is true, then the ordering of multi-indices in the k'th coordinate direction is reversed, and is left undisturbed otherwise.
More formally, we introduce the function
(define flip-multi-index
(let* ((domain (array-domain array))
(lowers (interval-lower-bounds->list domain))
(uppers (interval-upper-bounds->list domain)))
(lambda (multi-index)
(map (lambda (i_k flip?_k l_k u_k)
(if flip?
(- (+ l_k u_k -1) i_k)
i_k))
multi-index
(vector->list flip?)
lowers
uppers))))Then if array is specialized, then array-reverse returns
(specialized-array-share
array
domain
(lambda multi-index
(apply values
(flip-multi-index multi-index))))Otherwise, if array is mutable, then array-reverse returns
(make-array
domain
(lambda multi-index
(apply (array-getter array)
(flip-multi-index multi-index)))
(lambda (v . multi-index)
(apply (array-setter array)
v
(flip-multi-index multi-index)))))Finally, if array is immutable, then array-reverse returns
(make-array
domain
(lambda multi-index
(apply (array-getter array)
(flip-multi-index multi-index))))) It is an error if array and flip? don't satisfy these requirements.
Procedure: array-sample array scales
We assume that array is an array all of whose lower bounds are zero, and scales is a vector of positive exact integers whose length is the same as the dimension of array.
array-sample returns a new array that is specialized, mutable, or immutable according to whether array is specialized, mutable, or immutable, respectively. Informally, if we construct a new matrix $S$ with the entries of scales on the main diagonal, then the $\vec i$th element of (array-sample array scales) is the $S\vec i$th element of array.
More formally, if array is specialized, then array-sample returns
(specialized-array-share
array
(interval-scale (array-domain array)
scales)
(lambda multi-index
(apply values
(map * multi-index (vector->list scales)))))Otherwise, if array is mutable, then array-sample returns
(make-array
(interval-scale (array-domain array)
scales)
(lambda multi-index
(apply (array-getter array)
(map * multi-index (vector->list scales))))
(lambda (v . multi-index)
(apply (array-setter array)
v
(map * multi-index (vector->list scales)))))Finally, if array is immutable, then array-sample returns
(make-array
(interval-scale (array-domain array)
scales)
(lambda multi-index
(apply (array-getter array)
(map * multi-index (vector->list scales)))))It is an error if array and scales don't satisfy these requirements.
Procedure: array-map f array . arrays
If array, (car arrays), ... all have the same domain and f is a procedure, then array-map
returns a new array with the same domain and getter
(lambda multi-index
(apply f
(map (lambda (g)
(apply g multi-index))
(map array-getter
(cons array arrays)))))It is assumed that f is appropriately defined to be evaluated in this context.
It is an error to call array-map if its arguments do not satisfy these conditions.
Procedure: array-for-each f array . arrays
If array, (car arrays), ... all have the same domain and f is an appropriate procedure, then array-for-each
calls
(interval-for-each
(lambda multi-index
(apply f
(map (lambda (g)
(apply g multi-index))
(map array-getter
(cons array
arrays)))))
(array-domain array))In particular, array-for-each always walks the indices of the arrays in lexicographical order.
It is expected that array-map and array-for-each will specialize the construction of
(lambda multi-index
(apply f
(map (lambda (g)
(apply g multi-index))
(map array-getter
(cons array
arrays)))))It is an error to call array-for-each if its arguments do not satisfy these conditions.
Procedure: array-fold kons knil array
If we use the defining relations for fold over lists from SRFI-1:
(fold kons knil lis)
= (fold kons (kons (car lis) knil) (cdr lis))
(fold kons knil '())
= knil
then (array-fold kons knil array) returns
(fold kons knil (array->list array))It is an error if array is not an array, or if kons is not a procedure.
Procedure: array-fold-right kons knil array
If we use the defining relations for fold-right over lists from SRFI-1:
(fold-right kons knil lis)
= (kons (car lis) (fold-right kons knil (cdr lis)))
(fold-right kons knil '())
= knil
then (array-fold-right kons knil array) returns
(fold-right kons knil (array->list array))It is an error if array is not an array, or if kons is not a procedure.
Procedure: array-any pred array1 array2 ...
Assumes that array1, array2, etc., are arrays, all with the same domain, which we'll call interval. Also assumes that pred is a procedure that takes as many arguments as there are arrays and returns a single value.
array-any first applies (array-getter array1), etc., to the first element of interval in lexicographical order, to which values it then applies pred.
If the result of pred is not #f, then that result is returned by array-any. If the result of pred is #f, then array-any continues with the second element of interval, etc., returning the first nonfalse value of pred.
If pred always returns #f, then array-any returns #f.
If it happens that pred is applied to the results of applying (array-getter array1), etc., to the last element of interval, then this last call to pred is in tail position.
The functions (array-getter array1), etc., are applied only to those values of interval necessary to determine the result of array-any.
It is an error if the arguments do not satisfy these assumptions.
Procedure: array-every pred array1 array2 ...
Assumes that array1, array2, etc., are arrays, all with the same domain, which we'll call interval. Also assumes that pred is a procedure that takes as many arguments as there are arrays and returns a single value.
array-every first applies (array-getter array1), etc., to the first element of interval in lexicographical order, to which values it then applies pred.
If the result of pred is #f, then that result is returned by array-every. If the result of pred is nonfalse, then array-every continues with the second element of interval, etc., returning the first value of pred that is #f.
If pred always returns a nonfalse value, then the last nonfalse value returned by pred is also returned by array-every.
If it happens that pred is applied to the results of applying (array-getter array1), etc., to the last element of interval, then this last call to pred is in tail position.
The functions (array-getter array1), etc., are applied only to those values of interval necessary to determine the result of array-every.
It is an error if the arguments do not satisfy these assumptions.
Procedure: array->list array
Stores the elements of array into a newly-allocated list in lexicographical order. It is an error if array is not an array.
Procedure: list->specialized-array l interval [ result-storage-class generic-storage-class ] [ safe? (specialized-array-default-safe?) ]
Returns a specialized-array with domain interval whose elements are the elements of the list l stored in lexicographical order. It is an error if l is not a list, if interval is not an interval, if the length of l is not the same as the volume of interval, if result-storage-class (when given) is not a storage class, if safe? (when given) is not a boolean, or if any element of l cannot be stored in the body of result-storage-class, and this last error shall be detected and raised if safe is #t.
We provide an implementation in Gambit-C; the nonstandard techniques used
in the implementation are: DSSSL-style optional and keyword arguments; a
unique object to indicate absent arguments; define-structure;
and define-macro.
Final SRFIs 25, 47, 58, and 63 deal with "Multi-dimensional Array Primitives", "Array", "Array Notation", and "Homogeneous and Heterogeneous Arrays", respectively. Each of these previous SRFIs deal with what we call in this SRFI specialized-arrays. Many of the functions in these previous SRFIs have corresponding forms in this SRFI. For example, from SRFI 63, we can translate:
(array? obj)(array? obj)(array-rank a)(array-dimension obj)(make-array prototype k1 ...)(make-specialized-array (make-interval (vector 0 ...) (vector k1 ...)) storage-class).(make-shared-array array mapper k1 ...)(specialized-array-share array (make-interval (vector 0 ...) (vector k1 ...)) mapper)(array-in-bounds? array index1 ...)(interval-contains-multi-index? (array-domain array) index1 ...)(array-ref array k1 ...)((array-getter array) k1 ...)(array-set! array obj k1 ...)((array-setter array) obj k1 ...)At the same time, this SRFI has some special features:
Image processing applications provided significant motivation for this SRFI.
Reading an image file in PGM format. On a system with eight-bit chars, one can write a function to read greyscale images in the PGM format of the netpbm package as follows. The lexicographical order in array->specialized-array guarantees the the correct order of execution of the input procedures:
(define make-pgm cons)
(define pgm-greys car)
(define pgm-pixels cdr)
(define (read-pgm file)
(define (read-pgm-object port)
(skip-white-space port)
(let ((o (read port)))
;; to skip the newline or next whitespace
(read-char port)
(if (eof-object? o)
(error "reached end of pgm file")
o)))
(define (skip-to-end-of-line port)
(let loop ((ch (read-char port)))
(if (not (eq? ch #\newline))
(loop (read-char port)))))
(define (white-space? ch)
(case ch
((#\newline #\space #\tab) #t)
(else #f)))
(define (skip-white-space port)
(let ((ch (peek-char port)))
(cond ((white-space? ch)
(read-char port)
(skip-white-space port))
((eq? ch #\#)
(skip-to-end-of-line port)
(skip-white-space port))
(else #f))))
;; The image file formats defined in netpbm
;; are problematical, because they read the data
;; in the header as variable-length ISO-8859-1 text,
;; including arbitrary whitespace and comments,
;; and then they may read the rest of the file
;; as binary data.
;; So we give here a solution of how to deal
;; with these subtleties in Gambit Scheme.
(call-with-input-file
(list path: file
char-encoding: 'ISO-8859-1
eol-encoding: 'lf)
(lambda (port)
;; We're going to read text for a while,
;; then switch to binary.
;; So we need to turn off buffering until
;; we switch to binary.
(port-settings-set! port '(buffering: #f))
(let* ((header (read-pgm-object port))
(columns (read-pgm-object port))
(rows (read-pgm-object port))
(greys (read-pgm-object port)))
;; now we switch back to buffering
;; to speed things up
(port-settings-set! port '(buffering: #t))
(make-pgm
greys
(array->specialized-array
(make-array
(make-interval '#(0 0)
(vector rows columns))
(cond ((or (eq? header 'p5)
(eq? header 'P5))
;; pgm binary
(if (< greys 256)
;; one byte/pixel
(lambda (i j)
(char->integer
(read-char port)))
;; two bytes/pixel,
;;little-endian
(lambda (i j)
(let* ((first-byte
(char->integer
(read-char port)))
(second-byte
(char->integer
(read-char port))))
(+ (* second-byte 256)
first-byte)))))
;; pgm ascii
((or (eq? header 'p2)
(eq? header 'P2))
(lambda (i j)
(read port)))
(else
(error "not a pgm file"))))
(if (< greys 256)
u8-storage-class
u16-storage-class)))))))Viewing two-dimensional slices of three-dimensional data. One example might be viewing two-dimensional slices of three-dimensional data in different ways. If one has a $1024 \times 512\times 512$ 3D image of the body stored as a variable body, then one could get 1024 axial views, each $512\times512$, of this 3D body by (array-curry body 2); or 512 median views, each $1024\times512$, by (array-curry (array-permute body '#(1 0 2)) 2); or finally 512 frontal views, each again $1024\times512$ pixels, by (array-curry (array-permute body '#(2 0 1)) 2); see Anatomical plane.
Calculating second differences of images. For another example, if a real-valued function is defined on a two-dimensional interval $I$, its second difference in the direction $d$ at the point $x$ is defined as $\Delta^2_df(x)=f(x+2d)-2f(x+d)+f(x)$, and this function is defined only for those $x$ for which $x$, $x+d$, and $x+2d$ are all in $I$. See the beginning of the section on "Moduli of smoothness" in these notes on wavelets and approximation theory for more details.
Using this definition, the following code computes all second-order forward differences of an image in the directions $d,2 d,3 d,\ldots$, defined only on the domains where this makes sense:
(define (all-second-differences image direction)
(let ((image-domain (array-domain image)))
(let loop ((i 1)
(result '()))
(let ((negative-scaled-direction
(vector-map (lambda (j)
(* -1 j i))
direction))
(twice-negative-scaled-direction
(vector-map (lambda (j)
(* -2 j i))
direction)))
(cond ((interval-intersect
image-domain
(interval-translate
image-domain
negative-scaled-direction)
(interval-translate
image-domain
twice-negative-scaled-direction))
=>
(lambda (subdomain)
(loop
(+ i 1)
(cons
(array->specialized-array
(array-map
(lambda (f_i f_i+d f_i+2d)
(+ f_i+2d
(* -2. f_i+d)
f_i))
(array-extract
image
subdomain)
(array-extract
(array-translate
image
negative-scaled-direction)
subdomain)
(array-extract
(array-translate
image
twice-negative-scaled-direction)
subdomain)))
result))))
(else
(reverse result)))))))
We can define a small synthetic image of size 8x8 pixels and compute its second differences in various directions:
(define image
(array->specialized-array
(make-array (make-interval '#(0 0) '#(8 8))
(lambda (i j)
(exact->inexact (+ (* i i) (* j j)))))))
(define (expose difference-images)
(pretty-print (map (lambda (difference-image)
(list (array-domain difference-image)
(array->list difference-image)))
difference-images)))
(begin
(display
"\nSecond-differences in the direction $k\times (1,0)$:\n")
(expose (all-second-differences image '#(1 0)))
(display
"\nSecond-differences in the direction $k\times (1,1)$:\n")
(expose (all-second-differences image '#(1 1)))
(display
"\nSecond-differences in the direction $k\times (1,-1)$:\n")
(expose (all-second-differences image '#(1 -1))))
On Gambit 4.8.5, this yields (after some hand editing):
Second-differences in the direction $k\times (1,0)$: ((#<##interval #2 lower-bounds: #(0 0) upper-bounds: #(6 8)> (2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2. 2.)) (#<##interval #3 lower-bounds: #(0 0) upper-bounds: #(4 8)> (8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8. 8.)) (#<##interval #4 lower-bounds: #(0 0) upper-bounds: #(2 8)> (18. 18. 18. 18. 18. 18. 18. 18. 18. 18. 18. 18. 18. 18. 18. 18.))) Second-differences in the direction $k\times (1,1)$: ((#<##interval #5 lower-bounds: #(0 0) upper-bounds: #(6 6)> (4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.)) (#<##interval #6 lower-bounds: #(0 0) upper-bounds: #(4 4)> (16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16.)) (#<##interval #7 lower-bounds: #(0 0) upper-bounds: #(2 2)> (36. 36. 36. 36.))) Second-differences in the direction $k\times (1,-1)$: ((#<##interval #8 lower-bounds: #(0 2) upper-bounds: #(6 8)> (4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4. 4.)) (#<##interval #9 lower-bounds: #(0 4) upper-bounds: #(4 8)> (16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16. 16.)) (#<##interval #10 lower-bounds: #(0 6) upper-bounds: #(2 8)> (36. 36. 36. 36.)))
You can see that with differences in the direction of only the first coordinate, the domains of the difference arrays get smaller in the first coordinate while staying the same in the second coordinate, and with differences in the diagonal directions, the domains of the difference arrays get smaller in both coordinates.
Separable operators. Many multi-dimensional transforms in signal processing are separable, in that that the multi-dimensional transform can be computed by applying one-dimensional transforms in each of the coordinate directions. Examples of such transforms include the Fast Fourier Transform and the Fast Wavelet Transform. Each one-dimensional subdomain of the complete domain is called a pencil, and the same one-dimensional transform is applied to all pencils in a given direction. Given the one-dimensional array transform, one can compute the multidimensional transform as follows:
(define (make-separable-transform 1D-transform)
(lambda (array)
;; Works on arrays of any dimension.
(let* ((n
(array-dimension array))
(permutation
;; we start with the identity permutation
(let ((result (make-vector n)))
(do ((i 0 (fx+ i 1)))
((fx= i n) result)
(vector-set! result i i)))))
;; We apply the one-dimensional transform
;; in each coordinate direction.
(do ((d 0 (fx+ d 1)))
((fx= d n))
;; Swap the d'th and n-1'st coordinates
(vector-set! permutation (fx- n 1) d)
(vector-set! permutation d (fx- n 1))
;; Apply the transform in the d'th coordinate
;; direction to all "pencils" in that direction.
;; array-permute re-orders the coordinates to
;; put the d'th coordinate at the end, array-curry
;; returns an $n-1$-dimensional array of
;; one-dimensional subarrays, and 1D-transform
;; is applied to each of those sub-arrays.
(array-for-each
1D-transform
(array-curry (array-permute array permutation)
1))
;; return the permutation to the identity
(vector-set! permutation d d)
(vector-set! permutation (fx- n 1) (fx- n 1))))))
We can test this by turning a one-dimensional Haar wavelet transform into a multi-dimensional Haar transform:
(define (1D-Haar-loop a)
(let ((getter (array-getter a))
(setter (array-setter a))
(n (interval-upper-bound (array-domain a) 0)))
(do ((i 0 (fx+ i 2)))
((fx= i n))
(let* ((a_i
(getter i))
(a_i+1
(getter (fx+ i 1)))
(scaled-sum
(fl/ (fl+ a_i a_i+1) (flsqrt 2.0)))
(scaled-difference
(fl/ (fl- a_i a_i+1) (flsqrt 2.0))))
(setter scaled-sum i)
(setter scaled-difference (fx+ i 1))))))
(define (1D-Haar-transform a)
;; works only on mutable arrays with domains
;; $[0, 2^k)$ for some $k$
(let ((n (interval-upper-bound (array-domain a) 0)))
(if (fx< 1 n)
(begin
;; calculate the scaled sums and differences
(1D-Haar-loop a)
;; Apply the transform to the
;; sub-array of scaled sums
(1D-Haar-transform (array-sample a '#(2)))))))
(define (1D-Haar-inverse-transform a)
;; works only on mutable arrays with domains
;; $[0, 2^k)$ for some $k$
(let* ((n (interval-upper-bound (array-domain a) 0)))
(if (fx< 1 n)
(begin
;; Apply the inverse transform to
;; get the array of scaled sums
(1D-Haar-inverse-transform
(array-sample a '#(2)))
;; reconstruct the array values from
;; the scaled sums and differences
(1D-Haar-loop a)))))
(define Haar-transform
(make-separable-transform 1D-Haar-transform))
(define Haar-inverse-transform
(make-separable-transform 1D-Haar-inverse-transform))
We then define an image that is a multiple of a single, two-dimensional Haar wavelet, compute its transform (which should be nonzero for only a single Haar coefficient), and then the inverse transform:
(let ((image
(array->specialized-array
(make-array (make-interval '#(0 0) '#(4 4))
(lambda (i j)
(if (fx< i 2) 1. -1.))))))
(display "\nInitial image: \n")
(pretty-print (list (array-domain image)
(array->list image)))
(Haar-transform image)
(display "\nArray of Haar wavelet coefficients: \n")
(pretty-print (list (array-domain image)
(array->list image)))
(Haar-inverse-transform image)
(display "\nReconstructed image: \n")
(pretty-print (list (array-domain image)
(array->list image))))
This yields:
Initial image: (#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)> (1. 1. 1. 1. 1. 1. 1. 1. -1. -1. -1. -1. -1. -1. -1. -1.)) Array of Haar wavelet coefficients: (#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)> (0. 0. 0. 0. 0. 0. 0. 0. 3.9999999999999987 0. 0. 0. 0. 0. 0. 0.)) Reconstructed image: (#<##interval #11 lower-bounds: #(0 0) upper-bounds: #(4 4)> (.9999999999999993 .9999999999999993 .9999999999999993 .9999999999999993 .9999999999999993 .9999999999999993 .9999999999999993 .9999999999999993 -.9999999999999993 -.9999999999999993 -.9999999999999993 -.9999999999999993 -.9999999999999993 -.9999999999999993 -.9999999999999993 -.9999999999999993))
In perfect arithmetic, this Haar transform is orthonormal, in that the sum of the squares of the elements of the image is the same as the sum of the squares of the Haar coefficients of the image. We can see that this is approximately true here.
The SRFI author thanks Edinah K Gnang, John Cowan, Sudarshan S Chawathe, Jamison Hope, and Per Bothner for their comments and suggestions, and Arthur A Gleckler, SRFI Editor, for his guidance and patience.
© 2016 Bradley J Lucier. All Rights Reserved.
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.