数字图像处理之图像复原

爱被打了一巴掌 2022-05-31 12:54 1105阅读 0赞

# 一、图像复原与图像增强的区别 #

图像的增强是一个主观的过程，其目的是改善图片的质量，对感兴趣的部分加以增强，对不感兴趣的部分予以抑制。而图像复原是一个客观的过程，针对质量降低或失真的图像，试图恢复其原始的内容或质量。复原技术是面向退化模型的，并且采用相反的过程进行处理，以便恢复出原图像。在进行图像复原之前要先建立起其退化模型，根据该模型进行图像复原。

课本中的图像退化过程建模为一个退化函数和一个加性噪声。

![70][]

假设H是一个线性时不变的过程，则我们可以得到![70 1][]，式子中的“\*”表示卷积。其频率域的表示为：![70 2][]。

对于只有加性噪声的情况，我们可以通过一些噪声模型（例如高斯噪声，瑞利噪声，椒盐噪声等等）以及对这些噪声参数的估计，来选择合适的空间滤波器（如均值滤波器，中值滤波器）或者频率滤波器（带阻/带通滤波器，低通/高通滤波器，陷波滤波器）来进行滤波。

# 二、退化函数的估计 #

在图像复原中，我们需要对退化函数进行估计。主要有观察法，实验法，数学建模法。

观察法通过选择噪声较小的子图像（减少噪声的影响）来得到H(u,v)，然后根据此信息来构建全图的H(u,v)，之后利用后面的复原方法来复原。实验法是指我们可以使用或设计一个与图像退化过程相似的装置（过程），使其成像一个脉冲，可得到退化系统的冲激响应 H(u,v) = G(u,v) / A。而建模估计则是从引起图像退化的基本原理进行推导，进而对原始图像进行模拟，在模拟过程中调整模型参数以获得尽可能精确的退化模型。课本中有两种模型，大气湍流模型和运动模糊模型。

## 大气湍流模型： ##

通用形式为![70 3][]。其中，k是与湍流性质有关的常数，k越大，图像越模糊，与高斯低通滤波器有着相同的形式。在实现该滤波器的过程中，由于中心化，要注意u,v应该分别替换为各自与频率中心之差，假设频率中心为(M/2,N/2)，则替换为u-M/2和v-N/2。与我们之前频率滤波器的实现相同。

其傅里叶频谱如下：

![70 4][]

注意：实现过程与我们之前的频率滤波器的实现一样，都需要中心化。

我在对图片施加大气湍流退化模型的时候，采用与之前相同的过程，但是得到的结果图像比课本的结果图更加地模糊。k = 0.0025（较模糊，比较容易看出区别）的情况如下：

![70 5][]

这个问题困扰了我很久，我尝试了各种不同的方法，在不中心变换的情况下进行滤波，发现情况更加糟糕。后来在做后面运动模糊模型的时候，得到的结果也与课本的结果相差甚远。之后，在尝试了不对图像进行填充（即不对原图像填充0值至大小为\[2m, 2n\]），发现得到的结果与课本的一致。对得到的结果图像直接取整（若是在最开始将值归一化到\[0,1\],自然不需要取整这一步），没有归一化到\[0, 255\]则与课本的结果一致。代码执行结果如下（同样是k = 0.0025的情况）：

![70 6][]

代码如下：

大气湍流实现函数为atmoTurbulence，输入为初始图像和大气湍流模型的系数k，输出为进行滤波之后的图像，下面给出的代码为了与课本的结果保持一致，没有在开始的时候对原图像进行填充，但结果仍然归一化到了\[0, 255\]。

