1.

Introduction

The image resampling is the process of geometrically transforming a digital image into a new image. It is an important field in computer graphics and has many applications with the satellite image processing and the medical image processing, in which the images need to be projected or co-registered.¹ The challenges of the image resampling are twofold. The first is how to preserve the information in the input images without introducing artifacts such as aliasing, blurring, or ringing. The second is how to make the algorithm more efficient, and this has been particularly important for the satellite image processing due to strict time requirements. The image resampling is particularly important in the satellite image processing, in which the images in different visible and infrared (IR) channels are required to be co-registered. The image resampling enables the image in one channel to be aligned with another channel, so that the pixels at the same coordinate in the images in different channels correspond to the same geolocation on Earth. This is crucial to the weather and environmental products that require the combinations of image data in different channels.

In this article, we present the fast Fourier transformation resampling (FFTR) algorithm for the image resampling. Instead of the traditional localized approach, the new algorithm represents an image as a global continuous function in the Fourier expansion form, and the Fourier spectrum is obtained through the Fourier transformation of the original image. The image resampling is achieved by introducing a noninteger phase shift in the Fourier expansion. This approach can be easily implemented through the fast Fourier transformation algorithm. The FFTR algorithm is accurate, since both the resampled and the original images share the same Fourier spectrum so that no information is lost during the resampling process. The main characteristics of the resampled image, such as its histogram, remain mostly the same as that of its original image. The FFTR is also generally reversible, and the reversibility has been one of the criteria for the performance of resampling algorithms. The Fourier spectrum obtained from the original image can also be obtained from the resampled image as well, and the resampling can then be performed on the resampled image to obtain the original image. The FFTR is efficient enough for the real-time operations in satellite image processing. The simulation shows that the performance of the FFTR algorithm meets the latency requirements of the Geostationary Operational Environmental Satellite (GOES) instrument data processing.

The FFTR algorithm generally works very well for most of the satellite images. However, the FFTR algorithm does generate artifacts for the images with the hot spots. The hot spots are the pixels on an image whose second-order derivatives are an order of magnitude larger than the average value. The hot spots could be the fire pixels in the IR channel images, of which the temperature could be tens of degrees hotter than the nearby pixels. The hot spots are present in the GOES images in the 3.9-μm IR channel because of the channel’s increased sensitivity to temperature. The hot spots can be regarded as the discontinuous points in the data so that the Fourier expansion alone for the global function is no-longer sufficient for the resampling purpose. In addition to the Fourier expansion, a local Gaussian function is introduced to model the contribution from a hot spot to the global function, so that the contributions from the hot spots to the Fourier expansion would be continuous. Section 5 discusses the procedure on how to estimate the contributions from Gaussian function and Fourier expansion to the hot spot pixels. Once the global function is determined from the existing pixels, the resampling operation becomes straightforward.

The FFT in the image processing is not new. There have been extensive discussions in the literature on the application of FFT in image resampling.¹ However, the FFTR algorithm presented here is significantly different from the existing approaches. Unlike the existing algorithm for the image reconstruction that uses the complex function, the FFTR here uses the real sine or cosine functions in the Fourier expansion, so that the resampling output is real as well. Furthermore, the input data array for FFT with the size requirement of $2^n$ for the integer $n$ is constructed differently. The discrete Fourier transformation (DFT) with a fractional phase shift in the coordinate space becomes the so called $z$-transform, which has been discussed in the literature.² However, the application of FFT algorithm in the image resampling here focuses on the real sine and cosine functions, whereas the application of the FFT in the $z$-transform³ in the signal processing uses the complex functions.

This article is organized as follows. Section 2 provides the detailed discussions on the FFTR algorithm. Section 3 presents the results of the FFTR algorithm and the comparison with the cubic spline (CS) approach, which is commonly used in the interpolation. Section 4 discusses the reversibility of the FFTR algorithm. The image resampling with the hot spots is presented in Sec. 5. Section 6 discusses the application of the FFTR algorithm in the satellite image processing, in particular, the IR channel images of GOES Imager. Finally, the summary is provided in Sec. 7.

