实习背景与目的意义¶

本次实习旨在通过数字图像处理的若干基础操作，让学生深入理解图像滤波的相关理论与实践内容，以提升对数字图像处理基础知识的掌握。本次实验内容涵盖了图像的彩色合成、噪声模拟、滤波器实现、二维傅里叶变换及频率域滤波操作，着重训练学生对图像基本操作的程序实现，并为后续复杂图像处理打下基础。具体目标包括掌握不同滤波技术（如中值滤波、平滑滤波、拉普拉斯滤波）的应用效果，以及通过频域操作（如高通和低通滤波）分析和处理图像特性，进一步提升对图像特性与滤波方法的理解。

实习内容与步骤¶

  1. 假彩色合成与保存：通过读取影像数据，调整波段顺序以实现标准假彩色合成。将近红外波段（NIR）、红波段（Red）和绿波段（Green）分别赋予RGB三个通道，合成为彩色影像并保存。

  2. 图像加噪：加入椒盐噪声和高斯噪声两种随机噪声，模拟含噪图像。其中，椒盐噪声通过随机极端值替换图像像素，高斯噪声通过高斯分布上的随机值叠加像素值，实现不同噪声效果的影像生成。

  3. 空间滤波：针对不同噪声图像进行滤波处理：椒盐噪声图像采用中值滤波降低噪声干扰，并保留边缘细节。高斯噪声图像使用均值滤波实现平滑处理。原始灰度图像通过拉普拉斯滤波进行锐化处理，增强边缘特征。

  4. 傅里叶变换：对原始影像和椒盐噪声图像进行频域变换，输出影像的频谱特性，便于分析图像的高低频信息分布。

  5. 高通与低通滤波：利用傅里叶变换设计高通滤波器和低通滤波器，分别用于突出高频边缘特征和降低平滑噪声干扰。 滤波后通过逆变换还原图像，验证频域滤波效果并保存滤波结果。

  6. 去噪结果定量评价，不同的去噪方法横向评价

（1）彩色合成与保存¶

任务：读取data目录下GF2.tif数据，将标准假彩色合成结果保存到./result/GF2-false.tif

In [19]:
import rasterio
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")

with rasterio.open('data/GF2.tiff') as src:
    meta = src.meta.copy()
    img = src.read().astype(np.float32)
    meta.update({'dtype': 'float32'})

def normized(data):
    min_val = np.min(data)
    max_val = np.max(data)
    normed_data = (data - min_val) / (max_val - min_val)
    return normed_data

for i in range(img.shape[0]):
    img[i, :, :] = normized(img[i, :, :])

blue = img[0, :, :]
green = img[1, :, :]
red = img[2, :, :]
nir = img[3, :, :]

true_rgb = np.stack([red, green, blue])
false_rgb = np.stack([nir, red, green])
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.title('True Color Composite')
plt.imshow(np.transpose(true_rgb, (1, 2, 0)))
plt.axis('off')
plt.subplot(1, 2, 2)
plt.title('False Color Composite')
plt.imshow(np.transpose(false_rgb, (1, 2, 0)))
plt.axis('off')


false_rgb_meta = meta.copy()
false_rgb_meta.update({'count': 3})

with rasterio.open('result/GF2-false.tif', 'w', **false_rgb_meta) as src:
    src.write(false_rgb)

# free memory
del img
del blue
del green
del red
del nir
del true_rgb
del false_rgb
del meta
del false_rgb_meta
No description has been provided for this image

（2）图像加噪¶

任务：读取data目录下的JAX_gray.tif数据，对影像分别添加椒盐噪声和高斯噪声，保存到data/JAX_noise_xxx.tif

In [20]:
with rasterio.open('data/JAX_gray.tif') as src:
    meta = src.meta.copy()
    img = src.read(1).astype(np.float32)
    meta.update({'dtype': 'float32'})

img = normized(img)

def add_salt_and_pepper_noise(image, salt_prob, pepper_prob):
    noisy_image = np.copy(image)
    total_pixels = image.size
    num_salt = int(total_pixels * salt_prob)
    num_pepper = int(total_pixels * pepper_prob)

    salt_coords = [np.random.randint(0, i - 1, num_salt) for i in image.shape]
    noisy_image[salt_coords[0], salt_coords[1]] = 1

    pepper_coords = [np.random.randint(0, i - 1, num_pepper) for i in image.shape]
    noisy_image[pepper_coords[0], pepper_coords[1]] = 0

    return noisy_image

def add_gaussian_noise(image, mean=0, var=0.01):
    sigma = var ** 0.5
    gaussian_noise = np.random.normal(mean, sigma, image.shape)
    noisy_image = image + gaussian_noise
    noisy_image = np.clip(noisy_image, 0, 1)
    return noisy_image

