基于Radon变换改进的运动模糊图像PSF参数估计算法

2020年
软 件
2020, V ol. 41, No. 6
基金项目: 国家国际科技合作专项项目(2014DFA00670)，贵州省研究生教育教学改革重点课题(黔教研合JG 字[2016]15)
作者简介: 陈健(1996–)，男，研究生，主要研究方向：数字图像处理；张欣(1976–)，男，副教授，主要研究方向：信号与信息处理、下一代无线通信及应用；陈忠仁(1992–)，男，研究生，主要研究方向：数字图像处理。

基于Radon 变换改进的运动模糊
图像PSF 参数估计算法
陈 健，张 欣，陈忠仁
（贵州大学 大数据与信息工程学院，贵州 贵阳 550025）
摘 要: 为了提高运动模糊图像的点扩展函数中模糊尺度，模糊角度估记的精准性，本文提出了一种改进的模
糊图像PSF 的参数估计算法。

首先将运动模糊图像进行3＊3分块，找出最能代表原始模糊图像模糊特征的子图像块，这样可以很好的估计运动参数。

然后将特征图像块的频谱图二值化后，进行形态学闭运算消除十字型亮线，有效的提高参数估计的精度。

其次，对频谱图进行腐蚀运算，用Radon 变换计算出模糊角度。

最后根据求出的角度将图像旋转到水平方向后，使用微分法求出模糊尺度。

并将本文算法与传统Radon 变换算法进行对比，表明该算法对提高模糊图像PSF 参数估计是有效的。

关键词: Radon 变换；运动模糊；图像分块；参数估计
中图分类号: TP391.41 文献标识码: A DOI ：10.3969/j.issn.1003-6970.2020.06.001
本文著录格式：陈健，张欣，陈忠仁. 基于Radon 变换改进的运动模糊图像PSF 参数估计算法[J]. 软件，2020，41（06）：01 06
Improved PSF Parameter Estimation Algorithm of Motion Blurred
Image Based on Radon Transform
CHEN Jian, ZHANG Xin, CHEN Zhong-ren
(College of big data and information engineering,Guizhou University, Guiyang 550025, China )
【Abstract 】: In order to improve the accuracy of blurred scale and blurred angle estimation in the point spread function of motion blurred image, an improved PSF parameter estimation algorithm is proposed. Firstly, the motion blurred image is divided into 3 * 3 blocks, and the sub image block which can best represent the blurred features of the original blurred image is found, so that the motion parameters can be estimated well. Then, after binarization of the spectrum of the feature image block, morphological close operation is carried out to eliminate the cross bright lines, effectively improving the accuracy of parameter estimation. Secondly, the spectrum was corroded and the am-biguity angle was calculated by Radon transform. Finally, according to the angle, the image is rotated to the hori-zontal direction, and the blurred scale is obtained by differential method. Compared with the traditional Radon transform algorithm, this algorithm is effective to improve the PSF parameter estimation of blurred image. 【Key words 】: Radon transform; Motion blurred; Image blocking; Parameter estimation
0 引言
人类从自然界接收的各种信息，超过80%是通过视觉获得的，图像作为一种关键的信息源，是人们感知世界、捕获信息、传递信息的重要手段。

图像在获取的过程中，许多因素会使得图像退化，具体表现为图像变得不清晰，根据退化原因主要有离焦模糊、高斯模糊以及运动模糊[1]三类。

其中，由
拍摄物体和成像设备发生相对运动引起的运动模糊在日常生活中最为常见，本文主要以运动模糊为研究对象。

在实际应用中我们需要的是清晰的图像，故通过一定的技术将图像变得清晰尤为重要。

通过一定方法计算出原始清晰图像退化时所形成的一些先验知识，包括点扩散函数[1]（point spread function ，PSF ），是实现运动模糊图像复原的关键一步。

PSF 模糊尺度L 和模糊角度θ这两个参数的估记是实
第41卷 第6期 软 件
2
《软件》杂志欢迎推荐投稿：**************.com
现图像复原的关键，模糊尺度L 和模糊角度θ能否准确估算出来直接影响图像复原的质量。

因此深入研究点扩散函数PSF 的参数估记有重大实际价值。

