CN104067312B

CN104067312B - 对噪声有鲁棒性的图像配准方法和***

Info

Publication number: CN104067312B
Application number: CN201280064221.4A
Authority: CN
Inventors: G·陈; S·库隆布
Original assignee: ECOLE TECHNOLOGIE SUPERIEURE
Current assignee: ECOLE TECHNOLOGIE SUPERIEURE
Priority date: 2011-12-24
Filing date: 2012-12-24
Publication date: 2017-10-24
Anticipated expiration: 2032-12-24
Also published as: IL232189A0; CA2850854A1; US8942512B2; KR20140109439A; EP2795579A1; IL232189A; JP2015503778A; US8958659B2; EP2795579B1; JP6147273B2; US20130272627A1; US20130163896A1; WO2013091088A1; CN104067312A; EP2795579A4

Abstract

公开了一种用于将失真图像处理成参考原始图像对准的配准图像的图像配准方法。对原始图像的失真可以包括缩放、旋转和噪声。所述方法基于关联两个图像的拉东变换来确定旋转角，以及在任何添加的噪声将会抵消的前提下通过除以每个图像的整体亮度的平均值来确定缩放因数。使用快速傅里叶变换(FFT)估算全局空间移位。在一个实施方式中，失真图像在旋转之前首先被缩放成原始图像的大小。在另一个实施方式中，在旋转失真图像之前原始图像首先被缩放成失真图像的大小，最后将失真图像进行缩放以匹配原始图像。还提供了用于图像配准的相应***。

Description

对噪声有鲁棒性的图像配准方法和***

技术领域

本发明总体上涉及数字图像处理，特别地，涉及结合参考图像来处理图像。

本发明涉及图像配准，包括估算参考图像与失真图像之间的平移、旋转和缩放参数。特别地，本发明涉及对噪声有鲁棒性的图像配准 (registration)方法和***。

本发明涉及评估图像的视觉质量，以及特别涉及对仿射变换具有鲁棒性的评估图像的视觉质量的方法和***。

背景技术

图像配准是诸如图像质量评估、医学成像、自动目标识别等应用中的重要任务。现有的图像配准方法包括例如：A.Averbuch、R.R.Coifman、 D.L.Donoho、M.Israeli、Y.Shkolnisky和I.Sedelnikov在SIAM期刊 Scientific Computing2008年1月第30卷第二本785-803页的论文“A Framework for Discrete Integral Transformations II—The2DDiscrete Radon Transform”，他们使用了对图像的大的平移、旋转和缩放的基于伪极坐标的估算；G.Wolberg和S.Zokai在图像处理IEEE国际会议论文集(ICIP)的2000年493-496页的论文“Robust Image Registration Using Log-Polar Transform”，他们还通过使用对数极坐标变换研究了鲁棒图像配准；G.Varghese和Z.Wang在视频技术电路和***IEEE学报2010 年第20卷第七本1032-1040页的论文“Video denoising based on aspatiotemporal Gaussian scale mixture model”，他们使用了傅里叶变换来估算全局空间转移；B.S.Reddy and B.N.Chatterji在图像处理IEEE学报1996年第5卷第八本1266-1271页的论文“An FFT-based technique for translation,rotation and scale-invariant image registration”，他们提出了基于FTT的平移、旋转和缩放不变的图像配准技术；K.Jafari-Khouzani 和H.Soltanian-Zadeh在模式分析与机器智能IEEE学报2005年第27卷第六本1004-1008页的论文“Radon transform orientation estimationfor rotation invariant texture analysis”，他们将拉东变换运用到旋转不变纹理分析的方向估算上；E.De Castro和C.Morandi在模式分析与机器智能 IEEE学报1987年第PAMI-95卷700-703页的论文“Registration of translated and rotated images usingfinite Fourier transforms”，他们使用有限傅里叶变换开发了平移和旋转图像的图像配准方法，但是该方法不能够呈现缩放变化；W.Wei、S.Wang,X.Zhang和Z.Tang在信息取证与安全IEEE学报2010年第5卷第三本507-517页的论文“Estimation of image rotationangle using interpolation-related spectral signatures with application toblind detection of image forgery”，他们通过使用插值相关的光谱特征来应用到盲检测伪造图像估算了图像旋转角；D.G.Lowe在国际期刊Computer Vision2004年第60卷第2本91-110页的论文“Distinctive Image Features from Scale-Invariant Keypoints”，他提出了缩放不变关键点提取和配准(SIFT)的独特的图像特征方法；以及H.Bay、A.Ess、T.Tuytelaars和L.Van Gool在计算视觉和图像理解(CVIU)2008年第110 卷第三本346-359页的论文“SURF:Speeded Up Robust Features”，他们开发了SURF(快速鲁棒特征)，其自称比SIFT更快。

在全参考图像质量评估中，通常将失真图像与参考图像比较来获得视觉质量度量，通过比较两个图像来计算该视觉质量度量。存在多个用于生成这种质量度量的方法，包括：Z.Wang、A.C.Bovik、H.R.Sheikh和 E.P.Simoncelli在图像处理IEEE学报2004年第13卷第四本600-612页的论文“Image quality assessment:From error visibility tostructural similarity”中提出的简单的峰值信噪比(PSNR)测量的结构相似度(SSIM) 指数、H.R.Sheikh和A.C.Bovik在图像处理IEEE学报2006年第15卷第二本430-444页的论文“Image information and visual quality”中提出的视觉信息保真(VIF)指数。而且，Rezazadeh和Coulombe在信号、图像和视频处理2011年9月1-15页的论文“Novel discretewavelet transform framework for full reference image quality assessment”中和在Procedia 工程(Procedia Engineering)(爱思唯尔期刊)中发表的“A novel discretewavelet domain error-based image quality metric with enhanced perceptualperformance”˙中提出了全参考图像质量评价的离散小波变换框架，而Qian和Chen在奥地利维也纳的摄影测量、遥感和空间信息科学ISPRS国际期刊2010年7月5-7日第204-208页的论文“Four reduced-reference metrics for measuring hyperspectral images afterspatial resolution enhancement”中研发了用于测量空间分辨率增强后高光谱图像的视觉质量的4种减少的参考度量。

然而，还需要在工业中研发用参考图像调整图像的简单而快速并且准确的方法，该方法可以避免和减轻现有技术的缺点，能够改善对视觉图像质量的客观测量，以及可以提供对噪声有鲁棒性的图像配准。

发明内容

本发明的目的是提供一种用于将第一图像Y处理成与第二图像X对准的配准图像Y*的图像配准方法和***。

根据本发明的一个方面，提供了一种将第一图像Y处理成与第二图像 X对准的配准图像Y*的图像配准方法，所述方法包括：

(a)确定参考图像X与失真图像Y之间的缩放因数“a”；

(b)通过缩放因数“a”改变失真图像Y的大小，从而生成规范化的失真图像Y2；

(c)确定参考图像X和规范化的失真图像Y2之间的旋转角“θ₀”；

(d)通过旋转角“θ₀”旋转规范化的失真图像Y2，从而获得旋转的失真图像Y3；以及

(e)通过将旋转的失真图像Y3平移成与参考图像X水平和垂直对准来获得配准的图像Y*。

在上述方法中，步骤(a)包括通过计算参考图像X的像素值的和与失真图像Y的像素值的和之比来确定缩放因数“a”。

步骤(a)还包括用零值像素填充参考图像X和失真图像Y中较小者，以使得参考图像X和失真图像Y的水平和垂直尺寸相等(equalise)。

步骤(b)还包括使失真图像Y居中从而形成居中的失真图像Y1，以及通过缩放因数“a”改变居中的失真图像Y1的大小来生成规范化的失真图像Y2。

步骤(c)包括分别形成参考图像X和规范化的失真图像Y2的拉东变换R1和R2，以及计算拉东变换R1和R2的行的互相关来确定旋转角“θ₀”。

计算互相关包括：

(i)计算拉东变换R1的每行与拉东变换R2的每行之间的循环互相关集，每个循环互相关定义各行之间的旋转偏移角“θ”；

(ii)为每行选择具有最高值的循环互相关；以及

(iii)为每行确定每个选择的循环互相关定义的旋转偏移“θ”，以

及将旋转角“θ₀”设置为等于确定的旋转偏移“θ”的中值。

计算互相关还包括：

(iv)计算拉东变换R1和拉东变换R2的一维(1D)快速傅里叶变换 (FFT)；

(v)将拉东变换R1和拉东变换R2的1D FFT合并成中间结果R3；

(vi)生成包括循环互相关的第二中间结果R4，第二中间结果R4基于中间结果R3的元素；

(vii)在中间结果R3的每行中选择具有最高值的循环互相关；以及

(viii)确定每个选择的循环互相关中定义的对应的旋转偏移“θ”，以及将旋转角“θ₀”设置成等于确定的旋转偏移“θ”的中值。

合并步骤包括形成拉东变换R1的1D FFT与拉东变换R2的1D FFT 的“共轭”函数的按分量逐个作出的乘积(component-wise product)。

在上述方法中，步骤(vi)还包括：对应于各个角“θ”为每个分量对第二中间结果R4中的行确定最大值。该方法还包括取最大值的中值，所述中值对应于旋转角“θ₀”。

该方法还包括对预处理图像Y*执行视觉质量评估，例如执行以下视觉质量评估中的任一者：确定峰值信噪比(PSNR)；计算结构相似性(SSIM) 指数；以及计算视觉信息保真(VIF)指数。

根据本发明的另一方面，提供了一种将第一图像处理成与第二图像对准的配准图像Y*的图像配准方法，所述方法包括：

(a)确定第一图像与第二图像之间的缩放因数“a”；

(b)通过用缩放因数“a”改变第一图像的大小来使得第一图像和第二图像的大小相等，从而生成大小改变版本的第一图像；

(c)使用拉东变换来确定大小改变版本的第一图像与第二图像之间的旋转角“θ₀”；

(d)通过旋转角“θ₀”旋转大小改变版本的第一图像，并在将其平移成与第二图像水平和垂直对准后将其输出为配准的图像 Y*。

根据本发明的再一个方面，提供了一种将失真图像Y处理成与参考图像X对准的配准图像Y*的***，该***包括：

处理器；

计算机可读存储介质，具有存储于其上的用于由处理器执行的计算机可读指令，形成：

缩放因数估算模块，用于确定缩放因数“a”；

图像大小改变模块，用于使用缩放因数“a”将失真图像Y的大小改变成大小改变的失真图像Y2；

旋转角确定单元，用于估算旋转角“θ₀”；

图像旋转模块，用于将大小改变的失真图像Y2旋转旋转角“-θ₀”形成修正的图像Y3；以及

图像平移模块，用于平移修正的图像Y3，以生成配准图像Y*。

在上述***中，所述旋转角确定单元还包括：

拉东变换模块，用于分别生成大小改变的失真图像Y2和参考图像X 的拉东变换R1和R2；以及

拉东变换修正模块，用于根据拉东变换R2和R1估算旋转角“θ₀”。

该***还包括图像估算模块处理模块，用于对配准图像Y*进行视觉质量评估。

在上述***中，图像估算模块处理模块通过以下任一种方式执行视觉质量评估：

确定峰值信噪比(PSNR)；