2.

FFTR Algorithm

An image can be characterized by a two-dimensional array with the element $P_{i,j}$ for the pixel value at the location $\{i,j\}$. The image resampling is to find the array with the element $P_{i+{\delta }_i,j+{\delta }_j}$ at the location $\{i+{\delta }_i,j+{\delta }_j\}$ based on the input array $P_{i,j}$, where both $i$ and $j$ are integers and ${\delta }_i$ and ${\delta }_j$ are generally fractional numbers. Generally, the resampling along ${\delta }_i$ and ${\delta }_j$ directions can be performed separately. Therefore, we will focus on our discussion to the resampling in one dimension, which can be easily extended to the resampling in two dimensions.

For one-dimensional resampling, one can define a global continuous function $G(x)$ with the condition

Equations (6) and (7) are the Fourier transforms for an input array $P_i$. The resampled array $P_i^{\delta }$ can be obtained by adding the phase shift $\delta$ to the Fourier expansion

The resampling procedure using Eqs. (7) and (8) is straightforward. One can obtain the Fourier amplitude $g(k)$ using Eq. (7) first and then calculating the resampled array $P_i^{\delta }$ using Eq. (8).

The efficiency of the resampling using Eqs. (7) and (8) is $O(N^2)$ for a one-dimensional array, which should be sufficient if the size of an image is small. For images with larger sizes, i.e., generally the case for satellite images, the procedure could not meet the latency requirements for real-time systems. The efficiency of the resampling procedure can be improved to $O(N{\log }_{2}N)$ by using the FFT algorithm.⁴ The FFT algorithm has the size requirement $M=2^m$ for the integer $m$. The global function $G(x)$ in Eq. (5) for FFT becomes

Therefore, one can use Eq. (4) to obtain the amplitude $g(k)$. The resampled array $P_i$ can be obtained by adding the phase shift $\delta$ to the Fourier expansion in Eq. (3):

For an integer $\delta$, Eq. (11) becomes

The elements at $i$ of the resampled array are related to the element at $i+\delta$ of the input array, which should be generally true for any resampling algorithm.

For a noninteger $\delta$, Eq. (11) cannot be used in FFT directly because of the phase shift $\delta$, and direct calculation of Eq. (11) has the efficiency of $O(N^2)$ for a one-dimensional array. In order to use FFT, Eq. (11) is rewritten as

This global function with a phase shift can be evaluated with two separate FFTs. The steps for FFTR are the following:

The $c$ routines for sine and cosine function FFTs come from Ref. 4. Equation (12) is a form of fractional DFT transformations that are discussed in the literature² as the $z$-transform, and the application of the FFT on the $z$-transform has been investigated.³ It should be pointed out that the $z$-transform is different from the fractional Fourier transform (FRFT), which is a rotation in the frequency-time domain. The application of the FFT on FRFT is discussed in Refs. 56.–7.

Figure 1 shows an example of the global function $G(x)$ obtained from an input array using FFTR. The data are from a GOES Imager image in the 10.7-μm IR channel, which is a typical scan line for GOES images in IR channels. The pixel values for GOES images are integers and have the range from 0 to 1023. The result shows that the global function $G(x)$ reproduces exactly the input array data. This can also be treated as the special case of the resampling with phase shift $\delta =0$. The other special cases are $\delta =1$ and $\delta =-1$, so that the resampled array is related to the input array by $P_i^1=P_{i+1}$ and $P_i^{-1}=P_{i-1}$, respectively, as shown in Eq. (12). The numerical results show that these relations are reproduced by the FFTR algorithm.

Fig. 1

Global function $G(x)$ obtained by FFT. The blue line represents the function $G(x)$, and the red dots are the values of the input array with the size of 2126.