本文将在传统图像Radon 变换来估计PSF 参数的算法上做出改进，提高PSF 参数估计的精确性。

1 运动模糊图像退化模型和频谱图
在现实情况下，我们可以把由拍摄物体和成像
设备发生相对运动引起运动模糊的过程用一个线性模型来表示，该模型在时域上用公式表示如下：
(,)(,)*(,)(,)g x y f x y h x y n x y =+
（1） 在频域上用公式表示如下：
(,)(,)(,)(,)G u v F u v H u v N u v =+
（2）
根据造成图像模糊的运动方式，可将运动模糊图像划分为匀速运动造成的模糊、变速运动造成的模糊、无规则运动造成的模糊[2]等。

图像在积分时间内的模糊模型如下所示：
000
(,)((),())(,)T
g x y f x x t y y t dt n x y =--+⎰ （3）
其中，T 表示曝光时间，0()x t 和0()y t 分别表示图像在二维XOY 平面上两个方向上的运动距离，(,)n x y 表示加性噪声，若将噪声忽略，则图像在积分时间内的模糊模型变为如下所示：
000
(,)((),())T
g x y f x x t y y t dt =--⎰
（4）
该模糊模型在频域上的公式表示如下：
000
(,)((),())T
G u v f x x t y y t dtdxdy +∞+∞-∞-∞=--=⎰
⎰⎰
002(()())
0(,)T
j ux t vy t F u v e
dt π-+⎰
（5）
由公式（4）和公式（5）可得：
002(()())0
(,)T j ux t vy t H u v e dt π-+=⎰
（6）
为了我们研究方便，不妨将图像的运动看作匀速的，且运动轨迹是直线。

不妨设物体曝光时间为T ，被拍摄物体在水平方向和竖直方向的运动距离分别为a 和b ，则物体在二维XOY 平面上两个方向的运动分量分别为0()/x t at T =和0()/y t bt T =。

将这两个分量代入公式（6）中可以得到：
()(,)sin(())()j ua vb T
H u v ua vb e ua vb πππ-+=
++（7） 由于计算机处理的图像都是数字图像，不妨设图片的尺寸为N *M ，我们对公式（7）离散化后，便能得到公式8如下：
2(,)sin ua vb j M N T ua vb H u v e ua vb M N M N πππ⎛⎫
-+ ⎪⎝⎭⎡⎤⎛⎫=
+⎢⎥ ⎪⎛⎫⎝⎭⎣⎦+ ⎪⎝⎭
（8） 其中，(,)H u v 被称为运动模糊图像的PSF ，由公式（7）中可以看出(,)H u v 是一个正弦函数，存
在着周期性的零值，由于运动方向与退化图像的频带图中存在的暗色条纹的方向是互相垂直的，且这
些暗色条纹是彼此平行的，退化图像的频带图中的暗色条纹也具有周期性[3]。

下图1为原始Lena 图像与其灰度图像，图2、
3为运动退化图像与其频带图像：
图1 原始图像（左）和灰度图像（右）
Fig.1 Original image (left) and grayscale image (right)
图2 运动退化图像（左）和频带谱像（右） Fig.2 Motion degraded image (left) and
spectrum image (right)
图3 运动退化图像（左）和频谱图像（右） Fig.3 Motion degraded image (left) and
spectrum image (right)
以上图1实验采用原始真彩色Lena 图像的尺寸为512*512，环境为MATLAB R2015a 。

图2对Lena 灰度图像人为添加模糊角度为15度，模糊尺度为20PX 后使其退化。

图3是Cameraman 的退化图像和对应的频谱图像，其中Cameraman 原始灰度图像的尺寸为256*256，模人为添加糊角度为45度，模糊尺度为10PX 后使其退化。

通过频谱图像我们可以发现频谱图中暗条纹与水平方向的夹角的角度与真实模糊角度的差接近90度，也即是说暗条纹的方向与模糊角度是可以近似当作是垂直的。

在图2的
陈 健等：基于Radon 变换改进的运动模糊图像PSF 参数估计算法
3
《软件》杂志欢迎推荐投稿：**************.com
频带图中我们可以明显的看出，频带图中存在许多距离相等的暗条纹，而这些暗条纹分别对应着PSF 函数在频域中的零值。

