孫 冬，高清維，盧一相，竺 德

隨著公共安全、軍事遙感、醫(yī)學成像、消費電子領域中圖像技術的廣泛應用，人們對成像分辨率提出越來越高的要求.圖像采集系統(tǒng)受光電傳感器有限的分辨率、透鏡點擴散函數(shù)、光學衍射和噪聲等因素的影響，往往只能獲得降質(zhì)的低分辨率圖像，影響了后繼圖像分析與識別技術的進一步應用.從頻域看，圖像采集系統(tǒng)等效于一個低通濾波器，其傳遞函數(shù)在截止頻率外的值為0.為了恢復丟失的高頻信息，可以對低分辨率圖像進行插值.這種基于軟件的分辨率增強不需要對現(xiàn)有圖像采集系統(tǒng)進行改造升級，不增加硬件成本，因此研究高精度的圖像插值算法具有重要的理論意義和應用價值.

根據(jù)Shannon采樣定理，對于一個帶寬有限的圖像，當采樣頻率不低于其Nyquist頻率時，用sinc函數(shù)進行內(nèi)插可以恢復原始的連續(xù)圖像.由于sinc函數(shù)物理不可實現(xiàn)，一些學者使用最近鄰、雙線性和B樣條［1－3］等一系列的經(jīng)典算法進行近似.這類算法實現(xiàn)簡單、速度快，但經(jīng)常會在插值圖像中形成模糊或鋸齒效應，插值質(zhì)量不高，其根本原因在于所采用的是一個過度簡化的分段光滑圖像模型.

基于局部自相似結(jié)構(gòu)的圖像描述方法是近年來提出的一種分形描述方法［4］，它將圖像定義為一個帶映射的局部迭代函數(shù)系統(tǒng)（local iterated f unction system with gray layer maps，簡稱LIFSM）的吸引子，并基于拼貼定理進行重建［5－8］.LIFSM通常由“父塊－子塊對”和相應的相似變換關系表示.因為塊之間的相似關系與實際的數(shù)字圖像分辨率無關，所以由相似關系所決定的吸引子能夠以任意的精度重建.根據(jù)這一原理，作者對分形圖像插值算法進行研究.

1 分形編碼原理

記（M，d）為數(shù)字圖像構(gòu)成的測度空間，d為測度，對于定義在（M，d）上的壓縮映射ω＝｛ω1，ω2，…，ωn｝及圖像μ∈（M，d），若

則μ是由ω唯一決定的吸引子.對于任意給定的初始圖像ζ∈（M，d），可利用下式對μ進行迭代求解

第i次迭代后的圖像ωi（ζ）與μ之間的距離滿足

其中：s＝max（s1，s2，…，sn），si為映射ωi的壓縮因子.

在μ可分塊的情況下，（1）式可寫為

其中：Ri和Di均為μ中的圖像塊，構(gòu)成μ的一個覆蓋，且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，此時有

其中：參數(shù)ai和bi可通過下式計算

其中：cov、var和exp分別代表求協(xié)方差、方差和均值.

根據(jù)（5）式，（4）式可寫為

（7）式給出了對圖像進行分形編碼的方法：首先將圖像分割成互不重疊的子塊｛Rii＝1，2，…，n｝，對于每個Ri，在圖像中搜索一個與其最相似的父塊Di，使得拼貼誤差

取最小值.此時近似有

相應的參數(shù)Fi＝（ai，bi，gi，u，v）稱為子塊Ri 的分形碼，其中（u，v）為 Di 的左上角坐標.當所有子塊的分形碼｛Fii＝1，2，…，n｝獲取完畢時，圖像編碼結(jié)束.

對于分形圖像解碼，任意給定一幅與μ尺寸相等的初始圖像ζ，根據(jù)（2）式進行L次迭代后，有

當L足夠大時，μ可作為原始圖像μ的近似.

