使用带阻滤波器去除周期性噪声

原图

结果图

处理过程

①傅里叶变换：对原图的三个通道进行傅里叶变换

②带阻滤波器：设计一个椭圆形的带阻滤波器

③滤波：对傅里叶变换的结果进行滤波

④逆傅里叶变换：对滤波后的结果进行逆傅里叶变换

⑤图像增强：使用二阶全变分算法对逆傅里叶变换的结构进行去噪，二阶全变分算法能够有效的平滑平坦区域，保留图像边缘细节。

MATLAB代码：

clear;
close all;
clc;
img_rgb = imread('figure.jpg');
figure,imshow(img_rgb,[]);

W=50;
n=5;
d0=112;
[l1,l2,l3]=size(img_rgb);
img_out = zeros(l1,l2,l3);

for k=1:3
    img=im2double(img_rgb(:,:,k));
    img_fft=fft2(img);
    g=fftshift(img_fft);
    
    [N1,N2]=size(g);
    Home = zeros(N1, N2);
    n1=round(N1/2);
    n2=round(N2/2);
    for i=1:N1 
        for j=1:N2
            d=sqrt(((i-n1)*N2/N1)^2+((j-n2))^2);
            Home(i,j)=1/(1+((d*W)/(d^2-d0^2))^(2*n));
        end
    end
    y=g.*Home;  
    if k==1
        %figure,imshow(g,[]);
        figure,imshow(log(abs(g)+1),[]);
        figure,imshow(Home,[]);
        figure,imshow(log(abs(y)+1),[]);
    end
    y=ifftshift(y);
    E1=ifft2(y);
    img_out(:,:,k)=real(E1);
end    
figure,imshow(img_out,[]), title('处理后图像');

%用一个TV_L2算子对处理结果进行去噪
img_out_filter = TV_L2_filter(img_out);
figure,imshow(img_out_filter,[]), title('处理后图像'); 
imwrite(img_out_filter,'.\result.jpg')

TV_L2_filter.m

function [S] = TV_L2_filter(Im, lambda, kappa, betamax)

if ~exist('kappa','var')
    kappa = 2.0;
end
if ~exist('lambda','var')
    lambda = 2e-2;
end
if ~exist('betamax','var')
    betamax = 1e5;
end
Im = im2double(Im);
S = Im;


fx = [1, -1];
fy = [1; -1];
[N,M,D] = size(Im);
sizeI2D = [N,M];
otfFx = psf2otf(fx,sizeI2D);
otfFy = psf2otf(fy,sizeI2D);
Normin1 = fft2(S);
Denormin2 = abs(otfFx).^2 + abs(otfFy ).^2;
if D>1
    Denormin2 = repmat(Denormin2,[1,1,D]);
end
beta = 2*lambda;
while beta < betamax
    lambeta = lambda/beta;
    Denormin   = 1 + beta*Denormin2;
    % h-v subproblem
    u = [diff(S,1,2), S(:,1,:) - S(:,end,:)];
    v = [diff(S,1,1); S(1,:,:) - S(end,:,:)];
    u = max(abs(u)-lambeta,0).*sign(u);
    v = max(abs(v)-lambeta,0).*sign(v);
    % S subproblem
    Normin2 = [u(:,end,:) - u(:, 1,:), -diff(u,1,2)];
    Normin2 = Normin2 + [v(end,:,:) - v(1, :,:); -diff(v,1,1)];
    FS = (Normin1 + beta*fft2(Normin2))./Denormin;
    S = real(ifft2(FS));
    beta = beta*kappa;
end
end