观察图3我们不难发现其频带图中还呈现了十字形的亮线。

2 传统Radon 变换参数估计算法理论
二维XOY 平面上的Radon 变换是把面xy 上的某条直线ρ在Radon 内映射为点(,)ρθ，也即是在一个平面内沿不同的直线对(,)f x y 做线积分，得到的像函数就是函数(,)f x y 的Radon 变换。

ρ及连续图像的Radon 变换表达式分别如下：
cos sin x y ρθθ=+
（9）
(,)(,)(cos sin )D
R f x y x y dxdy ρθδρθθ=--⎰⎰（10） 其中，ρ表示原点到退化后图像中的任意一点
(,)x y 的间距，θ是水平方向与原点和直线垂直线的
夹角，D 代表图像所在的平面，(,)f x y 是点(,)x y 对应像素点的一个灰度值，δ代表单位冲击函数。

将公式（10）变形后我们可以得到下式：
()(cos sin ,sin s R x f x y x y co dy ξξξξξξ
θθθθθ)+∞
-∞
=--⎰
（11）
在公式（10）中，x ξ和y ξ的表示如下矩阵所示：
cos sin sin cos x x y y ξξθθθθ⎛⎫⎛⎫⎛⎫
= ⎪ ⎪⎪ ⎪-⎝⎭⎝⎭⎝⎭
（12） 在求得退化图像的频谱图后，我们采取Radon 变换来对运动退化图像的模糊角度进行估量，其估量原理为通过对图像采取Radon 变换来检测频带图中的直线并求出其角度[4]，由于运动退化方向与暗色条纹是互相垂直的，当求出暗色条纹的角度时，模糊角度的值自然能够求出。

另外一个参数模糊尺度指的是在很短的曝光时间内，物体所运动距离的大小，单位为像素（PX ）。

在得到运动模糊图像的频谱图的基础上，通过分析能够知中央位置暗条纹间的跨距是其他暗条纹跨距的两倍[5]，且暗色条纹的数目和跨距直接反应出模糊尺度的像素值大小。

当模糊尺度的像素值逐渐增加时，暗色条纹之间的跨距会逐渐减小，与此同时数目也逐渐增多。

设暗色条纹的跨距大小为d ，即
为点(,)(0,0)u v =到直线1ua vb
M N
+=的距离，其数学
表达式如下所示：
d =
（13） 假设中间位置条纹的跨距为D ，模糊尺度像素
值为L ，则有2D d =，cos a L θ=，
sin b L θ=，为方便计算，设在退化图像中
()0M N λλ=>。

将公式（13）做变换如下所示：
D =
（14）
也即是模糊尺度L 的计算公式为：
L =
（15）
当1λ=时，有M N =，此时退化图像的宽以及
高是相等的，则有2/L M D =，此时模糊尺度L 与模糊图像频带图中暗条纹的跨距D 成反比例关系，模糊尺度L 的大小便可计算出来。

综上所诉，我们可以把传统运动模糊图像PSF 参数估计算法用下面图4流程图表示。

图4 传统运动模糊图像PSF 参数估计算法流程图 Fig.4 Flow chart of traditional motion blur image PSF
parameter estimation algorithm
3 改进的运动模糊图像的PSF 参数估计算法
上述算法虽然能实现对PSF 参数的估计，但是发现该算法遇到某类图像时，估计出来的参数值误差不小，这类图像有一个共同的特点就是图像中存在一大块区域是色调很单一，在前文所述的图像退化模型中，我们考虑的是像素点扩散方式是相同的。

在色调单一区域内的像素点受到PSF 的影响后，因为领域像素点的灰度值差异小，受到模糊核[6]卷积运算后造成的像素叠加从视觉上看并不明显，色调单一区域会对真实值产生较大影响，故需要将色调单一的区域去掉后再进行PSF 参数估计。

基于此我们先对图像进行3*3的分块[7]处理，找出最能代表原始模糊图像模糊特征的部分区域图像，对具有替代性的子区域进行PSF 参数估量。