计算结构相似度(SSIM)指数；以及

计算视觉信息保真(VIF)指数。

该***还包括：

图像居中模块，用于将失真图像Y居中，以及生成居中图像Y1；

其中所述图像大小改变模块被配置为使用缩放因数“a”对居中图像Y1改变大小成大小改变的失真图像Y2。

该***还包括：

图像遮蔽模块，用于从失真图像Y中提取遮蔽图像Y0；以及

质心(centroid)计算模块，用于根据遮蔽图像Y0确定中心参数“c_r”和“c_c”；

其中图像居中模块被配置为使用中心参数“c_r”和“c_c”根据失真图像 Y生成居中图像Y1。

该***还包括：

平移估算模块，用于确定居中图像Y3与参考图像X之间的偏移向量“TV”；

其中所述图像平移模块被配置为使用偏移向量“TV”平移修正图像 Y3以生成配准图像Y*。

下面还描述了本发明另一个实施方式的***和方法。

根据本发明的另一个方面，提供了一种将失真Y处理成与参考图像X 对准的配准图像Y*的图像配准方法，该方法包括：

(a)确定参考图像X与失真图像Y之间的缩放因数“a”；

(b)通过缩放因数“a”的倒数(inverse)改变参考图像X的大小，从而生成大小改变的参考图像X2；

(c)确定大小改变的参考图像X2和失真图像Y之间的旋转角“θ₀”；

(d)将失真图像Y通过旋转旋转角“-θ₀”，从而获得旋转的失真图像Y3；以及

(e)通过缩放因数“a”改变旋转的失真图像Y3的大小，从而获得配准的图像Y*。

在上述方法中，步骤(a)包括通过计算参考图像x的像素值的和与失真图像Y的像素值的和的比例来确定缩放因数“a”。

步骤(a)还包括用零值像素填充参考图像X和失真图像Y中较小者，以使得参考图像X和失真图像Y的水平和垂直尺寸相等。

步骤(c)还包括使失真图像Y居中从而形成居中的失真图像Y1，以及进行大小改变的参考图像X2与居中的失真图像Y1之间的旋转角“θ₀”的确定。

步骤(c)包括分别形成大小改变的参考图像X2和失真图像Y的拉东变换R1和R2，以及计算拉东变换R1和R2的行的互相关来确定旋转角“θ₀”。

步骤(c)包括分别形成大小改变的参考图像X2和居中的失真图像 Y1的拉东变换R1和R2，以及计算拉东变换R1和R2的行的互相关来确定旋转角“θ₀”。

计算互相关包括：

(i)计算拉东变换R1的每行与拉东变换R2的每行之间的循环互相关集，每个循环互相关定义了各行之间的旋转偏移角“θ”；

(ii)为每行选择具有最高值的循环互相关；以及

(iii)为每行确定每个选择的循环互相关定义的旋转偏移“θ”，以及将旋转角“θ₀”设置为等于确定的旋转偏移“θ”的中值。

计算互相关还包括：

(v)将拉东变换R1和拉东变换R2的1D FFT合并成中间结果R3；

合并步骤的步骤(v)包括形成拉东变换R1的1D FFT与拉东变换 R2的1D FFT的“共轭”函数的按分量逐个作出的乘积。

步骤(vi)包括对应于各个角“θ”为第二中间结果R4中的每行确定最大值。

该方法还包括取最大值的中值，所述旋转角“θ₀”对应于中值指示的互相关所指的角。

该方法还包括对配准图像Y*执行视觉质量评估。

根据权利要求12所示的方法，其中执行以下视觉质量评估中的任一者：确定峰值信噪比(PSNR)；计算结构相似性(SSIM)指数；以及计算视觉信息保真(VIF)指数。

根据本发明的再一个方面，提供了一种将第一图像处理成与第二图像对准的配准图像Y*的图像配准方法，所述方法包括：

(a)确定第一图像与第二图像之间的缩放因数“a”；

(b)通过用缩放因数“a”改变图像中一个图像的大小来使得第一图像和第二图像的大小相等，从而生成各个大小改变的版本；

(c)使用拉东变换来确定第一图像与第二图像或第一图像与第二图像的各个大小改变版本之间的旋转角“θ₀”；以及

(d)将第一图像旋转旋转角“-θ₀”，并在将其平移成与大小改变版本的第二图像水平和垂直对准，以及在用缩放因数“a”第一次改变其大小后将其输出为配准的图像Y*。

根据本发明的再一个方面，提供了一种通过将失真图像Y与参考图像 X对其来将失真图像Y处理成配准图像Y*的***，该***包括：

处理器；

缩放因数估算模块，用于确定参考图像X与失真图像Y之间的缩放因数“a”；

旋转角确定单元，用于估算参考图像X与失真图像Y之间的旋转角“θ₀”；

图像旋转模块，用于将失真图像Y旋转旋转角“-θ₀”以形成旋转的失真图像Y3；以及

图像缩放模块，用于改变旋转的失真图像Y3的大小，以生成配准图像Y*。

该***还包括可选的居中模块，该可选的居中模块用于生成居中图像 Y1，其中所述旋转角确定单元被配置为估算参考图像X与居中图像Y1之间的旋转角“θ₀”；

该***还包括图像预缩放模块，该图像预缩放模块用于将参考图像X 缩放缩放因数“a”的倒数以成为大小改变的参考图像X2，其中所述旋转角确定单元被配置为估算大小改变的参考图像X2与居中图像Y1之间的旋转角“θ₀”。

该***还包括平移估算模块，该平移估算模块用于确定修正的图像Y3 与参考图像X之间的偏移向量“TV”，所述图像平移模块被配置为使用偏移向量“TV”平移旋转的失真图像Y3以生成配准图像Y*。

图像估算模块被配置为通过以下任一种方式执行视觉质量评估：确定峰值信噪比(PSNR)；计算结构相似度(SSIM)指数；以及计算视觉信息保真(VIF)指数。

由此，提供了改进的图像配准方法和***。

附图说明

根据下面对优选实施方式的详细描述，将会更全面地理解本发明的实施方式，应该参考附图来阅读所述优选选实施方式。所述附图不应该被理解为限制本发明，而应该仅仅是示例性的。

图1示出了图像对准***100的总框图；

图2a示出了基于仿射变换的图像调整***200.1的框图；

图2b示出了本发明第二实施方式的基于可替换的放射变换的图像调整***200.0的框图；

图3示出了根据本发明第一实施方式的基于仿射变换的重新对准方法 300的流程图；

图4包括图4a、4b和4c，示出了图3的基于仿射变换的重新对准方法的图形结果；

图5-13示出了详细描绘了用图3的重新对准方法获得的实验结果的图形化结果；

图5a、6a和7a示出了“原始”测量的DMOS(纵轴)相对于PSNR、 MSSIM和VIF分数(横轴)的散布图；

图5b、6b和7b示出了对图像进行仿射预处理之后的测量的相应散布图；

图8a、9a和10a示出了对于分别带有另外的仿射失真的LIVE图像数据库中所有失真图像的“原始”测量的DMOS(纵轴)相对于PSNR、 MSSIM和VIF分数(横轴)的散布图；

图8b、9b和10b示出了对图像进行仿射预处理之后的测量的相应散布图；

图11a、11b和11c示出了分别针对度量PSNR、MSSIM和VIF进行平移来使LIVE数据库中的选择的图像失真的实验结果的曲线图；

图12a、12b和12c示出了分别针对度量PSNR、MSSIM和VIF在 1-15度范围进行旋转来使LIVE数据库中的选择的图像失真的实验结果的曲线图；

图13a、13b和13c示出了分别针对度量PSNR、MSSIM和VIF在 0.5-1.0范围进行缩放来使LIVE数据库中的选择的图像失真的实验结果的曲线图；

图14示出了基于拉东变换的对准***的一组框图的示意图1400；

图15是本发明第三实施方式的基于第一拉东变换的图像对准*** 1500的框图；

图16示出了针对图15的***的基于第一拉东变换的方法1600；

图17示出了对图16的步骤1660“相关变换”的解释；

图18-22示出了用于表示图23中描述的本发明第四实施方式的执行的实验结果，其中图18-22中每个示出了原始图像、其失真版本和由本发明第三实施方式的方法和由现有技术方法配准的对应的失真图像；

图23示出了本发明第四实施方式的基于增强的拉东变换的图像对准***2300的框图；

图24示出了图23的***的基于第二拉东变换的方法2400；

图25示出了对图24的步骤2450“估算旋转角”的解释；以及

图26示出了图像对准***100的实施方式的结合的框图2600。

具体实施方式

本发明提出了用于对失真图像进行预处理的***和方法，通过该***和方法，可以在计算质量度量之前将失真图像与参考图像对准。用于将失真图像与参考图像对准的过程还可以不同方式被称为“进行配准”或“配准”。

对失真图像进行图像质量评估的全参考方法基于将失真图像与参考图像进行比较的特定形式。如果失真包括改变几何性质，诸如两个图像之间的尺寸差异、失真图像在其包络线内的水平和/或垂直位移或旋转、或其他变形，现有的全参考质量评估方法不能够生成与视觉评估一致的质量度量。

本发明实施方式中推荐的图像配准方法的其他应用为将两个图像进行自动比较，所述两个图像在外部观察上看上去(几乎)一样但在尺寸、旋转上不同，或者其中一个图像有噪声。

从而图像配准可以应用在遥感(地图更新)和计算机视觉中。

医学图像配准(对同一患者在诸如变换检测或肿瘤检测的不同时间点拍摄的图像数据)可以用于弥补由于例如呼吸或解剖改变造成的弹性(也称为非刚性)变形。医学图像的非刚性配准也可能用于将患者的图像数据配准成解剖图谱，诸如用于神经成像的塔莱拉什图谱(Talairach atlas)。

图像配准的还一个用处可以是在天体摄影术中对准为太空拍摄的图像。使用控制点(自动或手动输入)，可以在一个图像上进行转换从而将主特征与第二图像对准。

图像配准还可以是全景图创建的基本部分。本发明推荐的该技术可以实时地实施，并运行于嵌入式设备上，诸如相机或可拍照手机。

图1示出了图像对准***100的总框图，该图像对准***100包括提供参考图像X的参考图像源110、提供失真图像Y的失真图像源120和包括对准参数估算过程140、图像对准过程150和图像估计过程160的图像预处理***130。

参考图像X和失真图像Y两者输入到图像预处理***130中，在所述图像预处理***130中通过参考参考图像X对失真图像Y进行预处理以生成预处理的图像Y*，该预处理的图像Y*也被称为配准图像Y*。

通常，图像估计过程160将参考图像X和预处理图像Y*作为输入接收，以及生成估计结果Q0。

失真图像Y可以追溯到参考图像X，其中失真图像的失真可以是以下任意项的单独的或者结合的结果：压缩之后对图像解压缩；以旋转、平移为特征的几何变换或者图像的线性缩放；裁剪图像造成部分图像损失；由于在噪声信道上传输造成的噪声增加。

对准参数估算过程140的目的是通过参考参考图像X来分析失真图像 Y，以及生成一组对准参数AP，所述对准参数AP提供对失真的几何方面的估算。

