Rotation of 3D volumes by Fourier-interpolated shears

Graphical Models 68 (2006) 356–370 www.elsevier.com/locate/gmod

Rotation of 3D volumes by Fourier-interpolated shears Joel S. Welling a

a,b,*,1 ,

William F. Eddy

c,2

, Terence K. Young

d,3

Department of Statistics, Carnegie Mellon University, Pittsburgh, PA 15213-3890, USA b Pittsburgh Supercomputing Center, Pittsburgh, PA 15213, USA c Carnegie Mellon University, Pittsburgh, PA 15213-3890, USA d Department of Geophysics, Colorado School of Mines, Golden, CO 80401, USA

Received 16 February 2005; received in revised form 26 August 2005; accepted 30 November 2005

Abstract Algorithms for rotation of 2D images by multiple shearing transformations are well known; algorithms which use as few as four shears to perform an arbitrary 3D rotation on a 3D volume have also been described. By using Fourier transform methods to implement these shears, rotations can be performed completely reversibly and without loss of information. This is of great utility when the rotated data are used as input to statistical calculations, for example in human brain imaging. In general, six different patterns of four shears can be used to implement a given 3D rotation. We examine the mathematical and implementation details of these rotation algorithms. This paper provides a classification of these patterns, demonstrates that a consistent scheme must be used to select shear patterns for various rotations, and presents several such consistent schemes.  2005 Elsevier Inc. All rights reserved. Keywords: Volume rotation; Fourier; Shear; Reversibility; Image registration

1. Introduction The discipline of computer graphics provides a number of algorithms for rotating raster images. One particularly simple and efficient algorithm for *

Corresponding author. Fax: +1 412 268 5832. E-mail address: [email protected] (J.S. Welling). 1 Partially supported by NICHHD Grants P01-HD35490 and U19-HD035469-08. 2 Partially supported by NINDS Grant P01-NS35949, NIMH Grant RO1 MH67166-03, and NICHHD Grants P01-HD35490 and RO1-HD042386-03. 3 Portions of this work were performed while this author was Visiting Scholar at Carnegie Mellon University from Mobil Exploration & Producing Technical Center, Dallas, TX 752650232.

2D rotation is to factor the 2D rotation matrix into the product of three in-plane shearing matrices (by Paeth [18,19] and also by Tanaka et al. [22]). Eddy, Fitzgerald, and Noll [9] have demonstrated an implementation of this algorithm which uses Fourier transforms and phase shifts to implement the shears. This algorithm produces unbiased estimates when used for motion estimation in 2D. Eddy and Young [10] have further shown the algorithm to be information-preserving when applied to complex-valued data, in the sense that applying a rotation followed by its inverse recovers the original data exactly. This is because the three shear factors of the inverse rotation exactly cancel the original three factors.

1524-0703/$ - see front matter  2005 Elsevier Inc. All rights reserved. doi:10.1016/j.gmod.2005.11.004

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

In Eddy et al. [9] it was shown that the method generalized to 3D, factoring a 3D rotation matrix into a product of three extended 2D matrices and thus into a product of nine shearing matrices. Schroeder and Salem [21] collapsed two adjacent shears to produce an eight shear factorization. Cox and Jesmanovicz [8] showed that there is a four shear decomposition of the 3D rotation matrix. They implement these shears using Fourier transforms and phase shifts as in [9], or alternatively using polynomial interpolation along the shearing direction. Note that Wittenbrink and Somani [25] and Toffoli and Quick [23] have proposed a three shear factorization, but their definition of a single shear differs from that of [9,21] and [8] (see Section 3). In the literature of direct volume rendering, Chen and Kaufman [4,5] have examined 3D rotation using four shears of a variety of types; the results reported here are a generalization and formalization of their beam-shear method. Their implementation uses polynomial resampling for speed in the implementation of volume rotation; this paper uses Fourier methods to produce results which are statistically unbiased and information-preserving. A generalized shearing transformation called a beam-slice shear with five free parameters is introduced by [5]. This transformation allows a 3D rotation to be decomposed into only two such shears, but because of the extra free parameters the shears cannot be implemented by Fourier methods. Fourier methods are slow for large rasters, but many disciplines must rotate small rasters. For example, fMRI or PET rasters are typically of the order 64 · 64 · 32, for which Fourier methods are competitive with polynomial methods. Because the Fourier methods do not induce unnecessary correlation between adjacent voxels [16], they avoid some biases which can occur when optimizing the alignment of translated rasters. When used in conjunction with sum-of-squares error metrics, unbiased methods provide better sub-voxel estimates of registration [14]. This level of accuracy is important given the large voxels in these rasters. This paper examines 3D rotation by shearing, particularly decompositions using four shears. In Section 2, we consider the quaternion representation of 3D rotations and provide some algebraic details. In Section 3, we define shearing transformations and describe their implementation via Fourier transforms; Section 4 describes the algebraic properties of shear transformations. Section 5 considers the effect of rotation using Fourier methods on a

357

discrete image of finite extent. Section 6 studies the factorization of a rotation into a product of shears, gives a general set of solutions to the decomposition of a 3D rotation into four shears and examines the limits and singularities of those decompositions. Section 7 considers the interaction of small rotations with the boundaries of the sampled region of Fourier transform space. Section 8 deals with algorithms for choosing an appropriate decomposition for any given 3D rotation. Section 9 gives computational results showing the importance of correctly selecting shear decompositions to implement lossless volume rotation. We conclude with a brief discussion of these results. 2. Quaternion representation of rotations To analyze rotation algorithms, we must first choose an algebraic representation for the rotations themselves. The most commonly used representation is that of Euler angles, in which sequential rotation angles about the coordinate axes are given. There are a number of variants of Euler angles, but the most common in computer graphics is to give a rotation about the X axis, followed by a rotation about the (rotated) Y axis, followed by a rotation about the (twice rotated) Z axis. Euler angles are conceptually simple, but they suffer several severe disadvantages. One is the phenomenon of ‘‘gimbal lock’’, a coordinate degeneracy for values in which two of the rotation axes align. This name comes from a failure mode of a real object called a gimbal made of concentric rings. Each ring is mounted to its outer neighbor so as to rotate about one axis, and to its inner neighbor so as to rotate about a perpendicular axis. In general, this allows the innermost ring to rotate freely in any direction; gyroscopes are often mounted this way. However, if two axes of rotation come to line up, it becomes impossible to rotate the innermost ring about one axis- hence ‘‘gimbal lock’’. More importantly for our purposes, the different coordinate axes are not treated equivalently in Euler angle representations. Thus, the physical symmetry between different axes is lost when algebra is done in terms of Euler angles. For this reason, we favor a quaternion representation of rotations. This representation is widely used in computer graphics; see for example [17]. Conceptually, quaternions are much like the angleaxis representations of rotations. They contain four real values called X, Y, Z, and W, such that their

