Interpolation in Fourier Space

⇐ ПредыдущаяСтр 52 из 85Следующая ⇒

Interpolation reduces to a simple operation in the Fourier domain. As shown by Eq. (10.45), the transfer function of an ideal interpolation ker- nel is a rectangular (box) function which is zero outside the wave num- bers that can be represented. This basic fact suggests the following interpolation procedure in Fourier space:

Enlarge the matrix of the Fourier transformed image. If an M M matrix is increased to an M' M' matrix, the image in the spatial domain is also increased to an M' M' image. Because of the reci- procity of the Fourier transform, the image size remains unchanged. Only the spacing between pixels is decreased, resulting in a higher spatial resolution:

M∆ k → M'∆ k ◦ • ∆ x = 1 → ∆ x' = 1

(10.48)

M∆ k M'∆ k

2. Fill the padded area in the Fourier space with zeroes and compute an inverse Fourier transform.

Theoretically, this procedure results in a perfectly interpolated im- age. Unfortunately, it has three serious drawbacks.

1. The Fourier transform of a ﬁ nite image implies a cyclic repetition of the image both in the spatial and Fourier domain. Thus, the convo- lution performed by the Fourier transform is also cyclic. This means that at the right or left edge of the image, convolution continues with the image at the opposite side. As the real world is not periodic and interpolation masks are large, this may lead to signiﬁ cant distortions of the interpolation even at quite large distances from the edges of the image.

The Fourier transform can be computed eﬃ ciently only for a spec- iﬁ ed number of values for M'. Best known are the fast radix-2 al- gorithms that can be applied only to images of the size M' 2N' (Section 2.5.2). Therefore, the Fourier transform-based interpolation is limited to scaling factors of powers of two.

3. As the Fourier transform is a global transform, it can be applied only to scaling. In principle, it could also be applied to rotation and aﬃ ne

272 10 Pixel Processing

1/2(g1/2+g-1/2)

g-1/2

g1/2

linear interpolation kernel

g3/2

g-2/3

g-1/2

g1/2

g3/2

-3/2

-1/2

1/2

3/2

-3/2

-1/2

1/2

3/2

Figure 10.17: Illustration of linear interpolation: a at x = 0, the mean of g_1/2

and g₋_1/2 is taken, b at x = 1/2, g_1/2 is replicated.

transforms. But then the interpolation problem is only shifted from the spatial domain to the wave number domain.

Linear Interpolation

Linear interpolation is the classic approach to interpolation. The inter- polated points lie on pieces of straight lines connecting neighboring grid points. In order to simplify the expression, we use in the following nor- malized spatial coordinates x˜ = x/∆ x. We locate the two grid points at

− 1/2 and 1/2. This yields the interpolation equation

g(x˜ ) = g 1 / 2 + g − 1 / 2 +

.g1/2 − g− 1/2Σ x˜

for |x˜ |≤ 1/2. (10.49)

By comparison of Eq. (10.49) with Eq. (10.43), we can conclude that the continuous interpolation mask for linear interpolation is

0 otherwise

h (x˜ ) = . 1 − |x˜ | |x˜ |≤ 1

. (10.50)

Its interpolatory nature is illustrated in Fig. 10.17. The transfer function of the interpolation mask for linear interpolation, the triangle function h₁(x) Eq. (10.50), is the squared sinc function (± R5)

hˆ 1

(k˜ ) =

sin² π k˜ /2. (π k˜ /2)2

(10.51)

A comparison with the ideal transfer function for interpolation Eq. (10.45) shows that two distortions are introduced by linear interpolation:

While low wave numbers (and especially the mean value k˜ 0) are interpolated correctly, high wave numbers are slightly reduced in am- plitude, resulting in some degree of smoothing. At k˜ = 1, the transfer function is reduced to about 40 %: hˆ 1(1) = (2/π )2 ≈ 0.4.

10.6 Interpolation† 273

2. Since hˆ 1(k˜ ) is not zero at wave numbers k˜ > 1, some spurious high wave numbers are introduced. If the continuously interpolated image is resampled, this yields moderate aliasing. The ﬁ rst sidelobe has an amplitude of (2/3π )2 ≈ 0.045.

