《数字图像处理》第五章 图像复原与重建学习笔记（5.1~5.3）

请注意:笔记内容片面粗浅，请读者批判着阅读！

一、 图像复原模型与噪声模型

核心概念

图像复原旨在通过建立退化模型，从退化的观测图像中恢复原始图像。其核心公式为：
$$g(x,y)=h(x,y)\ast f(x,y)+\eta (x,y)$$
其中：

• $$g(x,y)$$：退化图像
• $$h(x,y)$$：点扩散函数（PSF，退化模型）
• $$f(x,y)$$：原始图像
• $$\eta (x,y)$$：加性噪声

与图像增强的区别：复原基于退化模型的逆过程，需要先验知识（如噪声类型、PSF），而增强无需建模。

常见噪声模型

1. 高斯噪声：概率密度函数为钟形曲线，常见于传感器噪声。
   $$p(z)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(z-\mu {)}^2}{2{\sigma }^2}}$$

2. 椒盐噪声：随机出现黑白像素点，模拟传感器失效。

3. 瑞利噪声：非对称分布，多用于雷达图像。

4. 均匀噪声：概率密度在区间内恒定。

高斯噪声例子：

import cv2
import numpy as np

def add_gaussian_noise(image, mean=0, var=0.01):
    image = image / 255.0
    noise = np.random.normal(mean, var**0.5, image.shape)
    noisy = np.clip(image + noise, 0, 1.0)
    return np.uint8(noisy * 255)

img = cv2.imread('input.jpg', 0)
noisy_img = add_gaussian_noise(img, var=0.04)
cv2.imwrite('noisy.jpg', noisy_img)

二、空间域滤波复原

常用滤波器

1.均值滤波器：

算术均值滤波：简单平均邻域像素，适用于高斯噪声。

原理：用邻域像素的算术平均值替代中心像素，适用于高斯噪声和均匀噪声。
公式：
$$\hat{f}(x,y)=\frac{1}{m\times n}\sum _{(i,j)\in S_{xy}}g(i,j)$$

import cv2

img = cv2.imread('1noisy_image.png', 0)
# 使用5x5算术均值滤波
mean_filtered = cv2.blur(img, (5, 5))
# cv2.imwrite('arithmetic_mean.jpg', mean_filtered)
cv2.imshow('arithmetic_mean', mean_filtered)
cv2.waitKey(0)

几何均值滤波：对噪声抑制更平滑，但计算复杂。

原理：对邻域像素的乘积取几何平均，抑制噪声更平滑，但计算复杂度高。
公式：
$$\hat{f}(x,y)={\left(\prod _{(i,j)\in S_{xy}}g(i,j)\right)}^{1/(m\times n)}$$

import cv2
import numpy as np


def geometric_mean_filter(image, kernel_size=3):
    pad = kernel_size // 2
    padded = cv2.copyMakeBorder(image, pad, pad, pad, pad, cv2.BORDER_REFLECT)
    filtered = np.zeros_like(image, dtype=np.float32)
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            region = padded[i:i+kernel_size, j:j+kernel_size]
            product = np.prod(region.astype(np.float64) + 1e-6)  # 避免零值
            filtered[i, j] = product ** (1/(kernel_size**2))
    return np.uint8(filtered)


img = cv2.imread('1noisy_image.png', 0)
geo_result = geometric_mean_filter(img)
# cv2.imwrite('geometric_mean.jpg', geo_result)
cv2.imshow('geometric_mean', geo_result)
cv2.waitKey(0)

谐波均值滤波

原理：对噪声的调和平均，适用于盐粒噪声（正脉冲）。
公式：
$$\hat{f}(x,y)=\frac{m\times n}{{\sum }_{(i,j)\in S_{xy}}\frac{1}{g(i,j)}}$$

import cv2
import numpy as np


def harmonic_mean_filter(image, kernel_size=3):
    pad = kernel_size // 2
    padded = cv2.copyMakeBorder(image, pad, pad, pad, pad, cv2.BORDER_REFLECT)
    filtered = np.zeros_like(image)
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            region = padded[i:i+kernel_size, j:j+kernel_size].astype(np.float32)
            filtered[i,j] = (kernel_size**2) / np.sum(1 / (region + 1e-6))  # 避免除以零
    return np.uint8(filtered)


img = cv2.imread('1noisy_image.png', 0)
harmonic_result = harmonic_mean_filter(img)
# cv2.imwrite('harmonic_mean.jpg', harmonic_result)
cv2.imshow('harmonic_mean', harmonic_result)
cv2.waitKey(0)

逆谐波均值滤波

原理：通过指数参数 $$Q$$ 控制对正/负脉冲噪声的响应，$$Q>0$$ 抑制盐粒噪声，$$Q<0$$ 抑制胡椒噪声。
公式：
$$\hat{f}(x,y)=\frac{{\sum }_{(i,j)\in S_{xy}}g(i,j{)}^{Q+1}}{{\sum }_{(i,j)\in S_{xy}}g(i,j{)}^Q}$$

import cv2
import numpy as np


def inverse_harmonic_mean_filter(image, Q=1.5, kernel_size=3):
    pad = kernel_size // 2
    padded = cv2.copyMakeBorder(image, pad, pad, pad, pad, cv2.BORDER_REFLECT)
    filtered = np.zeros_like(image)
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            region = padded[i:i+kernel_size, j:j+kernel_size].astype(np.float32)
            numerator = np.sum(region ** (Q + 1))
            denominator = np.sum(region ** Q) + 1e-6
            filtered[i,j] = numerator / denominator
    return np.uint8(filtered)


