数字图像处理高级应用:基于MATLAB与CUDA的实现(第2版) (精通MATLAB)
上QQ阅读APP看书,第一时间看更新

1.4 基于运动估计的视频倍频插帧

对于一段视频图像,通常情况下,人眼所能感知的没有明显闪烁的最小帧率是每秒24帧。但是,在具体应用中,受到图像的采集装置、存储介质和传输硬件等条件的影响,所获得的片源有可能出现帧率不高的情况。比如,早期的医学X射线影像透视设备,受到装置热容量的限制,或者为了减少对病人的射线辐射剂量,往往选择15帧拍摄,但是明显的闪烁感很容易造成视觉疲劳,给医生诊断带来很大的负担;又如,某些电视为了达到更为精细的显示效果或满足3D显示需要,原始的片源帧率不能满足需求。

为了解决上述问题,提供更好的视觉感受,需要在原始的视频图像之间插入一帧或者多帧过渡图像,其中的关键问题就是应该插入怎样的图像。如果插入一帧与前一帧或者后一帧完全一样的图像,与没有插入图像几乎没有区别。所以,所插入的图像应该是一幅新的图像,该图像基于已有的图像计算而得到,但图像之间复杂的变化情况给这个问题带来不小的难度。对于这些问题,有没有比较成熟的解决思路和方法呢?目前,用来获得图像插帧最为常见的思路是运动估计和运动补偿。

1.4.1 运动估计简介

运动估计(Motion Estimation, ME)常常和运动补偿(Motion Compensation, MC)一起使用,MEMC是目前帧间编码最为常用的方法之一。运动估计的基本思路是:首先将一幅图像划分为若干细小的单元,然后将每个细小的单元与前一幅或后一幅图像对应位置的某个区域进行对比,找到匹配度最高的位置,将该位置和原来位置之间的坐标差记录下来,这个差就是一个向量,称为运动向量(Motion Vector, MV)。利用这些运动向量,可以帮助实现视频压缩、图像插帧倍频等应用。而运动估计研究的主要内容就是如何快速、有效地获得有足够精度的运动向量。反过来,通过MV计算出或者预测出完整图像的过程,就是运动补偿。

1.4.2 运动估计的应用领域

利用运动估计实现的应用大致分为两大方向:视频压缩和视频插帧。它们分别是怎么实现的呢?

先来看看视频压缩,上面介绍过,求运动估计(ME)可以得到运动向量(MV)。这样,在保存视频的时候,在保存一张完整的图像之后,第二幅图像甚至后面若干幅图像都不用保存原始图,而只需要保存MV即可,显示的时候通过这个完整的图像和所得到的MV就可以计算还原出来。显然,一般情况下,相比起每个图像小块所要存储的像素,这个由只包含两个整数的一组MV所占用的空间远远小于一幅完整的图像,因此,在保存为MV的同时,视频也就被压缩了。事实上,常见的视频编码标准H.261、H.263、H.264以及MPEG-1/2/4等主流标编码标准中都采用了基于运动估计的压缩方法。

上述过程中,一段视频经过了求运动估计(ME)得到运动向量(MV),再经过运动补偿(MC)还原得到一幅完整的图像。那么,如果现在进行这样一个实验:对第一幅图和第三幅图计算ME,并将所得到的MV1,3与刚才第一幅图和第二幅图的MV1,2对比。通过统计,可以发现得到的MV1,3里面每个向量方向几乎都与MV1,2一样,但是长度却变为两倍。这一规律常被用来作图像插帧的MC:求出两幅图像之间的MV,并将该MV长度减小为一半MV′,再用原图和MV′计算得到新的图像,该图像即为插入帧。

以上就是运动估计应用的简单介绍,具体的方法需要根据不同的应用背景来设计。

图1.4.1为运动估计和运动补偿的流程示意图,运动估计(ME)的过程需要参考帧和当前帧参与,计算得到运动向量(MV),这些MV构成压缩编码所需的必要信息;而解码图像的过程,也即是运动补偿(MC)的过程,需要参考帧和运动向量(MV)。

图1.4.1 MEMC图像编解码处理流程

1.4.3 运动估计方法分类

20世纪70年代,运动估计方法刚刚被提出的时候,分为两个大的方向:基于像素的方法和基于块的方法。基于像素的方法主要是像素递归法,而基于块的方法主要是块匹配方法。

像素递归法根据像素间亮度的变化和梯度,通过递归修正的方法来估计每个像素的运动向量,每个像素都有一个运动向量与之对应。该方法的优点是估计精度高,缺点是解码端复杂,不利于一发多收。

相比之下,块匹配方法通用性好,改进空间大,很快便成为研究的热点和主流。而随着时间的推移,常见应用中,基于像素的方法几乎已难觅踪影。现在用到的各类常见的运动估计方法都是基于块匹配的思想,各种方法的不同之处,主要包括搜索的起始方向、匹配方法和搜索算法。