| |≥ =

= − −

If we interpolate only the intermediate grid points at x˜ 0, the con- tinuous interpolation function Eq. (10.50) reduces to a discrete convolu- tion mask with values at x˜ [... 3/2 1/2 1/2 3/2...]. As Eq. (10.50) is zero for x˜ 1, we obtain the simple interpolation mask H 1/2[11] with the transfer function

Hˆ 1(k˜ ) = cos π k˜ /2. (10.52)

− +

∈ −

The transfer function is real, so no phase shifts occur. The signiﬁ cant amplitude damping at higher wave numbers, however, shows that struc- tures with high wave numbers are not correctly interpolated. Phase shifts do occur at all other values except for the intermediate grid points at x˜ 0. We investigate the phase shift and amplitude attenuation of linear interpolation at arbitrary fractional integer shifts H [ 1/2, 1/2]. The interpolation mask for a point at H is then [1/2 H, 1/2 H]. The mask contains a symmetric part [1/2, 1/2] and an antisymmetric part [− H, H]. Therefore, the transfer function is complex and has the form

hˆ 1(H, k˜ ) = cos π k˜ /2 + 2iH sin π k˜ /2. (10.53)

In order to estimate the error in the phase shift, it is useful to compen- sate for the linear phase shift ∆ ϕ = Hπ k˜ caused by the displacement H. According to the shift theorem (Theorem 3, p. 52, ± R4), it is required to multiply Eq. (10.53) by exp(− iHπ k˜ ): Then we obtain

= = =

hˆ 1(H, k˜ ) = (cos π k˜ /2 + 2iH sin π k˜ /2) exp(− iHπ k˜ ). (10.54) Only for H 0 and H 1/2 is the transfer function real: hˆ 1(0, k˜ )

cos π k˜ /2, h₁(1/2, k˜ ) 1; but at all other fractional shifts, a non-zero

phase shift remains, as illustrated in Fig. 10.18. The phase shift ∆ ϕ is expressed as the position shift ∆ x of the corresponding periodic struc- ture, i. e., ∆ x = ∆ ϕ λ /2π = ∆ ϕ /(π k˜ ).

Polynomial Interpolation

Given the signiﬁ cant limitations of linear interpolation as discussed in Section 10.6.3, we ask whether higher-order interpolation schemes per- form better. The basic principle of linear interpolation was that a straight line was drawn to pass through two neighboring points. In the same way, we can use a polynomial of degree P that must pass through P + 1 points with P + 1 unknown coeﬃ cients a_p:

g_r (x˜ ) = a_px˜ ^p. (10.55)

p=0

274 10 Pixel Processing

0.8

0.6

0.4

0.2

0.95

0.9

0.85

0.8

0.75

0.7

0.1

0.05

-0.05

-0.1

0.06

0.04

0.02

-0.02

-0.04

-0.06

-0.4 -0.2 0 0.2 0.4 e

= =

Figure 10.18: Amplitude attenuation (left column) and phase shift expressed as a position shift ∆ x ∆ ϕ λ /2π (right column) in radians for wave numbers k˜ 1/4, 1/2, 3/4, as indicated, displayed as a function of the fractional position from

−

1/2 to 1/2 for linear interpolation ( a and b ) and cubic B-spline interpolation ( c and d ).

For symmetry reasons, in case of an even number of grid points, we set their position at half-integer values

x˜ _p = 2p − P . (10.56)

+ +

From the interpolation condition at the grid points g_r (x˜ _p) g_p, we obtain a linear equation system with P 1 equations and P 1 unknowns a_P of the following form when P is odd:

 g0

  1 − P /2 P 2/4 − P 3/8 ··· 

 

 ..  ..

a0 

 g   1 − 1/2 1/4 − 1/8 ···  . 

 (P− 1)/2  

  

(10.57)



g_(P_+1)/2  =  1 1/2 1/4 1/8 ···  . 

 .  .

  

  gP

   

