基于超分辨率分形解码和误差补偿的图像插值算法

孙 冬，高清维，卢一相，竺 德

随着公共安全、军事遥感、医学成像、消费电子领域中图像技术的广泛应用，人们对成像分辨率提出越来越高的要求.图像采集系统受光电传感器有限的分辨率、透镜点扩散函数、光学衍射和噪声等因素的影响，往往只能获得降质的低分辨率图像，影响了后继图像分析与识别技术的进一步应用.从频域看，图像采集系统等效于一个低通滤波器，其传递函数在截止频率外的值为0.为了恢复丢失的高频信息，可以对低分辨率图像进行插值.这种基于软件的分辨率增强不需要对现有图像采集系统进行改造升级，不增加硬件成本，因此研究高精度的图像插值算法具有重要的理论意义和应用价值.

根据Shannon采样定理，对于一个带宽有限的图像，当采样频率不低于其Nyquist频率时，用sinc函数进行内插可以恢复原始的连续图像.由于sinc函数物理不可实现，一些学者使用最近邻、双线性和B样条［1－3］等一系列的经典算法进行近似.这类算法实现简单、速度快，但经常会在插值图像中形成模糊或锯齿效应，插值质量不高，其根本原因在于所采用的是一个过度简化的分段光滑图像模型.

基于局部自相似结构的图像描述方法是近年来提出的一种分形描述方法［4］，它将图像定义为一个带映射的局部迭代函数系统（local iterated f unction system with gray layer maps，简称LIFSM）的吸引子，并基于拼贴定理进行重建［5－8］.LIFSM通常由“父块－子块对”和相应的相似变换关系表示.因为块之间的相似关系与实际的数字图像分辨率无关，所以由相似关系所决定的吸引子能够以任意的精度重建.根据这一原理，作者对分形图像插值算法进行研究.

1 分形编码原理

记（M，d）为数字图像构成的测度空间，d为测度，对于定义在（M，d）上的压缩映射ω＝｛ω1，ω2，…，ωn｝及图像μ∈（M，d），若

则μ是由ω唯一决定的吸引子.对于任意给定的初始图像ζ∈（M，d），可利用下式对μ进行迭代求解

第i次迭代后的图像ωi（ζ）与μ之间的距离满足

其中：s＝max（s1，s2，…，sn），si为映射ωi的压缩因子.

在μ可分块的情况下，（1）式可写为

其中：Ri和Di均为μ中的图像块，构成μ的一个覆盖，且Ri∩Rj＝∅.一般情况下，Ri与Di分别可取为M×M及N×N像素的方块（N＞M）；ωi为收缩变换c、几何变换g、灰度变换f的组合：ωi＝fi◦gi◦c.收缩变换c将父块Di的尺寸从N×N像素缩小到与子块Ri相同的M×M像素，如图1所示.

几何变换ｇ 对块进行变形，并且限定为图２所示的８种特殊变换之一.

灰度变换f对父块的像素灰度值进行变换，通常使用线性变换f（x）＝a·x＋b，此时有

其中：参数ai和bi可通过下式计算

其中：cov、var和exp分别代表求协方差、方差和均值.

根据（5）式，（4）式可写为

（7）式给出了对图像进行分形编码的方法：首先将图像分割成互不重叠的子块｛Rii＝1，2，…，n｝，对于每个Ri，在图像中搜索一个与其最相似的父块Di，使得拼贴误差

取最小值.此时近似有

相应的参数Fi＝（ai，bi，gi，u，v）称为子块Ri 的分形码，其中（u，v）为 Di 的左上角坐标.当所有子块的分形码｛Fii＝1，2，…，n｝获取完毕时，图像编码结束.

对于分形图像解码，任意给定一幅与μ尺寸相等的初始图像ζ，根据（2）式进行L次迭代后，有

当L足够大时，μ可作为原始图像μ的近似.

μ与μ误差主要由（8）式给出的拼贴误差ε构成.为了减少误差，可采取块尺寸可变的四叉树分形编码策略：以一个较大的子块尺寸Mstart×Mstart为起点，为每个子块R搜索其最优父块D.如果R对应的拼贴误差ε大于某个给定的阈值，则将该子块分解为4个较小的块Rq、Ru、Ra和Rd（见图3）.图3a中的子块R在图3b中分解为4个较小的块Rq、Ru、Ra、Rd（以不同灰度标识），它们分别与相应的父块Dq、Du、Da、Dd存在变换关系.为每个块重新搜索最优父块并编码，该过程以递归方式进行，直到图像中所有的块都编码完毕.

2 插值算法

分形码给出了图像中具有相似性的块Ri与Di之间的变换关系.因为块之间的相似性与数字化图像的分辨率无关，所以当μ在行、列方向各自增强K倍分辨率（图像扩大K倍）后，Ri与Di仍是相似的，并且这种相似关系保持不变.

记μK为分辨率增强K 倍后的图像，由（4）式知

其中：l为插值算子，R′i和D′i为Ri和Di在高分辨率图像μK中的对应块.记Ri和Di的左上角坐标分别为（x，y）和（u，v），则R′i和D′i的左上角坐标分别为（Kx，Ky）和（Ku，Kv），尺寸分别为K M×K M和K N×K N.

与分形图像解码类似，任意给定一幅与μK尺寸相等的初始图像ζK，可利用下式对μK进行迭代求解

L次迭代后可得到关于图像μ的高分辨率估值＾μK.将（12）式给出的解码过程称为超分辨率分形解码.

因为自然图像不具有严格的局部自相似结构，所以分形图像编码是一种有损编码，其重构图像＾μ与原始图像μ之间存在误差.记

显然，＾μ为μ中具有严格局部自相似的分形部分，e为残余部分（误差图像）.对（13）式两端应用插值算子l，得

