IDW interpolation method suitable for grid image resolution enhancement

Abstract

The application discloses an IDW interpolation method suitable for improving resolution of grid images and belongs to the technical field of remote sensing image processing. The method comprises the following steps: firstly, configuring interpolation parameters for inputted original grid images, and extracting rational coordinates and pixel width; then, calculating weights based on a distance matrix from target pixels to reference pixels, normalizing the weights, and reorganizing the weights to generate reusable three-dimensional weight templates; traversing the original two-dimensional numerical matrix by using a set window, and adopting boundary padding processing for the edges of the matrix; using broadcast multiplication to batch operate the three-dimensional weight templates and the window matrix, calculating target pixel values, splicing the target pixels obtained by interpolation of each window to generate high-resolution grid images. Block noise and bull's eye effect can be eliminated by using Gaussian filtering, spatial continuity is improved, geographical reference is restored, and high-resolution TIFF images are outputted. The application can be applied to grid image processing in the fields of meteorology, ecological environment and the like.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, and in particular to an IDW interpolation method suitable for improving the resolution of raster images. Background Technology

[0002] Improving image resolution is an important technique in remote sensing image processing. Currently, the mainstream methods can be divided into traditional interpolation methods, pixel- or neighborhood-based improved interpolation methods, and deep learning methods.

[0003] Traditional interpolation methods include nearest neighbor interpolation, bilinear interpolation, and bicubic interpolation. These methods are simple in principle, computationally efficient, and suitable for rapid resampling in low-computation scenarios. The second category of methods improves upon traditional interpolation methods by introducing distance weights, neighborhood features, and terrain constraints. These include inverse distance weighted interpolation, kriging interpolation, and spline interpolation. These methods retain the advantages of low computational cost and ease of implementation, while significantly improving accuracy compared to the first category. Deep learning methods are primarily applied in super-resolution reconstruction. By learning the spatial features and texture patterns of raster images through neural networks, they can further improve accuracy. However, this requires higher-performance computing power and a large training dataset.

[0004] Inverse distance weighted interpolation (IDW), as a compromise between traditional interpolation methods and deep learning, is characterized by its simplicity, ease of implementation, strong controllability, and high accuracy, making it a superior choice for downscaling single-band raster images when deep learning frameworks are unavailable. Current quantitative remote sensing products and reanalysis meteorological data such as ERA5 generally exhibit characteristics of long-term series, high temporal resolution, and wide spatial coverage. Performing batch, long-term spatial downscaling on such data involves a massive computational load and demands high processing efficiency. However, current mainstream applications such as ArcGIS's built-in IDW interpolation cannot directly process raster images; it requires converting raster data into point data before interpolation, a cumbersome process with low computational efficiency and strong platform dependence, failing to meet the practical needs of efficient downscaling of long-term raster data. Therefore, this algorithm aims to improve the IDW interpolation method by enhancing computational efficiency. Summary of the Invention

[0005] The purpose of this invention is to address the problems of redundant steps, low computational efficiency, and strong platform dependence in traditional IDW interpolation methods when processing raster data, and to propose an IDW interpolation method suitable for improving the resolution of raster images, so as to achieve direct, efficient, and high-precision downscaling of raster images.

[0006] The technical solution adopted in this invention is: an IDW interpolation method suitable for improving the resolution of raster images, which includes the following steps:

[0007] Step 1, Data Preprocessing

[0008] Input the raw raster image data and configure the difference parameters, including: sliding window size, image magnification, and IDW power parameter;

[0009] The original raster image data is converted into a floating-point two-dimensional numerical matrix (i.e., the original two-dimensional numerical matrix), and the geographic coordinates and pixel width are extracted for subsequent accurate calculation of the new pixel size and spatial location of the high-resolution image. This ensures that the output image has correct georeferenced information and avoids spatial offset and resolution distortion.

[0010] Step 2, pre-calculation of weight template

[0011] The coordinate matrix of the reference pixel is generated according to the sliding window size, and then each reference pixel is subdivided according to the image magnification to obtain multiple target pixels, and the coordinate matrix of the target pixels is generated.

[0012] Calculate the distance matrix from the target pixel to the reference pixel, set a minimum distance threshold δ to replace the minimum distance in the distance matrix, thereby avoiding a denominator of 0; that is, for each distance element in the distance matrix, if its value is less than δ, then set it to δ; then calculate the weights according to the IDW formula and normalize them, and then reorganize the weights to generate a reusable 3D weight template.