其中,对搜索算法的研究十分活跃,最近十多年,涌现了非常多的新算法和改进方法,比一开始的遍历式的全搜索方法有了长足进步。在后续小节中,将会对一些经典的搜索算法进行介绍。

1.4.4 基于块匹配方法的运动估计

基于块匹配方法的运动估计,基本思想是:先将图像划分为许多小子块(Block),然后对当前帧中的每一小块根据一定的匹配准则在相邻帧中找出当前块的匹配块,由此得到两者的相对位移,即当前块的运动向量(MV),该过程即为运动估计(ME)计算。进行运动补偿(MC)计算时,将得到的MV和原始图的对应小子块进行计算,还原出新图像。

由此得到基于块匹配方法的运动估计的一般流程为:

(1)获得前后两帧图像,前一帧习惯上称为参考帧(Referenced Image),后一帧称为当前帧(Current Image);

(2)为了方便计算,常常将图像中要做运动估计计算的部分裁剪为长宽各为8或16的倍数;

(3)将参考帧图像划分为很多边长为MBSize的小子块,比如16×16像素的小块,设定搜索区域range的大小,比如令range=8,即上下左右各8个像素;

(4)对于每一块得到的小子块,在当前帧与其对应的特定搜索区域(如图1.4.2所示)中进行搜索,通过价值函数找到匹配程度最好的那个子块位置,该位置相对于参考帧中对应子块的位置之差,即为运动向量;

图1.4.2 子块搜索区域

(5)如果要进一步还原出图像,则利用运动向量(MV)和参考帧对应子块,将子块按照MV所指的位置移动过去,所得图像即为还原出来的当前帧图像,该过程也即成为运动补偿。

基于块匹配的运动估计,在实际应用中面临下面的一些难点:

首先,准确度问题。运动估计(ME)和运动补偿(MC)之后所得到的图像应当足够准确,也就是说计算得到的运动向量(MV)应当足够准确。这里说的“准确”,包含两方面的意义:每个小块的MV是否准确以及所有小块的MV的准确率是否足够高。由于搜索区域大小的有限性、小块之间匹配算法的局限性,对于运动速度较大、亮度变化较大或是图像中图形不连续的部分,在计算MV的时候,有一定出错的可能性。

其次,计算效率问题。将一幅图像划分为很多小块,每个小块都需要进行各自独立的搜索,这对搜索过程和匹配算法的效率提出了很高的要求。最初的遍历式搜索方法能够保证找到搜索区域中匹配度最高的点,但是动辄每帧图像十几甚至几十秒的运算速度,也让实时应用难以实现。

为了又快又准地进行运动估计计算,研究人员主要从搜索算法寻求突破,研究的方向主要有搜索起始点的方向、能够尽可能减少搜索点的搜索路径、MV在时间和空间上的关联预测等,来对算法进行改进。

对于经过MEMC计算和还原的图像,评价其品质好坏的方法是计算其峰值信噪比(PSNR)。

PSNR即到达噪声比率的峰值信号,它是原图像与被处理图像之间的均方误差相对于(2 n-1)2的对数值(信号最大值的平方,n是每个采样值的比特数),它的单位是dB。PSNR是一种评价图像的客观标准,PSNR值越大,就代表失真越小。

式(1.4.1)中的MSE是最小均方差(Mean Squared Error),其计算公式如下:

这里N是子块的边长CijRij是各个当前子块和参考子块中对比的像素值。

在前面介绍了运动估计的难点,归结起来就是怎样才能又准又快地计算出MV,而目前的研究热点主要是针对这些难点进行改进的。这些改进主要包括搜索起始点的位置和方向、匹配算法和搜索算法,其中搜索算法是近几年来研究进展最快的,严格地说,起始点位置和方向也属于搜索算法中的一部分。

1.4.5 相关概念

起始点位置和方向是每个小子块在搜索一开始可以做的一个选择,选择搜索区域中从哪个位置开始搜索,朝哪个方向开始搜索。最初的遍历式的全局方法中并不考虑这个问题,只是单纯地将搜索区域中所有点遍历一遍,然后找出匹配度最高的位置。

但是随着研究的深入以及终止条件的引入,人们发现,对于同一个小子块的运动总有着一定的规律:这些运动方向往往与当前帧中周围位置的子块(如图1.4.3所示)或者前面几帧中相同位置的子块(如图1.4.4所示)有着类似的运动趋势。于是,就利用这些条件来预测搜索起始点的位置和方向,从而提高搜索的效率和速度。

图1.4.3 利用周围运动向量信息预测搜索起始点

图1.4.4 利用前面若干帧向量信息预测搜索起始点

1.4.6 匹配方法:价值函数

