图像复原
6.1 图像复原及退化模型
6.2 噪声模型
6.3 空间域滤波复原
6.4 频率域滤波复原
6.5 估计退化函数
6.6 逆滤波和维纳滤波

一、图像复原与图像增强的区别
图像的增强是一个主观的过程，其目的是改善图片的质量，对感兴趣的部分加以增强，对不感兴趣的部分予以抑制。而图像复原是一个客观的过程，针对质量降低或失真的图像，试图恢复其原始的内容或质量。复原技术是面向退化模型的，并且采用相反的过程进行处理，以便恢复出原图像。在进行图像复原之前要先建立起其退化模型，根据该模型进行图像复原。
课本中的图像退化过程建模为一个退化函数和一个加性噪声。

假设H是一个线性时不变的过程，则我们可以得到，式子中的“*”表示卷积。其频率域的表示为：。
对于只有加性噪声的情况，我们可以通过一些噪声模型（例如高斯噪声，瑞利噪声，椒盐噪声等等）以及对这些噪声参数的估计，来选择合适的空间滤波器（如均值滤波器，中值滤波器）或者频率滤波器（带阻/带通滤波器，低通/高通滤波器，陷波滤波器）来进行滤波。
二、退化函数的估计
在图像复原中，我们需要对退化函数进行估计。主要有观察法，实验法，数学建模法。
观察法通过选择噪声较小的子图像（减少噪声的影响）来得到H(u,v)，然后根据此信息来构建全图的H(u,v)，之后利用后面的复原方法来复原。实验法是指我们可以使用或设计一个与图像退化过程相似的装置（过程），使其成像一个脉冲，可得到退化系统的冲激响应 H(u,v) = G(u,v) / A。而建模估计则是从引起图像退化的基本原理进行推导，进而对原始图像进行模拟，在模拟过程中调整模型参数以获得尽可能精确的退化模型。课本中有两种模型，大气湍流模型和运动模糊模型。
大气湍流模型：
通用形式为。其中，k是与湍流性质有关的常数，k越大，图像越模糊，与高斯低通滤波器有着相同的形式。在实现该滤波器的过程中，由于中心化，要注意u,v应该分别替换为各自与频率中心之差，假设频率中心为(M/2,N/2)，则替换为u-M/2和v-N/2。与我们之前频率滤波器的实现相同。
其傅里叶频谱如下：

注意：实现过程与我们之前的频率滤波器的实现一样，都需要中心化。
我在对图片施加大气湍流退化模型的时候，采用与之前相同的过程，但是得到的结果图像比课本的结果图更加地模糊。k = 0.0025（较模糊，比较容易看出区别）的情况如下：

这个问题困扰了我很久，我尝试了各种不同的方法，在不中心变换的情况下进行滤波，发现情况更加糟糕。后来在做后面运动模糊模型的时候，得到的结果也与课本的结果相差甚远。之后，在尝试了不对图像进行填充（即不对原图像填充0值至大小为[2m, 2n]），发现得到的结果与课本的一致。对得到的结果图像直接取整（若是在最开始将值归一化到[0,1],自然不需要取整这一步），没有归一化到[0, 255]则与课本的结果一致。代码执行结果如下（同样是k = 0.0025的情况）：

代码如下：
大气湍流实现函数为atmoTurbulence，输入为初始图像和大气湍流模型的系数k，输出为进行滤波之后的图像，下面给出的代码为了与课本的结果保持一致，没有在开始的时候对原图像进行填充，但结果仍然归一化到了[0, 255]。
function [image_out, g] = atmoTurbulence(image_in, k)
% 模拟受大气湍流影响得到的图像
% 输入为初始灰度图像，大气湍流模型的系数k
% 输出为受到大气湍流影响的结果图像
% 该实现方法没有对原图像进行填充。

[m, n] = size(image_in);
fp = zeros(m, n);
% 中心化
for i = 1 : m
    for j = 1 : n
        fp(i, j) = double(image_in(i,j)) * (-1)^(i+j);
    end
