数字图像学笔记——13. 图像退化与复原（退化函数的评估方法：观察法、实验法、数学建模法与湍流导致的退化）

在对受到多种原因影响的图像进行复原时，我们经常需要先行评估对图像质量产生影响的退化函数，有时甚至需要尝试建模。通过这些手段，能够最大程度上恢复图像上的噪音，并重建高清的图像细节。

我们在前面的两篇文章里，已经分别介绍过对于空间噪音、频域噪音对于图像的影响。

数字图像学笔记——11. 图像退化与复原（空间噪音函数对图像的影响）
数字图像学笔记——12. 图像退化与复原（频域噪音函数对图像的影响、以及什么是陷滤波器）

在这里我们来详细讨论一下都有哪些传统的对退化的图像的评估方法。

线性位置不变的退化

在原书5.1节，曾提到过这样一个概念。图像的退化过程，可以表示为一个退化函数和一个加性噪音。

空间退化

其空间退化可以表示如下

$$g(x,y)=d(x,y)\odot f(x,y)+\eta (x,y)$$

其中 $\odot$ 表示卷积运算。

频域退化

对于频域来说也遵从类似形式：

$$G(u,v)=D(u,v)\ast F(u,v)+\mathrm{H}(u,v)$$

而在这里的 $\ast$ 表示矩阵的元素乘积。

在以上公式中，$h(x,y)$ 表示为二维噪音退化函数，$f(x,y)$ 表示为二维原始数据（这里指的是还未退化前的图像），自然 $\eta (x,y)$ 就是这里的加性噪音函数了。

所以，从上述公式中，我们可以看书，对于某个点 $(x_i,y_j)$ 它在退化前 $f(x_i,y_j)$ 的位置和退化后的位置 $g(x_i,y_j)$ 的位置并为发生改变，发生改变的只是点 $P(x_i,y_j)\to I_{(x_i,y_j)}$ 对应的像素值的变化， 这就是本章要特别强调的 线性位置不变 这一重要的概念与假设前提。

另外，除了由于采集设备因为电磁、温度干扰而导致的退化外，大部分采集得到的数据属于空间上的退化，所以以下内容将重点讨论在空间上的线性位置不变的退化。

一般评估方法

观察法

假设我们得到了一张退化后的图片 $g(x,y)$，基于线性位置不变的假设。我们通过一些手段（例如锐化、均值、中值滤波等）得到一张处理后的图片 $\hat{f}(x,y)$ （要求它比较接近还原效果），尽管这张图片可能不是特别理想的结果，但是我们可以假设噪音影响为0的前提下，通过这张处理后的图片来评估退化函数。

$$D(u,v)=\frac{G(u,v)}{\hat{F}(u,v)}$$

注意，为了更好的评估退化函数，我们通常会在这时把图片从空间转化到频域，如果存在退化函数，我们大概率能从频域图上看到带有某种光滑趋势曲线图，然后可以通过均值或平滑算法，得到退化函数在频域的估计$H(x,y)$，然后再通过卷积公式可以计算得到最终恢复的图像的频域值。

$$F(u,v)=D^{-1}(u,v)\odot G(u,v)\to f(x,y)$$

很显然，这种方法很繁琐，而且有一定的局限性，例如通常用来还原老照片。

试验法

现在，我们再来讨论第二种方法。

如果我们有条件可以得到和退化效果有相似的实现装置，从理论上讲，存在可以无限接近退化估计的可能。当我们的输入图像，经过这种退化实验装置后，得到的结果是我们想要恢复的退化图像的效果时。

就可以采用退化冲激响应来对退化函数进行采样了，对于冲激响应不是很理解的同学可以看看我写的这篇文章介绍 冲激函数——信号采样的利器。

在图像处理领域，使用的冲激响应模拟方法是激光，也就是拿激光照射实验设备，然后得到底片，通过分析底片的退化程度来分析采样。

上面左图是激光照射后应该形成的图片（未经过退化的结果），右侧是通过退化装置后照射得到的结果（经过退化的结果）。

冲激函数在数学上是一个脉冲点（在某一点上可能无限大，在其他位置上为0的脉冲），而如果在离散的傅立叶空间里，冲激所对应的就是一个常量A。所以我们就可以很容易的对退化函数使用以下的公式：

$$D(u,v)=\frac{G(u,v)}{A}$$

需要注意，这里的 $G(u,v)$ 已经不再是需要恢复的原始图像的FFT，而是这里实验的冲激退化后的FFT，而通常我们习惯上令 A = 1 （离散冲激函数的性质或特点），所以上式又可以简化为：

$$D(u,v)=G(u,v)$$

而且事实上也可以通过数学证明，对于矩阵来说，是否有倍数不影响矩阵的性质和结果。

数学建模法

如果没有条件搭建这样的实验设备和环境，那么我们可以尝试使用数学建模法来评估退化函数。以教材里提到的例子为例，Hufnagel and Stanley [1964] ¹ 最终给出的对大气湍流导致的图像退化其通式为：

$$D(u,v)=e^{-k(u^2+v^2{)}^{5/6}}$$

然后我们来实现这个代码看看。

def degradation_kernel(dft, k):
    # derive width, height, channel
    width, height, _ = dft.shape

    # center pointer
    p = width / 2 + 1.0
    q = height / 2 + 1.0

    # generate an empty kernel
    kernel = np.zeros((width, height), dtype=np.float32)

    # generate turbulence kernel
    for u in range(width):
        for v in range(height):
            power = -k * np.power((u - p) ** 2 + (v - q) ** 2, 5 / 6)
            kernel[u, v] = np.power(np.e, power)

    return kernel

这样我们就可以得到一个湍流退化函数了，接下来我们看看怎么让这个退化函数影响到我们的图片。

def update_dft_with_degradation(dft, kernel):

    # derive width, height, channel
    width, height, _ = dft.shape

    # shift dft
    dft_backup = np.fft.fftshift(dft)

    # apply the kernel
    dft_backup[:, :, 0] = dft_backup[:, :, 0] * kernel
    dft_backup[:, :, 1] = dft_backup[:, :, 1] * kernel

    # shift back
    dft_backup = np.fft.fftshift(dft_backup)

    return dft_backup


def turbulence_analysis(img, k):

    # convert byte to float
    dft = cv2.dft(np.float32(img), flags=cv2.DFT_COMPLEX_OUTPUT)

    # generate turbulence degradation
    kernel = degradation_kernel(dft, k)

    # apply kernel
    final_dft = update_dft_with_degradation(dft, kernel)

    # convert dft image back
    final_img = cv2.idft(final_dft, flags=cv2.DFT_COMPLEX_INPUT | cv2.DFT_SCALE | cv2.DFT_REAL_OUTPUT)

    # return
    return final_img

最终，我们输出如下：

---

源码可见：Github

---

1. Modulation Transfer Function Associated with Image Transmission through Turbulent Media https://www.osapublishing.org/josa/abstract.cfm?uri=josa-54-1-52 ↩︎