358

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

squares sum to unity. If the X, Y, and Z values are taken to be the components of a 3-vector, that vector lies along the rotation axis and the ratio of its magnitude to that of the W value is the tangent of half the total angle of rotation. Any given rotation can be represented by two distinct quaternions, representing rotations of opposite senses about antiparallel axes. Formally, a quaternion q is defined as

Thus, R must represent a rotation about the (X, Y, Z)T direction. For a rotation purely about the X axis (Y = Z = 0), 0 1 0 1 1 0 0 1 0 0 B C B C Rx ¼ @ 0 1  2X 2 2XW A  @ 0 cos h  sin h A:

q  W þ X i þ Y j þ Zk;

Using the half-angle relation tanðh=2Þ ¼ ð1  cos hÞ= sin h, it is easy to show that X ¼ sinðh=2Þ and W ¼ cosðh=2Þ. Extending this to the general case of rotation about any axis, one finds that a rotation by angle h about the unit vector (X, Y, Z)T can be represented by the quaternion

where W, X, Y, and Z are real and i, j, and k do not commute and obey i2 ¼ j2 ¼ k2 ¼ 1 and ij ¼ k ¼ ji;

jk ¼ i ¼ kj;

ki ¼ j ¼ ik:

Two types of quaternions of particular interest are the unit quaternions and pure quaternions. Unit quaternions are those for which X2 + Y2 + Z2 + W2 = 1, while for pure quaternions W = 0. Thus, one can write a pure quaternion v as follows: v  ai þ bj þ ck: To demonstrate the relationship between quaternions and rotations, we will consider the expression v 0 = q v q1 where v and v 0 are unnormalized pure quaternions, q and q1 are unit quaterions, and q

1

¼ W  X i  Y j  Zk:

We will demonstrate that this can be rewritten as v 0 = Rv where v 0 and v are geometrical 3-vectors with the same coefficients as those of the corresponding pure quaternions and R has the properties of a rotation matrix. The i, j, and k components of v 0 = qvq1 must be separately equated. This provides three simultaneous equations, which can be put into matrix form as follows: 1 0 a0 1  2ðY 2 þ Z 2 Þ 2ðXY  ZW Þ B 0C B @ b A ¼ @ 2ðXY þ ZW Þ 1  2ðX 2 þ Z 2 Þ c0 2ðXZ  YW Þ 2ðYZ þ XW Þ 0

2

2

10 1 a 2ðXZ þ YW Þ CB C 2ðYZ  XW Þ A@ b A: c 1  2ðX 2 þ Y 2 Þ 2

2

Remembering that X + Y + Z + W = 1, it is easy to show that this matrix (call it R) is orthogonal and has determinant 1. Thus, it is a pure rotation matrix. R leaves the vector (X, Y, Z)T invariant: 0 1 0 1 X X B C B C R@ Y A ¼ @ Y A : Z Z

0

2XW

1  2X 2

0 sin h

cos h

h h q ¼ cos þ ðX i þ Y j þ ZkÞ sin : 2 2 Note that a rotation by h about an axis v is equivalent to a rotation by h about the axis v, so there exist two quaternions for every rotation. Quaternions are free of gimbal lock and other coordinate singularities, and have a number of other pleasant algebraic properties. Most importantly, the X, Y, and Z axis directions are treated equivalently, so permutations of axes are easy to implement. Henceforth, we write the quaternion q = q(X, Y, Z, W) in the following matrix representation: 0

1  2ðY 2 þ Z 2 Þ

B 2ðXY þ ZW Þ B RðqÞ ¼ B @ 2ðXZ  YW Þ 0

2ðXY  ZW Þ 1  2ðZ 2 þ X 2 Þ 2ðYZ þ XW Þ 0

2ðXZ þ YW Þ

0

1

0C C C 1  2ðX þ Y Þ 0 A 0 1 2ðYZ  XW Þ 2

2

ð1Þ

This and later matrices are written in homogeneous coordinates [13] for convenience when we later deal with translations. For small rotations X, Y, and Z become small while W approaches 1.0. Specifically, for a rotation of , X ¼ rx ; Y ¼ ry ; Z ¼ rz ;

ð2Þ

W ¼ 1  Oð2 Þ where r = (rx, ry, rz)T is the unit vector in the direction of the rotation axis. For a 1 rotation,  is approximately 0.0013, so a small- approximation is useful for motions of the scale found in, for example, functional magnetic resonance imaging (fMRI) (see, for example [9]).

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

3. Shearing transformations A shearing transformation is one in which displacement in one coordinate direction scales linearly with location in another direction. A shearing transformation performed on a deck of cards, parallel to the card surfaces, skews it in such a way that the cards remain flat and do not rotate with respect to each other, but that each card is displaced a horizontal distance proportional to its vertical distance from the center card of the deck producing a parallelogram in side view. In 3D, there are three distinct shears, given by 0

1 a

1

b

Dx

1

0

0 0 0

1 0 0

0 C C C; 0 A

1 0

a 1

0 0 1 0

0 0

B0 1 B S z ða; b; Dz Þ ¼B @a b

0 1

1 1 0 0C C C: Dz A

0 0

0

1

B0 B S x ða; b; Dx Þ ¼B @0 0 0 1 Bb B S y ða; b; Dy Þ ¼B @0 0

1 1 0 Dy C C C; 0 A

ð3Þ

For brevity we will refer to Sx (a, b, 0) as Sx (a, b) etc. Note that this convention differs from that of [8] in that our Sy (a, b) corresponds to their S2 (b, a). We use this form so that Sx, Sy, and Sz are related by cyclic permutation of the X, Y, and Z axes. The D terms represent translations in the given coordinate directions. Where D terms are not given they are assumed to be 0. The important feature of shear transformations is that they can be implemented using simple linear shifts [24]. In particular, since linear shifts can be implemented using phase shifts in Fourier transform space, these transformations can be efficiently and accurately implemented using the Fast Fourier Transform. Specifically, to apply Sx (a, b,Dx) to data sampled on an Nx by Ny by Nz grid using Fourier transforms, one would: 1. Perform an FFT in the Y and Z directions. aL n z nz 2. Apply a phase shift by 2pðDx þ LxyN yy þ bL Þnx . Lx N z 3. Perform an inverse FFT in the Y and Z directions.

359

