2026/1/10 9:07:50
网站建设
项目流程
合肥网站建设首选众龙,外贸网站外链怎么做,品牌网站建设四川,做环球资源网站有没有效果目录 3.3 均值滤波#xff08;续#xff09;
频域理解与信号处理视角
边界效应分析
计算复杂度优化
变体#xff1a;加权均值滤波
均值滤波完整Python程序
3.4 高斯滤波
核心思想与直观理解
数学原理#xff1a;二维高斯函数
高斯核的离散化与创建
高斯滤波的性质…目录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是归一化频率。边界效应分析当滤波器核移动到图像边界时需要特殊处理。常见的边界处理策略有零填充Zero-padding在图像外围填充0。优点简单缺点在边界处引入暗边效应复制填充Replicate复制最近的边界像素值。优点避免暗边缺点可能造成边界模糊镜像填充Reflect镜像边界外的像素。优点保持连续性缺点计算稍复杂循环填充Wrap假设图像是周期性的。优点数学上优美缺点实际图像通常不是周期性的计算复杂度优化对于大尺寸均值滤波直接卷积的计算复杂度是O(m×n×M×N)其中m×n是核大小M×N是图像大小。我们可以使用积分图像或滑动窗口累加技术将复杂度降至O(M×N)。积分图像算法计算积分图像I_ΣI_Σ(x,y) I(x,y) I_Σ(x-1,y) I_Σ(x,y-1) - I_Σ(x-1,y-1)任意矩形区域的像素和可以通过积分图像的4个值快速计算sum I_Σ(x2,y2) - I_Σ(x1-1,y2) - I_Σ(x2,y1-1) I_Σ(x1-1,y1-1)均值 sum / (矩形面积)变体加权均值滤波标准均值滤波对所有邻居一视同仁。但直觉上离中心越近的像素应该越相关。因此出现了加权均值滤波给予不同位置不同的权重。常用的加权方案距离加权权重与到中心的距离成反比高斯加权权重服从高斯分布这过渡到了高斯滤波均值滤波完整Python程序import cv2 import numpy as np import matplotlib.pyplot as plt from scipy import ndimage import time def mean_filter_custom(image, kernel_size3, padding_modereflect): 自定义均值滤波函数 :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), dtypenp.float32) / (kernel_size * kernel_size) # 获取填充大小 pad kernel_size // 2 # 根据填充模式处理图像 if padding_mode zero: padded np.pad(image.astype(np.float32), pad, modeconstant, constant_values0) elif padding_mode replicate: padded np.pad(image.astype(np.float32), pad, modeedge) elif padding_mode reflect: padded np.pad(image.astype(np.float32), pad, modereflect) elif padding_mode wrap: padded np.pad(image.astype(np.float32), pad, modewrap) else: raise ValueError(f不支持的填充模式: {padding_mode}) # 初始化输出图像 output np.zeros_like(image, dtypenp.float32) # 手动实现卷积为了教学目的 height, width image.shape for i in range(pad, height pad): for j in range(pad, width pad): # 提取邻域 neighborhood padded[i-pad:ipad1, j-pad:jpad1] # 计算加权和 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_size3): 使用积分图像优化的均值滤波 :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((height1, width1), dtypenp.float32) # 计算积分图像 for i in range(1, height1): row_sum 0 for j in range(1, width1): 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[x21, y21] - integral[x1, y21] - integral[x21, 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_size15): 比较不同实现方式的性能 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(fOpenCV 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, modereflect) result_scipy np.clip(result_scipy, 0, 255).astype(np.uint8) scipy_time time.time() - start_time print(fSciPy 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, ksize3) sobel_y cv2.Sobel(original, cv2.CV_64F, 0, 1, ksize3) gradient_orig np.sqrt(sobel_x**2 sobel_y**2) sobel_x_f cv2.Sobel(filtered, cv2.CV_64F, 1, 0, ksize3) sobel_y_f cv2.Sobel(filtered, cv2.CV_64F, 0, 1, ksize3) 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), dtypenp.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, cmapgray, vmin0, vmax255) axes[0, 0].set_title(原始图像) axes[0, 0].axis(off) axes[1, 0].imshow(img_gaussian, cmapgray, vmin0, vmax255) axes[1, 0].set_title(高斯噪声图像\n(σ25)) axes[1, 0].axis(off) axes[2, 0].imshow(img_salt_pepper, cmapgray, vmin0, vmax255) 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, idx1].imshow(filtered_gaussian, cmapgray, vmin0, vmax255) axes[0, idx1].set_title(f高斯噪声均值滤波\n核大小{ksize}×{ksize}) axes[0, idx1].axis(off) # 处理椒盐噪声图像 filtered_sp mean_filter_custom(img_salt_pepper, ksize, reflect) axes[1, idx1].imshow(filtered_sp, cmapgray, vmin0, vmax255) axes[1, idx1].set_title(f椒盐噪声均值滤波\n核大小{ksize}×{ksize}) axes[1, idx1].axis(off) # 分析边缘保留 grad_orig, grad_filt, retention analyze_edge_preservation(img, filtered_gaussian) axes[2, idx1].imshow(grad_filt, cmaphot) axes[2, idx1].set_title(f梯度图\n边缘保留率: {retention:.2%}) axes[2, idx1].axis(off) plt.tight_layout() plt.savefig(mean_filter_comparison.png, dpi300, bbox_inchestight) # 4. 性能比较 print(\n *60) print(均值滤波性能与效果分析) print(*60) # 测试不同实现方式的性能 perf_results compare_performance(img, kernel_size15) # 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), dtypenp.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, cmapgray) axes2[0, i].set_title(f边界处理: {mode}) axes2[0, i].axis(off) # 显示中心行剖面 center_row filtered_small[25, :] axes2[1, i].plot(center_row, b-, linewidth2) axes2[1, i].set_title(f中心行灰度剖面) axes2[1, i].set_xlabel(列坐标) axes2[1, i].set_ylabel(灰度值) axes2[1, i].grid(True, alpha0.3) axes2[1, i].set_ylim([0, 255]) plt.tight_layout() plt.savefig(boundary_effects.png, dpi300, bbox_inchestight) # 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高斯核的离散化与创建在实际应用中我们需要将连续的高斯函数离散化为一个有限大小的卷积核。创建高斯核的步骤确定核大小通常根据σ选择经验法则是核大小 ≈ 6σ1确保覆盖99.7%的能量计算离散坐标对于核中的每个位置(i,j)计算其到中心的距离应用高斯公式计算每个位置的权重归一化确保所有权重之和为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]]高斯滤波的性质与优势旋转对称性二维高斯函数是旋转对称的这意味着高斯滤波在各个方向上的平滑效果相同。可分离性这是高斯滤波最重要的性质之一。二维高斯函数可以分解为两个一维高斯函数的乘积G(x,y) G(x) * G(y) (1/(√(2π)σ))exp(-x²/(2σ²)) * (1/(√(2π)σ))exp(-y²/(2σ²))这意味着二维高斯滤波可以分解为两次一维高斯滤波先水平后垂直将计算复杂度从O(m×n)降低到O(mn)。傅里叶变换性质高斯函数的傅里叶变换仍是高斯函数。在频域中高斯滤波器是理想的低通滤波器没有振铃效应。尺度空间理论高斯滤波是构建图像尺度空间的基础这在SIFT等特征检测算法中至关重要。σ参数的影响标准差σ是高斯滤波的关键参数小σ如0.5-1.0高斯核尖锐主要平滑高频噪声保留大部分细节中σ如1.5-2.5平衡噪声抑制和细节保留大σ如3.0以上高斯核宽平强平滑效果大量细节丢失与均值滤波的对比特性均值滤波高斯滤波权重分配均匀高斯分布中心重边缘轻边缘保持差较好计算复杂度O(m×n)O(m×n) 或 O(mn)可分离时频域响应sinc函数有旁瓣高斯函数无旁瓣振铃效应可能有无参数控制仅核大小核大小和σ高斯滤波完整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_modereflect): 可分离高斯滤波实现先水平后垂直 :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)), modeconstant) 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:jsize] * kernel_1d) else: # 使用SciPy的convolve1d支持多种边界模式 temp ndimage.convolve1d(img_float, kernel_1d, axis1, modepadding_mode) # 垂直方向卷积 if padding_mode zero: padded np.pad(temp, ((pad_width, pad_width), (0, 0)), modeconstant) 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:isize, j] * kernel_1d) else: output ndimage.convolve1d(temp, kernel_1d, axis0, modepadding_mode) # 裁剪并转换类型 output np.clip(output, 0, 255).astype(np.uint8) return output def analyze_frequency_response(sigma_values, kernel_size31): 分析不同σ值的高斯滤波器频率响应 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)), modeconstant) 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), idx1) im1 ax1.imshow(kernel, cmaphot, interpolationnearest) ax1.set_title(fσ{sigma}\n高斯核) ax1.axis(off) plt.colorbar(im1, axax1, fraction0.046, pad0.04) # 绘制3D表面图 ax2 plt.subplot(3, len(sigma_values), idx1len(sigma_values), projection3d) 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, cmapviridis, alpha0.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), idx12*len(sigma_values)) center magnitude.shape[0] // 2 profile magnitude[center, center-50:center50] ax3.plot(profile, b-, linewidth2) ax3.set_title(f频率响应剖面) ax3.set_xlabel(空间频率) ax3.set_ylabel(幅度(dB)) ax3.grid(True, alpha0.3) plt.tight_layout() return fig def compare_gaussian_mean(image, sigma2.0, kernel_size15): 比较高斯滤波和均值滤波的效果 # 添加高斯噪声 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_levels4, sigma1.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), dtypenp.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, dpi300, bbox_inchestight) # 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, sigma2.0, kernel_size15) print(f\n 质量指标对比:) print(f 高斯滤波 - PSNR: {comparison[metrics][gaussian][psnr]:.2f} dB, fSSIM: {comparison[metrics][gaussian][ssim]:.4f}) print(f 均值滤波 - PSNR: {comparison[metrics][mean][psnr]:.2f} dB, fSSIM: {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_levels5, sigma1.6) # 7. 可视化所有结果 print(\n7. 生成可视化结果...) # 创建综合可视化图 fig2, axes2 plt.subplots(3, 4, figsize(16, 12)) # 原始图像和噪声图像 axes2[0, 0].imshow(img, cmapgray) axes2[0, 0].set_title(原始图像) axes2[0, 0].axis(off) axes2[0, 1].imshow(comparison[noisy], cmapgray) 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, cmapgray) axes2[row, col*2].set_title(f高斯滤波 σ{sigma}) axes2[row, col*2].axis(off) # 显示剖面线 center_line filtered[150, :] axes2[row, col*21].plot(center_line, b-, linewidth2) axes2[row, col*21].set_title(f中心行剖面 (σ{sigma})) axes2[row, col*21].set_xlabel(列坐标) axes2[row, col*21].set_ylabel(灰度值) axes2[row, col*21].grid(True, alpha0.3) axes2[row, col*21].set_ylim([0, 255]) plt.tight_layout() plt.savefig(gaussian_filter_results.png, dpi300, bbox_inchestight) # 8. 可视化高斯金字塔 fig3, axes3 plt.subplots(1, len(pyramid), figsize(15, 4)) for i, level_img in enumerate(pyramid): axes3[i].imshow(level_img, cmapgray) axes3[i].set_title(f金字塔层级 {i}\n大小: {level_img.shape}) axes3[i].axis(off) plt.tight_layout() plt.savefig(gaussian_pyramid.png, dpi300, bbox_inchestight) # 9. 实际应用示例文档图像预处理 print(\n8. 实际应用文档图像预处理示例...) # 模拟文档图像文字噪声 doc_img np.ones((200, 300), dtypenp.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, cmapgray) ax.set_title(title) ax.axis(off) plt.tight_layout() plt.savefig(document_processing.png, dpi300, bbox_inchestight) # 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(mn)通过可分离性优化)