灰度图像加性噪声污染和运动模糊图像复原

关键字

代码

I=imread('6.bmp');
noise=0.18*randn(size(I));
psf=fspecial('motion',21,11);

blurred=imfilter(I,psf,'circular');
>> blurrednoisy=im2uint8(blurred+noise);
??? Error using ==> plus
Integers can only be combined with integers of the same class, or scalar doubles.

>>

figure,imshow(noise);


figure,imshow(blurred);
>>

实验和结果情况:

image

image

>> size(I)

ans =

   512   512     3

>> info=imfinfo('6.bmp');
>> info

info =

              Filename: '6.bmp'
           FileModDate: '13-May-2009 15:26:56'
              FileSize: 786486
                Format: 'bmp'
         FormatVersion: [1x33 char]
                 Width: 512
                Height: 512
              BitDepth: 24
             ColorType: 'truecolor'
       FormatSignature: 'BM'
    NumColormapEntries: 0
              Colormap: []
               RedMask: []
             GreenMask: []
              BlueMask: []
       ImageDataOffset: 54
      BitmapHeaderSize: 40
             NumPlanes: 1
       CompressionType: 'none'
            BitmapSize: 786432
        HorzResolution: 0
        VertResolution: 0
         NumColorsUsed: 0
    NumImportantColors: 0

E1=mat2gray(I(:,:,1));
noise=0.1*randn(size(E1));
psf1=fspecial('motion',21,11);
blurred1=imfilter(E1,psf1,'circular');
burrednoisy1=im2uint8(blurred1+noise);
 E2=mat2gray(I(:,:,2));
blurred2=imfilter(E2,psf1,'circular');
burrednoisy2=im2uint8(blurred2+noise);
E3=mat2gray(I(:,:,3));
blurred3=imfilter(E3,psf1,'circular');
burrednoisy3=im2uint8(blurred3+noise);
burrednoisy(:,:,1)=burrednoisy1;
burrednoisy(:,:,2)=burrednoisy2;
burrednoisy(:,:,3)=burrednoisy3;

figure,imshow(burrednoisy);

                                              这样就实现效果如下:

image

NP=abs(fftn(noise)).^2;                         %噪声功率
NPOW=sum(NP(:))/prod(size(noise));    %噪声自相关
NCORR=fftshift(real(ifftn(NP)));

IP1=abs(fftn(E1)).^2;                           %R分量图像功率
IPOW1=sum(IP1(:))/prod(size(E1));      %R分量图像自相关
IPORR1=fftshift(real(ifftn(IP1)));

IPORRC1=IPORR1( :,ceil(size(E1,1)/2));

IP2=abs(fftn(E2)).^2;                                         %G分量图像功率
IPOW2=sum(IP2(:))/prod(size(E2));                  %G分量图像自相关
IPORR2=fftshift(real(ifftn(IP2)));

IPORRC2=IPORR2( :,ceil(size(E2,1)/2));

IP3=abs(fftn(E3)).^2;                                          %B分量图像功率
IPOW1=sum(IP3(:))/prod(size(E3));                   %B分量图像自相关
IPORR3=fftshift(real(ifftn(IP3)));

IPORRC3=IPORR3( :,ceil(size(E3,1)/2));

q(:,:,1)=deconvwnr(burrednoisy1,psf,NCORR,IPORR1);

q(:,:,2)=deconvwnr(burrednoisy2,psf,NCORR,IPORR2);

q(:,:,3)=deconvwnr(burrednoisy3,psf,NCORR,IPORR3);

figure,imshow(q);

title('deconvwnr(burrednoisy,psf,NCORR,IPORR)')

v(:,:,1)=deconvwnr(burrednoisy1,psf,NCORR,IPORRC1);

v(:,:,2)=deconvwnr(burrednoisy2,psf,NCORR,IPORRC2);

v(:,:,3)=deconvwnr(burrednoisy3,psf,NCORR,IPORRC3);