Here, Lx is the overall physical length of the grid in the X direction, and similarly for Ly and Lz. The a and b are given in terms of fractional shears, and have values between 1 and 1 inclusive. Dx is given as a fraction of Lx and thus also has a value between 1 and 1. We take nx to be the grid index in the X direction, centered so as to range from N/2 to (N/2)  1, and similarly for ny and nz. This corresponds to band limited sinc filtering as described in [24]. Note that we have imposed no requirement that the grid have matching dimensions or voxel sizes in the X, Y, and Z directions. 3.1. Fourier shearing of real data There is an important complication to the application of the phase shift described here when the input data are purely real. The N grid locations span a range from bN/2c to d(N/2)e  1. If N is even, all grid points except the first will have matching points equidistant on the other side of the origin. The first point, at N/2, contains data sampled at the Nyquist frequency limit [3]. For purely real imagespace data, the properties of the discrete Fourier transform guarantee that the corresponding k-space data will have Hermitian symmetry about the origin. For even N the first grid location has no symmetrical partner, however, and so it is not surprising that the value in k-space must be purely real. This can be seen directly by examination of the image-space terms which sum to give this value. When a shear operation is performed on purely real data, one expects that the resulting sheared data will also be purely real. For this to be true, the kspace sample at the Nyquist frequency must be real both before and after all phase shifts have been applied. If the ‘‘energy’’ associated with this highest frequency sample is to be preserved, the magnitude of the sample must also be unchanged. Thus, it is clear that no phase shift is applied to the N/2 sample under these circumstances. The physical interpretation of this effect is that there is not enough information in a discrete real signal to determine the phase of the Nyquist frequency component of its Fourier transform. Only a minimum energy (associated with a phase of zero) can be calculated. Assuming that one wants the data before and after shearing to have the same minimum energy, the Nyquist frequency component must not be phase shifted if the input data are real. Failure to follow this rule produces clear information loss, as demonstrated in Section 9.

360

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

4. Algebra of shears We give here some algebraic rules for shear matrices as defined by Eq. (3). These shear matrices are related by cyclic permutation of the coordinate axes, so their algebraic rules can also be generated by cyclic permutation. Rules for one shearing direction are given. Shear transformations along one axis can be composed with shears along the same axis and, consequently, they commute. Specifically, S x ða; b; D1 ÞS x ðc; d; D2 Þ ¼ S x ða þ c; b þ d; D1 þ D2 Þ ¼ S x ðc; b; D1 ÞS x ða; d; D2 Þ ¼ S x ða; d; D1 ÞS x ðc; b; D2 Þ ¼ S x ðc; d; D1 ÞS x ða; b; D2 Þ ¼ S x ða; b; D2 ÞS x ðc; d; D1 Þ ¼ S x ðc; b; D2 ÞS x ða; d; D1 Þ ¼ S x ða; d; D2 ÞS x ðc; b; D1 Þ ¼ S x ðc; d; D2 ÞS x ða; b; D1 Þ

ð4Þ

and similarly for Sy and Sz. An obvious corollary is S 1 x ða; b; DÞ ¼ S x ða; b; DÞ and similarly for Sy and Sz. The general commutation relation for two shears along different axes is 0

ad B0 B S x ða; bÞS y ðc; dÞ  S y ðc; dÞS x ða; bÞ ¼ B @0 0

0 ad

ac bd

0 0

0 0

1 0 0C C C: 0A 0

Two useful special cases are S x ða; bÞS y ðc; 0Þ  S y ðc; 0ÞS x ða; bÞ ¼ S x ð0; acÞ  I; S x ð0; bÞS y ðc; dÞ  S y ðc; dÞS x ð0; bÞ ¼ S y ðbd; 0Þ  I; where I is the identity matrix. 5. Rotation and the k-space unit cell If discrete image space (i-space) data of finite extent are translated to k-space using a discrete Fourier transform, the discrete k-space data will occupy a finite region whose boundaries are defined by the maximum Nyquist sampling frequency acquired. If the image space data consist of N samples separated by the distance l, the kspace data will consist of N samples separated by steps of 1/l in spatial frequency. We refer to this N/l region in k-space as the unit cell. The remainder of k-space consists of periodic repeti-

tions of this unit cell. This becomes clear if one uses phase shifts to shear a region of k-space as described in the previous section; data shifted outside the unit cell through one edge reenter the cell through the opposite edge. Thus, when a rotation is to be performed, it is important to consider the interaction of the rotation with the boundaries of the unit cell. A true global rotation would rotate all replications of the cell in k-space about their single common center, but our goal is to rotate the unit cell at the origin about the center of its own cell. Note that all correct implementations of rotation must carry the same volume through the cell boundaries, as that volume is a function only of the overall rotation of the cell. One approach is to require that a rotation and its inverse be implemented in such a way that the inverse repairs the effect on the unit cell caused by the original rotation. Implementations of rotation with this property are information-conserving in the sense of [10]. We will refer to this property as reversibility. An implementation of rotation can be irreversible either because the fundamental operators involved destroy information or because the sequence of operations involved in a rotation is not exactly undone by the sequence used to implement the inverse rotation. In general, neither rotations nor shears commute, so a requirement of reversibility places strong constraints on the order of operators used to implement a rotation and its inverse. If the data are zero near the boundaries of the image, well-chosen implementations of small rotations are not likely to carry non-zero data through the boundaries. In fMRI data acquired as axial slices it is typically the case that the data are compact in the slice plane. Care must be taken to insure that the data also have compact support in the throughplane direction, typically by padding with zeros. If the data to be rotated are in fact k-space data (so that the data are actually transformed to i-space in the process of the rotation) the situation is more difficult. In this case, the data are not compact, and in fact the highest-frequency components of the data lie on the boundaries. Translation operators do not in general commute with rotations or shears either, so an algorithm for constructing operators implementing combined translations and rotations must be very carefully chosen if the operators it produces are to be reversible. We do not implement these combined opera-

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

