2018年认证杯SPSSPRO杯数学建模
动态模糊图像复原
B题 动态模糊图像
原题再现：
人眼由于存在视觉暂留效应，所以看运动的物体时，看到的每一帧画面都包含了一段时间内 (大约 1/24 秒) 的运动过程，所以这帧画面事实上是模糊的。对电影的截图来说，动态画面的每一帧也都是模糊的，例如图 1 为某部电影截图，展现的是在高速飞行中的拍摄效果，所以俯拍到的路面字迹是模糊的。但是一般来说，电脑游戏的每一帧画面都是以清晰的静态方式绘制出的，所以需要较高的帧率才能感觉到平滑，否则感觉会不够流畅。为了以较低的帧率能够取得较流畅的感受，在计算机视觉技术中，人们开发出了能够模拟动态模糊效果的算法。
第二阶段问题： 在拍摄到的某一帧画面中，往往只有部分场景处于运动状态。例如图 2 是伦敦的公交车经过电话亭时拍摄到的画面，图中的主要景观中，只有车辆是运动的，所以产生了运动模糊的效果。请你设计有效的数学模型，将只有部分景观存在动态模糊的帧恢复出尽可能清晰的画面。在论文中可以使用图 2 ，也可以另行选用动态模糊的图像作为算例。
整体求解过程概述(摘要)
图像的模糊退化是数字图像处理和视觉系统中存在的一个重要问题，导致这种现象出现的原因有很多。根据产生模糊的原因不同，图像的模糊可以分为很多种，生活中常见的是运动模糊，它的实质大多是由于原目标图像与相机摄像头之间的相对运动造成了运动方向上的平滑滤波，即原始清晰图像与退化函数的卷积过程。在本题目中所给的图像，就是属于运动模糊，处理这种类型的图像，可以分为几个步骤。
首先，对图像进行预处理，即对图像进行格式转换、加窗处理。相机在视觉范围之内进行成像时，成像范围之外的边缘信息并未被卷积关系牵涉进来，因而在图像复原时由于缺少这些边缘信息而在图像边缘处产生截断现象，即产生边缘振铃效应，为了有效地抑制复原图像的振铃效应，提出了一种基于边缘恢复的抑制运动模糊图像振铃效应的方法——最优窗算法。其可在保留原图像边缘细节，且不增加运算量的前提下，有效抑制图像恢复过程中的振铃效应。然后对预处理后的图像进行傅里叶变换，和灰度均衡，得到图像的频谱，以便进行下一步的图像处理。
其次，为了得到点扩散函数的计算模型，需要对图像进行二值化和 Radon 变换。图像二值化就是将图像上的像素点的灰度值设置为 0 或 255，就是将整个图像呈现出明显的黑白效果的过程，从而得到图像的轮廓信息，图像的二值化有很多方法，有固定门限二值化、Dither 矩阵二值化、基本自适应门限二值化等，这些方法的运用需要选定一个合适的阈值，阈值可以手动选择也可以自动生成，按图片的需要进行选择。然后对图像进行 Radon 变换，在数学上，Radon 是一种积分变换，这个变换将二维平面函数变换成一个定义在二维空间上的一个线性函数,即相当于得到函数的投影讯号，对图像做 1 到180 的 Radon 变换，得到一个变换矩阵，找出这个矩阵中的最大值，此最大值所对应的列数即为图像的模糊角度。模糊的图像的频谱会出现黑色带条，通过分析频谱图，得到一个量化的反比关系公式和黑带的精确位置，通过微分自适应可以估算出运动模糊图像的模糊长度。
通过以上的分析，用 MATLAB 和 OPENCV 建立相应的数学模型，处理题目中和自行拍摄的图像，得出模拟结果，并对结果进行分析，指出缺点和不足。
问题分析：
观察所给图片，图中所示为公交车在路上行驶的画面，其中静止的行人和电话亭，以及商店的招牌都清晰可见，公交车因为运动而造成了在图像中的模糊。分析可知，图像的拍摄与公交车的平面有一定的角度，本题中的相机拍摄时是静止的，只有公交车是运动的，所以在求解时，需要对动态模糊的角度和动态模糊的长度的估算。在处理这种部分模糊的图像时，可以分为以下步骤：
1.首先需要对部分动态模糊图像进行分离，将图像中模糊部分和非模糊部分分离。
2.对分离后的模糊图像进行傅里叶变换得出其频谱，从而观察动态模糊的方向。
3.利用灰度均衡和消除噪声技术，对图像进行预处理，方便后面的处理过程。按照不同的方式分类，噪声分为白噪声、椒盐噪声、冲击噪声和量化噪声等，对噪声的消除有很多种方法，比如均值滤波，自适应维纳滤波，中值滤波等方法都可以很好的消除图片的噪声；
4.运用二值化和 Radon 变换求出图片的轮廓信息。
5.综合以上的处理过程，求出模糊图像处理中的两个重要参数：模糊角度和模糊长度。得到模糊图像的点扩散函数。
6.对模糊图像进行去模糊处理。其方法主要有逆滤波恢复、有约束最小二乘法恢复、维纳滤波恢复和 Richardson-Lucy 滤波恢复等方法。
在本文中我们基于 MATLAB 和 OPENCV 建立模糊图像的模型，利用改进后的维纳滤波方法处理模糊图像，并取得了较好的效果。因为带最优窗维纳滤波法可以有效的抑制振铃效应，获得高质量的恢复图像。
模型假设：
1、假设镜头已经准确的对焦
2、假设所有的噪声因素是白噪声
3、不考虑其他因素对图片造成的运动模糊
4、拍照的过程是一个瞬间动作。
论文缩略图：
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
部分程序代码：(代码和文档not free)
close all;
clear all;
%% 读入并显示图像
filename = '图 2.jpg';
I = imread(filename);
figure(1), imshow(uint8(I)); title('原图');
%% 对动态模糊图像进行灰度化，并进行二维快速傅里叶变换，生成其频谱图
gb = rgb2gray(I);
figure(4), imshow(uint8(gb));
[p,q]=size(gb);
gb=rot90(gb,2);
figure(5), imshow(uint8(gb));
PQ = paddedsize(size(gb));
F = fft2(gb, PQ(1), PQ(2));
figure(5), imshow(uint8(F));
%% 作出倒频谱
F1 = log(1+abs(F));
F2 = abs(F1).^2;
F3 = real(ifft2(F2));
figure(6), imshow(uint8(F3));
%% 将倒频谱压缩，居中
H = log(1+abs(F3));% 将倒频谱动态范围进行压缩
Hc = fftshift(H); % 将压缩结果进行循环移位，使低频成分居中
figure(9), imshow(uint8(Hc));
%% 通过阈值处理，边缘检测“canny”算子二值化倒频谱
T = graythresh(Hc);
bw=edge(Hc, 'canny', T);
figure(7), imshow(bw);
%% 对倒频谱从 1°到 180°作 radon 变换，以求出模糊角度
theta = 1:180;
R = radon(bw, theta);
figure(8), imshow(R);
%% 计算出通过倒频谱 radon 变换估计出的模糊角度
MAX = max(max(R));
[m, n] = find(R == MAX);
if 90 < n <= 180
beita = n - 90;
elseif 0 < n < 90
beita = n + 90;
elseif n == [90;90] | n == [180;180]
beita = n(1)
end
beita
len=5;
nsr=0.005;
psf=fspecial('motion',len, beita);
wnr1=deconvwnr(gb,psf,nsr);
figure(8),imshow(wnr1,[0 255]),title('维纳滤波复原图像(NSR=0.015)');
