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


    关键字

    代码

    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

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


  • 相关阅读:
    排序
    阿里实习 电面,面试
    外部排序
    error: LNK 2019 无法解析的外部符号
    c++ 变量定义 的初始化规则
    Spring 让 LOB 数据操作变得简单易行
    深入浅出JMS(四)--Spring和ActiveMQ整合的完整实例
    Spring jdbc call oralce procedure or function
    JSON 简单封装
    extjs3 用户管理 页面
  • 原文地址:https://www.cnblogs.com/fleetwgx/p/1489161.html
Copyright © 2020-2023  润新知