tors reversibly, but rather assume that translation always follows rotation as described above. 6. Rotations via shearing We wish to represent full 3D rotation matrices as products of some number of shears, extending the method of [18,19] to 3D. A (non-homogeneous) 3D rotation has nine matrix elements, but is orthonormal. This provides six constraints, leaving three free parameters. Intuitively these correspond to the three Euler angles necessary to specify a 3D rotation. A 3D rotation cannot be represented as a product of fewer than four shears. Shear matrices have determinant 1, so a product of shears will also have determinant 1. The product of four shears will have eight free parameters and a unit determinant. This provides sufficient freedom to set the matrix elements of the product equal to the matrix elements of the desired rotation matrix, so four shears are sufficient to represent a 3D rotation. Three shears are not sufficient, as the six degrees of freedom are over-constrained. Decompositions using more shears are also possible, of course. For example, one can represent a 3D rotation as a product of three Euler angle rotations, and each of these can be represented with three shears after the manner of [18]. This leads to a nine-shear representation as described in [9]. In the following subsections we develop a classification scheme for the possible four-shear rotations and examine their properties. Section 6.1 enumerates the decompositions while Section 6.2 examines their symmetries and singularities. Section 6.3 examines small-angle expansions. Finally, Section 6.4 describes the effects of incorporating a translation into a set of four shears. 6.1. Classification of four-shear rotations Successive shears in a decomposition should be along different axes because of the composition rule in Eq. (4). Thus, there are 3 · 2 · 1 · 1 = 6 possible four-shear decompositions of a 3D rotation. Three of these decompositions can be reached from the decomposition SySzSxSy by cyclic permutation of the coordinate axes; the other three can be reached by cyclic permutations starting from SySxSzSy. We will arbitrarily call the first set forward decompositions, and the second set backward decompositions. Each element of the first

361

set will be shown to have a corresponding complement in the second set. We will take one possible decomposition as a prototype. Given the expression RðqðX ; Y ; Z; W ÞÞ ¼ S y ða; bÞS z ðc; dÞS x ðe; f ÞS y ðg; hÞ ð5aÞ (which expresses the rotation parameterized by a quaternion as the product of four shears) and Eqs. (1) and (3) one can solve for the real shear coefficients a, . . . , h. The result is X2 þ Y2 ; YZ  XW Y 2  Z2 bðX ; Y ; Z; W Þ ¼ XY  ZW YZðX 2 þ Y 2 Þ ; 2 ðXY  ZW ÞðYZ  XW Þ YZ ; cðX ; Y ; Z; W Þ ¼ 2 XY  ZW dðX ; Y ; Z; W Þ ¼ 2ðXW  YZÞ; eðX ; Y ; Z; W Þ ¼ 2ðXY  ZW Þ; XY f ðX ; Y ; Z; W Þ ¼  2 ; YZ  XW X2  Y2 gðX ; Y ; Z; W Þ ¼ YZ  XW XY ðY 2 þ Z 2 Þ ; þ2 ðXY  ZW ÞðYZ  XW Þ aðX ; Y ; Z; W Þ ¼

hðX ; Y ; Z; W Þ ¼

ð5bÞ

1 þ X 2 þ W 2 XY  ZW

Note that if X, Y, Z, and W are all negated in the above solution, the values of a, . . ., h are unchanged. This is expected, because both quaternions represent the same rotation. Given this decomposition, two other decompositions can be generated by cyclic permutation of X, Y, and Z. For example, the equation RðqðX ; Y ; Z; W ÞÞ ¼ S z ða0 ; b0 ÞS x ðc0 ; d 0 ÞS y ðe0 ; f 0 ÞS z ðg0 ; h0 Þ has the solution p0 ðX ; Y ; Z; W Þ ¼ pðY ; Z; X ; W Þ where p is any of a, . . ., h in Eq. (5b). Consider ~a; . . . ; ~h such that the rotation R can be decomposed as follows: ~ y ð~a; ~bÞ: RðqðX ; Y ; Z; W ÞÞ ¼ S y ð~ g; ~hÞS x ð~e; f~ ÞS z ð~c; dÞS A rotation by an angle h about a given axis is the inverse of the rotation by h about the same axis,

362

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

so R (q(X, Y, Z, W)) = R1 (q(X, Y, Z, W)). Writing the right-hand side in terms of Eq. (5a) and expanding the left-hand side as SySxSzSy as above, one immediately finds ~pðX ; Y ; Z; W Þ ¼ pðX ; Y ; Z; W Þ

ð6Þ

in where again p is any of a, . . ., h in Eq. (5b). We refer to this solution for the SySxSzSy decomposition as the complement of the original SySzSxSy decomposition. Cyclic permutation of X, Y, and Z produce two more solutions related to SySxSzSy; these decompositions are the complements of those produced by cyclic permutations of SySzSxSy.

cannot be decomposed in this way. This is no surprise, as the 2D shear expansion for a rotation about Y (as in [18]) cannot be put in the form of Eq. (5a). In this manner each decomposition and its complement will be unusable for rotations about one coordinate axis. Substituting Y = 0 in Eq. (5b) produces the following simplified form, valid for Eq. (5a) when the axis of rotation lies in the XZ plane. aðX ; 0; Z; W Þ ¼  X =W ; bðX ; 0; Z; W Þ ¼Z=W ; cðX ; 0; Z; W Þ ¼0; dðX ; 0; Z; W Þ ¼2XW ;

6.2. Symmetry and singularity of the shear decompositions

eðX ; 0; Z; W Þ ¼  2ZW ; f ðX ; 0; Z; W Þ ¼0;

Consider the coefficients for the SySzSxSy decomposition given in Eq. (5b). This decomposition cannot be used for rotations which would make any coefficient undefined. Reversing the sign of any pair of X, Y, Z, and W in these coefficients leaves the magnitude of a, . . ., h unchanged (although some coefficients change sign). Thus, if Eq. (5b) gives a valid decomposition for any quaternion q, reflections of that quaternion through any coordinate axis will also have valid decompositions. By inspection of Eq. (5b) where W 5 0 and either X 5 0 or Z 5 0, if the decomposition of some quaternion q(X, Y, Z, W) is singular, the decomposition of q(X, Y, Z, W) will be nonsingular. However, by Eq. (6), these coefficients are the negatives of the coefficients of the complement decomposition of the original quaternion. Thus, at any point where a given decomposition is singular, its complement decomposition will be nonsingular. One of the two quaternions representing any given rotation will have W P 0. Consider Eq. (5b) subject to this choice of W. (No solutions will be lost, as we have already noted that the shear coefficients are unchanged if X, Y, Z, and W all change sign). By inspection, no forward decomposition has any singularity within the octant where X < 0, Y > 0, Z > 0. By Eq. (6), in the opposite octant no backward decomposition will have any singularity. We have already noted that reflection through any coordinate axis leaves singularities invariant, so either the forward or backward decompositions will be non-singular in any given octant. Setting X = Z = 0 in Eq. (5b) produces singular coefficients. Thus, a rotation about the Y axis

