图像复原和重建技术

1、常见的图像噪声模型

 图像复原是将图像退化的过程加以估计,并补偿退化过程造成的失真,以便获得未经干扰退化的原始图像或原始图像的最优估值，从而改善图像质量的一种方法。
 典型的图像复原方法是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用滤波等手段进行处理，使得复原后的图像符合一定的准则，达到改善图像质量的目的。

 图像复原
沿着质量降低的逆过程来重现真实的原始图像。通过去模糊函数而去除图像模糊。
 特点
图像增强：实际上是一种对比度拉伸，评价的标准是观察者主观感受；
图像复原：则是通过定量的去模糊函数而去除图像的模糊成分。

 图像噪声
图像获取过程：图像获取的环境条件和传感器质量。
图像传输过程：传输信道受到干扰。
 示例
成像系统镜头聚焦不准产生的散焦； 
相机与景物之间的相对运动； 
成像系统存在的各种非线性因素以及系统本身的性能；
射线辐射大气湍流等因素造成的照片畸变；
成像系统的像差、畸变、有限带宽等； 
底片感光图像显示时会造成记录显示失真； 
成像系统中存在的各种随机噪声 ；

 图像噪声（Probability density function简称PDF）

高斯噪声：

瑞利噪声：

伽马噪声：

指数分布噪声：

均匀分布噪声：

脉冲噪声（椒盐噪声）：

2、空间滤波图像复原

clc; clear all; close all;
% 加载图像
I = imread('cameraman.tif');
% 椒盐噪声
J = imnoise(I, 'salt & pepper');
J = double(J);
figure;
subplot(1, 2, 1); imshow(I, []); title('原始图像');
subplot(1, 2, 2); imshow(J, []); title('椒盐噪声的图像');
S1 = medfilt2(J);

figure;
subplot(1, 3, 1); imshow(I, []); title('原始图像');
subplot(1, 3, 2); imshow(J, []); title('椒盐噪声图像');
subplot(1, 3, 3); imshow(S1, []); title('中值滤波');

% 最大值滤波
S2 = ordfilt2(J, 25, true(5));
figure;
subplot(1, 3, 1); imshow(I, []); title('原始图像');
subplot(1, 3, 2); imshow(J, []); title('椒盐噪声图像');
subplot(1, 3, 3); imshow(S2, []); title('最大值滤波');

S3 = ordfilt2(J, 1, true(5));
figure;
subplot(1, 3, 1); imshow(I, []); title('原始图像');
subplot(1, 3, 2); imshow(J, []); title('椒盐噪声图像');
subplot(1, 3, 3); imshow(S3, []); title('最小值滤波');

S4 = ordfilt2(J, 5, true(3));
figure;
subplot(1, 3, 1); imshow(I, []); title('原始图像');
subplot(1, 3, 2); imshow(J, []); title('椒盐噪声图像');
subplot(1, 3, 3); imshow(S4, []); title('中值滤波');

3、频域滤波图像复原

clc; close all; clear all;
% 加载图像
I = imread('origin.jpg');
if ndims(I) == 3
    I = rgb2gray(I);
end
% 三次高斯带阻滤波器
J1 = GsFilter(I, 60, 30);
J2 = GsFilter(J1, 103, 30);
J3 = GsFilter(J2, 150, 20);
figure;
subplot(1,2,1); imshow(I, []); title('滤波前');
subplot(1,2,2); imshow(J3, []); title('滤波后');

4、逆滤波图像复原

5、维纳滤波（最小均方误差滤波）图像复原

●寻找到一个真实图像的估计值，使它们之间的均方误差最小。

clc; close all; clear all;
% 加载图像
I = im2double(imread('cameraman.tif'));
figure; imshow(I); title('原始图像');
% 运动模糊滤波器
LEN = 21;
THETA = 11;
PSF = fspecial('motion', LEN, THETA);
blurred = imfilter(I, PSF, 'conv', 'circular');
figure; imshow(blurred); title('退化图像');
% 高斯噪声
noise_mean = 0;
noise_var = 0.0001;
blurred_noisy = imnoise(blurred, 'gaussian', ...
    noise_mean, noise_var);
figure; imshow(blurred_noisy); title('退化+噪声图像');
% 逆滤波
estimated_nsr = 0;
wnr2 = deconvwnr(blurred_noisy, PSF, estimated_nsr);figure; imshow(wnr2); title('直接逆滤波图像');
estimated_nsr = noise_var / var(I(:));
wnr3 = deconvwnr(blurred_noisy, PSF, estimated_nsr);figure; imshow(wnr3); title('预估逆滤波图像');

6、约束的最小二乘滤波图像复原

在最小二乘复原处理中，常常需要附加某种约束条件。例如令Q为f的线性算子，那么最小二乘方复原的问题可以看成使形式为

最小化问题，这种有附加条件的极值问题可以用拉格朗日乘数法来处理。寻找一个  ，使下述准则函数为最小：                                                                                    

clc; clear all; close all;
I = checkerboard(8);
PSF = fspecial('gaussian',7,10);
V = .01;
BlurredNoisy = imnoise(imfilter(I,PSF),'gaussian',0,V);
NOISEPOWER = V*numel(I);
[J, LAGRA] = deconvreg(BlurredNoisy,PSF,NOISEPOWER);
subplot(221); imshow(BlurredNoisy);title('模糊及加噪图像');
subplot(222); imshow(J);title('直接处理');
subplot(223); imshow(deconvreg(BlurredNoisy,PSF,[],LAGRA/10));title('加入0.1*LAGRA');
subplot(224); imshow(deconvreg(BlurredNoisy,PSF,[],LAGRA*10));title('加入10*LAGRA');

7、Lucy-Richardson图像复原

       L-R算法是一种迭代非线性复原算法，它是从最大似然公式引出来的，图像用泊松分布加以模型化的。当下面这个迭代收敛时模型的最大似然
      函数就可以得到一个令人满意的方程：

clc; clear all; close all;
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('模糊及加噪图像');
subplot(222);imshow(J1);title('直接处理');
subplot(223);imshow(J2);title('加入NI,DP');
subplot(224);imshow(J3);title('加入NI,DP,WT');

8、案例演示：盲卷积图像复原

       最后，让我们演示基于盲卷积图像复原的编程实战。在图像复原过程中，最困难的问题之一是，如何获得PSF的恰当估计。那些
不以PSF为基础的图像复原方法统称为盲区卷积。
      它以MLE为基础的，即一种用被随机噪声所干扰的量进行估计的最优化策略。工具箱通过函数deconvblind来执行盲区卷积。

clc; clear all; close all;
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('模糊及加噪图像');
subplot(222);imshow(PSF,[]);title('实际的PSF');
subplot(223);imshow(J);title('复原图像');
subplot(224);imshow(P,[]);title('重构的PSF');