精华内容
下载资源
问答
  • %%======================= 1.3 谐波均值滤波器=============================== %定义子窗口的尺寸 m=5; n=5; %确定要扩展的行列数 len_m=floor(m/2); len_n=floor(n/2); %将原始图像进行扩展,这里采用了镜像扩展...
    %%%%%%------------------图像复原之空间滤波---------------------------------
    clc;
    clear;
    %读入图像,并转换为double型
    I=imread('D:\Gray Files\5-13.tif');
    I_D=im2double(I);
    [MM,NN]=size(I_D);
    
    %%%%%----------------------1、均值滤波器-----------------------------------
    %%======================= 1.3 谐波均值滤波器===============================
    %定义子窗口的尺寸
    m=5;
    n=5;
    %确定要扩展的行列数
    len_m=floor(m/2);
    len_n=floor(n/2);
    %将原始图像进行扩展,这里采用了镜像扩展,以进行图像边缘计算
    I_D_pad=padarray(I_D,[len_m,len_n],'symmetric');
    %获得扩展后的图像尺寸
    [M,N]=size(I_D_pad);
    J_Harmonic=zeros(MM,NN);
    %逐点计算子窗口的谐波平均
    for i=1+len_m:M-len_m
        for j=1+len_n:N-len_n
            %从扩展图像中取出子图像
            Block=I_D_pad(i-len_m:i+len_m,j-len_n:j+len_n);
            %求子窗口的谐波平均
            s=sum(sum(1./Block));
            J_Harmonic(i-len_m,j-len_n)=m*n/s;
        end
    end
    imshow(J_Harmonic);
    title('谐波均值滤波器');

     

    展开全文
  • %%======================= 1.4 逆谐波均值滤波器============================= %定义子窗口的尺寸 m=3; n=3; %确定要扩展的行列数 len_m=floor(m/2); len_n=floor(n/2); %将原始图像进行扩展,这里采用了镜像扩展...
    %%%%%%------------------图像复原之空间滤波---------------------------------
    clc;
    clear;
    %读入图像,并转换为double型
    I=imread('D:\Gray Files\5-13.tif');
    I_D=im2double(I);
    [MM,NN]=size(I_D);
    
    %%%%%----------------------1、均值滤波器-----------------------------------
    %%======================= 1.4 逆谐波均值滤波器=============================
    %定义子窗口的尺寸
    m=3;
    n=3;
    %确定要扩展的行列数
    len_m=floor(m/2);
    len_n=floor(n/2);
    %将原始图像进行扩展,这里采用了镜像扩展,以进行图像边缘计算
    I_D_pad=padarray(I_D,[len_m,len_n],'symmetric');
    %获得扩展后的图像尺寸
    [M,N]=size(I_D_pad);
    %滤波器阶数
    Q=1.5;
    J_Contraharmonic=zeros(MM,NN);
    %逐点计算子窗口的谐波平均
    for i=1+len_m:M-len_m
        for j=1+len_n:N-len_n
            %从扩展图像中取出子图像
            Block=I_D_pad(i-len_m:i+len_m,j-len_n:j+len_n);
            %求子窗口的谐波平均
            s1=sum(sum(Block.^(Q+1)));
            s2=sum(sum(Block.^Q));
            J_Contraharmonic(i-len_m,j-len_n)=s1/s2;
        end
    end
    imshow(J_Contraharmonic);
    title('逆谐波均值滤波器');

     

    展开全文
  • 算术均值 close all clear allf=imread('D:/testData/filtering.tif');[w,h]=size(f); image= f(:,:); fsize1=3; fsize2=5; fsize3=9;fssize1=(fsize1-1)/2; fssize2=(fsize2-1)/2; fssize3=(fsize3-1)/2;figure();...

    算术均值

    close all 
    clear all
    
    f=imread('D:/testData/filtering.tif');
    
    [w,h]=size(f);
    image= f(:,:);
    fsize1=3;
    fsize2=5;
    fsize3=9;
    
    fssize1=(fsize1-1)/2;
    fssize2=(fsize2-1)/2;
    fssize3=(fsize3-1)/2;
    
    figure();
    subplot(1,2,1);
    imshow(image);
    xlabel('原图像');
    
    resultImage = image;
    
    for x=1+fssize1:1:w-fssize1
        for y=1+fssize1:1:w-fssize1
            is=f(x-fssize1:1:x+fssize1,y-fssize1:1:y+fssize1);
            resultImage(x,y)=sum(is(:))/numel(is); 
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('3*3算术均值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize2:1:w-fssize2
        for y=1+fssize2:1:w-fssize2
            is=f(x-fssize2:1:x+fssize2,y-fssize2:1:y+fssize2);
            resultImage(x,y)=sum(is(:))/numel(is); 
    
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('5*5算术均值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize3:1:w-fssize3
        for y=1+fssize3:1:w-fssize3
            is=f(x-fssize3:1:x+fssize3,y-fssize3:1:y+fssize3);
            resultImage(x,y)=sum(is(:))/numel(is); 
    
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('9*9算术均值');

    几何均值

    close all 
    clear all
    
    f=imread('D:/testData/filtering.tif');
    
    [w,h]=size(f);
    image= f(:,:);
    fsize1=3;
    fsize2=5;
    fsize3=9;
    
    fssize1=(fsize1-1)/2;
    fssize2=(fsize2-1)/2;
    fssize3=(fsize3-1)/2;
    
    figure();
    subplot(1,2,1);
    imshow(image);
    xlabel('原图像');
    
    resultImage = image;
    
    for x=1+fssize1:1:w-fssize1
        for y=1+fssize1:1:w-fssize1
            is=f(x-fssize1:1:x+fssize1,y-fssize1:1:y+fssize1);
            resultImage(x,y)=prod(prod(is(:)))^(1/numel(is)); 
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('3*3几何均值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize2:1:w-fssize2
        for y=1+fssize2:1:w-fssize2
            %遍历每个点的四周
            is=f(x-fssize2:1:x+fssize2,y-fssize2:1:y+fssize2);
            resultImage(x,y)=prod(prod(is(:)))^(1/numel(is)); 
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('5*5几何均值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize3:1:w-fssize3
        for y=1+fssize3:1:w-fssize3
            %遍历每个点的四周
            is=f(x-fssize3:1:x+fssize3,y-fssize3:1:y+fssize3);
            resultImage(x,y)=prod(prod(is(:)))^(1/numel(is)); 
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('9*9几何均值');

    谐波均值

    close all 
    clear all
    
    f=imread('D:/testData/filtering.tif');
    
    [w,h]=size(f);
    image= f(:,:);
    fsize1=3;
    fsize2=5;
    fsize3=9;
    
    fssize1=(fsize1-1)/2;
    fssize2=(fsize2-1)/2;
    fssize3=(fsize3-1)/2;
    
    figure();
    subplot(1,2,1);
    imshow(image);
    xlabel('原图像');
    
    resultImage = image;
    
    for x=1+fssize1:1:w-fssize1
        for y=1+fssize1:1:w-fssize1
            is=f(x-fssize1:1:x+fssize1,y-fssize1:1:y+fssize1);
            is=1./is;
            resultImage(x,y)=numel(is)/sum(is(:)); 
    
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('3*3谐波均值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize2:1:w-fssize2
        for y=1+fssize2:1:w-fssize2
            %遍历每个点的四周
            is=f(x-fssize2:1:x+fssize2,y-fssize2:1:y+fssize2);
            is=1./is;
            resultImage(x,y)=numel(is)/sum(is(:));  
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('5*5谐波均值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize3:1:w-fssize3
        for y=1+fssize3:1:w-fssize3
            %遍历每个点的四周
            is=f(x-fssize3:1:x+fssize3,y-fssize3:1:y+fssize3);
            is=1./is;
            resultImage(x,y)=numel(is)/sum(is(:));   
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('9*9谐波均值');

    逆谐波

    close all 
    clear all
    
    f=imread('D:/testData/filtering.tif');
    
    [w,h]=size(f);
    image= f(:,:);
    Q1 = 1.5;
    fsize1=3;
    fsize2=5;
    fsize3=9;
    
    fssize1=(fsize1-1)/2;
    fssize2=(fsize2-1)/2;
    fssize3=(fsize3-1)/2;
    
    figure();
    subplot(1,2,1);
    imshow(image);
    xlabel('原图像');
    
    resultImage = image;
    
    for x=1+fssize1:1:w-fssize1
        for y=1+fssize1:1:w-fssize1
            is=f(x-fssize1:1:x+fssize1,y-fssize1:1:y+fssize1);
            resultImage(x,y) = sum(is(:).^(Q1+1))/sum(is(:).^(Q1));
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('3*3逆谐波均值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize2:1:w-fssize2
        for y=1+fssize2:1:w-fssize2
            %遍历每个点的四周
            is=f(x-fssize2:1:x+fssize2,y-fssize2:1:y+fssize2);
            resultImage(x,y) = sum(is(:).^(Q1+1))/sum(is(:).^(Q1));
    
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('5*5逆谐波均值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize3:1:w-fssize3
        for y=1+fssize3:1:w-fssize3
            %遍历每个点的四周
            is=f(x-fssize3:1:x+fssize3,y-fssize3:1:y+fssize3);
            resultImage(x,y) = sum(is(:).^(Q1+1))/sum(is(:).^(Q1));
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('9*9逆谐波均值');

    中值

    close all 
    clear all
    
    f=imread('D:/testData/filtering.tif');
    
    [w,h]=size(f);
    image= f(:,:);
    
    fsize1=3;
    fsize2=5;
    fsize3=9;
    
    fssize1=(fsize1-1)/2;
    fssize2=(fsize2-1)/2;
    fssize3=(fsize3-1)/2;
    
    figure();
    subplot(1,2,1);
    imshow(image);
    xlabel('原图像');
    
    resultImage = image;
    
    for x=1+fssize1:1:w-fssize1
        for y=1+fssize1:1:w-fssize1
            is=f(x-fssize1:1:x+fssize1,y-fssize1:1:y+fssize1);
            temp = sort(is(:));
            resultImage(x,y)= temp((numel(temp)-1)/2);
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('3*3中值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize2:1:w-fssize2
        for y=1+fssize2:1:w-fssize2
            %遍历每个点的四周
            is=f(x-fssize2:1:x+fssize2,y-fssize2:1:y+fssize2);
            temp = sort(is(:));
            resultImage(x,y)= temp((numel(temp)-1)/2);
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('5*5中值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize3:1:w-fssize3
        for y=1+fssize3:1:w-fssize3
            %遍历每个点的四周
            is=f(x-fssize3:1:x+fssize3,y-fssize3:1:y+fssize3);
            temp = sort(is(:));
            resultImage(x,y)= temp((numel(temp)-1)/2);
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('9*9中值');

    最大值

    close all 
    clear all
    
    f=imread('D:/testData/filtering.tif');
    
    [w,h]=size(f);
    image= f(:,:);
    
    fsize1=3;
    fsize2=5;
    fsize3=9;
    
    fssize1=(fsize1-1)/2;
    fssize2=(fsize2-1)/2;
    fssize3=(fsize3-1)/2;
    
    figure();
    subplot(1,2,1);
    imshow(image);
    xlabel('原图像');
    
    resultImage = image;
    
    for x=1+fssize1:1:w-fssize1
        for y=1+fssize1:1:w-fssize1
            is=f(x-fssize1:1:x+fssize1,y-fssize1:1:y+fssize1);
            temp = is(:);
            resultImage(x,y)= max(temp);
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('3*3最大值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize2:1:w-fssize2
        for y=1+fssize2:1:w-fssize2
            %遍历每个点的四周
            is=f(x-fssize2:1:x+fssize2,y-fssize2:1:y+fssize2);
            temp = is(:);
            resultImage(x,y)= max(temp);
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('5*5最大值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize3:1:w-fssize3
        for y=1+fssize3:1:w-fssize3
            %遍历每个点的四周
            is=f(x-fssize3:1:x+fssize3,y-fssize3:1:y+fssize3);
            temp = is(:);
            resultImage(x,y)= max(temp);
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('9*9最大值');

    最小值

    close all 
    clear all
    
    f=imread('D:/testData/filtering.tif');
    
    [w,h]=size(f);
    image= f(:,:);
    
    fsize1=3;
    fsize2=5;
    fsize3=9;
    
    fssize1=(fsize1-1)/2;
    fssize2=(fsize2-1)/2;
    fssize3=(fsize3-1)/2;
    
    figure();
    subplot(1,2,1);
    imshow(image);
    xlabel('原图像');
    
    resultImage = image;
    
    for x=1+fssize1:1:w-fssize1
        for y=1+fssize1:1:w-fssize1
            is=f(x-fssize1:1:x+fssize1,y-fssize1:1:y+fssize1);
            temp = is(:);
            resultImage(x,y)= min(temp);
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('3*3最小值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize2:1:w-fssize2
        for y=1+fssize2:1:w-fssize2
            %遍历每个点的四周
            is=f(x-fssize2:1:x+fssize2,y-fssize2:1:y+fssize2);
            temp = is(:);
            resultImage(x,y)= min(temp);
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('5*5最小值');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize3:1:w-fssize3
        for y=1+fssize3:1:w-fssize3
            %遍历每个点的四周
            is=f(x-fssize3:1:x+fssize3,y-fssize3:1:y+fssize3);
            temp = is(:);
            resultImage(x,y)= min(temp);
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('9*9最小值');

    中点

    close all 
    clear all
    
    f=imread('D:/testData/filtering.tif');
    
    [w,h]=size(f);
    image= f(:,:);
    
    fsize1=3;
    fsize2=5;
    fsize3=9;
    
    fssize1=(fsize1-1)/2;
    fssize2=(fsize2-1)/2;
    fssize3=(fsize3-1)/2;
    
    figure();
    subplot(1,2,1);
    imshow(image);
    xlabel('原图像');
    
    resultImage = image;
    
    for x=1+fssize1:1:w-fssize1
        for y=1+fssize1:1:w-fssize1
            is=f(x-fssize1:1:x+fssize1,y-fssize1:1:y+fssize1);
            temp = is(:);
            resultImage(x,y)= (max(temp) + min(temp))/2;
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('3*3中点');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize2:1:w-fssize2
        for y=1+fssize2:1:w-fssize2
            %遍历每个点的四周
            is=f(x-fssize2:1:x+fssize2,y-fssize2:1:y+fssize2);
            temp = is(:);
            resultImage(x,y)= (max(temp) + min(temp))/2;
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('5*5中点');
    
    resultImage= f(:,:);
    figure();
    subplot(1,2,1);
    imshow(f);
    xlabel('原图像');
    
    for x=1+fssize3:1:w-fssize3
        for y=1+fssize3:1:w-fssize3
            %遍历每个点的四周
            is=f(x-fssize3:1:x+fssize3,y-fssize3:1:y+fssize3);
            temp = is(:);
            resultImage(x,y)= (max(temp) + min(temp))/2;
        end
    end
    
    
    subplot(1,2,2);
    imshow(resultImage);
    xlabel('9*9中点');
    展开全文
  • 标题只存在噪声的复原 - 空间滤波均值滤波器算术平均滤波器几何均值滤波器谐波平均滤波器反(逆)谐波平均滤波器 只存在噪声的复原 - 空间滤波 仅被加性噪声退化 g(x,y)=f(x,y)+η(x,y)(5.21)g(x, y) = f(x, y) + \...

    只存在噪声的复原 - 空间滤波

    仅被加性噪声退化
    g(x,y)=f(x,y)+η(x,y)(5.21)g(x, y) = f(x, y) + \eta(x, y) \tag{5.21}

    G(u,v)=F(u,v)+N(u,v)(5.22)G(u, v) = F(u, v) + N(u, v) \tag{5.22}

    噪声项通常是未知的,因此不能直接从退化图像中减去噪声项来得到复原的图像。

    对于周期噪声,只是个例外而不是规律,我们可以用谱G(u,v)G(u, v)来估计N(u,v)N(u, v),从G(u,v)G(u, v)中减去N(u,v)N(u, v)能够得到原图像的一个估计。

    均值滤波器

    算术平均滤波器

    算术平均滤波器是最简单的均值滤波器。
    f^(x,y)=1mn(r,c)Sxyg(r,c)(5.23)\hat{f}(x, y) = \frac{1}{mn} \sum_{(r,c)\in S_{xy}} g(r,c) \tag{5.23}

    SxyS_{xy}表示中心为(x,y)(x, y)、大小为m×nm\times{n}的矩形子图像窗口(邻域)的一组坐标。算术平均滤波器在由SxyS_{xy}定义的区域中,计算被污染图像g(x,y)g(x, y)的平均值。复原的图像f^\hat{f}(x,y)(x, y)处的值,是使用SxyS_{xy}定义的邻域中的像素算出的算术平均值。

    均值滤波平滑图像中的局部变化,它会降低图像中的噪声,但会模糊图像

    def arithmentic_mean(image, kernel):
        """
        define arithmentic mean filter, math: $$\hat{f}(x, y) = \frac{1}{mn} \sum_{(r,c)\in S_{xy}} g(r,c)$$
        param image: input image
        param kerne: input kernel, actually use kernel shape
        return: image after arithmentic mean filter, 
        """
        img_h = image.shape[0]
        img_w = image.shape[1]
    
        m, n = kernel.shape[:2]
    
        padding_h = int((m -1)/2)
        padding_w = int((n -1)/2)
        
        # 这样的填充方式,可以奇数核或者偶数核都能正确填充
        image_pad = np.pad(image, ((padding_h, m - 1 - padding_h), \
                                      (padding_w, n - 1 - padding_w)), mode="edge")
    
        image_mean = image.copy()
        for i in range(padding_h, img_h + padding_h):
            for j in range(padding_w, img_w + padding_w):
                temp = np.sum(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1])
                image_mean[i - padding_h][j - padding_w] = 1/(m * n) * temp
        return image_mean
    
    # 算术平均滤波器
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接读为灰度图像
    
    mean_kernal = np.ones([3, 3])
    mean_kernal = mean_kernal / mean_kernal.size
    
    img_arithmentic = arithmentic_mean(img_ori, kernel=mean_kernal)
    
    img_cv2_mean = cv2.filter2D(img_ori, ddepth= -1, kernel=mean_kernal)
    
    plt.figure(figsize=(18, 6))
    plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(img_arithmentic, 'gray'), plt.title('Self Mean'), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(img_cv2_mean, 'gray'), plt.title('CV2 mean'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    几何均值滤波器

    f^(x,y)=[(r,c)Sxyg(r,c)]1mn(5.24)\hat{f}(x, y) = \Bigg[\prod_{(r,c)\in S_{xy}} g(r,c) \Bigg]^{\frac{1}{mn}} \tag{5.24}

    几何均值滤波器实现的平滑可与算术平均滤波器的相比,但损失的图像细节更少

    def geometric_mean(image, kernel):
        """
        define geometric mean filter, math: $$\hat{f}(x, y) = \Bigg[\prod_{(r,c)\in S_{xy}} g(r,c) \Bigg]^{\frac{1}{mn}}$$
        param image:  input image
        param kerne:  input kernel, actually use kernel shape
        return: image after geometric mean filter, 
        """
        img_h = image.shape[0]
        img_w = image.shape[1]
    
        m, n = kernel.shape[:2]
        order = 1 / (kernel.size)
    
        padding_h = int((m -1)/2)
        padding_w = int((n -1)/2)
        
        # 这样的填充方式,可以奇数核或者偶数核都能正确填充
        image_pad = np.pad(image.copy(), ((padding_h, m - 1 - padding_h), \
                                      (padding_w, n - 1 - padding_w)), mode="edge")
    
        image_mean = image.copy()
        # 这里要指定数据类型,但指定是uint64或者float64,但结果不正确,反而乘以1.0,也是float64,但却让结果正确
        for i in range(padding_h, img_h + padding_h):
            for j in range(padding_w, img_w + padding_w):
                prod = np.prod(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1]*1.0)
                image_mean[i - padding_h][j - padding_w] = np.power(prod, order)
    
        return image_mean
    
    def geometric_mean(image, kernel):
        """
        :param image: input image
        :param kernel: input kernel
        :return: image after convolution
        """
        
        img_h = image.shape[0]
        img_w = image.shape[1]
    
        kernel_h = kernel.shape[0]
        kernel_w = kernel.shape[1]
    
        # padding
        padding_h = int((kernel_h -1)/2)
        padding_w = int((kernel_w -1)/2)
    
        image_pad = np.pad(image.copy(), (padding_h, padding_w), mode="constant", constant_values=1)
    
        image_convol = image.copy()
        for i in range(padding_h, img_h + padding_h):
            for j in range(padding_w, img_w + padding_w):
                temp = np.prod(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1] * kernel)
                image_convol[i - padding_h][j - padding_w] = np.power(temp, 1/kernel.size)
    
        return image_convol
    
    # 几何均值滤波器
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接读为灰度图像
    
    mean_kernal = np.ones([3, 3])
    arithmetic_kernel = mean_kernal / mean_kernal.size
    
    img_geometric = geometric_mean(img_ori, kernel=mean_kernal)
    
    img_arithmentic = arithmentic_mean(img_ori, kernel=arithmetic_kernel)
    
    plt.figure(figsize=(18, 6))
    plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(img_arithmentic, 'gray'), plt.title('Arithmetic mean'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    # plt.show()
    

    在这里插入图片描述

    谐波平均滤波器

    f^(x,y)=mn(r,c)Sxy1g(r,c)(5.25)\hat{f}(x, y) = \cfrac{mn}{\sum_{(r,c)\in S_{xy}} \cfrac{1}{g(r,c)}} \tag{5.25}

    可以处理盐粒噪声,又能处理类似于高斯噪声的其他噪声,但不能处理胡椒噪声

    def harmonic_mean(image, kernel):
        """
        define harmonic mean filter, math: $$\hat{f}(x, y) = \Bigg[\prod_{(r,c)\in S_{xy}} g(r,c) \Bigg]^{\frac{1}{mn}}$$
        param image:  input image
        param kerne:  input kernel, actually use kernel shape
        return: image after harmonic mean filter, 
        """
        epsilon = 1e-8
        img_h = image.shape[0]
        img_w = image.shape[1]
    
        m, n = kernel.shape[:2]
        order = kernel.size
    
        padding_h = int((m -1)/2)
        padding_w = int((n -1)/2)
        
        # 这样的填充方式,可以奇数核或者偶数核都能正确填充
        image_pad = np.pad(image.copy(), ((padding_h, m - 1 - padding_h), \
                                      (padding_w, n - 1 - padding_w)), mode="edge")
    
        image_mean = image.copy()
        # 这里要指定数据类型,但指定是uint64或者float64,但结果不正确,反而乘以1.0,也是float64,但却让结果正确
        # 要加上epsilon,防止除0
        for i in range(padding_h, img_h + padding_h):
            for j in range(padding_w, img_w + padding_w):
                temp = np.sum(1 / (image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1]*1.0 + epsilon))
                image_mean[i - padding_h][j - padding_w] = order / temp
    
        return image_mean
    
    # 谐波滤波处理高斯噪声
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接读为灰度图像
    
    mean_kernal = np.ones([3, 3])
    arithmetic_kernel = mean_kernal / mean_kernal.size
    
    img_geometric = geometric_mean(img_ori, kernel=mean_kernal)
    
    img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernal)
    
    plt.figure(figsize=(18, 6))
    plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    # 谐波滤波处理盐粒噪声
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(b)(circuit-board-salt-prob-pt1).tif', 0) #直接读为灰度图像
    
    mean_kernal = np.ones([3, 3])
    arithmetic_kernel = mean_kernal / mean_kernal.size
    
    img_geometric = geometric_mean(img_ori, kernel=mean_kernal)
    
    img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernal)
    
    plt.figure(figsize=(18, 6))
    plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    # 谐波滤波处理胡椒噪声
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) #直接读为灰度图像
    
    mean_kernal = np.ones([3, 3])
    arithmetic_kernel = mean_kernal / mean_kernal.size
    
    img_geometric = geometric_mean(img_ori, kernel=mean_kernal)
    
    img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernal)
    
    plt.figure(figsize=(18, 6))
    plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    反(逆)谐波平均滤波器

    f^(x,y)=(r,c)Sxyg(r,c)Q+1(r,c)Sxyg(r,c)Q(5.26)\hat{f}(x, y) = \frac{\sum_{(r,c)\in S_{xy}} g(r,c)^{Q+1}}{\sum_{(r,c)\in S_{xy}} g(r,c)^Q} \tag{5.26}

    Q称为滤波器的阶数。这种滤波器适用于降低或消除椒盐噪声。Q值为正时,该滤波器消除胡椒噪声;Q值为负时,该滤波器消除盐粒噪声。然而,该滤波器不能同时消除这两种噪声。注意当Q=0Q=0时,简化为算术平均滤波器;当Q=1Q=-1时,简化为谐波平均滤波器。

    def inverse_harmonic_mean(image, kernel, Q=0):
        """
        define inverse harmonic mean filter, 
        math: $$\hat{f}(x, y) = \frac{\sum_{(r,c)\in S_{xy}} g(r,c)^{Q+1}}{\sum_{(r,c)\in S_{xy}} g(r,c)^Q}$$
        param image : input image
        param kernel: input kernel, actually use kernel shape
        param Q     : input order of the filter, default is 0, which equal to arithmentic mean filter, while is -1 is harmonic mean filter
        return: image after inverse harmonic mean filter, 
        """
        epsilon = 1e-8
        img_h = image.shape[0]
        img_w = image.shape[1]
    
        m, n = kernel.shape[:2]
    
        padding_h = int((m -1)/2)
        padding_w = int((n -1)/2)
        
        # 这样的填充方式,可以奇数核或者偶数核都能正确填充
        image_pad = np.pad(image.copy(), ((padding_h, m - 1 - padding_h), \
                                      (padding_w, n - 1 - padding_w)), mode="edge")
    
        image_mean = image.copy()
        # 这里要指定数据类型,但指定是uint64或者float64,但结果不正确,反而乘以1.0,也是float64,但却让结果正确
        # 要加上epsilon,防止除0
        for i in range(padding_h, img_h + padding_h):
            for j in range(padding_w, img_w + padding_w):
                temp = image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1] * 1.0 + epsilon
                # image_mean[i - padding_h][j - padding_w] = np.sum(temp**(Q+1)) / np.sum(temp**Q + epsilon)
                image_mean[i - padding_h][j - padding_w] = np.sum(np.power(temp, (Q+1))) / np.sum(np.power(temp, Q) + epsilon)
    
        return image_mean
    
    # 反谐波滤波处理胡椒噪声
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) #直接读为灰度图像
    
    mean_kernel = np.ones([3, 3])
    arithmetic_kernel = mean_kernel / mean_kernel.size
    
    img_inverse_harmonic = inverse_harmonic_mean(img_ori, kernel=mean_kernel, Q=1.5)
    
    img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernel)
    
    plt.figure(figsize=(18, 6))
    plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(img_inverse_harmonic, 'gray'), plt.title('Inverse Harmonic Mean'), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    # 反谐波滤波处理椒盐噪声
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0510(a)(ckt-board-saltpep-prob.pt05).tif', 0) #直接读为灰度图像
    
    mean_kernel = np.ones([3, 3])
    arithmetic_kernel = mean_kernel / mean_kernel.size
    
    img_inverse_harmonic = inverse_harmonic_mean(img_ori, kernel=mean_kernel, Q=1.5)
    img_arithmentic_mean = inverse_harmonic_mean(img_ori, kernel=mean_kernel, Q=0)
    # img_arithmentic_mean = arithmentic_mean(img_ori, kernel=mean_kernel)
    
    plt.figure(figsize=(18, 6))
    plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(img_inverse_harmonic, 'gray'), plt.title('Inverse Harmonic Mean'), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(img_arithmentic_mean, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    下面是各种滤波器的对比与总结

    # 算术平均滤波器和几何均值滤波器对比
    img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(a)(ckt-board-orig).tif', 0) #直接读为灰度图像
    img_noise = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接读为灰度图像
    
    mean_kernel = np.ones([3, 3])
    arithmetic_kernel = mean_kernel / mean_kernel.size
    
    img_arithmentic = arithmentic_mean(img_noise, kernel=mean_kernel)
    img_geometric   = geometric_mean(img_noise, kernel=mean_kernel)
    
    plt.figure(figsize=(12, 12))
    plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(222), plt.imshow(img_noise, 'gray'), plt.title('Noisy'), plt.xticks([]), plt.yticks([])
    plt.subplot(223), plt.imshow(img_arithmentic, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]), plt.yticks([])
    plt.subplot(224), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric mean'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    反谐波滤波器可以处理胡椒与盐粒噪声。在处理胡椒噪声的时候,Q值一般为正值;在处理盐粒噪声时,Q值一般为负值。错误的选择Q值,不仅不能得到好的滤波效果,反而带来灾难性的结果。

    # 反谐波滤波器分别使用不同的Q值对胡椒与盐粒噪声的滤波器
    img_pepper = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) 
    img_salt = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(b)(circuit-board-salt-prob-pt1).tif', 0)
    
    mean_kernel = np.ones([3, 3])
    arithmetic_kernel = mean_kernel / mean_kernel.size
    
    inverse_harmonic_15 = inverse_harmonic_mean(img_pepper, kernel=mean_kernel, Q=1.5)
    inverse_harmonic_15_ = inverse_harmonic_mean(img_salt, kernel=mean_kernel, Q=-1.5)
    
    plt.figure(figsize=(12, 12))
    plt.subplot(221), plt.imshow(img_pepper, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(222), plt.imshow(img_salt, 'gray'), plt.title('Noisy'), plt.xticks([]), plt.yticks([])
    plt.subplot(223), plt.imshow(inverse_harmonic_15, 'gray'), plt.title('Pepper Noise Q = 1.5'), plt.xticks([]), plt.yticks([])
    plt.subplot(224), plt.imshow(inverse_harmonic_15_, 'gray'), plt.title('Salt Noise Q = -1.5'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    # 错误的使用Q值,反谐波滤波器得到的是严重的后果
    img_pepper = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) 
    img_salt = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(b)(circuit-board-salt-prob-pt1).tif', 0)
    
    mean_kernel = np.ones([3, 3])
    arithmetic_kernel = mean_kernel / mean_kernel.size
    
    inverse_harmonic_15 = inverse_harmonic_mean(img_pepper, kernel=mean_kernel, Q=-1.5)
    inverse_harmonic_15_ = inverse_harmonic_mean(img_salt, kernel=mean_kernel, Q=1.5)
    
    plt.figure(figsize=(12, 12))
    plt.subplot(221), plt.imshow(img_pepper, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
    plt.subplot(222), plt.imshow(img_salt, 'gray'), plt.title('Noisy'), plt.xticks([]), plt.yticks([])
    plt.subplot(223), plt.imshow(inverse_harmonic_15, 'gray'), plt.title('Pepper Noise Q = -1.5'), plt.xticks([]), plt.yticks([])
    plt.subplot(224), plt.imshow(inverse_harmonic_15_, 'gray'), plt.title('Salt Noise Q = 1.5'), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    展开全文
  • 《数字图像处理》课后编程习题答案5.2-5.5:使用几何均值滤波器、谐波均值滤波器以及逆谐波均值滤波器重做习题5.1;资源包括四个问题的编程源码,以及获取问题所需的原始图像
  • 均值滤波器介绍均值滤波器属于低通滤波器;输出为模板内领域像素的简单平均值;主要用于图像的模糊和降噪,去除尖锐部分,比滤波器模板尺寸小的像素区域将会过滤掉;与此同时,边缘也会被平滑、模糊。线性均值滤波器...

空空如也

空空如也

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

谐波均值滤波器