数字图像处理第二次实验

实验三

技术点分析

根据实验要求，需要实现以下图像空间域滤波技术：

1. 噪声生成：

   • 高斯噪声

   • 椒盐噪声

2. 空间域滤波：

   • 均值滤波（3×3, 5×5, 7×7）

   • 中值滤波（3×3, 5×5, 7×7）

   • 最大值滤波

   • 最小值滤波

3. 图像处理流程：

   • 读取原始图像

   • 添加噪声（高斯/椒盐）

   • 应用各种滤波器

   • 可视化对比结果

完整示例代码

import cv2
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
import random

# 设置中文字体支持
font = FontProperties(fname=r"c:\windows\fonts\simhei.ttf", size=10)  # Windows
# font = FontProperties(fname='/System/Library/Fonts/PingFang.ttc', size=10)  # macOS

def add_gaussian_noise(image, mean=0, sigma=25):
    """添加高斯噪声"""
    noise = np.random.normal(mean, sigma, image.shape).astype(np.uint8)
    noisy = cv2.add(image, noise)
    return np.clip(noisy, 0, 255)

def add_salt_pepper_noise(image, salt_prob=0.01, pepper_prob=0.01):
    """添加椒盐噪声"""
    noisy = np.copy(image)
    # 添加盐噪声（白点）
    salt_mask = np.random.random(image.shape) < salt_prob
    noisy[salt_mask] = 255
    # 添加椒噪声（黑点）
    pepper_mask = np.random.random(image.shape) < pepper_prob
    noisy[pepper_mask] = 0
    return noisy

def apply_filters(noisy_img):
    """应用各种滤波器"""
    # 均值滤波
    mean_3x3 = cv2.blur(noisy_img, (3, 3))
    mean_5x5 = cv2.blur(noisy_img, (5, 5))
    mean_7x7 = cv2.blur(noisy_img, (7, 7))
    
    # 中值滤波
    median_3x3 = cv2.medianBlur(noisy_img, 3)
    median_5x5 = cv2.medianBlur(noisy_img, 5)
    median_7x7 = cv2.medianBlur(noisy_img, 7)
    
    # 最大值滤波
    kernel = np.ones((3, 3), np.uint8)
    max_filter = cv2.dilate(noisy_img, kernel)
    
    # 最小值滤波
    min_filter = cv2.erode(noisy_img, kernel)
    
    return {
        'mean_3x3': mean_3x3,
        'mean_5x5': mean_5x5,
        'mean_7x7': mean_7x7,
        'median_3x3': median_3x3,
        'median_5x5': median_5x5,
        'median_7x7': median_7x7,
        'max_filter': max_filter,
        'min_filter': min_filter
    }

def calculate_psnr(original, processed):
    """计算PSNR评估图像质量"""
    mse = np.mean((original - processed) ** 2)
    if mse == 0:
        return float('inf')
    max_pixel = 255.0
    psnr = 20 * np.log10(max_pixel / np.sqrt(mse))
    return psnr

