第5章 Python 数字图像处理(DIP) - 图像复原与重建6

第5章 Python 数字图像处理(DIP) - 图像复原与重建6 - 椒盐噪声

标题

- - - 椒盐噪声

椒盐噪声

如果 $k$ 是一幅数字图像中表示灰度值的比特数，则灰度值可能是 $0, 2^k -1]$ 。椒盐噪声的PDF为：

$\begin{cases} P_s, & z = 2^k -1 \\ P_p, & z=0 \\ 1-(P_s + P_p), & z=V \end{cases} \tag{5.16}$

$0<V<2^k-1 $的任意整数

均值和方差为
$zˉ=(0)Pp+K(1−Ps−Pp)+(2k−1)Ps(5.17)\bar z =(0)P_p + K(1 -P_s - P_p) + (2^k -1) P_s \tag{5.17}$
$σ2=(0−zˉ)2Pp+(K−zˉ)2(1−Ps−Pp)+(2k−1)2Ps(5.18)\sigma^2 =(0-\bar z)^2 P_p + (K - \bar z)^2 ( 1 -P_s -P_p) + (2^k -1)^2 P_s \tag{5.18}$

添加椒盐噪声的方法
令 $η(x,y)\eta(x, y)$ 表示是一幅椒盐噪声图像，其密度值满足式(5.16)。我们使用椒盐噪声污染图像 $f (x, y)$ 的方法为：
在 $f$ 中 $η\eta$ 为0的所有位置赋0，在 $f$ 中 $η\eta$ 为 $2^k-1$ 的所有位置赋 $2^k-1$ 值，保留 $f$ 中 $η\eta$ 为 $V$ 的所有位置的值不变。
噪声密度
像素被盐粒或胡椒噪声污染的概率 $P$ 为 $P=P_s + P_p$ 。 $P$ 称为噪声密度。

def salt_pepper_pdf(ps=0.1, pp=0.1):"""create salt and pepper PDF, math $$P(z) = \begin{cases} P_s, & z = 2^k -1 \\ P_p, & z=0 \\ 1-(P_s + P_p), & z=V \end{cases}$$param: z: input grayscale value of iamgeparam: v:  float, while z = v, probability of other valuesparam: ps: float, probability of the saltparam: pp: float, probability of the pepper"""salt_pepper = np.zeros([3])salt_pepper[0] = pssalt_pepper[2] = ppsalt_pepper[1] = 1 - (ps + pp)return salt_pepper

def add_salt_pepper(img, ps=0.01, pp=0.01):"""add salt pepper noise to imageparam: img: input image, uint8 [0, 255]param: ps:  probability of salt noise, which is white noise, default is 0.01param: pp:  probability of peper noise, which is black noise, default is 0.01return image with salt pepper noise, [0, 255]"""h, w = img.shape[:2]mask = np.random.choice((0, 0.5, 1), size=(h, w), p=[pp, (1-ps-pp), ps])img_out = img.copy()img_out[mask==1] = 255img_out[mask==0] = 0return img_out

def add_salt_pepper_1(img, prob):output = np.zeros(img.shape, np.uint8)thres = 1 - probfor i in range(img.shape[0]):for j in range(img.shape[1]):rdn = np.random.random()if rdn < prob:output[i][j] = 0elif rdn > thres:output[i][j] = 255else:output[i][j] = img[i][j]return output

# 椒盐噪声
z = np.array([0, 128, 255])
ps = 0.1
pp = 0.3
k = 8
z_ = 0 * pp + k * (1 - ps - pp) + (2**k - 1) * ps
sigma = (0 - z_) ** 2 * pp + (k - z_) ** 2 * (1 - ps - pp) + (2**k -1) ** 2 * psprint(f"z_ -> {z_}, sigma -> {sigma}")salt_pepper = salt_pepper_pdf(ps=ps, pp=pp)plt.figure(figsize=(9, 6))
plt.bar(z, salt_pepper), plt.xticks([]), plt.yticks([])
plt.show()

z_ -> 30.3, sigma -> 7076.301

在这里插入图片描述