1 P /2 P 2/4 P 3/8 ···

   aP 

from which we can determine the coeﬃ cients of the polynomial. For a cubic polynomial (P = 3), the equations system is

10.6 Interpolation† 275

A b

1 1

0.8

0.6

0.4

0.99

1 3

0.98

5 7

0.97

0.2

0.96

0.95 0 0.2 0.4 0.6 0.8 ~ 1

Figure 10.19: Transfer function of discrete polynomial interpolation ﬁ lters to interpolate the value between two grid points. The degree of the polynomial (1 = linear, 3 = cubic, etc.) is marked on the graph. The dashed line marks the transfer function for cubic B-spline interpolation (Section 10.6.5). a The full range, b a 5 % margin below the ideal response hˆ (k˜ ) = 1.

 g0 

 1 − 3/2 9/4 − 27/8   a₀ 



    =

g₁ 1 − 1/2 1/4 − 1/8 a₁ 

(10.58)

g₂ g3 

1 1/2 1/4 1/8 a₂

 

  



 1 3/2 9/4 27/8   a₃ 

with the solution

 a0

  − 3 27 27 − 3   g₀ 



 





a₁ = 2 − 54 54 − 2 g₁

 . (10.59)



 

 a2  48

12 − 12 − 12 12 g₂ 

 a3 

 − 8 24 − 24 8   g₃ 

− −

= = − + + −

From this solution, we can infer, for example, that the point at x˜ 0 is interpolated by g_r (0) a₀ 1/16g₀ 9/16g₁ 9/16g₂ 1/16g₃ corresponding to the interpolation mask 1/16[ 1, 9, 9, 1].

Figure 10.19 shows the transfer functions for a polynomial interpo- lation of various degrees. With increasing degree P of the interpolating polynomial, the transfer function approaches the box function better. However, convergence is slow. For an accurate interpolation, we must take a large interpolation mask.

10.6.5 Spline-Based Interpolation‡

Besides its limited accuracy, polynomial interpolation has another signiﬁ cant disadvantage. The interpolated curve is not continuous at the grid points al- ready in its ﬁ rst derivative. This is due to the fact that for each interval between grid points another polynomial is taken. Thus, only the interpolated function is continuous at the grid points but not the derivatives.

Splines avoid this disadvantage by additional constraints for the continuity of derivatives at the grid points. From the wide classes of splines, we will here dis- cuss only one class, the B-splines. As B-splines are separable, it is suﬃ cient to

276 10 Pixel Processing

0.8

0.6

0.4

0.2

-0.2 -2 -1 0 1 2

0.8

0.6

0.4

0.2

-0.2

-3 -2 -1 0 1 2 3

Figure 10.20: a B-spline interpolation kernels generated by cascaded convolu- tion of the box kernel of order 0 (nearest neighbor), 1 (linear interpolation), 2 (quadratic B-spline), and 3 (cubic B-spline); b corresponding transfer functions.

discuss the properties of 1-DB-splines. From the background of image process- ing, the easiest access to B-splines is their convolution property. The kernel of a P-order B-spline curve is generated by convolving the box function P 1 times with itself (Fig. 10.20a):

β _P (x˜ ) = Π ( x ˜ ) ∗. . . ∗ Π ( x ˜ ). (10.60)

(P+1˛ )¸ times r

The transfer function of the box function is the sinc function ( R5). Therefore, the transfer function of the P-order B-spline is

(π k˜ /2)

˜ P+1

β ˆ P (kˆ ) = . si n π k/ 2 Σ

. (10.61)

Figure 10.20b shows that the B-spline function does not make a suitable inter- polation function. The transfer function decreases too early, indicating that B-spline interpolation performs too much averaging. Moreover, the B-spline kernel does not meet the interpolation condition Eq. (10.46) for P > 1.

B-splines can only be used for interpolation if ﬁ rst the discrete grid points are transformed in such a way that a following convolution with the B-spline kernel restores the original values at the grid points. This transformation is known as the B-spline transformation and constructed from the following condition:

− = =