function [image_out, g] = atmoTurbulence(image_in, k)
    % 模拟受大气湍流影响得到的图像
    % 输入为初始灰度图像，大气湍流模型的系数k
    % 输出为受到大气湍流影响的结果图像
    % 该实现方法没有对原图像进行填充。
    
    [m, n] = size(image_in);
    fp = zeros(m, n);
    % 中心化
    for i = 1 : m
        for j = 1 : n
            fp(i, j) = double(image_in(i,j)) * (-1)^(i+j);
        end
    end
    % 傅里叶变换
    F1 = fft2(fp);
    % 变换中心
    p = m / 2 + 1.0;
    q = n / 2 + 1.0;
    % 生成大气湍流模型函数
    H = zeros(m, n);
    for u = 1 : m
        for v = 1 : n
            temp1 = (u-p)^2+(v-q)^2;
            H(u,v) = exp((-k)*(temp1^(5/6)));
        end
    end
    % 滤波
    G = H .* F1;
    % 反傅里叶变换，并取实部
    gp = real(ifft2(G));
    % 反中心化
    
    g = zeros(m, n);
    for i = 1 : m
        for j = 1 : n
            g(i, j) = gp(i, j) * (-1)^(i+j);
        end
    end
    % 归一化输出图像到[0, 255],g取uint8的结果与课本一致
    image_out = zeros(m, n, 'uint8');
    mmax = max(g(:));
    mmin = min(g(:));
    range = mmax-mmin;
    for i = 1 : m
        for j = 1 : n
            image_out(i,j) = uint8(255 * (g(i, j)-mmin) / range);
        end
    end
    
    end

测试代码就不给出了，与之前的代码对函数的调用方法一致。

## 运动模糊模型： ##

退化函数为

![70 7][]

其中T表示曝光时间，a和b分别表示水平和垂直方向上的移动量。注意有![70 8][]，当π(ua+vb) = 0时，H(u,v) = T。同理，在实现该滤波器的过程中，由于中心化，要用u-M/2和v-N/2分别替换u和v。

当a = b  = 0.1，T = 1时，其傅里叶频谱如下：

imshow(log(abs(H) + 1), [])

![70 9][]

在对图像进行运动模糊退化之后，由于图像的运动，图像的尺寸与位置会发生一些的变化。在取a = b = 0.1，T = 1的情况下，不知道应该裁剪的数值，我直接显示了不进行裁剪的图像（由于对原图像填充0至图像大小为\[2m, 2n\]，在最后的结果需要裁剪），如下所示。但是，发现对图像进行裁剪之后，与课本给出的图像相差很大。

![70 10][]

课本给出的结果如下：

![70 11][]

经过对课本图片的仔细观察，发现图片的左边和上边应该输入右边和下边移动过来的。因此猜测课本的图片应该是没有对原图像进行填0扩充得到的。修改代码，得到的结果如下：

![70 12][]

代码如下：

clear all;
    clc;
    close all;
    
    image = imread('7.tif'); % 原图
    
    [m, n] = size(image);
    P = m / 2 + 1.0;
    Q = n / 2 + 1.0;
    fp = zeros(m, n);
    %对图像乘以(-1)^(x+y) 以移到变换中心
    for i = 1 : m
        for j = 1 : n
            % fp(i, j) = double(image(i, j));
            fp(i, j) = double(image(i, j)) * (-1)^(i+j);
        end
    end
    % 对填充后的图像进行傅里叶变换
    F1 = fft2(fp);
    
    a = 0.1;
    b = 0.1;
    T = 1;
    
    % 生成运动模型函数，中心在(P,Q)
    Mo = zeros(m, n);
    for u = 1 : m
        for v = 1 : n
            temp1 = pi * ((u-P)*a + (v-Q)*b);
            if (temp1 == 0)
                Mo(u,v) = T;
            else
                Mo(u, v) =  T * sin(temp1) / temp1 * exp(-1i * (temp1));
            end
        end
    end
    %进行滤波
    G = F1 .* Mo;
    
    % 反傅里叶变换
    gp = ifft2(G);
    
    % 处理得到的图像
    image_out = zeros(m, n, 'uint8');
    gp = real(gp);
    g = zeros(m, n);
    for i = 1 : m
        for j = 1 : n
            g(i, j) = gp(i, j) * (-1)^(i+j);
        end
    end
    mmax = max(g(:));
    mmin = min(g(:));
    range = mmax-mmin;
    for i = 1 : m
        for j = 1 : n
            image_out(i,j) = uint8(255 * (g(i, j)-mmin) / range);
        end
    end
    imshow(image_out);

