企业官网建设流程全解析

目录

3.3 均值滤波（续）

频域理解与信号处理视角

边界效应分析

计算复杂度优化

变体：加权均值滤波

均值滤波完整Python程序

3.4 高斯滤波

核心思想与直观理解

数学原理：二维高斯函数

高斯核的离散化与创建

高斯滤波的性质与优势

σ参数的影响

与均值滤波的对比

高斯滤波完整Python程序

3.3 均值滤波（续）

频域理解与信号处理视角

从信号处理的角度深入理解均值滤波，我们需要引入傅里叶变换和频率域的概念。

空域与频域：

• 空域：我们通常看到的图像是空域表示，像素值随空间位置变化。

• 频域：将图像视为不同频率的波的叠加。低频对应大面积平滑区域，高频对应细节、边缘和噪声。

均值滤波的频域响应：
均值滤波核的傅里叶变换是一个sinc函数。在频域中，它像一个低通滤波器：

• 对于低频分量（变化缓慢的区域），sinc函数值接近1，几乎完全通过。

• 对于高频分量（变化剧烈的区域），sinc函数值衰减很快，被显著抑制。

数学上，N×N均值滤波器的频率响应为：

H(u,v) = (1/N²) * sin(πuN)/sin(πu) * sin(πvN)/sin(πv)

其中u,v是归一化频率。

边界效应分析

当滤波器核移动到图像边界时，需要特殊处理。常见的边界处理策略有：

1. 零填充（Zero-padding）：在图像外围填充0。

   • 优点：简单

   • 缺点：在边界处引入暗边效应

2. 复制填充（Replicate）：复制最近的边界像素值。

   • 优点：避免暗边

   • 缺点：可能造成边界模糊

3. 镜像填充（Reflect）：镜像边界外的像素。

   • 优点：保持连续性

   • 缺点：计算稍复杂

4. 循环填充（Wrap）：假设图像是周期性的。

   • 优点：数学上优美

   • 缺点：实际图像通常不是周期性的

计算复杂度优化

对于大尺寸均值滤波，直接卷积的计算复杂度是O(m×n×M×N)，其中m×n是核大小，M×N是图像大小。我们可以使用积分图像或滑动窗口累加技术将复杂度降至O(M×N)。

积分图像算法：

1. 计算积分图像I_Σ：

   I_Σ(x,y) = I(x,y) + I_Σ(x-1,y) + I_Σ(x,y-1) - I_Σ(x-1,y-1)

2. 任意矩形区域的像素和可以通过积分图像的4个值快速计算：

   sum = I_Σ(x2,y2) - I_Σ(x1-1,y2) - I_Σ(x2,y1-1) + I_Σ(x1-1,y1-1)

3. 均值 = sum / (矩形面积)

变体：加权均值滤波

标准均值滤波对所有邻居一视同仁。但直觉上，离中心越近的像素应该越相关。因此出现了加权均值滤波，给予不同位置不同的权重。

常用的加权方案：

1. 距离加权：权重与到中心的距离成反比

2. 高斯加权：权重服从高斯分布（这过渡到了高斯滤波）

均值滤波完整Python程序