Figure 2 shows that Fourier spectrum $g(k)$ used in Eq. (3), where the $x$-axis is the normalized frequency variable $r=k/M$ with the range from 0 to 1. The Fourier spectrum is dominated by the low-frequency amplitudes, and it converges very quickly to zero as the frequency increases. This is generally true for the satellite images with little noise. The Fourier spectrum is in fact a data profile in the frequency space for a given scan line in the image. Since the resampled array is based on the same Fourier spectrum as the input array, the noises in the resampled array should be in the exact same level as in the original array.

Fig. 2

The Fourier spectrum for the data shown in Fig. 1.

3.

Resampling with Noninteger Shifts and Comparison with Cubic Spline Algorithm

For a noninteger shift $\delta$, the image resampling is basically an interpolation of the equal-spaced input arrays. For an image with little noise, there is an expectation of the smoothness behavior between the two neighboring pixels for an interpolation algorithm. The existing interpolation approaches are generally local, which are based on the neighboring pixels. The most often used one is the CS algorithm. (See Ref. 4 for the discussions and routines of the CS interpolation algorithm.) The CS interpolation assumes the continuous first- and second-order derivatives of the input array, and the second-order derivative coefficients used in the interpolation are determined in a slightly nonlocal manner, which makes it more stable than the simple polynomial interpolations. The CS interpolation is more efficient, which is order of $O(N)$ for a one-dimensional array with the size $N$.

Figure 3 shows the resampled lines from the CS and FFTR (FFT), and the original input data are also shown in Fig. 3. The input data array comes from a full disk image of GOES Imager in the 3.9-μm IR channel shown in Fig. 4. The image covers the western hemisphere with 5208 pixel width and 2704 pixel height. A particular interesting local area in Fig. 4 is shown in Fig. 5 that covers Mexico, which contains the dark pixels corresponding to the hot spots. The resampled line in Fig. 3 contains one of the hot pixels shown in Fig. 5. The constant data on the left and right sides are the space data, and the data in the middle are the scan across the Earth along the east-west direction. The data size of the input data array is 5208, which is a typical size for the full disk Earth image for the IR channels. The value of the shift $\delta$ is 0.5, which represents the interpolation point with the largest uncertainty.

Fig. 3

The resampled data from cubic spline (CS) and FTTR algorithms. The original data (Raw) is also shown.

Fig. 4

A full disk image from GOES Imager in the 3.9-μm infrared (IR) channel.

Fig. 5

The enlarged local area in Fig. 4 shows the dark pixels corresponding to the hot spots.

The data shown in Figs. 3–5 are typical for satellite images. The data on the space part are generally constant with small variations with averages of one or two counts. The data on the Earth segment generally change significantly from one place to another with many local minimums and maximums due to water, land, and clouds. For some hot spots, the difference between the values of the two neighboring pixels could be as large as more than 100 counts, while the average difference between the neighboring pixels for the whole image is generally less than 10 counts. Mathematically, both first- and the second-order derivatives are very large on the sharper edges in the Earth segment, and the contributions from the higher order derivatives could also be very significant. This suggests that the CS approach assuming the continuous second-order derivative is no-longer sufficient for the interpolations of the satellite images.

It is very difficult to see the visual effects of the image resampling algorithm on the resampled images, as the changes of the pixel counts are only a few percent of their original counts. To highlight the difference of different resampling algorithms, we show the one-dimensional plot in Fig. 6, which are two different local areas in Fig. 3 where the edge effects are more significant in its neighboring area. These differences are significant, as they leads to the large temperature differences for the IR channel images that could impact the weather forecast. In the region where the changes between the neighboring pixels are small, both FFTR and the CS results agree with each other. The CS interpolation does not produce the smoothness in the area near the sharper edge, whereas FFTR produces more smoothness in the resampling data. The main reason for such difference is that the global function used in FFTR algorithm is continuous in every order derivatives. This shows that the FFTR algorithm is better for the resampling of the satellite images.

Fig. 6