# 三、逆滤波 #

## 1. 直接逆滤波 ##

知道了图像的退化模型之后，最简单的复原方法就是采用逆滤波的方法。计算如下：

![70 13][]

但是，当H(u,v)为0或者接近0的时候，N(u,v)/H(u,v)会变得很大，成为支配整个图像的主要部分。即使N(u,v)很小或者为0，当H(u,v)接近0的时候，在计算的时候，F(u,v)也会变得非常大，从而使复原之后的图像没有任何信息。

因此，可以对H(u,v)进行修改，①可以对G(u,v)/H(u,v)应用一个低通滤波器，虑去其中病态的高频成分（即虑去H(u,v)中接近0的部分）,②或者规定一个值，当|H(u,v)| ≤ δ 时，1/H(u,v) = 0。这两种方法都是去掉或者说削弱H(u,v)接近0时的影响。低通滤波器的方法对于上面提到的运动模糊模型几乎没有用，因为运动模糊模型的傅里叶频谱并不是从中间向四周赋值逐渐减小的，但是方法②仍然有用。

我在对收到了大气湍流模型之后的图像做直接逆滤波的时候，得到的最终图片一片模糊，得不到任何信息，与课本的结果相差较大，这可能是在对计算G(u,v)/H(u,v)的时候，在H(u,v)接近0的时候所做的处理不一致所导致的。这个直接造成了在之后对G(u,v)/H(u,v)使用butterworth低通滤波器的时候，最佳效果的截止频率与课本的不一致。

## 2. 维纳滤波（最小均方差误差滤波） ##

一个维纳滤波的计算公式如下，其中K是一个特定的常数，与噪声和未退化图像之间的信噪比有关。不同情况下取值不同。在平均意义下是最优的。

![70 14][]

## 3. 约束最小二乘方滤波器  ##

![70 15][]，其中r是一个参数，P(u,y)是矩阵p(x,y) = \[0, -1, 0; -1, 4,-1; 0, -1, 0\]的傅里叶变换。该算法相比于维纳滤波，对其应用的每幅图像都能产生最优的结果。

## 4. 代码实现 ##

下面给出在不填充图片的情况下，利用运动模型对图片进行模糊(a = b = 0.1, T = 1)，添加高斯噪声(均值为0，方差为500)，对有噪声的运动模糊图像直接使用逆滤波，使用维纳滤波（经过试验，k取0.02较好），使用约束最小二乘方滤波（经过试验，r取0.5较好）得到的结果。

代码如下（下面的代码再进行傅里叶变换之前没有对图像进行填充）：