end
% 傅里叶变换
F1 = fft2(fp);
% 变换中心
p = m / 2 + 1.0;
q = n / 2 + 1.0;
% 生成大气湍流模型函数
H = zeros(m, n);
for u = 1 : m
    for v = 1 : n
        temp1 = (u-p)^2+(v-q)^2;
        H(u,v) = exp((-k)*(temp1^(5/6)));
    end
end
% 滤波
G = H .* F1;
% 反傅里叶变换，并取实部
gp = real(ifft2(G));
% 反中心化

g = zeros(m, n);
for i = 1 : m
    for j = 1 : n
        g(i, j) = gp(i, j) * (-1)^(i+j);
    end
end
% 归一化输出图像到[0, 255],g取uint8的结果与课本一致
image_out = zeros(m, n, 'uint8');
mmax = max(g(:));
mmin = min(g(:));
range = mmax-mmin;
for i = 1 : m
    for j = 1 : n
        image_out(i,j) = uint8(255 * (g(i, j)-mmin) / range);
    end
end

end
测试代码就不给出了，与之前的代码对函数的调用方法一致。
运动模糊模型：
退化函数为
其中T表示曝光时间，a和b分别表示水平和垂直方向上的移动量。注意有，当π(ua+vb) = 0时，H(u,v) = T。同理，在实现该滤波器的过程中，由于中心化，要用u-M/2和v-N/2分别替换u和v。
当a = b  = 0.1，T = 1时，其傅里叶频谱如下：
imshow(log(abs(H) + 1), [])

在对图像进行运动模糊退化之后，由于图像的运动，图像的尺寸与位置会发生一些的变化。在取a = b = 0.1，T = 1的情况下，不知道应该裁剪的数值，我直接显示了不进行裁剪的图像（由于对原图像填充0至图像大小为[2m, 2n]，在最后的结果需要裁剪），如下所示。但是，发现对图像进行裁剪之后，与课本给出的图像相差很大。

课本给出的结果如下：

经过对课本图片的仔细观察，发现图片的左边和上边应该输入右边和下边移动过来的。因此猜测课本的图片应该是没有对原图像进行填0扩充得到的。修改代码，得到的结果如下：

代码如下：
clear all;
clc;
close all;

image = imread('7.tif'); % 原图

[m, n] = size(image);
P = m / 2 + 1.0;
Q = n / 2 + 1.0;
fp = zeros(m, n);
%对图像乘以(-1)^(x+y) 以移到变换中心
for i = 1 : m
    for j = 1 : n
        % fp(i, j) = double(image(i, j));
        fp(i, j) = double(image(i, j)) * (-1)^(i+j);
    end
end
% 对填充后的图像进行傅里叶变换
F1 = fft2(fp);

a = 0.1;
b = 0.1;
T = 1;

% 生成运动模型函数，中心在(P,Q)
Mo = zeros(m, n);
for u = 1 : m
    for v = 1 : n
        temp1 = pi * ((u-P)*a + (v-Q)*b);
        if (temp1 == 0)
            Mo(u,v) = T;
        else
            Mo(u, v) =  T * sin(temp1) / temp1 * exp(-1i * (temp1));
        end
    end
end
%进行滤波
G = F1 .* Mo;

% 反傅里叶变换
gp = ifft2(G);

% 处理得到的图像
image_out = zeros(m, n, 'uint8');
gp = real(gp);
g = zeros(m, n);
for i = 1 : m
    for j = 1 : n
        g(i, j) = gp(i, j) * (-1)^(i+j);
    end
end
mmax = max(g(:));
mmin = min(g(:));
range = mmax-mmin;
for i = 1 : m
    for j = 1 : n
        image_out(i,j) = uint8(255 * (g(i, j)-mmin) / range);
    end
end
imshow(image_out);
三、逆滤波
1. 直接逆滤波
知道了图像的退化模型之后，最简单的复原方法就是采用逆滤波的方法。计算如下：
                                         
