多维滤波

概述：

噪声对图像处理的影响很大，它影响图像处理的输入、采集和处理等各个环节以及输出结果。因此，在进行其它的图像处理前，需要对图像进行去噪处理。尤其在医学图像中可能有大量的3d数据，本文将从陆续实现2d，3d常用滤波。

1.均值滤波

均值滤波，是图像处理中最常用的手段，从频率域观点来看均值滤波是一种低通滤波器，高频信号将会去掉，因此可以帮助消除图像尖锐噪声，实现图像平滑，模糊等功能。理想的均值滤波是用每个像素和它周围像素计算出来的平均值替换图像中每个像素。采样Kernel数据通常是3X3的矩阵，如下表示：

红色为中心像素，计算九个像素的平均值，替换中心像素。

优点：算法简单，计算速度快；

缺点：降低噪声的同时使图像产生模糊，特别是景物的边缘和细节部分。

2.中值滤波

中值滤波也是消除图像噪声最常见的手段之一，特别是消除椒盐噪声，中值滤波的效果要比均值滤波更好。中值滤波是将周围像素和中心像素排序以后，取中值，一个3X3大小的中值滤波如下：

优点：抑制效果很好，画面的清析度基本保持；

缺点：对高斯噪声的抑制效果不是很好。

3.高斯滤波

高斯滤波器是一种线性滤波器，能够有效的抑制噪声，平滑图像。其作用原理和均值滤波器类似，都是取滤波器窗口内的像素的均值作为输出。其窗口模板的系数和均值滤波器不同，均值滤波器的模板系数都是相同的为1；而高斯滤波器的模板系数，则随着距离模板中心的增大而系数减小。所以，高斯滤波器相比于均值滤波器对图像个模糊程度较小。二维高斯函数如下：

一个3X3大小的高斯滤波如下：

快速将图像分成一块块的卷积块：

def conv2d_data2(data,k_size):
    padd_h = k_size[0]//2
    padd_w = k_size[1]//2
    x_pad = np.pad(data,((padd_h,padd_h),(padd_w,padd_w)))
    h_num = data.shape[0]
    w_num = data.shape[1]
    temp_h = np.array([x_pad[i:i+k_size[0],:] for i in range(h_num)])
    temp_w = np.array([temp_h]*k_size[1]).transpose([1,0,2,3])
    for i in range(1,k_size[1]):
        temp_w[:,i,:,:-i] = temp_w[:,i,:,i:]
    result = temp_w[:,:,:,:-(k-1)].transpose(0,3,2,1)
    return result

高斯2d快速滤波：

def GaussianFilter2d_speed(img,K_size,sigma):#单通道，K_size[h,w]
    #生成卷积块
    data = conv2d_data(img,K_size)
    
    #生成卷积核
    padd_h = K_size[0]//2
    padd_w = K_size[1]//2
    K = np.zeros((K_size[0],K_size[1]),dtype=np.float)
    for x in range(-padd_w,-padd_w+K_size[1]):
        for y in range(-padd_h,-padd_h+K_size[0]):
            K[y+padd_h,x+padd_w] = np.exp(-(x**2+y**2)/(2*(sigma**2)))
    K /= (sigma**2*2*np.pi)
    K /=  K.sum()
    
    #滤波
    result = np.squeeze(np.einsum("abef,fg->abeg",data.reshape(data.shape[0], data.shape[1], 1, -1) ,K.reshape(-1, 1)))
    return result

上述代码只需替换卷积和部分，就可以完成大部分的滤波操作。

三维高斯函数如下：

快速划分3d卷积块：

def conv3d_data2(data,k_size):#stride只能是1,pad是旁边填0,k_size[c,h,w]
    padd_c = k_size[0]//2
    padd_h = k_size[1]//2
    padd_w = k_size[2]//2
    x_pad = np.pad(data,((padd_c,padd_c),(padd_h,padd_h),(padd_w,padd_w)))
    c_num = data.shape[0]
    h_num = data.shape[1]
    w_num = data.shape[2]
    temp_c = np.array([x_pad[i:i+k_size[0],:,:] for i in range(c_num)])
    temp_h = np.array([temp_c[:,:,i:i+k_size[1],:] for i in range(h_num)])
    temp_w = np.array([temp_h]*k_size[2]).transpose([2,1,0,3,4,5])
    for i in range(1,k_size[2]):
        temp_w[:,:,i,:,:,:-i] = temp_w[:,:,i,:,:,i:]
    result = temp_w[:,:,:,:,:,:-(k_size[2]-1)].transpose([0,1,5,3,4,2])
    return result

高斯3d快速滤波：

def GaussianFilter3d_speed(img,K_size,sigma):#单通道，K_size[h,w]
    #生成卷积块
    data = conv3d_data2(img,K_size)
    
    #生成卷积核
    padd_c = K_size[0]//2
    padd_h = K_size[1]//2
    padd_w = K_size[2]//2
    K = np.zeros((K_size[0],K_size[1],K_size[2]),dtype=np.float)
    for x in range(-padd_w,-padd_w+K_size[2]):
        for y in range(-padd_h,-padd_h+K_size[1]):
            for z in range(-padd_c,-padd_c+K_size[0]):
                K[z+padd_c,y+padd_h,x+padd_w] = np.exp(-(x**2+y**2+z**2)/(2*(sigma**2)))
    K /= (sigma*np.sqrt(2*np.pi))**3
    K /=  K.sum()
    
    #滤波
    result = np.squeeze(np.einsum("abcdef,dfg->abcdeg",data.reshape(data.shape[0], data.shape[1], data.shape[2], 1, 1, -1) ,K.reshape(1, -1, 1)))
    return result

只需替换上述的卷积和部分，也就能完成大部分的3d滤波操作。[300,512,512]的数据大概需要10s。

上述代码是我经过多次优化后的结果，要是有大神有更加好的优化方式，欢迎私信我或者给我留言。