matlab在空域與頻域中去除週期噪聲、椒鹽噪聲的簡單應用

時間 2020-07-06

原文原文鏈接

前言

去年圖像處理的DLL，有學弟問我作的思路，便放到博客裏
github地址，歡迎star
圖像加強處理：設計一套空間域與頻率域結合的圖像加強算法，處理如下任一組圖片中的帶噪聲圖像，去除噪聲，提升圖像質量。
（1）已知：噪聲爲隨機噪聲和週期噪聲混合噪聲；
（2）要求：
a）去噪處理後，計算均方偏差評估去噪處理後圖像的去噪效果
b）撰寫完整的科技報告（形式相似科技論文）表述本身的算法設計，算法實現與算法評估過程。git

對 swanNoise 圖像去噪

空域去噪

根據老師講解，swanNoise.bmp 所包含的噪聲爲椒鹽噪聲與週期噪聲的混合。
對於傳統圖像中的椒鹽噪聲，適合使用 k 近鄰濾波、中值濾波(二維統計濾波)、自適應中值濾波來去除噪聲。k 近鄰濾波能保留圖像細節，使圖像保持必定的清晰度，但椒鹽噪聲仍有少量干擾。中值濾波能徹底去除椒鹽噪聲，但圖像細節信息也損失了許多。
對於本圖像，選取默認的參數調用三個濾波器去噪：

結果發現 k 近鄰濾波圖像細節損失很多。
再將二維統計濾波結果與自適應中值濾波結果比較：

結果發現適應中值濾波去噪效果最好。github

頻域去噪

將空域去噪的結果頻譜圖進行對比：

結果發現週期噪聲突出的頻率在每一個區域均勻分佈在各個點中。打開 photoshop 確認各點座標。
對於頻域濾波來講，能夠採用高通濾波器、帶阻濾波器、陷波濾波器、小波濾波器。下面對這幾種濾波器進行比較：算法

高通濾波器對於週期噪聲的濾波效果並很差
帶阻濾波器對噪聲之外的成分也有衰減
陷波濾波器對某個點進行衰減，對其他的成分不形成損失
小波濾波器對不一樣的頻率成分分解到互不重疊的頻帶，對其他成分不形成損失

*因爲試用版 matlab 沒法安裝 wavetoolbox ，在機房中實驗的小波濾波器對本圖去噪效果並很差，故採用陷波濾波器。
陷波濾波器結果以下：
設計

小結

通過實驗，發現先對頻域去噪再對空域去噪效果最好，在這種狀況下，自適應中值濾波比之原圖的均方偏差最小，結果以下：
3d

二維統計濾波：
function y = TwostaticFilter(imageWithNoise,k,boxSize) 
% iamgeWithNoise:噪聲圖像
% k:k值
% boxSize:模板尺寸
% 二維統計濾波

y = imageWithNoise;
[rows,cols]=size(y);
template = zeros(boxSize);
for i = 1:rows-boxSize+1
    for j = 1:cols-boxSize+1
        % 取模板內元素
        template = imageWithNoise(i:i+(boxSize-1),j:j+(boxSize-1));
        % 用排序後第k個值替換模板中心點像素值
        v = sort(template(:));
        y(i+(boxSize-1)/2,j+(boxSize-1)/2) = v(k);
    end
end
非局部均值濾波：
function DenoisedImg=NLmeans(I,ds,Ds,h)
 [m,n]=size(I);
DenoisedImg=zeros(m,n);
% 擴展圖像邊界
PaddedImg = padarray(I,[ds,ds],'symmetric','both');
% 定義d值
kernel=ones(2*ds+1,2*ds+1);
kernel=kernel./((2*ds+1)*(2*ds+1));
% 定義噪聲功率
h2=h*h;
for i=1:m
    for j=1:n
        % 原圖像對應擴展圖像的偏移量
        i1=i+ds;
        j1=j+ds;
        % 在擴展圖像中以(i1,j1)爲中心的鄰域窗口1
        W1=PaddedImg(i1-ds:i1+ds,j1-ds:j1+ds);
        average=0; % 加權和
        sweight=0; % 歸一化係數
        % 搜索窗口
        rmin = max(i1-Ds,ds+1); % 搜索窗口上邊界最低限制到原圖像上邊界
        rmax = min(i1+Ds,m+ds); % 搜索窗口下邊界最高限制到原圖像下邊界
        smin = max(j1-Ds,ds+1); % 搜索窗口左邊界最低限制到原圖像左邊界
        smax = min(j1+Ds,n+ds); % 搜索窗口右邊界最高限制到原圖像右邊界
        % r與s爲搜索窗口內像素點的座標,對搜索窗口內的每一個像素點求類似度
        for r=rmin:rmax
            for s=smin:smax
                % 不能與本身比較類似度
                if(r==i1&&s==j1) 
                    continue;
                end
                % 以搜索窗口內的像素點爲中心的鄰域窗口2
                W2=PaddedImg(r-ds:r+ds,s-ds:s+ds);
                % 計算鄰域間距離
                Dist2=sum(sum(kernel.*(W1-W2).*(W1-W2)));
                % 計算權值w(x,y)
                w=exp(-Dist2/h2);
                sweight=sweight+w;
                average=average+w*PaddedImg(r,s);
            end
        end
        % 將加權和歸一化並替換原像素點
        DenoisedImg(i,j)=average/sweight;
    end
end
三維塊匹配濾波：
function [y_est] = BM3D(y, z, sigma)
elseif strcmp(transform_type, 'dct') == 1,
    Tforward    = dct(eye(N));
elseif strcmp(transform_type, 'dst') == 1,
    Tforward    = dst(eye(N));
elseif strcmp(transform_type, 'DCrand') == 1,
    x = randn(N); x(1:end,1) = 1; [Q,R] = qr(x); 
    if (Q(1) < 0), 
        Q = -Q; 
    end;
    Tforward = Q';
else
    dwtmode('per','nodisp');
    Tforward = zeros(N,N);
    for i = 1:N
        Tforward(:,i)=wavedec(circshift([1 zeros(1,N-1)],[dec_levels i-1]), log2(N), transform_type);  %% construct transform matrix
    end
end
Tforward = (Tforward' * diag(sqrt(1./sum(Tforward.^2,2))))'; 
Tinverse = inv(Tforward);
return;

陷波濾波：
function [im,fftim] = swannotchfilter(Image,D)  
f = Image;
f = mat2gray(f,[0 255]);  
[M,N] = size(f);  
P = 2*M;  
Q = 2*N;  
fc = zeros(M,N);  
for x = 1:1:M  
    for y = 1:1:N  
        fc(x,y) = f(x,y) * (-1)^(x+y);  
    end  
end  
F = fft2(fc,P,Q);  
H_NF = ones(P,Q);  
for x = (-P/2):1:(P/2)-1  
     for y = (-Q/2):1:(Q/2)-1    
        v_k = 200; u_k = 145;  
        D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);  
        H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4);  
     end  
end  
G_1 = H_NF .* F;  
g_1 = real(ifft2(G_1));  
g_1 = g_1(1:1:M,1:1:N);       
for x = 1:1:M  
    for y = 1:1:N  
        g_1(x,y) = g_1(x,y) * (-1)^(x+y);  
    end  
end  
im=g_1;
fftim=G_1;