第5章 Python 数字图像处理(DIP) - 图像复原与重建11 - 空间滤波 - 自适应滤波器 - 自适应局部降噪、自适应中值滤波器

自适应滤波器

自适应局部降噪滤波器

均值是计算平均值的区域上的平均灰度，方差是该区域上的图像对比度

$g(x,y)$噪声图像在$(x,y)$处的值
${\sigma }_{\eta }^2$ 为噪声的方差，为常数，需要通过估计得到
${\bar{z}}_{S_{xy}}$ 为局部平均灰度
${\sigma }_{S_{xy}}^2$ 为局部方差

$$\hat{f}(x,y)=g(x,y)-\frac{{\sigma }_{\eta }^2}{{\sigma }_{S_{xy}}^2}[g(x,y)-{\bar{z}}_{S_{xy}}]\text{\hspace{1em}(5.32)}$$

当${\sigma }_{\eta }^2>{\sigma }_{\eta }^2$ 时，将比率设置为1，这样做会使得滤波器是非线性的，但可阻止因缺少图像噪声方差的知识而产生无意义的结果（即负灰度级）。

若${\sigma }_{\eta }^2$估计值太低时，算法会因校正量小于就有的值而返回与原图像接近的图像。估计值太高会使得方差的比率在1.0处被水平，与正常情况相比，算法会更频繁地从图像中减去平均值。若允许为负值，且最后重新标定图像，则如前所述，结果 将损失图像的动态范围。

整个图像的相关值：
${\mu }_n={\sum }_{i=0}^{L-1}(r_i-m{)}^{n}p(r_i)$ 为灰度值r相对于其均值同的第n阶矩
$m={\sum }_{i=0}^{L-1}{r}_{i}p(r_i)$ 均值
${\sigma }^2={\mu }_2={\sum }_{i=0}^{L-1}(r_i-m{)}^{2}p(r_i)$ 方差

邻域内的相关值：
$m_{S_{xy}}={\sum }_{i=0}^{L-1}{r}_i{P}_{S_{xy}}(r_i)$ 均值
${\sigma }_{S_{xy}}^2={\mu }_2={\sum }_{i=0}^{L-1}(r_i-m_{S_{xy}}{)}^2{P}_{S_{xy}}(r_i)$ 方差

def calculate_sigma_eta(image):
    """
    :param image: input image
    :return: sigma eta of image
    """
    hist, bins = np.histogram(image.flatten(), bins=256, range=[0, 256], density=True)
    
    r = bins[:-1]
    m = np.sum(r * hist)
    sigma = np.sum((r - m)**2 * hist)
    
    return sigma

def calculate_sigma_sxy(block):
    """
    param block: input block
    return: sigma sxy of block
    """
    #==========公式法==========
    hist, bins = np.histogram(block.flatten(), bins=256, range=[0, 256], density=True)
    
    r = bins[:-1]
    m = (r * hist).sum()
    sigma = ((r - m)**2 * hist).sum()
    
    #===========等价于========
#     m = np.mean(block)
#     sigma = np.var(block)
    
    return sigma, m

def adaptive_local_denoise(image, kernel, sigma_eta=1):
    """
    adaptive local denoising 
    math: $$\hat{f}(x,y) = g(x,y) - \frac{\sigma_{\eta} ^ 2 }{\sigma_{S_{xy}}^2} \big[ g(x,y) - \bar z_{S_{xy}}\big]$$
    param: image:  input image for denoising
    param: kernel: input kernel, actually only use kernel shape, just want to keep the format as mean filter
    return: image after adaptive local denoising 
    """
    epsilon = 1e-8
    height, width = image.shape[:2]
    m, n = kernel.shape[:2]

    padding_h = int((m -1)/2)
    padding_w = int((n -1)/2)
    
    # 这样的填充方式，可以奇数核或者偶数核都能正确填充
    image_pad = np.pad(image, ((padding_h, m - 1 - padding_h), \
                                  (padding_w, n - 1 - padding_w)), mode="edge")
    
    img_result = np.zeros(image.shape)
    for i in range(height):
        for j in range(width):
            block = image_pad[i:i + m, j:j + n]
            gxy = image[i, j]
            z_sxy = np.mean(block)
            sigma_sxy = np.var(block)
            rate = sigma_eta / (sigma_sxy + epsilon)
            if rate >= 1:
                rate = 1
            img_result[i, j] = gxy - rate * (gxy - z_sxy)
    return img_result

# 自适应局部降噪滤波器处理高斯噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0513(a)(ckt_gaussian_var_1000_mean_0).tif', 0) #直接读为灰度图像

kernel = np.ones([7, 7])
img_arithmentic_mean = arithmentic_mean(img_ori, kernel=kernel)
img_geometric_mean = geometric_mean(img_ori, kernel=kernel)
img_adaptive_local = adaptive_local_denoise(img_ori, kernel=kernel, sigma_eta=1000)

plt.figure(figsize=(10, 10))

plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('With Gaussian noise'), plt.xticks([]),plt.yticks([])
plt.subplot(222), plt.imshow(img_arithmentic_mean, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]),plt.yticks([])
plt.subplot(223), plt.imshow(img_geometric_mean, 'gray'), plt.title('Geomentric mean'), plt.xticks([]),plt.yticks([])
plt.subplot(224), plt.imshow(img_adaptive_local, 'gray'), plt.title('Adaptive local denoise'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()

自适应中值滤波器

$z_{\text{min}}是S_{xy}$邻域中的最小灰度值；$z_{\text{max}}是S_{xy}$邻域中的最大灰度值；$z_{\text{med}}是S_{xy}$邻域中的中值， $z_{xy}$是坐标$(x,y)$处的灰度值；$S_{\text{max}}$是$S_{xy}$允许的最大尺寸，是大于1的奇正整数。

自适应中值滤波算法在点$(x,y)$处使用两个处理层次，分别表示为层次$A$和层次$B$:

层次$A$:
若$z_{\text{min}}<z_{\text{med}}<z_{\text{max}}$ 则转到层次B
否则，增$S_{xy}$的尺寸，
若$S_{xy}\le S_{\text{max}}$, 则重复层次A
否则，输出$z_{\text{med}}$

层次$B$：
若$z_{\text{min}}<z_{xy}<z_{\text{max}}$，则输出$z_{xy}$
否则，输出 $z_{\text{med}}$

这算法有3个主要的目的：去除椒盐（冲激）噪声，平滑其他非常冲激噪声，减少失真（如目标边界的过度细化）。该算法统计上认为$z_{\text{min}}$和$z_{\text{max}}$是区域$S_{xy}$的“类冲激”噪声分量，即使它们不是图像中的最小像素和最大像素值。

能很好的处理椒盐噪声，但对高斯噪声不敏感

def adaptive_median_denoise(image, sxy=3, smax=7):
    """
    adaptive median denoising 
    param: image: input image for denoising
    param: sxy  : minimum kernel size
    param: smax : maximum kernel size
    return: image after adaptive median denoising 
    """
    epsilon = 1e-8
    height, width = image.shape[:2]
    m, n = smax, smax

    padding_h = int((m -1)/2)
    padding_w = int((n -1)/2)
    
    # 这样的填充方式，可以奇数核或者偶数核都能正确填充
    image_pad = np.pad(image, ((padding_h, m - 1 - padding_h), \
                                  (padding_w, n - 1 - padding_w)), mode="edge")
    
    img_new = np.zeros(image.shape)
    
    for i in range(padding_h, height + padding_h):
        for j in range(padding_w, width + padding_w):
            sxy = 3     #每一轮都重置  
            k = int(sxy/2)
            block = image_pad[i-k:i+k+1, j-k:j+k+1]
            zxy = image[i - padding_h][j - padding_w]
            zmin = np.min(block)
            zmed = np.median(block)
            zmax = np.max(block)
            
            if zmin < zmed < zmax:
                if zmin < zxy < zmax:
                    img_new[i - padding_h, j - padding_w] = zxy
                else:
                    img_new[i - padding_h, j - padding_w] = zmed
            else:
                while True:
                    sxy = sxy + 2
                    k = int(sxy / 2)
                    
                    if zmin < zmed < zmax or sxy > smax:
                        break
                        
                    block = image_pad[i-k:i+k+1, j-k:j+k+1]
                    zmed = np.median(block)
                    zmin = np.min(block)
                    zmax = np.max(block)
                        
                if zmin < zmed < zmax or sxy > smax:
                    if zmin < zxy < zmax:
                        img_new[i - padding_h, j - padding_w] = zxy
                    else:
                        img_new[i - padding_h, j - padding_w] = zmed
                        
    return img_new

# 自适中值滤波器处理椒盐噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0514(a)(ckt_saltpep_prob_pt25).tif', 0) #直接读为灰度图像

kernel = np.ones([7, 7])

img_arithmentic_mean = median_filter(img_ori, kernel=kernel)
img_adaptive_median = adaptive_median_denoise(img_ori)

plt.figure(figsize=(15, 10))
plt.subplot(231), plt.imshow(img_ori, 'gray'), plt.title('Salt pepper noise'), plt.xticks([]),plt.yticks([])
plt.subplot(232), plt.imshow(img_arithmentic_mean, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]),plt.yticks([])
plt.subplot(233), plt.imshow(img_adaptive_median, 'gray'), plt.title('Adaptive median'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()

# 自适中值滤波器处理高斯噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0513(a)(ckt_gaussian_var_1000_mean_0).tif', 0) #直接读为灰度图像

kernel = np.ones([7, 7])
img_arithmentic_mean = median_filter(img_ori, kernel=kernel)
img_adaptive_median = adaptive_median_denoise(img_ori)

plt.figure(figsize=(15, 10))

plt.subplot(231), plt.imshow(img_ori, 'gray'), plt.title('With Gaussian noise'), plt.xticks([]),plt.yticks([])
plt.subplot(232), plt.imshow(img_arithmentic_mean, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]),plt.yticks([])
plt.subplot(233), plt.imshow(img_adaptive_median, 'gray'), plt.title('Adaptive median'), plt.xticks([]),plt.yticks([])
plt.tight_layout()
plt.show()