其中：eK为高分辨率图像μK与其估值＾μK的误差，使用下式对eK进行估算

为了减少插值误差，可将＾eK作为补偿项累加至＾μK中，并以此作为μK的一个更佳的估计.式中lbicubic为双立方插值算子.

综上所述，分形图像插值具体步骤如下（流程框图见图4）：

（1）构造图像的覆盖集｛Ri｝，其中每个子块Ri的尺寸均为Mstart×Mstart，且Ri∩Rj＝∅.

（2）将所有的子块都标记为“未编码”，并加入编码任务队列Qcode.

（3）从Qcode取出一个标记为“未编码”的子块R.

（4）记R的尺寸为M×M.以穷举的方式，在图像中所有尺寸为2 M×2 M的块中，寻找一个最优父块D，使其与R在（8）式的近似下的误差ε最小.

（5）若ε≤εquad或R的尺寸小于Mstop×Mstop：

（a）记录R的尺寸及分形码F；

（b）将R标记为“已编码”.

（6）否则：

（a）把R分解为4棵较小的子块Rq、Ru、Ra和Rd（如图3b所示），将它们标记为“未编码”并加入到Qcode；

（b）从Qcode中删除原有的子块R.

（7）返回第（3）步，对Qcode中下一棵子块进行编码，直到Qcode中所有的子块都标记为“已编码”.

（8）根据第（2）～（7）步得到的分形码｛Fi｝，对图像进行解码，并计算与原始图像μ的误差e.

（9）将Qcode中所有的子块R及相应父块D 的左上角坐标及尺寸分别乘以K 倍，得到新的覆盖集｛R′i｝以及覆盖集中每个子块R′i所对应的父块D′i.

（10）任选一幅尺寸K倍于μ的初始图像，利用 （12）式进行超分辨率解码，得到关于μ的高分辨率估值＾μK.

（11）对第（8）步所得误差图像e进行K 倍双立方插值，得到误差补偿项＾eK.

（12）将～μK＝＾μK＋＾eK作为μ的最终插值输出.

3 实 验

为验证作者所提出的分形插值算法的正确性和有效性，对一组标准测试图像进行插值实验，从主观视觉质量和PSNR值对插值的效果进行考查.所使用的测试图像全部取自标准测试图像库，图像尺寸为512×512像素，8位灰度量化.分形插值算法中的参数为：Mstart＝16，Mstop＝4，εquad＝100.

图5～6给出使用分形插值算法分别对Lena和Peppers进行2×2倍插值的结果.从帽子的边缘、眼睛的轮廓和辣椒的边缘等其他细节可以看出，由分形插值算法得到的图像具有清晰的边缘结构和纹理，不存在锯齿和模糊现象，主观视觉质量较高.

作为客观质量评价，表1对双线性、分形插值、边缘导向插值［9］（new edge direct interpolation，简称NEDI）和基于方向滤波器及数据融合插值［10］（directional filtering and coefficients f usion，简称DFCF）算法的插值结果进行PSNR对比，其中低分辨率图像的尺寸为256×256，由原始的标准测试图像作2×2局部平均后下采样得到.

表1 插值算法PSNR值对比Tab.1 Inter polation algorith ms comprised with PSNR

由表1可见，分形插值可获得比其他3种算法更高的PSNR值，因此插值精度更高，这与图像主观视觉质量的感受是一致的.

4 结束语

综上所述，基于分形的图像插值算法能够对结构细节实现准确有效的恢复，不会造成边缘模糊和锯齿效应，具有较高的插值精度和图像质量.下一步的研究目标应主要集中在进一步提高插值精度上，可以从两方面开展研究：一是改进图像分块策略，使用形状更自由的图像子块进行分形编码，减小拼贴误差；二是研究变换域中的分形图像插值算法.针对分形图像插值过长的计算时间，可以考虑利用GPU对编码进行加速，缩短编码时间.

［1］ Leh mann T M，G¨onner C，Spitzer K.Survey：interpolation methods in medical i mage processing［J］.IEEE Transactions on Medical Imaging，1999，18（11）：1049－1075.

［2］ Unser M.Splines：a perfect fit f or signal and i mage processing［J］.IEEE Transactions on Signal Pr ocessing Magazine，1999，16（6）：22－38.

［3］ Leh mann T M，Gonner C，Spitzer K.Addendu m：B－spline inter polation in medical i mage pr ocessing［J］.IEEE Transactions on Medical Imaging，2001，20（7）：660－665.

［4］ Bar nsley M F，Jacquin A.Application of recurrent iterated f unction systems to i mages［J］.Mathematics and Statistics，1989，5（1）：3－31.

［5］ Jacquin A.Image coding based on a fractal theor y of iterated contractive i mage transfor mations［J］.IEEE Transactions on Image Pr ocessing，1992，1（1）：18－30.

［6］ Davis G.A wavelet－based analysis of fractal i mage compression［J］.IEEE Transactions on Image Processing，1998，7（2）：141－154.

［7］ Hyun J，Lickho S.A spectrum－based searching technique f or the most favorable section of digital music［J］.IEEE Transactions on Consu mer Electr onics，2009，55（4）：2122－2126.

［8］ Lu J，Ye Z X，Zou Y R.Huber fractal i mage coding based on a fitting plane［J］.IEEE Transactions on Image Pr ocessing，2013，22（1）：134－145.

［9］ Li X，Orchard M T.New edge directed inter polation［J］.IEEE Transactions on Image Processing，2001，10（10）：1521－1527.

［10］ Zhang L.An edge guided i mage interpolation algorith m via directional filtering and coefficients f usion［J］.IEEE Transactions on Image Processing，2006，15（8）：2226－2238.