现对MATLAB R2015a 自带的标准测试图像Cameraman 做3*3分块处理，测试图像尺寸为256*256，通过MATLAB 软件分别对Cameraman 原图像、运动退化图像、运动退化图像的边缘检测效果如下图5所示。

第41卷 第6期 软 件
4
《软件》杂志欢迎推荐投稿：**************.com
图5 Cameraman 灰度图像（左）、运动模糊
图像（中）和边缘图像（右）
Fig.5 Cameraman grayscale image (left)、motion degraded image (middle) and edge image (right)
其中，对Cameraman 原始图像手动添加的模糊尺度为15PX ，模糊角度为15度后使其退化。

本文使用Canny 算子对模糊图像进行边缘检测，主要是因为Canny 算子较其它几个算子不仅可以充分的识别出模糊图像的边沿信息，而且还能抑制噪声干扰[8]。

同时，这会使后续检测更加精确。

对Cameraman 原始灰度图像做3*3分块处理效果如下图6所示。

图6 Cameraman 3*3分块图像 Fig.6 Cameraman 3*3 block image
图6分块图像边缘检测结果如下图7所示。

图7 3*3分块图像边缘检测图 Fig.7 3*3 Block image edge detection
从图7中找出边缘信息最多的4个图像块，如下图8所示。

图8 边缘信息最多的4个图像快
Fig.8 Four images with the most edge information
图8所对应的特征图像块如下图9所示。

图9 特征图像块
Fig.9 Feature block image
我们将模糊图像分块处理后再对其边缘检测，找到对应边缘信息最多几个区域，原因是在边缘处的像素值灰度变换会比其它部分的像素值的灰度变更剧烈，然后像素叠加后在视觉上的模糊效果更明显，对这些区域进行PSF 参数估计后的结果会更精准。

为了减小某些运动退化图像频带中存在的十字亮线对模糊尺度估计的千扰，我们将模糊图像的频带图后做一次二值化[9]后的闭运算操作，最后对图像不同程度腐蚀。

针对不同的模糊图像中心条纹的宽度不同，对进行闭运算后的图像先设置阈值[9]进行形态学腐蚀操作。

先测中心条纹的宽与高，得到宽和高的值中小的那个值，设置阈值Ｔ，如果该值大于Ｔ则继续腐蚀，反之则对图像进行Radon 变换求得模糊角度。

运动模糊图像频谱图的二值化、闭运算和腐蚀后如图10所示。

图10 二值化频谱（左）、闭运算后图像（中）
和腐蚀后的图像（右）
Fig.10 Binary spectrum (left), image after closed operation (middle), and image after corrosion (right)
以上实验原图为Lena ，尺寸512*512，人为设置的模糊角度为30度，模糊尺度为10PX 。

陈 健等：基于Radon 变换改进的运动模糊图像PSF 参数估计算法
5
《软件》杂志欢迎推荐投稿：**************.com
由于利用 Radon 变换直接求取模糊尺度[10]存在一定的缺点，当模糊图像中含有的噪声强度较大时，其频谱图的暗条纹会变得模糊不清。

不利于估量这些暗条纹之间的间隔，且当模糊尺度比较小时，频带图经过阈值化后会呈现出长椭圆形状，导致在该方向上的 Radon 变换得到的投影图发生振荡，当模糊尺度比较大时，估算出来的结果产生的误差值较大。

故使用微分法[11]来提高估计精度。

当我们用Radon 变换把模糊角度计算出来之后，把图像的运动方向回转到水平方向，然后将图像在竖直方向上求微分运算，微分表达式如下：
(,)(,1)(,)y g i j g i j g i j =+- （16） 将旋转到水平方向的图像做一阶微分运算，以达到去除图像相邻像素之间的相关性，而且增强图像的高频信息，抑制低频信息。

一阶微分表达式如下：
(,)(1,)(,)g i j g i j g i j =+- （17）
微分图像中的每一行的自相关量如下所示：
()()()i i i P n g m n g m +∞
-∞=+∑
（18）
其中，()i P n 为微分图像(,)g i j 中第i 行的自相关向量，()i g m 表示图像中第i 行m 列的像素值。