g_p(x) =. c_nβ _P (x − x_n) with g_p(x_n) = g(x_n). (10.62)

If centered around a grid point, the B-spline interpolation kernel is unequal to zero only for three grid points. The coeﬃ cients β ₃( 1) β ₋₁, β ₃(0) β ₀, and β ₃(1) β ₁ are 1/6, 2/3, and 1/6. The convolution of this kernel with the unknown B-spline transform values c_n should result in the original values g_n at the grid points. Therefore,

g = c ∗ β ₃ or g_n =

n'.=− 1

c_n₊_n' β _n'. (10.63)

10.6 Interpolation† 277

Equation (10.63) constitutes the sparse linear equation system

 g0

 4 1 0 ... 0 1 

...

 1 4 1 0 0 



 = 6 

 



  0 1 4 1 0 ... 

c0 

 1  ..

....   



(10.64)

  gN− 1  

 . . .. . . 1 4 1 0    cN  

− 1

  0 0 1 4 1

... 



 1 0 ... 0 1 4 

using cyclic boundary conditions. The determination of the B-spline transfor- mation thus requires the solution of a linear equation system with N unknowns. The special form of the equation system as a convolution operation, however, allows for a more eﬃ cient solution. In Fourier space, Eq. (10.63) reduces to

g ˆ = β ˆ ˆ c. (10.65)

= +

The transfer function of β ₃ is β ˆ 3(k˜ ) 2/3 1/3 cos(π k˜ ). As this function has no zeroes, we can compute c by inverse ﬁ ltering (Section 4.2.10), i. e., convoluting g with a mask that has the transfer function

β ˆ − 1(k˜ ) = β ˆ T (k˜ ) = 1 . (10.66)

3 2/3 + 1/3 cos π k˜

Such a transfer function is a kind of recursive ﬁ lter (Section 4.2.10) that is ap- plied ﬁ rst in the forward and then in the backward direction with the following

recursion [187]:

√ '

gn'

= g_n − (2 − √ 3)(g_n₋₁ − g_n)

(10.67)