然后图像对准过程150使用对准参数AP来将失真图像Y转换成预处理图像Y*,该预处理图像Y*也可以被称为配准图像Y*。使用对准参数AP，图像对准过程150就这样对失真图像Y的一些或所有几何失真进行补偿，但是不改变图像的视觉性质。结果，通过图像对准过程150将失真图像Y 转换成配准图像Y*，即将失真图像Y配准和与参考图像X对准。

图像估计过程160可以基于传统的视觉质量测量，其中估计结果Q0 表示诸如结构相似性(SSIM)指数或视觉信息保真(VIF)指数的质量度量。在这种情况下，当对预处理图像Y*计算估计结果Q0时,估计结果Q0 假定为失真图像Y的真实视觉质量度量，这是因为预处理动作仅修正几何失真而不改变图像的视觉内容。例如，压缩和噪声效应保持不变。

在许多其他应用中，图像估计过程160的输出Q0可以简单的为失真图像Y与参考图像X之间的视觉相似性指标，例如当失真图像X和参考图像Y两个图像为同一场景的独立的照片时。

图像预处理***130的第一实施方式基于：

-从参考图像X和失真图像Y两者提取多组关键点；

-从X和Y中匹配选定数量的多对最靠近的关键点；

-对第一仿射变换的参数求解；

-生成为第一仿射变换的倒数的第二仿射变换；

-在失真图像Y上应用第二仿射变换，从而生成预处理图像Y*。

关键点就是所谓的D.G.Lowe在计算机视觉国际期刊2004年第60卷第二本91-110页上的“Distinctive Image Features from Scale-Invariant Keypoints”中提出的尺度不变特征变换(SIFT)参数集，其用于表征图像中的特有区域。提出了对准参数估算过程140的特定实现和图像对准过程150的对应实现。

图2a示出了图像预处理***130的第一实施方式的基于仿射变换图像对准***200.1的框图，该图像预处理***130包括图1的对准参数估算过程140的实施方式的仿射变换参数估算子***210和图1的图像对准过程150的实施方式的仿射图像变换模块220。

仿射变换参数估算子***210包括像素填充模块230、X-SIFT提取模块240、Y-SIFT提取模块250、关键点匹配模块260、仿射变换生成模块 270和反向变换模块280。

在像素填充模块230中，以不改变的方式传递两个图像X和Y，除非两个图像按照其以像素测量的水平和垂直尺寸“M”和“N”来说大小不一样。如果他们大小不相等，则通过适当地用零值像素填充一个或两个图像来调整他们以及使他们的水平和垂直尺寸大小相等(对于每个尺寸，我们用其尺寸中最小的分辨率填充图像)。尺寸相等的图像将被称为均等参考图像X’和均等失真图像Y’，均等参考图像X’和均等失真图像Y’两者都是M x N像素。

为了简单起见，在下面的描述中，我们理解在任何给定情况中像素填充实际上可以是不必须的，由此均等参考图像X’和均等参考图像Y’将会持续被简称为参考图像X和失真图像Y。对图像处理概念熟悉的人很容易理解如果完成了像素填充，则可以执行处理步骤(在实施方式的描述中未示出)从配准图像Y*中移除多行或多列像素从而保证配准图像Y*的尺寸与(真实)参考图像X的尺寸匹配。

在X-SIFT提取模块240和Y-SIFT提取模块250中，分别以关键点集Kx和Ky的形式从参考图像X和失真图像Y中提取D.G.Lowe在计算机视觉国际期刊2004年第60卷第二本91-110页上的“Distinctive Image Features from Scale-Invariant Keypoints”中描述的尺度不变特征变换 (SIFT)特征。

将关键点集Kx和Ky输入到关键点匹配模块260中。

在关键点匹配模块260中，从关键点集“Kx”和“Ky”中选择预定数目的具有最接近匹配特征的向量的关键点作为匹配的关键点子集Kx’和 Ky’。选择的关键点子集Kx’和Ky’被发送到仿射变换生成模块270，用于计算正向仿射变换“FAT”的参数，所述正向仿射变换几乎非常接近参考图像X与失真图像Y之间的几何关系。将FAT发送到反向变换模块来被反向，从而生成反向仿射变换IAT，该反向仿射变换表示图1的对准参数 AP的集的第一实施方式。

在仿射变换参数估算子***210中，应用反向仿射变换IAT来将失真图像Y变换成配准图像Y*。

图2b示出了图像预处理***130的第二实施方式的基于变换的仿射变换的图像对准***200.2的框图。图2b的基于变换的仿射变换的图像对准***200.2提供了和图2a的基于仿射变换的图像对准***200.1相同的功能，除了反向变换模块280和仿射图像变换模块220被结合成单个的、更有效的仿射补偿图像生成模块290。

图3是根据本发明第一实施方式的基于仿射变换的重新对准方法300 的流程图，包括以下步骤：

310“零填充图像”；

320“提取关键他正”；

330“匹配关键点”；

340“匹配像素的数目<NK？”

350“估算仿射参数”；

360“反向变换”；

370“执行仿射变换”；

380“生成仿射补偿图像”；以及

390“填充空白像素”。

参考图像X被写作x(i,j)，其中i＝1到m,而j＝1到n，且失真图像Y 被写作y(i,j)，其中i＝1到m’,而j＝1到n’。假设参考图像X受到各种失真 (包括压缩、噪声污染和可能的仿射变换)生成了失真图像Y。预计两个图像(X和Y)不会是完全不同的图像。

在图2a的像素填充模块230中执行的步骤310“零填充图像”中，参考图像X的尺寸m和n与失真图像Y的对应尺寸m’和n’进行比较。

如果参考图像X和失真图像Y不具有相同数目的像素，则将均等图像 X和Y的尺寸设置为M和N，其中M比m和m’大，N比n和n’大。然后通过在较短的尺寸上在原始图像周围添加零像素来填充图像，以根据具体情况创建替代参考图像X和失真图像Y的均等参考图像X’或均等失真图像Y’。为了简化描述，被填充或没有被填充的图像将被称为参考图像X 和失真图像Y。

在图2a的X SIFT提取模块240和Y SIFT提取模块250中执行的步骤320“提取关键特征”中，计算缩放不变特征变换(SIFT)来分别检测和描述参考图像X和失真图像Y中的本地特征，David Lowe在计算机视觉国际期刊2004年第60卷第二本91-110页上的“Distinctive Image Features from Scale-Invariant Keypoints”中提出了详细的方法。

从参考图像X和失真图像Y两者中提取SIFT关键点。根据David Lowe的“SIFT”方法，获得关键点列表，其中每个关键点包括本地描述符的向量(即诸如颜色、纹理等对内容的视觉特征的描述)和其在图像中的位置。对于参考图像X和失真图像Y，通过该方法分别获得多个(多个集合的)Kx和Ky关键点。每个关键点由128个值的描述符向量表示，并耦合到包括4个值的位置记录上，所述4个值包括关键点的行指数、列指数、缩放因数和方向。

然后将关键点描述符标准化为表示N维欧几里得空间中的点，其中 N＝128，对应于每个描述符的128个值。

在图2a的关键点匹配模块260中执行的步骤330“匹配关键点”中，计算参考图像X的关键点与失真图像X的关键点之间的NK个接近匹配，其中通过128维空间中的距离值确定匹配的接近度。如果在关键点特征方面其距离小于到第二接近匹配距离的距离比(distRatio)倍数，则只能接受匹配。可以从0.2到0.6的范围中选择参数距离比，在下面所述的实验中使用了值0.4作为距离比。

为了有效进行计算，计算单位向量之间的点积比欧式距离可以更有效。值得注意的是，由单位矢量的点积的反余弦(acos)计算的角度之比近似于小角度的欧氏距离的比例。因此，在从最近到第二近的邻域的向量角的比例小于距离比时仅保持匹配。

在步骤340“匹配像素的数目<NK？”中，确定是否在先前步骤330中找到了至少NK个匹配的关键点对。NK的值优选设置为4。如果找到了小于4个匹配的关键点对(从步骤340输出“是”)，则不能够执行失真图像Y和仿射处理，并且功能必须停止，否则(从步骤340输出“否”)，则执行进行到步骤350。

在图2a的放射变换生成模块260中执行的步骤350“估算仿射参数”中，可以使用线性最小二乘法来确定匹配的关键点对之间的全局放射变换：

具体地，如下确定仿射变换中的参数。