# 椒盐噪声，可随意控制胡椒与盐粒的概率大小
# img_ori = np.ones((512, 512)) * 128
img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0503 (original_pattern).tif", 0)img_salt_pepper = add_salt_pepper(img_ori, ps=0.05, pp=0.02)plt.figure(figsize=(9, 6))
plt.subplot(121), plt.imshow(img_ori, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(img_salt_pepper, 'gray', vmin=0, vmax=255), plt.xticks([]), plt.yticks([])plt.tight_layout()
plt.show()

在这里插入图片描述

hist, bins = np.histogram(img_salt_pepper.flatten(), bins=255, range=[0, 255], density=True)
bar = plt.bar(bins[:-1], hist[:])

在这里插入图片描述

# 一些重要的概率密度函数
a = 2
b = 8z = np.linspace(0, 10, 200)gaussian = gauss_pdf(z)
rayleigh = rayleigh_pdf(z, a=a, b=b)
ireland  = ireland_pdf(z, a=a, b=b-6)
exponent = exponential_pdf(z, a=a)
average  = average_pdf(z, a=a, b=b)#=============椒盐=============
ps = 0.2
pp = 0.3
salt_pepper = salt_pepper_pdf(ps=ps, pp=pp)show_list = ['gaussian', 'rayleigh', 'ireland', 'exponent', 'average', 'salt_pepper']fig = plt.figure(figsize=(16, 8))for i in range(len(show_list)):ax = fig.add_subplot(2, 3, i+1)if i == 5:z = np.array([0, 100, 255])ax.bar(z, eval(show_list[i]), width=2), ax.set_xticks([0, 100, 255]), ax.set_yticks([0, eval(show_list[i]).max() ])ax.set_xlim(-1, 257), ax.set_ylim([0, eval(show_list[i]).max() + 0.05]), ax.set_title(show_list[i])else:ax.plot(z, eval(show_list[i])), ax.set_xticks([0, 10]), ax.set_yticks([0, eval(show_list[i]).max() ])ax.set_xlim(0, 10), ax.set_ylim([0, eval(show_list[i]).max() + 0.05]), ax.set_title(show_list[i])
plt.tight_layout()
plt.show()

在这里插入图片描述

一般重要噪声的例子

高斯噪声
- 由电子电路及（光照不足和/或高温引起的）传感器噪声等因素导致的
瑞利噪声
- 有助于表征距离成像中的噪声现象
指数和伽马密度
- 在激光成像中广泛应用
冲激噪声
- 出现在成像期间的快速瞬变（如开关故障）
均匀密度
- 是对实际情况最起码的描述

# 一些重要的噪声
img_ori = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0503 (original_pattern).tif", 0)img_gauss       = add_gaussian_noise(img_ori, mu=0, sigma=0.05)
img_rayleigh    = add_rayleigh_noise(img_ori, a=1)
img_gamma       = add_gamma_noise(img_ori, scale=2)
img_exponent    = add_exponent_noise(img_ori, scale=3)
img_average     = add_average_noise(img_ori, mean=10, sigma=1.5)
ps = 0.05
pp = 0.02
img_salt_pepper = add_salt_pepper(img_ori, ps=ps, pp=pp)show_list = ['img_gauss', 'img_rayleigh', 'img_gamma', 'img_exponent', 'img_average', 'img_salt_pepper']fig = plt.figure(figsize=(14.5, 20))for i in range(len(show_list)):if i >= 3:# 显示图像ax = fig.add_subplot(4, 3, i + 3 + 1)ax.imshow(eval(show_list[i]), 'gray'), ax.set_xticks([]), ax.set_yticks([]), ax.set_title(show_list[i].split('_')[-1])# 对应图像的直方图ax = fig.add_subplot(4, 3, i + 1 + 6)hist, bins = np.histogram(eval(show_list[i]).flatten(), bins=255, range=[0, 255], density=True)bar = ax.bar(bins[:-1], hist[:]), ax.set_xticks([]), ax.set_yticks([]),else:# 显示图像ax = fig.add_subplot(4, 3, i + 1)ax.imshow(eval(show_list[i]), 'gray'), ax.set_xticks([]), ax.set_yticks([]), ax.set_title(show_list[i].split('_')[-1])# 对应图像的直方图ax = fig.add_subplot(4, 3, i + 1 + 3)hist, bins = np.histogram(eval(show_list[i]).flatten(), bins=255, range=[0, 255], density=True)bar = ax.bar(bins[:-1], hist[:]), ax.set_xticks([]), ax.set_yticks([]),plt.tight_layout()
plt.show()