在R中计算分辨率不同的两幅图像的平均值

核心前提 #

由于两个Z矩阵的坐标分辨率差异极大（x1/y1步长1，x2/y2步长0.01），首先需要统一坐标网格，将低分辨率数据插值到高分辨率网格上，才能进行后续的区域匹配与计算。

1. 准备工作与网格统一 #

使用akima包的插值函数将低分辨率的z1转换为与z2同维度的矩阵：

# 安装并加载依赖包
if (!require(akima)) install.packages("akima")
library(akima)

# 原始数据定义
x1 = seq(1,10, by =1)
y1 = seq(1,10, by =1)
z1 = outer(x1,y1)

x2 = seq(7,12, by =.01)
y2 = seq(5,12, by =.01)
z2 = outer(x2,y2, FUN = "*")

# 插值z1到x2-y2的高分辨率网格
z1_interp = interp(x1, y1, z1, xo = x2, yo = y2)$z
# 填充插值后的NA值（z1未覆盖的区域设为0，不影响后续非重叠区域取值）
z1_interp[is.na(z1_interp)] = 0

2. 定位重叠区域 #

计算并转换重叠区域的坐标为网格索引，方便后续操作：

# 重叠区域的坐标边界
x_overlap_min = max(min(x1), min(x2))
x_overlap_max = min(max(x1), max(x2))
y_overlap_min = max(min(y1), min(y2))
y_overlap_max = min(max(y1), max(y2))

# 转换为网格索引
x_idx_overlap = which(x2 >= x_overlap_min & x2 <= x_overlap_max)
y_idx_overlap = which(y2 >= y_overlap_min & y2 <= y_overlap_max)

3. 无缓冲区的合并矩阵 #

直接在重叠区域取两个矩阵的算术平均值，非重叠区域保留各自原始有效数据：

# 初始化合并矩阵
z_combined_no_buffer = matrix(NA, nrow = length(x2), ncol = length(y2))

# 填充z1覆盖但z2未覆盖的区域
z1_mask = (x2 >= min(x1) & x2 <= max(x1)) & (y2 >= min(y1) & y2 <= max(y1))
z_combined_no_buffer[z1_mask] = z1_interp[z1_mask]

# 填充z2覆盖但z1未覆盖的区域
z2_mask = !z1_mask
z_combined_no_buffer[z2_mask] = z2[z2_mask]

# 重叠区域取均值
z_combined_no_buffer[x_idx_overlap, y_idx_overlap] = (z1_interp[x_idx_overlap, y_idx_overlap] + z2[x_idx_overlap, y_idx_overlap]) / 2

4. 带缓冲区的平滑合并 #

通过滑动窗口对重叠区域的两个矩阵做平滑处理，再取均值，实现缓冲区的过渡效果。这里采用3x3窗口的简单平滑：

# 定义滑动窗口平滑函数
smooth_matrix = function(mat, window_size = 3) {
  pad_width = floor(window_size / 2)
  # 给矩阵边缘填充NA，避免边界计算出错
  mat_padded = matrix(NA, nrow = nrow(mat) + 2*pad_width, ncol = ncol(mat) + 2*pad_width)
  mat_padded[(pad_width+1):(nrow(mat)+pad_width), (pad_width+1):(ncol(mat)+pad_width)] = mat
  
  smoothed_mat = matrix(NA, nrow = nrow(mat), ncol = ncol(mat))
  # 遍历每个窗口计算均值
  for (i in 1:nrow(mat)) {
    for (j in 1:ncol(mat)) {
      window = mat_padded[i:(i+window_size-1), j:(j+window_size-1)]
      smoothed_mat[i,j] = mean(window, na.rm = TRUE)
    }
  }
  return(smoothed_mat)
}

# 对插值后的z1和原始z2分别做平滑
z1_smoothed = smooth_matrix(z1_interp)
z2_smoothed = smooth_matrix(z2)

# 初始化带缓冲区的合并矩阵（继承非重叠区域的值）
z_combined_with_buffer = z_combined_no_buffer

# 重叠区域用平滑后的值取均值
z_combined_with_buffer[x_idx_overlap, y_idx_overlap] = (z1_smoothed[x_idx_overlap, y_idx_overlap] + z2_smoothed[x_idx_overlap, y_idx_overlap]) / 2

5. 结果可视化 #

对比两种合并结果的差异：

par(mfrow = c(1,2), mar = c(3,3,2,1))
image(x2, y2, z_combined_no_buffer, main = "无缓冲区合并", asp = 1, xlab = "", ylab = "")
image(x2, y2, z_combined_with_buffer, main = "带缓冲区平滑合并", asp = 1, xlab = "", ylab = "")

内容的提问来源于stack exchange，提问作者M. Beausoleil