参考图像中的点[i_kj_kl^T到失真图像中的点[u_kv_k]^T的仿射变换可以写成：

其中k为匹配的关键点的指数。未知变量为变换参数a₁,a₂,a₃,a₄,t₁和t₂。等式(1)可以写成：

其中，每对匹配的关键点向最左侧的矩阵和最右侧的向量提供两行。由此最右侧的向量为[u₁v₁u₂v₂…u_NKv_NK]^T。可以通过线性最小二乘法来解该等式。我们将等式(2)简写为PX＝Q。然后，我们获得：

X＝(P^TP)^-1P^TQ (3)

所以，变换参数为a₁＝X(1),a₂＝X(2),a₃＝X(3),a₄＝X(4),t₁＝X(5),和 t₂＝X(6)。

由此，步骤350“估算仿射参数”的输出为如上所示获取的放射变换参数a₁、a₂、a₃、a₄、t₁和t₂的集，其表示参考图像X与失真图像Y之间的前向(forward)仿射变换(FAT)。

通过向失真图像应用反向FAT，现在可以几何处理失真图像，即缩放、旋转、移位(平移)甚至剪切，以类似参考图像，从而将失真图像Y转换成配准图像Y*，即仿射补偿图像y^*(i^*,j^*)，该仿射补偿图像将具有与参考图像X相同的分辨率。

对于仿射补偿图像y^*(i^*,j^*)，其中i^*∈[1,m]和j^*∈[1,n]，我们可以以两种方式计算[i^*j^*]^T。

在图2a的反向变换模块260中执行的步骤360“反向变换”中，反向仿射变换执行如下：

之后是步骤370“执行仿射变换”，该步骤在仿射图像中执行。

通过图2a的变换模块220获得预处理图像如下：

y^*(i^*,j^*)＝y(u,v) (5)

但是因为计算的值i*和j*可以具有非整数值，需要解决找仿射补偿图像Y*的y^*(i^*,j^*)在整数位置(存在的，但在步骤370中未清楚地示出)的元素的问题。

在图2b的放射补偿图像生成模块290中执行的优选步骤380“生成放射补偿图像”中，图中 360和370的上述函数结合起来直接用预计在整数指数位置的像素生成仿射补偿图像y^*(i^*,j^*)。用双线性插值计算像素值如下：

对于放射补偿图像y^*(i^*,j^*)中的每个像素(i^*,j^*)，我们有

由于变换的指数值(u,v)^T为实值，可以使用双线性插值来计算值 y^*(i^*,j^*)。

当指数值u和v在1≤u≤m和1≤v≤n范围内时，下面的等式(7)提供对步骤380的过程的简短说明：

其中和分别是实值u的顶函数和地板函数。如果另一方面u和 v的值不在该范围内，即u<1或m<u或v<1或n<v，则我们分配相应的像素值为0：y^*(i^*,j^*)＝0.

在下面的MATLAB代码示例1中还示出了步骤380，即等式(6)和 (7)的示例性执行方式。注意：MATLAB(矩阵实验室)是一个数值计算环境和***编程语言，由MathWorks公司开发 (http://www.mathworks.com)。MATLAB代码示例1。

在步骤390“可选地填充空白像素”中，任何接近仿射补偿图像y*的图像边界的空白像素可以保留为空白，或者可以可选地用图像边界处非空白图像像素的镜像值填充。

例如，如果我们有以下图像：

0 0 0 0 0 0 0 0

0 0 1 2 3 4 0 0

0 0 5 6 7 8 0 0

0 0 9 1 2 3 0 0

0 0 0 0 0 0 0 0

那么，通过在空白与非空白值之间的边界的非空白像素值进行镜像操作，通过镜像值来填充缺失的像素(零值)。其提供了：

2 1 1 2 3 4 4 3

6 5 5 6 7 8 8 7

1 9 9 1 2 3 3 2

为了评价本发明第一实施方式的目的，所述像素保留为空白。

总之，基于仿射变换的重新对准方法300接收参考图象X和失真图像 y作为输入，使用结合仿射变换技术的SIFT技术，输出对准的配准图像 Y*。

图4a、4b和4c示出了图3的基于仿射变换的重新对准方法300的图形结果。图4a中示出了用于替换参考图像X的原始图像。为了论证基于仿射变换的重新对准方法300，通过对图4a的原始图像的仿射变换来创建失真图像以用作失真图像Y，并且显示了缩放、旋转、平移和剪切失真的结合。图4b中示出了失真图像Y。用图3的基于仿射变换的重新对准方法300预处理图4b的失真图像Y，来生成预处理图像，即图4c中示出的配准图像Y*。通过目测，可以看到预处理方法能够在用参考图像对准失真图像中有出色表现。

为了测试基于仿射变换的重新对准方法300的性能，使用来自H.R. Sheikh、Z.Wang,L.Cormack和A.C.Bovik发布的LIVE图像质量评估数据库发行版本2的图像执行多个实验，所述多个实验可以在互联网上从网址http://live.ece.utexas.edu/research/quality获得。该数据库包括使用五种失真从29个原始图像获得的779个失真图像。失真包括JPEG压缩、 JPEG2000压缩、高斯白噪声(GWN)、高斯模糊(GBlur)和瑞利(Rayleigh) 快衰落(FF)信道模型。另外，还通过缩放、旋转和在水平和垂直两个方向上空间移位从原始图像生成失真图像。

为了测试，与原始(参考)图像对准的预处理配准图像Y*受到视觉质量评估，如图1所示(参见图像估计过程160)。使用标准矩阵来测量参考图像X的视觉质量和失真图像Y的预处理版本，也就是配准图像Y*。还可以为测量图像视觉质量考虑三个全参考矩阵、峰值信噪比(PSNR)、平均结构相似性指数测量(MSSIM)和视觉信息保真(VIF)测量。

考虑对性能的三种测量，即差值平均意见分(DMOS)与非线性回归后的客观模型输出之间的相关系数(CC)、均方根误差(RMSE)和史比尔曼等级相关系数(SROCC)。所述SROCC定义如下：

其中N为进行评价的图像中像素的数目，以及为第i个图像在主观和客观评价之间的差异范围。SROCC为基于非参数排序的相关度量，其与任何 Z.Wang和Q.Li在图像处理IEEE学报2011年第20卷第五本1185-1198 页中描述的主观和客观分之间的单调非线性映射无关。

下面的表1-5示出了根据对没有预处理的失真图像进行视觉质量测量 (确定为“原始”方法)和对根据本发明第一实施方式进行仿射预处理的失真图像进行视觉质量测量(确定为“推荐的”方法)之间的比较，得到的概括的实验结果。

表1示出了没有加任何失真的该数据库的三个全参考度量(PSNR、 MSSIM和VIF)的实验结果。最好的执行办法用粗体强调。表1中示出的结果证明对失真图像预处理没有比原始数据库的原始方法对图像质量评估的性能产生更消极的影响，实际上其效果更好。

表2以表格形式示出了当用缩放因数0.9、旋转角0.1×180/π度、水平方向空间移位4个像素和垂直方向空间移位4个像素对LIVE数据库中每个图像进行失真时的结果。在表中，标准全参考度量PSNR、MSSIM和 VIF得到非常低的质量分，而预处理方法给出了非常好的视觉质量分，即对预处理图像执行的客观质量评估比没有预处理的图像与主观评估要一致的多。最好的执行方法用粗体强调。实际上，我们可以看到在没有预处理的情况下，现有的质量度量将会给出非常差的性能。

表2示出了当失真包括缩放、旋转和移位时对原始方法进行改进后的性能。通过图8、9和10中所示的推荐的度量和人给出的分数之间的较好的相关，可以看到这种改善的性能。图11、12、13中示出的结果证明本发明推荐的方法比原始方法受原始失真的影响较小。

下面给出了推荐的对失真图像的预处理的目的和理由。我们想要补偿放射变换，诸如不会影响整个视觉质量的旋转和缩放。实际上，参见图11 和12，可以看到推荐的方法对于小的旋转和移位是很平坦的(不变的)。但是缩放要更复杂，因为如果图像被明显地缩放从而丢失了图像细节的话我们不能指望质量保持不变。但是如图13c所示为了VIF需要方法给出合理的分数，该分数通过缩放因数的改变而有效减小。在该情况下，我们推荐的方法分数有效减小并保持合理(例如，如果缩放为0.5，分数为0.7) 而原始方法(即，没有预处理)对于0.5的缩放来说分数低于0.1，因为如果将图像缩放到原始分辨率的一半时质量不会那么差，所以在直觉上来说不太说得通。

表1.对LIVE图像数据库中所有失真测试图像进行图像质量评估的整体性能比较。

表2.对LIVE图像数据库中每个图像进行图像质量评估的整体性能比较，其中所述图像通过缩放、旋转和在水平和垂直两个方向上空间移位而失真。

表3-5示出了分别对于LIVE图像数据中的五种失真的PSNR、 MSSIM和VIF的CC、SROCC和RMSE值。有仿射预处理的分数比没有仿射预处理的分数要好。对于VIF，有仿射预处理的分数与没有仿射预处理的分数相当。

表3.对于LIVE图像数据库在非线性回归之后的CC值

表4.对于LIVE图像数据在非线性回归后的SROCC值

表5.对于LIVE图像数据在非线性回归后的RMSE值

图5-13更详细地示出了表示实验结果的图形结果。

图5a、6a和7a示出了分别对于没有仿射变换的LIVE图像数据库中所有失真图像的“原始”测量(即，没有仿射变换参数估算和仿射补偿) 的DMOS(纵轴)相对于PSNR、MSSIM和VIF分数(横轴)的散布图，而图5b、6b和7b示出了在对图像进行仿射预处理之后对测量的相应散布图。

图8a、9a和10a示出了分别对于分别有另外的仿射失真(缩放因数 0.9、旋转角0.1×180/π度、水平方向空间移位4个像素和垂直方向空间移位4个像素)的LIVE图像数据库中所有失真图像的“原始”测量的DMOS (纵轴)相对于PSNR、MSSIM和VIF分数(横轴)的散布图，而图8b、 9b和10b示出了对图像进行仿射预处理之后的测量的相应散布图；

图8-10的散布图示出了对有另外的失真的图像的相关，显示了和原始测量相比由于预处理而得到的显著改进。通过不同的横轴标度这变得更清楚。对于没有预处理的图像，VIF分数为0.004到0.2这么低。然而，在仿射预处理后，VIF分数范围可以从0到1。

图11a、11b和11c示出了分别针对度量PSNR、MSSIM和VIF进行平移来使LIVE数据库中的选择的图像失真的实验结果的曲线图。

类似的，图12a、12b和12c示出了分别针对度量PSNR、MSSIM和 VIF在1-15度范围进行旋转来使LIVE数据库中的选择的图像失真的实验结果的曲线图，以及图13a、13b和13c示出了分别针对度量PSNR、 MSSIM和VIF在0.5-1.0范围进行缩放来使LIVE数据库中的选择的图像失真的实验结果的曲线图。

图11和12中的横轴标度不包括(0,0)点，以用于更清楚地显示比较关系。实际上，在(0，0)处没有失真，所以原始图像和仿射预处理图像将会有相同的度量分数。

如在平移失真的曲线图中所示，推荐的对失真图像进行仿射预处理的方法对度量PSNR、MSSIM和VIF生成接近恒定的分数。由于小的平移不会影响视觉质量，所以当执行好的预处理时预计会得到接近恒定的分数。对于旋转失真，推荐的方法对于MSSIM和VIF生成接近恒定的分数。 PSNR值随着空间移位增加而减小，但是推荐的方法仍然比没有预处理的标准PSNR更好。再一次，由于小的旋转不会影响视觉质量，所以当执行好的预处理时预计会得到接近恒定的分数。在平移和旋转时度量分数基本保持恒定的情况显示了在采用度量之前进行推荐的预处理的效果。对于缩放失真，推荐的方法比标准方法对于PSNR、MSSIM和VIF能生成高得多的分数。因为图像质量会随着其变小而降低(变得更模糊)，所以随着缩放因数的减小，预计质量也会降低。然而，从直觉来讲图像质量不能如图所示用原始方法那样显著下降。如图所示用推荐的方法质量逐渐下降会更合适。

如图7a、7b和7c所示，其中度量值随着缩放因数减小而减小，但是用本发明的预处理方法，改变是逐渐进行的，这表示关于缩放图像质量是逐渐损失的。我们发现，当图像被仿射失真时，MSSIM和VIF是比PSNR 更好的度量。不过，用本发明的预处理步骤生成的所有三个度量比没有预处理的那些标准度量可以生成更高的分数。很容易看到，推荐的预处理方法对不同种类的失真有很好的鲁棒性。

M.P.Sampat、Z.Wang、S.Gupta、A.C.Bovik和M.K.Markey (Sampat等)在图像处理IEEE学报2009年第18卷第11本2385-2401 页的“Complex wavelet structuresimilarity:A new image similarity index”中提出了复小波域中的新的图像相似性指数(CW-SSIM)，该图像相似性指数仅对小的平移和小的旋转是不变的。在另一方面，在本发明应用中推荐的预处理方法对任何程度的平移、旋转和缩放因数都有更好的鲁棒性。这表明本发明的推荐的方法对于测量不同种类失真条件下的图像视觉质量更可行。

图像预处理***130的第一和第二实施方式为使用SIFT关键点匹配和仿射变换的“基于仿射变换的”实施方式，而下面将描述基于拉东变换的另外的实施方式。

基于仿射变换的方法和基于拉东变换的方法两者都服务于将失真图像处理成与参考图像对准的同一目的，但在性能上有一些差异。例如，基于仿射变换的方法能够修正许多几何失真，包括倾斜和翻转(颠倒)，但基于拉东变换的方法不能够纠正倾斜和翻转失真。另一方面，基于拉东变换的方法的鲁棒性在噪声很大的图像的情况下可以超过基于仿射变换的方法的性能。

基于拉东变换的图像预处理***130的实施方式是基于用于确定将失真图像Y与参考图像X“配准”的缩放因数和旋转角的方法。这些另外的实施方式有利于有效的实施，并且在存在噪声的情况下有更好的鲁棒性，但是针对某些种类的失真有局限性，对于这些种类的失真前述第一和第二实施方式可能更好(例如剪切失真，参见图4b)。这些实施方式还基于水平和垂直应用同样的缩放因数的假设。然而，在噪声环境中推荐的方法可以更好地起作用。

图14示出了基于拉东变换的对准***的一组框图的示意图1400，所述一组框图包括用于执行以下功能的多个模块：

尺寸相等功能1410，用于使得输入图像(即，参考图像X和失真图像 Y)的水平尺寸相等；

居中功能1420，用于将图像关于其质心居中；

缩放因数估算功能1430，用于确定缩放因数“a”作为参考图像X的大小与失真图像Y的大小之比；

角度估算功能1440，用于确定参考图像X与失真图像Y的方向之间存在的旋转角θ₀；

缩放功能1450，用于通过缩放因数“a”对失真图像Y或其处理后的版本进行缩放，从而使其大小大体上等于参考图像X的大小；

预缩放功能1460，用于通过缩放因数“a”的倒数对参考图像X或其处理后的版本进行预缩放，从而使其大小整体上等于失真图像Y或其处理后的版本的大小；

旋转功能1470，用于通过旋转角θ₀旋转失真图像Y或其处理后的版本；

平移功能1480，用于执行对失真图像Y的处理后的版本的平移，以将其水平和垂直与参考图像X或其处理后的版本对准。

下面描述的本发明的实施方式以各种结合和顺序使用所述一组框图 1400，以通过处理失真图像Y和参考图像X来生成配准图像Y*。

图15是图1的本发明第三实施方式的基于第一拉东变换的图像对准***1500的框图，包括以下模块：

缩放因数估算模块1502，用于确定缩放因数“a”；

图像遮蔽模块1504，用于从失真图像Y中提取遮蔽图像Y₀；

质心计算模块1506，用于根据遮蔽图像Y₀确定中心参数“c_r”和“c_c”；

图像居中模块1508，用于使用中心参数“c_r”和“c_c”对失真图像Y 进行居中并生成居中图像Y₁；以及

图像大小改变模块1510，用于使用缩放因数“a”将居中图像Y₁的大小改变成大小改变的失真图像Y₂。

大小改变的失真图像Y₂与参考图像X有相同的大小，并且已经被居中，但是还需要被旋转。

基于第一拉东变换的图像对准***1500还包括：

旋转角确定单元1511，包括拉东变换模块1512和拉东变换相关模块 1516，所述拉东变换模块1512用于分别生成大小改变的失真图像Y₂和参考图像X的拉东变换R₁和R₂，所述拉东变换相关模块1516用于根据拉东变换R₂和R₁估算旋转角θ₀。

图像旋转模块1518，其中通过旋转角θ₀来旋转大小改变的失真图像 Y₂来形成修正的图像Y₃，然而该修正的图像Y₃可能还从参考图像X横向偏离；

平移估算模块1520，其中确定修正的图像Y₂与参考图像X之间的偏移向量“TV”；

以及图像平移模块1522，其中将偏移向量“TV”应用到修正的图像 Y₃来生成配准图像Y*。

对于第三实施方式的方法，假设在失真图像Y内的实际视觉内容周围有黑色背景。还很重要的是需要注意，由于旋转角明显远离0、90、180 或270度的旋转角，所以用该方法会丢失更多的像素。离所述角度值越远，丢失的像素数越多。所以，如果对准方法之后是图像评价方法，我们离所述角度值越远，质量评估将会越不准确。

图16示出了可以实例化为基于第一拉东变换的图像对准***1500的模块的基于第一拉东变换的方法1600，其包括以下步骤：

310 “零填充图像”；

1610 “估算缩放因数a”；

1620 “计算图像Y的质心”；

1630 “将图像Y居中以生成Y₁”；

1640 “改变图像Y₁的大小以生成Y₃”；

1650 “执行拉东变换R₁和R₂”；

1660 “相关变换以生成θ₀”；

1670 “旋转图像Y₁以生成Y₃”；以及

1680 “平移修正的图像Y₃”。

在如之前所述的在步骤310“零填充图像”中使得参考图像X与失真图像Y的大小相等之后执行以下步骤序列，通过以下步骤来缩放、旋转和平移失真图像Y以变成配准图像Y*。

在缩放因数估算模块1502中执行的步骤1610“估算缩放因数”中，通过对每个图像的亮度进行平均来确定缩放因数“a”，可以通过使用傅里叶变换并划分结果来完成对每个图像的亮度的平均。这假设在创建失真图像时，没有发生改变平均亮度(诸如对比度或亮度操作)的操作：

将失真图像Y考虑为参考图像X的旋转或缩放版本。在等式(9)和 (10)中给出了每个大小为M×N的图像X和Y的傅里叶变换F₁和F2：

在极坐标***中，在等式(11)中给出了其大小关系，参见B.S.Reddy 和B.N.Chatterji在IEEE学报图像处理1996年第5卷第八本1266-1271 页中的“An FFT-basedtechnique for translation,rotation and scale-invariant image registration”，具体为：

其中a为缩放因数(水平和垂直应用同一因数)，以及θ₀为两个图像之间的旋转角，以及

θ＝tan^-1(n/m) (13)

所以，缩放因数可以获得如下：

还可以在不计算图像X和Y的2D FFT的情况下更容易地获得缩放因数a，因为

以及

现在，考虑失真图像Y还可能被平移、旋转和缩放。只要移位操作相对于失真图像Y的像素的总数移位非常少的像素，就还保持用于计算缩放因数a的方法。

在图像遮蔽模块1504执行的步骤1620“创建遮蔽图像”中，然后通过以下规则来创建遮蔽图像Y0：如果Y(m,n)>τ，Y₀(m,n)＝1；否则Y₀(m,n)＝0。在本发明的该实施方式中，发明人选择了τ＝40。

在步骤1630“将图像居中”中，使图像Y居中成为居中的图像Y₁。

遮蔽图像Y₀的质心(c_r，c_c)可以计算如下：

其中，参数c_r和c_c分别表示被称为“质心”的像素的横坐标和纵坐标。在质心计算模块1506中执行该操作。

然后将对应于Y(c_r，c_c)的位置移位到图像Y的中心，在图像居中模块1508中执行平移。如下，用等式19计算居中的图像Y₁：

Y₁(m,n)＝Y(m-c_r+M/2,n-c_c+N/2). (19)

其中m∈[1,M]且n∈[1,N]。如果(m-c_r+M/2,n-c_c+N/2)在指数范围[1,M]×[1,N]之外，则Y₁(m,n)被设置为等于0。

在图像大小改变模块1510中执行的步骤1640“缩放图像”中，然后使用缩放因数“a”来将居中的图像Y₁的大小改变为标准化图像Y₂，从而使得其与参考图像X具有相同的比例尺度。在准备估算旋转角时完成这个操作。通过Y₂(m,n)＝Y₁(m/a,n/a)来定义标准化图像Y₂。

标准化图像Y₂的指数m和n必须为整数值，但是在大小改变计算之后，m和h通常会为浮点数，然后需要将所述浮点数删节或四舍五入变成整数值。对于大小改变的图像中的像素值，然后可以使用最近的像素值，或者使用双线性插值来从多个邻近像素中确定像素值。

在以下两个步骤中，根据相应的拉东变换来获得参考图像X与标准化图像Y₂之间的旋转角θ₀。

在拉东变换模块1512中执行步骤1650“执行拉东变换”，以计算参考图像X和标准化图像Y₂的相应的拉东变换。

参考图像X的拉东变换R₁定义为：

R₁(r,θ)＝∫∫X(m,n)δ)δ-mcosθ-nsinθ)dmdn (20)