The local area of Fig. 3 with the sharper edge. (a) is the local area with the positive sharper edge and (b) is one of the local areas with the negative sharper edge in Fig. 3.

4.

Reversibility of FFTR Algorithm

Whether a resampling algorithm is reversible is one of the important tests of the accuracy of resampling algorithms. Since the resampled array in FFTR and original data array are represented by the same set of the Fourier spectrum, one can perform the reverse resampling on the resampled array to obtain the original array. For a resampled image with the shift $\delta$, the reverse resampling on the resampled image has the shift $-\delta$. The reversibility of a resampling algorithm should be generally true if the algorithm satisfies Eq. (12) for integer shifts. For example, the reverse resampling on the resampled array with $\delta =1$ has the shift value $\delta =-1$. The resampled array with the integer shift $\delta$ should be reversible for the resampling algorithms that satisfy Eq. (12) for the integer shifts.

For noninteger shift $\delta$, the FFTR output is generally floating numbers, while the pixel counts for satellite images are generally integers. Thus, the resampled array has to be converted from floating numbers of the FFTR outputs to integers, which leads to the loss of precision. The reverse resampling for the images with the integer pixel counts cannot exactly reproduce the original image because of the loss of precision in the conversion process. The numerical evaluation shows that the difference between the two images for a given pixel is 1 pixel count. Figure 7 shows the typical difference between the original array shown in Fig. 3 and the reversed resampled array in the area with the sharper edges.

Fig. 7

The difference between the reverse resampled and original arrays for the data in Fig. 6(a) in the area with the sharp edge.

In practice, there are special cases that the difference between the reverse resampled array and original array could be larger than 1 pixel count. For satellite images, there is an upper limit for the pixel count. For example, the pixel counts for the GOES images have the range from 0 to 1023. If the output of a resampled pixel value is larger than the maximum count, then it is required to change the output value to the maximum count that corresponds to the pixel saturation. If this happens, then the Fourier spectrum is changed accordingly, and the difference between the reverse resampled array and the original array could be larger. The data simulation shows that this is limited to the local area where the saturation occurs.

5.

Resampling Algorithm for the Images with Hot Spots

The hot spots in the satellite images in IR channels could represent the fire pixels in the images, in which the temperature on the hot spot could be tens of degrees higher than that for the neighboring pixels. If the edges on hot spots in satellite images are significantly sharper, then the FFTR resampling also creates the artifacts in the neighboring pixels around the hot spots. Figure 8 shows an example of the artifacts, in which the shift value $\delta$ is 0.5 that represents the maximum resampling uncertainty in the FFTR algorithm. The data in Fig. 8 are the local area of the full Earth scan with the size of 5208, and the hot spot corresponds to the dark pixel near the center in Fig. 6. In this case, the pixel count difference between the neighboring pixels could be as high as 250 counts for the hot spot pixels, while the average difference between the neighboring pixels counts is less than 10 counts. The hot spot pixels in a scan line on a satellite image, in fact, could be regarded as a discontinuity point in the data, and the Fourier expansion does not produce the smooth behavior between the two neighboring pixels around hot spots. The output of the FFTR algorithm shows the sinc function, $\sin (kx)/x$, behavior in the neighboring pixels around the discontinuous point for the noninteger shifts. This behavior is most significant with the shift value $\delta =0.5$. The sinc function behavior for discontinuity points is generally true for the FFTR algorithm, and the discontinuity points may correspond to the hot pixels in the satellite images or the images with sharper edges that need to be treated differently.

Fig. 8

The local area with the hot spots in a data array with the size of 5208. The resampled output is in blue with the shift value 0.5, and the original data are in red.

The discontinuity points, such as the hot spots, can be determined by the second-order derivative. The second-order derivative for an input data array at the element $i$ is defined as

If $D_i^2$ value is larger than a threshold value for a pixel, then it becomes discontinuous that the simple FFTR algorithm will generate the sinc behavior in the neighboring pixels. The $D_i^2$ value is used in the algorithm for the detection of the discontinuity points. The actual threshold value depends on the average $D_i^2$ value of an image.