[0013] Step 3, sliding window interpolation

[0014] Traverse the two-dimensional numerical matrix using the sliding window size, and perform boundary padding on the edges of the current window (i.e., the local window matrix) to obtain the padded local window matrix, and then flatten it into a one-dimensional array.

[0015] Broadcast multiplication is used to perform batch operations on the three-dimensional weight template and the one-dimensional array corresponding to the window matrix to calculate the target pixel value. Then, each target pixel is stitched together to generate a high-resolution raster image. Based on the extracted geographic coordinates and pixel width, georeference restoration is performed on the generated high-resolution raster image.

[0016] Furthermore, in step 2, the coordinate matrix of the reference pixel is specifically as follows:

[0017] Using the top-left corner of the first reference cell in the sliding window as the origin, the coordinates of the center point of each reference cell are generated to represent the position of the reference cell;

[0018] Based on the image magnification factor r, the center grid of the current sliding window is subdivided into r×r target pixels, and the center coordinates of each target pixel represent the position of the target pixel.

[0019] Furthermore, the location of the target pixel is:

[0020]

[0021]

[0022] in, , The first subdivision of the mesh Line number The x and y coordinates of the center point of the target cell in the column. Sliding window size.

[0023] Furthermore, in step 2, the calculation and normalization of weights according to the IDW formula are specifically as follows:

[0024] Based on the distance from i reference pixels to the j-th target pixel Calculate initial weights ;in, The preset IDW exponentiation parameter;

[0025] Then adjust the initial weights Normalization is performed to obtain normalized weights. , where n is the number of target pixels.

[0026] Furthermore, in step 2, the reusable 3D weight template is: normalized weights The current sliding window is organized into a three-dimensional structure based on the number of rows, columns, and reference pixels of the center grid of the current sliding window, generating a reusable three-dimensional weight template.

[0027] Furthermore, in step 3, the target pixel value is:

[0028]

[0029] in, Let be the pixel value of the i-th target pixel. The pixel is the j-th reference pixel.

[0030] Furthermore, step 3 also includes: after eliminating the generated high-resolution raster image using Gaussian filtering, georeference restoration is performed to output a high-resolution TIFF image.

[0031] The technical solution provided by this invention brings at least the following beneficial effects:

[0032] (1) Significantly improved computational efficiency: The redundant steps of raster to point conversion are removed, the weight template is calculated once and reused globally, and the broadcast multiplication accelerates the operation. The time complexity is reduced from O(n²) to O(n), and the computation time is shortened by 4 to 100 times compared with ArcGIS built-in IDW.

[0033] (2) No loss of interpolation accuracy: The interpolation results are highly consistent with the spatial distribution of the original image, with a maximum deviation of 0.008%, a minimum deviation of 0.004%, no deviation of the mean, and a standard deviation of 0.54%, which is better than the built-in ArcGIS algorithm.

[0034] (3) Result quality optimization: edge completion processing improves edge pixel accuracy, Gaussian filtering eliminates block effect and bullseye effect, and the spatial continuity of the image is better.

[0035] (4) Strong versatility and applicability: It is independent of the ArcGIS platform, directly processes raster images, and is suitable for downscaling large-scale, long-term series, and multi-magnification raster data. Attached Figure Description

[0036] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0037] Figure 1 This invention provides an overall technical flowchart of an IDW interpolation method for improving the resolution of raster images.

[0038] Figure 2 A schematic diagram showing the coordinates of the center point of the reference pixel within a 3×3 sliding window;

[0039] Figure 3 A schematic diagram of the target pixel distribution after subdivision of the center grid of the sliding window;

[0040] Figure 4 This is a comparison chart of interpolation results;

[0041] Figure 5 A detailed comparison of the two interpolation algorithms and the original image results;

[0042] Figure 6 Comparison of histograms for core indicators. Detailed Implementation

[0043] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be described in detail and completely below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. Generally, the components of the embodiments of the present invention described and shown in the accompanying drawings can be arranged and designed using different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed application, but merely represents selected embodiments of the present invention.

[0044] This invention provides an IDW interpolation method suitable for improving the resolution of raster images, focusing on the downscaling of remote sensing images. It can be applied to raster image processing in fields such as meteorology and ecological environment. This embodiment uses ERA5-Land2m temperature TIFF raster imagery as experimental data, with an original spatial resolution of 0.1°×0.1°. Three experimental regions, region1 to region3, with progressively increasing areas, are selected: Chengdu, Sichuan, and Southwest China, respectively. The number of raster cells in the three regions are 651, 15561, and 61341, respectively. Under the same hardware environment and parameter configuration, the proposed method (i.e., the improved IDW interpolation algorithm) is compared with the built-in IDW interpolation algorithm in ArcGIS. The computational efficiency differences of five sliding window sizes (3×3, 5×5, 7×7, 9×9, and 11×11) at magnification ratios of 5x, 10x, and 20x are verified.

