目录
1.绪
图像复原目的：获得关于原始图像的近似估计f(x,y)。
由于种种原因(拍摄、传输、保存等)，图像变得模糊，或者被噪声污染，图像退化，退化的模型可以被模型化为一个退化函数和一个加性噪声项，这个时候想得到一个比较清晰完整的图像，就需要图像复原。
主要处理以噪声引起的退化。
2.噪声
2.1 一些重要噪声的概率密度函数
=====================================================================
=========================================================================
=========================================================================
2.2 添加噪声(实例1)
图1.1
对图1.1分别加入上面各种噪声，并显示出其直方图。
=============================MATLAB代码============================
注意将不是灰度图的转为灰度图，如果已经是灰度图就要将rgb2gray函数去掉，后面的代码同理。
如果想显示出比较高清的图片，可以使用imwrite(fig,'a.jpg');fig为要保存的图片，a.jpg为保存的名字和格式，可以直接保存。
直方图可以点击编辑窗口，进行参数范围，字体大小等的调整，使得呈现的图片更加饱满。
噪声添加，也可以使用imnoise(f,type,parameters);其中，f是输入图像，type选择噪声类型，parameters为可调整的参数(噪声的均值和方差等)。
clc;
clear;
%读入图像，并转换为double型
B=imread('huade.jpg');
I = rgb2gray(B); %转化为灰度图
I_D=im2double(I);
%获得图像大小
[M,N]=size(I_D);
%========================生成高斯白噪声==================================
a=0;
b=0.001;
N_Gau=a+sqrt(b)*randn(M,N);
%将噪声叠加到图像上
J_Gaussian=I_D+N_Gau;
%保存噪声数据
%======================生成瑞利噪声======================================
a=0;
b=0.08;
B=1;
N_Ray1=a+b*raylrnd(B,M,N);
%将噪声叠加到图像上
J_rayl=I_D+N_Ray1;
%=====================叠加伽马噪声=======================================
a=0;
b=0.08;
A=1;
B=2;
N_Gam=a+b*gamrnd(A,B,[M,N]);
%将噪声叠加到图像上
J_Gamma=I_D+N_Gam;
%====================叠加指数噪声=======================================
a=0;
b=0.08;
mu=2;
N_exp=a+b*exprnd(mu,[M,N]);
J_exp=I_D+N_exp;
%=====================叠加均匀分布噪声=======================================
a=0;
b=0.08;
A=0;
B=2;
N_unif=a+b*unifrnd(A,B,[M,N]);
J_unif=I_D+N_unif;
%=====================叠加椒盐分布噪声=======================================
a=0.25;
J_salt = imnoise(I_D,'salt & pepper',a);
subplot(2,3,1);
imshow(J_Gaussian);
title('高斯噪声');
subplot(2,3,2);
imshow(J_rayl);
title('瑞利噪声');
subplot(2,3,3);
imshow(J_Gamma);
title('伽马噪声');
subplot(2,3,4);
imshow(J_exp);
title('指数噪声');
subplot(2,3,5);
imshow(J_unif);
title('均匀分布噪声');
subplot(2,3,6);
imshow(J_salt);
title('椒盐噪声');
figure;
subplot(2,3,1);
imhist(J_Gaussian);
title('高斯噪声');
subplot(2,3,2);
imhist(J_rayl);
title('瑞利噪声');
subplot(2,3,3);
imhist(J_Gamma);
title('伽马噪声');
subplot(2,3,4);
imhist(J_exp);
title('指数噪声');
subplot(2,3,5);
imhist(J_unif);
title('均匀分布噪声');
subplot(2,3,6);
imhist(J_salt);
title('椒盐噪声');运行结果：
2.3 周期噪声
周期噪声是在图像获取时从电力或机电干扰中产生的。这是唯一的一种空间依赖型噪声。通常有正弦噪声和余弦噪声等。周期噪声通过频域滤波可以显著地减少。
2.4 添加正弦噪声(实例2)
图2
向图2加入正弦噪声，并画出其频谱图。
=============================MATLAB代码============================
clear;
clc;
%=======加入正弦噪声========
t=imread('einstein.jpg');
[m,n]=size(t);
t_1=t;
for i=1:m
for j=1:n
t_1(i,j)=t(i,j)+40*sin(40*i)+40*sin(40*j);
end
end
imshow(t),title('原图');
imwrite(t_1,'einstein1.jpg');
%===========谱=============
data = imread('einstein1.jpg');
data = rgb2gray(data);
F = fft2(data);
FC = fftshift(F);
S2 = log(1+abs(FC));
imshow(S2,[]);运行结果：
3.均值滤波
3.1 算术均值滤波
最简单的均值滤波器，m*n的矩形子图像窗口的各个值相加的和除以m*n。均值简单地平滑了一幅图像的局部变化，在模糊图像的同时减少了噪声。
3.2 几何均值滤波
每个被复原像素由子图像窗口中像素点的乘积并自乘到1/mn次幂。与算术均值相比，在滤波过程中会更少丢失图像细节。
3.3 均值滤波1(实例3)
图3
对图3加入均值为0，方差为0.02(归一化后)的高斯噪声，再分别用同样大小的窗口进行算术均值滤波和几何均值滤波。
=============================MATLAB代码============================
%高斯噪声部分
clc;clear;
I=imread('test1.jpg');
I=im2double(I);
%添加高斯噪声
I_noise=double(imnoise(I,'gaussian',0,0.02));%均值为0
imwrite(I_noise,'a.jpg');
%算术均值滤波
%定义子窗口的尺寸
m=3;
n=3;
%确定要扩展的行列数
len_m=floor(m/2);
len_n=floor(n/2);
%将原始图像进行扩展，这里采用了镜像扩展，以进行图像边缘计算
I_2=im2double(I_noise);
[M2,N2]=size(I_2);
I_2_pad=padarray(I_2,[len_m,len_n],'symmetric');
I_3=fspecial('average',[3,3]);%3*3均值滤波 建立预定义的滤波算子
I_3=imfilter(I_2_pad,I_3);%(待处理矩阵，滤波器)
imwrite(I_3,'b.jpg');
%几何均值滤波
I_D=im2double(I_noise);
[MM,NN]=size(I_D);
%定义子窗口的尺寸
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);
J_Geometric=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);
%求子窗口的几何均值
J_Geometric(i-len_m,j-len_n)=(prod(prod(Block)))^(1/(m*n));
end
end
imwrite(J_Geometric,'c.jpg');运行结果：
(a)原图 (b)高斯噪声污染
(c)3*3均值滤波 (d)3*3几何均值滤波
3.4 谐波均值滤波
谐波均值滤波器对“盐”噪声效果更好，但是不适用于“胡椒”噪声，它善于处理像高斯噪声那样的其他噪声。
3.5 逆谐波均值滤波
适合减少或在实际中消除椒盐噪声的影响。
Q称为滤波器的阶数。
当Q值为正数时，滤波器用于消除“胡椒”噪声。
当Q值为负数时，滤波器用于消除“盐粒”噪声。
当Q=0时，为算术均值滤波器。
当Q为-1时，为谐波均值滤波器。
3.6 均值滤波2(实例4)
对lena图分别加入概率为0.2的胡椒噪声和盐粒噪声，再用不同Q值得逆谐波滤波器分别对图像进行滤波。
=============================MATLAB代码============================
%椒盐噪声部分
clc;clear;
I=imread('Lena.jpg');
I=rgb2gray(I);
%imshow(I);
I=im2double(I);
%添加胡椒信号
pepper_ind = randperm(numel(I),0.1*numel(I));
tif_pepper = I;
tif_pepper(pepper_ind) = 0;
tif_pepper = im2double(tif_pepper);
figure;imshow(tif_pepper);
%Q=1.5逆谐波滤波器滤波对胡椒
Q=1.5;
I_mean=imfilter(tif_pepper.^(Q+1),fspecial('average',3))./imfilter(tif_pepper.^Q,fspecial('average',3));
figure;imshow(I_mean);
% Salt Noise
salt_ind = randperm(numel(I),0.1*numel(I));
tif_salt = I;
tif_salt(salt_ind) = 255;
tif_salt = im2double(tif_salt);
figure;imshow(tif_salt);
%Q=-1.5逆谐波滤波器滤波对盐粒
Q=-1.5;
I_mean=imfilter(tif_salt.^(Q+1),fspecial('average',3))./imfilter(tif_salt.^Q,fspecial('average',3));
figure;imshow(I_mean);运行结果：
(a)被概率为0.1的胡椒噪声污染; (b)被概率为0.1的盐粒噪声污染;
(c)使用大小为3*3，阶数为1.5的逆谐波均值对(a)滤波;
(d)使用大小为3*3，阶数为-1.5的逆谐波均值对(b)滤波;
如果Q值选取错误，即用正值的Q对“盐粒”进行滤波，用负值的Q对“胡椒”进行滤波，会导致更不好的效果。
=============================MATLAB代码============================
%椒盐噪声部分
clc;clear;
I=imread('Lena.jpg');
I=rgb2gray(I);
%imshow(I);
I=im2double(I);
%添加胡椒信号
pepper_ind = randperm(numel(I),0.1*numel(I));
tif_pepper = I;
tif_pepper(pepper_ind) = 0;
tif_pepper = im2double(tif_pepper);
figure;imshow(tif_pepper);
%Q=1.5逆谐波滤波器滤波对胡椒
Q=1.5;
I_mean=imfilter(tif_pepper.^(Q+1),fspecial('average',3))./imfilter(tif_pepper.^Q,fspecial('average',3));
figure;imshow(I_mean);
% Salt Noise
salt_ind = randperm(numel(I),0.1*numel(I));
tif_salt = I;
tif_salt(salt_ind) = 255;
tif_salt = im2double(tif_salt);
figure;imshow(tif_salt);
%Q=-1.5逆谐波滤波器滤波对盐粒
Q=-1.5;
I_mean=imfilter(tif_salt.^(Q+1),fspecial('average',3))./imfilter(tif_salt.^Q,fspecial('average',3));
figure;imshow(I_mean);运行结果：
(a)Q=-1.5对“胡椒”进行滤波 (b)Q=1.5对“盐粒”进行滤波
4.统计排序滤波
4.1 中值滤波
最著名的统计排序滤波器就是中值滤波器，就是用该像素的相邻像素的灰度中值来代替该像素的值。相同尺寸下，比线性平滑滤波器引起的模糊更少。对单极或双极脉冲噪声非常有效。
4.2 多次中值滤波(实例5)
对lena图加入0.2的椒盐噪声，再进行中值滤波。
=============================MATLAB代码============================
%椒盐噪声部分
clc;clear;
I=imread('test1.jpg');
I=im2double(I);
I_noise=double(imnoise(I,'salt & pepper',0.2));%salt & pepper注意中间的空格 无空格报错
imshow(I_noise);title('椒盐噪声');
%中值滤波
f1 = medfilt2(I_noise,'symmetric');%symmetric指出图像按照镜像反射方式对称地沿边界扩展
figure;imshow(f1);
%再次滤波
f2 = medfilt2(f1,'symmetric');%symmetric指出图像按照镜像反射方式对称地沿边界扩展
figure;imshow(f2);
%再一次滤波
f3 = medfilt2(f2,'symmetric');%symmetric指出图像按照镜像反射方式对称地沿边界扩展
figure;imshow(f3);运行结果：
(a)被概率为0.2的椒盐噪声污染 (b)使用3*3的中值滤波对(a)处理
(c)使用3*3的中值滤波对(b)处理 (d)使用3*3的中值滤波对(c)处理
4.3 最大值和最小值滤波
最大值：用该像素的相邻像素的灰度最大值来代替该像素的值，滤除“胡椒”比较好。
最小值：用该像素的相邻像素的灰度最小值来代替该像素的值，滤除“盐粒”比较好。
4.4 最大最小滤波(实例6)
向lena图分别加入0.5的胡椒噪声和盐粒噪声，然后进行最大值最小值滤波。
=============================MATLAB代码============================
%椒盐噪声部分
clc;clear;
I=imread('Lena.jpg');
I=rgb2gray(I);
I=im2double(I);
%添加胡椒信号
pepper_ind = randperm(numel(I),0.2*numel(I));
tif_pepper = I;
tif_pepper(pepper_ind) = 0;
tif_pepper = im2double(tif_pepper);
figure;imshow(tif_pepper);
% Salt Noise
salt_ind = randperm(numel(I),0.2*numel(I));
tif_salt = I;
tif_salt(salt_ind) = 255;
tif_salt = im2double(tif_salt);
figure;imshow(tif_salt);
%使用最大值
f3 = ordfilt2(tif_pepper,9,ones(3,3));
figure;imshow(f3);
%使用最小值
f4 = ordfilt2(tif_salt,1,ones(3,3));
figure;imshow(f4);运行结果：
4.5 修正后的阿尔法均值滤波器(实例7)
=============================MATLAB代码============================
clc;clear;
I=imread('eee.jpg');
I=rgb2gray(I);
I=im2double(I);
%获得图像大小
[M,N]=size(I);
% a=0;
% b=0.08;
% A=-1;
% B=1;
% N_unif=a+b*unifrnd(A,B,[M,N]);
% J_unif=I_D+N_unif;
% imshow(J_unif);
sigma = 0.1;
a = sqrt(3)*sigma;
N_unif = I + 2*(rand(M,N)-.5)*a;
imshow(N_unif);
%0.01的椒盐噪声
J_salt = imnoise(N_unif,'salt & pepper',0.01);
figure;imshow(J_salt);
%算术均值滤波
I_1=fspecial('average',[7,7]);
I_1=imfilter(J_salt,I_1);
figure;imshow(I_1);
%几何均值滤波
I_D=im2double(J_salt);
[MM,NN]=size(I_D);
%定义子窗口的尺寸
m=7;
n=7;
%确定要扩展的行列数
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_Geometric=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);
%求子窗口的几何均值
J_Geometric(i-len_m,j-len_n)=(prod(prod(Block)))^(1/(m*n));
end
end
figure;imshow(J_Geometric);
%中值滤波
I_3 = medfilt2(J_salt,'symmetric');%symmetric指出图像按照镜像反射方式对称地沿边界扩展
figure;imshow(I_3);
%阿尔法
d=7;
J_Alpha=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);
%计算矩阵的阿尔法均值
J_Alpha(i-len_m,j-len_n)=sum(sum(Block))/(m*n-d);
end
end
figure;imshow(J_Alpha);运行结果：
(a)被0.1均匀噪声污染 (b)再被0.01椒盐噪声污染
(c)5*5均值滤波 (d)5*5几何均值滤波
(e)5*5中值滤波 (f)5*5阿尔法均值滤波
4.6 自适应降噪滤波(实例8)
=============================MATLAB代码============================
clc;
clear;
%读入图像，并转换为double型
I=imread('I_noise.jpg');
imshow('I_noise.jpg');
%I=rgb2gray(I);
I_D=im2double(I);
[MM,NN]=size(I_D);
%======================= 自适应噪声削减滤波器=========================
%定义子窗口的尺寸
m=7;
n=7;
%确定要扩展的行列数
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);
%加性高斯白噪声的方差
var_noise=0.015;
J_AdaptiveReduction=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);
%将多维矩阵转换为一维数组
Block=Block(:);
%求局部图像的方差
var_block=var(Block);
%求局部图像的均值
mb=mean(Block);
if var_noise>var_block
v=1;
else
v=var_noise/var_block;
end
%叠加到原始图像上
J_AdaptiveReduction(i-len_m,j-len_n)=I_D(i-len_m,j-len_n)-v*(I_D(i-len_m,j-len_n)-mb);
end
end
figure;imshow(J_AdaptiveReduction);运行结果：
(a)被0.015高斯噪声污染 (b)算术均值滤波
(c)几何均值滤波 (d)自适应降噪滤波
4.7 自适应中值滤波(实例9)
=============================MATLAB代码============================
clc;
clear;
f=imread('hong.jpg');
image_gray=rgb2gray(f);%得到灰度图像
ff =image_gray;
f1=imnoise(image_gray,'salt & pepper',0.3);%添加椒盐噪声后的图像
alreadyProcessed = false(size(image_gray));%生成逻辑非的矩阵
% 迭代.
Smax=7;
for k = 3:2:Smax
zmin = ordfilt2(f1, 1, ones(k, k), 'symmetric');
zmax = ordfilt2(f1, k * k, ones(k, k), 'symmetric');
zmed = medfilt2(f1, [k k], 'symmetric');
processUsingLevelB = (zmed > zmin) & (zmax > zmed) & ...
~alreadyProcessed; %需要转到B步骤的像素
zB = (f1 > zmin) & (zmax > f1);
outputZxy = processUsingLevelB & zB;%满足步骤AB的输出原值对应像素的位置
outputZmed = processUsingLevelB & ~zB;%满足A不满足B的输出中值对应的像素位置
ff(outputZxy) = f1(outputZxy);
ff(outputZmed) = zmed(outputZmed);
alreadyProcessed = alreadyProcessed | processUsingLevelB;
if all(alreadyProcessed(:))
break;
end
end
ff(~alreadyProcessed) = zmed(~alreadyProcessed);
f2=medfilt2(f1,[7,7]);%中值滤波后的图像
imshow(image_gray);
figure;imshow(f1);
figure;imshow(f2);
figure;imshow(ff);运行结果：
(a)0.3的椒盐噪声污染 (b)7*7中值滤波 (c)Smax=7的自适应中值
1074
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