参考帧中的子块在当前帧中的搜索区域搜索,要找到最佳匹配位置,每一次搜索会计算出一个匹配值,而最佳区域就是满足条件的峰值点。这里,如何计算每次匹配的值,就要用到价值函数,也就是所谓的匹配方法。在图像匹配中,可以用作价值函数的方法包括:绝对平均误差(MAD)、绝对差和(SAD)、最小均方误差(MSE)和归一化互相关(NCC)。下面对这几种匹配方法进行简单介绍。

1.绝对平均误差(MAD)

绝对平均误差也称平均绝对失真,是对两个子块间对应像素误差的描述。其计算公式见式(1.4.3):

其中,(i,j)是位移向量,(m,n)为当前子块左上角坐标,fkfk-1分别为当前帧和上一帧的灰度值,M×N为子块大小。若在某一点(i0,j0), MAD(i0,j0)达到最小,则该点为要找的最优匹配点,对应的块为最优匹配块。

2.绝对差和(SAD)

绝对差和是对两个子块间对应像素误差之和的描述,与MAD的区别在于绝对差和不求平均值,因此,比MAD少一个计算步骤,而所得到的匹配结果与MAD是一样的。其计算公式见式(1.4.4):

与MAD类似,其中,(i,j)是位移向量,(m,n)为当前子块左上角坐标,fkfk-1分别为当前帧和上一帧的灰度值,M×N为子块大小。若在某一点(i0,j0), SAD(i0,j0)达到最小,则该点为要找的最优匹配点,对应的块为最优匹配块。

3.最小均方差(MSE)

最小均方差也是对两个子块间对应像素误差之和的描述,与SAD相比,在计算中多了相乘的步骤,使得计算量大大增加。其计算公式见式(1.4.5):

其中,(i,j)是位移向量,(m,n)为当前子块左上角坐标,fkfk-1分别为当前帧和上一帧的灰度值,M×N为子块大小。若在某一点(i0,j0), MAE(i0,j0)达到最小,则该点为要找的最优匹配点,对应的块为最优匹配块。

4.归一化互相关(NCC)

归一化互相关是模板匹配中常用的方法,用到能量函数的思想,该方法有很多改进算法,包括扩展到频率域等。常用于尺寸较大的图像的匹配,匹配准确率较高。

其中,是当前帧子块的平均值,是参考帧中当前帧对应位置子块的平均值。若在某一点(i0,j0), NCC(i0,j0)达到最大,则该点为要找的最优匹配点,对应的块为最优匹配块。

对比上述几种常见匹配方法,计算最为复杂、准确率最高的是NCC,最简单的是SAD,其匹配准确率在图像很小的时候也比较不错。考虑到资源占用和性能要求,往往选择SAD作为匹配算法。实际使用中,经典的非对称十字多层次六边形搜索算法UMHexagonS中,就以SAD作为匹配算法。

考虑到图像之间亮度可能发生变化,因此,具有较强鲁棒性且运算速度较快的价值函数也成为研究的需要。

1.4.7 搜索算法

怎样尽可能快地找到匹配点?除了匹配方法,更重要的就是寻找最佳匹配点的方法。这一部分将介绍若干经典的搜索算法。

1.穷尽搜索法(ES)/全搜索法(FS)

穷尽搜索法(Exhaustive Search, ES),也称全搜索法(Full Search, FS),顾名思义,就是将待搜索区域每个可能位置遍历一遍,如图1.4.5所示,找到最佳匹配位置。这样的结果就是它不会遗漏并且能够发现搜索区域中最好的可能的匹配,并且获得所有块匹配算法中最高的PSNR值。

图1.4.5 全搜索法搜索的点覆盖整个搜索区域

后面提出的各类快速块匹配算法都是试图在尽可能少的计算中做到相同的PSNR。ES的优点是找到的MV精度高,而其缺点也很明显,就是搜索窗越大,越耗费资源。

2.三步搜索法(TSS/3SS)

三步搜索法(Three Step Search, TSS/3SS)可以追溯到20世纪80年代中期,是最早尝试快速块匹配的方法之一,大致思想如图1.4.6所示。假设对于一个普通的设为range=7的搜索区域参数,该算法的起点是搜索区域的中心,并且设置“搜索步长”为S=4;然后搜索8个位置,环绕位置(0,0)±S个像素;从这9个位置搜索,找到最小值的那一点;然后以此点作为新的搜索原点;然后设置新的搜索步长SS/2,重复类似搜索两个周期,直到S=1;在该点发现价值函数的最小值,子块在这点位置被认为是最好匹配。

图1.4.6 三步搜索法流程,运动向量是(5, -3)

图1.4.6中,●为第一步,▲为第二步,■为第三步的检索点。

TSS背后的思想是:在每一个子块的运动中产生的错误构成的曲面是单峰值的。一个单峰值的曲面是一个碗形曲面,这样由价值函数产生的值在全局极小值中是单调递增的。