[0045] See Figure 1 The specific implementation steps of the method proposed in this invention include:

[0046] 1. Data Preprocessing

[0047] 1.1 Parameter Configuration

[0048] Pre-set the input / output file paths, sliding window size, image magnification, IDW power parameter, and Gaussian filter smoothing parameter.

[0049] 1.2 Data Parsing and Format Conversion

[0050] The input raw raster image is directly converted into a floating-point two-dimensional numerical matrix by reading the TIFF raster image. Simultaneously, the spatial resolution and georeferenced information of the raster image are analyzed to extract the original pixel resolution and spatial coordinate parameters, providing a benchmark for subsequent pixel width calculation and georeferenced reconstruction after downscaling. The floating-point two-dimensional numerical matrix provides the carrier for subsequent vector operations and broadcast multiplication, enabling direct processing of the raster image.

[0051] 2. Weight template pre-calculation

[0052] 2.1 Generating the coordinate matrix

[0053] Using the top-left corner of the first reference pixel in the sliding window as the origin, the coordinates of the center point of each reference pixel are generated to represent the position of the reference pixel. Based on the sliding window size and the set magnification, target pixels are generated by subdividing the window center grid by magnification multiplied by the magnification factor. The position of the target pixel is then represented by the center coordinates of each target pixel.

[0054] The coordinates of the target pixel are calculated as shown in formulas (2-1) and (2-2), where x ij The x-coordinate of the center point of the target cell in the i-th row and j-th column of the subdivided grid is represented by y. ij The column coordinates are represented by t, the sliding window size is represented by t, and the magnification factor is represented by r.

[0055] (2-1)

[0056] (2-2)

[0057] Taking a 3x3 window as an example, Figure 2 The coordinates of the center point of each window grid are displayed. Figure 3 The results of subdividing the window's center grid are shown, where the original pixels within the window are reference pixels, and the smaller pixels after subdivision of the center grid are target pixels. The coordinates of the center point of each pixel are manually set to represent its position, allowing for the calculation of the distance from each target pixel to each reference pixel.

[0058] 2.2 Distance Matrix Calculation

[0059] Based on the distance formula, the distance matrix from all target pixels to the reference pixels is calculated. A minimum distance threshold δ is set, and the minimum values ​​in the distance matrix that are less than the threshold are replaced with δ to avoid the denominator being 0 in subsequent calculations. As shown in formula (2-3), where d ij This represents the distance from the i-th target pixel to the j-th reference pixel.

[0060] (2-3)

[0061] 2.3 Weight Calculation and Normalization

[0062] The initial weight matrix is ​​calculated using formula (2-4), where w ij d represents the inverse distance weight of the i-th reference pixel to the j-th target pixel. ij Let W represent the distance from the i-th reference pixel to the j-th target pixel, and p represent the IDW power parameter. Then, the initial weight matrix is ​​normalized to ensure the sum of the weights is 1, guaranteeing the rationality of the weighted summation, as shown in formula (2-5), where W... ijThese are the normalized weight values. The normalized weights corresponding to all target pixels are organized into a three-dimensional structure in the form of (number of rows in the center grid, number of columns in the center grid, and number of reference pixels), forming a three-dimensional weight template that can be pre-computed and reused globally.

[0063] (2-4)

[0064] (2-5)

[0065] 3. Sliding window interpolation

[0066] 3.1 Sliding Window Traversal

[0067] The two-dimensional numerical matrix of the original raster image is traversed using a set sliding window size, and a local window matrix centered on each original pixel is cropped.

[0068] 3.2 Edge Reflection Filling

[0069] For local window matrices where the window size is insufficient after edge pixel clipping, a reflection fill mode is used to fill the blank areas, making the edge window matrix consistent with the set sliding window size, thus solving the technical problem of low edge pixel interpolation accuracy in traditional IDW interpolation.

[0070] 3.3 Broadcast Multiplication Interpolation Calculation