To remove the artifacts related to the hot spots, the hot spot pixels have to be treated separately from the rest. The global function is modified with the addition of a local function

The $s$ and $e$ in Eq. (18) are the start and end indices for the hot spot pixels in a one-dimension data array. The variable $m$ in Eq. (18) is the scale factor determined by the number of pixels within the hot spot boundary, which is expressed as

Then Eq. (11) becomes

The local function $h(x)$ has a Gaussian form, which can be easily estimated analytically. Thus, the global function within the boundary of the hot pixels comes from both local function $h(x)$ and the Fourier expansion. To separate the contributions between the Fourier expansion and the Gaussian function, the contributions to the Fourier expansion are defined as the linear interpolation based on the pixel count outside the hot spot pixel boundary, which can be written as

Therefore the local function $h(x)$ within the hot spot boundary can be written as

The steps for the image resampling with the local function $h(x)$ are

How these three parameters are estimated depends on how many pixels are there within the hot spot boundary. Generally, the number of hot spot pixels in a satellite image has the range from 1 to 4. The general expressions for these three parameters are

The steps for the FFTR with the local function are the following:

Although the algorithm here is relatively more computational intensive, it is still an order of $O(N\log N)$. In practice, there are only few hot spot pixels in a satellite image; the overhead to the image resampling is generally very small for the entire image.

Figure 9 shows the FFTR with the local function results of same local area as that in Fig. 6. The number of the pixels within the hot spot boundary is 2. The parameter $\beta$ is found to be 270.43, which is approximately the difference between the maximum of the hot spot and the average value of the neighboring pixels. The threshold value used to determine the discontinuity point is 150 pixel counts. The fixed value used to determine the start and end points is 50. The result in Fig. 9 shows a smooth behavior between the neighboring pixels, and the artifact shown in Fig. 8 has been removed. The simulation on the GOES images in 3.9-μm IR channel shows that the introduction of the Gaussian function is working very well in resampling of the images with the hot spots.

Fig. 9

The same local area with the hot spot. The resampled output (in blue) comes from the modified FFTR algorithm.

The threshold values in determining the discontinuity points and start/end points depend on the characteristics of images. For GOES images, the hot spot pixels are concentrated in the 3.9-μm IR channel images. Although the image contents change, the average values of the first- and second-order derivatives are nearly constant. Therefore, the fixed threshold value in determining the discontinuity points and start/end points for the hot spot pixels can be used.

It should also be pointed out that not all discontinuity points are hot pixels with the Gaussian-like distribution. There are also discontinuity points attributed to the sharp edges due to the saturations where the pixel counts reach to upper limit, which also causes the oscillation behavior near the discontinuity points. The saturation of the image pixels is mostly due to the contamination from sun light during the certain time each day, which could lead to the sudden drop of the pixels counts from the constant saturation points to the nearby pixels. The approach to the discontinuity points is the same, using Eq. (20) for the Fourier expansion, so that it is continuous and a local function for the sharp edge. One can certainly use the same approach for the hot pixels by defining a local function for the sharp edge around the saturation points. This is left for the future investigation.

The introduction of the local function does have an impact on the reversibility of the resampling algorithm, particularly around the start and end points, since the start and end points are not the same during the reverse resampling and resampling processes. The difference between the reverse resampled array and the original array at the start and end points could be as high as 4 pixel counts. However, the larger differences are only limited to the pixels at the start and end points, which is localized.

6.

Resampling on GOES Images

The resampling in the satellite image processing is mostly used in correcting the co-registration errors in the IR channels. The Imager for the current generation of GOES has wavelengths centered at 3.9, 13.3, 10.7, and 6.5 μm, which are identified as channels 2, 6, 4, and 3, respectively.⁸ The larger co-registration errors among the IR channels certainly affect the reliability of the products that rely on the combinations of the different IR channels. The samples of these products are the cloud mask, fog, and the hot spot detection. Studies⁹ have found that the co-registration errors are indeed more than 1 IR pixel for some GOES during a certain time each day.

