数字图像处理 - 实验作业四 - Python

第五章 图像复原与重建

我们所看到的事物并不是事物的本身……物体本身是什么对我们而言完全是未知的，并且远离我们的感知。除了对物体的感知，我们什么也不知道——伊曼努尔·康德

特别说明：

第三个和第四个实验的代码大部分都不是原创，因为自己理解不是很深，借助了网上的资料，文末给出了参考来源

图像退化/复原过程的模型：

退化图像：
$$g(x,y)=f(x,y)\ast h(x,y)+\eta (x,y)$$

$$G(u,v)=F(u,v)H(u,v)+N(u,v)$$

我们研究的就是在退化之后如何进行复原

1、实现自适应均值滤波，并和算术均值滤波的结果做对比

算术均值滤波器的表达式为
$$\hat{f}(x,y)=\frac{1}{mn}\sum _{(s,t)\in S_{mn}}g(s,t)$$

其中的S_mn代表中心在点（x，y）处，大小为m×n的矩形子图像窗口（邻域）的一组坐标

算术均值滤波会平滑一幅图像中的局部变化。虽然模糊了结果，但是降低了噪声

自适应均值滤波的表达式为
$$\hat{f}(x,y)=g(x,y)-\frac{{\sigma }_{\eta }^2}{{\sigma }_L^2}\left[g(x,y)-m_L\right]$$

参数以及处理的原则为

其中需要估计的就是噪声的方差

实现的结果如下所示：

可见自适应均值滤波在算术均值滤波的基础上，对于边缘的处理更加好了。消除了原有的模糊信息，保留的更多的边缘的部分

实现的代码如下所示

# 自适应均值滤波和算术均值滤波
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image


def mean_filter(src_image, kernel_size=3):
    """The arithmetic average filter"""
    edge_size = int(kernel_size / 2)
    padded_img = np.pad(src_image, edge_size, 'edge')

    height, width = src_image.shape

    src_after_filter = np.zeros_like(src_image, dtype="float64")

    for i in range(height):
        for j in range(width):
            box = padded_img[i:i + kernel_size, j:j + kernel_size]
            src_after_filter[i, j] = np.mean(box)

    return src_after_filter


def self_adaptaion_mean_filter(src_image, kernel_size=3, variance=1000):
    """self adaptaion mean filter"""
    edge_size = int(kernel_size / 2)
    padded_img = np.pad(src_image, edge_size, 'edge')

    height, width = src_image.shape

    src_after_filter = np.zeros_like(src_image, dtype="float64")

    for i in range(height):
        for j in range(width):
            box = padded_img[i:i + kernel_size, j:j + kernel_size]
            box_mean = np.mean(box)
            box_var = np.var(box)
            src_after_filter[i, j] = src_image[i, j] - (variance / box_var) * (src_image[i, j] - box_mean)

    return src_after_filter


def main(img_path):
    src = np.array(Image.open(img_path).convert("L"))
    noise = np.random.normal(0, 15, size=src.shape)
    noise_src = src + noise

    img_mean_filter = mean_filter(noise_src, 7)
    img_self_adap_mean_filter = self_adaptaion_mean_filter(noise_src, 7, 225)

    img_list = [src, noise_src, img_mean_filter, img_self_adap_mean_filter]
    img_name = ["原图像", "被加性高斯噪声污染的图像", "7×7算术均值滤波", "7×7自适应均值滤波"]

    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False

    _, axs = plt.subplots(2, 2, figsize=(12, 12))

    for i in range(2):
        for j in range(2):
            axs[i][j].imshow(img_list[i * 2 + j], cmap='gray')
            axs[i][j].set_title(img_name[i * 2 + j])
            axs[i][j].axes.get_xaxis().set_visible(False)
            axs[i][j].axes.get_yaxis().set_visible(False)

    plt.show()


if __name__ == '__main__':
    main('man.jpg')

---

2、实现自适应中值滤波，并和中值滤波的结果做对比