整幅图像的微分自相关函数如下所示：
1i P P M
= （19）
通过计算曲线中2个负共轭相关峰[12]的间距我们可以得到模糊尺度。

综上，本文使用的PSF 参数估计算法的流程图如下图11所示。

图11 改进的运动模糊图像PSF 参数估计算法流程图 Fig.11 Flow chart of improved motion blur image PSF
parameter estimation algorithm
4 实验仿真及结果分析
4.1 仿真环境
本实验利用MATLAB R2015a 来编程实现，使用的图像数据为MATLAB R2015a 自带的4张不一样的测试图像，其详细参数如下表1所示。

表1 实验图像参数表
Tab.1 Experimental image
图像名称 图像尺寸图像格式 灰度图像/RGB 图像
Cameraman 256*256.tif 灰度图像 Lena
512*512
.bmp
RGB 图像
其中，Lena 图像为RGB 图像，在实验时我们需要将它转换成灰度图像。

4.2 实验过程与数据处理
使用MATLAB R2015a 软件读取上述四幅图像，按照图4传统模糊图像PSF 参数估计算法求取模糊图像的模糊角度θ和模糊尺度L 。

按照图11改进的运动模糊图像PSF 参数估计算法（本文算法）求取模糊图像的模糊角度θ和模糊尺度L ，当我们对上述四幅图像进行3*3分块处理后，统一选取4个小特征图像块，对该4个图像快进行模糊角度θ估计得到四组数据，去掉四组数据中的最大值与最小值，将剩下的两组数据取平均值后作为最终的θ参数估量值。

4.3 实验结果与分析
首先选取Cameraman 、Lena 两幅无噪测试图像，对两幅无噪图像手动添加模糊参数，再用本文PSF 参数估计算法与传统PSF 参数估计算法分别对人为模糊后的图像进行PSF 参数的估量，其中，对Cameraman 无噪图像模糊角度自15度起，以15度为间隔取6组角度值，模糊尺度自5PX 起，每个5PX 取6组像素值。

对Lena 无噪图像模糊角度自20度起，以20度为间隔取6组角度值，模糊尺度自15PX 起，每个5PX 取6组像素值。

得出的数据如下表2、表3所示。

表2 Cameraman 图像模糊参数估计结果 Tab.2 The result of Cameraman image
blurred parameter estimation
真实模糊角度 传统算法估计角度本文算法 估计角度 真实模糊 尺度 传统算法估计尺度本文算法估计尺度15°
13° 15° 5PX 6PX 6PX 30° 28° 30° 10PX 12PX 10PX 45° 44° 44° 15PX 16PX 16PX 60° 62° 61° 20PX 22PX 21PX 75° 78° 76° 25PX 28PX 26PX 90°
93°
92°
30PX
27PX
29PX
现在我们将Cameraman 、Lena 两幅测试图像手动添加均值为0，方差为0.005的高斯白噪声，使其 成为含噪图像，再用本文PSF 参数估计算法与传统PSF 参数估计算法分别对人为模糊后的图像进行模糊参数的估量，其中，对Cameraman 无噪图像的模
第41卷 第6期 软 件
6
《软件》杂志欢迎推荐投稿：**************.com
糊角度自15度起，以15度为间隔取6组值，模糊尺度自5PX 起，每个5PX 取6组像素值。

对Lena 无噪图像模糊角度自20度起，以20度为间隔取6组角度值，模糊尺度自15PX 起，每个5PX 取6组像素值。

得出的数据如下表4、表5所示。