figure,imshow(v);

title('deconvwnr(burrednoisy,psf,NCORR,IPORRC)')

实现效果如下:

image

image

遗憾的是,当移植到一幅彩色图像时候,出错,如下

错误点1

>> E1=mat2gray(I(:,:,1));
>> noise=0.1*randn(size(E1));
>> psf1=fspecial('motion',21,11);
>> blurred1=imfilter(E1,psf1,'circular');
>> burrednoisy1=im2uint8(blurred1+noise);
>> E2=mat2gray(I(:,:,2));
>> blurred2=imfilter(E2,psf1,'circular');
>> burrednoisy2=im2uint8(blurred2+noise);
>> E3=mat2gray(I(:,:,3));
>> blurred3=imfilter(E3,psf1,'circular');
>> burrednoisy3=im2uint8(blurred3+noise);
>> burrednoisy(:,:,1)=burrednoisy1;
??? Subscripted assignment dimension mismatch.

info=imfinfo('fly.bmp')

info =

              Filename: 'fly.bmp'
           FileModDate: '16-May-2009 10:16:42'
              FileSize: 856374
                Format: 'bmp'
         FormatVersion: [1x33 char]
                 Width: 640
                Height: 446
              BitDepth: 24
             ColorType: 'truecolor'
       FormatSignature: 'BM'
    NumColormapEntries: 0
              Colormap: []
               RedMask: []
             GreenMask: []
              BlueMask: []
       ImageDataOffset: 54
      BitmapHeaderSize: 40
             NumPlanes: 1
       CompressionType: 'none'
            BitmapSize: 856320
        HorzResolution: 0
        VertResolution: 0
         NumColorsUsed: 0
    NumImportantColors: 0

可以导致错误2

>> figure,imshow(burrednoisy1);
>> burrednoisy1=mat2gray(burrednoisy1);
>> burrednoisy(:,:,1)=burrednoisy1;
??? Subscripted assignment dimension mismatch. 所以,这里问题真是很大很大,搞不明白为什么这里会出错

可以导致错误3:

>>  I=imread('fly.bmp');
>> E1=I(:,:,1);
>> figure,imshow(E1);
>> noise=0.1*randn(size(E1));
>> psf1=fspecial('motion',21,11);
>> blurred1=imfilter(E1,psf1,'circular');
>> burrednoisy1=im2uint8(blurred1+noise);
??? Error using ==> plus
Integers can only be combined with integers of the same class, or scalar doubles.

测试分析:

>> size(burrednoisy2)

ans =

   446   640

>> size(blurred2)

ans =

   446   640

>>

>> burrednoisy=zeros(size(I));
>> size(burrednoisy)

ans =

   446   640     3

>>

                                                                  结果是ok:

                                                                     >> burrednoisy(:,:,1)=burrednoisy1;
                                                                     >>

                                                                       恐怕上述问题在于缺少了这个burrednoisy=zeros(size(I));

可是,运行结果还是问题大得不可理喻:


I=imread('fly.bmp');
burrednoisy=zeros(size(I));
q=zeros(size(I));
E1=mat2gray(I(:,:,1));
noise=0.1*randn(size(E1));
psf1=fspecial('motion',21,11);
blurred1=imfilter(E1,psf1,'circular');
burrednoisy1=im2uint8(blurred1+noise);
E2=mat2gray(I(:,:,2));
blurred2=imfilter(E2,psf1,'circular');
burrednoisy2=im2uint8(blurred2+noise);
E3=mat2gray(I(:,:,3));
blurred3=imfilter(E3,psf1,'circular');
burrednoisy3=im2uint8(blurred3+noise);
burrednoisy(:,:,1)=burrednoisy1;
burrednoisy(:,:,2)=burrednoisy2;
burrednoisy(:,:,3)=burrednoisy3;

figure,imshow(burrednoisy);