中值滤波器的表达式为
f ^ ( x , y ) = m e d i a n { g ( s , t ) } ( s , t ) ∈ S x y \hat{f}(x, y)=median\{g(s,t)\} \quad (s,t)∈S_{xy} f^​(x,y)=median{g(s,t)}(s,t)∈Sxy​

不同于平均值，这里取的是中值，对于椒盐噪声类型的噪声，可以提供更好的去噪能力，且比相同尺寸的线性平滑滤波器引起的模糊更少

自适应中值滤波的参数为

• $S_{xy}$：滤波器的作用区域，该区域中心点为图像中第y行第x列个像素点；
• $Z_{min}$：$S_{xy}$中最小的灰度值；
• $Z_{max}$$S_{：}xy$中最大的灰度值；
• $Z_{med}$：$S_{xy}$中所有灰度值的中值；
• $Z_{xy}$：表示图像中第y行第x列个像素点的灰度值；
• $S_{max}$：$S_{xy}$所允许的最大窗口尺寸；

自适应中值滤波器分为以下两个过程，A和B：

Process A

1. $A1$ =$Z_{med}-Z_{min}$
2. $A2$ =$Z_{med}-Z_{max}$
3. 如果$A1>0$ 且$ A2<0$，则跳转到B
4. 否则，增大窗口的尺寸
5. 如果增大后的尺寸$\le S_{max}$，则重复A
6. 否则，直接输出$Z_{med}$

Process B

1. $B1$ =$Z_{xy}-Z_{min}$
2. $B2$ =$Z_{xy}-Z_{max}$
3. 如果$B1>0$ 且 $B2<0$，则输出$Z-xy$
4. 否则输出$Z_{med}$

实现结果如下所示：

可见自适应中值滤波对于被污染的图像的去噪有着更好的效果，而一般的中值滤波残留的噪声还是比较多的。

实现代码如下所示：

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
import sys


def salt(image, prob):
    """prob：Noise Scale"""
    output = np.zeros(image.shape, np.float64)
    prob_other = 1 - prob
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            rdn = np.random.random()
            if rdn < prob:
                output[i][j] = 0
            elif rdn > prob_other:
                output[i][j] = 255
            else:
                output[i][j] = image[i][j]
    return output


def median_filter(src_image, kernel_size):
    """Median filter"""
    edge_size = int(kernel_size / 2)
    padded_img = np.pad(src_image, edge_size, 'edge')

    height, width = src_image.shape

    src_after_filter = np.zeros_like(src_image, dtype="float64")

    for i in range(height):
        for j in range(width):
            box = padded_img[i:i + kernel_size, j:j + kernel_size]
            src_after_filter[i, j] = np.median(box)

    return src_after_filter


def self_adaptation_median_filter(src_image, kernel_size, output):
    """self-adaptation median filter"""

    height, width = src_image.shape

    num = height * width
    loc = 0

    src_after_filter = np.zeros_like(src_image, dtype="float64")

    for i in range(height):
        for j in range(width):
            rst = filter_process(src_image, kernel_size, i, j, 7)
            src_after_filter[i][j] = rst
            loc = loc + 1
            count = loc / num * 100
            output.write(f'\rcomplete percent:{count:.0f}')

    return src_after_filter


def filter_process(src_image, kernel_size, i, j, s_max=7):
    edge_size = int(kernel_size / 2)
    padded_img = np.pad(src_image, edge_size, 'edge')

    box = padded_img[i:i + kernel_size, j:j + kernel_size]
    z_min = np.min(box)
    z_max = np.max(box)
    z_med = np.median(box)
    z_xy = src_image[i, j]

    rst = a_process(src_image, kernel_size, i, j, s_max, z_xy, z_min, z_max, z_med)

    return rst


def a_process(src_image, kernel_size, i, j, s_max, z_xy, z_min, z_max, z_med):
    A1 = z_med - z_min
    A2 = z_med - z_max

    if A1 > 0 and A2 < 0:
        return b_process(z_xy, z_min, z_max, z_med)
    elif (kernel_size + 2) <= s_max:
        return filter_process(src_image, kernel_size + 2, i, j, s_max)
    else:
        return z_med