# 主程序
def main():
    # 1. 读取原始图像
    original_img = cv2.imread('lena.jpg', cv2.IMREAD_GRAYSCALE)
    if original_img is None:
        print("错误：无法读取图像，请检查路径")
        return
    
    # 2. 添加噪声
    gaussian_noisy = add_gaussian_noise(original_img, sigma=30)
    saltpepper_noisy = add_salt_pepper_noise(original_img, salt_prob=0.03, pepper_prob=0.03)
    
    # 3. 应用滤波器
    gaussian_filters = apply_filters(gaussian_noisy)
    saltpepper_filters = apply_filters(saltpepper_noisy)
    
    # 4. 创建画布
    plt.figure(figsize=(18, 8))
    
    # 第一行：高斯噪声处理结果
    images_row1 = [
        original_img,
        gaussian_noisy,
        gaussian_filters['mean_3x3'],
        gaussian_filters['mean_5x5'],
        gaussian_filters['mean_7x7'],
        gaussian_filters['max_filter'],
        gaussian_filters['min_filter']
    ]
    
    titles_row1 = [
        '原始图像',
        '高斯噪声图像',
        '均值滤波(3×3)\nPSNR: {:.2f}'.format(calculate_psnr(original_img, gaussian_filters['mean_3x3'])),
        '均值滤波(5×5)\nPSNR: {:.2f}'.format(calculate_psnr(original_img, gaussian_filters['mean_5x5'])),
        '均值滤波(7×7)\nPSNR: {:.2f}'.format(calculate_psnr(original_img, gaussian_filters['mean_7x7'])),
        '最大值滤波\nPSNR: {:.2f}'.format(calculate_psnr(original_img, gaussian_filters['max_filter'])),
        '最小值滤波\nPSNR: {:.2f}'.format(calculate_psnr(original_img, gaussian_filters['min_filter']))
    ]
    
    # 第二行：椒盐噪声处理结果
    images_row2 = [
        original_img,
        saltpepper_noisy,
        saltpepper_filters['median_3x3'],
        saltpepper_filters['median_5x5'],
        saltpepper_filters['median_7x7'],
        saltpepper_filters['max_filter'],
        saltpepper_filters['min_filter']
    ]
    
    titles_row2 = [
        '原始图像',
        '椒盐噪声图像',
        '中值滤波(3×3)\nPSNR: {:.2f}'.format(calculate_psnr(original_img, saltpepper_filters['median_3x3'])),
        '中值滤波(5×5)\nPSNR: {:.2f}'.format(calculate_psnr(original_img, saltpepper_filters['median_5x5'])),
        '中值滤波(7×7)\nPSNR: {:.2f}'.format(calculate_psnr(original_img, saltpepper_filters['median_7x7'])),
        '最大值滤波\nPSNR: {:.2f}'.format(calculate_psnr(original_img, saltpepper_filters['max_filter'])),
        '最小值滤波\nPSNR: {:.2f}'.format(calculate_psnr(original_img, saltpepper_filters['min_filter']))
    ]
    
    # 5. 显示结果
    for i in range(7):
        # 第一行
        plt.subplot(2, 7, i+1)
        plt.imshow(images_row1[i], cmap='gray')
        plt.title(titles_row1[i], fontproperties=font)
        plt.axis('off')
        
        # 第二行
        plt.subplot(2, 7, i+8)
        plt.imshow(images_row2[i], cmap='gray')
        plt.title(titles_row2[i], fontproperties=font)
        plt.axis('off')
    
    plt.tight_layout()
    plt.savefig('spatial_filtering_results.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    # 6. 分析报告
    print("="*80)
    print("空间域滤波增强实验结果分析")
    print("="*80)
    
    # 高斯噪声滤波效果分析
    print("\n高斯噪声图像滤波效果 (PSNR值越高越好):")
    print("-"*70)
    print(f"均值滤波 3×3: {calculate_psnr(original_img, gaussian_filters['mean_3x3']):.2f} dB")
    print(f"均值滤波 5×5: {calculate_psnr(original_img, gaussian_filters['mean_5x5']):.2f} dB")
    print(f"均值滤波 7×7: {calculate_psnr(original_img, gaussian_filters['mean_7x7']):.2f} dB")
    print(f"最大值滤波: {calculate_psnr(original_img, gaussian_filters['max_filter']):.2f} dB")
    print(f"最小值滤波: {calculate_psnr(original_img, gaussian_filters['min_filter']):.2f} dB")
    
    # 椒盐噪声滤波效果分析
    print("\n椒盐噪声图像滤波效果 (PSNR值越高越好):")
    print("-"*70)
    print(f"中值滤波 3×3: {calculate_psnr(original_img, saltpepper_filters['median_3x3']):.2f} dB")
    print(f"中值滤波 5×5: {calculate_psnr(original_img, saltpepper_filters['median_5x5']):.2f} dB")
    print(f"中值滤波 7×7: {calculate_psnr(original_img, saltpepper_filters['median_7x7']):.2f} dB")
    print(f"最大值滤波: {calculate_psnr(original_img, saltpepper_filters['max_filter']):.2f} dB")
    print(f"最小值滤波: {calculate_psnr(original_img, saltpepper_filters['min_filter']):.2f} dB")
    
    # 滤波效果对比分析
    print("\n滤波效果对比分析:")
    print("-"*70)
    print("1. 高斯噪声处理:")
    print("   - 均值滤波效果随窗口增大而提升，但会导致图像模糊")
    print("   - 最大值滤波会增强亮噪声，最小值滤波会增强暗噪声")
    
    print("\n2. 椒盐噪声处理:")
    print("   - 中值滤波效果最佳，能有效消除孤立的椒盐噪声点")
    print("   - 3×3窗口对细节保留较好，7×7窗口去噪更彻底但模糊更明显")
    print("   - 最大值滤波消除椒噪声(黑点)，最小值滤波消除盐噪声(白点)")
    
    print("\n3. 综合建议:")
    print("   - 高斯噪声: 使用5×5或7×7均值滤波")
    print("   - 椒盐噪声: 使用3×3或5×5中值滤波")
    print("   - 混合噪声: 可结合使用中值滤波和均值滤波")

if __name__ == "__main__":
    main()

代码说明

1. 噪声生成函数：

   • add_gaussian_noise()：添加高斯噪声

   • add_salt_pepper_noise()：添加椒盐噪声

2. 滤波处理函数：

   • apply_filters()：应用所有滤波器并返回结果

   • calculate_psnr()：计算PSNR评估滤波效果

3. 核心滤波方法：

   • 均值滤波：cv2.blur()

   • 中值滤波：cv2.medianBlur()

   • 最大值滤波：cv2.dilate()（膨胀操作）

   • 最小值滤波：cv2.erode()（腐蚀操作）

4. 可视化与分析：

   • 2行×7列对比展示所有结果

   • 标题显示滤波类型和PSNR值

   • 控制台输出详细分析报告

实验结果分析

1. 高斯噪声处理：

   • 均值滤波窗口越大，去噪效果越好，但图像越模糊

   • 最大值/最小值滤波会增强某些噪声成分

2. 椒盐噪声处理：

   • 中值滤波效果最佳，能有效消除孤立噪声点

   • 3×3窗口保留细节，7×7窗口去噪更彻底

   • 最大值滤波消除椒噪声(黑点)，最小值滤波消除盐噪声(白点)

3. PSNR指标：

   • 量化评估滤波效果（值越高越好）

   • 均值滤波对高斯噪声PSNR提升明显

   • 中值滤波对椒盐噪声PSNR提升显著

使用说明

1. 安装依赖库：

   pip install opencv-python matplotlib numpy

2. 准备测试图像（如lena.jpg）放在代码同级目录

3. 运行代码后将显示：

   • 第一行：原始图像、高斯噪声图像及各种滤波结果

   • 第二行：原始图像、椒盐噪声图像及各种滤波结果

4. 控制台输出详细的分析报告和PSNR指标

5. 如果提示没有对应的字体请先下载字体，或者用英文表示图表标签

实验原理

空间域滤波直接在图像像素上操作：

1. 均值滤波：取邻域像素平均值，公式：
   g(x,y)=1mn∑(i,j)∈Sxyf(i,j)g(x,y)=mn1​∑(i,j)∈Sxy​​f(i,j)

   • 有效平滑高斯噪声

   • 窗口越大，模糊效果越明显

2. 中值滤波：取邻域像素中值

   • 有效消除椒盐噪声

   • 保留边缘信息

3. 最大值滤波：取邻域最大值
   g(x,y)=max⁡(i,j)∈Sxyf(i,j)g(x,y)=max(i,j)∈Sxy​​f(i,j)

   • 消除暗噪声点（椒噪声）

4. 最小值滤波：取邻域最小值
   g(x,y)=min⁡(i,j)∈Sxyf(i,j)g(x,y)=min(i,j)∈Sxy​​f(i,j)

   • 消除亮噪声点（盐噪声

实验四

实验主要涉及以下技术点：

1. ​​傅里叶变换与频谱分析​​
2. ​​频域滤波器实现​​（理想/巴特沃斯/高斯滤波器）
3. ​​相位/幅度分离重构​​
4. ​​图像去噪处理​​（椒盐/高斯噪声）
5. ​​高频增强与直方图均衡化​​
6. ​​振铃效应分析​​

示例代码： 

import numpy as np
import cv2
import matplotlib.pyplot as plt
from skimage import data, util

# 1. 傅里叶变换与频谱可视化
def show_fourier_transform(img_path):
    img = cv2.imread(img_path, 0)
    if img is None:
        raise ValueError(f"Image {img_path} not found")
        
    # 傅里叶变换
    dft = np.fft.fft2(img)
    dft_shift = np.fft.fftshift(dft)
    magnitude = np.log(np.abs(dft_shift) + 1)
    
    plt.figure(figsize=(12, 6))
    plt.subplot(121), plt.imshow(img, cmap='gray'), plt.title('Original Image')
    plt.subplot(122), plt.imshow(magnitude, cmap='magma'), plt.title('Fourier Magnitude')
    plt.show()
    return dft_shift

# 2. 频域滤波器实现
def butterworth_lowpass(shape, D0=30, n=2):
    rows, cols = shape
    center = (rows//2, cols//2)
    filter = np.zeros(shape)
    for i in range(rows):
        for j in range(cols):
            dist = np.sqrt((i - center[0])**2 + (j - center[1])**2)
            filter[i, j] = 1 / (1 + (dist/D0)**(2*n))
    return filter

def apply_filter(dft_shift, filter):
    filtered = dft_shift * filter
    inv_shift = np.fft.ifftshift(filtered)
    img_back = np.fft.ifft2(inv_shift)
    return np.abs(img_back)

# 3. 相位/幅度分离重构实验
def phase_vs_magnitude(dft_shift):
    magnitude = np.abs(dft_shift)
    phase = np.angle(dft_shift)
    
    # 仅幅度重构
    recon_mag = np.fft.ifft2(np.fft.ifftshift(magnitude))
    
    # 仅相位重构
    recon_phase = np.fft.ifft2(np.fft.ifftshift(np.exp(1j * phase)))
    
    plt.figure(figsize=(15, 5))
    plt.subplot(131), plt.imshow(np.abs(recon_mag), cmap='gray'), plt.title('Magnitude Only')
    plt.subplot(132), plt.imshow(np.abs(recon_phase), cmap='gray'), plt.title('Phase Only')
    plt.show()

# 4. 噪声处理与去噪
def noise_processing(img):
    # 添加椒盐噪声
    noisy_sp = util.random_noise(img, mode='s&p', amount=0.05)
    # 添加高斯噪声
    noisy_gauss = util.random_noise(img, mode='gaussian', var=0.01)
    
    # 去噪处理（高斯滤波）
    denoised_sp = cv2.GaussianBlur((noisy_sp*255).astype(np.uint8), (5,5), 0)
    denoised_gauss = cv2.GaussianBlur((noisy_gauss*255).astype(np.uint8), (5,5), 0)
    
    plt.figure(figsize=(12, 8))
    plt.subplot(231), plt.imshow(img, cmap='gray'), plt.title('Original')
    plt.subplot(232), plt.imshow(noisy_sp, cmap='gray'), plt.title('Salt & Pepper')
    plt.subplot(233), plt.imshow(noisy_gauss, cmap='gray'), plt.title('Gaussian Noise')
    plt.subplot(235), plt.imshow(denoised_sp, cmap='gray'), plt.title('Denoised S&P')
    plt.subplot(236), plt.imshow(denoised_gauss, cmap='gray'), plt.title('Denoised Gauss')
    plt.tight_layout()
    plt.show()

# 主程序
if __name__ == "__main__":
    try:
        # 示例图像处理
        img_path = 'rect1.bmp'  # 替换为你的图像路径
        
        # 1. 傅里叶变换
        dft = show_fourier_transform(img_path)
        
        # 2. 低通滤波
        img = cv2.imread(img_path, 0)
        bw_filter = butterworth_lowpass(img.shape, D0=30, n=2)
        filtered_img = apply_filter(dft, bw_filter)
        
        plt.figure(figsize=(12, 5))
        plt.subplot(121), plt.imshow(img, cmap='gray'), plt.title('Original')
        plt.subplot(122), plt.imshow(filtered_img, cmap='gray'), plt.title('Butterworth Filtered')
        plt.show()
        
        # 3. 相位/幅度实验
        phase_vs_magnitude(dft)
        
        # 4. 噪声处理（使用自带图像示例）
        noise_processing(data.camera())
        
    except Exception as e:
        print(f"Error: {e}")

关键技术说明：

​​傅里叶变换核心​​

dft = np.fft.fft2(img)
dft_shift = np.fft.fftshift(dft)  # 中心化低频

​​巴特沃斯滤波器实现​​

def butterworth_lowpass(shape, D0=30, n=2):
    filter = np.zeros(shape)
    for i in range(rows):
        for j in range(cols):
            dist = np.sqrt((i - center[0])**2 + (j - center[1])**2)
            filter[i, j] = 1 / (1 + (dist/D0)**(2*n))  # 巴特沃斯公式

​​相位/幅度分离原理​

magnitude = np.abs(dft_shift)     # 获取幅度谱
phase = np.angle(dft_shift)       # 获取相位谱
recon_phase = np.exp(1j * phase)  # 仅相位重建

​

​​振铃效应控制​​

• 通过调整巴特沃斯滤波器的阶数n控制振铃强度
• 高斯滤波器可完全避免振铃（平滑过渡）

​​噪声处理流程​

# 添加椒盐噪声
noisy_sp = util.random_noise(img, mode='s&p', amount=0.05)

# 高斯去噪
denoised = cv2.GaussianBlur(noisy_img, (5,5), 0)

实验建议：

1. 尝试不同截止频率(D0=10/30/50)观察效果
2. 对比理想/巴特沃斯/高斯滤波器的振铃差异
3. 修改巴特沃斯阶数(n=1/2/10)观察振铃变化
4. 比较相位/幅度重构图像的视觉差异
5. 尝试高频增强+直方图均衡化不同顺序的效果