salt_pepper_noisy_img = add_salt_and_pepper_noise(img, salt_prob=0.02, pepper_prob=0.02)
gaussian_noisy_img = add_gaussian_noise(img, mean=0, var=0.01)

# 显示三个图像
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
plt.title('Original Image')
plt.imshow(img, cmap='gray')
plt.axis('off')
plt.subplot(1, 3, 2)
plt.title('Salt and Pepper Noise')
plt.imshow(salt_pepper_noisy_img, cmap='gray')
plt.axis('off')
plt.subplot(1, 3, 3)
plt.title('Gaussian Noise')
plt.imshow(gaussian_noisy_img, cmap='gray')
plt.axis('off')

with rasterio.open('result/JAX_noise_salt_pepper.tif', 'w', **meta) as src:
    src.write(salt_pepper_noisy_img, 1)

with rasterio.open('result/JAX_noise_gaussian.tif', 'w', **meta) as src:
    src.write(gaussian_noisy_img, 1)

del img
del salt_pepper_noisy_img
del gaussian_noisy_img
del meta
No description has been provided for this image

（3）空间滤波¶

任务：对(2)中生成的椒盐噪声图像进行中值滤波，对(2)中生成的高斯噪声图像进行均值滤波，对原始灰度图像进行拉普拉斯滤波，将滤波后的结果保存到data/JAX_filter_xx.tif。

In [21]:
def median_filter(image):
    padded_image = np.pad(image, pad_width=1, mode='edge')
    filtered_image = np.zeros_like(image)

    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            neighborhood = padded_image[i:i+3, j:j+3]
            filtered_image[i, j] = np.median(neighborhood)

    return filtered_image

def mean_filter(image):
    padded_image = np.pad(image, pad_width=1, mode='edge')
    filtered_image = np.zeros_like(image)

    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            neighborhood = padded_image[i:i+3, j:j+3]
            filtered_image[i, j] = np.mean(neighborhood)

    return filtered_image

def laplacian_filter(image):
    kernel = np.array([[0, -1, 0],
                       [-1, 4, -1],
                       [0, -1, 0]])
    padded_image = np.pad(image, pad_width=1, mode='edge')
    filtered_image = np.zeros_like(image)

    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            neighborhood = padded_image[i:i+3, j:j+3]
            filtered_image[i, j] = np.sum(neighborhood * kernel)

    filtered_image = np.clip(filtered_image, 0, 1)
    return filtered_image

with rasterio.open('data/JAX_gray.tif') as src:
    img = src.read(1).astype(np.float32)
    meta = src.meta.copy()
    meta.update({'dtype': 'float32'})
    img = normized(img)
    laplacian_filtered_img = laplacian_filter(img)


with rasterio.open('result/JAX_noise_salt_pepper.tif') as src:
    salt_pepper_noisy_img = src.read(1).astype(np.float32)
    salt_pepper_filtered_img = median_filter(salt_pepper_noisy_img)

with rasterio.open('result/JAX_noise_gaussian.tif') as src:
    gaussian_noisy_img = src.read(1).astype(np.float32)
    gaussian_filtered_img = mean_filter(gaussian_noisy_img)

plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
plt.title('Median Filter on Salt and Pepper Noise')
plt.imshow(salt_pepper_filtered_img, cmap='gray')
plt.axis('off')
plt.subplot(1, 3, 2)
plt.title('Mean Filter on Gaussian Noise')
plt.imshow(gaussian_filtered_img, cmap='gray')
plt.axis('off')
plt.subplot(1, 3, 3)
plt.title('Laplacian Filter on Original Image')
plt.imshow(laplacian_filtered_img, cmap='gray')
plt.axis('off')

with rasterio.open('result/JAX_filter_median.tif', 'w', **meta) as src:
    src.write(salt_pepper_filtered_img, 1)

with rasterio.open('result/JAX_filter_mean.tif', 'w', **meta) as src:
    src.write(gaussian_filtered_img, 1)

with rasterio.open('result/JAX_filter_laplacian.tif', 'w', **meta) as src:
    src.write(laplacian_filtered_img, 1)

del img
del salt_pepper_filtered_img
del gaussian_filtered_img
del laplacian_filtered_img
del meta
No description has been provided for this image

（4）傅里叶变换¶

任务：对原始影像和椒盐噪声图像进行傅里叶变换，输出频率谱。

In [23]:
def compute_fft(image):
    f = np.fft.fft2(image)
    fshift = np.fft.fftshift(f)
    magnitude_spectrum = 20 * np.log(np.abs(fshift) + 1)
    return magnitude_spectrum

with rasterio.open('data/JAX_gray.tif') as src:
    img = src.read(1).astype(np.float32)
    meta = src.meta.copy()
    meta.update({'dtype': 'float32'})
    img = normized(img)