但是，当H(u,v)为0或者接近0的时候，N(u,v)/H(u,v)会变得很大，成为支配整个图像的主要部分。即使N(u,v)很小或者为0，当H(u,v)接近0的时候，在计算的时候，F(u,v)也会变得非常大，从而使复原之后的图像没有任何信息。
因此，可以对H(u,v)进行修改，①可以对G(u,v)/H(u,v)应用一个低通滤波器，虑去其中病态的高频成分（即虑去H(u,v)中接近0的部分）,②或者规定一个值，当|H(u,v)| ≤ δ 时，1/H(u,v) = 0。这两种方法都是去掉或者说削弱H(u,v)接近0时的影响。低通滤波器的方法对于上面提到的运动模糊模型几乎没有用，因为运动模糊模型的傅里叶频谱并不是从中间向四周赋值逐渐减小的，但是方法②仍然有用。
我在对收到了大气湍流模型之后的图像做直接逆滤波的时候，得到的最终图片一片模糊，得不到任何信息，与课本的结果相差较大，这可能是在对计算G(u,v)/H(u,v)的时候，在H(u,v)接近0的时候所做的处理不一致所导致的。这个直接造成了在之后对G(u,v)/H(u,v)使用butterworth低通滤波器的时候，最佳效果的截止频率与课本的不一致。
2. 维纳滤波（最小均方差误差滤波）
一个维纳滤波的计算公式如下，其中K是一个特定的常数，与噪声和未退化图像之间的信噪比有关。不同情况下取值不同。在平均意义下是最优的。

3. 约束最小二乘方滤波器 
，其中r是一个参数，P(u,y)是矩阵p(x,y) = [0, -1, 0; -1, 4,-1; 0, -1, 0]的傅里叶变换。该算法相比于维纳滤波，对其应用的每幅图像都能产生最优的结果。
4. 代码实现
下面给出在不填充图片的情况下，利用运动模型对图片进行模糊(a = b = 0.1, T = 1)，添加高斯噪声(均值为0，方差为500)，对有噪声的运动模糊图像直接使用逆滤波，使用维纳滤波（经过试验，k取0.02较好），使用约束最小二乘方滤波（经过试验，r取0.5较好）得到的结果。
代码如下（下面的代码再进行傅里叶变换之前没有对图像进行填充）：
clear all;
close all;
clc;

image = imread('7.tif');
[m, n] = size(image);
% 参数如下：
% p,q为频率中心，a,b,T为运动模糊参数，
% k为维纳滤波参数， r为约束最小二乘方参数
% M,V 分别为高斯噪声的均值和方差
p = m / 2 + 1.0;
q = n / 2 + 1.0;
a = 0.1;
b = 0.1;
T = 1;
k = 0.05;
r = 0.5;
M = 0;
V = 500;

% 对原图像进行傅里叶变换以及中心化
fp = double(image);
F = fftshift(fft2(fp));

% 生成运动模型的傅里叶变换，中心在(p,q)
Mo = zeros(m, n);
for u = 1 : m
    for v = 1 : n
        temp = pi * ((u-p)*a + (v-q)*b);
        if (temp == 0)
            Mo(u,v) = T;
        else
            Mo(u, v) =  T * sin(temp) / temp * exp(-1i * (temp));
        end
    end
end

% 生成维纳滤波的傅里叶变换
Wiener = (abs(Mo).^2) ./ (abs(Mo).^2 + k) ./ Mo;

% 生成约束最小二乘方滤波的傅里叶变换
lp = [0, -1, 0; -1, 4, -1; 0, -1, 0];
Flp = fftshift(fft2(lp, m, n));
Hw = conj(Mo) ./ (abs(Mo).^2 + r * abs(Flp).^2);

% 生成高斯噪声的傅里叶变换
noise = 500^0.5 * randn([m, n]);
Fn = fftshift(fft2(noise));

% 运动模糊图像，并且加上高斯噪声
image1 = zeros(m, n, 'uint8');
G1 = F .* Mo + Fn;
gp1 = ifft2(fftshift(G1));
g1 = real(gp1);
% 归一化图像到 [0, 255];
mmax = max(g1(:));
mmin = min(g1(:));
range = mmax-mmin;
for i = 1 : m
    for j = 1 : n
        image1(i,j) = uint8(255 * (g1(i, j)-mmin) / range);
    end
end

% 为了接近真实情况，对归一化之后的加噪图像进行逆滤波
F2 = fftshift(fft2(image1));
% 直接逆滤波
G2 = F2 ./ Mo;
gp2 = ifftshift(G2);
g2 = real(gp2);

