一、章节概述

图像复原与重建是数字图像处理的核心内容之一，旨在将退化的图像恢复到原始状态。与图像增强不同，图像复原是基于数学模型的客观恢复过程，而增强是主观的视觉改善过程。本文将从退化模型、噪声分析、各类复原算法到图像重建，全方位讲解图像复原技术，并通过 Python 代码实现所有关键算法的可视化对比。

二、核心内容讲解

5.1 图像退化 / 复原过程模型

5.2 噪声模型

5.2.1 噪声的空间域与频域特性

• 空间域：噪声表现为像素值的随机波动
• 频域：不同噪声在频域有不同的分布（如高斯噪声是均匀分布，周期性噪声是离散的峰值）

5.2.2 常用噪声概率密度函数

5.2.3 周期性噪声

由电力线、机械振动等周期性干扰引起，在频域表现为离散的冲激点，可通过频域滤波消除。

5.2.4 噪声参数估计

常用方法：

1. 从图像的平坦区域（背景）估计
2. 使用统计方法（均值、方差、直方图拟合）

5.3 仅含噪声的图像复原 —— 空间域滤波

当图像仅受噪声退化时，可通过空间域滤波实现复原。

5.3.1 均值滤波器

• 算术均值滤波：用邻域内像素的平均值替换中心像素
• 几何均值滤波：用邻域内像素的几何平均值替换中心像素
• 谐波均值滤波：适合消除盐粒噪声
• 逆谐波均值滤波：可根据 Q 值消除盐粒 / 胡椒噪声

5.3.2 统计排序滤波器

• 中值滤波：最常用，对椒盐噪声效果极佳
• 最大值滤波：消除胡椒噪声
• 最小值滤波：消除盐粒噪声
• 中点滤波：结合排序和平均，适合高斯 / 均匀噪声

5.3.3 自适应滤波器

根据局部区域的统计特性（均值、方差）自适应调整滤波参数，效果优于固定滤波器。

5.4 基于频域滤波的周期性噪声抑制

5.4.1 带阻滤波器

抑制特定频率范围的噪声，分为理想、巴特沃斯、高斯带阻滤波器。

5.4.2 带通滤波器

允许特定频率范围通过，可用于提取周期性噪声。

5.4.3 陷波滤波器

抑制 / 通过特定频率点（而非频率范围），是消除周期性噪声最常用的滤波器。

5.4.4 最优陷波滤波

结合多个陷波滤波器，自适应消除周期性噪声。

5.5 线性移不变退化

5.6 退化函数的估计

5.6.1 基于图像观测的估计方法

从退化图像中选取平坦区域，估计 PSF。

5.6.2 基于实验的估计方法

使用与退化图像相同的成像系统，对已知目标成像，直接测量 PSF。

5.6.3 基于建模的估计方法

根据物理过程（如运动模糊、大气湍流）建立数学模型，推导 PSF。

5.7 逆滤波

5.8 最小均方误差（维纳）滤波

5.9 约束最小二乘滤波

5.10 几何均值滤波

5.11 基于投影的图像重建

5.11.1 引言

基于投影的重建是 CT、MRI 等医学成像的核心技术，通过多角度投影数据重建断层图像。

5.11.2 计算机断层扫描（CT）原理

CT 通过 X 射线源和探测器阵列绕人体旋转，采集不同角度的投影数据，再通过数学方法重建断层图像。

5.11.3 投影与拉东变换

5.11.4 傅里叶切片定理

5.11.5 基于平行束滤波反投影的重建

步骤：

1. 对每个角度的投影进行滤波
2. 将滤波后的投影反投影到图像空间
3. 对所有角度的反投影结果求和

5.11.6 基于扇形束滤波反投影的重建

针对扇形束扫描（实际 CT 常用）的重建方法，需要进行扇形束到平行束的转换。

三、完整 Python 代码实现

