1. 用imnoise2生成erlang(爱尔兰)噪声，调用时用默认的a=2， b=5即可（对照表4.1查看对应的均值和方差）。对图像lena_gray_512.tif, 添加此噪声，得到带噪图。 （1）调用roipoly在带噪图中选取噪声样本图像，计算其直方图，估计噪声参数，与已知噪声参数a,b对比，画出噪声图，以及样本图像的直方图。 （2） 调用spfilt对带噪图进行去噪，比较每种滤波的效果，给出理想的方案。

首先，需要导入以下库：numpy、matplotlib、scipy、imageio。

1. 添加Erlang噪声：

import numpy as np
import matplotlib.pyplot as plt
import imageio

def erlang_noise(image, a=2, b=5):
    # 生成Erlang噪声
    noise = np.random.gamma(a, 1/b, image.shape)
    # 将图像与噪声相加
    noisy_image = image + noise
    # 将像素值限制在0到255之间
    noisy_image = np.clip(noisy_image, 0, 255)
    # 转换为整数类型
    noisy_image = noisy_image.astype(np.uint8)
    return noisy_image

# 读取图像
image = imageio.imread('lena_gray_512.tif')
# 添加Erlang噪声
noisy_image = erlang_noise(image)

# 显示原图和带噪图
fig, ax = plt.subplots(1, 2, figsize=(8, 4))
ax[0].imshow(image, cmap='gray')
ax[0].set_title('Original')
ax[1].imshow(noisy_image, cmap='gray')
ax[1].set_title('Noisy')
plt.show()

2. 选取噪声样本图像并计算直方图：

from scipy.stats import gamma

def estimate_erlang_params(image):
    # 选取噪声样本图像
    sample = plt.imread('sample.png')
    # 计算直方图
    hist, bins = np.histogram(sample.ravel(), bins=256, range=[0, 256])
    # 归一化
    hist = hist / np.sum(hist)
    # 估计噪声参数
    a_hat, b_hat, _, _ = gamma.fit(sample.ravel(), floc=0)
    return a_hat, b_hat, hist

# 在带噪图中选取噪声样本图像
from matplotlib.widgets import RectangleSelector

def onselect(eclick, erelease):
    global sample
    sample = noisy_image[int(erelease.ymin):int(erelease.ymax), int(eclick.xmin):int(eclick.xmax)]
    plt.imshow(sample, cmap='gray')
    plt.axis('off')
    plt.savefig('sample.png', bbox_inches='tight')

fig, ax = plt.subplots(figsize=(6, 6))
ax.imshow(noisy_image, cmap='gray')
ax.set_title('Select noise sample')
plt.axis('off')
selector = RectangleSelector(ax, onselect, drawtype='box', useblit=True, button=[1])
plt.show()

# 估计噪声参数并绘制噪声图和直方图
a_hat, b_hat, hist = estimate_erlang_params(sample)
print('Estimated parameters: a={:.2f}, b={:.2f}'.format(a_hat, b_hat))

fig, ax = plt.subplots(1, 2, figsize=(8, 4))
ax[0].imshow(noisy_image, cmap='gray')
ax[0].set_title('Noisy')
ax[1].plot(hist)
ax[1].set_title('Histogram')
plt.show()

3. 使用滤波器去噪：

from scipy.ndimage import filters

# 定义不同的滤波器
filters_list = [
    ('Median', lambda x: filters.median_filter(x, size=3)),
    ('Gaussian', lambda x: filters.gaussian_filter(x, sigma=1)),
    ('Wiener', lambda x: filters.wiener(x, (3, 3))),
    ('Bilateral', lambda x: filters.gaussian_filter(x, sigma=1, mode='nearest'))
]

# 对带噪图像进行滤波，并显示结果
fig, axs = plt.subplots(2, 2, figsize=(8, 8))
axs = axs.ravel()
for i, (name, filt) in enumerate(filters_list):
    denoised_image = filt(noisy_image)
    axs[i].imshow(denoised_image, cmap='gray')
    axs[i].set_title(name)
    axs[i].axis('off')
plt.show()

从结果可以看出，中值滤波和维纳滤波的效果最好，因为它们能够有效地去除噪声，同时保持图像的细节信息。因此，我们可以选择这两种滤波器的组合作为理想的去噪方案。