def b_process(z_xy, z_min, z_max, z_med):
    B1 = z_xy - z_min
    B2 = z_xy - z_max

    if B1 > 0 and B2 < 0:
        return z_xy
    else:
        return z_med


def add_noise(src_image, noise_method, prob):
    return noise_method(src_image, prob)


def main(img_path):
    src = np.array(Image.open(img_path).convert("L"))
    noise_src = add_noise(src, salt, 0.1)

    output = sys.stdout

    img_median_filter = median_filter(noise_src, 3)
    img_self_adap_median_filter = self_adaptation_median_filter(noise_src, 3, output)

    output.flush()

    img_list = [src, noise_src, img_median_filter, img_self_adap_median_filter]
    img_name = ["原图像", "被椒盐噪声污染的图像", "3×3中值滤波", "自适中值滤波"]

    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False

    _, axs = plt.subplots(2, 2, figsize=(12, 12))

    for i in range(2):
        for j in range(2):
            axs[i][j].imshow(img_list[i * 2 + j], cmap='gray')
            axs[i][j].set_title(img_name[i * 2 + j])
            axs[i][j].axes.get_xaxis().set_visible(False)
            axs[i][j].axes.get_yaxis().set_visible(False)

    plt.show()


if __name__ == '__main__':
    main('man.jpg')

---

3、逆滤波与维纳滤波

原题：用下式对图像进行模糊处理，然后加高斯白噪声，得到降质图像。用逆滤波和维纳滤波恢复图像，对结果进行分析

题目中所要求的模糊处理的式子为
$$H(u,v)=e^{-k{\left(u^2+v^2\right)}^{5/6}}$$

模糊处理可以理解为我们在摄像的过程中，由于抖动导致图像在某个方向上变得模糊了，这个模糊反应在我们的式子中就是$H(u,v)$

而逆滤波的原理非常简单，怎么变过来就这么回去，F代表原图，G代表模糊后的图，N为噪声

可以看出来的是，在没有噪声的时候，逆滤波就可以很好的复原，但是如果有噪声，对原图像的影响还是很大的。为了应对这种问题，我们引入了维纳滤波

维纳滤波用来去除含有噪声的模糊图像，其目标是找到未污染图像的一个估计，使它们之间的均方差最小，可以去除噪声，同时清晰化模糊图像

具体原理可以参考《数字图像处理》或文末给出的连接

其对应的公式为
$$G(f)=\frac{1}{H(f)}\left[\frac{|H(f){|}^2}{|H(f){|}^2+\frac{N(f)}{S(f)}}\right]=\frac{1}{H(f)}\left[\frac{|H(f){|}^2}{|H(f){|}^2+\frac{1}{SNR(f)}}\right]$$

可见噪声为0的时候，等同与逆滤波

实现的效果如下所示

可见逆滤波在有噪声的时候不能很好的复原图像，但是维纳滤波可以，虽然在画质和边缘上有一些损失，但是还是很好的复原了图像

实现代码如下所示：

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image


def motion_process(height, weight, k=0.001):
    h_uv = np.zeros((height, weight), dtype='complex128')
    for u in range(height):
        for v in range(weight):
            q = np.power((u ** 2 + v ** 2), (5.0 / 6.0))
            h_uv[u][v] = np.exp(-(k * q))
    return h_uv


def blur_noise(src, h_uv):
    noise = np.random.normal(0, 1, size=src.shape)
    f_uv = np.fft.fft2(src)
    arr = np.multiply(f_uv, h_uv) + noise
    return arr, noise


def winner(src, h_uv, k=0.001):
    p_uv = np.conj(h_uv) / (np.abs(h_uv) ** 2 + 0.00005*k)
    rst = np.multiply(src, p_uv)
    return np.abs(np.fft.ifft2(rst))


def inv_filter(src, h_uv):
    p_uv = 1 / h_uv
    rst = np.multiply(src, p_uv)
    return np.abs(np.fft.ifft2(rst))


