matlab自带函数-盲卷积-加噪-卷积-滤波-小结

        总结自网上、matlab帮助文档等，都是图像复原中经常用到的基础函数或操作。可以模拟图像降质过程和用一些经典方法盲解卷积复原的过程。

一、卷积：conv2、convn、convmtx2

卷积的计算步骤：

（1）卷积核绕自己的核心元素顺时针旋转180度

（2）移动卷积核的中心元素，使它位于输入图像待处理像素的正上方

（3）在旋转后的卷积核中，将输入图像的像素值作为权重相乘

（4）第三步各结果的和做为该输入像素对应的输出像素

相关的计算步骤：

（1）移动相关核的中心元素，使它位于输入图像待处理像素的正上方

（2）将输入图像的像素值作为权重，乘以相关核

（3）将上面各步得到的结果相加做为输出

        可以看出他们的主要区别在于计算卷积的时候，卷积核要先做旋转。而计算相关过程中不需要旋转相关核。

例如： magic(3) =[8 1 6;3 5 7;4 9 2]，旋转180度后就成了[2 9 4;7 5 3;6 1 8]。

A = rand(3); 

B = rand(4); 

C = conv2(A,B)  % C is 6-by-6

C =

   0.1838  0.2374  0.9727  1.2644  0.7890  0.3750

   0.6929  1.2019  1.5499  2.1733  1.3325  0.3096

   0.5627  1.5150  2.3576  3.1553  2.5373  1.0602

   0.9986  2.3811  3.4302  3.5128  2.4489  0.8462

   0.3089  1.1419  1.8229  2.1561  1.6364  0.6841

   0.3287  0.9347  1.6464  1.7928  1.2422  0.5423

Cs = conv2(A,B,'same')   % Cs is the same size as A: 3-by-3

Cs =

   2.3576  3.1553  2.5373

   3.4302  3.5128  2.4489

   1.8229  2.1561  1.6364

二、滤波：filter2、fspecial、imfilter、spfilt、admedin

1、imfilter 可进行多维图像（RGB等）进行空间滤波，且可选参数较多
2、filter2 只能对二维图像（灰度图）进行空间滤波

3、Fspecial创建滤波算子，使用fpecial对图像模糊、退化函数建模，使用imfilter来创建一个已知PSF计算得到PSF的。

fspecial函数用于建立预定义的滤波算子，其语法格式为：

        h = fspecial(type)

        h = fspecial(type，para)

        其中type指定算子的类型，para指定相应的参数；

        type的类型有：

(1)'average'

        averaging filter

        为均值滤波，参数为hsize代表模板尺寸，默认值为【3，3】。

        H = FSPECIAL('average',HSIZE);

(2)'disk'

        为圆形区域均值滤波，参数为radius代表区域半径，默认值为5.

        H = FSPECIAL('disk',RADIUS) 

(3)'gaussian'

        为高斯低通滤波，有两个参数，hsize表示模板尺寸，默认值为【3 3】，sigma为滤波器的标准值，单位为像素，默认值为0.5.

        H = FSPECIAL('gaussian',HSIZE,SIGMA)

(4)'laplacian' 

        为拉普拉斯算子，参数alpha用于控制算子形状，取值范围为【0，1】，默认值为0.2.

       H = FSPECIAL('laplacian',ALPHA) 

(5)'log'

      为拉普拉斯高斯算子，有两个参数，hsize表示模板尺寸，默认值为【3 3】，sigma为滤波器的标准差，单位为像素，默认值为0.5.

      H = FSPECIAL('log',HSIZE,SIGMA)

(6)'motion'

      为运动模糊算子，有两个参数，表示摄像物体逆时针方向以theta角度运动了len个像素，len的默认值为9，theta的默认值为0；

      H = FSPECIAL('motion',LEN,THETA)

(7)'prewitt'

    用于边缘增强，大小为【3 3】，无参数

    H = FSPECIAL('prewitt')

(8)'sobel'

    用于边缘提取，无参数

    H = FSPECIAL('sobel') 

(9)'unsharp'

    为对比度增强滤波器。参数alpha用于控制滤波器的形状，范围为【0，1】，默认值为0.2。

    H = FSPECIAL('unsharp',ALPHA)

