7. 图像复原（下）

7.1 图像退化与复原模型

公式表示：

7.2 估计退化函数

我们在处理图像复原问题时，往往都是假设退化函数已知，那么这个退化函数该如何去知道呢？这就需要事先估计退化函数，主要方法有以下三种：
1）图像观察估计法
设Gs(u，v)为观测到图像的子图像，Fs(u,v)为原始子图像的估计值，假设由于选择我们选择的是强信号区域，因此噪声可以忽略，有:

然后我们可以从Hs(u,v)推导出完整的函数H(u,v).
2）实验估计方法
假设我们已知图像的成像设备，拿这个设备去拍摄一幅已知频域的图像，比如一面白墙，经过过傅里叶变换，求得H(u,v):

3）模型估计法
Hufnagel等[1964]基于大气湍流的物理特征提出：

大气湍流模型的示例：

7.3 建模运动带来的图像模糊退化

1）现在，假设一幅图像由于相机的匀速直线运动变得成像模糊了，该过程可以建模为：

我们对它进行傅里叶变换，有：

特别地，若 ：

那么：

若：

那么：

2）利用matlab对运动产生的图像模糊进行建模
运动产生的图像模糊可以通过 matlab 图像处理工具相中的函数 fspecial 进行建模：

%len 表示移动几个像素 theta为移动的角度
PSF = fspecial(’motion’,len,theta)
g= imfilter(f,PSF,'circular');
g = g+noise;

测试图像可以通过函数 checkerboard 来生成，它的语法是：

%NP为每个正方形一边的像素数 M为棋盘行数，N为列数
C = checkerboard(NP,M,N)
f=checkerboard(8);
imshow(f)

利用matlab对运动产生的图像模糊进行建模举例:

f=checkerboard(8);
subplot(2,2,1);imshow(f);title('原图');
PSF=fspecial('motion',7,45);
gb=imfilter(f,PSF,'circular');
noise=imnoise(zeros(size(f)),'gaussian',0,0.001);
g=gb+noise;
subplot(2,2,2);imshow(gb,[]);title('添加运动模糊');
subplot(2,2,3);imshow(noise,[]);title('噪声图');
subplot(2,2,4);imshow(pixeldup(g,8),[]);title('添加噪声图');

7.4 逆滤波

图像退化过程：

对于一幅被退化函数 H 降质的图像，最简单的复原方法就是直接采用逆滤波，即：

加上噪声项：

如果退化函数在某一个频率点的值为零或很小的值，那么比值 可以很容易压制信号，造成 的估计存在较大误差。一种解决方法是只利用原点附近退化函数的值来进行逆滤波的计算。
逆滤波实验：

逆滤波实际情况很少用。

7.5 维纳滤波

1）维纳滤波器寻求能最小化下面的统计误差函数的估计值：

其中 E 为求取数学期望值的算子，f为没有退化的图像。这个表达式在频域中的解为：

其中H(u,v)为退化函数

通常情况下，我们使用维纳滤波的变型：

其中K为噪信功率比，我们通常取一个噪信功率的均值之比。如果K=0，就变成了逆滤波器。
在IPT中，维纳滤波是使用deconvwnr函数实现的，具体形式如下：

fr=deconvwnr(g,PSF)
fr=deconvwnr(g,PSF,NSPR)
fr=deconvwnr(g,PSF,NACORR,FACORR)

其中fr表示复原图像，g表示退化图像，PSF表示图像模糊，NSPR表示噪信功率比K，NACORR和FACORR表示噪声和未退化图像的自相关函数。

2）Matlab中的维纳滤波

f=checkerboard(8);
subplot(2,2,1);imshow(f);title('原图');
PSF=fspecial('motion',7,45);
gb=imfilter(f,PSF,'circular');
noise=imnoise(zeros(size(f)),'gaussian',0,0.001);
g=gb+noise;

fr1=deconvwnr(g,PSF);
subplot(2,2,2);imshow(fr1);title('K为0的维纳滤波');
Sn=abs(fft2(noise)).^2;
nA=sum(Sn(:))/prod(size(noise));
Sf=abs(fft2(f)).^2;
fA=sum(Sf(:))/prod(size(f));
R=nA/fA;
fr2=deconvwnr(g,PSF,R);
subplot(2,2,3);imshow(fr2);title('K为均值之比的维纳滤波');
NCORR=fftshift(real(ifft2(Sn)));
ICORR=fftshift(real(ifft2(Sf)));
fr3=deconvwnr(g,PSF,NCORR,ICORR);
subplot(2,2,4);imshow(fr3);title('使用相关系数维纳滤波');

通常情况下，噪声自相关系数和原图像自相关系数我们是不知道的，所以就要去调节K的大小，来获得较好复原。
如果复原的图像显示振铃现象，有时在调用 deconvwnr 之前使用 edgetaper函数会有所帮助：

J=edgetaper(I,PSF)

7.6 约束最小二乘法图像复原算法

1）对于本章讨论的所有方法来说，普遍都需要已知退化函数 的所有信息。
2）使用维纳滤波方法存在的困难:必须知道未退化图像和噪声的功率谱。
3）对功率谱比的恒定估计有时可以取得很好的结果，但并不意味着它总是一种很好的解决方法。
4）下面将要讨论的约束最小二乘滤波方法只需要知道噪声的均值和方差。

5）如果我们用矩阵来建模退化过程，有：

可以将卷积过程变为矩阵和向量乘积。
该方法的核心是H对噪声的敏感性问题。缓解这个问题的一种方法是将复原计算建立在一个图像平滑度的度量上。拉普拉斯算子(图像的二阶导数)似乎是一个不错的选择。

该问题在频域中的解为：

令：

并且：

可以证明关于r的单调递增函数 所以我们可通过调整 ,使得：

问题是：我们如何计算||n||2呢？

计算过程为已知||n||2，给定一个初始gamma，因为gama关于||r||2是递增的，通过迭代使得如下式子成立：

就可以解除一个gamma值，通过反傅里叶变换就可以恢复图像。

6）约束最小二乘在IPT使用deconvreg函数实现，语法为：

fr = deconvreg(g,PSF,NOISEPOWER,RANGE)

g是被污染的图像，fr是复原的图像，NOISEPOWER与||n||2成正比，也就是gamma的初始值，RANGE为gamma的取值范围。NOISEPOWER较好的估计值为 MN* [均值方 + 方差]。

7）示例：

f=checkerboard(8);
subplot(2,2,1);imshow(f);title('原图');
PSF=fspecial('motion',7,45);
gb=imfilter(f,PSF,'circular');
noise=imnoise(zeros(size(f)),'gaussian',0,0.001);
g=gb+noise;
fr = deconvreg(g,PSF,4);
subplot(1,2,1);imshow(fr,[]);title('默认gamma取值范围');
%调整gamma取值范围
fr = deconvreg(g,PSF,0.4,[1e-7, 1e7]);
subplot(1,2,2);imshow(fr);title('改变gamma取值范围并换初值');