clear all;
    close all;
    clc;
    
    image = imread('7.tif');
    [m, n] = size(image);
    % 参数如下：
    % p,q为频率中心，a,b,T为运动模糊参数，
    % k为维纳滤波参数， r为约束最小二乘方参数
    % M,V 分别为高斯噪声的均值和方差
    p = m / 2 + 1.0;
    q = n / 2 + 1.0;
    a = 0.1;
    b = 0.1;
    T = 1;
    k = 0.05;
    r = 0.5;
    M = 0;
    V = 500;
    
    % 对原图像进行傅里叶变换以及中心化
    fp = double(image);
    F = fftshift(fft2(fp));
    
    % 生成运动模型的傅里叶变换，中心在(p,q)
    Mo = zeros(m, n);
    for u = 1 : m
        for v = 1 : n
            temp = pi * ((u-p)*a + (v-q)*b);
            if (temp == 0)
                Mo(u,v) = T;
            else
                Mo(u, v) =  T * sin(temp) / temp * exp(-1i * (temp));
            end
        end
    end
    
    % 生成维纳滤波的傅里叶变换
    Wiener = (abs(Mo).^2) ./ (abs(Mo).^2 + k) ./ Mo;
    
    % 生成约束最小二乘方滤波的傅里叶变换
    lp = [0, -1, 0; -1, 4, -1; 0, -1, 0];
    Flp = fftshift(fft2(lp, m, n));
    Hw = conj(Mo) ./ (abs(Mo).^2 + r * abs(Flp).^2);
    
    % 生成高斯噪声的傅里叶变换
    noise = 500^0.5 * randn([m, n]);
    Fn = fftshift(fft2(noise));
    
    % 运动模糊图像，并且加上高斯噪声
    image1 = zeros(m, n, 'uint8');
    G1 = F .* Mo + Fn;
    gp1 = ifft2(fftshift(G1));
    g1 = real(gp1);
    % 归一化图像到 [0, 255];
    mmax = max(g1(:));
    mmin = min(g1(:));
    range = mmax-mmin;
    for i = 1 : m
        for j = 1 : n
            image1(i,j) = uint8(255 * (g1(i, j)-mmin) / range);
        end
    end
    
    % 为了接近真实情况，对归一化之后的加噪图像进行逆滤波
    F2 = fftshift(fft2(image1));
    % 直接逆滤波
    G2 = F2 ./ Mo;
    gp2 = ifftshift(G2);
    g2 = real(gp2);
    
    % 维纳滤波
    G3 = F2 .* Wiener;
    gp3 = ifft2(fftshift(G3));
    g3 = real(gp3);
    
    % 约束最小二乘方滤波
    G4 = F2 .* Hw;
    gp4 = ifft2(fftshift(G4));
    g4 = real(gp4);
    subplot(2,3,1), imshow(image), title('原图');
    subplot(2,3,2), imshow(image1), title('运动加噪图像');
    subplot(2,3,3), imshow(g2, []), title('直接逆滤波');
    subplot(2,3,4), imshow(g3, []), title('维纳滤波');
    subplot(2,3,5), imshow(g4, []), title('约束最小二乘方滤波');

执行结果如下：

![70 16][]

[70]: /images/20220531/3800d70953c143aab471cc683552e851.png
[70 1]: https://img-blog.csdn.net/20180225163316900?watermark/2/text/aHR0cDovL2Jsb2cuY3Nkbi5uZXQvY2pzaF8xMjM0NTY=/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70
[70 2]: /images/20220531/51c7e48c1e1e40e09a5aafb17fb1003e.png
[70 3]: /images/20220531/76ab0356060949039a59a7b78a4cd05f.png
[70 4]: /images/20220531/6620da39b35a4291aaa4f105c2b92709.png
[70 5]: /images/20220531/f000e9c0bd574cf988f4b2b3e443cfda.png
[70 6]: /images/20220531/74800e79b65946b286291b073640fc0c.png
[70 7]: /images/20220531/2197a1440c2c4dc584e5d55e90cb66e6.png
[70 8]: /images/20220531/b7196f36f866482b81d844149769a331.png
[70 9]: /images/20220531/d0401272c6e84bc59bc3a767f91c8400.png
[70 10]: /images/20220531/72c7fb8fc07e442096bd25e791419db4.png
[70 11]: /images/20220531/f6c1535b6a63492f8ef5e30ea459b9ae.png
[70 12]: /images/20220531/7396b473da4a424f8a23c0cf0f1d7286.png
[70 13]: /images/20220531/182580ee11cd4ffb910e02ba82fc34c8.png
[70 14]: /images/20220531/7af8b676797f4e90a788db16073dc786.png
[70 15]: /images/20220531/1a9bb6bb49ea43899694000c4d633779.png
[70 16]: /images/20220531/291285796b1448ccbe0e32c0b1483566.png