import cv2 import numpy as np import matplotlib.pyplot as plt from scipy import ndimage from scipy.fftpack import fft2, ifft2, fftshift, ifftshift import warnings warnings.filterwarnings('ignore') # ====================== 全局设置 ====================== # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans'] plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题 plt.rcParams['figure.figsize'] = (16, 12) # ====================== 1. 噪声生成函数 ====================== def add_gaussian_noise(image, mean=0, var=0.001): """添加高斯噪声""" image = np.float32(image) / 255.0 sigma = var ** 0.5 noise = np.random.normal(mean, sigma, image.shape) noisy_image = image + noise noisy_image = np.clip(noisy_image, 0, 1) return np.uint8(noisy_image * 255) def add_salt_pepper_noise(image, prob=0.05): """添加椒盐噪声""" noisy_image = np.copy(image) # 盐噪声（白色） salt = np.random.rand(*image.shape) < prob/2 noisy_image[salt] = 255 # 椒噪声（黑色） pepper = np.random.rand(*image.shape) < prob/2 noisy_image[pepper] = 0 return noisy_image def add_poisson_noise(image): """添加泊松噪声""" vals = len(np.unique(image)) vals = 2 ** np.ceil(np.log2(vals)) noisy_image = np.random.poisson(image * vals) / float(vals) return np.uint8(np.clip(noisy_image, 0, 255)) # ====================== 2. 空间域滤波函数 ====================== def arithmetic_mean_filter(image, kernel_size=3): """算术均值滤波""" return cv2.blur(image, (kernel_size, kernel_size)) def geometric_mean_filter(image, kernel_size=3): """几何均值滤波""" image = np.float32(image) + 1e-8 # 避免0值 kernel = np.ones((kernel_size, kernel_size)) log_image = np.log(image) filtered = ndimage.convolve(log_image, kernel) / (kernel_size * kernel_size) filtered = np.exp(filtered) return np.uint8(np.clip(filtered, 0, 255)) def median_filter(image, kernel_size=3): """中值滤波""" return cv2.medianBlur(image, kernel_size) def adaptive_median_filter(image, kernel_size=3, max_kernel_size=7): """自适应中值滤波""" h, w = image.shape pad_size = max_kernel_size // 2 padded = np.pad(image, pad_size, mode='reflect') result = np.copy(image) for i in range(h): for j in range(w): current_size = kernel_size while current_size <= max_kernel_size: # 提取邻域 y1, y2 = i, i + current_size x1, x2 = j, j + current_size neighborhood = padded[y1:y2, x1:x2] # 计算统计量 z_min = np.min(neighborhood) z_max = np.max(neighborhood) z_med = np.median(neighborhood) z_xy = image[i, j] # 阶段A A1 = z_med - z_min A2 = z_med - z_max if A1 > 0 and A2 < 0: # 阶段B B1 = z_xy - z_min B2 = z_xy - z_max if B1 > 0 and B2 < 0: result[i, j] = z_xy else: result[i, j] = z_med break else: current_size += 2 if current_size > max_kernel_size: result[i, j] = z_med break return result # ====================== 3. 频域滤波（周期性噪声抑制） ====================== def create_notch_filter(shape, center, radius=5, notch_type='reject'): """创建陷波滤波器""" rows, cols = shape crow, ccol = center mask = np.ones((rows, cols), np.float32) # 创建圆形掩码 x, y = np.ogrid[:rows, :cols] mask_area = (x - crow)**2 + (y - ccol)**2 <= radius**2 if notch_type == 'reject': mask[mask_area] = 0 else: # pass mask[mask_area] = 1 mask[~mask_area] = 0 # 对称陷波（频域共轭对称） mask[rows-crow, cols-ccol] = mask[crow, ccol] return mask def notch_filtering(image, noise_freqs, radius=5): """陷波滤波消除周期性噪声""" # 转换为浮点型 img_float = np.float32(image) / 255.0 # 傅里叶变换 f = fft2(img_float) f_shift = fftshift(f) # 创建陷波滤波器 notch_filter = np.ones_like(f_shift) for freq in noise_freqs: notch_filter *= create_notch_filter(img_float.shape, freq, radius) # 应用滤波器 f_shift_filtered = f_shift * notch_filter f_ishift = ifftshift(f_shift_filtered) img_filtered = ifft2(f_ishift) img_filtered = np.abs(img_filtered) # 归一化 img_filtered = np.uint8(np.clip(img_filtered * 255, 0, 255)) return img_filtered, f_shift, f_shift_filtered # ====================== 4. 图像复原（逆滤波、维纳滤波） ====================== def motion_blur_kernel(shape, angle=45, len=15): """生成运动模糊核（PSF）""" kernel = np.zeros(shape) angle_rad = np.deg2rad(angle) # 计算模糊核的中心和方向 center = (shape[0]//2, shape[1]//2) x_len = len * np.cos(angle_rad) y_len = len * np.sin(angle_rad) # 生成线条 x = np.linspace(center[0] - x_len/2, center[0] + x_len/2, len) y = np.linspace(center[1] - y_len/2, center[1] + y_len/2, len) x = np.round(x).astype(int) y = np.round(y).astype(int) # 确保坐标在范围内 mask = (x >= 0) & (x < shape[0]) & (y >= 0) & (y < shape[1]) x = x[mask] y = y[mask] kernel[x, y] = 1 kernel = kernel / np.sum(kernel) # 归一化 return kernel def inverse_filtering(blurred_image, psf, epsilon=1e-6): """逆滤波""" # 转换为浮点型 img_float = np.float32(blurred_image) / 255.0 # 傅里叶变换 f = fft2(img_float) psf_padded = np.zeros_like(img_float) psf_padded[:psf.shape[0], :psf.shape[1]] = psf psf_padded = np.roll(psf_padded, (-psf.shape[0]//2, -psf.shape[1]//2), (0, 1)) # 退化函数的傅里叶变换 H = fft2(psf_padded) # 逆滤波（添加epsilon避免除以0） F_hat = f / (H + epsilon) # 逆傅里叶变换 img_restored = ifft2(F_hat) img_restored = np.abs(img_restored) # 归一化 img_restored = np.uint8(np.clip(img_restored * 255, 0, 255)) return img_restored def wiener_filtering(blurred_image, psf, K=0.01): """维纳滤波（简化版，K为噪声功率与信号功率比）""" # 转换为浮点型 img_float = np.float32(blurred_image) / 255.0 # 傅里叶变换 f = fft2(img_float) psf_padded = np.zeros_like(img_float) psf_padded[:psf.shape[0], :psf.shape[1]] = psf psf_padded = np.roll(psf_padded, (-psf.shape[0]//2, -psf.shape[1]//2), (0, 1)) # 退化函数的傅里叶变换 H = fft2(psf_padded) H_conj = np.conj(H) # 维纳滤波 F_hat = (H_conj / (np.abs(H)**2 + K)) * f # 逆傅里叶变换 img_restored = ifft2(F_hat) img_restored = np.abs(img_restored) # 归一化 img_restored = np.uint8(np.clip(img_restored * 255, 0, 255)) return img_restored # ====================== 5. 主函数：演示所有关键算法 ====================== def main(): # 1. 读取图像（转为灰度图） img = cv2.imread('lena.jpg', 0) # 替换为你的图像路径 if img is None: print("无法读取图像，请检查路径！") print("使用内置测试图像...") # 创建测试图像（512x512 Lena替代） img = np.uint8(ndimage.imread('https://upload.wikimedia.org/wikipedia/en/7/7d/Lenna_%28test_image%29.png', mode='L')) # ====================== 实验1：噪声生成与空间域滤波 ====================== print("=== 实验1：噪声生成与空间域滤波 ===") # 添加不同噪声 img_gaussian = add_gaussian_noise(img, var=0.01) img_salt_pepper = add_salt_pepper_noise(img, prob=0.05) img_poisson = add_poisson_noise(img) # 空间域滤波 # 高斯噪声滤波 img_gaussian_mean = arithmetic_mean_filter(img_gaussian, 3) img_gaussian_geo = geometric_mean_filter(img_gaussian, 3) # 椒盐噪声滤波 img_salt_pepper_median = median_filter(img_salt_pepper, 3) img_salt_pepper_adaptive = adaptive_median_filter(img_salt_pepper, 3, 7) # 可视化对比 plt.figure(figsize=(20, 15)) # 原始图像 plt.subplot(3, 4, 1) plt.imshow(img, cmap='gray') plt.title('原始图像', fontsize=12) plt.axis('off') # 高斯噪声及滤波 plt.subplot(3, 4, 2) plt.imshow(img_gaussian, cmap='gray') plt.title('添加高斯噪声', fontsize=12) plt.axis('off') plt.subplot(3, 4, 3) plt.imshow(img_gaussian_mean, cmap='gray') plt.title('算术均值滤波 (3x3)', fontsize=12) plt.axis('off') plt.subplot(3, 4, 4) plt.imshow(img_gaussian_geo, cmap='gray') plt.title('几何均值滤波 (3x3)', fontsize=12) plt.axis('off') # 椒盐噪声及滤波 plt.subplot(3, 4, 5) plt.imshow(img_salt_pepper, cmap='gray') plt.title('添加椒盐噪声', fontsize=12) plt.axis('off') plt.subplot(3, 4, 6) plt.imshow(img_salt_pepper_median, cmap='gray') plt.title('中值滤波 (3x3)', fontsize=12) plt.axis('off') plt.subplot(3, 4, 7) plt.imshow(img_salt_pepper_adaptive, cmap='gray') plt.title('自适应中值滤波', fontsize=12) plt.axis('off') # 泊松噪声 plt.subplot(3, 4, 8) plt.imshow(img_poisson, cmap='gray') plt.title('添加泊松噪声', fontsize=12) plt.axis('off') plt.tight_layout() plt.show() # ====================== 实验2：周期性噪声与陷波滤波 ====================== print("=== 实验2：周期性噪声与陷波滤波 ===") # 生成周期性噪声 rows, cols = img.shape x, y = np.meshgrid(np.arange(cols), np.arange(rows)) # 添加正弦噪声（模拟周期性噪声） noise = 30 * np.sin(2 * np.pi * 10 * x / cols) + 30 * np.sin(2 * np.pi * 10 * y / rows) img_periodic_noise = img + noise img_periodic_noise = np.uint8(np.clip(img_periodic_noise, 0, 255)) # 陷波滤波 # 噪声频率位置（频域中心附近） noise_freqs = [(rows//2, cols//2 + cols//20), (rows//2 + rows//20, cols//2)] img_notch_filtered, f_shift, f_shift_filtered = notch_filtering(img_periodic_noise, noise_freqs) # 可视化 plt.figure(figsize=(18, 6)) plt.subplot(1, 3, 1) plt.imshow(img_periodic_noise, cmap='gray') plt.title('含周期性噪声的图像', fontsize=12) plt.axis('off') plt.subplot(1, 3, 2) plt.imshow(20 * np.log1p(np.abs(f_shift)), cmap='gray') plt.title('噪声图像的频域', fontsize=12) plt.axis('off') plt.subplot(1, 3, 3) plt.imshow(img_notch_filtered, cmap='gray') plt.title('陷波滤波复原图像', fontsize=12) plt.axis('off') plt.tight_layout() plt.show() # ====================== 实验3：运动模糊与图像复原 ====================== print("=== 实验3：运动模糊与图像复原 ===") # 生成运动模糊核 psf = motion_blur_kernel((15, 15), angle=45, len=15) # 添加运动模糊 img_blurred = ndimage.convolve(img, psf, mode='reflect') # 添加少量高斯噪声 img_blurred_noisy = add_gaussian_noise(img_blurred, var=0.001) # 图像复原 img_inverse = inverse_filtering(img_blurred_noisy, psf) img_wiener = wiener_filtering(img_blurred_noisy, psf, K=0.01) # 可视化 plt.figure(figsize=(20, 8)) plt.subplot(1, 5, 1) plt.imshow(img, cmap='gray') plt.title('原始图像', fontsize=12) plt.axis('off') plt.subplot(1, 5, 2) plt.imshow(img_blurred, cmap='gray') plt.title('运动模糊图像', fontsize=12) plt.axis('off') plt.subplot(1, 5, 3) plt.imshow(img_blurred_noisy, cmap='gray') plt.title('模糊+噪声图像', fontsize=12) plt.axis('off') plt.subplot(1, 5, 4) plt.imshow(img_inverse, cmap='gray') plt.title('逆滤波复原', fontsize=12) plt.axis('off') plt.subplot(1, 5, 5) plt.imshow(img_wiener, cmap='gray') plt.title('维纳滤波复原', fontsize=12) plt.axis('off') plt.tight_layout() plt.show() # ====================== 运行主函数 ====================== if __name__ == "__main__": # 安装依赖（如果需要） # !pip install opencv-python numpy matplotlib scipy main()