以及标准化图像Y₂的拉东变换R₂定义为：

R₂(r,θ)＝∫∫Y₂(m,n)δ)δ-mcosθ-nsinθ)dmdn (21)

其中δ()为狄拉克函数。

A.Averbuch、R.R.Coifman、D.L.Donoho、M.Israeli和Y. Shkolnisky在科学计算SIAM期刊2008年30(2)第764-784页的“A framework for discrete integraltransformations I-the pseudo-polar Fourier transform”中指出了，对于现代应用，进行离散2D拉东变换是很重要的，在最近20多年其一直是许多作者关注的对象。直到最近，拉东变换缺乏对2D离散图像的在代数上准确的、可逆的且可快速计算的连贯离散定义。所以，Averbuch等定义了用于离散2D图像的2D离散拉东变换的概念，这是基于沿着绝对值斜率小于1的线的求和。在零填充网格上使用三角插值来定义非网格位置处的值。这证明其定义提供了对连续统的准确描述，因为其随着离散化步步趋近于零而收敛成连续拉东变换。可以在Averbuch等的文章中找到关于离散拉东变换的更多细节。根据 Averbuch的文章，我们可以执行X和Y₂的离散拉东变换，分别表示为 R₁(r,θ)和R₂(r,θ),所述两个离散拉东变换的大小都为K×L。

基于两个图像(参考图像X和标准化失真图像Y₂)(除了在两个图像之间的旋转和可能的噪声以外)是相同的或者至少非常相似的前提，之后是其拉东变换R₁(r,θ)和R₂(r,θ)大体上或完全相同，除了沿着θ方向上做循环移位，即R₂(r,θ)＝R₁(r,θ+θ₀)。由参数K和L表征R₁和R₂的大小，基于对于半径“r”和角度theta所需的精度来选择参数K和L。在基于拉东变换的***的实施方式中，从使用的软件(MATLAB，见上)的“拉东”函数来选择参数K作为默认值，所述参数K与图像大小相关，大约K＝ sqrt(M*M+N*N)+5，即稍微超过用像素数表示的图像对角线的长度。为了角度的精度，L的范围选择为L＝0到179，该范围可以提供大约1度的精度。

在MATLAB文档中描述的拉东函数如下：

“R＝拉东(I，theta)返回角度为theta度的强度图像I的拉东变换R。拉东变换是沿以特定角度为方向的径向线的图像强度的投影。如果theta 为a标量，R为包括theta度的拉东变换的列向量。如果theta为向量，R 为其中每列为theta中一个角的拉东变换的矩阵”。

根据理论，如果一个图像是另一个图像的旋转版本，则两个图像将会具有同样的拉东变换，除了沿着角度尺寸的循环移位以外。在当前情况下，参考图像X和标准化失真图像Y₂有更多的差异，将不能够准确保持变换的等效性。为此，为多个r值计算角度以及取角度的中值，将会在估算角度上具有更好的鲁棒性。

在拉东变换相关模块1516执行的步骤1660“相关变换”中，通过相关方法***地比较两个拉东变换，从而相对于参考图像X来估算失真图像的最可能的取向角。

可以使用循环互相关来计算θ₀，但是这很耗时，其计算复杂度为 O(KL²)，因为对于R1和R2的每行需要循环对R2移位然后计算两行之间的互相关，这就具有O(L2)的计算复杂度。因为总共有K行，所以总的计算复杂度就为O(KL²)。

作为更有效的替换方式，可以改为使用一维(1D)FFT，该1D FFT 具有O(KLlogL)的复杂度。

