- Email: [email protected]
Fast Fourier method for the accurate rotation of sampled images
15June
1997
OPTICS COMMUNICATIONS ELSEVIER
Optics Communications
139 (1997) 99-106
Full length article
Fast Fourier method for the accurate rotation of sampled images Kieran G. Larkin a,‘, Michael A. Oldfield b, Hanno Klemm a ADepartment of Physical Optics, School of Physics, The Uniuersin; of Sydney, NSW 2006, Australia ’ Australian Key Centre for Microscopy and Microanalysis, The Unioersity of Sydney, Sydney, NSW 2006, Australia Received 30 August 1996; accepted
12 February
1997
Abstract At present the best methods for rotation of discrete sampled images use a combination of (fast) Fourier interpolation followed by cubic interpolation onto a rotated grid. A method is presented which uses only Fourier interpolation. The new method has a similar computational complexity to the old, and is exactly reversible. The method uses the well-known decomposition of rotation into three pure shears. Each shear is performed using a 2D extension of the 1D Fourier shift theorem. This allows the fast Fourier transform (FFT) to be used. With appropriate data padding (such as zero padding) in both the real and Fourier domains, the procedure gives near perfect results and minimal loss of information in multiple rotation tests. Keywords: Fourier interpolation;
Digital image rotation; Lossless image transformation:
1. Introduction The literature regarding rotation of a two-dimensional sampled image (which nominally satisfies the criteria for a band-limited function) has often attempted to resolve some persistently inconsistent definitions. For example, the underlying test image may be a continuous sinusoidal grating truncated at a non-integer period [ 11. Such an image is not strictly band-limited and cannot, in general, be reconstructed using the sampling theorem. Rotation by multiples of 90” are straightforward enough, being merely row/column transpositions, but other angles present a number of difficulties. One useful approach to the problem is to express it as a remapping (or resampling) of one rectangular sampling grid onto a similar, but rotated, grid. Unfortunately this ignores the inevitable non-overlapping comer regions (see Fig. 1). This pervasive difficulty is at the heart of the doubtful definition of rotation for 2-D band-limited functions. If the preceding difficulty is ignored or, rather, elimi-
’ Present address: Canon Information Systems tralia. E-mail: [email protected]. 0030-4018/97/$17.00 Copyright PI1 SOO30-4018(97)00097-7
Research
Aus-
Image warping;
Image processing
algorithm
nated by appropriate (zero) padding of the image in the problematic comer regions, then some progress can be made towards an operational definition of rotation. Significant progress in overcoming various artifacts has been made by practitioners in the field [2,3]. Even so, a number of useful properties have not been exhibited by methods developed previously. The most important of these being the preservation of the image spectrum; it is well known that simplistic approaches to the problem, such as polynomial interpolation and cubic convolution attenuate high frequency components of the image [4]. This means that possible requirements such as (exact) reversibility and commutativity cannot be met. In short, some information is lost in the process. By taking a rather different approach which is firmly based on the Fourier information content of an image we have developed a procedure which is exactly reversible (no information loss) if not quite commutative. By formulating the problem as a series of one-dimensional shifts it has been possible to utilize the Fast Fourier Transform (FFI’) at all stages thus ensuring high computational efficiency. The procedure is useful for the rotation of datasets which contain significant information up to the sampling limit. Comparing the new method with other (3 stage shearing) rotation techniques we see that the
0 1997 Elsevier Science B.V. All rights reserved.
100
K.G. Lurkin rt d/Optics
Conmunicutions
139 (I9971 99-106
considered here because most sampled images are tacitly assumed to be band-limited. A two-dimensional sampled image is, implicitly, a doubly periodic function (in x and v) with values given only at discrete locations determined by a rectangular grid (square in all cases that we consider). Some implications of this periodicity for the representation of strictly non-periodic functions have been discussed by Fraser [7]. The sampling theorem (some readable accounts are available [8-l 11) prescribes the procedure for recovering values at any location in between grid points. Essentially this procedure is a convolution with a sine function. Conveniently this procedure can be applied in the Fourier domain; more of which will be explained in the next section. The operational definition of rotation (or any other resampling such as shifting, zooming, or generalised image warping [12]) can be expressed simply as the evaluation of underlying function values on the new grid by use of the sampling theorem. This definition is convenient and simple but has a few defects if the full sampling theorem is applied to the resampled dataset. Fig. l(a) shows the full implications for the most extreme case of interest here which is a 4.5” rotation. Using the full discrete Fourier transform (DFT) rather than the FFT it is possible to evaluate exactly all the underlying values on the new 45” grid but the comer regions are inevitably scrambled and the computation is infeasibly long for all but the smallest of images. If the new dataset is itself expected to satisfy all the sampling theorem requirements then the periodically repeating image in Fig. I(b) is the implied result. Clearly some comer information is missing and other information is duplicated. Thus, unless another not
(al
03
Fig. 1. (a) A periodic band-limited image showing the repeating structure. The central image is also shown in relation to a 45” rotated co-ordinate system (the black diamond shaped box). (b) After perfect resampling onto a 49 rotated grid the original image displays unwanted corner artifacts,
band-limited significant
spectrum computational
can
now
be preserved
without
a
penalty.
2. Possible definitions of discrete image rotation We need to be absolutely clear about what a sampled image represents before a useful definition of discrete image rotation can be obtained. For any meaningful interpretation a discrete image must represent an underlying continuous function; that is a function which exists at all possible intermediate values of the discretely sampled function. From hereon we assume that the underlying continuous image function satisfies the Shannon-Nyquist [5.6] sampling theorem. Other choices are possible but are
-Smallest square enclosing circular zone
Fig. The The that
2. The smallest meaningful rotation zones for a square image. circular zone encloses the image regardless of rotation angle. large square encloses the circular zone, and has an area twice of the original square.
K.G. Larkin et cd /Optics
Cornrnuni~ntior~.s 13Y f IY971 YY-106
suitable definition can be conceived, we are prompted to remove or otherwise ameliorate this serious defect. The problem is due to the obvious fact that a rotated square region does not exactly overlap the original square region unless the rotation is a multiple of 90”. A circular region, however, has the desired property of exact overlap for any angle of rotation. Unfortunately a circular region does not easily fit in with the original definition of a bandpass sampled image on a square grid. A possible remedy is to surround the single period image by a region of zeroes which extends outside an essentially circular region to the limits of an enlarged square region as shown in Fig. 2. In this paper we assume that the image has been padded by such a border of zeroes. There are a number of problems and solutions related to this zero padding which are discussed in Section 5. Our definition of rotation, then. is for a band-limited function, sampled on a square grid with a border of zeroes in the problematic comer and edge regions.
3. Image shear and the Fourier shift theorem The Fourier transform of a (2D) function which has undergone the geometric distortion known as the affine transform is itself distorted and scaled in a well defined manner [IO, 131. Rotation and shear are just special cases of the affine transformation. The decomposition of rotation into two separate scanline operations was first utilised for computer graphics by Catmull and Smith [14]. More recently the decomposition into three simple shears (without resealing) has been noted in a number of papers [ lS.l6]. As the component shears can be chosen to be purely in the x- and purely in the v-directions, a rotation can be conveniently split into three one-dimensional operations. This property has been put to good use in the fast (parallel) implementation of image rotation for computer graphics. The process is most clearly defined in terms of a matrix equation as follows:
0 ( .r’ ?.’ =
case -sin@
;~;j.(;)=M,-(;j.
(2)
A shear in the .x-direction is represented by M, or M’,. and similarly for a y-shear. The actual matrices can be shown to be of the following form [15,16]
’
1 tan(o/2)
i0
M’y=(_tan;0,2)
]
)’ Mvz( -ijn(jy)’
:).
(3)
Curiously, a simple application of the Fourier shift theorem to sheared images appears to have been overlooked until now, although a recent publication by Unser et al. [I71 develops a Fourier-based shift as the (infinite order) limit of a one-dimensional B-spline interpolation. The following analysis is initially presented for the case of continuous functions. Some important details necessary for a discrete formulation are presented in the Appendix. Consider a two-dimensional function f( x, y) which has a Fourier transform (or FT for short) F(u,t,) conventionally defined in the following equation
Typically the 2D FT is performed first by performing the ID FT on all the rows of an image, then by another ID Ff on all the resulting columns. We can define these intermediate transforms using a concise, yet informative, operator formalism. Firstly the Fourier transform along the x-direction (rows) is defined as the operator U( 1. F,(u,~)=U(f.(.r,v)}=/xf(.r,~)exp(-2~~~~~-)d.r. --z (5) Secondly the transform along the y-direction similarly defined as V{ }, F2(_~.~~)=V{f(x,~))
=/=
(columns)
is
,f(.x,y)exp(-2rripv)dv. -Zr. (6)
The following inverses
relations
apply to the transforms
V’{F(u,t,)]
=l_*IF(n,r)exp(?~i~,~)d~’ = u{.+..v)}.
M,=M;M;M,=M’;M’;M’~.
=
si;“)
and their
(1)
The original image co-ordinates are (x.,v). The new image co-ordinates are (x’,!‘) where the image has been rotated clockwise by an angle 0. which is equivalent to the new co-ordinate system at an angle 0 anticlockwise to the original system. The three component shears are as follows:
M
M’,=(:)
101
U’{F(u,L.)}
(7)
=JX F(u,L’)exp(2~iltx)drr -Z = V{f( _Y,Y)].
(8)
The simple effect of a shear on the ID Ff can now be demonstrated using the well-known Fourier shift theorem. Consider a .x-sheared function g( X, v 1 defined as g(.r.!‘)
=f(.r+av,v).
(9)
where the constant a controls the shear angle. The Ff in the direction of shear is as follows: G,(m,v)=U{g(x,p))=F,(u,v)exp(-2rriua?,) = exp(-2niUay)U{f(
x,.v)}.
(IO)
102
K.G. Lokin et d/Optics
In the operator salient features
formalism
U{j(x+ay,.v))
Communications
it is much easier to see the
=exp(-2~iuuy)U(f(x,y)).
(11)
Clearly the magnitude of the 1D FT of a sheared function is unchanged, whereas the phase is altered by a bilinear term (actually hyperboloidal in the two variables u and y). If a further transform (perpendicular to the first) is performed to obtain the full 2D fl then the following relation results from application of the shift theorem again, =F(u,v-au).
G(u,~~)=VU{f(x+uy,y)}
(12)
Eq. (12) simply illustrates the well-known result that a shear in the spatial domain is equivalent to a perpendicular shear in the Fourier domain.
4. Implementation
of the Fourier shift theorem
I39 (1997) 99-106
5. Edge effects and a&wing
Each shearing operation produces some rather disconcerting edge effects and aliasing problems in the discrete formulation, which are not present in the original continuous function analysis. In essence the problem arises from the implicit limited bandwidth and periodicity of the DFT. Fig. 3 shows the regions of a periodic function which are folded after the application of the 3 shears process. The regions affected in the worst case of 45” rotation are only marginally larger than the comer regions predicted for a single step rotation (compare Figs. 1 and 3). Interestingly these defects are completely reversible (unfolded) when the process is followed by the reverse (i.e. a - 0 rotation following a + 0 rotation). As far as we know, no other rotation method for discrete images has this reversible property. This reversibility was first noted by Unser et al. L171.
From the operator formalism developed in the preceding section we can assemble a full procedure for Fourier interpolated rotation for any angle 13(restricted in practice to a maximum of 45”, other angles being reached with additional transpositions of multiples of 90”). X-shear operation g.,( x,y)
= U-‘{exp(
-2 7rizLuy)~U{f(x,y)]}.
(a
(13)
Y-shear operation gJ
x.y)
= V-l{ exp( -27riubx)
.V{ g,( x,Y)}). (14)
X-shear again gs(x,.v)
=grJrv) =I_-‘{exp(-2rriuay).U{g,.,(x,y)}}.
(15)
The above procedure defines the three-step process in its entirety. So a rotation through an angle 0 can be accomplished by three 1D forward FFTs and three 1D inverse FFTs performed on 2D images. Suitable values of the shear parameters a and b are given by Eq. (3), specifically a = tan(0/2), and b = -sin0 for the XYX shear process. The process also requires three 2D complex multiplications by the phase factors exp(- 2n-iuay) and exp( - 2 rr iuk). These factors can be efficiently generated by a simple recursion often used for trigonometric evaluation in FFT algorithms:
(b
exp(iizo)=cosno+isinncr=c,+is,. c,+ I = c,c - s,s,
s”+I = S”C+ c,s.
(16)
The conversion of Eqs. (13) to (161 into their discrete form is trivial, bearing in mind the Nyquist constraint detailed in the Appendix. Comments relating to the computational complexity of the full procedure are delayed in the next section which deals with the practical solution of some edge effect problems raised in Section 2.
Fig. 3. (a) A simple image before the three shear rotation process. (b) After the three shear process the comer regions are incorrectly transformed. For the worst case of 45’ rotation the “scrambled” comers are shown shaded.
K.G. Larkin et al. / Optics Communications Aesthetically
and
practically
we may
wish
problem
6. Edge blending Although the edge effects we discuss here are well known to many researchers we believe they are worth mentioning because they are all too often ignored in the
103
Edge Effects with zero padding
to remove
regions. In the real space domain, for rotations in the range -45” < 0 5 45”, the wraparound comer regions are such that a zero padding (border regions around the original image) ensures that the wrapped regions are near zero intensity and consequently have little bearing of the overall rotated image. The wrapped regions are not exactly zero because the zero padding has the well known effect of introducing the Gibbs phenomenon [ 18,191 in the edge regions of the equivalent continuous image function (which is, after all, the function that this interpolation procedure selects). In situations where the discrete sampling theorem results in “ringing” boundary regions and this is counter to the intuitive assessment of what the edge regions should be, then appropriate edge blending methods may be introduced (see Section 6). For simplicity of implementation it is recommended that the zero padding increases the number of samples by a factor of 2 in each direction. Marginal gains in efficiency may be attained by utilising non-radix-2 FFTs but again these fine details will be omitted here. So, an N X N image zero padded to size 2 N X 2 N is our nominal starting point. Looking in the Fourier domain, this padding increases the number of sampled frequencies by a factor of 4. It can, therefore, be considered as a type of frequency interpolation. The maximum frequency sampled is unaltered; only the sampling interval is halved in both frequency coordinates. A shear of the padded image results in an orthogonal shear of the FT. This means that frequency wraparound (aliasing in other words) still occurs. To circumvent this wraparound we recommend zero padding in the Fourier domain as well. Again, a doubling in sample numbers in both axes is quite sufficient to circumvent the potential aliasing. The overall image size is now 4N X 4N. Edge blending in the Fourier domain is not recommended because it alters the original sample values. Any residual errors in this method of image rotation are due to subtle differences in overlap and reshuffling of these near zero zones both in real space and in frequency space. The computational requirements for this method are comparable to the method of Fraser [l] which requires of the order of (pNl’log,(pN)’ floating point operations. Typically p = 4 in the case of four times interpolation, however, p = 8 would be required for images with a zero padded border region. The Fourier shear method outlined in this paper requires order 3(4 N)‘log ,(4 N )z. If the initial image has N = 64. then floating point operation counts are in the range of one million to five million. these
139 (1997) 99-106
O.Ob
I
/
1-7
Fig. 4. A cross-section through the edge region of a uniform image that has been zero padded to twice its original width, then Fourier interpolated to reveal the underlying continuous function. The ripple structure decays to about 0.01 of the initial value once it reaches the boundary.
rotation literature. It should be stated that the zero padding in both domains is subject to the ubiquitous uncertainty principle. To paraphrase: a band-limited function must have infinite extent in the real space domain. Similarly a function with finite spatial support cannot be band-limited. The principle manifests itself rather subtly in this instance. Seemingly we have a function zero padded in both domains, but the zeroes only occur at our initial sampling points in the boundary regions. If we consider the continuous function that these samples actually represent in both domains we can see an oscillating function which just happens to pass through these zeros at certain well defined sample points. The continuous function in these regions is a doubly oscillating surface not unlike that found in the design of egg boxes! Fig. 4 shows the typical ripple structure on the edge of a zero padded region. The overshot/undershoot is about + 13% for a 64 pixel image padded up to 128 pixels. The Edge Effects with Simple (Mean ) Blending
Fig. 5. A cross-section through the edge region of a uniform image that has been zero padded to twice its original width. However, in this instance the first zero from the edge was replaced by a value midway between the edge value and zero. Such blending of the edge region reduces the ripple in the underlying continuous function. The ripple structure decays to about O.COO3 of the initial value once it reaches the boundary. The effect for 2D images is even more marked because the ripple in the difficult comer regions is now squared (O.OOfKrOOl).
KG. Larkin et al. /Optics
104
Communications
ripple structure dies out rather quickly reaching less than 1% at the extremes. By incorporating a smooth transition region from image to zero border the ripples can be substantially reduced. The process can be termed “blending” and is unlike (convolutional) smoothing because original image pixel values are unchanged; only the values in the padding region are modified. Blending can be considered as a polynomial fitting process. High order blending requires continuity of derivatives, however good results can be obtained with the very simplest blending method which inserts a mean pixel value between the last image pixel and the first zero pixel in the padding region. Fig. 5 shows the ripple structure for this form of blended edge. Note that the overshoot and undershoot is now only +3% with a quick decay down to less than 0.03% at the extremes.
7. Practical
example of rotation procedure
The starting point in this example is a 64 X 64 X 1 byte (8 bits) greyscale image of the Mona Lisa corresponding to the central zone of Fig. 1. First the image is simply zero padded in both real and Fourier space. Then the image is rotated by 45”. and then 45” again. Edge ringing artifacts are present in both images as expected from the sharp edge discontinuities in a band-limited function. As mentioned earlier the rotation procedure is exactly reversible and no is lost in this case except for information truncation/rounding errors due to the finite precision floating point arithmetic and, more importantly, quantization errors when the floating point value is converted back into a 1 byte per pixel representation. To assess errors due to multiple rotations the above case of + 45” rotation applied twice in succession was considered. The intermediate image is converted to 1 byte per pixel so we can expect even more quantization errors. The resultant image was compared with a 90” transposed image. One image was subtracted from the other and the resultant difference image was found to be identically zero in the unpadded (i.e. original) image region. The expectation even for perfect rotation would be a + 1 grey level difference due to two successive quantizations. Actually such an expectation is only likely for an image containing a uniform frequency spectrum.
8. Conclusion A simple and accurate method of rotating band-limited images has been presented. The method has a computational complexity comparable to other current procedures, but has the useful property that it is exactly reversible. Exact bounds on the accuracy have not yet been determined, however, a typical image has been rotated through two successive 45” with no measurable degradation.
I39 (1997) 99-106
The method can be extended to 3-D rotations using orthogonal shears in three directions, and has already been used for the accurate rotation of 3-D tomographic images
ml.
9. Note added in revision Since the manuscript was originally submitted our attention has been drawn to the excellent work of Unser et al. [ 171 by several colleagues and one of the reviewers. As mentioned in our revised text the work of Unser et al. develops a rotation algorithm using the Fourier interpolated shear. However, the emphasis of their work is upon the one-dimensional B-spline interpolation of Nth order. The Fourier method then appears as the limit as N tends to infinity. Whereas our approach is purely Fourier based and considers some of the broader issues regarding the awkward and inevitable edge effects implicit in discrete image rotation. Our initial simulations of the Unser method indicate that edge effects cause degradation in all parts of the image which is clearly visible in multiple rotation tests. With appropriate attention to edge effects our method shows no such degradation even after one hundred rotations. Clearly there is the usual trade-off between computational efficiency and accuracy, and this may be explored in further work. We would like to thank Jason Jones, Rebecca Powles and Matt Sheumack for their careful work simulating multiple rotations using these rotation algorithms.
Acknowledgements The authors would like to thank The Science Foundation in the School of Physics and the Australian Research Council for funding this work.
Appendix A The discrete Fourier shifr theorem Using discrete Fourier transforms for interpolation of functions has become a common technique for resampiing discrete data [21]. The interpolation takes N samples of a function, and produces a block of NM interpolants. Although it is often not explicitly stated, the resampling criterion is that the interpolated function is a bandwidth limited function passing through the sampled data points 1221. For interpolation at arbitrary points between the original sample points, either very large values for M are required, or further interpolations must be made between the resampled points. The Fourier shift theorem allows interpolation at arbitrary points between the sampled points
K.G. Larkin
et al./
Optics Convttunication.s
by translating the function by an arbitrary amount. The interpolation can be done in place, the interpolants at shifted positions replacing the original N sample values. More importantly, the interpolation can be carried out without loss of information, ensuring that interpolations can be performed serially in any order, while retaining all information contained in the original N sample points. The shifr theorem: Suppose that a function g,, is obtained from a function ,t, translated by an amount m to the right:
(Al)
R,, = .I;, ,,, If N’ is defined a\; if N is even, if N is odd.
(42)
the Fourier shift theorem for the discrete Fourier transform (DFT) states [8] that FL (the DFT of .fl and G, (the DFI of R) satisfy the following equation: e
G, =
-‘n,,r,h
/‘V
!i = 0.. . . N’/2 - I.
FL
e-la,lllo!L.V,,!v
i
F
I
!i = N’/2..
. N - I.
(A31
In general. nt can have an integer part and a fractional part. which we will write as tn, and /.Lrespectively so that rn = m,, + p. An equivalent way of interpreting the translated function g,, in Eq. (All is to consider g, as a set of interpolants on ,f;,. Eq. (A31 provides a means of calculating the DFI of R,, from the DFT of f,,. If ,f,, is a real valued function. its DFI will satisfy FL = FN-k.
k = 12..
.N’/2
- I,
(A41
where * signifies the complex conjugate. F, is real. and if N is even. FN,: is also real (if N is odd, FN,? does not exist). If G, is calculated using Eq. (A31 then its inverse transform will be
13
(IYY7)
105
YY-106
In order to obtain a real valued interpolating define the function efnpG ,,,,?
if N is even and k = N/2.
G,
otherwise.
G; =
function.
(A’)
With this definition, g:, the inverse DFT of G;, is a real valued interpolating function for .f;,. The phase of the N/2-harmonic differs from the “correct” interpolating function by up to 7r/2 radians (if ni,, is chosen to be the nearest integer to nr). but does retain the full power component of that harmonic. This choice for G; (and hence ~1,) also means that the interpolation is fully reversible: although some phase information has been discarded, it can be reconstructed by an appropriate interpolation on the newly constructed g:,. An alternative, and perhaps simpler way to view the complications arising at the Nyquist sampling frequency is as follows. Only real components occur at the Nyquist frequency for a real sampled function. All imaginary components are lost at the Nyquist frequency. We require a real, band-limited interpolation function for a meaningful definition of the discrete shift theorem. But a naive application of the shift theorem to this interpolation problem produces an imaginary component at the Nyquist frequency. Setting this imaginary component to zero avoids the problem. Interestingly the problem does not arise if the DFT is defined over a symmetric range of summation indices. It should be noted that if p= I/M, where M is an integer, the algorithm fully reproduces the interpolation values obtained using Fraser’s [21] algorithm. Numerical experiments show that interpolation using the Fourier shift theorem takes comparable time to interpolation with M = 2 using the Fraser algorithm, so if interpolation at midpoints between sampled data is required, there is little to choose between the two algorithms. As M increases, the Fraser algorithm requires increased storage space and longer computing time, whereas the shift algorithm described above can be done in place, and computation time is independent of the choice of interpolation position.
References I + z
I
G,, + x en”’ G_ V,/1’
(A9
From Eqs. tA3) and (A41 it is clear that G, = G,,mI.
X = 1.2 . . . . . N’/Z - I.
(A61
If N is odd. g,, is the sum of a set of complex conjugate pairs and a single real value. so g,, is real. The final term listed in Eq. (A5). G,,,. only exists if N is even. In that case, R,, is the sum of a set of complex conjugate pairs, a single real value and the value N-‘e”“‘G_,,, = N-1 eV,,?I-+II,,-H) F _v,2. This final term is complex if p # 0 and N is even. so that g,, is a complex function.
valued
[I] D. Fraser, Computer Vision Graphics Image Processmy -L6 (1989) 267. [2] D. Fraser, Aliasing in multipass image rotations. in: DICTA9 I. Digital Image Computing: Techniques and Applications. I99 I. Melbourne. [3] D. Fraser, R.A. Schowengerdt. IEEE Trans. Image Processing 3 (1994) 721. [4] S.K. Park, R.A. Schowengerdt, Computer Vision Graphic Image Processing 23 (1983) 258. [5] C.E. Shannon, Proc. Inst. Radio Eng. 37 t 1949) 10. [6] H. Nyquist, Trans. Amer. Inst. Electrical Eng. 2 (191-8) 617. [7] D. Fraser, Australian Computer J. 19 (1987) I 19. [8] R.N. Bracewell, The Fourier transform and its applications. McGraw-Hill Electrical and Electronic Engineering Series. Ed. F.E. Terman (McGraw Hill, New York. 1978).
106
K.G. Lnrkin et al. / Optics Communications 139 119971 99-106
[91 0. Brigham, The fast Fourier transform. 2nd Ed. (Prentice Hall, New Jersey, 1988). [lOIJ.D. Gaskill, Linear systems. Fourier transforms. and Optics (Wiley, New York, 1978). IllI A. Kohlenberg, J. Appl. Phys. 21 (1953) 1432. [I21 G. Wolberg, Digital Image Warping (Los Alamitos, California, IEEE Computer Press, 1992). Imaging, Signal ProcessiI31R.N. Bracewell, Two-Dimensional ing, Ed. A.V. Oppenheim (Prentice Hall. Englewood Cliffs, NJ, 199.5). iI41E. Catmull and A.R. Smith, Proc. SIGGRAPH’SO (1980) p. 279. [151A.W. Paeth, Graphics Interface ‘86 (1986) p. 77. [161A. Tanaka, M. Kameyama, S. Kazana and 0. Watanabe, IEEE Conf. on Computer Vision and Pattern Recognition (1986) p. 272.
[ 171 M. Unser. P. Thevenaz, L. Yaroslavsky, IEEE Trans. Image Processing 4 (1995) 1371. [18] R.W. Hamming, Digital Filters, Ed. A.V. Oppenheim &entice Hall, Englewood Cliffs, NJ, 1983). [19] G. Helmberg, J. Approximation 78 (1994) 41. [20] C.J. Cogswell, K.G. Larkin and H.U. Klemm, in: Three-Dimensional Microscopy: Image Acquisition and Processing III. IS and T/SPIE Symposium on Electronic Imaging Science and Technology, 1996, San Jose, CA. SPIE Proc. 7.655. p. 109. [21] D. Fraser, IEEE Trans. Acoustics Speech Signal Processing 37 (1989) 665. [22] R.W. Schafer and L.R. Rabiner. Proc. IEEE 61 (1973) 692.
1997
OPTICS COMMUNICATIONS ELSEVIER
Optics Communications
139 (1997) 99-106
Full length article
Fast Fourier method for the accurate rotation of sampled images Kieran G. Larkin a,‘, Michael A. Oldfield b, Hanno Klemm a ADepartment of Physical Optics, School of Physics, The Uniuersin; of Sydney, NSW 2006, Australia ’ Australian Key Centre for Microscopy and Microanalysis, The Unioersity of Sydney, Sydney, NSW 2006, Australia Received 30 August 1996; accepted
12 February
1997
Abstract At present the best methods for rotation of discrete sampled images use a combination of (fast) Fourier interpolation followed by cubic interpolation onto a rotated grid. A method is presented which uses only Fourier interpolation. The new method has a similar computational complexity to the old, and is exactly reversible. The method uses the well-known decomposition of rotation into three pure shears. Each shear is performed using a 2D extension of the 1D Fourier shift theorem. This allows the fast Fourier transform (FFT) to be used. With appropriate data padding (such as zero padding) in both the real and Fourier domains, the procedure gives near perfect results and minimal loss of information in multiple rotation tests. Keywords: Fourier interpolation;
Digital image rotation; Lossless image transformation:
1. Introduction The literature regarding rotation of a two-dimensional sampled image (which nominally satisfies the criteria for a band-limited function) has often attempted to resolve some persistently inconsistent definitions. For example, the underlying test image may be a continuous sinusoidal grating truncated at a non-integer period [ 11. Such an image is not strictly band-limited and cannot, in general, be reconstructed using the sampling theorem. Rotation by multiples of 90” are straightforward enough, being merely row/column transpositions, but other angles present a number of difficulties. One useful approach to the problem is to express it as a remapping (or resampling) of one rectangular sampling grid onto a similar, but rotated, grid. Unfortunately this ignores the inevitable non-overlapping comer regions (see Fig. 1). This pervasive difficulty is at the heart of the doubtful definition of rotation for 2-D band-limited functions. If the preceding difficulty is ignored or, rather, elimi-
’ Present address: Canon Information Systems tralia. E-mail: [email protected]. 0030-4018/97/$17.00 Copyright PI1 SOO30-4018(97)00097-7
Research
Aus-
Image warping;
Image processing
algorithm
nated by appropriate (zero) padding of the image in the problematic comer regions, then some progress can be made towards an operational definition of rotation. Significant progress in overcoming various artifacts has been made by practitioners in the field [2,3]. Even so, a number of useful properties have not been exhibited by methods developed previously. The most important of these being the preservation of the image spectrum; it is well known that simplistic approaches to the problem, such as polynomial interpolation and cubic convolution attenuate high frequency components of the image [4]. This means that possible requirements such as (exact) reversibility and commutativity cannot be met. In short, some information is lost in the process. By taking a rather different approach which is firmly based on the Fourier information content of an image we have developed a procedure which is exactly reversible (no information loss) if not quite commutative. By formulating the problem as a series of one-dimensional shifts it has been possible to utilize the Fast Fourier Transform (FFI’) at all stages thus ensuring high computational efficiency. The procedure is useful for the rotation of datasets which contain significant information up to the sampling limit. Comparing the new method with other (3 stage shearing) rotation techniques we see that the
0 1997 Elsevier Science B.V. All rights reserved.
100
K.G. Lurkin rt d/Optics
Conmunicutions
139 (I9971 99-106
considered here because most sampled images are tacitly assumed to be band-limited. A two-dimensional sampled image is, implicitly, a doubly periodic function (in x and v) with values given only at discrete locations determined by a rectangular grid (square in all cases that we consider). Some implications of this periodicity for the representation of strictly non-periodic functions have been discussed by Fraser [7]. The sampling theorem (some readable accounts are available [8-l 11) prescribes the procedure for recovering values at any location in between grid points. Essentially this procedure is a convolution with a sine function. Conveniently this procedure can be applied in the Fourier domain; more of which will be explained in the next section. The operational definition of rotation (or any other resampling such as shifting, zooming, or generalised image warping [12]) can be expressed simply as the evaluation of underlying function values on the new grid by use of the sampling theorem. This definition is convenient and simple but has a few defects if the full sampling theorem is applied to the resampled dataset. Fig. l(a) shows the full implications for the most extreme case of interest here which is a 4.5” rotation. Using the full discrete Fourier transform (DFT) rather than the FFT it is possible to evaluate exactly all the underlying values on the new 45” grid but the comer regions are inevitably scrambled and the computation is infeasibly long for all but the smallest of images. If the new dataset is itself expected to satisfy all the sampling theorem requirements then the periodically repeating image in Fig. I(b) is the implied result. Clearly some comer information is missing and other information is duplicated. Thus, unless another not
(al
03
Fig. 1. (a) A periodic band-limited image showing the repeating structure. The central image is also shown in relation to a 45” rotated co-ordinate system (the black diamond shaped box). (b) After perfect resampling onto a 49 rotated grid the original image displays unwanted corner artifacts,
band-limited significant
spectrum computational
can
now
be preserved
without
a
penalty.
2. Possible definitions of discrete image rotation We need to be absolutely clear about what a sampled image represents before a useful definition of discrete image rotation can be obtained. For any meaningful interpretation a discrete image must represent an underlying continuous function; that is a function which exists at all possible intermediate values of the discretely sampled function. From hereon we assume that the underlying continuous image function satisfies the Shannon-Nyquist [5.6] sampling theorem. Other choices are possible but are
-Smallest square enclosing circular zone
Fig. The The that
2. The smallest meaningful rotation zones for a square image. circular zone encloses the image regardless of rotation angle. large square encloses the circular zone, and has an area twice of the original square.
K.G. Larkin et cd /Optics
Cornrnuni~ntior~.s 13Y f IY971 YY-106
suitable definition can be conceived, we are prompted to remove or otherwise ameliorate this serious defect. The problem is due to the obvious fact that a rotated square region does not exactly overlap the original square region unless the rotation is a multiple of 90”. A circular region, however, has the desired property of exact overlap for any angle of rotation. Unfortunately a circular region does not easily fit in with the original definition of a bandpass sampled image on a square grid. A possible remedy is to surround the single period image by a region of zeroes which extends outside an essentially circular region to the limits of an enlarged square region as shown in Fig. 2. In this paper we assume that the image has been padded by such a border of zeroes. There are a number of problems and solutions related to this zero padding which are discussed in Section 5. Our definition of rotation, then. is for a band-limited function, sampled on a square grid with a border of zeroes in the problematic comer and edge regions.
3. Image shear and the Fourier shift theorem The Fourier transform of a (2D) function which has undergone the geometric distortion known as the affine transform is itself distorted and scaled in a well defined manner [IO, 131. Rotation and shear are just special cases of the affine transformation. The decomposition of rotation into two separate scanline operations was first utilised for computer graphics by Catmull and Smith [14]. More recently the decomposition into three simple shears (without resealing) has been noted in a number of papers [ lS.l6]. As the component shears can be chosen to be purely in the x- and purely in the v-directions, a rotation can be conveniently split into three one-dimensional operations. This property has been put to good use in the fast (parallel) implementation of image rotation for computer graphics. The process is most clearly defined in terms of a matrix equation as follows:
0 ( .r’ ?.’ =
case -sin@
;~;j.(;)=M,-(;j.
(2)
A shear in the .x-direction is represented by M, or M’,. and similarly for a y-shear. The actual matrices can be shown to be of the following form [15,16]
’
1 tan(o/2)
i0
M’y=(_tan;0,2)
]
)’ Mvz( -ijn(jy)’
:).
(3)
Curiously, a simple application of the Fourier shift theorem to sheared images appears to have been overlooked until now, although a recent publication by Unser et al. [I71 develops a Fourier-based shift as the (infinite order) limit of a one-dimensional B-spline interpolation. The following analysis is initially presented for the case of continuous functions. Some important details necessary for a discrete formulation are presented in the Appendix. Consider a two-dimensional function f( x, y) which has a Fourier transform (or FT for short) F(u,t,) conventionally defined in the following equation
Typically the 2D FT is performed first by performing the ID FT on all the rows of an image, then by another ID Ff on all the resulting columns. We can define these intermediate transforms using a concise, yet informative, operator formalism. Firstly the Fourier transform along the x-direction (rows) is defined as the operator U( 1. F,(u,~)=U(f.(.r,v)}=/xf(.r,~)exp(-2~~~~~-)d.r. --z (5) Secondly the transform along the y-direction similarly defined as V{ }, F2(_~.~~)=V{f(x,~))
=/=
(columns)
is
,f(.x,y)exp(-2rripv)dv. -Zr. (6)
The following inverses
relations
apply to the transforms
V’{F(u,t,)]
=l_*IF(n,r)exp(?~i~,~)d~’ = u{.+..v)}.
M,=M;M;M,=M’;M’;M’~.
=
si;“)
and their
(1)
The original image co-ordinates are (x.,v). The new image co-ordinates are (x’,!‘) where the image has been rotated clockwise by an angle 0. which is equivalent to the new co-ordinate system at an angle 0 anticlockwise to the original system. The three component shears are as follows:
M
M’,=(:)
101
U’{F(u,L.)}
(7)
=JX F(u,L’)exp(2~iltx)drr -Z = V{f( _Y,Y)].
(8)
The simple effect of a shear on the ID Ff can now be demonstrated using the well-known Fourier shift theorem. Consider a .x-sheared function g( X, v 1 defined as g(.r.!‘)
=f(.r+av,v).
(9)
where the constant a controls the shear angle. The Ff in the direction of shear is as follows: G,(m,v)=U{g(x,p))=F,(u,v)exp(-2rriua?,) = exp(-2niUay)U{f(
x,.v)}.
(IO)
102
K.G. Lokin et d/Optics
In the operator salient features
formalism
U{j(x+ay,.v))
Communications
it is much easier to see the
=exp(-2~iuuy)U(f(x,y)).
(11)
Clearly the magnitude of the 1D FT of a sheared function is unchanged, whereas the phase is altered by a bilinear term (actually hyperboloidal in the two variables u and y). If a further transform (perpendicular to the first) is performed to obtain the full 2D fl then the following relation results from application of the shift theorem again, =F(u,v-au).
G(u,~~)=VU{f(x+uy,y)}
(12)
Eq. (12) simply illustrates the well-known result that a shear in the spatial domain is equivalent to a perpendicular shear in the Fourier domain.
4. Implementation
of the Fourier shift theorem
I39 (1997) 99-106
5. Edge effects and a&wing
Each shearing operation produces some rather disconcerting edge effects and aliasing problems in the discrete formulation, which are not present in the original continuous function analysis. In essence the problem arises from the implicit limited bandwidth and periodicity of the DFT. Fig. 3 shows the regions of a periodic function which are folded after the application of the 3 shears process. The regions affected in the worst case of 45” rotation are only marginally larger than the comer regions predicted for a single step rotation (compare Figs. 1 and 3). Interestingly these defects are completely reversible (unfolded) when the process is followed by the reverse (i.e. a - 0 rotation following a + 0 rotation). As far as we know, no other rotation method for discrete images has this reversible property. This reversibility was first noted by Unser et al. L171.
From the operator formalism developed in the preceding section we can assemble a full procedure for Fourier interpolated rotation for any angle 13(restricted in practice to a maximum of 45”, other angles being reached with additional transpositions of multiples of 90”). X-shear operation g.,( x,y)
= U-‘{exp(
-2 7rizLuy)~U{f(x,y)]}.
(a
(13)
Y-shear operation gJ
x.y)
= V-l{ exp( -27riubx)
.V{ g,( x,Y)}). (14)
X-shear again gs(x,.v)
=grJrv) =I_-‘{exp(-2rriuay).U{g,.,(x,y)}}.
(15)
The above procedure defines the three-step process in its entirety. So a rotation through an angle 0 can be accomplished by three 1D forward FFTs and three 1D inverse FFTs performed on 2D images. Suitable values of the shear parameters a and b are given by Eq. (3), specifically a = tan(0/2), and b = -sin0 for the XYX shear process. The process also requires three 2D complex multiplications by the phase factors exp(- 2n-iuay) and exp( - 2 rr iuk). These factors can be efficiently generated by a simple recursion often used for trigonometric evaluation in FFT algorithms:
(b
exp(iizo)=cosno+isinncr=c,+is,. c,+ I = c,c - s,s,
s”+I = S”C+ c,s.
(16)
The conversion of Eqs. (13) to (161 into their discrete form is trivial, bearing in mind the Nyquist constraint detailed in the Appendix. Comments relating to the computational complexity of the full procedure are delayed in the next section which deals with the practical solution of some edge effect problems raised in Section 2.
Fig. 3. (a) A simple image before the three shear rotation process. (b) After the three shear process the comer regions are incorrectly transformed. For the worst case of 45’ rotation the “scrambled” comers are shown shaded.
K.G. Larkin et al. / Optics Communications Aesthetically
and
practically
we may
wish
problem
6. Edge blending Although the edge effects we discuss here are well known to many researchers we believe they are worth mentioning because they are all too often ignored in the
103
Edge Effects with zero padding
to remove
regions. In the real space domain, for rotations in the range -45” < 0 5 45”, the wraparound comer regions are such that a zero padding (border regions around the original image) ensures that the wrapped regions are near zero intensity and consequently have little bearing of the overall rotated image. The wrapped regions are not exactly zero because the zero padding has the well known effect of introducing the Gibbs phenomenon [ 18,191 in the edge regions of the equivalent continuous image function (which is, after all, the function that this interpolation procedure selects). In situations where the discrete sampling theorem results in “ringing” boundary regions and this is counter to the intuitive assessment of what the edge regions should be, then appropriate edge blending methods may be introduced (see Section 6). For simplicity of implementation it is recommended that the zero padding increases the number of samples by a factor of 2 in each direction. Marginal gains in efficiency may be attained by utilising non-radix-2 FFTs but again these fine details will be omitted here. So, an N X N image zero padded to size 2 N X 2 N is our nominal starting point. Looking in the Fourier domain, this padding increases the number of sampled frequencies by a factor of 4. It can, therefore, be considered as a type of frequency interpolation. The maximum frequency sampled is unaltered; only the sampling interval is halved in both frequency coordinates. A shear of the padded image results in an orthogonal shear of the FT. This means that frequency wraparound (aliasing in other words) still occurs. To circumvent this wraparound we recommend zero padding in the Fourier domain as well. Again, a doubling in sample numbers in both axes is quite sufficient to circumvent the potential aliasing. The overall image size is now 4N X 4N. Edge blending in the Fourier domain is not recommended because it alters the original sample values. Any residual errors in this method of image rotation are due to subtle differences in overlap and reshuffling of these near zero zones both in real space and in frequency space. The computational requirements for this method are comparable to the method of Fraser [l] which requires of the order of (pNl’log,(pN)’ floating point operations. Typically p = 4 in the case of four times interpolation, however, p = 8 would be required for images with a zero padded border region. The Fourier shear method outlined in this paper requires order 3(4 N)‘log ,(4 N )z. If the initial image has N = 64. then floating point operation counts are in the range of one million to five million. these
139 (1997) 99-106
O.Ob
I
/
1-7
Fig. 4. A cross-section through the edge region of a uniform image that has been zero padded to twice its original width, then Fourier interpolated to reveal the underlying continuous function. The ripple structure decays to about 0.01 of the initial value once it reaches the boundary.
rotation literature. It should be stated that the zero padding in both domains is subject to the ubiquitous uncertainty principle. To paraphrase: a band-limited function must have infinite extent in the real space domain. Similarly a function with finite spatial support cannot be band-limited. The principle manifests itself rather subtly in this instance. Seemingly we have a function zero padded in both domains, but the zeroes only occur at our initial sampling points in the boundary regions. If we consider the continuous function that these samples actually represent in both domains we can see an oscillating function which just happens to pass through these zeros at certain well defined sample points. The continuous function in these regions is a doubly oscillating surface not unlike that found in the design of egg boxes! Fig. 4 shows the typical ripple structure on the edge of a zero padded region. The overshot/undershoot is about + 13% for a 64 pixel image padded up to 128 pixels. The Edge Effects with Simple (Mean ) Blending
Fig. 5. A cross-section through the edge region of a uniform image that has been zero padded to twice its original width. However, in this instance the first zero from the edge was replaced by a value midway between the edge value and zero. Such blending of the edge region reduces the ripple in the underlying continuous function. The ripple structure decays to about O.COO3 of the initial value once it reaches the boundary. The effect for 2D images is even more marked because the ripple in the difficult comer regions is now squared (O.OOfKrOOl).
KG. Larkin et al. /Optics
104
Communications
ripple structure dies out rather quickly reaching less than 1% at the extremes. By incorporating a smooth transition region from image to zero border the ripples can be substantially reduced. The process can be termed “blending” and is unlike (convolutional) smoothing because original image pixel values are unchanged; only the values in the padding region are modified. Blending can be considered as a polynomial fitting process. High order blending requires continuity of derivatives, however good results can be obtained with the very simplest blending method which inserts a mean pixel value between the last image pixel and the first zero pixel in the padding region. Fig. 5 shows the ripple structure for this form of blended edge. Note that the overshoot and undershoot is now only +3% with a quick decay down to less than 0.03% at the extremes.
7. Practical
example of rotation procedure
The starting point in this example is a 64 X 64 X 1 byte (8 bits) greyscale image of the Mona Lisa corresponding to the central zone of Fig. 1. First the image is simply zero padded in both real and Fourier space. Then the image is rotated by 45”. and then 45” again. Edge ringing artifacts are present in both images as expected from the sharp edge discontinuities in a band-limited function. As mentioned earlier the rotation procedure is exactly reversible and no is lost in this case except for information truncation/rounding errors due to the finite precision floating point arithmetic and, more importantly, quantization errors when the floating point value is converted back into a 1 byte per pixel representation. To assess errors due to multiple rotations the above case of + 45” rotation applied twice in succession was considered. The intermediate image is converted to 1 byte per pixel so we can expect even more quantization errors. The resultant image was compared with a 90” transposed image. One image was subtracted from the other and the resultant difference image was found to be identically zero in the unpadded (i.e. original) image region. The expectation even for perfect rotation would be a + 1 grey level difference due to two successive quantizations. Actually such an expectation is only likely for an image containing a uniform frequency spectrum.
8. Conclusion A simple and accurate method of rotating band-limited images has been presented. The method has a computational complexity comparable to other current procedures, but has the useful property that it is exactly reversible. Exact bounds on the accuracy have not yet been determined, however, a typical image has been rotated through two successive 45” with no measurable degradation.
I39 (1997) 99-106
The method can be extended to 3-D rotations using orthogonal shears in three directions, and has already been used for the accurate rotation of 3-D tomographic images
ml.
9. Note added in revision Since the manuscript was originally submitted our attention has been drawn to the excellent work of Unser et al. [ 171 by several colleagues and one of the reviewers. As mentioned in our revised text the work of Unser et al. develops a rotation algorithm using the Fourier interpolated shear. However, the emphasis of their work is upon the one-dimensional B-spline interpolation of Nth order. The Fourier method then appears as the limit as N tends to infinity. Whereas our approach is purely Fourier based and considers some of the broader issues regarding the awkward and inevitable edge effects implicit in discrete image rotation. Our initial simulations of the Unser method indicate that edge effects cause degradation in all parts of the image which is clearly visible in multiple rotation tests. With appropriate attention to edge effects our method shows no such degradation even after one hundred rotations. Clearly there is the usual trade-off between computational efficiency and accuracy, and this may be explored in further work. We would like to thank Jason Jones, Rebecca Powles and Matt Sheumack for their careful work simulating multiple rotations using these rotation algorithms.
Acknowledgements The authors would like to thank The Science Foundation in the School of Physics and the Australian Research Council for funding this work.
Appendix A The discrete Fourier shifr theorem Using discrete Fourier transforms for interpolation of functions has become a common technique for resampiing discrete data [21]. The interpolation takes N samples of a function, and produces a block of NM interpolants. Although it is often not explicitly stated, the resampling criterion is that the interpolated function is a bandwidth limited function passing through the sampled data points 1221. For interpolation at arbitrary points between the original sample points, either very large values for M are required, or further interpolations must be made between the resampled points. The Fourier shift theorem allows interpolation at arbitrary points between the sampled points
K.G. Larkin
et al./
Optics Convttunication.s
by translating the function by an arbitrary amount. The interpolation can be done in place, the interpolants at shifted positions replacing the original N sample values. More importantly, the interpolation can be carried out without loss of information, ensuring that interpolations can be performed serially in any order, while retaining all information contained in the original N sample points. The shifr theorem: Suppose that a function g,, is obtained from a function ,t, translated by an amount m to the right:
(Al)
R,, = .I;, ,,, If N’ is defined a\; if N is even, if N is odd.
(42)
the Fourier shift theorem for the discrete Fourier transform (DFT) states [8] that FL (the DFT of .fl and G, (the DFI of R) satisfy the following equation: e
G, =
-‘n,,r,h
/‘V
!i = 0.. . . N’/2 - I.
FL
e-la,lllo!L.V,,!v
i
F
I
!i = N’/2..
. N - I.
(A31
In general. nt can have an integer part and a fractional part. which we will write as tn, and /.Lrespectively so that rn = m,, + p. An equivalent way of interpreting the translated function g,, in Eq. (All is to consider g, as a set of interpolants on ,f;,. Eq. (A31 provides a means of calculating the DFI of R,, from the DFT of f,,. If ,f,, is a real valued function. its DFI will satisfy FL = FN-k.
k = 12..
.N’/2
- I,
(A41
where * signifies the complex conjugate. F, is real. and if N is even. FN,: is also real (if N is odd, FN,? does not exist). If G, is calculated using Eq. (A31 then its inverse transform will be
13
(IYY7)
105
YY-106
In order to obtain a real valued interpolating define the function efnpG ,,,,?
if N is even and k = N/2.
G,
otherwise.
G; =
function.
(A’)
With this definition, g:, the inverse DFT of G;, is a real valued interpolating function for .f;,. The phase of the N/2-harmonic differs from the “correct” interpolating function by up to 7r/2 radians (if ni,, is chosen to be the nearest integer to nr). but does retain the full power component of that harmonic. This choice for G; (and hence ~1,) also means that the interpolation is fully reversible: although some phase information has been discarded, it can be reconstructed by an appropriate interpolation on the newly constructed g:,. An alternative, and perhaps simpler way to view the complications arising at the Nyquist sampling frequency is as follows. Only real components occur at the Nyquist frequency for a real sampled function. All imaginary components are lost at the Nyquist frequency. We require a real, band-limited interpolation function for a meaningful definition of the discrete shift theorem. But a naive application of the shift theorem to this interpolation problem produces an imaginary component at the Nyquist frequency. Setting this imaginary component to zero avoids the problem. Interestingly the problem does not arise if the DFT is defined over a symmetric range of summation indices. It should be noted that if p= I/M, where M is an integer, the algorithm fully reproduces the interpolation values obtained using Fraser’s [21] algorithm. Numerical experiments show that interpolation using the Fourier shift theorem takes comparable time to interpolation with M = 2 using the Fraser algorithm, so if interpolation at midpoints between sampled data is required, there is little to choose between the two algorithms. As M increases, the Fraser algorithm requires increased storage space and longer computing time, whereas the shift algorithm described above can be done in place, and computation time is independent of the choice of interpolation position.
References I + z
I
G,, + x en”’ G_ V,/1’
(A9
From Eqs. tA3) and (A41 it is clear that G, = G,,mI.
X = 1.2 . . . . . N’/Z - I.
(A61
If N is odd. g,, is the sum of a set of complex conjugate pairs and a single real value. so g,, is real. The final term listed in Eq. (A5). G,,,. only exists if N is even. In that case, R,, is the sum of a set of complex conjugate pairs, a single real value and the value N-‘e”“‘G_,,, = N-1 eV,,?I-+II,,-H) F _v,2. This final term is complex if p # 0 and N is even. so that g,, is a complex function.
valued
[I] D. Fraser, Computer Vision Graphics Image Processmy -L6 (1989) 267. [2] D. Fraser, Aliasing in multipass image rotations. in: DICTA9 I. Digital Image Computing: Techniques and Applications. I99 I. Melbourne. [3] D. Fraser, R.A. Schowengerdt. IEEE Trans. Image Processing 3 (1994) 721. [4] S.K. Park, R.A. Schowengerdt, Computer Vision Graphic Image Processing 23 (1983) 258. [5] C.E. Shannon, Proc. Inst. Radio Eng. 37 t 1949) 10. [6] H. Nyquist, Trans. Amer. Inst. Electrical Eng. 2 (191-8) 617. [7] D. Fraser, Australian Computer J. 19 (1987) I 19. [8] R.N. Bracewell, The Fourier transform and its applications. McGraw-Hill Electrical and Electronic Engineering Series. Ed. F.E. Terman (McGraw Hill, New York. 1978).
106
K.G. Lnrkin et al. / Optics Communications 139 119971 99-106
[91 0. Brigham, The fast Fourier transform. 2nd Ed. (Prentice Hall, New Jersey, 1988). [lOIJ.D. Gaskill, Linear systems. Fourier transforms. and Optics (Wiley, New York, 1978). IllI A. Kohlenberg, J. Appl. Phys. 21 (1953) 1432. [I21 G. Wolberg, Digital Image Warping (Los Alamitos, California, IEEE Computer Press, 1992). Imaging, Signal ProcessiI31R.N. Bracewell, Two-Dimensional ing, Ed. A.V. Oppenheim (Prentice Hall. Englewood Cliffs, NJ, 199.5). iI41E. Catmull and A.R. Smith, Proc. SIGGRAPH’SO (1980) p. 279. [151A.W. Paeth, Graphics Interface ‘86 (1986) p. 77. [161A. Tanaka, M. Kameyama, S. Kazana and 0. Watanabe, IEEE Conf. on Computer Vision and Pattern Recognition (1986) p. 272.
[ 171 M. Unser. P. Thevenaz, L. Yaroslavsky, IEEE Trans. Image Processing 4 (1995) 1371. [18] R.W. Hamming, Digital Filters, Ed. A.V. Oppenheim &entice Hall, Englewood Cliffs, NJ, 1983). [19] G. Helmberg, J. Approximation 78 (1994) 41. [20] C.J. Cogswell, K.G. Larkin and H.U. Klemm, in: Three-Dimensional Microscopy: Image Acquisition and Processing III. IS and T/SPIE Symposium on Electronic Imaging Science and Technology, 1996, San Jose, CA. SPIE Proc. 7.655. p. 109. [21] D. Fraser, IEEE Trans. Acoustics Speech Signal Processing 37 (1989) 665. [22] R.W. Schafer and L.R. Rabiner. Proc. IEEE 61 (1973) 692.