gðX ; 0; Z; W Þ ¼  X =W ; hðX ; 0; Z; W Þ ¼Z=W : This solution specializes easily to the cases of rotations about the X or Z axes. Substituting W = 0 and additionally setting any of X, Y or Z to zero in Eq. (5b) produces singular coefficients, and no symmetry or cyclic permutation of coordinates can remove the singularities. Thus, no four-shear decomposition exists for any 180 rotation about an axis in the XY, YZ or ZX planes. With the exception of these 180 rotations, any rotation can thus be analytically represented as an appropriate product of four shears. Numerical implementation requires some additional care, as the coefficients given in Eq. (5b) are subject to large numerical error in some regions. For example, the identity rotation (X = Y = Z = 0, W = 1) has ´ Hospital’s rule. An implementor a, . . . , h = 0 by L must take care to avoid numerical errors in this region. Some of Eqs. (5b) will be singular if XY  ZW = 0 or if YZ  XW = 0. Considering only the three forward decompositions, there is only one solution to these expressions for which all cyclic permutations will also be zero. This occurs if X = Y = Z = W. Interestingly, this quaternion corresponds to a 120 rotation about the (1, 1, 1) axis, which is the rotation which carries the x coordinate axis into the y coordinate direction, y into z, and z into x. That is, there is no forward decomposition for the rotation which cyclically permutes the coordinate axes. (This quaternion has a backward decomposition, however.)

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

6.3. Small angle expansions of the decompositions Poor numerical accuracy can result when evaluating the expressions of Eq. (5b) for rotation angles near zero. This difficulty can be avoided by expanding these results for small rotations as described in Eq. (2). The shear coefficients of Eq. (5a) become 2 aðr; Þ ¼ r1 x ð1  rz Þ;

bðr; Þ ¼

2 r1 z ðrz



D0x ¼Dx ; D0y ¼Dy  ðbDx þ aDz Þ;

ð8Þ

D0z ¼Dz  cDx : Results for other decompositions of R follow similarly. The results in Eq. (8) are identical to those of [8]. 7. Four-shear rotations and the k-space unit cell

r2y Þ;

cðr; Þ ¼ 2ry ; dðr; Þ ¼ 2rx ;

ð7Þ

eðr; Þ ¼ 2rz ; f ðr; Þ ¼ 2ry ; 2 2 gðr; Þ ¼ r1 x ðry  rx Þ; 2 hðr; Þ ¼ r1 z ð1  rx Þ

neglecting lower order terms (of order o()). Note that this expansion must be modified if any of the components of r are small. 6.4. Adding a translation Each of the shear transformations in Eq. (3) can include a translation in the shearing direction. This allows an opportunity to carry out an arbitrary translation in the course of doing a series of shears by Fourier transform without requiring any additional Fourier transforms. Specifically, consider a translation T(Dx, Dy, Dz) given by the following homogeneous matrix: 0

363

1

1

0

0

Dx

B0 B T ðDx ; Dy ; Dz Þ ¼ B @0

1 0

0 1

Dy C C C: Dz A

0

0

0

1

Assuming that the translation is to be applied after the rotation and that the rotation is to be decomposed as a series of shears R(q(X, Y, Z, W)) = SySzSxSy, one can write T ðDx ; Dy ; Dz ÞRðqðX ; Y ; Z; W ÞÞ ¼ S y ða; b; D0y ÞS z ðc; d; D0z ÞS x ðe; f ; D0x ÞS y ðg; hÞ: Solving the matrix equation, one finds that the values of a, . . . , h are unchanged from those given by Eq. (5b) and that the appropriate D’s are

We now examine in detail the interaction of fourshear rotations with the k-space unit cell. The origin corresponds to the zero frequency point. We adopt the convention that the cell occupies k-space coordinates from (1/2, 1/2] in each direction, so that the boundaries of the cell are open in the negative-facing directions. The boundary conditions of the cell are periodic, so that any data which are shifted through the face at coordinate 1/2 re-emerge from the 1/2 face and vice versa. Assuming that all six decompositions of some given rotation are non-singular, they must all carry the same volume through the cell boundaries, since that volume is a function only of the overall rotation. It is easy to compute this volume in the small angle limit of Eq. (7), as the volumes contributed by shears in different directions become effectively independent. A volume displaced out of a unit cell by the initial Sy(g, h) of Eq. (5a) can be returned to the cell by the final Sy(a, b), so the pairs of coefficients in those shears must be treated together. Thus, it is easy to see that the displaced volume Vd scales as follows: V d / ja þ gj þ jb þ hj þ jcj þ jdj þ jej þ jf j

ð9Þ

for small rotations . This result can be seen to extend to the other shear decompositions as well by considering the symmetries by which they were generated from the solution to Eq. (5a). It is interesting to examine the excursion of a point originally corresponding to a corner of the unit cell as the series of shears proceeds. Successive shears carry this point through cell boundaries and back again in a complicated pattern which depends on the initial corner, the sequence of shears chosen and the relative weights of the shear coefficients (which in turn depend on the rotation and/or translation being implemented). A small volume of data originally just within the corner are transported with this point. This subvolume is very small, scaling as 3 for small angles, so it has little numerical

364

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

impact. The geometry of the excursion is informative, however. Consider the effect of rotating the point at the (1/ 2, 1/2, 1/2)T corner of the cell through a small angle, using our canonical forward decomposition SySzSxSy (see Fig. 1). We will choose coordinate axes such that all components of the rotation axis r are positive, and the ry component is largest. Assuming that all of the components of r are large relative to the angle scale factor , the shear coefficients are well-behaved and have the values given by Eq. (7). The first shear applied is Sy(g, h); for this rotation axis g and h are both positive. This shear will carry the point of interest (point 1 of Fig. 1) through the plane given by the equation y = 1/2, specifically to the point (1/2, 1/2 + (g + h)/2, 1/2)T. This point is equivalent to the point (1/2, 1/2 + (g + h)/2, 1/2)T; the point has reentered the unit cell through the opposite face (point 2). If one carries out the remaining shears in the sequence one finds that for this choice of rotation axis Sx(e, f) carries the point through the x = 1/2 plane to point 3 of Fig. 1, and that Sz(c, d) carries the point through the z = 1/2 plane to point 4. Because of the choice of signs and relative magnitudes for the components of r, Sy(a, b) does not carry the point through a boundary plane (point 5). To leading order the final location of the point is 011 2 B C S y ða; bÞS z ðc; dÞS x ðe; f ÞS y ðg; hÞ@ 12 A 1 2

0

1  12 þ ðry þ rz Þ B 1 r2y r2y  C C ¼B @  2 þ rx þ rz  A :  12 þ ðry  rx Þ ð10Þ

