高斯模糊图像的盲复原附带matlab源代码
- 格式:pdf
- 大小:1.50 MB
- 文档页数:39
matlab模糊函数代码在数学和图像处理领域中,模糊函数是一种常用的工具,用于对图像进行模糊处理以达到一定的效果。
Matlab提供了一些内置函数来实现图像的模糊处理,本文将介绍如何使用Matlab编写模糊函数代码。
Matlab中有多种不同类型的模糊函数,例如高斯模糊、均值模糊和运动模糊等。
下面将逐一介绍这些模糊函数代码的实现方式。
1. 高斯模糊:高斯模糊是最常用的模糊算法之一,它通过卷积图像与一个高斯核来实现。
以下是Matlab中实现高斯模糊的代码示例:```matlabfunction blurredImage = gaussianBlur(image, sigma)kernelSize = 2 * ceil(3 * sigma) + 1; % 根据sigma计算高斯核大小kernel = fspecial('gaussian', [kernelSize kernelSize], sigma); % 生成高斯核blurredImage = imfilter(image, kernel, 'conv'); % 对图像进行卷积操作end```2. 均值模糊:均值模糊是一种简单但常用的模糊算法,它通过计算邻域像素的平均值来实现。
以下是Matlab中实现均值模糊的代码示例:```matlabfunction blurredImage = meanBlur(image, kernelSize)kernel = ones(kernelSize) / (kernelSize^2); % 生成均值核blurredImage = imfilter(image, kernel, 'conv'); % 对图像进行卷积操作end```3. 运动模糊:运动模糊是一种模糊算法,它通过模拟相机快门打开时的移动效果来实现。
以下是Matlab中实现运动模糊的代码示例:```matlabfunction blurredImage = motionBlur(image, angle, distance)PSF = fspecial('motion', distance, angle); % 生成运动模糊核blurredImage = imfilter(image, PSF, 'conv'); % 对图像进行卷积操作end```以上是几种常见的模糊函数的Matlab代码实现。
% *-*--*-*-*-*-*-*-*-*-*-*-*图像处理*-*-*-*-*-*-*-*-*-*-*-*%{% (一)图像文件的读/写A=imread('drum.jpg'); % 读入图像imshow(A); % 显示图像imwrite(A,'drum.jpg');info=imfinfo('drum.jpg') % 查询图像文件信息% 用colorbar函数将颜色条添加到坐标轴对象中RGB=imread('drum.jpg');I=rgb2gray(RGB); % 把RGB图像转换成灰度图像h=[1 2 1;0 0 0;-1 -2 -1];I2=filter2(h,I);imshow(I2,[]);colorbar('vert') % 将颜色条添加到坐标轴对象中% wrap函数将图像作为纹理进行映射A=imread('4.jpg');imshow(A);I=rgb2gray(RGB);[x,y,z]=sphere;warp(x,y,z,I); % 用warp函数将图像作为纹理进行映射%}% subimage函数实现一个图形窗口中显示多幅图像RGB=imread('drum.jpg');I=rgb2gray(RGB);subplot(1,2,1);subimage(RGB); % subimage函数实现一个图形窗口中显示多幅图像subplot(1,2,2),subimage(I);% *-*--*-*-*-*-*-*-*-*-*-*-*图像处理*-*-*-*-*-*-*-*-*-*-*-*% (二)图像处理的基本操作% ----------------图像代数运算------------------%{% imadd函数实现两幅图像的相加或给一幅图像加上一个常数% 给图像每个像素都增加亮度I=imread('4.jpg');J=imadd(I,100); % 给图像增加亮度subplot(1,2,1),imshow(I);title('原图');subplot(1,2,2),imshow(J);title('增加亮度图');%% imsubtract函数实现将一幅图像从另一个图像中减去或减去一个常数I=imread('drum.jpg');J=imsubtract(I,100); % 给图像减去亮度subplot(1,2,1),imshow(I);%% immultiply实现两幅图像的相乘或者一幅图像的亮度缩放I=imread('drum.jpg');J=immultiply(I,2); % 进行亮度缩放subplot(1,2,1),imshow(I);subplot(1,2,2),imshow(J);%% imdivide函数实现两幅图像的除法或一幅图像的亮度缩放I=imread('4.jpg');J=imdivide(I,0.5); % 图像的亮度缩放subplot(1,2,1),imshow(I);subplot(1,2,2),imshow(J);%}% ----------------图像的空间域操作------------------%{% imresize函数实现图像的缩放J=imread('4.jpg');subplot(1,2,1),imshow(J);title('原图');X1=imresize(J,0.2); % 对图像进行缩放subplot(1,2,2),imshow(X1);title('缩放图');%% imrotate函数实现图像的旋转I=imread('drum.jpg');J=imrotate(I,50,'bilinear'); % 对图像进行旋转subplot(1,2,1),imshow(I);subplot(1,2,2),imshow(J);%% imcrop函数实现图像的剪切I=imread('drum.jpg');I2=imcrop(I,[1 100 130 112]); % 对图像进行剪切subplot(1,2,1),imshow(I);subplot(1,2,2),imshow(I2);%}% ----------------特定区域处理------------------%{% roipoly函数用于选择图像中的多边形区域I=imread('4.jpg');c=[200 250 278 248 199 172];r=[21 21 75 121 121 75];BW=roipoly(I,c,r); % roipoly函数选择图像中的多边形区域subplot(1,2,1),imshow(I);subplot(1,2,2),imshow(BW);%% roicolor函数式对RGB图像和灰度图像实现按灰度或亮度值选择区域进行处理a=imread('4.jpg');subplot(2,2,1),imshow(a);I=rgb2gray(a);BW=roicolor(I,128,225); % 按灰度值选择的区域subplot(2,2,4),imshow(BW);%% ploy2mask 函数转化指定的多边形区域为二值掩模x=[63 186 54 190 63];y=[60 60 209 204 601];bw=poly2mask(x,y,256,256); % 转化指定的多边形区域为二值掩模imshow(bw);hold onplot(x,y,'r','LineWidth',2);hold off%% roifilt2函数实现区域滤波a=imread('4.jpg');I=rgb2gray(a);c=[200 250 278 248 199 172];r=[21 21 75 121 121 75];BW=roipoly(I,c,r); % roipoly函数选择图像中的多边形区域h=fspecial('unsharp');J=roifilt2(h,I,BW); % 区域滤波subplot(1,2,1),imshow(I);subplot(1,2,2),imshow(J);%% roifill函数实现对特定区域进行填充a=imread('4.jpg');I=rgb2gray(a);c=[200 250 278 248 199 172];r=[21 21 75 121 121 75];J=roifill(I,c,r); % 对特定区域进行填充subplot(1,2,1),imshow(I);subplot(1,2,2),imshow(J);%}% ----------------图像变换------------------%{% fft2 和ifft2函数分别是计算二维的快速傅里叶变换和反变换f=zeros(100,100);subplot(1,2,1);imshow(f);f(20:70,40:60)=1;subplot(1,2,2);imshow(f);F=fft2(f); % 计算二维的快速傅里叶变换F2=log(abs(F));% 对幅值对对数figure;subplot(1,2,1),imshow(F),colorbar;subplot(1,2,2),imshow(F2),colorbar;%% fftsshift 函数实现了补零操作和改变图像显示象限f=zeros(100,100);subplot(2,2,1),imshow(f);title('f')f(10:70,40:60)=1;subplot(2,2,2),imshow(f);title('f取后')F=fft2(f,256,256);subplot(2,2,3),imshow(F);title('F')F2=fftshift(F); % 实现补零操作subplot(2,2,4),imshow(F2);title('F2')figure,imshow(log(abs(F2)));title('log(|F2|)')%% dct2 函数采用基于快速傅里叶变换的算法,用于实现较大输入矩阵的离散余弦变换% idct2 函数实现图像的二维逆离散余弦变换RGB=imread('drum.jpg');I=rgb2gray(RGB);J=dct2(I); % 对I进行离散余弦变换imshow(log(abs(J))),title('对原图离散后取对数'),colorbar;J(abs(J)<10)=0;K=idct2(J); % 图像的二维逆离散余弦变换figure,imshow(I),title('原灰度图')figure,imshow(K,[0,255]);title('逆离散变换');%% dctmtx 函数用于实现较小输入矩阵的离散余弦变figure;RGB=imread('4.jpg');I=rgb2gray(RGB);subplot(3,2,1),imshow(I),title('原灰度图');I=im2double(I);subplot(3,2,2),imshow(I),title('取双精度后');T=dctmtx(8); % 离散余弦变换subplot(3,2,3),imshow(I),title('离散余弦变换后');B=blkproc(I,[8,8],'P1*x*P2',T,T');subplot(3,2,4),imshow(B),title('blkproc作用I后的B');mask=[ 1 1 1 1 0 0 0 01 1 1 0 0 0 0 01 1 0 0 0 0 0 01 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 ];B2=blkproc(B,[8,8],'P1.*x',mask);subplot(3,2,5),imshow(B2),title('blkproc作用B后的B2');I2=blkproc(B2,[8,8],'P1*x*P2',T',T);subplot(3,2,6),imshow(I2),title('blkproc作用B2后的I2');%% edge函数用于提取图像的边缘RGB=imread('4.jpg');I=rgb2gray(RGB);BW=edge(I);imshow(I);figure,imshow(BW);%% radon 函数用来计算指定方向上图像矩阵的投影RGB=imread('4.jpg');I=rgb2gray(RGB);BW=edge(I);theta=0:179;[R,XP]=radon(BW,theta); % 图像矩阵的投影figure,imagesc(theta,XP,R);colormap(hot);xlabel('\theta(degrees)');ylabel('x\prime');title('R_{\theta}(x\prime)');colorbar;%}% ----------------图像增强、分割和编码------------------%{% imhist 函数产生图像的直方图A=imread('4.jpg');B=rgb2gray(A);subplot(2,1,1),imshow(B);subplot(2,1,2),imhist(B);%% histeq 函数用于对图像的直方图均衡化A=imread('4.jpg');B=rgb2gray(A);subplot(2,1,1),imshow(B);subplot(2,1,2),imhist(B);C=histeq(B); % 对图像B进行均衡化figure;subplot(2,1,1),imshow(C);subplot(2,1,2),imhist(C);%% filter2 函数实现均值滤波a=imread('4.jpg');I=rgb2gray(a);subplot(2,2,1),imshow(I);K1=filter2(fspecial('average',3),I)/255; % 3*3的均值滤波K2=filter2(fspecial('average',5),I)/255; % 5*5的均值滤波K3=filter2(fspecial('average',7),I)/255; % 7*7的均值滤波subplot(2,2,2),imshow(K1);subplot(2,2,3),imshow(K2);subplot(2,2,4),imshow(K3);%% wiener2 函数实现Wiener(维纳)滤波a=imread('4.jpg');I=rgb2gray(a);subplot(2,2,1),imshow(I);K1=wiener2(I,[3,3]); % 3*3 wiener滤波K2=wiener2(I,[5,5]); % 5*5 wiener滤波K3=wiener2(I,[7,7]); % 7*7 wiener滤波subplot(2,2,2),imshow(K1);subplot(2,2,3),imshow(K2);subplot(2,2,4),imshow(K3);%% medfilt2 函数实现中值滤波a=imread('4.jpg');I=rgb2gray(a);subplot(2,2,1),imshow(I);K1=medfilt2(I,[3,3]); % 3*3 中值滤波K2=medfilt2(I,[5,5]); % 5*5 中值滤波K3=medfilt2(I,[7,7]); % 7*7 中值滤波subplot(2,2,2),imshow(K1);subplot(2,2,3),imshow(K2);subplot(2,2,4),imshow(K3);%}% ----------------图像模糊及复原------------------%{% deconvwnr 函数:使用维纳滤波器I=imread('qier.jpg');imshow(I);% 对图像进行模糊处理LEN=31;THETA=11;PSF1=fspecial('motion',LEN,THETA); % 运动模糊PSF2=fspecial('gaussian',10,5); % 高斯模糊Blurred1=imfilter(I,PSF1,'circular','conv'); % 得到运动模糊图像Blurred2=imfilter(I,PSF2,'conv'); % 得到高斯噪声模糊图像figure;subplot(1,2,1);imshow(Blurred1);title('Blurred1--"motion"'); subplot(1,2,2);imshow(Blurred2);title('Blurred2--"gaussian"');% 对模糊图像加噪声V=0.002;BlurredNoisy1=imnoise(Blurred1,'gaussian',0,V); % 加高斯噪声BlurredNoisy2=imnoise(Blurred2,'gaussian',0,V); % 加高斯噪声figure;subplot(1,2,1);imshow(BlurredNoisy1);title('BlurredNoisy1'); subplot(1,2,2);imshow(BlurredNoisy2);title('BlurredNoisy2');% 进行维纳滤波wnr1=deconvwnr(Blurred1,PSF1); % 维纳滤波wnr2=deconvwnr(Blurred2,PSF2); % 维纳滤波figure;subplot(1,2,1);imshow(wnr1);title('Restored1,True PSF'); subplot(1,2,2);imshow(wnr2);title('Restored2,True PSF');%% deconvreg函数:使用约束最小二乘滤波器I=imread('qier.jpg');imshow(I);% 对图像进行模糊处理LEN=31;THETA=11;PSF1=fspecial('motion',LEN,THETA); % 运动模糊PSF2=fspecial('gaussian',10,5); % 高斯模糊Blurred1=imfilter(I,PSF1,'circular','conv'); % 得到运动模糊图像Blurred2=imfilter(I,PSF2,'conv'); % 得到高斯噪声模糊图像figure;subplot(1,2,1);imshow(Blurred1);title('Blurred1--"motion"');subplot(1,2,2);imshow(Blurred2);title('Blurred2--"gaussian"');% 对模糊图像加噪声V=0.002;BlurredNoisy1=imnoise(Blurred1,'gaussian',0,V); % 加高斯噪声BlurredNoisy2=imnoise(Blurred2,'gaussian',0,V); % 加高斯噪声figure;subplot(1,2,1);imshow(BlurredNoisy1);title('BlurredNoisy1');subplot(1,2,2);imshow(BlurredNoisy2);title('BlurredNoisy2');NP=V*prod(size(I));reg1=deconvreg(BlurredNoisy1,PSF1,NP); % 约束最小二乘滤波reg2=deconvreg(BlurredNoisy2,PSF2,NP); % 约束最小二乘滤波figure;subplot(1,2,1);imshow(reg1);title('Restored1 with NP');subplot(1,2,2);imshow(reg2);title('Restored2 with NP');%% deconvlucy函数:使用Lucy-Richardson滤波器I=imread('qier.jpg');imshow(I);% 对图像进行模糊处理LEN=31;THETA=11;PSF1=fspecial('motion',LEN,THETA); % 运动模糊PSF2=fspecial('gaussian',10,5); % 高斯模糊Blurred1=imfilter(I,PSF1,'circular','conv'); % 得到运动模糊图像Blurred2=imfilter(I,PSF2,'conv'); % 得到高斯噪声模糊图像figure;subplot(1,2,1);imshow(Blurred1);title('Blurred1--"motion"');subplot(1,2,2);imshow(Blurred2);title('Blurred2--"gaussian"');% 对模糊图像加噪声V=0.002;BlurredNoisy1=imnoise(Blurred1,'gaussian',0,V); % 加高斯噪声BlurredNoisy2=imnoise(Blurred2,'gaussian',0,V); % 加高斯噪声figure;subplot(1,2,1);imshow(BlurredNoisy1);title('BlurredNoisy1');subplot(1,2,2);imshow(BlurredNoisy2);title('BlurredNoisy2');luc1=deconvlucy(BlurredNoisy1,PSF1,5); % 使用Lucy-Richardson滤波luc2=deconvlucy(BlurredNoisy1,PSF1,15); % 使用Lucy-Richardson滤波figure;subplot(1,2,1);imshow(luc1);title('Restored Image,NUMIT=5'); subplot(1,2,2);imshow(luc2);title('Restored Image,NUMIT=15');%}% deconvblind 函数:使用盲卷积算法a=imread('4.jpg');I=rgb2gray(a);figure;imshow(I);title('Original Image');PSF=fspecial('motion',13,45); % 运动模糊figure;imshow(PSF);Blurred=imfilter(I,PSF,'circ','conv'); % 得到运动模糊图像figure;imshow(Blurred);title('Blurred Image');INITPSF=ones(size(PSF));[J,P]=deconvblind(Blurred,INITPSF,30); % 使用盲卷积figure;imshow(J);figure;imshow(P,[],'notruesize');% *-*--*-*-*-*-*-*-*-*-*-*-*图像处理*-*-*-*-*-*-*-*-*-*-*-* %{% 对图像进行减采样a=imread('lena.jpg');%subplot(1,4,1);figure;imshow(a);title('原图');b=rgb2gray(a);%subplot(1,4,2);figure;imshow(b);title('原图的灰度图');[wid,hei]=size(b);%---4倍减采样----quartimg=zeros(wid/2+1,hei/2+1);i1=1;j1=1;for i=1:2:widfor j=1:2:heiquartimg(i1,j1)=b(i,j);j1=j1+1;endi1=i1+1;j1=1;end%subplot(1,4,3);figure;imshow(uint8(quartimg));title('4倍减采样')% ---16倍减采样---quanrtimg=zeros(wid/4+1,hei/4+1);i1=1;j1=1;for i=1:4:widfor j=1:4:heiquanrtimg(i1,j1)=b(i,j);j1=j1+1;endi1=i1+1;j1=1;end%subplot(1,4,4);.figure;imshow(uint8(quanrtimg));title('16倍减采样');%}% 图像类型% 将图像转换为256级灰度图像,64级灰度图像,32级灰度图像,8级灰度图像,2级灰度图像a=imread('4.jpg');%figure;subplot(2,3,1);imshow(a);title('原图');b=rgb2gray(a); % 这是256灰度级的图像%figure;subplot(2,3,2);imshow(b);title('原图的灰度图像');[wid,hei]=size(b);img64=zeros(wid,hei);img32=zeros(wid,hei);img8=zeros(wid,hei);img2=zeros(wid,hei);for i=1:widfor j=j:heiimg64(i,j)=floor(b(i,j)/4); % 转化为64灰度级endend%figure;subplot(2,3,3);imshow(uint8(img64),[0,63]);title('64级灰度图像');for i=1:widfor j=1:heiimg32(i,j)=floor(b(i,j)/8);% 转化为32灰度级endend%figure;subplot(2,3,4);imshow(uint8(img32),[0,31]);title('32级灰度图像');for i=1:widfor j=1:heiimg8(i,j)=floor(b(i,j)/32);% 转化为8灰度级endend%figure;subplot(2,3,5);imshow(uint8(img8),[0,7]);title('8级灰度图像');for i=1:widfor j=1:heiimg2(i,j)=floor(b(i,j)/128);% 转化为2灰度级endend%figure;subplot(2,3,6);imshow(uint8(img2),[0,1]);title('2级灰度图像');% *-*--*-*-*-*-*-*-*-*-*-*-*图像处理*-*-*-*-*-*-*-*-*-*-*-* %{% ------------------ 图像的点运算------------------I=imread('lena.jpg');figure;subplot(1,3,1);imshow(I);title('原图的灰度图');J=imadjust(I,[0.3;0.6],[0.1;0.9]); % 设置灰度变换的范围subplot(1,3,2);imshow(J);title('线性扩展');I1=double(I); % 将图像转换为double类型I2=I1/255; % 归一化此图像C=2; % 非线性扩展函数的参数K=C*log(1+I2); % 对图像的对数变换subplot(1,3,3);imshow(K);title('非线性扩展');M=255-I;figure;subplot(1,3,1);imshow(M);title('灰度倒置');N1=im2bw(I,0.4); % 将此图像二值化,阈值为0.4N2=im2bw(I,0.7); % 将此图像二值化,阈值为0.7 subplot(1,3,2);imshow(N1);title('二值化阈值0.4');subplot(1,3,3);imshow(N2);title('二值化阈值0.7');%}%{% ------------------ 图像的代数运算------------------% 将两幅图像进行加法运算I=imread('lena.jpg');I=rgb2gray(I);J=imread('rice.png');% 以下把两幅图转化为大小一样for i=1:size(I)for j=size(J):size(I)J(i,j)=0;endendI=im2double(I); % 将图像转化为double型J=im2double(J);% imshow(I);figure;imshow(J);K=I+0.3*J; % 将两幅图像相加subplot(1,3,1);imshow(I);title('人物图');subplot(1,3,2);imshow(J);title('背景图');subplot(1,3,3);imshow(K);title('相加后的图');imwrite(K,'i_lena1.jpg');%%% 将两幅图像做减运算,分离背景与原图A=imread('i_lena1.jpg');B=imread('rice.png');% 以下把两幅图转化为大小一样for i=1:size(A)for j=size(B):size(A)B(i,j)=0;endendC=A-0.3*B;a=imread('lena.jpg');subplot(2,2,1);imshow(a);title('原图图');subplot(2,2,2);imshow(A);title('混合图');subplot(2,2,3);imshow(B);title('背景图');subplot(2,2,4);imshow(C);title('分离后的图');%% 设置掩模,需要保留下来的区域,掩模图像的值为1,否则为0 A=imread('drum.jpg');A=rgb2gray(A);A=im2double(A);sizeA=size(A);subplot(1,2,1);imshow(A);title('原图');B=zeros(sizeA(1),sizeA(2)); % 设置模板B(100:400,100:500)=1;K=A.*B; % 两幅图像相乘subplot(1,2,2);imshow(K);title('局部图');%}%{% ------------------ 图像的缩放------------------A=imread('drum.jpg');B1=imresize(A,1.5); % 比例放大1.5杯,默认采用的是最近邻法进行线性插值B2=imresize(A,[420 384]); % 非比例放大到420:384C1=imresize(A,0.7); % 比例缩小0.7倍C2=imresize(A,[150 180]); % 非比例缩小到150:180figure;imshow(B1);title('比例放大图');figure;imshow(B2);title('非比例放大图');figure;imshow(C1);title('比例缩小图');figure;imshow(C2);title('非比例缩小图');% 检测非比例缩放得到的图片是否能还原到原图a=size(A)d=imresize(C2,[a(1),a(2)]);figure;imshow(d);%}% ------------------ 图像的旋转------------------I=imread('drum.jpg');J=imrotate(I,45); % 图像进行逆时针旋转,默认采用最近邻插值法进行插值处理K=imrotate(I,90); % 默认旋转出界的部分不被截出subplot(1,3,1);imshow(I);subplot(1,3,2);imshow(J);subplot(1,3,3);imshow(K);% 检测旋转后的图像是否失真P=imrotate(K,270);figure;imshow(P);。
Matlab中的模糊图像恢复与图像重建技术详解引言:随着数码相机、移动设备以及各种图像处理软件的普及,人们对图像质量要求越来越高。
然而,在图像获取和传输过程中,由于种种原因,图像可能会变得模糊,失真或损坏。
为了解决这些问题,图像恢复和重建技术应运而生。
本文将详细介绍基于Matlab的模糊图像恢复与图像重建技术。
一、图像模糊恢复技术1. 模糊图像的概念和原因模糊图像是指由于摄像机移动、图像采集设备问题、环境光线等因素而导致图像失真的现象。
图像模糊会降低图像的细节和清晰度,使得图像难以辨认和识别。
常见的模糊原因有运动模糊、焦距模糊、镜头畸变等。
2. 模糊图像恢复方法为了恢复模糊图像的清晰度和细节,研究人员提出了各种方法。
其中,基于傅里叶变换的频域滤波是最常用的方法之一。
该方法通过将模糊图像转换到频域,应用适当的频域滤波器来消除模糊效果。
Matlab提供了丰富的函数和工具箱来实现这些滤波方法,比如利用低通滤波器恢复运动模糊图像。
另外,基于对图像恢复的数学建模和优化算法也是常用的方法。
例如,最小二乘法、最小化总变差等。
3. Matlab中的模糊图像恢复函数Matlab提供了多种函数用于模糊图像恢复。
其中,`deconvwnr`函数可以用于模糊图像的逆滤波处理。
该函数通过对图像进行频域滤波,去除模糊效果。
另外,`deconvblind`函数可以用于盲去卷积处理,即对图像进行反卷积操作以恢复图像细节。
二、图像重建技术1. 图像重建的意义和应用图像重建指的是利用已有的图像信息来还原、修复或生成新的图像。
与图像恢复类似,图像重建技术对于改善图像质量、还原损坏图像、生成虚拟图像等方面有着重要的应用。
图像重建技术在医学影像、图像压缩和增强、虚拟现实等领域都有广泛的应用。
2. 图像重建算法在Matlab中,图像重建可以通过多种算法实现。
其中一种常用的算法是基于插值的图像重建方法。
该方法通过对已有图像的像素进行插值来生成新的图像。
实验六 模糊图像恢复一、实验目的本实验是一个综合性实验,要求学生巩固学习多个知识点和内容,主要有: 1、理解掌握运动图像的退化模型; 2、掌握维纳滤波法的原理和实现方法;3、在不同的噪声和点扩散函数参数下进行恢复,并比较结果;4、通过分析和实验得出相应的结论。
二、实验准备1、运动模糊退化模型:运动模糊是图像退化的一种,可以用数学表达式刻画出来。
对线性移(空)不变系统,退化模型可表示为:g(x,y)=h(x,y)*f(x,y)+n(x,y)。
对匀速直线运动而言,退化图像为:()()()[]⎰--=Tdt t y y t x x f y x g 000,,其中x 0(t)和y 0(t)分别表示x 和y 方向的运动分量。
并假设退化系统是线性移不变的,光学成像过程是完善的,快门开关是瞬间完成的。
对上式进行傅立叶变换,则得频域表达式为()()()[]()()[]()[]()()()[]{}),(),(2exp ,2exp ,2exp ,,000000v u H v u F dt t vy t ux j v u F dtdxdy vy ux j t y y t x x f dxdy vy ux j y x g v u G T T=+-=⎥⎥⎦⎤⎢⎢⎣⎡+---=+-=⎰⎰⎰⎰⎰⎰+∞∞-+∞∞-+∞∞-+∞∞-πππ 其中()()()[]{}dt t vy t ux j v u H T⎰+-=0002exp ,π假设景物只在x 方向匀速运动,在T 时间内共移动距离是a ,即x 0(t)=at/T ,y 0(t)=0,则()()[]ua j ua ua T dt T at uj v u H Tππππ-=⎥⎦⎤⎢⎣⎡-=⎰exp sin 2exp ,0 在Matlab 中可用滤波器卷积的方法仿真出运动模糊图像。
h=fspecial(‘motion ’,len,theta),表示在theta 方向移动len 长度,产生运动模糊的点扩散函数h 。
图像复原MATLAB实现前⾔:本篇博客先介绍滤波器滤除噪声,再介绍滤波器复原,侧重于程序的实现。
⼀:三种常见的噪声⼆:空间域滤波空间域滤波复原是在已知噪声模型的基础上,对噪声的空间域进⾏滤波。
空间域滤波复原⽅法主要包括: 均值滤波器 算术均值滤波器 ⼏何均值滤波器 谐波均值滤波器 逆谐波均值滤波器 顺序统计滤波器 中值滤波器 最⼤值/最⼩值滤波器2.1算数均值滤波器1 img=imread('D:/picture/ZiXia.jpg');2 img=rgb2gray(img);3 figure,imshow(img);//原图4 img_noise=double(imnoise(img,'gaussian',0.06));5 figure,imshow(img_noise,[]);//含有⾼斯噪声的图6 img_mean=imfilter(img_noise,fspecial('average',3));//滤波后的图7 figure;imshow(img_mean,[]);2.2⼏何均值滤波器1 img=imread('cameraman.tif');2 img=rgb2gray(img);3 figure,imshow(img);4 img_noise=double(imnoise(img,'gaussian',0.06));5 figure,imshow(img_noise,[]);6 img_mean=exp(imfilter(log(img_noise+1),fspecial('average',3)));7 figure;imshow(img_mean,[]);2.3谐波均值滤波器2.4逆谐波均值滤波器采⽤逆谐波均值滤波器对附加胡椒噪声图像进⾏滤波的matlab程序如下:1 img=imread('cameraman.tif'); figure,imshow(img);2 [M,N]=size(img);R=imnoise2('salt & pepper',M,N,0.1,0);3 img_noise=img;img_noise(R==0)=0;4 img_noise=double(img_noise); figure,imshow(img_noise,[]);5 Q=1.5;6 img_mean=imfilter(img_noise.^(Q+1),fspecial('average',3))./imfilter(img_noise.^Q,fspecial('average',3));7 figure;imshow(img_mean,[]);采⽤逆谐波均值滤波器对附加盐噪声图像进⾏滤波的matlab程序如下:1 img=imread('csboard.tif');figure,imshow(img);2 [M,N]=size(img);R=imnoise2('salt & pepper',M,N,0,0.1);3 img_noise=img;img_noise(R==1)=255;4 img_noise=double(img_noise); figure,imshow(img_noise,[]);5 Q=-1.5;6 img_mean=imfilter(img_noise.^(Q+1),fspecial('average',3))./imfilter(img_noise.^Q,fspecial('average',3));7 figure;imshow(img_mean,[]);2.5中值滤波器1 img=imread('cameraman.tif');2 img_noise=double(imnoise(img,'salt & pepper',0.06));3 img_mean=imfilter(img_noise,fspecial('average',5));4 img_median=medfilt2(img_noise);%⼀次中值滤波5 img_median2=medfilt2(img_median);%⼆次中值滤波2.6最⼤值,最⼩值滤波器利⽤最⼤值滤波器消除胡椒噪声污染图像的matlab程序如下。
如何利用Matlab进行图像恢复图像恢复是数字图像处理中的一个重要的研究领域。
Matlab作为一种功能强大的工具,被广泛应用于图像处理领域。
本文将介绍如何利用Matlab进行图像恢复,并探讨其中的原理和算法。
首先,图像恢复是一种通过消除或减小图像失真、模糊或噪声等问题,使图像更加清晰和还原的过程。
在实际应用中,图像常常受到噪声污染、运动模糊、光照变化等影响,导致图像质量下降。
利用图像恢复技术,可以提高图像的视觉质量和辨识度,对于图像处理、计算机视觉等领域具有重要意义。
Matlab作为一款高级的数学计算工具,提供了丰富的函数库和灵活的编程接口,能够方便地进行图像处理和分析。
在图像恢复中,Matlab提供了多种处理图像的函数和算法,可以帮助我们实现各种图像恢复的方法。
一种常用的图像恢复方法是基于空域滤波的处理。
在Matlab中,可以使用imfilter函数来实现各种空域滤波算法,如均值滤波、中值滤波、高斯滤波等。
这些滤波算法通过在图像像素之间进行加权平均或统计操作,可以消除图像中的噪声和模糊。
另一种常用的图像恢复方法是基于频域滤波的处理。
在Matlab中,可以使用fft2函数和ifft2函数来实现图像的傅里叶变换和反傅里叶变换。
通过将图像从空域转换到频域,可以利用频域滤波算法对图像进行处理,如理想低通滤波、巴特沃斯低通滤波、维纳滤波等。
这些滤波算法可以根据图像的频域特征,有选择地增强或抑制图像中的某些频率分量,从而实现图像的恢复。
此外,Matlab还提供了一些专门用于图像恢复的函数,如wiener2函数、deconvwnr函数等。
wiener2函数实现了维纳滤波算法,可以用于消除运动模糊或加性噪声的图像恢复。
deconvwnr函数实现了维纳滤波的变种算法,可以根据图像和模糊函数的噪声特性,自适应地调整滤波参数,从而实现更好的图像恢复效果。
除了上述方法,Matlab还提供了其他一些高级的图像恢复算法,如超分辨率恢复、图像拼接等。
高斯模糊实现(matlab)高斯模糊是一种图像模糊滤波器,它用正态分布计算图像中每个像素的变换。
N 维空间正态分布方程为(1)在二维空间定义为(2)其中r 是模糊半径,指模板元素到模板中心的距离。
σ 是正态分布的标准偏差,。
在二维空间中,这个公式生成的曲面的等高线是从中心开始呈正态分布的同心圆。
分布不为零的像素组成的卷积矩阵与原始图像做变换。
每个像素的值都是周围相邻像素值的加权平均。
原始像素的值有最大的高斯分布值,所以有最大的权重,相邻像素随着距离原始像素越来越远,其权重也越来越小。
这样进行模糊处理比其它的均衡模糊滤波器更高地保留了边缘效果。
1. 使用给定高斯模板平滑图像维基百科的实例高斯模糊矩阵:0.00000067 0.00002292 0.00019117 0.00038771 0.00019117 0.00002292 0.00000067 0.00002292 0.00078633 0.00655965 0.01330373 0.00655965 0.00078633 0.00002292 0.00019117 0.00655965 0.05472157 0.11098164 0.05472157 0.00655965 0.00019117 0.00038771 0.01330373 0.11098164 0.22508352 0.11098164 0.01330373 0.00038771 0.00019117 0.00655965 0.05472157 0.11098164 0.05472157 0.00655965 0.00019117 0.00002292 0.00078633 0.00655965 0.01330373 0.00655965 0.00078633 0.000022920.00000067 0.00002292 0.00019117 0.00038771 0.00019117 0.00002292 0.00000067用该矩阵进行高斯模糊的结果如下:使用代码如下:guass=[0.00000067 0.00002292 0.00019117 0.00038771 0.000191 17 0.00002292 0.00000067;0.00002292 0.00078633 0.00655965 0.01330373 0.006559650.00078633 0.00002292;0.00019117 0.00655965 0.05472157 0.11098164 0.054721570.00655965 0.00019117;0.00038771 0.01330373 0.11098164 0.22508352 0.110981640.01330373 0.00038771;0.00019117 0.00655965 0.05472157 0.11098164 0.054721570.00655965 0.00019117;0.00002292 0.00078633 0.00655965 0.01330373 0.006559650.00078633 0.00002292;0.00000067 0.00002292 0.00019117 0.00038771 0.00019117 0.00002292 0.00000067];TestImg=imread('Lena1.jpg');FuzzyImg=conv2(TestImg,guass,'full');subplot(121);imshow(TestImg);subplot(122);imshow(FuzzyImg/256);编程注意事项:在matlab中,我们常使用imshow()函数来显示图像,而此时的图像矩阵可能经过了某种运算。
matlab盲去模糊算法-回复matlab盲去模糊算法是一种常用于图像处理领域的技术。
在拍摄或传输过程中,图像往往会受到模糊的影响,导致细节失真或不清晰。
盲去模糊算法可以有效地恢复原始图像的清晰度和细节。
本文将一步一步地介绍matlab盲去模糊算法的原理和实现过程。
I. 模糊图像的生成在进行盲去模糊算法之前,首先需要生成一个模糊图像。
可以使用matlab 内置的图像模糊函数,如imfilter或imgaussfilt,对原始图像进行模糊处理。
将生成的模糊图像保存为一个矩阵,作为算法的输入。
II. 图像模型建立在盲去模糊算法中,需要建立一个数学模型来描述图像模糊的过程。
一种常用的模型是卷积模型,即假设模糊过程可以由一个矩阵与原始图像的卷积得到。
设原始图像为f,模糊图像为g,模糊核为h,那么模糊图像可以表示为g = f ∗h,其中∗表示卷积运算。
III. 反卷积算法盲去模糊的核心在于反卷积算法的实现。
基本思想是通过估计模糊核h的逆矩阵h_inv,将模糊图像g进行逆卷积恢复得到清晰的图像f_hat。
由于逆卷积是一个不稳定的问题,常常会引入一些正则化项来增加稳定性,如Tikhonov正则化或L1正则化等。
IV. 选择正则化项正则化项是为了约束逆卷积问题的解,防止过拟合和噪声放大。
在选择正则化项时,可以根据实际情况和需求进行选择。
常用的正则化项有Tikhonov正则化、最小二乘正则化和L1正则化等。
这些正则化项可以用于优化问题的目标函数中,以获得更好的逆卷积结果。
V. 正则化参数选择正则化参数是用来平衡模型复杂度和拟合程度的重要参数。
在选择正则化参数时,可以使用交叉验证或模型选择方法来确定最佳参数。
交叉验证方法将数据集分为训练集和验证集,通过在训练集上进行模型训练并在验证集上进行测试,选择表现最好的参数。
VI. 算法实现和结果评估在matlab中,可以使用函数deconvblind来实现盲去模糊算法。
该函数的输入参数包括模糊图像、模糊核的大小等。
附录:源代码:1、高斯噪声的添加以及滤波处理I= imread('E: \lena.bmp','bmp');J=imnoise(I,'gaussian',0,0.01);figure;subplot(1,2,1);imshow(I);title('源图像lena.bmp');subplot(1,2,2);imshow(J);title('加入gaussian噪声后的lena.bmp');n1=7;sigma1=1.5;n2=3;sigma2=1.5;theta=0;r=[cos(theta) -sin(theta); sin(theta) cos(theta)];for i = 1 : n2for j = 1 : n1u = r*[j-(n1+1)/2 i-(n2+1)/2]';h(i,j)=exp(-u(1)^2/(2*sigma1^2))/(sigma1*sqrt(2*pi))*exp(-u(2)^2/(2*sigma2^2))/(sigma2*sqrt(2*pi)); endendh = h / sqrt(sum(sum(h.*h)));f1=conv2(J,h,'same');subplot(1,2,2);figure;imagesc(f1);title('高斯平滑后的lena.bmp(7x7)');colormap(gray);2、图像模糊及添加噪声①图像的运动模糊I= imread('E: \lena.bmp','bmp');figure;subplot(1,2,1);imshow(I);title('源图像lena.bmp');f=double(I); % 数据类型转换,MATLAB不支持图像的无符号整型的计算g=fft2(f); % 傅立叶变换g=fftshift(g); % 转换数据矩阵[M,N]=size(g);a=0.5;b=0.5;T=0.1;m=fix(M/2);n=fix(N/2);j=sqrt(-1);for i=1:Mfor k=1:Nh=(T/(pi*(i*a+k*b)))*sin(pi*(i*a+k*b))*exp(-j*pi*(i*a+k*b)); endresult(i,k)=h*g(i,k);endresult=ifftshift(result);J1=ifft(result);J2=uint8(real(J1));subplot(1,2,2);imshow(J2);title('模糊化lena.bmp');②图像加噪figure;subplot(1,2,1);imshow(J2);title('运动模糊后的lena.bmp(角度为45)');J3=imnoise(J2,'gaussian',0,0.01);subplot(1,2,2);imshow(J3);title('加噪并模糊的lena.bmp');3、图像恢复①维纳滤波恢复图像I= imread('E: \lena.bmp','bmp');H=fspecial('motion',50,45);J=imfilter(I,H,'circular','conv');figure;subplot(1,2,1);imshow(J);title('运动模糊后的lena.bmp(角度为45)');J1=imnoise(J,'gaussian',0,0.01);subplot(1,2,2);imshow(J1);title('加噪并模糊的lena.bmp');%figure;J2=deconvwnr(J1,H),[]);imshow(J2);title('模糊噪声图像的维纳滤波复原');figure;noise=imnoise(zeros(size(I)),'gaussian',0,0.01);NSR=sum(noise(:).^2)/sum(im2double(I(:)).^2);J3=deconvwnr(J1,H,NSR);imshow(J3);title('引入SNR的维纳滤波复原');②约束最小二乘方法恢复图像I= imread('E:\lena.bmp','bmp');I1=checkerboard(8);PSF=fspecial('motion',50,45);V=0.0001;J=imfilter(I,PSF,'circular','conv'); J1=imnoise(J,'gaussian',0,0.01); figure;subplot(1,2,1);imshow(J1);title('模糊加噪图像');NoisePower=V*prod(size(I));[G,LAGRA]=deconvreg(J,PSF,NoisePower); subplot(1,2,2);imshow(G);title('约束最小二乘滤波复原');。
msr颜⾊复原matlab程序,Retinex图像去模糊(含MATLAB代码)Retinex 图像去模糊Retinex 理论:Retinex是由 Retina + Cortex 组成,被称为视⽹膜⽪层理论。
该理论指出 物体能够被观察到的信息是由两个因素来决定的:物体本⾝的反射性质和物体周围的光照强度。
其中,光照强度决定了原始图像中所有像素点的动态范围⼤⼩,⽽原始图像的固有属性(颜⾊)是由物体本⾝的反射系数决定的。
因此该⽅法是将⼀幅图像分成两个不同的图像:反射图像和亮度图像,然后去除光照的影响,保留物体的固有属性。
该⽅法⽬前主要有三个分⽀:单尺度Retinex(SSR),多尺度Retinex(MSR)和带颜⾊恢复的多尺度Retinex(MSRCR)。
Retinex理论的基本思路:在原始图像中,通过某种⽅法去除或者降低⼊射图像的影响,从⽽保留物体本质的反射属性图像。
单尺度Retinex(SSR):算法思路:1 读取原始图像,分别取出RGB三个通道的数值矩阵,转化成 double 型。
2 利⽤⾼斯模板对各个通道的进⾏卷积操作处理。
3 将RGB图像以及卷积后的图像转化到对数域,并相减,然后利⽤反对数转化到 实数域。
4 对获得的各个通道的图像进⾏线性拉伸,合并RGB就可的处理后的图像。
算法实现(MATLAB):close all; clear all; clcI = imread('E:\⽔下⽬标识别\chinamm2019uw_train\jpegimage\000000.jpg');I_r = double(I(:,:,1));I_g = double(I(:,:,2));I_b = double(I(:,:,3));I_r_log = log(I_r+1);I_g_log = log(I_g+1);I_b_log = log(I_b+1);Rfft1 = fft2(I_r);Gfft1 = fft2(I_g);Bfft1 = fft2(I_b);% SSR算法[m,n] = size(I_r);sigma = 200;f = fspecial('gaussian', [m, n], sigma); efft1 = fft2(double(f));D_r = ifft2(Rfft1.*efft1);D_g = ifft2(Gfft1.*efft1);D_b = ifft2(Bfft1.*efft1);D_r_log = log(D_r + 1);D_g_log = log(D_g + 1);D_b_log = log(D_b + 1);R = I_r_log - D_r_log;G = I_g_log - D_g_log;B = I_b_log - D_b_log;R = exp(R);MIN = min(min(R));MAX = max(max(R));R = (R - MIN)/(MAX - MIN);R = adapthisteq(R);G = exp(G);MIN = min(min(G));MAX = max(max(G));G = (G - MIN)/(MAX - MIN);G = adapthisteq(G);B = exp(B);MIN = min(min(B));MAX = max(max(B));B = (B - MIN)/(MAX - MIN);B = adapthisteq(B);J = cat(3, R, G, B);figure;subplot(121);imshow(I);subplot(122);imshow(J);figure;imshow(J)结果:多尺度Retinex(MSR):由于SSR需要在颜⾊保真度和细节保持上追求⼀个完美的平衡,⽽这个平衡不宜实现。
二维高斯滤波matlab -回复二维高斯滤波在数字图像处理中扮演着重要的角色,它能够有效地去除图像中的噪声,并使图像看起来更加清晰和平滑。
在本文中,我们将详细介绍二维高斯滤波的原理、实现步骤以及如何在MATLAB中进行应用。
1. 二维高斯滤波的原理二维高斯滤波是一种线性平滑滤波器,它利用高斯函数对输入图像进行卷积运算。
高斯函数是一种钟形曲线,具有中心对称性,其形状由两个参数决定,即均值μ和方差σ。
对于二维高斯分布,其数学表达式为:G(x, y) = [1 / (2πσ^2)] * exp(-(x^2 + y^2) / (2σ^2))其中,(x, y)代表图像中的像素位置,σ表示高斯函数的标准差,决定了滤波器的平滑程度。
值得注意的是,当σ=0时,高斯函数变为一个脉冲响应,即不起到平滑的作用。
2. 二维高斯滤波的实现步骤为了将二维高斯滤波应用于图像处理,我们需要按照以下步骤进行实现:步骤1: 设定滤波器的大小和标准差在开始滤波之前,我们需要指定滤波器的大小和标准差。
滤波器的大小决定了滤波器的模板,一般为奇数,例如3x3、5x5等。
标准差σ决定了滤波器的平滑程度,较小的σ会产生较弱的平滑效果,而较大的σ则会产生更强的平滑效果。
步骤2: 生成二维高斯滤波核接下来,我们需要根据指定的滤波器大小和标准差生成二维高斯滤波核。
我们可以使用高斯函数的离散近似来实现,即将连续高斯函数转化为离散形式。
离散二维高斯滤波核的计算公式如下:G(x, y) = (1 / (2πσ^2)) * exp(-(x^2 + y^2) / (2σ^2))步骤3: 对输入图像进行卷积运算在生成二维高斯滤波核之后,我们可以将其应用于输入图像。
以图像中的每个像素为中心,将滤波器与图像进行卷积运算。
对于位于图像边缘的像素,我们可以采用镜像填充(mirror padding)等方法来解决边缘效应的问题。
步骤4: 输出滤波后的图像最后,经过卷积运算后,我们将得到滤波后的图像。
逆滤波复原matlab -回复逆滤波复原是一种信号处理方法,用于恢复被模糊处理的图像或信号。
在计算机视觉和图像处理中,图像模糊可能是由于摄像头或其他传感器的运动模糊,或者是由于图像传输过程中的噪声引起的。
通过逆滤波复原,我们可以尝试从模糊的图像中恢复原始清晰的图像。
在本文中,我将一步一步解释逆滤波复原的原理,并提供一些在MATLAB中实现逆滤波复原的示例代码。
首先,让我们了解模糊是如何发生的。
在图像处理中,模糊通常是由某种系统的传输函数引起的。
这个传输函数描述了输入信号通过系统时所发生的改变。
在传输函数的作用下,输入信号的高频部分被抑制,而低频部分被保留下来。
这样就导致了图像的模糊。
逆滤波复原的目标是根据已知的传输函数,将模糊的图像恢复为清晰的原始图像。
逆滤波复原的方法是通过对模糊图像进行频域处理来实现的。
具体来说,我们需要将模糊图像和传输函数都转换到频域,然后将它们相除,得到的结果就是恢复过程中所需的滤波器。
在MATLAB中实现逆滤波复原的第一步是将模糊图像和传输函数转换到频域。
这可以通过应用二维离散傅里叶变换(DFT)来实现。
在MATLAB 中,可以使用fft2函数将图像和传输函数转换为频域表示形式。
假设我们的模糊图像为I,传输函数为H,则可以执行以下操作:matlabI_freq = fft2(I);H_freq = fft2(H);接下来,我们需要将传输函数的频域表示进行逆滤波。
逆滤波的目标是通过将传输函数的频域表示的每个像素值除以相应位置的模糊图像的频域表示,得到恢复过程中所需的滤波器。
然而,在实际应用中,这样的直接除法运算可能会导致过度放大高频噪声的问题。
为了解决这个问题,我们可以引入一个稳定因子(也称为正则化常数)来控制恢复的强度。
代码示例如下:matlabepsilon = 0.01; 正则化常数,用于控制恢复的强度F = conj(H_freq) ./ (abs(H_freq).^2 + epsilon);这里的`F` 是逆滤波复原所需的滤波器。
1、退化程序clc;clear all;close all;I=imread('F:\mmw\B1\图2.jpg');%读图figure;subplot(4,3,1);imshow(I);title('原图像');LEN=30;%运动长度30THETA=30;%运动角度30% LEN=60;% THETA=60;n=2;for i=1:3for j=1:3PSF=fspecial('motion',LEN*i,THETA*j);%退化并研究运动角度和长度对图片模糊程度的影响PSF=fspecial('motion',LEN,THETA);Blurred=imfilter(I,PSF,'circular','conv');subplot(4,3,n);imshow(uint8(Blurred));title('模糊化');hold onn=n+1;endend%imwrite(Blurred,'模糊∠60长60.png');%保存图2、运动角度的求解%求解模糊运动角度matlab代码close all;clc;clear all;im=imread('F:\mmw\B1\模糊∠60长60.png');img_gray=rgb2gray(im);%灰度化img_fft=fftshift(fft2(img_gray));N=abs(img_fft);P=(N-min(min(N)))/(max(max(N))-min(min(N)))*225;figure;imshow(P);title('频谱图(运动角度与光斑方向垂直)');len=35;theta=0;PSF=fspecial('motion',len,theta);B=imfilter(img_gray,PSF,'circular','conv');subplot(121);imshow(B);%模糊图像B1=fft2(double(B));B2=mat2gray(log(abs(B1)));subplot(122);imshow(B2);%模糊图像的频谱图C=sum(B2,1);%对频谱图求列和[m,n]=size(C); x=0:1:n-1; y=C;figure,plot(x,y);title('频谱列和曲线图1') %绘制频谱列和曲线图3、运动长度算法%求解模糊运动长度matlab代码:im=imread('F:\mmw\B1\模糊∠60长60.png'); img_gray=rgb2gray(im);%灰度化h=fspecial('sobel');%sobel边缘检测img_double=double(img_gray);J=conv2(img_double,h,'same');IP=abs(fft2(J));S=fftshift(real(ifft2(IP)));figure;plot(S);title('模糊运动长度');4、噪声分析%噪声分析clc;clear allim=imread('F:\mmw\B1\图1.png');[m,n,h]=size(im);f11=ones(192,162,3);f22=ones(130,130,3);f33=ones(100,100,3);f44=ones(70,70,3);for i=1:190for j=1:162for k=1:3f11(i,j,k)=im(i,j,k);endendendfor i=1:130for j=501:630fork=1:3;f22(i,j-500,k)=im(i,j,k);endendfor i=721:870for j=11:170for k=1:3f33(i-720,j-10,k)=im(i,j,k);endendendfor i=761:830for j=561:630for k=1:3f33(i-760,j-560,k)=im(i,j,k);endendendfigure;subplot(221),hist(f11,100);subplot(222),hist(f22,100);subplot(223),hist(f33,100);subplot(224),hist(f44,100);title('噪声分析2');5、去噪还原clc;clear all;close all;I=imread('F:\mmw\B1\图1.png');%读图Len=60;Theta=60;PSF=fspecial('motion',Len,Theta); %模糊化BlurredA=imfilter(I,PSF,'circular','conv');wnr1=deconvwnr(BlurredA,PSF);%维纳滤波BlurredD=imfilter(I,PSF,'circ','conv');INITPSF=ones(size(PSF));[K DePSF]=deconvblind(BlurredD,INITPSF,30);%盲去卷积法BlurredB=imfilter(I,PSF,'conv');v=0.02;Blurred_I_Noisy=imnoise(BlurredB,'gaussian',0,v);NP=v*prod(size(I));J=deconvreg(Blurred_I_Noisy,PSF,NP);%最小二乘法BlurredC=imfilter(I,PSF,'symmetric','conv');v=0.002;BlurredNoisy=imnoise(BlurredC,'gaussian',0,v);Luc=deconvlucy(BlurredNoisy,PSF,5);%L_Rl滤波subplot(221);imshow(I);title('原图');subplot(222);imshow(BlurredA);title('模糊化');%subplot(233);imshow(wnr1);title('维纳滤波');subplot(223);imshow(J);title('最小二乘法');imwrite(J,'min_recover1.png');subplot(224);imshow(Luc);title('L_R法');imwrite(Luc,'LR_recover1.png');6、截取部分图片进行对比程序clear all;clc;a=imread('F:\mmw\B1\模糊∠60长60.png');%未处理质量较差图像b=a([64:120],[67:126]);a=imread('F:\mmw\min_recover1.png');%算法处理后质量较好图象c=a([64:120],[67:126]);%%从eyechart3中截取测试参考图象,截取部分需要进行缩放------------------%%使之与eyechart1,eyechart2截取部分大小匹配-----------------------a=imread('F:\mmw\B1\图2.jpg');%高清晰参考图象d=a([64:120],[67:126]);e=imresize(d,[length(b(:,1)),length(b(1,:))],'bicubic');%调整imwrite(b,'area_模糊∠60长60.png');imwrite(c,'area_最小二乘法复原图.png');imwrite(e,'area_图2.png');subplot(1,3,1);imshow(e);title('模糊∠60长60截取参考');hold on;subplot(1,3,2);imshow(b);title('eyechart1截取部分');hold on;subplot(1,3,3);imshow(c);title('eyechart2截取部分');7、模糊系数、质量指数、PSNR的计算(评价)clc;clear;PSNRenable=1; %PSNR计算使能,为0不计算,为1,计算KBlurenable=1; %模糊系数KBlur计算使能,为0不计算,为1,计算Qenable=1; %质量指数Q计算使能,为0不计算,为1,计算for m=1:2imsrcnamehead='area_模糊∠60长60'; %源图象文件名头imsrcnameext='png'; %源图象文件名扩展if m==1 %以area_eyechart1.bmp为测试图象imdstname=strcat('area_图2','.',imsrcnameext);%污染图象文件名,可修改elseif m==2%以area_eyechart2.bmp为测试图象imdstname=strcat('area_最小二乘法复原图','.',imsrcnameext);%污染图象文件名,可修改end%--------------------------------------------------------------------------iminfo=imfinfo(strcat(imsrcnamehead,'.',imsrcnameext));%源图象信息读取imsrc=imread(strcat(imsrcnamehead,'.',imsrcnameext)); %源图象读取imdst=imread(imdstname,imsrcnameext); %污染图象读取doubleimsrc=double(imsrc); %转换为浮点类型doubleimdst=double(imdst); %转换为浮点类型%----------------------------------------------------源图象和污染图象读取W=iminfo.Width; %图象宽H=iminfo.Height; %图象高%----------------------------PSNR计算--------------------------------------if PSNRenable==1PSNR=0.0; %PSNR赋初值for j=1:Hfor i=1:WPSNR=PSNR+double((doubleimsrc(j,i)-doubleimdst(j,i))*(doubleimsrc(j,i )-doubleimdst(j,i)));endendPSNR=PSNR/W/H;PSNR=10*log10(255*255/PSNR);%---------------------------PSNR计算完毕-----------------------------------end%-------------------------模糊系数KBlur计算--------------------------------if KBlurenable==1Sin=0.0; %Sin赋初值Sout=0.0;for i=2:W-1t=doubleimsrc(j-1,i+1)+doubleimsrc(j+1,i-1)-doubleimsrc(j-1,i-1)-doub leimsrc(j+1,i+1);if t<0t=-t;endSin=Sin+t; %源图象邻域边缘能量计算t=doubleimdst(j-1,i+1)+doubleimdst(j+1,i-1)-doubleimdst(j-1,i-1)-doub leimdst(j+1,i+1);if t<0t=-t;endendendSout=Sout+t; %污染图象邻域边缘能量计算KBlur=Sout/Sin;end%-------------------------------KBlur计算完毕------------------------------%-------------------------------质量指数Q计算------------------------------if Qenable==1Q=0.0; %Q赋初值Qnum=0; %图象以7X7块大小计算每块的Q,逐象素的移动块窗口,这里Qnum为块数量的计数for j=4:H-3for i=4:W-3midsrc=0.0;middst=0.0;varsrc=0.0;vardst=0.0; %源图象和污染图象块内的平均值和方差赋初值varsrcdst=0.0;%源图象和污染图象块内的协方差赋初值for n=-3:3for m=-3:3midsrc=midsrc+doubleimsrc(j+n,i+m);middst=middst+doubleimdst(j+n,i+m);endendmidsrc=midsrc/49;middst=middst/49;%%源图象和污染图象块内的平均值计算----------------------------------for n=-3:3varsrc=varsrc+(doubleimsrc(j+n,i+m)-midsrc)*(doubleimsrc(j+n,i+m)-mid src);vardst=vardst+(doubleimdst(j+n,i+m)-middst)*(doubleimdst(j+n,i+m)-mid dst);varsrcdst=varsrcdst+(doubleimsrc(j+n,i+m)-midsrc)*(doubleimdst(j+n,i+ m)-middst);endendvarsrc=varsrc/48;vardst=vardst/48;varsrcdst=varsrcdst/48;if((varsrc+vardst)*(midsrc*midsrc+middst*middst))~=0%分母不为零的块才计算质量指数QQ=Q+4*varsrcdst*midsrc*middst/((varsrc+vardst)*(midsrc*midsrc+middst* middst));%%源图象和污染图象块内Q计算完毕------------------------------------Qnum=Qnum+1; %块计数加1endendendQ=Q/Qnum;endendQPSNRKBlur。
clear all% source = aviread('C:\Video\Source\traffic\san_fran_traffic_30sec_QVGA');source = mmreader('SampleVideo.avi');frameQYT=get(source,'NumberOfFrames');% ----------------------- frame size variables -----------------------fr = read(source,1); % 读取第一帧作为背景fr_bw = rgb2gray(fr); % 将背景转换为灰度图像fr_size = size(fr); %取帧大小width = fr_size(2);height = fr_size(1);fg = zeros(height, width);bg_bw = zeros(height, width);% --------------------- mog variables -----------------------------------C = 4; % 组成混合高斯的单高斯数目(一般3-5)M = 0; % 组成背景的数目D = 2.5; % 阈值(一般2.5个标准差)alpha = 0.01; % learning rate 学习率决定更新速度(between 0 and 1) (from paper 0.01)thresh = 0.75; % foreground threshold 前景阈值(0.25 or 0.75 in paper)sd_init = 6; % initial standard deviation 初始化标准差(for new components) var = 36 in paperw = zeros(height,width,C); % initialize weights array 初始化权值数组mean = zeros(height,width,C); % pixel means 像素均值sd = zeros(height,width,C); % pixel standard deviations 像素标准差u_diff = zeros(height,width,C); % difference of each pixel from mean 与均值的差p = alpha/(1/C); % initial p variable 参数学习率(used to update mean and sd)rank = zeros(1,C); % rank of components (w/sd)% ------initialize component means and weights 初始化均值和权值----------pixel_depth = 8; % 8-bit resolution 像素深度为8位pixel_range = 2^pixel_depth -1; % pixel range 像素范围2的7次方0—255(# of possible values)for i=1:heightfor j=1:widthfor k=1:Cmean(i,j,k) = rand*pixel_range; % means random (0-255之间的随机数)w(i,j,k) = 1/C; % weights uniformly distsd(i,j,k) = sd_init; % initialize to sd_initendendend%----- process frames -处理帧--,这里去第八帧n = 8;fr = read(source,n); % read in frame 读取帧fr_bw = rgb2gray(fr); % convert frame to grayscale 转换为灰度图像% calculate difference of pixel values from mean 计算像素差值for m=1:Cu_diff(:,:,m) = abs(double(fr_bw) - double(mean(:,:,m)));end% update gaussian components for each pixel 更新每个像素的背景模型for i=1:heightfor j=1:widthmatch = 0;for k=1:Cif (abs(u_diff(i,j,k)) <= D*sd(i,j,k)) % pixel matches component像素匹配了模型match = 1; % variable to signal component match 设置匹配记号% update weights, mean, sd, p 更新权值,均值,标准差和参数学习率w(i,j,k) = (1-alpha)*w(i,j,k) + alpha;p = alpha/w(i,j,k);mean(i,j,k) = (1-p)*mean(i,j,k) + p*double(fr_bw(i,j));sd(i,j,k) = sqrt((1-p)*(sd(i,j,k)^2) + p*((double(fr_bw(i,j)) - mean(i,j,k)))^2);else % pixel doesn't match component 几个模型中都没有匹配的w(i,j,k) = (1-alpha)*w(i,j,k); % weight slighly decreases 权值减小endendbg_bw(i,j)=0;for k=1:Cbg_bw(i,j) = bg_bw(i,j)+ mean(i,j,k)*w(i,j,k); %更新背景if(bg_bw(i,j)>thresh)k=k-1;M=k;end% 这里有问题,背景权值和大于阈值时,背景建模的数目M取k-1,end% if no components match, create new component 如果没有匹配的模型则创建新模型if (match == 0)[min_w, min_w_index] = min(w(i,j,:));mean(i,j,min_w_index) = double(fr_bw(i,j));sd(i,j,min_w_index) = sd_init;endrank = w(i,j,:)./sd(i,j,:); % calculate component rank 计算模型优先级rank_ind = [1:1:C];% calculate foreground 计算前景while ((match == 0)&&(k>M))% 这里用于前景计算的高斯模型应该是C-M,所以这里k>Mif (abs(u_diff(i,j,rank_ind(k))) <= D*sd(i,j,rank_ind(k)))fg(i,j) = 0; %black = 0elsefg(i,j) = fr_bw(i,j);endk = k+1;if(k==5)k=k-1;breakendendendendfigure(1),subplot(3,1,1),imshow(fr) %显示输入图像subplot(3,1,2),imshow(uint8(bg_bw)) %显示背景图像subplot(3,1,3),imshow(uint8(fg)) %显示前景图像。
%%通过BP神经网络训练将模糊化的图像进行还原处理%%下面是主程序image=imread('image.bmp');image1=rgb2gray(image);figure,subplot(121),imshow(image);title('RGB原图'); subplot(122),imshow(image1);title('灰度图');image2=imresize(image1,[128,128]);%读入图像,并改变成90*90image2=double(image2)./256;len=5;theta=10;PSF=fspecial('motion',len,theta);%加入运动模糊Blurredmotion=imfilter(image2,PSF,'circular','conv' );w=fspecial('gaussian',[9 9],1);%加入高斯模糊blurred_image=imfilter(Blurredmotion,w);P_Matrix=zeros(32,(size(image2,1))*(size(image2,2)) /16);%输入矩阵T_Matrix=zeros(16,(size(image2,1))*(size(image2,2))/16);%目标矩阵blurred_BW = edge(blurred_image,'sobel');%提取模糊图像边缘BW = edge(image2,'sobel');%提取清晰图像边缘image_edge=double(BW).*image2;blurred_image_edge=double(blurred_BW).*blurred_imag e;%blurred_image_flat=blurred_image-blurred_image_edg e;%得到模糊图像非边缘区域,即平坦区域%target_image=zeros(size(image2,1),size(image2,2));%目标矩阵t=1;%初始化输入矩阵P_Matrix,假设神经网络是32-30-16的结构,故按如下方式初始化for i=1:4:size(image2,1)for j=1:4:size(image2,2)P_Matrix(1,t)=blurred_image(i,j);%前16位为degraded image的block,大小为4x4P_Matrix(2,t)=blurred_image(i,j+1);P_Matrix(3,t)=blurred_image(i,j+2);P_Matrix(4,t)=blurred_image(i,j+3);P_Matrix(5,t)=blurred_image(i+1,j);P_Matrix(6,t)=blurred_image(i+1,j+1);P_Matrix(7,t)=blurred_image(i+1,j+2);P_Matrix(8,t)=blurred_image(i+1,j+3);P_Matrix(9,t)=blurred_image(i+2,j);P_Matrix(10,t)=blurred_image(i+2,j+1);P_Matrix(11,t)=blurred_image(i+2,j+2);P_Matrix(12,t)=blurred_image(i+2,j+3);P_Matrix(13,t)=blurred_image(i+3,j);P_Matrix(14,t)=blurred_image(i+3,j+1);P_Matrix(15,t)=blurred_image(i+3,j+2);P_Matrix(16,t)=blurred_image(i+3,j+3);P_Matrix(17,t)=blurred_image_edge(i,j);%后16位为original image's edge 的block,大小为4x4P_Matrix(18,t)=blurred_image_edge(i,j+1);P_Matrix(19,t)=blurred_image_edge(i,j+2);P_Matrix(20,t)=blurred_image_edge(i,j+3);P_Matrix(21,t)=blurred_image_edge(i+1,j);P_Matrix(22,t)=blurred_image_edge(i+1,j+1);P_Matrix(23,t)=blurred_image_edge(i+1,j+2);P_Matrix(24,t)=blurred_image_edge(i+1,j+3);P_Matrix(25,t)=blurred_image_edge(i+2,j);P_Matrix(26,t)=blurred_image_edge(i+2,j+1);P_Matrix(27,t)=blurred_image_edge(i+2,j+2);P_Matrix(28,t)=blurred_image_edge(i+2,j+3);P_Matrix(29,t)=blurred_image_edge(i+3,j);P_Matrix(30,t)=blurred_image_edge(i+3,j+1);P_Matrix(31,t)=blurred_image_edge(i+3,j+2);P_Matrix(32,t)=blurred_image_edge(i+3,j+3);t=t+1;endendt=1;%初始化目标矩阵for i=1:4:size(image2,1)for j=1:4:size(image2,2)T_Matrix(1,t)=image2(i,j);T_Matrix(2,t)=image2(i,j+1);T_Matrix(3,t)=image2(i,j+2);T_Matrix(4,t)=image2(i,j+3);T_Matrix(5,t)=image2(i+1,j);T_Matrix(6,t)=image2(i+1,j+1);T_Matrix(7,t)=image2(i+1,j+2);T_Matrix(8,t)=image2(i+1,j+3);T_Matrix(9,t)=image2(i+2,j);T_Matrix(10,t)=image2(i+2,j+1);T_Matrix(11,t)=image2(i+2,j+2);T_Matrix(12,t)=image2(i+2,j+3);T_Matrix(13,t)=image2(i+3,j);T_Matrix(14,t)=image2(i+3,j+1);T_Matrix(15,t)=image2(i+3,j+2);T_Matrix(16,t)=image2(i+3,j+3);t=t+1;endendnet=newff(P_Matrix,T_Matrix,30);%得到神经网络,假设使用的是BP神经网络net=train(net,P_Matrix,T_Matrix);Y=sim(net,P_Matrix);%Y得到复原矩阵restored_image=zeros(size(image2,1),size(image2,2)) ;%复原图像t=1;%通过复原矩阵得到复原图像for i=1:4:size(image2,1)for j=1:4:size(image2,2)restored_image(i,j)=Y(1,t);restored_image(i,j+1)=Y(2,t);restored_image(i,j+2)=Y(3,t);restored_image(i,j+3)=Y(4,t);restored_image(i+1,j)=Y(5,t);restored_image(i+1,j+1)=Y(6,t);restored_image(i+1,j+2)=Y(7,t);restored_image(i+1,j+3)=Y(8,t);restored_image(i+2,j)=Y(9,t);restored_image(i+2,j+1)=Y(10,t);restored_image(i+2,j+2)=Y(11,t);restored_image(i+2,j+3)=Y(12,t);restored_image(i+3,j)=Y(13,t);restored_image(i+3,j+1)=Y(14,t);restored_image(i+3,j+2)=Y(15,t);restored_image(i+3,j+3)=Y(16,t);t=t+1;endendPSNR_edge_BP=Cal_PSNR(image2*255,restored_image*255 );SSIM_edge_BP=ssim_index(image2*255,restored_image*2 55);PSNR_edge_BP_blurred=Cal_PSNR(image2*255,blurred_im age*255);SSIM_edge_BP_blurred=ssim_index(image2*255,blurred_ image*255);figure,%显示结果subplot(221),imshow(image2);title('原图像');subplot(222),imshow(blurred_image);title('运动模糊后的图像');subplot(223),imshow(restored_image);title('BP神经网络复原图像');imwrite(image2,'image.bmp');imwrite(blurred_image,'blurred_image.bmp');imwrite(restored_image,'restored_image.bmp');image=imread('image.bmp');image1=rgb2gray(image);figure,subplot(121),imshow(image);title('RGB原图'); subplot(122),imshow(image1);title('灰度图');image2=imresize(image1,[128,128]);%读入图像,并改变成90*90image2=double(image2)./256;len=5;theta=10;PSF=fspecial('motion',len,theta);%加入运动模糊Blurredmotion=imfilter(image2,PSF,'circular','conv' );w=fspecial('gaussian',[9 9],1);%加入高斯模糊blurred_image=imfilter(Blurredmotion,w);P_Matrix=zeros(32,(size(image2,1))*(size(image2,2)) /16);%输入矩阵T_Matrix=zeros(16,(size(image2,1))*(size(image2,2)) /16);%目标矩阵blurred_BW = edge(blurred_image,'sobel');%提取模糊图像边缘BW = edge(image2,'sobel');%提取清晰图像边缘image_edge=double(BW).*image2;blurred_image_edge=double(blurred_BW).*blurred_imag e;%blurred_image_flat=blurred_image-blurred_image_edg e;%得到模糊图像非边缘区域,即平坦区域%target_image=zeros(size(image2,1),size(image2,2));%目标矩阵t=1;%初始化输入矩阵P_Matrix,假设神经网络是32-30-16的结构,故按如下方式初始化for i=1:4:size(image2,1)for j=1:4:size(image2,2)P_Matrix(1,t)=blurred_image(i,j);%前16位为degraded image的block,大小为4x4P_Matrix(2,t)=blurred_image(i,j+1);P_Matrix(3,t)=blurred_image(i,j+2);P_Matrix(4,t)=blurred_image(i,j+3);P_Matrix(5,t)=blurred_image(i+1,j);P_Matrix(6,t)=blurred_image(i+1,j+1);P_Matrix(7,t)=blurred_image(i+1,j+2);P_Matrix(8,t)=blurred_image(i+1,j+3);P_Matrix(9,t)=blurred_image(i+2,j);P_Matrix(10,t)=blurred_image(i+2,j+1);P_Matrix(11,t)=blurred_image(i+2,j+2);P_Matrix(12,t)=blurred_image(i+2,j+3);P_Matrix(13,t)=blurred_image(i+3,j);P_Matrix(14,t)=blurred_image(i+3,j+1);P_Matrix(15,t)=blurred_image(i+3,j+2);P_Matrix(16,t)=blurred_image(i+3,j+3);P_Matrix(17,t)=blurred_image_edge(i,j);%后16位为original image's edge 的block,大小为4x4P_Matrix(18,t)=blurred_image_edge(i,j+1);P_Matrix(19,t)=blurred_image_edge(i,j+2);P_Matrix(20,t)=blurred_image_edge(i,j+3);P_Matrix(21,t)=blurred_image_edge(i+1,j);P_Matrix(22,t)=blurred_image_edge(i+1,j+1);P_Matrix(23,t)=blurred_image_edge(i+1,j+2);P_Matrix(24,t)=blurred_image_edge(i+1,j+3);P_Matrix(25,t)=blurred_image_edge(i+2,j);P_Matrix(26,t)=blurred_image_edge(i+2,j+1);P_Matrix(27,t)=blurred_image_edge(i+2,j+2);P_Matrix(28,t)=blurred_image_edge(i+2,j+3);P_Matrix(29,t)=blurred_image_edge(i+3,j);P_Matrix(30,t)=blurred_image_edge(i+3,j+1);P_Matrix(31,t)=blurred_image_edge(i+3,j+2);P_Matrix(32,t)=blurred_image_edge(i+3,j+3);t=t+1;endendt=1;%初始化目标矩阵for i=1:4:size(image2,1)for j=1:4:size(image2,2)T_Matrix(1,t)=image2(i,j);T_Matrix(2,t)=image2(i,j+1);T_Matrix(3,t)=image2(i,j+2);T_Matrix(4,t)=image2(i,j+3);T_Matrix(5,t)=image2(i+1,j);T_Matrix(6,t)=image2(i+1,j+1);T_Matrix(7,t)=image2(i+1,j+2);T_Matrix(8,t)=image2(i+1,j+3);T_Matrix(9,t)=image2(i+2,j);T_Matrix(10,t)=image2(i+2,j+1);T_Matrix(11,t)=image2(i+2,j+2);T_Matrix(12,t)=image2(i+2,j+3);T_Matrix(13,t)=image2(i+3,j);T_Matrix(14,t)=image2(i+3,j+1);T_Matrix(15,t)=image2(i+3,j+2);T_Matrix(16,t)=image2(i+3,j+3);t=t+1;endendnet=newff(P_Matrix,T_Matrix,30);%得到神经网络,假设使用的是BP神经网络net=train(net,P_Matrix,T_Matrix);Y=sim(net,P_Matrix);%Y得到复原矩阵restored_image=zeros(size(image2,1),size(image2,2)) ;%复原图像t=1;%通过复原矩阵得到复原图像for i=1:4:size(image2,1)for j=1:4:size(image2,2)restored_image(i,j)=Y(1,t);restored_image(i,j+1)=Y(2,t);restored_image(i,j+2)=Y(3,t);restored_image(i,j+3)=Y(4,t);restored_image(i+1,j)=Y(5,t);restored_image(i+1,j+1)=Y(6,t);restored_image(i+1,j+2)=Y(7,t);restored_image(i+1,j+3)=Y(8,t);restored_image(i+2,j)=Y(9,t);restored_image(i+2,j+1)=Y(10,t);restored_image(i+2,j+2)=Y(11,t);restored_image(i+2,j+3)=Y(12,t);restored_image(i+3,j)=Y(13,t);restored_image(i+3,j+1)=Y(14,t);restored_image(i+3,j+2)=Y(15,t);restored_image(i+3,j+3)=Y(16,t);t=t+1;endendPSNR_edge_BP=Cal_PSNR(image2*255,restored_image*255 );SSIM_edge_BP=ssim_index(image2*255,restored_image*2 55);PSNR_edge_BP_blurred=Cal_PSNR(image2*255,blurred_im age*255);SSIM_edge_BP_blurred=ssim_index(image2*255,blurred_ image*255);figure,%显示结果subplot(221),imshow(image2);title('原图像');subplot(222),imshow(blurred_image);title('运动模糊后的图像');subplot(223),imshow(restored_image);title('BP神经网络复原图像');imwrite(image2,'image.bmp');imwrite(blurred_image,'blurred_image.bmp');imwrite(restored_image,'restored_image.bmp');function[mssim, ssim_map] = ssim_index(img1, img2, K, window, L)% % 结构相似度(SSIM)用于图像复原(去噪,去模糊,盲复原,修复)后的客观评价!% % 该评价标准比较适合于纹理结构多的图像质量评价if (nargin < 2 | nargin > 5) %判断输入变量的个数mssim = -Inf; %设定为负的无穷大ssim_map = -Inf;return;endif (size(img1) ~= size(img2)) %两个图像的大小不相等mssim = -Inf;ssim_map = -Inf;return;end[M N] = size(img1);if (nargin == 2)if ((M < 11) | (N < 11))mssim = -Inf;ssim_map = -Inf;returnendwindow = fspecial('gaussian',11,1.5); %统计的本地窗口默认是高斯窗口,模板尺寸为11,滤波器标准值为1.5像素K(1) = 0.01; %默认的设置K(2) = 0.03;L = 255; %图像的动态范围endif (nargin == 3)if ((M < 11) | (N < 11))mssim = -Inf;ssim_map = -Inf;returnendwindow = fspecial('gaussian',11,1.5);L = 255;if(length(K) == 2) %SSIM指标公式常数的长度为2if (K(1) < 0 | K(2) < 0)mssim = -Inf;ssim_map = -Inf;return;endelsemssim = -Inf;ssim_map = -Inf;return;endendif (nargin == 4)[H W] = size(window);if ((H*W) < 4 | (H > M) | (W > N)) mssim = -Inf;ssim_map = -Inf;returnendL = 255;if (length(K) == 2)if (K(1) < 0 | K(2) < 0)mssim = -Inf;ssim_map = -Inf;return;endelsemssim = -Inf;ssim_map = -Inf;return;endendif (nargin == 5)[H W] = size(window);if ((H*W) < 4 | (H > M) | (W > N)) mssim = -Inf;ssim_map = -Inf;returnendif (length(K) == 2)if (K(1) < 0 | K(2) < 0)mssim = -Inf;ssim_map = -Inf;return;endelsemssim = -Inf;ssim_map = -Inf;return;endendC1 = (K(1)*L)^2;C2 = (K(2)*L)^2;window = window/sum(sum(window));img1 = double(img1); %将图像矩阵uint8型数据转换成double型img2 = double(img2);% % 计算图片的平均值和方差--------------------------------------------------mu1= filter2(window, img1, 'valid'); %使用矩阵h中的二维FIR(有限长单位冲%激响应)过滤器过滤变量X中的数据,使用二维相关性计算结果Y,并且返回与X的维数相%同的相关项的中心部分。
高斯模糊图像的盲复原
顾亚芳
【期刊名称】《科教文汇》
【年(卷),期】2008(000)005
【摘要】介绍了高斯模糊图像的盲复原方法.通过实验比较了6种去噪方法:采用最大似然估计方法来寻找最相似于退化图像的点扩展函数;利用已估计出的点扩展函数进行约束最小二乘图像复原,并分析了正则化参数的选取方法.最后比较了EM算法和本文算法.
【总页数】1页(P74-74)
【作者】顾亚芳
【作者单位】昆山第一职业高级中学,江苏·昆山,215316
【正文语种】中文
【中图分类】G632
【相关文献】
1.盲复原高斯模糊图像
2.逆主元法盲目复原高斯模糊图像
3.高斯模糊图像的复原处理与研究
4.复原处理大噪声高斯型模糊图像的前置滤波器的设计研究
5.基于低秩矩阵与稀疏约束的运动模糊图像盲复原
因版权原因,仅展示原文概要,查看原文内容请购买。
东南大学
硕士学位论文
高斯模糊图像的盲复原姓名:顾亚芳
申请学位级别:硕士专业:信号与信息处理指导教师:吴乐南
20051225
东南大学硕士学位论文
3.3.2实验结果
取噪声方差占2=O.01,则各种去噪方法所复原的图像示于图3.1。
可以看出,维纳去噪的视觉效果最好,二次均值其次,而小波变换去噪的效果对于高斯噪声而言并不是很好。
图3.1去噪图像的主观对比
表3.1去噪图像均方误差的比较
62=0.00162=0.Ol62=O.162=1.0含噪图像457.23581.0103e+0034.7299e+0031.3682e_卜004
小波变换1.1296e+0091.1316e+0091.1201e+0091.0779e+009
维纳滤波421.9284479.0220515.48232.5666e+003
维纳+均值1.1098e-卜0091.1117叶0091.0894e+0091.0284e+009
~次均值1.1161e+0091.1225e+0091.1250e+0091.1044e+009
二次均值1.1064e+0091.1109e+0091.0999e+0091.0546e+009
三次均值1.0986e+0091.1025e+0091.0878e+0091.0357e_卜009各种方法去嵘后的图像均方误差如表3.1所示。
比较来看,对于高斯噪声,无论噪声大小,维纳滤波的去噪效果总是最好的
3.4小结
本章通过实验比较了6种去除图像噪声的方法:一次、二次和三次均值滤波,一次维纳去噪,一次维纳组合一次均值去噪,以及基于独立自适应阈值的小波去噪。
主观评判和基于复员图像均方误差的客观比较都表明,对于高斯噪声,维纳滤波的去噪效果在这6种方法中效果最好。
东南大学硕士学位论文
4.2.4实验结果
取噪声方差万2=O.01,图4.2给出了预处理分别采用3.3节的6种去噪方法后所得到的点扩展函数的主观对比;而表4.1则给出其均方误差的对比。
图4.26种去噪方法采用EM算法获取的点扩展函数对比
表4.16种方法左噪后图像所获得的PSF均方误筹的比较
62=0.00162=0.0162=0.162=1.0含噪图像1.7876e.0040.01100.31200.9588
小波变换1.1590e-0041.2369e.0041.225le一0045.9282e.004
维纳滤波1.2009e-0041.1220e.0041.9583e一0041.9226e-004
维纳+均值8.7668e一0058.2955e.0055.2122e.0053.7933e.005
一次均值7.8996e.0055.7468e.0056.9729e.0040.0046
二次均值7.0210e.0055.4030e.0057.7613e一0055.4512e.004
三次均值0.00470.00460.02790.5049
从复原后的PSF图像和PSF的均方误差可以看出:在噪卢较大时,采用维纳加均值组合去噪滤波后所得到的点扩展函数最接近于真实的PSF;但当噪声较小时,采用二次均值滤波所获得的点扩展函数最接近于真实的PSF。
本节只是估计点扩散函数,而不需要追求复原图像效果,所以迭代次数只设定为lO次,这样可以减少复原时间。
(4.39)其中MED表示带噪图像信号小波变换_,=l尺度上的小波系数对角矩阵的中值。
(2)以肛旷一gll2=11.112为约束条件,最小化8j矿一刎2+∥I阿112,求得正则化参数∥,进一步求得正则解.厂。
(3)通常白动确定复原滤波器比人为调整滤波器参数的复原结果要差,特别是约束展小二乘滤波器完全由单一标量参数来决定时更是如此。
所以在初步确定Ⅳ以后,再根据视觉效果人为地进行调整。
4.3.3复原实验结果
仍取噪声方差万2=0.01,分别采用3.3节6种去噪方法获取点扩展函数后进行复原所得图像如图4.2。
复原图像均方误差的比较则如表4.2所示。
图4.26种去噪方法获得点扩展函数后的复原图像对比
表4.2采用6种之噪方法获得的点扩展函数所复原图像的均方误若比较
62=0.00162=o.Ol62=0.162=1.0含噪图像457.608310lO.35089.413534
小波变换301.6989393.4022769.07453956.7
维纳滤波301.9000393.4436768.90363956.7
维纳+均值301.1769393.37lO768.49233956.7
一次均值301.0759393.2766771.49523956.7
二次均值301.0759393.2766771.49523956.7
三次均值32J.561l405.40159】J.85】03956.7
从图4.2最后所获得的复原图像和表4.2中所列出的均方误差对比可以看到,对于高斯模糊图像当噪声比较小时,采用一次均值滤波去噪后所获得的图像最好,随着噪卢增加,维纳加均值组合去噪后的复原效果最好,但当噪声大到己淹没图像时,任何一种去噪方法的效果都很差,即对图像复原已无能为力了。
第4章高斯模糊图像的盲复原
4.4总体性能比较
4.4.1EM算法存在的问题
4.2节介绍了EM算法的原理,本节分析一下EM算法目前还存在的问题和需要改进之处:
1)首先可以从以上的算法分析中看出,此算法非常繁琐,要得到复原结果,计算时间非常长;
2)在此复原算法中,由于点扩展函数的系数之间没有固定的联系,因此需要估计大量的独立参数来获得点扩展函数,图像尺寸太大会使计算过程复杂,计算时间延长,因此,算法只适用于小尺寸的图像;
3)在此复原算法的分析过程中,为了推导相关公式,需要假设支持域的尺寸已知,实际情况是两者往往未知,在计算过程中,需假定支持域的范围,但这可能比实际的点扩展函数偏大或偏小,会带来不可避免的图像复原误差;
4)EM算法可能收敛到局部最优点,而不是全局最优点。
这是由于初始猜测点的不同,使得算法的收敛值有所不同。
4.4.2性能对比
仍取噪声方差艿2=0.001,图4.3是EM算法和本文算法所复原图像主观效果的对比;而表4.3则是
EM算法和本文算法所复原图像均方误差的对比。
图4.3EM算法和本文算法复原图像的对比
表4.3两种算法复原图像均方误差对比
62=O.00152=O.01
含噪图像455.87641010.O
EM算法1.1564e+0091.1839e+009
本文算法304.5428398.2778
从视觉图像和均方误差可看出,在噪声较小时,EM算法对图像的去模糊效果也很好,但也由于迭代带来了噪声的放大。
而在噪声比较大时,EM算法的复原效果已很差。
但采用本文算法无论是从视觉I型像评价还是从均方误差对比,都能获得较好的复原效果。
4.5小结
本章在第3章去噪预处理的基础上,利用EM算法获取点扩展函数后进行约束最小二乘法图像复原,最后再调节正则化参数以达到更佳的复原效果。
经过实验对比得出结论:
高斯模糊图像的盲复原
作者:顾亚芳
学位授予单位:东南大学
1.期刊论文顾亚芳高斯模糊图像的盲复原-科教文汇2008,""(5)
介绍了高斯模糊图像的盲复原方法.通过实验比较了6种去噪方法:采用最大似然估计方法来寻找最相似于退化图像的点扩展函数;利用已估计出的点扩展函数进行约束最小二乘图像复原,并分析了正则化参数的选取方法.最后比较了EM算法和本文算法.
本文链接:/Thesis_Y1040222.aspx
授权使用:武汉理工大学(whlgdx),授权号:3e6de5ec-a33a-4e85-811a-9da200d4cca2,下载时间:2010年6月
27日。