μ與μ誤差主要由（8）式給出的拼貼誤差ε構(gòu)成.為了減少誤差，可采取塊尺寸可變的四叉樹分形編碼策略：以一個較大的子塊尺寸Mstart×Mstart為起點，為每個子塊R搜索其最優(yōu)父塊D.如果R對應的拼貼誤差ε大于某個給定的閾值，則將該子塊分解為4個較小的塊Rq、Ru、Ra和Rd（見圖3）.圖3a中的子塊R在圖3b中分解為4個較小的塊Rq、Ru、Ra、Rd（以不同灰度標識），它們分別與相應的父塊Dq、Du、Da、Dd存在變換關系.為每個塊重新搜索最優(yōu)父塊并編碼，該過程以遞歸方式進行，直到圖像中所有的塊都編碼完畢.

2 插值算法

分形碼給出了圖像中具有相似性的塊Ri與Di之間的變換關系.因為塊之間的相似性與數(shù)字化圖像的分辨率無關，所以當μ在行、列方向各自增強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）式給出的解碼過程稱為超分辨率分形解碼.

因為自然圖像不具有嚴格的局部自相似結(jié)構(gòu)，所以分形圖像編碼是一種有損編碼，其重構(gòu)圖像＾μ與原始圖像μ之間存在誤差.記

顯然，＾μ為μ中具有嚴格局部自相似的分形部分，e為殘余部分（誤差圖像）.對（13）式兩端應用插值算子l，得

其中：eK為高分辨率圖像μK與其估值＾μK的誤差，使用下式對eK進行估算

為了減少插值誤差，可將＾eK作為補償項累加至＾μK中，并以此作為μK的一個更佳的估計.式中l(wèi)bicubic為雙立方插值算子.

綜上所述，分形圖像插值具體步驟如下（流程框圖見圖4）：

（1）構(gòu)造圖像的覆蓋集｛Ri｝，其中每個子塊Ri的尺寸均為Mstart×Mstart，且Ri∩Rj＝?.

（2）將所有的子塊都標記為“未編碼”，并加入編碼任務隊列Qcode.

（3）從Qcode取出一個標記為“未編碼”的子塊R.

（4）記R的尺寸為M×M.以窮舉的方式，在圖像中所有尺寸為2 M×2 M的塊中，尋找一個最優(yōu)父塊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）根據(jù)第（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 實 驗

為驗證作者所提出的分形插值算法的正確性和有效性，對一組標準測試圖像進行插值實驗，從主觀視覺質(zhì)量和PSNR值對插值的效果進行考查.所使用的測試圖像全部取自標準測試圖像庫，圖像尺寸為512×512像素，8位灰度量化.分形插值算法中的參數(shù)為：Mstart＝16，Mstop＝4，εquad＝100.

圖5～6給出使用分形插值算法分別對Lena和Peppers進行2×2倍插值的結(jié)果.從帽子的邊緣、眼睛的輪廓和辣椒的邊緣等其他細節(jié)可以看出，由分形插值算法得到的圖像具有清晰的邊緣結(jié)構(gòu)和紋理，不存在鋸齒和模糊現(xiàn)象，主觀視覺質(zhì)量較高.

作為客觀質(zhì)量評價，表1對雙線性、分形插值、邊緣導向插值［9］（new edge direct interpolation，簡稱NEDI）和基于方向濾波器及數(shù)據(jù)融合插值［10］（directional filtering and coefficients f usion，簡稱DFCF）算法的插值結(jié)果進行PSNR對比，其中低分辨率圖像的尺寸為256×256，由原始的標準測試圖像作2×2局部平均后下采樣得到.

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

由表1可見，分形插值可獲得比其他3種算法更高的PSNR值，因此插值精度更高，這與圖像主觀視覺質(zhì)量的感受是一致的.

4 結(jié)束語

綜上所述，基于分形的圖像插值算法能夠?qū)Y(jié)構(gòu)細節(jié)實現(xiàn)準確有效的恢復，不會造成邊緣模糊和鋸齒效應，具有較高的插值精度和圖像質(zhì)量.下一步的研究目標應主要集中在進一步提高插值精度上，可以從兩方面開展研究：一是改進圖像分塊策略，使用形狀更自由的圖像子塊進行分形編碼，減小拼貼誤差；二是研究變換域中的分形圖像插值算法.針對分形圖像插值過長的計算時間，可以考慮利用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.