对于range=7的情况,ES将要搜索225个子块,而TSS只需要计算25个子块,可见,该方法可以显著减少计算量。

3.新三步搜索法(NTSS)

三步搜索法(TSS/3SS)对于运动估计使用了固定的检索模式,容易错过小的运动。相比之下,新三步搜索法(New Three Step Search, NTSS)的步骤如图1.4.7所示。在第一步中,除了搜索原点之外的16个点(8个●和8个■),找到价值函数的最小值点。这些多出来的搜索位置,8个点是距离原点S=4的点(与TSS类似),另外8个点是距离原点S=1的点。如果最小值在原点,那么搜索就停在这里,并且运动向量就为(0,0);如果最小值为S=1是8个点中的任意一个点,则改变搜索原点到该点,并检查它邻近的值。按照不同的点,或许能够检查5个◆点(见图1.4.8(b))或3个▲点就结束(见图1.4.8(c))。给出最小值的位置就是最接近的匹配,运动向量就被设置为那个位置。从另一方面来说,如果最小值在第一步之后,是S=4中的8个点其中之一,那么就按照TSS的标准方法。因此,尽管这个方法在理想情况下需要对每个子块最少检测17个点,但是也有可能在最差情况下检索33个点。

图1.4.7 新三步搜索法(NTSS)

如图1.4.7所示,●是三步搜索法(TSS/3SS)中第一步的检索点,■是NTSS中加入的属于第一步的多余的8个点。▲和◆是NTSS的第二步,当第一步中最小值在8个相邻点中一个的中心时,后面要检索的点数分两种情况:5个◆点或者是3个▲点。

新三步搜索法改进了三步搜索法的结果,通过增加一个有偏向性的中心搜索方案和一个中途停止条件来减少计算量。这是最早被广泛接受的快速算法之一,被频繁用于实施早期的编码标准,比如MPEG 1和H.261。

4.四步搜索法(FSS/4SS)

与NTSS类似,四步搜索法(Four Step Search, FSS/4SS)同样采用中心偏置搜索,并且有中途停止条件。如图1.4.8所示,FSS在第一步中设置了一个固定搜索步长S=2,而不考虑搜索区域参数r的值是多少。因此,它在5×5的窗口中搜索9个位置。如果最小值在搜索窗中心被发现,搜索直接跳到第四步;如果最小值为8个位置中的一个,则将该点作为搜索原点并进入到第二步。搜索窗仍然维持在5×5。也许可以在检索到3个点或5个点时结束,这取决于最小值的位置。

图1.4.8 搜索子过程可能出现的各种情况

四步搜索法的搜索模式如图1.4.9所示。此外,如果最小值的位置在5×5搜索窗的中心,可跳到第四步或者继续第三步。第三步和第二步一模一样,在第四步,窗口尺寸缩小到3×3,比如S=1。最小值的位置就是最佳匹配宏块,运动向量就设置成指向该位置。

图1.4.9 四步搜索法

图1.4.8中显示了一个处理例子,●是第一步检索的点,■是第二步检索的点,▲是第三步检索的点,◆是第四步检索的点。这个搜索算法最好情况下检索17个点,最差情况下检索27个点。

5.菱形搜索法(DS)

菱形搜索法(Diamond Search, DS)和FSS几乎一样,但是搜索点模式从方形换成了菱形,并且对于算法的搜索步骤数没有限制。DS使用了两种不同的混合的模式,一种是大菱形搜索(LDSP),另一种是小菱形搜索(SDSP)。

这两种模式DS的执行过程如图1.4.10所示,黑色点为大菱形搜索模式,灰色点为小菱形搜索模式。该例中得到运动向量(-4, -2)总共5步:4步使用了LDSP,1步使用了SDSP。

因为搜索模式既不太小也不太大,事实上对于搜索的步数并没有限制,该算法可以比较精确地找到全局最小值。

6.非对称十字多层次六边形搜索法(UMHexagonS)

非对称十字多层次六边形搜索法(Unsymmetrical-Cross Multi-Hexagon Search,UMHexagonS)是目前比较主流的搜索算法,因其快速和高精度的特点,被H.264编码标准所采用。需要说明的是,UMHexagonS方法并非最快速的方法,相比其他简单算法,计算复杂度反而有所上升。

图1.4.10 菱形搜索(DS)过程

图1.4.11显示了UMHexagonS算法搜索过程,总体上分为四步,其中第一步为初始判断,之后根据所得结果与阈值进行判断,按照SAD值的满意程度选择进入Step2~4中的任意一步。

图1.4.11 UMHexagonS算法搜索过程示例