c_n' = g_n' − (2 − 3)(c_n₊₁ − g_n' ).

The whole operation takes only two multiplications and four additions.

The B-spline interpolation is applied after the B-spline transformation. In the continuous case, using Eqs. (10.61) and (10.66), this yields the eﬀ ective transfer function

β ˆ I (k˜ ) =

sin⁴(π k˜ /2)/(π k˜ /2)4

(2/3 + 1/3 cos π k˜ ) (10.68)

Essentially, the B-spline transformation performs an ampliﬁ cation of high wave numbers (at k˜ 1 by a factor 3) which compensates the smoothing of the B- spline interpolation to a large extent.

We investigate this compensation at the grid points and at the intermediate points. From the equation of the cubic B-spline interpolating kernel Eq. (10.60) (see also Fig. 10.20a) the interpolation coeﬃ cients for the grid points and inter- mediate grid points are

1/6 [14 1] and

1/48 [1 23 23 1],

(10.69)

278 10 Pixel Processing

respectively. Therefore, the transfer functions are

2/3 + 1/3 cos π k˜ and

23/24 cos(π k˜ /2) + 1/24 cos(3π k˜ /2), (10.70)

respectively. At the grid points, the transfer function exactly compensates — as expected — the application of the B-spline transformation Eq. (10.66). Thus, the interpolation curve goes through the values at the grid points. At the interme- diate points the eﬀ ective transfer function for the cubic B-spline interpolation

is then

β ˆ (1/2, k˜ ) = 23/24 cos(π k˜ /2) + 1/24 cos(3π k˜ /2). (10.71)

I 2/3 + 1/3 cos π k˜

−

The amplitude attenuation and the phase shifts expressed as a position shift in pixel distances are shown in Fig. 10.18c, d. Note that the shift is related to the intermediate grid. The shift and amplitude damping is zero at the grid points [ 0.5, 0.5]^T. While the amplitude damping is maximal for the interme-

diate point, the position shift is also zero at the intermediate point because of symmetry reasons. Also, at the wave number k˜ 3/4 the phase shift is unfor- tunately only about two times smaller than for linear interpolation (Fig. 10.18b).

It is still signiﬁ cant with a maximum of about 0.13.

This value is much too high for algorithms that ought to be accurate in the 1/100 pixel range. If no better interpolation technique can be applied, this means that the maximum wave number should be lower than 0.5. Then, the maximum shift is lower than 0.01 and the amplitude damping less than 3 %.

Note that these comments on phase shifts only apply for arbitrary fractional shifts. For pixels on the intermediate grid, no position shift occurs at all. In this special case — which often occurs in image processing, for instance for pyramid computations (Chapter 5) — optimization of interpolation ﬁ lters is quite easy because only the amplitude damping must be minimized over the wave number range of interest.

10.6.6 Optimized Interpolation‡

Filter design for interpolation — like any ﬁ lter design problem — can be treated in a mathematically more rigorous way as an optimization problem. The general idea is to vary the ﬁ lter coeﬃ cients in such a way that the derivation from the ideal transfer function reaches a minimum. For non-recursive ﬁ lters, the transfer function is linear in the coﬃ cients h_r:

hˆ (k˜ ) = h_r fˆ r (k˜ ). (10.72)

r=1

Let the ideal transfer function be hˆ I (k˜ ). Then the optimization procedure should minimize the integral

∫ . R  .

w(k˜ ). . h_r fˆ r (k˜ ) − hˆ I (k˜ ).

dk˜. (10.73)

. r=1 .

In this expression, a weighting function w(k˜ ) has been introduced which allows control over the optimization for a certain wave number range. In equation

10.6 Interpolation† 279

0.8

0.6

0.4

0.2

0 0.2 0.4 0.6 0.8 1

1.01

1.005

0.995

0.99

Figure 10.21: Transfer function of interpolation kernels optimized with the weighted least squares technique of Eq. (10.76) and Eq. (10.77) with R = 3 (solid line) and of Eq. (10.78) for R = 2 (dashed line). The weighting function used for the optimization is shown in a as a thin solid line; b shows a narrow sector of the plot in a for a better estimation of small deviations from ideal values.

Eq. (10.73) an arbitrary L_n-norm is included. Mostly the L₂-norm is taken, which minimizes the sum of squares.

For the L₂-norm, the minimization problem results in a linear equation system for the R coeﬃ cients of the ﬁ lter which can readily be solved:

Mh = d (10.74)



with

 hI fˆ 1 

ˆ 2 fˆ 1fˆ 2 ···

fˆ 1fˆ R 

d = 

h_Ifˆ 2

. 

and M = 

fˆ 1fˆ 2

ˆ 2 ···

fˆ 2fˆ R



. 

   





. .



h_Ifˆ R

where the abbreviation

fˆ 1fˆ R

fˆ 2fˆ R ··· ˆ 2

eˆ (k˜ ) = ∫

w(k˜ ) · eˆ (k˜ )dk˜

(10.75)

for an arbitrary function e(k˜ ) has been used.

The ﬂ exibility of the least squares optimization technique for ﬁ lter design is given by the free choice of the weighting function w(k˜ ) and the careful consid- eration of the symmetry properties and other features of a ﬁ lter by the choice of the transfer function in Eq. (10.72). For illustration, we discuss the following

two approaches:

hˆ (k˜ ) =. h_r cos. 2 r − 1 π k˜ Σ (10.76)

and

r=1 2

hˆ (k˜ ) = cos. 1 π k˜ Σ +. h_r  cos. 2 r − 3 π k˜ Σ − cos. 1 π k˜ Σ  . (10.77)

2 r=2 2 2

Both ﬁ lters result in a symmetric mask by the choice of the cosine function. Equation Eq. (10.77) ensures that hˆ (0) 1, i. e., the mean gray values are pre- served by the interpolation. This is done by forcing the ﬁ rst coeﬃ cient, h₁, to

280 10 Pixel Processing

be one minus the sum of all others. Equation Eq. (10.76) does not apply this constraint. Figure 10.21 compares the optimal transfer functions with both approaches for R = 3.

The ﬁ lters are signiﬁ cantly better than those obtained by polynomial and cu- bic B-spline interpolation (Fig. 10.19). The additional degree of freedom for Eq. (10.76) leads to signiﬁ cantly better solutions for the wave number range where the weighting function is maximal.

Even better interpolation masks can be obtained by using a combination of non- recursive and recursive ﬁ lters, as with the cubic B-spline interpolation:

cos.1/2 π k˜ Σ +. h_r Σ cos.(2r − 3)/2 π k˜ Σ − cos.1/2 π k˜ Σ Σ

hˆ (k˜ ) =

r=2

1 − α + α cos.π k˜ Σ

. (10.78)

With recursive ﬁ lters, the least squares optimization becomes nonlinear be- cause hˆ (k˜ ) in Eq. (10.78) is nonlinear in the parameter α of the recursive ﬁ lter. Then, iterative techniques are required to solve the optimization problem. Fig- ure 10.21c, d shows the transfer functions for R = 2. A more detailed discussion of interpolation ﬁ lters including tables with optimized ﬁ lters can be found in Jä hne [81].

10.6.7 Fast Algorithms for Geometric Transforms‡

With the extensive discussion on interpolation we are well equipped to devise fast algorithms for the diﬀ erent geometric transforms. Basically, all fast in- terpolation algorithms use the following two principles: eﬃ cient computation employing the interpolation coeﬃ cients and partition into 1-Dgeometric trans- forms.

First, many computations are required to compute the interpolation coeﬃ cients for fractional shifts. For each shift, diﬀ erent interpolation coeﬃ cients are re- quired. Thus we must devise the transforms in such a way that we need only constant shifts for a certain pass of the transform. If this is not possible, it might still be eﬃ cient to precompute the interpolation coeﬃ cients for diﬀ erent fractional shifts and to save them for later usage.

Second, we learnt in Section 10.6.1 that interpolation is a separable procedure. Taking advantage of this basic fact considerably reduces the number of opera- tions. In most cases it is possible to separate the two- and higher dimensional geometric transforms into a series of 1-Dtransforms.

10.7 Further Readings‡

Holst [70, 72] and Biberman [7] deal with radiometric calibration of sensors and cameras in the visible and infrared. A detailed discussion of interpola- tion ﬁ lters including tables with ﬁ lter coeﬃ cients and eﬃ cient algorithms for geometric transforms can be found in Jä hne [81, Chap. 8]. Readers interested in the mathematical background of interpolation are referred to Davis [25] and Lancaster and Salkauskas [103]. Wolberg [200] expounds geometric transforms.

Part III

Feature Extraction

Averaging

Introduction

In this chapter we will discuss neighborhood operations for performing the elementary task of averaging. This operation is of central importance for low-level image processing. It is one of the building blocks for the more complex feature extraction operators discussed in Chapters 13–15. In the simplest case, objects are identiﬁ ed as regions of constant radiance, i. e., gray values. Then, averaging gives adequate mean values of the gray values within the object. This approach, of course, implies a simple model of the image content. The objects of interest must indeed be characterized by constant gray values that are clearly diﬀ erent from

the background and/or other objects.

However, this assumption is seldom met in real-world applications. The intensities will generally show some variations. These variations may be an inherent feature of the object or could be caused by the image formation process. Typical cases are noise, a non-uniform illumination, or inhomogeneous background.

In complex cases, it is not possible to distinguish objects from the background with just one feature. Then it may be a valid approach to compute more than one feature image from one and the same image. This results in a multicomponent or vectorial feature image.

The same situation arises when more than one image is taken from a scene as with color images or any type of multispectral image. There- fore, the task of averaging must also be applied to vectorial images. In image sequences, averaging is extended into the time coordinate to a spatiotemporal averaging.

⇐ Предыдущая 47 48 49 50 515253 54 55 56 Следующая ⇒

Последнее изменение этой страницы: 2019-05-04; Просмотров: 217; Нарушение авторского права страницы