The form of the Y component is different from that of the other two components because the decomposition used contains two shears along Y. By the same method, it is possible to construct a rotation as a product of shears so that small neighborhoods of all corners return to their initial positions. The somewhat surprising fact that the implementation of a single rotation can have this property is a consequence of the degeneracy of k-space in the neighborhood of the cell corners. One implementation with this property can be produced by composing the rotation from three axis-aligned rotations in the manner of Euler angles. If each of these three rotations is implemented as two half-rotations, each of which is implemented as per [19], 18 shears result. After combining adjacent shears which lie along the same axis, the rotation can be implemented using 13 shears. This is a very high computational cost given the small volume of the corner regions, but it is the simplest known method which returns all points in the volume to their original positions. 8. Selecting a four-shear decomposition Theoretically, a rotation R can be implemented using any shear decomposition with non-singular coefficients a, . . . , h. However, decompositions with small coefficients are to be preferred, as any loss or error in the implementation of shearing will grow as the shear coefficients grow. Thus, we define a figure of merit fm(a, . . . , h) P 0 and select the decomposition for which fm is smallest. Since presumably any error or loss in the implementation of a shear will be symmetrical in the shear coefficients, fm should be symmetrical in a, . . . , h and in the signs of a, . . . , h as well. Some obvious candidates are fm ða; . . . ; hÞ ¼ maxðjaj; jbj; . . . ; jhjÞ; fm ða; . . . ; hÞ ¼jaj þ jbj þ . . . þ jhj; 2

2

ð11Þ

2

fm ða; . . . ; hÞ ¼a þ b þ . . . þ h ; where the first expression given is used by [8]. Another possibility is fm ða; . . . ; hÞ ¼ ja þ gj þ jb þ hj þ jcj þ jdj þ jej þ jf j

Fig. 1. Excursion of a point near the boundary of the k-space unit cell under the action of a four-shear rotation. See the text for further details.

after Eq. (9). If fm is chosen to be symmetrical in the shear coefficients, a given rotation and its inverse will produce decompositions which are complements.

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

For example, if some quaternion q = (X, Y, Z, W) has as its best decomposition S y S z S x S 0y for some particular a, . . . , h,  q ¼ ðX ; Y ; Z; W Þ will have 1 1 1 the best decomposition S 10 by Eq. (6). y Sx Sz Sy Thus, a rotation and its inverse will exactly cancel, since the individual shear terms pair up and cancel. Assuming that the individual shear operations conserve information, pairs of rotations and their inverses will be reversible in the sense of Section 5. Because of the symmetries of the various shear decompositions, basing fm on Eq. (9) shares this feature. Figs. 2–7 show the shear decomposition which produces the smallest fm for each of the formulae above. In these figures the possible values of a quaternion q = (X, Y, Z, W) are shown as a sphere of unit radius. Each point within the sphere has X, Y, and Z values corresponding to those quaternion

Fig. 2. Optimal forward four-shear decompositions using fm(a, . . . , h) = max(jaj, . . . , jhj). See the text for further details.

Fig. 3. A pair of complementary four-shear decompositions using fm (a, . . . , h) = max(jaj, . . . , jhj). See the text for further details.

365

components; the corresponding positive W is chosen. Since every possible rotation R(q(X, Y, Z, W)) can also be represented as R(q(X, Y, Z, W)), all possible rotations can be represented as points within this sphere. Figs. 2 and 3 show three-dimensional views of this spherical volume. The center of the sphere is located at the origin and the bounding box spans the values between 1 and 1 along each axis. The X, Y, and Z coordinate axes are shown in red, green, and blue respectively; the crossing point of the axes is (1, 1, 1). Fig. 2 shows the regions where fm is minimal for each of the forward shear decompositions S y S z S x S 0y , S z S x S y S 0z and S x S y S z S 0x in red, green, and blue, respectively. Animated views of this and the following figure are available at http://www.stat.cmu.edu/ ~welling/rotations_by_shearing. Regions where a backward decomposition minimizes fm are transparent and appear as holes in the unit sphere. Note that the shape is symmetrical under reflection through any coordinate axis. Roughness in the outer surface of the sphere is a rendering artifact, as is the slight rounding at points of contact between different decompositions. Fig. 3 shows the regions where fm is minimal for the forward decomposition S y S z S x S 0y (in red) or for its complement S 0y S x S z S y (in cyan). All other decompositions are transparent. Note that if the forward decomposition is optimal for some point, the backward decomposition will be optimal for the reflection of that point through any coordinate plane. This implies that points reflected through the origin have the same property, so that any rotation and its inverse will be optimally represented by a complementary pair of decompositions. Figs. 4–7 show a series of constant-Z slices through this spherical volume for the various fm above. In each figure the Z location of the slice is most negative at the top left and increases left to right and top to bottom. Slices are chosen in equal steps and symmetrically to include pairs of slices at positive and negative Z. Positions where the three forward decompositions are optimal are shown in red, green, and blue respectively, while their complementary decompositions are shown in cyan, magenta, and yellow. (Slight asymmetries between slices at positive and negative Z are due to slight asymmetries in the numerical sampling process which produced the images.)

366

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

Fig. 4. Slices through the volume of optimal four-shear decompositions using fm(a, . . . , h) = max(jaj, . . . , jhj). See the text for further details.

Fig. 5. Slices through the volume of optimal four-shear decompositions using fm (a, . . ., h) = jaj + jbj +. . .+ jhj. See the text for further details.

Fig. 6. Slices through the volume of optimal four-shear decompositions using fm (a, . . ., h) = a2 + b2 + . . . + h2. See the text for further details.

Fig. 7. Slices through the volume of optimal four-shear decompositions using fm (a, . . . , h) = ja + gj + jb + hj + jcj + jdj + jej + jfj. See the text for further details.

9. Results A series of numerical experiments were performed comparing different methods of implementing rotation using shears, and comparing different methods of implementing the shears themselves. Translations are well understood and hence were

not tested. The software package FIASCO [15,12] was used to implement these methods. Each 3D rotation tested was decomposed into shears in several ways. The four shear decomposition used was chosen to minimize the sum-of-squares. A 7-shear decomposition was formed by pairing

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