NP=abs(fftn(noise)).^2;                         %噪声功率
NPOW=sum(NP(:))/prod(size(noise));    %噪声自相关
NCORR=fftshift(real(ifftn(NP)));

IP1=abs(fftn(E1)).^2;                           %R分量图像功率
IPOW1=sum(IP1(:))/prod(size(E1));      %R分量图像自相关
IPORR1=fftshift(real(ifftn(IP1)));

IPORRC1=IPORR1( :,ceil(size(E1,1)/2));

IP2=abs(fftn(E2)).^2;                                         %G分量图像功率
IPOW2=sum(IP2(:))/prod(size(E2));                  %G分量图像自相关
IPORR2=fftshift(real(ifftn(IP2)));

IPORRC2=IPORR2( :,ceil(size(E2,1)/2));

IP3=abs(fftn(E3)).^2;                                          %B分量图像功率
IPOW1=sum(IP3(:))/prod(size(E3));                   %B分量图像自相关
IPORR3=fftshift(real(ifftn(IP3)));

IPORRC3=IPORR3( :,ceil(size(E3,1)/2));

q(:,:,1)=deconvwnr(burrednoisy1,psf,NCORR,IPORR1);

q(:,:,2)=deconvwnr(burrednoisy2,psf,NCORR,IPORR2);

q(:,:,3)=deconvwnr(burrednoisy3,psf,NCORR,IPORR3);

figure,imshow(q);

title('deconvwnr(burrednoisy,psf,NCORR,IPORR)')

实验2

>> clear
>> I=imread('fly.bmp');
>> E1=mat2gray(I(:,:,1));
>> burrednoisy(:,:,1)=E1;
>> E2=mat2gray(I(:,:,2));
>> burrednoisy(:,:,2)=E2;
>> E3=mat2gray(I(:,:,3));
>> burrednoisy(:,:,3)=E3;
>> figure,imshow(burrednoisy);

image

实验2 :

>> I=imread('fly.bmp');
>> E1=mat2gray(I(:,:,1));
>> noise=0.1*randn(size(E1));
>>  psf1=fspecial('motion',21,11);
blurred1=imfilter(E1,psf1,'circular');
burrednoisy1=im2uint8(blurred1+noise);
>>  burrednoisy(:,:,1)=burrednoisy1;

成功了

天不负有心人啊,终于搞出个名堂;

I=imread('fly.bmp');
E1=mat2gray(I(:,:,1));
noise=0.1*randn(size(E1));
psf1=fspecial('motion',21,11);
blurred1=imfilter(E1,psf1,'circular');
burrednoisy1=im2uint8(blurred1+noise);
burrednoisy(:,:,1)=burrednoisy1;
E2=mat2gray(I(:,:,2));
blurred2=imfilter(E2,psf1,'circular');
burrednoisy2=im2uint8(blurred2+noise);
burrednoisy(:,:,2)=burrednoisy2;
E3=mat2gray(I(:,:,3));
blurred3=imfilter(E3,psf1,'circular');
burrednoisy3=im2uint8(blurred3+noise);
burrednoisy(:,:,3)=burrednoisy3; 
figure,imshow(burrednoisy);

区别仅仅在于

burrednoisy(:,:,x)=burrednoisyx位置

image

>> NP=abs(fftn(noise)).^2;                         %噪声功率
NPOW=sum(NP(:))/prod(size(noise));    %噪声自相关
NCORR=fftshift(real(ifftn(NP)));

IP1=abs(fftn(E1)).^2;                           %R分量图像功率
IPOW1=sum(IP1(:))/prod(size(E1));      %R分量图像自相关
IPORR1=fftshift(real(ifftn(IP1)));

>> q(:,:,1)=deconvwnr(burrednoisy1,psf,NCORR,IPORR1);
??? Undefined function or variable 'psf'.