Step1(图1.4.11中的●):搜索起点预测。具体预测模式包括中值预测(MP)、上层预测(UP)、时间域邻近参考帧预测(NRP)、上一参考帧对应块(CP)、原点预测(OP)这五种方法来预测当前块的运动向量(MV),将该点作为搜索起点;此外,还可以利用SAD的相关性来进行预测。

当上述几种预测完成之后,找到最佳预测搜索起点,根据所得SAD值的满意程度,选择跳转进入Step2~4中的一个。

Step2:不满意搜索模式。

Step2.1(图1.4.11中的△):以预测得到的搜索起点为中心,用非对称十字形搜索模板进行搜索;获得当前最佳点,判断此处SAD值是否属于满意或很满意,跳到相应的步骤3或步骤4或继续搜索。

Step2.2(图1.4.11中的〇):以当前最佳点为中心,在边长为5的方形区域中进行ES搜索;获得当前最佳点,判断此处SAD值是否属于满意或很满意,跳到相应的Step3或Step4或继续搜索。

Step2.3(图1.4.11中的□):用每次扩大最小六边形直径一倍的大六边形模板进行搜索,直至搜索到的SAD能符合相应阈值而进入Step3或Step4的点为止;否则结束Step2的搜索进入Step3。

Step3(图1.4.11中的▽):满意搜索模式。以当前最佳点为中心进行小六边形搜索,如果当前最佳点在六边形边上,则将搜索中心移动到该点,如此反复,直至最佳点为六边形的中心。

Step4(图1.4.11中的■):很满意搜索模式。以当前最佳点为中心进行菱形搜索,直至最佳点为菱形的中心。

1.4.8 实际应用举例

本小节将对医学X射线透视视频使用运动估计的方法进行插帧,看看是否能让视频效果有所提升。

实验1所用图像为一段由某型号X射线机采集的每秒30帧的视频,为了让结果更为准确直观,以一张分辨率检测卡为对象进行实验(如图1.4.12所示),获得的原始视频包含56幅灰度值为0~255的灰度图像。

图1.4.12 分辨率检测卡实验

为了更好地理解运动估计方法的局限性,提供了实验2,该实验所用图像为X射线机采集的每秒15帧的视频,如图1.4.13所示,采集对象为人体胸正位透视图像,原始视频包含107幅灰度值为0~255的灰度图像。

实验3所用视频与实验2相同,但是在运动估计计算之前,对亮度进行了补偿,用来与实验2效果进行对比。

实验1和实验2流程框图如图1.4.14所示,所用方法为全搜索法(ES)(参看MATLAB程序main_MEMC_ES.m),具体包括以下环节:

(1)设置初始参数:子块大小MBSize,搜索范围range。

(2)读入名为XRayVideo.avi的原始图像视频。

图1.4.13 人体胸正位实验

(3)从视频中取出两相邻帧图像,前一帧为参考帧ImageRef,后一帧为当前帧ImageCur。

(4)将两帧图像分割出边长为MBSize的子块,进行运动估计计算。对于当前帧中的每个子块,通过全搜索法,计算其在参考帧中对应位置的子块的价值函数(MAD)的值,存放在一个和搜索区域对应大小的矩阵Costs中,再找到Costs矩阵中的最小值点以及该点所在的位置dy和dx,这个dy和dx就为该子块的运动向量(MV)。

(5)重复步骤(4),直到计算完当前图像中的所有MV,将所得结果存放在MotionVect中。

(6)计算插入帧图像:将所得MV的值变为原来的一半MV′,再利用新的MV′进行运动补偿MC计算,所得的图像即为插入帧图像。

(7)所得结果保存在result文件夹中。

图1.4.14 实验1和实验2流程框图