% 维纳滤波
G3 = F2 .* Wiener;
gp3 = ifft2(fftshift(G3));
g3 = real(gp3);

% 约束最小二乘方滤波
G4 = F2 .* Hw;
gp4 = ifft2(fftshift(G4));
g4 = real(gp4);
subplot(2,3,1), imshow(image), title('原图');
subplot(2,3,2), imshow(image1), title('运动加噪图像');
subplot(2,3,3), imshow(g2, []), title('直接逆滤波');
subplot(2,3,4), imshow(g3, []), title('维纳滤波');
subplot(2,3,5), imshow(g4, []), title('约束最小二乘方滤波');
执行结果如下：

破损图像复原
简介
实验工具介绍-MATLAB
实验步骤
实验及结论

简介
因摄像机与物体相对运动，以及系统误差、畸变、噪声等因素的影响，使图像往往不是真是景物的完善映像，因而图像复原技术就此产生；在图像恢复中，需要建立造成图像质量下降的退化模型，然后运用相反过程来恢复原来图像，并运用一定准则来判定是否得到图像的最佳恢复。图像复原的原理：使用可以精确描述失真的点扩散函数（PSF）对模糊图像进行去卷积计算。图像恢复的方法主要有逆滤波复原、未纳滤波复原、约束的最小二乘方滤波复原。Lucy-Richardson滤波复原及盲去卷积滤波复原。
实验工具介绍-MATLAB
本次实验选取MATLAB作为实验工具，MATLAB是一款数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言(如C、Fortran)的编辑模式，代表了当今国际科学计算软件的先进水平。
实验步骤
选择如下图1破损图片进行复原：





图1
实验及结论
一、灰度变换增强——用调整图像的灰度增加图像的对比度

灰度变换增强是把图像的对比度从弱变强的过程，所以灰度变换增强也通常被称为对比度增强。它主要利用点运算来修正像素灰度，由输入像素灰度点的值确定相应输出点的灰度值，是一种基于图像变换的操作。


使用imadjust函数调整灰度
   J=imadjust(I);
   J=imadjust(I,[low_in;high_in],[low_out;high_out]);
   J=imadjust(I,[low_in;high_in],[low_out;high_out],gamma);



注：I为输入图像，J是返回调整后的图像，该函数把在[low_in;high_in]的像素值调整到[low_out;high_out]，对于3中的gamma，当gamma为1时，为线性变换；当gamma大于1时，此映射偏重更低数值（灰暗）输出；当gamma小于1时，此映射偏重更高数值（明亮）输出。

a)  设置J=imadjust(I,[],[0.4 0.6]); 调整图像的灰度到0.4-0.6范围内；效果图如下图2所示：



图2
b)  设置J=imadjust(I,[],[0.2 0.6]); 调整图像的灰度到0.2-0.6范围内；效果图如下图3所示：

图3
c)  设置J=imadjust(I,[],[0.4 0.8]); 调整图像的灰度到0.4-0.8范围内；效果图如下图4所示：



图4
d)  结论：调整图像的灰度到0.4-0.6范围内时效果最好，能获得较好的灰度图像。



直方图均衡化
         J_histeq=histeq(I); 使用histeq函数对图像进行增强
         J_adapthisteq=adapthisteq(I); 使用adaptjisteq函数对图像进行增强,即自适应直方图均衡化

使用各函数效果对比图如下图5所示：



图5

结论：观察效果图可知，使用直方图均衡化能更清楚得看到破损图像的轮廓。


空域滤波增强

使用空域模板滤波进行的图像处理被称为图像的空域滤波增强，模板本身被称为空域滤波器。空域滤波增强的机理就是在待处理的图像中逐点移动模板，滤波器在改点的响应通过事先定义的滤波器系数与滤波模板扫过区域的相应像素值的关系来计算。空域滤波器分为平滑滤波器、中值滤波器、自适应除噪滤波器和锐化滤波器，但在实验过程中锐化滤波器并不能达到本实验的图像复原效果，故不采用。




（1）平滑滤波器滤波