import cv2 import numpy as np import matplotlib.pyplot as plt from scipy import ndimage import time def mean_filter_custom(image, kernel_size=3, padding_mode='reflect'): """ 自定义均值滤波函数 :param image: 输入灰度图像 :param kernel_size: 滤波器大小，必须是奇数 :param padding_mode: 边界填充模式 :return: 滤波后的图像 """ if kernel_size % 2 == 0: raise ValueError("kernel_size必须是奇数") # 创建均值核 kernel = np.ones((kernel_size, kernel_size), dtype=np.float32) / (kernel_size * kernel_size) # 获取填充大小 pad = kernel_size // 2 # 根据填充模式处理图像 if padding_mode == 'zero': padded = np.pad(image.astype(np.float32), pad, mode='constant', constant_values=0) elif padding_mode == 'replicate': padded = np.pad(image.astype(np.float32), pad, mode='edge') elif padding_mode == 'reflect': padded = np.pad(image.astype(np.float32), pad, mode='reflect') elif padding_mode == 'wrap': padded = np.pad(image.astype(np.float32), pad, mode='wrap') else: raise ValueError(f"不支持的填充模式: {padding_mode}") # 初始化输出图像 output = np.zeros_like(image, dtype=np.float32) # 手动实现卷积（为了教学目的） height, width = image.shape for i in range(pad, height + pad): for j in range(pad, width + pad): # 提取邻域 neighborhood = padded[i-pad:i+pad+1, j-pad:j+pad+1] # 计算加权和 output[i-pad, j-pad] = np.sum(neighborhood * kernel) # 裁剪到[0, 255]范围并转换为uint8 output = np.clip(output, 0, 255).astype(np.uint8) return output def mean_filter_integral(image, kernel_size=3): """ 使用积分图像优化的均值滤波 :param image: 输入灰度图像 :param kernel_size: 滤波器大小，必须是奇数 :return: 滤波后的图像 """ if kernel_size % 2 == 0: raise ValueError("kernel_size必须是奇数") # 转换为float32以便计算 img_float = image.astype(np.float32) height, width = img_float.shape pad = kernel_size // 2 # 1. 计算积分图像 integral = np.zeros((height+1, width+1), dtype=np.float32) # 计算积分图像 for i in range(1, height+1): row_sum = 0 for j in range(1, width+1): row_sum += img_float[i-1, j-1] integral[i, j] = integral[i-1, j] + row_sum # 2. 初始化输出图像 output = np.zeros_like(img_float) # 3. 使用积分图像快速计算区域和 for i in range(height): for j in range(width): # 计算矩形边界（考虑边界情况） x1 = max(0, i - pad) y1 = max(0, j - pad) x2 = min(height-1, i + pad) y2 = min(width-1, j + pad) # 使用积分图像快速计算区域和 # 注意：积分图像的索引比原图像大1 area_sum = (integral[x2+1, y2+1] - integral[x1, y2+1] - integral[x2+1, y1] + integral[x1, y1]) # 计算实际区域面积 area = (x2 - x1 + 1) * (y2 - y1 + 1) # 计算平均值 output[i, j] = area_sum / area # 裁剪到[0, 255]范围并转换为uint8 output = np.clip(output, 0, 255).astype(np.uint8) return output def compare_performance(image, kernel_size=15): """ 比较不同实现方式的性能 """ print(f"\n性能测试 - 核大小: {kernel_size}×{kernel_size}") print("-" * 50) # 测试自定义实现 start_time = time.time() result_custom = mean_filter_custom(image, kernel_size, 'reflect') custom_time = time.time() - start_time print(f"自定义卷积实现: {custom_time:.4f} 秒") # 测试积分图像实现 start_time = time.time() result_integral = mean_filter_integral(image, kernel_size) integral_time = time.time() - start_time print(f"积分图像实现: {integral_time:.4f} 秒") # 测试OpenCV实现 start_time = time.time() result_opencv = cv2.blur(image, (kernel_size, kernel_size)) opencv_time = time.time() - start_time print(f"OpenCV blur函数: {opencv_time:.4f} 秒") # 测试SciPy实现 start_time = time.time() kernel = np.ones((kernel_size, kernel_size)) / (kernel_size * kernel_size) result_scipy = ndimage.convolve(image.astype(float), kernel, mode='reflect') result_scipy = np.clip(result_scipy, 0, 255).astype(np.uint8) scipy_time = time.time() - start_time print(f"SciPy convolve函数: {scipy_time:.4f} 秒") return { 'custom': result_custom, 'integral': result_integral, 'opencv': result_opencv, 'scipy': result_scipy, 'times': { 'custom': custom_time, 'integral': integral_time, 'opencv': opencv_time, 'scipy': scipy_time } } def analyze_edge_preservation(original, filtered): """ 分析滤波器的边缘保持能力 """ # 计算梯度幅度（使用Sobel算子） sobel_x = cv2.Sobel(original, cv2.CV_64F, 1, 0, ksize=3) sobel_y = cv2.Sobel(original, cv2.CV_64F, 0, 1, ksize=3) gradient_orig = np.sqrt(sobel_x**2 + sobel_y**2) sobel_x_f = cv2.Sobel(filtered, cv2.CV_64F, 1, 0, ksize=3) sobel_y_f = cv2.Sobel(filtered, cv2.CV_64F, 0, 1, ksize=3) gradient_filt = np.sqrt(sobel_x_f**2 + sobel_y_f**2) # 计算梯度保留率 mask = gradient_orig > 10 # 只考虑明显的边缘 if np.sum(mask) > 0: retention_ratio = np.mean(gradient_filt[mask] / gradient_orig[mask]) else: retention_ratio = 0 return gradient_orig, gradient_filt, retention_ratio # 主程序 if __name__ == "__main__": # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False # 1. 读取或创建测试图像 img = cv2.imread('test_image.jpg', cv2.IMREAD_GRAYSCALE) if img is None: print("未找到test_image.jpg，创建合成测试图像...") # 创建包含边缘、纹理和平滑区域的合成图像 img = np.zeros((256, 256), dtype=np.uint8) # 添加垂直边缘 img[:, 100:110] = 255 # 添加斜边缘 for i in range(256): for j in range(256): if abs(i - j) < 5: img[i, j] = 128 # 添加纹理区域 for i in range(50, 150): for j in range(150, 250): img[i, j] = 100 + 50 * np.sin(i/10) * np.cos(j/10) # 添加平滑渐变区域 for i in range(200, 256): img[i, :] = i - 200 # 2. 添加混合噪声 np.random.seed(42) # 添加高斯噪声 gaussian_noise = np.random.normal(0, 25, img.shape) img_gaussian = img.astype(float) + gaussian_noise img_gaussian = np.clip(img_gaussian, 0, 255).astype(np.uint8) # 添加椒盐噪声 img_salt_pepper = img.copy() salt_pepper_mask = np.random.random(img.shape) img_salt_pepper[salt_pepper_mask < 0.01] = 0 img_salt_pepper[salt_pepper_mask > 0.99] = 255 # 3. 应用不同大小的均值滤波器 kernel_sizes = [3, 7, 15, 31] # 准备可视化 fig, axes = plt.subplots(3, len(kernel_sizes)+1, figsize=(20, 15)) # 显示原始图像和噪声图像 axes[0, 0].imshow(img, cmap='gray', vmin=0, vmax=255) axes[0, 0].set_title('原始图像') axes[0, 0].axis('off') axes[1, 0].imshow(img_gaussian, cmap='gray', vmin=0, vmax=255) axes[1, 0].set_title('高斯噪声图像\n(σ=25)') axes[1, 0].axis('off') axes[2, 0].imshow(img_salt_pepper, cmap='gray', vmin=0, vmax=255) axes[2, 0].set_title('椒盐噪声图像\n(密度=2%)') axes[2, 0].axis('off') # 对每个核大小进行处理 for idx, ksize in enumerate(kernel_sizes): # 处理高斯噪声图像 filtered_gaussian = mean_filter_custom(img_gaussian, ksize, 'reflect') axes[0, idx+1].imshow(filtered_gaussian, cmap='gray', vmin=0, vmax=255) axes[0, idx+1].set_title(f'高斯噪声+均值滤波\n核大小={ksize}×{ksize}') axes[0, idx+1].axis('off') # 处理椒盐噪声图像 filtered_sp = mean_filter_custom(img_salt_pepper, ksize, 'reflect') axes[1, idx+1].imshow(filtered_sp, cmap='gray', vmin=0, vmax=255) axes[1, idx+1].set_title(f'椒盐噪声+均值滤波\n核大小={ksize}×{ksize}') axes[1, idx+1].axis('off') # 分析边缘保留 grad_orig, grad_filt, retention = analyze_edge_preservation(img, filtered_gaussian) axes[2, idx+1].imshow(grad_filt, cmap='hot') axes[2, idx+1].set_title(f'梯度图\n边缘保留率: {retention:.2%}') axes[2, idx+1].axis('off') plt.tight_layout() plt.savefig('mean_filter_comparison.png', dpi=300, bbox_inches='tight') # 4. 性能比较 print("\n" + "="*60) print("均值滤波性能与效果分析") print("="*60) # 测试不同实现方式的性能 perf_results = compare_performance(img, kernel_size=15) # 5. 分析不同核大小对噪声抑制和边缘模糊的影响 print("\n" + "-"*60) print("核大小对滤波效果的影响分析") print("-"*60) for ksize in [3, 5, 9, 15, 25]: filtered = mean_filter_custom(img_gaussian, ksize, 'reflect') # 计算噪声抑制效果 noise_original = np.std(img_gaussian - img) noise_filtered = np.std(filtered - img) noise_reduction = (1 - noise_filtered/noise_original) * 100 # 计算边缘模糊程度 _, _, retention = analyze_edge_preservation(img, filtered) print(f"核大小 {ksize:2d}×{ksize:<2d} | " f"噪声抑制: {noise_reduction:6.2f}% | " f"边缘保留: {retention:.2%}") # 6. 可视化边界处理效果 print("\n" + "-"*60) print("不同边界处理策略比较") print("-"*60) padding_modes = ['zero', 'replicate', 'reflect', 'wrap'] fig2, axes2 = plt.subplots(2, 4, figsize=(16, 8)) # 创建一个小图像以突出边界效果 small_img = np.zeros((50, 50), dtype=np.uint8) small_img[10:40, 10:40] = 200 small_img[20:30, 20:30] = 50 for i, mode in enumerate(padding_modes): filtered_small = mean_filter_custom(small_img, 15, mode) axes2[0, i].imshow(filtered_small, cmap='gray') axes2[0, i].set_title(f'边界处理: {mode}') axes2[0, i].axis('off') # 显示中心行剖面 center_row = filtered_small[25, :] axes2[1, i].plot(center_row, 'b-', linewidth=2) axes2[1, i].set_title(f'中心行灰度剖面') axes2[1, i].set_xlabel('列坐标') axes2[1, i].set_ylabel('灰度值') axes2[1, i].grid(True, alpha=0.3) axes2[1, i].set_ylim([0, 255]) plt.tight_layout() plt.savefig('boundary_effects.png', dpi=300, bbox_inches='tight') # 7. 显示所有结果 plt.show() print("\n" + "="*60) print("实验总结") print("="*60) print("1. 核大小的影响:") print(" - 小核(3×3): 轻微降噪，边缘保留好，但椒盐噪声去除不彻底") print(" - 大核(15×15以上): 强降噪，但边缘严重模糊") print("") print("2. 边界处理策略:") print(" - zero: 产生暗边，不推荐") print(" - replicate: 常用，避免暗边但可能模糊边界") print(" - reflect: 推荐，保持连续性") print(" - wrap: 仅适用于周期性图像") print("") print("3. 均值滤波的局限性:") print(" - 对高斯噪声: 有效但边缘会模糊") print(" - 对椒盐噪声: 效果差，噪声被扩散而非去除") print(" - 计算效率: 直接实现慢，积分图像法对大核更高效") print("") print("4. 实用建议:") print(" - 轻度降噪: 使用3×3或5×5小核") print(" - 优先选择reflect边界处理") print(" - 对于椒盐噪声，考虑使用中值滤波而非均值滤波")

