自定义高斯滤波代码与scipy同参数处理结果不一致排查

问题现象 #

以下代码用于实现高斯滤波功能，但将其运行结果与scipy库自带的gaussian_filter函数输出对比时，二者结果存在明显差异：使用相同标准差参数的情况下，自定义代码处理得到的图像比scipy处理结果更偏白。

• 自定义高斯滤波处理结果
  
• scipy高斯滤波处理结果
  
• 含噪声原图
  
• 两处理结果差值图
  

复现代码 #

from scipy import ndimage as nd
import numpy as np
from skimage import io, img_as_float
image = img_as_float(io.imread(r'C:\Users\Y\Downloads\MRI_clean.tif'))
gauss = np.random.normal(0,0.05,(940,934))
noisy = image + gauss
def convolution(oldimage, kernel):

    kernel_h = kernel.shape[0]

    if (len(oldimage.shape) == 3):
        image_pad = np.pad(oldimage, pad_width=(
        (kernel_h // 2, kernel_h // 2), (kernel_h // 2,kernel_h // 2), (0, 0)), mode = 'constant', constant_values = 0).astype(np.float32)
    elif (len(oldimage.shape) == 2):
        image_pad = np.pad(oldimage, pad_width=(
        (kernel_h // 2, kernel_h // 2), (kernel_h // 2, kernel_h // 2)), mode = 'constant', constant_values = 0).astype(np.float32)

    h = kernel_h // 2

    image_conv = np.zeros(image_pad.shape)

    for i in range(h, image_pad.shape[0] - h):
        for j in range(h, image_pad.shape[1] - h):
            x = image_pad[i - h:i - h + kernel_h, j - h:j - h + kernel_h]
            x = x * kernel
            image_conv[i][j] = x.sum()
    h_end = -h

    if (h == 0):
         return image_conv[h:, h:h_end]

    return image_conv[h:h_end, h:h_end]

def f(x,y,sigma):
    return (1/sqrt(2*math.pi*(sigma**2))) * np.exp(-(x**2+y**2)/(2*sigma**2))

def gaussian_kernel(size,sigma):
    kernel = np.zeros((size,size))
    n = (size - 1)/2
    for i in range(size):
        for j in range(size):
            kernel[i, j] = f(i-n, j-n, sigma)
    return kernel/(np.sum(kernel))

gaussian = nd.gaussian_filter(noisy, sigma = 5)
filtred =convolution(noisy,gaussian_kernel(31, 5))

问题原因 #

结果整体偏白的核心原因是卷积函数对多通道图像的求和逻辑错误，同时代码存在其他几个影响结果一致性的问题：

1. 核心错误：多通道卷积时错误聚合通道维度
   读取的MRI图像为3通道RGB格式（即使三个通道像素值完全相同），img_as_float读入后数组形状为(940,934,3)。生成噪声时使用了2D形状(940,934)，通过numpy广播机制可以直接和3通道图像相加，不会抛出维度不匹配错误，因此该问题不容易被发现。在卷积计算时，x.sum()没有指定求和轴，会把31×31窗口内所有元素（包括3个通道的所有值）全部求和得到一个标量，再赋值给输出位置的3个通道，最终输出值约为正确值的3倍，直接导致图像整体偏白过曝。
2. 缺失依赖导入
   高斯核计算函数f中用到了math.sqrt，但代码没有导入math模块，直接运行会抛出NameError。
3. 边界填充模式不匹配
   scipy的gaussian_filter默认边界填充模式为reflect（镜像反射边缘像素），自定义卷积使用constant模式填充0，二者在边缘区域的计算结果会存在差异。
4. 高斯核尺寸不匹配
   scipy默认高斯核的截断半径为4*sigma，sigma=5时自动生成的核半径为20、尺寸为41；自定义使用的是尺寸31、半径15的核（截断半径3*sigma），模糊程度和scipy默认结果有差异。

修复方案 #

针对以上问题修改代码即可得到和scipy一致的结果：

1. 补上import math导入；
2. 卷积求和时，对多通道图像指定仅在空间维度（0、1轴）求和，把x.sum()改为x.sum(axis=(0,1))；
3. 如果要完全对齐scipy默认结果，把填充模式改为reflect，同时将高斯核尺寸改为41（对应4*sigma截断半径）。

修复后的核心代码示例：

# 补上缺失的导入
import math

def convolution(oldimage, kernel):
    kernel_h = kernel.shape[0]
    if (len(oldimage.shape) == 3):
        # 改为reflect填充对齐scipy默认逻辑
        image_pad = np.pad(oldimage, pad_width=(
        (kernel_h // 2, kernel_h // 2), (kernel_h // 2,kernel_h // 2), (0, 0)), mode = 'reflect').astype(np.float32)
    elif (len(oldimage.shape) == 2):
        image_pad = np.pad(oldimage, pad_width=(
        (kernel_h // 2, kernel_h // 2), (kernel_h // 2, kernel_h // 2)), mode = 'reflect').astype(np.float32)

    h = kernel_h // 2
    image_conv = np.zeros(image_pad.shape)
    for i in range(h, image_pad.shape[0] - h):
        for j in range(h, image_pad.shape[1] - h):
            x = image_pad[i - h:i - h + kernel_h, j - h:j - h + kernel_h]
            x = x * kernel
            # 多通道场景下指定仅在空间维度求和
            image_conv[i][j] = x.sum(axis=(0,1))
    h_end = -h
    if (h == 0):
         return image_conv[h:, h:h_end]
    return image_conv[h:h_end, h:h_end]

# 调用时核大小改为41，对齐scipy默认4*sigma截断规则
filtred = convolution(noisy,gaussian_kernel(41, 5))

内容的提问来源于stack exchange，提问作者Parados Hiver