【例1.4.1】图像插帧程序。

    *******************************************************************************
    % 图像插帧程序
    % 说明:
    %     本程序读入根目录下的一段名为XRayVideo的医学图像视频,使用MEMC中的ES
    %     方法进行图像插帧,得到的新视频XRayVideoNew帧数和帧率为原来两倍,结果存放
    %     在./result文件夹下
    % 注意:
    %       本程序运算速度很慢,需耐心等待
    close all;
    clear all;
    % 预设参数
    MBSize = 16;                                                 %子块大小
    Range  = 8; %搜索范围
    % 要读写的文件名
    % ReadFileName  = './XRayVideoCard';                       %分辨率卡透视
    ReadFileName  = './XRayVideoMan';                          %人体胸正位透视
    WriteFileName =[ReadFileName, 'New'];
    SaveDirectory= './result/';
    % 读出视频信息
    ReadVideoObject = VideoReader([ReadFileName, '.avi']);
    FrameRateOrigin= ReadVideoObject.FrameRate;                  %帧率
    FrameNums = ReadVideoObject.NumberOfFrames;                  %帧数
    % 创建要写入的视频对象
    WriteVideoObject = VideoWriter([SaveDirectory, WriteFileName, '.avi']);
    WriteVideoObject.FrameRate = FrameRateOrigin * 2;
    open(WriteVideoObject);
    % 循环计算插入帧,保存为新视频文件
    for i = 1 :FrameNums - 1
        imgNumO= 2 * i - 1; %存放原始图序号
        imgNumI = 2 * i;                                          %存放插入图序号
        % 从视频中读出要处理的图像
        ImageRef0 = read(ReadVideoObject, i);
        ImageCur0 = read(ReadVideoObject, i + 1);
        ImageRef  = double(rgb2gray(ImageRef0));
        ImageCur  = double(rgb2gray(ImageCur0));
        % 生成并保存插入帧
        [MotionVect]=ME_ES(ImageCur, ImageRef, MBSize, Range);   %计算运动向量
        ImageComp    =MC(ImageRef, MotionVect, MBSize);           %计算插入帧
        % 保存JPG图像
        imwrite(uint8(ImageRef), [SaveDirectory, ReadFileName, int2str(imgNumO), '.jpg']);
        imwrite(uint8(ImageComp), [SaveDirectory, ReadFileName, int2str(imgNumI), '.jpg']);
        % 将图像写入AVI视频
        writeVideo(WriteVideoObject, ImageRef/255);
        writeVideo(WriteVideoObject, ImageComp/255);
   end
   close(WriteVideoObject);
   *******************************************************************************

经过以上步骤,实验1得到一段110帧图像的分辨率检测卡的视频,实验2得到一段212帧图像的人体胸正位视频,两个视频中奇数帧为原始图像,偶数帧为插入图像。

实验3(参见例1.4.2)的流程框图如图1.4.15所示,其步骤与实验1、2相比,首先将搜索范围range扩大到10像素,然后在步骤(3)后面新增加了两帧图像之间的亮度倍数计算,粗略地将两幅图像亮度相除,得到一个比例系数,再将该比例系数乘以参考帧图像的亮度,所得新图像再同当前帧图像进行运动估计计算得到MV。而最后插入帧图像步骤与上述步骤(6)、(7)相同。

图1.4.15 实验3流程框图

【例1.4.2】 修改后的图像插帧程序。

    *******************************************************************************
    % 图像插帧程序
    % 说明:
    %       本程序读入根目录下的一段名为XRayVideo的医学图像视频,使用MEMC中的ES
    %       方法进行图像插帧,得到的新视频XRayVideoNew帧数和帧率为原来两倍,结果存
    %       放在./result文件夹下
    % 注意:
    %       本程序运算速度很慢,需耐心等待
    close all;
    clear all;
    %%%%%%%%  预设参数  %%%%%%%%
    MBSize = 16;                                               %子块大小
    Range= 10; %搜索范围
    %%%%%%%%  要读写的文件名  %%%%%%%%
    % ReadFileName  = './XRayVideoCard';                     %分辨率卡透视
    ReadFileName  = './XRayVideoMan';                        %人体胸正位透视
    WriteFileName =[ReadFileName, 'New'];
    SaveDirectory= './result/';
    %%%%%%%%  读出视频信息  %%%%%%%%
    ReadVideoObject = VideoReader([ReadFileName, '.avi']);
    FrameRateOrigin= ReadVideoObject.FrameRate;                %帧率
    FrameNums = ReadVideoObject.NumberOfFrames;                %帧数
    %%%%%%%%  创建要写入的视频对象  %%%%%%%%
    WriteVideoObject = VideoWriter([SaveDirectory, WriteFileName, '.avi']);
    WriteVideoObject.FrameRate = FrameRateOrigin * 2;
    open(WriteVideoObject);
    %%%%%%%%  循环计算插入帧,保存为新视频文件  %%%%%%%%
    for i = 1 :FrameNums - 1
        imgNumO= 2 * i - 1; %存放原始图序号
        imgNumI = 2 * i;                                      %存放插入图序号
        %%%%%%%%  从视频中读出要处理的图像  %%%%%%%%
        ImageRef0 = read(ReadVideoObject, i);
        ImageCur0 = read(ReadVideoObject, i + 1);
        ImageRef  = double(rgb2gray(ImageRef0));
        ImageCur  = double(rgb2gray(ImageCur0));
        %%%%%%%%  亮度补偿  %%%%%%%%
        [Row, Col]= size(ImageRef);
        ImageRefCut = ImageRef(Row/4 + 1 :3 * Row/4, Col/4 + 1 :3 * Col/4);
        ImageCurCut = ImageCur(Row/4 + 1 :3 * Row/4, Col/4 + 1 :3 * Col/4);
        SumPixelRef = sum(ImageRefCut(:));
        SumPixelCur = sum(ImageCurCut(:));
        PixelTimes  = SumPixelCur/SumPixelRef;
        ImageRef1 = ImageRef * PixelTimes;
        %%%%%%%%  生成并保存插入帧  %%%%%%%%
        [MotionVect]=ME_ES(ImageCur, ImageRef1, MBSize, Range); %计算运动向量
        ImageComp    =MC(ImageRef, MotionVect, MBSize);        %计算插入帧
        %%%%%%%%  保存JPG图像  %%%%%%%%
        imwrite(uint8(ImageRef), [SaveDirectory, ReadFileName, int2str(imgNumO), '.jpg']);
        imwrite(uint8(ImageComp), [SaveDirectory, ReadFileName, int2str(imgNumI), '.jpg']);
        %%%%%%%%  将图像写入AVI视频  %%%%%%%%
        writeVideo(WriteVideoObject, ImageRef/256);
        writeVideo(WriteVideoObject, ImageComp/256);
    end
    close(WriteVideoObject);
    *******************************************************************************