4、spfilt空间噪声滤波器，空间滤波器的名称有算术平均，几何平均，调和平均，反调和平均，中值，最大，最小，中点。顺序-平衡均值他们都被封装在spfilt中可以直接实用此函数来实现滤波器。

5、admedin自适应空间滤波器，是基于自适应中值滤波算法的基础上用admedin函数能够实现。

    如果输出像素是输入像素邻域像素的线性组合则称为线性滤波（例如最常见的均值滤波和高斯滤波），否则为非线性滤波（中值滤波、边缘保持滤波等）。

    线性平滑滤波器去除高斯噪声的效果很好，且在大多数情况下，对其它类型的噪声也有很好的效果。线性滤波器使用连续窗函数内像素加权和来实现滤波。特别典型的是，同一模式的权重因子可以作用在每一个窗口内，也就意味着线性滤波器是空间不变的，这样就可以使用卷积模板来实现滤波。如果图像的不同部分使用不同的滤波权重因子，且仍然可以用滤波器完成加权运算，那么线性滤波器就是空间可变的。任何不是像素加权运算的滤波器都属于非线性滤波器．非线性滤波器也可以是空间不变的，也就是说，在图像的任何位置上可以进行相同的运算而不考虑图像位置或空间的变化。

    滤波通常是用卷积或者相关来描述，而线性滤波一般是通过卷积来描述的。

三、加噪：imnoise、 imnoise2、imnoise3 

    Imnoise加噪，共有六类噪声可以使用。使用指定的分布产生空间随机噪声。

    空间噪声值是随机数，以PDF和CDF来表征。随机变量分为均匀，高斯，椒盐，对数正态，瑞利，指数，厄兰。他们都被封装在了imnoise2中的函数中。

    图像的周期噪声通常是通过频域滤波来处理的它的噪声，模型是二维正弦波它和它的M*N ，DET被封装在了imnoise3函数中，周期噪声的参数一般是通过分析图像的傅立叶频谱来估计的。通过样本图像来估计均值和方差，然后利用这些估计来求解噪声PDF的参数。噪声参数通常直接由带噪声的图像和一组带噪声的图像来估计。所选择的方法是在一幅图像中选择一个尽可能与背景一样无特色的区域，以便使该区域的亮度值的可变性主要有噪声产生。用roipoly函数寻则一个感兴趣的区域，然后再通过绘制所选区域的直方图与使用imnoise2产生的数据进行比较，被覆盖区域的均值和方差是通过statmoments函数计算得到。

    通过频域滤波降低周期噪声，滤除这些成分的主要途径是通过陷滤波。n阶巴特沃兹陷滤波器函数。

四、盲复原：deconvwnr、deconvreg、deconvlucy、deconvblind  

1、deconvwnr：维纳滤波线性图像复原方法，维纳滤波器寻找一个统计误差函数。

%读图

I = im2double(imread('cameraman.tif'));

imshow(I);

title('Original Image (courtesy of MIT)');



%模糊

LEN = 21;

THETA = 11;

PSF = fspecial('motion', LEN, THETA);

blurred = imfilter(I, PSF, 'conv', 'circular');

figure, imshow(blurred)

%加噪

noise_mean = 0;

noise_var = 0.0001;

blurred_noisy = imnoise(blurred, 'gaussian', ...noise_mean, noise_var);

figure, imshow(blurred_noisy)

title('Simulate Blur and Noise')



%无噪声复原

estimated_nsr = 0;

wnr2 = deconvwnr(blurred_noisy, PSF, estimated_nsr);

figure, imshow(wnr2)

title('Restoration of Blurred, Noisy Image Using NSR = 0')



%加入性噪比复原

estimated_nsr = noise_var / var(I(:));

wnr3 = deconvwnr(blurred_noisy, PSF, estimated_nsr);

figure, imshow(wnr3)

title('Restoration of Blurred, Noisy Image Using Estimated NSR');



2、deconvreg：约束的最小二乘方滤波，也是实现线性复原的方法。deblur image using regularized filter（利用正则化滤波复原图像 ）

%读图

I = checkerboard(8);

%滤波加噪

PSF = fspecial('gaussian',7,10);

V = .01;