3.4 高斯滤波

核心思想与直观理解

如果说均值滤波是对邻域内所有像素"一视同仁"，那么高斯滤波则是"区别对待"。它的核心思想是：离中心像素越近的像素，对中心像素的影响应该越大；离得越远，影响应该越小。这种思想符合我们对自然图像的认知——相邻像素通常高度相关，而距离越远相关性越弱。

从概率论角度看，高斯滤波假设像素值在空间上的分布服从高斯分布（正态分布）。这种滤波器以其发明者卡尔·弗里德里希·高斯命名，是计算机视觉中最重要、最常用的线性滤波器之一。

数学原理：二维高斯函数

一维高斯函数的公式为：

G(x) = (1/(σ√(2π))) * exp(-x²/(2σ²))

扩展到二维，得到二维高斯函数：

G(x,y) = (1/(2πσ²)) * exp(-(x²+y²)/(2σ²))

其中：

• (x,y)是到中心点的距离

• σ是标准差，控制分布的宽度

• 系数1/(2πσ²)是归一化常数，确保函数积分为1

高斯核的离散化与创建

在实际应用中，我们需要将连续的高斯函数离散化为一个有限大小的卷积核。创建高斯核的步骤：

1. 确定核大小：通常根据σ选择，经验法则是核大小 ≈ 6σ+1（确保覆盖99.7%的能量）

