匀速直线运动的模糊图像的盲复原的复现

前言

很久没写博客了，学习的效率很低，还是觉得要写博客作为学习的驱动力。。。
之前看了几篇图像复原论文，做个笔记。
主要参考：
1–南昌航空大学硕士论文–运动模糊图像复原研究
2–中国科学院大学硕士论文–图像的运动模糊复原算法研究

一 名称解释

运动模糊：当图像采集设备图像场景发生相对位移，不论是采集设备运动还是背景中目标在运动，都会使采集的图像发生降质，产生图像模糊。
前者是整个图像场景发生模糊，是空间不变型，只有一个模糊核，本文主要研究该类型。后者是目标发生移动，是空间可变型，可能有多个模糊核。
空间不变型：退化函数h(x,y)与x,y无关且具有可加性
盲复原：不知道PSF（点扩散函数）情况进行图像的复原

二 图像复原理论基础

1退化模型

原图像 f 在点扩散函数h的作用下，加上加性噪声n，最后形成了退化图像g
公式表示为：g(x, y) =f(x, y) ∗h(x, y) +n(x, y)

2加性噪声模型

不同的噪声模型有不同的概率密度函数，常用的噪声模型有高斯噪声，瑞丽噪声，伽马噪声，指数噪声，均匀噪声，椒盐噪声等。只处理该噪声的方法常常就是空间滤波。均值滤波，统计排序滤波，自适应滤波等。详情见冈萨雷斯第五章。

3傅里叶变换

冈萨雷斯第四章

4图像的退化模型推导

（1）
一般退化模型：
g(x, y) =f(x, y) ∗h(x, y) +n(x, y)
由傅里叶变换的卷积定理得，空间域得卷积等于频率域的乘积：
G(u,v)=F(u,v) H(u,v)+N(u,v)

（2）
运动模糊可以认为是图像在沿着模糊方向进行了灰度的平均。
假设图像在记录的过程中物体在水平面上发生匀速直线运动，在x和Y方向上的运动分别用时间函数Xo(t)和Yo(t)表示。曝光时间为T时.

带入，进行傅里叶变换

交换积分次序可以得到

对比一般退化模型，可以得到H的公式如下

（3）
将H转为空间域可以得到直线均匀点扩散函数的定义为

由公式可知，点扩散函数只和模糊角度和模糊长度有关。所以我们只要求到这两个参数就可以得到点扩散函数从而推出原图像。

三图像质量评价

1，图像质量主观评价
2，图像质量客观评价
客观评价又分为无参考图像质量评价和有参考图像质量评价。
区别在于有参考评价需要知道原图像。

1有参考评价

均方误差MSE
复原后图像与原图像每个像素点的差值的平方。
MSE越小代表评价越好

信 噪 比 参 数SNR和 峰 值 信 噪 比 参 数PSNR
同理

基于结构相似性的图像质量评价方法
结构相似度是从三个方面：局部图像亮度、局部结构、局部图像对比度，这三个方面用于得到一个总的相似性度量

2 无参考图像质量评价参数

灰度平均梯度法GMG
求梯度，梯度越大质量越好

拉普拉斯算子
求二阶差分，一样的道理

四点扩散函数离散化

五 运动模糊方向估计

频谱图中的条纹变化具有一定的规律性，除中心条纹外其它的条纹间距相等并且以中心条纹对称分布，对于这一特征，我们对一次傅立叶变换频谱图再次进行一次傅里叶变换得到的频谱图象并对其进行二值化处理，在运动模糊方向会出现一条亮条纹，当频谱图象中条纹倾斜角发生变化时，条纹也随之变化并且方向总是保持垂直，图象该亮条纹方向就是反映了运动模糊点扩展函数 h(x, y) 的方向性

通过雷登变换求方向
雷登变换公式：
p代表与直线与原点的距离，theta代表直线的角度。
由公式可知，只有当冲击函数的参数满足直线时候，冲击函数才为1，才进行积分，否则不积分。


由图可知：
当theta等于0时候，射线为90度射下，得到的结果广而强度低。
当theta等于45度时候，射线为135度射下，得到的结果集中且强度高。
可以知道，每列最大值中的最大值所在的方向为直线的方向。

六 运动模糊长度估计

D 为频谱中极小条纹间间隔，即投影方向的 radon 变换所得数组中的极大值旁边两个极小值间间隔的 0.5 倍.M，N 为图像的尺寸

七 图像复原模型

1逆滤波

H(u,v) 的幅度会随着u ，v 平面原点距离的增加而迅速减小，而噪声 N(u,v) 一般情况下变换比较缓慢，在远离u ，v 平面的原点时 N(u,v) / H(u,v) 的值就会变得很大即噪声放大。这时候恢复只能在与原点较近（接近频域中心）的范围内进行

2 维纳滤波

八 匀速直线运动的模糊图像的盲复原的复现

I=imread('cameraman.tif');
I=im2double(I);
psf = fspecial('motion',20,45);%45度移动20像素
blurred = imfilter(I, psf, 'circular', 'conv');%创建模糊图片
imshow(blurred);

F=fft2(blurred);
F=fftshift(F);
F=abs(F);
F=log(F+1);%傅里叶图像
F=mapminmax(F,0,1);
F = im2uint8(F);
F=medfilt2(F,[3,3]);
thresh = graythresh(F);     %自动确定二值化阈值
B = im2bw(F,thresh);       %对图像二值化
figure(),imshow(B);
theta= 0:179;
[R,xp]=radon(B,theta);
figure,imagesc(theta,xp,R); 
title('R_theta X'); 
xlabel('theta(degree)'); ylabel('X\prime'); 
colormap(hot);  colorbar;
col_max=max(R);%求每列最大值
[temp,re_theta]=max(col_max);%最大值和角度
[x,y]=size(R);
col_re_theta=R(:,re_theta);
figure();plot(1:1:x,col_re_theta);
[temp,total_max_x]=max(col_re_theta);
for i=total_max_x:x-1
	if col_re_theta(i)<col_re_theta(i+1)
		break;
	end
end
D=i-total_max_x;
[x,y]=size(I);
t=x/y;
L= x/D*sqrt(power(sin(re_theta),2)+power(t*cos(re_theta),2));
re_psf = fspecial('motion', L, re_theta);
wnr2 = deconvwnr(blurred , re_psf, 0.001);
figure, imshow(wnr2);

九 结论