BlurredNoisy = imnoise(imfilter(I,PSF),'gaussian',0,V);

%噪声功率

NOISEPOWER = V*prod(size(I));

[J LAGRA] = deconvreg(BlurredNoisy,PSF,NOISEPOWER);

 

subplot(221); imshow(BlurredNoisy);

title('A = Blurred and Noisy');

subplot(222); imshow(J);

title('[J LAGRA] = deconvreg(A,PSF,NP)');

subplot(223); imshow(deconvreg(BlurredNoisy,PSF,[],LAGRA/10));

title('deconvreg(A,PSF,[],0.1*LAGRA)');

subplot(224); imshow(deconvreg(BlurredNoisy,PSF,[],LAGRA*10));

title('deconvreg(A,PSF,[],10*LAGRA)');

3、deconvlucy：使用LR算法的迭代非线性复原，基于模型的最大似然数方程

%读图

I = checkerboard(8);

%滤波加噪

PSF = fspecial('gaussian',7,10);

V = .0001;

BlurredNoisy = imnoise(imfilter(I,PSF),'gaussian',0,V);

 

WT = zeros(size(I));

WT(5:end-4,5:end-4) = 1;

J1 = deconvlucy(BlurredNoisy,PSF);

J2 = deconvlucy(BlurredNoisy,PSF,20,sqrt(V));

J3 = deconvlucy(BlurredNoisy,PSF,20,sqrt(V),WT);

 

subplot(221);imshow(BlurredNoisy);

title('A = Blurred and Noisy');

subplot(222);imshow(J1);

title('deconvlucy(A,PSF)');

subplot(223);imshow(J2);

title('deconvlucy(A,PSF,NI,DP)');

subplot(224);imshow(J3);

title('deconvlucy(A,PSF,NI,DP,WT)');

 

4、deconvblind盲去卷积，不以PSF知识为基础的图像复原方法统称为盲去卷积算法。

%读图

I = checkerboard(8);

%滤波加噪

PSF = fspecial('gaussian',7,10);

V = .0001;

BlurredNoisy = imnoise(imfilter(I,PSF),'gaussian',0,V);

 

WT = zeros(size(I));

WT(5:end-4,5:end-4) = 1;

INITPSF = ones(size(PSF));

[J P] = deconvblind(BlurredNoisy,INITPSF,20,10*sqrt(V),WT);

 

subplot(221);imshow(BlurredNoisy);

title('A = Blurred and Noisy');

subplot(222);imshow(PSF,[]);

title('True PSF');

subplot(223);imshow(J);

title('Deblurred Image');

subplot(224);imshow(P,[]);

title('Recovered PSF');

五、实例

1、模糊加噪滤波

%加入高斯噪声

noise_mean = 0;

noise_var = 0.001;

I1=imnoise(I,'gaussian',noise_mean, noise_var);

noise=imnoise(zeros(size(I)),'gaussian',0,0.001);

%点扩展函数，sigma=0.5，模版为3*3

%PSF = fspecial('gaussian', 3, 0.5);

PSF=fspecial('motion',5,5); 

%信噪比

%I2=im2double(I);

%estimated_nsr = noise_var / var(I2(:));

 

sn=abs(fft2(noise)).^2;% 噪声功率谱 

na=sum(sn(:))/prod(size(noise));% 噪声的平均功率 

sf=abs(fft2(I)).^2;%输入图像的功率谱

fa=sum(sf(:))/prod(size(I)); %输入图像的平均功率

estimated_nsr=na/fa;%噪信比

%维纳滤波

I3=deconvwnr(I1,PSF,estimated_nsr);

2、第二个例子

subplot(2,2,2); 

INITPSF = ones(size(PSF));

[J P]= deconvblind(Blurred,INITPSF,30);

imshow(J); 

title('Restored Image');

title('Restored PSF');

 

hsize=[3,3];

H = fspecial('average',hsize);

smoothed = imfilter(I,H,'replicate');

subplot(2,2,3); 

imshow(smoothed); 

title('3*3 Smoothed Image');

 

hsize=[5,5];

H = fspecial('average',hsize);

smoothed = imfilter(I,H,'replicate');

subplot(2,2,4); 

imshow(smoothed); 

title('5*5 Smoothed Image');