2. 计算离散坐标：对于核中的每个位置(i,j)，计算其到中心的距离

3. 应用高斯公式：计算每个位置的权重

4. 归一化：确保所有权重之和为1

例如，当σ=1.0时，一个5×5高斯核近似为：

[[0.003, 0.013, 0.022, 0.013, 0.003],
[0.013, 0.059, 0.097, 0.059, 0.013],
[0.022, 0.097, 0.159, 0.097, 0.022],
[0.013, 0.059, 0.097, 0.059, 0.013],
[0.003, 0.013, 0.022, 0.013, 0.003]]

高斯滤波的性质与优势

1. 旋转对称性：二维高斯函数是旋转对称的，这意味着高斯滤波在各个方向上的平滑效果相同。

2. 可分离性：这是高斯滤波最重要的性质之一。二维高斯函数可以分解为两个一维高斯函数的乘积：

G(x,y) = G(x) * G(y) = (1/(√(2π)σ))exp(-x²/(2σ²)) * (1/(√(2π)σ))exp(-y²/(2σ²))

1. 这意味着二维高斯滤波可以分解为两次一维高斯滤波（先水平后垂直），将计算复杂度从O(m×n)降低到O(m+n)。

2. 傅里叶变换性质：高斯函数的傅里叶变换仍是高斯函数。在频域中，高斯滤波器是理想的低通滤波器，没有振铃效应。