平滑滤波器也称均值滤波器，它是用领域的均值代替像素值，减小了图像灰度的尖锐变化，
         J=filter2(fspecial('average',x),I)/255; %    x值可变，x不同值对应不同模板，默认值为3



3*3滤波模板：
      J=filter2(fspecial('average',x),I)/255;

效果图如下图6所示：

图6


5*5滤波模板：
J=filter2(fspecial('average',5),I)/255;

效果图如下图7所示：

图7


各个滤波模板间的比较如下图8所示：

图8


结论：由以上效果图可知，5*5滤波模板效果最好




（2）中值滤波器

中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术，其基本原理是把数字图像或数字序列中一点的值用该点的一个领域中各点值的中值代替，让周围的像素值接近真实值，从而消除孤立的噪声点。
        J=medfilt2(I);   %默认值时为3*3模板
        J=medfilt3(I,[m,n]);  % m*n 模板



3*3滤波模板
   J=medfilt2(I);

效果图如下图9所示：

图9


5*5滤波模板
    J=medfilt3(I,[5,5]);

效果图如下图10所示：

图10


各个滤波模板的比较效果图如下图11所示：


图11

结论：由以上图可以看出，随着模板的增大，噪声会越来越弱，但是图像会逐渐不清晰，而在5*5时效果是最好的。



（3）自适应滤波器

自适应滤波是指根据环境的改变，使用自适应算法来改变滤波器的参数和结构，使用wiener2函数用于对图像进行自适应除噪滤波，其可以估计每个像素的局部均值与方差。

J=wiener2(I,[M,N],noise)；%使用M×N大小领域局部图像均值与偏差，采用像素式自适应滤波器对图像I进行滤波。

各个滤波模板对比结果如下图12所示：

图12




频域滤波增强-低通滤波器

图像的频域增强是利用图像变换方法将原来的图像空间中的图像以某种形式转换到其他空间中，然后利用该空间的特有性质再进行图像处理，最后转换回原来的图像空间中，从而得到处理后的图像。频域滤波器中又分为低通滤波器和高通滤波器，二者之间的作用相反，又由于噪声主要集中在高频部分，所以本实验中采用低通滤波器才能达到一定的实验效果，故不考虑高通滤波器。低通滤波滤波除去了高频成分，而低频信息基本无损失地通过。



（1）理想低通滤波器

设傅里叶平面上理想低通滤波器离开原点的截止频率为，则理想低通滤波器的传递函数为：

