插值以提高空间分辨率——反距离插值与双三次插值

前言

因为一些工作需求，进行了反距离插值与双三次插值的尝试，试图将图像分辨率体高20倍，这是一个比较大的倍率。现将两段代码稍加记录。

反距离插值

matlab代码

代码如下：

clear;
close all;
clc;

% 设置变量
data_folder = "E:\反距离权重\cut_new";
output_folder = "E:\反距离权重\M\"; 
num_nearest_neighbors = 8; 
p = 2; % 幂指数


% 循环遍历每个数据块
for i = 1:5
    for j = 1:4
        % 生成当前块的文件名
        data_filename = sprintf('data_block_%d-%d.mat', i, j);
        full_data_path = fullfile(data_folder, data_filename);
        
        % 加载当前块的数据
        loaded_data = load(full_data_path);
        data = loaded_data.data;
        lat = loaded_data.lat;
        lon = loaded_data.lon;

        % 获取当前块的行列数
        [rows, cols] = size(data);
        
        % 计算新的分辨率
        new_rows = ceil(rows * 20);
        new_cols = ceil(cols * 20);

        % 生成原始数据的网格
        [X, Y] = meshgrid(1:cols, 1:rows);

        % 生成新的插值网格
        [Xq, Yq] = meshgrid(linspace(1, cols, new_cols), linspace(1, rows, new_rows));

        % IDW插值
        Zq = idw_interpolation_optimized(X(:), Y(:), data(:), Xq, Yq, num_nearest_neighbors, p);

        % 重塑插值结果以匹配新网格
        Zq_reshaped = reshape(Zq, size(Xq));
        
        % 生成插值数据的保存路径
        interpolated_filename = sprintf('interpolated_data_block_%d-%d.mat', i, j);
        full_interpolated_path = fullfile(output_folder, interpolated_filename);

        % 保存插值结果
        save(full_interpolated_path, 'Zq_reshaped');
        
        % 打印进度
        fprintf('Block %d-%d processed and saved.\n', i, j);
    end
end

%定义函数
function Zq = idw_interpolation_optimized(X, Y, Z, Xq, Yq, num_nearest_neighbors, p)
    % 创建KD树以优化搜索
    ns = createns([X, Y], 'nsmethod', 'kdtree');
    eps = 1e-6;  % 防止除以零

    % 初始化输出矩阵
    Zq = zeros(size(Xq));

    % 标记NaN值的位置
    nan_mask = isnan(Z);

    % 遍历每个查询点进行插值
    for ix = 1:length(Xq(:))
        % 查询点
        queryPoint = [Xq(ix), Yq(ix)];
        
        % 使用KD树找到最近的num_nearest_neighbors个点
        [idx, dist] = knnsearch(ns, queryPoint, 'k', num_nearest_neighbors);

        % 排除NaN值
        valid_idx = idx(~nan_mask(idx));
        valid_dist = dist(~nan_mask(idx));

        if isempty(valid_idx)
            % 如果所有最近邻点都是NaN，则设置插值结果为NaN
            Zq(ix) = NaN;
        else
            % 获取这些点的Z值
            Z_values = Z(valid_idx);

            % 防止距离为零，这里使用eps来替换零距离
            valid_dist = max(valid_dist, eps);

            % 计算权重
            weights = 1 ./ (valid_dist.^p);

            % 计算加权平均
            Zq(ix) = sum(weights .* Z_values) / sum(weights);
        end
    end
end

matlab有内嵌的插值函数，内置四种插值方法，但是是不包括反距离插值的。这段代码中使用for循环嵌套进行插值点数据的计算，可能还有更好的方法。该代码在一开始设置了最近邻数量，即插值时考虑的最近邻数目，以及权重，通过调节这两个参数可以对插值效果以及运行速率进行调控。

除此之外，这段代码其实还可以添加限制搜索范围的语句，即在计算某网格点的数值时，只考虑与该网格点距离在搜索范围之内的数据影响，从而避免每次计算一个网格点都需要比较大量网点之间的距离。在这段代码中，因为没有设置搜索范围，在绘图结果中会看到在插值前表现为全部nan值的区域，在插值后出现很稀疏的数据点。因为这类点的数据来源来自距离较远的区域，我认为这种点并不包含有效数据。

其中，在python代码中（以下给出），距离在不同数据中具有不同含义，假如数据中没有提供其他信息，那么距离指代的就是网格的数目，假如提供了距离、经纬度等信息，那么距离指的可能就是m、km、°等。

另外，这段代码在nan值的处理上存在一定缺陷。该代码在插值时忽略nan值，在插值结束后，将插值得到数据的网格点之外的网格全部填充为nan，这实际上造成了在插值得到的网格中，有数据与无数据的区域的界限并没有移动，在绘制图像的视觉效果上表现为，数据区域边缘变得更加平滑，但是区域并没有扩大。有很大一部分的无效值点并不是因为没有值而体现为nan值，而是因为云遮档、耀斑等原因被掩掉了。换句话说，这种nan值的处理属于非常保守的一种。

python代码

以下给出python代码

import os
import numpy as np
import scipy.io
from scipy.spatial import cKDTree
from multiprocessing import Pool, cpu_count
import time
#IDW插值
# 设置特殊数值来替代NaN和插值参数
special_value = -9999.0  # 使用float类型的特殊值
num_nearest_neighbors = 8  # 最近邻数量
p = 2  # 反距离的幂次
scale_factor = 20  # 分辨率扩大的倍数

# 设置变量
data_folder = "E:\\反距离权重\\cut_new"
output_folder = "E:\\反距离权重\\M\\"