[0071] The local window matrix is ​​flattened into a one-dimensional array. Based on the characteristics of vectorized broadcast multiplication, the pre-calculated three-dimensional weight template is batch-operated with the one-dimensional window array, and the weighted sum is calculated to obtain the values ​​of all target pixels, as shown in formula (3-1), where v i u represents the value of the i-th target pixel. j W represents the value of the j-th reference cell. ij This represents the normalized weight value of the j-th reference pixel to the i-th target pixel. Broadcast multiplication eliminates nested loops through vectorized operations, reducing the computational time complexity from O(n²) to O(n).

[0072] (3-1)

[0073] 3.4 Generating High-Resolution Matrices

[0074] The target pixel value after subdividing the center grid of each original pixel is written into the corresponding position of the high-resolution two-dimensional numerical matrix, and finally a high-resolution raster image with a size of original number of rows × magnification factor and original number of columns × magnification factor is generated.

[0075] 4. Post-processing of results

[0076] 4.1 Gaussian Filtering Smoothing

[0077] Gaussian filtering is applied to the generated high-resolution raster image to eliminate local block noise and bullseye effect caused by window shape and center grid subdivision, thereby improving the spatial continuity of the raster image.

[0078] 4.2 Results Output

[0079] Based on the original pixel resolution and downscaling factor, the spatial reference information of the output image is updated, the high-resolution raster data is saved in the universal TIFF format, and the algorithm configuration parameters and running indicators are recorded in a document to achieve complete preservation of the results.

[0080] 5. Test Results and Analysis

[0081] 5.1 Comparison of computational efficiency

[0082] Table 1. Magnification 5x (unit: seconds)

[0083]

[0084] Table 2 Magnification 10x (Unit: seconds)

[0085]

[0086] Table 3 Magnification 20x (Unit: seconds)

[0087]

[0088] The three tables above show the test time results. My_IDW represents the improved algorithm, and ArcGIS_IDW represents the ArcGIS built-in IDW interpolation algorithm. As can be seen from the time data in Tables 1, 2, and 3, the computation time of the improved IDW interpolation algorithm is significantly lower than that of the ArcGIS built-in IDW interpolation algorithm in all test scenarios. Moreover, the time change trends of the two are significantly different. The time of the improved algorithm increases steadily and linearly with the increase of area, magnification, and window size, while the time of the ArcGIS built-in algorithm increases exponentially and linearly, fully demonstrating the efficiency advantage of the improved algorithm.

[0089] Increasing the window size means increasing the number of reference pixels. ArcGIS's built-in algorithm needs to recalculate the distance and weight of more reference pixels for each target pixel, resulting in an exponential increase in repetitive calculations. In contrast, the improved algorithm achieves one-time calculation and global reuse through weight template pre-calculation. The window size only affects the time taken for template pre-calculation and has no significant impact on subsequent interpolation calculations. Therefore, the increase in time is much lower than that of ArcGIS's built-in algorithm.

[0090] Based on three sets of test data, the improved IDW interpolation algorithm improves computational efficiency by 4 to 100 times compared to ArcGIS's built-in IDW interpolation algorithm. Specifically, the efficiency improvement is more than 50 times for small areas, high magnification, and large window sizes, and remains stable at more than 40 times for large areas, high magnification, and large window sizes. The efficiency improvement is still more than 4 times for low magnification and small window sizes, fully meeting the raster image downscaling efficiency requirements in different scenarios.

[0091] 5.2 Comparison of Interpolation Effects

[0092] To verify that the improved IDW interpolation algorithm can improve efficiency without sacrificing interpolation accuracy and optimize the spatial quality of the results, this test uses ERA5-Land 2m temperature raster imagery of Sichuan Province as the core object, employing 20x magnification and a 5×5 window. The downscaling effect of the improved algorithm is compared with that of ArcGIS's built-in IDW interpolation algorithm. The interpolation results are compared as follows: Figure 4 As shown, the original image serves as the benchmark for the interpolation results, visually reflecting the downscaling and restoration capabilities and optimization effects of the two algorithms. Visually, both the ArcGIS built-in algorithm interpolation results and the improved algorithm interpolation results highly match the spatial distribution characteristics of temperature in the original image. Figure 5 The three images can be clearly distinguished in terms of detail differences, and the improved algorithm interpolation results are better in terms of edge transition and texture continuity.