$H(u,v)
\left\{
\begin{aligned} 
1     ,D(u,v)≤D_0  \\ 
\\ 
0     ,D(u,v)＞D_0
\end{aligned} 
\right.$
式中，$D(u,v)=(u^2+v^2 )^(1^/2^)$表示点(u,v)到原点的距离，D0表示截止频率点到原点的距离。
主要代码：
 I=imread('F:\inner课程作业\图像与视频处理\图片1.png');
 [f1,f2]=freqspace(size(I),'meshgrid');  %生成频率阵序列矩阵
 Hd=ones(size(I));
 r=sqrt(f1.^2+f2.^2);
 Hd(r>D_0)=0;        %构造滤波器, D_0可变
 Y=fft2(double(I));
 Y=fftshift(Y);
 Ya=Y.*Hd;                           %滤波
 Ya=ifftshift(Ya);
 Ia0X=ifft2(Ya);

不同D0值时对应效果图如下图13所示：

图13

结论：由上图可知，当r为0.1时效果最好


（2）	巴特沃斯低通滤波器

n阶巴特沃斯滤波器的传递函数：

$H(u,v)=1/(1+[D(u,v)/D_0 ]^2n )$

其特性时连续衰减，不像理想滤波器陡峭变化。

主要代码：
 I=imread('F:\inner课程作业\图像与视频处理\图片1.png');
 [f1,f2]=freqspace(size(I),'meshgrid');
 D=0.4;            %截止频率
 n=x;              %    x可变
 Hd=ones(size(I));
 r=sqrt(f1.^2+f2.^2);
 for i=1:size(I,1)
    for j=1:size(I,2)
      t=r(i,j)/(D*D);
      Hd(i,j)=1/(t^n+1);      %构造滤波函数
    end
 end
B=fft2(double(I));
B=fftshift(B);
Ba=B.*Hd;
Ba=ifftshift(Ba);
Ia1=ifft2(Ba);

不同n值时对应效果图如下图14所示：

图14

结论：由上图可以看出，n为10时效果最好。


（3）	指数低通滤波器

指数低通滤波器的转移函数为
$H(u,v)=exp^{-[D(u,v)/D_0  ]^n}$

主要代码：
I=imread('F:\inner课程作业\图像与视频处理\图片1.png');
[f1,f2]=freqspace(size(I),'meshgrid');
D=X/size(I,1);                         %  D为X，其中X可变
Hd=ones(size(I));
r=f1.^2+f2.^2;
for i=1:size(I,1)
   	for j=1:size(I,2)
 			    t=r(i,j)/(D*D);
  				Hd(i,j)=exp(-t);
   	end
end
E=fft2(double(I));
E=fftshift(E);
Ea=E.*Hd;
Ea=ifftshift(Ea);
IaX=ifft2(Ea);          %   X 可变

不同X值时对应效果图如下图15、16所示：

图15
图16


结论：当x 为100时效果最好

emmm第一次使用CSDN啦啦啦啦啦以后继续努力！！！！！!

图像复原技术的目的是使退化了的图像尽可能恢复到原来的真实面貌。
图像复原与图像增强的联系与区别：
图像复原和图像增强相似，两者都是要得到某种意义上改进的图像，或者说，希望改进输入图像的质量。两者不同之处是图像增强技术一般要借助人的视觉系统的特性，以取得看起来比较好的视觉效果。而图像复原则认为图像在某种情况下退化或者恶化啦，现在需要根据相应的退化模型和知识重建或者复原原始图像。尽管两者有相交叉的邻域，但图像增强主要是一个主观的过程，而图像复原的大部分过程是一个客观的过程。也就是说图像复原技术将图像退化的过程模型化，并据此采取相反的过程以得到原始的图像。
图像复原的前提是图像退化。图像退化是指图像在形成、记录、处理、传输过程中由于成像系统、记录设备、处理方法和传输介质的不完善，导致的图像质量下降。具体的说，常见的退化原因大致有成像系统的像差或有限孔径或存在衍射、成像系统的离焦、成像系统与景物的相对运动、底片感光特性曲线的非线性、显示器显示时失真、遥感成像中大气散射和大气扰动、遥感摄像机的运动和扫描速度不稳定、系统各个环节的噪声干扰、模拟图像数字化引入的误差等。
对退化图像的复原，一般采用两种方法。（1）适用于对图像缺乏已知信息的情况，此时可以对退化过程（模糊和噪声）建立模型，进行描述，并进而寻找一种去除或削弱其影响过程，由于这种方法试图估计图像被一些相对良性的退化过程影响以前的情况，故事一种估计方法
（2）若对于原始图像有足够的已知信息，则对原始图像建立一个数学模型并根据他对退化图像进行拟合会更有效。例如：已知图像中仅含有确定大小的圆形物体（如星辰、颗粒）这样由于仅是原始图像很少的几个参数（数目、位置、幅度等）未知，因此这是一种检测问题。
 
如：用约束的最小二乘方滤波复原模糊噪声图像
f=checkerboard(8);%生成标准棋盘测试图

PSF=fspecial('motion',7,45);%产生空间滤波器

gb=imfilter(f,PSF,'circular');

subplot(1,3,1);

imshow(f);

title('原图')

%产生噪声

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

g=gb+noise;

subplot(1,3,2);

%将原图像放大到512*512 ，并显示

%         imshow(pixeldup(f,8),[]);

%         subplot(1,3,2);

% %将噪声图像放大到512*512 ，并显示

% umshow(pixeldup(g,8),[]);

imshow(g);

title('加入噪声图像')

%利用函数deconvreg() 复原模糊图像

fr=deconvreg(g,PSF,0.4,[1e-7 1e7]);

subplot(1,3,3);

imshow(fr);

title('约束最小二乘方滤波复原图像');
结果如下：