表3 Lena 图像模糊参数估计结果 Tab.3 The result of Lena image blurred
parameter estimation
真实模糊 角度 传统算法 估计角度 本文算法 估计角度 真实模糊 尺度 传统算法 估计尺度 本文算法估计尺度20° 19° 20° 15PX 13PX 15PX 40° 38° 40° 20PX 18PX 19PX 60° 58° 60° 25PX 24PX 25PX 80° 82° 81° 30PX 28PX 30PX 100° 103° 98° 35PX 32PX 34PX 120°
123°
120°
40PX
38PX
39PX
表4 Cameraman 含噪图像模糊参数估计结果 Tab.4 The result of Cameraman image with noise
blurred parameter estimation
真实模糊角度 传统算法估计角度 本文算法估计角度 真实模糊尺度 传统算法估计尺度 本文算法估计尺度15° 12° 15° 5PX 6PX 5PX 30° 28° 30° 10PX 12PX 10PX 45° 42° 44° 15PX 13PX 17PX 60° 63° 61° 20PX 16PX 19PX 75° 79° 76° 25PX 26PX 31PX 90°
94°
92°
30PX
32PX
34PX
表5 Lena 含噪图像模糊参数估计结果 Tab.5 The result of Lena image with noise
blurred parameter estimation
真实模糊角度 传统算法估计角度 本文算法估计角度 真实模糊尺度 传统算法估计尺度 本文算法估计尺度20° 21° 20° 15PX 13PX 14PX 40° 37° 41° 20PX 17PX 19PX 60° 57° 59° 25PX 22PX 23PX 80° 83° 81° 30PX 33PX 31PX 100° 104° 98° 35PX 31PX 32PX 120°
123°
122°
40PX
43PX
42PX
通过这两组实验数据可以看出，无论是有噪声图像还是无噪声图像，总的来说本文算法在求模糊图像的PSF 参数中的模糊角度与模糊尺度的问题上都要比传统PSF 参数求解方法更加精确。

对比表2与表4可以知道，当Cameraman 图像的模糊尺度为
25PX 、30PX 时并且含噪声时，本文算法估计效果较差。

对比表3与表5可以看出，本文算法在估计运动模糊图像PSF 参数时都更加精确些。

5 结论
本文在在原始模糊图像的PSF 参数的求解中，结合了图像分块理论和一些形态学图像处理操作，并在MATLAB R2015a 仿真软件上进行了实验，对实验结果进行分析表明该算法能提高运动模糊图像PSF 参数估计的准确率。

在下一步的学习过程中，我们可以充分考虑不同强度的高斯噪声以及其他类型的噪声对PSF 参数估计的影响，并且继续研究运动模糊图像在某些特定模糊参数下出现的PSF 参数估计精准性较差的问题。

参考文献
[1]
王玉全, 隋宗宾. 运动模糊图像复原算法综述[J]. 微型机与应用, 2014, 33(19): 54-57. [2]
杨丹, 赵海滨, 龙哲, 等.MATLAB 图像处理实例详解[M].北京: 清华大学出版社, 2013. [3]
盛文正.运动模糊图像复原技术的研究[D]. 长沙: 长沙理工大学电气与信息工程学院, 2015. [4]
吴海洋, 周仲礼, 曾凡霞, 等. 盲小波算法在遥感图像去噪中的应用[J]. 软件, 2012, 33(1): 9-11. [5]
J. Jeon, I. Yoon, and others, "Fully Digital Auto-Focusing System with Automatic Focusing Region Selection and PSF Estimation," IEEE Trans.on Consumer Electronics, Aug. 2015, PP. 1204-10.
[6] 张书豪. 单幅运动模糊图像复原技术的研究与实现[D].湘潭: 湘潭大学物理与光电工程学院, 2017.
[7]
陈俊杰, 鲁康. 基于直方图平移的图像二次分块可逆数据隐藏方法[J]. 软件, 2015, 36(7): 114-117.
[8] 曹妍, 陈伟, 徐森. 图像去噪方法研究与仿真[J]. 软件, 2015, 36(4): 33-36.
[9] 方奋奇, 边学军. 一种基于二值图像的信息隐藏技术. [J]. 软件, 2012,33(1):27-28
[10] Y T. Hwang, H R. Wen and B C. Tsai, “Efficient Block
Adaptive Point Spread Function Estimation and Out-of-Focus Image Restoration Scheme,” in Proc. of IEEE ICCE-TW, 2015.
[11] 廖永忠, 蔡自兴, 何湘华. 基于点扩散函数参数辨识的运
动模糊图像的盲恢复研究[J]. 光学技术, 2014, 40(3): 235-239.
[12] 郭璇, 郑菲, 赵若晗, 等. 基于阈值的医学图像分割技术
的计算机模拟及应用[J]. 软件, 2018, 39(3): 12-15.。