a forward four-shear decomposition and its complement, each rotating through half the total angle. For small angles this sequence of shears is left-right symmetrical. A longer decomposition was formed by factoring the total rotation into a series of rotations about the Z, Y, X, Y, and Z axes in that order, such that pairs of rotations about a given axis were identical. (This decomposition can be found by a simple iterative algorithm). Appropriately chosing the set of three shears for each of these five rotations leads to a 13-shear decomposition which is completely leftright symmetrical even for large overall rotation angles. Timings for these three decompositions for several matrix sizes are given in Table 1. These timings are made using code based on FFTW [11], an efficient, portable Fast Fourier Transform package. Note that this implementation of the code is not optimal, in that it does not make good use of the special properties of the Fourier transform of real data. The expected scaling rules (proportionality to number of shears and to N log(N) in the matrix size) are roughly followed, but some differences occur, presumably due to details of the Fourier Transform implementation. These rotation algorithms were applied to two datasets. The first was constructed from a set of brain images acquired by EPI fMRI at 64 · 64 (3.125 mm2 in-slice voxel size) by 7-slice resolution (7 mm slices with 6 mm slice separation), averaged over an experiment consisting of roughly 300 acquisitions. The averaged volume was padded with nine slices consisting of zero values and linearly interpolated (‘‘resliced’’) between slices to produce a 64 · 64 by 64 volume consisting of cubic 3.125 mm voxels. The rotation algorithm does not require cubic voxels, since the procedure described in Section 3 does not require them. The choice of cubic voxels here is made purely for ease in comparing the performance of the algorithms with this dataset and the second (described below). In image space the brain appears at the Table 1 Times on a 2.2 GHz Pentium IV single processor workstation with a 512 KB L2 cache to perform a single rotation of the given type. Times are in seconds Matrix Size

32 · 32 · 32

64 · 64 · 64

128 · 128 · 128

4 shear 7 shear 13 shear

0.055 0.098 0.19

0.49 0.88 1.7

4.7 8.5 17.0

367

Fig. 8. Central 25 slices of brain dataset for 3D rotation tests. Padding of the dataset gives the impression of missing slices.

center of this volume with a wide empty field on all sides. Fig. 8 shows the central 25 slices of this dataset.4 The second dataset was also a 64 · 64 · 64 cube of voxels. In this case, the dataset consists of a dot at the center point surrounded by a series of concentric cubical 3D shells. The center point here is the point to which a 3D Fourier transform maps a constant input, the zero-frequency point, rather than the geometrical center of the grid (which is offset a half voxel from this point). This dataset is extremely high contrast and thus very demanding to rotate accurately. Note that some other software packages (for example, AIR [26,2] and AFNI [7,1]) rotate about the geometrical center and hence would be at a disadvantage in rotating this particular dataset. Fig. 9 shows the central 25 slices of this dataset. A modified version of this dataset was constructed by zeroing out the Nyquist frequency component of the image in k-space in all three directions. This ‘‘filtered shells’’ dataset was used in some numerical experiments which follow to demonstrate the dependency of the results on the Nyquist frequency component. For each dataset, the following numerical experiments were performed. Datasets were rotat4

This data are not distorted; this is the shape of the subject’s brain anatomy.

368

J.S. Welling et al. / Graphical Models 68 (2006) 356–370 Table 2 Summed squared errors for one 1 forward rotation about (1, 1, 1) followed by its inverse

Fig. 9. Central 25 slices of shells dataset for 3D rotation tests. Padding of the dataset gives the impression of missing slices.

ed in 1 increments about the (1, 1, 1) axis. In the first set of experiments the datasets were rotated one step and then rotated back to their original position (a forward rotation and its backward inverse). In the second set the datasets underwent 360 consecutive 1 rotations. In the final set the datasets underwent 360 positive 1 rotations followed by 360 negative 1 rotations. In all three cases the final position should be identical to the initial position, so direct measurement of the effect of the rotations should be possible. Note that this may not be exactly true for the 360-steps-forward experiment, because even small numerical errors in setting the step size will result in a slight angular displacement between the initial and final positions. In the other two experiments we presume that these errors cancel; the backward rotation step is assumed to be of exactly the same magnitude as the forward step even if that magnitude is not precisely 1. Each of these numerical experiments was performed twice, once with the data in image space and once in k-space. In the former case the (multiple) rotation operations were simply applied to the original image-space data, producing image space output. In the latter case, the image-space input data were first transformed to k-space. All rotations were then sequentially applied to this k-space data, treating it as fully complex rather than as the Fourier transform of a real image. The final result was transformed to image space, taking the modulus in the process. The image space method does not rotate the Nyquist (highest) spatial frequency component, as discussed in Section 3.1.

Method

Brain

Shells

Image

k-space

Image

k-space

4 shear 7 shear 13 shear

8.15e  9 1.48e  8 4.11e  8

3.55e  9 6.14e  9 2.67e  8

2.32e  10 4.37e  10 1.09e  9

9.28e  11 1.68e  10 5.70e  10

4 shear fixed cubic quintic heptic

0.317 161.5 122.6 106.7

0.296 — — —

0.099 — — —

0.111 — — —

Table 3 Summed squared errors for 360 1 forward rotations about (1, 1, 1) Method

4 shear 7 shear 13 shear

Brain

Shells

Filtered shells

Image

k-space

Image

k-space

Image

k-space

320.06 320.17 321.03

240.85 239.33 238.91

135.96 134.28 134.68

223.88 223.82 223.73

14.108 14.295 14.281

15.107 15.141 15.003

For each experiment, the voxelwise mean squared difference between the initial and final datasets was calculated over the central 32 · 32 · 32 cube of these 64 · 64 · 64 datasets. This difference is calculated in image space, whether the data were rotated in image space or k-space; clipping to the central region of the image avoids direct contamination by edge effects. Tables 2–4 give the results of these experiments. Tables 2–4 show summed squared errors for these test datasets after various patterns of repetitive rotations. In Table 4 and in all rows of Table 2 except the ‘‘4 shear fixed’’ case, each forward rotation is paired with an inverse rotation, so the SSE arises purely from numerical round-off error. The k-space rotation method effectively assigns more bits for the phase information, reducing this source of error, so for these tests k-space methods produce better results. For Table 2 equivalent values have been computed using various orders of polynomial interpolation for comparison5. Because all finite impulse response filters have a limited passband, one expects that higher order interpolating filters or non-polynomial interpolating filters such as those derived using the Parks–McClellan algorithm [20,6] would produce results intermediate between those

5

There values were computed with AFNI’s 3drotate utility.

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

369

Table 4 Summed squared errors for 360 1 forward rotations about (1, 1, 1) followed by 360 inverse rotations Method

4 shear 7 shear 13 shear

Brain

Shells

Filtered shells

Image

k-space

Image

k-space

Image

k-space

1.308e  5 4.080e  5 2.211e  4

1.251e  5 7.188e  5 1.756e  4

4.776e  7 8.857e  7 4.663e  6

2.749e  7 4.587e  7 4.227e  6

1.696e  9 3.145e  9 3.173e  8

1.575e  9 3.793e  9 2.766e  8