下面的图表提供了通过FFT计算互相关的一般原理，第三实施方式的优选方法基于该原理：

在此我们简略地叙述互相关的定义及其实施。根据D.Nikolic、R.C. Muresan、W.,Feng和W.Singer在欧洲期刊神经学2012年1-21页的“Scaled correlation analysis:abetter way to compute a cross-correlogram”，f[n]和g[n]的离散实值函数之间的互相关为：

V.Sobolev在斯普林格(Springer)出版的作者为Michiel Hazewinkel 的数学百科全书ISBN978-1-55608-010-4中给出了f[n]和g[n]之间的卷积：

函数f[n]和g[n]的互相关等于h[n]＝f[-n]与g[n]的卷积：

类似于卷积定理，互相关满足：

FFT{f*g}＝FFT{h⊕g}

＝FFT{h}·FFT{g}

＝conj(FFT(f))·FFT{g}

其中FFT表示快速傅里叶变换，conj为复数共轭，以及点表示逐个分量相乘。通过采用反向FFT，我们可以获得其快速实施。

图17示出了步骤1660“相关变换”的扩展，包括步骤：

1710 “计算R₁和R₂的1D正向FFT来生成R₃”；

1720 “对R₃/|R₃|应用1D反向FFT来生成R₄”；

1730 “逐行确定最高互相关系数”；以及

1740 “根据阵列指数最大互相关系数的中值来估算θ₀”。

我们将R₁和R₂的第m行分别表示为r1(n)和r2(n)，其中n＝1到L，表示角度的个数，即，R₁和R₂的循环移位。

在步骤1710“计算R₁和R₂的1D正向FFT来生成R₃”中，在r₁(n) 和r₂(n)每个上面执行1D正向FFT，来生成第一中间结果R3＝r₃(u)：

r₃(u)＝FFT(r₁)·conj(FFT(r₂))

其中“conj”代表复数共轭，是按分量逐个作出的乘积。

在步骤1720“对R₃/|R₃|应用1D反向FFT来生成R₄”中，将1D反向FFT应用到r₃(u)/|r₃(u)|，即，第二中间结果R4＝r₄(n)计算为：

r₄(n)＝IFFT(r₃(u)/||r₃(u)||)。

由于r₂(n)＝r₁(n+θ₀)，我们得到

r₄(n)＝δ(n+θ₀) (24)

所以，r₄(n)在-θ₀处具有峰值，其余所有的值都为0。对于两个拉东变换之间的给定偏移角，R4＝r₄(n)的每个值表示R₁和R₂之间的互相关系数的向量(即，阵列)。

在步骤1730“逐行确定最高互相关系数”中，我们找到R4＝r₄(n)的每行上最大值的位置指数θ*(m)，即

在步骤1740“根据阵列指数最大互相关系数的中值来估算θ₀”中，通过取指向所述最大值的各个阵列指数的中值，来获得两个图像(参考图像 X与居中图像Y₁)之间的旋转角的估值为θ₀＝median(θ*(m))。

下面的MATLAB代码示例2实现了对旋转角的估算。

现在回到对图16的描述。

在图像旋转模块1518执行的步骤1670“旋转图像Y₂”中，将大小改变的图像Y₂旋转-θ₀度来补偿方向差，从而获得修正的图像Y₃。

在图像平移模块1522执行的步骤1680“平移修正的图像Y₃”中，将旋转的失真图像Y₃平移到最后位置，从而生成配准图像Y*。

对于平移，可以选择G.Varghese和Z.Wang在用于视频技术的电路和***IEEE学报2010年第20卷第七本1032-1040页的“Video denoising based on a spatiotemporalGaussian scale mixture model”中提出的全局运动补偿(MC)方法，该方法是简单、快速和可靠的方法，提供了整数像素的精度。让Y₃(m,n)为已经为缩放和旋转做过补偿的修正的图像，以及 X(m,n)为参考图像。并且，我们对X和Y₃执行2D FFT：

然后，互相关(CC)函数定义为：

k_CC(m,n)＝IFFT₂(F₁(u,v)·conj(F₂(u,v))) (27)

其中，IFFT₂为反向2D傅里叶变换，为按分量逐个作出的乘积，以及 conj表示复数共轭。估算的运动向量给出为：

在平移向量估算模块1520中计算的该估算的运动向量解释为需要用于将旋转的失真图像Y₃移动到最终位置成为配准图像Y*(即， Y*(m,n)＝Y₃(m-m_opt,n-n_opt))的平移向量。

图18-22显示了用于示出下面所述的图23中描述的本发明第四实施方式中推荐的配准方法的性能的实验结果。

图18(a-e)到图22(a-e)中每个图示出了：

-原始图像(a)；

-其失真图像(b)；以及

-通过以下方法配准的失真图像：

本发明第三实施方式的方法(c)；

通过SURF方法(d)；以及

分别通过SIFT方法(e)。

特别地，图18a、18b、18c、18d和18e显示了原始图像、其失真版本 (缩放＝0.8，旋转＝0.2弧度，平移＝(40,40)，噪声标准偏差σ_n＝10)和分别使用本发明实施方式的推荐的方法、SURF和SIFT配准的图像。

图19a、19b、19c、19d和19e显示了原始图像、其失真版本(缩放＝0.8，旋转＝0.2弧度，平移＝(4，4)，噪声标准偏差σ_n＝50)和分别使用本发明实施方式的推荐的方法、SURF和SIFT配准的图像。

在图18(a-e)和19(a-e)中，可以看到所有三个方法对于失真在从大的到中等的空间移位(分别是40,40和4,4)以及结合相对较小到中等程度的噪声(分别是σ_n＝10和σ_n＝50)方面都很有效。

图20a、20b、20c、20d和20e显示了原始图像、其失真版本(缩放＝0.8，旋转＝0.2弧度，平移＝(4，4)，噪声标准偏差σ_n＝100)和分别使用本发明实施方式的推荐的方法、SURF和SIFT配准的图像。可以看到， SURF和SIFT都没有效，但是在这种情况下推荐的方法是有效的。

图21a、21b、21c、21d和21e显示了原始图像、其失真版本(缩放＝0.5，旋转＝0.2弧度，平移＝(4，4)，噪声标准偏差σ_n＝10)和分别使用本发明实施方式的推荐的方法、SURF和SIFT配准的图像。在缩放因数为0.5的情况下，所有三个方法都有效。

图22a、22b、22c、22d和22e显示了原始图像、其失真版本(缩放＝0.1，旋转＝0.2弧度，平移＝(4，4)，噪声标准偏差σ_n＝10)和分别使用本发明实施方式的推荐的方法、SURF和SIFT配准的图像。在缩放因数为0.1的情况下，本发明实施方式的我们推荐的方法非常有效。然而，SURF 和SIFT两者没有效。从而，在输入图像缩放得非常小时本发明实施方式的推荐的方法比SIFT和SURF要好。

总之，可以看到，在σ_n＝100时，SURF和SIFT方法两者没有效，但是在这种情况下我们建议的方法非常有效。另外，在缩放因数a＝0.5时，所有三个方法都能生成好的配准(参见图21)。然而，在缩放因数a＝0.1 时，SURF和SIFT两者没有效，而我们推荐的方法非常有效(参见图22c)。将本发明第三实施方式的推荐的方法的计算复杂度与两种现有方法进行了比较。在现有技术文献中，B.S.Reddy和B.N.Chatterji在图像处理IEEE 学报1996年第五卷第八本1266-1271页的“An FFT-based technique for translation,rotation andscale-invariant image registration”中提出了一种用于估算平移、旋转和缩放参数的技术。其方法使用了6个前向2D FFT 和6个反向2D FFT，这比较耗时并且容易出错。在本发明的第三实施方式中，我们使用平均值之比来估算缩放因数，以及使用拉东变换来估算参考图像与失真图像之间的旋转。可以用恒定时间O(1)来在可重构网格上执行图像的拉东变换，如Y.Pan,K.Li和M.Hamdi在***、人和控制论A 部分-***和人类IEEE学报1999年第29卷第四本417-421页的“An improved constant-time algorithm for computing theradon and hough transforms on a reconfigurable mesh”中所示。最多，可以以2D FFT的复杂度执行拉东变换，如K.Jafari-Khouzani和H.Soltanian-Zadeh在模式分析和机器智能IEEE学报2005年第27卷第六本1004-1008页的“Radon transform orientationestimation for rotation invariant texture analysis”中所示。然后我们可以利用2DFFT来获得全局空间移位，参见 G.Varghese和Z.Wang在用于视频技术的电路和***IEEE学报2010年第20卷第七本1032-1040页的“Video denoising based on a spatiotemporalGaussian scale mixture model”。本发明第三实施方式的方法值是有3个前向2D FFT,2个反向2D FFT和2个拉东变换。所以，我们推荐的方法比上面引用的B.S.Reddy和B.N.Chatterji的文献中所述的现有技术的方法更快。在E.De Castro和C.Morandi在模式分析和机器智能IEEE学报 1987年PAMI-95卷700-703页的“Registration of translatedand rotated images using finite Fourier transforms”中，作者提出了图像配准方法，该图像配准方法需要180个反向2D FFT。而且，如果存在缩放变化，则该方法没有效。反过来，在存在缩放变换的情况下，本发明第三实施方式的方法有效，且其计算复杂度比E.DeCastro和C.Morandi描述的方法小得多。

本发明第三实施方式的方法和***比现有技术方法提供了以下优点。本发明第三实施方式的方法在噪声环境执行得非常好，而诸如D.G.Lowe 在计算机视觉国际期刊2004年第60卷第二本91-110页的“Distinctive Image Features from Scale-InvariantKeypoints”中和H.Bay、A.Ess、T. Tuytelaars和L.Van Gool在计算机视觉和图像理解(CVIU)2008年第 110卷第三本346-359页的“SURF:Speeded Up Robust Features”中提出的现有方法在噪声水平太高时是无效的。另外，本发明实施方式的方法在计算复杂度O(MNlog(MN)+KL log(L))方面很快，其中输入图像大小为 M×N且拉东图像大小为K×L。而且，本发明实施方式的方法可以估算所有三个参数(平移、旋转和缩放)，而大多数现有的方法仅计算一个或两个参数，例如K.Jafari-Khouzani和H.Soltanian-Zadeh在模式识别和机器智能IEEE学报2005年第27卷第六本1004-1008页的“Radon transform orientationestimation for rotation invariant texture analysis”和W.Wei、 S.Wang、X.Zhang和Z.Tang在信息取证与安全IEEE学报2010年第5 卷第三本507-517页的“Estimation ofimage rotation angle using interpolation-related spectral signatures withapplication to blind detection of image forgery”中提到的方法。在下面的部分得到的实验结果显示了我们推荐的方法尤其是在噪声环境中配准图像是可行的。

本发明第三实施方式的试验结果