3. 尺度空间理论：高斯滤波是构建图像尺度空间的基础，这在SIFT等特征检测算法中至关重要。

σ参数的影响

标准差σ是高斯滤波的关键参数：

• 小σ（如0.5-1.0）：高斯核尖锐，主要平滑高频噪声，保留大部分细节

• 中σ（如1.5-2.5）：平衡噪声抑制和细节保留

• 大σ（如3.0以上）：高斯核宽平，强平滑效果，大量细节丢失

与均值滤波的对比

高斯滤波完整Python程序

import cv2 import numpy as np import matplotlib.pyplot as plt from scipy import ndimage import time from mpl_toolkits.mplot3d import Axes3D def gaussian_kernel_2d(size, sigma): """ 生成二维高斯卷积核 :param size: 核大小（奇数） :param sigma: 标准差 :return: 二维高斯核 """ if size % 2 == 0: raise ValueError("核大小必须是奇数") # 创建坐标网格 ax = np.arange(-size // 2 + 1., size // 2 + 1.) xx, yy = np.meshgrid(ax, ax) # 计算二维高斯函数 kernel = np.exp(-(xx**2 + yy**2) / (2. * sigma**2)) # 归一化 kernel /= np.sum(kernel) return kernel def gaussian_kernel_1d(size, sigma): """ 生成一维高斯卷积核（用于可分离滤波） :param size: 核大小（奇数） :param sigma: 标准差 :return: 一维高斯核 """ if size % 2 == 0: raise ValueError("核大小必须是奇数") ax = np.arange(-size // 2 + 1., size // 2 + 1.) kernel = np.exp(-ax**2 / (2. * sigma**2)) kernel /= np.sum(kernel) return kernel def gaussian_filter_separable(image, size, sigma, padding_mode='reflect'): """ 可分离高斯滤波实现（先水平后垂直） :param image: 输入图像 :param size: 核大小 :param sigma: 标准差 :param padding_mode: 边界处理模式 :return: 滤波后的图像 """ # 生成一维高斯核 kernel_1d = gaussian_kernel_1d(size, sigma) # 转换为float32以便计算 img_float = image.astype(np.float32) # 水平方向卷积 if padding_mode == 'zero': pad_width = size // 2 padded = np.pad(img_float, ((0, 0), (pad_width, pad_width)), mode='constant') temp = np.zeros_like(img_float) for i in range(img_float.shape[0]): for j in range(img_float.shape[1]): temp[i, j] = np.sum(padded[i, j:j+size] * kernel_1d) else: # 使用SciPy的convolve1d（支持多种边界模式） temp = ndimage.convolve1d(img_float, kernel_1d, axis=1, mode=padding_mode) # 垂直方向卷积 if padding_mode == 'zero': padded = np.pad(temp, ((pad_width, pad_width), (0, 0)), mode='constant') output = np.zeros_like(img_float) for i in range(img_float.shape[0]): for j in range(img_float.shape[1]): output[i, j] = np.sum(padded[i:i+size, j] * kernel_1d) else: output = ndimage.convolve1d(temp, kernel_1d, axis=0, mode=padding_mode) # 裁剪并转换类型 output = np.clip(output, 0, 255).astype(np.uint8) return output def analyze_frequency_response(sigma_values, kernel_size=31): """ 分析不同σ值的高斯滤波器频率响应 """ fig = plt.figure(figsize=(15, 10)) for idx, sigma in enumerate(sigma_values): # 生成高斯核 kernel = gaussian_kernel_2d(kernel_size, sigma) # 计算频率响应（2D FFT） kernel_padded = np.pad(kernel, ((100, 100), (100, 100)), mode='constant') fft_kernel = np.fft.fft2(kernel_padded) fft_shifted = np.fft.fftshift(fft_kernel) magnitude = 20 * np.log(np.abs(fft_shifted) + 1e-10) # 绘制核 ax1 = plt.subplot(3, len(sigma_values), idx+1) im1 = ax1.imshow(kernel, cmap='hot', interpolation='nearest') ax1.set_title(f'σ={sigma}\n高斯核') ax1.axis('off') plt.colorbar(im1, ax=ax1, fraction=0.046, pad=0.04) # 绘制3D表面图 ax2 = plt.subplot(3, len(sigma_values), idx+1+len(sigma_values), projection='3d') x = np.arange(kernel_size) - kernel_size//2 y = np.arange(kernel_size) - kernel_size//2 X, Y = np.meshgrid(x, y) ax2.plot_surface(X, Y, kernel, cmap='viridis', alpha=0.8) ax2.set_title(f'σ={sigma}的3D视图') ax2.set_xlabel('X') ax2.set_ylabel('Y') ax2.set_zlabel('权重') # 绘制频率响应 ax3 = plt.subplot(3, len(sigma_values), idx+1+2*len(sigma_values)) center = magnitude.shape[0] // 2 profile = magnitude[center, center-50:center+50] ax3.plot(profile, 'b-', linewidth=2) ax3.set_title(f'频率响应剖面') ax3.set_xlabel('空间频率') ax3.set_ylabel('幅度(dB)') ax3.grid(True, alpha=0.3) plt.tight_layout() return fig def compare_gaussian_mean(image, sigma=2.0, kernel_size=15): """ 比较高斯滤波和均值滤波的效果 """ # 添加高斯噪声 noisy = image.astype(float) + np.random.normal(0, 25, image.shape) noisy = np.clip(noisy, 0, 255).astype(np.uint8) # 应用高斯滤波 gaussian_filtered = gaussian_filter_separable(noisy, kernel_size, sigma, 'reflect') # 应用均值滤波（相同核大小） mean_filtered = cv2.blur(noisy, (kernel_size, kernel_size)) # 计算性能 start_time = time.time() _ = gaussian_filter_separable(noisy, kernel_size, sigma, 'reflect') gaussian_time = time.time() - start_time start_time = time.time() _ = cv2.blur(noisy, (kernel_size, kernel_size)) mean_time = time.time() - start_time # 计算质量指标 def calculate_metrics(original, filtered): # PSNR (峰值信噪比) mse = np.mean((original.astype(float) - filtered.astype(float)) ** 2) if mse == 0: return float('inf'), 0 psnr = 20 * np.log10(255.0 / np.sqrt(mse)) # SSIM (结构相似性) 的简化版本 C1 = (0.01 * 255) ** 2 C2 = (0.03 * 255) ** 2 mu_x = np.mean(original) mu_y = np.mean(filtered) sigma_x = np.std(original) sigma_y = np.std(filtered) sigma_xy = np.mean((original - mu_x) * (filtered - mu_y)) ssim = ((2 * mu_x * mu_y + C1) * (2 * sigma_xy + C2)) / \ ((mu_x**2 + mu_y**2 + C1) * (sigma_x**2 + sigma_y**2 + C2)) return psnr, ssim gaussian_psnr, gaussian_ssim = calculate_metrics(image, gaussian_filtered) mean_psnr, mean_ssim = calculate_metrics(image, mean_filtered) return { 'noisy': noisy, 'gaussian': gaussian_filtered, 'mean': mean_filtered, 'times': {'gaussian': gaussian_time, 'mean': mean_time}, 'metrics': { 'gaussian': {'psnr': gaussian_psnr, 'ssim': gaussian_ssim}, 'mean': {'psnr': mean_psnr, 'ssim': mean_ssim} } } def multi_scale_gaussian_pyramid(image, num_levels=4, sigma=1.6): """ 构建高斯金字塔（尺度空间） """ pyramid = [image] current = image.astype(float) for level in range(1, num_levels): # 应用高斯滤波 filtered = gaussian_filter_separable(current.astype(np.uint8), 5, sigma, 'reflect') # 下采样 downsampled = filtered[::2, ::2] pyramid.append(downsampled.astype(np.uint8)) current = downsampled return pyramid # 主程序 if __name__ == "__main__": # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False print("="*70) print("高斯滤波深入分析与实验") print("="*70) # 1. 创建或读取测试图像 print("\n1. 准备测试图像...") img = cv2.imread('test_image.jpg', cv2.IMREAD_GRAYSCALE) if img is None: print("创建合成测试图像...") # 创建包含多种特征的测试图像 img = np.zeros((300, 300), dtype=np.uint8) # 添加不同频率的内容 for i in range(300): for j in range(300): # 低频区域（平滑渐变） if i < 100 and j < 100: img[i, j] = i + j # 中频区域（正弦纹理） elif i < 200 and j < 200: img[i, j] = 128 + 50 * np.sin(i/20) * np.cos(j/20) # 高频区域（棋盘格） else: img[i, j] = 255 if ((i//20) + (j//20)) % 2 == 0 else 0 # 添加边缘 img[150:160, :] = 255 # 水平边缘 img[:, 150:160] = 255 # 垂直边缘 # 2. 分析不同σ值的影响 print("\n2. 分析不同σ值的高斯核特性...") sigma_values = [0.5, 1.0, 2.0, 4.0] fig1 = analyze_frequency_response(sigma_values) plt.savefig('gaussian_frequency_response.png', dpi=300, bbox_inches='tight') # 3. 可分离性验证 print("\n3. 验证高斯滤波的可分离性...") sigma = 2.0 kernel_size = 15 # 生成二维核 kernel_2d = gaussian_kernel_2d(kernel_size, sigma) # 生成一维核 kernel_1d = gaussian_kernel_1d(kernel_size, sigma) # 验证可分离性：kernel_2d ≈ kernel_1d * kernel_1d.T kernel_separable = np.outer(kernel_1d, kernel_1d) separation_error = np.mean(np.abs(kernel_2d - kernel_separable)) print(f" 可分离性验证误差: {separation_error:.6e}") print(f" 二维核计算量: O({kernel_size}²) = {kernel_size**2} 次乘法/像素") print(f" 可分离计算量: O(2×{kernel_size}) = {2*kernel_size} 次乘法/像素") print(f" 加速比: {(kernel_size**2)/(2*kernel_size):.1f}倍") # 4. 性能测试：可分离 vs 不可分离 print("\n4. 性能测试：可分离实现 vs OpenCV实现...") # 添加噪声 noisy_img = img.astype(float) + np.random.normal(0, 30, img.shape) noisy_img = np.clip(noisy_img, 0, 255).astype(np.uint8) # 测试不同实现 implementations = [ ("OpenCV GaussianBlur", lambda: cv2.GaussianBlur(noisy_img, (kernel_size, kernel_size), sigma)), ("可分离实现（自定义）", lambda: gaussian_filter_separable(noisy_img, kernel_size, sigma, 'reflect')), ("SciPy gaussian_filter", lambda: ndimage.gaussian_filter(noisy_img.astype(float), sigma)) ] results = {} for name, func in implementations: start_time = time.perf_counter() result = func() elapsed = time.perf_counter() - start_time results[name] = {'time': elapsed, 'result': result} print(f" {name:30s}: {elapsed:.4f} 秒") # 5. 与均值滤波的对比 print("\n5. 高斯滤波 vs 均值滤波对比分析...") comparison = compare_gaussian_mean(img, sigma=2.0, kernel_size=15) print(f"\n 质量指标对比:") print(f" 高斯滤波 - PSNR: {comparison['metrics']['gaussian']['psnr']:.2f} dB, " f"SSIM: {comparison['metrics']['gaussian']['ssim']:.4f}") print(f" 均值滤波 - PSNR: {comparison['metrics']['mean']['psnr']:.2f} dB, " f"SSIM: {comparison['metrics']['mean']['ssim']:.4f}") print(f"\n 计算时间对比:") print(f" 高斯滤波: {comparison['times']['gaussian']:.4f} 秒") print(f" 均值滤波: {comparison['times']['mean']:.4f} 秒") # 6. 构建高斯金字塔（尺度空间） print("\n6. 构建高斯金字塔（尺度空间）...") pyramid = multi_scale_gaussian_pyramid(img, num_levels=5, sigma=1.6) # 7. 可视化所有结果 print("\n7. 生成可视化结果...") # 创建综合可视化图 fig2, axes2 = plt.subplots(3, 4, figsize=(16, 12)) # 原始图像和噪声图像 axes2[0, 0].imshow(img, cmap='gray') axes2[0, 0].set_title('原始图像') axes2[0, 0].axis('off') axes2[0, 1].imshow(comparison['noisy'], cmap='gray') axes2[0, 1].set_title('加噪图像 (σ=30)') axes2[0, 1].axis('off') # 不同σ值的滤波结果 for idx, sigma in enumerate([0.5, 1.0, 2.0, 4.0]): row = idx // 2 + 1 col = idx % 2 filtered = gaussian_filter_separable(comparison['noisy'], 15, sigma, 'reflect') axes2[row, col*2].imshow(filtered, cmap='gray') axes2[row, col*2].set_title(f'高斯滤波 σ={sigma}') axes2[row, col*2].axis('off') # 显示剖面线 center_line = filtered[150, :] axes2[row, col*2+1].plot(center_line, 'b-', linewidth=2) axes2[row, col*2+1].set_title(f'中心行剖面 (σ={sigma})') axes2[row, col*2+1].set_xlabel('列坐标') axes2[row, col*2+1].set_ylabel('灰度值') axes2[row, col*2+1].grid(True, alpha=0.3) axes2[row, col*2+1].set_ylim([0, 255]) plt.tight_layout() plt.savefig('gaussian_filter_results.png', dpi=300, bbox_inches='tight') # 8. 可视化高斯金字塔 fig3, axes3 = plt.subplots(1, len(pyramid), figsize=(15, 4)) for i, level_img in enumerate(pyramid): axes3[i].imshow(level_img, cmap='gray') axes3[i].set_title(f'金字塔层级 {i}\n大小: {level_img.shape}') axes3[i].axis('off') plt.tight_layout() plt.savefig('gaussian_pyramid.png', dpi=300, bbox_inches='tight') # 9. 实际应用示例：文档图像预处理 print("\n8. 实际应用：文档图像预处理示例...") # 模拟文档图像（文字+噪声） doc_img = np.ones((200, 300), dtype=np.uint8) * 200 # 添加文字（高频信息） cv2.putText(doc_img, 'Computer Vision', (30, 80), cv2.FONT_HERSHEY_SIMPLEX, 1.5, 50, 3) cv2.putText(doc_img, 'Gaussian Filter Demo', (30, 140), cv2.FONT_HERSHEY_SIMPLEX, 1.2, 30, 2) # 添加噪声 doc_noisy = doc_img.astype(float) + np.random.normal(0, 40, doc_img.shape) doc_noisy = np.clip(doc_noisy, 0, 255).astype(np.uint8) # 应用不同滤波 doc_gaussian = gaussian_filter_separable(doc_noisy, 5, 1.0, 'reflect') doc_mean = cv2.blur(doc_noisy, (5, 5)) # 可视化文档处理结果 fig4, axes4 = plt.subplots(1, 4, figsize=(16, 5)) images = [doc_img, doc_noisy, doc_gaussian, doc_mean] titles = ['原始文档', '加噪文档', '高斯滤波(σ=1.0)', '均值滤波(5×5)'] for i, (ax, img, title) in enumerate(zip(axes4, images, titles)): ax.imshow(img, cmap='gray') ax.set_title(title) ax.axis('off') plt.tight_layout() plt.savefig('document_processing.png', dpi=300, bbox_inches='tight') # 10. 显示所有图表 plt.show() print("\n" + "="*70) print("实验总结与关键发现") print("="*70) print("\n1. 高斯核的性质:") print(" - 可分离性：二维高斯滤波可分解为两次一维滤波，大幅降低计算复杂度") print(" - 旋转对称性：在各方向平滑效果一致，适合处理各向同性噪声") print(" - 尺度空间：通过不同σ构建图像多尺度表示，是特征检测的基础") print("\n2. σ参数的影响:") print(" - σ小(0.5-1.0)：轻微平滑，保留细节，适合预处理") print(" - σ中(1.0-2.0)：平衡噪声抑制和细节保留，通用选择") print(" - σ大(>2.0)：强平滑，损失细节，用于显著降噪") print("\n3. 与均值滤波对比:") print(" - 边缘保持：高斯滤波明显优于均值滤波（权重分布合理）") print(" - 计算效率：可分离实现使高斯滤波效率接近均值滤波") print(" - 理论性质：高斯滤波有完美的数学性质（可分离、傅里叶变换等）") print("\n4. 实际应用建议:") print(" - 通用降噪：σ=1.0-2.0，核大小=3×3到7×7") print(" - 预处理：轻度高斯滤波(σ=0.5-1.0)可减少后续算法对噪声的敏感度") print(" - 尺度空间：使用多个σ值构建图像金字塔，用于多尺度特征检测") print(" - 实时应用：务必使用可分离实现或硬件优化的库函数") print("\n5. 局限性与注意事项:") print(" - 高斯滤波是线性滤波器，对脉冲噪声（椒盐）效果有限") print(" - 过度平滑会损失重要细节和纹理信息") print(" - 边界处理需要谨慎选择，reflect模式通常最安全") print("\n" + "="*70) print("关键公式回顾") print("="*70) print("1. 二维高斯函数: G(x,y) = (1/(2πσ²)) * exp(-(x²+y²)/(2σ²))") print("2. 可分离性: G(x,y) = G(x) * G(y)") print("3. 离散核大小经验公式: size ≈ 6σ + 1 (覆盖99.7%能量)") print("4. 计算复杂度: O(m×n) → O(m+n)（通过可分离性优化）")