of the given polynomial interpolators and the Fourier methods. In the case of Table 3 all rotations are forward and should sum to precisely 360. Rounding error in the angle of the individual rotations can cause a slight mis-alignment between the rotated and unrotated datasets, introducing another contribution to the SSE. In the case of the brain dataset numerical rounding error dominates as it does in the Table 2, and again the k-space method produces slightly better results than the image space method. The shells dataset has a very strong component at the Nyquist spatial frequency, however. As discussed in Section 3.1, that component is not transported under the image-space rotation method. Thus, that component of the output signal cannot become mis-aligned with the matching component of the input signal, reducing the SSE. This causes the ispace rotation method to produce better results for the shell dataset under this rotation pattern. This can be seen by comparison with the results produced by performing these rotations on the ‘‘filtered shells’’ version of this dataset, which differs from the ‘‘shells’’ dataset in that the Nyquest spatial frequency components have been set to zero. In this case, the image space and k-space versions of the rotation produce very similar SSEs. In these cases it is desirable to perform the rotations in k-space, because of the additional phase information preserved by this method. In the case of 360 successive forward rotations, that additional information results in reduced overall error. In the cases where rotations are followed by their inverse rotations, the accuracy of the result is ultimately limited to round-off error. Table 2 includes entries for a rotation method in which a single shear decomposition is always used, identified as ‘‘4 shear fixed,’’ specifically the SySzSxSy decomposition. In this case, the inverse rotation does not undo the changes introduced by the forward rotation. This is a demonstration that, for four-shear rotations, an appropriate method of picking decompositions is necessary to produce reversibility.

10. Discussion Decomposition of 3D rotations into sequences of shears provides an efficient and easy-to-implement algorithm for carrying out 3D rotations. The polynomial interpolation methods of [24] can be used to implement these shears, but for problems requiring preservation of information an FFT-based implementation is preferable. These methods can be quite efficient and extraordinarily accurate if care is taken to respect subtle details. Appropriate choice of a sequence of shears is one such detail. When the input data are purely real the phase of the Nyquist frequency component is unknown; failure to respect this fact leads to loss of energy in that frequency band. Any image rotation carries some fraction of the input image outside of the output region, and fills some part of the output image with information from outside the input. In traditional image-space rotations the information carried into the image is assumed to be some arbitrary ‘‘background color’’, a rather unsatisfactory solution. Carrying out the rotation by k-space shears has the feature that this information comes from the adjacent k-space unit cells, which contain replicas of the input image. This provides a natural and sensible implementation of periodic boundary conditions. References [1] AFNI (Analysis of Functional Neuroimages): This software package is available from . [2] AIR (Automatic Image Registration): This software package is available from . [3] R.N. Bracewell, The Fourier Transform and its Applications, McGraw-Hill, New York, 1978. [4] Baoquan Chen, Arie Kaufman, 3D Volume Rotation Using Shear Transformations, Graphical Models (formerly Graphical Models and Image Processing, GMIP) 62 (2000) 308–322. [5] Baoquan Chen, Arie Kaufman, Two-pass image and volume rotation, IEEE Workshop on Volume Graphics 2001 (2001). [6] Ingrid Carlbom, Optimal filter design for volume reconstruction and visualization, IEEE Visualization, 1993 (1993) 54–61.

370

J.S. Welling et al. / Graphical Models 68 (2006) 356–370

[7] R.W. Cox, AFNI: Software for analysis and visualization of functional magnetic resonance neuroimages, Computers and Biomedical Research 29 (1996) 162–173. [8] Robert W. Cox, Andrzej Jesmanowicz, Real-time 3D image registration for functional MRI, Magnetic Resonance in Medicine 42 (1999) 1014–1018. [9] William F. Eddy, Mark Fitzgerald, Douglas C. Noll, Improved image registration by using Fourier interpolation, Magnetic Resonance in Medicine 36 (6) (1996) 923–931. [10] William F. Eddy, Terence K. Young, Optimizing the resampling of registered images, in: I.N. Bankman (Ed.), Handbook of Medical Image Processing, Processing and Analysis, Academic Press, London, 2000, pp. 603–612. [11] Fastest Fourier Transform in the West, Available from: . [12] Fiasco Functional Image Analysis Software - Computational Olio, Available from: . [13] Foley, van Dam, Feiner, Hughes, Computer Graphics Principles and Practice, second ed., Addison-Wesley, Uppersadle River, NJ, 1997. [14] J.V. Hajnal, N. Saeed, E.J. Soar, A. Oatridge, I.R. Young, G.M. Byddr, A registration and interpolation procedure for subvoxel matching of serially acquired MR images, Computer Assisted Tomography 19(2) (1995) 289–296. [15] N.A. Lazar, W.F. Eddy, C.R. Genovese, J.S. Welling, Statistical Issues in fMRI for Brain Imaging, International Statistical Review 69 (1) (2001) 105–127. [16] Thomas M. Lehmann, Claudia Go¨nner, Klaus Spitzer, Survey: interpolation methods in medical image processing, in: IEEE Transactions on Medical Imaging 18(11) (1999) 1049–1075.

[17] Patrick-Gilles Maillot, Using quaternions for coding 3D transformations, in: Andrew Glassner (Ed.), Graphics Gems, 1990. [18] Alan W. Paeth, A Fast Algorithm for General Raster Rotation, Proceedings, Graphics Interface ’86, Canadian Information Processing Society, Vancouver, 1986, pp. 77–81. [19] Alan W. Paeth, A fast algorithm for general raster rotation, in: Andrew Glassner (Ed.), Graphics Gems, 1990. [20] T.W. Parks, J.H. McClellan, Chebyshev approximation for nonrecursive digital filter design with linear phase, IEEE Transactions on Circuit Theory CT-19 (1972) 189–194. [21] P. Schro¨der, J.B. Salem, Fast rotation of volume data on parallel architectures, IEEE Visualization ’91 Proceedings, (1991) 50–57. [22] A. Tanaka, M. Kameyama, S. Kazama, O. Watanabe, A rotation method for raster image using skew transformation, Proceedings, IEEE Conference on Computer Vision and Pattern Recognition (1986) 272–277. [23] T. Toffoli, J. Quick, Three-dimensional rotations by three shears, Graphical Models and Image Processing 59(2) (1997) 89–95. [24] M. Unser, P. Thevenaz, L. Yaroslavsky, Convolution-based interpolation for fast, high-quality rotation of images, IEEE Transactions on Image Processing 4(10) (1995) 1371–1381. [25] C.M. Wittenbrink, A.K. Somani, Permutation warping for data parallel volume rendering, ACM SIGGRAPH Symposium on Parallel Rendering (1993) 57–60. [26] R.P. Woods, S.R. Cherry, J.C. Mazziotta, Rapid automated algorithm for aligning and reslicing PET images, Journal of Computer Assisted Tomography 16 (1992) 633.