[0093] To further quantify and verify the interpolation accuracy of the improved algorithm, four core statistical indicators—maximum value, minimum value, average value, and standard deviation—were extracted from the original image, the interpolation results of the ArcGIS built-in algorithm, and the interpolation results of the improved algorithm. The numerical deviations of each algorithm's results with the original image were compared, and the quantitative results are shown in Table 4. As can be seen from the quantitative data in Table 4, the statistical indicators of the improved algorithm's interpolation results are highly consistent with those of the original image, with a standard deviation of 0.54%. The average value is completely consistent with the original image, and both indicators are lower than the numerical deviations of the ArcGIS built-in algorithm, with maximum and minimum value deviations of only 0.008% and 0.004%, respectively. Figure 6 The comparison of the core indicator histograms can more intuitively show the numerical differences of the three in the four key statistical indicators of maximum value, minimum value, average value and standard deviation, which fully verifies the high-precision restoration capability of the improved algorithm at the numerical level.

[0094] Table 4 Comparison of core statistical indicators of interpolation results (unit: K)

[0095]

[0096] Overall, the improved IDW interpolation algorithm not only enhances computational efficiency but also ensures consistency with the statistical features and spatial distribution of the original image. Furthermore, through dual optimization of edge reflection filling and Gaussian filtering, it improves the spatial continuity of the interpolation results. In terms of interpolation accuracy and result quality, it achieves the same effect as the built-in IDW interpolation algorithm in ArcGIS, thus improving computational efficiency while maintaining accuracy.

[0097] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the foregoing embodiments have described the present invention in detail, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the embodiments of the present invention.

[0098] The above descriptions are merely some embodiments of the present invention. Those skilled in the art can make various modifications and improvements without departing from the inventive concept of the present invention, and these all fall within the scope of protection of the present invention.

Claims

1. An IDW interpolation method suitable for improving the resolution of a raster image, characterized in that, Includes the following steps: Step 1, Data Preprocessing Input the raw raster image data and configure the difference parameters, including: sliding window size, image magnification, and IDW power parameter; The original raster image data is converted into a floating-point two-dimensional numerical matrix, and geographic coordinates and cell widths are extracted. Step 2, pre-calculation of weight template The coordinate matrix of the reference pixel is generated according to the sliding window size, and then each reference pixel is subdivided according to the image magnification to obtain multiple target pixels, and the coordinate matrix of the target pixels is generated. Calculate the distance matrix from the target pixel to the reference pixel, set the minimum distance threshold δ, reset the distance elements in the distance matrix with values ​​less than δ to δ, then calculate the weights according to the IDW formula and normalize them, and then reorganize the weights to generate a reusable 3D weight template. Step 3, sliding window interpolation Traverse the two-dimensional numerical matrix by sliding window size, and perform boundary padding on the edge of the current window to obtain the padded local window matrix, and then flatten it into a one-dimensional array. Broadcast multiplication is used to perform batch operations on the three-dimensional weight template and the one-dimensional array corresponding to the window matrix to calculate the target pixel value. Then, each target pixel is stitched together to generate a high-resolution raster image. Based on the extracted geographic coordinates and pixel width, georeference restoration is performed on the generated high-resolution raster image.

2. The method of claim 1, wherein, In step 2, the coordinate matrix of the reference pixel is specifically as follows: Using the top-left corner of the first reference cell in the sliding window as the origin, the coordinates of the center point of each reference cell are generated to represent the position of the reference cell; Based on the image magnification factor r, the center grid of the current sliding window is subdivided into r×r target pixels, and the center coordinates of each target pixel represent the position of the target pixel.

3. The method of claim 2, wherein, The location of the target pixel is: in, , The first subdivision of the mesh Line number The x and y coordinates of the center point of the target cell in the column. Sliding window size.

4. The method as described in claim 1, characterized in that, In step 2, the calculation and normalization of weights according to the IDW formula are as follows: Based on the distance from i reference pixels to the j-th target pixel Calculate initial weights , Among them, among them, The preset IDW exponentiation parameter; Then adjust the initial weights Normalization is performed to obtain normalized weights. , where n is the number of target pixels.

5. The method as described in claim 4, characterized in that, In step 2, the reusable 3D weight template is: normalized weights The current sliding window is organized into a three-dimensional structure based on the number of rows, columns, and reference pixels of the center grid of the current sliding window, generating a reusable three-dimensional weight template.

6. The method as described in claim 5, characterized in that, In step 3, the target pixel value is: in, Let i be the pixel value of the i-th target pixel. The pixel is the j-th reference pixel.

7. The method as described in claim 1, characterized in that, Step 3 also includes: after eliminating the generated high-resolution raster image using Gaussian filtering, georeferenced restoration is then performed.

Images