There are two components in correcting the co-registration errors in GOES IR channel images.¹⁰ The first component is the GOES Imager IR Channel to Channel Co-Registration Characterization (GII4C) algorithm to determine the co-registration error among the IR channels based on the existing GOES images, and the second one is to correct the co-registration by the image resampling algorithm. The GII4C algorithm uses the statistical correlations between the brightness temperatures for the images between the different channels to determine the co-registration errors. The statistical correlation is a function of the phase shift factor in the resampling algorithm, and the co-registration error between the images of two IR channels is the phase shift that corresponds to the maximum correlation between the brightness temperatures for the images between the two IR channels. Thus, the resampling algorithm also plays the key role in determining the co-registration errors. The co-registration error is used as the input for the resampling algorithm in the GOES ground system to correct the error in the images in the IR channels. Because the differences between the original and resampled images are too small to see the visual effects between images, we show a resampled image in the 3.9-μm IR channel for GOES-N (or GOES13) in Fig. 10 which is a typical GOES image that covers the continental United States (CONUS). The image is resampled along the west/east direction, in which the co-registration error is found to be large. Figure 11 shows that the same image in channel 4, which is used as the reference channel, used to calculate the co-registration error by GII4C algorithm. The shift value $\delta$ in this case is $-1.1706$.

Fig. 10

Channel 2 Image for GOES 13, the image covers the continental United States (CONUS).

Fig. 11

Channel 4 Image for GOES 13 is used as the reference image to calculate the co-registration errors between the images in channels 2 and 4.

Figure 12 shows the histograms for the resampled image shown in Fig. 10 and its original image. The histogram is a statistical profile of an image, and the results in Fig. 10 shows that the statistical profile remains mostly the same as that for the original image. The main difference between the resampled and the original images is that the histogram for the resampled image has less fluctuation, which is smoother than the histogram for the original image. The other characteristics of the satellite images, such as the mean value of the image and the gradient histogram, have also not been changed in the resampled image.

Fig. 12

The histograms of the resampled and the original images.

7.

Summary

The FFTR algorithm uses a global continuous function in the Fourier expansion form to represent an image. The Fourier amplitudes are obtained from the Fourier transformation of the original image, and the image resampling is achieved by introducing a phase shift in the Fourier expansion. This algorithm is accurate and reversible, as the original and the resampled images belong to the same global function based on the same set of Fourier spectrum. Our study finds that the local interpolation approach is not as accurate as the FFTR algorithm for the satellite image resampling.

The hot spots in the satellite images are treated as the local Gaussian function, so that the global function for an image with hot spots consists of the Fourier expansion and local Gaussian function. By including the Gaussian function for the hot spot, the contribution of the hot spot pixels to the Fourier becomes more continuous, so that the artifacts associated with the Fourier expansion alone are removed. This modified FFTR algorithm is mostly useful for the resampling of the 3.9-μm IR channel images, in which the hot spot pixels are present. For other longer wave IR channels, the original FFTR algorithm would be sufficient. For other types of images, the Gaussian function alone is no-longer sufficient to model the local discontinuity points. For example, the discontinuity points could be a step function, which has to be treated differently. Thus, FFTR could not be used directly on the images with the simple objects, in which the discontinuity points are the main characteristics of the images.

The histograms are nearly identical between the resampled and the original images. The further evaluations also show that the mean value and mean gradient value are nearly identical as well. These results show that the FFTR algorithm maintains the image quality in the resampling process.

The FFTR algorithm is used in the GOES Imager channel to channel co-registration correction program, in which the co-registration errors are extracted and corrected by FFTR. The co-registration error correction images improve the reliability of the weather products that rely on the combinations of different IR channels significantly.¹⁰ It is expected that this program to be implemented and be operational in the near future.