img = cv2.imread('1noisy_image.png', 0)
ihm_result = inverse_harmonic_mean_filter(img, Q=1.5)
# cv2.imwrite('inverse_harmonic.jpg', ihm_result)
cv2.imshow('inverse_harmonic', ihm_result)
cv2.waitKey(0)

2.统计排序滤波器：

中值滤波：取邻域中值，有效去除椒盐噪声。

from scipy.ndimage import median_filter
import cv2


img = cv2.imread('1noisy_image.png', 0)
median_result = median_filter(img, size=3)
cv2.imwrite('median_filtered.jpg', median_result)
cv2.imshow('median_filtered', median_result)
cv2.waitKey(0)

最大值和最小值滤波

最大值滤波：提取邻域最大值，去除胡椒噪声（暗点）。
最小值滤波：提取邻域最小值，去除盐粒噪声（亮点）。

# 最大值滤波
import cv2
import numpy as np


img = cv2.imread('1noisy_image.png', 0)
max_result = cv2.dilate(img, np.ones((3,3), np.uint8))
cv2.imwrite('max_filtered.jpg', max_result)
cv2.imshow('max_filtered', max_result)

# 最小值滤波
min_result = cv2.erode(img, np.ones((3,3), np.uint8))
cv2.imwrite('min_filtered.jpg', min_result)
cv2.imshow('min_filtered', min_result)
cv2.waitKey(0)

中点滤波

原理：取邻域最大值和最小值的均值，平衡噪声抑制与细节保留。
// TODO: 例子后期补上

修正阿尔法均值滤波：剔除极端值后取平均，适合混合噪声。

原理：剔除邻域中最大和最小的 $$d/2$$ 个值后取平均，适合混合噪声。
公式：
$$\hat{f}(x,y)=\frac{1}{m\times n-d}\sum _{(i,j)\in S_{xy}^{\prime }}g(i,j)$$

import cv2
import numpy as np


def alpha_trimmed_mean_filter(image, d=4, kernel_size=5):
    pad = kernel_size // 2
    padded = cv2.copyMakeBorder(image, pad, pad, pad, pad, cv2.BORDER_REFLECT)
    filtered = np.zeros_like(image)
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            region = padded[i:i+kernel_size, j:j+kernel_size].flatten()
            sorted_region = np.sort(region)
            trimmed = sorted_region[d//2 : len(sorted_region)-d//2]
            filtered[i,j] = np.mean(trimmed)
    return np.uint8(filtered)


img = cv2.imread('1noisy_image.png', 0)
alpha_result = alpha_trimmed_mean_filter(img, d=4)
# cv2.imwrite('alpha_trimmed.jpg', alpha_result)
cv2.imshow('alpha_trimmed', alpha_result)
cv2.waitKey(0)

3.自适应滤波器：

自适应局部降低噪声滤波

原理：根据局部噪声方差动态调整滤波强度。
公式：
$$\hat{f}(x,y)=g(x,y)-\frac{{\sigma }_{\eta }^2}{{\sigma }_L^2}(g(x,y)-{\mu }_L)$$
其中 $${\sigma }_L^2$$ 为局部方差，$${\mu }_L$$ 为局部均值。

import cv2
import numpy as np


def adaptive_local_filter(image, noise_var=100, kernel_size=7):
    mean = cv2.blur(image, (kernel_size, kernel_size))
    mean_sq = cv2.blur(image**2, (kernel_size, kernel_size))
    local_var = mean_sq - mean**2
    local_var = np.clip(local_var, 1e-6, None)  # 避免除以零
    result = image - (noise_var / local_var) * (image - mean)
    return np.uint8(np.clip(result, 0, 255))


img = cv2.imread('1noisy_image.png', 0)
adaptive_local_result = adaptive_local_filter(img, noise_var=100)
# cv2.imwrite('adaptive_local.jpg', adaptive_local_result)
cv2.imshow('adaptive_local', adaptive_local_result)
cv2.waitKey(0)

自适应中值滤波

原理：动态扩展窗口大小，区分噪声与细节，处理高密度椒盐噪声。

import cv2
import numpy as np


def adaptive_median_filter(image, max_window=7):
    output = image.copy()
    padded = cv2.copyMakeBorder(image, max_window//2, max_window//2,
                               max_window//2, max_window//2, cv2.BORDER_REFLECT)
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            window_size = 3
            while window_size <= max_window:
                half = window_size // 2
                region = padded[i:i+window_size, j:j+window_size]
                z_min, z_med, z_max = np.min(region), np.median(region), np.max(region)
                if z_min < z_med < z_max:
                    if z_min < image[i,j] < z_max:
                        output[i,j] = image[i,j]
                    else:
                        output[i,j] = z_med
                    break
                else:
                    window_size += 2  # 增大窗口
            else:
                output[i,j] = z_med  # 窗口已达最大值
    return output.astype(np.uint8)


img = cv2.imread('1noisy_image.png', 0)
adaptive_median_result = adaptive_median_filter(img)
# cv2.imwrite('adaptive_median.jpg', adaptive_median_result)
cv2.imshow('adaptive_median', adaptive_median_result)
cv2.waitKey(0)

总结