精华内容
下载资源
问答
  • 数字图像处理图像复原

    万次阅读 多人点赞 2018-02-26 17:58:04
    一、图像复原与图像增强的区别图像的增强是一个主观的过程,其目的是改善图片的质量,对感兴趣的部分加以增强,对不感兴趣的部分予以抑制。而图像复原是一个客观的过程,针对质量降低或失真的图像,试图恢复其原始的...

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

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

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


    假设H是一个线性时不变的过程,则我们可以得到,式子中的“*”表示卷积。其频率域的表示为:

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

    二、退化函数的估计

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

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

    大气湍流模型:

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

    其傅里叶频谱如下:


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

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


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


    代码如下:

    大气湍流实现函数为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

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

    运动模糊模型:

    退化函数为

    其中T表示曝光时间,a和b分别表示水平和垂直方向上的移动量。注意有,当π(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), [])


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


    课本给出的结果如下:


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


    代码如下:

    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. 直接逆滤波

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

                                             

    但是,当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是一个特定的常数,与噪声和未退化图像之间的信噪比有关。不同情况下取值不同。在平均意义下是最优的。


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

    ,其中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('约束最小二乘方滤波');
    执行结果如下:


    展开全文
  • 数字图像处理第四章 图像复原图像复原与重建1 图像退化/复原过程的模型2 噪声模型2.1 用imnoise函数为图像添加噪声2.2 用给定分布产生空间随机噪声2.3 周期噪声3 仅有噪声的复原——空间滤波3.1 空间噪声滤波器3.2 ...

    图像复原与重建

    图像复原技术主要是以预先确定的目标来改善图像,与之前我们学习的图像增强相比,虽然有重叠之处,但是图像增强主要是一个主观过程,而图像复原大部分是一个客观过程。图像复原试图利用退化现象的某种先验知识来复原被退化的图像。因此,复原技术是面向退化模型的,并且采用相反的过程处理,以便恢复出原图像。

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

    在本章节中,退化过程被建模为一个退化函数和一个加性噪声项,对一幅输入图像f(x,y)进行处理,产生一幅退化后的图像g(x,y)。给定g(x,y)和关于退化函数H的一些知识以及关于退化函数H的一些知识以及关于加性噪声项η(x,y)的一些知识后,图像复原的目的就是获得原始图像的一个估计f*(x,y)。g(x,y)=H[f(x,y)]+η(x,y)g(x,y)=H[f(x,y)]+η(x,y)在这里插入图片描述
    本章节中使用的大部分复原方法都是以不同类型的图像复原滤波器为基础的。如果H是一个线性的、位置不变的过程,那么空间域中的退化图像可由下式给出:g(x,y)=h(x,y)f(x,y)+η(x,y)g(x,y)=h(x,y)★f(x,y)+η(x,y)式中,h(x,y)是退化函数的空间表示,被称为点扩散函数;★表示空间卷积,同时空间卷积等于频域乘积,那么我们把上式写成等价的频率域表示:G(u,v)=H(u,v)F(u,v)+N(u,v)G(u,v)=H(u,v)F(u,v)+N(u,v)式中,H(u,v)有时被称为光传递函数。

    2 噪声模型

    能够模仿噪声的行为和效果的能力是图像复原的核心。本章中,我们对两种基本噪声模型感兴趣;空间域中的噪声(用噪声概率密度函数来描述)以及频域中的噪声(用噪声的各种傅里叶特性来描述)。

    2.1 用imnoise函数为图像添加噪声

    图像处理工具箱采用imnoise函数,使噪声污染一幅图像,基本语法

    g = imnoise(f, type, parameters)
    

    函数imnoise在为图像加噪声之前,需要转换为范围在[0,1]的double类。

    下面让我们将均值为m、方差为var的高斯噪声加到图像上。

    >> f = imread('jimei2.jpg');
    >> d = im2double(f);
    >> g = imnoise(d, 'gaussian', 0.1, 0.05); //这里均值为0.1,方差为0.05
    >> subplot(1,2,1),imshow(d),title('Original Image');
    >> subplot(1,2,2),imshow(g,[]),title('Polluted By Gaussian');
    

    在这里插入图片描述
    除了高斯噪声,imnoise函数还可以添加椒盐噪声与泊松噪声。

    2.2 用给定分布产生空间随机噪声

    在函数imnoise中,能够产生可用类型和参数的噪声是很有必要的。空间噪声值是随机数,以概率密度函数(PDF)或者等价的、相应的累积分布函数(CDF)为特征。我们感兴趣的分布类型的随机数的产生要符合概率论的一些最简单规则。许多随机数产生器以区间(0,1)内具有均匀CDF的随机数产生问题为基础,我们就需要用其构造自己想要的随机变量。

    如果w是在区间(0,1)内均匀分布的随机变量,那么我们可以通过求解下面这个方程来得到具有规定的CDF和F的随机变量z:z=F1(w)z = F^{-1}(w)这个简单但却强有力的结果也可以等价地通过寻找方程F(z)=w的解来说明。下面是常用的随机变量以及他们的PDF:在这里插入图片描述
    接下来让我们用M-函数编程去定义自己的imnoise2函数可以生成上图中我们想要的随机变量

    function R=imnoise2(type,M,N,a,b)%type是函数类型,M*N是噪声数组的大小,a,b为两个参数
    %设置默认值
    if nargin==1%如果函数的输入参数为1,则默认a=0;b=1;M=1;N=1
        a=0;b=1;
        M=1;N=1;
    elseif nargin==3%如果函数的输入参数为3,则默认a=0;b=1
        a=0;b=1;
    end
    %开始运行程序
    switch lower(type)
        case 'gaussian'%如果是高斯类型,执行下面方程
            R=a+b*randn(M,N);
        case 'salt & pepper'%如果是焦盐类型,当输入参数小于等于3,a=0.05,b=0.05
            if nargin<=3
                a=0.05;b=0.05;
            end
            %检验Pa+Pb是否大于1
            if(a+b)>1
                error('The sum Pa+Pb must be not exceed >1')
            end
            R(1:M,1:N)=0.5;
            X=rand(M,N);%(0,1)范围内产生一个M*N大小的均匀随机数组
            c=find(X<=a);%寻找X中小于等于a的数,并赋值为0
            R(c)=0;
            u=a+b;
            c=find(X>a & X<=u);%寻找X中大于a并小于等于u的数,并赋值为1
            R(c)=1;
        case 'lognormal'%对数正态类型,当输入参数小于等于3,a=1,b=0.25,执行下面方程
            if nargin<=3
                a=1;b=0.25;
            end
            R=a*exp(b*randn(M,N));
        case 'rayleigh'%瑞利类型,执行下面方程
            R=a+(-b*log(1-rand(M,N))).^0.5;
        case 'exponential'%指数类型,执行下面方程
            if nargin<=3%如果输入参数小于等于3,a=1
                a=1;
            end
            if a<=0%如果a=0,错误类型
                error('Parameter a must be positive for exponential type.')
            end
            k=-1/a;
            R=k*log(1-rand(M,N));
        case 'erlang'%厄兰类型,如果输入参数小于等于3,a=2,b=5
            if nargin<=3
                a=2;b=5;
            end
          if(b~=round(b)|| b<=0)%如果b=0,错误类型
              error('Param b must a positive integer for Erlang.')
          end 
          k=-1/a;
          R=zeros(M,N);
          for j=1:b
              R=R+k*log(1-rand(M,N));
          end
        otherwise%如果不是以上类型,输出未知分配类型
            error('Unknown distribution type.')
    end
    

    有了imnoise2函数我们就可以得到我们自己想要类型的随机变量,让我们来生成一下他们的直方图观察一下分布:

    >> r = imnoise2('gaussian',10000,1,0,1);
    >> p = rand(10000,1);
    >> q = imnoise2('lognormal',10000,1,1,0.25);
    >> a = imnoise2('rayleigh',10000,1,0,1);
    >> b = imnoise2('exponential',10000,1,1);
    >> c = imnoise2('erlang',10000,1,2,5);
    >> subplot(2, 3, 1), hist(r,50), title('Gaussian');
    >> subplot(2, 3, 2), hist(p,50), title('Uniform');
    >> subplot(2, 3, 3), hist(q,50), title('Lognormal');
    >> subplot(2, 3, 4), hist(a,50), title('Rayleigh');
    >> subplot(2, 3, 5), hist(b,50), title('Exponential');
    >> subplot(2, 3, 6), hist(c,50), title('Erlang');
    

    在这里插入图片描述
    这里我们选取10000的数量级呈现出来的直方图,明显可以看出来每一个随机变量都分布都是符合其概率密度函数变化的。

    2.3 周期噪声

    一幅图像的周期噪声典型地产生于图像获取过程中的电气和电动机械的干扰。图像的周期噪声通常通过频率滤波来处理。周期噪声的模型是2D正弦波,他有如下所示的方程式:r(x,y)=Asin[2πu0(x+Bx)/M+2πv0(y+By)/N]r(x,y)=Asin[2\pi u_0(x+B_x)/M+2\pi v_0(y+B_y)/N]A是振幅,u0和v0分别确定了关于x轴和y轴的正弦频率。Bx和By是关于原点的相移。这个方程式的DFT是:R(u,v)=jAMN2[ej2π(u0Bx/M+v0By/N)δ(u+u0,v+v0)ej2π(u0Bx/M+v0By/N)δ(uu0,vv0)]R(u,v)=j\frac{AMN}{2}[e^{-j2\pi (u_0B_x/M+v_0B_y/N)}\delta(u+u_0,v+v_0)-e^{-j2\pi (u_0B_x/M+v_0B_y/N)}\delta(u-u_0,v-v_0)]我们可以根据DFT来编写imnoise3函数生成周期噪声

    function [r, R, S] = imnoise3(M, N, C, A, B)
    
    % 根据周期噪声的DFT公式得到,A是振幅,u0,v0是关于x轴和y轴确定正弦频率,这里是C表示
    % Bx和By是关于原点的相移。
    % [r,R,S] = IMNOISE3(M,N,C,A,B),产生大小为M-by-N的空间正弦噪声模式r,
    % 其傅里叶变换R和频谱S。 其余参数如下:
    % C是K-by-2矩阵,包含K对频域坐标(u,v), 表示频域中脉冲的位置。 
    % 这些位置相对于频率矩形中心(M / 2 + 1,N / 2 + 1)。 每个脉冲只需要一个坐标。 
    % 程序自动生成共轭对称脉冲的位置。 这些脉冲信号决定了r的频率成分。
    % A是1×k矢量,包含每个K脉冲对的幅度。 如果参数中不包含A,则使用的默认值为A = ONES(1,K)。 
    % 然后B自动设置为其默认值(参见下一段)。 为A(j)指定的值与C(j,12)中的坐标相关联。 
    % B是K-by-2矩阵,包含每个脉冲对的Bx和By相位分量。 B的默认值为B(1:K,12= 0% Process input parameters.
    [K, n] = size(C);
    if nargin ==3
        A(1:K) = 1.0;
        B(1:K, 1:2) = 0;
    elseif nargin == 4
        B(1:K, 1:2) = 0;
    end
    
    % Generate R
    R = zeros(M, N);
    for j = 1:K
        u1 = M/2 + 1 + C(j, 1);
        v1 = N/2 + 1 + C(j, 2);
        R(u1, v1) = i * (A(j)/2) * exp(i * 2 * pi * C(j, 1)* B(j,1)/M);  % 位于(u + u0, v + v0) 的冲击处的R
        % Complex conjugate.
        u2 = M/2 + 1 - C(j, 1);
        v2 = N/2 + 1 - C(j, 2);
        R(u2, v2) = -i * (A(j)/2) * exp(i * 2 * pi * C(j, 2)* B(j,2)/M);% 位于(u - u0, v - v0) 的冲击处的R
    end
    % Compute spectrum and spatial sinusoidal pattern.
    S = abs(R);
    r = real(ifft2(ifftshift(R)));
    

    我们可以指定傅里叶谱中的脉冲坐标,将其反变换为周期噪声。

    >> C = [0 64; 0 128; 32 32; 64 0; 128 0; -32 32];
    >> [r, R, S] = imnoise3(512,512,C);
    >> subplot(1,2,1),imshow(S,[]),title('Specified Pulse');
    >> subplot(1,2,2),imshow(r,[]),title('Spatial sinusoidal noise');
    

    在这里插入图片描述
    周期噪声在频域傅里叶谱中只是一个明亮的点,选择好滤波器在频域内很容易除去周期噪声。

    3 仅有噪声的复原——空间滤波

    3.1 空间噪声滤波器

    在这一节我们假设图像出现的退化仅仅是噪声,即退化遵循:g(x,y)=f(x,y)+η(x,y)g(x,y)=f(x,y)+\eta(x,y)在这种情况下,选择的降低噪声的方法是空间滤波。

    我们用M-函数去定义一个可以生成不同空间滤波器的函数:

    function f=spfilt(g,type,m,n,parameter)
    %spfilt执行线性和非线性的空间滤波器,g为原图像,type为滤波器类型,M*N为滤波器模板大小
    %处理输入参数
    if nargin==2
        m=3;n=3;Q=1.5;d=2;
    elseif nargin==5
        Q=parameter;d=parameter;
    elseif nargin==4
        Q=1.5;d=2;
    else
        error('Wrong number of inputs.');
    end
    %开始执行滤波
    switch type
        case 'amean'%算数平均滤波
            w=fspecial('average',[m n]);
            f=imfilter(g,w,'replicate');
        case 'gmean'%几何平均滤波
            [g, revertClass] = tofloat(g);
            f = exp(imfilter(log(g), ones(m,n), 'replicate')).^(1/m/n);
            f = revertClass(f);
        case 'hmean'%调和平均滤波
            [g, revertClass] = tofloat(g);
            f = m * n ./ imfilter(1./(g + eps), ones(m,n), 'replicate');
            f = revertClass(f);
        case 'median'%中值滤波
            f=medfilt2(g,[m n],'symmetric');
        case 'max'%最大值滤波
            f=ordfilt2(g,m*n,ones(m,n),'symmetric');
        case 'min'%最小值滤波
            f=ordfilt2(g,1,ones(m,n),'symmetric');
        case 'midpoint'%中值滤波
            f=ordfilt2(g,1,ones(m,n),'symmetric');
            f=ordfilt2(g,m*n,ones(m,n),'symmetric');
            f=imlincomb(0.5,f1,0.5,f2);
        otherwise
            error('Unkown filter type.')
    end
    

    下面我们使用spfilt去生成一个中值滤波器处理被椒盐噪声污染的图像。

    >> f = imread('jimei2.jpg');
    >> [M,N] = size(f);
    >> R = imnoise2('salt & pepper', M, N, 0.1, 0.1);
    >> g = f;
    >> g(R ==0) = 0;
    >> g(R ==1) = 255;
    >> g1 = rgb2gray(g);
    >> fs = spfilt(g1, 'median', 3, 3);
    >> subplot(1,3,1),imshow(g1),title('Polluted By Salt And Pepper');
    >> subplot(1,3,2),imshow(fs),title('Restore By Median');
    >> fs1 = spfilt(fs, 'median', 3, 3);
    >> subplot(1,3,3),imshow(fs1),title('Restore By Median Again');
    

    在这里插入图片描述
    这里结果就是简单的用中值滤波器处理了两次椒盐噪声污染的图像,反复处理最后结果会越来越清晰,空间滤波器应该已经很熟悉了。

    3.2 自适应空间滤波器

    上一节的空间滤波器我们在第二章已经很详细的学过了,但是在第二章中讨论用于图像处理的滤波器,并未考虑图像中的一点对其他点的特征变化,下面我们来学习两个简单的自适应滤波器,滤波器的特性变换以m×n矩形窗口Sxy定义的滤波器区域内图像的统计特性为基础。

    自适应局部降低噪声滤波器

    随机变量最简单的统计度量是其均值和方差。这也是我们在拿到一张图像以后最容易得到的两个属性,均值给出了在其上计算均值的区域中的平均灰度的度量,而方差则给出了该区域的对比度的度量。

    滤波器作用于局部区域Sxy。滤波器在该区域中心任意一点(x,y)上的响应基于以下4个量: (a) g(x,y),带噪图像在点(x,y)的值;
    (b) δ2η,污染f(x,y)以形成g(x,y)的噪声的方差;
    (c)mL,Sxy中像素的局部均值;
    (d)δ2L,Sxy中像素的局部方差。

    我们希望滤波器性能如下:
    1.如果 δ2η为零,则滤波器应该简单地返回g(x,y)的值。这无关紧要,在零噪声情况下g(x,y)等于f(x,y)。
    2.如果局部方差与δ2η是高度相关的,则滤波器返回g(x,y)的一个近似值。高局部方差通常与边缘相关,并且应该保护这些边缘。
    3.如果两个方差相等,我们则希望滤波器返回Sxy中的像素的算术均值。这种情况发生在局部区域与整个图像有相同特性的条件下,并且局部噪声将通过简单地求平均值来降低。

    基于这些假设得到的f*(x,y)的自适应表达式可以写成f^(x,y)=g(x,y)δη2δL2[g(x,y)mL]\hat{f}(x,y) = g(x,y)-\frac{\delta_\eta^2}{\delta_L^2}[g(x,y)-m_L]

    自适应中值滤波器
    对于上一小节使用的中值滤波器,只要脉冲噪声的空间密度不大,性能就会很好,这里我们要介绍的自适应中值滤波器可以平滑非脉冲噪声时会试图保留细节,这是传统中值滤波器所做不到的。我们令:zmin=Sxyzmax=Sxyzmed=Sxyzxy=(x,y)Smax=Sxyz_{min}=S_{xy}中的最小灰度值 \\ z_{max}=S_{xy}中的最大灰度值 \\ z_{med}=S_{xy}中的灰度值的中值 \\ z_{xy}=坐标(x,y)处的灰度值 \\ S_{max} = S_{xy}允许的最大尺寸自适应中值滤波算法以两个进程工作,分别表示为进程A和进程B,如下所示:
    进程AA1=zmedzminA2=zmedzmaxA1>0A2<0,BSmaxAzmedA_1=z_{med}-z_{min} \\ A_2 = z_{med}-z_{max} \\ 如果A_1>0且A_2<0,则转到进程B,否则增大窗口尺寸 \\ 如果窗口尺寸≤S_{max},则重复进程A,否则输出z_{med}进程BB1=zxyzminB2=zxyzmaxB1>0B2<0zxyzmedB_1 = z_{xy}-z_{min} \\ B_2 = z_{xy}-z_{max} \\ 如果B_1>0且B_2<0,则输出z_{xy},否则输出z_{med}理解该算法的关键在于,要记住它有3个主要目的:去除椒盐(脉冲)噪声,平滑其他非脉冲噪声,并减少诸如物体边界细化或粗化等失真。值zmin和zmax在算法统计上认为是类脉冲噪声成分,即使它们在图像中并不是最低最高的可能像素值,进程A的目的是确定中值滤波器的输出zmed是否是一个脉冲,进程B的目的是确定zmed非脉冲以后进行中值排序。下面我们来实现一下自适应中值滤波:

    function f = adpmedian(g,Smax)
    if (Smax<=1)||(Smax/2==round(Smax/2))||(Smax~=round(Smax))
        error('Smax must be an odd integer >1.')
    end
    f=g;
    f(:)=0;
    alreadyProcessed=false(size(g));
    for k=3:2:Smax
        zmin=ordfilt2(g,1,ones(k,k),'symmetric');
        zmax=ordfilt2(g,k*k,ones(k,k),'symmetric');
        zmed=medfilt2(g,[k k],'symmetric');
        
        processUsingLevelB=(zmed>zmin)&(zmax>zmed)&...
            ~alreadyProcessed;
        zB=(g>zmin)&(zmax>g);
        outputZxy=processUsingLevelB & zB;
        outputZmed=processUsingLevelB&~zB;
        f(outputZxy)=g(outputZxy);
        f(outputZmed)=zmed(outputZmed);
        
        alreadyProcessed=alreadyProcessed | processUsingLevelB;
        if all(alreadyProcessed(:))
            break;
        end
    end
    f(~alreadyProcessed)=zmed(~alreadyProcessed);
    end
    

    我们来对比一下传统滤波器和自适应滤波器的区别:

    >> w = imread('jimei2.jpg');
    >> f = im2double(rgb2gray(w));
    >> g = imnoise(f, 'salt & pepper', 0.25);
    >> f1 = medfilt2(g, [7,7], 'symmetric');
    >> f2 = adpmedian(g,7);
    >> subplot(1,3,1),imshow(g,[]),title('Polluted By Salt And Pepper');
    >> subplot(1,3,2),imshow(f1,[]),title('Restore By Traditional median filter');
    >> subplot(1,3,3),imshow(f2,[]),title('Restore By adaptive median filter');
    

    在这里插入图片描述
    我们可以看出来,传统的中值滤波器和自适应中值滤波器都很好的实现了去噪声的功能,但是明显可以看出传统滤波器处理后图像细节也被抹除了不少,边缘变得有点模糊,自适应中值滤波器就很好的保留了细节又消除了噪声。

    4 退化函数建模

    在图像复原问题中,遇到的主要退化是图像模糊。场景和传感器产生的模糊可以用空间域或频域的低通滤波器来建模。另一个重要的退化模型是由于在图像获取时传感器和场景之间的均匀线性运动而产生的图像模糊,下面我们对运动退化函数建模:

    >> f = checkerboard(8);
    >> PSF = fspecial('motion', 7, 45);
    >> gb = imfilter(f, PSF, 'circular');
    >> noise = imnoise(zeros(size(f)),'gaussian',0,0.001);
    >> g = gb + noise;
    >> subplot(2, 2, 1), imshow(pixeldup(f,8),[ ]);title('Chessboard');
    >> subplot(2, 2, 2), imshow(gb);title('Motion blurred image');
    >> subplot(2, 2, 3), imshow(noise);title('Noise');
    >> subplot(2, 2, 4), imshow(g);title('Noise And Motion');
    

    在这里插入图片描述
    我们将在下面两节去复原图像。

    5 维纳滤波

    维纳滤波是一种最早、也是最为熟知的线性图像复原方法。维纳滤波器寻找统计误差函数最小的估计f^:e2=E{(ff^)2}e^2=E\{(f-\hat{f})^2\}其中,E是期望值算子,f是未退化图像。这个表达式在频域中的表达为:F^(u,v)=[1H(u,v)H(u,v)2H(u,v)2+Sη(u,v)/Sf(u,v)]G(u,v)\hat{F}(u,v)=[\frac{1}{H(u,v)}\frac{|H(u,v)|^2}{|H(u,v)|^2+S_{\eta}(u,v)/S_f(u,v)}]G(u,v)比例式Sη(u,v)/Sf(u,v)被称为信噪功率比。我们感兴趣的两个量为平均噪声功率和平均图像功率,定义为:ηA=1MNuvSη(u,v)\eta_A=\frac{1}{MN}\sum_u\sum_vS_\eta(u,v)fA=1MNuvSf(u,v)f_A=\frac{1}{MN}\sum_u\sum_vS_f(u,v)其中,M和N分别表示图像和噪声数组的垂直和水平大小。这些量都是标量,它们的比例:R=ηAfAR=\frac{\eta_A}{f_A}也是标量,有时被用来代替函数Sη(u,v)/Sf(u,v)以产生常量数组。我们用维纳滤波去复原上一小节运动模糊的图像:

    >> fr1 = deconvwnr(g,PSF);
    >> Sn = abs(fft(noise)).^2;
    >> nA = sum(Sn(:))/prod(size(noise));
    >> Sf = abs(fft2(f)).^2;
    >> fA = sum(Sf(:))/prod(size(f));
    >> NSPR = nA/fA;
    >> fr2 = deconvwnr(g,PSF,NSPR);
    >> NACORR = fftshift(real(ifft2(Sn)));
    >> FACORR = fftshift(real(ifft2(Sf)));
    >> fr3 = deconvwnr(g,PSF,NACORR,FACORR);
    >> subplot(2,2,1),imshow(g,[]),title('Polluted Image');
    >> subplot(2,2,2),imshow(fr1,[]),title('Inverse filtering');
    >> subplot(2,2,3),imshow(fr2,[]),title('Wiener filter with constant ratio');
    >> subplot(2,2,4),imshow(fr3,[]),title('Wiener filtering of autocorrelation function');
    

    在这里插入图片描述
    从结果可以看出,简单的维纳滤波(逆滤波)后的结果仍然模糊,含有噪声,使用常数比率的维纳滤波后的结果虽然比逆滤波好多了但是仍有模糊,而使用自相关函数的维纳滤波后的结果虽然仍有噪声存在,但是最接近原始图像,实质上由于噪声的存在逆滤波就不能起到很好的效果,如果不加噪声光使用逆滤波就可以很好的复原了。

    6 约束的最小二乘法(规则化)滤波

    关于退化函数H必须具有一些先验知识对于本章讨论都所有方法都很常见,维纳滤波也存在一定问题:未退化图像和噪声的功率谱必须是已知的。本节给出的算法仅要求噪声方差和均值的知识,并且应用到的每幅图像都能产生最优结果。

    我们之前讨论的线性退化模型:g(x,y)=h(x,y)f(x,y)+η(x,y)g(x,y) = h(x,y)★f(x,y)+\eta(x,y)用向量矩阵的形式表达就是:g=Hf+ηg = Hf + \eta由此得出复原问题可以简化为简单的矩阵运算,这样的结论看起来是顺理成章的,遗憾的是我们在处理较大数量级的矩阵并不是一件很轻松的事情,这里不推导最小二乘法,该方法的核心是H对噪声的敏感性问题。减少噪声敏感性问题的一种方法是将复原的最优化性建立在平滑度量上,要想有意义,必须用手边问题的参数来约束复原,因此我们期望找到的是一个标准函数C的最小值,定义如下:C=x=0M1y=0N1[2f(x,y)]2C=\sum^{M-1}_{x=0}\sum^{N-1}_{y=0}[▽^2f(x,y)]^2其约束为gHf^2=η2||g-H\hat{f}||^2=||\eta||^2式中,||w||2=wTw是欧几里得范数,hat{f}是未退化图像的估计。拉普拉斯算子▽2在空域滤波器学过。我们下面用约束最小二乘方滤波去复原图像:

    >> f = checkerboard(8);
    >> PSF = fspecial('motion',7,45);
    >> gb = imfilter(f, PSF, 'circular');
    >> noise = imnoise(zeros(size(f)),'gaussian',0,0.001);
    >> g = gb + noise;
    >> fr1 = deconvreg(g,PSF,4);
    >> fr2 = deconvreg(g,PSF,0.4,[1e-7 1e7]);
    >> subplot(1,3,1),imshow(g,[]),title('Polluted Image');
    >> subplot(1,3,2),imshow(fr1,[]),title('Restored By NOISEPOWER');
    >> subplot(1,3,3),imshow(fr2,[]),title('Restored By NOISEPOWER Range [1e-7 1e7]');
    

    在这里插入图片描述
    结果图像与原始图像相比,稍微有些改善,还是没有使用自相关函数的维纳滤波后的复原效果好,但是该方法使用的先验知识比维纳滤波要少,可见我们掌握越多的先验知识就可以将结果做到越好。

    展开全文
  • 数字图像处理实验5图像复原

    千次阅读 2018-07-23 21:07:35
    一、实验目的 (1)了解图像复原的目的及意义,加深对图像复原理论的认识。 (2)掌握维纳滤波复原基本原理。 (3)掌握约束最小二乘方复原方法。...I=imread('E:\大三课件\大三下\数字图像处理\实验\...

    一、实验目的

    (1)了解图像复原的目的及意义,加深对图像复原理论的认识。

    (2)掌握维纳滤波复原基本原理。

    (3)掌握约束最小二乘方复原方法。

    (4)掌握盲解卷积复原方法

    二、实验内容

     (1)维纳滤波复原。

     (2)约束最小二乘方复原

     (3)盲解卷积复原

    三、实验代码及结果、分析

    (1)维纳滤波复原

    • 代码:

    I=imread('E:\大三课件\大三下\数字图像处理\实验\实验5\lena.jpg');

    I=rgb2gray(I);

    noise=5*randn(size(I));      

    noise = noise - min(min(noise));

    J  = double(I) +noise;

    R1=wiener2(J,[10 10]);    %未知噪声

    R2=wiener2(J,[10 10],noise);  %已知噪声分布

    figure      

    subplot(2,2,1),imshow(uint8(I));title('原始图像');      

    subplot(2,2,2),imshow(uint8(J));title('退化图像');      

    subplot(2,2,3),imshow(uint8(R1));title('盲复原');

    subplot(2,2,4),imshow(uint8(R2));title('非盲复原');

    • 结果:

    • 分析:
    1. 维纳滤波又叫最小均方差滤波,它的目标是找到一个原图像f的估计图像f',使得它们之间的均方误差最小。  e^2 = E((f-f')^2)
    2. randa()函数是产生随机噪声,nosize应该是一个矩阵,min(min(noise))表示先取每列中的最小值,再取最小值中的最小值
    3. 关于维纳滤波函数有wiener2和deconvwnr,分别适用于灰度图像和彩色图像
    4. R1为未知噪声的分布,R2为已知噪声分布

     

    (2)约束最小二乘方复原方法

    • 代码:

    I=imread('E:\大三课件\大三下\数字图像处理\实验\实验5\lena.jpg');

    I=rgb2gray(I);I = im2double(I);[hei,wid,~] = size(I);

    PSF = fspecial('motion', 21, 11);%摄像物体逆时针方向以11角度运动了21个像素

    blurred = imfilter(I, PSF, 'conv', 'circular');%运动模糊图像

    % 逆滤波结果

    If = fft2(blurred);Pf = psf2otf(PSF,[hei,wid]);deblurred =ifft2(If./Pf);

    % 运动模糊+高斯噪声

    noise_mean = 0;noise_var = 0.00001;

    blurred_noisy = imnoise(blurred, 'gaussian',noise_mean, noise_var);

    % 约束最小二乘法

    p = [0 -1 0;-1 4 -1;0 -1 0];P = psf2otf(p,[hei,wid]);

    gama = 0.002;%γ设置成0.001左右会有比较好的效果If = fft2(blurred_noisy);numerator = conj(Pf);

    denominator = Pf.^2 + gama*(P.^2);

    deblurred2 = ifft2( numerator.*If./ denominator );

    subplot(2,3,1),imshow(I);title('原始图像');

    subplot(2,3,2), imshow(blurred); title('运动模糊图像');

    subplot(2,3,3), imshow(deblurred); title('运动模糊复原')

    subplot(2,3,4), imshow(blurred_noisy),title('运动模糊+噪声图像')

    subplot(2,3,5), imshow(deblurred2),title('无约束最小二乘法复原');

    subplot(2,3,6); imshow(deconvreg(blurred_noisy, PSF,0)); title('deconvreg函数复原');%自带的去模糊deconvreg函数

    • 结果:

    • 分析:

    1、fspecial('type', len,theta)函数,用于预定于滤波算子,type=motion为运动算子,theta参数表示逆时针方向的角度,len表示运动了len个像素

    2、约束最小二乘方滤波效果比自带的去模糊deconvreg函数好

    3、约束最小二乘方滤波对高噪声和中等噪声产生结果要好于维纳滤波,对于低噪声两种滤波产生结果基本相等

     

    (3)盲解卷积复原方

    • 代码:

    I=checkerboard(8);

    PSF=fspecial('gaussian',7,10);

    V=.0001;

    BlurredNoisy=imnoise(imfilter(I,PSF),'gaussian',0,V);

    WT=zeros(size(I));

    WT(5:end-4,5:end-4)=1;

    INITPSF=ones(size(PSF));

    FUN=inline('PSF+P1','PSF','P1');

    [J P]=deconvblind(BlurredNoisy,INITPSF,20,10*sqrt(V),WT,FUN,0);

    subplot(131);imshow(I);title('原图像');

    subplot(132);imshow(BlurredNoisy);title('运动模糊+噪声图像');

    subplot(133);imshow(J);title('盲去卷积');

    • 结果:

    • 分析:
    1. 盲解卷积复原是指在没有图像退化先验知识、对退化系统了解不足的条件下,通过观察退化图像的多个图像以某种方式抽出退化信息,进行图像复原
    2. 函数imnoise 是表示添加噪声污染一幅图像,叫做噪声污染图像函数,g=imnoise(f,'localvar',V)将均值为0,局部方差为V的高斯噪声添加到图像f上,其中V是与f大小相同的一个数组,它包含了每一个点的理想方差值
    3. Inline()函数:定义一个内置函数,本质上说跟function干的是一样的事,只不过它可以直接内嵌在命令行里,不用另外单独定义function
    展开全文
  • 图像复原技术 图像在采集、传送、转换过程中会混入噪声,造成图像模糊、失真、有噪声,图像复原技术的目的是使图像尽可能恢复到原本的样子。 一、常见的噪声类型 高斯噪声 椒盐噪声 均匀分布噪声 指数分布噪声 伽马...

    图像复原技术

    图像在采集、传送、转换过程中会混入噪声,造成图像模糊、失真、有噪声,图像复原技术的目的是使图像尽可能恢复到原本的样子。

    一、常见的噪声类型

    • 高斯噪声
    • 椒盐噪声
    • 均匀分布噪声
    • 指数分布噪声
    • 伽马分布噪声

    在Matlab中,函数 imnoise 可以给图像加入高斯噪声、椒盐噪声、泊松分布噪声和乘性噪声

    二、空域内的滤波复原方法

    1、均值滤波

    (1)、算术均值滤波

    令 Sxy 表示中心在点(x, y)处、大小为 m * n的矩形子窗口区域(邻域)的一组坐标,然后计算这一组坐标所在像素值的平均值,将其赋值给点 (x, y)。

    通俗点讲,就是选图像上的一块区域,把这块区域的所有像素点相加求平均值再赋值给这块区域中心的那个像素点,后边的几种方法只是算法不同而已。

    公式为:f^(x,y)=1mn(s,t)Sxyg(s,t)\hat{f}(x, y) = \frac{1}{mn}\sum_{(s, t) \in S_{xy}}g(s, t)

    (2)、几何均值滤波

    算术均值滤波是把 Sxy 所代表点的像素值求和再取平均值,而几何均值滤波是把像素值求乘再求像素的数量次幂。几何均值滤波实现的平滑可与算术均值滤波相比,且丢失的图像细节更少
    公式为:f^(x,y)=[(s,t)Sxyg(s,t)]1mn\hat{f}(x, y) = \left[ \prod_{(s, t) \in S_{xy}}g(s, t)\right]^{\frac{1}{mn}}

    (3)、谐波均值滤波器

    谐波均值滤波器对于盐粒噪声效果很好,但不适用于胡椒噪声,擅于处理高斯噪声
    公式为:f^(x,y)=mn(s,t)Sxy1g(s,t)\hat{f}(x, y) = \frac{mn}{\sum_{(s, t) \in S_{xy}} \frac{1}{g(s, t)}}

    (4)、逆谐波均值滤波器

    Q 称为滤波器的阶数,适用于减少或在实际中消除椒盐噪声干扰的影响。当 Q 为正时,该滤波器可消除胡椒噪声,当 Q 为负值时该滤波器可消除盐粒噪声,当 Q 为 0 时,公式可化简为算术均值滤波,当 Q = -1 时则为谐波均值滤波器。
    公式为:f^(x,y)=(s,t)Sxyg(s,t)Q+1(s,t)Sxyg(s,t)Q\hat{f}(x, y) = \frac{\sum_{(s, t) \in S_{xy}}g(s, t)^{Q+1}}{\sum_{(s, t) \in S_{xy}}g(s, t)^{Q}}

    下列代码为算术均值滤波和几何均值滤波的实例代码及结果:

    I = imread('火影9.jpg');
    I = im2double(I);
    I = imnoise(I, 'gaussian', 0.05);   %加入高斯噪声
    %算术均值滤波
    PSF = fspecial('average', 3);   %创建预定义的二维滤波器
    J = imfilter(I, PSF);
    %几何均值滤波
    K = exp(imfilter(log(I), PSF));
    
    figure;
    subplot(131);imshow(I);title('加入高斯噪声');
    subplot(132);imshow(J);title('算术均值滤波');
    subplot(133);imshow(K);title('几何均值滤波');
    

    冈萨雷斯大佬的《数字图像处理》那本书上确实写了几何均值滤波比算术均值滤波要保留更多细节,我tm咋感觉是假的呢 o(一︿一+)o,结果如下:
    在这里插入图片描述

    2、统计排序滤波

    其实我觉得和均值滤波的原理差不多,都是选取一块区域,把这块区域里的像素们(一长串数字)进行各种运算操作。

    (1)、中值滤波器

    中值滤波函数:medfilt2

    正如其名字一样,使用一个像素邻域中的灰度级的中值来替代该像素的值。该滤波器针对单极或双极脉冲噪声尤其有效
    公式为:f^(x,y)=median(s,t)Sxy{g(s,t)}\hat{f}(x, y) = median_{(s, t) \in S_{xy}} \left\{ g(s, t) \right\}

    I = imread('火影9.jpg');
    I = im2double(I);
    I = imnoise(I, 'salt & pepper', 0.05);   %加入椒盐噪声
    
    J = medfilt2(I, [3, 3]);    %中值滤波
    
    figure;
    subplot(121);imshow(I);title('加入椒盐噪声');
    subplot(122);imshow(J);title('中值滤波');
    

    从结果可以看出,图像很好的去除了椒盐噪声,但是图像有些模糊了:
    在这里插入图片描述

    (2)、最大值滤波器

    就是对邻域中的所有像素点取最大值赋值给该像素点,最大值滤波器对于发现图像中的最亮点非常有用!另外它会从黑色物体边缘去除一些黑色像素
    公式为:f^(x,y)=max(s,t)Sxy{g(s,t)}\hat{f}(x, y) = max_{(s, t) \in S_{xy}} \left\{ g(s, t) \right\}

    (3)、最小值滤波器

    就是对邻域中的所有像素点取最小值赋值给该像素点,最小值滤波器对于发现图像中的最暗点非常有用!另外它会从白色物体边缘去除一些白色像素
    公式为:f^(x,y)=min(s,t)Sxy{g(s,t)}\hat{f}(x, y) = min_{(s, t) \in S_{xy}} \left\{ g(s, t) \right\}

    (4)、中点滤波器

    选取最大值与最小值的中点(就是最大最小值的平均)替换该点的像素值,结合了统计排序和求平均,对于处理随机噪声表现很好,比如高斯噪声或均匀噪声。
    公式为:f^(x,y)=12[max(s,t)Sxy{g(s,t)}+min(s,t)Sxy{g(s,t)}]\hat{f}(x, y) = \frac{1}{2} [ max_{(s, t) \in S_{xy}} \left\{ g(s, t) \right\} + min_{(s, t) \in S_{xy}} \left\{ g(s, t) \right\} ]

    • 二维顺序统计滤波函数:ordfilt2,可以实现中值、最大最小值滤波
    I = imread('火影9.jpg');
    I = rgb2gray(I);
    I = im2double(I);
    I = imnoise(I, 'salt & pepper', 0.1);   %加入椒盐噪声
    
    domain = [1 1 1 1; 1 1 1 1; 1 1 1 1; 1 1 1 1];
    J1 = ordfilt2(I, 1, domain);
    J2 = ordfilt2(I, 8, domain);
    J3 = ordfilt2(I, 16, domain);
    
    figure;
    subplot(221);imshow(I);title('加入椒盐噪声');
    subplot(222);imshow(J1);title('最小值滤波');
    subplot(223);imshow(J2);title('中值滤波');
    subplot(224);imshow(J3);title('最大值滤波');
    

    结果展示:
    在这里插入图片描述

    3、自适应滤波

    二维自适用去噪滤波:wiener2,自适应滤波器的性能比目前所有滤波器的性能都好

    I = imread('火影9.jpg');
    I = rgb2gray(I);
    
    %I = imcrop(I, [100, 100, 200, 200]);    %图像裁剪
    
    I = imnoise(I, 'gaussian', 0, 0.003);   %加入高斯噪声
    [K, noise] = wiener2(I, [5, 5]);
    
    figure;
    subplot(121);imshow(I);title('加入高斯噪声');
    subplot(122);imshow(K);title('自适应滤波');
    

    可以看出,效果还不错:
    在这里插入图片描述

    三、图像复原方法

    1、逆滤波

    大意就是我们知道了原图受到了什么算法(叫做退化函数 H )带来的噪声,那么我们从噪声图反着用这个算法就可以推算出原图。
    不过即使知道了退化函数,也还是不能完全恢复到原图。

    I = imread('cameraman.tif');
    I = im2double(I);
    [m, n] = size(I);
    M = 2 * m;
    N = 2 * n;
    
    u = -m / 2:(m / 2 - 1);
    v = -n / 2:(n / 2 - 1);
    [U, V] = meshgrid(u, v);           %基于向量 u 和 v 中包含的坐标返回二维网格坐标
    D = sqrt(U.^2 + V.^2);
    D0 = 130;                           %滤波器截止频率为20
    H = exp(-(D.^2)./(2 * (D0^2)));
    N = 0.01 * ones(size(I, 1), size(I, 2));
    N = imnoise(N, 'gaussian', 0, 0.001);
    J = fftfilter(I, H) + N;
    
    HC = zeros(m, n);
    M1 = H > 0.1;
    HC(M1) = 1 ./ H(M1);
    K = fftfilter(J, HC);
    HC = zeros(m, n);
    M2 = H > 0.01;
    HC(M2) = 1 ./ H(M2);
    L = fftfilter(J, HC);
    
    figure;
    subplot(221);imshow(J);
    subplot(222);imshow(J, []);
    subplot(223);imshow(K, []);
    subplot(224);imshow(L, []);
    

    其中,fftfilter函数如下:

    function Z = fftfilter(X, H)
    F = fft2(X, size(H, 1), size(H, 2));
    Z = H .* F;
    Z = ifftshift(Z);
    Z = abs(ifft2(Z));
    Z = Z(1:size(X, 1), 1:size(X, 2));
    end
    

    在这里插入图片描述

    完整目录

    Matlab数字图像处理——图像处理工具箱Image Processing Toolbox
    Matlab数字图像处理——图像类型的转换
    Matlab数字图像处理——图像文件的读取
    Matlab数字图像处理——图像文件的显示
    Matlab数字图像处理——视频文件的读写
    Matlab数字图像处理——图像的像素运算(灰度变换)
    Matlab数字图像处理——图像的空间变换
    Matlab数字图像处理——图像的平移、邻域操作、区域选取
    Matlab数字图像处理——图像增强
    Matlab数字图像处理——图像复原

    展开全文
  • 数字图像处理之傅里叶变换 by方阳 版权声明:本文为博主原创文章,转载请指明转载地址 http://www.cnblogs.com/fydeblog/p/7070055.html 1. 前言 这篇博客主要介绍常见的噪声及其概率密度函数,并用MATLAB...
  • 第五章 图像复原 一. 图像退化/复原的模型 图像增强:利用人的视觉系统特性,取得较好的视觉效果,不一定逼近原始图像 图像复原:针对图像退化原因进行补偿,要有一定的先验知识,利用逆过程复原图像,尽可能...
  • 图像复原与重建图像退化/复原过程的模型噪声模型噪声的空间和频率特性一些重要噪声的概率密度函数周期噪声噪声参数的估计仅存在噪声时的复原 图像退化/复原过程的模型 退化过程被建模为一个退化函数和一个加性噪声项...
  • 数字图像处理第四章数字图像处理---图像复原(一)图像退化/复原处理的模型(二)噪声模型2.1 用imnoise函数为图像添加噪声4.2 用给定分布产生空间随机噪声用均匀随机数来产生指定分布的随机数2.3 周期噪声2.4 估计...
  • 文章目录图像复原1. 图像退化/复原过程的模型2. 噪声模型2.1 周期噪声2.2 噪声参数估计3 只有噪声存在的空间域图像恢复 图像复原 与图像增强相似,图像复原的目的也是改善给定的图像,但图像增强主要是一个主观的...
  • 数字图像处理——图像退化与复原

    千次阅读 2019-09-28 18:49:56
    图像退化与复原的原理1.1 图像退化的数学模型1.2 图像退化的原理1.3 图像复原的原理2. 图像去噪2.1 噪声模型2.2 噪声参数的估计2.3 针对噪声的图像复原——空间、频域滤波 内容简介   图像复原和图像增强两者有...
  • matlab数字图像处理(六)图像复原

    千次阅读 2020-06-09 02:43:55
    图像复原 图像退化/复原过程的模型 图像复原与图像增强不同,图像复原必须逼近原始图像。 图像在形成、记录、处理和传输过程中,由于成像 系统、记录设备、传输介质和处理方法的不完善, 会导致图像质量下降。这一...
  • 数字图像处理图像复原

    千次阅读 2013-04-10 16:36:39
    图像增强是主观过程,图像复原大部分是一个客观过程。根据退化模型进行相反的处理,恢复出原图像。 1. 图像退化/复原的模型 空间域退化图像:线性移不变系统g(x) = h(x) * f(x) + n(x),h(x)为退化函数的空间描述...
  • 数字图像处理第四章数字图像处理---图像复原(三)仅有噪声的复原——空间滤波3.1 空间噪声滤波器4.3.2 自适应空间滤波 数字图像处理图像复原 (三)仅有噪声的复原——空间滤波 如果出现的退化仅仅是噪声,...
  • 数字图像处理(五) 图像复原

    千次阅读 2019-11-23 18:57:36
    本文首发于公众微信号-AI研究订阅号。   本节主要目的是介绍本书所用到的数字图像处理的一些基本概念。来源于东北大学 魏颖教授的模式识别课程笔记。 ...
  • 图像复原:根据图像畸变或退化的原因,进行模型化处理,将质量退化的图像重建或恢复到原始图像。即恢复退化图像的本来面目,忠于原图像,因此必须根据一定的图像退化模型来进行图像复原图像复原方法思路:...
  • 数字图像处理-MATLAB》运动模糊图像复原 图像复原技术也常被称为图像恢复技术图像复原技术能够去除或减轻在获取数字图像过程中发生的图像质量下降(退化)问题,从而使图像尽可能地接近于真实场景
  • 5.1 图像退化/复原过程的模型 退化:由于成像系统各种 因素的影响,使得图像质量降低 5.2 噪声模型 一些重要的噪声 高斯噪声(正态噪声)(gaussian) 瑞利噪声 伽马(爱尔兰)噪声 指数分布噪声 均匀分布噪声 脉冲...
  • 数字图像处理笔记—图像复原(一)0 概述定义:图像复原,即利用退化现象的某些先验知识来重建或复原被退化的图像。 与图像增强的区别:图像增强是一个主观的过程,即突出所关心的内容,满足人的视觉系统;图像复原...
  • 八、数字图像处理图像复原(1)

    千次阅读 2019-04-16 19:39:14
    图像复原是客观的处理,图像增强是主观的处理 (一)图像退化、复原处理的模型 如图所示,用退化函数对退化过程建模,它和附加噪声选项一起,作用于输入图像f(x,y),产生一幅退化的图像g(x,y): 如果H是线性的...
  • 第五章:图像复原和重建 1、 图像复原的目的是以预先确定的目标来改善图像; 2、 如果图像退化是因为有噪声,那么用空间域滤波就可以;如果是图像模糊,用频域滤波比较合适。 3、 噪声的主要来源是图像...
  • 数字图像处理(Digital Image Processing)是通过计算机对图像进行去除噪声、增强、复原、分割、提取特征等处理的方法和技术。本专栏将以学习笔记形式对数字图像处理的重点基础知识进行总结整理,欢迎大家一起学习交流...
  • 图像退化/复原过程的模型 退化:由于成像系统各种 因素的影响,使得图像质量降低 噪声模型 一些重要的噪声 高斯噪声 瑞利噪声 伽马(爱尔兰)噪声 指数分布噪声 均匀分布噪声 脉冲噪声(椒盐噪声) 几种噪声的...
  • 本文为武汉理工大学数字图像处理课程设计:彩色图像复原的实现过程,使用的语言为Python,如有不足之处欢迎指出。 先上效果图,对于jpg图片,在一分钟内处理完成,以下分别是原始图片,彩色图像,裁剪等优化后的...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 558
精华内容 223
关键字:

数字图像处理图像复原