我们在512×512的Barbara和Lena图像上进行一些实验，并且使用了LIVE图像质量评估数据库发布版本2(H.R.Sheikh、Z.Wang、L. Cormack和A.C.Bovik，“LIVE imagequality assessment database release 2,”http://live.ece.utexas.edu/research/quality.)，其包括使用五种失真从原始图像得到的779个失真图像。失真包括JPEG压缩、JPEG2000压缩、高斯白噪声(GWN)、高斯模糊(GBlur)和瑞利(Rayleigh)快衰落(FF) 信道模型。我们使用MATLAB代码来用于本发明第三实施方式的我们推荐的方法。

表6和7给出了带有平移、旋转和缩放失真以及具有添加的噪声的 Lena和Barbara图像的实验结果。在两个表中，我们估算的参数与造成失真图像的输入失真参数非常接近。另外，我们推荐的方法在噪声环境中非常有效，而在σ_n＝100时，SIFT不能找到匹配关键点，以及SURF不能为Lena图像生成精确结果。

应该指出在本发明实施方式推荐的方法中的旋转角的精度取决于点L 的数目。在我们的模拟中，我们在拉东变换的旋转角中使用L＝360的实例点。所以，精度为1度(等于0.0174弧度)。我们可以通过增加实例点的数目来增加精度，但是这会增加计算复杂度。所以，在参数精度和计算复杂度之间有折中。

表8以表格形式列出了结果，当LIVE数据库中的每个失真图像还被缩放因数0.9、旋转角0.1×180/π度、水平方向空间平移4个像素和垂直方向平移4个像素失真。对于该数据库，所有三个(推荐的、SIFT和SURF) 方法在没有添加噪声时允许比较好的估算参数。我们估计如果我们向该数据库中的图像添加噪声，则本发明实施方式推荐的方法比SIFT和SURF 效果更好，就如表6和7所示。

表9列出了对于三种方法配准一对512×512大小的图像的以秒为单位的执行时间。可以看到，本发明实施方式我们推荐的方法比SURF花的时间少，但比SIFT花的时间多。注意，用C语言写SIFT，所以比MATLAB 实现快。根据H.Bay、A.Ess、T.Tuytelaars和L.VanGool在计算机视觉与图像理解(CVIU)2008年第110卷第三本346-359页中的“SURF:Speeded Up Robust Features”，SIFT和SURF如果都以MATLAB些，则SIFT没有SURF快。所以，我们总结出本发明实施方式的我们推荐的方法比SURF和SIFT都快。另外，如我们的实验中所证明的，我们推荐的方法在噪声环境中比SURF和SIFT都表现好。

表6.具有不同噪声水平的平移(4，4)、旋转(0.1弧度)和缩放(0.9) 的Lena图像的匹配结果

表7.具有不同噪声水平的平移(4，4)、旋转(0.1弧度)和缩放(0.9) 的Barbara图像的匹配结果

表8.LIVE图像数据库中每个已经失真的图像的总体估算参数，所述图像还被缩放因数0.9、旋转角0.1×180/π度、水平方向空间平移4个像素和垂直方向空间平移4个像素

表9.用于配准一对512×512大小的图像的执行时间。SURF和推荐的方法使用MATLAB代码，而SIFT使用C代码。根据H.Bay、A.Ess、T. Tuytelaars和L.Van Gool在计算机视觉与图像理解(CVIU)2008年第 110卷第三本346-359页中的“SURF:Speeded Up RobustFeatures”，如

果SIFT和SURF都使用MATLAB，则SIFT应该比SURF慢。

公开了作为本发明第三实施方式的替换的第四实施方式，下面对第四实施方式进行描述。

图像预处理***130的第四实施方式是第三实施方式的增强版，其特征是以不同顺序来旋转和改变失真图像的大小。概述第三实施方式的操作：失真图像Y第一次居中形成居中的图像Y₁，然后改变大小(缩放)以匹配参考图像的大小生成标准化图像Y₂。通过将参考图像X和标准化图像 Y₂的拉东变换的角度相关来确定旋转角。然后旋转标准化图像Y₂来将其方向与参考图像X的方向对准，生成旋转的失真图像Y₃，以及最后水平和垂直平移旋转的失真图像Y₃已将其与参考图像X对准，从而产生配准图像Y*。

重要的是，第四实施方式的不同之处在于，不再是在计算拉东变换之前缩放失真图像Y，在第四实施方式中在计算拉东变换之前缩放参考图像 X来匹配失真图像的大小，以及在完成旋转和平移运动后在最后实现对失真图像的所需大小改变。

注意，为了便于阅读和描述，在第四实施方式中再次使用一些第三实施方式的变量名。

图23示出了是图1的图像预处理***130的第四实施方式的基于增强的拉东变换的图像对准***2300的框图，包括以下模块：

图15的缩放因数估算模块1502，用于确定缩放因数“a”；

可选的居中模块2305，用于生成居中的图像Y₁，包括图15的图像遮蔽模块1504、质心计算模块1506和图像居中模块1508；

图像预缩放模块2310，用于将参考图像X缩放缩放因数“a”的倒数以成为大小改变的参考图像X₂；

第二旋转角确定单元，包括：

拉东变换模块，用于分别生成大小改变的参考图像X₂和居中图像Y₁的拉东变换R₁和R₂；和旋转角估算块2320，用于从拉东变换R₁和R₂提取角度，消除来自角度的平移影响，以及确定旋转角“θ₀”；

图15的图像旋转模块1518，其中将居中图像Y₁旋转旋转角“θ₀”以形成旋转的失真图像Y₃，然而该旋转的失真图像Y₃可能还相对于参考图像X横向偏离；

图15的平移估算模块1520，其中确定旋转的失真图像Y₃与大小改变的参考图像X₂之间的偏移向量“TV”；

图15的图像平移模块1522，其中应用偏移向量“TV”到旋转的失真图像Y₃来生成补偿图像Y₄；

图像缩放模块2330，用于通过缩放因数“a”将补偿图像Y₄缩放成参考图像X的大小，从而生成配准图像Y*。

基于增强的拉东变换的对准***2300所用的模块1502、1512、1518、 1520和1522具有和图15的基于第一拉东变换的图像对准***1500相同的功能，但是其输入和输出数据可以改变。

基于增强的拉东变换的图像对准***2330的负责预缩放图像X的模块2310和最后将图像Y4缩放生成配准图像Y*的模块2330适于改变的拓扑结构。

包括图15的图像遮蔽模块1504、质心计算模块1506和图像居中模块 1508并具有相同功能的模块2305在基于增强的拉东变换的图像对准*** 2300中是可选的。图像居中模块的功能还可以用于在预缩放X之前将图像 X居中。

在基于增强的拉东变换的图像对准***2300的变换方式中，图14的构建模块也可以以不同方式组合，例如能够对居中图像Y₁进行预缩放而不是预缩放X，以用于确定旋转角。

图24示出了基于第二拉东变换的方法2400，其可在基于增强的拉东变换的图像对准***2300的模块中举例说明，该方法包括：

2410“估算缩放因数”；

2420“创建遮蔽图像”

2430“将图像居中”；

2440“改变参考图像X的大小来生成大小改变的参考图像X₂”；

2450“估算旋转角”；

2460“旋转失真图像”；

2470“估算平移”；

2480“平移图像”；以及

2490“改变失真图像大小”。

如果参考图像X和失真图像Y不具有相同数目的像素，则我们在先前所述的步骤310“零填充图像”中在较小的图像周围填充零。从现在开始， X和Y指的是X和Y的零填充版本。

在步骤2410“估算缩放因数”中，这等同于图16的步骤1610，如先前所述的确定缩放因数“a”。

步骤2420“创建遮蔽图像”和用于将失真图像Y居中的步骤2430“将图像居中”(分别等同于基于第一拉东变换的方法1600的图16的步骤1620 和1630)在基于第二拉东变换的方法2400中是可选的，并且在这里可能不需要，因为使用了更具有鲁棒性的拉东变换操作。然而，如果对图像的最终旋转将会造成许多像素丢失(位于帧外)，然后参数估计就不会那么准确。所以，为了增强鲁棒性，可以还选择计算质心和对准图像(步骤1620 和1630)来将失真图像Y转换成居中的图像Y₁。类似的，在此没有示出，也可以可选地对参考图像X进行居中。

在步骤2440“改变参考图像X的大小以生成大小改变的参考图像X2”中，使用计算的缩放因数“a”的倒数1/a来从参考图像X生成大小改变的参考图像X₂，从而使得大小改变的参考图像X₂与失真图像Y具有同样的尺寸。

该步骤在估算旋转角方面是有利的。我们将大小改变的参考图像表示为X₂(m,n)＝X(m*a,n*a)。

在步骤2450“估算旋转角”中，分别通过大小改变的参考图像X2与失真图像Y1的拉东变换R1和R2的输出的相关得到参考图像X与失真图像Y1之间的旋转角“θ₀”。

图25示出了步骤2450“估算旋转角”的扩展，包括步骤：

2510“计算具有鲁棒性的拉东变换”；

2520“计算R₁和R₂的1D正向FFT来生成R₃”；

2530“向R₃应用1D反向FFT来生成R₄”；

2540“计算R₄的每行中的最大值的指数”；以及

2550“根据最大值的中值计算θ₀”。

在步骤2510“计算具有鲁棒性的拉东变换R1和R2”中，分别为大小改变的参考图像X2和居中的图像Y1(或者如果跳过步骤1620和1630中的对图像居中，则用失真图像Y而不是居中图像Y₁)计算拉东变换R1和R2。通过A.Averbuch、R.R.Coifman、D.L.Donoho、M.Israeli、Y. Shkolnisky和I.Sedelnikov在科学计算SIAM期刊2008年30(2)：785-803 的“A framework for discrete integral transformations II-the2D discrete Radontransform”中提到的下面的等式(29)定义2D离散函数A(x,y) 的拉东变换：

其中δ()是狄拉克delta函数。

根据Averbuch等，我们可以执行X₂和Y(或者如果没有跳过居中步骤则为Y1)的离散拉东变换，分别表示为R₁(r,θ)和R₂(r,θ)，其中两者的大小都为K×L。

然后沿着R₁(r,θ)和R₂(r,θ)两者的r方向进行1D FFT，并取产生的矩阵的大小。这会消除沿着角度为θ的线的移位。例如，如果旋转的位置为 (x0，y0)，则在角度为θ的线上的空间移位为：

b＝x0cos(θ)+y0sin(θ).

沿着该线进行FFT和获得大小，我们随后可以消除移位b。这会改善旋转角和全局空间移位两者的精度。

基于两个图像(大小改变的参考图像X₂和失真图像Y或其居中的版本Y₁)(除了两个图像之间的旋转和可能的噪声之外)相同或者至少非常相似的前提，将得到其拉东变换R₁(r,θ)和R₂(r,θ)大体上是相同的，除了沿θ方向的旋转角之外，即R₂(r,θ)＝R₁(r,θ+θ₀)。循环互相关可以用于计算θ₀，但是这很耗时，其计算复杂度为O(KL²)，因为R₁和R₂的每行需要循环地对R₂移位并且再计算两行之间的互相关。这有O(L²)的计算复杂度。由于总共有K行，所以总的计算复杂度为O(KL²)。

作为更有效的替换方式，可以使用一维(1D)FFT作为替代，该1D FFT 具有O(KLlogL)的复杂度。

第四实施方式的优选方法的用于通过FFT计算互相关的方法和第三实施方式中所用的方法一样，请参见上面步骤1660的一般描述，但是其在消除上述移位效应方面上有所增强。

基于上述讨论，我们使用快速互相关来找到θ₀，其具有O(KL log L)的计算复杂度。

我们分别将R₁和R₂的第m行表示为r₁(n)和r₂(n)，其中n∈[1,L]。我们寻找互相关为最大值的偏移，即，找到两行r₁(n)和r₂(n)最匹配时的偏移。

在步骤2520“计算R₁和R₂的1D正向FFT来生成R₃”中，对r₁(n) 和r₂(n)每者执行1D正向FFT，来生成第一中间结果R3＝r₃(u)：

r₃(u)＝conj(FFT(r₁))·FFT(r₂),

其中，“conj”表示复数共轭，而“·”表示按分量逐个作出的乘积。

在步骤2530“向R₃应用1D反向FFT来生成R₄”中，将1D反向FFT 应用到r₃(u)，即计算第二中间结果R4＝r₄(n)＝IFFT(r₃(u))。第二中间结果R4组成拉东变换R1和R2的行的组合之间的互相关值的向量。

在步骤2540“计算R₄的每行中的最大值的指数”中，我们在获得的第二中间结果R4＝r₄(n)的每行上找到最大值的位置指数θ*(m)，即

在步骤2550“根据最大值的中值计算θ₀”中，我们取这些值的中值作为θ₀＝median(θ*(m))。

下面的MATLAB代码示例3实现了对旋转角的估算。

回到对图24的描述。

在步骤2460“旋转失真图像”中，失真图像Y被旋转-θ度以补偿方向差。旋转的结果是旋转的失真图像Y₃。

在步骤2470“估算平移”中，由X₂和Y₃的函数确定对旋转的失真图像Y₃进行的用于和大小改变的图像X₂对准的必要平移。

对于平移，可以选择G.Varghese和Z.Wang在用于视频技术的电路和***IEEE学报2010年第20卷第七本1032-1040页的“Video denoising based on a spatiotemporalGaussian scale mixture model”中提出的全局运动补偿(MC)方法，该方法是简单、快速和可靠的方法，提供了整数像素的精度。例如图像已经被做过旋转补偿，X₂(m,n)为缩放的参考图像。而且，

然后，我们可以在等式(32)中将快速互相关(FC)定义为：

k_FC(m,n)＝IFFT₂(conj(F₁(u,v))·F₂(u,v))

(32)

其中IFFT₂是反向2D傅里叶变换，conj是复数共轭，以及符号表示按分量逐个作出的乘积。用等式(33)给出估算的运动向量(平移向量)为：

在步骤2480“平移图像”中，将计算的平移向量应用到旋转的失真图像Y3，以补偿平移并生成平移的图像Y4。我们将该函数表示为 Y₄(m,n)＝Y₃(m-m_opt,n-n_opt)。

在步骤2490“改变失真图像的大小”中，使用计算的缩放因数a来改变平移图像Y₄的大小，从而使其具有和参考图像X同样的尺寸，以获得最终图像，即配准图像Y*，现在对该配准图像Y*进行了平移、旋转和缩放补偿。我们将该函数正式表示为Y*(m,n)＝Y₄(m/a,n/a)。如图23所示，本发明第四实施方式的***部分由于用于拉东变换的过程而具有很好的鲁棒性，以及与第三实施方式的***对比的是，对于任何旋转角都很有效。

图26示出了为了验证目的实施的图像对准***100的实施方式的结合的框图2600，包括：处理器2602，包括CPU2604、网络I/O***和命令接口2608；非临时性计算机可读存储介质形式的计算机存储器2610，包括存储于其上的由处理器执行的计算机可读指令；以及图像数据库2612，包括用于图像对准***100处理的参考图像和测试图像。存储在计算机存储器2610中的计算机可读指令包括软件模块，用于执行至少一个本发明的实施方式，具体为：基于仿射变换的图像对准***200.1；基于可替换的仿射变换的图像对准***200.2；基于第一拉东变换的图像对准***1500；和基于增强的拉东变换的图像对准***2300。计算机存储器2610还可以包括图像估计程序2614(例如图像估计过程160)的存储器，和软件库2614，例如在实施本发明实施方式中使用的MATLAB环境2618。

图26示出了本发明不同实施方式作为交叉使用的模块组的个体模块，例如基于仿射变换的图像对准***200.1和可替换的基于仿射变换的图像对准***200.2共享模块230、240、250、260和270，而基于仿射变换的图像对准***200.1还包括模块220和280；以及可替换的基于仿射变换的图像对准***200.2除了共享的模块外包括模块290。类似的，基于第一拉东变换的图像对准***1500和基于增强的拉东变换的图像对准***2300 共享一组模块，即1502、1512、1518、1520和1522，而基于第一拉东变换的图像对准***1500还使用模块1504、1506、1508、1510和1516；并且基于增强的拉东变换的图像对准***2300除了共享的模块外还包括模块2305、2310、2320和2330。

本发明的实施方式包括具有CPU(可以为多核CPU)和计算机可读介质(例如存储器、DVD、CD-ROM、软盘、磁带或其他存储介质)的通用和专用计算机，所述计算机可读介质具有存储于其上的计算机可读指令，当CPU执行该计算机可读指令时形成上述***的模块。可替换的，图2a、 2b、15和23的***可以包括专门的专用硬件或固件和具有计算机可读存储介质的专门的计算机的结合，所述计算机可读存储介质具有存储于其上的由CPU执行以形成上述那些***的模块的计算机可读指令。图2a、2b、 15和23的每个模块可以包括固件，或者可替换的包括存储于计算机可读存储介质中由图26中所示的处理器执行的计算机可读指令。

最后，本发明提供了用于对失真图像预处理以将其与参考图像对准的改善的***和方法，可以用于提供更精确的视觉质量评估以及用于将失真图像与参考图像配准。

虽然详细描述了本发明的实施方式，但是对于本领域技术人员来说显而易见的是，可以在所附权利要求的范围内作出改变和修改。

Claims

1.一种将失真图像Y处理成与参考图像X对准的配准图像Y*的图像配准方法，该方法包括：

(a)确定用于将所述参考图像X变换为所述失真图像Y的缩放因数a；

(b)通过所述缩放因数a的倒数来改变所述参考图像X的大小，从而生成大小改变的参考图像X2；

(c)确定从所述大小改变的参考图像X2到所述失真图像Y的旋转角θ₀；

(d)通过旋转角-θ₀旋转所述失真图像Y，从而生成旋转的失真图像Y3；以及

(e)通过所述缩放因数a改变所述旋转的失真图像Y3的大小，从而获得所述配准图像Y*。

2.根据权利要求1所述的方法，其中所述步骤(a)包括通过计算所述参考图像X的像素值的和与所述失真图像Y的像素值的和之比来确定所述缩放因数a。

3.根据权利要求1所述的方法，其中所述步骤(a)还包括用零值像素填充所述参考图像X和所述失真图像Y中较小者，以使得所述参考图像X和所述失真图像Y的水平和垂直尺寸相等。

4.根据权利要求1-3中任一项所述的方法，其中所述步骤(c)还包括相对于相对应的遮蔽图像Y0使所述失真图像Y居中从而形成居中的失真图像Y1，以及确定从所述大小改变的参考图像X2到所述居中的失真图像Y1的旋转角θ₀，遮蔽图像Y0是二进制图像，具有值1和0，1对应于失真图像Y中具有高于预定水平的像素，0对应于失真图像Y中具有不高于预定水平的像素。

5.根据权利要求1-3中任一项所述的方法，其中所述步骤(c)包括分别计算所述大小改变的参考图像X2和所述失真图像Y的拉东变换R1 和R2，以及使用所述拉东变换R1和R2来确定旋转角θ₀。

6.根据权利要求4所述的方法，其中所述步骤(c)还包括分别计算所述大小改变的参考图像X2和所述居中的失真图像Y1的拉东变换R1和R2，以及计算所述拉东变换R1和R2的行的互相关来确定所述旋转角θ₀。

7.根据权利要求5所述的方法，该方法还包括计算所述拉东变换R1和R2的行的互相关以确定所述旋转角θ₀，包括：

(i)计算所述拉东变换R1的每行与所述拉东变换R2的每行之间的循环互相关集，每个循环互相关定义了各行之间的旋转偏移角θ；

(ii)为每行选择具有最高值的循环互相关；以及

(iii)为每行确定由每个选择的循环互相关定义的旋转偏移θ，以及

将所述旋转角θ₀设置为等于确定的旋转偏移θ的中值。

8.根据权利要求1-3和6-7中任一项所述的方法，该方法还包括对所述配准图像Y*执行视觉质量评估。

9.一种通过将失真图像Y与参考图像X配准来将所述失真图像Y处理成配准图像Y*的***，该***包括：

处理器；

存储装置，具有存储于其上的用于由所述处理器执行的计算机可读指令，形成：

缩放因数估算模块，用于确定用于将所述参考图像X变换为所述失真图像Y的缩放因数a；

图像预缩放模块，用于通过缩放因数a的倒数将参考图像X缩放为大小改变的参考图像X2；

旋转角确定模块，用于估算从所述大小改变的参考图像X2到所述失真图像Y的旋转角θ₀；

图像旋转模块，用于将所述失真图像Y旋转旋转角-θ₀以形成旋转的失真图像Y3；以及

图像缩放模块，用于改变所述旋转的失真图像Y3的大小，以生成所述配准图像Y*。

10.根据权利要求9所述的***，该***还包括可选的居中模块，该可选的居中模块用于相对于相对应的遮蔽图像使所述失真图像Y居中从而生成居中的失真图像Y1，其中所述旋转角确定单元被配置为估算从所述大小改变的参考图像X2到所述居中的失真图像Y1之间的旋转角θ₀。

11.根据权利要求9所述的***，其中所述旋转角确定模块还被配置为分别计算所述大小改变的参考图像X2和所述失真图像Y的拉东变换R1和R2，以及使用所述拉东变换R1和R2来确定旋转角θ₀。

12.根据权利要求11所述的***，其中所述旋转角确定模块还被配置为计算所述拉东变换R1和R2的行的互相关以确定所述旋转角θ₀，包括：

(i)用于计算所述拉东变换R1的每行与所述拉东变换R2的每行之间的循环互相关集的部件，每个循环互相关定义了各行之间的旋转偏移角θ；

(ii)用于为每行选择具有最高值的循环互相关的部件；以及

(iii)用于为每行确定由每个选择的循环互相关定义的旋转偏移θ，

以及将所述旋转角θ₀设置为等于确定的旋转偏移θ的中值的部件。

13.根据权利要求9-11中任一项所述的***，该***还包括：

平移估算模块，用于确定所述旋转的失真图像Y3与所述参考图像X之间的偏移向量TV；以及

图像平移模块，被配置为使用所述偏移向量TV平移所述旋转的失真图像Y3来生成配准图像Y*。

14.根据权利要求9-11中任一项所述的***，该***还包括：图像估算模块处理模块，该图像估算模块处理模块用于对所述配准图像Y*执行视觉质量评估。

15.根据权利要求14所述的***，其中所述图像估算模块处理模块被配置为通过以下任意方式执行所述视觉质量评估：

确定峰值信噪比PSNR；计算结构相似性SSIM指数；以及计算视觉信息保真VIF指数。