original_fft = compute_fft(img)
salt_pepper_fft = compute_fft(salt_pepper_noisy_img)

plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.title('FFT of Original Image')
plt.imshow(original_fft, cmap='gray')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.title('FFT of Salt and Pepper Noisy Image')
plt.imshow(salt_pepper_fft, cmap='gray')
plt.axis('off')

with rasterio.open('result/JAX_fft_original.tif', 'w', **meta) as src:
    src.write(original_fft, 1)
with rasterio.open('result/JAX_fft_salt_pepper.tif', 'w', **meta) as src:
    src.write(salt_pepper_fft, 1)

del img
del salt_pepper_noisy_img
del original_fft
del salt_pepper_fft
del meta
No description has been provided for this image

（5）高通滤波¶

任务：对data目录下的GF2.tiff进行高通滤波(基于傅里叶变换)，保存高通滤波后的结果到data/high_pass_result.png。

（6）低通滤波¶

任务：对data目录下的GF2.tiff进行低通滤波(基于傅里叶变换)，保存低通滤波后的结果到data/low_pass_result.png。

In [26]:
with rasterio.open('data/GF2.tiff') as src:
    meta = src.meta.copy()
    img = src.read(1).astype(np.float32)
    meta.update({'dtype': 'float32'})
    img = normized(img)

def high_pass_filter(image, cutoff_frequency):
    f = np.fft.fft2(image)
    fshift = np.fft.fftshift(f)

    rows, cols = image.shape
    crow, ccol = rows // 2 , cols // 2

    mask = np.ones((rows, cols), np.uint8)
    r = cutoff_frequency
    center = (crow, ccol)
    y, x = np.ogrid[:rows, :cols]
    mask_area = (x - center[1])**2 + (y - center[0])**2 <= r*r
    mask[mask_area] = 0

    fshift_filtered = fshift * mask
    f_ishift = np.fft.ifftshift(fshift_filtered)
    img_back = np.fft.ifft2(f_ishift)
    img_back = np.abs(img_back)

    img_back = normized(img_back)
    return img_back

def low_pass_filter(image, cutoff_frequency):
    f = np.fft.fft2(image)
    fshift = np.fft.fftshift(f)

    rows, cols = image.shape
    crow, ccol = rows // 2 , cols // 2

    mask = np.zeros((rows, cols), np.uint8)
    r = cutoff_frequency
    center = (crow, ccol)
    y, x = np.ogrid[:rows, :cols]
    mask_area = (x - center[1])**2 + (y - center[0])**2 <= r*r
    mask[mask_area] = 1

    fshift_filtered = fshift * mask
    f_ishift = np.fft.ifftshift(fshift_filtered)
    img_back = np.fft.ifft2(f_ishift)
    img_back = np.abs(img_back)

    img_back = normized(img_back)
    return img_back

high_pass_filtered_img = high_pass_filter(img, cutoff_frequency=30)
high_pass_filtered_img = normized(high_pass_filtered_img)
low_pass_filtered_img = low_pass_filter(img, cutoff_frequency=30)
low_pass_filtered_img = normized(low_pass_filtered_img)

plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.title('High Pass Filtered Image')
plt.imshow(high_pass_filtered_img, cmap='gray')
plt.axis('off')
plt.subplot(1, 2, 2)
plt.title('Low Pass Filtered Image')
plt.imshow(low_pass_filtered_img, cmap='gray')
plt.axis('off')

with rasterio.open('data/high_pass_result.png', 'w', **meta) as src:
    src.write(high_pass_filtered_img, 1)

with rasterio.open('data/low_pass_result.png', 'w', **meta) as src:
    src.write(low_pass_filtered_img, 1)

del img
del high_pass_filtered_img
del meta
No description has been provided for this image

遇到的问题与解决策略¶

循环处理效率低：原先逐像素双重 for 循环计算中值滤波耗时过长，改用 NumPy 的 stride tricks 或 OpenCV 中位数函数显著提升速度。

频谱可视化问题：直接对幅值取对数后显示对比度更佳

边缘效应与零填充：空间滤波需对图像边界作对称或零填充，避免卷积时尺寸不匹配或边缘失真。

实验收获与心得¶

理论与实践结合

通过编程实现空间域与频域滤波，深刻理解低频对应图像的整体亮度与缓变部分，高频对应边缘与细节； 卷积核、掩模设计与信号处理理论相互印证，增强了对数字图像处理原理的掌握。 工具与优化

熟练掌握 NumPy、Matplotlib 等常用 Python 库； 体会到向量化与库函数调用在大规模数据处理中的重要性。 问题解决能力

在实验过程中遇到性能瓶颈与调试难题，提高了查阅文档、测试与调优的综合能力； 学会分步骤分析问题，将复杂任务拆分为子模块逐个攻克。 团队协作与时间管理