例1.4.1和例1.4.2所用到的几个子函数如下:

    *******************************************************************************
    function Cost = CostFuncMAD(CurBlk, RefBlk, MBSize)
    % 计算两个子块的平均绝对误差(MAD)
    % hehao20131216
    % 输入
    % CurBlk:当前帧子块
    % RefBlk:参考帧子块
    % MBSize :子块边长
    % 输出
    % Cost :两个子块之间的MAD值
    Err = 0;
    for i = 1 :MBSize
        for j = 1 :MBSize
            Err = Err + abs((CurBlk(i, j)- RefBlk(i, j)));
        end
    end
    Cost = Err/(MBSize *MBSize);
    *******************************************************************************
    *******************************************************************************
    function ImageComp =MC(Image, MotionVect, MBSize)
    % 生成运动补偿图像
    % 输入
    %   Image:参考图
    %  MotionVect :运动向量
    %  MBSize :子块边长
    %
    % 输出
    %   ImageComp:运动补偿图像
    [Row, Col]= size(Image);
    ImageC = Image;
    MBCount = 1;
    for i = 1 :MBSize :Row -MBSize + 1
        for j = 1 :MBSize :Col -MBSize + 1
            dy=MotionVect(1, MBCount);
            dx =MotionVect(2, MBCount);
            RefBlockY = i + round(dy/2);
            RefBlockX = j + round(dx/2);
            ImageC(i :i +MBSize - 1, j :j +MBSize - 1)= Image(RefBlockY:RefBlockY+
    MBSize - 1, RefBlockX :RefBlockX +MBSize - 1);
            MBCount =MBCount + 1;
        end
    end
    ImageComp = ImageC;
    *******************************************************************************
    *******************************************************************************
    function[MotionVect]=ME_ES(ImageCur, ImageRef, MBSize, Range)
    % 全搜索法搜索运动向量
    % 输入
    %   ImageCur:当前帧
    %   ImageRef :参考帧
    %  MBSize :子块边长
    %   Range:搜索范围
    % 输出
    %  MotionVect :运动向量MV
    [Row, Col]= size(ImageRef);
    Vectors = zeros(2, Row* Col/MBSize 2);
    Costs   = ones(2 * Range + 1,2 * Range + 1)* 256;
    MBCount = 1; %MV计数
    for i = 1 :MBSize :Row-MBSize + 1                        %分割出子块
        for j = 1 :MBSize :Col -MBSize + 1
            for m= -Range:Range                             %处理每个子块搜索区
                for n = -Range :Range
                    RefBlockY = i + m;
                    RefBlockX = j + n;
                    if(RefBlockY   1 ‖ RefBlockY + MBSize - 1   Row ‖ RefBlockX   1 ‖
    RefBlockX +MBSize - 1  Col)
                        continue;
                    end
                    Costs(m + Range + 1, n + Range + 1)= CostFuncMAD(ImageCur(i :i +
    MBSize - 1, j :j +MBSize- 1), ImageRef(RefBlockY:RefBlockY+MBSize- 1, RefBlockX:
    RefBlockX+MBSize- 1), MBSize); %计算MAD矩阵
                end
            end
            [dy, dx]=MinCost(Costs); %找到MAD最小值点作为运动向量的值
            Vectors(1, MBCount)= dy- Range- 1; %对于每个子块的运动向量
            Vectors(2, MBCount)= dx - Range - 1;
            MBCount =MBCount + 1;
            Costs = ones(2 * Range + 1,2 * Range + 1)* 256;
         end
     end
     MotionVect = Vectors;
     *******************************************************************************
     *******************************************************************************
     function[dy, dx]=MinCost(Costs)
     % 寻找最小值位置
     % hehao20131216
     % 输入
     %   Costs :包含一个子块在搜索区中各个位置价值函数值的矩阵
     % 输出
     %   dy:运动向量竖直方向的值
     %   dx:运动向量水平方向的值
     [Row, Col]= size(Costs);
     Min = 256;
     for i = 1 :Row
         for j = 1 :Col
             if(Costs(i, j)  Min)
                 Min = Costs(i, j);
                 dy= i;
                 dx = j;
             end
         end
     end
     *******************************************************************************

