案例7夜间图像去噪

崔丽（1）.噪声的扫描变换 现在图像系统的输入光电变换都是先把二维图像信号扫描变换成一维电信号再进行处理加工。最后再将一维电信号变成二维图像信号。噪声也存在着同样的变换方式。 （2）.噪声与图像的相关性 使用光导摄象管的摄像机，可以认为，信号幅度和噪声幅度无关。而使用超正析摄像机的信号和噪声相关，黑暗部分噪声大 ，明亮部分噪声小，在数字图像处理技术中量化噪声是肯定存在的，它和图像相位有关，如图像内容接近平坦时，量化噪声呈现伪轮廓，但在此时图像信号中的随机噪声就会因为颤噪效应反而使量化噪声变得不那么明显。 （3）.噪声的迭加性 在串联图像传输系统中，各部分窜入噪声若是同类噪声可以进行功率相加，依次信噪比要下降。若不是同类噪声应区别对待，而且要考虑视觉检出特性的影响。但是因为视觉检出特性中的许多问题还没有研究情趣，所以也只能进行一些主观的评价试验。如空间频率特性不同的噪声迭加要考虑到视觉空间频谱的带通特性。而时间特性不同的噪声迭加就要考虑视觉滞留和其闪烁的特性等等。亮度和色度噪声的迭加一定要清楚视觉的彩色特性。而以上的这些都因为视觉特性的未获解决而无法进行分析。从图像中找到一块平滑区域，根据统计数字来获得噪声方差，计算其均值和方差。Matlab命令：选取小块roipoly(f)灰度直方图imhist(f)自编程序_冈萨雷斯数字图像处理： 几何中心矩的计算statement 噪声分析imnoise2 定义:中心矩 中心坐标： 中心矩：( , )pqpqmx y f x y dxdy 01000/xmm00100/ymm向量的几何矩向量的几何矩00() ()( , )pqpqxxyyf x y dxdy0()( )ppmE xxf xStatement函数解读：function v, unv = statmoments(p, n)%STATMOMENTS 计算图像直方图的中心矩，输入：图像直方图p，n为矩的阶数返回值：v为中心矩，unv为没有平移中心的矩Imnoise2函数解读function R = imnoise2(type, M, N, a, b)输入：type表示什么样的噪声，M，N为产生噪声的矩阵大小，a，b为相应的噪声参数；输出：R为大小为M*N矩阵Type型噪声，该噪声带有参数为a，b。estimate_noise.m文件估计出噪声强度，通过矩来计算用imnoise2，把噪声重现查看属于哪一类噪声。通过imnoise,imnoise2函数实现。具体程序见fnoise.m模糊化处理程序EstimateH，如何找到这个模糊核？Matlab命令deconvreg J = deconvreg(I,PSF)J = deconvreg(I,PSF,NOISEPOWER)J = deconvreg(I,PSF,NOISEPOWER,LRANGE)程序解读：I是结果图像，PSF是滤波器（比原图像要小），NOISEPOWER是噪声强度（实数）去噪函数见Constrained Least squares Filtering（CLSF.m）模型是 g = hf+还有待研究。在共享的文件夹里面，还有若干函数，可以供大家学习各种滤波。如spatialfilter.m(空间滤波)课上可完成学习inversefilter.m(反滤波和Wiener滤波)Inversewiener.m (Wiener 滤波)Notchfilter.m(Notch 滤波)OptimalNotch.m(最优Notch滤波)推荐有兴趣同学看书冈萨雷斯数字图像处理（Matlab版）Matlab命令CA,CH,CV,CD = DWT2(X,wname)或者CA,CH,CV,CD = DWT2(X,Lo_D,Hi_D)；X = IDWT2(CA,CH,CV,CD,wname)或者X = IDWT2(CA,CH,CV,CD,Lo_R,Hi_R)；单独重构可用代替不用的部分。Wavemenuwavelet 2D对信号进行小波分解针对小波系数进行处理方法： 模极大值去噪 阈值法去噪 相关性去噪重构信号衡量标准缺点：硬阈值函数在w=T时是不连续的，用它重构信号时会产生振荡；软阈值函数虽然连续性好，但|w|T时，它与w存在恒定的偏差，直接影响重构信号的性质。 为噪声方差，N为信号长度，T值估计应考虑平稳性与信噪比（自适应的选取）ddencmp，thselect，wbmpen，wdcbm2实现小波去噪的函数有wden，wdencmp，wpdencmp，wthresh，wpthcoef和wthcoef2Wavemenu工具箱一维二维小波变换对应下的denoise。nblocr1.mat；nblocr2.mat（带噪的方波信号）nbumpr1.mat； nbumpr2.mat； nbumpr3.mat；ndoppr1.mat；nelec.matnoischir.mat;noisbloc.mat；noisbump.mat；noisdopp.mat；noismima.mat；noissin.mat（高斯噪声）noiswom.mat；noissi2d.mat; nbarb1.mat1.读入信号或者图像X2.加入随机噪声， 随机噪声 e*rand（M，N）;3.带有噪声的图像为 X随机噪声Matlab命令：不同阈值选择方法下的阈值 THR = THSELECT(X,TPTR) 用TPTR字符串定义的阈值选取方法选择与输入信号X有关的阈值。 输入参数：X为带噪信号 TPTR为阈值选择准则的字符串： rigrsure，用 Steins 无偏估计的原理进行的自适应阈值； heursure是 启发式方法； sqtwolog 通用阈值方法，阈值为； minimaxi 最大最小阈值法。 阈值选取方法的噪声模型为y = f(t) + e 其中e 是满足N(0,1)的白噪声 处理非白噪声或者非校正的噪声可以用SCAL参数来重新调整，见wden命令。 返回值：THR为阈值。进行软硬阈值函数运算命令 Y = WTHRESH(X,SORH,T) 输入参数：X为输入矩阵，SORH为s(软阈值函数)，h(硬阈值函数) ，T是阈值T 返回值：用软(硬)阈值处理后的矩阵。XD,CXD,LXD = WDEN(X,TPTR,SORH,SCAL,N,wname) 1-D信号去噪函数 输入参数： X为原始带噪信号， TPTR为阈值选择准则的字符串：见 THSELECT命令 SORH : 选择s 或者h ，h是硬阈值，s为软阈值。见WTHRESH命令。 SCAL定义多种方法的阈值重新调节: one ：没有重新调节。 sln 基于第一层系数进行噪声估计来重新调节 mln层层独立进行噪声估计来重新调节 N 为小波分解层数， wname为小波函数 返回值： XD为去噪后的信号；CXD,LXD去噪后信号的小波分解的系数。另一种表达： XD,CXD,LXD = WDEN(C,L,TPTR,SORH,SCAL,N,wname) 与第一种不同之处在于用分解系数第N层的小波分解系数C,L代替直接的带噪信号X% 在原始图像中加入Gaussian白噪声，椒盐噪声和乘性噪声。 f=imread(original_pattern.tif); m=64/255; var = 400/2552; g_gauss=imnoise(f,gaussian,m,var); d = 0.05;%d表示噪声强度 g_salt = imnoise(f,salt & pepper,d); v = 0.06; g_speckle = imnoise(f,speckle, v);% 显示加噪图像 figure (1)subplot(2,2,1), imshow(f),title(原始图像);subplot(2,2,2), imshow(g_gauss), title(Guassian 噪声图像);subplot(2,2,3), imshow(g_salt), title(椒盐噪声图像);subplot(2,2,4), imshow(g_speckle), title( speckle 噪声图像); figure(2)title(直方图对比);subplot(2,2,1), imhist(f);%原始图像的灰度直方图subplot(2,2,2), imhist(g_gauss);%Guassian 噪声图像的直方图subplot(2,2,3), imhist(g_salt);%椒盐噪声图像的直方图subplot(2,2,4), imhist(g_speckle);%speckle噪声图像的直方图function v, unv = statmoments(p, n)% STATMOMENTS计算图像灰度直方图的统计中心矩。% 参数说明：输入p为图像直方图，n为矩的阶数；% 返回值V(1)为均值, V(2)为方差,V(3)为三阶中心矩，以此类推V(n)为第n阶中心矩。% 返回值UNV为像素值取值再0,255非归一化的时候的矩的值。Lp = length(p);G = Lp - 1;% 直方图归一化p = p/sum(p);% 灰度值随机变量及其归一化取值范围到0, 1.z = 0:G;z = z./G;% 均值m = z*p;% 平移中心z = z - m;% 计算中心矩v = zeros(1, n);v(1) = m;for j = 2:n v(j) = (z.j)*p;endif nargout 1 % 计算非中心矩 unv = zeros(1, n); unv(1)=m.*G; for j = 2:n unv(j) = (z*G).j)*p; endend* g_Gaussian=imread(gaussian-noise.tif);% 取出图像的一块 B,c,r=roipoly(g_Gaussian);% 画出原始图像f在B这块的直方图p p=imhist(g_Gaussian(B);% 求和 npix = sum(B(:);%计算2阶中心矩 v, unv = statmoments(p, 2) vv =0.7868 0.0033 unvunv = 200.6315 214.4707%生成Gaussian噪声，均值为unv(1)，方差为unv(2) X =unv(1)+sqrt(unv(2)*randn(npix,1); M = max(p)M =86 subplot(2,2,1),imshow(g_Gaussian),title(噪声图像);subplot(2,2,2),imshow(B),title(交互模式选取的区域);subplot(2,2,3),bar(p,1);subplot(2,2,4),hist(noise,130),axis(0,300,0,150)*%计算两个信号f和g的峰值信噪比tfunction t = PSNR(f,g)m,n = size(f);m1,n1 = size(g);if(m1=m |n1=n) error(Signals are not the same size.,msg)endt = 10*log(2552/(1/(m*n)*sum(f-g).2 );*%计算两个信号f和g的均方根误差tfunction t = mse(f,g)m,n = size(f);m1,n1 = size(g);if(m1=m |n1=n) error(Signals are not the same size.,msg)endt = 1/(m*n)*sum(f-g).2 );* load noisbloc;X=noisbloc; Y=fft(X); plot(abs(Y); T=abs(Y);tt=find(T1000);Y(tt)=0;IY=ifft(Y);subplot(2,1,1), plot(X),title(原始信号);subplot(2,1,2),plot(IY),title(Fourier变换去噪后信号);子函数：用模极大值重构信号，采用交错投影法，见李世雄小波变换域应用function int=Py(int,len);%Py投影function int=Pg(int,lev,sr);%Pg投影function w2=Pyg(w1,wp, sr); Pg和Py投影去噪函数function s=mden(f,lev,n,wf)%单区间Py投影：对区间进行裁减即Py投影，返回裁剪后的区间信号%输入参数：int为单区间的点、len为区间中点的个数%返回值int为Py投影后的区间function int=Py(int,len);if sign(int(1)=sign(int(len)%若区间左右端点同号 int=int.*(sign(int)=sign(int(1);%区间上，只保留本身符号与左端点符号相同的那些inte=interp1(1,len,int(1),int(len),(1:len),linear);%只用区间端点处的值，在区间上进行线性插值，得到区间的差值后的点，与原来区间点个数相同int=sign(int(1)*(abs(inte)-(abs(inte)-abs(int).*(abs(inte)-abs(int)0);%只保留区间上的点的模值比线性插值对应点模值小的那些点，其余的为0，符号与左端点一致。else sgn=sign(int(len)-int(1);%若区间端点不同号，则取相减的符号。两极值点异号，中间有单调性 intem=max(int(1),int(len);将端点中最大最小点分别赋值inten=min(int(1),int(len);%从区间端点开始，循环找寻该区间的极大值点和极小值点 for i=1:len-2 if sign(int(i+1)-int(i)=sgn%若差商符号不是sgn int(i+1)=int(i);%令二者相等，保证两个极值点之间的单调性 end if int(i+1)intem%比端点值大 int(i+1)=intem;%则是端点值 end if int(i+1)2 w1(m,fr(j):fr(j+1)=Py(int,len);%区间内多于两个点，即中间有点，则进行Py投影 end end % 再逐一区间进行Pgama投影 for j=1:num_int int=err(m,fr(j):fr(j+1);%在两极值点之间的原来与现在的差 err(m,fr(j):fr(j+1)=Pg(int,m,sr); end w2(m,:)=w1(m,:)+err(m,:);%投影最终结果，对应参考文献【3】中的（5.28）式。end*%模极大值法去噪：%参数C，ep，sr可以调节。%输入参数为：f为带去噪信号，lev为小波分解层数，n去噪重构算法的迭代次数，wf为小波函数名称function s=mden(f,lev,n,wf)pp=size(f);pp=pp(2);%所处理数据的长度 sr=360;%抽样率Lo_D,Hi_D,Lo_R,Hi_R=wfilters(wf);%二进小波变换swa,swd = swt(f,lev,Lo_D,Hi_D);%求二进小波变换的模极大值及其位置%初始结果%dw为局部极大位置，wp为局部极大序列。dw=zeros(size(swd);zdw=dw;fdw=dw;%要找模极大值，把小波系数中大于0，和小于0的分别考虑。%小波系数大于0的赋值给zwzw=swd.*(swd0);%留下左边元素比右边元素小的值的位置记为1，存入zdw中zdw=(zw(:,1:pp-1)-zw(:,2:pp)0);%小波系数小于0的赋值给fwfw=swd.*(swd0);%再计算1,0 间隔的点，即找到模极大值的点的位置fdw(:,2:pp-1)=(fdw(:,1:pp-2)-fdw(:,2:pp-1)0);%将zdw和fdw中的模极大值点位置合并dw=zdw|fdw;%第一列最后一列保留dw(:,1)=1;dw(:,pp)=1;%wp存放模极大值点的值wp=dw.*swd;%处理模极大值：最高层的模极大值点开始Dwp(:,lev)=wp(lev,:);M=max(Dwp(:,lev);%模极大值的阈值用去噪函数ddencmp来获得。Thr,sorh,keepapp = ddencmp(den,wv,f);C=0.8;%阈值参数Thr=C*M/lev; %阈值计算，见参考文献”夏敏华，3mm波段脉冲雷达系统研究和小波去噪分析”。%将大于阈值Thr的最后一层的模极大值留下Dwp(:,lev)=Dwp(:,lev).*(abs(Dwp(:,lev)Thr);%处理最后一层%模极大值的处理方式：%在尺度j上极大值点位置，构造一个搜索区域，%在尺度j-1中，极大值点落在该区域的点保留，其他的置0；ep=3;%该参数确定邻域for j=1:lev-2 Dwp(:,lev-j)=wp(lev-j,:); Dp(:,lev-j+1)=(Dwp(:,lev-j+1)=0); DD=Dp(:,lev-j+1); for Pd=ep:(length(Dp(:,lev-j+1)-ep); if DD(Pd)=1 for i=1:ep-1; DD(Pd-i)=1; end end end Dp(:,lev-j+1)=DD; Dwp(:,lev-j)=Dwp(:,lev-j).*Dp(:,lev-j+1); Dp(:,lev-j)=(Dwp(:,lev-j)=0);end%第一层单独处理，在第二层极大值点位置上，保留第一层相应极大值点：Dwp(:,1)=wp(1,:);Dwp(:,1)=Dwp(:,1).*Dp(:,2);wp=Dwp;wp=wp;%重构信号：s=swa(lev,:); % s: 为待重建的信号w0=zeros(1,pp);a,d=swt(w0,lev,Lo_D,Hi_D);%从0开始w2=d;%重建高频for j=1:n w2=Pyg(d,wp,sr);%Py投影和Pgama投影 w0=iswt(s,w2,Lo_R,Hi_R);%Pv投影 a,d=swt(w0,lev,Lo_D,Hi_D);ends=iswt(swa(lev,:),w2,Lo_R,Hi_R); % 计算重建信号figure,subplot(2,1,1)plot(f,r);grid on;title(原始信号);subplot(2,1,2);plot(s);grid on;title(模极值去噪后信号)* load sumsin;%读入带噪信号 f = sumsin zeros(1,24);%补零使得信号变成2的幂次，这里是信号长度变为1024。 lev = 4;%模极大值提取时候的分解层数 n = 20;%模极大值重构信号时候的迭代次数 wf = db3;%采用的小波函数 s =mden(f,lev,n,wf);%模极大值去噪算法阈值法去噪法load sumsin;f = sumsin zeros(1,24);lev = 4;n = 5;wf = db3;%软硬阈值去噪s_s=wden(f,heursure,s,one,lev,wf);s_h=wden(f,heursure,h,one,lev,wf);%计算PSNR值P1 = PSNR(f,s_s)P2 = PSNR(f,s_h)m1 = mse(f,s_s)m2 = mse(f,s_h)figure,subplot(3,1,1)plot(f,r);grid on;title(原始信号);subplot(3,1,2);plot(s_s);grid on;title(软阈值去噪后信号)subplot(3,1,3);plot(s_h);grid on;title(硬阈值去噪后信号)P1 = 116.1074P2 = 116.1912m1 = 0.5897m2 =%输入参数：f为带噪信号，n为分解层数，wf为小波函数名称%返回值f_n为去噪信号function f_n=xden(f,n,wf)Lo_D,Hi_D,Lo_R,Hi_R=wfilters(wf);swa,swd = swt(f,n,Lo_D,Hi_D);%swd是细节系数，swa是近似系数ss=swd; s1=swd;T=zeros(size(ss);%先把系数处理矩阵设置为全0。%在1:n-1分解层次上对高频系数处理，最后一层无法求相关系数，所以不作处理。for j=1:n-1 %nd=s1(j,:); %nd=nd(1:80); %nv=var(nd); %以信号的前80个只含有噪声的点估计噪声在各层的方差。 %pwv=var(s1(j,:); cw=s1(j,:).*s1(j+1,:); %定义相关系数为相邻两层的乘积。 %thr,sorh,keepapp = ddencmp(den,wv,f); %_用以设定停止迭代的 噪声能量阈值，需要根据情况调节。_% %cc=thr; %while pwvcc*nv pw=sum(abs(swd(j,:).2); %计算小波能量 pcw=sum(abs(cw).2); %计算相关系数能量 w_n=cw.*(pw/pcw)0.5); %归一化 acw=abs(w_n); aw=abs(swd(j,:); ss=swd(j,:).*(acwaw);%(acwaw)返回0或者1,1为信号，0为噪声 ss1=(ss=0); T(j,:)=T(j,:)+ss1; %将选出的点赋给系数处理矩阵相应位置。 %ss0=ones(size(ss1); %ss0=ss0-ss1; %swd(j,:)=swd(j,:).*ss0; %将高频系数选出大值后的地方置0。 %pwv=var(swd(j,:); %w_n=w_n.*ss0; %将相关系数选出大值后的地方置0。 %cw=w_n; %endend mask_max=ones(1,length(T);T=T(1:(n-1),:);mask_max; %最后一层系数处理矩阵全置1。s1=s1.*T;f_n=iswt(swa,s1,wf);%画图：figure; %信号滤波前后比较。subplot(2,1,1);plot(f,r); grid on;title(原始信号);subplot(2,1,2);plot(f_n); grid on;title(相关性去噪后信号)*注意：相关性去噪也可以采用先估计噪声方差，进而阈值处理一下，再进行相关性去噪另外，相关性去噪的改进算法，还可以继续研究。针对实际图片的去噪实验夜间模式下的照片去噪和复原，借助同一时间段的傻瓜照片