# Digital image warping

Download 2.54 Mb.
 Page 15/30 Date 18.10.2016 Size 2.54 Mb.
 f (x) fo f l f Xo X 1 X 2 X 3 X 4 X 5 X 6 Figure 5.9: A spline consisting of 6 piecewise cubic polynomials. The four coefficients offt can be defined in terms of the data points and their first (or second) derivatives. Assuming that the data samples are on the integer lattice, each 134 IMAGE RESAMPLING spaced one unit apart, then the coefficients, defined in terms of the data samples and their first derivatives, are given below. a0 = Y (5.4.16a) a = y (5.4.16b) a2 = 3Ay, - 2y -Y,+I (5.4.16c) a3 = -2Ay +y +Y+I (5.4.16d) where Ay = y+ - y. Although the derivatives are not supplied with the data, they are derived by solving the following system of linear equations. 41 141 y y 1 4 _ 4 ,_ -Sy0 + 4yl +Y2 30'2 -Y0) 30'3 -Yl) 30'n-t - Yn-3) -Y-3 - 4yn-2 + 5yn- (5.4.17) The not-a-knot boundary condition [de Boor 78] was used above, as reflected in the first and last rows of the matrices. It is superior to the artificial boundary conditions com- monly reported in the literature, such as the natural or cyclic end conditions, which have no relevance in our application. Note that the need to solve a linear system of equations arises from global dependencies introduced by the constxaints for continuous first and second derivatives at the knots. A complete derivation is given in Appendix 2. In order to compare interpolating cubic splines with other methods, we must analyze the interpolation kernel. Thus far, however, the piecewise interpolating polyno- mials have been derived without any reference to an interpolation kernel. We seek to express the interpolating cubic spline as a convolution in a manner similar to the previous algorithms. This can be done with the use of cubic B-splines as interpolation kernels [Hou 78]. 5.4.5.1. B-Splines A B-spline of degree n is derived through n convolutions of the box filter, B 0. Thus, B t =B0*B0 denotes a B-spline of degree 1, yielding the familiar triangle filter shown in Fig. 5.7a. Interpolation by B  consists of a sequence of stxaight lines joined at the knots continuously. This is equivalent to linear interpolation. S.4 INTERPOLATION KERNELS 13 The second-degree B-spline B2 is produced by convolving Bo*B 1. Using B2 to interpolate data yields a sequence of parabolas that join at the knots continuously together with their slopes. The span orB2 is limited to three points. The cubic B-spline B 3 is generated from convolving Bo*B2. That is, B 3 = Bo*Bo*Bo*B o. The interpolation with B 3 is composed of a series of cubic poly- nomials that join at the knots continuously together with their slopes and curvatures, i.e., their first and second derivatives. Figure 5.10 summarizes the shapes of these low-order B-splines. -l.5 -.5 .5 1.5 -2 -1 0 1 2 Figure 5.10: Low-order B-splines are derived from repeated box filters. Denoting the cubic B-spline interpolation kernel as h, we have the following piece- wise cubic polynomials defining the kemel. [31xl3-6lxl 2+4 0-< Ixl < l h(x) = -}lx13+61x12-121xl+8 l 2_< Ixl This kemel is sometimes called the Parzen window. There are several properties of cubic B-splines worth noting. As in the cubic con- volution method, the extent of the cubic B-spline is over four points. This allows two points on each side of the centxal interpolated region to be used in the convolution. Con- sequently, the cubic B-spline is shift-invariant as well. Unlike cubic convolution, however, the cubic B-spline kernel is not interpolatory since it does not satisfy the necessary consmint that h (0)= 1 and h(1)= h(2)= 0. Instead, it is an approximating function that passes near the points but not necessarily through them. This is due to the fact that the kernel is strictly positive. 136 IMAGE RESAMPLING The posifivity of the cubic B-spline kernel is actually attractive for our image pro- cessing application. When using kernels with negative lobes, (e.g., the cubic convolution and windowed sinc functions), it is possible to generate negative values while interpolat- ing positive data. Since negative intensity values are meaningless for display, it is desir- able to use strictly positive interpolation kernels to guarantee the positivity of the interpo- lated image. There are problems, however, in dkectly interpolating the data with kernel h, as given in Eq. (5.4.18). Due to the low-pass (blur) characteristics of h, the image under- goes considerable smoothing. This is evident by examining its frequency response where the stopband is effectively suppressed at the expense of additional attenuation in the passband. This leads us to the development of an interpolation method built upon the local support of the cubic B-spline. 5.4.5.2. Interpolating B-Splines Interpolating with cubic B-splines requires that at data point x/, we again satisfy Eq. (5.4.7). Namely, j+2 f(xj) =  ch(xj-x,) (5.4.19) k=j-2 From Eq. (5.4.18), we have h(0)=4/6, h(-1)= h(1)= 1/6, and h(-2)= h(2)= 0. This yields 1 f (xj) = (cj_ 1 q- 4cj + Cj+l) (5.4.20) Since this must be true for all data points, we have a chain of global dependencies for the ck coefficients. The resulting linear system of equations is similar to that obtained for the derivatives of the cubic interpolating spline algorithm. We thus have, f0 [4 1 fl i41 f2 141 C2 (5.4.21) Labeling the three matrices above as F, K, and C, respectively, we have F = K C (5.4.22) The coefficients in C may be evaluated by multiplying the known data points F with the inve[se of the tridiagonal matrix K. C = K - F (5.4.23) 5.4 INTERPOLATION KERNELS 137 The inversion of tridiagonal matrix K has an efficient'algorithm that is solvable in linear time [Press 88]. In [Lee 83], the matrix inversion step is modified to introduce high-frequency emphasis. This serves to compensate for the undesirable low-pass filter imposed by the point-spread function of the imaging system. In all the previous methods, the coefficients c were taken to be the data samples themselves. In the cubic spline interpolation algorithm, however, the coefficients must be determined by solving a tridiagonal matrix problem. After the interpolation coefficients have been computed, cubic spline interpolation has the same computational cost as cubic convolution. 5.4.6. Windowed Sinc Function Sampling theory establishes that the sine function is the ideal interpolation kernel. Although this interpolation filter is exact, it is not practical since it is an IIR filter defined by a slowly converging infinite sum. Nevertheless, it is perfectly reasonable to consider the effects of using a trancated, and therefore finite, sinc function as the interpolation ker- nel. The results of this operation are predicted by sampling theory, which demonstrates that huncation in one domain leads to ringing in the other domain. This is due to the fact that truncating a signal is equivalent to multiplying it with a rectangle function Rect(x), defined as Rect(x) = .5 -< Ix l (5.4.24) Since multiplication in one domain is convolution in the other, lynncation amounts to convolving the signal's spectram with a sinc function, the transform pair ofRect (x). We have already seen an example of this in Fig. 4.7. Since the stopband is no longer elim- inated, but rather attenuated by a ringing filter (i.e., a sinc), the input is not bandlimited and aliasing artifacts are introduced. The most typical problems occur at step edges, where the Gibbs phenomena becomes noticeable in the form of undershoots, overshoots, and ringing in the vicinity of edges. In [Ratzel 80], the author found this method to per- form poorly. The Rect function above served as a window, or kemel, that weighs the input signal. In Fig. 5.11a, we see the Rect window extended over three pixels on each side of its center, i.e., Rect(6x) is plotted. The corresponding windowed sinc function h(x) is shown in Fig. 5.1lb. This is simply the product of the sine function with the window function, i.e., sinc(x)Rect(6x). Its spectrum, shown in Fig. 5.11c, is nearly an ideal low-pass filter. Although it has a fairly sharp mmsifion from the passband to the stop- band, it is plagued by ringing. In order to more clearly see the values in the spectrum, we use a logarithmic scale for the vertical axis of the spectram in Fig. 5.11 d. The next few figures will be illustrated by using this same four-part format. Ringing can be mitigated by using a different windowing function exhibiting smoother fall-off than the rectangle. The resulting windowed sine function can yield 138 IMAGE REVAMPLING S.4 INTERPOLATION KERNELS 139 Rect (x) -.2s L--.i..-...!-....--i...-.-i.-.-..i.-.-..!...-...::......i.....-:: ] 4-3-2-101234 4-33-101234 (a) ( -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 (C) (d) Figure 5.11: (a) Rectangular window; (b) Windowed sinc; (c) Spectrum; (d) Log plot. better results. However, since slow fall-off requires larger windows, the computation remains costly. Aside from the rectangular window mentioned above, the most frequently used win- dow functions are: Harm, t Hamming, Blackman, and Kaiser [Antoniou 79]. These filters identify a quantity known as the ripple ratio, defined as the ratio of the maximum side- lobe amplitu.d.e to the main-lobe amplitude. Good filters will have small ripple ratios to achieve effective attenuation in the stopband. A tradeoff exists, however, between ripple ratio and main-lobe width. Therefore, as the ripple ratio is decreased, the main-lobe width is increased. This is consistent with the reciprocal relationship between the spatial and frequency domains, i.e., narrow bandwidths correspond to wide spatial functions. In general, though, each of these smooth window functions is defined over a small finite extent. This is tantamount to multiplying the smooth window with a rectangle function, While this is better than the Rect function alone, there will inevitably be some form of aliasing. Nevertheless, the window functions described below offer a good compromise between tinging and blurring. ' Due to Julius yon Harm. It is often mistakenly referred to as the Htinning window. 5.4.6.1. Hann and Hamming Windows The Hann and Hamming windows are defined as { 2for N- 1 +(1-a)cos2 i- Ixl < 2 HannlHamming(x) = - (5.4.25) otherwise where N is the number of samples in the windowing function. The two windowing func- tions differ in the choice of x. In the Hann window x=0.5, and in the Hamming window 0=0.54. Since they both amount to a scaled and shifted cosine function, they are also known as the raised cosine window. The spectra for the Hann and Hamming windows can be shown to be the sum of a sinc, the spectrum of Rect(x), with two shifted counterparts: a sinc shifted to the tight by 2x/(N - 1), as well as one shifted to the left by the same amount. This serves to cancel the right and left side lobes in the specmm of Rect(x). As a result, the Hann and Ham- ming windows have reduced side lobes in their spectra as compared to those of the rec- tangular window. The Hann window is illustrated in Fig. 5.12. -4 -3 -2 -I 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 (a) (b) IH(f)l IHCf)l 0.0011 le-04 I le-05 le-06 -4-3-2-1 0 1 2 3 4 -4-3-2-1 0 1 2 3 4 (C) (d) Figure 5.12: (a) Hann window; (b) Windowed sinc; (c) Spectram; (d) Log plot. IIq] I  140 IMAGE RESAMPL1NG Notice that the passband is only slightly attenuated, but the stopband continues to retain high frequency components in the stopband, albeit less than that of Rect(x). It performs somewhat better in the stopband than the Hamming window, as shown in Fig. 5.13. This is partially due to the fact that the Hamming window is discontinuous at its ends, giving rise to "kinks" in the spectrum. -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -I 0 1 2 3 4 (a) (b) IHCf)[ .7s ...h....;..-..L....h., ,..i.......i......L....i-. IHCf)} 0.0011 le-04 I le05 le-06 -4-3-2-101234 -4-3-2-101234 ½) Figure 5.13: (a) Hamming window; (b) Windowed sinc; (c) Spectrum; (d) Log plot. 5.4.6.2. Blackman Window The Blackman window is similar to the Hann and Hamming windows. It is defined as i 2rx 4rx N - 1 ß 42+0.5cos----T+0.08cos---ZT Ix{ < 2 Blackman (x) = - - (5.4.26) otherwise The purpose of the additional cosine term is to further reduce the ripple ratio. This win- dow function is shown in Fig. 5.14. 5.4 INTERPOLATION KERNELS 141 Blackman (x ) -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 (a) (b) IHtf)l 0.1 0.01 o.ol ] 1-04 [ 1½-0 ] le-06 L -4-3-2-1 0 1 2 3 4 -4-3-2-1 0 1 2 3 4 (C) (d) Figure 5.14: (a) Blackman window; (b) Windowed sinc; (c) Spectrum; (d) Log plot. 5.4.6.3. Kaiser Window The Kaiser window is defined as f 1o(b) N- 1 Kaiser(x) =l: (0 'xl< 2 (5.4.27) otherwise f f 211t' I0 is the zeroth-order Bessel function of the first kind. This can be evaluated to any desired degree of accuracy by using the rapidly converging series Io(n) = 1 +  (5.4.29) k=l 142 IMAGE RESAMPLING The Kaiser window leaves the filter designer much flexibility in controlling the rip- ple ratio by adjusting the parameter at. As at is incremented, the level of sophistication of the window function grows as well. Therefore, the rectangular window corresponds to a Kaiser window with at = 0, while more sophisticated windows such as the Hamming win- dow correspond to o = 5. This formulation facilitates a tradeoff between ringing and edge softening. 5.4.6.4. Lanczos Window Windowed sinc functions are notorious for producing ringing artifacts near edges. Although they are an improvement over truncated sinc functions, they retain a fairly sharp transition from passband to stopband. Superior filters can be designed by imposing further constraints on the filter response in the frequency domain. Reasonable constraints to impose on the kernel include: unity gain in the low-pass region with cut-off at frequency fi, zero gain at high frequencies beyond f2, and linear fall-off in the transition range between fi and f2. This frequency response can be expressed as the convolution of two boxes. In the spatial domain, this corresponds to the multiplication of two sinc functions, yielding a function known as the Lanczos window. The widths of the two sinc functions determine the extent of the transition range. The two-lobed Lanczos window function is defined as sin(;c/2) 0 -< Ix I < 2 Lanczos 2(x) = :x/2 (5.4.30) 0 2- The Lanczos2 window function is the central lobe of a sinc function. It is wide enough to extend over two lobes of the ideal low-pass filter, i.e., a second sinc function. The windowed sinc function is therefore given by the product sinc (x)Latczos2(x). This can be rewritten as sinc (x) sinc (x/2) Rect(x/4), where the first term is the ideal low-pass filter, the second term is Lanczos 2(x), and Rect(x/4) is the rectangular function that tan- cares Lanczos2 past x=2. Note that its abscissa is x/4 because Rect is defined over -.5 < x < .5. The spectrum of this product is Rect(f)*Rect(2f)*sinc (4f), where * is convolution. The Lanczos 2(x) window function is shown in Fig. 5.15. This formulation can be generalized to an N-lobed window function by replacing the value 2 in Eq. (5.4.30) to the value N. For instance, the 3-lobed Lanczos window is defined as sin(mr/3) 0 -< I x I < 3 Lanczos3(x) = rx/3 (5.4.31) 0 3- The Lanczos3(x) window function is shown in Fig. 5.16. As we let more more lobes pass under the Lanczos window, then the spectrum of the windowed sinc function becomes Rect(f)*Rect(Nf)*sinc(2Nf). This proves to be a superior frequency 5.4 INTERPOLATION KERNELS 143 Irtf)l 0.1 0.01 [ le-04] le,-05 [ -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 (C) (d) Figure 5.15: (a) Lanczos2 window; (b) Windowed sinc; (c) Spectrum; (d) Log plot. response than that of the 2-lobed Lanczos window because the Rect(Nf) term causes faster fall-off in the transition region and sinc (2N f) is a narrower sinc function that pro- duces less deleterious ringing artifacts. 5.4.6.5. Gaussian Window The Gaussian function is defined as 1 e_X2/2o2 (5.4.32) Gauss(x) = 2q'o where o is the standard deviation. Ganssians have the nice property that their spectrum is also a Gaussian. They can be used to directly smooth the input for prefihering pur- poses or to smooth the sinc function for windowed sinc reconstruction filters. Further- more, since the tails of a Gaussian diminish rapidly, they may be truncated and still pro- duce results that are not plagued by excessive ringing. The rate of fall-off is determined by o, with low values of  resulting in faster decay. The general form of the Gaussian may be expressed more conveniently as Gausso(x) = 2 -(x)2 (5.4.33) 144 IMAGE RESAMPLING -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 (a) (b) IHf)l 0.1 0.01 0.e01 1½-04 lo-05 lo-06 -4-3-2-1 0 1 2 3 4 -4-3-2-1 0 1 2 3 4 (C) (d) Figure 5.16: (a) Lanczos3 window; (b) Windowed sinc; (c) Spectram; (d) Log plot. Two plots of this function are shown in Fig. 5.17 with = 1/2 and 1/'. The latter is a wider Gaussian than the first, and its magnitude doesn't become negligible until two sam- ples away from the center. The benefit of these choices of ct is that many of their coefficients for commonly used resampling ratios are scaled powers of two, which makes way for fast computation. In [Turkowskl 88a], these two functions have been examined for use in image resampling. -2 -1 0 1 2 Figure 5.17: Two Gaussian functions. 5.4 INTERPOLATION KERNELS 145 5.4.7. Exponential Filters A superior class of reconstruction filters can be derived using exponential functions. Consider, for instance, the hyberbolic tangent function tanh defined in Eq. (5.4.34). e x -- e-x tanh (x) = -- (5.4.34) e x + e -x This function has several desirable properties. First, it converges quickly to + 1. Second, its transition from -1 to 1 is sharp. We can sharpen the transition even further by scaling the domain, i.e., use tanh(kx) for k Z 1. In addition, this function is infinitely differenti- able everywhere, i.e., it satisfies an important smoothness constraint. These properties are readily apparent in Fig. 5.18, which illustrates tanh (kx) for k = 1, 4, and 10. Notice that the function quickly approximates Rect for larger values of k. -2 -1 0 1 2 Figure 5.18: Scaled hyperbolic tangent function. Given tanh (kx) as our starting point, we can define a new function that resembles the ideal low-pass filter Rect(f), i.e., a box in the frequency domain. This is done by treating tanh as one half of Rect, and then merely compositing that with a mirror image of itself. Since tanh lies between -1 and 1, some care must be taken to normalize the expression so that it yields a box of unity height. The resulting function is given as Hk(f)=[tanh(k(5+fc))+l'l[tanh(k(-+fc))+l- 1 (5.4.35) where fc is the cut-off frequency. In our examples, we shall use f½ = .5 to conform to the Nyquist rate. The purpose of the addition and division operations is to normalize H,(f) so that 0 < Hk(f) < 1. Function H,t is treated as the desired spectrum of our reconstruction filter. By vary- ing k, we can control the shape of the spectram. For low values of k, H is smooth and resembles a Gaussian function. As k is made larger, Hk will have increasingly sharper comers, eventually approximating a Rect function. Figure 5.19 shows Hk(f ) for k = 1, 4, and 10. Having established H, to be the desired spectrum of our interpolation kemel, the actual kemel is derived by computing the inverse Fourier transform of Eq. (5.4.35). This 146 IMAGE nESAMPLING -2 -1 0 1 2 Figure 5.19: Spectrum HkO e ) is a function of tanh (kx). gives us hk(x), as shown in Fig. 5.20. Not surprisingly, it has infinite extent. However, unlike the sine function that decays at a rate of l/x, hk decays exponentially fast. This is readily verified by inspecting the log plots in Fig. 5.20? This means that we may truncate it with negligible penalty in reconstruction quality. The truncation is, in effect, implicit in the decay of the filter. In practice, a 7-point kernel (with 3 points on each side of the center) yields excellent results [Massalin 90]. n4(x) -4 -3 -2 -1 0 1 2 3 4 (a) h4(x) 0.1 0.01 0.001 [ lc-04 J le-05 I tL -4 -3 -2 -1 0 1 2 3 4 %) h10(x) h10(x) o.oot I le-05 re-06 -4-34-101234 -4-34-101234 (C)  Figure 5.20: Interpolation kernels derived from H4(f) and Ht0Oe). t Note that a linear fall-off in log scale conesponds to an exponential function. COMPARISON OF INTERPOLATION METHODS 147 5.5. COMPARISON OF INTERPOLATION METHODS The quality of the popular interpolation kernels are ranked in ascending order as fol- lows: nearest neighbor, linear, cubic convolution, cubic spline, and sine function. These interpolation metheds are compared in many sources, including [Andrews 76, Parker 83, Maeland 88, Ward 89]. Below we give some examples of these techniques for the magnification of the Star and Madonna images. The Star image helps show the response Directory: filedownloadDownload 2.54 Mb.Share with your friends:

The database is protected by copyright ©ininet.org 2020
send message

Main page