>> q(:,:,1)=deconvwnr(burrednoisy1,psf1,NCORR,IPORR1);
>> IP2=abs(fftn(E2)).^2;                                         %G分量图像功率
IPOW2=sum(IP2(:))/prod(size(E2));                  %G分量图像自相关
IPORR2=fftshift(real(ifftn(IP2)));
>> q(:,:,2)=deconvwnr(burrednoisy2,psf1,NCORR,IPORR2);
>> IP3=abs(fftn(E3)).^2;                                          %B分量图像功率
IPOW1=sum(IP3(:))/prod(size(E3));                   %B分量图像自相关
IPORR3=fftshift(real(ifftn(IP3)));
>> q(:,:,3)=deconvwnr(burrednoisy3,psf1,NCORR,IPORR3);
>> figure,imshow(q);

title('deconvwnr(burrednoisy,psf,NCORR,IPORR)')

这个结果有点怪,不知对不对.

image

对第三个图像处理,悟出刚才一些问题:为什么会出现 ??? Subscripted assignment dimension mismatch.

解决方法就是一句一句写入,或者用m文件

实验2  比较                   fspecial('motion',x,y);

>> clear
>> I=imread('audi.jpg');
>> E1=mat2gray(I(:,:,1));
>> noise=0.1*randn(size(E1));
>> psf1=fspecial('motion',21,11);
>> blurred1=imfilter(E1,psf1,'circular');
>> burrednoisy1=im2uint8(blurred1+noise);
>> burrednoisy(:,:,1)=burrednoisy1;
>> E2=mat2gray(I(:,:,2));
>> blurred2=imfilter(E2,psf1,'circular');
>> burrednoisy(:,:,2)=burrednoisy2;
??? Undefined function or variable 'burrednoisy2'.

>> burrednoisy2=im2uint8(blurred2+noise);
>> burrednoisy(:,:,2)=burrednoisy2;
>> E3=mat2gray(I(:,:,3));
>> blurred3=imfilter(E3,psf1,'circular');
>> burrednoisy3=im2uint8(blurred3+noise);
>> burrednoisy(:,:,3)=burrednoisy3;
>> figure,imshow(burrednoisy);
>> psf1=fspecial('motion',41,31);
>> blurred1=imfilter(E1,psf1,'circular');
>> burrednoisy1=im2uint8(blurred1+noise);
>> burrednoisy(:,:,1)=burrednoisy1;
>> E2=mat2gray(I(:,:,2));
>> blurred2=imfilter(E2,psf1,'circular');
>> burrednoisy2=im2uint8(blurred2+noise);
>> burrednoisy(:,:,2)=burrednoisy2;
>> E3=mat2gray(I(:,:,3));
>> blurred3=imfilter(E3,psf1,'circular');
>> burrednoisy3=im2uint8(blurred3+noise);
>> burrednoisy(:,:,3)=burrednoisy3;
>> figure,imshow(burrednoisy);

audi

 image

image

image

非常遗憾,仍然有大量细节无法修复出来.

幸运的是,及时上网得到一幅非常不错的图片

http://www-rohan.sdsu.edu/doc/matlab/toolbox/images/deblurr6.html 一个不错的网站,效果好是因为没有加噪声

deblur4a

>> I = imread('audi.jpg');
>> figure;imshow(I);title('Original Image');
>> % create PSF
LEN = 31;
THETA = 11;
PSF = fspecial('motion',LEN,THETA);

% blur the image
Blurred = imfilter(I,PSF,'circular','conv');
figure; imshow(Blurred);title('Blurred Image');

% deblur the image
wnr1 = deconvwnr(Blurred,PSF);
figure;imshow(wnr1);
title('Restored, True PSF');
>>

image

image

image

最后,为了验证噪声的作用到底有多大,我们可以用灰度图像做两个实验,一个用后一种没有噪声进行,一个用书上的,看看效果怎样,比较就知道了


【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/fleetwgx/p/1489161.html