def main(img_path):
    src = np.array(Image.open(img_path).convert("L"))
    height, weight = src.shape
    h_uv = motion_process(height, weight, k=0.001)
    noise_src, noise = blur_noise(src, h_uv)
    mov_noi_src = np.abs(np.fft.ifft2(noise_src))
    k = np.power(np.abs(np.fft.fft2(noise)), 2) / np.power(np.abs(np.fft.fft2(src)), 2)

    src_processed_inv_filter = inv_filter(noise_src, h_uv)
    src_processed_winner = winner(noise_src, h_uv, k)

    img_list = [src, mov_noi_src, src_processed_inv_filter, src_processed_winner]
    img_name = ["原图像", "运动模糊+高斯噪声", "逆滤波", "维纳滤波"]

    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False

    _, axs = plt.subplots(2, 2, figsize=(12, 12))

    for i in range(2):
        for j in range(2):
            axs[i][j].imshow(img_list[i * 2 + j], cmap='gray')
            axs[i][j].set_title(img_name[i * 2 + j])
            axs[i][j].axes.get_xaxis().set_visible(False)
            axs[i][j].axes.get_yaxis().set_visible(False)

    plt.show()


if __name__ == '__main__':
    main('man.jpg')

---

4、维纳滤波和有约束最小二乘方滤波

原题：用运动模糊对图像进行模糊处理，然后加白高斯噪声，得到降质图像。用维纳滤波和有约束最小二乘方滤波恢复图像，约束线性算子用拉普拉斯算子。对实验结果进行对比分析

维纳滤波基本原理不赘述了

约束最小二乘方滤波核心是H对噪声的敏感性问题。减少噪声敏感新问题的一种方法是以平滑度量的最佳复原为基础的，建立约束式
$$C=\sum _0^{M-1}\sum _0^{N-1}{\left[{\nabla }^{2}f(x,y)\right]}^2$$

约束条件为
$$\Vert G-H\hat{F}{\Vert }_2^2=\Vert N{\Vert }_2^2$$

将上式表示成矩阵形式，同时将约束项转换成拉格朗日乘子项，并最小化代价函数

具体过程可以见《数字信号处理》或文末给出的连接

最后可得到
$$\hat{F}=\frac{\lambda H^{\ast }G}{\lambda H^{\ast }H+P^{\ast }P}=\left[\frac{H^{\ast }}{\Vert H{\Vert }_2^2+\lambda \Vert P{\Vert }_2^2}\right]G$$

P是函数
$$p=\left[\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right]$$

的傅里叶变换

实现的效果如下所示

可以看出在没有噪声的时候还是逆滤波的效果最好，但是在有噪声的时候，维纳滤波和约束最小二乘方滤波从两个角度给出了解答，虽然都具有一定的模糊，但是对图像都能够成功的修复，在我所选的因子下，维纳滤波的效果好一点，这个和调参也有关系

实现的代码如下所示：

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from numpy import fft


def motion_process(image_size, motion_angle):
    PSF = np.zeros(image_size)
    center_position = (image_size[0] - 1) / 2

    slope_tan = np.tan(motion_angle * np.pi / 180)
    slope_cot = 1 / slope_tan
    if slope_tan <= 1:
        for i in range(15):
            offset = round(i * slope_tan)  # ((center_position-i)*slope_tan)
            PSF[int(center_position + offset), int(center_position - offset)] = 1
        return PSF / PSF.sum()  # 对点扩散函数进行归一化亮度
    else:
        for i in range(15):
            offset = round(i * slope_cot)
            PSF[int(center_position - offset), int(center_position + offset)] = 1
        return PSF / PSF.sum()


def psf2otf(psf, outSize):
    psfSize = np.array(psf.shape)
    outSize = np.array(outSize)
    padSize = outSize - psfSize
    psf = np.pad(psf, ((0, padSize[0]), (0, padSize[1])), 'constant')
    for i in range(len(psfSize)):
        psf = np.roll(psf, -int(psfSize[i] / 2), i)
    otf = np.fft.fftn(psf)
    nElem = np.prod(psfSize)
    nOps = 0
    for k in range(len(psfSize)):
        nffts = nElem / psfSize[k]
        nOps = nOps + psfSize[k] * np.log2(psfSize[k]) * nffts
    if np.max(np.abs(np.imag(otf))) / np.max(np.abs(otf)) <= nOps * np.finfo(np.float32).eps:
        otf = np.real(otf)
    return otf