四、代码运行说明

4.1 环境要求

• Python 3.7+
• 依赖库：

  pip install opencv-python numpy matplotlib scipy

4.2 运行步骤

1. 将代码保存为image_restoration.py
2. 准备测试图像（如 lena.jpg），或使用代码中的内置图像
3. 运行代码：python image_restoration.py
4. 查看自动弹出的效果对比图

4.3 关键参数调整

• 噪声强度：调整add_gaussian_noise中的var参数、add_salt_pepper_noise中的prob参数
• 滤波器大小：调整均值 / 中值滤波的kernel_size参数
• 维纳滤波参数：调整wiener_filtering中的K参数（噪声功率比）

五、小结

核心知识点回顾

1. 图像退化模型：核心公式 g(x,y)=H[f(x,y)]+η(x,y)，退化函数 + 噪声共同导致图像质量下降
2. 噪声处理：空间域滤波（均值、中值、自适应滤波）适合随机噪声，频域滤波（陷波滤波）适合周期性噪声
3. 图像复原：逆滤波简单但对噪声敏感，维纳滤波（最小均方误差）是更鲁棒的选择，通过权衡退化逆过程和噪声抑制实现最优复原
4. 图像重建：基于投影和拉东变换，是 CT 等医学成像的核心技术，滤波反投影是最常用的重建算法

工程应用要点

1. 实际应用中需先分析图像退化类型（噪声 / 模糊 / 混合），选择合适的复原算法
2. 参数调优是关键（如滤波器大小、维纳滤波的 K 值），需结合具体场景调整
3. 自适应滤波器在保留细节和去除噪声之间取得了更好的平衡，是实际项目中的首选

📌 本文所有代码均可直接运行，建议你修改参数、更换测试图像，亲身体验不同算法的效果差异，加深对图像复原技术的理解