这里所用为Intel i7-2860QM 2.5GHz处理器,4GB内存,第1、2两组实验中每帧图像的计算时间都是每幅图约13.6s,第3组实验中每帧图像计算时间约为22s。

实验1:成功的分辨率检测卡插帧实验XRayVideoCard range=8。

对比原始图像和插帧之后的图像,可以发现,经过插帧的图像人眼感知到的分辨率线条更为精细,同时白色十字的跳跃感更小。说明插帧之后的图像能够改善人眼视觉感受。

实验2:失败的人体胸正位插帧实验XRayVideoMan range=8。

该视频原始图像为15帧,在人体进行平移的时候闪烁感强烈,而经过插帧之后的图像闪烁感减少,但是会有明显的错误。所产生的错误都是由运动估计方法本身的一些局限性造成的,这些局限性将在下节进行介绍。

实验3:效果改善的人体胸正位插帧实验XRayVideoMan range=10 LumComp。

经过亮度补偿之后,亮度变化时抖动有所改善。

实验2的结果很不理想。那么,算法中的哪些局限性造成实验的失败呢?对比实验1和实验2,以及在其他实验中出现的问题,可以发现运动估计包括以下一些局限性:

(1)整体亮度变化时容易失败。实验1的灰度一直保持整体均匀,而实验2的整体灰度在不断变化,每隔特定时间就会从亮到暗跳变几帧。这是由X射线机的成像原理决定的,射线发射装置会依据图像整体亮度进行调节,于是,在人体移动的过程中就会发生照射剂量的改变,图像整体灰度也因此改变。对于前后两张图像,物体不动时,亮度改变,运动估计时也常常计算失败。

如图1.4.16所示,每当亮度发生变化时,运动估计往往出现大范围失败。

图1.4.16 实验2匹配失败位置

(2)图形运动速度大于搜索范围会失败。为了节约计算量,搜索范围range不能设置得太大,否则计算速度将非常缓慢。但是,如图1.4.17所示,如果两帧图像间的物体或者物体中的某一部分运动速度非常快,超过了range的范围时,那可以肯定图像是无论如何也匹配不上。因此,运动速度或者局部运动速度过快,都会导致匹配失败。

(3)计算速度慢。由于需要计算的小块很多,计算速度过于缓慢也是运动估计应用的一大问题。

(4)马赛克现象,也就是说子块之间出现非连续性,严重影响视觉质量。

(5)对于两帧图像之间图像不连续会失败,即某种图形在前后两帧图像中突然出现或者突然消失时,插入图像也会出现错误。

图1.4.17 移动超过搜索范围造成匹配失败

对于运动估计的研究一直处于不断进步之中,研究的热点主要集中在如何又准确又快速地获得运动向量,各类更快速、搜索范围更大的搜索算法层出不穷。这些算法的提出和计算机硬件的升级,正逐渐让运动估计的实时应用变为可能。

但是,在实际应用中,一些问题仍然是运动估计方法所急需解决的。这其中比较常见的有:

(1)两帧图像亮度变化。目前算法的价值函数在子块亮度发生变化之后,很有可能搜索到错误的位置。因此,改进的方向可以是改进匹配方法的计算规则,或者是对图像亮度进行补偿,又或是对图像进行归一化计算。

(2)局部运动速度过大。两帧图像之间,如果全部或者部分区域的运动位移超过搜索范围,那么可以断定不可能找准匹配位置。扩大搜索范围还是改进搜索方法,也是未来可以考虑研究的方向。

(3)计算速度过慢。除了从算法上改进之外,还可以考虑引入并行计算,Nvidia的显卡所附带的CUDA并行计算功能,可以让传统算法速度有很大的提升,对于实时应用要求较高的读者不妨尝试一下。

(4)马赛克现象。由于人为将图像划分成为小子块,那么子块边缘的非连续性不可避免,或多或少会带来一定影响,如何解决这个问题也可进行研究。

(5)图像不连续的情况。对于两帧图像不连续,比如血管造影图像,当造影剂打入血管的一瞬间,前后两帧图像很可能大不相同,这就是不连续的情况。要从图形凭空生成或者从消失的图像中产生新的图像,是一个十分困难但又非常有意义的课题。

从医学影像的角度来说,如果能够获得估计精确的插入图像,不但可以让医生看着更加舒适,更重要的是能够让病人和医生接受更少的X射线照射。

对于其他领域,运动估计和运动补偿也发挥着重要作用。需要注意的是,读者在应用过程中应当结合自己的项目背景和成像原理进行计算策略的调整,以达到更好的成像效果。