• 图像处理 -- 空域高斯滤波与频域高斯滤波


    卷积定理

    函数空间域的卷积的傅里叶变换是函数傅里叶变换的乘积。对应地,频率域的卷积与空间域的乘积存在对应关系。
    即:

        

        

    由卷积定理可知所有频域的滤波理论上都可以转化为空域的卷积操作。

    给定频率域滤波器,可对其进行傅里叶逆变换得到对应的空域滤波器;滤波在频域更为直观,但空域适合使用更

    小的滤波模板以提高滤波速度。因为相同尺寸下,频域滤波器效率高于空域滤波器,故空域滤波需要一个更小尺寸的模板近似得到需要的滤波结果。

    空域卷积

    将模板在图像中逐像素移动,将卷积核的每个元素分别和图像矩阵对应位置元素相乘并将结果累加,累加和作为模板中心对应像素点的卷积结果。通俗的讲,卷积就是对整幅图像进行加权平均的过程,每一个像素点的值,都由其本身和邻域内的其他像素值经过加权平均后得到。在像素的处理上,是先将结果暂存在于一个副本,最后统一拷贝,故不会出现处理顺序不同而结果不同的情况。

    二维连续卷积的数学定义:

     离散形式:

     频域滤波:

    频率域是由傅里叶变换和频率变量 (u,v)定义的空间,频域滤波处理过程:先对图像进行傅里叶变换,转换至频率域,在频域使用滤波函数进行滤波,最后将结果反变换至空间域。即:

     高斯函数:

     形状:

    高斯函数的特殊性:高斯函数傅里叶变换仍是高斯函数,但标准差已经变化,频域标准差越大(高斯函数越宽),变换后空域标准差越小(高斯函数越窄)。

    如频域函数对应的空域函数为

    空域高斯平滑滤波:

    高斯模板的生成

    因为图像是离散存储的,故我们需要一个高斯函数的离散近似。具体地,对高斯函数进行离散化,以离散点上的高斯函数值作为权值,组成一定尺寸的模板,用此模板对图像进行卷积。由于高斯分布在任意点处都非零,故理论上需要一个无穷大的模板,但根据"准则",即数据分布在的概率是0.9974,距离函数中心超过数据所占权重可以忽略,因此只需要计算的矩阵就可以保证对高斯函数的近似了。

    假设二维模板大小,则模板上元素 处的值为:

    前面的系数在实际应用中常被忽略,因为是离散取样,不能使取样和为1,最后还要做归一化操作。
    程序:

    function filt=mygaussian(varargin)
    %参数初始化,使用varargin处理可变参数情况
    siz=varargin{1};%模板尺寸
    if(numel(siz)==1)
        siz=[siz,siz];
    end
    std=varargin{2};%方差
    centa = (siz(1)+1)/2;%此处不要取整
    centb = (siz(1)+1)/2;
     
    filt = zeros(siz(1),siz(2));
    summ=0;
    for i=1:siz(1)
        for j=1:siz(2)
            radius = ((i-centa)^2+(j-centb)^2);
            filt(i,j) = exp(-(radius/(2*std^2)));
            summ=summ+filt(i,j);
        end
    end
    filt=filt/summ;%归一化

    测试:

    执行mygaussian(4,1)得:

    0.0181   0.0492    0.0492    0.0181

    0.0492   0.1336    0.1336    0.0492

    0.0492   0.1336    0.1336    0.0492

    0.0181   0.0492    0.0492    0.0181

    执行fspecial('gaussian',4,1)得:

    0.0181   0.0492    0.0492    0.0181

    0.0492   0.1336    0.1336    0.0492

    0.0492   0.1336    0.1336    0.0492

    0.0181   0.0492    0.0492    0.0181

    可以看出与Matlab结果相同。
    查看fspecial的Matlab源码,将gaussian部分提取出来:

    function h = hisfspecial(varargin)
    siz=varargin{1};
    if(numel(siz)==1)
        siz=[siz,siz];
    end
    std   = varargin{2};
    siz   = (siz-1)/2;
    [x,y] = meshgrid(-siz(2):siz(2),-siz(1):siz(1));
    arg   = -(x.*x + y.*y)/(2*std*std);
    h     = exp(arg);
    h(h<eps*max(h(:))) = 0;
    sumh = sum(h(:));
    if sumh ~= 0,
        h  = h/sumh;
    end;

    与mygaussian类似,只是更具Matlab特色。

    模板与图像卷积

    直接利用Matlab提供的fspecial、imfilter。

    I=imread('lena.bmp');
    h=fspecial('gaussian',7,8);
    O=imfilter(I,h,'replicate','same','conv');
    subplot(1,2,1);imshow(I,[]);title('input');
    subplot(1,2,2);imshow(O,[]);title('output');

    具体的实现细节还有很多值得研究的,如边界处理、将二维高斯模板分解为两个一维模板等。可参考文后链接。

    结果:

    频域高斯低通滤波

    在频率域,高斯函数表示为:

    或写成:

    其中为截止频率, 为图像频率空间点(u,v)处的值,;可以看出距离频率空间的(0,0)处越远,越大,也就是频率越高。而对高斯函数距离中心点 越远值越小,故利用可以滤除高频,保留低频。

    程序:

    I=imread('lena.bmp');
    I=double(I);
    [M,N]=size(I);
    m=(M+1)/2;
    n=(N+1)/2; 
    d0=30;                       %截止频率
    F=fftshift(fft2(I));               %空域转频域,平移中心
    H=zeros(M,N);
    G=zeros(M,N);
    for i=1:M
       for j=1:N
           d=(i-m)^2+(j-n)^2;
           H(i,j)=exp(-d/(2*d0^2));
           G(i,j)=H(i,j)*F(i,j);
        end
    end
    O=ifft2(ifftshift(G));                %平移中心,频域转空域
    O=real(O);                        %取实数部分
     
    subplot(1,2,1);imshow(I,[]);title('inpit');
    subplot(1,2,2);imshow(O,[]);title('output');
    结果:

    对比

    方差:空域高斯函数的方差越大,高斯函数越宽,模板尺寸越大,处理的图像越模糊;

               频域高斯函数的方差越小,高斯函数越窄,滤除的低频成分越多,图像越模糊;

    计算量:空域高斯滤波的计算花费随着模板的规模的增大而增大;频域高斯滤波的计算量独立于滤波函数。

    注:Matlab中提供了两个类似的空域滤波操作:卷积、相关,它们的不同之处是:卷积操作在计算前会先将模板旋转180度,下面会给出验证(对于我们的高斯模板没有影响,因为高斯模板是对称的)。

    验证:

    f=[0 0 0 0 0;

        0 0 0 0 0 ;

        0 0 1 0 0;

        0 0 0 0 0;

        0 0 0 0 0];

    g=[1 2 3 ;

         4 5 6 ;

         7 8 9];

    w1=imfilter(f,g,'conv');

    w2=imfilter(f,g,'corr');

    w1 =

        0     0     0    0     0

        0     1     2    3     0

        0     4     5    6     0

        0     7     8    9     0

        0     0     0    0     0

    w2 =

        0     0     0    0     0

        0     9     8    7     0

        0     6     5    4     0

        0     3     2    1     0

         0    0     0     0    0

    参考

    [1] http://blog.csdn.net/likezhaobin/article/details/6835049

    [2]http://blog.csdn.net/zddblog/article/details/7450033

    [3]http://homepages.inf.ed.ac.uk/rbf/HIPR2/freqfilt.htm

    [4]http://blog.csdn.net/zddblog/article/details/7521424


    原文链接:https://blog.csdn.net/u010839382/article/details/41908541

     

     

     

  • 相关阅读:
    树莓派4B踩坑指南
    树莓派4B踩坑指南
    树莓派4B踩坑指南
    从java反射到注解再到动态代理,一锅端!
    Lc_551学生出勤记录I
    aop日志测试类
    转载-使用@Async异步注解导致该Bean在循环依赖时启动报BeanCurrentlyInCreationException异常的根本原因分析,以及提供解决方案【享学Spring】
    linux 启动和停止脚本
    mysql 2003
    项目启动连不上mysql, Communications link failure
  • 原文地址:https://www.cnblogs.com/zzzsj/p/14648436.html
Copyright © 2020-2023  润新知