def load_and_preprocess_data(i, j):
    # 生成当前块的文件名
    data_filename = f'data_block_{i}-{j}.mat'
    full_data_path = os.path.join(data_folder, data_filename)

    # 加载当前块的数据
    loaded_data = scipy.io.loadmat(full_data_path)
    data = loaded_data['data']

    # 将NaN值替换为特殊数值，并更改数据类型为float32
    nan_mask = np.isnan(data)
    data = data.astype(np.float32)
    data[nan_mask] = special_value

    return data

def interpolate_block(i, j):
    data = load_and_preprocess_data(i, j)

    # 获取行列数并计算扩大后的大小
    rows, cols = data.shape
    new_rows, new_cols = rows * scale_factor, cols * scale_factor

    # 创建扩大后的数据数组
    expanded_data = np.full((new_rows, new_cols), special_value, dtype=np.float32)

    # 创建KD树
    X, Y = np.meshgrid(np.arange(cols), np.arange(rows), indexing='ij')
    valid_points = np.column_stack((X.ravel(), Y.ravel()))
    valid_values = data.ravel()
    tree = cKDTree(valid_points)

    # 添加时间监视
    t1 = time.time()

    # 对扩大后的每个点进行插值
    for new_x in range(new_cols):
        for new_y in range(new_rows):
            # 找到最近的原始数据点
            original_x, original_y = new_x // scale_factor, new_y // scale_factor

            # 如果原始数据点是特殊值，执行插值
            if data[original_y, original_x] == special_value:
                distances, indices = tree.query((original_x, original_y), k=num_nearest_neighbors)

                # 检查距离中是否包含零，并适当处理
                if np.any(distances == 0):
                    # 如果查询点与一个数据点重合，直接使用该点的值
                    zero_distance_index = np.where(distances == 0)[0][0]
                    interpolated_value = valid_values[indices[zero_distance_index]]
                else:
                    weights = 1 / np.power(distances, p)
                    interpolated_value = np.sum(weights * valid_values[indices]) / np.sum(weights)

                # 如果插值结果为负，赋值为NaN
                expanded_data[new_y, new_x] = interpolated_value if interpolated_value >= 0 else np.nan
            else:
                # 直接复制原始数据点
                expanded_data[new_y, new_x] = data[original_y, original_x]

    # 保存插值结果
    interpolated_filename = f'data_block_{i}-{j}_interpolated_expanded.npz'
    full_interpolated_path = os.path.join(output_folder, interpolated_filename)
    np.savez(full_interpolated_path, interpolated_data=expanded_data)

    # 计算处理耗时并打印
    t2 = time.time()
    print(f'Block {i}-{j} expanded and interpolated. Time taken: {t2 - t1} seconds')

def main():
    # 使用多进程并行处理每个数据块
    num_processes = cpu_count()  # 使用所有可用的CPU核心
    with Pool(processes=num_processes) as pool:
        tasks = [(i, j) for i in range(1, 6) for j in range(1, 5)]
        pool.starmap(interpolate_block, tasks)

if __name__ == '__main__':
    main()

如上所示，该代码中采用了更新的对于nan值的处理方法，即将nan赋值为-9999，在插值之后将所有0以下的数据全部赋值为nan。但是另一方面，该代码依然没有设置最大搜索范围，因而在绘图时会出现噪点。图像已经删除，无法展示，简单来说就是在下图左侧的大片空白区中，插值过后在很靠左的区域出现零星不规则分布的噪点。

双三次插值

这就没什么好说的了，matlab有内置的三次插值算法，只需要选择函数正确就能进行一维、二维、三维的三次插值。但是这段函数中没有对nan值点进行任何处理，那么实际上就是使用的插值算法的默认处理方式，即不考虑在内，效果与反距离插值的matlab代码一致，边缘平滑但是视觉效果上无效值占比增大。

clear; close all; clc; 

% 插值倍数
scale_factor = 20;

% 指定输入文件夹名称和输出文件夹名称
input_folder = "E:\反距离权重\cut_new";  % 这里替换为你的输入文件夹路径
output_folder = 'Cubic';

% 如果输出文件夹不存在，则创建它
if ~exist(output_folder, 'dir')
    mkdir(output_folder);
end

% 获取指定文件夹中所有.mat文件
files = dir(fullfile(input_folder, 'data_block_*.mat'));

% 循环处理每一个文件
for k = 1:length(files)
    file = files(k);
    
    % 获取完整的文件路径
    filename = fullfile(file.folder, file.name);
    
    % 加载数据
    loaded_data = load(filename);
    
    % 获取原始数据
    data = loaded_data.data;
    lat = loaded_data.lat;
    lon = loaded_data.lon;

    % 原始数据网格
    [X, Y] = meshgrid(lon, lat); % 创建二维网格
    
%     % 新的查询点网格
%     new_lat = linspace(min(lat), max(lat), scale_factor*size(lat,1));
%     new_lon = linspace(min(lon), max(lon), scale_factor*size(lon,2));
    % 新的查询点网格
new_lat = linspace(min(lat), max(lat), numel(lat) * scale_factor);
new_lon = linspace(min(lon), max(lon), numel(lon) * scale_factor);

    [Xq, Yq] = meshgrid(new_lon, new_lat); % 创建新的二维网格
    
    % 进行三次插值
    interpolated_data = interp2(X, Y, data, Xq, Yq, 'cubic');
    
    % 保存数据的新变量
    data_new = interpolated_data;  % 保持数据的新命名
    lat_new = new_lat;  % 新的纬度数组
    lon_new = new_lon;  % 新的经度数组
    
    % 构造保存文件的完整路径
    save_filename = fullfile(output_folder, file.name);

    
    % 保存插值后的数据
    save(save_filename, 'data_new', 'lat_new', 'lon_new', '-v7.3');

    
    % 输出状态
    fprintf('File %s has been processed and saved.\n', save_filename);
end