dct椒盐噪声去噪效果不好_matlab在空域与频域中去除周期噪声、椒盐噪声的简单应用...

前言

去年图像处理的DLL，有学弟问我做的思路，便放到博客里 github地址，欢迎star 图像增强处理：设计一套空间域与频率域结合的图像增强算法，处理以下任一组图片中的带噪声图像，去除噪声，提高图像质量。 (1)已知：噪声为随机噪声和周期噪声混合噪声； (2)要求： a)去噪处理后，计算均方误差评估去噪处理后图像的去噪效果 b)撰写完整的科技报告(形式类似科技论文)表述自己的算法设计，算法实现与算法评估过程。

对 swanNoise 图像去噪

空域去噪

根据老师讲解，swanNoise.bmp 所包含的噪声为椒盐噪声与周期噪声的混合。

对于传统图像中的椒盐噪声，适合使用 k 近邻滤波、中值滤波(二维统计滤波)、自适应中值滤波来去除噪声。k 近邻滤波能保留图像细节，使图像保持一定的清晰度，但椒盐噪声仍有少许干扰。中值滤波能完全去除椒盐噪声，但图像细节信息也损失了许多。

对于本图像，选取默认的参数调用三个滤波器去噪：

结果发现 k 近邻滤波图像细节损失不少。

再将二维统计滤波结果与自适应中值滤波结果比较：

结果发现适应中值滤波去噪效果最好。

频域去噪

将空域去噪的结果频谱图进行对比：

结果发现周期噪声突出的频率在每个区域均匀分布在各个点中。打开 photoshop 确认各点坐标。

对于频域滤波来说，可以采用高通滤波器、带阻滤波器、陷波滤波器、小波滤波器。下面对这几种滤波器进行比较：

高通滤波器对于周期噪声的滤波效果并不好

带阻滤波器对噪声以外的成分也有衰减

陷波滤波器对某个点进行衰减，对其余的成分不造成损失

小波滤波器对不同的频率成分分解到互不重叠的频带，对其余成分不造成损失

*由于试用版 matlab 无法安装 wavetoolbox ，在机房中实验的小波滤波器对本图去噪效果并不好，故采用陷波滤波器。

陷波滤波器结果如下：

小结

经过实验，发现先对频域去噪再对空域去噪效果最好，在这种情况下，自适应中值滤波比之原图的均方误差最小，结果如下：

二维统计滤波：

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;

对 dogDistorted 图像去噪

空域去噪

根据老师讲解，dogDistorted.bmp 所包含的噪声为高斯噪声与周期噪声的混合。

对于传统图像中的高斯噪声，适合使用非局部均值滤波、中值滤波(二维统计滤波)、三维块匹配滤波 非局部均值滤波均方误差更小，但三维块匹配滤波对于细节的保留程度更高。

频域去噪

同样的，采用陷波滤波器去噪，查看噪声图像频谱图：

采用一样的方法陷波滤波：

小结：