def CLSF(blurred, PSF, gamma=0.05):
    PF = psf2otf(PSF,blurred.shape)
    size = blurred.shape
    kernel = np.array([[0, -1, 0],
                       [-1, 4, -1],
                       [0, -1, 0]])

    PF_kernel = psf2otf(kernel, size)
    IF_noisy = fft.fft2(blurred)

    numerator = np.conj(PF)
    denominator = PF ** 2 + gamma * (PF_kernel ** 2)
    CLSF_deblurred = fft.ifft2(numerator * IF_noisy / denominator)
    CLSF_deblurred = np.real(CLSF_deblurred)
    return CLSF_deblurred


def inverse(blurred_image, PSF, eps):  # 逆滤波
    input_fft = fft.fft2(blurred_image)
    PSF_fft = fft.fft2(PSF) + eps  # 噪声功率，这是已知的，考虑epsilon
    result = fft.ifft2(input_fft / PSF_fft)  # 计算F(u,v)的傅里叶反变换
    result = np.abs(fft.fftshift(result))
    return result


def wiener(blurred_img, PSF, eps, k=0.01):  # 维纳滤波，K=0.01
    input_fft = fft.fft2(blurred_img)
    PSF_fft = fft.fft2(PSF) + eps
    PSF_fft_1 = np.conj(PSF_fft) / (np.abs(PSF_fft) ** 2 + k)
    result = fft.ifft2(input_fft * PSF_fft_1)
    result = np.abs(fft.fftshift(result))
    return result


def make_blurred(image, PSF, eps):
    input_fft = fft.fft2(image)  # 进行二维数组的傅里叶变换
    PSF_fft = fft.fft2(PSF) + eps
    blurred = fft.ifft2(input_fft * PSF_fft)
    blurred = np.abs(fft.fftshift(blurred))
    return blurred


def main(img_path):
    img = np.array(Image.open(img_path).convert("L"))

    img_size = img.shape
    PSF = motion_process(img_size, 60)

    blurred = make_blurred(img, PSF, 1e-3)
    blurred_inv = inverse(blurred, PSF, 1e-3)
    blurred_wiener = wiener(blurred, PSF, 1e-3, 0.01)
    blurred_CLSF = CLSF(blurred, PSF)

    blurred_noise = make_blurred(img, PSF, 1e-3) + np.random.normal(0, 5, img_size)
    blurred_noise_inv = inverse(blurred_noise, PSF, 1e-3)
    blurred_noise_wiener = wiener(blurred_noise, PSF, 1e-3, 0.01)
    blurred_noise_CLSF = CLSF(blurred_noise, PSF,0.03)

    img_list = [img, blurred, blurred_inv, blurred_wiener, blurred_CLSF,
                blurred_noise, blurred_noise_inv,blurred_noise_wiener, blurred_noise_CLSF]
    img_name = ["原图像", "运动模糊", "逆滤波", "维纳滤波", "约束最小二乘方",
                "运动模糊(噪声)", "逆滤波(噪声)", "维纳滤波(噪声)", "约束最小二乘方(噪声)"]

    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.rcParams['axes.unicode_minus'] = False

    _, axs = plt.subplots(3, 3, figsize=(12, 12))

    for i in range(3):
        for j in range(3):
            axs[i][j].imshow(img_list[i * 3 + j], cmap='gray')
            axs[i][j].set_title(img_name[i * 3 + j])
            axs[i][j].axes.get_xaxis().set_visible(False)
            axs[i][j].axes.get_yaxis().set_visible(False)

    plt.show()


if __name__ == '__main__':
    main('man.jpg')

参考链接

Python实现psf2otf_python

图像去模糊（约束最小二乘方滤波）

python实现逆滤波与维纳滤波客

Matlab 函数circShift、psf2otf 的 python（Numpy）实现

OpenCV—Python 图像去模糊（维纳滤波，约束最小二乘方滤波）